導彈彈道有理Bezier曲線造型與優化方法

張 松 侯明善 唐成師

1.西北工業大學自動化學院，西安710072

2.上海機電工程研究所，上海201109

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

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

1 彈道優化問題描述

1.1 導彈縱向運動描述

導彈力平衡及本文所采用的符號如圖1 所示。若不考慮地球自轉，認為重力場是平行力場，且導彈取為質點模型，其運動限制在鉛垂平面內，那么導彈運動方程組可寫為:

其中，

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

圖1 導彈力平衡與符號

1.2 彈道優化問題

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

末端時間tf(給定或未定)時狀態滿足約束:

狀態變量和控制變量約束條件為:

設彈道優化指標為

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

2 彈道Bezier 曲線造型與優化

2.1 有理Bezier 型彈道

有理n 次Bezier 曲線可表示為

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

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

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

設彈道邊界條件為:

則控制多邊形頂點b0，b1，bn－1以及bn應滿足下列條件:

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

則由于xs和xf為已知量，自由參數簡化為

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

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

2.2 基于彈道造型的優化問題描述

經彈道造型后，彈道優化問題就歸結為:確定自由參數Φ，使得與之對應的彈道在某項性能指標下最優。這就需要根據自由參數和導彈運動方程組確定性能指標。在本文方法中，有理Bezier 型彈道的位置變量x，h 及彈道傾角變量γ 完全由自由參數Φ確定，還需確定速度變量V，攻角α 以及時間變量t以計算性能指標。

記

則彈道傾角γ 可表示為

上式兩邊對參數ω 求導，可得

引入變換關系

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

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

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

因此，確定性能指標所需的各變量，彈道優化問題(9)～(12)就歸結為:確定自由參數Φ，使得性能指標

最小，并滿足狀態和控制約束

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

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

本文通過引入罰函數將多約束軌跡優化問題轉換為無約束軌跡優化問題，并保證在給定采樣節點處狀態和控制滿足約束(27)，具體步驟如下:

首先，令

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

其次，根據簡單約束(28)構造罰函數

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

從而約束優化問題(25)和(26)最終簡化為:確定自由參數Φ，使得指標函數

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

由有理Bezier 彈道的特性可知，優化問題(30)中優化變量數量少，所得彈道完全光滑，而且因狀態和控制變量都是通過等價轉換和積分來確定，故各變量精確滿足導彈運動方程組，即解的精度高。

3 彈道優化算例

本節以某型空空導彈為例，研究其最短時間飛行彈道。對于空空導彈而言，飛行時間短有利于快速打擊敵方目標，減少載機受攻擊的機會，擴大載機執行其它任務的靈活性，對于采用慣導形式中制導的空空導彈，更能有效減小誤差。因而要提升空空導彈的打擊能力必須把短的飛行時間作為一個基本的性能指標。本節取優化指標為

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

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

表1 邊界條件

仿真所用導彈的氣動參數來源于文獻［12］，其中CLα=7，CD0=0.4，K =0.01，S =0.02m2。空氣密度采用標準大氣模型計算，發動機推力T 是時間的函數，關系式為:

導彈質量變化規律˙m 滿足:

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

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

表2 自由參數優化搜索結果

表3 自由參數優化搜索結果(續)

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

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

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

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

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

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

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

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

4 結論

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

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

［2］雍恩米，陳磊，唐國金. 飛行器軌跡優化數值方法綜述［J］. 宇航學報，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 偽譜法的高超聲速飛行器多約束三維再入軌跡優化［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］洪蓓，辛萬青.基于hp 自適應偽譜法的固體運載火箭軌 跡優化［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］閆循良，廖守億，張金生，等.基于節點改善策略的偽譜軌跡優化［J］. 宇航學報，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］施法中.計算機輔助幾何設計與非均勻有理B 樣條［M］.北京:高等教育出版社，2001:55-62. (Shi F Z.CAGD＆NURBS［M］. Beijing:Higher Education Press，2001.)

［12］Indig N，Ben-Asher J Z，Farber 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.