樊朋飞,郭云鹤2,凡永华,闫 杰

(1.西北工业大学 航天学院，西安 710072; 2.上海机电工程研究所，上海 201109)

0 引言

高超声速滑翔飞行器具有高升阻比气动外形，能够在大气层内实现高速、远距离滑翔和大范围横向机动。因其落点覆盖范围大、机动能力强，难以进行轨迹预测和拦截，因而具有重大的军事应用前景。平衡滑翔式轨迹是指滑翔飞行器飞行中升力在纵向的投影与离心力、重力近似满足平衡关系的一类轨迹。该类轨迹具有良好的热流率曲线以及稳定的操控特性和机动能力，在轨迹优化算法的研究和工程应用中受到广泛关注[1]。

再入可达区域体现了飞行器的打击覆盖能力和机动性，是衡量其作战性能的重要指标，对该指标进行精确高效的计算对于打击任务规划、目标和航路点选取具有实用参考价值。文献[2]总结了计算可达区域的四种常用方法：轨迹优化法、剖面规划法、常值倾侧角法和椭圆近似法。轨迹优化法中，常采用伪谱法[3-4]、粒子群算法[5]、遗传算法[6]等优化算法求解可达域的边界轨迹。该类算法由于最优性能指标的保证，具有较高计算精度，但其缺点在于求解再入轨迹这种具有复杂状态约束和动力学模型的问题时，效率不高。同时，由于高速、高升阻比的特点，再入滑翔轨迹具有天然的“振荡特性”，平衡滑翔轨迹在数值轨迹优化算法中难以得到。文献[7]和[8]通过在动力学方程中加入振荡抑制反馈项以提升轨迹的阻尼特性，并采用伪谱法对改进后的动力学模型进行轨迹优化。文献[9]通过引入无损松弛技术，得到了具有仿射结构的动力学模型，以减弱控制量对动力学方程的耦合作用。在此基础上，文献[10]设计了光滑的标称航迹倾角剖面，并将其与实际倾角的偏差引入最优指标中，保证求解轨迹的光滑特性。

本文研究采用连续凸优化算法求解平衡滑翔式轨迹的再入可达区域问题。为提高优化求解的效率，首先引入“准平衡滑翔假设”条件对动力学方程进行降阶，进而以速度变量替代时间作为自变量，使得动力学方程的阶次减小一倍。通过这种方法，不仅大大降低了最优轨迹的求解难度，而且消除了高升阻比飞行器浮沉特性对于控制输入敏感的问题，此时，轨迹的平滑度仅与控制量的平滑程度相关，而无需引入额外的阻尼项。随后，对优化轨迹的求解问题进行了转化，通过将降阶的动力学方程进行线性化、离散化处理，并引入二阶锥约束条件对控制量的平滑度进行约束，将非线性最优控制问题转化为可被高效求解的二阶锥规划(SOCP)问题，通过连续求解一系列凸优化子问题获得最优的可达域边界轨迹。以CAV-H高超声速滑翔飞行器为模型的可达区域求解算例验证了方法的有效性。

1 再入飞行器动力学建模

忽略地球自转的影响，将地球视为静止的球体，建立再入滑翔飞行器的质心动力学方程如下：

(1)

(2)

(3)

(4)

(5)

(6)

式中,状态量V，γ，ψ，r，θ和φ分别表示飞行器的速度、航迹倾角、航向角、地心距、经度和纬度；m为飞行器质量；g为重力加速度；σ为倾侧角。D和L分别为阻力和升力，其计算公式为:

(7)

(8)

再入过程中，飞行器受到再入走廊的约束，即满足：

(9)

(10)

(11)

式(9)～(11)分别表示热流率约束、法向过载约束和动压约束。α为攻角；CQ为热流常数。

为使轨迹末端满足一定的交接班条件，需对末端高度和速度进行约束，其形式如下：

r(tf)-R0=Hf

V(tf)=Vf

(12)

式中,R0为地球半径。

2 平衡滑翔轨迹可达区域计算

采用轨迹优化方法计算飞行器的可达区域时，通常将问题转化一系列纵程固定，横程指标最大的最优控制问题。首先，计算可达区域的远界和近界，得到纵程的取值范围；然后，在纵程区间内依次取值作为末端状态约束，并求解最优控制使末端横程指标最大。不失一般性，本文取0°经度、纬度和90°航向角作为初始位置航向条件，则飞行器的经度、纬度值可直接作为纵程和横程使用，且可达区域具有关于纵轴对称的特性。该条件计算得到的可达区域可经球面坐标变换[11]后投影至其他初始位置航向条件的情况。

2.1 运动方程简化

(13)

根据升阻比关系，计算阻力：

(14)

则简化后的动力学方程为：

(15～18)

考虑到可达区域的计算不关心飞行时间，而速度V的初始和末端状态固定，将上述方程组对V求导，消去速度方程和时间变量，得到以V为自变量的运动方程：

(19～21)

式中，下标V表示状态量相对速度的导数。此时，方程中状态量仅剩ψ、θ和φ；升阻比CL/CD为马赫数Ma和攻角α的函数，由于再入飞行器的攻角通常采用固定剖面，若忽略高空声速的变化，则CL/CD仅与飞行速度相关。在方程中的状态量和控制量σ确定后，方程中未包含的状态信息H和γ也可计算得出。其中H可通过联立式(7)和式(13)求解：

(22)

将上式对V求导：

(23)

根据式(1)与式(4)，相除后又可得：

(24)

联立上两式得：

(25)

式中,Dm=D/m为阻力方向的加速度。

2.2 约束条件转化

再入轨迹的优化中，再入走廊通常以V-H或者V-Dm走廊的形式出现。而本文中，H或Dm并不显含于运动方程(19～21)中，基于 “准平衡滑翔假设”，上述形式的约束可通过式(13)或式(14)转化为V-σ形式的约束，即满足：

|σ|≤σr(V)

(26)

式中,σr表示平衡滑翔状态下不超出再入走廊约束的倾侧角最大值。同理，初始和末端高度约束可转化为初始、末端倾侧角约束：

σ0=σ(H0,V0)

σf=σ(Hf,Vf)

(27)

2.3 优化问题求解

本节将最优轨迹求解问题由非线性规划问题转化为连续求解的一系列SOCP[12]问题。SOCP是凸优化问题的一种特例，其要求性能指标函数为优化变量的线性组合，受到线性约束和二阶锥约束。由于该类问题可运用对偶内点法[13](primal-dual interior-point algorithm)高效求解，因此，将文中所述优化问题转化为可被SOCP方法求解的标准形式成为问题的关键。将非线性运动方程写为如下一般形式：

(28)

首先，对非线性动力学模型进行线性化、离散化处理，使其满足线性等式约束形式。令{xk(V);σk(V)}表示连续求解过程中的第k次解，将动力学方程在其附近一阶展开可得：

(29)

其中：

c(xk,σk,V)=f(xk,σk,V)-A(xk,σk,V)xk-

B(xk,σk,V)σk

(30)

随后，将速度区间[V0,Vf]等分为N+1个子区间，则其步长ΔV=(Vf-V0)/N，离散点可表示为{V0,V1,V2,…,VN}，其中，Vi=V0+iΔV。则状态量x和控制量σ被离散化为xi=x(Vi)，σi=σ(Vi)。状态方程可进行数值积分转化为如下形式：

(31)

(32)

Mz=F

(33)

其中：

(34)

(35)

经过线性化和离散化后的最优问题可描述为：在离散化的运动方程(33)约束下，求取最优控制，使横程指标J=-φ(Vf)最小，并满足，

1)初始状态约束:

x(V0)=x0

(36)

2)末端经度约束:

θ(Vf)=θf

(37)

3)再入走廊约束:

|σi|≤σr(Vi)，i=0,…,N

(38)

以及式(27)中的初始和末端倾侧角约束：

σ(V0)=σ0

σ(Vf)=σf

式(37)中，θf∈[θfmin,θfmax]。

可以看出，此时最优问题只包含线性的约束条件和性能指标，符合SOCP问题的求解要求。为改善最优问题的可解性，并且保证最优轨迹的平滑特性，对上述最优问题的指标和约束项进行改进。其中最优指标变为：

J=-φ(Vf)+kθθa+kσσvar

(39)

式中,θa、σvar为新增优化变量，kθ和kσ为其权重系数。θa表示对末端经度指标的满足程度，有如下不等式关系：

|θ(Vf)-θf|≤θa

(40)

σvar用于增加控制量的平滑度，其与离散控制量满足不等式关系：

(41)

考虑到线性化处理后，最优解的可信域问题，应对每次求解结果的范围进行约束，即需要满足：

|xk+1(Vi)-xk(Vi)|≤δ，i=0,…,N

(42)

式中,取δ=[10 5 5]T。

(43)

最优轨迹求解时，首先给出一组初始猜测值{x0(V);σ0(V)}，并反复将求解结果作为下一次求解的输入值，当满足：

max|xk+1(Vi)-xk(Vi)|≤ε，i=0,…,N

(44)

时，则认为连续凸优化过程已收敛至最优解。式中：ε为一小量。

3 仿真分析

采用CAV-H的气动、总体参数建立再入动力学模型。仿真算例的初始条件为：H0=56 km，V0=6 000 m/s。

计算程序在MATLAB环境中编写实现，运用YALMIP软件[14]进行问题建模并调用MOSEK软件包[13]中的SOCP问题求解器进行求解。

图1 再入可达区域

图2 倾侧角曲线

图3 高度曲线

再入可达区域的计算结果如图1～4所示。其中，图1为可达域边界轨迹的经纬度信息，从中可以看出，可达域形状呈近似椭圆形分布，在上述条件约束下，可达域纵向远边界达到69.6°，侧向边界在经度45°处取得最大，为25.9°。

图2～4分别为轨迹的倾侧角、高度和阻力加速度信息随速度的变化曲线。从中可以看出，随着纵程的减小，轨迹的倾侧角逐渐增大，且轨迹的初段倾向于以最大倾侧角变化率增加至某一较大值，随后逐渐减小，至末端时迅速收敛至满足末端倾侧角约束。从高度曲线图可以看出，最优轨迹倾向于在初段以较低的高度飞行，以尽快调整轨迹偏角，而在中、末端以较高高度飞行以增大航程并满足末端状态约束，这与倾侧角曲线图的表现是一致的。

图4 阻力加速度曲线

为研究可达域范围随末端速度的变化规律，分别计算了末速Vf=2 000 m/s和Vf=1 000 m/s时的可达区域，并与图1中Vf=1 500 m/s时的结果进行了对比，其结果如图5所示。从中可以看出，可达区域的范围随末速减小而显著增大。其最大纵、横程比较结果如表1所示。可以看出，末速度变化相同数值时，可达域的最大纵、横程值增加量随末速度的减小而减小，这是由于随着速度减小，相同速度增量带来的航程收益逐渐变小。

图5 可达区域随末速度变化图

末速/(m/s)纵程/(°)横程/(°)100070.9726.82150069.5825.88200067.3024.19

为验证本文方法的有效性和优越性，基于末速Vf=1 500 m/s的算例，将采用本文方法与文献[2]中叙述的工程常用的剖面规划法进行了对比，其结果如图6所示。从中可以看出，两种方法得到的可达区域形状基本相符。在可达远边界和最大横程的求解方面，本文方法都优于剖面规划法。需要指出的是，由于本文采用了“准平衡滑翔”假设简化运动方程，在求取较近射程的最优轨迹时，所得横程边界小于剖面规划法，且可达区域的近界亦较远。这是因为对于较近射程的轨迹，其状态变化比较剧烈，往往不满足平衡假设条件，此时引入该约束则限制了轨迹变化的幅度，使得轨迹近界的计算不如采用完整动力学模型精确。但是由于远边界的计算相对极限近界而言更受关注，因此本文中的假设条件仍是合理的。

图6 两种方法的可达区域对比

4 结束语

本文研究了再入滑翔飞行器平衡滑翔轨迹的可达区域计算问题。基于“准平衡滑翔”假设简化了运动模型，并采用连续凸优化方法对再入走廊约束下的可达区域进行了求解。仿真结果表明：本文方法能够降低最优轨迹的求解难度，获得光滑无振荡的平衡滑翔轨迹；所计算得到的可达区域远界和横向边界都优于常规工程方法，体现了该方法的有效性和优越性。