近空间可变翼飞行器(near-space vehicle，NMV)是一种翼形可变的工作于近空间的飞行器，其外形结构呈三角形，采用翼身融合构型。为了适应不同飞行环境，飞行器通过伸缩翼来改变气动外形，完成不同的飞行任务，具有重要的研究意义与军事应用价值，受到了国内外众多学者和研究机构的重视^[1-3]。

近空间飞行器轨迹优化一般是非线性、带有控制约束及状态约束的最优控制问题^[4]。由于本文研究对象的特殊性，在控制变量中加入小翼伸缩量，所以该问题可以转变为受多约束条件下的高度非线性最优控制问题。该类最优控制问题根据求解方式分成两大类。1)间接法：将问题转化为Hamilton边值问题而后求解，此类方法在满足最优性的一阶必要条件同时求解的精度高，但存在收敛域小x(t)、难以估算共轭变量初值等缺点；2)直接法：将时序连续的问题转化为离散的NLP问题来获取最优解^[5]，存在求解精度较低等缺点。

近年来，在解决轨迹优化问题上，与其他方法相比，直接伪谱法表现出全局收敛及高效率、高精度，因此在航空航天领域得到了广泛的应用^[6-7]。文献[8-9]采用Gauss伪谱法解决在多约束下的滑翔轨迹优化问题；文献[10]应用Gauss伪谱法解决了临近空间飞行器上升段的轨迹优化问题。文献[11]应用一种改进的点位经纬度变换条件GPM的方法，提升了航迹规划。文献[12]采用Gauss伪谱法针对可变翼飞行器小翼伸缩与最省燃油轨迹进行优化。

上述优化方法存在一定的局限性：1)当目标函数复杂多变时，须增加优化函数的配点数，才可在一定程度上提升优化精度；2)随着配点数的增多，状态微分矩阵的维数也随即增大，导致计算量大、效率低；3)收敛速度会因控制变量的不连续而降低。针对上述局限性，本文提出一种改进的多段整合优化的高斯伪谱法优化策略，首先在时间区间上进行逐一分段，而后采用低阶多项式插值逼近状态变量，通过逐步增加分段的个数，提高逼近精度。

1 非线性最优控制及多段轨迹优化问题描述

考虑多约束条件的非线性最优控制问题，即Bolza问题^[13]：

$ \begin{array}{l} {\rm{min}}:J = \varPhi (x({t_0}),{t_0},x({t_f}),{t_f}) + \int_{{t_0}}^{{t_f}} g (x(t),u(x),t){\rm{d}}t\\ {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\rm{s}}{\rm{.t}}{\rm{.}}\left\{ {\begin{array}{*{20}{l}} {\dot x(t) = f(x(t),u(t),t)}\\ {E(x(t),u(t),t) = 0}\\ {C(x(t),u(t),t) \le 0,t \in [{t_0},{t_f}]} \end{array}} \right. \end{array} $

(1)

式中: x(t)、u(t)分别为状态变量和控制变量(优化变量)；t₀、t_f为起始时间和末端时间；E为等式约束；C表示不等式约束。

1.1 一般Gauss伪谱法

针对Bolza问题通过引入时间变量τ，将时间区间[t₀，t_f]投射到[-1, 1]^[14]，采用非等间距的Gauss伪谱法求解：

$ t = \frac{{{t_f} - {t_0}}}{2}\tau + \frac{{{t_f} + {t_0}}}{2},\;\:\tau \in [ - 1,1] $

(2)

将时序连续的无限维最优控制问题转换得到有限维的NLP问题：

$ \begin{array}{l} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} \begin{array}{*{20}{l}} {{\rm{min}}:J = \varPhi (x( - 1),{t_0},x(1),{t_f}) + }\\ {\frac{{{t_f} - {t_0}}}{2}\int_{ - 1}^1 g (x(\tau ),u(\tau ),\tau ;{t_0},{t_f}){\rm{d}}\tau } \end{array}\\ {\rm{s}}{\rm{.t}}{\rm{.}}\left\{ \begin{array}{l} \dot x(t) = \frac{{{t_f} - {t_0}}}{2}f(x(\tau ),u(\tau ),\tau ;{t_0},{t_f})\\ E(x( - 1),u( - 1),{t_0};\\ {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} x(1),u(1),{t_f};x(\tau ),u(\tau ),\tau ) = 0\\ C(x(\tau ),u(\tau ),\tau ;{t_0},{t_f}) \le 0 \end{array} \right. \end{array} $

(3)

采用全局Lagrange多项式插值，动力学微分方程和控制变量离散化，路径约束、边界约束、终端状态约束与目标函数的精确近似。可以建立NLP问题：

$ \begin{array}{l} \begin{array}{*{20}{l}} { {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\rm{min}} :J = \varPhi (x( - 1),{t_0},x(1),{t_f}) + }\\ {{\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} \frac{{{t_f} - {t_0}}}{2}\sum\limits_{k = 1}^N {{{\hat w}_k}} g({X_k},{U_k},{\tau _k};{t_0},{t_f})} \end{array}\\ {\rm{s}}{\rm{.t}}{\rm{.}}\left\{ {\begin{array}{*{20}{l}} {\sum\limits_{i = 0}^N {{D_{k,i}}} (\tau ){X_k} - \frac{{{t_f} - {t_0}}}{2}f({X_k},{U_k},{\tau _k};{t_0},{t_f}) = 0}\\ {C({X_k},{U_k},{\tau _k};{t_0},{t_f}) \le 0,k = 1,2 \cdots ,N}\\ {{E_1}(x( - 1),u( - 1),{t_0};) = 0}\\ {{E_2}(x(1),u(1),{t_f}) = 0}\\ {{X_f} - {X_0} + \frac{{{t_f} - {t_0}}}{2}\sum\limits_{k = 1}^N {{w_k}} f({X_k},{U_k},{\tau _k};{t_0},{t_f}) = 0} \end{array}} \right. \end{array} $

(4)

1.2 分段优化问题描述

由于采用一般的高斯伪谱法进行优化时，是将整个时间区间内的状态变量用一个全局Lagrange函数逼近，一般需要较多的配点才能达到所需的精度要求。为了减少配点数，降低优化时间，提高优化精度，提出采用分段优化，即将时间段划分p个区间进行优化。记为：

$ [{t_0},{t_1}],[{t_1},{t_2}], \cdots ,[{t_{p - 1}},{t_p}] $

(5)

式中：t_p=t_f，[t_p-1, t_p](p=1，2，…，P)表示第p个阶段的时间段，记为阶段p，其中t_p(p=1，2，…，p-1)称为分割点。由于对时间段划分了p个区间后，其存在p-1个间断点，在每个间断点上需要考虑其系统状态的连续性问题，因此分段优化需要满足以下连接约束^[15-16]：

$ {x^{(p + 1)}}({t_p}) - {x^{(p)}}({t_p}) = 0,\;\:p = 1,2, \cdots ,p - 1 $

(6)

分段示例如图 1所示。图 1中可以看出存在3个阶段，其中包含2个连接约束点。每个连接点即为前一段的结束时刻和下一段的初始时刻。可以在每一个优化阶段选择合适的LG配点数n_p，(p=1，2，…，p-1)，在每段上，根据LG配点拟合状态量与控制量，给出每个子段间约束条件并建立目标函数等，最后然后采用SQP算法求解。分段后的Bolza型问题表示为：

性能指标函数：

$ \begin{array}{*{20}{l}} {{\rm{min}} :{J_p} = \varPhi ({x^{(1)}}({t_0}),{t_0},{x^{(p)}}({t_p}),{t_p}) + }\\ {{\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} \sum\limits_{p = 1}^p {\int_{{t_{p - 1}}}^{{t_p}} g } ({x^{(p)}}(t),{u^{(p)}}(t),t){\rm{d}}t} \end{array} $

(7)

状态方程：

$ \begin{array}{*{20}{c}} {{{\dot x}^{(p)}}(t) = f({x^{(p)}},{u^{(p)}},t)t \in [{t_{p - 1}},{t_p}]}\\ {p = 1,2, \cdots ,P} \end{array} $

(8)

边界约束：

$ E({x^{(1)}}({t_0}),{t_0},{x^{(p)}}({t_p}),{t_p}) = 0 $

(9)

路径约束：

$ \begin{array}{*{20}{c}} {C({x^{(p)}}(t),{u^{(p)}}(t),t) \le 0,t \in [{t_{p - 1}},{t_p}],}\\ {p = 1,2, \cdots ,P} \end{array} $

(10)

连接约束：

$ {x^{(p + 1)}}({t_p}) - {x^{(p)}}({t_p}) = 0 $

(11)

式中:x^(p)∈Rⁿ, u^(p)∈R^m分别是阶段p=1, 2, …, P上的状态变量和待优化的控制变量。

2 爬升多子段优化模型与优化策略建立

针对爬升过程中小翼伸缩最省燃料轨迹优化，综合考虑爬升过程中声速、密度、引力及推力变化等因素的影响，根据发动机工作状态及飞行状态特性对爬升段合理划分，建立多子段的优化模型。

2.1 近空间可变翼飞行器构型参数

本文的研究对象相关外形参数如下：飞行器总质量1.002×10⁵ kg，机身总长60.69 m，机翼可伸缩面积30 m²，翼弦长30 m，俯仰力矩惯性积8.47×10⁶ kg·m²。与一般的飞行器相比，机翼翼梢装有可伸缩小翼。根据飞行需求，通过伸缩小翼改变翼型，引起不同的气动效应，适应不同的飞行环境，完成各种飞行任务。

2.2 可变翼飞行器模型建立

1) 动力学模型：可变翼飞行器爬升段数学模型^[17]为：

$ \left\{ {\begin{array}{*{20}{l}} {\dot V = \frac{{T{\rm{cos}}{\kern 1pt} {\kern 1pt} {\kern 1pt} \alpha - D}}{m} - \frac{\mu }{{{r^2}}}{\rm{sin}}{\kern 1pt} {\kern 1pt} {\kern 1pt} \gamma }\\ {\dot h = V{\rm{sin}}{\kern 1pt} {\kern 1pt} {\kern 1pt} \gamma }\\ {\dot \gamma = \frac{{L + T{\rm{sin}}{\kern 1pt} {\kern 1pt} {\kern 1pt} \alpha }}{{mV}} - \frac{{(\mu - {V^2}r){\rm{cos}}{\kern 1pt} {\kern 1pt} {\kern 1pt} \gamma }}{{V{r^2}}}}\\ {\dot m = - \frac{T}{{{I_{sp}}{g_0}}}} \end{array}} \right. $

(12)

$ {I_{sp}} = \left\{ {\begin{array}{*{20}{l}} {4{\kern 1pt} {\kern 1pt} {\kern 1pt} 700 - \frac{h}{{100}},}&{Ma < 4}\\ {5{\kern 1pt} {\kern 1pt} {\kern 1pt} 680 - 245Ma - \frac{h}{{100}},}&{Ma \ge 4} \end{array}} \right. $

(13)

式中：V、h、γ、m、α分别为飞行器速度、高度、航迹角、质量、迎角；L、D、T分别为升力、阻力、发动机推力；I_sp、g₀、r、μ分别为燃料比冲、海平面重力加速度、地球半径、引力参数；Ma为马赫数。

2) 发动机模型：

$ T = 0.5\rho {V^2}(S + {S_{{\rm{ xiaoyi }}}}){C_T} $

(14)

式中：S为飞行器有效参考面积；S_xiaoyi小翼的伸缩面积。C_T为推力系数，且：

$ {C_T} = \left\{ {\begin{array}{*{20}{l}} {0.025{\kern 1pt} {\kern 1pt} {\kern 1pt} 76\beta ,}&{\beta \ge 1}\\ {0.022{\kern 1pt} {\kern 1pt} {\kern 1pt} 4 + 0.003{\kern 1pt} {\kern 1pt} {\kern 1pt} 36\beta ,}&{\beta < 1} \end{array}} \right. $

(15)

式中β为发动机节流阀调定值状态。

3) 大气环境模型^[18]：

$ {g = {g_0}{{\left( {\frac{{{r_0}}}{{{r_0} + h}}} \right)}^2}} $

(16)

$ {\rho = {\rho _0}{{\rm{e}}^{ - h/7315.2}}} $

(17)

式中：r₀为地球半径；g₀即标准重力加速度；ρ₀=1.226 6 kg/m³，即地球表面的大气密度。本文数值仿真验证中采用的标准重力加速度g₀=9.806 65 m/s²，地球半径r₀=6 356 766 m。

4) 气动力模型：

飞行器飞行过程中受到阻力D和升力L的计算公式为：

$ \begin{array}{*{20}{l}} {L = 0.5\rho {V^2}(S + {S_{{\rm{ xiaoyi }}}}){C_L}}\\ {D = 0.5\rho {V^2}(S + {S_{{\rm{ xiaoyi }}}}){C_D}} \end{array} $

(18)

式中：C_L和C_D分别为升力系数和阻力系数，其中小翼伸缩时的升阻力系数通过插值法得到。图 2给出小翼伸、缩状态下气动参数变化与α、Ma的三维关系曲面。

结合近空间飞行环境模型，通过高精度插值法获得小翼伸缩时的气动参数，建立气动参数拟合数据库，根据建立的气动参数拟合数据库查询各项气动数据，进行飞行器非线性数学模型的后续计算，具体获取过程如图 3所示。

综上，通过气动数据库查询获得各项气动数据之后，代入到各气动力和气动力矩中，即可进行近空间飞行器非线性数学模型的后续计算。

2.3 爬升各子段划分

近空间可变翼飞行器在爬升模态下需要经历平流层、中间层和部分热层。随着自身速度和高度的增加，飞行状态发生显著变化，飞行约束条件愈加苛刻。针对爬升模态下可变翼飞行器的参数状态变化、发动机推力变化等特性，将爬升模态划分3个爬升子段进行研究。

近空间可变翼飞行器从起飞段加速达到一定速度后开始拉升进入爬升段阶段。在飞行爬升过程中，飞行器受到的发动机推力较大，燃料消耗严重，飞行器的质量、重心、惯性力矩等都会产生明显的变化。根据爬升段发动机工作状态变化及飞行器高度、速度变化将其划分成3个子段，表 2给出了飞行器各分段的状态及对飞行器的飞行要求。

依据飞行器在爬升时的状态变化及发动机的推力特性，将上升轨迹合理分段：将原最优控制问题转化为多段最优控制问题后，采用改进的多段高斯伪谱法进行整合并进行优化计算。

2.4 爬升段约束条件 2.4.1 边界约束

近空间飞行器一般是通过助推器助推到预定的高度和速度后，启动发动机，给定飞行器初始与末端的状态：

$ \left\{ {\begin{array}{*{20}{l}} {h({t_0}) = {h_0},}&{h({t_f}) = {h_f}}\\ {V({t_0}) = {V_0},}&{V({t_f}) = {V_f}}\\ {\gamma ({t_0}) = {\gamma _0},}&{\gamma ({t_f}) = {\gamma _f}}\\ {m({t_0}) = {m_0}}&{} \end{array}} \right. $

(19)

2.4.2 路径约束

飞行器在加速爬升过程中，受到大气压强作用，会破坏机体产生形变并摩擦生热。因此，飞行过程中动压、热流率及过载须满足一定的约束^[18]。

1) 热流率约束：

$ \dot Q = C{\rho ^p}{V^q} $

(20)

式中：C=7.968 6×10^-5; p=0.5;q=3.15。

2) 动压约束：

$ {q_{{\rm{min}}}} \le q = \frac{1}{2}\rho {v^2} \le {q_{{\rm{max}}}} $

(21)

3) 过载约束：

$ n = \frac{{qS\sqrt {(C_L^2 + C_D^2)} }}{{mg}} \le {n_{{\rm{max}}}} $

(22)

2.4.3 控制量与状态量约束

飞行器执行机构在爬升过程中受限，因而控制量也受到一定的范围内约束，且爬升中满足机动要求状态量相应受到约束。对于可变翼飞行器，存在迎角、发动机节流阀、小翼伸缩面积，u=[α  β  S_xiaoyi]，3个控制量和速度、高度、航迹角、质量，x=[V  h  γ  m]，4个状态量。

$ \left\{ {\begin{array}{*{20}{l}} {{\alpha _{{\rm{min}}}} \le \alpha \le {\alpha _{{\rm{max}}}},\;\:{V_{{\rm{min}}}} \le V \le {V_{{\rm{max}}}}}\\ {{\beta _{{\rm{min}}}} \le \beta \le {\beta _{{\rm{max}}}},\;\:{h_{{\rm{min}}}} \le h \le {h_{{\rm{max}}}}}\\ {{S_{{\rm{ xiaoyi}}{{\rm{ }}_{{\rm{min}}}}}} \le {S_{{\rm{ xiaoyi }}}} \le {S_{{\rm{ xiaoyi}}{{\rm{ }}_{{\rm{max}}}}}},\;\:{\gamma _{{\rm{min}}}} \le \gamma \le {\gamma _{{\rm{max}}}}} \end{array}} \right. $

(23)

2.4.4 间断点约束

针对划分的3个子阶段，出现的2个间断点，需要满足如下间断点约束条件：

$ \left\{ {\begin{array}{*{20}{l}} {V_{{t_1}}^{{p_2}} - V_{{t_1}}^{{p_1}} = 0,}&{V_{{t_2}}^{{p_3}} - V_{{t_2}}^{{p_2}} = 0}\\ {h_{{t_1}}^{{p_2}} - h_{{t_1}}^{{p_1}} = 0,}&{h_{{t_2}}^{{p_3}} - h_{{t_2}}^{{p_2}} = 0}\\ {\gamma _{{t_1}}^{{p_2}} - \gamma _{{t_1}}^{{p_1}} = 0,}&{\gamma _{{t_2}}^{{p_3}} - \gamma _{{t_2}}^{{p_2}} = 0}\\ {m_{{t_1}}^{{p_2}} - m_{{t_1}}^{{p_1}} = 0,}&{m_{{t_2}}^{{p_3}} - m_{{t_2}}^{{p_2}} = 0} \end{array}} \right. $

(24)

2.5 优化指标函数建立

飞行器在爬升飞行过程中, 动能和势能增长迅速，需要为其提供大量的能量。基于本文所研究的近空间可变翼飞行器，其在各个子段爬升过程中质量的减少近似为燃料的消耗。由于飞行器执行任务过程中携带的燃油有限，为了能节省更多燃料用于后续任务，因此建立以最省燃料为优化指标进行爬升，即飞行器在爬升末端质量保持最大：

$ {J_m} = - {m_{{t_f}}} = \mathop {{\rm{max}}}\limits_t (m({t_0}) - \sum\limits_{p = 1}^{p = 3} ( \int_{{t_{p - 1}}}^{{t_p}} {{{\dot m}_p}} (\tau ){\rm{d}}\tau )) $

(25)

2.6 爬升多子段小翼伸缩优化模型建立

首先针对近空间可变翼飞行器划分的P个子段，确定每个子段的初始配点数n_p，然后在每个子段上进行离散化处理，将时间段[t_p-1, t_p]变换到[-1, 1]上：

$ t = \frac{{{t_p} - {t_{p - 1}}}}{2}\tau + \frac{{{t_p} + {t_{p - 1}}}}{2} $

(26)

式中:τ表示时间t在[t_p-1, t_p]上转化后的时间变量，p=1, 2, …, P。

分别将x=[V  h  γ  m]在每子段LG配点上进行离散化：

$ \begin{array}{*{20}{c}} {{\mathit{\boldsymbol{x}}^{(p)}}(\tau ) = {\mathit{\boldsymbol{X}}^{(p)}}(\tau ) = \sum\limits_{i = 0}^{{n_p}} {L_i^{(p)}} (\tau ){\mathit{\boldsymbol{X}}^{(p)}}(\tau _i^{(p)}) = }\\ {\sum\limits_{i = 0}^{{n_p}} {L_i^{(p)}} (\tau )\mathit{\boldsymbol{X}}_i^{(p)}} \end{array} $

(27)

$ L_i^{(p)}(\tau ) = \prod\limits_{j = 0,j \ne i}^{{n_p}} {\frac{{\tau - \tau _j^{(p)}}}{{\tau _i^{(p)} - \tau _j^{(p)}}}} ,\;\:i = 0,1, \cdots ,{n_p} $

(28)

对式子(27)状态变量x两边分别求导可得：

$ {\mathit{\boldsymbol{\dot X}}^{(p)}}(\tau ) = \sum\limits_{i = 0}^{{n_p}} {\dot L_i^{(p)}} (\tau )\mathit{\boldsymbol{X}}_i^{(p)} $

(29)

令$ D_{k,i}^{(p)} = \dot L_i^{(p)}(\tau )$，代入式(29), 可得到状态变量在LG配点上的离散状态方程：

$ {\mathit{\boldsymbol{\dot X}}^{(p)}}(\tau ) = \sum\limits_{i = 0}^{{n_p}} {\mathit{\boldsymbol{D}}_{k,i}^{(p)}} \mathit{\boldsymbol{X}}_i^{(p)} $

(30)

式中：D_{k, i}^(p)称为阶段p的状态微分矩阵，是一个n_p×(n_p+1)的非稀疏性矩阵。

其次同样将控制变量u=[α  β_c  S_xiaoyi]在LG配点上进行离散化：

$ {{\mathit{\boldsymbol{u}}^{(p)}}(\tau ) \approx {\mathit{\boldsymbol{U}}^{(p)}}(\tau ) = \sum\limits_{i = 1}^{{n_p}} {\tilde L_i^{(p)}} (\tau )\mathit{\boldsymbol{U}}_i^{(p)}} $

(31)

$ {\tilde L_i^{(p)}(\tau ) = \prod\limits_{j = 0,j \ne i}^{{n_p}} {\frac{{\tau - \tau _j^{(p)}}}{{\tau _i^{(p)} - \tau _j^{(p)}}}} ,\;\:i = 0,1, \cdots ,{n_p}} $

(32)

最后同理将爬升段各段的约束条件、边界点连接约束及性能指标表示为^[20-21]：

$ {\mathit{\boldsymbol{X}}_0^{(p)} - \mathit{\boldsymbol{X}}_0^{(p)} - \sum\limits_{i = 1}^{{n_p}} {\mathit{\boldsymbol{x}}_i^{(p)}} \sum\limits_{k = 1}^{{n_p}} {w_k^{(p)}} D_{k,i}^{(p)} = 0} $

(33)

$ {\mathit{\boldsymbol{X}}_{{N^{(p)}} + 1}^{(p)} - \mathit{\boldsymbol{X}}_0^{(p)} - \sum\limits_{i = 1}^{{n_p}} {\mathit{\boldsymbol{x}}_i^{(p)}} \sum\limits_{k = 1}^{{n_p}} {w_k^{(p)}} D_{k,i}^{(p)} = 0} $

(34)

$ \begin{array}{*{20}{c}} {\int_{{t_{p - 1}}}^{{t_p}} g ({\mathit{\boldsymbol{x}}^{(p)}}(t),{\mathit{\boldsymbol{u}}^{(p)}}(t),t){\rm{d}}t \approx }\\ {\frac{{{t_p} - {t_{p - 1}}}}{2}\sum\limits_{k = 1}^{{n_p}} {w_k^{(p)}} g(\mathit{\boldsymbol{X}}_k^{(p)},\mathit{\boldsymbol{U}}_k^{(p)},\tau _k^{(p)})} \end{array} $

(35)

根据离散结果最后得到改进的多段高斯伪谱法优化模型：

$ \begin{array}{l} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} \begin{array}{*{20}{c}} {{\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\rm{min}}:{J_p} = \varPhi (\mathit{\boldsymbol{x}}_0^{(1)}({t_0}),{t_0},\mathit{\boldsymbol{x}}_{{N^{(p)}} + 1}^{(p)}({t_p}),{t_p}) + }\\ {\sum\limits_{p = 1}^p {\frac{{{t_p} - {t_{p - 1}}}}{2}} \sum\limits_{k = 1}^{{N^{(p)}}} {w_k^{(p)}} g(\mathit{\boldsymbol{X}}_k^{(p)},\mathit{\boldsymbol{U}}_k^{(p)},\tau _k^{(p)};{t_{p - 1}},{t_p})} \end{array}\\ {\rm{s}}{\rm{.t}}{\rm{.}}\left\{ \begin{array}{l} \sum\limits_{i = 0}^{{N^{(p)}}} {D_{k,i}^{(p)}} (\tau )\mathit{\boldsymbol{X}}_i^{(p)} - \frac{{{t_p} - {t_{p - 1}}}}{2}f(\mathit{\boldsymbol{X}}_k^{(p)},\mathit{\boldsymbol{U}}_k^{(p)},\tau _k^{(p)};\\ {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {t_{p - 1}},{t_p}) = 0\\ C({\mathit{\boldsymbol{x}}^{(p)}}(t),{\mathit{\boldsymbol{u}}^{(p)}}(t),t) \le 0,k = 1,2 \cdots ,N,\\ {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} p = 1,2, \cdots ,p\\ \begin{array}{*{20}{l}} {{E_1}(\mathit{\boldsymbol{x}}_0^1({t_0}),{t_0}) = 0}\\ {{E_2}({\mathit{\boldsymbol{x}}^p}({t_p}),{t_p}) = 0}\\ {\mathit{\boldsymbol{X}}_{{N^{(p)}} + 1}^{(p)} - \mathit{\boldsymbol{X}}_0^{(p)} - \sum\limits_{i = 1}^{{N^{(p)}}} {\mathit{\boldsymbol{x}}_i^{(p)}} \sum\limits_{k = 1}^{{N^{(p)}}} {w_k^{(p)}} D_{k,i}^{(p)} = 0} \end{array} \end{array} \right. \end{array} $

(36)

式中：由于选取的优化指标为最省燃油爬升飞行，所以J_p可由式(25)J_m变换得到，C为不等式形式下的路径约束包括热流、动压、过载等，E₁、E₂为等式形式下的初始和终端的边界约束。

2.7 多段优化策略建立

改进多段高斯维谱法优化策略流程如图 4。

利用上述改进的多段高斯维谱法优化策略进行飞行器爬升段优化。根据明确的物理特性对时间分段后，通过选取少量的LG配点进行优化，一方面有利于提高计算效率，另一方面保证优化的精度与收敛性。

3 仿真分析 3.1 初始条件及参数设定

针对上述改进的优化策略，对小翼伸缩量及最省燃油轨迹进行仿真实验，每段参数设计如下：

1) 各爬升子段初值与终值设置如下：

爬升速度：V₀=590 m/s，V_tf=3 100 m/s；

爬升高度：h₀=4 500 m，h_f=36 000 m；

航迹角：γ₀=0.12 rad，γ_tf=0 rad；

飞行器总质量：m₀=100 200 kg，max(m_tf)；

爬升时间：t₀=0 s，t_f=310 s。

2) 控制变量u=[α, β_c, S_xiaoyi]约束如下：

迎角：α_min=-3°，α_max=8°；

油门开合状态：β_min=0，β_max=1；

小翼伸缩面积：S_xiaoyimin=0 m²，S_xiaoyimax=30 m²。

3) 状态变量限制范围如下：

速度：V_min=590 m/s，V_max=3 150 m/s；

高度：h_min=4 500 m，h_max=36 500 m；

航迹角：γ₀=0 rad，γ_tf=0.785 3 rad；

总质量：m_max=100 200 kg。

4) 飞行器在各个子段爬升过程中，每个子段的初始与前一子段的终端值必须保持平滑连接，因此各子段连接点的值域约束范围如下：

第1、2阶段连接点：V_tfmin^p₁=1 100 m/s，V_tfmax^p₁=1 300 m/s，h_tfmin^p₁=16 000 m，h_tfmax^p₁=18 000 m

第2、3阶段连接点：V_tfmin^p₂=2 700 m/s，V_tfmax^p₂=2 900 m/s，h_tfmin^p₂=30 000 m，h_tfmax^p₂=32 000 m

5) 路径约束参数如下：热流约束：$ {{\dot Q}_{{\text{max}}}} = 8 \times {10^5}{\text{W}}/{{\text{m}}^2}$

动压约束：q_min=4×10⁴ Pa，q_max=1.6×10⁴ kPa；

过载约束：n_max=3.5。

6) 其他参数：飞行器参考面积S=360 m²，地表标准重力加速度g₀=9.806 65 m/s²，地表大气密度值ρ₀=1.226 6 kg/m³，地球半径r₀=6 356 766 m，每个子段初始LG配点数均为5。

3.2 仿真结果分析

为了获得爬升高精度优化结果，当进行单段优化时，为了实现误差低于10%，配点数N必须大于80，增大优化时间，影响爬升优化效率。采用多段优化，将爬升时间划分为3个子段[0, 60]、[60, 260]、[260, 310]、LG最终收敛配点个数分别为n₁=20，n₂=35，n₃=20，总优化时长为178.9 s，最大相对误差为6.73%，在同等配点数情况下，多段优化逼近比未分段的逼近精度提高了3.27%。仿真结果如图 5。

由图 5 (a)、(b)仿真速度、高度可知，飞行器在每个子段末都能稳定达到最终收敛值或处于约束条件内。随着速度的不断增加，初始段高度值先稳定而后迅速增加。这种先加速后爬升的飞行轨迹有2大优势：1)先加速可使飞行器快速达到安全速度; 2)先加速后爬升可以满足最省燃油优化指标。图 5(c)为航迹角曲线，随着爬升进行，每子段内航迹角基本都是先增大后减小，这也说明速度矢量先是向垂直方向增加，然后往水平方向移动，先加速向上爬升而后爬升趋于缓慢，与每个阶段速度、高度曲线相对应。最终航迹角收敛为零，飞行器达到巡航速度高度，进入平稳巡航飞行阶段。

图 6是3个优化控制变量变化曲线。随着飞行高度、速度的不断增加，迎角在不断减少，最后稳定在一定的幅值范围内。每个阶段中发动机节流阀状态都相对稳定在一定幅度范围内，但在阶段连接点处存在跳变情况，跳变原因：1)由于在阶段的连接处，为了飞行状态的平滑连接，油门存在突变; 2)小翼在阶段连接处伸缩变化引起油门大幅度的跳变。小翼在初始爬升段的时候，由于初始爬升速度较小，此时小翼伸缩对飞行器的升阻比影响较大，小翼缩回会更大程度上减少阻力的产生，优化表明，爬升初始段，在速度较小的情况下，小翼选择收回，更易于飞行器加速爬升减少飞行燃油消耗。切换过后，油门状态保持较小的恒定值，对于加速爬升油量控制具有可观的影响。飞行器在爬升至150 s左右时，加速到高超声速阶段，由于发动机推力在做转换，且小翼伸缩量对升阻比影响较大，伸出小翼的同时气动焦点后移，使得稳定性与机动特性得到增强。因此小翼在150 s处从伸出到缓慢收回，此过程对飞行器油量节省也具有一定的影响。在最后一个子段，由于小翼的伸出可以增大阻力特性，为了降低油量消耗，采用伸出的小翼，此策略意在配合阻力使得加速度逐渐减少，最后结合油门使得飞行器达到最终收敛巡航速度和高度。

图 5 (d)为质量曲线，与固定翼状态和伸缩翼对比可知，小翼伸缩可节省燃油爬升，在高超声速段，小翼变化对阻力变化明显，找到合适的小翼状态会缩减燃料消耗。整个爬升过程中耗油约12 t，与小翼收回相比节省燃料约800~1 000 kg，相比完全展开小翼节省燃料约1 600~2 000 kg，相当于爬升段耗油量的10%，为飞行完成下个任务节省了足够的油量。

4 结论

1) 为了能模拟真实飞行环境与气动数据，本文综合考虑近空间区域声速、密度、引力及推力等众多变化因素对飞行器的影响，通过气动参数拟合数据库，使得建模更为准确。

2) 设计并改进多段高斯伪谱法优化策略，有效解决爬升模态下的多段小翼伸缩轨迹优化，在减少计算耗时的同时提高了优化精度，进一步提升多配点问题的优化效率。

3) 针对可伸缩小翼这一特性，对飞行器爬升模态下进行最省燃油优化，应用伸缩翼达到节省大量的燃油，得到高效飞行收益，为今后可变翼飞行研究提供了参考，具有一定的工程应用价值。