探討基于改進(jìn)模擬退火算法的復(fù)合材料層合板屈曲優(yōu)化

1、探討基于改進(jìn)模擬退火算法的復(fù)合材料層合板屈曲優(yōu)化0引言復(fù)合材料層合板在航空、汽車、船舶等工程領(lǐng)域被廣泛用于板、殼等結(jié)構(gòu)中，屈曲失穩(wěn)成為其不容忽視的失效形式之一，層合板的屈曲性能表現(xiàn)因此倍受關(guān)注。然而，與傳統(tǒng)金屬板殼相比，各向異性特性、鋪層層數(shù)和鋪層角度的離散化特征，使得層合板的屈曲分析與優(yōu)化更加復(fù)雜困難。較多的層合板屈曲優(yōu)化設(shè)計(jì)以正交各向異性板和0、30、45、60、90等幾種離散鋪層角為研究對象，從優(yōu)化算法和屈曲響應(yīng)分析兩方面研究問題的解決途徑。通過鋪層角度連續(xù)化，利用序列線性規(guī)劃、可行方向法等數(shù)學(xué)規(guī)劃法實(shí)現(xiàn)層合板屈曲優(yōu)化，但連續(xù)解的離散化會導(dǎo)致設(shè)計(jì)結(jié)果非最優(yōu)化或不滿足約束，且優(yōu)化易陷入局部

2、極值。遺傳算法(GA)、模擬退火(SA)、粒子群(PSO)、蟻群(ACO)等隨機(jī)算法能較好地求解離散型層合板屈曲優(yōu)化問題，其目標(biāo)函數(shù)評價(jià)成本遠(yuǎn)高于數(shù)學(xué)規(guī)劃法。Erdal等9采用具有記憶功能的直接搜索模擬退火(DSA)算法開展了層合板屈曲優(yōu)化研究，考慮了角度增量為10、15、30時(shí)不同離散角度對優(yōu)化結(jié)果的影響，但DSA算法易陷入局部解。Karakay等則比較了GA、SA和ACO算法求解層合板屈曲和頻率優(yōu)化問題時(shí)的性能表現(xiàn)。其他如分層優(yōu)化方法、分散搜索算法、進(jìn)化算法及兩層次優(yōu)化方法等也相繼被用于層合板屈曲問題的鋪層順序設(shè)計(jì)。另一方面，準(zhǔn)確求解層合板屈曲響應(yīng)，實(shí)現(xiàn)目標(biāo)函數(shù)的精確計(jì)算，可以提高優(yōu)化效率

3、、減小誤差，但因存在彎扭耦合故很難獲得屈曲控制方程的精確解。一些研究通過忽略彎扭耦合以簡化屈曲響應(yīng)分析來優(yōu)化鋪層順序，并結(jié)合限定彎扭耦合相對大小、鋪層鋪設(shè)形式等附加約束來減小計(jì)算偏差，這必然會對最優(yōu)解的獲得造成阻礙。屈曲問題的近似求解方法如里茲法和有限元法可以考慮彎扭耦合，提供滿足精度要求的計(jì)算結(jié)果。有限元法能求解復(fù)雜結(jié)構(gòu)的屈曲響應(yīng)，但與隨機(jī)算法結(jié)合的計(jì)算成本不容忽視，研究者常采用近似模型替代原問題模型來提高計(jì)算效率，而里茲法能方便有效地求解規(guī)則形狀和常規(guī)邊界條件的層合板屈曲響應(yīng)問題。針對層合板屈曲優(yōu)化的鋪層順序設(shè)計(jì)問題，本文以規(guī)則簡支板為研究對象，以離散鋪層角度為設(shè)計(jì)變量，通過改進(jìn)DSA算法

4、的新解產(chǎn)生方式提高算法全局收斂性和穩(wěn)定性，采用里茲法進(jìn)行層合板屈曲響應(yīng)分析來考慮彎扭耦合對優(yōu)化結(jié)果的影響，從而實(shí)現(xiàn)了層合板屈曲載荷系數(shù)最大化的鋪層順序優(yōu)化。此外，本文還研究了不同角度增量下具有更大離散設(shè)計(jì)空間時(shí)改進(jìn)算法的性能表現(xiàn)和優(yōu)化解的變化規(guī)律，并比較層合板不同長寬比、載荷比以及不同鋪層數(shù)對優(yōu)化結(jié)果的影響。1優(yōu)化問題描述對稱層合板由N =2n層鋪層組成，長寬為ab，每層鋪層具有相同厚度t0，第k層的纖維鋪層角度為k(k =1，2，n)，則層合板鋪層順序可表示為1/2/k/ns。圖1中，Nx、Ny分別為x 和y 方向作用在板中面的壓力載荷。軸壓載荷作用下對稱層合板的屈曲控制方程為D114wx4

5、 +2(D12+2D66) 4wx2y2+D224wy4 +4D16 4wx3y+4D26 4wxy3 =Nx2wx2 +Ny2wy2 (1)Dij = 23 N/2k=1Q(k)ij (z3k-z3k-1)i，j=1，2，6式中，w 為層合板中面位移;Dij為層合板彎曲剛度;Q(k)ij為第k 層的材料偏軸剛度。Q(k)ij由材料正軸剛度Qij按下式計(jì)算得到：Q11Q12Q16Q22Q26Q熿燀燄66燅=c42c2s2s44c2s2c2s2c4+s4c2s2 -4c2s2c3scs3-c3scs32(cs3-c3s)s42c2s2c44c2s2cs3c3s-cs3 -c3sc3s-cs3c2

6、s2 -2c2s2 -c2s2c2-s熿燀燄2 燅Q11Q12Q22Q熿燀燄66燅(2)c=coss=sinQ11 =E1/(1-1221)Q12 =12E2/(1-1221)Q66 =G12式中，E1、E2分別為單層材料在1、2主軸方向的彈性模量;G12為1-2平面內(nèi)的剪切模量;12、21分別為1-2平面內(nèi)的縱向和橫向泊松比。從式(2)可以看到，由于鋪層偏軸剛度與鋪層角具有復(fù)雜的函數(shù)關(guān)系，造成屈曲優(yōu)化問題具有多極值特征，故使得以鋪層角為設(shè)計(jì)變量的鋪層順序設(shè)計(jì)較為困難。四邊簡支板的邊界條件為x =0，a;w =0;Mx =0y=0，b;w =0;My = 0(3)式中，Mx、My分別為x 軸和

7、y 軸方向的力矩。上述邊界條件中，對于彎扭耦合項(xiàng)D16、D26為零或其值相對較小的幾種層合板，可采用封閉形式求解屈曲載荷系數(shù)b：b(m，l)=2D11 (ma) 4+2(D12+2D66)(malb) 2+D22 (lb) 4(ma) 2Nx +(lb) 2Ny(4)其中，m、l為x、y軸向的半波數(shù)(m，l=1，2，)，不同m、l組合對應(yīng)不同b(m，l)，相應(yīng)屈曲載荷為bNx和bNy。而D16和D26不為零時(shí)，因法向位移不能通過變量分離技術(shù)獲得封閉形式解析解，故常采用里茲法近似求解。根據(jù)里茲法，取層合板中面z向位移函數(shù)：w(x，y)=Mm=1Ll=1cmlsinmxa sinlyb(5)其中，

8、cml為待定系數(shù)，M =L 為事先給定的正整數(shù)，一般在1，10內(nèi)取值。層合板的彎曲應(yīng)變能V 為V = 1 2b0a0dxdy (6)=D11(2wx2)2+D22(2wy2)2+4D66(2wxy)2+2D122wx22wy2 +4(D162wx2 +D262wy2)2wxy中面載荷做功為W = 1 2b0a0Nx(wx)2+Ny(wy)2dxdy (7)基于最小勢能原理有：(V +W)cml=0m =1，2，M;l=1，2，L(8)由此建立關(guān)于cml的M L個(gè)線性方程組，問題表現(xiàn)為方程組系數(shù)矩陣的特征值問題，通過子空間迭代法等標(biāo)準(zhǔn)特征值求解方法獲得一系列特征值，即屈曲載荷系數(shù)b(m，l)。而

9、臨界屈曲載荷系數(shù)cb對應(yīng)層合板能抵抗屈曲的最小臨界屈曲載荷，是一系列b(m，l)值中的最小值，即cb = minb(m，l) (9)最小臨界屈曲載荷則通過公式Nx，cr=cbNx和Ny，cr=cbNy求得。層合板屈曲優(yōu)化問題就是要最大化臨界屈曲載荷系數(shù)cb，以層合板鋪層角= (1，2，k，n)T(k=1，2，n)為設(shè)計(jì)變量，臨界屈曲載荷系數(shù)最大化問題的數(shù)學(xué)表達(dá)為min-cb()Lkk Uk= (1，2，k，n)Tk=1，2，烍烌n 烎(10)Uk和Lk分別為鋪層角k的上下限，構(gòu)成設(shè)計(jì)變量的搜索空間S; 為離散設(shè)計(jì)變量向量，其取值為Lk+、Lk+2、Uk -、Uk ，為角度增量。如鋪層角k -9

10、0，90，=45時(shí)，k取值為-45、0、45、90。因此越小，鋪層角的設(shè)計(jì)空間越大。2模擬退火算法的改進(jìn)SA算法是一種模擬高溫金屬冷卻過程的全局隨機(jī)搜索算法，基于單點(diǎn)串行搜索并以一定概率接受差點(diǎn)，因此會丟失優(yōu)化過程中獲得的優(yōu)化點(diǎn)。DSA算法區(qū)別于SA算法之處在于DSA算法擁有一個(gè)點(diǎn)集合，按兩種機(jī)制產(chǎn)生的新點(diǎn)被接受后僅替換該集合中的最差點(diǎn)，實(shí)現(xiàn)了優(yōu)化過程最佳點(diǎn)的記憶保存。DSA算法的新點(diǎn)產(chǎn)生混合機(jī)制覆蓋了全局粗搜索與局部細(xì)搜索，但降低了算法搜索效率和優(yōu)化穩(wěn)定性，易導(dǎo)致算法陷于局部極值。本文針對層合板屈曲優(yōu)化問題，通過改進(jìn)DSA算法的初始溫度確定方法和新點(diǎn)修正方法，增加動態(tài)調(diào)整的新點(diǎn)產(chǎn)生方式，提高

11、算法求解層合板鋪層順序優(yōu)化的穩(wěn)定性和計(jì)算效率。初始集合和初始溫度確定首先在搜索空間S 內(nèi)按均勻分布隨機(jī)產(chǎn)生20n(n為鋪層角度個(gè)數(shù))個(gè)點(diǎn)的集合，若集合中包含可行解則初始溫度T0 =500n，否則T0 =20XXn，使得T0與變量規(guī)模和問題難易程度相關(guān)，從而確定一個(gè)相對高且計(jì)算經(jīng)濟(jì)的初始溫度值。因式(10)為無約束問題，因此T0 =500n。繼續(xù)從該集合中選出目標(biāo)函數(shù)值較小的7(n+1)個(gè)點(diǎn)構(gòu)成初始點(diǎn)集合A，并標(biāo)記A 中對應(yīng)的最大目標(biāo)函數(shù)值fH和最小目標(biāo)函數(shù)值fL。這樣避免了DSA算法以概率1抽樣確定高初始溫度所耗費(fèi)的目標(biāo)評價(jià)時(shí)間和因初始溫度高而增加的迭代收斂過程，從而有助于提高算法計(jì)算效率。新

12、點(diǎn)產(chǎn)生方式在DSA算法中新點(diǎn)按混合機(jī)制產(chǎn)生：=2-n+1PU(S) U烅烄烆i(13)式中，i為修正后的第i個(gè)分量;mod()為取余數(shù)函數(shù)。測試數(shù)據(jù)表明式(13)可以提高計(jì)算效率且不影響優(yōu)化結(jié)果。新點(diǎn)接受概率新點(diǎn)的接受概率Pa計(jì)算式為Pa =1f()fHe(f()-fH)/Tf()f H(14)式中，f()為新點(diǎn)的目標(biāo)函數(shù)值;T 為當(dāng)前溫度。如果新點(diǎn) 被接受，則替代A 中最差點(diǎn)(fH所對應(yīng)點(diǎn))，并更新fH和fL。降溫策略DSA算法采用指數(shù)形式降溫策略：Tk+1 =k+1Tkk0 (15)降溫系數(shù)k+1min，max(k1)按下式計(jì)算：k+1 =maxLk Lkmin+(k -min)Lk/Lk

13、+1 Lk =Lk，Lk+1 Lkmax-(max-k)Lk+1/Lk Lk =Lk，Lk+1 L烅烄烆k(16)Lk =10n+10n(1-e(fL-fH)式中，Lk為溫度Tk時(shí)的計(jì)算馬氏鏈長度，隨A中fL與fH的差值的減小而減小;Lk為Tk時(shí)實(shí)際馬氏鏈搜索長度。顯然，在Lk內(nèi)找到更好的新點(diǎn)時(shí)，應(yīng)緩慢降溫以強(qiáng)化該鄰域的搜索，因此有k+1 =max;而在Lk內(nèi)沒搜尋到好于當(dāng)前最佳點(diǎn)的新點(diǎn)時(shí)，有Lk=Lk，再結(jié)合Lk+1計(jì)算k+1，從而形成自適應(yīng)降溫策略。初始降溫系數(shù)1 =(min+max)，本文取max =和min =。收斂準(zhǔn)則DSA算法采取雙迭代收斂條件：Tk 1fH -fL 2(17)其中

14、，1和2為較小實(shí)數(shù)，本文取1 =2 =，即迭代溫度足夠低和A中的點(diǎn)足夠集中時(shí)算法終止。3算例分析基于MATLAB平臺編程實(shí)現(xiàn)了DSA 算法和改進(jìn)DSA算法，并分別采用這兩種算法進(jìn)行臨界屈曲載荷系數(shù)最大化的鋪層順序優(yōu)化設(shè)計(jì)，比較算法改進(jìn)效果。鋪層材料為Graphite/Epoxy6，E1=，E2 = GPa，G12=，12 =，厚度t0 =。每個(gè)問題獨(dú)立運(yùn)行50次。算例1考慮四邊簡支、兩軸向承壓的64層對稱層合板，長a=，寬b=a/2，厚度t= mm，層合板承受面內(nèi)壓力Nx =/mm，Ny=l1Nx，l1為兩軸向載荷大小比例。假設(shè)該層合板為對稱均衡層合板，由鋪層角為02、45或902(即=45)

15、的雙層鋪層組鋪設(shè)而成。l1=1時(shí)DSA算法和改進(jìn)DSA算法基于式(4)和里茲法獲得的優(yōu)化結(jié)果，Nfa為平均目標(biāo)函數(shù)評估次數(shù)，Nit為平均迭代次數(shù)，成功率為所獲優(yōu)化解與理論最優(yōu)值的絕對誤差小于設(shè)定精度(本文設(shè)為)的次數(shù)與50次統(tǒng)計(jì)次數(shù)之比。DSA 算法、改進(jìn)DSA 算法以及SA算法的優(yōu)化結(jié)果為獨(dú)立運(yùn)行50次統(tǒng)計(jì)所得，而GA算法的優(yōu)化數(shù)據(jù)是獨(dú)立運(yùn)行200次統(tǒng)計(jì)結(jié)果。結(jié)果顯示改進(jìn)DSA 算法獲得了與SA 算法、GA算法相同的7個(gè)全局最優(yōu)解，而DSA算法僅獲得了其中3個(gè)，且GA 算法獲得的優(yōu)化解和實(shí)為局部極值，說明本文改進(jìn)措施提高了DSA算法的全局搜索能力，能有效求解層合板屈曲優(yōu)化問題。數(shù)據(jù)顯示DSA

16、算法與改進(jìn)DSA算法雖然均獲得3個(gè)優(yōu)化鋪層，但改進(jìn)DSA 僅用了DSA所需Nfa和Nit的約1%的計(jì)算成本，獲得了高于DSA 的收斂成功率和優(yōu)于DSA 算法的cb優(yōu)化結(jié)果，說明本文改進(jìn)措施有效避免了DSA算法易陷入局部極值和迭代后期收斂緩慢的不足，提高了DSA算法的計(jì)算效率和穩(wěn)定性。另外，基于式(4)所獲優(yōu)化鋪層按里茲法計(jì)算獲得的cb基本各不相同，且均小于基于里茲法獲得的cb，說明層合板屈曲優(yōu)化結(jié)果受響應(yīng)分析方法計(jì)算精度的影響明顯，考慮彎扭耦合的里茲法提高了計(jì)算精度，能獲得更好的優(yōu)化結(jié)果。假設(shè)上述64層對稱層合板由單層鋪層鋪設(shè)而成，且k 0，90， 分別為30、15、10、5，l1分別為、時(shí)改

17、進(jìn)DSA算法求解式(10)的優(yōu)化結(jié)果。算例2四邊簡支兩軸向承壓的對稱層合板，長a = mm，長寬比為a/b，Nx =/mm，Ny =l1Nx。鋪層角度取值范圍為-90，90，總鋪層數(shù)為N?；诶锲澐ê透倪M(jìn)DSA算法求解式(10)，獲得單軸(l1 =0)和雙軸(l1 =1)向壓力作用下，不同N、a/b及時(shí)的優(yōu)化結(jié)果，表中Nopt、Iopt分別為50次獨(dú)立運(yùn)算中獲得全局優(yōu)化解的個(gè)數(shù)和次數(shù)?；谑?4)的優(yōu)化結(jié)果與最優(yōu)鋪層完全一致，但其對應(yīng)的cb與按里茲法計(jì)算的cb差異較大，與l1=時(shí)的優(yōu)化結(jié)果相比，可發(fā)現(xiàn)基于里茲法獲得的cb均高于基于式(4)所獲最優(yōu)鋪層的里茲法解，說明對于彎扭耦合不可忽視的層合板

18、，式(4)因忽略彎扭耦合產(chǎn)生較大計(jì)算偏差而不能獲得優(yōu)化解，考慮彎扭耦合的里茲法能有效獲得問題優(yōu)化解。l1一定時(shí)cb隨減小而增大，說明增大設(shè)計(jì)空間能獲得更好的優(yōu)化結(jié)果。隨著l1增大，cb逐漸減小且受的影響減弱，層合板趨于用較大的鋪層角來增強(qiáng)窄邊承載能力以提高屈曲穩(wěn)定性。l1=0的優(yōu)化結(jié)果，可以看到，給定N和a/b時(shí)，分別取45、15、5均獲得相同的優(yōu)化鋪層，獲得的cb隨N 減小而減小，隨a/b增大而增大，當(dāng)a/b=時(shí)0鋪層角最優(yōu)，而a/b1時(shí)優(yōu)化鋪層角為45和-45的不同組合。l1 =1時(shí)N 分別取32、64和128時(shí)的優(yōu)化結(jié)果，數(shù)據(jù)顯示，N 一定時(shí)，ab隨a/b增大而增大，鋪層角隨a/b增大逐漸從接近x軸向分布趨于接近y軸向分布，當(dāng)a/b=1時(shí)的優(yōu)化鋪層角為45和-45的不同組合。另外，從數(shù)據(jù)發(fā)現(xiàn)，給定N、a/b和下的優(yōu)化問題50次獨(dú)立運(yùn)算所獲得的cb差異均較小，相對波動幅值小于1%，說明改進(jìn)DSA算法具有較好的計(jì)算穩(wěn)定性，能穩(wěn)定地進(jìn)行層合板屈曲問題的鋪層順序優(yōu)化;而給定N 和a/b時(shí)，隨著減小，Iopt逐步減小，cb及其波動幅值逐漸增大，說明小值增大了