基于伪谱法的滑翔轨迹多级迭代优化策略

(西南科技大学 计算机科学与技术学院，四川 绵阳 621010)

0 引言

自中国进入“超级航空2018”，就预示着在未来几十年中我国将在航空领域迎来更大的突破，而高超声速技术则是突破航空领域的重要研究点。飞行器最优轨迹规划是高超声速技术研究的关键点之一，其求解方法主要有两大类：间接法和直接法[1]。伪谱法以较少的参数、较高的计算精度、指数级的收敛速度，同时兼具了直接法和间接法的优势，被认为是目前最具潜力的轨迹最优控制求解方法[2]。

直接求解飞行器轨迹最优控制问题难度较大，利用伪谱法将轨迹最优控制问题进行离散，对离散化得到的非线性规划(Nonlinear Programming，NLP)问题进行插值获得原问题的解，可将求解最优轨迹控制问题上的难度转移到NLP问题的求解上，转化得到的NLP问题表现出优化空间大、约束控制量多等特点[3-4]。近几年，在选择求解该类NLP问题的策略上，很多研究者偏向使用序列二次规划(Sequential Quadratic Programming,SQP)算法[5]来求解。由于SQP算法可以将一些约束添加到变量中去，对于强约束问题的求解有很大优势。近两年有研究者使用SQP算法对高超声速飞行器周期巡航轨迹[6]、爬升轨迹[7]进行求解，取得了较好的结果。

作为求解强约束NLP问题的主流方法，SQP算法具有计算速度快，求解精度高等优势，但当NLP问题规模过大时，使用SQP算法求解存在初值敏感，收敛速度减慢的情况。因此求解较大规模的NLP问题时，选择一个可靠的初值是十分必要的。K.Okuda等人[8]在求解带翼火箭弹道优化时，提出使用具有全局最优性的遗传算法为SQP算法提供初始值，有效的解决初值选取问题。但遗传算法等全局优化算法在求解约束性较强的NLP问题时，存在精度低，收敛速度慢，结果随机抖动等问题。

本文结合伪谱法与SQP算法各自的优势，提出基于伪谱法的滑翔轨迹多级迭代优化策略。以X-51A[9]相似飞行器模型为研究对象，利用增量法与查表插值建立飞行器纵向气动力模型。使用伪谱法与SQP算法对滑翔段轨迹进行求解，提出多级迭代优化策略弥补轨迹规划过程中，SQP算法求解大规模NLP问题的不足之处，并与传统方法求解出的状态量与控制量仿真飞行状态进行对比。仿真结果表明本文提出的多级迭代优化策略稳定性好，计算精度高，收敛速度快，在实际工程应用中起到了很好的支撑作用。

1 滑翔段建模

1.1 气动模型

X-51A相似飞行器是气动推进一体化飞行器，如图1所示，模型总长4.26 m，最大宽度0.58 m，最大高度0.566 m，翼展最大0.927 m，两舵面与水平面成45°，内流道宽0.225 m，进气道总收缩比为4.9，内收缩比为2.4。

图1 一体化飞行器数模

本文借鉴意大利航空航天中心(CIRA)[10]利用增量法实现PRORA-USV1高超声速飞行器风洞试验数据与数值计算数据的融合建模思路，使用数值方法模拟X-51A相似飞行器的数据集进行建模，为仿真实验提供数据基础。结合增量建模理论，在基准数据集与舵效数据集的基础上构建增量模型，流程图如图2所示。

图2 增量模型建模流程图

(1)

(2)

式中，ρ∞为大气密度，L为升力，D为阻力，m表示飞行器质量，Sref=1m2,bref=1m2，V∞表示飞行速度。

1.2 滑翔段轨迹描述

滑翔段轨迹规划是X-51A相似飞行器轨迹优化的重要阶段，滑翔飞行能延长飞行器整体航程，实现远距离投送或打击。由于高超声速飞行器的动力学系统较为复杂，其建模所涉及的因素繁多，为研究需要，将地球视作一个均匀球体，认为滑翔飞行过程中飞行器是一个质点，质量保持不变，并忽略外界环境的干扰。滑翔段的运动学方程形式含多个状态变量与非线性等式约束，且方程中隐含有控制变量。本文以最大滑翔距离作为优化目标，其余指标作为约束条件，构建最优轨迹控制问题，以分析X-51A相似飞行器的滑翔弹道轨迹，运动模型表达式如下：

(3)

式中，x=[H,V,γ]为轨迹状态变量，分别表示飞行高度、飞行速度、飞行航迹角；u=α为轨迹的控制变量，α为攻角；tf为终端时间；D为阻力方向的合力，L为升力，m为飞行器质量，g为重力加速度，δ为最优控制变量α对应飞行姿态下保持纵向平衡的俯仰舵偏。J为目标函数，表示滑翔段的最大飞行距离。

飞行器滑翔段为无动力飞行，为保证平衡飞行，在滑翔段轨迹设计过程中要考虑相应的约束条件，本文的参数约束如表1所示。

表1 滑翔段约束条件

2 伪谱法参数化

伪谱法源于求解微分方程的谱方法，其思想是在离散的状态变量与控制变量基础上，构造逼近式，求得微分矩阵与积分算子，将最优控制问题转换为NLP问题,其核心技术包括：离散网格划分策略的选择；离散点插值函数对连续函数的逼近；插值函数的微积分对连续函数的微积分的近似。伪谱法参数化过程细化为以下几步。

1)时域变换：

本文使用Legendre多项式作为插值基，其系数{Ln(τ)}是τ∈[-1,1]上的正交多项式。当n≥1时满足以下递推关系：

(4)

由于插值多项式的作用域为τ∈[-1,1]，所以需要将t∈[0,tf]时域变换到τ∈[-1,1]：

(5)

2)变量离散化：

通过离散网格划分策略选取合适的配点，能有效避免逼近可能出现的Runge现象，通过Gauss-Lobato求积，取n个节点将问题离散，获得离散时间：Tn=[τ1,τ2…τn+1]。离散时间对应的状态变量xn=[Hn,Vn,γn]，Hn=[h1,h2…hn+1]，Vn=[v1,v2…vn+1]，γn=[γ1,γ2…γn+1]，mn=[m1,m2…mn+1]；控制变量un=[αn]，αn=[α1,α2…αn+1]。

3)LGL配点近似逼近状态与控制变量：

(6)

其中，Xi∈Rnx,Ui∈Rnu，φi(τ)为n阶Lagrange插值基数：

(7)

4)最优控制转化为NLP问题：

通过前几个步骤，可得到参数化后的NLP问题，具体表达式如下所示：

(8)

其中: [Hn,Vn,γn]为轨迹状态变量的离散向量，n为伪谱配点数；D,L分别为各离散状态量下的阻力与升力；D为微分矩阵；wi为积分算子在对应第i个状态变量下的积分权重：

(9)

3 非线性规划问题求解

3.1 多级迭代优化策略

初始值的选取对于NLP问题的求解来说具有很大影响。不可靠的初值，不仅会影响收敛速度，甚至有可能导致无法收敛，找不到可行解或最优解的情况。并且随着NLP问题规模的增大，约束性的强增，初值对求解结果的影响也会不断增大。伪谱法求解飞行器最优轨迹的精度是由配点数量的多少所决定的，一般情况下配点数越多，求解精度越高。但是转化得到的NLP问题的规模也会变得更大，在使用SQP算法求解时，会出现初值敏感、收敛速度减慢、陷入局部最优的情况。而配点数过少，则会导致伪谱法求解精度无法满足要求。

本文在滑翔轨迹问题离散过程中取n个配点时，会形成n+1个离散变量，对应决策向量的维度为4n+5，其中状态变量与控制变量构成4(n+1)维决策变量，最后一维决策变量为终端时间tf，采用一个多级迭代优化策略，对含高维度决策变量的NLP问题求解初值进行选取。

首先设置较少伪谱配点，保证问题规模不至于影响SQP算法的求解，对于得到的结果进行插值，然后将插值后的结果作为增加配点后的大规模NLP问题的初始值。详细过程如下：

1)在变量上下边界条件之间进行初始猜测插值，获得一个初始值x0；

2)设置少量伪谱配点，通过伪谱法参数化构建中小规模NLP问题，以x0作为初始值，使用SQP算法求解NLP问题得到粗略解x1；

3)增加配点数，对获得的粗略解x1进行插值，将插值结果作为增大规模后的NLP问题的初始值，求解NLP问题，重复，直到收敛。

多级迭代优化策略流程图如图3所示。

图3 多级迭代优化策略流程图

在初始猜测插值过程中，由于控制变量不光滑，如果与状态变量一样选择使用Lagrange多项式插值，其变化过大时会出现Gipps现象，导致控制变量超出边界。分段三次Hermite插值多项式不仅要求在节点上的值要等于目标函数值，而且要求在节点上的导数值也要与目标函数的导数值相等，光滑性更好。选择使用分段三次Hermite插值对离散的控制变量进行插值，可有效避免出现Gipps现象。

3.2 NLP问题求解

SQP算法是求解强约束NLP问题最有效的方法之一，其思想是将原NLP问题拆分为多个二次规划子问题，通过逐步求解二次规划子问题的解，为主迭代过程提供搜索方向，保证主迭代逐渐收敛于NLP问题的最优解。SQP算法的简要计算步骤如下：

1)初值设置：选择初始点x0，初始矩阵H1是n×n的正定矩阵，令k=1，计算精度0<ε<Δ；

2)求解二次规划子问题，确定下次迭代搜索方向d：

(10)

得到满足原问题Kuhn-Tucker条件的近似解(d(k),λ(k));

3)若‖d(k)‖≤ε，则x(k+1)=x(k)+d(k)为求得的最优解，停止计算；

4)如果P(x(k)+d(k),r)

5)一维搜索：对minP(x(k)+λd(k),r)进行求解，得到最优步长λ(k)，令x(k+1)=x(k)+λd(k)；

6)更新矩阵H(k)：根据拟牛顿公式更新Hessen矩阵H(k)，确定新的正定矩阵H(k+1)，令k=k+1，转到第2)步。

4 实验结果与分析

通过MATLAB进行仿真验证，仿真实验过程中设置X-51A相似飞行器质量m=1 600 kg，考虑到滑翔飞行的终点与投射、打击目标之间还存在一定的距离，因此需要对高度与速度设置相应的终端条件。设置初始速度V0=6 mach，初始高度H0=27 000 m，终端速度Vf=3 mach，终端高度Hf=1 000 m，为得到最远滑翔航程，不对终端时间做限制。

经过反复计算，整理结果如表2～3所示。表2中n为配点数，z为对应的决策变量数，t(s)为使用传统方法求解对应配点数规模的NLP问题计算时间，titer(s)为使用多级迭代优化策略下的计算时间。通过对比可以看出随着伪谱法配点个数的增加，计算所需要的时间也在快速增加，多级迭代优化策略可有效的减少计算时间。配点数越多，节省的计算时间越多，体现出的效果越好。

仿真实验过程中设置的最大迭代次数为1 000，表3中Kopt为最优迭代次数，L为计算求得的目标值。通过分析比较，本文提出的多级迭代优化策略能大幅度减少迭代次数，加快收敛速度，获得更精确的解。图4、图5为配点数n=40时，传统方法与使用多级迭代优化策略计算下，目标值的迭代情况。可以看到，使用多级迭代优化策略时，SQP算法能更快的收敛得到最优值，减少计算时间。

表2 计算时间对比情况

表3 目标值与迭代次数

图4 传统方法目标值迭代情况

图5 多级迭代优化策略目标值迭代情况

SQP算法对伪谱法的离散配点数量敏感，配点数越多决策变量维度越高，算法收敛时间越长。但配点数过少，则会影响结果精度。对不同配点数求得的结果，进行反解常微分方程，结果表明配点数过少时，求得的结果会出现很大偏差，情况如图6所示。

图6 配点数对解的影响

由于控制变量的不光滑性，使用Lagrange多项式插值，在变量变化较大的情况下会出现Gipps现象，造成控制变量超出边界，情况如图7所示。采用分段三次Hermite插值方法对控制变量进行插值，能保证光滑性，有效避免Gipps现象与Runge现象。

图7 PCHIP插值克服Gipps现象

通过仿真实验，对比本文方法与传统方法求解的状态变量与俯仰舵偏的飞行状态。结果表明多级迭代优化策略相比于传统方法，在求解精度和收敛速度上有很大的优势，图8～11对比了两种方法下配点数n=40时，状态变量和俯仰舵偏角的变化情况，虚线为使用传统方法的变化情况。传统方法求解下，飞行高度与终端条件设置的高度相差大约1818 m，速度相差大约0.3 mach。多级迭代优化策略方法求解下，飞行高度与终端条件设置的高度相差大约86 m，速度相差大约0.01 mach。本文方法与传统方法相比，求解偏差降低了很多，所得结果更加接近最优轨迹。并且在多级迭代优化策略求解下，飞行器滑翔时间更长，得到的滑翔水平距离比传统方法远大约220 m，能获得更远航程，实现更远距离的投送或打击。通过以上各项指标，证明了本文提出的多级迭代优化策略的有效性与高效性。

图8 滑翔轨迹对比图

图9 速度变化对比图

图10 航迹角变化对比图

图11 俯仰舵偏变化对比图

5 结束语

本文提出的多级迭代优化策略有效的结合了伪谱法参数化最优轨迹控制问题与SQP算法求解强约束非线性规划问题的优势。伪谱法参数化后的NLP问题具有大规模、强约束等特征，直接使用SQP算法求解时存在对初值敏感、收敛速度减慢等问题。多级迭代优化策略，为求解大规模NLP问题提供一个较为有效的初始值，可有效解决SQP算法在求解过程中存在的不足。仿真验证表明，在多级迭代优化策略下，NLP求解算法求解精度更高，计算时间更少。且该方法在实际工程应用中取得了良好效果。

随着高超声速技术的发展，全轨迹优化的研究也在火热进行中，所面对的问题规模将更大，约束性将更强，多级迭代策略体现的价值也将更大。