基于改進遺傳算法的雙向漸進結構優化方法研究

吳貝尼，夏利娟

（1.上海交通大學海洋工程國家重點實驗室，上海200240；2.高新船舶與深海開發裝備協同創新中心，上海200240）

0 引 言

結構拓撲優化是在給定的設計空間里確定滿足給定約束的最優材料分布以求獲得最佳結構性能。在過去幾十年中，各種優化方法已被成功運用到各類工程結構拓撲優化設計中，包括均勻化方法[1]、變密度法[2-4]、雙向漸進結構優化方法（BESO）[5-7]、獨立連續映射方法[8]和水平集方法[9]等。其中，BESO方法原理簡單，收斂性好，通用性強，易于編程實現，因此備受青睞。

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

遺傳算法基于適者生存思想，有很強的全局尋優能力，而且遺傳算法中優勝劣汰的思想十分契合結構拓撲優化中高效單元保留、低效單元淘汰的原則。因此，將遺傳算法結合到BESO方法中，能有效克服BESO方法的缺陷，進一步提高優化求解質量。

基于此，本文提出了一種基于改進遺傳算法的軟刪除雙向漸進結構優化法（G-BESO），采用一種考慮單元密度歷史信息的材料插值模型以增強恢復誤刪高效單元的能力，利用遺傳算法中的交叉和變異操作啟發式更新結構狀態以提高全局尋優能力。結合有限元平衡方程建立了連續體結構優化模型并給出具體求解流程，最后用經典算例驗證了方法的有效性和穩健性，同時對比標準BESO 方法表明了其優越性。

1 基于改進遺傳算法的雙向漸進結構優化模型構建

1.1 算法思想

首先利用有限單元將連續體結構離散化，設計區域為種群，離散單元為個體，以離散的單元密度值作為設計變量（0 或1），其中0 代表孔洞單元，1 表示實體單元；對于個體（單元）使用二進制編碼，編碼中1的數目代表著其重要程度即適應度。適應度決定個體編碼的更新方式，個體之間通過交叉、變異更新其編碼，編碼決定個體保留與淘汰。通過這種方式，實現種群的進化即結構的優化。

標準BESO方法中每一步迭代優化通常需要刪除由進化率決定的固定體積的單元，完全根據靈敏度大小機械地增添單元。本文提出的G-BESO 方法基于優勝劣汰原則，根據由適應度決定的單元編碼中1的數量啟發式增添不定體積的單元。具體而言，若單元為實單元，且其編碼值均為0，則從結構中刪除；若為空單元，且其編碼中1的數量達到一定值，則添加到結構中。

1.2 數學模型

本文研究連續體結構在給定體積約束下的最大剛度設計問題，拓撲優化模型可以描述為

確定設計變量：

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

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

1.3 靈敏度計算與收斂準則

單元靈敏度值代表其對整個結構剛度貢獻的大小，關系著單元的刪除與保留，決定最終結構拓撲形狀。對于剛度最大優化問題，基于所提出的材料插值模型的單元應變能靈敏度為

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

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

式中，N為過濾范圍內節點總數，rmin為過濾半徑，rij為單元i中心與節點j之間的距離。

同時，為穩定優化過程，當前單元靈敏度需要根據上一步所得值進行修正[11]：

結構拓撲優化是一個不斷更新單元狀態的迭代過程，直到滿足體積約束和收斂準則，迭代終止。本文采用文獻[11]中推薦的收斂終止公式：

式中：τ為收斂因子即迭代精度，本文取0.001；k為迭代步；N為常數，一般取為5。表征結構的應變能在連續10次以上的迭代中變化量足夠小即認為得到穩定結構。

2 基于改進遺傳算法的雙向漸進結構優化求解流程

2.1 技術手段

為提高G-BESO 方法的實用性，使之能用于優化復雜的工程結構，本文以Python 為編程手段，通過調用HyperMesh軟件更新有限元模型，調用Nastran求解器進行有限元計算，可實現對實際工程結構的拓撲優化設計。

2.2 求解步驟

基于上一章構建的拓撲優化模型，給出G-BESO方法的具體求解步驟如下：

（1）使用HyperMesh對連續體結構進行網格劃分，定義設計區域，施加邊界條件和載荷。

（2）利用Nastran進行結構有限元分析，計算過濾后的單元靈敏度。

（3）根據靈敏度將單元劃分為三個集合：高效單元集、過渡單元集和低效單元集。若為第一次進化，對設計區域單元進行編碼：對所有單元賦予少量0、較多1 的編碼。計算進度指標Prg，交叉概率Pc，變異概率Pm。

式中，Vcur，V*分別為結構當前體積分數和約束體積分數。

為解決遺傳算法隨機性引發的跳躍解、非可行解（結構不協調）等問題，參照文獻[14]，引入懲罰因子如下，使得Pc與Pm隨結構優化逐漸增大：

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

（4）個體（單元）之間實施交叉、變異以更新單元編碼。

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

變異操作只作用于高效單元集和低效單元集，直接改變單元二進制編碼中的數位值[14]。對于每個數位來說，其變異的概率設置為Pm，如果產生的隨機數值小于Pm，高效單元編碼中0 變為1，低效單元中的1變為0。這樣，變異可以使高效單元編碼中含有更多的1，而低效單元編碼中含有更多的0，結構逐步優化。

（5）根據單元編碼中1的數量確定單元狀態，更新設計變量，調用HyperMesh修改結構。

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

應力離散度表征結構受力的均勻程度。隨著優化的進行，結構中低效單元逐漸刪除，高效單元得以保留，應力離散度值越來越低。因此每次迭代選擇應力離散度最小的種群為當前種群參與后續迭代優化。式中，N表示實單元數量，σVMi為實單元i的相當應力，σVMmax為結構最大相當應力。

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

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

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

3 算例分析

3.1 L型支架優化

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

圖2 L型支架初始結構Fig.2 L-bracket initial design

由于標準BESO 方法中的兩種方法（硬刪除方法和軟刪除方法）的計算結果類似，為簡潔計算，以下均僅列出前者的結果。為便于描述，標準BESO 方法中的硬刪除BESO 方法簡稱為H-BESO，采用本文提出的材料插值模型的軟刪除BESO 方法簡稱為S-BESO，本文提出的基于改進遺傳算法的雙向漸進結構優化方法簡稱為G-BESO。

三種方法所得優化構型和結果數據如表1和表2所示。

表1 不同優化方法對L型支架優化所得拓撲構型Tab.1 Topological layouts of L-bracket based on different methods

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

對比S-BESO和H-BESO方法，隨著進化率由1%增大到5%，前者均能得到最優拓撲構型，而后者在進化率增大到5%時優化失敗（表1），這體現出本文提出的考慮材料歷史信息的插值模型的優勢。該插值模型基于單元密度-單元權重系數-材料彈性模量插值，單元的密度是逐步減小而非（在標準BESO 方法中）從1 直接減到0 或某極小值（0.001），且部分誤刪的高效單元也能通過密度正向生長而重新被添加到結構中，因此單元誤刪率更低。

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

3.2 短懸臂梁優化

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

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

圖3 短懸臂梁初始結構Fig.3 Short cantilever initial design

表3 不同優化方法對短懸臂梁優化所得拓撲構型Tab.3 Topological layouts of short cantilever based on different methods

續表3

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

4 結 論

本文提出的G-BESO 方法綜合了漸進結構優化法，能得到簡單清晰的拓撲構型以及遺傳算法全局尋優能力強的優勢，計算結果表明：

（1）建立的基于單元密度-單元權重系數-材料彈性模量插值的材料插值模型，考慮了單元密度的歷史信息，相較標準BESO方法單元誤刪率更低，優化結果更好。

（2）對基本遺傳算法所做的改進及與BESO 方法的結合，賦予G-BESO 方法更強的全局尋優能力、更高的計算效率。

（3）G-BESO方法能穩定得到最優拓撲構型而且計算效率高，具有應用到實際工程結構優化設計之中的潛力和價值。