彭 科，胡 凡，張為華

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

基于序列近似優(yōu)化算法的固體運載火箭彈道設計

彭 科，胡 凡，張為華

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

采用打靶法研究固體運載火箭彈道優(yōu)化設計問題。針對打靶法中智能優(yōu)化算法效率、優(yōu)化問題中約束條件提法等方面的不足，引入序列近似優(yōu)化算法，提出算法中徑向基函數(shù)高斯核寬度的高效確定方法。數(shù)值仿真實驗表明，提出的方法可在基本不帶來計算量增加的前提下，顯著提高代理模型精度，使序列近似優(yōu)化算法得到改進，提高優(yōu)化效率。設計固體運載火箭飛行程序。建立彈道優(yōu)化設計問題數(shù)學模型，并將模型中相關等式約束合理轉化為不等式約束，降低優(yōu)化問題求解難度；基于改進后的序列近似優(yōu)化算法完成某固體運載火箭彈道優(yōu)化，原始計算模型調用308次之后便搜索到最優(yōu)解，較傳統(tǒng)智能優(yōu)化算法顯著提高了優(yōu)化效率，優(yōu)化方案末助推級液體推進劑消耗比原方案減少25.7%。

固體運載火箭；彈道優(yōu)化；打靶法；序列近似優(yōu)化；徑向基函數(shù)；高斯核寬度

0 引言

固體運載火箭具有響應速度快、機動性強、成本低、可靠性高的特點，是各國研制的熱點[1]。運載能力是火箭發(fā)射最大有效載荷的能力，是描述運載火箭性能的關鍵指標[1]。固體火箭運載能力相對較小，運載能力受飛行彈道影響較大，開展彈道優(yōu)化設計研究對提高其運載能力和降低發(fā)射成本具有重要意義[2]。

運載火箭彈道優(yōu)化是一類終端時刻自由、終端狀態(tài)固定且?guī)в新窂郊s束的多階段、非線性最優(yōu)控制問題[3]，求解該問題的數(shù)值方法一般分為間接法和直接法[4]。間接法存在收斂域小、難以估計共軛變量初值等不足，且推導過程較為復雜[3]；直接法采用參數(shù)化方法，將連續(xù)空間的最優(yōu)控制問題轉化為非線性規(guī)劃問題進行數(shù)值求解，在飛行器彈道優(yōu)化領域得到廣泛應用[4]。根據(jù)參數(shù)化方法不同，直接法可分為同時離散控制變量和狀態(tài)變量的配點法[5]和僅離散控制變量的打靶法，配點法難以保證獲得的解是原最優(yōu)控制問題的解[4]。本文采用打靶法完成固體運載火箭彈道優(yōu)化。

固體運載火箭彈道優(yōu)化問題存在多個不等式約束和等式約束，對打靶法中優(yōu)化算法的全局尋優(yōu)能力和精度提出嚴峻挑戰(zhàn)。已有研究中主要采用智能優(yōu)化算法[6-7]、智能優(yōu)化與局部尋優(yōu)結合的算法[8-9]、智能優(yōu)化與非線性方程組迭代求解方法結合的算法[10-12]?，F(xiàn)有優(yōu)化方法在減少計算量、提高全局優(yōu)化能力、減輕對優(yōu)化初值依賴等方面尚有較大改進空間。序列近似優(yōu)化算法(Sequential Approximate Optimization，SAO)基于少量初始采樣點構造初始代理模型，采取一定的加點策略更新采樣點，逐步提高代理模型對最優(yōu)解的近似精度，根據(jù)收斂準則終止算法并輸出最優(yōu)解[13-15]。相比智能優(yōu)化算法，SAO算法可在保證全局最優(yōu)的前提下大幅降低原模型調用次數(shù)，顯著提高優(yōu)化效率[14]。代理模型構造是SAO算法的關鍵技術，提高其精度對改善SAO算法效率具有顯著作用[16]。代理模型構造基本思想是在一定的精度條件下，用一個相對簡單的逼近函數(shù)近似替代復雜高精度性能分析模型，利用已知點(樣本)的響應信息來預測未知點響應值，國內外學者對代理模型相關方法進行了大量研究[13, 17-19]。徑向基函數(shù)在精度和魯棒性方面皆較為可靠[20]，是SAO算法中被廣泛使用的代理模型，其核寬度的合理確定對近似精度有決定性影響。文獻[13]提出了一種采樣點個數(shù)趨于無窮情況下的非均勻核寬度確定方法；文獻[14]提出了基于采樣點局部密度的核寬度確定方法，但算法中總影響體積這一關鍵參數(shù)依賴經(jīng)驗確定，難以實現(xiàn)最優(yōu)的近似精度；文獻[18]解決了采樣點較為均勻條件下的徑向基函數(shù)核寬度確定問題。序列近似優(yōu)化過程中樣本點數(shù)量由少逐漸增多，樣本分布較不均勻，建立不同規(guī)模、非均勻樣本的徑向基函數(shù)(radial basis function, RBF)核寬度確定方法，對提高序列近似優(yōu)化過程代理模型精度、提高優(yōu)化效率具有重要意義。此外，運載火箭彈道優(yōu)化問題中等式約束的存在使得可行域范圍較小，往往需要對等式約束進行特殊處理，處理策略尚無一般規(guī)則，是一個有待研究的難題[6]。將等式約束轉化為不等式約束，從而改進優(yōu)化問題描述方式，可顯著利于運載火箭彈道優(yōu)化問題的求解。

針對現(xiàn)有彈道優(yōu)化方法的不足，引入SAO算法，提出算法中徑向基函數(shù)高斯核寬度的高效確定方法，提高代理模型精度；設計固體運載火箭飛行程序；構建彈道優(yōu)化設計問題數(shù)學模型，并將其中相關等式約束合理轉化為不等式約束，降低優(yōu)化問題求解難度；基于SAO算法完成某固體運載火箭彈道優(yōu)化，并分析結果。

1 序列近似優(yōu)化算法及其改進

為方便討論，給出一般優(yōu)化問題的數(shù)學表述：

(1)

本文SAO算法的流程如圖1所示[14]，主要步驟：(1)采用式(2)將設計空間線性映射到n維單位立方體內，采用拉丁超立方法得到初始采樣點，并計算初始采樣點目標函數(shù)與約束值；(2)基于徑向基函數(shù)構建目標函數(shù)與各約束值的代理模型；(3)基于代理模型，采用自適應采樣策略更新樣本點；(4)收斂判斷。

(2)

式中XiU、XiL分別為第i維設計變量的上下界；Xi、xi分別為原設計空間與映射后單位立方體中第i維設計變量的取值。

樣本點更新的自適應采樣策略取文獻[21]中的方法，該策略通過求解以下優(yōu)化問題實現(xiàn)，將最優(yōu)解作為新采樣點加入樣本中。

(3)

(4)

(i,j=1,2,…,N,i≠j)

(5)

式中N為采樣點個數(shù)；xi、xj為樣本點的坐標。

采樣初期樣本點較少，δ數(shù)值較大，使得新采樣點向未知區(qū)域探索，隨著樣本點數(shù)量的逐漸增多，δ數(shù)值逐漸減小，逐漸趨向最優(yōu)區(qū)域搜索，即采用式(5)所示的最小距離約束可使得根據(jù)優(yōu)化問題(3)得到的新采樣點具有自適應特性[21]。

本文采用粒子群算法求解式(3)所示優(yōu)化問題，收斂判斷準則采用文獻[14]中的方法。提出新的徑向基函數(shù)核寬度確定方法，提高代理模型精度，改進SAO算法，提高優(yōu)化效率。

給定N個訓練樣本S:[xi,yi](i=1,2,…,N)，基本徑向基函數(shù)的數(shù)學模型為[22]：

(6)

式中x為設計變量向量；xi是第i個樣本點的位置向量；ri=x-xi為歐式距離；wi為第i個基函數(shù)的權系數(shù)；φ為基函數(shù)。

本文取基函數(shù)為Gauss函數(shù)，即：

(7)

式中σi為第i個基函數(shù)的核寬度。

將N個訓練樣本代入式(6)，即得包含N個未知數(shù)wi(i=1,2,…,N)的N維線性方程組，求解該方程組即可求得基函數(shù)的權系數(shù)wi，本文采用QR分解法求解。

式(7)中基函數(shù)核寬度σi的確定對代理模型預測精度有決定性影響，文獻[14]提出了基于樣本點局部密度并估計總影響體積的核寬度σi計算方法，局部密度定義如下：

(8)

其中

(9)

文獻[14]中，總影響體積這一關鍵參數(shù)較難確定，難以實現(xiàn)較優(yōu)的核寬度估計。本文在樣本點局部密度計算的基礎上，提出新的核寬度確定方法。

核寬度σi表征了第i個樣本點影響區(qū)域的大小[14]，σi的n次方應與密度函數(shù)ρ(xi)成反比，即

(10)

式(10)包含N-1個獨立的方程，待解的核寬度σi為N個，求得任意樣本點核寬度后即可求解其余樣本點的核寬度。

考慮局部密度最小的樣本點xs，設其核寬度為σs，其基函數(shù)在距xs最近的樣本點的數(shù)值為

(11)

若σs數(shù)值過小，其余樣本點的核寬度σi也將過小，將導致徑向基函數(shù)代理模型不光滑；反之若σs數(shù)值過大，將導致龍格現(xiàn)象。合理的σs取值應保證樣本點xs的基函數(shù)的影響域應到達離其最近的樣本點的位置，即φ(ds,min)的數(shù)值應足夠大，φ(ds,min)與σs/ds,min的數(shù)值關系如圖2所示。

當σs/ds,min<0.659 2時，φ(ds,min)<0.1，數(shù)值偏小，即樣本點xs的影響偏弱；當σs/ds,min>1時，φ(ds,min)>0.376 9。本文取σs/ds,min=1，以保證樣本點xs的影響域到達離其最近的樣本點，確保代理模型光滑，即

σs=ds,min

(12)

由式(10)與式(12)得，各樣本點核寬度σi(i=1,2,…,N)計算如下：

(13)

采用下式準則比較本文提出的核寬度確定方法和文獻[13]與文獻[23]提出的方法對應徑向基函數(shù)的近似精度：

(14)

R2≤1，越接近于1表明近似模型精度越高;R2=1時表明近似模型在驗證樣本處的誤差為0。

取如下所示5個典型測試函數(shù)。

函數(shù)I(一維函數(shù))：

f(x)=xsin(1.5x)，0≤x≤10

函數(shù)II(低維低階函數(shù))：

f(x)= sin(x1+x2)+(x1-x2)2-1.5x1+

2.5x2+1，|xi|≤5

f(x)= (10+x1cosx1)×[3+exp(-x22)]，

|xi|≤5

函數(shù)IV(高維低階函數(shù))：

f(x)= 2(x1-1)2+x22-x1x2+(4-x3)2+3(x4-6)2+

0.5(x5-2)2+(x6-10)2+4(x7-9)2+

2(x8-5)2+33，|xi|≤20

函數(shù)V(高維高階函數(shù))：

|xi|≤4

采用文獻[24]的方法，驗證本文提出的核寬度確定方法的性能，函數(shù)I取10個隨機樣本點，函數(shù)II、III取50個隨機樣本點，函數(shù)IV、V取200個隨機樣本點，根據(jù)3種核寬度確定方法建立近似模型，為保證近似精度R2計算的準確性，取足夠多的驗證樣本點，數(shù)量為1 000，進行20次獨立的隨機數(shù)值仿真實驗，計算R2與R2的平均值MR2。3種核寬度確定方法的MR2值比較如圖3所示。結果表明，本文方法近似精度顯著優(yōu)于其他2種方法，可改善SAO算法效率和精度。

2 固體運載火箭總體參數(shù)與飛行程序設計

2.1 總體參數(shù)

適宜氣候：溫暖濕潤；年均氣溫10～20 ℃，1月份平均氣溫3～9 ℃，7月份平均氣溫24～28 ℃，極端最高氣溫低于35 ℃，極端最低氣溫高于0 ℃，年均降水量600～1 200 mm，年平均日照600～1 200 h[7]。

本文研究的固體運載火箭由三級固體發(fā)動機與液體末助推級組成，各級固體發(fā)動機參數(shù)性能參數(shù)如表1所示。末助推級液體推進劑質量為230 kg，真空比沖Isp=2 850 N·s/kg，真空推力Tb=2 000 N。全箭最大直徑為1.4 m，目標軌道為300 km太陽同步圓軌道，入軌質量約450 kg。

2.2 飛行程序設計

本文固體運載火箭采用直接入軌方式，飛行過程一級動力飛行段、一級滑行段、二級動力飛行段、二級滑行段、三級動力飛行段、三級滑行段、末助推段等7個階段[8]，飛行程序如下：

參數(shù)名稱一子級二子級三子級直徑/mm140014001200裝藥量/kg1520078002900平均比沖/(N·s/kg)2350(地面)2810(真空)2850(真空)質量比0．9130．9060．923工作時間/s66．462．257．4

(1)一級飛行段

運載火箭一級在稠密大氣層中飛行，火箭垂直起飛，而后按攻角程序轉彎，飛行速度進入跨聲速前結束轉彎[25]，而后保持零攻角飛行，飛行程序為

(15)

其中

α(t)=-αmsin2f(t)，t1

(16)

(17)

為便于處理，引入λm=(tm-t1)/(t2-t1)。

(2)二級飛行段

運載火箭真空飛行動力段最優(yōu)俯仰角程序近似為線性[25]：

(18)

運載火箭二級及以上飛行段大氣已較稀薄，飛行程序按真空飛行段設計[25]：

(19)

(3)三級飛行段

運載火箭三級飛行段飛行程序設計：

(20)

(4)末助推段

末助推段飛行程序設計：

末助推關機時刻應達到入軌條件。

3 彈道優(yōu)化設計問題數(shù)學模型

3.1 設計變量

根據(jù)飛行程序，取固體運載火箭彈道優(yōu)化設計變量為

x=[x1x2x3x4x5x6x7x8x9x10x11x12x13]T

(22)

式中A0為發(fā)射方位角。

本文取設計變量上下限如表2所示。

表2 固體運載火箭彈道優(yōu)化設計變量上下限

3.2 目標函數(shù)

本文彈道優(yōu)化設計問題目標函數(shù)的實質是燃料消耗最省。由于助推系統(tǒng)三級固體發(fā)動機均采用耗盡關機工作模式，選取末助推級消耗液體推進劑質量mp最小作為目標函數(shù)[3]，即

(23)

3.3 約束條件

本文彈道優(yōu)化設計問題約束條件包括不等式約束和等式約束。為保證火箭安全飛行，對法向過載ny、一級程序轉彎結束時馬赫數(shù)Mat2與一二級分離時動壓qsep進行約束，得如下所示不等式約束：

(24)

式中nymax為法向過載上限值，nymax=0.05；Mat2max為一級程序轉彎結束時馬赫數(shù)上限值，Mat2max=0.75；qsepmax為一二級分離時動壓上限值，qsepmax=10 kPa。

目標軌道為300 km太陽同步圓軌道，得等式約束：

(25)

式中h、v、e、i分別為實際入軌高度、入軌速度、偏心率和軌道傾角；hobj=300 km、vobj=7 725.84 m/s、eobj=0、iobj=96.67°分別為目標軌道高度、速度、偏心率和軌道傾角。

式(25)中軌道傾角i主要受設計變量中發(fā)射方位角A0影響，優(yōu)化過程中Ceq4(x)較易滿足。Ceq1(x)、Ceq2(x)、Ceq3(x)對大多數(shù)設計變量敏感，優(yōu)化過程難以滿足，約束條件的處理難度較大，優(yōu)化問題難以求得全局最優(yōu)解。

本文建立的運載火箭彈道優(yōu)化問題不考慮式(25)中Ceq1(x)、Ceq2(x)、Ceq3(x)等3個等式約束。根據(jù)末助推關機點參數(shù)計算出軌道近地點高度hp，以hp≥hobj=300 km為約束，因目標函數(shù)為燃料消耗最省，優(yōu)化問題的最優(yōu)解必然為300 km圓軌道，即滿足式(25)中的Ceq1(x)、Ceq2(x)、Ceq3(x)等3個等式約束。得本文運載火箭彈道優(yōu)化問題約束條件如下：

(26)

4 優(yōu)化結果與分析

采用改進后的SAO算法，優(yōu)化過程彈道計算采用發(fā)射系中空間三自由度彈道方程，具體計算模型見文獻[25]。取初始樣本點為200個，優(yōu)化過程目標函數(shù)收斂曲線如圖4所示，迭代至第108步，原始計算模型調用308次之后，便搜索到全局最優(yōu)解，較傳統(tǒng)智能優(yōu)化算法上萬次的原始計算模型調用次數(shù)[6]，顯著提高了優(yōu)化效率。

優(yōu)化后彈道與優(yōu)化前的原始彈道設計變量與性能參數(shù)的對比如表3、表4。其中，原始彈道為人工調節(jié)得到的彈道方案。優(yōu)化方案末助推級液體推進劑消耗量為147.62 kg，比原始方案減少25.7%，各項設計指標與約束滿足要求。優(yōu)化彈道高度、速度、動壓、俯仰角、當?shù)厮俣葍A角隨時間變化曲線如圖5、圖6所示。

表3 固體運載火箭彈道優(yōu)化方案與原方案設計變量數(shù)值

表4 固體運載火箭彈道優(yōu)化方案與原方案對應特征參數(shù)

5 結論

(1)提出了巧妙、高效、可靠的徑向基函數(shù)高斯核寬度的確定方法?；跇颖军c局部密度計算結果，合理確定局部密度最小的樣本點基函數(shù)核寬度，進而巧妙地確定其他樣本點的核寬度。提出的核寬度確定方法在基本不帶來計算量增加的前提下，顯著提高了代理模型精度。

(2)設計了固體運載火箭飛行程序。建立了彈道優(yōu)化設計問題數(shù)學模型，并將模型中與入軌點高度、速度、當?shù)厮俣葍A角關聯(lián)的3個等式約束合理轉化為與軌道近地點高度關聯(lián)的1個不等式約束，降低了優(yōu)化問題求解難度。

(3)基于改進后的SAO算法，完成了某固體運載火箭彈道優(yōu)化，原始計算模型調用308次之后，便搜索到最優(yōu)解，較傳統(tǒng)智能優(yōu)化算法顯著提高了優(yōu)化效率。優(yōu)化方案各項指標與約束滿足設計要求，末助推級液體推進劑消耗比原方案減少25.7%。

[1] 洪蓓, 梁欣欣, 辛萬青. 固體運載火箭多約束彈道優(yōu)化[J]. 導彈與航天運載技術, 2012(3): 1-5.

[2] 龍樂豪.總體設計(上)[M].北京:宇航出版社,1989.

[3] 楊希祥, 張為華. 基于Gauss偽譜法的固體運載火箭上升段軌跡快速優(yōu)化研究[J]. 宇航學報, 2011, 31(1): 15-21.[4] 雍恩米, 陳磊, 唐國金. 飛行器軌跡優(yōu)化數(shù)值方法綜述[J]. 宇航學報, 2008, 29(2): 397-406.

[5] 龔春林, 韓璐. RBCC可重復使用運載器上升段軌跡優(yōu)化設計[J]. 固體火箭技術, 2012, 35(3): 290-295.

[6] 楊希祥, 江振宇, 張為華. 基于粒子群算法的固體運載火箭上升段彈道優(yōu)化設計研究[J]. 宇航學報, 2010, 31(5): 1304-1309.

[7] 李振華, 鮮勇, 雷剛, 等. 基于混合粒子群算法的上升段交會彈道快速優(yōu)化設計[J]. 航空動力學報, 2015, 30(12): 3029-3034.

[8] 楊希祥, 張為華, 肖飛, 等. 小型固體運載火箭運載能力分析[J]. 固體火箭技術, 2009, 32(4): 355-359.

[9] Volker Maiwald. About combining tisserand graph gravity-assist sequencing with low-thrust trajectory optimization[C]//International Conference on Astrodynamics Tools and Techniques, Darmstadt: 2016.

[10] 崔乃剛, 黃盤興, 路菲, 等. 基于混合優(yōu)化的運載器大氣層內上升段軌跡快速規(guī)劃方法[J]. 航空學報, 2015, 36(6): 1915-1923.

[11] 黃文博, 張強, 肖飛, 等. 空間快速響應航天器軌道/彈道一體化規(guī)劃[J]. 固體火箭技術, 2012, 35(1): 11-16.

[12] Abolfazl Shirazi. Trajectory optimization of spacecraft high-thrust orbit transfer using a modified evolutionary algorithm[J]. Engineering Optimization,2015: 1-19. http://dx.doi.org/10.1080/0305215X.2015.1115026.

[13] Kitayama S,Arakawa M,Yamazaki K.Sequential approximate optimization using radial basis function network for engineering optimization[J]. Optimization and Engineering,2011, 12: 535-557.

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

[15] Kitayama S, Arakawa M, Yamazaki K. Sequential approximate optimization for discrete design variable problems using radial basis function network[J]. Applied Mathematics and Computation,2012, 219(8): 4143-4156.

[16] Jacobs J H,Etman L F P,Keulen F Van,et al. Framework for sequential approximate optimization[J]. Structural And Multidisciplinary Optimization,2004, 27: 384-400.

[17] Deng Y M,Lam I C,Tor S B. A CAD-CAE integrated injection molding design system[J]. Eng. Computation,2002, 18: 80-92.

[18] Kitayama S,Yamazaki K. Simple estimate of the width in Gaussian kernel with adaptive scaling technique[J]. Applied Soft Computing,2011, 11: 4726-4737.

[19] Luo C,Zhang S,Wang C. A meta model-assisted evolutionary algorithm for expensive optimization[J]. Journal of Computational and Applied Mathematics,2011, 236: 759-764.

[20] Jin R,Chen W,Simpson W T. Comparative studies of metamodeling techniques under multiple modeling criteria[J]. Journal of Structural and Multidisciplinary Optimization,2001, 23: 1-13.

[21] 武澤平, 王東輝, 楊希祥，等. 基于徑向基代理模型的序列優(yōu)化中自適應再采樣策略[J]. 國防科技大學學報, 2014, 36(6): 18-24.

[22] Hardy R L. Multiquadratic equations of topography and other irregular surfaces[J]. Journal of Geophysical Research,1971, 76: 1905-1915.

[23] Nakayama H,Arakawa M,Sasaki R. Simulation-based optimization using computational intelligence[J]. Optimization and Engineering,2002(3): 201-214.

[24] Wang Dong-hui,Hu Fan,Ma Zhen-yu, et al. A CAD/CAE integrated framework for structural design optimization using sequential approximation optimization[J]. Advances in Engineering Software,2014, 76: 56-68.

[25] 陳克俊, 劉魯華, 孟云鶴. 遠程火箭飛行動力學與制導[M]. 北京: 國防工業(yè)出版社, 2014.

(編輯：呂耀輝)

Solid launch vehicles trajectory design based on sequential approximate optimization algorithm

PENG Ke, HU Fan, ZHANG Wei-hua

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

Solid launch vehicles trajectory optimization problem was researched using shooting method. Current shooting methods have shortcomings on the efficiency of intelligent optimization algorithms and constraints formulation. The sequential approximate optimization (SAO) algorithm was introduced to solve solid launch vehicles trajectory optimization problem. A high efficiency Gaussian kernel width determining method for radial basis function (RBF) used in SAO was proposed. Numerical experiments show that the proposed kernel width determining method improves the approximate accuracy of RBF remarkably and generates almost no excessive calculation cost. Flight program and trajectory optimization problem for solid launch vehicles were established. In order to reduce the difficulty of solving solid launch vehicles trajectory optimization problem, equality constraints in optimization problem were transformed into an inequality constraint ingeniously. Based on improved SAO algorithm, trajectory optimization for a launch vehicle was accomplished, and the global optimal solution was obtained after calculating original model 308 times. Propellant consumption of the upstage in optimization scheme is 25.7% less than primary scheme.

solid launch vehicle；trajectory optimization；shooting method；sequential approximate optimization；radial basis function；Gaussian kernel width

2016-01-27；

2016-04-20。

彭科(1989—)，男，博士生，研究方向為飛行器總體設計優(yōu)化。E-mail:pengke_pk@163.com

張為華，男，教授，博士生導師，研究方向為飛行器總體設計。E-mail:zhangweihua@nudt.edu.cn

V412.1

A

1006-2793(2017)02-0250-07

10.7673/j.issn.1006-2793.2017.02.021