吳貝尼，夏利娟

（1.上海交通大學(xué)海洋工程國(guó)家重點(diǎn)實(shí)驗(yàn)室，上海200240；2.高新船舶與深海開發(fā)裝備協(xié)同創(chuàng)新中心，上海200240）

0 引 言

結(jié)構(gòu)拓?fù)鋬?yōu)化是在給定的設(shè)計(jì)空間里確定滿足給定約束的最優(yōu)材料分布以求獲得最佳結(jié)構(gòu)性能。在過去幾十年中，各種優(yōu)化方法已被成功運(yùn)用到各類工程結(jié)構(gòu)拓?fù)鋬?yōu)化設(shè)計(jì)中，包括均勻化方法[1]、變密度法[2-4]、雙向漸進(jìn)結(jié)構(gòu)優(yōu)化方法（BESO）[5-7]、獨(dú)立連續(xù)映射方法[8]和水平集方法[9]等。其中，BESO方法原理簡(jiǎn)單，收斂性好，通用性強(qiáng)，易于編程實(shí)現(xiàn)，因此備受青睞。

BESO 方法分為硬刪除（hard-kill）方法[10]和軟刪除（soft-kill）方法[11-12]，前者直接刪除低效單元，因此計(jì)算效率高，而后者引入帶懲罰項(xiàng)人工材料插值模型，并非完全移除單元，而是賦予其非常低的材料屬性。兩種BESO 方法（本文中稱為標(biāo)準(zhǔn)BESO 方法）在優(yōu)化求解的過程中同時(shí)考慮刪去低效單元和增添高效單元，有很大幾率獲得最優(yōu)結(jié)構(gòu)，但是二者增刪單元是依據(jù)某種優(yōu)化準(zhǔn)則的，每次迭代優(yōu)化受控于進(jìn)化率，從而可能導(dǎo)致誤刪最優(yōu)拓?fù)浣Y(jié)構(gòu)中的高效單元而無法恢復(fù)，最終得到的只是局部最優(yōu)解。

遺傳算法基于適者生存思想，有很強(qiáng)的全局尋優(yōu)能力，而且遺傳算法中優(yōu)勝劣汰的思想十分契合結(jié)構(gòu)拓?fù)鋬?yōu)化中高效單元保留、低效單元淘汰的原則。因此，將遺傳算法結(jié)合到BESO方法中，能有效克服BESO方法的缺陷，進(jìn)一步提高優(yōu)化求解質(zhì)量。

基于此，本文提出了一種基于改進(jìn)遺傳算法的軟刪除雙向漸進(jìn)結(jié)構(gòu)優(yōu)化法（G-BESO），采用一種考慮單元密度歷史信息的材料插值模型以增強(qiáng)恢復(fù)誤刪高效單元的能力，利用遺傳算法中的交叉和變異操作啟發(fā)式更新結(jié)構(gòu)狀態(tài)以提高全局尋優(yōu)能力。結(jié)合有限元平衡方程建立了連續(xù)體結(jié)構(gòu)優(yōu)化模型并給出具體求解流程，最后用經(jīng)典算例驗(yàn)證了方法的有效性和穩(wěn)健性，同時(shí)對(duì)比標(biāo)準(zhǔn)BESO 方法表明了其優(yōu)越性。

1 基于改進(jìn)遺傳算法的雙向漸進(jìn)結(jié)構(gòu)優(yōu)化模型構(gòu)建

1.1 算法思想

首先利用有限單元將連續(xù)體結(jié)構(gòu)離散化，設(shè)計(jì)區(qū)域?yàn)榉N群，離散單元為個(gè)體，以離散的單元密度值作為設(shè)計(jì)變量（0 或1），其中0 代表孔洞單元，1 表示實(shí)體單元；對(duì)于個(gè)體（單元）使用二進(jìn)制編碼，編碼中1的數(shù)目代表著其重要程度即適應(yīng)度。適應(yīng)度決定個(gè)體編碼的更新方式，個(gè)體之間通過交叉、變異更新其編碼，編碼決定個(gè)體保留與淘汰。通過這種方式，實(shí)現(xiàn)種群的進(jìn)化即結(jié)構(gòu)的優(yōu)化。

標(biāo)準(zhǔn)BESO方法中每一步迭代優(yōu)化通常需要?jiǎng)h除由進(jìn)化率決定的固定體積的單元，完全根據(jù)靈敏度大小機(jī)械地增添單元。本文提出的G-BESO 方法基于優(yōu)勝劣汰原則，根據(jù)由適應(yīng)度決定的單元編碼中1的數(shù)量啟發(fā)式增添不定體積的單元。具體而言，若單元為實(shí)單元，且其編碼值均為0，則從結(jié)構(gòu)中刪除；若為空單元，且其編碼中1的數(shù)量達(dá)到一定值，則添加到結(jié)構(gòu)中。

1.2 數(shù)學(xué)模型

本文研究連續(xù)體結(jié)構(gòu)在給定體積約束下的最大剛度設(shè)計(jì)問題，拓?fù)鋬?yōu)化模型可以描述為

確定設(shè)計(jì)變量：

式中：F、u分別是載荷向量和位移向量；K是剛度矩陣；V*是給定的體積約束；xi是設(shè)計(jì)變量，其值為0或1，表征空/實(shí)單元。

本文提出的G-BESO 方法是一種軟刪除方法，采用了一種新的考慮單元密度歷史信息的材料插值模型。引入文獻(xiàn)[13]中單元權(quán)重系數(shù)的概念，結(jié)合SIMP 材料插值模型，建立了單元權(quán)重系數(shù)I與單元密度X、單元材料屬性E的關(guān)系式：

1.3 靈敏度計(jì)算與收斂準(zhǔn)則

單元靈敏度值代表其對(duì)整個(gè)結(jié)構(gòu)剛度貢獻(xiàn)的大小，關(guān)系著單元的刪除與保留，決定最終結(jié)構(gòu)拓?fù)湫螤?。?duì)于剛度最大優(yōu)化問題，基于所提出的材料插值模型的單元應(yīng)變能靈敏度為

式中，ui，Ki，Ii分別是單元i的位移向量、剛度矩陣和權(quán)重系數(shù)。

為解決棋盤格問題和網(wǎng)格依賴性，本文采用文獻(xiàn)[10]中的靈敏度過濾方法：

式中，N為過濾范圍內(nèi)節(jié)點(diǎn)總數(shù)，rmin為過濾半徑，rij為單元i中心與節(jié)點(diǎn)j之間的距離。

同時(shí)，為穩(wěn)定優(yōu)化過程，當(dāng)前單元靈敏度需要根據(jù)上一步所得值進(jìn)行修正[11]：

結(jié)構(gòu)拓?fù)鋬?yōu)化是一個(gè)不斷更新單元狀態(tài)的迭代過程，直到滿足體積約束和收斂準(zhǔn)則，迭代終止。本文采用文獻(xiàn)[11]中推薦的收斂終止公式：

式中：τ為收斂因子即迭代精度，本文取0.001；k為迭代步；N為常數(shù)，一般取為5。表征結(jié)構(gòu)的應(yīng)變能在連續(xù)10次以上的迭代中變化量足夠小即認(rèn)為得到穩(wěn)定結(jié)構(gòu)。

2 基于改進(jìn)遺傳算法的雙向漸進(jìn)結(jié)構(gòu)優(yōu)化求解流程

2.1 技術(shù)手段

為提高G-BESO 方法的實(shí)用性，使之能用于優(yōu)化復(fù)雜的工程結(jié)構(gòu)，本文以Python 為編程手段，通過調(diào)用HyperMesh軟件更新有限元模型，調(diào)用Nastran求解器進(jìn)行有限元計(jì)算，可實(shí)現(xiàn)對(duì)實(shí)際工程結(jié)構(gòu)的拓?fù)鋬?yōu)化設(shè)計(jì)。

2.2 求解步驟

基于上一章構(gòu)建的拓?fù)鋬?yōu)化模型，給出G-BESO方法的具體求解步驟如下：

（1）使用HyperMesh對(duì)連續(xù)體結(jié)構(gòu)進(jìn)行網(wǎng)格劃分，定義設(shè)計(jì)區(qū)域，施加邊界條件和載荷。

（2）利用Nastran進(jìn)行結(jié)構(gòu)有限元分析，計(jì)算過濾后的單元靈敏度。

（3）根據(jù)靈敏度將單元?jiǎng)澐譃槿齻€(gè)集合：高效單元集、過渡單元集和低效單元集。若為第一次進(jìn)化，對(duì)設(shè)計(jì)區(qū)域單元進(jìn)行編碼：對(duì)所有單元賦予少量0、較多1 的編碼。計(jì)算進(jìn)度指標(biāo)Prg，交叉概率Pc，變異概率Pm。

式中，Vcur，V*分別為結(jié)構(gòu)當(dāng)前體積分?jǐn)?shù)和約束體積分?jǐn)?shù)。

為解決遺傳算法隨機(jī)性引發(fā)的跳躍解、非可行解（結(jié)構(gòu)不協(xié)調(diào)）等問題，參照文獻(xiàn)[14]，引入懲罰因子如下，使得Pc與Pm隨結(jié)構(gòu)優(yōu)化逐漸增大：

式中，Penc，Penm分別是預(yù)先選定的交叉和變異懲罰因子。

（4）個(gè)體（單元）之間實(shí)施交叉、變異以更新單元編碼。

類似于精英保留策略，高效單元不作為父代（但作為母代）參與交叉直接進(jìn)入下一代，以避免高效單元意外喪失導(dǎo)致優(yōu)化陷入局部最優(yōu)解。其他單元作為父代隨機(jī)選擇單元（母代）進(jìn)行一次配對(duì)交叉，產(chǎn)生的子代之一進(jìn)入下一代，方式為兩點(diǎn)交叉。選擇同一集合中的單元配對(duì)交叉的概率為Pc，選擇其他集合中單元的概率為( 1 -Pc)/2。

變異操作只作用于高效單元集和低效單元集，直接改變單元二進(jìn)制編碼中的數(shù)位值[14]。對(duì)于每個(gè)數(shù)位來說，其變異的概率設(shè)置為Pm，如果產(chǎn)生的隨機(jī)數(shù)值小于Pm，高效單元編碼中0 變?yōu)?，低效單元中的1變?yōu)?。這樣，變異可以使高效單元編碼中含有更多的1，而低效單元編碼中含有更多的0，結(jié)構(gòu)逐步優(yōu)化。

（5）根據(jù)單元編碼中1的數(shù)量確定單元狀態(tài)，更新設(shè)計(jì)變量，調(diào)用HyperMesh修改結(jié)構(gòu)。

盡管遺傳算法能以較大概率找到全局最優(yōu)解，由于其固有的隨機(jī)性，仍可能限于局部最優(yōu)。為控制每一次進(jìn)化結(jié)果，每一步交叉變異更新單元編碼的過程進(jìn)行若干次（10次），根據(jù)應(yīng)力離散度[15]進(jìn)行評(píng)估，從中預(yù)選出最優(yōu)解，更新設(shè)計(jì)變量。需要指出的是，此處不需要進(jìn)行有限元計(jì)算，只需進(jìn)行的簡(jiǎn)單代數(shù)運(yùn)算并不耗時(shí)。

應(yīng)力離散度表征結(jié)構(gòu)受力的均勻程度。隨著優(yōu)化的進(jìn)行，結(jié)構(gòu)中低效單元逐漸刪除，高效單元得以保留，應(yīng)力離散度值越來越低。因此每次迭代選擇應(yīng)力離散度最小的種群為當(dāng)前種群參與后續(xù)迭代優(yōu)化。式中，N表示實(shí)單元數(shù)量，σVMi為實(shí)單元i的相當(dāng)應(yīng)力，σVMmax為結(jié)構(gòu)最大相當(dāng)應(yīng)力。

（6）判斷體積約束和收斂條件是否滿足，若滿足則終止優(yōu)化，否則返回至第（2）步。

G-BESO方法流程如圖1所示。

圖1 G-BESO方法流程圖Fig.1 G-BESO flow chart

3 算例分析

3.1 L型支架優(yōu)化

圖2所示為經(jīng)典L 型支架優(yōu)化問題。支架的上端固支，右側(cè)頂端受豎直向下的集中力F=4 N。材料彈性模量E=206 GPa，泊松比ν=0.3。體積分?jǐn)?shù)約束為0.4。將結(jié)構(gòu)劃分為1 600個(gè)4×4 mm 四邊形單元。過濾半徑rmin=12 mm，交叉概率Pc?min=0.6，Pc?max=1.0，變異概率Pm?min=0.6，Pm?max=1.0，權(quán)重因子δ=0.5。為便于與標(biāo)準(zhǔn)BESO 法進(jìn)行對(duì)比，相同參數(shù)設(shè)為一致，BESO法的其他參數(shù)設(shè)置為：進(jìn)化率取1%～5%，最大添加率6%。為驗(yàn)證G-BESO方法的穩(wěn)健性，連續(xù)計(jì)算5次。

圖2 L型支架初始結(jié)構(gòu)Fig.2 L-bracket initial design

由于標(biāo)準(zhǔn)BESO 方法中的兩種方法（硬刪除方法和軟刪除方法）的計(jì)算結(jié)果類似，為簡(jiǎn)潔計(jì)算，以下均僅列出前者的結(jié)果。為便于描述，標(biāo)準(zhǔn)BESO 方法中的硬刪除BESO 方法簡(jiǎn)稱為H-BESO，采用本文提出的材料插值模型的軟刪除BESO 方法簡(jiǎn)稱為S-BESO，本文提出的基于改進(jìn)遺傳算法的雙向漸進(jìn)結(jié)構(gòu)優(yōu)化方法簡(jiǎn)稱為G-BESO。

三種方法所得優(yōu)化構(gòu)型和結(jié)果數(shù)據(jù)如表1和表2所示。

表1 不同優(yōu)化方法對(duì)L型支架優(yōu)化所得拓?fù)錁?gòu)型Tab.1 Topological layouts of L-bracket based on different methods

表2 不同優(yōu)化方法對(duì)L型支架優(yōu)化結(jié)果數(shù)據(jù)（柔順度單位：×10-3 N?mm）Tab.2 Optimization data of L-bracket based on different methods(Compliance unit:×10-3 N?mm)

對(duì)比S-BESO和H-BESO方法，隨著進(jìn)化率由1%增大到5%，前者均能得到最優(yōu)拓?fù)錁?gòu)型，而后者在進(jìn)化率增大到5%時(shí)優(yōu)化失?。ū?），這體現(xiàn)出本文提出的考慮材料歷史信息的插值模型的優(yōu)勢(shì)。該插值模型基于單元密度-單元權(quán)重系數(shù)-材料彈性模量插值，單元的密度是逐步減小而非（在標(biāo)準(zhǔn)BESO 方法中）從1 直接減到0 或某極小值（0.001），且部分誤刪的高效單元也能通過密度正向生長(zhǎng)而重新被添加到結(jié)構(gòu)中，因此單元誤刪率更低。

此外，本文提出的基于改進(jìn)遺傳算法的G-BESO 方法不僅在5 次運(yùn)行中均收斂到最優(yōu)結(jié)構(gòu)，而且平均柔順度8.395 9×10-3N?mm 也優(yōu)于其他兩種方法所得值，證明了其穩(wěn)健性和有效性；平均迭代步數(shù)44.4明顯少于另外兩種方法所需步數(shù)，表明其優(yōu)化效率也更高。

3.2 短懸臂梁優(yōu)化

圖3所示為經(jīng)典短懸臂梁算例。梁的左端固支、右端中點(diǎn)承受豎直向下集中力F=100 N。梁的尺寸為長(zhǎng)L=120 mm，寬H=60 mm，厚度t=1 mm；材料彈性模量E=206 GPa，泊松比ν=0.3。體積分?jǐn)?shù)約束為0.5。將結(jié)構(gòu)劃分為60×30 個(gè)四邊形單元。過濾半徑rmin=6 mm，其他參數(shù)均與上例一致。

從最終優(yōu)化得到的拓?fù)錁?gòu)型（表3）上看，H-BESO 方法受進(jìn)化率的影響較大，5次優(yōu)化下只有2次（3%、4%）得到最優(yōu)構(gòu)型，S-BESO 方法受進(jìn)化率的影響相對(duì)較小，3次（2%、4%、5%）得到最優(yōu)構(gòu)型，再次表明本文建立的材料插值模型的優(yōu)勢(shì)；受益于遺傳算法的全局尋優(yōu)能力，G-BESO 方法5 次運(yùn)行均得到最優(yōu)構(gòu)型，表明G-BESO 方法的穩(wěn)健性和有效性；從計(jì)算效率和優(yōu)化目標(biāo)值（表4）上看，G-BESO方法平均迭代步數(shù)38.6、平均柔順度1.5516 N?mm均優(yōu)于其他兩種方法。

圖3 短懸臂梁初始結(jié)構(gòu)Fig.3 Short cantilever initial design

表3 不同優(yōu)化方法對(duì)短懸臂梁優(yōu)化所得拓?fù)錁?gòu)型Tab.3 Topological layouts of short cantilever based on different methods

續(xù)表3

表4 不同優(yōu)化方法對(duì)短懸臂梁優(yōu)化結(jié)果數(shù)據(jù)（柔順度單位：N?mm）Tab.4 Optimization data of short cantilever based on different methods(Compliance unit:N?mm)

4 結(jié) 論

本文提出的G-BESO 方法綜合了漸進(jìn)結(jié)構(gòu)優(yōu)化法，能得到簡(jiǎn)單清晰的拓?fù)錁?gòu)型以及遺傳算法全局尋優(yōu)能力強(qiáng)的優(yōu)勢(shì)，計(jì)算結(jié)果表明：

（1）建立的基于單元密度-單元權(quán)重系數(shù)-材料彈性模量插值的材料插值模型，考慮了單元密度的歷史信息，相較標(biāo)準(zhǔn)BESO方法單元誤刪率更低，優(yōu)化結(jié)果更好。

（2）對(duì)基本遺傳算法所做的改進(jìn)及與BESO 方法的結(jié)合，賦予G-BESO 方法更強(qiáng)的全局尋優(yōu)能力、更高的計(jì)算效率。

（3）G-BESO方法能穩(wěn)定得到最優(yōu)拓?fù)錁?gòu)型而且計(jì)算效率高，具有應(yīng)用到實(shí)際工程結(jié)構(gòu)優(yōu)化設(shè)計(jì)之中的潛力和價(jià)值。