助推-滑翔飞行器弹道最优控制研究*

范文锋，许波，郝昀

(北京机电工程总体设计部，北京 100854)

0 引言

近年来，随着临近空间的利用逐渐受到重视，一种助推-滑翔式飞行器(boost-glide vehicle, BGV)得到国内外学者广泛关注[1-6]。其特点是通过主动段的快速助推使飞行器达到较大速度和高度，随后在临近空间以无动力跳跃滑行的方式进行长时间机动飞行，具有可实现远程快速精确攻击、覆盖区域大、机动性能好及突防能力较强等优点。本文拟对助推-滑翔飞行器弹道最优控制问题进行研究。

目前求解此类弹道优化问题主要有间接法和直接法两大类[7]。间接法是根据Pontryagin极小值原理将约束最优控制问题转化为Hamiltonian边值问题(Hamiltonian boundary value problem, HBVP)进行求解，其特点是解的精度较高且满足一阶最优性必要条件，然而收敛半径小且协态变量初值难以给定等问题成为此类方法广泛应用的主要困难。直接法采用参数化方法将连续空间的最优控制问题离散转化为非线性规划(non-linear programming, NLP)问题，进而使用成熟的NLP求解方法得到最优弹道，由于直接法不必推导一阶最优必要条件，收敛半径大，无需协态变量的初始值，因此直接法比间接法应用更为广泛。

本文以助推-滑翔飞行器为研究对象，建立无量纲化的弹道动力学模型并提炼出总体设计参数，考虑飞行中的实际约束，以射程最优为目标建立BGV弹道最优控制模型；采用Radau伪谱方法将约束最优控制问题转化为非线性规划问题，通过引入连接点概念处理多阶段不连续问题，并使用序列二次规划方法(sequential quadratic programming, SQP)进行求解；最后给出数值优化算例，说明弹道最优控制设计特点以及文中方法的实用性和有效性。

1 弹道动力学模型

1.1 基本假设

本文以助推-滑翔飞行器为研究对象，仅限于均匀重力场中的纵向飞行弹道，考虑地球形状，忽略地球自转及扁率的影响。

1.2 气动力模型

借鉴文献[8]中拟合得到的高精度气动模型，其结果如下：

CD=0.000 68Ma4-0.014 1Ma3+0.110 21Ma2-

0.408 38Ma+0.822 73α+0.174 92，

(1)

CL=2.423 9×10-6Ma5-1.895×10-5Ma4+

1.863 910-5Ma3-0.000 48Ma2+0.003 94Ma-

0.000 93α2+α+69.656.

(2)

1.3 大气密度模型

假设大气密度遵循指数变化规律，其表达式为

ρ=ρ0exp-βH,

(3)

式中：ρ0为海平面大气密度，1.225 kg/m3；β=1/7 200。

1.4 无量纲化动力学模型

当发动机工作结束后开始被动段飞行，此时仅考虑动力学方程组中前4项即为被动段动力学方程组。

综上可知：助推-滑翔飞行器的主要设计参数可归纳为tb，μt，Kfg，Kgs，CD和CL，由于tb=Ispμt/Kfg，因此最大射程可表示为

Lmax=fIsp,μt,Kfg,Kgs,CD,CL.

(5)

2 弹道最优控制模型

2.1 约束条件

助推-滑翔飞行器飞行弹道经历严酷的力、热环境，需要考虑气动加热、结构承载以及姿态稳定的需求，提出如下过程约束以确保飞行正常。

驻点热流密度约束：

(6)

总过载约束：

(7)

动压约束：

(8)

此外，为保证对目标的打击效果，对落点状态提出如下终端约束：

Θtf=Θf,

(9)

(10)

2.2 BGV弹道最优控制问题

在以上工作基础上，本文研究的助推-滑翔飞行器弹道最优控制问题可以描述为求解飞行攻角α，使得如下目标函数最大：

(11)

且满足状态方程组(4)、过程约束(6)～(8)以及终端约束(9)～(10)。

3 Radau伪谱方法

伪谱法也称正交配点法，主要包括Legendre伪谱法(Legendre pseudo-spectral method, LPM)、Gauss伪谱法(Gauss pseudo-spectral method, GPM)以及Radau伪谱法(Radau pseudo-spectral method, RPM)。这些方法相同之处都采用Lagrange多项式对状态变量和控制变量进行全局逼近，区别在于配点的选取不同。文献[9]对以上3种伪谱方法特点进行了详细讨论，指出采用GPM 和RPM得到的NLP问题KKT(Karush-Kuhn-tucker)条件与原HBVP一阶最优性必要条件等价。由于RPM配点包含初始点，容易得到协态变量高精度初始值，本文采用RPM对助推-滑翔飞行器弹道最优控制问题进行求解。

3.1 基本原理

Radau伪谱法将约束最优控制问题的状态变量和控制变量在Legendre-Gauss-Radau(LGR)配点处进行离散，采用全局插值的Lagrange多项式对状态变量和控制变量进行全局逼近，通过对多项式求导来近似状态方程中的导数，且在一系列配点上满足控制方程中右函数的约束，从而将状态方程转化为代数形式的等式约束；目标函数中的积分项由Gauss-Radau积分近似计算；至此可将约束最优控制问题转化为非线性规划问题[10]。

3.2 连续最优控制问题的离散转换

3.2.1 连续最优控制问题的一般格式

不失一般性，考虑如下Bolza最优控制问题。求解控制函数ut及相应的状态函数xt，使得如下目标函数最小：

(12)

满足状态方程：

(13)

过程约束：

gxt,ut,t≤0，

(14)

边界约束：

ψxt0,t0,xtf,tf=0.

(15)

3.2.2 离散转换的主要内容

(1) 时域转换

Radau伪谱法的求解时域为 [-1,1]，因此需要将实际求解时域[t0,tf]进行转化：

(16)

(2) 状态变量与控制变量的近似

Radau伪谱法采用全局正交的Lagrange插值多项式分别对状态变量与控制变量进行逼近。

状态变量的近似公式为

(17)

(18)

式中：τjj=1,…,N为LGR配点，对应LN-1τ+LNτ的根，LNτ表示N阶Legendre多项式；τN+1不属于LGR配点但参与状态变量离散近似过程，其取值为1。

控制变量只在LGR配点处进行离散近似，其近似公式为

(19)

(20)

(3) 状态方程转换

将式(17)对τ进行微分，可得到LGR配点处的状态变量导数为

(21)

式中：Dkj为N×N+1阶Radau伪谱微分矩阵分量。

(22)

在此基础上，可将LGR配点处状态方程转化为代数形式的等式约束：

k=1,2，…,N.

(23)

(4) 积分项的近似

最优控制问题中的积分项采用Gauss-Radau积分进行处理。不失一般性，其转换格式为

(24)

式中：ωk为Gauss-Radau积分权重。

3.2.3 转换后的非线性规划问题

基于以上离散处理可将最优控制问题式(12)～(15)转换为如下形式：

J=HxNτ1,τ1,xNτN+1,τN+1+

(25)

满足约束：

k=1,2，…,N.

显然式(25)，(26)构成一个非线性规划问题，即求解离散点处的状态变量xNτk(k=1,2，…,N+1)、LGR配点处的控制变量uNτk(k=1,2，…,N)以及终端时刻tf(若未给定)，使式(25)在满足式(26)的约束下最小，可采用成熟的SQP方法求解。

3.3 多阶段不连续问题处理

对于助推-滑翔飞行器，主动段及被动段动力学模型存在差别，其弹道优化本质上是一个多阶段不连续最优控制问题。通过引入连接点的概念形成多阶段伪谱方法[11-12]，即根据不同阶段动力学模型及约束条件变化特点将最优控制问题划分成若干子区间，在每个子区间上分别对状态变量和控制变量进行离散配点，同时对相邻子区间连接点加入约束条件，以处理多阶段不连续的最优控制问题。

假设在某一点τ1∈[τ0,τ2]存在不连续，记

则连接点约束条件为

φ(x-(τ1),x+(τ1))≤0.

(27)

当τ1为不定时，可将其作为优化变量，同时加入不等式约束条件：

τ1-τ0τ1-τ2≤0.

(28)

4 数值优化算例

优化过程中根据计算精度需求及弹道特点选取配点数，其中主动助推段LGR配点数取75个，被动滑翔段LGR配点数取105个。总体设计参数取值如表1所示，过程约束及终端约束指标如表2所示。以助推-滑翔飞行器总射程最大为目标进行优化，结果如图1～7所示。

表1 总体设计参数值Table 1 Parameters for overall design

表2 约束指标Table 2 Specifications for constraints

图1 纵向平面弹道Fig.1 Planar trajectory

图2 当地弹道倾角Fig.2 Flight path angle

图3 飞行攻角Fig.3 Angle of attack

图4 总速度与动压Fig.4 Total velocity and dynamic pressure

优化结果表明：各状态及控制变量曲线过渡光滑，且满足过程及终端约束。由图3飞行攻角曲线可知：采用2次攻角转弯的策略进行主动段飞行程序设计，有利于精细调节主动段能量损失；被动滑翔段采用波浪式攻角变化规律，而非采用最大升阻比状态定常攻角飞行策略，更有利于提高射程。

由图6可知与飞行距离相关的协态变量值恒为-1，由图7可知Hamiltonian函数在主动助推段和被动滑翔段均保持为常值，基于助推-滑翔飞行器弹道最优控制模型特点，可知图6～7结果从侧面说明了采用的多阶段Radau伪谱方法得到的KKT条件与原HBVP一阶最优性必要条件等价，计算结果与理论最优解一致。

图5 驻点热流密度与总过载Fig.5 Heat flux at the stagnation point and total overload

图6 飞行距离相关的协态变量Fig.6 Co-state about flight range

图7 哈密顿函数Fig.7 Time history for Hamiltonian

5 结论

本文从多阶段多约束最优控制的角度对助推-滑翔飞行器弹道最优控制问题进行研究，得到结论如下：

(1) 通过无量纲化方法得到助推-滑翔飞行器无量纲化动力学方程组，可提取出Isp，μt，Kfg，Kgs，CD和CL作为BGV总体设计的主要参数。

(2) 建立了助推-滑翔飞行器弹道最优控制模型，采用Radau伪谱方法将最优控制问题转化为NLP问题并使用SQP方法求解；通过数值优化算例验证了模型及方法实用有效，优化结果满足一阶最优性必要条件。

(3) 助推-滑翔飞行器弹道包括主动助推段和被动滑翔段2部分，其弹道优化是典型的多阶段不连续最优控制问题；通过引入连接点的概念处理多阶段不连续问题，可避免以往弹道分段优化的缺陷，以全射程最优直接处理助推-滑翔飞行器弹道最优控制问题。

参考文献：

[1] 刘欣, 杨涛, 张青斌. 助推-滑翔导弹弹道优化与总体参数分析[J].弹道学报, 2012, 24(3): 43-48.

LIU Xin, YANG Tao, ZHANG Qing-bin. Trajectory Optimization and Parameter Analysis for Boost-Glide Missile[J]. Journal of Ballistics, 2012, 24(3): 43-48.

[2] 李瑜, 杨志红, 崔乃刚. 助推-滑翔导弹最大射程优化[J]. 弹道学报, 2008, 20(4): 53-56.

LI Yu, YANG Zhi-hong, CUI Nai-gang. Optimization of Maximum Range for Boost-Glide Missile[J]. Journal of Ballistics, 2008, 20(4): 53-56.

[3] 李珂, 聂万胜, 冯必鸣. 助推-滑翔飞行器弹道分段研究[J]. 指挥控制与仿真, 2012, 34(5): 21-25.

LI Ke, NIE Wan-sheng, FENG Bi-ming. Research on Multi-phase Trajectory Optimization for Boost-Glide Vehicle[J]. Command Control & Simulation, 2012, 34(5): 21-25.

[4] 刘欣, 杨涛. 滑翔导弹再入拉起段弹道优化与制导[J].国防科技大学学报, 2012, 34(1): 67-71.

LIU Xin, YANG Tao. Trajectory Optimization and Guidance in Reentry Phase for Glide Missile[J]. Journal of National University of Defense Technoloty, 2012, 34(1): 67-71.

[5] 王晨曦, 李新国. 助推-滑翔导弹射程管理技术研究[J]. 固体火箭技术, 35(2): 143-147.

WANG Chen-xi, LI Xin-guo.Range Management Technology Research of Boost Glide Missiles[J]. Journal of Solid Rocket Technology, 35(2): 143-147.

[6] 刘君, 陈克俊, 谢愈, 等. 助推滑翔飞行器发射诸元计算方法研究[J]. 弹箭与制导学报, 2012, 32(6): 33-36.

LIU Jun, CHEN Ke-jun, XIE Yu, et al. The Research on Computing Method for Firing Data of Boost-Glide Aerocraft[J]. Journal of Projectiles, Rocket, Missiles and Guidance, 2012, 32(6): 33-36.

[7] 杨秀霞, 张毅, 施建洪, 等. 助推-滑翔飞行器轨迹设计研究综述[J]. 海军航空工程学院学报, 2012, 27(3): 245-252.

YANG Xiu-xia, ZHANG Yi, SHI Jian-hong, et al. A Survey of Boost-Glide Aircraft Trajectory Design[J]. Journal of Naval Aeronautical and Astronautical University, 2012, 27(3): 245-252.

[8] 宗群, 田柏苓, 窦立谦. 基于Gauss伪谱法的临近空间飞行器上升段轨迹优化[J]. 宇航学报, 2010, 31(7): 1775-1781.

ZONG Qun, TIAN Bai-ling, DOU Li-qian. Ascent Phase Trajectory Optimization for Near Space Vehicle Based on Gauss Pseudo-Spectral Method[J]. Journal of Astronautics，2010, 31(7): 1775-1781.

[9] Divya Garg. Advances in Global Pseudo-Spectral Methods for Optimal Control [D].Gainesville,USA: University of Florida, 2011.

[10] HAN Peng, SHAN Jia-yuan, MENG Xiu-yun. Reentry Trajectory Optimization Using a Multiple-Interval Radau Pseudo-Spectral Method[J]. Journal of Beijing Institute of Technology, 2013, 22(1): 20-27.

[11] 杨希祥, 张为华. 基于Gauss伪谱法的固体火箭上升段轨迹快速优化研究[J]. 宇航学报, 2011, 32(1): 15-21.

YANG Xi-xiang, ZHANG Wei-hua. Rapid Optimization of Ascent Trajectory for Solid Launch Vehicles Based on Gauss Pseudo-Spectral Method[J]. Journal of Astronautics,2011, 32(1):15-21.

[12] 关成启, 陈聪. 基于Gauss伪谱法的助推-滑翔飞行器多阶段约束轨迹优化[J]. 宇航学报, 2010,31(11): 2512-2518.

GUAN Cheng-qi, CHEN Cong. Multiphase Path-Constrained Trajectory Optimization for the Boost-Glide Vehicle Via the Gauss Pseudo-Spectral Method[J]. Journal of Astronautics， 2010, 31(11): 2512-2518.