基于自适应伪谱法的高超声速飞行器再入轨迹优化

2019-12-02任鹏飞王洪波周国峰

北京航空航天大学学报 2019年11期

任鹏飞，王洪波，周国峰

(中国运载火箭技术研究院, 北京 100076)

高超声速飞行器飞行速度高，机动范围大，具备入轨和再入大气的能力，可实现远程快速物资运送或精确打击，具有较高的经济军事价值，近年来已成为各国的研究热点[1]，其中具有代表性的为美国通用空天飞行器(Common Aero Vehicle, CAV)。再入飞行段作为此类飞行器重要的飞行阶段，由于飞行器的力热载荷环境严酷，且对于控制量高度敏感，因此再入轨迹优化设计也面临着一定挑战。

此类问题可以描述为具有多种复杂约束的非线性连续时间最优控制问题。在近20年里，直接配点法已经成为数值求解非线性最优控制问题的有力工具之一。直接配点法的基本思想是在所研究时域内合理的选取一系列不同的点，对状态变量和控制变量进行离散并建立一定的约束关系，将连续时间最优控制问题转化为有限维的非线性规划 (Nonlinear Programming,NLP)问题，从而通过著名的NLP求解器，例如SNOPT[2]、IPOPT[3]、 KNITRO[4]等进行有效求解。相比于间接法，直接配点法无需推导Pontryagin极小值原理的最优控制一阶必要条件[5]，降低了初值敏感性，减小了求解的难度，在飞行器轨迹优化等领域得到了广泛应用。

起初，直接配点法中发展了h方法，例如Runge-Kutta方法等，其通过将时间区间划分为多个子区间，在每个子区间上利用固定阶数的多项式来近似状态变量和控制，通过增加子区间数量来实现h方法的收敛[6]。

近年来，大量的研究都聚焦于采用一系列Gauss正交配点的方法[7-10]，称为伪谱法或正交配点法。由于高精度Gauss积分规则，LG(Legendre- Gauss)、LGR(Legendre-Gauss-Radau)以及LGL(Legendre-Gauss-Lobatto)这3种配点类型被广泛采用，它们分别是Legendre多项式的根、Legendre多项式线性组合的根以及Legendre多项式导数的根。伪谱法属于p方法，通常基于Gauss正交规则，在整个时间区间选取一系列配点，通常采用Lagrange全局插值多项式近似状态变量和控制变量，通过配置微分代数方程来近似微分约束，将连续时间最优控制问题转换为NLP问题。通过增加多项式的阶数来实现p方法的收敛[11]。

自适应伪谱法[12-14]结合了h方法和p方法的特点，通过估计当前解的相对误差，按照一定规则调整区间数量和多项式阶数，相比单独使用h方法或p方法，提高了求解的精度及效率，并能够有效应对状态控制等不连续的问题。

本文针对高超声速飞行器再入轨迹优化问题，建立了三自由度再入运动方程，以CAV-H飞行器为对象，考虑主要约束条件，基于Legendre多项式近似理论，建立了LGR伪谱法误差估计关系式，采用自适应网格重构技术，实现了3种典型性能指标再入轨迹优化问题的高效求解。本文算法相对传统自适应伪谱法具有配点和区间分配合理，收敛速度快及对人工参数不敏感等优点。

1 再入轨迹优化问题

1.1 三自由度再入运动方程

本文考虑地球为旋转圆球，飞行器无侧滑情况下，建立无动力高超声速飞行器再入运动方程：

(1)

式中：h、V、θ、φ、γ、ψ、m、ωe、σ和r分别为高度、速度、经度、纬度、航迹角、航向角、飞行器质量、地球自转角速度、倾侧角和地心距；引力常数μ=3.986×1014m3/s2；地球半径Re=6 378 137 m；L和D分别为升力和阻力，其表达式为

(2)

其中：ρ为大气密度；海平面大气密度ρ0=1.2 2 56kg/m3；归一化高度H0=7 110 m；升力系数CL和阻力系数CD由飞行器外形和飞行状态决定；S为参考面积。本文以CAV-H为例，相关总体和气动参数来自文献[15]，将气动系数拟合为马赫数Ma和迎角α的函数。

(3)

1.2 约束条件

1) 端点约束

为满足再入段末端能量管理要求，对末端高度hf、速度Vf及航迹角γf提出限制：

(4)

2) 过程约束

(5)

式中：KQ为与飞行器外形和材料相关的常数，本文取KQ=2.582 0×10-7;g0为海平面重力加速度。

3) 控制约束

为减轻姿态控制系统压力，要求迎角和倾侧角变化较为平缓，同时对其大小提出限制：

(6)

1.3 目标函数

考虑3种典型再入段轨迹优化目标，即最大纵程，最大横程及最小航迹角变化。通过引入权重系数表征设计偏好，则一般性目标函数可以表示为

(7)

式中：w1、w2、w3为权重系数；t0为初始时刻；tf为末端时刻。

2 Bolza型连续时间最优控制问题

考虑如下Bolza型连续时间最优控制问题：

minJ=Φ(x(t0),t0,x(tf),tf)+

(8)

式中：x(t)∈RNx为状态变量，Nx为状态变量维数；u(t)∈RNu为控制变量，Nu为控制变量维数。Mayer型性能指标Φ，Lagrange型性能指标L，动力学微分约束f，路径约束C及端点约束φ定义为

(9)

其中：NC和Nφ分别为路径约束和端点约束维数。

3 多区间LGR伪谱法

多区间伪谱法的思想是将整个时间区间划分为多个小的子区间，在每个子区间内采用较少数量的配点来离散最优控制问题。全局伪谱法则可认为是只有一个区间的特例。不失一般性，将t∈[t0,tf]分为s(s≥1)个子区间Sk=[tk-1,tk]，其中:k=1,2,…,s;t0

(10)

LGR伪谱法的配点采用τ∈[-1,1]上的N个LGR点或反转的LGR点，本文选用前者。LGR点为(τ+1)(PN(τ)+PN-1(τ))的根，其中:PN(τ)为N阶Legendre多项式。因此对t进行区间转化，可使每个子区间满足Legendre正交多项式的定义区间。

LGR配点和Gauss-Radau积分权重的计算方法主要为特征值方法和牛顿迭代法[16]。本文采用特征值方法进行求解，N-1阶Jacobi矩阵为

(11)

式中：

(12)

AN-1的第i个特征值为λi，对应单位特征向量的首个元素为v1i。根据Radau点与Gauss点的对应关系[16]，可得到LGR配点τi及积分权值ωi为

(13)

式中：i=1,2,…,Nk-1。

在每个时间子区间内采用Lagrange基函数近似状态变量和控制变量。状态变量的插值节点为Nk个LGR点和τNk=1点，即插值节点数比配点数多1个；控制变量的插值节点为Nk个LGR点，末端控制点通过外插得到。因此第k个子区间的状态变量和控制变量可近似表示为

(14)

式中：

(15)

状态变量在配点处对τ的微分可表示为

j= 0,1,…,Nk-1

(16)

式中：DNk×(Nk+1)为Radau微分矩阵，表示各Lagrange基函数在各LGR配点处的微分值，其表达式为[16]

(17)

Z(k)(τ)=(1-τ)(PNk(τ)+PNk-1(τ))

(18)

至此最优控制问题(8)可以离散为微分形式的NLP问题：

minJ=Φ(X(1)(-1),t0,X(s)(1),tf)+

(19)

4 自适应网格重构技术

自适应伪谱法的基本思想：基于当前解的最大相对误差来评估解的质量，当不满足求解精度要求时，按照一定规则调整子区间内多项式阶数或调整子区间数量，以达到改善解精度的目的。

4.1 误差评估

文献[13]给出了一种近似相对误差的算法，其基本思想是通过对比2种状态近似解，从而构建相对误差的近似表达式。

j=1,2,…,Nk+1

(20)

(21)

此时可以构建第i维状态变量绝对误差和相对误差的近似表达式为

(22)

则Sk内最大相对误差为

(23)

4.2 Legendre多项式近似

根据收敛性理论[17]，对于光滑问题，全局伪谱法以指数形式收敛。文献[17]指出，对于分段光滑函数，采用Legendre多项式的近似误差与多项式系数具有相同的衰减率，且衰减率与多项式的阶数相关。因此可以基于Legendre多项式系数构建多项式阶数与近似误差的关系。

对于分段光滑有界函数y(τ)可以展开为Legendre级数：

(24)

式中：pi为Legendre系数；Pi(τ)为Legendre多项式。取N阶近似为

(25)

则截断误差为

(26)

根据Legendre多项式的正交性：

(27)

式中：δmn为Kronecker函数。则式(26)化简为

(28)

根据文献[17]分段光滑有界函数y(τ)的Legendre系数近似为

(29)

(30)

因此可根据离散Legendre近似形式：

(31)

采用最小二乘法求解c和υ。

同时考虑到级数关系式：

(32)

因此可构建近似误差关系式：

(33)

4.3 光滑性判定

4.4 多项式阶数更新

当υ>υε时，根据式(34)可计算当前网格的求解误差：

(34)

(35)

(36)

为了严格增加配点数量，采用向上取整算子:

(37)

4.5 区间更新

当υ≤υε时，考虑到υ接近于0时，由式(33)确定区间分段数会产生很大的区间数量，从而导致NLP问题的规模过大，使得求解效率降低。因此采用υε代替υ，根据式(36)估计所需的多项式阶数：

(38)

令每个新子区间与当前子区间的配点数相同，则区间分段数Hk为

(39)

4.6 网格缩减

减少区间配点数的基本思想是将Lagrange插值多项式展开为τ的幂级数形式，并计算各阶次系数。为了减少该多项式阶数同时保证相对误差仍满足要求，考虑到τ∈[-1,1]，令τ=τmax=1，因此仅需从最高阶次系数向最低阶次系数依次累加，直至累加值超过容许误差ε，则可减少的配点数不超过累加次数。

合并相邻子区间的基本思想与上述算法类似，只是将相邻子区间的值在当前区间对应多项式上进行外插，然后类似地判断容许误差是否满足。

网格缩减算法的具体推导过程见文献[14]，本文不再赘述。

4.7 算法流程

本文自适应伪谱法的计算流程如下：

步骤1给定求解容许误差ε和初始网格。

步骤2求解NLP问题(19)。

步骤5重复步骤2。

综上，本文算法的计算流程图见图1。

图1 自适应伪谱法流程图Fig.1 Flowchart of adaptive pseudospectral method

5 算例分析

再入轨迹优化算例设置见表1，边界条件设置见表2。计算环境为Windows XP 2.83 GHz，3.25 GB内存，软件环境为MATLAB R2014a，NLP求解器采用SNOPT v7.0，梯度计算采用有限差分算法。

为了验证本文自适应伪谱法求解的有效性与快速性，采用本文算法(记为hpL)求解最优控制变量，并采用变步长Runge-Kutta-Fehlberg法(记为RKF45)进行积分验证，积分相对误差为10-13。对比算法分别采用文献[12]算法(记为hp)以及文献[13]算法(记为ph)。本文算法轨迹优化结果见图2～图6，不同算法对比结果见表3。

表1 算例设置Table 1 Setup of calculation examples

表2 边界条件设置Table 2 Setup of boundary conditions

由图2～图4可以发现，对于3种典型算例，采用hpL算法的求解结果与采用变步长RKF45算法积分的结果均保持一致。优化获得的飞行器再入轨迹末端均严格满足末端高度约束，速度约束以及航迹角约束。由图4可以发现，h-V平面内,3种算例的最优轨迹均位于再入走廊下边界的上方，即全程满足飞行过程中的热流率约束、法向过载约束以及动压约束要求。可以注意到，对于本文CAV-H飞行器，法向过载在再入段过程中均不是主导约束，再入初始阶段，飞行器飞行高度高，大气非常稀薄，即使飞行速度较高，相比于动压，热流率为主导约束，由图4可以发现，此时飞行器采用较大迎角以获得较大升力，从而有效减缓飞行高度的快速下降，以满足热流率要求。当飞行器能量下降至一定时，动压为主导约束。对于最小航迹角变化率问题，整个飞行轨迹平滑无跳跃，初始以较大迎角飞行，飞行速度快速下降以匹配平稳飞行的能量要求。由图5和图6可以发现，整个飞行过程中，3种算例的迎角和倾侧角均变化缓慢，其变化率均不超过1(°)/s。从而验证了本文算法的有效性。

图2 3D轨迹Fig.2 Three-dimensional trajectory

图3 航迹角变化Fig.3 Path angle versus time

图4 再入走廊边界Fig.4 Boundary of reentry corridor

图5 迎角与倾侧角变化Fig.5 Angle of attack and heeling angle versus time

图6 迎角与倾侧角变化率变化Fig.6 Change rate of angle of attack and heeling angle versus time

表3 不同算法轨迹优化结果Table 3 Trajectory optimization results by different methods

由表3可以发现，本文hpL算法具有较好的稳健性，人工参数衰减容限υε一般取值为(0,1]之间，本文设置为0.25和0.5时，3种算例的求解时间和网格迭代次数差别不大，原因是采用Legendre衰减系数算法进行光滑性判定时，衰减系数υ的区分度较大，使得人工参数在一定区间内取值时，对增加配点数和分割区间数的选择影响较小，因此轨迹优化效率对于人工参数的选择不敏感。对于hpL算法，增加的配点数取决于衰减系数υ，即与状态平滑性相关，而ph算法可认为是υ≡lgN的算法，当N<10时，ph算法估计的多项式阶数偏大，其求解效率受到区间最大配点数的影响较大，对于最大射程问题，最大配点数为6的求解时间仅为最大配点数为10的0.1倍，此时ph算法更倾向于分段，从而避免稠密的约束Jacobian矩阵。类似地，对于hp算法可认为是υ≡1的算法。本文phL算法对于3种算例均展现了较好的求解快速性，并且相比于hp算法和ph算法，hpL算法的网格迭代次数更少，收敛速度更快，同时采用网格缩减技术，使得配点和区间总数均相对更少。验证了基于Legendre衰减系数的误差估计的有效性以及本文算法的快速性。