飞行器轨迹优化的造型降维方法

张松，侯明善（.中航工业第一飞机设计研究院，陕西西安70089；.西北工业大学自动化学院，陕西西安7007）



飞行器轨迹优化的造型降维方法

张松1，侯明善2
（1.中航工业第一飞机设计研究院，陕西西安710089；2.西北工业大学自动化学院，陕西西安710072）

摘要：提出一种基于Bezier曲线造型的降维轨迹优化方法，用于提高轨迹优化的精度。该方法利用Bezier曲线优良的形状刻画能力描述最优轨迹，将优化问题的边界条件化为Bezier曲线的造型参数约束，从而使原问题表示为较低维含造型参数的优化问题，减少了计算量。轨迹造型法可有效协调飞行器动力学模型中不同时间尺度变量动态特性差异对非线性规划求解条件数的不利影响，更易于得到光滑最优解。按照常规高斯伪普法和网格自适应伪谱法分别对一种轨迹优化的原问题和造型降维问题进行了仿真验证，证明了所提出方法的可靠性。

关键词：飞行器控制、导航技术；轨迹优化；轨迹造型；Bezier型轨迹；伪谱法

0　引言

轨迹优化指确定一条满足给定动力学方程和相关约束条件的空间运动轨迹，使飞行器沿该轨迹飞行时指定的性能指标最优。求解轨迹优化问题的数值方法可分为间接法和直接法［1-2］。间接法［3-5］根据最优控制原理将优化问题转换为Hamilton边值问题进行求解，解的最优性能够得到保障，但难以适应复杂的路径约束情况，对初始猜测要求苛刻。直接法［6-9］根据一定准则离散化状态、控制，从而将连续最优控制问题转化为有限维非线性规划问题，进行迭代搜索求解。相对间接法，直接法收敛域更宽，鲁棒性更强。

近年来，直接法中的伪谱法受到越来越多的关注［2，7-9］。该方法采用全局插值多项式的有限基在一系列离散点上逼近状态和控制变量，并使近似状态在给定配点处满足动力学方程约束。配点通常选为正交多项式的根，因而此类方法能以较少的离散点获得较高精度的解，且对光滑问题是指数收敛的。然而，弹道优化问题具有高度非线性、多约束和多尺度等特性，通常是非光滑的。对于非光滑问题，伪谱法需要大量配点才满足精度要求，而随着配点的增加，对应非线性规划问题中约束雅可比矩阵和海塞矩阵的密度和规模将急剧增大，算法收敛性和精度将变差［10］。而且用同一离散化格式处理动力学模型中不同时间尺度的变量也不尽合理。此外，由于伪谱法中正交配点主要集中于区间两端，配点的增加将导致区间两端配点过度集中，进一步恶化非线性规划问题的条件数。为此，文献［11-12］等提出了不同的节点细化方案，将求解区间分段，并在每个区间使用低阶全局多项式逼近状态，以提高伪谱法解决非光滑问题的性能，但上述问题依然存在。提高算法鲁棒性和求解精度的一种有效方法是降低轨迹优化问题的维度，减少优化变量数。文献［13］利用非线性控制理论工具，将原始优化问题映射到低维空间进行求解，但是映射中的坐标变换并不容易实现。文献［14］根据逆动态方法，通过定义一组期望输出将控制变量消除，以达到降维的目的。但是该方法在优化过程中需要求解一组微分代数方程以确定控制变量。对于复杂的轨迹优化问题，这一微分代数方程并不容易求解，而且由此产生的计算误差会被引入到优化问题中。

从几何本质上看，飞行器飞行轨迹总是一条满足边界条件的光滑曲线。我们知道，对给定光滑曲线，曲线的切矢量和法矢量可通过方程计算确定。等价地说，如果轨迹方程已知，则轨迹的航迹角、航迹角速度就是确定的。研究发现，几乎所有的光滑飞行轨迹均可用Bezier曲线构造。受此启发，本文提出了一种能够有效改善直接轨迹优化方法性能的优化策略，根据边界条件将待求轨迹表示成光滑且只含少量造型变量的Bezier曲线形式，并结合动力学方程利用参数变换方法，将原始轨迹优化问题转换为低维最优控制问题，利用直接法进行求解。本文使用的直接法为伪谱法。由于能够有效改善直接法中非线性规划问题的条件数，本文方法解算轨迹优化问题的鲁棒性更强，精度更高。

1　轨迹优化问题描述

1.1动力学方程

惯性坐标系飞行器在铅垂平面内的质点动力学模型为

式中：x表示飞行器水平飞行距离；h表示飞行器飞行高度；v表示飞行器飞行速度；γ表示轨迹倾角；T表示飞行器发动机推力；α表示攻角；m表示飞行器质量；g表示重力加速度；Isp表示发动机比冲；D表示阻力；L表示升力；ρ表示大气密度；S表示参考面积；CL表示升力系数；CLα表示升力的攻角气动导数；α0是平衡点攻角；CD为阻力系数；CD0零升阻力系数；K为升阻极曲线系数。

1.2轨迹优化问题

取状态向量z=［x，h，v，γ］，控制变量uc=［α，T］，其容许集为U.设初始状态为

末端时间tf（给定或未定）时状态满足约束：

状态变量和控制变量约束条件为

设轨迹优化指标为

式中：Ψ（·）为Mayer型性能指标函数；L［·］为Lagrange型性能指标函数。

轨迹优化问题可表述为：在给定边界条件（8）式和（9）式情况下，确定满足约束条件（10）式和动力学方程（1）式～（5）式的控制变量uc（t）∈U使性能指标（11）式最优。下文中若不特别说明，则（·）s和（·）f分别表示变量（·）的初始值和末端值。

2　轨迹Bezier曲线造型与优化

2.1Bezier曲线

Bezier曲线是利用一组控制多边形定义的曲线，可通过合理选取控制多边形对曲线形状进行调整［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曲线具有诸多利于轨迹造型的优良性质，例如：Bezier曲线是n次连续的，首末端点分别是控制多边形的首末顶点，在首末端点处分别与控制多边形首末边相切，而且由同一组控制多边形定义的Bezier曲线唯一等。

2.2Bezier型轨迹

由前面的讨论可知，Bezier曲线是关于控制多边形顶点bi=［xi，hi］T，i=0，1，…，n的函数。若称xi和hi为造型变量，那么通过改变造型变量的值可以得到不同形状的Bezier曲线。使用Bezier曲线进行轨迹造型时，首先需要确定待造型轨迹的边界条件，即轨迹的起点、终点与相应的轨迹倾角，一般写为

为了满足轨迹的边界条件（16）式，控制多边形顶点b0、b1、bn-1以及bn应满足

用于调节Bezier曲线形状的造型变量成为x1，…，xn-1，h2，…，hn-2.通常飞行器飞行高度可以表示为射程的函数，因而选择造型变量x2，…，xn-2可取为区间［xs，xf］上具有特定分布特征的点，比如等距分布的点那么，由于xs和xf为已知量，造型变量简化为x1，xn-1和h2，…，hn-2.令Φ=［x1，xn-1，h2，…，hn-2］T，则Bezier曲线可根据（13）式表示为关于Φ的函数，即不妨称P（Φ，ω）为Bezier型轨迹。

需要强调的是，造型变量的多少完全由边界条件和Bezier曲线的次数确定。若不考虑末端轨迹倾角约束（18）式，选择造型变量 x2，…，xn-1为区间［xs，xf］上具有特定分布特征的点，则x1和hn-1也成为造型变量。为简化计算，选择造型变量 x2，…，xn-1为区间［xs，xf］上已知的等距分布点：（i-1），这时最终的造型变量为Φ=［x1，h2，…，hn-1］T.

2.3优化问题降维

现结合飞行器动力学方程与Bezier轨迹对原轨迹优化问题进行降维处理。

因而轨迹倾角γ可表示为

（21）式两边对参数ω求导，得到轨迹倾角关于变量ω的导数为

引入变换关系

对（1）式作变换得

因而有

对（3）式和（5）式作（23）式变换，并将（24）式代入可得速度和质量关于参数 ω∈［0，1］的微分方程组：

同理，对（4）式作变换，得到

因此，轨迹优化问题（8）式～（11）式就转换为如下关于造型变量Φ的降维优化问题：

可见，需要离散的变量由原来的［x，h，γ，v，m，α，T］变为［v，m，α，T］.

由于采用光滑的Bezier曲线对轨迹进行表示降低了优化问题的维数，利用直接法求解降维优化问题（27）式将比求解原优化问题鲁棒性更强，所得最优解精度也更高。

3　仿真研究

本节将仿真研究某型远程空空导弹的轨迹优化问题。远程空空导弹通常采用中制导与末制导相结合的复合制导模式，其中，中制导段的作用是将导弹由起始点引导至某一预定拦截点，并以有利态势进入末制导段。中制导段是整个制导过程中持续时间最长的阶段，而且该阶段结束时的速度决定了导弹在末制导阶段的能量状态和机动能力。因此，有必要以最大末速为指标研究导弹的中制导轨迹优化问题。

仿真使用的导弹模型参数取值［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］分别求解原轨迹优化问题（8）式～（11）式（下称原问题）和本文提出的造型降维优化问题（27）式（下称降维问题）。伪谱法是一种将连续轨迹优化问题转换为非线性规划问题求解的数值优化方法。本文采用序列二次规划算法（SQP）求解非线性规划问题。

本文采用的两种伪谱法分别是常规高斯伪普法［17］和网格自适应伪谱法［18］。为方便，称原问题高斯伪普法求解法为FDTO1，原问题网格自适应伪谱法求解为FDTO2；称降维问题高斯伪普法求解法为RDTO1，降维问题网格自适应伪谱法求解为RDTO2.

优化结果如图1～图6所示。图1～图6中FDTO1和FDTO2分别表示采用常规高斯伪普法［15］取50个配点以及网格自适应伪谱法［16］迭代3次对原问题的计算结果，而RDTO1和RDTO2则为相同条件下对本文降维问题的计算结果。图1表明，求解原问题与降维问题所得的最大末速轨迹基本一致。由图5和图6中的推力曲线可见，伪谱法在求解原问题时推力（虚线）存在剧烈的震荡，结果不可靠，而且引入网格自适应细化后，随着节点数的增加，震荡变得更加严重。类似震荡在速度（见图2）和攻角（见图4）曲线中同样存在。相同方法求解本文降维问题所得的推力（实线）及相应状态变量则更为平滑。这表明，与求解原轨迹优化问题相比，直接法求解造型降维后优化问题所形成非线性规划问题的规模大幅降低，因而非线性规划求解器的计算精度能够得到保证，所得最优解也更为可靠。

图1　最优轨迹比较Fig.1 Comparison of optimal trajectories

图2　速度-时间历程Fig.2 Velocity vs.time

下面对本文方法的收敛性和速度进行对比分析。本文用于求解非线性规划问题的SQP算法包含主迭代和子迭代两层结构：主迭代生成满足线性约束的迭代序列｛yk｝，且该序列最终收敛到满足一阶最优性xf条件的最优解；子迭代则通过解一系列二次规划（QP）子问题来产生主迭代中从点yk到点yk+1的搜索方向，其自身也存在迭代计算。因此，可用SQP算法的主迭代次数和子迭代总次数来表征本文方法的收敛性和速度。

图3　质量-时间历程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个配点条件下，利用常规高斯伪普法分别求解本文方法生成的降维问题与原问题时，SQP算法的迭代次数，其中，降维问题中造型曲线阶数分别取10、12、14和16.由表1可知，在优化指标值基本相同的情况下，与解算原问题相比，常规高斯伪普法解算降维问题所需的主迭代次数和子迭代总次数均较少。这表明对原问题进行降维处理后，求解的收敛性和速度显著提高。表1同时表明，增加造型曲线的阶数有利于提高解的最优性，但SQP算法将通过增加主迭代次数，来达到给定精度。

表1　迭代次数比较Tab.1 Comparison of iterations

4　结论

轨迹优化问题因具有高度非线性、多约束和多尺度等特性，通常是非光滑的，采用常规直接法求解存在最优解不光滑现象。为此，本文提出了基于Bezier曲线轨迹造型的降维轨迹优化策略，根据边界条件将待求轨迹表示成光滑且只含少量造型变量的Bezier曲线形式，并结合导弹动力学方程，利用参数变换方法将原轨迹优化问题转换为低维最优控制问题求解。仿真实验表明，由于条件数得到改善，直接法求解此类问题的精度更高，所得最优解更为光滑。

参考文献（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］ 雍恩米，陈磊，唐国金.飞行器轨迹优化数值方法综述［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伪谱法的制导炸弹最优滑翔弹道研究［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伪谱最优制导律［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