改進算法在導彈結構熱防護一體化中優化研究

齊名軍 , 王志寶, 吳 凱, 鄧 紅, 谷海紅

(1. 鶴壁職業技術學院, 河南 鶴壁 458030; 2. 東北石油大學, 黑龍江 大慶 163318; 3. 天津農學院, 天津 300191; 4. 黑龍江工程學院, 哈爾濱 150050; 5. 河南理工大學, 河南 鶴壁 458030)

0 引 言

導彈進入地球大氣層超高速飛行時, 其外部由摩擦產生的熱量在300 ℃到1 500 ℃之間, 給結構抗變形能力及熱防護體系帶來極大的考驗[1], 而結構熱防護一體化技術的出現, 為導彈在大氣層超高速飛行所面臨的困境提供了新思路[2]。 結構熱防護一體化技術的實質是將承載負荷和隔熱雙重功能集于一身的技術[3], 該技術還有較高的比強度、 比剛度及便于修理卸載等優點, 在各類高超聲速飛行器中得到廣泛應用[4]。 文獻[5]把結構熱防護一體化技術應用到航天器中, 并構建航天器一體化結構綜合性能評價方法與體系。 文獻[6]把該技術應用到超高速飛行導彈中, 實驗表明該技術在承載和隔熱兩個方面具有良好的效果。 仿生學的出現加快新型承載/防熱一體化結構的發展, 如蓮藕結構等, 雖然雙重功能良好但制造工藝較為困難[7]。 目前該技術在超高聲速飛行器應用研究中的最大難點是在壁寬毫米量級超高速飛行器上難以實現[8], 原因在于隨著導彈隔熱層壁寬減小, 腹板熱短路效應增強, 熱防護功能失靈, 因此腹板鏤空設計被淘汰。 文獻[9]研究證實了波紋夾芯結構制造工藝相對簡單, 未來應用優勢明顯。

隨著群體智能優化算法的出現, 國內外學者為提高一體化結構性能, 把智能優化算法嘗試性應用到該體系中。 文獻[10]采用SA算法且僅從結構方面降低波紋夾芯一體化結構質量; 文獻[11]利用GA算法僅從隔熱功能改善波紋夾芯板熱防護效率; 文獻[7]利用GA算法改善兼顧承載/熱防護結構一體化的性能, 效果一般。 對算法深入研究發現, 群體智能優化算法在諸多領域應用中取得很大成功, 尤其在低維函數應用表現較佳, 但在解決復雜高維函數優化方面存在如下缺點[12]: 算法在種群個體進化期間多樣性差, 收斂速度慢且易陷入“早熟收斂”, 無法獲得高精度解。 樽海鞘算法(Salp Swarm Algorithm, SSA)也屬于群體智能算法, 同樣具有上述不足。 本文針對該算法存在的不足, 做出以下改進: (1)引入Sobol 序列對初始種群進行均勻化處理, 保持種群多樣性, 提升算法全局勘探能力; (2)使用動態權重調節機制策略, 提高算法收斂速度; (3)融合變異因子和反向解策略, 有利于算法跳出局部最優, 提高算法精度; (4)多用一個臨時儲備庫可以從單目標理論上最佳解集合獲得滿足實際需要的Pareto最優解。 最后, 把改進算法應用到導彈波紋夾芯結構熱防護一體化中的熱防護效率、 結構質量等多目標優化函數中, 利用罰函數法把多目標函數轉化為單目標最小化進行求解計算。

1 導彈結構熱防護一體化數學模型

導彈在太空大氣層超高聲速飛行時, 氣動熱與氣動力載荷都要考慮到, 以文獻[7]為例, 如果導彈熱防護目標為300 ℃, 其艙體模型大小如圖1所示, 氣動熱載荷轉化為溫度載荷如圖2所示, 且作用在結構外壁, 氣動力載荷等效處理為1 000 N·m。 若實驗模型無輻射, 結構內壁底部絕緣, 仿真將會推出一種理想的結論。

圖1 導彈艙體模型Fig.1 Missile cabin model

圖2 溫度載荷Fig.2 Temperature load

通過文獻[7]中對導彈的熱-結構耦合實驗, 得到導彈艙體的溫度分布, 如圖3所示。 分析可知, 其溫度分布呈“U型”軸對稱狀態, 隨著位置的增加, 兩端高溫區逐漸對稱地過渡到中部低溫區, 位置在區間38～160 mm時溫度保持不變, 距離艙體兩端16 mm范圍內, 內表面溫度大于300 ℃。 為防止彈體損壞, 在兩端位置安裝連接環, 如圖4所示。 此時, 范圍總寬度增加, 熱容增大, 使艙體內表溫度下降。

圖3 溫度分布Fig.3 Temperature distribution

圖4 艙體結構設計Fig.4 Cabin structure design

熱-結構耦合分析的應力分布結果如圖5所示。 由圖5可知, 該結構真實分布特點反映出其具有較好的承載能力。 艙體中部應力強度大, 艙體結構極易發生破壞, 屬于應力薄弱范圍區域。 由文獻[7] 可知, 通過熱-結構仿真計算時, 對網格進行了收斂性檢驗, 所得到的距離艙體頂端16 mm的溫度與結構最大應力如表1所示。 其中, 綜合計算精度與效率, 網格尺度取0.1 mm。

表1 溫度與結構最大應力Table 1 Temperature and maximum structural stress

圖5 應力云圖Fig.5 Stress cloud

由表1可知, 距離艙體端部16 mm時, 溫度越低越好, 并且導彈結構所對應的最大von-Mises應力越小越好, 這樣才能達到熱防護目標。

文獻[7]研究了外壁、 內壁及隔熱層寬度對熱防護性能的影響, 得出如下結論: 外壁寬度對熱防護性能作用較小, 擴大隔熱層與內壁寬度確實能夠增強結構熱防護性能。 但局限于艙體壁寬大小約束, 隔熱層與內壁寬度不能無限增大, 且外壁寬度會隨之降低, 結構承載能力變弱而發生變形。 因此, 在導彈結構熱防護一體化優化研究中, 既要考慮結構熱防護性能, 又要考慮彈體結構輕質化。

由以上模型分析可知, 進行彈體優化時, 需要研討的要點如下: (1)保證導彈壁寬為3 mm; (2)結構承載能力要強, 在兼顧力-熱雙重作用下, 優化后彈體結構應力不能超出材料極限承載力; (3)考慮目前工廠加工精確程度, 將導彈壁寬精確到千分位。 實驗隔熱層填充材料采用 SY1000, 其余結構部分均采用鈦合金材料, 具體材料參數見文獻[7]。

本文案例是多目標優化問題, 即以結構熱防護效率最大和結構質量最小為最終求解目標, 具體定義如下:

miny={y1,y2}

(1)

(2)

式中:x1為結構隔熱層寬度;x2為結構內壁寬度;η(x1,x2)為結構熱防護效率;m(x1,x2)為結構質量。

導彈一體化熱防護結構的優化數學模型約束條件為

(3)

式中:σ(x1,x2)為結構最大應力。

本優化案例是多目標函數問題, 但考慮到以提高結構熱防護效率優先, 采用單目標優化, 因此將式(2)中的結構質量和式(3)中的不等式約束條件都采用懲罰函數法, 轉化為如下函數:

μ(m(x1,x2))=max{m(x1,x2)-m0, 0}

φ(σ(x1,x2))=max{σ(x1,x2)-σ0, 0}

ψ(x1,x2)=max{(x1+x2)-3, 0}

δ(x1)=max{x1-3, 0}

θ(x1)=max{x1-3, 0}

當以上函數均滿足約束條件時, 其值為0; 反之其值為正。 因此, 本文目標優化函數轉化為

F(x1,x2)=min{y1}+PF(μ(m(x1,x2)+

φ(σ(x1,x2))+ψ(x1,x2)+δ(x1)+

θ(x2))

(4)

式中:PF為懲罰因子。

根據文獻[7, 13]中熱-結構耦合模型有限元分析得到結構應力函數為

(5)

參照文獻[13], 對于結構應力von Mises仿真計算方法為: 先對結構進行瞬態溫度歷程計算, 再選取某時刻的溫度值作為載荷, 然后利用式(5)求解。 本文同樣采用文獻[7]中的熱防護效率仿真公式:

η=100×(tout-tin)/tout

(6)

式中:tout為測點外壁溫度;tin為測點內壁溫度。

本文數學建模采用邊界元法, 其中用到有限元分析軟件中的應力單元和溫度單元, 模型準確的邊界條件為: 保證導彈壁寬為3 mm, 距離艙體端部16 mm; 距離艙體端部16 mm處的溫度載荷最大不超過300 ℃, 導彈結構在這種溫度下的最大von-Mises應力不超過396 MPa。

2 基本樽海鞘算法

設M×K維為種群覓食空間,M為覓食點數,K為種群數量。 覓食范圍內食物為A={A1,A2, …,AK}T, 其相應位置為xn={xn1,xn2, …,xnk}T,n=1, 2, …,M, 其解上限為bh, 下限為bl,xM×K=rand(M,K))(bh-bl)+bl, 種群領導者位置更新為

(7)

(8)

3 個體擾動多策略樽海鞘算法

針對樽海鞘算法在優化方面的不足, 提出了個體擾動多策略樽海鞘算法(Individual Disturbance Multi-Strate-gy Salp Swarm Algorithm, IDMSSA)。 新算法包括Sobol 序列, 動態權重調節機制, 早熟收斂判斷機制、 變異因子和反向解策略, 以及針對本文多目標優化問題所用到的Pareto最佳解儲備庫, 其核心思想是提高算法性能。

3.1 初始種群均勻化處理

在算法進化過程中, 群體智能算法優化性能與初始種群分布有關, 初始種群均勻分布能夠保證算法多樣性, 增強算法的全局搜索能力, 提高算法效率[14]。 受文獻[15]的啟發, 結合初始種群分布均勻化的思想, 考慮到樽海鞘算法初始化是采取隨機序列, 例如在[0, 1]之間生成維度為9、 種群數為200的隨機序列初始化(如圖6(a)所示), 其個體均勻分布比較差。 為加強算法全局勘探能力, 提高算法的多樣性, 充分發揮初始解的作用, 引進Sobol 序列來初始化種群。 該序列在求解區間呈均勻化分布, 周期小, 采樣快, 在處理高維復雜性非線性函數效率更高[16], 其在同樣范圍內種群分布如圖6(b)所示。 對比可知, Sobol序列明顯優于隨機序列, 其種群分布更均勻。

圖6 種群初始化Fig.6 Population initialization

Sobol 序列公式如下:

yn=ymin+εn·(ymax-ymin)

(9)

式中:εn∈rand(0, 1);ymin,ymax為當前代個體適應值的最小值和最大值。

3.2 動態權重調節機制

通過研究基本樽海鞘算法機理及位置更新式(8)可知, 追隨者只受前一個樽海鞘個體的影響, 新一代個體位置更新是由前一代相鄰個體位置所決定, 相鄰不同優劣個體所起作用相同, 種群沒有把優勢個體的優勢展現出來, 優勢個體對整個種群中其他個體的影響逐代下減, 使得種群群體協作能力發揮效用不佳, 降低算法的收斂速度。 受文獻[17]的啟發, 結合動態權重能夠引領整個群體朝優化方向進化的思想, 本文引入動態權重, 充分發揮其在進化過程中優勢個體的引導作用, 減少劣勢個體帶來的不良影響, 擺脫追隨者對前一個樽海鞘位置的完全依賴。 因此, 種群中樽海鞘個體能夠根據優勢權值動態變化, 不僅有利于提高搜索的靈活性, 還能增強整個算法勘探和開發能力, 從而進一步提升算法速度。 假設種群規模有N個個體, 其適應度分別為f1,f2, …,fN, 相對應的權重為w1,w2, …,wN, 則新一代位置更新按照式(10)進行:

(10)

動態權重wi為

(11)

3.3 早熟收斂判斷機制、 變異因子和反向解策略

由式(10)可知, 動態權重調節機制確實能夠使算法朝著最優解的方向搜索, 其變化能強化進化過程中優勢個體的影響, 同時也弱化了優勢個體進入早熟時, 算法難以跳出局部最優的能力, 從而降低算法收斂速度。 為此, 引入早熟收斂判斷及早熟處理機制。 該機制能夠有效地根據所處狀態實時脫離局部最優, 提高算法效率。 局部最優早熟判斷與處理如圖7所示。

圖7 早熟判斷預處理機制Fig.7 Preprocessing mechanism of precocious judgment

早熟收斂判斷是“早熟處理”的前提[18]。 研究發現, 樽海鞘算法無論是“局部收斂”還是“全局收斂”, 樽海鞘種群個體都會出現“聚堆”現象。 本文將樽海鞘個體最佳適應值的變化狀態作為“早熟收斂”判斷的條件, 設xavg為樽海鞘當代的平均位置,α2為樽海鞘群體位置方差[19], 定義為

(12)

式中:x為歸一化因子, 主要約束α2的大小, 即

(13)

式(12)能說明當代種群個體的“堆集”程度:α越小, 說明樽海鞘種群的“堆集”程度越大。 如果算法不滿足預設結束條件, 將使樽海鞘缺乏多樣性而呈現“早熟收斂狀態”。 當α2≤c時(c為預設常數), 則認為算法陷入“停滯”, 即算法進入早熟狀態。 受文獻[19]的啟發, 高斯變異、 柯西變異兩種函數同范圍同緯度內的不同分布如圖8所示, 考慮到柯西變異比高斯變異更具有優勢(其分布區間較大而且緊湊), 充分利用其分布函數“橫軸具有較大或較小的值, 縱軸也有一定概率得到相應的值”的思想, 即對種群個體的最優位置施加一個變異擾動, 使其發揮強大的局部調節能力, 進而幫助算法跳出局部最優, 克服早熟的不足, 提高了算法的求解精度。

圖8 高斯、 柯西變異分布Fig.8 Gaussian and Cauchy variation distribution

對局部早熟最優解進行柯西變異:

x(t+1)=Cauchy(0, 1)⊕xbest(t)+xbest(t)

(14)

式中:xbest為當代最優解; Cauchy(0, 1)為標準柯西分布。

圖9 C點情況圖Fig.9 C point situation diagram

設Y=(y1,y2, …,yd)為d維空間的一點,yj∈[aj,bj],j∈1, 2, …,d, 則yj的反向解為

(15)

3.4 多目標Pareto最佳解替換機制

與單目標函數優化不同, 多目標函數優化的最優解不一定是滿足實際需要的最佳解。 因此, 受文獻[22]的啟發, 結合本文多目標函數優化目標的特點, 引入兩個存儲種群最佳解集合: 當代種群進化集合X(t)和Pareto最優解臨時儲備庫M, 其集合數目都為n。 二者的不同之處在于X(t)僅是理論上的最佳解, 而Pareto最優解臨時儲備庫M不僅是理論上的最佳解還是符合實際需要的最佳解。 隨代逐步進化替換, 其中Pareto最優解臨時儲備庫M進行更換的公式[23]如下:

(16)

式中:x,y為兩個不同種群個體;ε為調節Pareto最優解間隔常數。 編程實現時, 在常數基礎上逐漸減小ε的大小, 直到種群數量達到規定的群體規模數, 這樣得到當前種群中Pareto最優解, 即Pareto最優陣面。

3.5 算法流程圖

個體擾動多策略樽海鞘算法從多方面增強樽海鞘算法的尋優性能, 其算法流程如圖10所示。

圖10 算法流程圖Fig.10 Algorithm flow chart

4 仿真實驗及其分析

為了充分驗證每個策略對算法性能的影響, 與文獻[7]中的遺傳算法(GA)進行仿真實驗對比。 設基本樽海鞘算法為SSA, 單增Sobol 序列策略樽海鞘算法為SSA1, 單增動態權重調節機制樽海鞘算法為SSA2, 單增早熟收斂判斷機制、 變異因子和反向解策略樽海鞘算法為SSA3。 所有對比算法均在相同仿真平臺上進行演示, 進化種群數量均設置為30, 最大進化代數Tmax= 100,ε=0.001,c=0.01,n=10, 其他參數與文獻[7]相同。 遺傳算法中, 交叉因子為0.7, 變異算子為0.1。 實驗仿真結果輸出后, 在Pareto最優解臨時儲備庫M中輸出實際需要的最佳解, 如表2所示。

表2 不同算法最優解對比Table 2 Comparison of optimal solutions of different algorithms

由表2可知, 在熱防護效率方面, SSA3>SSA2>SSA1>SSA。 這說明, 早熟收斂判斷機制、 變異因子和反向解策略對算法性能作用較大, 動態權重調節機制對算法性能作用次之, Sobol 序列策略初始均勻化對算法性能的作用較小。 主要原因在于SSA3能夠通過早熟判斷機制進行位置變異而跳出“局部最優”, 提高了算法精度。 3種改進算法所得結果都優于基本樽海鞘算法的優化結果, 說明不同的改進算法都是有效的, 只是不同策略對改進算法性能的貢獻大小不同。

整體來看, 6種不同優化算法都優于初始模型所得熱防護效率和結構質量這兩項指標的優化結果, IDMSSA模型最優解中熱防護效率為73%, 明顯優于其他對比算法, 而且其所求結構質量也是所有算法中最小的。

從導彈結構熱防護的優化目標函數角度可知: (1)要滿足結構熱防護效率越大越好; (2)所求導彈結構質量越小越好; (3)要滿足距離艙體端部16 mm處的溫度不超過300 ℃, 并且導彈結構在該溫度下的最大von-Mises應力越小越好, 最大不超過396 MPa(與300 ℃相對應), 其最理想的結構最大應力為392 MPa。

根據式(5)和文獻[6]中熱-結構耦合分析方法流程圖, 可計算出本文改進算法求得的溫度為294 ℃, 導彈結構最大von-Mises應力為379 MPa; SSA所得溫度為297 ℃, 導彈結構最大von-Mises應力為388 MPa。 與文獻[7]中的遺傳算法(GA)相比(所得溫度為296.2 ℃, 導彈結構最大von-Mises應力為384 MPa), 本文改進算法所得的熱防護效率、 導彈結構質量、 溫度及對應的最大von-Mises應力都是最優的, 均優于初始模型算法和遺傳算法。

本文改進算法所得的質量、 溫度、 應力水平和熱防護效率上均比遺傳算法、 初始模型算法具有優勢。 原因在于: 在過去, 基于嚴格機理模型所得到的優化命題通常具有方程多、 非線性強、 變量維度高、 因數多、 影響廣、 難度高和規模大等特點, 這使相關量的存儲及命題的求解都相當困難, 這些問題必須由一個相當有效的優化工具來進行求解。 面對大型問題, 常規的優化方法已經無能為力, 計算速度、 收斂性、 初值敏感性等都不能滿足需要。 群體智能優化算法的興起, 豐富了現在優化技術的內涵, 促進優化技術的快速發展, 成為求解現代優化問題最優值的有力工具, 為解決具有非線性、 多極值等以前難以處理的傳統高維復雜優化函數提供了新的思路。 而與遺傳算法相比, IDMSSA雖然也屬于群體智能優化算法, 但其具有的最大優點是: 該算法引入Sobol序列策略對初始種群均勻化處理, 提高解的質量, 保持種群多樣性; 使用動態權重調節機制, 更好地發揮優勢個體的引導作用, 提高算法速度; 利用早熟收斂判斷機制、 變異因子和反向解策略, 增強種群探索效率, 有利于算法跳出局部最優, 提高求解精度; 而遺傳算法具有所求精度不高, 易陷入“局部最優”的缺點。 因此, 本文算法能夠克服其不足而同時實現其質量、 溫度、 應力水平和熱防護效率的優勢。

5 結 束 語

針對樽海鞘算法在解決高維復雜多目標函數優化方面存在多樣性差、 早熟收斂及無法求得精確解等缺點, 引入Sobol序列, 動態權重調節機制, 早熟收斂判斷機制、 變異因子和反向解策略來提高算法性能, 并對導彈基于波紋夾芯結構熱防護一體化熱防護效率、 結構質量進行仿真實驗。 通過對比3種算法的優化結果, 驗證了改進算法求得熱防護效率、 結構質量兩種參數的精度高于其他算法, 驗證了改進算法的有效性和優越性。