賴(lài)文星， 鄧忠民， 張?chǎng)谓?/p>

（北京航空航天大學(xué) 宇航學(xué)院，北京100191）

1 引 言

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

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

2 基于支配強(qiáng)度的NSGA2改進(jìn)算法

2.1 快速支配強(qiáng)度排序法

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

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

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

2.2 分布性保持方法

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

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

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

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

2.3 自適應(yīng)精英保留策略

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

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

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

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

2.4 進(jìn)化操作

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

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

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

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

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

2.5 算法流程

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

（1）隨機(jī)初始化種群。

（2）利用快速支配強(qiáng)度排序法對(duì)種群中全部個(gè)體進(jìn)行分層排序。

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

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

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

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

3 結(jié)構(gòu)動(dòng)力學(xué)模型確認(rèn)方法

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

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

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

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

采用Delta有導(dǎo)師學(xué)習(xí)算法，根據(jù)實(shí)際輸出與期望輸出差別來(lái)調(diào)整連接權(quán)，

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

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

選取固有頻率作為動(dòng)力學(xué)響應(yīng)的對(duì)比參數(shù)。假設(shè)某結(jié)構(gòu)前k階真實(shí)固有頻率為Freal＝［f＊1，f＊2，…，f＊k］T，模態(tài)分析所得固有頻率為F＝［f1，f2，…，fk］T。為評(píng)價(jià)F 與Freal的真實(shí)差異，選取馬氏距離作為評(píng)判指標(biāo)，

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

為計(jì)算各階固有頻率的偏離程度，選取相對(duì)誤差的方差作為模型確認(rèn)的魯棒性指標(biāo)。F與Freal的平均相對(duì)誤差公式為

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

因此，考慮魯棒性的結(jié)構(gòu)動(dòng)力學(xué)模型確認(rèn)方法的多目標(biāo)優(yōu)化模型為

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

4 GARTEUR飛機(jī)模型算例

GARTEUR飛機(jī)模型是由法國(guó)航空航天研究院1995年設(shè)計(jì)制造的標(biāo)準(zhǔn)飛機(jī)模型，用來(lái)驗(yàn)證振動(dòng)試驗(yàn)的分析技術(shù)［15］。本文用PATRAN建立GARTEUR的有限元模型，如圖1所示。該模型由1704個(gè)節(jié)點(diǎn)及846個(gè)六面體單元組成，材料密度D ＝2.7g／cm3，彈性模量E ＝70GPa，泊松比μ＝0.3。

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

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

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

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

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

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

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

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

表1 GARTEUR飛機(jī)模型確認(rèn)結(jié)果Tab.1 GARTEUR aircraft model validation result

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

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

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

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

從圖6可以看出，在目標(biāo)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 結(jié) 論

本文提出了一種基于支配強(qiáng)度的NSGA2改進(jìn)算法，建立以馬氏距離和魯棒性作為優(yōu)化目標(biāo)的結(jié)構(gòu)動(dòng)力學(xué)模型確認(rèn)方法，并將INSGA2-DS改進(jìn)算法用于求解該多目標(biāo)優(yōu)化問(wèn)題。研究結(jié)果表明，

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

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

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

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