范文鋒，許波，郝昀

(北京機(jī)電工程總體設(shè)計(jì)部，北京 100854)

0 引言

近年來(lái)，隨著臨近空間的利用逐漸受到重視，一種助推-滑翔式飛行器(boost-glide vehicle, BGV)得到國(guó)內(nèi)外學(xué)者廣泛關(guān)注[1-6]。其特點(diǎn)是通過(guò)主動(dòng)段的快速助推使飛行器達(dá)到較大速度和高度，隨后在臨近空間以無(wú)動(dòng)力跳躍滑行的方式進(jìn)行長(zhǎng)時(shí)間機(jī)動(dòng)飛行，具有可實(shí)現(xiàn)遠(yuǎn)程快速精確攻擊、覆蓋區(qū)域大、機(jī)動(dòng)性能好及突防能力較強(qiáng)等優(yōu)點(diǎn)。本文擬對(duì)助推-滑翔飛行器彈道最優(yōu)控制問(wèn)題進(jìn)行研究。

目前求解此類彈道優(yōu)化問(wèn)題主要有間接法和直接法兩大類[7]。間接法是根據(jù)Pontryagin極小值原理將約束最優(yōu)控制問(wèn)題轉(zhuǎn)化為Hamiltonian邊值問(wèn)題(Hamiltonian boundary value problem, HBVP)進(jìn)行求解，其特點(diǎn)是解的精度較高且滿足一階最優(yōu)性必要條件，然而收斂半徑小且協(xié)態(tài)變量初值難以給定等問(wèn)題成為此類方法廣泛應(yīng)用的主要困難。直接法采用參數(shù)化方法將連續(xù)空間的最優(yōu)控制問(wèn)題離散轉(zhuǎn)化為非線性規(guī)劃(non-linear programming, NLP)問(wèn)題，進(jìn)而使用成熟的NLP求解方法得到最優(yōu)彈道，由于直接法不必推導(dǎo)一階最優(yōu)必要條件，收斂半徑大，無(wú)需協(xié)態(tài)變量的初始值，因此直接法比間接法應(yīng)用更為廣泛。

本文以助推-滑翔飛行器為研究對(duì)象，建立無(wú)量綱化的彈道動(dòng)力學(xué)模型并提煉出總體設(shè)計(jì)參數(shù)，考慮飛行中的實(shí)際約束，以射程最優(yōu)為目標(biāo)建立BGV彈道最優(yōu)控制模型；采用Radau偽譜方法將約束最優(yōu)控制問(wèn)題轉(zhuǎn)化為非線性規(guī)劃問(wèn)題，通過(guò)引入連接點(diǎn)概念處理多階段不連續(xù)問(wèn)題，并使用序列二次規(guī)劃方法(sequential quadratic programming, SQP)進(jìn)行求解；最后給出數(shù)值優(yōu)化算例，說(shuō)明彈道最優(yōu)控制設(shè)計(jì)特點(diǎn)以及文中方法的實(shí)用性和有效性。

1 彈道動(dòng)力學(xué)模型

1.1 基本假設(shè)

本文以助推-滑翔飛行器為研究對(duì)象，僅限于均勻重力場(chǎng)中的縱向飛行彈道，考慮地球形狀，忽略地球自轉(zhuǎn)及扁率的影響。

1.2 氣動(dòng)力模型

借鑒文獻(xiàn)[8]中擬合得到的高精度氣動(dòng)模型，其結(jié)果如下：

CD=0.000 68Ma4-0.014 1Ma3+0.110 21Ma2-

0.408 38Ma+0.822 73α+0.174 92，

(1)

CL=2.423 9×10-6Ma5-1.895×10-5Ma4+

1.863 910-5Ma3-0.000 48Ma2+0.003 94Ma-

0.000 93α2+α+69.656.

(2)

1.3 大氣密度模型

假設(shè)大氣密度遵循指數(shù)變化規(guī)律，其表達(dá)式為

ρ=ρ0exp-βH,

(3)

式中：ρ0為海平面大氣密度，1.225 kg/m3；β=1/7 200。

1.4 無(wú)量綱化動(dòng)力學(xué)模型

當(dāng)發(fā)動(dòng)機(jī)工作結(jié)束后開(kāi)始被動(dòng)段飛行，此時(shí)僅考慮動(dòng)力學(xué)方程組中前4項(xiàng)即為被動(dòng)段動(dòng)力學(xué)方程組。

綜上可知：助推-滑翔飛行器的主要設(shè)計(jì)參數(shù)可歸納為tb，μt，Kfg，Kgs，CD和CL，由于tb=Ispμt/Kfg，因此最大射程可表示為

Lmax=fIsp,μt,Kfg,Kgs,CD,CL.

(5)

2 彈道最優(yōu)控制模型

2.1 約束條件

助推-滑翔飛行器飛行彈道經(jīng)歷嚴(yán)酷的力、熱環(huán)境，需要考慮氣動(dòng)加熱、結(jié)構(gòu)承載以及姿態(tài)穩(wěn)定的需求，提出如下過(guò)程約束以確保飛行正常。

駐點(diǎn)熱流密度約束：

(6)

總過(guò)載約束：

(7)

動(dòng)壓約束：

(8)

此外，為保證對(duì)目標(biāo)的打擊效果，對(duì)落點(diǎn)狀態(tài)提出如下終端約束：

Θtf=Θf,

(9)

(10)

2.2 BGV彈道最優(yōu)控制問(wèn)題

在以上工作基礎(chǔ)上，本文研究的助推-滑翔飛行器彈道最優(yōu)控制問(wèn)題可以描述為求解飛行攻角α，使得如下目標(biāo)函數(shù)最大：

(11)

且滿足狀態(tài)方程組(4)、過(guò)程約束(6)～(8)以及終端約束(9)～(10)。

3 Radau偽譜方法

偽譜法也稱正交配點(diǎn)法，主要包括Legendre偽譜法(Legendre pseudo-spectral method, LPM)、Gauss偽譜法(Gauss pseudo-spectral method, GPM)以及Radau偽譜法(Radau pseudo-spectral method, RPM)。這些方法相同之處都采用Lagrange多項(xiàng)式對(duì)狀態(tài)變量和控制變量進(jìn)行全局逼近，區(qū)別在于配點(diǎn)的選取不同。文獻(xiàn)[9]對(duì)以上3種偽譜方法特點(diǎn)進(jìn)行了詳細(xì)討論，指出采用GPM 和RPM得到的NLP問(wèn)題KKT(Karush-Kuhn-tucker)條件與原HBVP一階最優(yōu)性必要條件等價(jià)。由于RPM配點(diǎn)包含初始點(diǎn)，容易得到協(xié)態(tài)變量高精度初始值，本文采用RPM對(duì)助推-滑翔飛行器彈道最優(yōu)控制問(wèn)題進(jìn)行求解。

3.1 基本原理

Radau偽譜法將約束最優(yōu)控制問(wèn)題的狀態(tài)變量和控制變量在Legendre-Gauss-Radau(LGR)配點(diǎn)處進(jìn)行離散，采用全局插值的Lagrange多項(xiàng)式對(duì)狀態(tài)變量和控制變量進(jìn)行全局逼近，通過(guò)對(duì)多項(xiàng)式求導(dǎo)來(lái)近似狀態(tài)方程中的導(dǎo)數(shù)，且在一系列配點(diǎn)上滿足控制方程中右函數(shù)的約束，從而將狀態(tài)方程轉(zhuǎn)化為代數(shù)形式的等式約束；目標(biāo)函數(shù)中的積分項(xiàng)由Gauss-Radau積分近似計(jì)算；至此可將約束最優(yōu)控制問(wèn)題轉(zhuǎn)化為非線性規(guī)劃問(wèn)題[10]。

3.2 連續(xù)最優(yōu)控制問(wèn)題的離散轉(zhuǎn)換

3.2.1 連續(xù)最優(yōu)控制問(wèn)題的一般格式

不失一般性，考慮如下Bolza最優(yōu)控制問(wèn)題。求解控制函數(shù)ut及相應(yīng)的狀態(tài)函數(shù)xt，使得如下目標(biāo)函數(shù)最小：

(12)

滿足狀態(tài)方程：

(13)

過(guò)程約束：

gxt,ut,t≤0，

(14)

邊界約束：

ψxt0,t0,xtf,tf=0.

(15)

3.2.2 離散轉(zhuǎn)換的主要內(nèi)容

(1) 時(shí)域轉(zhuǎn)換

Radau偽譜法的求解時(shí)域?yàn)?[-1,1]，因此需要將實(shí)際求解時(shí)域[t0,tf]進(jìn)行轉(zhuǎn)化：

(16)

(2) 狀態(tài)變量與控制變量的近似

Radau偽譜法采用全局正交的Lagrange插值多項(xiàng)式分別對(duì)狀態(tài)變量與控制變量進(jìn)行逼近。

狀態(tài)變量的近似公式為

(17)

(18)

式中：τjj=1,…,N為L(zhǎng)GR配點(diǎn)，對(duì)應(yīng)LN-1τ+LNτ的根，LNτ表示N階Legendre多項(xiàng)式；τN+1不屬于LGR配點(diǎn)但參與狀態(tài)變量離散近似過(guò)程，其取值為1。

控制變量只在LGR配點(diǎn)處進(jìn)行離散近似，其近似公式為

(19)

(20)

(3) 狀態(tài)方程轉(zhuǎn)換

將式(17)對(duì)τ進(jìn)行微分，可得到LGR配點(diǎn)處的狀態(tài)變量導(dǎo)數(shù)為

(21)

式中：Dkj為N×N+1階Radau偽譜微分矩陣分量。

(22)

在此基礎(chǔ)上，可將LGR配點(diǎn)處狀態(tài)方程轉(zhuǎn)化為代數(shù)形式的等式約束：

k=1,2，…,N.

(23)

(4) 積分項(xiàng)的近似

最優(yōu)控制問(wèn)題中的積分項(xiàng)采用Gauss-Radau積分進(jìn)行處理。不失一般性，其轉(zhuǎn)換格式為

(24)

式中：ωk為Gauss-Radau積分權(quán)重。

3.2.3 轉(zhuǎn)換后的非線性規(guī)劃問(wèn)題

基于以上離散處理可將最優(yōu)控制問(wèn)題式(12)～(15)轉(zhuǎn)換為如下形式：

J=HxNτ1,τ1,xNτN+1,τN+1+

(25)

滿足約束：

k=1,2，…,N.

顯然式(25)，(26)構(gòu)成一個(gè)非線性規(guī)劃問(wèn)題，即求解離散點(diǎn)處的狀態(tài)變量xNτk(k=1,2，…,N+1)、LGR配點(diǎn)處的控制變量uNτk(k=1,2，…,N)以及終端時(shí)刻tf(若未給定)，使式(25)在滿足式(26)的約束下最小，可采用成熟的SQP方法求解。

3.3 多階段不連續(xù)問(wèn)題處理

對(duì)于助推-滑翔飛行器，主動(dòng)段及被動(dòng)段動(dòng)力學(xué)模型存在差別，其彈道優(yōu)化本質(zhì)上是一個(gè)多階段不連續(xù)最優(yōu)控制問(wèn)題。通過(guò)引入連接點(diǎn)的概念形成多階段偽譜方法[11-12]，即根據(jù)不同階段動(dòng)力學(xué)模型及約束條件變化特點(diǎn)將最優(yōu)控制問(wèn)題劃分成若干子區(qū)間，在每個(gè)子區(qū)間上分別對(duì)狀態(tài)變量和控制變量進(jìn)行離散配點(diǎn)，同時(shí)對(duì)相鄰子區(qū)間連接點(diǎn)加入約束條件，以處理多階段不連續(xù)的最優(yōu)控制問(wèn)題。

假設(shè)在某一點(diǎn)τ1∈[τ0,τ2]存在不連續(xù)，記

則連接點(diǎn)約束條件為

φ(x-(τ1),x+(τ1))≤0.

(27)

當(dāng)τ1為不定時(shí)，可將其作為優(yōu)化變量，同時(shí)加入不等式約束條件：

τ1-τ0τ1-τ2≤0.

(28)

4 數(shù)值優(yōu)化算例

優(yōu)化過(guò)程中根據(jù)計(jì)算精度需求及彈道特點(diǎn)選取配點(diǎn)數(shù)，其中主動(dòng)助推段LGR配點(diǎn)數(shù)取75個(gè)，被動(dòng)滑翔段LGR配點(diǎn)數(shù)取105個(gè)。總體設(shè)計(jì)參數(shù)取值如表1所示，過(guò)程約束及終端約束指標(biāo)如表2所示。以助推-滑翔飛行器總射程最大為目標(biāo)進(jìn)行優(yōu)化，結(jié)果如圖1～7所示。

表1 總體設(shè)計(jì)參數(shù)值Table 1 Parameters for overall design

表2 約束指標(biāo)Table 2 Specifications for constraints

圖1 縱向平面彈道Fig.1 Planar trajectory

圖2 當(dāng)?shù)貜椀纼A角Fig.2 Flight path angle

圖3 飛行攻角Fig.3 Angle of attack

圖4 總速度與動(dòng)壓Fig.4 Total velocity and dynamic pressure

優(yōu)化結(jié)果表明：各狀態(tài)及控制變量曲線過(guò)渡光滑，且滿足過(guò)程及終端約束。由圖3飛行攻角曲線可知：采用2次攻角轉(zhuǎn)彎的策略進(jìn)行主動(dòng)段飛行程序設(shè)計(jì)，有利于精細(xì)調(diào)節(jié)主動(dòng)段能量損失；被動(dòng)滑翔段采用波浪式攻角變化規(guī)律，而非采用最大升阻比狀態(tài)定常攻角飛行策略，更有利于提高射程。

由圖6可知與飛行距離相關(guān)的協(xié)態(tài)變量值恒為-1，由圖7可知Hamiltonian函數(shù)在主動(dòng)助推段和被動(dòng)滑翔段均保持為常值，基于助推-滑翔飛行器彈道最優(yōu)控制模型特點(diǎn)，可知圖6～7結(jié)果從側(cè)面說(shuō)明了采用的多階段Radau偽譜方法得到的KKT條件與原HBVP一階最優(yōu)性必要條件等價(jià)，計(jì)算結(jié)果與理論最優(yōu)解一致。

圖5 駐點(diǎn)熱流密度與總過(guò)載Fig.5 Heat flux at the stagnation point and total overload

圖6 飛行距離相關(guān)的協(xié)態(tài)變量Fig.6 Co-state about flight range

圖7 哈密頓函數(shù)Fig.7 Time history for Hamiltonian

5 結(jié)論

本文從多階段多約束最優(yōu)控制的角度對(duì)助推-滑翔飛行器彈道最優(yōu)控制問(wèn)題進(jìn)行研究，得到結(jié)論如下：

(1) 通過(guò)無(wú)量綱化方法得到助推-滑翔飛行器無(wú)量綱化動(dòng)力學(xué)方程組，可提取出Isp，μt，Kfg，Kgs，CD和CL作為BGV總體設(shè)計(jì)的主要參數(shù)。

(2) 建立了助推-滑翔飛行器彈道最優(yōu)控制模型，采用Radau偽譜方法將最優(yōu)控制問(wèn)題轉(zhuǎn)化為NLP問(wèn)題并使用SQP方法求解；通過(guò)數(shù)值優(yōu)化算例驗(yàn)證了模型及方法實(shí)用有效，優(yōu)化結(jié)果滿足一階最優(yōu)性必要條件。

(3) 助推-滑翔飛行器彈道包括主動(dòng)助推段和被動(dòng)滑翔段2部分，其彈道優(yōu)化是典型的多階段不連續(xù)最優(yōu)控制問(wèn)題；通過(guò)引入連接點(diǎn)的概念處理多階段不連續(xù)問(wèn)題，可避免以往彈道分段優(yōu)化的缺陷，以全射程最優(yōu)直接處理助推-滑翔飛行器彈道最優(yōu)控制問(wèn)題。

參考文獻(xiàn)：

[1] 劉欣, 楊濤, 張青斌. 助推-滑翔導(dǎo)彈彈道優(yōu)化與總體參數(shù)分析[J].彈道學(xué)報(bào), 2012, 24(3): 43-48.

LIU Xin, YANG Tao, ZHANG Qing-bin. Trajectory Optimization and Parameter Analysis for Boost-Glide Missile[J]. Journal of Ballistics, 2012, 24(3): 43-48.

[2] 李瑜, 楊志紅, 崔乃剛. 助推-滑翔導(dǎo)彈最大射程優(yōu)化[J]. 彈道學(xué)報(bào), 2008, 20(4): 53-56.

LI Yu, YANG Zhi-hong, CUI Nai-gang. Optimization of Maximum Range for Boost-Glide Missile[J]. Journal of Ballistics, 2008, 20(4): 53-56.

[3] 李珂, 聶萬(wàn)勝, 馮必鳴. 助推-滑翔飛行器彈道分段研究[J]. 指揮控制與仿真, 2012, 34(5): 21-25.

LI Ke, NIE Wan-sheng, FENG Bi-ming. Research on Multi-phase Trajectory Optimization for Boost-Glide Vehicle[J]. Command Control & Simulation, 2012, 34(5): 21-25.

[4] 劉欣, 楊濤. 滑翔導(dǎo)彈再入拉起段彈道優(yōu)化與制導(dǎo)[J].國(guó)防科技大學(xué)學(xué)報(bào), 2012, 34(1): 67-71.

LIU Xin, YANG Tao. Trajectory Optimization and Guidance in Reentry Phase for Glide Missile[J]. Journal of National University of Defense Technoloty, 2012, 34(1): 67-71.

[5] 王晨曦, 李新國(guó). 助推-滑翔導(dǎo)彈射程管理技術(shù)研究[J]. 固體火箭技術(shù), 35(2): 143-147.

WANG Chen-xi, LI Xin-guo.Range Management Technology Research of Boost Glide Missiles[J]. Journal of Solid Rocket Technology, 35(2): 143-147.

[6] 劉君, 陳克俊, 謝愈, 等. 助推滑翔飛行器發(fā)射諸元計(jì)算方法研究[J]. 彈箭與制導(dǎo)學(xué)報(bào), 2012, 32(6): 33-36.

LIU Jun, CHEN Ke-jun, XIE Yu, et al. The Research on Computing Method for Firing Data of Boost-Glide Aerocraft[J]. Journal of Projectiles, Rocket, Missiles and Guidance, 2012, 32(6): 33-36.

[7] 楊秀霞, 張毅, 施建洪, 等. 助推-滑翔飛行器軌跡設(shè)計(jì)研究綜述[J]. 海軍航空工程學(xué)院學(xué)報(bào), 2012, 27(3): 245-252.

YANG Xiu-xia, ZHANG Yi, SHI Jian-hong, et al. A Survey of Boost-Glide Aircraft Trajectory Design[J]. Journal of Naval Aeronautical and Astronautical University, 2012, 27(3): 245-252.

[8] 宗群, 田柏苓, 竇立謙. 基于Gauss偽譜法的臨近空間飛行器上升段軌跡優(yōu)化[J]. 宇航學(xué)報(bào), 2010, 31(7): 1775-1781.

ZONG Qun, TIAN Bai-ling, DOU Li-qian. Ascent Phase Trajectory Optimization for Near Space Vehicle Based on Gauss Pseudo-Spectral Method[J]. Journal of Astronautics，2010, 31(7): 1775-1781.

[9] Divya Garg. Advances in Global Pseudo-Spectral Methods for Optimal Control [D].Gainesville,USA: University of Florida, 2011.

[10] HAN Peng, SHAN Jia-yuan, MENG Xiu-yun. Reentry Trajectory Optimization Using a Multiple-Interval Radau Pseudo-Spectral Method[J]. Journal of Beijing Institute of Technology, 2013, 22(1): 20-27.

[11] 楊希祥, 張為華. 基于Gauss偽譜法的固體火箭上升段軌跡快速優(yōu)化研究[J]. 宇航學(xué)報(bào), 2011, 32(1): 15-21.

YANG Xi-xiang, ZHANG Wei-hua. Rapid Optimization of Ascent Trajectory for Solid Launch Vehicles Based on Gauss Pseudo-Spectral Method[J]. Journal of Astronautics,2011, 32(1):15-21.

[12] 關(guān)成啟, 陳聰. 基于Gauss偽譜法的助推-滑翔飛行器多階段約束軌跡優(yōu)化[J]. 宇航學(xué)報(bào), 2010,31(11): 2512-2518.

GUAN Cheng-qi, CHEN Cong. Multiphase Path-Constrained Trajectory Optimization for the Boost-Glide Vehicle Via the Gauss Pseudo-Spectral Method[J]. Journal of Astronautics， 2010, 31(11): 2512-2518.