基于比例微分優(yōu)化準則的拓撲優(yōu)化方法

趙龍彪 高 亮 陳志敏 邱浩波

1．華中科技大學(xué)數(shù)字制造裝備與技術(shù)國家重點實驗室，武漢，430074 2．中國艦船研究設(shè)計中心，武漢，430064

0 引言

拓撲優(yōu)化是在一個確定的區(qū)域內(nèi)尋求滿足設(shè)計約束的最優(yōu)拓撲結(jié)構(gòu)的優(yōu)化方法。設(shè)計目標和約束可以為重量、頻率、強度等。其主要求解思路是將尋求結(jié)構(gòu)的最優(yōu)拓撲結(jié)構(gòu)問題轉(zhuǎn)化為在給定的設(shè)計空間內(nèi)尋求最優(yōu)材料的分布問題［1］。它是一種能幫助設(shè)計者選擇合適的初始拓撲結(jié)構(gòu)的有效方法，在概念設(shè)計階段具有十分重要的指導(dǎo)意義。

在拓撲優(yōu)化的眾多模型求解方法中，優(yōu)化準則法（optimality criteria method，OC）［2］是根據(jù)滿足各種約束條件（應(yīng)力、位移、頻率等）的最佳準則，從可行的設(shè)計中找出最佳方案的方法，以充分發(fā)揮材料的剛度、強度和穩(wěn)定性的潛力，實現(xiàn)等強度、等應(yīng)變能的最佳傳導(dǎo)應(yīng)力路徑。最早的優(yōu)化準則法包括互應(yīng)變能法（constant mutual energy design）、應(yīng)力比法（stress ratio method）和滿應(yīng)力法（fully stressed design）［3］。目前，學(xué)者們已經(jīng)對OC法有了比較深入的探討，Meske等［4］對OC法的健壯性和快速收斂性進行了研究，并應(yīng)用在結(jié)構(gòu)的固有頻率優(yōu)化上。Logo等［5］把OC法應(yīng)用在隨機拓撲優(yōu)化設(shè)計問題上。Chiandussi等［6］把遺傳路徑加入到OC法中。

優(yōu)化準則法的每一次迭代都是根據(jù)上一次迭代中各單元的材料密度和敏度信息，按照固定法則來計算各單元新的材料密度。在求解過程中，搜索效率不高，容易產(chǎn)生局部震蕩和難以收斂的現(xiàn)象，會降低求解的速度。鑒于此，本文在傳統(tǒng)OC法的基礎(chǔ)上，基于SIMP密度函數(shù)插值模型，以結(jié)構(gòu)的柔度最小化為目標函數(shù)，借助比例微分控制的思想，提出了比例微分優(yōu)化準則法（proportional and differential optimality criterion method，PDOC），改進了迭代算子，給出了設(shè)計變量的迭代更新方案。在提高搜索效率的同時，加快了優(yōu)化速度，并使得求解更容易收斂。

1 結(jié)構(gòu)拓撲優(yōu)化設(shè)計基礎(chǔ)

1.1 結(jié)構(gòu)拓撲優(yōu)化模型

結(jié)構(gòu)拓撲優(yōu)化的目的是確定設(shè)計域中材料的最優(yōu)分布。為實現(xiàn)材料的最優(yōu)分布，需要進行材料密度的連續(xù)化插值。變密度法是人為假定單元的密度和材料物理屬性（如許用應(yīng)力、彈性模量）之間的某種對應(yīng)關(guān)系，以連續(xù)變量的密度函數(shù)形式顯式地表達這種對應(yīng)關(guān)系。變密度法主要的密度－剛度插值方式有帶懲罰指數(shù)的固體各向同性微結(jié)構(gòu)模型SIMP［7］和材料屬性的理性近似模型RAMP［8］兩種，本文以SIMP模型為例討論連續(xù)體結(jié)構(gòu)拓撲優(yōu)化模型。

SIMP模型主要通過引入懲罰因子，在材料的彈性模量和單元相對密度之間建立起一種顯式的非線性對應(yīng)關(guān)系。其目的是當設(shè)計變量的值在（0，1）之間時，對中間密度值進行懲罰，使其逐漸向0／1兩端聚集，這樣可以使連續(xù)變量的拓撲優(yōu)化模型能很好地逼近原來0－1離散變量的優(yōu)化模型。

SIMP材料模型的數(shù)學(xué)表達形式：

式中，xj為一般設(shè)計變量，表示離散單元的相對密度；p為SIMP模型中對中間密度材料的懲罰因子；E0、Emin分別為固體和空洞部分材料的彈性模量；Ep為插值后的彈性模量；K為插值以后的剛度矩陣；Kj為第j個單元固體材料的剛度矩陣。

懲罰因子的作用是當設(shè)計變量的值在［0，1］之間時，通過逐漸增加p的值對設(shè)計變量的中間值進行懲罰，隨著p值的增大，設(shè)計逐漸接近0／1設(shè)計。為有效壓縮中間密度材料，要求p＞2［9］。結(jié)構(gòu)拓撲優(yōu)化過程中，通過不斷地優(yōu)化迭代計算，來保留對結(jié)構(gòu)傳力路徑有利的單元，刪除對結(jié)構(gòu)傳力路徑作用不大的單元，從本質(zhì)上講，結(jié)構(gòu)拓撲優(yōu)化問題是一個單元集合的增減問題。對于連續(xù)體拓撲優(yōu)化模型，對每個單元j的增減操作變成對xj（j＝1，2，…，n）值的加減操作，并且使其值在［0，1］之間變化。

以結(jié)構(gòu)拓撲優(yōu)化中比較典型的最小柔度問題為例，在SIMP材料插值方法的基礎(chǔ)上，拓撲優(yōu)化模型為

式中，C為目標函數(shù)，定義C為結(jié)構(gòu)的總體柔度；F為單元載荷矢量；U為單元位移列陣；K為結(jié)構(gòu)總體剛度矩陣；V1為每個單元的設(shè)計域的初始體積；V0為設(shè)計域的初始體積；V為優(yōu)化后的結(jié)構(gòu)體積；f為優(yōu)化體積比；uj為單元位移列向量。

1.2 優(yōu)化準則法

優(yōu)化準則法（OC）是根據(jù)工程經(jīng)驗、力學(xué)概念以及數(shù)學(xué)規(guī)劃的最優(yōu)條件，預(yù)先建立某種準則，通過相應(yīng)的迭代方法，獲得滿足這一準則的設(shè)計方案，作為問題的最優(yōu)解的一種優(yōu)化方法。在實際運用中，OC法可以并行處理設(shè)計空間的所有單元，快速地發(fā)現(xiàn)應(yīng)力傳遞的主要路徑，能夠求解超大規(guī)模設(shè)計變量的拓撲優(yōu)化問題。

優(yōu)化準則法中應(yīng)用最成功的是Kuhn－Tucker條件。下式所示的是不等式約束的多元函數(shù)極值問題，其中X＝ （x1，x2，…，xn）T為設(shè)計變量，受到m個不等式約束。

拓撲優(yōu)化問題是典型的不等式約束的多元函數(shù)極值問題，它通過Kuhn－Tucker條件構(gòu)造多元函數(shù)的極值條件函數(shù)，通過求極值條件函數(shù)間接求多元函數(shù)的極值點。優(yōu)化準則法從一個空間的初始設(shè)計點x（k）出發(fā)，著眼于每次迭代應(yīng)滿足的優(yōu)化條件，依據(jù)迭代式：

式中，ζ為阻尼因子（0＜ζ＜1）；k為迭代步數(shù)；Vj為體積比。得到一個改進的設(shè)計x（k＋1）。式（6）中常常引入一些經(jīng)驗系數(shù)來調(diào)整優(yōu)化過程的收斂性和穩(wěn)定性，如步長因子、阻尼因子等。

總體來講，OC法簡單明了、易于理解，但其設(shè)計變量更新來源于一種啟發(fā)式的迭代方式，在實際的應(yīng)用中還存在著一些缺陷：

（1）設(shè)計變量的更新通過預(yù)定的迭代控制策略進行搜索，總是按照預(yù)先規(guī)定的路線進行，這種搜索具有盲目性，效率不高，會降低求解速度。

（2）求解通用性較差，對于復(fù)雜的目標函數(shù)，不容易構(gòu)造設(shè)計變量的迭代公式，搜索過程中的指導(dǎo)信息只有當前結(jié)果和目標結(jié)果的偏差，所以容易出現(xiàn)震蕩和難以收斂的現(xiàn)象。

（3）一般比較適用于單目標函數(shù)、單約束條件下的優(yōu)化問題求解，多目標問題需要推導(dǎo)不同的優(yōu)化準則，難以適用于復(fù)雜問題的求解。

針對這些情況，我們在OC法的基礎(chǔ)上，通過引入比例微分控制（PD）的思想，改進其迭代算子，并給出設(shè)計變量的迭代更新方案，提出了PDOC方法，并將其用于拓撲優(yōu)化中，以求在傳統(tǒng)OC法優(yōu)勢基礎(chǔ)上，進一步加快收斂，提高求解速度，改善優(yōu)化效果。

2 PDOC方法及算法實現(xiàn)

2.1 PD控制基礎(chǔ)

其思想來源于工程控制中廣泛應(yīng)用的比例、積分和微分控制方法，簡稱PID控制。PID控制在實際中又可以根據(jù)情況簡化為PI控制和PD控制。

（1）比例控制（proportional control）。比例控制是一種最簡單的控制方式。其控制器的輸出與輸入誤差信號成比例關(guān)系，可適當縮放信號的幅值?？刂破髦?，比例項的引入，能夠放大誤差的幅值。

（2）微分控制（differential control）。在微分控制中，控制器的輸出與輸入誤差信號的微分（即誤差的變化率）成正比關(guān)系?？刂破髦校⒎猪椀囊?，能夠預(yù)測誤差變化的趨勢。

（3）比例微分控制（proportional and differential control）。單一的比例控制或者微分控制都不能完全解決問題，只有同時具有比例、微分控制的控制器，才能夠提前使抑制誤差的控制作用等于零，甚至為負值，從而避免被控量的嚴重超調(diào)。

因此，對有較大慣性或滯后的被控對象，比例微分（PD）控制器能改善系統(tǒng)在調(diào)節(jié)過程中的動態(tài)特性，通過引入比例項，放大誤差的增幅；增加微分項，便能預(yù)測誤差變化的趨勢。這樣，具有比例、微分環(huán)節(jié)的控制器，就能夠提前使抑制誤差的控制作用等于零，甚至為負值，從而避免了被控量的嚴重超調(diào)。

2.2 PDOC方法介紹

傳統(tǒng)OC法迭代公式計算簡單、易于數(shù)值實現(xiàn)，但其收斂速度不夠。在實際的數(shù)值試驗中，為加快迭代法的收斂速度，可以通過多級定常迭代的思想構(gòu)造迭代公式，如：

即在計算x（k＋1）時考慮前面l個迭代值，這樣的迭代法稱為l級定常迭代法。如果φk與k無關(guān)，迭代公式可表達成以下形式：

多級迭代法在數(shù)值計算過程中比單級迭代法需要保存更多的信息，增加了存儲量。所以一般l不宜太大，常在計算x（k＋1）時考慮前面2個迭代式。迭代式如下：

在實際優(yōu)化問題中，求解的非線性方程組結(jié)構(gòu)復(fù)雜，很難通過確切的數(shù)學(xué)推導(dǎo)過程構(gòu)造二級定常迭代法。這與控制系統(tǒng)中的很多情況類似，我們可以把x（k＋1）看作是控制系統(tǒng)φ（x（k），x（k－1））的輸出變量，而x（k），x（k＋1），…，x（k－l＋1）為控制系統(tǒng)的輸入變量?？刂葡到y(tǒng)φ（x（k），x（k－1））是一個黑匣子，無法得知其內(nèi)在的控制機理。在控制工程中，當被控對象的結(jié)構(gòu)和參數(shù)不能完全明確，或得不到精確的數(shù)學(xué)模型，而控制理論的其他技術(shù)難以采用時，系統(tǒng)控制器的結(jié)構(gòu)和參數(shù)必須依靠經(jīng)驗和現(xiàn)場調(diào)試來確定，這時最適合應(yīng)用PD控制方法。PD控制器就是根據(jù)系統(tǒng)的誤差，利用比例、微分計算出控制量進行控制的。輸入誤差信號可以用差分形式表達如下：

結(jié)合優(yōu)化準則迭代方法和微分控制的思想，可以構(gòu)造比例微分迭代公式如下：

式中，α為微分項的影響因子，簡稱微分參數(shù)。

稱這種結(jié)合優(yōu)化準則迭代方法和微分控制思想的二級定常迭代法為PDOC方法。

考慮迭代算子的移動極限和收斂效果，可以把迭代算子改進為

式中，m為移動極限（0＜m＜1）；ζ為阻尼因子（0＜ζ＜1）；xmin為材料密度的下限值，xmin＝0．0001；Λ 為拉格朗日乘子。

在每一步迭代中，確保體積約束滿足的拉格朗日乘子Λ是變化的，在第k步，可以采用二分法求解Λ：

（4）重復(fù)步驟 （2）和 步驟 （3），直到滿 足｜V（k）－V＊｜≤δ（δ＝0．0001）。

2.3 PDOC算法實現(xiàn)

圖1為基于PDOC法拓撲優(yōu)化設(shè)計的算法流程圖，可以采用 MATLAB實現(xiàn)算法，步驟如下：

（1）有限元模型的前處理，主要包括網(wǎng)格劃分、定義約束和載荷。

（2）初始化設(shè)計變量，定義單元的初始材料密度，默認的值均為0．5。

（3）有限元分析求解，計算出各單元的敏度和剛度，并進行敏度分析與過濾。

（4）根據(jù)上一步迭代更新的各單元材料密度，采用PDOC法計算各單元新的材料密度和結(jié)構(gòu)柔度。

（5）判斷是否達到優(yōu)化設(shè)計目標。如果未達到優(yōu)化設(shè)計目標，則轉(zhuǎn)步驟（3）繼續(xù)優(yōu)化；如果達到優(yōu)化設(shè)計目標，則轉(zhuǎn)步驟（6）。一般終止設(shè)計目標有兩種情況：一是材料密度總和達到最小約束界限；二是結(jié)構(gòu)的整體柔度值的改變量達到預(yù)定界限。

（6）輸出結(jié)果，結(jié)束。

圖1 PDOC法拓撲優(yōu)化求解流程圖

3 實例分析

3.1 懸臂梁剛度優(yōu)化

圖2所示為懸臂梁剛度優(yōu)化問題的設(shè)計域示意圖，其設(shè)計空間長A＝60mm，高B＝20mm。材料的彈性模量為207GPa，泊松比為0．3，體積比為0．5。結(jié)構(gòu)左側(cè)的豎邊受X和Y方向的自由度約束，右側(cè)中間位置受垂直向下的10kN載荷作用。

圖2 懸臂梁柔度最小化問題的設(shè)計域示意圖

首先通過試驗分析PDOC法中的微分參數(shù)對求解效果的影響。試驗電腦CPU主頻率為2．0GHz，內(nèi)存為2GB。采用MATLAB語言實現(xiàn)算法。

網(wǎng)格劃分為120×40，單元數(shù)為4800。分析式（12）可知，當α≠0時，即為PDOC方法，當α＝0時，PDOC方法即退化為OC法。經(jīng)過測試發(fā)現(xiàn)，實際求解過程中，一般取0．5＜α≤0．9時，求解效果較好。為比較PDOC法與OC法的求解效果，我們僅取α＝0和α＝0．8作為典型代表進行對比分析。求解結(jié)果如圖3所示，解的柔度進化曲線如圖4所示。從求解結(jié)果可以看出，當采用PDOC法求解時，求解過程均比較穩(wěn)定，經(jīng)過10步迭代就能接近最優(yōu)值173．68，每步迭代耗13．673s。當采用OC法求解時，大概需經(jīng)過20步迭代才能接近最優(yōu)值，每步迭代耗13．452s。從優(yōu)化過程看，PDOC法比OC法有更快的收斂速度，更容易收斂；從優(yōu)化結(jié)果看，PDOC法比OC法能獲得剛度更大的結(jié)構(gòu)，在加快收斂的同時，提高了結(jié)構(gòu)剛度。

圖3 網(wǎng)格劃分120×40懸臂梁的解

圖4 網(wǎng)格劃分120×40的懸臂梁柔度進化曲線比較圖

3.2 簡支梁剛度優(yōu)化

如圖5所示的簡支梁剛度優(yōu)化問題的設(shè)計域示意圖，其設(shè)計空間長A＝240mm，高B＝40mm。材料的彈性模量為207GPa，泊松比為0．3，體積比為0．5。中間位置受垂直向下的10kN載荷作用。試驗電腦CPU主頻率為2．0GHz，內(nèi)存為2GB，采用 MATLAB語言實現(xiàn)算法。

圖5 簡支梁柔度最小化問題的設(shè)計域示意圖

設(shè)定網(wǎng)格劃分為240×40，單元數(shù)為9600。懲罰因子p＝3，體積比V＊＝0．5，阻尼因子ζ＝0．5，移動極限m＝0．3，最大迭代步數(shù)為40。

求解結(jié)果如圖6所示，解的柔度進化曲線如圖7所示。從求解結(jié)果可以看出，當采用PDOC法求解時，求解過程均比較穩(wěn)定，經(jīng)過15步迭代就能接近最優(yōu)值192．15，每步迭代耗14．792s。當采用OC法求解時，大概需經(jīng)過25步迭代才能接近最優(yōu)值，每步迭代耗14．358s。從優(yōu)化過程看，PDOC法比OC法有更快的收斂速度，更容易收斂；從優(yōu)化結(jié)果看，PDOC法比OC法能獲得剛度更大的結(jié)構(gòu)，在加快收斂的同時，提高了結(jié)構(gòu)剛度。

圖7 網(wǎng)格劃分240×40的簡支梁柔度進化曲線比較圖

3.3 柔性位移反向器

柔性機構(gòu)是通過其部分或全部具有柔性的構(gòu)件變形而產(chǎn)生位移的機械結(jié)構(gòu)，它的設(shè)計主要要求在輸入端施加輸入載荷以后，在輸出端產(chǎn)生位移運動。柔性機構(gòu)是多目標設(shè)計問題，優(yōu)化的兩個目標函數(shù)分別為機構(gòu)的共有應(yīng)變能（mutual strain energy，MSE）和結(jié)構(gòu)的應(yīng)變能（strain energy，SE）［10］。

本文以位移反相器的拓撲優(yōu)化為例，驗證PDOC方法的求解性能。拓撲優(yōu)化的兩個目標函數(shù)分別是最大化輸出端的MSE以及最小化結(jié)構(gòu)的SE。設(shè)計的體積比為30%。設(shè)計域尺寸為80μm×80μm，如圖8所示。利用對稱性，選取設(shè)計域的上半部分進行分析。在上半部分均布80×40個網(wǎng)格。材料彈性模量為1MPa，泊松比為0．3，輸入驅(qū)動力為fin＝0．5mN，uout為輸出位移。在輸入端口固定一個剛度為Kin＝1的彈簧，輸入能量通過輸入驅(qū)動力和彈簧的剛性來表示。在輸出端固定一個剛度為Kout＝1的彈簧來模擬工件對機構(gòu)的作用力。

圖8 位移反相器的設(shè)計域示意圖

求解結(jié)果如圖9所示，解的MSE進化曲線如圖10所示。從求解結(jié)果可以看出，當采用PDOC法求解時，求解過程均比較穩(wěn)定，經(jīng)過15步迭代就能接近最優(yōu)值MSE＝0．821 45，每步迭代耗9．876s。當采用OC法求解時，大概需經(jīng)過40步迭代才能接近最優(yōu)值MSE＝0．789 46，每步迭代耗9．435s。從優(yōu)化過程看，PDOC法比OC法有更快的收斂速度，更容易收斂；從優(yōu)化結(jié)果看，PDOC法比OC法能獲得剛度更大的結(jié)構(gòu)，因此，在多目標問題的拓撲優(yōu)化求解中，PDOC方法仍然具有較好的實際應(yīng)用效果，在加快收斂的同時，提高了結(jié)構(gòu)剛度。

圖9 位移反相器簡支梁的解

圖10 位移反相器MSE進化曲線比較圖

4 結(jié)語

本文在PD控制理論的基礎(chǔ)上，提出了基于PDOC的拓撲優(yōu)化方法。PDOC法通過引入比例微分控制的思想來改進迭代算子，通過構(gòu)造更合理的數(shù)值迭代公式以加快收斂，提高計算速度。引入比例、微分控制的求解算法，以預(yù)測誤差的變化方向，提前抑制誤差的作用，從而避免了被控量的超調(diào)現(xiàn)象和震蕩現(xiàn)象，經(jīng)過實例測試，采用改進后的PDOC方法進行拓撲優(yōu)化，能夠使優(yōu)化求解過程經(jīng)過10次左右迭代就能收斂，明顯比OC法的搜索速度要快，同時還能夠求出柔度更小且結(jié)構(gòu)更清晰的解。另外，通過一個多目標優(yōu)化實例測試，驗證了PDOC方法在多目標拓撲優(yōu)化問題求解中的有效性。綜上所述，PDOC方法是一種比OC方法求解速度更快、求解結(jié)果更好的方法。

［1］ 羅震，陳立平，黃玉盈，等．連續(xù)體結(jié)構(gòu)的拓撲優(yōu)化設(shè)計［J］．力學(xué)進展，2004，34（4）：463－476．

［2］ Zhou M，Rozvany G I N．The COC Algorithm，Part II：Topological Geometrical and Generalized Shape Optimization［J］．Computer Methods in Applied Mechanics and Engineering，1991，89：197－224．

［3］ Levy R，Lavan O．Fully Stressed Design of Passive Controllers in Framed Structures for Seismic Loadings［J］．Structural and Multidisciplinary Optimization，2006，32（6）：485－498．

［4］ Meske R，Lauber B，Schnack E．A New Optimality Criteria Method for Shape Optimization of Natural Frequency Problems［J］．Structural and Multidisciplinary Optimization，2006，31（4）：135－140．

［5］ Logo J．New Type of Optimality Criteria Method in Case of Probabilistic Loading Conditions［J］．Structural and Multidisciplinary Optimization，2007，35（2）：147－162．

［6］ Chiandussi G，Codegone M，F(xiàn)errero S．Topology Optimization with Optimality Criteria and Transmissible Loads［J］．Computers and Mathematics with Applications，2009，57：772－788．

［7］ Rietz A．Sufficiency of a Finite Exponent in SIMP（Power Law）Method［J］．Structural and Multidiscipline Optimization，2001，21：159－163．

［8］ Stolpe M，Svanberg K．An Alternative Interpolation Scheme for Minimum Compliance Topology Optimization［J］．Structural and Multidiscipline Optimization，2001，22：116－124．

［9］ 羅震，陳立平，黃玉盈，等．基于RAMP密度－剛度插值格式的結(jié)構(gòu)拓撲優(yōu)化［J］．計算力學(xué)學(xué)報，2005，22（5）：585－590．

［10］ Frecker M，Ananthasuresh G K，Nishiwaki S，et al．Topological Synthesis of Compliant Mechanisms Using Multi－criteria Optimization［J］．Mech．Des．Trans．ASME，1997，119（2）：238－245