基于多目標優(yōu)化NSGA2改進算法的結構動力學模型確認

賴文星， 鄧忠民， 張鑫杰

（北京航空航天大學 宇航學院，北京100191）

1 引 言

多目標進化算法從20世紀90年代開始迅速發(fā)展，Deb等［1］提出第二代帶精英保留策略的快速非支配排序算法NSGA2。NSGA2采用快速非支配排序方法，基于擁擠距離的分布性方法和精英保留策略，憑借簡單及高效等優(yōu)點，廣泛應用于科學計算和工程設計等領域。Kollat等［2］將Epsilon支配概念引入 NSGA2，提出Epsilon-NSGA2算法；Zhang等［3］提出了基于分解的多目標進化算法MOEA／D，MOEA／D將多個目標分為若干組，再并行優(yōu)化求解；Elhossini等［4］提出粒子群算法和進化算法的混合算法；Deb等［5］提出一種基于參考點的NSGA2算法，以提高高維優(yōu)化能力；Shim等［6］將非支配排序與目標分解結合，以提高算法優(yōu)化性能；Qiu等［7］提出用于多目標優(yōu)化的自適應交叉差分演化算子。以上改進一定程度上改善了NSGA2的優(yōu)化性能，但NSGA2還存在較多的設計缺陷，如計算復雜度較高、無法識別偽非支配解和擁擠距離公式不合理等。針對上述不足，本文提出一種基于支配強度的NSGA2改進算法，在非支配排序、分布性保持和精英保留策略三個方面對NSGA2進行改進，致力于提高該算法的計算效率、收斂性和分布性。

結構動力學模型確認是工程結構設計的主要組成部分之一［8］。2002年國際振動模態(tài)分析會議上，美國Los Alamos國家實驗室首次提出關于有限元模型不確定性描述的討論，并且率先開展了此項研究。Calvi［9］從靈敏度方法研究、有限元建模和貝葉斯估計方法等方面展開了研究。文獻［10］利用模型確認方法對氣動不確定建模問題進行了驗證，證明模型確認方法的可行性。文獻［11］提出基于蒙特卡洛法的結構動力學模型確認方法，將蒙特卡洛法與多元回歸分析相結合，并分析了不同距離準則下的確認精度。但是，傳統(tǒng)模型確認方法普遍采用單目標描述多個動力學特性與真實值之間的差異，并未考慮設計變量的干擾和忽視魯棒性的影響［12］，對于復雜結構難以得到理想的確認結果。本文提出一種考慮魯棒性的結構動力學模型確認方法，建立以馬氏距離和魯棒性作為優(yōu)化目標的數(shù)學模型，并將INSGA2-DS改進算法應用于求解GARTEUR飛機模型算例。

2 基于支配強度的NSGA2改進算法

2.1 快速支配強度排序法

NSGA2采用快速非支配排序法，而傳統(tǒng)的支配關系容易造成非支配集NDSet中存在大量偽非支配解，降低收斂效率和解集質量。為去除NDSet中的偽非支配解，INSGA2-DS引入支配強度：

式中 M為子目標數(shù)量，fmaxk和fmink分別代表NDSet中第k個子目標的最大值和最小值。

通過引入支配強度去除偽非支配解，保證性能優(yōu)異的個體優(yōu)先獲得遺傳，提高解集質量。

2.2 分布性保持方法

分布性是指與真實Pareto邊界的相似程度和分布的廣泛程度，反映解集分布的均勻性和多樣性。NSGA2僅考慮個體與相鄰個體之間的距離，在不同子目標上擁擠距離差異程度較大的個體無法獲得更多的遺傳機會，不利于保持分布性。本文提出一種考慮方差的新型擁擠距離公式為

經(jīng)過上述計算后，種群個體都具有非支配集等級rank、支配強度ξ和擁擠距離iD三個參數(shù)。在INSGA2-DS中，先比較rank，其次比較ξ，最后比較iD。個體i優(yōu)于個體j的比較算子如下，

ranki＜rankj或ranki＝rankj∩ξi＜ξj

或ranki＝rankj∩ξi＝ξj∩iD＞jD

2.3 自適應精英保留策略

NSGA2的精英保留規(guī)模是一個固定值，不利于解集收斂。INSGA2-DS采用自適應精英保留策略，進化前期，限制精英保留規(guī)模，增加解集多樣性；進化后期，增加精英保留規(guī)模，提高算法收斂性。第k代種群精英保留規(guī)模Elitistk為

式中αk為第k代精英保留規(guī)模的影響因子，Pop為種群規(guī)模。αk的自適應迭代公式為

式中αk＋1為第k＋1代精英保留規(guī)模的影響因子，ρ為第k代種群的NDSet與Pop 的比值。

為了防止精英保留規(guī)模走向極端，本文設定，當αk＞0.8時，則αk＝0.8；當αk＜0.2時，則αk＝0.2。

2.4 進化操作

INSGA2-DS采用模擬二進制交叉和多項式變異產(chǎn)生子代。若進入交叉操作，從父輩中隨機選取兩個個體P1和P2，子代個體Q1和Q2計算如下，

式中β與隨機數(shù)μ∈［0，1］有關，

式中ηc為交叉分布指數(shù)。

若進入變異操作，從父輩種群中隨機選取一個個體P3，產(chǎn)生一個子代個體Q3，

式中 隨機數(shù)r∈［0，1］，ηm為變異分布指數(shù)。

2.5 算法流程

INSGA2-DS的算法流程具體如下。

（1）隨機初始化種群。

（2）利用快速支配強度排序法對種群中全部個體進行分層排序。

（3）根據(jù)種群個體的非支配等級、支配強度和擁擠距離，通過自適應精英保留策略，選擇精英個體直接進入下一代種群。

（4）采用模擬二進制交叉和多項式變異方法，由父代種群產(chǎn)生子代種群。

（5）將父代種群中精英個體與子代種群混合組成下一代種群。

（6）若達到進化終止條件，則停止循環(huán)，輸出結果；否則轉至步驟（2）。

3 結構動力學模型確認方法

結構動力學模型確認通常采用蒙特卡羅模擬，對有限元模型進行隨機性分析［13］。蒙特卡羅模擬存在計算量過大的缺點，而基于代理模型的蒙特卡羅模擬，可以極大降低計算時間［14］。首先，本文對結構有限元模型進行隨機性分析，獲得一組不確定性參數(shù)輸入和動力學響應輸出；然后，采用神經(jīng)網(wǎng)絡作為代理模型，得到模型不確定性參數(shù)與動力學響應的映射關系，再建立考慮魯棒性的多目標優(yōu)化模型；最后，采用INSGA2-DS算法進行求解。

3.1 神經(jīng)網(wǎng)絡代理模型

首先由輸入樣本X＝［x1，x2，…，xn］T和輸出樣本Y＝［y1，y2，…，ym］T建立神經(jīng)元，

式中 wij為神經(jīng)元j到神經(jīng)元i的連接權值，θ為閾值，neti為神經(jīng)元i的凈激活，f為激活函數(shù)。采用雙極S形激活函數(shù)，

采用Delta有導師學習算法，根據(jù)實際輸出與期望輸出差別來調整連接權，

式中α為學習速度，yi為神經(jīng)元的實際輸出，di為神經(jīng)元i的期望輸出。

3.2 考慮魯棒性的多目標優(yōu)化模型

選取固有頻率作為動力學響應的對比參數(shù)。假設某結構前k階真實固有頻率為Freal＝［f＊1，f＊2，…，f＊k］T，模態(tài)分析所得固有頻率為F＝［f1，f2，…，fk］T。為評價F 與Freal的真實差異，選取馬氏距離作為評判指標，

式中 Σ－1為協(xié)方差矩陣的逆矩陣。

為計算各階固有頻率的偏離程度，選取相對誤差的方差作為模型確認的魯棒性指標。F與Freal的平均相對誤差公式為

式中ωi為第i階固有頻率的權重系數(shù)。固有頻率相對誤差的方差計算公式為

因此，考慮魯棒性的結構動力學模型確認方法的多目標優(yōu)化模型為

式中 不確定參數(shù)X＝［x1，x2，…，xn］T，和分別為第i個不確定參數(shù)的最小值和最大值。X與固有頻率F的映射關系由神經(jīng)網(wǎng)絡代理模型提供。

4 GARTEUR飛機模型算例

GARTEUR飛機模型是由法國航空航天研究院1995年設計制造的標準飛機模型，用來驗證振動試驗的分析技術［15］。本文用PATRAN建立GARTEUR的有限元模型，如圖1所示。該模型由1704個節(jié)點及846個六面體單元組成，材料密度D ＝2.7g／cm3，彈性模量E ＝70GPa，泊松比μ＝0.3。

對GARTEUR的D和E進行確認，選取第1～10階固有頻率的測量值作為模型確認的真實值。用MATLAB軟件產(chǎn)生56組D和E的隨機數(shù)作為輸入樣本，D和E均服從正態(tài)分布，均值為μD＝3.0g／cm3和μE＝7.5243×1010Pa，變異系數(shù)θE＝θD＝0.1。用NASTRAN軟件進行模態(tài)分析，取前10階固有頻率作為輸出樣本。通過神經(jīng)網(wǎng)絡代理模型，建立D和E與前10階固有頻率f1，f2，…，f10之間的映射關系，并建立考慮魯棒性的結構動力學模型確認的多目標優(yōu)化模型。

選取INSGA2-DS、NSGA2和MOPSO共3種算法，分別對上述多目標優(yōu)化模型進行優(yōu)化。算法參數(shù)設置為，種群大小50，進化代數(shù)100，交叉概率0.9，變異概率0.1，交叉參數(shù)和變異分布指數(shù)都為20。神經(jīng)網(wǎng)絡參數(shù)設置為，訓練網(wǎng)絡隱含層10，訓練次數(shù)1000，全局最小誤差0.0001，學習速率0.025。

從圖2可以看出，INSGA2-DS計算所得的Pareto邊界的收斂性和分布性均優(yōu)于NSGA2，可以獲得更多馬氏距離更小、魯棒性更好的解集；MOPSO的解集與INSGA2-DS的部分解集相近，但是很明顯MOPSO陷入局部最優(yōu)解集，沒有獲得更多魯棒性更好的解集。

為具體比較3種算法的優(yōu)化結果，選取每種算法Pareto邊界最左上端和最右下端的個體，即馬氏距離最小個體和魯棒性最好的個體，分別記為目標1和目標2。具體結果統(tǒng)計列入表1，表中D，E和固有頻率f的單位分別為g／cm3，GPa和Hz。

圖1 GARTEUR飛機有限元模型Fig.1 Finite element model of GARTEUR aircraft

圖2 GARTEUR模型確認優(yōu)化結果比較Fig.2 GARTEUR model validation optimization results comparison

根據(jù)表1數(shù)據(jù)，計算3種算法模型確認結果與真實值的相對誤差。馬氏距離最?。繕?）和魯棒性最好（目標2）兩種情況下的相對誤差如圖3和圖4所示。

從圖3和圖4可以看出，INSGA2-DS計算結果除第2階固有頻率f2和第7階固有頻率f7相對誤差超過2%，其他參數(shù)的相對誤差均低于2%，并且總體上優(yōu)于NSGA2和MOPSO的計算結果。馬氏距離最小時，INSGA2-DS，NSGA2和 MOPSO的相對誤差的平均數(shù)分別為1.28%，3.70%和2.06%；魯棒性最好時，INSGA2-DS，NSGA2 和MOPSO的相對誤差的平均數(shù)分別為1.01%，3.40%和2.19%。NSGA2計算結果比INSGA2-DS差，但除了f3和f9，其他參數(shù)相對誤差均小于5%；MOPSO 的計算結果介于INSGA2-DS和NSGA2之間?？梢?種算法中，不論馬氏距離還是魯棒性，INSGA2-DS的相對誤差均最小，優(yōu)化性能最佳。

表1 GARTEUR飛機模型確認結果Tab.1 GARTEUR aircraft model validation result

圖3 馬氏距離最小的相對誤差比較Fig.3 Relative error comparison on minimal Mahalanobis distance

圖4 魯棒性最好的相對誤差比較Fig.4 Relative error comparison on the best robustness

為進一步比較3種算法的收斂性能，對GARTEUR模型進行了500次優(yōu)化分析，馬氏距離最小值和魯棒性最小值的迭代結果如圖5和圖6所示。

從圖5可以看出，在目標1馬氏距離上，INSGA2-DS收斂最快，收斂至4.9附近；MOPSO收斂速度次之，NSGA2最慢，但都收斂至5.0左右?？傮w上，3種算法在馬氏距離收斂性能上大致相當。

從圖6可以看出，在目標2魯棒性上，3種算法收斂速度差異不大，但最終INSGA2-DS，NSGA2和 MOPSO分別收斂至1.8×10－4，2.4×10－4和2.98×10－4；MOPSO陷入局部最優(yōu)化，不確定參數(shù)的總體魯棒性較差。與NSGA2和MOPSO相比，INSGA2-DS在馬氏距離和魯棒性上都具有更好的收斂性能。

圖5 馬氏距離迭代曲線Fig.5 Iteration curve of Mahalanobis distance

圖6 魯棒性迭代曲線Fig.6 Iteration curve of robustness

5 結 論

本文提出了一種基于支配強度的NSGA2改進算法，建立以馬氏距離和魯棒性作為優(yōu)化目標的結構動力學模型確認方法，并將INSGA2-DS改進算法用于求解該多目標優(yōu)化問題。研究結果表明，

（1）INSGA2-DS引入非支配強度，去除偽非支配解；改進擁擠距離計算公式，增加解集的多樣性；采用自適應精英保留策略，提高算法收斂性能。改進算法有效地改善了NSGA2的設計缺陷，提高其求解復雜和非線性多目標問題的優(yōu)化能力。

（2）神經(jīng)網(wǎng)絡代理模型通過對輸入輸出樣本的擬合逼近，可建立不確定參數(shù)與動力學響應之間的映射關系，大大減少蒙特卡羅模擬的計算量，具有很強的非線性映射能力、容錯能力和魯棒性優(yōu)勢。

（3）與傳統(tǒng)模型確認方法采用單目標優(yōu)化不同，考慮魯棒性的模型確認方法可獲得同時滿足多個不確定性量化指標的解集。GARTEUR算例中，INSGA2-DS解集的收斂性和分布性均優(yōu)于NSGA2，多樣性明顯優(yōu)于MOPSO，其收斂速度最快且解集質量最高，模型確認結果總體誤差小于2%，滿足實際工程需求。

考慮魯棒性的結構動力學模型確認方法可以獲得更多滿足不同設計目標的解集，為實際工程應用提供一種有益的參考。此外，本文只考慮了待確定參數(shù)，并未考慮其他參數(shù)可能存在微小擾動對模型確認精度的影響，還有待對魯棒性進行進一步研究。