基于C/GMRES模型预测控制的滑翔段轨迹跟踪研究

闫 瑞，泮斌峰

（西北工业大学航天学院，西安710072）

1 引言

滑翔式高超声速飞行器作为一种典型的再入飞行器，具有速度快、航程远、机动性强等特点，其再入制导技术是保证飞行器完成飞行任务的关键技术。而滑翔式飞行器的再入飞行过程具有高度和速度变化范围大、总体约束苛刻、飞行环境变化剧烈且气动参数不确定等特点［1-2］，增大了制导律设计的难度。

再入制导方法一般可以分为标称轨迹跟踪制导和预测校正制导。标称轨迹制导法在制导系统工作前已规划好了飞行参考轨迹，计算简单，对弹载计算机的存储容量和计算速度要求较低，更容易实现；预测校正制导按落点预测方法的不同可分为解析预测制导法和数值预测制导法。前者在线计算量小，但是落点的预测精度不高；后者精度较高，但是数值积分计算量大，对机载计算机要求高［3］。因此研究标称轨迹跟踪制导具有较大的应用价值。

标称轨迹跟踪制导中的轨迹跟踪方法有：反馈线性化、线性二次型调节器（Linear Quadratic Regulator，LQR）、模型预测控制［4］等。Bharadwaj等［5］采用近似反馈线性化方法，将三维轨迹跟踪转化为对横程和纵程的跟踪，反馈线性化方法在模型存在较大的不确定性，非线性项无法精确对消，控制性能将变差甚至失稳；Dukeman［6］基于LQR方法设计轨迹跟踪控制律，根据二次型性能指标离线计算最佳反馈系数，虽然LQR控制律对初始条件的拉偏具有良好的鲁棒性，但是控制效果和精度取决于参考轨迹的特征，当轨迹陡峭时，精确度将显著降低；Lu［7］首先提出采用非线性模型预测控制完成轨迹跟踪，与LQR方法相比，模型预测控制权重矩阵参数选择较不敏感，该方法可以准确地追踪三维轨迹并可保证闭环稳定性［8］。

非线性模型预测控制在每个采样点通过求解有限时域最优控制问题获得最优控制律。当采用间接法在线求解最优控制律时，需利用极小值原理将有限时域最优控制问题转化为两点边值问题并求解，但是Riccati矩阵微分方程的传统求解方法计算量较大，耗时较长［9］。

针对求解控制律时计算量大的问题，本文采用基于C/GMRES算法［10-12］的非线性模型预测控制方法设计合理的跟踪控制律对纵向参考轨迹进行跟踪控制：首先采用延拓法得到用于更新控制序列的微分方程组，使用广义极小残差法（Generalized Minimal Residual Methods，GMRES）高效解算近似线性方程，替代求解Riccati微分方程组，从而降低求解的复杂度，提高计算效率，保证了制导系统的实时性；然后基于航向角误差走廊进行横向轨迹控制；最后进行了数值仿真验证。

2 滑翔段纵向运动方程

设再入飞行器滑翔段纵向运动方程的状态变量X=[ r Vγ]T，控制量u=σ，以剩余航程为自变量的高超声速滑翔式再入飞行器纵向再入平面无量纲运动方程［13］（不考虑地球自转）可表示为式（1）：

式中：

其中，r为无量纲的飞行器到地球中心的距离，无量纲化参数为地球平均半径R0=6 378 137m；V为无量纲的飞行器相对地球的速度值，无量纲化参数为为海平面重力加速度值；γ为当地速度倾角；n为有量纲的飞行器质量；σ为倾侧角，L和D分别为无量纲的升力和阻力，表示为式（2）：

其中，ρ为有量纲的空气密度，Sref为有量纲的飞行器参考面积，CL和CD分别为升力系数和阻力系数。

3 滑翔段参考轨迹跟踪控制

滑翔段轨迹跟踪制导分为纵向制导和横向制导两部分，本节设计控制律来跟踪滑翔段纵向参考轨迹，由于纵向运动只与倾侧角σ的大小有关，因此将控制变量改写为并对其进行跟踪控制；横向运动方面，采用倾侧反转逻辑确定倾侧角的符号，将航向角误差控制在合理范围内，实现横向轨迹控制。

3.1 纵向轨迹跟踪控制律设计

采用基于C/GMRES算法的非线性模型预测控制方法来跟踪纵向参考轨迹。非线性模型预测控制问题本质上是一类有限时域最优控制问题，最小化性能指标如式（3）所示：

其中Δs为预测时域的长度，Xref为参考轨迹状态量，uref为参考轨迹倾侧角绝对值，Q和R为权重系数。在每一剩余航程点s处，最小化预测时域τ∈[ s，s+Δs]内的性能指标J，同时满足式（1）所示等式约束，以得到最优控制输入u*(τ)，只将初始值u*(s)作为剩余航程为s时的实际控制输入u( s)。在下一航程点处，系统重新以新的航程点为初始状态，求解最优控制输入，并继续将其初始值施加给受控对象，如此循环，直至完成整个控制过程。

定义哈密尔顿函数如式（4）所示：

其中，λ为协态向量。根据极小值原理，上述以剩余航程s处的实际状态量X( s) 为初始状态的最优控制问题可以转化为两点边值问题，其中正则方程表示为式（5）：

横截条件为式（6）：

极值条件为式（7）：

令采样周期为Δl，将上述连续两点边值问题离散化，得到离散两点边值问题，如式（8）所示：

定义式（9）

在每一剩余航程s处求解方程组（10）得到最优控制序列U( s)，实际控制输入u( s) 取最优控制序列的第一个元素( s)。

采用延拓法求解方程组（10），将其重写为式（11）：

其中As为使F=0的稳定矩阵。若FU非奇异，则可得到U( s) 的微分方程，如式（12）所示：

U′(s) 为控制量相对剩余航程的微分。由于式（12）中仍包含计算复杂的雅可比矩阵FX，FU和与相关的线性方程，为了减小求解雅可比矩阵和线性方程的计算量，采用雅可比矩阵和向量乘积的前向差分来近似F′，将式（11）近似为线性方程组，如式（13）所示：

其中，h为正实数。式（13）可以等价为式（14）：

前向差分近似和GMRES算法组合应用于求解关于U′的近似线性方程（14），称为FDGMRES算法［14］，表示为FDGMRES（U，X，X′，′，h，kmax），其中表示一个特定的初始猜 测，kmax表示最大迭代次数。表1中给出了算法的详细步骤

表1 FDGMRES算法Table 1 FDGMRES algorithm

通过FDGMRES算法计算得到U′，实时积分U′得到U( s)。用下缀n表示在剩余航程s0+nΔl处的相应变量，s0表示剩余总航程。将FDGMRES算法与上文中基于极小值原理的延拓法相结合，得到C/GMRES算法，其详细过程如表2所述。

应当注意的是，迭代法仅用于求解关于U′的线性方程式（14），并且通过积分U′，得到非线性方程组F( U，X) =0的解U，而无需任何线性搜索或牛顿迭代求解。C/GMRES算法在每个采样点仅求解一次线性方程式（14），因此与牛顿法等迭代方法多次求解线性方程以确定搜索方向相比，C/GMRES算法所需的计算量要少得多。此外，C/GMRES算法不涉及线性搜索，这也与标准优化方法有显著差异。

表2 C/GMRES算法Table 2 C/GMRES algorithm

3.2 基于航向角误差走廊的横向轨迹控制

飞行器飞行过程中航向误差角Δψ应一直在一个小范围内变化，使飞行器始终保持向目标区域飞行。本文根据吕瑞［15］研究所采用的倾侧反转策略和航向角误差走廊设计思路，离线规划航向角误差走廊以控制倾侧角反转，实现对飞行器飞行方向的控制。

定义航向角误差为导弹当前航向角与导弹相对目标的视线角之差，如式（15）所示：

式中，导弹当前位置相对目标点 ( θT，φT) 的视线角 ψLOS为式（16）：

采用漏斗形航向角误差走廊的设计思想［15］，在再入的初始阶段将误差控制在较小范围内，来减小制导系统的工作负荷；中间段适当放大误差走廊，将其设计为速度的分段线性函数；在终端将误差缩小到终端约束要求范围内。这种设计思路既可保证终端航向角的制导精度，也可使倾侧角的反转不过于频繁。

对于特定飞行器，可通过多次仿真实验获得适当的航向角误差走廊边界值。最终的航向角误差走廊设计为图1所示分段函数。

基于航向角误差走廊的反转逻辑的数学表达式为式（17）：

式中，sign( σi-1)为制导系统前一工作周期倾侧角的符号。ΔψdH为航向角误差走廊边界的绝对值。该反转策略表示当导弹的航向角误差的绝对值大于ΔψdH时，改变倾侧角σ的符号；在边界区域内时，σ保持原来的符号，不进行反转；在滑翔段结束时，航向角误差将符合终端约束。

图1 航向角误差走廊Fig.1 Heading angle error corridor

4 仿真验证与分析

4.1 无扰动的轨迹跟踪控制

以滑翔式再入飞行器X-33为例进行仿真验证，该飞行器质量为35 828 kg，参考面积为149.388 1m2，参考轨迹来自基于凸优化方法的三维再入轨迹规划［16］。滑翔段初始状态为：轨道地心距为6 452 735.38 m，速度为7528.77 m/s，速度倾角为-1.328 4°，经度为237.58°，纬度为-9.18°，航向角为2°；终端条件：轨道地心距为6 408 761.57 m，速度为1133.37 m/s，经度为238.25°，纬度为10.01°，航向角误差的绝对值小于等于6°。剩余总航程为2138.589 km。对状态量进行无量纲化，相应的标准值如第2节所述。设置仿真所需参数：采样周期Δl为50 m，预测区间长度 Δs为250Δl，h为1E-7，kmax为1000。考虑以跟踪高度为主，依据Bryson准则选取Q、R分别为diag（［3E10，0，0］）、10，选取制导指令更新周期为1 s，仿真步长为0.1 s。

图2～图5分别给出了滑翔段纵向轨迹跟踪图像、高度误差、倾侧角以及航向角误差随时间变化过程。

图2 纵向轨迹跟踪图Fig.2 Tracking effect of longitudinal reference trajectory

图3 高度误差变化曲线Fig.3 A ltitude error history

图4 速度倾角变化曲线Fig.4 Velocity path angle history

由图2可知，基于C/GMRES算法的非线性模型预测控制律跟踪效果较好；图3中高度误差随时间在约（-80 m，60 m）范围内变化，并且误差逐渐趋近于0，高度跟踪效果较好；图4给出了速度倾角随时间变化曲线，速度倾角在约（-8°，1°）范围内变化，再入飞行器滑翔段倾侧角反转次数为3～6次比较合理［1，4，17］；图5中倾侧角反转5次，反转次数较少，保证了制导系统的稳定性；图6中航向角误差在误差走廊范围内合理变化，滑翔段结束时航向角误差Δψf为-5.75°，符合终端约束要求，横向制导精度较高。

图5 倾侧角变化曲线Fig.5 Bank angle history

图6 航向角误差变化曲线Fig.6 Heading error history

4.2 有扰动的轨迹跟踪控制

进一步进行含有初始状态扰动和动态干扰的蒙特卡洛仿真，检验算法鲁棒性。滑翔式飞行器受到的控制力由气动力提供，所以气动力的扰动对再入制导的精度有重要的影响。气动力扰动主要可以分为大气密度扰动和气动参数扰动，且大气密度扰动可等效为气动参数发生扰动［3］，因此改变动力学方程中的气动系数进行仿真验证，采用的偏差量如表3所示。在含有初始状态扰动和动态干扰的条件下进行1000次蒙特卡洛打靶仿真，选取其中100次仿真结果得到高度和航向角偏差随时间变化图像如图7、8所示，1000次打靶的终端高度及航向角误差如图9、10所示。

表3 打靶加入的偏差Table 3 Perturbations for shoot

图7 高度误差变化曲线Fig.7 Altitude error histories

图8 航向角误差变化曲线Fig.8 Heading error histories

图9 终端高度误差Fig.9 Term inal altitude error

图10 终端航向角误差Fig.10 Term inal heading error

预先设定打靶落点高度误差合理范围为［-200 m，200 m］。图7中高度偏差随着时间逐渐减小，最后稳定在约（-200 m，150 m）范围内，符合终端误差要求；图8中航向角误差均在航向角误差走廊范围内变化，满足制导要求；图9中打靶落点在（-200 m，150 m）范围内变化，分布范围较小且均匀，满足终端误差要求；图10中终端航向角误差均在（-6°，-6°）范围内变化，符合终端航向角误差要求。以上仿真分析表明所设计的制导律对初始状态偏差、动态干扰鲁棒性较好，满足制导要求。

为了验证该控制律的快速性，本文在2 GB内存和CPU主频为2.8 GHz的计算机上以MATLAB R2018a为仿真平台，对制导律的计算时间进行统计。以1000次蒙特卡洛仿真为例，分别统计每次打靶中的最大制导计算时间和平均制导计算时间，如图11、图12所示。

图11 最大制导计算时间直方图Fig.11 Histogram of maximum guidance calculation time

图12 平均制导计算时间直方图Fig.12 Histogram of average guidance calculation tim e

图11中最大制导计算时间分布在（0.03 s，0.28 s）范围内，并且在（0.05 s，0.06 s）范围内分布最多，远小于制导周期1 s；图12中平均制导计算时间分布在（0.022 s，0.051 s）范围内，远小于制导周期1 s，且在（0.024 s，0.025 s）范围内分布最多。另外，采用信赖域算法求解上述非线性方程组（10）得到最优控制律的平均制导计算时间约为3.892 7 s。采用C/GMRES算法求解制导律的平均制导计算时间约为采用信赖域算法求解计算时间的1%左右。以上统计比较充分地证明了基于C/GMRES算法的非线性模型预测控制律的快速性和实时性。

5 结论

本文采用基于C/GMRES算法的非线性模型预测控制律对再入飞行器滑翔段参考轨迹进行跟踪控制。

1）通过仿真分析验证了该控制方法对初始状态扰动和动态干扰等不确定因素具有较强的鲁棒性。

2）1000 次蒙特卡洛仿真的最大制导计算时间在（0.05 s，0.06 s）范围内分布最多，平均制导计算时间都分布在（0.022 s，0.051 s）范围内，远小于制导周期1 s，表明求解该控制律所需计算量小，计算效率高且易于实现，能够满足制导实时性要求。