武澤平，王東輝，胡　凡，張為華

(國防科技大學(xué) 航天科學(xué)與工程學(xué)院，長沙　410073)

?

基于代理模型的固體火箭發(fā)動機裝藥設(shè)計①

武澤平，王東輝，胡凡，張為華

(國防科技大學(xué) 航天科學(xué)與工程學(xué)院，長沙410073)

提出一種基于代理模型的固體火箭發(fā)動機裝藥優(yōu)化設(shè)計方法。將代理模型技術(shù)用于基于Level Set的固體火箭發(fā)動機裝藥設(shè)計，并提出了基于Pareto前沿采樣的并行加點準(zhǔn)則，進一步提升算法性能。以后翼柱型裝藥為例，以推力平穩(wěn)為設(shè)計目標(biāo)，得到了優(yōu)良的推力性能，驗證了方法的有效性。

固體火箭發(fā)動機；代理模型；并行加點準(zhǔn)則；優(yōu)化設(shè)計方法；LevelSet；裝藥設(shè)計

0　引言

固體火箭發(fā)動機設(shè)計是火箭/導(dǎo)彈設(shè)計的重要內(nèi)容，發(fā)動機設(shè)計質(zhì)量好壞直接影響火箭/導(dǎo)彈性能。裝藥設(shè)計是固體火箭發(fā)動機設(shè)計的重要組成部分，對固體火箭發(fā)動機性能的優(yōu)劣起著決定性作用[1]。裝藥設(shè)計與發(fā)動機總體設(shè)計聯(lián)系密切，往往在發(fā)動機總體方案論證時，就已確定了推進劑的配方，甚至選定了藥型。在通常情況下，裝藥設(shè)計最主要的任務(wù)是確定藥型尺寸，使其滿足一定的內(nèi)彈道性能。在過去的幾十年里，形成了一系列的裝藥設(shè)計方法[2-4]，此類方法通常采用智能優(yōu)化方法驅(qū)動計算量較小的低精度燃面退移仿真方法進行裝藥設(shè)計，需要大量的迭代計算，且計算精度難以保證。

隨著高性能計算技術(shù)不斷發(fā)展和廣泛應(yīng)用，發(fā)動機仿真計算精度明顯提高，利用計算機技術(shù)進行發(fā)動機優(yōu)化設(shè)計，以尋求性能優(yōu)良、結(jié)構(gòu)可靠的發(fā)動機設(shè)計方案。然而，精確仿真模型往往具有計算復(fù)雜度高、計算耗時長的特點。此時，智能優(yōu)化方法的大量迭代計算通常難以接受[5]。代理模型技術(shù)可通過設(shè)計空間少量樣本點，對全空間輸出進行預(yù)測，在優(yōu)化過程中，能大量降低模型仿真次數(shù)，對基于大計算量高精度模型的優(yōu)化設(shè)計，能顯著提高優(yōu)化設(shè)計的效率。

本文將基于代理模型的高效優(yōu)化方法引入固體火箭發(fā)動機裝藥設(shè)計過程中，為進一步提升優(yōu)化性能，提出并行加點準(zhǔn)則。最后，以后翼柱型裝藥為例，驗證本文方法的有效性。

1　優(yōu)化計算模型

1.1基于LevelSet的燃面計算方法

燃面計算是固體火箭發(fā)動機設(shè)計的重要環(huán)節(jié)，其方法優(yōu)劣直接影響發(fā)動機內(nèi)彈道性能的預(yù)估及發(fā)動機裝藥設(shè)計品質(zhì)。燃面推移是一個復(fù)雜的過程，推移過程中伴有燃面的消失和新燃面的生成，導(dǎo)致燃面結(jié)構(gòu)劇烈變化[6]。

代理模型技術(shù)通過少量樣本點的輸入輸出特性，對整個設(shè)計空間響應(yīng)進行預(yù)測，將其引入優(yōu)化設(shè)計過程，可有效減少模型調(diào)用次數(shù)。對基于復(fù)雜耗時模型的優(yōu)化設(shè)計，可顯著提高設(shè)計效率。本文將代理模型技術(shù)引入基于LevelSet燃面計算的固體火箭發(fā)動機裝藥優(yōu)化設(shè)計，并提出高效的并行加點準(zhǔn)則，進一步提高優(yōu)化設(shè)計效率。

1.2彈道計算方法

零維內(nèi)彈道計算的基本微分方程組如下：

(1)

(2)

(4)

1.3優(yōu)化目標(biāo)與約束

以往的裝藥設(shè)計目標(biāo)通常是發(fā)動機性能最優(yōu)，而不關(guān)心推力隨時間的變化規(guī)律。因此，在裝藥燃燒過程中，往往會出現(xiàn)由于燃面的消失和新燃面的生成導(dǎo)致推力發(fā)生震蕩的現(xiàn)象，而這往往是彈/箭總體設(shè)計中應(yīng)盡量避免的。因此，本文將發(fā)動機工作過程中的推力平穩(wěn)性作為優(yōu)化目標(biāo)，同時將發(fā)動機的性能作為約束進行優(yōu)化設(shè)計，以期滿足彈/箭總體設(shè)計需求。具體實施過程中，本文將設(shè)計發(fā)動機推力與總體要求平均推力之間的偏差最小化作為優(yōu)化目標(biāo)，其數(shù)學(xué)表達(dá)式如式(5)所示。其中，T為發(fā)動機工作時間，考慮到發(fā)動機點火和耗盡時不可避免會出現(xiàn)推力迅速變化的情況，本文略去開始工作后5%和工作結(jié)束前10%的推力變化影響。

(5)

其中，x為設(shè)計變量，即裝藥的幾何構(gòu)型參數(shù)。設(shè)計約束為裝藥性能參數(shù)，包括馬赫數(shù)M、裝填系數(shù)J、喉通比η以及裝藥質(zhì)量mp。

上述參數(shù)均為裝藥幾何構(gòu)型參數(shù)的函數(shù)。因此，本文構(gòu)建的優(yōu)化模型如式(6)所示。其中，M0、J0、η0和mp0為總體設(shè)計給出的各參數(shù)臨界值。

min:f(x)

J(x)>J0

(6)

2　基于代理模型的優(yōu)化設(shè)計方法

LevelSet對復(fù)雜三維藥型燃面計算有較高精度。然而，由于計算相對較為耗時，目前對基于該方法的裝藥設(shè)計方法僅有少量研究。本文將代理模型技術(shù)引入裝藥設(shè)計過程中，可有效減少LevelSet燃面計算模型調(diào)用次數(shù)，實現(xiàn)基于LevelSet的固體發(fā)動機裝藥優(yōu)化設(shè)計?；诖砟Ｐ偷膬?yōu)化設(shè)計方法主要由數(shù)據(jù)與處理、近似建模和序列加點3部分構(gòu)成，算法流程如圖1所示。

數(shù)據(jù)預(yù)處理階段進行采樣點的選擇，不同設(shè)計變量尺度的統(tǒng)一，以及在采樣點運行模型得到模型輸出，為近似建模提供訓(xùn)練樣本。

近似建模階段分為近似模型的構(gòu)造和驗證。首先，通過插值或擬合技術(shù)，根據(jù)訓(xùn)練樣本對設(shè)計空間的輸出進行預(yù)測；其次，最小化驗證樣本點預(yù)測誤差進行參數(shù)選優(yōu)與模型驗證。

采樣點更新階段的核心是構(gòu)造用于加密采樣點的優(yōu)化問題，即加點準(zhǔn)則。當(dāng)加點準(zhǔn)則確定后，采用一種高效的內(nèi)優(yōu)化算法(通常為全局優(yōu)化能力較強的智能優(yōu)化算法)，根據(jù)加密準(zhǔn)則，選擇下一采樣點。

收斂準(zhǔn)則用來終止算法迭代，高效的收斂準(zhǔn)則可大量減少不必要的計算消耗。本文采用二步收斂判定方法[16]：新點的加入對代理模型是否起到更新作用；新點的加入能否夠提高最優(yōu)點處代理模型近似性能。

該收斂判定準(zhǔn)則能有效減少復(fù)雜模型調(diào)用次數(shù)，提高優(yōu)化設(shè)計效率。

圖1　基于代理模型的優(yōu)化設(shè)計流程

2.1最優(yōu)拉丁超立方實驗設(shè)計

基于代理模型的優(yōu)化方法中，實驗設(shè)計是其前處理步驟，合理的實驗設(shè)計手段可有效選擇初始采樣點，提高近似模型對樣本點信息的利用效率[7]。由于拉丁超立方設(shè)計(Latin Hypercube Sampling, LHS)任意設(shè)計變量個數(shù)均能自由設(shè)計采樣點個數(shù),因此在安排計算機仿真實驗中應(yīng)用最廣。

LHS在采樣空間中，將每一個設(shè)計變量等分成個間隔(水平)。首先，在每個變量的N個等間隔中，隨機選取一個不重復(fù)的間隔，按均勻分布隨機產(chǎn)生一個樣本點，循環(huán)N次，即可得到N個樣本點。為了保證采樣點的空間均勻分布特性，需要對拉丁超立方進行合理設(shè)計，即優(yōu)化拉丁超立方設(shè)計。近年來，有大量學(xué)者對此進行了研究[8-14]，本文不再贅述，僅將其成果用于構(gòu)造優(yōu)化拉丁超立方設(shè)計。

為了避免不同變量之間數(shù)量級差異對近似模型精度的影響，在實驗設(shè)計之前，對所有變量進行歸一化處理；之后，所有過程均在歸一化后的N維單位超立方中進行。

2.2代理模型構(gòu)造方法

代理模型的本質(zhì)是多變量函數(shù)逼近，徑向基函數(shù)插值方法本質(zhì)上采用一元函數(shù)來對多元函數(shù)進行描述[15]，在多變量函數(shù)逼近中發(fā)揮著重要作用。徑向基插值可表示為

(7)

其中，N為采樣點個數(shù)；wk為每個樣本點的權(quán)系數(shù)；φk(rk)為基函數(shù)。其中，rk表示點x距中心點xk的歐氏距離。將所有采樣點xk的響應(yīng)值yk代入式(7)，得方程組

(8)

其中，φk(rik)為基函數(shù)，本文采用式(9)所示的高斯函數(shù)。

(9)

c為基函數(shù)的形狀參數(shù)，本文采用Wang等[16]提出的方法進行確定。為了增加徑向基函數(shù)插值對帶隨機誤差數(shù)據(jù)點的逼近能力，在求解方程組(8)時，引入正則化項，即在系數(shù)矩陣加上充分小的對角陣Λ(通常取0.001I)，使其能夠處理帶有隨機誤差的數(shù)據(jù)，其物理意義為在考慮模型對樣本點的逼近程度的同時，增加模型的光滑性。此時，徑向基插值的模型為

(10)

通過求解方程組(10)，可得系數(shù)w，將其代入式(7)，可得輸出的代理模型。

2.3并行加點準(zhǔn)則

加點準(zhǔn)則的構(gòu)造是基于代理模型優(yōu)化設(shè)計中的關(guān)鍵技術(shù)，合理的加點準(zhǔn)則不僅能有效提高代理模型的全局精度，而且能對最優(yōu)解做出快速預(yù)測，已有的加點準(zhǔn)則多為單點加點準(zhǔn)則[17-19]，效率低下。本文針對此問題，提出一種基于Pareto前沿采樣的并行加點準(zhǔn)則。

目前，單點加點準(zhǔn)則主要有探索準(zhǔn)則、開發(fā)準(zhǔn)則和平衡準(zhǔn)則，探索準(zhǔn)則側(cè)重在采樣點稀疏區(qū)域進行加點，開發(fā)準(zhǔn)則側(cè)重在代理模型的最優(yōu)解處加點，平衡準(zhǔn)則是對二者的折中。本文提出的并行加點準(zhǔn)則使加點過程同時具備探索和開發(fā)能力，具體實施如下。

構(gòu)造多目標(biāo)優(yōu)化問題(11)，其解為圖2所示的Pareto前沿。

(11)

其中，s(x)為代理模型輸出，最大化s(x)導(dǎo)致采樣點向代理模型預(yù)測的最優(yōu)區(qū)域聚集，使該準(zhǔn)則具備開發(fā)能力，d(x)為下一采樣點與歷史采樣點之間的最小距離，最大化d(x)導(dǎo)致采樣點向未采樣區(qū)域探索，使該準(zhǔn)則具備探索能力。因此，通過求解問題(11)得到的最優(yōu)解(即圖2中Pareto前沿上的點)，同時具備探索和開發(fā)能力；之后，在Pareto最優(yōu)解中，選取n個點進行下一輪迭代計算，n為并行計算節(jié)點個數(shù)。

圖2　基于Pareto前沿加點示意圖

3　測試算例分析

3.1藥柱參數(shù)化建模方法

后翼柱型裝藥構(gòu)型如圖3所示。在不考慮加工過程中倒圓倒角的情況下，圖3中尺寸、前后封頭橢球比以及后翼個數(shù)與翼寬，可唯一確定裝藥的幾何構(gòu)型。

圖3　后翼柱型裝藥幾何構(gòu)型

由于實際應(yīng)用中，上述參數(shù)彼此不能獨立變化。因此，優(yōu)化過程中，獨立參數(shù)選擇L1、L13、L12、D1、R6、翼寬Bw，其余參數(shù)固定不變。

給定裝藥配方后，裝藥密度已知，L1、D1的范圍根據(jù)裝藥體積Vp、馬赫數(shù)M及裝填系數(shù)J等性能指標(biāo)確定，如式(12)所示。為防止尺寸欠約束導(dǎo)致建模干涉，本算例不以實際尺寸為設(shè)計變量，采用比例映射方式處理其余各尺寸參數(shù)，優(yōu)化過程中，除L1、D1外，各設(shè)計變量如(13)所示。

(12)

(13)

經(jīng)過式(13)處理后的各設(shè)計變量范圍均可在0～1之間變化，且有效避免了建模干涉問題，是對此類構(gòu)型進行參數(shù)化的有效方法。

根據(jù)本文提出的方法，裝藥優(yōu)化設(shè)計計算步驟如下：

(1)采用優(yōu)化拉丁超立方實驗設(shè)計選取采樣點，在樣本點出計算目標(biāo)函數(shù)和所有約束值；

(2)通過樣本點的輸出值，根據(jù)2.2節(jié)的方法構(gòu)造目標(biāo)函數(shù)和約束的代理模型；

(3)采用NSGA-II算法，運行2.3節(jié)的并行加點，采樣點個數(shù)選擇為6；

(4)計算終止判定，是，停止計算，輸出最優(yōu)解；否，返回(2)，重新構(gòu)造代理模型。

3.2優(yōu)化設(shè)計與分析

本算例采用8片翼的后翼柱裝藥構(gòu)型，設(shè)計輸入指標(biāo)如表1所示。其中，推力、總沖、工作時間為設(shè)計指標(biāo)，裝藥相關(guān)參數(shù)由裝藥配方給出，喉徑和膨脹比為總體優(yōu)化設(shè)計給出的最優(yōu)結(jié)果。

優(yōu)化設(shè)計初始采樣點個數(shù)為48，并行計算節(jié)點數(shù)為6，最大模型仿真次數(shù)設(shè)置為300，目標(biāo)函數(shù)迭代曲線如圖4所示，當(dāng)模型調(diào)用次數(shù)為198次時，算法收斂。最優(yōu)解如表2所示。

圖4　優(yōu)化目標(biāo)迭代曲線

表2　優(yōu)化設(shè)計結(jié)果

最優(yōu)裝藥構(gòu)型對應(yīng)的內(nèi)彈道曲線和燃面-時間曲線如圖5和圖6所示。

圖5　最優(yōu)推力曲線

圖6　最優(yōu)燃面曲線

由圖5和圖6可得，推力和燃面曲線較平穩(wěn)，無明顯波動，表明優(yōu)化算法性能良好。

根據(jù)最終得到的推力曲線，與設(shè)計輸入的總體指標(biāo)進行對比，結(jié)果如表3所示。由表3可得，各相對誤差均在1%以內(nèi)，二者吻合較好。表明該方法是對固體火箭發(fā)動機裝藥設(shè)計的有效方法。

表3　計算性能結(jié)果對比

4　結(jié)論

(1)將代理模型引入固體火箭發(fā)動機裝藥設(shè)計，使基于LevelSet高精度燃面計算的應(yīng)用成為現(xiàn)實，提高了固體火箭發(fā)動機裝藥設(shè)計精度。

(2)提出了基于Pareto前沿并行加點的序列加點準(zhǔn)則，將其用于基于代理模型的裝藥設(shè)計，僅一百多次高精度仿真，即可找到最優(yōu)解。相比常規(guī)智能算法，大大減少了迭代次數(shù)。

(3)基于推力匹配的思想，對后翼柱裝藥進行優(yōu)化，完成固體火箭發(fā)動機裝藥設(shè)計，設(shè)計推力與推力需求能很好地匹配，且總體指標(biāo)也能很好的吻合，驗證了算法的有效性，豐富了固體火箭發(fā)動機裝藥設(shè)計方法。

[1]王有元. 固體火箭發(fā)動機設(shè)計[M]. 北京: 國防工業(yè)出版社, 1984.

[2]王鵬，李旭昌，徐穎軍. 固體火箭發(fā)動機總體優(yōu)化設(shè)計[J]. 火箭推進, 2007, 33(4): 16-19.

[3]高波，葉定友，侯曉，等. 固體火箭發(fā)動機三維藥柱燃面自動化設(shè)計方法研究[J]. 推進技術(shù), 1998, 19(6): 24-26.

[4]楊涓，戎海武，朱燕堂，等. 固體發(fā)動機藥柱優(yōu)化設(shè)計[J]. 推進技術(shù), 1997, 18(2): 98-102.

[5]Alexander I J, Forrester, Keane A J. Recent advances in surrogate-based optimization[J]. Progress in Aerospace Sciences, 2009, 45(1-3): 50-79.

[6]費陽. 含裂紋固體火箭發(fā)動機性能分析[D]. 長沙: 國防科技大學(xué), 2010.

[7]武澤平. 序列近似方法及其應(yīng)用研究[D]. 長沙: 國防科技大學(xué), 2013.

[8]劉曉路，陳英武，荊顯榮，等. 優(yōu)化拉丁方試驗設(shè)計方法及其應(yīng)用[J]. 國防科技大學(xué)學(xué)報, 2011, 33(5): 73-77.

[9]孫法省. 模型未知下試驗設(shè)計的構(gòu)造[D]. 南開大學(xué), 2010.

[10]劉新亮，郭波. 基于改進ESE算法的多目標(biāo)優(yōu)化試驗設(shè)計方法[J]. 系統(tǒng)工程與電子技術(shù), 2010, 32(2): 410-414.

[11]劉新亮，郭波. 適用于復(fù)雜系統(tǒng)仿真試驗的試驗設(shè)計方法[J]. 國防科技大學(xué)學(xué)報, 2009, 31(6): 95-99.

[12]Sun F, Liu M, Lin D K J. Construction of orthogonal Latin hypercube designs[M]. Biometrika, 2009.

[13]Grosso A, Jamali A, Locatelli M. Finding maximin Latin hypercube designs by Iterated Local Search heuristics[J]. European Journal of Operational Research, 2009, 197(2): 541-547.

[14]Cioppa T M, Lucas T W. Efficient nearly orthogonal and space-filling latin hypercubes[J]. Technometrics, 2007, 49(1): 45-55.

[15] 穆雪峰. 多學(xué)科設(shè)計優(yōu)化代理模型技術(shù)的研究和應(yīng)用[D]. 南京航空航天大學(xué), 2004.

[16]Wang D, Hu F, Ma Z, et al. A CAD/CAE integrated framework for structural design optimization using sequential approximation optimization[J]. Advances in Engineering Software, 2014, 76: 56-68.

[17]Jones D R. A taxonomy of global optimization methods based on response surfaces[J]. Journal of Global Optimization, 2001, 23: 345-383.

[18] 武澤平，王東輝，楊希祥，等. 應(yīng)用徑向基代理模型實現(xiàn)序列自適應(yīng)在采樣優(yōu)化策略[J]. 國防科技大學(xué)學(xué)報, 2014, 36(6): 18-24.

[19]Wang D, Wu Z, Fei Y, et al. Structural design employing a sequential approximation optimization approach[J]. Computers & Structures, 2014, 134: 75-87.

(編輯：崔賢彬)

Surrogate based grain design optimization for solid rocket motor

WU Ze-ping, WANG Dong-hui, HU Fan, ZHANG Wei-hua

(College of Aerospace Science and Engineering, National University of Defense Technology, Changsha410073, China)

An optimization design method for solid rockot motor grain based on surrogate model was proposed. To enhance the efficiency of surrogate based SRM grain design, a Pareto front based parallel infill criteria was proposed. Taking a finocyl grain design as an example, the variance of the thrust was minimized and a smooth thrust-time curve was obtained, which proves the effectiveness of the method.

solid rocket motor；surrogate model；parallel infill criteria；design optimization method；LevelSet；grain design

2015-04-20；

2015-06-15。

武澤平(1990—)，男，博士生，研究方向為飛行器總體設(shè)計。E-mail:zeping90315@126.com

V438

A

1006-2793(2016)03-0321-06

10.7673/j.issn.1006-2793.2016.03.005