考虑气动参数扰动的再入轨迹快速优化方法

杨 奔 雷建长 王宇航

(中国运载火箭技术研究院，北京 100071)

引言

再入飞行器利用气动力在大气层内进行长时间、远距离、高马赫数的飞行，期间会受到极其复杂的非线性热流、动压、过载约束以及终端的高精度约束，使其轨迹规划问题成为了一个极富挑战性的问题，引起了国内外许多学者的研究兴趣[1].传统的再入弹道优化设计方法主要有基于极大值原理的间接法和基于非线性规划理论的直接法[2]等优化算法，两者或因巨大的计算量，或因不易收敛，都只适合离线的轨迹优化[3-6].虽然航天飞机所采用的基于阻力加速度剖面的轨迹规划[7]，曾取得过辉煌的成绩，但该方法的简化过程限制了飞行器的机动能力.总的来说，基于极大值原理的直接法推导最优解的过程较为繁琐，对于复杂问题几乎无法得到解析解；以伪谱法为代表的直接法对多阶段等复杂问题的建模过程比较麻烦；而以粒子群法等现代启发式算法对于计算量大、优化参数多的复杂模型优化问题，需反复交叉、变异、迭代，计算效率较低.

强耦合性、强非线性决定了再入轨迹规划问题难以直接求解，通常需要采用一系列易于求解的子问题逼近原问题.凸优化具有多项式复杂度，局部最优解即全局最优解等特点，加上独特的内点法，使得数学上所有易解的规划问题几乎都指向了凸优化[8-9].美国加州理工大学的Acikmese 等[10-15]最早开展了基于凸优化的轨迹优化研究.Acikmese 等[10-11]针对火星软着陆问题最先提出了无损凸化(lossless convexification)的概念，使用松弛的凸形式约束代替原非凸约束，将原问题转化为凸问题，并利用极大值原理对二者的等价性进行了证明，保证转化过程的“无损性”.针对具有非凸约束的一般形式线性系统最优控制问题，Acikmese 等[12-13]讨论了无损凸化的方法，对原问题与凸优化问题的等价性进行了证明.相关理论研究成果已经在NASA 的G-FOLD 飞行试验中得到了成功应用[14].

近年来，许多学者都在尝试将凸优化引入到再入轨迹规划问题中来[15-20].Liu 等[15]在给定攻角剖面的条件下，通过变量替换得到以能量为自变量的动力学方程，采用序列凸优化解决了再入轨迹优化问题.Wang 等[16]则以时间为自变量并以倾侧角的导数为控制量，求解了终端时间固定的轨迹规划.现有的基于凸优化的滑翔再入轨迹规划方法基本上都直接使用了小扰动线性化来处理微分方程的非线性同时相应的添加了信赖域约束来保证子问题向原问题逼近.虽然目前的序列凸优化技术简单可行，寻优时间较传统直接法有一定优势，但是在实际仿真计算中还会遇到以下问题:(1)反复迭代的序列凸优化技术本质上还是牛顿法，其收敛范围小，并且在迭代后期可能会出现类似与牛顿法的Maratos 效应，导致迭代寻优失败; (2) 当初值给的比较粗糙时，子问题的可行性将是一个严峻的问题，可能出现“伪不可行”[17](artificial infeasibility)，的情况使得迭代失败.

针对上述序列凸优化的不足之处，本文提出了几种改善其收敛性的方法.首先通过增加虚拟控制量，可避免迭代初期由于控制能力的不足带来的“伪不可行”问题，降低该算法对初始参考轨迹的依赖.其次通过B 样条曲线离散控制量，抑制由于离散带来的锯齿化现象.最后，引入一种“回溯搜索”的方法，改善迭代过程中状态量的振动，从而加快算法的收敛.

在高超声速飞行器实际的再入飞行过程中，大气参数扰动以及外界随机干扰都会对再入轨迹产生一定的影响.现有的基于序列凸优化的轨迹优化方法，均属于确定性轨迹优化方法，即在优化的过程中未考虑参数扰动和初值不确定性对优化轨迹的影响.在上述方法设计的最优输入的控制下，高超声速飞行器一旦受到不确定扰动或干扰的影响，可能会偏离预定优化轨迹，难以满足预期设定的过程约束、终端约束以及禁飞区等地理约束，甚至可能出现飞行器烧蚀坠落的事故.因此，在轨迹优化设计的过程中，除了确定性的动力学建模分析之外，还需要综合考虑参数扰动的影响，完成复杂不确定条件下的高超声速再入轨迹的优化设计，提高再入飞行轨迹的鲁棒性.

目前鲁棒轨迹优化方法在高超声速领域应用较少.Rangaraj 等[21]利用多项式混沌方法设计了一种改进粒子滤波器，与传统的粒子滤波器相比，计算效率大大提高.Xiong 等[22]针对一类具有不确定性的轨迹优化问题，将不确定的路径约束使用多项式混沌方法表示为均值和方差的约束，提高了设计效率.Li 等[23]设计了一种将不确定约束转化为确定性优化的非侵入的多项式混沌飞行轨迹优化方法.Michael 等[24]对直接法进行改进，提出了一种快速鲁棒轨迹优化方法.此外结合高斯伪谱法和多项式混沌理论的鲁棒轨迹优化方法也成为研究的热点[24-26].

本文为研究飞行器再入过程中的气动参数扰动问题，提出了采样点少、易于实现，计算效率高的广义混沌多项式理论研究方法，建立了基于广义混沌多项式和凸优化相结合的再入轨迹鲁棒优化模型，该模型可有效降低优化轨迹对气动参数扰动的敏感程度，在不确定条件的干扰下，展现出了较强的鲁棒稳定性.

1 确定性轨迹优化模型

1.1 再入动力学方程

考虑地球自转的影响，建立如下再入飞行器无量纲运动方程

其中，r,θ,φ,V,γ,ψ 分别为无量纲化的地心距、经度、纬度、对地速度、航迹角、航向角，地心距按地球平均半径归一化，对地速度按第一宇宙速度归一化.速度σ 为倾侧角，Ωg为地球自转角速度，其对归一化，无量纲化升力L、阻力D的表达式为

其中，ρ 是大气密度，为简化问题，采用指数形式的大气模型ρ=ρ0exp(−βh)，S为飞行器的气动参考面积，升力系数CL、阻力系数CD均认为是攻角和马赫数的函数.

1.2 基本约束

飞行器再入过程中，为保证飞行任务的成功，需满足热流、动压、过载等路径约束，可表示为

考虑伺服结构的约束，保证优化结果光滑可行，选择倾侧角的变化率为控制量，式(1)中需增加一个额外的状态方程=u，则状态变量为x=，该策略可避免对控制量进行凸化处理.

至此，得到该最优控制问题(P1)如下

式(7)中动力学方程可表示为

其中，第一项是与状态变量相关的项，第二项是与控制量相关的项，第三项是与地球自旋相关项.可以看出，选择倾侧角变化率为控制量，使得控制量与状态量自然解耦，更加有利于优化收敛.

1.2 凸化处理

由1.2 节可知，P1是一个高度的非线性最优控制问题，无法直接采用凸优化方法求解.目前针对这种复杂的多约束最优控制问题基本上均是通过线性化处理来求解.因此本文采用基于信赖域约束的序列凸优化方法.通过求解一系列的凸优化子问题，来逐渐逼近原有问题的解.

将式(10)一阶泰勒展开，由于地球自旋速度较慢，可认为fΩ(x)≈fΩ(xk)，则线性化的动力学方程为

式(3)∼式(5)中的基本约束也可通过上述凸化方法处理，序列凸优化的迭代寻优过程是在给定的参考轨迹附近线性化，为保证子问题向原问题逼近，需增加信赖域约束.

为避免迭代初期“伪不可行”的问题，降低对初值的依赖程度，同时弥补线性化带来的误差，在式(11)中增加一个虚拟控制量w，可得

可在迭代过程中，在目标函数中控制虚拟控制量的范围，使其在迭代后期为小量，从而保证问题的准确性,即

2 算法表述

2.1 问题离散

为求解无限维最优控制问题的数值解，必须对其进行离散.为了简单起见，状态变量像通常一样离散成等距的数值点，各个数值点处的控制量则通过B 样条曲线控制点的表达式表示，控制点个数的选取可在实际仿真中调整，优化过程中的积分项均采用梯形积分表示.仿真结果显示这种方法可以十分有效地抑制数值求解中控制量的锯齿化现象，更加有利于迭代收敛，并且减少了需要求解的变量个数.将时间N等份，可得离散后的动力学方程为

2.2 回溯搜索

大气层内各个状态变量与控制量相互耦合，具有极强的非线性，且各个约束集也均为非凸，导致无法直接运用内点法进行求解.受序列二次规划算法的的启发，国内外学者提出了基于信赖域约束的序列凸优化技术.其信赖域约束也就是为了保障了线性化的可靠程度，当然每次外层迭代所使用的参考轨迹的不同使得每次迭代过程并不是同一个问题，也就是说并不是严格的最优解.若取固定信赖域，且每次迭代间均采用全步长更新寻优结果，在迭代后期可能会出现如图1 和图2 的情况，各个状态变量在迭代后期产生剧烈振荡，难以符合停机准则，导致寻优失败.

图1 振荡现象Fig.1 Vibration phenomena

针对这种具有两个局部极小值，得益于牛顿法的启发，本文提出一种“回溯搜索”的收敛办法.通过设置阻尼步长和回归搜索来快速处理这种局部极小值问题.具体过程如下：

(1) 初值选取，并线性化，取两个回溯参数α ∈(0,1],β ∈(0,1] 以及信赖域∆0

(2)解决子SOCP 问题，得到∆x.

(3)判断停机准则.

(4)判断回溯条件：I f J(k+1)>J(k)，x(k+1)=x(k)+α∆x，∆0=β∆0，进入步骤一.

该回溯搜索的关键点在于信赖域的收缩权值以及阻尼步长权值的选取，调试初期可先设为黄金分割数，后期在具体的仿真过程中适当调整.一方面信赖域的收缩速度不宜过快，否则可能及早地陷入局部最优点；另外一方面阻尼步长是在两次轨迹迭代轨迹之间的权衡，该值的选择决定了迭代收敛速度的快慢，回溯系数实质上代表了迭代搜索的方向，可决定收敛速度的快慢及计算结果的最优性.

3 考虑气动参数扰动的随机轨迹优化方法

3.1 广义混沌多项式理论

混沌多项式理论的基本思想是利用一系列特定分布的随机变量的正交多项式对随机过程与扰动进行逼近[27-29].多项式混沌方法能够用于多种随机模型的求解过程，如集成电路的响应分析[30]、不同类型流体问题研究[31-32]以及其他问题的研究中.

这里采用符号ω 表示随机事件，Ω 为随机事件的集合.定义随机变量矢量为

其中每个元素δi(ω)均为随机变量.同时在空间定义其内积为

根据Cameron-Martin 理论，这些随机变量δi(ω)的多项式运算在L2(Ω)空间中是稠密的，也就是说任何随机过程都可以由L2(Ω)空间中随机变量的正交多项式逼近.基于上述结论，任意随机过程X(ω)都可以由混沌多项式近似为

采用不高于p阶的多项式进行有限项截断，并写成紧凑形式为

其中，Pj(∆(ω))为随机变量的第j阶基底多项式，xj为混沌多项式系数.在随后的分析中将随机变量∆(ω)简写为∆，混沌多项式Ψ(∆(ω)) 简写为Ψ.NPC是多项式展开式的项数，可由下式确定

对于分布函数或者概率密度已知的随机变量，式中的系数为确定值.多项式的选择依赖于随机变量的概率密度函数，根据Askey 法则，对于不同概率密度函数，存在不同的最优多项式，具体对应形式如表1 所示[28].

表1 正交多项式与变量分布类型对应表Table 1 Table of orthogonal polynomials with distribution types

根据正交多项式法则，式(20)中的多项式展开系数可以表示为

等式右侧的计算，可以使用配点法求解，将连续积分问题转化为离散点位置函数值的加权和.配点和权重需要依据相关变量的分布类型进行选择.对于均匀分布的随机变量，需要选择Legendre-Gauss 点和Gauss 权重.则相应的积分计算为

其中，Q为配点的总个数，相应的正交权重由τm表示.

应用混沌多项式可方便得到最优轨迹对应的随机变量的统计信息，这些信息可通过混沌多项式的系数得到，其均值可表示为

方差为

3.2 基于广义混沌多项式的鲁棒轨迹优化

在飞行器再入飞行阶段，存在气动参数扰动等不确定性的作用，若仍然利用基于确定性系统的优化所得到的倾侧角控制律，热流、动压等路径约束将存在违背的情况，从而可能导致发射任务的失败.因此考虑再入段的气动参数扰动不确定性扰动，本节提出再入段鲁棒轨迹优化方法，通过该方法可降低最优轨迹对扰动的敏感程度，确保优化结果满足所有规定约束.

考虑到非线性最优控制问题的复杂性，这里主要考虑采用状态的均值和标准差近似描述其分布特性.将传统轨迹优化问题转化为鲁棒轨迹优化问题,如下

为提高计算效率，本文采用回归法[30]求解式(20)中的广义混沌多项式展开系数，首先通过配点法对随机变量δ 进行采样，即δ={δ(1),δ(2),···,δ(NS)}.将优化过程中的7 个状态变量x=[r;θ;ϕ;V;γ;ψ;σ]通过混沌多项式展开，并将离散点代入式(20)可得

将式(28)代入动力学方程可表示为以下矩阵形式为

可将式(29)简写为=Fi，由于采样点的数目大于多项式混沌展开系数的个数，可通过最小二乘法求解

非线性约束函数g(x,u,t,a)和目标函数O(x,u,t,a)也可采用类似的方法表示为混沌多项式的展开形式，从而得到相应的展开式系数

各个系数向量计算完成后，可得到各个状态变量、非线性约束函数以及目标函数的均值、方差等统计特性.

4 仿真分析

本节以美国一可重复使用飞行器的再入轨迹优化为例，对第3 节和第4 节提出的确定性轨迹优化方法和考虑气动参数扰动条件的随机性轨迹优化方法进行仿真分析验证.其质量m=104 305 kg，气动参考面积S=391.22 m2，气动模型如下

借鉴航天飞机的程序攻角设计方法，采用典型的二次分段标称攻角剖面，其具体表达形式为

仿真中暂不考虑飞行时间对轨迹的影响，规定总再入时间tf=1600 s，离散点总数为200，为保证计算的效率，设定最高迭代次数为100 次.仿真所需参数见表2.仿真中，序列凸优化算法中信赖域约束和收敛条件如下

本文的仿真均是基于matlab 环境，调用Mosek8.0[34]的求解子凸问题的求解器，以初始状态和末端状态的直线为初始参考轨迹，大约迭代5∼7 次就基本收敛，说明该方法具有很强的鲁棒性.

表2 再入飞行参数Table 2 Reentry flight parameters

4.1 确定性再入轨迹优化结果分析

假设气动参数没有扰动，利用序列凸技术进行再入轨迹优化，由于凸优化具有多项式时间复杂度，且局部最优解即全局最优解等特点，因此该寻优时间比较短，大约为7 s.若在嵌入式C 语言环境下，计算时间大概减少数倍，具有在线轨迹计算的潜力.

图2(a)是优化结果中的速度--高度曲线和热流、动压、过载等路径约束的示意图，在提供再入点与终点连线的轨迹作为初始参考值的情况下，序列凸优化技术也很好地求解出满足各非线性约束的最优轨迹，表明该方法对初值的依赖度较低，鲁棒性较强.图3 是对控制量分别进行均匀离散和B 样条曲线离散的寻优结果，可以发现，本文采取的方法有效避免了数值计算中锯齿化现象，使得寻优结果更为光滑可行.

图2 序列凸优化优化结果Fig.2 Results of sequential convex optimization

图2 序列凸优化优化结果(续)Fig.2 Results of sequential convex optimization(continued)

图3 时间--控制量曲线Fig.3 Time-control curve

为验证序列凸优化的有效性和快速性，同时为计算状态量的的统计特性，比较广义混沌多项式与蒙特卡洛方法的特点，现对气动参数拉偏，进行400 次打靶测试，拉偏范围Ka=[0.7,1.3]，CL=，测试结果如图4 所示.

图4 打靶测试结果Fig.4 Results of target test

可以看出，序列凸优化可以适应各种拉偏情况，平均耗时7.14 s 左右，算法本身的鲁棒性较强.取8 个Gauss 配点的数据，通过最小二次回归求解广义混沌多项式的系数.利用广义多项式混沌GPC 和蒙特卡洛仿真MCS 分别计算满足均匀分布的气动力符合系数扰动条件下的优化轨迹各状态变量的统计特性.下图对比分析了不确定条件下状态变量的均值变化曲线与标准差变化曲线.可以看出，两种方法计算的各状态变量的均值变化曲线基本一致.

由图5(b)和图5(d)可知，在气动力扰动过程中，各状态变量的标准差均发生剧烈的变化.其中高度扰动峰值超过2000 m，倾侧角在飞行过程中也产生了超过10◦的扰动，说明气动参数的扰动对飞行状态产生了较大的影响.

图5 混沌多项式计算结果Fig.5 Results of chaotic polynomial calculation

图5 混沌多项式计算结果(续)Fig.5 Results of chaotic polynomial calculation(continued)

综合分析两种方法计算所得的统计特性，不难发现广义多项式混沌方法与传统蒙特卡洛仿真方法的计算结果基本一致.本节利用广义多项式混沌方法运用配点法计算统计特性的过程中，对于一维随机变量设置的配点个数为8，而蒙特卡洛仿真方法达到收敛精度的采样次数较大，说明广义多项式混沌方法对于计算气动参数扰动条件下状态变量的统计特性具有很大的计算优势.

4.2 考虑气动参数扰动的再入轨迹优化仿真分析

上一小节的仿真结果验证了序列凸优化的快速性和有效性，但是再入过程中如果存在气动参数扰动，若依然利用基于确定性模型的轨迹优化问题所得到的倾侧角控制量，路径约束如热流、动压等将存在违背的情况.因此本节主要研究在飞行器再入过程中考虑气动参数扰动的轨迹优化问题，力争降低最优轨迹对不确定扰动的敏感程度，提高寻优结果的鲁棒性，从而保证飞行任务的圆满成功.

由式(20)知，混沌多项式的复杂度随着展开项的增多而急剧增大，为保证计算效率，本节的优化过程中只考虑10%的气动参数扰动，并且基底取一阶Legendre 多项式，仅取3 个Gauss 点的样本值进行展开式的系数计算.考虑到鲁棒性的要求与控制能力有限的条件，现将最大热流密度调整至=16 kW/m2，式(27)中参数k=1.

图6 是利用确定性模型下所得到的控制量进行打靶测试所得到热流与动压的结果.可以发现，确定性模型的优化结果并不能满足所有扰动工况的路径约束，在实际飞行时，可能会使飞行器遭到损坏，导致任务失败.

图6 确定性优化的热流和动压约束Fig.6 and q of the deterministic optimization

为提高寻优结果的稳定性和鲁棒性，本节利用3.2 节提出的鲁棒轨迹优化方法进行仿真验证.在序列凸优化的技术上，充分考虑气动参数的干扰，利用混沌多项式计算得到非线性约束函数的统计特性，在寻优过程中将其包含在内，并将优化所得的控制量进行气动参数拉偏打靶测试，继而得到如图7 的打靶结果.

图7 确定性优化的热流和动压约束Fig.7 and q of the robust optimization

由图7 可知，经打靶测试优化所得的控制量可以满足所有既定工况的路径约束，通过广义混沌多项式所估计的非线性约束函数的统计特性基本符合事实，测试中所有的约束都在估计范围内，这也说明了混沌多项式在描述不确定度传播方面的优势.图7 中所估计的上下偏差虽然包含所有工况的约束，在一定程度上保证了满足非线性约束的裕度，但估计量仍与真实值存在偏差.这是由于计算能力的限制，仿真过程中只采用了一阶多项式且采样点仅有3 个Gauss 点，计算精度较为有限.同时需要说明的是，这里只考虑了10%的气动参数扰动，且认为该扰动服从均匀分布，虽然已验证了所提出的鲁棒轨迹优化方法的可行性，但若进一步增加其干扰程度，可能会导致无可行解的情况.另外本次仿真中，取k值为1，k值越大，考虑的干扰裕度越大即轨迹优化的鲁棒性越强，但最优性损失越多，在实际应用过程中，需要在最优性和鲁棒性之间作出权衡.

3 结论

本文以考虑气动参数扰动的再入鲁棒轨迹优化设计为主要研究内容，通过广义混沌多项式方法对随机变量的统计特性进行有效估计，同时结合一种改进的序列凸优化技术进行鲁棒轨迹优化设计策略.轨迹优化过程中以倾侧角的变化率作为控制量，为抑制数值优化过程中的锯齿化现象，采用B 样条曲线离散控制量，同时增加额外虚拟控制量，通过一种“回溯直线”搜索的方法，提高算法的稳定性、快速性和寻优结果的光滑性.在考虑气动参数扰动条件的寻优过程中，满足热流、动压等多重约束的同时，降低气动参数扰动对优化结果的影响，提出了基于广义混沌多项式与序列凸优化技术相结合的鲁棒轨迹优化方法.仿真结果表明，该策略可以有效降低优化结果对气动参数扰动的敏感程度，可以适应既定的所有工况，表现出较强的鲁棒稳定性，具有一定的工程应用价值.