0 引　言

当前飞行任务复杂化、多样化的发展趋势对飞行器动力系统提出了新的需求，其中火箭基组合循环（Rocket Based Combined Cycle，RBCC）动力系统以其优越的性能受到了各国研究人员的关注。相对传统的火箭动力，RBCC将高推重比、低比冲的火箭发动机和低推重比、高比冲的吸气式冲压发动机组合在一起，具有可在宽来流范围内工作、实现多模态一体化结构设计的特点，尤其适合执行宽空域、宽速域飞行任务^[1]。

与传统运载火箭不同，RBCC可充分利用大气层中的氧气从而减少对氧化剂的需求，综合比冲是现有火箭发动机的数倍，因而是一种能够有效降低运载任务费用的动力装置。但由于RBCC动力发动机核心是吸气式发动机，其工作性能与飞行状态相互影响，且工作范围较宽，推力变化较大，这也使得RBCC动力运载器的运载能力与轨迹设计关系更为密切，且轨迹设计约束更为严格和复杂，传统火箭助推运载器的轨迹设计经验和规律已难以胜任，因而有必要针对RBCC动力系统开展相应的轨迹优化设计工作。

作为RBCC动力系统的核心，吸气式发动机工作性能与飞行状态的耦合是轨迹优化中的一大难点。针对这一问题，Lu等^[2]指出，吸气式动力运载器的上升轨迹优化问题可转化为典型两点边值问题求解。Pan等^[3]通过简化协态方程降低了间接法的求解难度，同时借助多重打靶法实现了初值的快速生成。李惠峰等^[4]基于有限差分方法和改进牛顿法，设计了一种以参考面积为同伦参数的迭代方法，实现了吸气式运载器上升轨迹的快速优化。Derek等^[5]在使用代理模型进行吸气式动力运载器轨迹设计时，证明了加速度设计对于燃料消耗的影响。

同时，由于RBCC组合式发动机对飞行状态比较敏感，有许多学者将关注的重点转移至上升段动压剖面的设计。同时，POST等标准程序在处理RBCC动力飞行器轨迹优化问题时可能会出现不合理状态点。闫晓东等^[6]研究了RBCC飞行器上升段等动压飞行轨迹设计方法，虽然方法具有较好的可行性，但同样未考虑上升段轨迹的最优性。孙佩华等^[7]总结分析了不同定动压下的上升轨迹特性从而设计出上升段动压曲线，但设计的动压曲线较为简单，优化程度有限。

此外，一些数值优化算法同样在RBCC动力飞行器的轨迹优化问题中得到了应用。龚春林等^[8]利用Radau伪谱法建立了最优轨迹求解模型，并针对上升段燃料消耗最小问题进行优化。李响等^[9]对控制量进行参数化处理，利用遗传算法进行寻优，对飞行器射程最大轨迹进行了设计。周宏宇等^[10]采用三次多项式设计组合动力发动机各动力段的攻角剖面，并利用改进粒子群算法进行求解。这些数值算法均取得了较好的优化效果，但优化效率均较低。

凸优化方法也在飞行器上升段轨迹设计中得到了一些应用。Szumk等^[11]以推力方向作为控制量，利用无损凸化技术求解运载火箭上升段轨迹优化问题，但该模型忽略了运载器所受升力。Liu等^[12]将推力方向和气动系数作为控制量，但考虑的气动力模型较为简单。王嘉炜等^[13]在对固体火箭助推飞行器轨迹优化问题进行建模时，选择攻角作为唯一控制量，有效处理了飞行器气动力的计算。

文中以RBCC动力高超声速飞行器为研究对象，基于凸优化理论建立优化模型，对复杂约束下飞行器上升段轨迹优化设计问题进行研究，并通过算例验证该方法的有效性。文中的研究旨在结合RBCC飞行器的特点建立优化设计方法、优化建模思路并分析结论，为这类新型飞行器的轨迹设计提供部分参考。

1 轨迹优化问题描述

RBCC动力高超声速飞行器采用一体化设计，由载机水平投放后由自身的RBCC发动机提供动力，经过多种工作模态将飞行器助推至预定的速度和高度。参考文献[14]根据RBCC工作特性的不同将工作模态划分为以下四个阶段：

（1）引射模态：爬升初段，飞行状态尚未达到冲压发动机开机条件，动力系统主要依靠小比冲火箭发动机所提供的大推力实现加速爬升。

（2）冲压模态：该模态主要由RBCC动力系统中的大比冲冲压发动机提供推力，是整个动力系统的核心工作模态。根据进入冲压发动机燃烧室的空气速度不同，该模态可分为亚燃冲压和超燃冲压两种。该模态下，发动机推力性能受燃油当量比、动压与攻角影响较大，是轨迹设计中的一大难点。

（3）冲压火箭模态：由于大气密度随飞行高度的增大会逐渐减小，冲压发动机进气量也会随之减小，进而导致推力大小无法满足任务需求，此时火箭发动机开机，与冲压发动机协同工作。

（4）纯火箭模态：当飞行高度超过某个临界点时，空气过于稀薄，导致冲压发动机无法正常工作而关机，此时仅由火箭发动机提供推力。

可见，与传统运载火箭的控制方式不同，RBCC动力系统工作模态多，动力系统与飞行状态耦合程度高；工作性能最佳的冲压模态对飞行状态极为敏感，工作条件极为苛刻，这使得RBCC动力高超声速飞行器上升段轨迹优化难度更大，需要建立合理的优化求解模型。

2 优化问题模型

2.1 动力学模型

飞行器上升段作为主动飞行段，通常选择发射坐标系作为该阶段飞行运动参考坐标系。若不考虑横向机动飞行，且认为姿态控制系统理想工作，其纵向射面内质心运动方程可表示为^[7]：

$ \left\{ \begin{gathered} \dot h = V\sin \theta \\ \dot L = \dfrac{{Rv\cos \theta }}{r} \\ \dot V = - g\sin \theta + \dfrac{{P\cos \alpha - X}}{m} \\ \dot \theta = \left( {\dfrac{V}{r} - \dfrac{g}{V}} \right)\cos \theta + \dfrac{{Y + P\sin \theta }}{{mV}} \\ \dot m = - {m_s} \\ \end{gathered} \right. $ (1)

式中：h、L、V、$ \theta $、m、$ \alpha $、g、r、R、$ {m_s} $、P、X、Y分别为飞行器纵向平面内飞行高度、航程、速度、弹道倾角、飞行器质量、攻角、重力加速度、地心距、地球半径、组合发动机燃料秒消耗量、组合发动机推力、气动阻力以及气动升力。

由于航程项不影响其余状态量的计算，由此，为简化优化模型、提高计算效率，将纵向射面内质心运动方程表示为：

$ \left\{ \begin{gathered} \dot h = V\sin \theta \\ \dot V = - g\sin \theta + \dfrac{{P\cos \alpha - X}}{m} \\ \dot \theta = \left( {\dfrac{V}{r} - \dfrac{g}{V}} \right)\cos \theta + \dfrac{{Y + P\sin \theta }}{{mV}} \\ \dot m = - {m_s} \\ \end{gathered} \right. $ (2)

基于美国1976年公布的标准大气模型，参考文献[15]中利用拟合方法给出高度范围0~91 km间大气参数的计算公式。同时，重力加速度大小计算公式如下：

$ g = \dfrac{{{R^2}}}{{{r^2}}}{g_0} $ (3)

式中：R表示地球半径；r为地心距； $ {g_0} $为海平面处重力加速度。

文中所研究飞行器气动数据可通过数值模拟或地面风洞实验得到，并以数据表的方式进行存储。其所受气动力为：

$ \left\{ \begin{gathered} X = qS{C_x} \\ Y = qS{C_y} \\ {C_x} = {f_x}\left( {Ma,\alpha ,h} \right) \\ {C_y} = {f_y}\left( {Ma,\alpha ,h} \right) \\ \end{gathered} \right. $ (4)

式中： $ q $为飞行动压； $ \;\rho $为大气密度； $ S $为飞行器参考面积； $ {C_y} $、 $ {C_x} $分别为飞行器气动升、阻力系数，可通过气动数据插值得到。

RBCC发动机推力与燃料秒流量可表示为：

$ \begin{gathered} P = {P_R} + {P_A} \\ {{\dot m}_s} = {{\dot m}_{sr}} + {{\dot m}_{sa}} \\ \end{gathered} $ (5)

式中： $ {P_R} $和 $ {P_A} $分别代表火箭发动机推力与冲压发动机推力，可分别由下式计算得到：

$ \begin{gathered} {P_R}{\text{ = }}{{\dot m}_{sr}}{I_{spr}} \\ {P_A} = {{\dot m}_{sa}}{E_r}{I_{spa}} \\ {{\dot m}_{sa}} = \rho SV\phi \\ \end{gathered} $ (6)

式中： $ {\dot m_{sr}} $、 $ {I_{spr}} $分别为火箭发动机秒流量、火箭发动机比冲； $ {\dot m_{sa}} $、 $ \rho $、 $ S $、 $ \phi $、 $ {E_r} $、 $ {I_{spa}} $分别为冲压发动机秒流量、当地大气密度、冲压发动机进气道横截面积、流量比、当量比、冲压发动机比冲。

2.2 控制量选择

在不考虑姿态控制的轨迹优化设计过程中，若直接选取攻角 $\alpha $作为控制量，优化结果中可能出现控制量变化率过大甚至一阶不连续情况，这与飞行器实际飞行情况不符^[16]。针对这一问题，文中引入攻角变化率 $\dot \alpha $作为伪控制量，攻角 $\alpha $可视为附加的状态量，通过 $\dot \alpha $积分得到，即

$ \alpha {\text{ = }}{\alpha _0} + \int {\dot \alpha {\rm{d}}t} $ (7)

同时，由于文中研究对象尚处于初步设计阶段，难以针对不同模态给出准确的发动机工作模态划分准则和约束要求，所以文中暂不给定模态划分准则，而是将火箭发动机燃料秒耗量作为全程控制量，进而通过轨迹优化设计获得其取值规律，最终获得RBCC动力系统的工作模态转换准则和机理。因此，对于文中的上升段轨迹优化设计问题，最终的控制量选择为 $\dot \alpha $、 $ {\dot m_{sr}} $。

同理，在考虑攻角控制系统二阶滞后的情况下则选定控制量为 ${\dot \alpha _c}$、 $ {\dot m_{sr}} $。

2.3 约束条件

在RBCC动力高超声速飞行器的上升段轨迹优化设计中，文中主要考虑以下约束条件对飞行轨迹的限制：

（1）过程约束

为保证上升段飞行安全性，在飞行过程中过载、动压、热流均不能超过飞行器临界值，否则将对飞行器结构造成损害。由于驻点位置是飞行器受热严重区域，为保证整体受热满足约束要求，主要对驻点热流密度加以限制。动压、过载、驻点热流密度计算公式与约束条件如下：

$\left\{ \begin{array}{l} q = \dfrac{1}{2}\rho {V^2} < {q_{\max }}\\ n = \dfrac{{\left| {X\sin \alpha + Y\cos \alpha } \right|}}{{m{g_0}}} < {n_{\max }}\\ \dot Q = K\sqrt \rho {V^{3.15}} < {{\dot Q}_{\max }} \end{array} \right. $ (8)

式中： $K$为常数，其取值与飞行器气动外形相关。

（2）状态变量约束

通常选取能够描述飞行器速度或位置状态的参数量作为状态变量，例如高度、速度等。从实际可行性以及安全飞行角度考虑，将存在状态变量约束。若用s表示状态变量，则有：

$ {s_{\min }} \leqslant {{s}} \leqslant {s_{\max }} $ (9)

（3）控制量约束

综合考虑飞行性能、控制系统性能和发动机性能，可建立控制变量约束模型。用u表示控制变量，则有：

$ {u_{\min }} \leqslant u \leqslant {u_{\max }} $ (10)

（4）端点约束

文中研究问题中端点约束可分为初始约束与终端约束两种。初始约束代表飞行器初始飞行状态并影响后续飞行能力；终端约束主要考虑末端高度、速度以及速度倾角等，是判断飞行器能否完成飞行任务的重要条件。用s表示状态变量，则有：

$ \left\{ \begin{gathered} {s_{0\min }} \leqslant {\text{s}}({t_0}) \leqslant {s_{0\max }} \\ {s_{f\min }} \leqslant {\text{s}}({t_f}) \leqslant {s_{f\max }} \\ \end{gathered} \right. $ (11)

2.4 优化目标

在助推级轨迹优化中，为挖掘助推级运载潜力文中将优化目标选取为上升段末端机械能最大，该优化目标可以表示为：

$ J = \left( {\dfrac{1}{2}{m_f}V_f^2 + {m_f}{g_f}{h_f}} \right) \to \max $ (12)

3 优化设计方法

3.1 凸优化基础理论

凸优化问题如下式：

$ \begin{split} &\min\;\; {f_0}\left( x \right) \\ & {\text{subject to }}{f_i}\left( x \right) \leqslant 0,i = 1, \cdots ,m \\ & \qquad \qquad a_i^Tx = {b_i},i = 1, \cdots ,p \end{split} $ (13)

式中：目标函数 $ {f_0}\left( x \right) $和不等式约束函数 $ {f_i}\left( x \right) $都是凸函数，而等式约束函数 $ {h_i}\left( x \right) = a_i^Tx - {b_i} $是仿射函数。理论上只要能够将一个优化问题描述为凸形式，即可在多项式时间（代数运算次数为问题维度的有限次多项式函数）内得到物理可行问题的全局最优解，且不依赖任何初值。显然，上升段轨迹优化设计问题是一个高度的非线性最优控制问题，无法直接采用凸优化方法求解，因此需要先对优化问题模型进行一些处理。

3.2 时间自由问题处理

由于上升段总飞行时间未知，为了考虑该类时间自由问题，定义新的自变量 $ \tau \in \left[ {0,1} \right] $和控制量 $ {u_t} = {t_f} - {t_0} $，并且将原问题的时间区间映射到[0,1]上，表示为：

$ t = {t_0} + \left( {{t_f} - {t_0}} \right)\tau ,\tau \in \left[ {0,1} \right] $ (14)

将原时间自变量作为新的状态变量，则有：

$ \dfrac{{{\rm{d}}t}}{{{\rm{d}}\tau }} = {u_t},\tau \in \left[ {0,1} \right] $ (15)

以姿态控制系统理想情况为例，设原弹道状态量为 $ x = \left[ {h,V,\theta ,m,\alpha } \right] $，弹道方程记为 $\dfrac{{{\rm{d}}x}}{{{\rm{d}}t}} = f\left( {x,u} \right)$，则以 $ \tau $为自变量的弹道运动方程表达为：

$ \dfrac{{{\rm{d}}x}}{{{\rm{d}}t}} = f\left( {x,u} \right)\dfrac{{{\rm{d}}t}}{{{\rm{d}}\tau }} = f\left( {x,u} \right){u_\tau } $ (16)

3.3 模型凸化处理

3.3.1 微分动力学方程约束处理

为了应用凸优化求解弹道优化问题，因此需要对问题进行凸化。凸化处理的基本思想是将非线性的系统方程和非线性的约束条件在给定的轨迹点处进行线性化处理。以不考虑姿态控制的上升段末端机械能最大问题为例，模型处理过程如下，设存在一个优化的初始状态和控制量序列：

$ \begin{gathered} {{\bar x}_{(k)}} = {[{h_{(k)}},{V_{(k)}},{\theta _{(k)}},{m_{(k)}},{\alpha _{(k)}},{\beta _{(k)}},{\alpha _{c(k)}},{t_{(k)}}]^{\rm{T}}} \\ {u_{(k)}} = {\left[ {{{\dot \alpha }_{c(k)}},{m_{sr(k)}}} \right]^{\rm{T}}} \\ \end{gathered} $ (17)

采用一阶欧拉法对状态和控制量序列下的线性状态约束进行离散处理，得到一个线性递推公式（凸约束）。设离散采用周期为T，则原连续状态约束成为：

$ \begin{gathered} \bar x\left( {k + 1} \right) = A\left( k \right)\bar x\left( k \right) + B\left( k \right)u\left( k \right) + C\left( k \right) \\ {\text{ }}k = 1,2, \cdots ,m \\ \end{gathered} $ (18)

其中，m为采样个数。各系数表达式为：

$\begin{split} & {A_{\left( k \right)}} = {\left[ \begin{gathered} \dfrac{{\partial {{\bar f}_1}}}{{\partial h}}{\text{ }}\dfrac{{\partial {{\bar f}_1}}}{{\partial V}}{\text{ }}\dfrac{{\partial {{\bar f}_1}}}{{\partial \theta }}{\text{ }}\dfrac{{\partial {{\bar f}_1}}}{{\partial m}}{\text{ }}\dfrac{{\partial {{\bar f}_1}}}{{\partial \alpha }} \\ \dfrac{{\partial {{\bar f}_2}}}{{\partial h}}{\text{ }}\dfrac{{\partial {{\bar f}_2}}}{{\partial V}}{\text{ }}\dfrac{{\partial {{\bar f}_2}}}{{\partial \theta }}{\text{ }}\dfrac{{\partial {{\bar f}_2}}}{{\partial m}}{\text{ }}\dfrac{{\partial {{\bar f}_2}}}{{\partial \alpha }} \\ \dfrac{{\partial {{\bar f}_3}}}{{\partial h}}{\text{ }}\dfrac{{\partial {{\bar f}_3}}}{{\partial V}}{\text{ }}\dfrac{{\partial {{\bar f}_3}}}{{\partial \theta }}{\text{ }}\dfrac{{\partial {{\bar f}_3}}}{{\partial m}}{\text{ }}\dfrac{{\partial {{\bar f}_3}}}{{\partial \alpha }}{\text{ }}0 \\ \dfrac{{\partial {{\bar f}_4}}}{{\partial h}}{\text{ }}\dfrac{{\partial {{\bar f}_4}}}{{\partial V}}{\text{ }}\dfrac{{\partial {{\bar f}_4}}}{{\partial \theta }}{\text{ }}\dfrac{{\partial {{\bar f}_4}}}{{\partial m}}{\text{ }}\dfrac{{\partial {{\bar f}_4}}}{{\partial \alpha }} \\ \dfrac{{\partial {{\bar f}_5}}}{{\partial h}}{\text{ }}\dfrac{{\partial {{\bar f}_5}}}{{\partial V}}{\text{ }}\dfrac{{\partial {{\bar f}_5}}}{{\partial \theta }}{\text{ }}\dfrac{{\partial {{\bar f}_5}}}{{\partial m}}{\text{ }}\dfrac{{\partial {{\bar f}_5}}}{{\partial \alpha }} \\ {\text{ }}0{\text{ }}0 \\ \end{gathered} \right]_{{x_k},{u_k}}} \text{} \\ &{B_{1(k)}} = {\left[ \begin{gathered} \dfrac{{\partial {{\bar f}_1}}}{{\partial \dot \alpha }}{\text{ }}\dfrac{{\partial {{\bar f}_1}}}{{\partial {m_{sr}}}} \\ \dfrac{{\partial {{\bar f}_2}}}{{\partial \dot \alpha }}{\text{ }}\dfrac{{\partial {{\bar f}_2}}}{{\partial {m_{sr}}}} \\ \dfrac{{\partial {{\bar f}_3}}}{{\partial \dot \alpha }}{\text{ }}\dfrac{{\partial {{\bar f}_3}}}{{\partial {m_{sr}}}} \\ \dfrac{{\partial {{\bar f}_4}}}{{\partial \dot \alpha }}{\text{ }}\dfrac{{\partial {{\bar f}_4}}}{{\partial {m_{sr}}}} \\ \dfrac{{\partial {{\bar f}_5}}}{{\partial \dot \alpha }}{\text{ }}\dfrac{{\partial {{\bar f}_5}}}{{\partial {m_{sr}}}} \\ {\text{ 0 0 }} \\ \end{gathered} \right]_{{x_k},{u_k}}} \text{，}{B_{2(k)}} = {\left[ \begin{gathered} \dfrac{{\partial {{\bar f}_1}}}{{\partial {u_t}}} \\ \dfrac{{\partial {{\bar f}_2}}}{{\partial {u_t}}} \\ \dfrac{{\partial {{\bar f}_3}}}{{\partial {u_t}}} \\ \dfrac{{\partial {{\bar f}_4}}}{{\partial {u_t}}} \\ \dfrac{{\partial {{\bar f}_5}}}{{\partial {u_t}}} \\ {\text{ 1 }} \\ \end{gathered} \right]_{{x_k},{u_k}}} \text{}\\ &{C_{\left( k \right)}} = \bar f - {A_{\left( k \right)}}{\bar x_{\left( k \right)}} - {B_{1\left( k \right)}}{\alpha _{\left( k \right)}} - {B_{2\left( k \right)}}u \end{split} $ (19)

为保证线性化的有效性，需引入信赖域约束：

$ \left| {x - {x^k}} \right| \leqslant \delta $ (20)

式中： $ \delta $即为优化问题的信赖域，可依据具体的飞行器运动模型设定。

3.3.2 过程约束处理

飞行过程中的动压、过载、驻点热流密度等过程约束同样需要在参考点处进行线性化处理。记 ${f_c} = \left[ {{f_6}{\text{ }}{f_7}{\text{ }}{f_8}} \right] = {\left[ {\dot Q - {{\dot Q}_{\max }}{\text{ }}P - {P_{\max }}{\text{ }}n - {n_{\max }}} \right]^{\rm{T} }}$，过程约束可表示为

$ {f_c}\left( {x,u} \right) \leqslant 0 $ (21)

对其进行线性化可得：

$ {C_a}\left( k \right)x\left( k \right) + {C_b}\left( k \right) \leqslant 0 $ (22)

式中： $ {C_a}\left( k \right) $、 $ {C_b}\left( k \right) $的求解方法与公式(19)中 $ {A_{\left( k \right)}} $、 $ {C_{\left( k \right)}} $的求解方法类似，此处不做赘述。

3.3.3 状态变量约束及控制变量约处理

控制变量约束、状态变量约束的离散形式为：

$ \begin{gathered} {u_{\min }} \leqslant u\left( n \right) \leqslant {u_{\max }} \\ {x_{\min }} \leqslant x\left( n \right) \leqslant {x_{\max }} \\ \end{gathered} $ (23)

状态变量的端点条件表达为离散形式为： $h\left( 1 \right) = {h_0} ， V\left( 1 \right) = {V_0}$， $ \theta \left( 1 \right) = {\theta _0}，m\left( 1 \right) = {m_0} $， $ \alpha \left( 1 \right) = {\alpha _0} $， $ t\left( 1 \right) = {t_0} $， $ h\left( m \right) = {h_f} $， $ V\left( m \right) \geqslant {V_{f\min }} $， $m\left( m \right) \geqslant {m_{f\min }} $。

控制变量、状态变量约束及端点约束的离散形式已为凸约束，无需额外处理。

3.3.4 优化指标凸化处理

机械能公式为：

$ {f_9} = \dfrac{1}{2}m{V^2} + mgh $ (24)

对其进行线性化可得：

$ {f_9} = {C_A}\left( k \right)x\left( k \right) + {C_B}\left( k \right) $ (25)

其中，

$ \begin{split} &{C_{A(k)}} = {\left[ {\dfrac{{\partial {f_9}}}{{\partial h}}{\text{ }}\dfrac{{\partial {f_9}}}{{\partial V}}{\text{ }}\dfrac{{\partial {f_9}}}{{\partial \theta }}{\text{ }}\dfrac{{\partial {f_9}}}{{\partial m}}{\text{ }}\dfrac{{\partial {f_9}}}{{\partial \alpha }}{\text{ }}0} \right]_{{x_k},{u_k}}} \\ & {C_{B(k)}} = {f_9} - {C_{A(k)}}x\left( k \right) \end{split} $ (26)

综上所述，原轨迹优化问题可表示为：

$ \begin{gathered} \min {\text{ }}(25) \\ {\text{subject to }}(18),{\text{ }}(20),{\text{ }}(21),{\text{ }}(22),{\text{ }}(23) \\ \end{gathered} $ (27)

3.4 轨迹优化求解算法

总结上文，基于凸优化方法的轨迹优化求解步骤如下：

(1) 设k=0时， $ \left\{ {{{\bar x}_{\left( 0 \right)}};{u_{\left( 0 \right)}}} \right\} $为一条初始弹道；

(2) 在第k+1次迭代中，利用式参考弹道的相关参数，求解凸优化问题 $ P_0^{\left( {k + 1} \right)} $，得到序列最优解 $ \left\{ {{{\bar x}_{\left( {k + 1} \right)}};{u_{\left( {k + 1} \right)}}} \right\} $；

(3) 判断收敛判据是否满足：

$ \left| {{{\bar x}_{\left( {k + 1} \right)}}\left( n \right) - {{\bar x}_{\left( k \right)}}\left( n \right)} \right| \leqslant \varepsilon ,n = 1,2, \cdots ,m $ (28)

其中，收敛域 $ \varepsilon $可以根据具体问题设定。若判据满足则进入第(4)步，否则令k=k+1，判断回溯条件并进入第(2)步；

(4) $ \left\{ {{{\bar x}_{\left( {k + 1} \right)}};{u_{\left( {k + 1} \right)}}} \right\} $即是全局最优状态轨迹和最优控制变量。

4 算例仿真及分析

4.1 仿真条件

RBCC高超声速飞行器初始质量4 t，其中RBCC助推级共携带燃料3.2 t。飞行器由载机水平投放，初始飞行速度为388.9 m/s，飞行高度为10.5 km，攻角为3°。飞行全程过载不超过5 g, 动压不超过200 kPa，驻点热流密度不超过800 kW/m²。飞行过程中攻角应保持在−2°~10°之间，且攻角变化率的绝对值不超过0.15 (°)/s。

综合考虑计算效率与解的精度，文中仿真过程中共取61个离散点，信赖域设置为[20 km，1000 m/s，10°，200 kg，3°，20 s]，收敛域设置为[0.1 km，10 m/s，0.1°，0.1 kg，0.1°,0.1 s]。

4.2 上升段末端机械能最大算例

针对该算例，凸优化算法经过12轮迭代后收敛至最终优化结果，每轮迭代约耗时1.8 s，优化过程共耗时19.88 s，优化结果如图1~图3所示。其中，上升段飞行器共飞行135.66 s，终端质量为800 kg，终端高度为62.08 km、飞行速度为4866.99 km/s。最大飞行过程中动压、过载、驻点热流密度均满足给定的过程约束限制。

图 1. 状态量优化结果曲线

Fig. 1. State quantity optimization result curve

下载图片 查看所有图片

图 3. 控制量优化结果曲线

Fig. 3. Control quantity optimization result curve

下载图片 查看所有图片

图3中的虚线为凸优化所使用的初始猜测轨迹，这里使用的是初、末端状态猜测值的连线。显然，该初始猜测轨迹无法满足动力学方程约束。可以看出，在初始轨迹精度较低且不满足动力学方程约束的条件下，凸优化方法仍能有效、快速地收敛至优化结果。

图 2. 过程约束优化结果曲线

Fig. 2. Result curve of process constrained optimization

下载图片 查看所有图片

对于该问题而言，飞行器上升段飞行过程主要可分为三段。第一段为0~22.6 s，此时飞行器飞行马赫数较低，因此冲压发动机工作效率低，该段RBCC发动机主要以火箭模态工作，飞行高度变化不大，飞行马赫数逐渐增大；第二段为22.6~91 s，该段飞行高度、马赫数都较适宜冲压发动机工作，此时火箭发动机关机，RBCC发动机以冲压模态工作，飞行马赫数变化不大，飞行高度爬升至约30 km，燃料消耗显著降低；第三段为91 s至上升段结束，该段由于高度较高，冲压发动机效率再次降低，RBCC发动机主要以火箭模态工作，飞行器飞行高度和马赫数快速增加，并达到终端值。

5 结　论

文中开展了RBCC动力高超声速飞行器上升段轨迹优化设计研究，分别针对不考虑姿态控制及考虑控制系统二阶滞后的运动模型下末端机械能最大算例进行了仿真分析，结论如下：

（1）文中提出的轨迹优化方法可有效处理复杂工作模态下RBCC助推飞行器上升段轨迹优化问题，为未来关于这一类轨迹设计与优化的工作提供了一些新的思路。

（2） RBCC发动机火箭模态峰-谷-峰的工作过程能够最大限度的结合并发挥火箭发动机工作范围广和冲压发动机燃料消耗少的优势，从而使飞行器在燃料更省的情况下获得更大的终端高度、马赫数。