化可分離凸規劃為對偶規劃顯式模型的普適解法研究比較1)

彭細榮 隋允康

* (湖南城市學院土木工程學院,湖南益陽 413000)

? (北京工業大學材料與制造學部,北京 100022)

引言

Fleury 首先將非線性對偶規劃引到結構優化領域,進行建立模型和求解算法的研究[1-3].錢令希團隊用對偶序列二次規劃算法求解結構優化問題[4].結構優化中經典的優化求解算法如CONLIN (convex linearization method)[5]是基于對偶算法.Svanberg[6]提出的移動漸近線法(MMA)也是基于對偶算法.Beekers 等[7]和Hoppe 等[8]均是較早研究了利用原-對偶算法解決結構優化問題.隋允康等[9-10]結合累積迭代信息策略,研究了對偶算法在結構優化中的方法和應用.吳京等[11]和金鵬等[12]用對偶規劃研究了鋼管混凝土拱結構優化設計問題.Xuan 等[13]在優化層應用對偶內點法的線性規劃和在分析層應用二次規劃法求解非光滑結構優化問題.Etman 等[14]研究了利用對偶理論的近似對角二次規劃法.Cronje 等[15]探討了序列對偶近似規劃中約束的等效性.Liang等[16]通過序列整數規劃和典范對偶理論求解結構拓撲優化問題.

多數結構拓撲優化方法對應的模型歸結為非線性規劃,如均勻化方法[17],變密度方法[18]、ICM(獨立連續映射)方法[19-21]、水平集方法[22-23]、拓撲導數法[24-25]、MMC(可移動變形組件)方法[26]、CBS(close B-splines)方法[27]等,可以采用基于非線性規劃理論的各種求解算法.如果采用數學規劃對偶理論,借助原-對偶變量的顯式關系,將大規模的原問題轉化為小規模的對偶問題求解,可以大幅度地提高算法的效率[19-21].然而,對偶目標函數卻是一個含參數的極小化問題,難以求解出對偶目標函數的顯式表達式,只得轉而求其1 階導數和2 階導數,構造二次規劃逼近對偶規劃,迭代求解,最后由對偶問題最大解求出原問題的最小解.

隋允康等[28]從優化所建立的數學規劃模型的特點出發,針對一大類廣泛遇見的變量可分離的凸規劃問題,突破了只是停留在對偶目標二階近似的定勢,推導得出了顯式的對偶目標函數.在這一研究的基礎上,提出了便捷求解的對偶規劃顯式模型(dual programming-explicit model,DP-EM)方法,并將其應用于ICM 方法[19-21]求解連續體結構拓撲優化問題,將DP-EM 方法與基于DSQP 算法及MMA 算法[6]的方法進行計算效率對比,結果顯示DP-EM 方法具有更高的求解效率.

本文的工作旨在進一步提升DP-EM 方法的境界,隋允康等[28]提出的DP-EM 方法僅解決了優化模型中約束和目標函數的每一項均為冪函數形式的情況(對應于本文所列舉的情況3: 冪函數形式的解法),而本文把常見的一類變量可分離的凸規劃模型抽象為普適的可分離規劃列式,在需要滿足的一些條件下,轉換為DP-EM 解法.相關工作有4 個具體的成果: (1)對偶變量迭代逼近法;(2)指數函數形式的解法;(3)冪函數形式的解法;(4)基于變換的精確解法.這一工作有深度與廣度兩方面的意義: (1)深度方面,加深了結構優化對偶解法的研究,相關的結論在結構拓撲優化里得到了驗證;(2)廣度方面,對數學規劃對偶理論的發展做出了新貢獻.預期今后無論在數學規劃尋優算法和結構拓撲優化的求解中,都將會有進一步的應用.

1 普適的可分離凸規劃轉換為DP-EM 解法

可分離規劃通常可以表達為

其中fi和gji為任意的單變量函數,即目標函數與約束條件屬于可分離函數形式.

式(1)的拉格朗日函數為

其中 λj(j=1,2,···,M)為拉格朗日乘子.

原問題(1)的對偶問題為

式(4)的庫恩-塔克條件為

式(5)表明存在如下的函數關系

式(6)表達了抽象的原-對偶設計變量關系,當對偶規劃的最優解求出來之后,需根據具體的原-對偶設計變量關系式(5)由對偶規劃的最優解求出來對應的原規劃最優解.可以按如下兩步求解式(5).

先得到如下非線性方程的解

由式(7)的解得到

為了給推導對偶規劃的海賽陣做準備,先將式(6)代入式(5)得

對偶規劃有如下結論[2-3]

由式(10)可列出對偶目標函數的梯度向量,而其海賽陣的元素則可在式(10)繼續求偏導數給出

由式(12)得

式(13)是i∈Ia的結果;若i?Ia,則有

于是,若i?Ia時,將式(13)和式(15)代入式(11),就得到了

有了式(10)和式(16),就可以構造對偶規劃的近似2 階顯式

其中 λ(v)取上一次的最優解 λ?作為本次的初始點,以便構造第v次的對偶二次規劃模型.

當第v次尋優收斂后,需要代入式(6)求出對應的原變量最優解.不過式(6)只是理論上的表達,具體將第v+1 次的最優解 λ(v+1)代入式(7)求出xs,然后代入式(8)得到第 μ +1 次的原變量x(μ+1),這里的μ與v不同,它表示結構重分析的指標.經過在x(μ+1)的結構重分析之后,重新構造新的對偶規劃式(17),其中梯度向量和海賽陣分別按式(10)和式(16)構造,此時的式(17)中的v清零,取 λ(0)=λ(v+1),這里的 λ(v+1)是前一次的最優解,開始新一輪的對偶規劃尋優,式(3)是一個迭代的尋優序列,直到得到新一輪的最優解為止.然后,再求出新的原變量.

值得探討一下式(7)的算法實現.在求解式(7)時,于具體的 λ 值,除了可按變量xi非線性方程式求解,常常會遇到顯式解,更為方便.

由于式(16)中含有 λj(j=1,2,···,M),因而需對此進行處理,使對偶二階導數不顯含對偶變量.下面給出4 種處理方法,分別是: (1)對偶變量迭代逼近法;(2)指數函數形式的解法;(3)冪函數形式的解法;(4)基于變換的精確解法.

1.1 對偶變量迭代逼近法

在迭代尋優過程中,取式(13)中的fk(λj)=(ti)β為上一次的變量,于是海賽陣中就不顯含本次的對偶變量了.

至于由對偶規劃的最優解按式(7)和式(8)求出對應的原變量,一般要求解N個變量xi各自的非線性方程解,有時存在原-對偶變量之間的顯函數關系可以利用,從而勿需求解非線性方程.

1.2 指數函數形式的解法

若約束函數存在下述形式

其中aji和bi為實數值.

將式(18)代入式(16)得

為了將式(19)中含 λj(j=1,···,M)的項消除,現將式(18)代入到式(5),得

將式(20)代入式(19)得

為了求出原-對偶變量的函數關系,假定有

將式(22)代入到式(20)得

若上式右邊大于0 能夠保證,則得到原-對偶變量的如下顯式關系,避免了非線性方程的求解

1.3 冪函數形式的解法

若約束函數存在下述形式

其中aji和bi為實數值.

同上一情況類似,可以推導出

為了求出原-對偶變量的函數關系,假定有

將式(25)和式(27)代入式(5),得到原-對偶變量的如下顯式關系

上式開方時必須注意實根的存在.

1.4 基于變換的精確解法

第2 和第3 種處理針對兩種函數形式的情況,能否給出約束函數更一般情況的處理？本文就如下情況予以回答.

若

滿足凸函數條件,則將式(29)代入式(1)得

引入變換及其逆變換

代入式(30)得

仿前面的推導,得如下結果

其中表達原-對偶變量關系的式(36)并不需要解非線性方程,記

由式(37)解得后,可得規劃式(30)的解

以上的4 種處理盡量不采用第1 種,原因是多了一個對 λ 的迭代過程.當分離的每一項函數形式為指數函數或冪函數時分別按第2 種和第3 種處理,一般情況下可采用第4 種處理.

2 算例

本研究團隊以往用對偶規劃求解結構優化問題時,多遇到冪函數的情況,故多按第3 種情況進行了處理,當時的推導均為個案進行,現以算例2 至算例5 中4 類不同結構拓撲優化問題為例,按本文的結果進行了驗證.算例1 為特別設計的簡單數學算例,以比較按第1 種情況和第4 種情況處理的效率,其它算例均是按第3 種情況處理.

2.1 算例1: 普適處理的效果比較

取

計算如下可分離非線性優化問題

取初值 (x1,x2)=(0.01,0.01),按第1 種情況處理,迭代12 次,目標值1.922 6,最優點(x1,x2)?=(0.713 7,0.433 1).按第4 種情況處理,迭代9 次,目標值最優點值相同.第4 種處理方式的迭代次數更少.算例1 的目標函數等值線及約束可行域如圖1 所示,星號點為最優點.

圖1 算例1 的目標函數等值線及可行域Fig.1 Contour lines of the objective function and feasible region of example 1

2.2 算例2: 位移約束下連續體拓撲優化問題[19-21,28]

本算例以多工況位移約束下結構重量極小化結構拓撲優化問題為例,比較MMA 算法[6]與本文提出的DP-EM 算法的計算效率.算例數據如下.

如圖2 所示,基結構尺寸為240×120×1 的平面體,材料彈性模量為E=1,泊松比為0.3.左右邊界固定,工況1: 集中載荷F1=-1 豎直向下作用于上邊界中點;工況2: 集中載荷F2=1 豎直向上作用于下邊界中點(力以豎直向上為正).計算網格為240×120 個正方形單元.初始位移為工況1 沿集中載荷F1方向位移為-5.156,工況2 沿集中載荷F2方向位移為5.156 (位移方向豎直向上為正).位移約束為:工況1 時約束集中載荷F1方向位移大于等于-15,工況2 時約束集中載荷F2方向位移小于等于15.收斂精度取0.01.

圖2 算例2 結構力學簡圖Fig.2 Mechanical diagram of example 2

使用MMA 算法經過80 次迭代收斂,歷時247.904 s,最優點結構體積比為24.88%,上下結構邊中點處約束點位移分別為-15.001 和15.001,最優拓撲圖形如圖3(a)所示.使用DP-EM 算法經過75 次迭代收斂,歷時65.435 s,最優點結構體積比為23.76%,上下結構邊中點處約束點位移分別為-14.999 和14.999,最優拓撲圖形如圖3(b)所示.MMA 算法的用時是DP-EM 算法的3.789 倍.結果顯示DP-EM 算法具有更高的求解效率.

圖3 算例2 最優拓撲Fig.3 Optimized toplolgies of example 2

2.3 算例3: 應力約束下連續體拓撲優化問題[29-31]

算例3 的基結構及載荷與算例2 相同,只是為了避免集中力作用導致的應力集中,將集中載荷分散布置在相鄰的3 個節點上.Mises 應力約束小于或等于0.3.

取用本文算法,迭代66 次收斂,最優點體積比15.42%,最優拓撲圖形如圖4 所示.各工況的最大Mises 應力為0.295 7,兩個工況下的Mises 應力分布圖如圖5 所示,滿足應力約束條件.

圖4 算例3 最優拓撲Fig.4 Optimized toplolgy of example 3

圖5 算例3 Mises 應力分布Fig.5 Mises stress distributions of example 3

比較算例2 可知,位移約束與應力約束下的結構拓撲優化得到的最優拓撲是不相同的.

2.4 算例4: 疲勞約束下連續體拓撲優化問題[30-31]

如圖6 所示,材料彈性模量為E=2.1×105MPa,泊松比 ν=0.3,材料對應于循環次數1.0×106的疲勞極限為 σs=250 MPa,基結構為300 mm×100 mm×10 mm 的平面體,劃分為300×100 個正方形單元.基結構的左邊界固定,正弦形式循環載荷F=7500 N且均值為零,相位角為零,作用于基結構右邊界的中心點,為了避免應力集中將載荷施加在3 個節點上.疲勞壽命約束為大于或等于 1.0×105次.巴士昆公式σ=AL-m中,σ 為疲勞循環下的材料應力,L為對應的疲勞壽命,A與m為材料給定的常數量,此例中取值為m=0.10,A=995.3 MPa,S-N曲線如圖7 所示.

圖6 懸臂結構的分析優化模型Fig.6 Model of the cantilevel structure for analysis and optimization

圖7 m=0.10,A=995.3 MPa 材料的S-N 曲線Fig.7 S-N curve for the material with m=0.10,A=995.3 MPa

取用本文算法,經過62 次迭代收斂,最優體積比為39.73%,最大疲勞壽命為99 163 次.得到的最優拓撲如圖8 所示,對應的疲勞壽命以10 為底的對數分布如圖9 所示.

圖8 算例4 最優拓撲Fig.8 Optimized topology of example 4

圖9 算例4 最優時疲勞壽命分布Fig.9 Fatigue life distribution of of the optimized topology of example 4

2.5 算例5: 破損-安全拓撲優化問題[32]

如圖10 所示,基結構為160×100 的平面體,厚度為1,材料彈性模量為1,泊松比為0.3.左邊界固定,一集中載荷F=1 作用于右邊界中心位置.采用25×25 的正方形作為結構局部破損模式,如圖11中白色和灰色正方形所示,采用無縫平鋪的結構破損狀況預估分布模式,右邊界附近區域為保留不發生破損的區域(如圖11 所示方格子圖案填充部分).有限元網格為160×100 個正方形單元,位移約束條件為A點的豎直向下位移值小于120.

圖10 算例5 力學模型Fig.10 Mechanical model of example 5

圖11 算例5 的破損區域及其分布Fig.11 Fail regions and their distribution of example 5

不考慮破損-安全,按傳統安全理念進行拓撲優化,得到的最優點體積比為18.40%,最優拓撲如圖12(a)所示,最優點時約束點位移值為-119.9,滿足位移約束條件.考慮破損-安全的結構拓撲優化最優點體積比為38.12%,最優拓撲如圖12(b)所示.可以看出,考慮破損-安全之后得到的結構體積比更大,結構拓撲更復雜,但結構可滿足所有破損情況下的安全.圖13 為部分破損工況下的結構拓撲圖及位移,表1 為最優點時所有各破損工況位移約束點的位移值,表中各行列數據均與圖11 中的破損工況位置一一對應,例如表中的第1 行第1 列數據(-120.0)對應于圖11 中第1 行第1 列破損位置發生破損,如圖13(a)所示.可以看出有些破損工況是關鍵約束,有些破損工況不是關鍵約束,但所有破損工況均滿足位移約束條件.

表1 算例5 最優點時各破損工況對應的約束點位移值Table 1 Displacements at constraint point for all failure cases at optimization of example 5

圖12 算例5 最優拓撲Fig.12 Optimized toplolgies of example 5

圖13 算例5 帶破損區域的最優拓撲及對應最優點時約束位移值Fig.13 Optimized toplolgies with failure regions and their constraint displacements of example 5

3 結論

我們團隊以往經常把對偶規劃用于結構優化的建模和求解,尤其用于結構拓撲優化問題,以提高求解效率,不過那時建立對偶規劃的海賽陣,均為對個案的推導.

本文闡述了普適的處理方法,針對一類變量可分離凸規劃模型,給出了4 種處理方法: (1)對偶變量迭代逼近法;(2)指數函數形式的解法;(3)冪函數形式的解法;(4)基于變換的精確解法.其中方法(1)和方法(4)是普適的方法,由于方法4 不需要采用迭代逼近的近似求解,比方法1 的求解效率更高.而方法2 和方法3 則示例了當函數形式給定后,依據具體的函數形式,同樣可以給出直接解法,而不需要如方法1 中采取迭代解法.

不同于以往變密度法或ICM 方法中針對某一具體懲罰函數或過濾函數形式所建立的優化模型給出個案的求解公式及算法,本文針對一類變量可分離凸規劃模型所闡述的普適處理方法,對文獻[28]的DP-EM 方法進行了推廣,從原來的只處理冪函數形式推廣為可處理任意函數形式,因而保證今后懲罰函數或過濾函數采取不同函數形式時,優化模型的求解有一套普遍適用的可靠理論和算法.本文所提出方法相比廣泛使用的MMA 方法具有更高的求解效率.這說明理論高度的提升會導致更有效和廣泛的應用前景.相信本文的工作也會對于同行有借鑒作用.

附錄A

以下列舉基于ICM 方法建立的位移、應力、疲勞等約束下結構拓撲優化模型,以及考慮破損-安全的結構拓撲優化模型.

(1)位移約束連續體結構拓撲優化[19-21,28]

采用冪函數形式的過濾函數,單元重量及單元剛度用過濾函數識別

其中wi及ki為單元重量及單元剛度陣,及為單元固有重量及單元固有剛度陣.fw(ti)=(ti)α及fk(ti)=(ti)β分別為單元重量及單元剛度的冪函數形式過濾函數.

結構總重量表達為

其中N表示單元拓撲設計變量的總數.

位移利用莫爾定理表達為

用ICM 法建立的拓撲優化模型表示為

式中t=(t1,t2,···,tN)T是單元拓撲設計變量向量表示第j號位移約束的上限.

(2)應力或疲勞局部性能約束連續體結構拓撲優化[30-31]

對結構重量(或體積)、結構局部性能如應力或疲勞壽命等,均用過濾函數進行識別

式中,fw(ti)及fψ(ti)分別為結構重量(或體積)及結構性能如應力或疲勞壽命等的過濾函數,及分別為結構固有的重量(或體積)及結構性能如應力或疲勞壽命允許值等.

將N個約束應用K-S 函數進行集成,得到基于K-S 函數集成的局部性能約束優化模型

(3)考慮破損-安全的連續體結構拓撲優化[32]

應用ICM 方法,建立的結構拓撲優化模型為

其中,t為拓撲設計變量向量;V(t,Φ)為結構失效狀況集合Φ中所有結構失效狀況中體積最大值;ujls(t)為結構失效狀況集合 Φ 中第s號結構失效狀況 ?s在第l的載荷工況下的第j號位移約束函數,為第j號位移約束值;J為位移約束總數,L為載荷工況總數;為防止有限元分析時總剛度矩陣奇異而設置的拓撲變量下限值.

結構總體積及位移約束的顯式化與1 中的式(A1)～式(A4)類似.在此不再重復.