齊 峻，武 星

（長(zhǎng)安大學(xué)工程機(jī)械學(xué)院，陜西 西安 710064）

0 引言

拓?fù)鋬?yōu)化算法自20世紀(jì)被提出以來(lái)，一直在不斷更新。拓?fù)鋬?yōu)化的目標(biāo)是在滿足各種設(shè)計(jì)約束的情況下，找到材料的最優(yōu)分布[1]。雙向漸進(jìn)結(jié)構(gòu)優(yōu)化法（BESO）是一種將漸進(jìn)結(jié)構(gòu)優(yōu)化法（ESO）和材料添加策略算法（AESO）相結(jié)合的方法[2]。最初的ESO主要是去除最低應(yīng)力單元，經(jīng)過(guò)發(fā)展，BESO在ESO原有的基礎(chǔ)上可實(shí)現(xiàn)增加單元。原始的硬殺傷ESO/BESO方法消除了設(shè)計(jì)領(lǐng)域的元素，通常，作為設(shè)計(jì)變量的元件不能直接從設(shè)計(jì)領(lǐng)域中消除，除非軟元件完全等同于空心元件。換言之，對(duì)于一個(gè)新的拓?fù)鋬?yōu)化問(wèn)題，開發(fā)一種軟殺傷BESO方法，然后探索硬殺傷元件的可能性可能會(huì)更可靠[3-5]。

近年來(lái)，遺傳算法（GAs）[6]也不斷地發(fā)展，并逐漸應(yīng)用于拓?fù)鋬?yōu)化問(wèn)題[7]。由于GAs的隨機(jī)性，其直接應(yīng)用于拓?fù)鋬?yōu)化問(wèn)題時(shí)會(huì)產(chǎn)生棋盤格現(xiàn)象。Liu等[8]很好地意識(shí)到了GAs和ESO結(jié)合的問(wèn)題，但由于在設(shè)計(jì)過(guò)程中無(wú)法實(shí)現(xiàn)材料的增加，得到的結(jié)構(gòu)可能無(wú)法滿足使用要求。Zuo等[9]提出了GAs和BESO相結(jié)合的拓?fù)鋬?yōu)化方法，實(shí)現(xiàn)了迭代過(guò)程中材料的增減。NSGA-II是在GAs的基礎(chǔ)上發(fā)展起來(lái)的，通過(guò)快速的非支配排序和擁擠計(jì)算、競(jìng)爭(zhēng)選擇和精英策略來(lái)實(shí)現(xiàn)快速優(yōu)化過(guò)程。在全局優(yōu)化方面，NSGA-II顯現(xiàn)出比GAs更大的優(yōu)勢(shì)。

本文主要部分說(shuō)明如下：第1節(jié)，闡述了拓?fù)鋬?yōu)化的基本問(wèn)題；第2節(jié)，根據(jù)傳統(tǒng)的靈敏度分析，推導(dǎo)出了固體和軟單元的靈敏度數(shù)值；第3節(jié)，簡(jiǎn)要說(shuō)明了本文程序的原理和數(shù)值實(shí)現(xiàn)；第4節(jié)，通過(guò)具體的數(shù)值算例證明所提出的新方法的穩(wěn)定性和有效性。

1 問(wèn)題陳述

在拓?fù)鋬?yōu)化問(wèn)題中，往往追求給定體積分?jǐn)?shù)下的最大剛度。剛度通常用應(yīng)變能數(shù)值表示，對(duì)于具有n個(gè)離散單元的拓?fù)鋬?yōu)化問(wèn)題可以表示為：

式中，C為平均柔度，在設(shè)計(jì)中，通過(guò)計(jì)算平均柔度來(lái)達(dá)到最大剛度。u為全局位移向量，K為全局剛度矩陣。vi是單元i的體積，V*是目標(biāo)體積。

2 靈敏度數(shù)值計(jì)算和過(guò)濾方案

2.1 靈敏度數(shù)值計(jì)算

設(shè)計(jì)變量通過(guò)基于梯度的方法更新。對(duì)于單相固體材料-空洞優(yōu)化問(wèn)題，通過(guò)推導(dǎo)目標(biāo)函數(shù)可以得到單元靈敏度值：

式中，p表示懲罰指數(shù)，xi表示單相固體材料的單元密度值，當(dāng)某單元固相材料處于激發(fā)狀態(tài)時(shí)，xi=1，否則，關(guān)閉此激活狀態(tài)，單元最終轉(zhuǎn)化為空洞狀態(tài)，此時(shí)xi=xmin。ui表示第i個(gè)單元的位移矢量，ki表示第i個(gè)單元的剛度矩陣。

靈敏度值可用式（5）表示：

2.2 過(guò)濾方案

采用靈敏度濾波策略對(duì)靈敏度進(jìn)行修正。這不僅有效減少了棋盤格現(xiàn)象，還通過(guò)過(guò)濾半徑內(nèi)物理單元的靈敏度值推導(dǎo)出空單元的靈敏度，實(shí)現(xiàn)單元的增刪。本文的過(guò)濾策略是將與節(jié)點(diǎn)關(guān)聯(lián)的所有元素的靈敏度值均質(zhì)化為節(jié)點(diǎn)的元素靈敏度值，如式（6）所示：

式中，H為設(shè)計(jì)域的節(jié)點(diǎn)總數(shù)，w(rik)為權(quán)重因子，rik表示元素i的中心到節(jié)點(diǎn)k的距離。權(quán)重因子可以表示為：

式中，rmin表示過(guò)濾器半徑。為了提高方法優(yōu)化的收斂性，采用考慮歷史迭代信息的收斂策略來(lái)修改靈敏度值。具體操作是計(jì)算本次迭代和上次迭代的靈敏度的平均值，將其作為靈敏度的更新值。

式中， ,?igα 和 ,1?ig-α 分別為單元i第g次和g-1次迭代的靈敏度值。最后，根據(jù)該靈敏度值測(cè)量各單元對(duì)結(jié)構(gòu)性能的影響，然后將體積作為拓?fù)鋬?yōu)化設(shè)計(jì)的約束條件。

3 算法的提出與實(shí)現(xiàn)

3.1 算法提出

在本文提出的算法中，每個(gè)元素都被視為一個(gè)獨(dú)立的個(gè)體進(jìn)行分析，在迭代過(guò)程中保持種群維數(shù)穩(wěn)定，允許新種群代替原種群，但不允許合并。每個(gè)單元都由一定長(zhǎng)度的二進(jìn)制字符串表示。初始字符串在本文中定義為全“1”的數(shù)組，表示所有單元均為固相材料。一個(gè)全為“0”的數(shù)組意味著該單元將在迭代中被刪除。數(shù)字“1”代表材料屬性的強(qiáng)度。交叉和變異算子只有“0”和“1”之間的轉(zhuǎn)換。

通過(guò)初始化種群，實(shí)現(xiàn)設(shè)計(jì)域內(nèi)所有個(gè)體的編碼，單位編碼的初始值為長(zhǎng)度為6的矩陣。對(duì)模型進(jìn)行預(yù)處理后，按式（5）計(jì)算各元素的靈敏度值。此時(shí)，可以根據(jù)每個(gè)迭代步驟的目標(biāo)體積分?jǐn)?shù)設(shè)置靈敏度閾值，并根據(jù)元素靈敏度值進(jìn)行排序。對(duì)于具有相同靈敏度的單元，計(jì)算它們的擁擠程度，使同級(jí)個(gè)體之間存在清晰的層級(jí)關(guān)系。如果種群中目標(biāo)函數(shù)的最大值等于目標(biāo)函數(shù)的最小值，則將種群的擁擠程度設(shè)置為無(wú)窮大。一般來(lái)說(shuō)，擁擠的表達(dá)式如下：

3.2 算法實(shí)現(xiàn)

本文算法步驟如下：

1）通過(guò)有限元思想進(jìn)行預(yù)處理，主要是網(wǎng)格劃分、施加載荷和設(shè)置邊界條件[10]。

2）初始化種群，個(gè)體由一系列“1”碼組成，根據(jù)靈敏度計(jì)算公式加上懲罰因子和過(guò)濾方案的實(shí)現(xiàn)來(lái)計(jì)算靈敏度。指定懲罰指數(shù)（p）和體積進(jìn)化率（ER）等參數(shù)。

3）根據(jù)靈敏度值迭代進(jìn)化種群。主要進(jìn)行快速的非支配排序和交叉變異操作。

4）刪除設(shè)計(jì)域中單元基因串中僅包含數(shù)字“0”的元素。如果上次迭代中刪除的元素在當(dāng)前生成中具有高靈敏度值，則將該元素添加到設(shè)計(jì)域中。

5）根據(jù)靈敏度級(jí)別和主導(dǎo)關(guān)系將實(shí)體材料屬性分配給設(shè)計(jì)域中的其余元素。

6）檢查設(shè)計(jì)域的剩余體積和收斂標(biāo)準(zhǔn)。

7）重復(fù)步驟2至步驟6，直到達(dá)到目標(biāo)體積和收斂標(biāo)準(zhǔn)。

對(duì)于收斂準(zhǔn)則，當(dāng)滿足體積約束時(shí)，考慮了兩個(gè)連續(xù)優(yōu)化迭代之間最大柔順度的絕對(duì)變化，兩次連續(xù)優(yōu)化迭代的收斂性之差小于0.005。本文中S=5，即連續(xù)優(yōu)化迭代次數(shù)為5。

4 數(shù)值算例

本節(jié)通過(guò)數(shù)值算例來(lái)驗(yàn)證本文提出的算法的有效性、穩(wěn)定性和優(yōu)越性。考慮了一種虛擬材料，材料的物理性質(zhì)如表1所示。所涉及的算法均在Win10 Intel?CoreTM i5-9300 CPU、2.40 GHz、16 G運(yùn)行內(nèi)存的桌面上執(zhí)行。

表1 虛擬材料的物理性質(zhì)

如圖1所示，考慮了尺寸為50 mm×30 mm的二維懸臂梁結(jié)構(gòu)。懸臂梁結(jié)構(gòu)左側(cè)完全固定，右側(cè)中間施加1 N的集中荷載。將體積約束設(shè)置為初始體積的35%。網(wǎng)格大小設(shè)置為1，共1 500個(gè)單元，懲罰指數(shù)p=3，過(guò)濾器半徑rmin=1.5 mm，體積進(jìn)化率ER=0.02。

圖1 懸臂梁模型

4組懸臂梁結(jié)構(gòu)最終迭代結(jié)果如圖2所示。從圖2可以看出，模型會(huì)產(chǎn)生良好的布局效果，邊界更清晰。由于NSGA-II具有隨機(jī)性，每次迭代的結(jié)果可能會(huì)有所不同，但比較4次運(yùn)行的結(jié)果，結(jié)構(gòu)布局顯示出相似之處。懸臂梁模型在設(shè)置相同參數(shù)條件下通過(guò)BESO方法所迭代優(yōu)化出的最終結(jié)果如圖3所示。4組懸臂梁最終平均柔度結(jié)果如表2所示。通過(guò)比較表2中4個(gè)運(yùn)行結(jié)果的目標(biāo)函數(shù)，相對(duì)于平均值的誤差和最大誤差均不超過(guò)1%，4個(gè)運(yùn)行結(jié)果的最終體積分?jǐn)?shù)為35%，滿足約束條件。通過(guò)數(shù)值算例分析，證明該算法具有良好的穩(wěn)定性和收斂性。

圖2 4組懸臂梁結(jié)構(gòu)最終迭代結(jié)果

圖3 BESO方法下懸臂梁最終結(jié)果

表2 4組懸臂梁最終平均柔度結(jié)果

與BESO方法相比，本文提出的新方法具有更強(qiáng)的進(jìn)化能力。傳統(tǒng)的BESO方法在每次迭代過(guò)程中直接增刪滿足條件的單元，但不同的是本文提出的新方法引入了NSGA-II的演化機(jī)制。在迭代過(guò)程中，不會(huì)直接刪除敏感度閾值以下的單元格或增加單元格。相反，將“0”和“1”的數(shù)量用作添加和刪除的標(biāo)準(zhǔn)。不斷地通過(guò)“0”和“1”之間的演化，當(dāng)一個(gè)單元出現(xiàn)全“0”的字符串時(shí)，就意味著一次刪除操作。該算法避免了每次迭代過(guò)程中的多次增刪改查，有利于算法逐步優(yōu)化直至收斂。

本文提出的新算法與原始BESO方法的對(duì)比如表3所示，設(shè)置相同的約束、材料和參數(shù)。由于本文提出的新方法具有很強(qiáng)的全局搜索能力，保證了優(yōu)化的效果。本文算法與原始BESO方法相比，整體在效果上和原始BESO方法近乎相同。考慮到遺傳算法具有隨機(jī)性，其在靈敏度計(jì)算、快速非支配排序這兩個(gè)方面每次的運(yùn)行時(shí)間會(huì)較長(zhǎng)，但是通過(guò)對(duì)比本文算法與原始BESO方法的總運(yùn)行時(shí)間與總迭代次數(shù)，發(fā)現(xiàn)二者幾乎無(wú)差距，證明了本算法的效率和可靠性。

表3 懸臂梁最終平均柔度比較結(jié)果

5 結(jié)論

研究小組提出了一種新的結(jié)構(gòu)輕量化設(shè)計(jì)拓?fù)鋬?yōu)化方法。通過(guò)將NSGA-II算法融入BESO方法中，保證了算法的全局優(yōu)化效果，實(shí)現(xiàn)了結(jié)構(gòu)的剛度最大化，獲得了干凈的拓?fù)浣Y(jié)構(gòu)。通過(guò)數(shù)值算例分析，表明本文提出的新方法也能得到類似于BESO的結(jié)構(gòu)布局，證明了算法的穩(wěn)定性和有效性。本文算法通過(guò)結(jié)合遺傳算法達(dá)到了全局最優(yōu)搜索的效果，在相同約束條件下與傳統(tǒng)BESO算法比較，發(fā)現(xiàn)本文算法可以提高結(jié)構(gòu)剛度，證明了本文算法的優(yōu)越性。顯然，研究小組提出的方法可以更好地滿足復(fù)雜工況下的輕量化結(jié)構(gòu)設(shè)計(jì)，可以應(yīng)用于各個(gè)領(lǐng)域的工程問(wèn)題，提高結(jié)構(gòu)的性能，對(duì)工程結(jié)構(gòu)的優(yōu)化具有指導(dǎo)意義。