飛行器軌跡優(yōu)化的造型降維方法

張松，侯明善（.中航工業(yè)第一飛機設計研究院，陜西西安70089；.西北工業(yè)大學自動化學院，陜西西安7007）

?

飛行器軌跡優(yōu)化的造型降維方法

張松1，侯明善2
（1.中航工業(yè)第一飛機設計研究院，陜西西安710089；2.西北工業(yè)大學自動化學院，陜西西安710072）

摘要：提出一種基于Bezier曲線造型的降維軌跡優(yōu)化方法，用于提高軌跡優(yōu)化的精度。該方法利用Bezier曲線優(yōu)良的形狀刻畫能力描述最優(yōu)軌跡，將優(yōu)化問題的邊界條件化為Bezier曲線的造型參數(shù)約束，從而使原問題表示為較低維含造型參數(shù)的優(yōu)化問題，減少了計算量。軌跡造型法可有效協(xié)調(diào)飛行器動力學模型中不同時間尺度變量動態(tài)特性差異對非線性規(guī)劃求解條件數(shù)的不利影響，更易于得到光滑最優(yōu)解。按照常規(guī)高斯偽普法和網(wǎng)格自適應偽譜法分別對一種軌跡優(yōu)化的原問題和造型降維問題進行了仿真驗證，證明了所提出方法的可靠性。

關鍵詞：飛行器控制、導航技術；軌跡優(yōu)化；軌跡造型；Bezier型軌跡；偽譜法

0　引言

軌跡優(yōu)化指確定一條滿足給定動力學方程和相關約束條件的空間運動軌跡，使飛行器沿該軌跡飛行時指定的性能指標最優(yōu)。求解軌跡優(yōu)化問題的數(shù)值方法可分為間接法和直接法［1-2］。間接法［3-5］根據(jù)最優(yōu)控制原理將優(yōu)化問題轉(zhuǎn)換為Hamilton邊值問題進行求解，解的最優(yōu)性能夠得到保障，但難以適應復雜的路徑約束情況，對初始猜測要求苛刻。直接法［6-9］根據(jù)一定準則離散化狀態(tài)、控制，從而將連續(xù)最優(yōu)控制問題轉(zhuǎn)化為有限維非線性規(guī)劃問題，進行迭代搜索求解。相對間接法，直接法收斂域更寬，魯棒性更強。

近年來，直接法中的偽譜法受到越來越多的關注［2，7-9］。該方法采用全局插值多項式的有限基在一系列離散點上逼近狀態(tài)和控制變量，并使近似狀態(tài)在給定配點處滿足動力學方程約束。配點通常選為正交多項式的根，因而此類方法能以較少的離散點獲得較高精度的解，且對光滑問題是指數(shù)收斂的。然而，彈道優(yōu)化問題具有高度非線性、多約束和多尺度等特性，通常是非光滑的。對于非光滑問題，偽譜法需要大量配點才滿足精度要求，而隨著配點的增加，對應非線性規(guī)劃問題中約束雅可比矩陣和海塞矩陣的密度和規(guī)模將急劇增大，算法收斂性和精度將變差［10］。而且用同一離散化格式處理動力學模型中不同時間尺度的變量也不盡合理。此外，由于偽譜法中正交配點主要集中于區(qū)間兩端，配點的增加將導致區(qū)間兩端配點過度集中，進一步惡化非線性規(guī)劃問題的條件數(shù)。為此，文獻［11-12］等提出了不同的節(jié)點細化方案，將求解區(qū)間分段，并在每個區(qū)間使用低階全局多項式逼近狀態(tài)，以提高偽譜法解決非光滑問題的性能，但上述問題依然存在。提高算法魯棒性和求解精度的一種有效方法是降低軌跡優(yōu)化問題的維度，減少優(yōu)化變量數(shù)。文獻［13］利用非線性控制理論工具，將原始優(yōu)化問題映射到低維空間進行求解，但是映射中的坐標變換并不容易實現(xiàn)。文獻［14］根據(jù)逆動態(tài)方法，通過定義一組期望輸出將控制變量消除，以達到降維的目的。但是該方法在優(yōu)化過程中需要求解一組微分代數(shù)方程以確定控制變量。對于復雜的軌跡優(yōu)化問題，這一微分代數(shù)方程并不容易求解，而且由此產(chǎn)生的計算誤差會被引入到優(yōu)化問題中。

從幾何本質(zhì)上看，飛行器飛行軌跡總是一條滿足邊界條件的光滑曲線。我們知道，對給定光滑曲線，曲線的切矢量和法矢量可通過方程計算確定。等價地說，如果軌跡方程已知，則軌跡的航跡角、航跡角速度就是確定的。研究發(fā)現(xiàn)，幾乎所有的光滑飛行軌跡均可用Bezier曲線構造。受此啟發(fā)，本文提出了一種能夠有效改善直接軌跡優(yōu)化方法性能的優(yōu)化策略，根據(jù)邊界條件將待求軌跡表示成光滑且只含少量造型變量的Bezier曲線形式，并結合動力學方程利用參數(shù)變換方法，將原始軌跡優(yōu)化問題轉(zhuǎn)換為低維最優(yōu)控制問題，利用直接法進行求解。本文使用的直接法為偽譜法。由于能夠有效改善直接法中非線性規(guī)劃問題的條件數(shù)，本文方法解算軌跡優(yōu)化問題的魯棒性更強，精度更高。

1　軌跡優(yōu)化問題描述

1.1動力學方程

慣性坐標系飛行器在鉛垂平面內(nèi)的質(zhì)點動力學模型為

式中：x表示飛行器水平飛行距離；h表示飛行器飛行高度；v表示飛行器飛行速度；γ表示軌跡傾角；T表示飛行器發(fā)動機推力；α表示攻角；m表示飛行器質(zhì)量；g表示重力加速度；Isp表示發(fā)動機比沖；D表示阻力；L表示升力；ρ表示大氣密度；S表示參考面積；CL表示升力系數(shù)；CLα表示升力的攻角氣動導數(shù)；α0是平衡點攻角；CD為阻力系數(shù)；CD0零升阻力系數(shù)；K為升阻極曲線系數(shù)。

1.2軌跡優(yōu)化問題

取狀態(tài)向量z=［x，h，v，γ］，控制變量uc=［α，T］，其容許集為U.設初始狀態(tài)為

末端時間tf（給定或未定）時狀態(tài)滿足約束：

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

設軌跡優(yōu)化指標為

式中：Ψ（·）為Mayer型性能指標函數(shù)；L［·］為Lagrange型性能指標函數(shù)。

軌跡優(yōu)化問題可表述為：在給定邊界條件（8）式和（9）式情況下，確定滿足約束條件（10）式和動力學方程（1）式～（5）式的控制變量uc（t）∈U使性能指標（11）式最優(yōu)。下文中若不特別說明，則（·）s和（·）f分別表示變量（·）的初始值和末端值。

2　軌跡Bezier曲線造型與優(yōu)化

2.1Bezier曲線

Bezier曲線是利用一組控制多邊形定義的曲線，可通過合理選取控制多邊形對曲線形狀進行調(diào)整［15］。若給定控制多邊形的頂點bi=［xi，hi］T，i= 0，1，…，n，則n次Bezier曲線表達式為

式中：Bn，i（ω）為伯恩斯坦多項式，

使用（13）式德卡斯特里奧遞推算法求n次Bezier曲線：

n次Bezier曲線的導矢可寫為

式中：Δbi為向前的差分矢量，Δbi=bi+1-bi.對（14）式重復求導，可以得到Bezier曲線的高階導矢公式為

式中：高階向前差分矢量由低階向前差分矢量遞推地定義為Δkbi=Δk-1bi+1-Δk-1bi.

Bezier曲線具有諸多利于軌跡造型的優(yōu)良性質(zhì)，例如：Bezier曲線是n次連續(xù)的，首末端點分別是控制多邊形的首末頂點，在首末端點處分別與控制多邊形首末邊相切，而且由同一組控制多邊形定義的Bezier曲線唯一等。

2.2Bezier型軌跡

由前面的討論可知，Bezier曲線是關于控制多邊形頂點bi=［xi，hi］T，i=0，1，…，n的函數(shù)。若稱xi和hi為造型變量，那么通過改變造型變量的值可以得到不同形狀的Bezier曲線。使用Bezier曲線進行軌跡造型時，首先需要確定待造型軌跡的邊界條件，即軌跡的起點、終點與相應的軌跡傾角，一般寫為

為了滿足軌跡的邊界條件（16）式，控制多邊形頂點b0、b1、bn-1以及bn應滿足

用于調(diào)節(jié)Bezier曲線形狀的造型變量成為x1，…，xn-1，h2，…，hn-2.通常飛行器飛行高度可以表示為射程的函數(shù)，因而選擇造型變量x2，…，xn-2可取為區(qū)間［xs，xf］上具有特定分布特征的點，比如等距分布的點那么，由于xs和xf為已知量，造型變量簡化為x1，xn-1和h2，…，hn-2.令Φ=［x1，xn-1，h2，…，hn-2］T，則Bezier曲線可根據(jù)（13）式表示為關于Φ的函數(shù)，即不妨稱P（Φ，ω）為Bezier型軌跡。

需要強調(diào)的是，造型變量的多少完全由邊界條件和Bezier曲線的次數(shù)確定。若不考慮末端軌跡傾角約束（18）式，選擇造型變量 x2，…，xn-1為區(qū)間［xs，xf］上具有特定分布特征的點，則x1和hn-1也成為造型變量。為簡化計算，選擇造型變量 x2，…，xn-1為區(qū)間［xs，xf］上已知的等距分布點：（i-1），這時最終的造型變量為Φ=［x1，h2，…，hn-1］T.

2.3優(yōu)化問題降維

現(xiàn)結合飛行器動力學方程與Bezier軌跡對原軌跡優(yōu)化問題進行降維處理。

因而軌跡傾角γ可表示為

（21）式兩邊對參數(shù)ω求導，得到軌跡傾角關于變量ω的導數(shù)為

引入變換關系

對（1）式作變換得

因而有

對（3）式和（5）式作（23）式變換，并將（24）式代入可得速度和質(zhì)量關于參數(shù) ω∈［0，1］的微分方程組：

同理，對（4）式作變換，得到

因此，軌跡優(yōu)化問題（8）式～（11）式就轉(zhuǎn)換為如下關于造型變量Φ的降維優(yōu)化問題：

可見，需要離散的變量由原來的［x，h，γ，v，m，α，T］變?yōu)椋踲，m，α，T］.

由于采用光滑的Bezier曲線對軌跡進行表示降低了優(yōu)化問題的維數(shù)，利用直接法求解降維優(yōu)化問題（27）式將比求解原優(yōu)化問題魯棒性更強，所得最優(yōu)解精度也更高。

3　仿真研究

本節(jié)將仿真研究某型遠程空空導彈的軌跡優(yōu)化問題。遠程空空導彈通常采用中制導與末制導相結合的復合制導模式，其中，中制導段的作用是將導彈由起始點引導至某一預定攔截點，并以有利態(tài)勢進入末制導段。中制導段是整個制導過程中持續(xù)時間最長的階段，而且該階段結束時的速度決定了導彈在末制導階段的能量狀態(tài)和機動能力。因此，有必要以最大末速為指標研究導彈的中制導軌跡優(yōu)化問題。

仿真使用的導彈模型參數(shù)取值［16］為：CLα= 35.0，CD0=0.92，K=0.03，S=0.03 m2，Isp=200 s，m0=240 kg，mf=132 kg，α0=0°.

本文采用兩種偽譜法［17-18］分別求解原軌跡優(yōu)化問題（8）式～（11）式（下稱原問題）和本文提出的造型降維優(yōu)化問題（27）式（下稱降維問題）。偽譜法是一種將連續(xù)軌跡優(yōu)化問題轉(zhuǎn)換為非線性規(guī)劃問題求解的數(shù)值優(yōu)化方法。本文采用序列二次規(guī)劃算法（SQP）求解非線性規(guī)劃問題。

本文采用的兩種偽譜法分別是常規(guī)高斯偽普法［17］和網(wǎng)格自適應偽譜法［18］。為方便，稱原問題高斯偽普法求解法為FDTO1，原問題網(wǎng)格自適應偽譜法求解為FDTO2；稱降維問題高斯偽普法求解法為RDTO1，降維問題網(wǎng)格自適應偽譜法求解為RDTO2.

優(yōu)化結果如圖1～圖6所示。圖1～圖6中FDTO1和FDTO2分別表示采用常規(guī)高斯偽普法［15］取50個配點以及網(wǎng)格自適應偽譜法［16］迭代3次對原問題的計算結果，而RDTO1和RDTO2則為相同條件下對本文降維問題的計算結果。圖1表明，求解原問題與降維問題所得的最大末速軌跡基本一致。由圖5和圖6中的推力曲線可見，偽譜法在求解原問題時推力（虛線）存在劇烈的震蕩，結果不可靠，而且引入網(wǎng)格自適應細化后，隨著節(jié)點數(shù)的增加，震蕩變得更加嚴重。類似震蕩在速度（見圖2）和攻角（見圖4）曲線中同樣存在。相同方法求解本文降維問題所得的推力（實線）及相應狀態(tài)變量則更為平滑。這表明，與求解原軌跡優(yōu)化問題相比，直接法求解造型降維后優(yōu)化問題所形成非線性規(guī)劃問題的規(guī)模大幅降低，因而非線性規(guī)劃求解器的計算精度能夠得到保證，所得最優(yōu)解也更為可靠。

圖1　最優(yōu)軌跡比較Fig.1 Comparison of optimal trajectories

圖2　速度-時間歷程Fig.2 Velocity vs.time

下面對本文方法的收斂性和速度進行對比分析。本文用于求解非線性規(guī)劃問題的SQP算法包含主迭代和子迭代兩層結構：主迭代生成滿足線性約束的迭代序列｛yk｝，且該序列最終收斂到滿足一階最優(yōu)性xf條件的最優(yōu)解；子迭代則通過解一系列二次規(guī)劃（QP）子問題來產(chǎn)生主迭代中從點yk到點yk+1的搜索方向，其自身也存在迭代計算。因此，可用SQP算法的主迭代次數(shù)和子迭代總次數(shù)來表征本文方法的收斂性和速度。

圖3　質(zhì)量-時間歷程Fig.3 Mass vs.time

圖4　攻角-時間歷程Fig.4 Angle of attack vs.time

圖5　推力-時間歷程1Fig.5 Thrust vs.time 1

圖6　推力-時間歷程2Fig.6 Thrust vs.time 2

表1顯示的是取50個配點條件下，利用常規(guī)高斯偽普法分別求解本文方法生成的降維問題與原問題時，SQP算法的迭代次數(shù)，其中，降維問題中造型曲線階數(shù)分別取10、12、14和16.由表1可知，在優(yōu)化指標值基本相同的情況下，與解算原問題相比，常規(guī)高斯偽普法解算降維問題所需的主迭代次數(shù)和子迭代總次數(shù)均較少。這表明對原問題進行降維處理后，求解的收斂性和速度顯著提高。表1同時表明，增加造型曲線的階數(shù)有利于提高解的最優(yōu)性，但SQP算法將通過增加主迭代次數(shù)，來達到給定精度。

表1　迭代次數(shù)比較Tab.1 Comparison of iterations

4　結論

軌跡優(yōu)化問題因具有高度非線性、多約束和多尺度等特性，通常是非光滑的，采用常規(guī)直接法求解存在最優(yōu)解不光滑現(xiàn)象。為此，本文提出了基于Bezier曲線軌跡造型的降維軌跡優(yōu)化策略，根據(jù)邊界條件將待求軌跡表示成光滑且只含少量造型變量的Bezier曲線形式，并結合導彈動力學方程，利用參數(shù)變換方法將原軌跡優(yōu)化問題轉(zhuǎn)換為低維最優(yōu)控制問題求解。仿真實驗表明，由于條件數(shù)得到改善，直接法求解此類問題的精度更高，所得最優(yōu)解更為光滑。

參考文獻（References）

［1］ Betts J T.Survey of numerical methods for trajectory optimization ［J］.Journal of Guidance，Control and Dynamic，1998，21（2）：193-206.

［2］ 雍恩米，陳磊，唐國金.飛行器軌跡優(yōu)化數(shù)值方法綜述［J］.宇航學報，2008，29（2）：7-16. YONG En-mi，CHEN Lei，TANG Guo-jin.A survey of numerical methods for trajectory optimization of spacecraft［J］.Journal of Astronautics，2008，29（2）：7-16.（in Chinese）

［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］∥7th AIAA/ASME Joint Thermophysics and Heat Transfer Conference.Albuquerque，US：AIAA，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］ Benson A，Thorvaldsen T，Rao V.Direct trajectory optimization and costate estimation via an orthogonal collocation method［J］. Journal of Guidance，Control and Dynamics，2006，29（6）：1435-1440.

［8］ 袁宴波，張科，薛曉東.基于Radau偽譜法的制導炸彈最優(yōu)滑翔彈道研究［J］.兵工學報，2014，35（8）：1179-1186. YUAN Yan-bo，ZHANG Ke，XUE Xiao-dong.Optimization of glide trajectory of guided bombs using a Radau pseudo-spectral method［J］.Acta Armamentarii，2014，35（8）：1179-1186.（in Chinese）

［9］ 陳琦，王中原，常思江.帶落角約束的間接Gauss偽譜最優(yōu)制導律［J］.兵工學報，2015，36（7）：1203-1212. CHEN Qi，WANG Zhong-yuan，CHANG Si-jiang.Optimal guidance law with impact angle constraints based on indirect Gauss pseudospectral method［J］.Acta Armamentarii，2015，36（7）：1203-1212.（in Chinese）

［10］ Berrut J P，Baltensperger R，Mittelman H D.Recent developments in barycentric rational interpolation［J］.International Series of Numerical Mathematics，2005，151：27-51.

［11］ Darby C L，Rao A V.A state approximation-based mesh refinement algorithm for solving optimal control problems［C］∥AIAA Guidance，Navigation，and Control Conference.Chicago，US：AIAA，2009：10-13.

［12］ Darby C L，Hager W W，Rao A V.An hp-adaptive pseudospectral method for solving optimal control problems［J］.Optimal Control Application and Methods，2011，32（4）：476-502.

［13］ Milan M B，Mushumbi K，Murray R M.A new computational approach to real-time trajectory generation for constrained mechanical systems［C］∥IEEE Conference on Decision and Control.Sydney，NSW：IEEE，2000：845-851.

［14］ Lu P.Inverse dynamics approach to trajectory optimization for an aerospace plane［J］.Journal of Guidance，Control and Dynamics，1993，16（4）：726-732.

［15］ 施法中.計算機輔助幾何設計與非均勻有理B樣條［M］.北京：高等教育出版社，2001：55-62. SHI Fa-zhong.CAGD&NURBS［M］.Beijing：Higher Education Press，2001：55-62.（in Chinese）

［16］ Fumiaki I，Takeshi K，Susumu M.Optimal thrust control of a missile with a pulse motor［J］.Journal of Guidance，Control，and Dynamics，1991，14（2）：377-382.

［17］ HuntigtonGT.AdvancementandanalysisofaGauss pseudospectral transcription for optimal control［D］.Cambridage：Massachusetts Institute of Technology，2007.

［18］ Rao A V，Benson D A.Algorithm 902：GPOPS，a MATLAB software for solving multiple-phase optimal control problems using Gauss pseudo-spectral method［J］.ACM Transactions on Mathematical Software，2010，37（2）：22：1-22：39.

中圖分類號：V19

文獻標志碼：A

文章編號：1000-1093（2016）06-1125-06

DOI：10.3969/j.issn.1000-1093.2016.06.022

收稿日期：2015-03-10

作者簡介：張松（1985—），男，博士研究生。E-mail：zhangsong.gz@outlook.com；侯明善（1959—），男，教授，博士生導師。E-mail：mingshan@nwpu.edu.cn

Trajectory Optimization of Aerocraft Based on Shaping and Dimension Reduction

ZHANG Song1，HOU Ming-shan2
（1.The First Aircraft Institute，Aviation Industry Corporation of China，Xi'an 710089，Shaanxi，China；2.School of Automation，Northwestern Polytechnical University，Xi'an 710072，Shaanxi，China）

Abstract：A dimension reduction optimization strategy based on trajectory shaping is presented to improve the accuracy of trajectory optimization.The proposed method is able to improve the performance of the direct trajectory optimization method in solving non-smooth trajectory optimization problem.According to the boundary conditions，the unknown trajectory is shaped using the Bezier curves with a few shape variables.The concepts of inverse dynamic and parametric transform are introduced to reduce the original problem to a lower dimensional optimal control problem.The result shows that lower dimensional problem tends to have a better conditioning number when it is solved by direct optimization method，and the accuracy of solutions is improved dramatically.The validation of the proposed method over conventional direct optimization method is demonstrated by the simulation results.

Key words：control and navigation technology of aerocraft；trajectory optimization；trajectory shaping；Bezier trajectory；pseudo-spectral method