飞行器轨迹优化是指在特定的约束条件下, 寻找飞行器从初始点到目标点满足某种性能指标最优的运动轨迹。在数学本质上, 飞行器轨迹优化可以抽象为在包含微分方程、代数方程和不等式等约束下求解泛函极值的开环最优控制问题。一般将轨迹优化方法分为直接法和间接法, 间接法基于Pontryagin极大值原理将最优控制问题转换为Hamiltonian边值问题(HBVP), 其解的精度高, 满足一阶最优必要条件, 但HBVP问题的收敛域小, 需要估计协态变量的初值, 且难以求解存在路径约束的最优控制问题^[1-3]; 直接法无需求解最优必要条件, 而是将连续的最优控制问题离散并参数化, 用数值方法对性能指标寻优, 其解的收敛半径大, 近年来得到广泛应用^[4-7]。作为直接法中的一个分支, 全局伪谱法发展极为迅速, 宗群等^[7]将Gauss伪谱法与序列二次规划算法相结合, 对存在边值及加速度约束的轨迹优化问题进行了求解。相比于传统的间接法, 全局伪谱法收敛速度快, 然而在求解非光滑问题时, 其收敛速度变慢且无法捕捉状态突变点导致求解精度下降, 运用效果不理想^[8]。Jiang Zhao等^[9]通过分段优化以改善轨迹的非光滑特性, 然而在解出最优轨迹之前难以确定要分多少段以及在哪分段。赵吉松等^[10-11]基于局部配点法提出自适应网格细化算法, 在非光滑区域加密网格以提高求解精度, 能够有效控制问题规模, 但失去了全局快速收敛的优势。Darby等^[12-13]首次提出hp自适应伪谱法并将其应用到最优控制问题上, 该方法汲取了全局伪谱法和网格细分法的优点, 采用双层优化策略, 具有更合理的配点分布, 从而提高了算法的优化效率。文献[14-15]将hp自适应伪谱法分别应用于滑翔弹道快速优化问题和再入轨迹优化问题。

本文基于Radau伪谱法提出了一种改进的hp自适应伪谱法, 并将该算法应用于高超声速飞行器的最省燃料爬升问题。结果表明，与目前已有的伪谱法相比，改进的hp自适应伪谱法能够以较少的时间代价获得较高精度的解。

1 上升段轨迹优化数学模型

高超声速飞行器上升段轨迹优化问题的数学模型包括动力学模型、推力模型和气动力模型, 以及飞行器爬升过程中必须满足的路径约束条件。

1.1 动力学模型

仅考虑飞行器在纵向平面内的质心运动, 忽略地球曲率和自转的影响, 得到高超声速飞行器上升段的动力学方程如下

(1)

式中, [V, θ, x, y, m]^T为系统的状态变量, 分别表示飞行器的速度、航迹倾角、航程、飞行高度和质量; P为发动机推力, D为阻力, L为升力, α为攻角, I_sp为燃料比冲, μ为地球引力参数, R_e为地球半径。

1.2 推力模型

本文采用涡轮冲压发动机与超燃冲压发动机相结合的推力模型。当马赫数小于4时, 涡轮冲压发动机提供推力; 当马赫数大于4时, 超燃冲压发动机提供推力。推力系数和特征冲量的计算公式如下^[16]:

1) 推力系数C_p为

(2)

推力P为

(3)

2) 特征冲量I_sp为

(4)

式中, h为飞行高度, q为动压, A_e为进气口面积, M为飞行马赫数。

1.3 气动力模型

飞行器飞行过程中受到升力L和阻力D的作用, 计算公式如下

(5)

式中, S_ref为飞行器参考面积, ρ为大气密度, 计算公式为ρ=ρ₀e^-βh, 其中ρ₀为海平面大气密度, β为密度系数常数, h为飞行器高度。升力系数C_L和阻力系数C_D由攻角α和马赫数Ma表示

(6)

式中的系数C_Li和C_Di基于通用航空飞行器数据^[17]。

1.4 约束条件 1.4.1 热流率约束

爬升过程中, 为了满足严苛的热防护条件, 必须对热流率Q_d添加如下约束

(7)

式中, K_Q为热流率标准化常数。

1.4.2 动压约束

动压的上限由飞行器执行机构的铰链力矩确定, 下限值则是由超燃冲压发动机启动条件所限制。动压q的路径约束可表示为

(8)

1.4.3 过载约束

为保证飞行器的结构不被破坏, 需对过载进行约束, 即

(9)

式中, g₀为海平面重力加速度。

1.4.4 控制量和状态量约束

受机动性能等因素的影响, 飞行器在上升段飞行过程中, 要求满足一定的控制约束和状态约束。因此建立如下约束条件

(10)

1.4.5 终端约束

飞行器巡航段所需要的初始高度、速度和弹道倾角确定了其爬升段的终端约束, 即

(11)

1.5 目标函数

飞行器爬升过程中节省的燃料作为巡航段的燃油增加, 可以大大增加巡航段的飞行航程。据此, 选取上升段燃料最省为优化目标, 即

(12)

1.6 多区间轨迹优化问题描述

基于高超声速飞行器的动力学模型, 最省燃料轨迹优化问题的本质为寻找最优控制量α使得目标函数J最小, 同时满足多种约束条件((7)~(11)式)。为表述方便, 将该问题描述为一般形式的Bolza型最优控制问题:定义状态变量x=[h, V, θ, m]^T、控制变量u=α, 使代价函数最小化

(13)

约束条件为

(14)

式中:t₀为初始时间, t_f为终端时间; =f(·)表示动态约束; C(·)≤0表示路径约束; φ(·)=0表示边界约束。

多区间最优控制问题是指将上述Bolza型问题的时间区间[t₀, t_f]划分为K个子区间, 用t₀ < t₁ < … < t_K=t_f表示K+1个网点, N_k(k=1, …, K-1)表示第k个子区间[t_k-1, t_k]内的配点数, 并对每个区间进行下述时域变换

(15)

将时间区间t∈[t_k-1, t_k]变换为τ=[-1,1]。

2 Radau伪谱法

Radau伪谱法使用Legendre-Gauss-Radau(LGR)配点对最优问题进行离散化, 其NLP问题的最优解满足原连续时间最优控制问题的KKT一阶必要条件^[18], 且配点包含(-1)点, 更容易表示多区间最优控制问题的区间衔接, 即x(t_k^-)=x(t_k⁺)。因此本文选取Radau伪谱法作为hp自适应策略的伪谱方法。

在RPM中, 连续Bolza问题的状态变量近似为

(16)

(17)

式中, τ∈[-1, +1], L_j^(k) (τ), j=1, …N_k+1, 为拉格朗日基函数, (τ₁^(k), …, τ_Nk^(k))为LGR配点。代价函数即(13)式离散化为

(18)

式中, U_j^(k), j=1, …, N_k, 是控制量在LGR配点处的近似值, ω_j^(k)为高斯权重。同时将最优控制问题的微分约束、路径约束和边界约束在LGR配点处进行离散化, 可以得到如下微分-代数约束

(19)

式中，i=1, …, N_k, D_ij^(k)为Radau伪谱微分矩阵, 由拉格朗日基函数L_j^(k) (τ)对τ求导得到, 即

(20)

此外, 我们将X_Nk+1^(k)和X₁^(k+1)用同一变量表示, 以此来满足相邻区间的连续性。

3 改进的hp自适应伪谱法 3.1 误差评估准则

在数值求解过程中, 若使数值方法的解能够精确地近似原问题的确切解, 则需要在任意点处满足微分-代数约束方程。因此, 以微分-代数约束方程在配点间的满足程度作为解的近似误差评估准则。

已知第k个子区间[t_k-1, t_k]内的配点数为N_k(k=1, …, K-1), 且配点分布为N_k阶正交多项式的根, 选取M_k=N_k+1阶正交多项式的根(t₁^(k), …, t_Mk^(k))作为采样点, 则相邻配点之间均存在一个采样点。将M_k个采样点代入(7)式中, 得到采样点处的近似状态量X^(k)(t_l^(k)), l∈[1, …, M_k]。同时, 控制量也可通过构造插值多项式来近似, 即

(21)

将采样点t_l^(k))代入上式, 得到采样点处控制量的近似值U^(k)(t_l^(k)), l∈[1, …, M_k]。定义采样点t_l^(k))处状态量的绝对误差和相对误差分别为

(22)

(23)

式中, l=1, …, M_k+1, i=1, …, n_x, D_lj^(k)为M_k×M_k阶Radau微分矩阵。定义区间k内的最大相对误差为

(24)

3.2 多项式阶数的估计

依据3.1节的内容求得第k个子区间的最大相对误差e_max^(k), 设定容许误差ε, 若e_max^(k)≤ε, 则认为该区间内的解已达到预设精度; 否则, 需要通过细化区间或者增加配点数来提高求解精度。

对于光滑问题, 全局伪谱法以指数形式收敛, 其求解误差可近似为O(N^2.5-N), N为区间内的配点数^[19]。因此, 若将配点数增加为N+P, 解的误差相应减小为原来的N^-P。假定区间k的最大相对误差e_max^(k)>ε, 为满足所需精度, 增加的配点数P_k应满足, 即

(25)

为保证求得的配点数为整数, 对(25)式右边项向上取整, 即

(26)

3.3 hp自适应迭代策略

通过(26)式可以求得区间k所需的配点数为=N_k+P_k。设定子区间内配点数的最大值和最小值为N_min, N_max, 若≤N_max, 即所需配点数未超过设定的配点数上限, 增加区间k内的配点数为; 若>N_max, 则细化区间, 将该区间划分为更小的子区间。

为保证新子区间的配点数总和仍为, 设定细分子区间的个数为

(27)

每个区间内的配点数均设为配点数初始值N_min。

可以看出, 本文所提出的hp迭代策略优先使用p法来提高精度, 只有当配点数超过所设定的上限值时才进行区间划分, 划分后每个区间内继续采用p法, 以此规律不断迭代直到满足解的精度要求。

基于以上的迭代策略, hp自适应算法流程总结如下:

① 初始化。令迭代次数M=0, 设定区间的初始划分以及区间内配点数上下限N_min, N_max;

② 将最优控制问题离散化, 转化为NLP问题进行求解, 并按照3.1节内容计算解的相对最大误差e_max;

③ 若所有区间的e_max均小于ε或M>M_max, 停止迭代; 否则, 继续④;

④ 在e_max>ε的区间内采用3.3节的自适应策略更新网格, 直到所有的网格完成更新;

⑤ 令M=M+1, 转向步骤②进行下一次迭代。

4 仿真分析

本文以某高超声速飞行器为例, 飞行器初始质量为3 600 kg, 气动参考面积为0.6 m², 给定初始高度y₀=17 km, 初始速度v₀=1 100 m/s, 初始弹道倾角θ₀=4°, 终端高度y_f=27 km, 终端速度v_f=1 870 m/s, 终端航迹倾角θ_f=0°。爬升过程中状态变量、控制变量和路径约束的边界限制如下

优化计算过程中软件环境为:Windows7 32位操作系统和MATLAB R2013a;硬件环境为:Intel Core i5-2450处理器, 4G内存。在精度要求为10^-5时, 求得的爬升段最小燃料消耗为165.56 kg, 整个计算过程用时7.069 8 s, 优化结果如图 1~图 4中a)图的圆点所示。

为验证解的最优性, 将最优控制量通过三次样条插值代入到动力学方程中, 求得飞行器的实际飞行状态, 计算结果图 1~图 4中a)图的实线所示。离散点处的状态误差如图中星型线所示。图中可以看出两者的相对误差在10^-4～10^-6之间, 解的最优性得到验证。

在精度要求为ε的情况下, 采用3种不同的伪谱方法对同一航迹优化问题进行求解, 优化结果如表 2所示。

表中, p代表全局Radau伪谱法, h代表分段定阶Radau伪谱法, hp代表改进的hp自适应Radau伪谱法。分别表示初始子区间个数和区间内的初始配点数。K₀、N₀分别表示优化完成时的最大微分-代数约束的误差、计算时间区间个数、整个时间区间的配点总数、迭代次数以及目标函数。需要注意的是在采用全局伪谱法优化时, 由于求解问题的不光滑性, 精度要求较高时, 过多的配点导致优化时间过长, 无法得到优化结果。

分析表中数据, 我们可以得出如下结论:①精度要求较低时, 3种方法均可以优化出结果, 且在优化速度和配点总数上无明显差异; ②的精度要求下, hp伪谱法逐渐地显现出其优势, 不论是优化速度还是配点总数均优于其他2种方法; ③随着精度要求的提高, p方法仅通过增加配点数无法得到最优解, 而h方法和hp方法成功地克服了不光滑带来的问题, 得到最优解。两者相对而言, hp伪谱法在配点总数, 子区间划分和优化速度上均优于h方法。

图 5给出了10^-3精度下不同伪谱法的配点分布(图中纵坐标无实际物理意义, 仅为图示方便, 直线上点的横坐标分别表示p伪谱法、h伪谱法和hp伪谱法的配点位置), 从图中可以看出:①p方法为全局伪谱法, 配点分布具有两端密、中间疏的特点; ②h方法通过细化区间来提高精度, 每个区间内配点相同, 图中可以看出区间的划分; ③hp方法不具备统一的配点分布标准, 相比于h方法, 区间划分较少。

本文中所采用的hp自适应策略为:设定区间配点数上限, 当配点数未达到上限时, 通过增加区间内配点数(p法)提高精度, 当配点数达到上限仍未满足求解精度时, 再利用细化时间区间(h法)来提高精度。文献[12]提出另一种自适应策略:通过比较相对曲率r_k与曲率阈值r_max的大小关系来决定采取细分时间区间(h法)还是增加区间内的配点数(p法)来提高算法精度。表 3为这两种细化策略在不同的精度要求下的优化结果比较, 其中hp1表示本文细化策略, hp2为文献[12]中细化策略。

从表 3中可以看出, 精度要求较低时, hp2方法所需的配点总数小于hp1方法, 计算效率也高于hp1;精度要求提高到10^-5时, hp2方法逐渐失去了优势; 精度为10^-6时, hp1方法在计算效率上明显优于hp2方法。分析其可能的原因, 文献[12]中hp伪谱法在初始迭代中, NLP问题的解与原问题的解存在较大的误差, 此时根据曲率进行分段可能导致分段过多且分段点不准确, 且这种情况随着精度要求的提高表现得更为突出, 而本文中所采用的hp伪谱法优先采用p法来提高精度, 在一定程度上可以消除该问题对求解效率的影响。

5 结论

本文提出了一种基于hp自适应Radau伪谱法的hp自适应网格算法, 并将其应用于高超声速飞行器上升段的轨迹优化问题。仿真结果表明:①与全局伪谱法和固定阶伪谱法相比, hp伪谱法充分结合了p法的快速收敛性和p法的计算稀疏性, 具有更合理的配点分布, 能够以较少的时间代价获得更高精度的解; ②与以曲率为迭代判据的hp伪谱法相比, 本文所提出的伪谱法优先采用p方法, 在一定程度上能够避免分段过多且不准确的问题, 更适合求解高精度的最优控制问题, 具有一定的工程应用价值。