彭細(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]提出的移動漸近線法(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(可移動變形組件)方法[26]、CBS(close B-splines)方法[27]等,可以采用基于非線性規(guī)劃理論的各種求解算法.如果采用數(shù)學(xué)規(guī)劃對偶理論,借助原-對偶變量的顯式關(guān)系,將大規(guī)模的原問題轉(zhuǎn)化為小規(guī)模的對偶問題求解,可以大幅度地提高算法的效率[19-21].然而,對偶目標(biāo)函數(shù)卻是一個含參數(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 個具體的成果: (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)化的求解中,都將會有進(jìn)一步的應(yīng)用.

1 普適的可分離凸規(guī)劃轉(zhuǎn)換為DP-EM 解法

可分離規(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時,將式(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)造,此時的式(17)中的v清零,取 λ(0)=λ(v+1),這里的 λ(v+1)是前一次的最優(yōu)解,開始新一輪的對偶規(guī)劃尋優(yōu),式(3)是一個迭代的尋優(yōu)序列,直到得到新一輪的最優(yōu)解為止.然后,再求出新的原變量.

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

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

1.1 對偶變量迭代逼近法

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

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

1.2 指數(shù)函數(shù)形式的解法

若約束函數(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)系,避免了非線性方程的求解

1.3 冪函數(shù)形式的解法

若約束函數(shù)存在下述形式

其中aji和bi為實(shí)數(shù)值.

同上一情況類似,可以推導(dǎo)出

為了求出原-對偶變量的函數(shù)關(guān)系,假定有

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

上式開方時必須注意實(shí)根的存在.

1.4 基于變換的精確解法

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

若

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

引入變換及其逆變換

代入式(30)得

仿前面的推導(dǎo),得如下結(jié)果

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

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

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

2 算例

本研究團(tuán)隊(duì)以往用對偶規(guī)劃求解結(jié)構(gòu)優(yōu)化問題時,多遇到冪函數(shù)的情況,故多按第3 種情況進(jìn)行了處理,當(dāng)時的推導(dǎo)均為個案進(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 種情況處理.

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

取

計(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

2.2 算例2: 位移約束下連續(xù)體拓?fù)鋬?yōu)化問題[19-21,28]

本算例以多工況位移約束下結(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 個正方形單元.初始位移為工況1 沿集中載荷F1方向位移為-5.156,工況2 沿集中載荷F2方向位移為5.156 (位移方向豎直向上為正).位移約束為:工況1 時約束集中載荷F1方向位移大于等于-15,工況2 時約束集中載荷F2方向位移小于等于15.收斂精度取0.01.

圖2 算例2 結(jié)構(gòu)力學(xué)簡圖Fig.2 Mechanical diagram of example 2

使用MMA 算法經(jīng)過80 次迭代收斂,歷時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 次迭代收斂,歷時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 算法的用時是DP-EM 算法的3.789 倍.結(jié)果顯示DP-EM 算法具有更高的求解效率.

圖3 算例2 最優(yōu)拓?fù)銯ig.3 Optimized toplolgies of example 2

2.3 算例3: 應(yīng)力約束下連續(xù)體拓?fù)鋬?yōu)化問題[29-31]

算例3 的基結(jié)構(gòu)及載荷與算例2 相同,只是為了避免集中力作用導(dǎo)致的應(yīng)力集中,將集中載荷分散布置在相鄰的3 個節(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,兩個工況下的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ù)涫遣幌嗤?

2.4 算例4: 疲勞約束下連續(xù)體拓?fù)鋬?yōu)化問題[30-31]

如圖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 個正方形單元.基結(jié)構(gòu)的左邊界固定,正弦形式循環(huán)載荷F=7500 N且均值為零,相位角為零,作用于基結(jié)構(gòu)右邊界的中心點(diǎn),為了避免應(yīng)力集中將載荷施加在3 個節(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)時疲勞壽命分布Fig.9 Fatigue life distribution of of the optimized topology of example 4

2.5 算例5: 破損-安全拓?fù)鋬?yōu)化問題[32]

如圖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)楸Ａ舨话l(fā)生破損的區(qū)域(如圖11 所示方格子圖案填充部分).有限元網(wǎng)格為160×100 個正方形單元,位移約束條件為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)時約束點(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)時所有各破損工況位移約束點(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)時各破損工況對應(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)時約束位移值Fig.13 Optimized toplolgies with failure regions and their constraint displacements of example 5

3 結(jié)論

我們團(tuán)隊(duì)以往經(jīng)常把對偶規(guī)劃用于結(jié)構(gòu)優(yōu)化的建模和求解,尤其用于結(jié)構(gòu)拓?fù)鋬?yōu)化問題,以提高求解效率,不過那時建立對偶規(guī)劃的海賽陣,均為對個案的推導(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)化模型給出個案的求解公式及算法,本文針對一類變量可分離凸規(guī)劃模型所闡述的普適處理方法,對文獻(xiàn)[28]的DP-EM 方法進(jìn)行了推廣,從原來的只處理冪函數(shù)形式推廣為可處理任意函數(shù)形式,因而保證今后懲罰函數(shù)或過濾函數(shù)采取不同函數(shù)形式時,優(yōu)化模型的求解有一套普遍適用的可靠理論和算法.本文所提出方法相比廣泛使用的MMA 方法具有更高的求解效率.這說明理論高度的提升會導(dǎo)致更有效和廣泛的應(yīng)用前景.相信本文的工作也會對于同行有借鑒作用.

附錄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ù),單元重量及單元剛度用過濾函數(shù)識別

其中wi及ki為單元重量及單元剛度陣,及為單元固有重量及單元固有剛度陣.fw(ti)=(ti)α及fk(ti)=(ti)β分別為單元重量及單元剛度的冪函數(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)行識別

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

將N個約束應(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è)置的拓?fù)渥兞肯孪拗?

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