基于比例微分優化準則的拓撲優化方法

趙龍彪 高 亮 陳志敏 邱浩波

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

0 引言

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

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

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

1 結構拓撲優化設計基礎

1.1 結構拓撲優化模型

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

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

SIMP材料模型的數學表達形式：

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

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

以結構拓撲優化中比較典型的最小柔度問題為例，在SIMP材料插值方法的基礎上，拓撲優化模型為

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

1.2 優化準則法

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

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

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

式中，ζ為阻尼因子（0＜ζ＜1）；k為迭代步數；Vj為體積比。得到一個改進的設計x（k＋1）。式（6）中常常引入一些經驗系數來調整優化過程的收斂性和穩定性，如步長因子、阻尼因子等。

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

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

（2）求解通用性較差，對于復雜的目標函數，不容易構造設計變量的迭代公式，搜索過程中的指導信息只有當前結果和目標結果的偏差，所以容易出現震蕩和難以收斂的現象。

（3）一般比較適用于單目標函數、單約束條件下的優化問題求解，多目標問題需要推導不同的優化準則，難以適用于復雜問題的求解。

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

2 PDOC方法及算法實現

2.1 PD控制基礎

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

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

（2）微分控制（differential control）。在微分控制中，控制器的輸出與輸入誤差信號的微分（即誤差的變化率）成正比關系。控制器中，微分項的引入，能夠預測誤差變化的趨勢。

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

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

2.2 PDOC方法介紹

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

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

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

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

結合優化準則迭代方法和微分控制的思想，可以構造比例微分迭代公式如下：

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

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

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

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

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

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

2.3 PDOC算法實現

圖1為基于PDOC法拓撲優化設計的算法流程圖，可以采用 MATLAB實現算法，步驟如下：

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

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

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

（4）根據上一步迭代更新的各單元材料密度，采用PDOC法計算各單元新的材料密度和結構柔度。

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

（6）輸出結果，結束。

圖1 PDOC法拓撲優化求解流程圖

3 實例分析

3.1 懸臂梁剛度優化

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

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

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

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

圖3 網格劃分120×40懸臂梁的解

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

3.2 簡支梁剛度優化

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

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

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

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

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

3.3 柔性位移反向器

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

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

圖8 位移反相器的設計域示意圖

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

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

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

4 結語

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

［1］ 羅震，陳立平，黃玉盈，等．連續體結構的拓撲優化設計［J］．力學進展，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，Ferrero 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密度－剛度插值格式的結構拓撲優化［J］．計算力學學報，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