張 松 侯明善 唐成師

1.西北工業(yè)大學(xué)自動(dòng)化學(xué)院，西安710072

2.上海機(jī)電工程研究所，上海201109

彈道優(yōu)化是指確定一條滿足給定動(dòng)力學(xué)方程和相關(guān)約束條件的空間運(yùn)動(dòng)軌跡，使導(dǎo)彈沿該軌跡飛行時(shí)指定的性能指標(biāo)最優(yōu)。常規(guī)數(shù)值彈道優(yōu)化方法可分為間接法和直接法［1－2］。間接法［3－5］根據(jù)最優(yōu)控制原理將優(yōu)化問題轉(zhuǎn)換為Hamilton 邊值問題進(jìn)行求解，解的最優(yōu)性能夠得到保障，但難以適應(yīng)復(fù)雜的路徑約束情況，對(duì)初始猜測(cè)要求苛刻。直接法［6－9］將連續(xù)最優(yōu)控制問題轉(zhuǎn)化為有限維非線性規(guī)劃問題進(jìn)行迭代搜索求解。相對(duì)間接法，直接法收斂域更寬，但離散化配置點(diǎn)數(shù)量大，且依賴初始猜測(cè)。對(duì)此，文獻(xiàn)［7］設(shè)計(jì)了初值軌跡生成算法，彌補(bǔ)了Gauss 偽譜法對(duì)初值敏感的不足。文獻(xiàn)［8］采用雙層優(yōu)化策略對(duì)細(xì)化單元數(shù)和插值基函數(shù)的階次進(jìn)行自適應(yīng)調(diào)節(jié)，減少了優(yōu)化變量數(shù)量并能滿足精度要求。文獻(xiàn)［9］在利用偽譜法設(shè)計(jì)最優(yōu)軌跡時(shí)，采用二階微分方程技術(shù)消去部分狀態(tài)量以減少優(yōu)化變量。此外，文獻(xiàn)［10］將逆動(dòng)態(tài)方法與直接法相結(jié)合研究了彈道導(dǎo)彈上升段彈道優(yōu)化問題。該方法通過(guò)優(yōu)化期望輸出從而間接設(shè)計(jì)控制量，有效減少了需要的離散變量個(gè)數(shù)，但該方法中期望輸出的參數(shù)化方式不具一般性，且精度難以保證。

從幾何本質(zhì)上看，飛行彈道總是一條滿足邊界條件的光滑曲線。對(duì)給定光滑曲線，曲線的切矢量和法矢量可通過(guò)方程計(jì)算確定。等價(jià)地說(shuō)，如果彈道光滑且方程已知，則彈道傾角和彈道偏角及其變化規(guī)律就是確定的，結(jié)合導(dǎo)彈運(yùn)動(dòng)方程組、約束條件等通過(guò)數(shù)值積分計(jì)算就可確定彈道上的其它參量，包括速度、法向加速度、升力和阻力等。這從另一個(gè)方面提示:為了獲得最優(yōu)彈道也可采取先形成光滑幾何彈道，再計(jì)算彈道其它參量，并根據(jù)優(yōu)化指標(biāo)對(duì)彈道做出評(píng)價(jià)，依據(jù)評(píng)價(jià)結(jié)果改變彈道幾何形狀，最終獲取最優(yōu)彈道。顯然，這種方法首先需要解決光滑彈道的形成或造型問題。作者研究發(fā)現(xiàn)，幾乎所有的光滑飛行彈道均可用有理Bezier 曲線構(gòu)造。因此，提出了一種新的彈道優(yōu)化策略，基本思想為:根據(jù)邊界條件將待求彈道表示成光滑且只含少量自由參數(shù)的有理Bezier 曲線形式，根據(jù)導(dǎo)彈運(yùn)動(dòng)方程組導(dǎo)出狀態(tài)、控制量和最優(yōu)性指標(biāo)與自由參數(shù)間的函數(shù)關(guān)系，從而將連續(xù)的彈道優(yōu)化問題轉(zhuǎn)換為對(duì)自由參數(shù)的參數(shù)優(yōu)化問題。該方法優(yōu)化變量少且所得最優(yōu)彈道光滑性好，可用于解算全程控制導(dǎo)彈，如飛航式導(dǎo)彈、防空類導(dǎo)彈等彈道優(yōu)化問題，也可用于彈道導(dǎo)彈主動(dòng)段最優(yōu)彈道設(shè)計(jì)。

1 彈道優(yōu)化問題描述

1.1 導(dǎo)彈縱向運(yùn)動(dòng)描述

導(dǎo)彈力平衡及本文所采用的符號(hào)如圖1 所示。若不考慮地球自轉(zhuǎn)，認(rèn)為重力場(chǎng)是平行力場(chǎng)，且導(dǎo)彈取為質(zhì)點(diǎn)模型，其運(yùn)動(dòng)限制在鉛垂平面內(nèi)，那么導(dǎo)彈運(yùn)動(dòng)方程組可寫為:

其中，

式(1)～(8)中，x 為導(dǎo)彈水平飛行距離，h 為導(dǎo)彈飛行高度，V 為導(dǎo)彈飛行速度，γ 為彈道傾角，T 為導(dǎo)彈發(fā)動(dòng)機(jī)推力，α 為攻角，m 為導(dǎo)彈質(zhì)量，g 為重力加速度，D 為阻力，L 為升力，ρ 為大氣密度，S 為參考面積;CL為升力系數(shù)，CLα為升力的攻角氣動(dòng)導(dǎo)數(shù)，α0為平衡點(diǎn)攻角;CD為阻力系數(shù)，CD0為零升阻力系數(shù)，K 為升阻極曲線系數(shù)。這里假定導(dǎo)彈發(fā)動(dòng)機(jī)推力T 和質(zhì)量m 是時(shí)間t 的確定函數(shù)，即fT(t)，fm(t)是已知的。

圖1 導(dǎo)彈力平衡與符號(hào)

1.2 彈道優(yōu)化問題

取狀態(tài)向量x =［x，h，V，γ］，控制變量uc= α，其容許集為 。設(shè)初始狀態(tài)為:

末端時(shí)間tf(給定或未定)時(shí)狀態(tài)滿足約束:

狀態(tài)變量和控制變量約束條件為:

設(shè)彈道優(yōu)化指標(biāo)為

彈道優(yōu)化問題可表述為:在給定邊界條件(9)和(10)的情況下，確定滿足約束條件(11)和運(yùn)動(dòng)方程組(1)～(6)的控制變量uc(t)∈ 使性能指標(biāo)(12)最優(yōu)。

2 彈道Bezier 曲線造型與優(yōu)化

2.1 有理Bezier 型彈道

有理n 次Bezier 曲線可表示為

式中，bi，μi，i = 0，1，…，n 分別為控制多邊形頂點(diǎn)和權(quán)因子。平面內(nèi)控制多邊形頂點(diǎn)可寫為

有理Bezier 曲線計(jì)算簡(jiǎn)便［11］且具有諸多利于彈道造型的優(yōu)良性質(zhì)，例如:首末端點(diǎn)分別是控制多邊形的首末頂點(diǎn)，在首末端點(diǎn)處分別與控制多邊形首末邊相切，由同一組控制多邊形和權(quán)因子定義的有理Bezier 曲線唯一。

由式(13)可知，有理Bezier 曲線是關(guān)于控制多邊形頂點(diǎn)及其權(quán)因子的函數(shù)，若xi，hi和μi為自由參數(shù)，那么通過(guò)改變自由參數(shù)的值可以得到不同形狀的有理Bezier 曲線。使用有理Bezier 曲線形成參數(shù)化彈道時(shí)，為了滿足彈道邊界條件約束，控制多邊形部分頂點(diǎn)需要附加相應(yīng)的約束。

設(shè)彈道邊界條件為:

則控制多邊形頂點(diǎn)b0，b1，bn－1以及bn應(yīng)滿足下列條件:

這樣，用于調(diào)節(jié)有理Bezier 曲線形狀的自由參數(shù)成為x1，…，xn－1，h2，…，hn－2，μ0，…，μn。進(jìn)一步，若不考慮末端彈道傾角約束(16)，選擇自由參數(shù)x2，…，xn－2為區(qū)間［xs，xf］上具有特定分布特征的點(diǎn)，比如等距分布的點(diǎn)

則由于xs和xf為已知量，自由參數(shù)簡(jiǎn)化為

因此，有理Bezier 曲線可根據(jù)式(13)表示為關(guān)于Φ 的函數(shù)，即

不妨稱P(Φ，ω)為有理Bezier 型彈道。

2.2 基于彈道造型的優(yōu)化問題描述

經(jīng)彈道造型后，彈道優(yōu)化問題就歸結(jié)為:確定自由參數(shù)Φ，使得與之對(duì)應(yīng)的彈道在某項(xiàng)性能指標(biāo)下最優(yōu)。這就需要根據(jù)自由參數(shù)和導(dǎo)彈運(yùn)動(dòng)方程組確定性能指標(biāo)。在本文方法中，有理Bezier 型彈道的位置變量x，h 及彈道傾角變量γ 完全由自由參數(shù)Φ確定，還需確定速度變量V，攻角α 以及時(shí)間變量t以計(jì)算性能指標(biāo)。

記

則彈道傾角γ 可表示為

上式兩邊對(duì)參數(shù)ω 求導(dǎo)，可得

引入變換關(guān)系

對(duì)式(1)，(3)和(4)作變換得到

考慮式(7)和(19)，整理式(23)可得攻角α 關(guān)于速度、時(shí)間和自由參數(shù)Φ 的表達(dá)式。因此，聯(lián)立式(21)和(22)，就得到時(shí)間和速度關(guān)于參數(shù)ω ∈［0，1］的微分方程組

給定自由參數(shù)Φ，通過(guò)解微分方程組(24)可以得到時(shí)間和速度變量關(guān)于參數(shù)ω 的變化規(guī)律。

因此，確定性能指標(biāo)所需的各變量，彈道優(yōu)化問題(9)～(12)就歸結(jié)為:確定自由參數(shù)Φ，使得性能指標(biāo)

最小，并滿足狀態(tài)和控制約束

其中，Φ ∈UΦ，UΦ為自由參數(shù)Φ 的取值范圍。由于邊界條件約束(14)已通過(guò)彈道造型自然滿足，狀態(tài)和控制約束一般為高度、彈道傾角和攻角約束，即

通常，αmin，γmin＜0 且αmax，γmax，hmax＞0。

本文通過(guò)引入罰函數(shù)將多約束軌跡優(yōu)化問題轉(zhuǎn)換為無(wú)約束軌跡優(yōu)化問題，并保證在給定采樣節(jié)點(diǎn)處狀態(tài)和控制滿足約束(27)，具體步驟如下:

首先，令

式中，ωi，i = 1，…，m 為區(qū)間［0，1］上等距分布的采樣節(jié)點(diǎn)，m 為足正整數(shù)，本文取m =500。將約束(27)簡(jiǎn)化為:

其次，根據(jù)簡(jiǎn)單約束(28)構(gòu)造罰函數(shù)

其中，罰因子r ＞0 為常值。

從而約束優(yōu)化問題(25)和(26)最終簡(jiǎn)化為:確定自由參數(shù)Φ，使得指標(biāo)函數(shù)

最小。其中，Φ ∈UΦ，Φ為自由參數(shù)Φ 的取值范圍。

由有理Bezier 彈道的特性可知，優(yōu)化問題(30)中優(yōu)化變量數(shù)量少，所得彈道完全光滑，而且因狀態(tài)和控制變量都是通過(guò)等價(jià)轉(zhuǎn)換和積分來(lái)確定，故各變量精確滿足導(dǎo)彈運(yùn)動(dòng)方程組，即解的精度高。

3 彈道優(yōu)化算例

本節(jié)以某型空空導(dǎo)彈為例，研究其最短時(shí)間飛行彈道。對(duì)于空空導(dǎo)彈而言，飛行時(shí)間短有利于快速打擊敵方目標(biāo)，減少載機(jī)受攻擊的機(jī)會(huì)，擴(kuò)大載機(jī)執(zhí)行其它任務(wù)的靈活性，對(duì)于采用慣導(dǎo)形式中制導(dǎo)的空空導(dǎo)彈，更能有效減小誤差。因而要提升空空導(dǎo)彈的打擊能力必須把短的飛行時(shí)間作為一個(gè)基本的性能指標(biāo)。本節(jié)取優(yōu)化指標(biāo)為

邊界條件為x(t0)= xs，h(t0)= hs，γ(t0)= γs，x(tf)= xf，h(tf)= hf。

表1 所給出了2 組邊界條件。用于造型的有理Bezier 曲線的次數(shù)取為6，如2.2 節(jié)所述，自由參數(shù)(優(yōu)化變量)為優(yōu)化變量取值范圍設(shè)為0 ＜x1≤xf/5，μj＞0，j =0，1，…，6，8≤hi≤35(km)，其中，i=2，…，5。

表1 邊界條件

仿真所用導(dǎo)彈的氣動(dòng)參數(shù)來(lái)源于文獻(xiàn)［12］，其中CLα=7，CD0=0.4，K =0.01，S =0.02m2?？諝饷芏炔捎脴?biāo)準(zhǔn)大氣模型計(jì)算，發(fā)動(dòng)機(jī)推力T 是時(shí)間的函數(shù)，關(guān)系式為:

導(dǎo)彈質(zhì)量變化規(guī)律˙m 滿足:

攻角、彈道傾角和高度約束為:α ∈［－ 15°，15°］，γ ∈［－50°，50°］，hs＜h ≤25km。

本文采用Matlab 全局優(yōu)化工具箱中的遺傳算法工具進(jìn)行搜索，參數(shù)設(shè)置如下:種群規(guī)模25，最大繁衍代數(shù)120，其余參數(shù)采用默認(rèn)設(shè)置。對(duì)條件1和2，罰函數(shù)參數(shù)r 分別取為15 和50。為了驗(yàn)證算法的有效性，本文將所提出算法的優(yōu)化結(jié)果與基于hp 自適應(yīng)偽譜法的優(yōu)化工具GPOPS［13］的所得結(jié)果進(jìn)行了對(duì)比研究。GPOPS 網(wǎng)格迭代次數(shù)為10，其它參數(shù)為默認(rèn)設(shè)置。最終的參數(shù)優(yōu)化結(jié)果如表2 和表3 所示。

表2 自由參數(shù)優(yōu)化搜索結(jié)果

表3 自由參數(shù)優(yōu)化搜索結(jié)果(續(xù))

對(duì)條件1 和2，本文算法優(yōu)化彈道的飛行時(shí)間分別為101.3s 和151.5s，而GPOPS 軟件的優(yōu)化結(jié)果分別為108.7s 和159.4s。可見本文算法得到的優(yōu)化結(jié)果更逼近全局最優(yōu)解。圖2 ～5 為條件1 兩種方法得到的優(yōu)化彈道、速度、彈道傾角和攻角曲線;圖6 ～9 為條件2 相應(yīng)的特性曲線。圖中實(shí)線表示本文算法優(yōu)化結(jié)果(ORBT)，虛線為軌跡優(yōu)化工具GPOPS 的優(yōu)化結(jié)果(GPOPS)。由圖2 和6 的彈道曲線看，2 種優(yōu)化方法得到的彈道最大高度基本一致，但ORBT 彈道爬升段平緩，這一點(diǎn)在彈道傾角特性中也可以看到。分析圖3 和7 可知，ORBT 優(yōu)化彈道的飛行速度始終大于GPOPS 優(yōu)化彈道的速度，而彈道傾角差異不大，因而ORBT 優(yōu)化彈道的水平速度較大，飛行時(shí)間更短。由圖5 和9 比較可以看出，導(dǎo)彈沿GPOPS 優(yōu)化彈道飛行時(shí)攻角易處于飽和狀態(tài)，而ORBT 彈道則很好的避免了攻角飽和問題。

圖3 速度時(shí)間歷程(條件1)

圖4 彈道傾角時(shí)間歷程(條件1)

圖5 攻角時(shí)間歷程(條件1)

圖6 鉛垂面彈道(條件2)

圖7 速度時(shí)間歷程(條件2)

圖8 彈道傾角時(shí)間歷程(條件2)

圖9 攻角時(shí)間歷程(條件2)

4 結(jié)論

彈道優(yōu)化問題是具有高度非線性、多約束和多尺度特性的最優(yōu)控制問題，常規(guī)數(shù)值優(yōu)化方法處理可靠性和全局最優(yōu)性難以保證，為此，本文提出了基于有理Bezier 曲線的彈道造型和參數(shù)優(yōu)化方法，并結(jié)合遺傳算法與模式搜索來(lái)搜索最優(yōu)參數(shù)。這種方法將無(wú)限維彈道優(yōu)化問題轉(zhuǎn)換為關(guān)于少量自由參數(shù)的參數(shù)優(yōu)化問題，使得彈道優(yōu)化求解更方便，解的收斂性和全局最優(yōu)性好。與GPOPS 優(yōu)化工具箱計(jì)算結(jié)果的對(duì)比分析驗(yàn)證了方法的可行性和有效性。

［1］Betts J T. Survey of Numerical Methods for Trajectory Optimization［J］.Journal of Guidance，Control and Dynamic，1998，21(2):193-206.

［2］雍恩米，陳磊，唐國(guó)金. 飛行器軌跡優(yōu)化數(shù)值方法綜述［J］. 宇航學(xué)報(bào)，2008，29(2):7-16. (Yong E M，Chen L，Tang G J. A Survey of Numerical Methods for Trajectory Optimization of Spacecraft［J］. Journal of Astronautics，2008，29(2):7-16.)

［3］Vinh N X，Chern J S，Lin C F. Phugiod Osillations in Optimal Reentry Trajectories［J］. Acta Astronautica.1981，8(4):311-324.

［4］Istratie V. Optimal Skip Entry with Terminal Maximum Velocity and Heat Constraint［C］.The 7th AIAA/ASME Joint Thermophysics and Heat Transfer Conference，Albuquerque，USA，June 15-18，1998.

［5］Gath P F. CAMTOS-A Software Suite Combining Direct and Indirect Trajectory Optimization Methods［D］. Stutgart:University of Stuttgart，2002.

［6］Clarke K A. Performance Optimization Study of a Common Aero Vehicle Using a Legendre Pseudo-spectral Method［D］. Massachusetts:Massachusetts Institute of Technology，2003.

［7］姚寅偉，李華濱.基于Gauss 偽譜法的高超聲速飛行器多約束三維再入軌跡優(yōu)化［J］. 航天控制，2012，30(2):33-38. (YAO Yinwei，LI Huabin. The Generation of Three-dimensional Constrained Entry Trajectories for Hypersonic Vehicle Based on the Gauss Pseudospectral Method［J］. Aerospace Control，2012，30(2):33-38.)

［8］洪蓓，辛萬(wàn)青.基于hp 自適應(yīng)偽譜法的固體運(yùn)載火箭軌 跡優(yōu)化［J］. 航天控制，2012，30(4):18-31.(HONG Bei，XIN Wanqin. Trajectory Optimization of Solid Lauch Vehicle Based on Hp-adaptive Pseudospectral Method［J］. Aerospace Control，2012，30(4):18-31.)

［9］閆循良，廖守億，張金生，等.基于節(jié)點(diǎn)改善策略的偽譜軌跡優(yōu)化［J］. 宇航學(xué)報(bào)，2013，34(7):891-900.(Yan X L，Liao S Y，Zhang J S，et al. Trajectory Optimization Using Pseudospectral Method Based on A Grid Node Refinement Strategy［J］. Journal of Astronautics，2013，34(7):891-900.)

［10］Lu P. Inverse Dynamics Approach to Trajectory Optimization for an Aerospace Plane［J］.Journal of Guidance，Control and Dynamics，1993，16(4):726-732.

［11］施法中.計(jì)算機(jī)輔助幾何設(shè)計(jì)與非均勻有理B 樣條［M］.北京:高等教育出版社，2001:55-62. (Shi F Z.CAGD＆NURBS［M］. Beijing:Higher Education Press，2001.)

［12］Indig N，Ben-Asher J Z，F(xiàn)arber N. Near-optimal Spatial Mid-course Guidance Law with Angular Constraint［C］.AIAA Guidance，Navigation and Control Conference and Exhibit，Minneapolis:2012. AIAA-2012-4691.

［13］Rao A V，Benson D A. Algorithm 902:GPOPS，A Matlab Software for Solving Multiple-phase Optimal Control Problems Using Gauss Pseudo-spectral Method［EB/OL］.2010［2013］. http://vdol. mae. ufl. edu/Journal-Publications/tomsgpm.pdf.