基于模型预测控制的运载火箭轨迹优化方法

Model Predictive Control Based Launch Vehicle Trajectory Optimization Method

摘要: 本文针对运载火箭上升段的轨迹优化问题的高度不确定性和严格的飞行过程约束，利用了模型预测控控制(MPC)方法来构建最优控制问题，并将其采用多重打靶法转换为非线性规划(NLP)问题，最后采用内点法进行了在线求解。本文通过利用MPC的在线预测与优化能力减小不确定性影响，以及MPC的过程约束处理能力处理强约束条件，仿真结果表明该方法能有效适应不确定性及过程约束要求。

Abstract: A model predictive control (MPC) based trajectory optimization method is used in this paper for the high uncertainty and strict flight process constraints at the ascending stage of the launch vehicle. The model predictive control (MPC) method is used to construct the optimal control problem, and the multiple shooting method is used to convert it into a nonlinear planning (NLP) problem. At last, the interior point method was used to solve online. By using MPC's online prediction and optimi-zation capabilities to reduce the impact of uncertainty, and MPC's process constraint processing capabilities to deal with strong constraints, the simulation results show that the method can effec-tively adapt to uncertainties and process constraints.

1. 引言

传统的上升制导算法是在飞行的早期高动压力段采用“开环”制导，即第一阶段，从升空到固体火箭助推器(SRB)分离和之后根据预定的时间或事件，之后切换到闭环真空制导方案 [1]。开环模式通常使用预先加载的姿态角指令，制导指令一般在地面生成，并随将发射当天的风更新到弹道设计中。真空段的制导方法基于最优控制理论，最早源于土星火箭的迭代制导模式(IGM) [2] 以及航天飞机的动力显式制导(PEG) [3]。

上升段制导的挑战在于大气内飞行段，空气动力和风的存在极大地增加了最优控制问题的求解难度问题，使解决方法难以可靠，充分地收敛，目前这段还是会是开环的。但是，开环制导没有足够的鲁棒性来应付较大的系统建模不确定性，尤其是对于新型发射任务来说。这类问题只能通过上升段的轨迹优化来降低影响。

一般来说，对于模型火箭，推力水平是固定的，当前变推力发动机还较少见。上升段的轨迹是通过寻求最优的推力方向，构成以到达路径最小燃料为目标的最优控制问题，近些年来国外专家学者不断对这类问题进行分析与提升 [4] [5]，并在所涉及的多个方面进行了探索。在求解数值算法方面，一般分为直接法和间接法。直接方法 [6] [7] [8] 将非线性控制问题转化为最优控制问题，并通过离散化方法转化为非线性规划问题(NLP)进行求解，一般依赖非线性(NLP)软件(IPOPT、BOCOP、GESOP、SNOPT等)。直接法可以轻松处理任何动力学约束条件，但是计算代价较大，且不容易获取收敛解。间接法 [9] [10] 将最优控制问题转换为边值问题，利用Pontryaguin最大原理(PMP)显式或隐式确定的进行求解。间接法能够快速获得收敛解，但是对初始值灵敏度高。除此之外，近些年也有学者研究了混合方法，即努力将直接方法和间接方法结合起来 [11] [12]，首先使用直接方法来建立良好的初始猜测，然后使用间接方法来产生精确收敛。对于此问题，国内学者也进行了大量研究 [13] - [20]，主要围绕轨迹、总体参数优化设计及最新的基于伪谱法的最优控制求解问题。

上述的研究方法主要处理纯轨迹问题，初步设计研究旨在为参考任务寻找最佳的运载器配置。在工程实际中，要解决的任务一般为多学科问题，需要同时解决配置参数和轨迹优化问题。由于其整体复杂性，只能分步解决，包括在运载火箭构型和轨迹优化之间进行迭代，直到趋于可接受的设计。

这类问题的一个典型特点是高动态和不确定性，需要一定的预测能力和在线优化能力，同时对于各类约束条件有较好的适应性。基于此类特点，本文拟探索使用模型预测控制(MPC)来解决上升段轨迹的预测能力，并能够针对约束条件的变化，进行在线优化。首先，对动力学模型进行了推导，然后给出了模型预测控制的最有控制问题构建，并利用多重打靶法将其转换为非线性规划问题，再采用相应的求解器进行求解，最后将其应用到国外某型运载火箭队仿真结果进行了验证。结果表明该方法能够很好地处理约束条件，通过在线优化保证最终的最优控制目标的实现。

2. 控制模型

通用上升段运载火箭运动方程，在惯性系中的运动方程为：

$\begin{array}{l}\stackrel{\dot{}}{r}=v\\ \stackrel{\dot{}}{v}=g\left(r\right)+\frac{T}{m\left(t\right)}u+\frac{A}{m\left(t\right)}+\frac{N}{m\left(t\right)}\\ \stackrel{\dot{}}{m}=-\frac{T_{vac}}{g_0{I}_{sp}}\end{array}$ (1)

其中，r为惯性位置，v为速度矢量，g为重力加速度，T_vac为真空推力的幅度，T为当前的推力，A和N为轴向和法向的空气动力；u为推力矢量方向，m(t)为当前时刻的质量，I_sp和g₀为比冲和地球引力。

$\begin{array}{l}T=\left[\eta T_{vac}+\Delta T\left(r\right)\right]/m\left(t\right)g_0\\ V_r=V-{\varpi }_e\times r-V_w\\ A=A\left(h,\alpha ,M,v\right)=\frac{1}{2}\rho v_r^{2}S{C}_A\equiv qS{C}_A\\ N=N\left(h,\alpha ,M,v\right)=\frac{1}{2}\rho v_r^{2}S{C}_N\equiv qS{C}_N\end{array}$

$\begin{array}{l}\rho =\rho \left(h\right)={\rho }_0\exp \left({\rho }_{1}h+{\rho }_2{h}^2+{\rho }_3{h}^3\right)\\ a=a\left(h\right)=a_0+a_{1}h+a_2{h}^2+a_3{h}^3\\ C_A=C_{A0}\left(\alpha ,M\right)+C_{A1}\left(M\right)\alpha +C_{A2}\left(M\right){\alpha }^2\\ C_N=C_{N0}\left(\alpha ,M\right)+C_{N1}\left(M\right)\alpha +C_{N2}\left(M\right){\alpha }^2\end{array}$ (2)

引力模型采用J2的模型。

基于J2的模型：

$\begin{array}{l}g=\left[\begin{array}{c}{g}_x\\ g_y\\ g_z\end{array}\right]=\left[\begin{array}{c}{g}_r\frac{R_{0x}+X}{r}+g_{\omega }\cdot b_x\\ g_r\frac{R_{0y}+Y}{r}+g_{\omega }\cdot b_y\\ g_r\frac{R_{0z}+Z}{r}+g_{\omega }\cdot b_z\end{array}\right]\\ g_r=-\frac{\mu }{r^2}\left[1+\frac{3}{2}{J}_2{\left(\frac{a}{r}\right)}^2\cdot \left(1-5{\sin }^2{\phi }_x\right)\right]\\ g_{\omega }=-3J_2\left(\frac{\mu }{r^2}\right)\cdot {\left(\frac{a}{r}\right)}^2\cdot \sin {\phi }_x\end{array}$ (3)

$\begin{array}{l}{x}^{\prime }=x+R_{0x}\\ y^{\prime }=y+R_{0y}\\ z^{\prime }=z+R_{0z}\end{array}$

$\begin{array}{l}{\omega }_{ex}={\omega }_e\cos B_0\cos A_0\\ {\omega }_{ey}={\omega }_e\sin B_0\\ {\omega }_{ez}=-{\omega }_e\cos B_0\sin A_0\end{array}$

$\begin{array}{l}{V}_{ex}={\omega }_{ey}{z}^{\prime }-{\omega }_{ez}{y}^{\prime }\\ V_{ey}={\omega }_{ez}{x}^{\prime }-{\omega }_{ex}{z}^{\prime }\\ V_{ez}={\omega }_{ex}{y}^{\prime }-{\omega }_{ey}{x}^{\prime }\end{array}$

$\begin{array}{l}{\tilde{V}}_x=V_x-V_{ex}+f_x\\ {\tilde{V}}_y=V_y-V_{ey}+f_y\\ {\tilde{V}}_z=V_z-V_{ez}+f_z\end{array}$

$\begin{array}{l}{f}_x=-V_{wind}\cos \left(\left(A_{wind}-\pi \right)-A_0\right)\\ f_y=0\\ f_z=-V_{wind}\sin \left(\left(A_{wind}-\pi \right)-A_0\right)\\ {\tilde{V}}^2={\tilde{V}}_x^2+{\tilde{V}}_y^2+{\tilde{V}}_z^2\end{array}$ (4)

马赫数：

$\text{}M=\frac{\tilde{V}}{a_{SL}}\sqrt{\frac{T_0}{T_{SL}}}=\frac{\tilde{V}}{20.0468\sqrt{T_{SL}}}$

$q{S}_M=\frac{1}{2}{\rho }_{SL}\cdot \frac{\rho }{{\rho }_{SL}}\cdot {\tilde{V}}^2\cdot S_M=\frac{1}{2}\rho \cdot {\tilde{V}}^2\cdot S_M$

$\begin{array}{l}\sin \left(\alpha +{\alpha }_w\right)=-\frac{{\tilde{V}}_{y1}}{\tilde{V}}\\ \sin \left(\alpha +{\alpha }_w\right)=\frac{{\tilde{V}}_{z1}}{\sqrt{{\tilde{V}}_{x1}^2+{\tilde{V}}_{z1}^2}}\end{array}$ (5)

其中， $\left|\alpha +{\alpha }_w\right|\le \frac{\pi }{\text{2}}$， $\left|\beta +{\beta }_w\right|\le \frac{\pi }{\text{2}}$。

$\begin{array}{l}{\tilde{V}}_{x1}=a_{bi\_11}{\tilde{V}}_x+a_{bi\_21}{\tilde{V}}_y+a_{bi\_31}{\tilde{V}}_z\\ {\tilde{V}}_{y1}=a_{bi\_12}{\tilde{V}}_x+a_{bi\_22}{\tilde{V}}_y+a_{bi\_32}{\tilde{V}}_z\\ {\tilde{V}}_{z1}=a_{bi\_13}{\tilde{V}}_x+a_{bi\_23}{\tilde{V}}_y+a_{bi\_33}{\tilde{V}}_z\end{array}$

其终端条件包括位置幅度、飞行路径角及相应的轨道要素，可以写为：

$\begin{array}{l}r\left(t_f\right)=r_f\\ v\left(t_f\right)=v_f\\ r{\left(t_f\right)}^{\text{T}}v\left(t_f\right)-r_t{v}_t\sin \left({\gamma }_t\right)=0\end{array}$ (6)

飞行路径约束条件一般包括：动压、动压攻角乘积、动压侧滑角乘积，以及保证火箭的飞行高度和控制矢量的单位化。

$\begin{array}{l}\left|r\left(t_0\right)\right|\ge R_0\\ \left|u_i\left(t\right)\right|=1\\ \frac{1}{2}\rho v_r^2\le q_{\max }\\ \frac{\eta T_{vac}}{mg}\le a_{\max }\\ q\alpha \le {\left|q\alpha \right|}_{\max }\\ q\beta \le {\left|q\beta \right|}_{\max }\end{array}$ (7)

初始条件由发射点的坐标及火箭的初始条件确定的。

$\begin{array}{l}r\left(t_0\right)=r_0\\ v\left(t_0\right)=v_0\\ m\left(t_0\right)=m_0\end{array}$ (8)

终端条件：

终端条件保证速度、位置满足入轨条件。

$\begin{array}{l}a\left(t_f\right)=a_f\\ e\left(t_f\right)=e_f\\ i\left(t_f\right)=i_f\\ \Omega \left(t_f\right)={\Omega }_f\\ \omega \left(t_f\right)=\omega \end{array}$ (9)

以上的各参量分别为：半长轴、偏心率、轨道倾角、升交点经度及近地点幅。

3. 模型预测控制

模型预测控制(MPC)是一种实现反馈控制的方法，并提供了一种有效的通用适用于在线优化的设计方法，主要应用于多输入多输出(MIMO)系统的控制系统，其最大的问题时候计算成本较大。随着计算技术的进步，计算能力已经大大提高，计算硬件和软件的最新进展也大大降低了成本并改善了这些计算的速度和可靠性，从而使计算密集型控制系统的收益/成本比越来越大。以MPC为代表的控制系统应用也越来越普遍，大大增加了对复杂多输入多输出系统设计控制器的需求。典型的MPC系统如图1所示。

$x\left(k+1\right)=f\left(x\left(k\right),u\left(k\right)\right)$ (10)

在k时刻，获取x(k)，在此基础上，在预测区域N内计算最优控制序列u(x(k))，将其应用到采样周期 $\left[k,k+\text{1}\right]$，然后循环执行。

模型预测控制的三个主要过程：预测、在线优化和滚动时域执行。其本质上也是一类最优控制问题，其数学构成为：

控制目标：

$\mathcal{l}\left(x,u\right)={\Vert x_u-x^r\Vert }_Q^2+{\Vert u-u^r\Vert }_R^2$ (11)

代价函数：

$J_N\left(x,u\right)={\displaystyle \underset{k=0}{\overset{N-1}{\sum }}\mathcal{l}\left(x_u\left(k\right),u\left(k\right)\right)}$ (12)

最优控制构成：

$\begin{array}{l}\min J_N\left(x_0,u\right)={\displaystyle \underset{k=0}{\overset{N-1}{\sum }}\mathcal{l}\left(x_u\left(k\right),u\left(k\right)\right)}\\ x_u\left(k+1\right)=f\left(x_u\left(k\right),u\left(k\right)\right)\\ x_u\left(0\right)=x_0\\ u\left(k\right)\in U,x_u\left(k\right)\in X\end{array}$ (13)

Figure 1. MPC structure

图1. MPC结构图

4. NLP问题构建及求解方法

当前，最常用的将最优控制问题(OCP)转换为NLP问题的方法包括打靶法和配点法。打靶法包括单打靶法和多重打靶法。

4.1. 单打靶法

NLP决策变量：

$w=\left[\begin{array}{ccc}{u}_0& \cdots & u_{N-1}\end{array}\right]$ (14)

将x(u)作为w、x₀和t_k的函数：

$\begin{array}{l}{x}_u(.)=F\left(w,x_0,t_k\right)\\ F\left(w,x_0,t_0\right)=x_0\end{array}$ (15)

代价函数：

$\begin{array}{l}\min \Phi \left(F\left(w,x_0,t_k\right),w\right)\le 0\\ g_1\left(F\left(w,x_0,t_k\right),w\right)\le 0\end{array}$ (16)

4.2. 多重打靶法

直接多重打靶法是解决边值问题的一种数值方法，其算法结构图如图2所示。该方法将寻求解决方案的时间间隔分为几个较小的时间间隔，解决了每个较小时间间隔中的初始值问题，并施加了额外的匹配条件以在整个时间间隔上形成一个解决方案。与单打靶法相比，该方法在非线性分布和数值稳定性方面具有重大改进。

Figure 2. Structure of multiple shooting method

图2. 多重打靶法结构图

NLP决策变量：

$w=\left[\begin{array}{ccc}{u}_0& \cdots & u_{N-1},\begin{array}{ccc}{x}_0& \cdots & x_N\end{array}\end{array}\right]$ (17)

代价函数：

$\begin{array}{l}\min \Phi \left(F\left(w,x_0,t_k\right),w\right)\le 0\\ g_1\left(w\right)=\left[\begin{array}{c}{g}_1\left(x_0,u_0\right)\\ ⋮\\ g_1\left(x_{N-1},u_{N-1}\right)\\ g_1\left(x_N\right)\end{array}\right]\le 0\\ g_2\left(w\right)=\left[\begin{array}{c}{\bar{x}}_0-x_0\\ f\left(x_0,u_0\right)-x_1\\ ⋮\\ f\left(x_{N-1},u_{N-1}\right)-x_N\end{array}\right]=0\end{array}$ (18)

4.3. NLP求解算法

当前在处理NLP时，普遍采用了基于可行方向搜索的约束最优化方法，包括内部点信任区域线搜索算法和主动集信任区域线搜索算法。其中，具有代表性的是成熟于80年代中后期的序列二次规划算法(SQP)及内部点迭代算法(IPI)。

内部点技术实现了原始对偶内部点算法，使用障碍函数来确保算法可行。当优化问题包含许多不等式约束且大多数约束将作为严格的不等式满足时，内点法非常有用。如果只需要考虑很少的约束，就可以减小原始优化问题的大小和复杂性。

本文采用基于内点法的IPOPT对NLP问题进行求解。

5. 仿真校验

以国外某型两级运载火箭为模型对象，初始条件为：

初始高度H₀ = 10 m，初始速度V₀ = 0，初始质量M₀ = [20000 kg 1000 kg]，结构质量为[2000 kg 200 kg]，阻力系数为[0.3 0.2]，比冲为[300 350]。

仿真结果如图3~10：

Figure 3. MPC iterative solution process

图3. MPC迭代求解过程

Figure 4. Height result curve

图4. 高度结果曲线

Figure 5. Velocity result curve

图5. 速度结果曲线

Figure 6. Quality result curve

图6. 质量结果曲线

Figure 7. Acceleration result curve

图7. 加速度结果曲线

Figure 8. Thrust result curve

图8. 推力结果曲线

Figure 9. Voyage result curve

图9. 航程结果曲线

Figure 10. Dynamic pressure result curve

图10. 动压结果曲线

从图3可以看出，MPC在迭代过程中不断对目标函数进行优化，同时在线处理飞行约束条件，单条轨迹的收敛时间在9 s左右。从图7、图8可以看出，整条轨迹控制指令平滑，飞行最大过载较低。从图10可以看出，最大动压控制较好，飞行过程对约束条件的在线处理效果明显。

6. 结论

对于具有高度不确定性及严格飞行过程约束运载火箭上升段轨迹优化问题，本文利用MPC构建最优控制问题，该方法能够通过在线预测与优化能力，并且能够利用MPC的过程约束处理能力处理强约束条件。方针结果表明MPC在迭代过程中不断对目标函数进行优化，同时在线处理飞行约束条件，单条轨迹的收敛时间在9 s左右，飞行过程对约束条件的在线处理效果明显。

参考文献