彭細(xì)榮 隋允康
* (湖南城市學(xué)院土木工程學(xué)院,湖南益陽 413000)
? (北京工業(yè)大學(xué)材料與制造學(xué)部,北京 100022)
Fleury 首先將非線性對偶規(guī)劃引到結(jié)構(gòu)優(yōu)化領(lǐng)域,進(jìn)行建立模型和求解算法的研究[1-3].錢令希團(tuán)隊(duì)用對偶序列二次規(guī)劃算法求解結(jié)構(gòu)優(yōu)化問題[4].結(jié)構(gòu)優(yōu)化中經(jīng)典的優(yōu)化求解算法如CONLIN (convex linearization method)[5]是基于對偶算法.Svanberg[6]提出的移動(dòng)漸近線法(MMA)也是基于對偶算法.Beekers 等[7]和Hoppe 等[8]均是較早研究了利用原-對偶算法解決結(jié)構(gòu)優(yōu)化問題.隋允康等[9-10]結(jié)合累積迭代信息策略,研究了對偶算法在結(jié)構(gòu)優(yōu)化中的方法和應(yīng)用.吳京等[11]和金鵬等[12]用對偶規(guī)劃研究了鋼管混凝土拱結(jié)構(gòu)優(yōu)化設(shè)計(jì)問題.Xuan 等[13]在優(yōu)化層應(yīng)用對偶內(nèi)點(diǎn)法的線性規(guī)劃和在分析層應(yīng)用二次規(guī)劃法求解非光滑結(jié)構(gòu)優(yōu)化問題.Etman 等[14]研究了利用對偶理論的近似對角二次規(guī)劃法.Cronje 等[15]探討了序列對偶近似規(guī)劃中約束的等效性.Liang等[16]通過序列整數(shù)規(guī)劃和典范對偶理論求解結(jié)構(gòu)拓?fù)鋬?yōu)化問題.
多數(shù)結(jié)構(gòu)拓?fù)鋬?yōu)化方法對應(yīng)的模型歸結(jié)為非線性規(guī)劃,如均勻化方法[17],變密度方法[18]、ICM(獨(dú)立連續(xù)映射)方法[19-21]、水平集方法[22-23]、拓?fù)鋵?dǎo)數(shù)法[24-25]、MMC(可移動(dòng)變形組件)方法[26]、CBS(close B-splines)方法[27]等,可以采用基于非線性規(guī)劃理論的各種求解算法.如果采用數(shù)學(xué)規(guī)劃對偶理論,借助原-對偶變量的顯式關(guān)系,將大規(guī)模的原問題轉(zhuǎn)化為小規(guī)模的對偶問題求解,可以大幅度地提高算法的效率[19-21].然而,對偶目標(biāo)函數(shù)卻是一個(gè)含參數(shù)的極小化問題,難以求解出對偶目標(biāo)函數(shù)的顯式表達(dá)式,只得轉(zhuǎn)而求其1 階導(dǎo)數(shù)和2 階導(dǎo)數(shù),構(gòu)造二次規(guī)劃逼近對偶規(guī)劃,迭代求解,最后由對偶問題最大解求出原問題的最小解.
隋允康等[28]從優(yōu)化所建立的數(shù)學(xué)規(guī)劃模型的特點(diǎn)出發(fā),針對一大類廣泛遇見的變量可分離的凸規(guī)劃問題,突破了只是停留在對偶目標(biāo)二階近似的定勢,推導(dǎo)得出了顯式的對偶目標(biāo)函數(shù).在這一研究的基礎(chǔ)上,提出了便捷求解的對偶規(guī)劃顯式模型(dual programming-explicit model,DP-EM)方法,并將其應(yīng)用于ICM 方法[19-21]求解連續(xù)體結(jié)構(gòu)拓?fù)鋬?yōu)化問題,將DP-EM 方法與基于DSQP 算法及MMA 算法[6]的方法進(jìn)行計(jì)算效率對比,結(jié)果顯示DP-EM 方法具有更高的求解效率.
本文的工作旨在進(jìn)一步提升DP-EM 方法的境界,隋允康等[28]提出的DP-EM 方法僅解決了優(yōu)化模型中約束和目標(biāo)函數(shù)的每一項(xiàng)均為冪函數(shù)形式的情況(對應(yīng)于本文所列舉的情況3: 冪函數(shù)形式的解法),而本文把常見的一類變量可分離的凸規(guī)劃模型抽象為普適的可分離規(guī)劃列式,在需要滿足的一些條件下,轉(zhuǎn)換為DP-EM 解法.相關(guān)工作有4 個(gè)具體的成果: (1)對偶變量迭代逼近法;(2)指數(shù)函數(shù)形式的解法;(3)冪函數(shù)形式的解法;(4)基于變換的精確解法.這一工作有深度與廣度兩方面的意義: (1)深度方面,加深了結(jié)構(gòu)優(yōu)化對偶解法的研究,相關(guān)的結(jié)論在結(jié)構(gòu)拓?fù)鋬?yōu)化里得到了驗(yàn)證;(2)廣度方面,對數(shù)學(xué)規(guī)劃對偶理論的發(fā)展做出了新貢獻(xiàn).預(yù)期今后無論在數(shù)學(xué)規(guī)劃尋優(yōu)算法和結(jié)構(gòu)拓?fù)鋬?yōu)化的求解中,都將會(huì)有進(jìn)一步的應(yīng)用.
可分離規(guī)劃通??梢员磉_(dá)為
其中fi和gji為任意的單變量函數(shù),即目標(biāo)函數(shù)與約束條件屬于可分離函數(shù)形式.
式(1)的拉格朗日函數(shù)為
其中 λj(j=1,2,···,M)為拉格朗日乘子.
原問題(1)的對偶問題為
式(4)的庫恩-塔克條件為
式(5)表明存在如下的函數(shù)關(guān)系
式(6)表達(dá)了抽象的原-對偶設(shè)計(jì)變量關(guān)系,當(dāng)對偶規(guī)劃的最優(yōu)解求出來之后,需根據(jù)具體的原-對偶設(shè)計(jì)變量關(guān)系式(5)由對偶規(guī)劃的最優(yōu)解求出來對應(yīng)的原規(guī)劃最優(yōu)解.可以按如下兩步求解式(5).
先得到如下非線性方程的解
由式(7)的解得到
為了給推導(dǎo)對偶規(guī)劃的海賽陣做準(zhǔn)備,先將式(6)代入式(5)得
對偶規(guī)劃有如下結(jié)論[2-3]
由式(10)可列出對偶目標(biāo)函數(shù)的梯度向量,而其海賽陣的元素則可在式(10)繼續(xù)求偏導(dǎo)數(shù)給出
由式(12)得
式(13)是i∈Ia的結(jié)果;若i?Ia,則有
于是,若i?Ia時(shí),將式(13)和式(15)代入式(11),就得到了
有了式(10)和式(16),就可以構(gòu)造對偶規(guī)劃的近似2 階顯式
其中 λ(v)取上一次的最優(yōu)解 λ?作為本次的初始點(diǎn),以便構(gòu)造第v次的對偶二次規(guī)劃模型.
當(dāng)?shù)趘次尋優(yōu)收斂后,需要代入式(6)求出對應(yīng)的原變量最優(yōu)解.不過式(6)只是理論上的表達(dá),具體將第v+1 次的最優(yōu)解 λ(v+1)代入式(7)求出xs,然后代入式(8)得到第 μ +1 次的原變量x(μ+1),這里的μ與v不同,它表示結(jié)構(gòu)重分析的指標(biāo).經(jīng)過在x(μ+1)的結(jié)構(gòu)重分析之后,重新構(gòu)造新的對偶規(guī)劃式(17),其中梯度向量和海賽陣分別按式(10)和式(16)構(gòu)造,此時(shí)的式(17)中的v清零,取 λ(0)=λ(v+1),這里的 λ(v+1)是前一次的最優(yōu)解,開始新一輪的對偶規(guī)劃尋優(yōu),式(3)是一個(gè)迭代的尋優(yōu)序列,直到得到新一輪的最優(yōu)解為止.然后,再求出新的原變量.
值得探討一下式(7)的算法實(shí)現(xiàn).在求解式(7)時(shí),于具體的 λ 值,除了可按變量xi非線性方程式求解,常常會(huì)遇到顯式解,更為方便.
由于式(16)中含有 λj(j=1,2,···,M),因而需對此進(jìn)行處理,使對偶二階導(dǎo)數(shù)不顯含對偶變量.下面給出4 種處理方法,分別是: (1)對偶變量迭代逼近法;(2)指數(shù)函數(shù)形式的解法;(3)冪函數(shù)形式的解法;(4)基于變換的精確解法.
在迭代尋優(yōu)過程中,取式(13)中的fk(λj)=(ti)β為上一次的變量,于是海賽陣中就不顯含本次的對偶變量了.
至于由對偶規(guī)劃的最優(yōu)解按式(7)和式(8)求出對應(yīng)的原變量,一般要求解N個(gè)變量xi各自的非線性方程解,有時(shí)存在原-對偶變量之間的顯函數(shù)關(guān)系可以利用,從而勿需求解非線性方程.
若約束函數(shù)存在下述形式
其中aji和bi為實(shí)數(shù)值.
將式(18)代入式(16)得
為了將式(19)中含 λj(j=1,···,M)的項(xiàng)消除,現(xiàn)將式(18)代入到式(5),得
將式(20)代入式(19)得
為了求出原-對偶變量的函數(shù)關(guān)系,假定有
將式(22)代入到式(20)得
若上式右邊大于0 能夠保證,則得到原-對偶變量的如下顯式關(guān)系,避免了非線性方程的求解
若約束函數(shù)存在下述形式
其中aji和bi為實(shí)數(shù)值.
同上一情況類似,可以推導(dǎo)出
為了求出原-對偶變量的函數(shù)關(guān)系,假定有
將式(25)和式(27)代入式(5),得到原-對偶變量的如下顯式關(guān)系
上式開方時(shí)必須注意實(shí)根的存在.
第2 和第3 種處理針對兩種函數(shù)形式的情況,能否給出約束函數(shù)更一般情況的處理?本文就如下情況予以回答.
若
滿足凸函數(shù)條件,則將式(29)代入式(1)得
引入變換及其逆變換
代入式(30)得
仿前面的推導(dǎo),得如下結(jié)果
其中表達(dá)原-對偶變量關(guān)系的式(36)并不需要解非線性方程,記
由式(37)解得后,可得規(guī)劃式(30)的解
以上的4 種處理盡量不采用第1 種,原因是多了一個(gè)對 λ 的迭代過程.當(dāng)分離的每一項(xiàng)函數(shù)形式為指數(shù)函數(shù)或冪函數(shù)時(shí)分別按第2 種和第3 種處理,一般情況下可采用第4 種處理.
本研究團(tuán)隊(duì)以往用對偶規(guī)劃求解結(jié)構(gòu)優(yōu)化問題時(shí),多遇到冪函數(shù)的情況,故多按第3 種情況進(jìn)行了處理,當(dāng)時(shí)的推導(dǎo)均為個(gè)案進(jìn)行,現(xiàn)以算例2 至算例5 中4 類不同結(jié)構(gòu)拓?fù)鋬?yōu)化問題為例,按本文的結(jié)果進(jìn)行了驗(yàn)證.算例1 為特別設(shè)計(jì)的簡單數(shù)學(xué)算例,以比較按第1 種情況和第4 種情況處理的效率,其它算例均是按第3 種情況處理.
取
計(jì)算如下可分離非線性優(yōu)化問題
取初值 (x1,x2)=(0.01,0.01),按第1 種情況處理,迭代12 次,目標(biāo)值1.922 6,最優(yōu)點(diǎn)(x1,x2)?=(0.713 7,0.433 1).按第4 種情況處理,迭代9 次,目標(biāo)值最優(yōu)點(diǎn)值相同.第4 種處理方式的迭代次數(shù)更少.算例1 的目標(biāo)函數(shù)等值線及約束可行域如圖1 所示,星號點(diǎn)為最優(yōu)點(diǎn).
圖1 算例1 的目標(biāo)函數(shù)等值線及可行域Fig.1 Contour lines of the objective function and feasible region of example 1
本算例以多工況位移約束下結(jié)構(gòu)重量極小化結(jié)構(gòu)拓?fù)鋬?yōu)化問題為例,比較MMA 算法[6]與本文提出的DP-EM 算法的計(jì)算效率.算例數(shù)據(jù)如下.
如圖2 所示,基結(jié)構(gòu)尺寸為240×120×1 的平面體,材料彈性模量為E=1,泊松比為0.3.左右邊界固定,工況1: 集中載荷F1=-1 豎直向下作用于上邊界中點(diǎn);工況2: 集中載荷F2=1 豎直向上作用于下邊界中點(diǎn)(力以豎直向上為正).計(jì)算網(wǎng)格為240×120 個(gè)正方形單元.初始位移為工況1 沿集中載荷F1方向位移為-5.156,工況2 沿集中載荷F2方向位移為5.156 (位移方向豎直向上為正).位移約束為:工況1 時(shí)約束集中載荷F1方向位移大于等于-15,工況2 時(shí)約束集中載荷F2方向位移小于等于15.收斂精度取0.01.
圖2 算例2 結(jié)構(gòu)力學(xué)簡圖Fig.2 Mechanical diagram of example 2
使用MMA 算法經(jīng)過80 次迭代收斂,歷時(shí)247.904 s,最優(yōu)點(diǎn)結(jié)構(gòu)體積比為24.88%,上下結(jié)構(gòu)邊中點(diǎn)處約束點(diǎn)位移分別為-15.001 和15.001,最優(yōu)拓?fù)鋱D形如圖3(a)所示.使用DP-EM 算法經(jīng)過75 次迭代收斂,歷時(shí)65.435 s,最優(yōu)點(diǎn)結(jié)構(gòu)體積比為23.76%,上下結(jié)構(gòu)邊中點(diǎn)處約束點(diǎn)位移分別為-14.999 和14.999,最優(yōu)拓?fù)鋱D形如圖3(b)所示.MMA 算法的用時(shí)是DP-EM 算法的3.789 倍.結(jié)果顯示DP-EM 算法具有更高的求解效率.
圖3 算例2 最優(yōu)拓?fù)銯ig.3 Optimized toplolgies of example 2
算例3 的基結(jié)構(gòu)及載荷與算例2 相同,只是為了避免集中力作用導(dǎo)致的應(yīng)力集中,將集中載荷分散布置在相鄰的3 個(gè)節(jié)點(diǎn)上.Mises 應(yīng)力約束小于或等于0.3.
取用本文算法,迭代66 次收斂,最優(yōu)點(diǎn)體積比15.42%,最優(yōu)拓?fù)鋱D形如圖4 所示.各工況的最大Mises 應(yīng)力為0.295 7,兩個(gè)工況下的Mises 應(yīng)力分布圖如圖5 所示,滿足應(yīng)力約束條件.
圖4 算例3 最優(yōu)拓?fù)銯ig.4 Optimized toplolgy of example 3
圖5 算例3 Mises 應(yīng)力分布Fig.5 Mises stress distributions of example 3
比較算例2 可知,位移約束與應(yīng)力約束下的結(jié)構(gòu)拓?fù)鋬?yōu)化得到的最優(yōu)拓?fù)涫遣幌嗤?
如圖6 所示,材料彈性模量為E=2.1×105MPa,泊松比 ν=0.3,材料對應(yīng)于循環(huán)次數(shù)1.0×106的疲勞極限為 σs=250 MPa,基結(jié)構(gòu)為300 mm×100 mm×10 mm 的平面體,劃分為300×100 個(gè)正方形單元.基結(jié)構(gòu)的左邊界固定,正弦形式循環(huán)載荷F=7500 N且均值為零,相位角為零,作用于基結(jié)構(gòu)右邊界的中心點(diǎn),為了避免應(yīng)力集中將載荷施加在3 個(gè)節(jié)點(diǎn)上.疲勞壽命約束為大于或等于 1.0×105次.巴士昆公式σ=AL-m中,σ 為疲勞循環(huán)下的材料應(yīng)力,L為對應(yīng)的疲勞壽命,A與m為材料給定的常數(shù)量,此例中取值為m=0.10,A=995.3 MPa,S-N曲線如圖7 所示.
圖6 懸臂結(jié)構(gòu)的分析優(yōu)化模型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
取用本文算法,經(jīng)過62 次迭代收斂,最優(yōu)體積比為39.73%,最大疲勞壽命為99 163 次.得到的最優(yōu)拓?fù)淙鐖D8 所示,對應(yīng)的疲勞壽命以10 為底的對數(shù)分布如圖9 所示.
圖8 算例4 最優(yōu)拓?fù)銯ig.8 Optimized topology of example 4
圖9 算例4 最優(yōu)時(shí)疲勞壽命分布Fig.9 Fatigue life distribution of of the optimized topology of example 4
如圖10 所示,基結(jié)構(gòu)為160×100 的平面體,厚度為1,材料彈性模量為1,泊松比為0.3.左邊界固定,一集中載荷F=1 作用于右邊界中心位置.采用25×25 的正方形作為結(jié)構(gòu)局部破損模式,如圖11中白色和灰色正方形所示,采用無縫平鋪的結(jié)構(gòu)破損狀況預(yù)估分布模式,右邊界附近區(qū)域?yàn)楸A舨话l(fā)生破損的區(qū)域(如圖11 所示方格子圖案填充部分).有限元網(wǎng)格為160×100 個(gè)正方形單元,位移約束條件為A點(diǎn)的豎直向下位移值小于120.
圖10 算例5 力學(xué)模型Fig.10 Mechanical model of example 5
圖11 算例5 的破損區(qū)域及其分布Fig.11 Fail regions and their distribution of example 5
不考慮破損-安全,按傳統(tǒng)安全理念進(jìn)行拓?fù)鋬?yōu)化,得到的最優(yōu)點(diǎn)體積比為18.40%,最優(yōu)拓?fù)淙鐖D12(a)所示,最優(yōu)點(diǎn)時(shí)約束點(diǎn)位移值為-119.9,滿足位移約束條件.考慮破損-安全的結(jié)構(gòu)拓?fù)鋬?yōu)化最優(yōu)點(diǎn)體積比為38.12%,最優(yōu)拓?fù)淙鐖D12(b)所示.可以看出,考慮破損-安全之后得到的結(jié)構(gòu)體積比更大,結(jié)構(gòu)拓?fù)涓鼜?fù)雜,但結(jié)構(gòu)可滿足所有破損情況下的安全.圖13 為部分破損工況下的結(jié)構(gòu)拓?fù)鋱D及位移,表1 為最優(yōu)點(diǎn)時(shí)所有各破損工況位移約束點(diǎn)的位移值,表中各行列數(shù)據(jù)均與圖11 中的破損工況位置一一對應(yīng),例如表中的第1 行第1 列數(shù)據(jù)(-120.0)對應(yīng)于圖11 中第1 行第1 列破損位置發(fā)生破損,如圖13(a)所示.可以看出有些破損工況是關(guān)鍵約束,有些破損工況不是關(guān)鍵約束,但所有破損工況均滿足位移約束條件.
表1 算例5 最優(yōu)點(diǎn)時(shí)各破損工況對應(yīng)的約束點(diǎn)位移值Table 1 Displacements at constraint point for all failure cases at optimization of example 5
圖12 算例5 最優(yōu)拓?fù)銯ig.12 Optimized toplolgies of example 5
圖13 算例5 帶破損區(qū)域的最優(yōu)拓?fù)浼皩?yīng)最優(yōu)點(diǎn)時(shí)約束位移值Fig.13 Optimized toplolgies with failure regions and their constraint displacements of example 5
我們團(tuán)隊(duì)以往經(jīng)常把對偶規(guī)劃用于結(jié)構(gòu)優(yōu)化的建模和求解,尤其用于結(jié)構(gòu)拓?fù)鋬?yōu)化問題,以提高求解效率,不過那時(shí)建立對偶規(guī)劃的海賽陣,均為對個(gè)案的推導(dǎo).
本文闡述了普適的處理方法,針對一類變量可分離凸規(guī)劃模型,給出了4 種處理方法: (1)對偶變量迭代逼近法;(2)指數(shù)函數(shù)形式的解法;(3)冪函數(shù)形式的解法;(4)基于變換的精確解法.其中方法(1)和方法(4)是普適的方法,由于方法4 不需要采用迭代逼近的近似求解,比方法1 的求解效率更高.而方法2 和方法3 則示例了當(dāng)函數(shù)形式給定后,依據(jù)具體的函數(shù)形式,同樣可以給出直接解法,而不需要如方法1 中采取迭代解法.
不同于以往變密度法或ICM 方法中針對某一具體懲罰函數(shù)或過濾函數(shù)形式所建立的優(yōu)化模型給出個(gè)案的求解公式及算法,本文針對一類變量可分離凸規(guī)劃模型所闡述的普適處理方法,對文獻(xiàn)[28]的DP-EM 方法進(jìn)行了推廣,從原來的只處理冪函數(shù)形式推廣為可處理任意函數(shù)形式,因而保證今后懲罰函數(shù)或過濾函數(shù)采取不同函數(shù)形式時(shí),優(yōu)化模型的求解有一套普遍適用的可靠理論和算法.本文所提出方法相比廣泛使用的MMA 方法具有更高的求解效率.這說明理論高度的提升會(huì)導(dǎo)致更有效和廣泛的應(yīng)用前景.相信本文的工作也會(huì)對于同行有借鑒作用.
附錄A
以下列舉基于ICM 方法建立的位移、應(yīng)力、疲勞等約束下結(jié)構(gòu)拓?fù)鋬?yōu)化模型,以及考慮破損-安全的結(jié)構(gòu)拓?fù)鋬?yōu)化模型.
(1)位移約束連續(xù)體結(jié)構(gòu)拓?fù)鋬?yōu)化[19-21,28]
采用冪函數(shù)形式的過濾函數(shù),單元重量及單元?jiǎng)偠扔眠^濾函數(shù)識(shí)別
其中wi及ki為單元重量及單元?jiǎng)偠汝?及為單元固有重量及單元固有剛度陣.fw(ti)=(ti)α及fk(ti)=(ti)β分別為單元重量及單元?jiǎng)偠鹊膬绾瘮?shù)形式過濾函數(shù).
結(jié)構(gòu)總重量表達(dá)為
其中N表示單元拓?fù)湓O(shè)計(jì)變量的總數(shù).
位移利用莫爾定理表達(dá)為
用ICM 法建立的拓?fù)鋬?yōu)化模型表示為
式中t=(t1,t2,···,tN)T是單元拓?fù)湓O(shè)計(jì)變量向量表示第j號位移約束的上限.
(2)應(yīng)力或疲勞局部性能約束連續(xù)體結(jié)構(gòu)拓?fù)鋬?yōu)化[30-31]
對結(jié)構(gòu)重量(或體積)、結(jié)構(gòu)局部性能如應(yīng)力或疲勞壽命等,均用過濾函數(shù)進(jìn)行識(shí)別
式中,fw(ti)及fψ(ti)分別為結(jié)構(gòu)重量(或體積)及結(jié)構(gòu)性能如應(yīng)力或疲勞壽命等的過濾函數(shù),及分別為結(jié)構(gòu)固有的重量(或體積)及結(jié)構(gòu)性能如應(yīng)力或疲勞壽命允許值等.
將N個(gè)約束應(yīng)用K-S 函數(shù)進(jìn)行集成,得到基于K-S 函數(shù)集成的局部性能約束優(yōu)化模型
(3)考慮破損-安全的連續(xù)體結(jié)構(gòu)拓?fù)鋬?yōu)化[32]
應(yīng)用ICM 方法,建立的結(jié)構(gòu)拓?fù)鋬?yōu)化模型為
其中,t為拓?fù)湓O(shè)計(jì)變量向量;V(t,Φ)為結(jié)構(gòu)失效狀況集合Φ中所有結(jié)構(gòu)失效狀況中體積最大值;ujls(t)為結(jié)構(gòu)失效狀況集合 Φ 中第s號結(jié)構(gòu)失效狀況 ?s在第l的載荷工況下的第j號位移約束函數(shù),為第j號位移約束值;J為位移約束總數(shù),L為載荷工況總數(shù);為防止有限元分析時(shí)總剛度矩陣奇異而設(shè)置的拓?fù)渥兞肯孪拗?
結(jié)構(gòu)總體積及位移約束的顯式化與1 中的式(A1)~式(A4)類似.在此不再重復(fù).