基于改进的序列梯度–修复算法的飞行器上升段轨迹优化

An Improved Gradient-Restoration Sequential Algorithm for Endoatmospheric Ascent Trajectory Optimization

摘要: 针对飞行器大气层内上升段轨迹优化问题，提出了一种改进的序列梯度–修复算法。根据序列梯度–修复算法仅可以处理等式约束的特点，给出了不等式过程约束的转化过程，同时为了解决飞行器大气层内上升段运动数学模型的强非线性带来的状态量求解更新困难的难题，引入了状态积分，改进了算法的更新方法，最后对改进的算法在飞行器大气层内上升段轨迹优化中的实用性进行了仿真分析。仿真结果表明改进的序列梯度–修复算法能够较好的获得满足末端约束和过程约束的大气层内上升段的最优轨迹，具有较强的适应性。

Abstract: An improved sequential gradient-restoration algorithm is developed to solve endoatmospheric ascent trajectory optimization problem. First, the transformation process of trajectory constraints is introduced to accommodate the character that sequential gradient-restoration algorithm can deal with equation constraints only. Then the state integral method is introduced and the update method of the sequential gradient-restoration algorithm is improved to deal with the problem of state solve and update, due to the strong nonlinearity of endoatmospheric ascent dynamics. Finally, the simulation analysis of the improved algorithm ascent trajectory optimization in aircraft atmosphere of practicability is analyzed through simulation under five different scenarios. Simulation results demonstrate the validity and adaptability of the improved algorithm for trajectory optimization with terminal constraints and process constraints.

1. 引言

飞行器在大气层内飞行受到空气动力、推力和引力的作用，同时还需考虑动压、过载等约束，使得大气层内飞行器的运动数学模型较为复杂，求解大气层内的最优化飞行轨迹具有较大的计算量和难度。

在这种背景下，吸引了大批学者的研究兴趣，出现了一些飞行器大气层内上升段最优轨迹快速解算方法。综合现有方法可以发现，直接法求解速度较慢，对初值敏感，间接法求解速度较快，但对过程约束的处理较为复杂 [1] - [4] 。

序列梯度–修复算法是一种求解最优控制问题的间接法，由梯度和修复两大模块组成。采用统一的最优性条件，减少了推导伴随方程和横截条件等过程的复杂和繁琐，便于模块化实现。同时，由于修复模块的存在，其对初值不敏感 [5] - [10] 。鉴于序列–梯度修复算法的优点，本文研究了序列梯度–修复算法在飞行器上升段轨迹优化设计中的应用。给出了过程约束的转化过程，同时结合飞行器上升段运动数学模型的强非线性，引入了状态积分，改进了算法的更新方法，最后对改进的算法在飞行器上升段轨迹优化中的实用性进行了仿真分析。

2. 大气层内上升段无量纲化的运动数学模型

为了提高算法的快速性和收敛性，在发射惯性坐标系上建立飞行器大气层内上升段的无量纲化运动数学模型为：

(1)

式中：和分别为飞行器相对地心的无量纲化矢径以及其对应的值；为飞行器无量纲化的绝对速度矢量；为无量纲化的推力；和分别为无量纲化的气动轴向力和法向力；和分别为飞行器纵向和法向的单位矢量；为飞行器质量。

3. 序列梯度–修复算法求解的最优控制问题模型

序列梯度–修复算法求解的最优控制问题模型可描述为 [5] ：

(1) 性能指标

性能指标为：

(2)

(2) 状态方程

状态方程为：

(3)

(3) 初始条件

初始条件为：

(4)

(4) 终端约束

终端约束为：

(5)

(5) 过程约束

过程约束为：

(6)

上述式中：性能指标由积分型函数和末端型函数组成；为n维状态微分函数；为n维状态变量；为k维控制变量；为i维参数向量；为时间变量；为j维终端约束函数，且有；表示状态量的末端值；表示状态量的初值，式(4)表示状态量的初值为给定值；为维过程约束；对于任意积分区间，可将看成参量，转化成0~1区间。

4. 大气层内上升段轨迹优化模型

为了使飞行器射程远、机动能力强，其大气层内上升段最优轨迹的求解问题，实际上就是寻找最优的俯仰程序角，使飞行器在满足动压、过载等约束的同时以最大的末端速度到达交班点的最优控制问题。

同时从第2节中可以看出序列梯度–修复算法仅可以求解过程约束为等式约束的最优控制问题，对于非等式约束的最优控制问题需要进行相应的处理，变成等式约束。相应的处理过程和最终获得的飞行器大气层内上升段轨迹优化模型如下文所述。

4.1. 过程约束的处理

(1) 动压约束

动压约束为：

(7)

式中：为大气密度；为飞行器的相对速度；为给定的动压约束。引入辅助状态量，将式(7)变成等式约束，同时进行无量纲化处理有：

(8)

式中：为无量纲化的相对速度；为地球引力常量；为地球平均半径。由于上式(8)中的各参数不显含控制量，为了提高算法的收敛性，需对上式(8)求导有：

(9)

式中：为大气密度对无量纲化地心距的导数。再引入辅助控制量有：

(10)

(11)

(12)

为此，式(10)将作为过程约束的一个分量，式(11)将作为一个状态方程，状态量的初值由式(12)给出。

(2) 体系下法向过载约束

体系下法向过载约束为：

(13)

式中：为气动法向力；为法向过载约束；是标准海平面引力加速度。引入辅助控制量，将式(13)变成等式约束，同时进行无量纲化处理有：

(14)

由于上式(14)中无量纲化的气动法向力为攻角的函数，可认为其显含控制量，因此无量纲化的式(14)可作为过程约束的一个分量。

(3) 动压和攻角的乘积约束

动压和攻角的乘积约束为：

(15)

式中：为攻角；为动压和攻角的乘积约束。引入辅助控制量，将式(15)变成等式约束，同时进行无量纲化处理有：

(16)

由于上式(16)中显含控制量，因此可将式(16)作为过程约束的一个分量。

4.2. 大气层内上升段轨迹优化模型

结合上述过程约束的处理，以及飞行器上升段轨迹优化的目标可以得到大气层内上升段轨迹优化模型如下文所述。

(1) 性能指标

性能指标为：

(17)

式中：为主动段末端相对速度。对式(17)进行无量纲化处理得到性能指标为：

(18)

(2) 状态方程

结合飞行器上升段的数学模型和动压约束的辅助状态量可以得到轨迹优化模型的状态方程为：

(19)

(3) 过程约束

根据对过程约束的处理可以得到相应的过程约束为：

(20)

(4) 终端约束

主动段末端约束可取为：

(21)

(22)

式中：、为主动段末端高度和弹道倾角；、为主动段末端高度和弹道倾角约束。

(5) 控制量和参量

结合飞行器相应的辅助控制量可以得到，优化控制量为：

(23)

式中：为俯仰程序角。

参量为无量纲化的飞行时间，同时应对状态方程式(19)作相应的变换，将时间从区间化成0~1区间。

5. 改进的序列梯度–修复算法结构

以文献 [10] 中的序列梯度–修复算法为基础，考虑到飞行器上升段运动数学模型的强非线性，引入状态积分，改进算法的更新方法，相应的算法结构如下文所述。

5.1. 初值的给定

由于算法中修复模块的存在，算法对初值不敏感。可根据系统的相关知识，首先给出在初始点、末端点的值以及的值，然后根据积分节点数，由初始点值和末端点值线性插值得到中间各点的值，再积分状态方程式(3)得到，最终获得、、的初值。

5.2. 梯度阶段

(1) 若为第一次迭代则直接转入修复阶段，否则根据获得的、、的初值，在各积分节点上计算对状态变量、控制变量以及参量的导数、以及，状态微分函数对状态变量、控制变量以及参量的导数、以及，过程约束对状态变量、控制变量以及参量的导数、以及，在末端计算函数对状态变量、参量的导数、，约束对状态变量、参量的导数、。

(2) 取，，为脉冲函数，且有：

(24)

同时，则积分下式(25)和式(26)组成的方程次，即先由式(25)解得和，再由式(26)计算和，进而进行积分求解获得：、、、、，。、分别为第次积分给定的的初始值和的值；、、、分别为第次积分求解得到的、、、。

(25)

(26)

(3) 引入个待定系数，令

(27)

(28)

(29)

(30)

(31)

则有：

(32)

(33)

(34)

(35)

解式(32)~式(35)组成的方程组，可以得到个待定系数和，进而求解得到、、、、。

(4) 由下式(36)、式(37)求得、，而后积分状态方程式(3)获得。

(36)

(37)

式中：为梯度步长。

取增广函数和函数为：

(38)

(39)

式中：表示向量的模，即

采用二分法搜索梯度步长，即先给定，若搜索条件满足则结束，否则取，直至搜索条件满足。相应的搜索条件为

(40)

(41)

式中：为一小量。

(5) 将求得、、赋给、、准备下一次迭代。

5.3. 修复阶段

(1) 根据获得的、、，在各积分节点上计算状态微分函数对状态变量、控制变量以及参量的导数、以及，过程约束对状态变量、控制变量以及参量的导数、以及，在末端计算约束，及其对状态变量、参量的导数、。

(2) 取，，为脉冲函数，且有

(42)

同时，采用与梯度阶段相同的方法正向积分下式(43)和式(44)组成的方程组次得到：、、、、，。、分别为第次积分给定的的初始值和的值；、、、分别为第次积分所获得的、、、。

(43)

(44)

(3) 引入个待定系数，令

(45)

(46)

(47)

(48)

(49)

则有：

(50)

(51)

(52)

(53)

解式(50)~式(53)组成的方程组，可以得到个待定系数和，进而求解得到、、、、。

(4) 由下式(54)、式(55)求得、，而后积分状态方程获得。

(54)

(55)

式中：为修复步长。

采用二分法搜索修复步长，即先给定，若搜索条件满足则搜索结束，否则取，直至搜索条件满足。相应的搜索条件为：

(56)

(5) 判断式(57)是否成立，若成立则修复阶段结束，进入步(6)；否则，将、、的值赋给、、重复步(1)，重复进行修复阶段，直到式(57)成立。

(57)

式中：为小量，且有。

(6) 若算法为第一次迭代，也即对初始给定的初值进行修复，则将、、的值赋给、、进行梯度阶段；否则，判断是否成立，若成立则将、、的值赋给、、进行下一个梯度阶段，若则减少梯度阶段的步长，用式(36)、式(37)重新获得梯度阶段结束时的、、，重新进行修复阶段，直到成立。代表梯度阶段结束时刻的性能指标值，代表整个一次迭代开始时刻的性能指标值。

5.4. 收敛判断准则

整个序列梯度–修复算法结束的标志为式(58)中的两式同时成立，其中的、为梯度环节中的、。

(58)

式中：为小量。

6. 仿真分析

为了验证算法的适应性，取表1中五种不同的末端情况进行仿真分析。其中法向过载约束为1，动压约束为20 kPa，动压与攻角乘积约束为250 Pa·rad，仿真初始条件为：

(1) 相对发射性系的初始速度：288.0 m/s、49.0 m/s、198.0 m/s；

(2) 相对发射性系的初始位置：2594.0 m、222.0 m、1789.0 m。

最终的优化结果如图1~图4所示。

从图1~图4中可以看出，对于本文给出的五种上升情况，该方法均能够获得满足要求的大气层内上升段最优轨迹，对于不同的上升情况，其最终最优轨迹的末端约束和过程约束均能得到较好地满足，相应的动压均在20 kPa以内，动压和攻角的乘积均小于等于250 Pa·rad。

因此从上述仿真结果中可以看出，本文所给出的方法能够适应不同上升情况的优化分析，能够获得满足末端约束和过程约束的飞行器大气层内上升段的最优轨迹，具有较强的适应性和鲁棒性。

Table 1. Simulation examples

表1. 仿真算例

Figure 1. Pitch angle curve in Case 1 - Case 5

图1. Case 1~Case 5俯仰角曲线

Figure 2. Trajectory angle curve in Case 1 - Case 5

图2. Case 1~Case 5弹道倾角曲线

Figure 3. Altitude curve in Case 1 - Case 5

图3. Case 1~Case 5高度曲线

Figure 4. Process constraints curve in Case 1 - Case 5

图4. Case 1~Case 5过程约束曲线

7. 结论

本文研究了序列梯度–修复算法在飞行器上升段轨迹优化设计中的应用。给出了过程约束的转化过程，同时结合飞行器上升段运动数学模型的强非线性，引入了状态积分，改进了算法的更新方法，最后对改进的算法在飞行器上升段轨迹优化中的实用性进行了仿真分析，仿真结果表明：改进的序列梯度–修复算法均能够适应不同上升情况的优化分析，能够获得满足末端约束和过程约束的飞行器大气层内上升段的最优轨迹，具有较强的适应性和鲁棒性。