可重复使用运载器(RLV)是指能够自由往返于地球表面与空间轨道之间且可重复使用的多用途飞行器.其再入过程往往伴有参数不确定和外界干扰^[1-2]等问题.为了保证飞行器安全稳定地再入飞行，对于再入段轨迹的优化与制导律的设计尤为关键. RLV再入轨迹优化的目标是生成满足各种路径约束(动压、热流密度和过载约束)、状态和控制量约束以及边值约束条件并实现某个最优目标的飞行轨迹，引导飞行器从再入起始点安全到达要求的终端区域内.近年来求解轨迹优化问题运用较为广泛的是伪谱法^[3]，如Gauss伪谱法、Radua伪谱法等.该方法将最优控制问题离散化为非线性规划(nonlinear programming，NLP)问题，并通过数值方法求解NLP问题来获得最优轨迹.然而伪谱法由于计算量较大、求解时间不确定，实时性难以保证^[4]，难以应用到RLV的轨迹在线求解中.

与求解NLP问题相比，利用凸优化求解RLV再入轨迹优化问题在求解速度上具有较大的优势，且具备全局收敛性，使得其在RLV轨迹在线求解上的应用成为可能.文献[5]将凸优化方法中的二阶锥规划方法应用于多约束再入轨迹优化问题，并对其全局收敛性进行了证明.然而文章中对飞行速度采用近似化的能量函数进行逼近，轨迹求解的精确性可能因此受到影响.在文献[5]的基础上，文献[6]提出将序列凸规划(sequential convex programming，SCP)方法用于求解带有约束的行星再入轨迹优化问题，该方法具备良好的收敛速度，具有用于实时轨迹规划的潜力.序列凸规划化的主要思想是通过求解序列近似凸子问题，实现子问题的解向原问题的收敛.文献[6]提出的算法使用固定的信赖域约束，然而固定的信赖域可能导致算法失去提升收敛性能的机会，算法收敛速度因此受到约束.得益于其全局收敛性及优秀的收敛性能，序列凸规划方法不似伪谱法那般依赖于初始状态.但序列凸规划求解过程中仍存在初值猜测问题，若初始猜测轨迹与最优轨迹之间偏差较大，则需要较多的迭代次数，算法求解速度也相应下降.此外文中讨论的是终端时间固定的情况，对不确定的终端时间未做分析.文献[7]提出了一种信赖域更新策略的序列凸规划方法，设计了高精度返回快速轨迹优化算法并应用于火箭返回着陆问题中，进一步提高了算法的快速性，然而文章未给出具体的信赖域更新策略.本文提出了一种基于变信赖域策略的序列凸规划算法，在传统的序列凸规划算法的基础上，利用性能指标函数作为每次迭代后的判定条件，设计了信赖域更新策略，提高了算法的收敛性能.此外在序列凸规划求解过程中，利用预测校正算法解决序列凸规划的初值猜测问题，加快算法的收敛速度，与此同时获得了再入终端时间，解决了终端时间不确定问题.本文将所提出算法与SNOPT优化包下的Gauss伪谱法和传统的序列凸规划进行轨迹优化仿真对比，表明所提出算法在收敛性能和计算速度方面有着显著的提升.

考虑实际再入过程中，由于外界干扰导致飞行轨迹大幅度偏离参考轨迹以及在某些紧急情况下要求RLV临时更换着陆场时，飞行器可能无法继续按照原轨迹飞行，此时需要通过轨迹重构提供可行轨迹.文献[8]指出，突发情况下的轨迹重构需要解决两大问题：一是要有快速轨迹优化算法保证可以在线计算重构轨迹；二是需要量化扰动对模型和飞行器约束的影响.文献[9]提出实时轨迹重构策略以解决着陆点变更问题，并采用滚动时域策略以抑制由轨迹求解消耗时间引起的状态量跳变.然而文章由于采用实时轨迹重构替代了制导律的作用，因此轨迹重构是无条件的、不间断的，对于机载计算机负载过大.本文针对RLV再入过程中遭遇的突发事件，在线对轨迹进行重构，并考虑轨迹重构耗时，预测实际轨迹重构初始点，以抑制可能产生的状态量跳变.在重构轨迹求解完成后，为了保证实时性的要求，本文采用基于LQR的制导律对重构轨迹进行跟踪.文中通过仿真将所提出方法与近年来发展迅速的预测校正制导进行对比，表明本文提出的轨迹重构策略在保证安全性、鲁棒性的同时，具备良好的实时性.本文首先对RLV再入轨迹优化问题进行描述，并提出变信赖域序列凸规划轨迹快速求解方法，然后针对RLV再入返回中遭遇的突发事件进行轨迹重构并进行跟踪制导，最后给出仿真结果以及结论.

1 再入轨迹优化问题描述 1.1 动力学模型

在RLV再入段，假设飞行器为无动力飞行的质点，考虑地球为旋转椭球时，忽略再入过程中侧力以及地球自转的影响，并取侧滑角为零. RLV再入三自由度运动方程^[9]为

$ \left\{ \begin{array}{l} \dot r = v\sin \gamma ,\\ \dot \theta = \frac{{v\cos \gamma \sin \chi }}{{r\cos \varphi }},\\ \dot \varphi = \frac{{v\cos \gamma \cos \chi }}{r},\\ \dot v = - \frac{D}{m} - g\sin \gamma ,\\ \dot \gamma = \frac{1}{v}\left[ {\frac{{L\cos \sigma }}{m} - \left( {g - \frac{{{v^2}}}{r}} \right)\cos \gamma } \right],\\ \dot \chi = \frac{{L\sin \sigma }}{{mv\cos \gamma }} + \frac{v}{r}\cos \gamma \sin \chi \tan \varphi . \end{array} \right. $

(1)

式中：r，θ，φ，v，γ，χ分别为地心距、经度、纬度、飞行速度、航迹角和航向角；σ为倾侧角；m为飞行器质量；g为重力加速度，g=μ_g/r²，其中μ_g为引力参数；L为升力，L=q_dSC_L；D为阻力，D=q_dSC_D，其中S为RLV气动参考面积，q_d为动压，q_d=0.5ρv²；ρ为大气密度，$\rho=\rho_{0} \mathrm{e}^{-\beta\left(r-R_{\mathrm{e}}\right) / R_{\mathrm{e}}}$，其中ρ₀为海平面处的大气密度，R_e为地球半径，β为常值系数；升力系数C_L和阻力系数C_D表示为攻角α的函数，将在仿真中给出.

1.2 约束条件

在RLV再入过程中，为了保证安全稳定飞行，飞行器需要严格满足一些约束条件，主要包括边值约束、路径约束和状态量约束.边值约束规定了飞行器状态量$\boldsymbol{x}=[r, \theta, \varphi, v, r, \chi]^{T}$，在再入中起点与终点处的取值，定义再入起点约束和再入终点约束分别为

$ \mathit{\boldsymbol{x}}\left( {{t_0}} \right) = {\mathit{\boldsymbol{x}}_0}, $

(2)

$ \mathit{\boldsymbol{x}}\left( {{t_{\rm{f}}}} \right) = {\mathit{\boldsymbol{x}}_{\rm{f}}}. $

(3)

再入过程常见的路径约束包括热流密度约束、动压约束和过载约束, 计算公式分别为

$ \dot Q = {c_1}\left( \mathit{\boldsymbol{x}} \right) = {k_{\rm{Q}}}{\rho ^{0.5}}{V^{3.15}} \le {{\dot Q}_{\max }}, $

(4)

$ q = {c_2}\left( \mathit{\boldsymbol{x}} \right) = \frac{1}{2}\rho {v^2} \le {q_{\max }}, $

(5)

$ n = {c_3}\left( \mathit{\boldsymbol{x}} \right) = \sqrt {{L^2} + {D^2}} \le {n_{\max }}. $

(6)

式中：k_Q为与飞行器相关的常数；$\dot{Q}_{\max }$，q_max，n_max分别为最大允许的热流、动压和过载.

此外，受飞行器性能影响，在再入过程中，控制量u=[α, σ]^T和状态量x满足的约束为

$ \left\{ \begin{array}{l} {\mathit{\boldsymbol{u}}_{\min }} \le \mathit{\boldsymbol{u}} \le {\mathit{\boldsymbol{u}}_{\max }},\\ {\mathit{\boldsymbol{x}}_{\min }} \le \mathit{\boldsymbol{x}} \le {\mathit{\boldsymbol{x}}_{\max }}. \end{array} \right. $

(7)

1.3 优化问题描述

综合上述动力学模型和约束条件，考虑以终端状态量和再入过程热流、动压、过载积分函数为指标的目标函数为

$ {J_1} = \phi \left[ {\mathit{\boldsymbol{x}}\left( {{t_{\rm{f}}}} \right)} \right] + \int_0^{{t_{\rm{f}}}} {G\left( \mathit{\boldsymbol{x}} \right){\rm{d}}t} . $

(8)

式中:t_f为再入飞行总时间，ϕ[x(t_f)]为与终端状态量相关的函数，G(x)为与再入过程热流、动压、过载相关的函数. RLV再入轨迹优化问题P0可描述为：满足条件式(1)~(7)，求解目标函数J₁最小.

2 变信赖域序列凸规划轨迹求解 2.1 问题的凸化

经过上一节的描述，针对问题P0求解可以得到再入最优轨迹.但是问题P0是一个高度非线性的优化问题，其中动力学模型(1)、路径约束式(4)~(6)、以及目标函数式(8)都是非线性、非凸的^[5]，为了使用凸优化方法对问题进行求解，需要对问题P0中的非凸约束进行凸化处理.

2.1.1 控制量的重新选取

在凸化处理过程中，若仅针对现有控制量进行处理，会引入高频抖振并对问题的收敛性产生影响^[5]，因此需要引入新的控制变量.动力学模型(1)中的控制量为攻角α和倾侧角σ, 其中攻角α由攻角-马赫数剖面确定，剩下需要设计的唯一控制量为倾侧角σ.引入新的辅助控制变量，从而实现控制量从状态量中解耦，令新的控制量为

$ u' = \dot \sigma , $

(9)

则运动方程式(1)可改写为

$ \mathit{\boldsymbol{\dot x'}} = \mathit{\boldsymbol{f}}\left( {\mathit{\boldsymbol{x'}}} \right) + \mathit{\boldsymbol{B}}u'. $

(10)

式中：$\boldsymbol{x}^{\prime}=[r ; \theta ; \varphi ; v ; \gamma ; \mathcal{X} ; \sigma]^{\mathrm{T}}$为新的状态量，u′= $\dot{\boldsymbol{\sigma}}$为新的控制量，状态矩阵f(x′)和控制矩阵B分别为

$ \mathit{\boldsymbol{f}}\left( {\mathit{\boldsymbol{x'}}} \right) = \left[ {\begin{array}{*{20}{c}} {v\sin \gamma }\\ {\frac{{v\cos \gamma \sin \chi }}{{r\cos \varphi }}}\\ {\frac{{v\cos \gamma \cos \chi }}{r}}\\ { - \frac{D}{m} - g\sin \gamma }\\ {\frac{1}{v}\left[ {\frac{{L\cos \sigma }}{m} - \left( {g - \frac{{{v^2}}}{r}} \right)\cos \gamma } \right]}\\ {\frac{{L\sin \sigma }}{{mv\cos \gamma }} + \frac{v}{r}\cos \gamma \sin \chi \tan \varphi }\\ 0 \end{array}} \right], $

(11)

$ \mathit{\boldsymbol{B}} = {\left[ {\begin{array}{*{20}{c}} 0&0&0&0&0&0&1 \end{array}} \right]^{\rm{T}}}. $

(12)

从式(12)中可以看到控制矩阵B是与状态量x′无关的一个常数矩阵，如此便将u′从x′中完成了解耦.此外，对于新的控制量σ·的约束，有助于加强倾侧角σ的平滑性.经过上述变换，轨迹优化问题P0转变为新控制量下的优化问题P1：满足条件式(2)~(7)、式(10)，求解目标函数J₁最小.

值得注意的是，对于问题P1若存在最优解，那么该解是问题P0的一个可行解.这是由于问题P1是在问题P0上引入了一个新的辅助控制变量，而问题P0中对于倾侧角σ的控制量约束转化为问题P1中的状态量约束，问题P0的其他状态量约束也包含在问题P1中，由此得到问题P1中的最优解满足问题P0的所有约束，进而得出该解是问题P0的一个可行解.

2.1.2 线性化处理

为了使问题能够利用凸优化求解，对问题P1进行凸化处理，对P1中的非凸约束进行线性化处理以将其转换为线性的凸约束.对轨迹优化问题P1中的非线性约束，包括运动方程式(10)、路径约束式(4)~(6)以及目标函数式(8)，基于小扰动线性化理论进行线性化处理，在给定状态点x′^*(t)处利用一阶泰勒展开式对方程进行逼近，有：

$ \mathit{\boldsymbol{\dot x'}} \approx \mathit{\boldsymbol{f}}\left( {{{\mathit{\boldsymbol{x'}}}^ * }} \right) + \mathit{\boldsymbol{A}}\left( {{{\mathit{\boldsymbol{x'}}}^ * }} \right)\left( {\mathit{\boldsymbol{x'}} - {{\mathit{\boldsymbol{x'}}}^ * }} \right) + \mathit{\boldsymbol{B}}u', $

(13)

$ c\left( {\mathit{\boldsymbol{x'}}} \right) \approx c\left( {{{\mathit{\boldsymbol{x'}}}^ * }} \right) + \frac{{\partial c\left( {\mathit{\boldsymbol{x'}}} \right)}}{{\partial \mathit{\boldsymbol{x'}}}}\left| {_{\mathit{\boldsymbol{x'}} = {{\mathit{\boldsymbol{x'}}}^ * }}} \right.\left[ {\mathit{\boldsymbol{x'}} - {{\mathit{\boldsymbol{x'}}}^ * }} \right], $

(14)

$ \begin{array}{l} {J_2} \approx \phi \left[ {\mathit{\boldsymbol{x'}}\left( {{t_{\rm{f}}}} \right)} \right] + \int_0^{{t_{\rm{f}}}} {\left[ {G\left( {{{\mathit{\boldsymbol{x'}}}^ * }} \right) + } \right.} \\ \;\;\;\;\;\;\;\left. {\frac{{\partial G\left( {\mathit{\boldsymbol{x'}}} \right)}}{{\partial \mathit{\boldsymbol{x'}}}}\left| {_{\mathit{\boldsymbol{x'}} = {{\mathit{\boldsymbol{x'}}}^ * }}} \right.\left( {\mathit{\boldsymbol{x'}} - {{\mathit{\boldsymbol{x'}}}^ * }} \right)} \right]{\rm{d}}t. \end{array} $

(15)

式中A(x′^*)为f(x′)在x′^*(t)处的雅克比矩阵，$\boldsymbol{c}\left(\boldsymbol{x}^{\prime}\right)=\left[\boldsymbol{c}_{1}\left(x^{\prime}\right), \boldsymbol{c}_{2}\left(\boldsymbol{x}^{\prime}\right), \boldsymbol{c}_{3}\left(\boldsymbol{x}^{\prime}\right)\right]^{\mathrm{T}}$.

基于泰勒展开式的特点，只有优化变量在参考点附近取值时，线性化的运动方程式和约束式才是对原非线性问题的良好近似.因此为了尽可能地减小逼近误差，保证线性化约束合理逼近原约束，引入信赖域约束：

$ \left| {\mathit{\boldsymbol{x'}} - {{\mathit{\boldsymbol{x'}}}^ * }} \right| \le \mathit{\boldsymbol{\varepsilon }}, $

(16)

式中ε为信赖域的半径.

经过上述线性化处理，新控制量下的轨迹优化问题P1可转化为凸优化问题P2：满足条件式(2)、(3)、(7)、(13)、(14)、(16)，求解目标函数J₂最小.

2.2 离散化处理

由于凸优化问题P2描述的仍是一个连续的最优控制问题，为了求得问题P2的数值解，需要首先对问题进行离散化处理.本文中采用梯形离散化方法^[10]，将时域均匀划分为N段，产生N+1个节点，每一段的步长为Δt=(t_f－t₀)/N，终端时间t_f在后续序列凸规划初值猜测时由预测校正算法给出.离散点表示为{t₀, t₁, t₂, …, t_N}，其中t_i=t_i－1+Δt, i=1, 2, …, N.对应的状态量和控制量表示为{x′₀, x′₁, x′₂, …, x′_N}，{u′₀, u′₁, u′₂, …, u′_N}.则离散化之后的运动方程式(13)，路径约束式(14)，目标函数式(15)可分别表示为

$ \begin{array}{l} {\mathit{\boldsymbol{h}}_i} = \mathit{\boldsymbol{x}}_i^\prime - \mathit{\boldsymbol{x}}_{i - 1}^\prime - \frac{{\Delta t}}{2}\left[ {\left( {\mathit{\boldsymbol{A}}_{i - 1}^{k - 1}\mathit{\boldsymbol{x}}_{i - 1}^\prime + \mathit{\boldsymbol{B}}u_{i - 1}^\prime + \mathit{\boldsymbol{f}}_{i - 1}^{k - 1} - } \right.} \right.\\ \;\;\;\;\;\;\left. {\mathit{\boldsymbol{A}}_{i - 1}^{k - 1}\mathit{\boldsymbol{x}}_{i - 1}^{\prime k - 1}} \right) + \left( {\mathit{\boldsymbol{A}}_i^{k - 1}\mathit{\boldsymbol{x}}_i^\prime + \mathit{\boldsymbol{B}}u_i^\prime + \mathit{\boldsymbol{f}}_i^{k - 1} - } \right.\\ \;\;\;\;\;\;\left. {\left. {\mathit{\boldsymbol{A}}_i^{k - 1}\mathit{\boldsymbol{x}}_i^{\prime k - 1}} \right)} \right] = 0, \end{array} $

(17)

$ \begin{array}{l} {\mathit{\boldsymbol{C}}_i} = \mathit{\boldsymbol{c}}\left( {\mathit{\boldsymbol{x}}_i^{\prime k - 1}} \right) + {\mathit{\boldsymbol{c}}^\prime }\left( {\mathit{\boldsymbol{x}}_i^{\prime k - 1}} \right)\left( {\mathit{\boldsymbol{x}}_i^\prime - \mathit{\boldsymbol{x}}_i^{\prime k - 1}} \right) - \left[ {{{\dot Q}_{\max }},} \right.\\ \;\;\;\;\;\;\;{\left. {{q_{\max }},{n_{\max }}} \right]^{\rm{T}}} \le 0, \end{array} $

(18)

$ \begin{array}{l} {J_3} = \phi \left[ {{\mathit{\boldsymbol{x}}^\prime }\left( {{t_{\rm{f}}}} \right)} \right] + \sum\limits_{i = 1}^N {\left[ {G\left( {\mathit{\boldsymbol{x'}}_i^{k - 1}} \right) + } \right.} \\ \;\;\;\;\;\;\;\left. {G\left( {\mathit{\boldsymbol{x}}_i^{\prime k - 1}} \right)\left( {\mathit{\boldsymbol{x}}_i^\prime - \mathit{\boldsymbol{x}}_i^{\prime k - 1}} \right)} \right]\Delta t. \end{array} $

(19)

经过离散化处理，凸优化问题P2的最优解可通过求解离散序列凸规划问题P3得到，问题P3可表示为：满足条件式(2)、(3)、(7)、(16)~(18)，求解目标函数J₃最小. P3的收敛性证明在文献[6]中可以找到，本文不赘述.

2.3 变信赖域序列凸规划 2.3.1 初始轨迹求解

针对离散化后的凸优化问题P3，采用序列凸规划算法对其进行求解.为了得到满足线性化近似需求的泰勒展开点，对凸优化问题进行迭代求解，迭代过程中参考点选取为上一次迭代求得的最优解，直到满足收敛约束.在传统的序列凸规划算法中，初始轨迹是基于给定的控制量常值积分得到的^[6]，如此便存在初值猜测问题.因此本文采用预测校正算法^[11]对初始控制量进行求解，引入能量变量$e=\frac{\mu_{\mathrm{g}}}{r}-\frac{V^{2}}{2}$代替时间作为自变量，由式(1)可得

$ \dot e = \frac{{DV}}{m}. $

(20)

考虑终端能量$e_{\mathrm{f}}=\frac{\mu_{\mathrm{g}}}{r_{\mathrm{f}}}-\frac{V_{\mathrm{f}}^{2}}{2}$，其中r_f、V_f分别为终端高度和终端速度.故可将再入终端约束转换为关于航程的单一约束：

$ \left| {z\left( {{\sigma _0}} \right)} \right| = \left| {s\left( {{e_{\rm{f}}}} \right) - {s_{\rm{f}}}} \right| = 0. $

(21)

式中待飞航程s(e)=arccos[sin φsin φ_f+cos φcos φ_fcos(θ_f－θ)].由式(1)、(20)可得到以能量为自变量，关于航程、高度和航迹角的运动方程为

$ \left\{ \begin{array}{l} \frac{{{\rm{d}}s}}{{{\rm{d}}e}} = \frac{{ - m\cos \gamma }}{{rD}},\\ \frac{{{\rm{d}}r}}{{{\rm{d}}e}} = \frac{{m\sin \gamma }}{D},\\ \frac{{{\rm{d}}\gamma }}{{{\rm{d}}e}} = \frac{{m\left[ {L\cos \sigma - \left( {g - \frac{{{V^2}}}{r}} \right)\cos \gamma } \right]}}{{D{V^2}}}. \end{array} \right. $

(22)

在给定的攻角剖面下，由上一次迭代的倾侧角对式(22)积分得到预测航程差，基于预测航程差利用Gauss-Newton法式(23)对倾侧角进行迭代处理，直到预测航程差满足终端约束，从而得到最终可行的倾侧角幅值.并在侧向上则采用侧向反转逻辑以缩小航向角误差，从而确定倾侧角的符号.最终得到满足再入要求的倾侧角指令：

$ \sigma _0^{\left( {k + 1} \right)} = \sigma _0^{\left( k \right)} - \frac{{z\left( {\sigma _0^{\left( k \right)}} \right)\left( {\sigma _0^{\left( k \right)} - \sigma _0^{\left( {k - 1} \right)}} \right)}}{{z\left( {\sigma _0^{\left( k \right)}} \right) - z\left( {\sigma _0^{\left( {k - 1} \right)}} \right)}}. $

(23)

2.3.2 变信赖域策略

得到初始猜测轨迹后，在序列凸规划后续迭代过程中，考虑式(16)中的信赖域，信赖域半径ε的大小决定了序列凸规划的收敛性能：若ε过大，则问题P3可能大幅度偏离原问题从而难以收敛；若ε过小，则迭代步长受到限制从而导致算法收敛速度不够快.因此，为了提高序列凸规划的收敛性能，在传统的序列凸规划的基础上提出了变信赖域策略：在每次迭代求解后，对比相邻迭代的实际性能指标函数ψ′与预测性能指标函数ψ，并基于对比结果设计信赖域更新策略.实际性能与预测性能指标函数分别表示为

$ {\psi ^\prime } = {J^\prime } + {\mu _1}\sum\limits_{i = 1}^{N - 1} {\left\| {\mathit{\boldsymbol{h}}_i^\prime } \right\|} + {\mu _2}\sum\limits_{i = 1}^N {\left\| {\mathit{\boldsymbol{C}}_i^\prime } \right\|} , $

(24)

$ \psi = J + {\mu _1}\sum\limits_{i = 1}^{N - 1} {\left\| {{\mathit{\boldsymbol{h}}_i}} \right\|} + {\mu _2}\sum\limits_{i = 1}^N {\left\| {{\mathit{\boldsymbol{C}}_i}} \right\|} . $

(25)

式中μ₁、μ₂分别为违反运动方程约束和违反路径约束的惩罚因子. J′、h_i′、C_i′分别表示为

$ {J^\prime } = \varphi \left[ {{\mathit{\boldsymbol{x}}^\prime }\left( {{t_{\rm{f}}}} \right)} \right] + \sum\limits_{i = 1}^N G \left( {\mathit{\boldsymbol{x}}_i^{\prime k - 1}} \right)\Delta t, $

(26)

$ \mathit{\boldsymbol{h}}_i^\prime = \mathit{\boldsymbol{x}}_i^{' k} - \mathit{\boldsymbol{x}}{_{i - 1}' ^k} - \frac{{\Delta t}}{2}\left[ {\left( {\mathit{\boldsymbol{B}}u_{i - 1}^\prime + \mathit{\boldsymbol{f}}_{i - 1}^k} \right) + \left( {\mathit{\boldsymbol{Bu}}_i^{' k} + \mathit{\boldsymbol{f}}_i^k} \right)} \right], $

(27)

$ \mathit{\boldsymbol{C}}_i^\prime = \mathit{\boldsymbol{c}}\left( {\mathit{\boldsymbol{x}}_i^{\prime k}} \right) - {\left[ {{{\dot Q}_{\max }},{q_{\max }},{n_{\max }}} \right]^{\rm{T}}}. $

(28)

实际性能指标函数给出了实际的离散点处的性能指标，相邻迭代的差值$\Delta \psi^{\prime}=| \psi^{\prime}\left(\boldsymbol{x}^{\prime k}, u^{\prime k}\right)-$${{\psi }^{\prime }}\left( {{\bf{x}}^{\prime k-1}}, {{\bf{u}}^{\prime k-1}} \right)|$则给出了序列凸规划算法中每一次迭代对于轨迹性能的提升，将之与相邻迭代间基于线性化的预测性能指标函数差值$\Delta \psi =|\psi \left( {{\bf{x}}^{\prime k}} \right., $$\left.u^{\prime k}\right)-\psi\left(\boldsymbol{x}^{\prime k-1}, u^{\prime k-1}\right) |$进行对比，可以判断出当前信赖域是否有利于收敛性能的提升.

2.3.3 变信赖域序列凸规划求解

综合上述初始轨迹求解与信赖域更新策略，得到基于变信赖域序列凸规划的最优轨迹求解步骤如下，整体轨迹优化流程如图 1所示.

步骤1 令k=0，k代表序列凸规划迭代次数，设置初始状态x′(t₀)=x′₀，基于预测校正算法得到初始轨迹x⁰.

步骤2 对于k≥1，在第k次迭代中对离散序列凸规划问题P3利用前一次的轨迹x^k－1求解得到[x′^k, u′^k].

步骤3 检查序列收敛条件sup|x′^k－x′^k－1|≤δ是否满足，其中k≥2, δ为迭代收敛阈值.若条件满足，则转至步骤5，否则转至步骤4.

步骤4 计算预测性能指标函数ψ，和实际性能指标函数ψ′，并对$\Delta \psi=| \psi\left(\boldsymbol{x}^{\prime k}, u^{\prime k}\right)-\psi\left(\boldsymbol{x}^{\prime k-1}\right.$u'^k-1和$\Delta \psi^{\prime}=\left|\psi^{\prime}\left(\boldsymbol{x}^{\prime k}, \boldsymbol{u}^{\prime k}\right)-\boldsymbol{\psi}^{\prime}\left(\boldsymbol{x}^{\prime k-1}, \boldsymbol{u}^{\prime k-1}\right)\right|$进行对比.若Δψ≥ξΔψ′, ξ为给定的系数，则当次序列迭代过程中性能指标提升幅度相对较小，可以适当放大信赖域ε=β₁ε，以寻找更适合的收敛步长.反之，则缩小信赖域$\mathit{\boldsymbol{\varepsilon }} = {\mathit{\beta }_2}\mathit{\boldsymbol{\varepsilon }}\left( {0 < {\mathit{\beta }_{\rm{2}}}{\rm{ < 1 < }}{\mathit{\beta }_{\rm{1}}}} \right).$然后令k=k+1，转至步骤2.

步骤5 得到最优轨迹x^k，迭代停止.

3 突发事件下的轨迹在线重构与跟踪制导

在实际再入飞行过程中，在正常飞行状态下，通过制导实现对参考轨迹的跟踪以消除不确定带来的影响；当遭遇突发事件时，则需要在线重构轨迹作为飞行器新的跟踪目标.本节研究了突发事件下再入轨迹快速重构方法，分析突发事件对约束条件及飞行目标的影响，结合上一节的内容实现轨迹在线快速重构，并结合LQR方法实现对重构轨迹的跟踪制导.

3.1 轨迹快速重构

针对飞行过程中遇到的突发事件，首先考虑RLV偏离参考轨迹的情况，当外界突发干扰导致飞行器大幅度偏离参考轨迹，产生控制系统无法有效修正的大跟踪偏差，此时针对原定目标点进行轨迹重构.为了减少机载计算机的负荷，引入重构阈值走廊的概念.重构阈值走廊基于飞行走廊H_C进行设计.飞行走廊H_C下边界H_CL为飞行器满足热流、动压以及过载约束的最低飞行高度；H_C上边界H_CU由拟平衡滑翔约束式(29)得到，对于确定的倾侧角σ_QEG，拟平衡滑翔约束即为H－V走廊上边界.

$ \left( {g - \frac{{{V^2}}}{r}} \right)\cos \gamma - \frac{{L\cos {\sigma _{{\rm{QEC}}}}}}{m} \le 0. $

(29)

根据飞行走廊确定RLV的轨迹重构条件，设置阈值参数ζ₁, ζ₂∈[0, 1]，飞行器再入重构阈值可表示为

$ \left\{ {\begin{array}{*{20}{l}} {{H_{\max }}\left( v \right) = {H^*}\left( v \right) + {\zeta _1}\left( {{H_{{\rm{CU}}}}\left( v \right) - {H^*}\left( v \right)} \right),}\\ {{H_{\min }}\left( v \right) = {H^*}\left( v \right) - {\zeta _2}\left( {{H^*}\left( v \right) - {H_{{\rm{CL}}}}\left( v \right)} \right).} \end{array}} \right. $

(30)

式中H^*(v)为离线最优轨迹高度，H_max(v)为轨迹重构上边界阈值，H_min(v)为轨迹重构下边界阈值.

当飞行器大幅度偏离参考轨迹，超出重构阈值走廊时，此时基于传统的制导方法已无法消除偏差，需要进行轨迹在线重构求得新的参考，轨迹优化算法在上一节中已经给出.这种情况下，飞行器本身模型以及终端约束都没有发生改变，相较于上一节给出的轨迹优化算法，改变的仅仅是问题的初始条件x₀.假设轨迹重构的开始时间为t_c，相应的状态量以及控制量为x_c和u_c，若以x_c作为轨迹重构的初始状态，由于轨迹重构需要消耗一定的时间，轨迹重构完成时，实际的状态量与x_c将产生偏差.因此，在t_c时刻，对T_ave后的状态进行预测，T_ave为离线轨迹库中轨迹求解的平均时间. $t_{\mathrm{r}}=t_{\mathrm{c}}+T_{\mathrm{ave}}$机根据当前状态量x_c和u_c，采用预测校正制导方法通过积分预测得到.以此预测的状态量x_r作为轨迹重构的初始条件x₀从而消除轨迹重构求解时间带来的状态量偏差.

考虑飞行目标终点变更的情况，此时的首要目标是求得飞行器当前可达到的终点区域，以选择可行的飞行目标终点.计算飞行器的可达域实际上是计算可达域的边界^[12].求解完整的可达域边界可分为两步：1)求解初始状态下的最大经度θ_max、最小经度θ_min、最大纬度φ_max和最小纬度φ_min. 2)将目标函数式(8)中的关于终端状态量的部分选取为经度与纬度的加权组合，即$\phi\left[{\boldsymbol{x}}\left(t_{f}\right)\right]=w\left(\pm \theta_{f}\right)+(1-$$w)\left( \pm {{\varphi }_{\text{f}}} \right), $式中w∈[0, 1]为权重系数，并进行优化求解.其次，机载计算机便预测T_ave′之后，在预测校正制导的作用下，轨迹重构的初始状态x_r，其中T_ave′为包含可达区域求解的轨迹重构平均耗时.之后以x_r以及在可达区域内重新选择的目标终点作为轨迹优化问题的初始约束和终端约束, 基于变信赖域序列凸规划算法完成重构轨迹的求解.

3.2 轨迹跟踪制导

综合考虑对轨迹的跟踪性能与实时性的需求，本文采用LQR这一较为成熟的制导方法实现轨迹跟踪. LQR需要在制导采样点处进行小扰动线性化处理，得到如式(13)所示的用于求解反馈增益的线性时变系统，而线性化的泰勒展开点则可以从前文中得到的收敛性能较好的离散点中选择.

线性化后的轨迹跟踪制导问题可表示为：求最优控制律$\Delta \boldsymbol{u}(t)=-\boldsymbol{K}(t)\Delta \boldsymbol{x}(t), $其中K为反馈增益系数，Δx、Δu为实际状态与参考轨迹之间的状态量偏差和控制量变差.使性能指标函数J_c=$\int_{{t_0}}^{{t_f}} {\left[ {\Delta {\mathit{\boldsymbol{x}}^{\rm{T}}}(t)\mathit{\boldsymbol{Q}}\Delta \mathit{\boldsymbol{x}}(t) + \Delta {\mathit{\boldsymbol{u}}^{\rm{T}}}(t)\mathit{\boldsymbol{R}}\Delta \mathit{\boldsymbol{u}}(t)} \right]} {\rm{d}}t$最小，其中t_f为制导的终止时间，Q、R为加权矩阵，本文基于Bryson原则^[13]选取矩阵系数. LQR制导律通过求解Riccati方程$\boldsymbol{\dot{P}}={{\boldsymbol{A}}^{\text{T}}}\boldsymbol{P}-\boldsymbol{PA}-\boldsymbol{Q}+\boldsymbol{PB}{{\boldsymbol{R}}^{-1}}{{\boldsymbol{B}}^{\text{T}}}\boldsymbol{P}, $，求出P满足P^T=P≥0，从而求得K=R^－1B^TP，进而获得最优反馈控制律Δu.

4 仿真与分析

为了验证设计算法的性能，本文以文献[14]中的RLV模型为例，基于MATLAB 2016a环境实现算法的仿真实验. PC机配置为core i5-8500，主频3 GHz，8 GB内存.

仿真基础参数设置：飞行器质量m=104.305 kg；引力参数μ_g=3.986×10¹⁴m³/s²；飞行器面积S=391.22 m²；海平面大气密度ρ₀=0.002 38 kg/m³；常值系数β=1.387 5×10^－4，k_Q=9.44×10^－5.再入约束参数设置：$\dot{Q}_{\max }$_max=3×10⁶ W/m²，q_max=18 000 Pa，n_max=2.5.升力系数CL=－0.207 0+1.675 6α，阻力系数C_D=0.078 5－0.352 9α+2.040 0α².在变信赖域序列凸规划算法中，选择状态量约束x′_max= $\left[R_{\mathrm{e}}+90000, \frac{\pi}{2}, \frac{\pi}{2}, 9000, \frac{\pi}{2}, \pi, \pi\right]^{\mathrm{T}} ; {\boldsymbol {x}}_{\min }^{\prime}=$$\left[R_{e}, -\frac{\pi}{2}, -\frac{\pi}{2}, 1, -\frac{\pi}{2}, -\pi, -\pi\right]^{\mathrm{T}}$; 控制量约束$u_{\max }^{\prime}=\frac{10 \pi}{180}, \boldsymbol{u}_{\min }^{\prime}=-\frac{10 \pi}{180}$；信赖域ε= $\left[ 10000, \frac{40\pi }{180}, \frac{40\pi }{180}, 500 \right.{{\left. , \frac{40\pi }{180}, \frac{40\pi }{180}, \frac{40\pi }{180} \right]}^{\text{T}}}$；迭代收敛阈值δ= $\left[1 \quad 000, \frac{2 \pi}{180}, \frac{2 \pi}{180}, 50, \frac{2 \pi}{180}, \frac{2 \pi}{180}, \frac{2 \pi}{180}\right]^{\rm {T}}$；惩罚因子μ₁=μ₂=100；性能指标函数对比参数ξ=0.5；变信赖域系数β₁=1.2，β₂=0.5.在YALMIP^[15]工具箱下构建凸优化问题，并利用SeDuMi求解器^[16]针对问题P3进行求解.

为了验证本文提出的变信赖域序列凸规划算法，进行仿真实验并将结果与传统的序列凸规划算法以及基于非线性规划求解器SNOPT^[17]的Gauss伪谱法进行对比.仿真给定初始条件以及终端约束见表 1，考虑终端纬度最大的问题，性能指标函数设定为J=－φ(t_f).再入轨迹对比仿真结果如图 2所示. 图 3给出了变信赖域序列凸规划算法与传统序列凸规划收敛情况的对比，图中红线和绿线分别表示传统序列凸规划算法和变信赖域序列凸规划算法的迭代收敛过程.

从图 2中可以看出变信赖域序列凸规划、传统的序列凸规划和SNOPT得到的轨迹大体趋势是吻合的.变信赖域序列凸规划由于引入了预测校正初始轨迹以及变信赖域策略，收敛性能优于传统的序列凸规划.从倾侧角曲线可以看出，序列凸规划得到的控制曲线存在着抖动.总体来说，基于SNOPT的伪谱法得到的轨迹依然是三者中优化性能较好的.但是考虑求解时间方面，序列凸规划相较于SNOPT则具有显著的优势.仿真得出采用SNOPT方法的轨迹求解CPU时间为55.067 s，采用SCP方法的轨迹求解CPU时间为11.124 s，采用变信赖域SCP方法的轨迹求解CPU时间为5.892 s.此外，图 3给出了变信赖域SCP与SCP在收敛效率上的对比，图中红色实线代表SCP收敛迭代过程，蓝色点划线代表变信赖域SCP迭代过程，SCP求解迭代7次，变信赖域SCP求解迭代3次.通过上述仿真对比，验证了变信赖域序列凸规划算法在收敛性能和求解速度上的改进.

良好的求解速度以及较好的求解性能使得变信赖域序列凸规划算法具备应用于轨迹在线重构的潜力.分别针对大幅度偏离参考轨迹、飞行目标终点变更的情况进行轨迹重构仿真.

首先针对大幅度偏离参考轨迹情况，考虑沿离线最优轨迹飞行至某一状态点处x=[51 800 m, －7, 26°, 4 017 m/s，－0.17°, 60°]^T，给定外界扰动带来的偏差Δh=5 000 m，Δv=100 m/s，采用ζ₁=ζ₂=0.2的阈值走廊，轨迹重构结果如图 4所示，从图 4可看出在扰动作用下，轨迹超出阈值重构走廊，需要进行轨迹在线重构，基于本文提出的轨迹重构方法，可以在短时间内生成一条满足约束的最优轨迹，求解CPU时间为5.73 s. 图 5给出了基于LQR制导律对初始参考轨迹以及重构轨迹跟踪的结果以及扰动下若对原参考轨迹进行跟踪或者采用预测校正制导律的制导结果.从图 5中可以看出，过大的偏差导致LQR制导律不具备跟踪原参考轨迹的能力，若依然对原参考轨迹进行跟踪，产生的轨迹将逐渐偏离最优轨迹，如图中黑色虚线所示，而重构后的轨迹能很好地满足LQR制导律的跟踪条件，引导飞行器安全再入返回，如图中红色实线所示.此外，由于重构轨迹求解时间预测策略的存在，轨迹重构时初始状态量偏差处于很小的范围内，重构求解平均时间取T_ave=6 s.与预测校正制导进行对比，采用预测校正方法得到的轨迹仅仅是一条可行轨迹，而非本文中得到的最优轨迹.在求解制导指令速度方面，预测校正制导每次求解需要经过45次迭代，总用时0.265 s；而LQR制导律每次求解仅需0.04 s.综上所述，基于轨迹重构的LQR制导律在实时性和轨迹性能上都优于预测校正制导.

考虑与大幅度偏离参考轨迹情况下相同的状态点处变更目标终点，基于可达域选择新的目标终点[θ_f, φ_f]=[10°，30°]，三维轨迹重构与LQR跟踪制导结果如图 6所示，重构后的轨迹引导飞行器飞向新的目标终点.

5 结论

1) 针对RLV的轨迹再入轨迹重构问题，提出变信赖域序列凸规划算法对再入轨迹在线快速求解.

2) 在传统序列凸规划算法的基础上，采用预测校正算法求解序列凸规划的初始迭代轨迹，以此提升算法的收敛效率，此外为了加快算法收敛速度，设计了信赖域更新策略.

3) 基于提出的轨迹快速求解算法，针对RLV再入过程中可能发生的如大幅度偏离参考轨迹及目标点变更等突发事件，在线重构轨迹并考虑轨迹重构耗时，对实际轨迹重构初始点进行了预测以抑制可能出现的状态量跳变，同时基于LQR制导方法实现对参考轨迹的快速跟踪.

4) 仿真结果表明所提出的轨迹求解算法具备良好的收敛性能以及求解速度，重构轨迹能够有效地引导RLV再入返回.