高超末段机动突防/精确打击弹道建模与优化

李静琳， 陈万春， 闵昌万

(北京航空航天大学 宇航学院, 北京 100083)

高超声速飞行器是一种全球进攻性武器，能进行远距离、防区外作战，最终在末段实现高精度、大威力的定点打击。然而，多层次反导防御系统尤其是末段爱国者-3(PAC-3)低层防御系统的发展，对高超声速飞行器的生存造成很大的威胁。因此，有必要研究末段机动突防和精确打击问题。

再入弹头一般利用反导武器拦截时间短、过载易饱和的弱点，以程序机动的方式盲目突防，如蛇形机动、螺旋机动等[1-3]。但是其突防效果严重依赖于对反导拦截时机的准确把握，而且无法对末段实施精确打击。同时，由于高超声速飞行器末段具有极高的飞行速度，并处于极其恶劣的力、热环境，这些都严重限制了其大机动能力的发挥。为了获得一种最优机动规律，在精确打击的前提下，最大限度发挥导弹的机动性能、提高其突防能力，有必要从最优控制的角度研究该问题[4-10]。Imado和Uehara[8-9]对飞机的最优躲避机动进行了深入研究，采用梯度法求解获得了适用于低速飞行器的最优机动弹道。但是该方法不能求解复杂的多段问题，只能对弹目交战阶段进行优化，无法考虑包含命中精度在内的多种终端约束。

Radau多段伪谱法(Multiphase Radau Pseudospectral Method，MRPM)[11]是一种求解多段、多约束最优控制问题的配点法，广泛应用于弹道优化、轨道转移等航空航天实际问题[12-14]，具有收敛性好、求解精度高的优点。因此，基于MRPM，本文以高超声速飞行器CAV(Common Aero Vehicle)为研究对象，考虑末段机动突防和精确打击的需求，提出了一种最优机动突防弹道优化方法，得到一种适合于高超声速飞行器末段机动的规律，从而获得其最大机动能力。首先，同时对突防弹、拦截弹的运动进行建模，通过施加约束限制拦截弹按照比例导引律飞行；根据飞行任务和双方弹道特征的变化，将飞行阶段分为弹目交战段、拦截弹指令饱和段和对地打击段；结合各段的任务和特性，分别提出了突防性能指标和精确打击性能指标等，并以加权函数的形式将各段独立、矛盾的性能指标统一，最终形成多对象、多段、多约束的末段机动弹道优化问题。然后，采用MRPM求解，通过积分弹道作初值、约束渐强以及分段嵌套优化等一系列优化策略，解决了复杂高超弹道优化问题的初值敏感、可行域窄等问题，获得了高精度最优解。根据协态映射原理，通过NLP(Nonlinear Programming Problem)问题的KKT(Karush-Kuhn-Tucker)乘子准确估计了最优控制问题的协态，求解了哈密顿函数，证明了优化结果的最优性。最后，通过多次攻防对抗仿真，开展最优突防弹道相对拦截弹发射时间的灵敏度研究。

1 高超末段机动突防与对地打击

高超声速飞行器一般是指在大气层内飞行、速度超过5倍声速的飞行器。CAV是一种典型的高超声速飞行器，其末段速度高达5倍声速，可以实现高速、大威力对地打击。但是由于拦截弹的存在(如PAC-3)，高超声速飞行器必须进行机动突防，如图1所示。CAV采用升力体外形，升阻比较高，在大气层中飞行时具备良好的机动能力。因此，本文的研究旨在寻找一种最优的机动突防弹道，在保证终端对地打击精度的前提下，最大限度地发挥高超声速飞行器的机动能力。

图1 高超末段机动突防-对地打击示意图

1.1 动力学与运动学建模

建立攻防双方动力学模型是开展最优机动突防弹道优化研究的第一步。由于高超声速飞行器末段具有拦截时间短、拦截高度低等特点，可以忽略地球自转和地球曲率的影响；并采用“瞬时平衡”假设，不考虑转动惯量；无推力偏心，无侧滑。

于是，突防弹在地面坐标系下的质点动力学方程可以表示为

(1)

式中：xA、yA、zA为飞行器在地面坐标系下的位置；VA为飞行器的飞行速度；θA和ψA分别为飞行器的弹道倾角和航向角；σA为飞行器的倾侧角；mA为飞行器的质量；g为当地地球引力项；LA和DA分别为飞行器的升力和阻力。

同理，拦截弹在地面坐标系下的质点动力学方程可以表示为

(2)

式中：xD、yD、zD为拦截弹在地面坐标系下的位置；VD为拦截弹的飞行速度；θD和ψD分别为拦截弹的弹道倾角和航向角；σD为拦截弹的倾侧角；mD为拦截弹的质量；PD为拦截弹发动机提供的推力；LD和DD分别为拦截弹的升力和阻力；αD为拦截弹的攻角。

1.2 攻防对抗相对运动建模

在攻防双方动力学建模的基础上，可以得到二者相对位置的微分方程

(3)

拦截弹按照纯比例导引律对突防弹进行拦截。则根据纯比例导引律，指令加速度为

(4)

(5)

2 多对象、多段、多约束最优控制问题

2.1 最优控制问题分段

本文研究的是一个典型的目标-突防-拦截 (Target-Attack-Defend, TAD)问题[15]，包含地面目标、突防弹、拦截弹3个对象，根据突防弹飞行任务和飞行弹道特征，主要分为以下几段：

第1段(T0～T1)为弹目交战阶段，突防弹为了躲避拦截进行机动飞行，拦截弹的运动遵循纯比例导引律。

第2段(T1～T2)为拦截弹指令饱和阶段。由于目标快速机动且二者相对距离过小，造成视线角速率过大，拦截弹指令加速度饱和，拦截弹保持恒定的最大或最小攻角飞行。该阶段飞行时间极短，攻防双方运动不再受制导律约束。

第3段(T2～T3)为弹头成功突防后对地打击阶段，不再考虑拦截弹的运动。

弹道分段如图2所示。

通过分段，可以使得不同任务和性质的弹道相对独立，避免相互之间耦合影响，便于分别研究，摸清并尽可能发挥各段的能力。

图2 弹道分段示意图

2.2 变量选取

本文虽然研究的是突防弹道的运动形式，但为了充分考虑拦截弹运动的影响，准确反映攻防对抗情况，同时对攻防双方建模。需要说明的是，本文认为拦截弹的运动规律是确定的，严格按照比例导引律对突防弹进行拦截。但是为了方便利用MRPM直接对攻防双方的运动状态进行求解，避免数值积分求解拦截弹道，缩减计算规模，此处同时以突防弹和拦截弹的状态量

(6)

作为状态变量，突防弹和拦截弹的攻角、倾侧角变化率

(7)

作为控制变量。式中：αA为突防弹的攻角。

2.3 系统状态方程

系统的状态方程包含突防弹和拦截弹的动力学方程，在1.1节已进行了详细介绍，即式(1)、式(2)，此处不再赘述。

2.4 系统多约束设置

2.4.1 初始和终端约束

突防弹的初始状态为中、末交班状态，拦截弹的初始状态为地面发射点火时的状态，假设双方初始状态确定，该约束可以表示为

(8)

(9)

式中：xA0、yA0、zA0为突防弹的初始位置；VA0、θA0、ψA0、αA0和σA0分别为突防弹的初始速度、弹道倾角、航向角、攻角和倾侧角；xD0、yD0、zD0为拦截弹的初始位置；VD0、θD0、ψD0、αD0和σD0分别为拦截弹的初始速度、弹道倾角、航向角、攻角和倾侧角；上标“*”表示期望值。

为了实现对地面目标进行精确打击，对突防弹的终端位置进行约束

(10)

式中： (xA3f,yA3f,zA3f)为突防弹道终端位置坐标;(xT,yT,zT)为地面目标的位置。

为了保证对落点的杀伤效果，同时避免速度过大造成飞行器解体，对突防弹的终端落角、落速和攻角进行严格限制：

(11)

式中：θA3f、VA3f和αA3f分别为突防弹终端的落角、落速和攻角；θmax、Vmax和αmax分别为最大终端落角、落速和攻角；θmin、Vmin和αmin分别为最小终端落角、落速和攻角。

2.4.2 路径约束

高超声速飞行器再入末段飞行环境恶劣，典型的过程约束包括驻点热流密度约束、动压约束和过载约束，其具体表达式为

(12)

2.4.3 连接条件

对于多段最优控制问题，为了保证各段之间的连续性，段与段之间需要施加连接条件，即

Xif=X(i+1)0

(13)

式中：Xif=(xif,yif,zif,Vif,θif,ψif,αif,σif)表示第i段的终端状态；X(i+1)0=(x(i+1)0,y(i+1)0,z(i+1)0,V(i+1)0,θ(i+1)0,ψ(i+1)0,α(i+1)0,σ(i+1)0)表示第(i+1)段的初始状态。

2.4.4 预测命中点耦合约束

为了保证拦截弹的指令饱和段在弹目相对距离最近时刻结束，需要对指令饱和段的终端状态设置如下约束条件：

rAD·vAD=0

(14)

式中：rAD和vAD分别为该时刻拦截弹和突防弹的相对位置和速度矢量。

2.4.5 拦截弹运动规律约束

在拦截弹和突防弹正常交战过程中，假设拦截弹按照纯比例导引律对突防弹进行拦截，根据三维空间纯比例导引律可知，弹道坐标系下，拦截弹的制导指令表示为

(15)

式中：ncy和ncz为指令加速度；Δψ为拦截弹速度与弹目相对速度在地面投影的夹角。根据拦截弹的运动方程式(2)，可知其当前实际的加速度为

(16)

将指令加速度与实际加速度严格相等作为过程约束引入优化模型中，即

(17)

以限制拦截弹按照纯比例导引律运动。

2.5 综合性能指标选取

高超声速飞行器末段飞行的主要任务是机动突防和精确打击，二者相互矛盾却缺一不可。性能指标的选取，决定了最优弹道的能力和特性。

2.5.1 间接性能指标

在以往的突防弹道优化研究中，没有考虑拦截弹的动力学模型，无法获得脱靶量信息，因此只能采用间接的方式来衡量突防效果。一般认为，增大突防弹的速度和减小突防飞行时间是提高突防能力的有效手段，典型的性能指标如下：

JJ1=-VA2f

(18)

JJ2=t2f

(19)

式中：VA2f为拦截时刻突防弹头的飞行速度；t2f为拦截时刻。

2.5.2 突防性能指标

脱靶量是对突防/拦截效果最直接的评估指标。脱靶量越大，突防效果越好。本文同时建立了攻防双方的模型，根据指令饱和段的终端状态，可直接得到准确的拦截脱靶量，以此作为反映弹道突防性能的指标

JT=-DMiss

(20)

式中：

DMiss=[(xA2f-xD2f)2+(yA2f-yD2f)2+

(21)

其中：下标“2f”表示指令饱和段终端时刻即拦截时刻。

2.5.3 精确打击性能指标

为了实现对地精确打击，式(10)对末段弹道的终端位置进行了等式约束。但是如果为了满足打击精度要求，设置严格的等式约束，很难直接得到满足要求的可行解。因此，将终端位置约束转化为性能指标：

JQ1=(xA3f-xT)2+(yA3f-yT)2+(zA3f-zT)2

(22)

同理，为了控制终端落角、落速，将其转化为性能指标：

JQ2=(θA3f-θf)2

(23)

JQ3=(VA3f-Vf)2

(24)

式中：θf、Vf分别为期望的落角、落速。

2.5.4 控制平滑性能指标

虽然本文优化的目的是得到同时具备突防能力和精确打击能力的机动弹道，但是仅仅以上述性能指标直接优化，得到的最优控制多数是存在突变、振荡的，如图3所示。显然平滑的控制更具有实际意义。因此，为了使优化得到的控制在一定程度上尽量平滑，需要施加如下性能指标：

(25)

式中:k1、k2为权重系数。

图3 平滑性能指标对控制曲线的影响

2.5.5 综合性能指标

为了得到兼顾机动突防和精确打击的末段弹道，采用加权函数的形式对上述几种性能指标进行组合，得到如下几种综合性能指标：

(26)

(27)

(28)

k4[(xA3f-xT)2+(yA3f-yT)2+(zA3f-zT)2]

(29)

式中:k1、k2、k3、k4均为权重系数。

图4 弹道优化流程图

3 Radau多段伪谱法求解多对象、多段、多约束最优控制问题

为了得到末段最优机动弹道，采用MRPM求解该多对象、多段、多约束最优控制问题。MRPM是一种处理最优控制问题的离散优化方法。其主要思想是将状态和控制在一系列正交节点即Legendre-Gauss-Radau(LGR)节点上离散并满足状态微分方程，并以这些离散点为节点构造拉格朗日插值多项式来逼近状态和控制变量，从而将最优控制问题转化为参数优化问题，即NLP问题。而且MRPM适合求解多段问题，具有高精度、指数收敛性质。

3.1 弹道优化策略

3.1.1 离散节点初值按弹道积分赋值

高超声速弹道具有强非线性、快速变化的特点，且对控制变量高度敏感。GPOPS初值生成器是在给定初始、终端状态的初值后，通过线性插值获得离散节点上的初值，这很难满足动力学方程的约束，使得问题不能收敛到可行解，更无法获得最优解。

因此，考虑再入动力学模型，本文以常值控制下的积分弹道作为突防弹的状态初值，采用比例导引律对该积分弹道进行拦截，并以此拦截弹道作为拦截弹的状态初值，从而同时满足双方的动力学约束，以及二者之间拦截导引律的耦合约束。而且，常值攻角、倾侧角控制下的弹道较为平滑，容易满足热流密度、动压等过程约束。如图5和图6所示，在相同情况下，与GPOPS自动生成初值相比，以积分弹道作为初值，解的可行性精度和最优性精度更高，收敛速度更快。其中，可行性精度和最优性精度分别反映了优化结果对约束条件和KKT条件的满足情况。

图5 可行性精度对比

图6 最优性精度对比

3.1.2 约束渐强

高超声速飞行器飞行速度快，在末段面临严峻的热、力学环境，热流密度、动压、过载等各种复杂约束使得弹道的可行域限制在较为狭窄的范围内，尤其对于多对象问题，还应当满足耦合约束即拦截导引律的限制。然而优化初值一般很难满足全部约束，仅通过一次优化很难收敛到狭窄的可行域中，通常得不到高精度的可行解。因此，需要先在一个扩展的可行域内得到可行解，并由该可行解迭代替换为优化初值，再由少到多逐步加强约束，缩小允许的可行域，直到得到满足约束的高精度解。如图7所示，强动压约束下，只能得到低精度可行解。而图8中，逐渐加强动压约束，3次迭代优化即可将可行解的精度提高2个数量级。

图7 可行性精度收敛性

图8 逐次迭代可行性精度收敛性

3.1.3 分段嵌套优化

本文研究的是一个复杂的多段优化问题，每一段都需要满足动力学方程约束、过程约束，而且段与段之间需要通过连接条件即多个严格的等式约束来保证解的连续性，如式(13)所示。若直接对3段同时优化，很难得到同时满足复杂多约束的可行解。因此，采用分段优化策略。先对第1段即弹目交战段优化得到最优解后，将第1段的优化结果作为初值、第1段的终端状态作为第2段即指令饱和段的初始状态，对前2段同时进行优化，得到前2段的最优解后再以此为初值，对3段同时进行优化，获得满足多段约束的解。

综上所述，针对该问题的一系列初值生成策略可以由图9表示。

3.2 协态映射原理及最优性验证

图9 弹道分段嵌套优化流程

通过MRPM的求解原理可知，本文实际上是将最优控制问题转化为NLP问题求解。对于NLP问题，KKT条件是确保其最优解的充要条件，但是这不能保证该解为原最优控制问题的最优解。

本文研究的是一个终端时间自由的最优控制问题，根据最优控制理论，其哈密顿函数为

H[x1(t),u1(t),x2(t),u2(t),λ1(t),λ2(t),t]=

(30)

即

H[x1(t),u1(t),x2(t),u2(t),λ1(t),λ2(t),t]=

λ1yVAsinθA+λ1z(-VAcosθAsinψA)+

λ2yVDsinθD+λ2z(-VDcosθDsinψD)+

(31)

式中：λ1x、λ1y、λ1z、λ1V、λ1θ、λ1ψ、λ2x、λ2y、λ2z、λ2V、λ2θ和λ2ψ为协态变量。由于哈密顿函数不显含时间t，且终端时刻自由，则有

(32)

即沿最优曲线哈密顿函数恒等于常数。因此，通过哈密顿函数的取值，可以验证优化结果的最优性。

根据MRPM的协态映射原理[16]，通过NLP问题的KKT条件可以准确地估计最优控制问题的协态：

(33)

(34)

(35)

(36)

通过该映射，NLP问题的KKT条件和离散的最优控制问题的一阶必要条件是等价的。通过NLP问题的拉格朗日乘子，可以得到最优控制问题的协态变量，从而引入最优控制问题的哈密顿函数验证MRPM得到的优化结果的最优性。

4 算例仿真

4.1 优化结果

为了验证最优控制问题建模的可行性，研究不同性能指标下最优机动弹道，获得高超飞行器的最大机动能力，以CAV和PAC-3为研究对象进行末段攻防对抗仿真。假设拦截弹的制导系统是理想的，即制导系统产生的制导指令，弹体能迅速进行跟踪，不会造成控制系统的偏差。

飞行器的初始、终端状态约束如表1所示。同时突防弹以攻角10°，倾侧角-180°为基准飞行状态。拦截弹的初始状态为发射后进行程序转弯前的状态，如表2所示。

末段需要严格满足热流密度、动压和过载约束，如表3所示。

攻防双方均采用BTT控制，且攻角、倾侧角及其变化率需满足以下约束条件：

表1 突防弹初始及终端状态约束

表2 拦截弹初始状态约束

表3 突防弹常规路径约束

(37)

通过仿真，不同性能指标下的最优弹道如图10和图11所示。终端状态如表4所示。其中，在x、y、z3个方向的落点误差均在1 m以内，远远小于CAV的毁伤半径，满足精确打击的任务需要。

根据文献[1]显示，PAC-3的毁伤半径为8 m，表5为常值控制和几种不同性能指标下突防弹道的拦截脱靶量。从表5可知，当飞行器按照常值攻角和倾侧角飞行不进行规避机动时，不能成功突防，但是，其他优化所得的机动弹道均成功突防，证明了动力学建模和最优控制问题建模的可行性。而且直接以脱靶量即J3、J4为性能指标所得到的最优弹道，其脱靶量比其他弹道高出1～2个数量级。从控制曲线图12可知，直接以脱靶量为指标优化，得到了一种纵横向组合机动控制规律，发挥了飞行器的最大机动能力，尤其是横向机动能力。

因此，为了验证最优纵横向组合机动弹道的优越性，与垂直S型机动弹道对比。二者的突防效果如表6所示。可见最优机动的突防效果更好。

图13为二者的控制曲线。最优组合机动在攻角近似bang-bang机动的同时，不断进行倾侧翻转，纵横向同时机动以躲避拦截，最大限度利用了高超声速飞行器升阻比高、机动性强的优点。

图10 x-z曲线、x-y曲线和速度-时间曲线

图11 弹道倾角-时间曲线

后期攻角再次增大，以降低终端速度，而倾侧角变化则是为了返回落点。因此，通过纵横向组合机动，在满足落点、落角、落速要求的同时，增强了弹道的突防能力。而垂直S型机动是以最大加速度飞行，并通过180°倾侧翻转来改变加速度的符号[17]。程序机动后按照比例导引律打击落点。因此，为了保证终端打击精度，必须缩短机动时间，保留足够能力修正弹道返回落点。二者的三维弹道如图14所示。

表4 优化结果终端状态

表5 突防脱靶量对比

图12 突防弹攻角-时间曲线和倾侧角-时间曲线

机动形式最优机动程序机动脱靶量/m234.6119.7

如图15所示，由于精确打击指标J4中精确打击指标即式(22)～式(24)的应用，控制了最优组合机动弹道的落速，因此其热流密度、动压均很小，这对高超声速飞行器来说至关重要。经过大过载转弯后，过载在终端也控制到较小值，满足了过载约束。

图13 不同机动形式突防弹攻角-时间曲线和倾侧角-时间曲线

图14 三维弹道曲线

图15 热流密度-时间，动压-时间，过载-时间曲线

4.2 最优性验证

为了验证弹道优化策略的必要性和有效性，在相同条件下，将采用GPOPS初值生成器与弹道优化策略后的求解精度对比，如表7所示。可见，采用优化策略后，解的可行性精度和最优性精度提高了多个数量级。

为了验证MRPM优化结果是否是原最优控制问题的最优解，图16给出了最优弹道的哈密顿函数曲线。可见，哈密顿函数值在0附近，满足最优性一阶必要条件，充分说明了MRPM求得的结果是最优的。

表7 求解精度对比

图16 哈密顿函数曲线

4.3 灵敏度分析

对于末段低层拦截，由于攻防双方的相对速度很大，拦截弹的发射时机严重影响着拦截脱靶量。因此，针对前面的优化结果，改变拦截弹的发射时机，统计不同发射时机下的脱靶量。

观察图17的脱靶量统计结果可知，同一条机动弹道，其突防效果依赖于拦截弹的发射时机。但是对于以脱靶量(即J4)为性能指标的机动弹道，由于其脱靶量远远高于其他弹道，即使调整拦截弹的发射时机，脱靶量也在百米量级。因此，直接以脱靶量作为性能指标的最优机动弹道，飞行器的机动能力将得到充分发挥，突防能力最强。

图17 拦截弹发射时机影响

5 结 论

针对高超末段的机动突防与精确打击问题，本文从最优控制角度出发，提出了一种多对象、多段、多约束的末段机动弹道优化方法。通过MRPM求解得到了一种兼顾突防和精确打击的机动弹道，获得了高超声速飞行器的最大机动能力：

1) 通过多对象、多段建模，分别获得了独立、直接的突防性能指标和精确打击性能指标，通过加权处理，获得了综合性能最优的弹道，在保证终端打击精度1 m的前提下，将突防脱靶量提高到百米量级，充分发挥了高超声速飞行器的纵横向机动性能，相比于其他间接性能指标的优化弹道、程序机动弹道，突防效果显著提升。

2) 通过积分弹道作初值、约束渐强以及分段嵌套优化等一系列优化策略，解决了MRPM求解复杂多对象、多段、多约束弹道优化问题的初值敏感、可行域窄等问题，使得优化问题得以快速收敛，且解的可行性精度和最优性精度提高了多个量级。

3) 根据协态映射原理，通过NLP问题的KKT乘子准确估计最优控制问题的协态，求解了哈密顿函数，证明了优化结果的最优性。

4) 通过攻防对抗仿真，验证了最优机动弹道对拦截弹的发射时间不敏感。

参考文献(References)

[1] 陈万春,聂蓉梅,刘佳琪,等.PAC-3爱国者拦截弹末制导精度仿真研究[J].飞航导弹,1999,19(7):57-62.

CHEN W C,NIE R M,LIU J Q,et al.Simulation and research on the terminal guidance precision of PAC-3 patriot missile[J].Winged Missiles Journal,1999,19(7):57-62(in Chinese).

[2] 崔静,姜玉宪.拦截导弹动力学特性对摆动式机动策略突防效果的影响[J].宇航学报,2001,22(5):33-38.

CUI J,JIANG Y X.The effect of interceptor’s dynamic system order on the penetration efficiency of weaving maneuver strategy[J].Journal of Astronautics,2001,22(5):33-38(in Chinese).

[3] 姜玉宪,崔静.导弹摆动式突防策略的有效性[J].北京航空航天大学学报,2002,28(2):133-136.

JIANG Y X,CUI J.Effectiveness of weaving maneuver strategy of a missile[J].Journal of Beijing University of Aeronautics and Astronautics,2002,28(2):133-136(in Chinese).

[4] ZARCHAN P.Proportional navigation and weaving targets[J].Journal of Guidance,Control,and Dynamics,1995,18(5):969-974.

[5] SHINAR J,STEIBERG D.Analysis of optimal evasive maneuvers based on a linearized two-dimensional kinematic model[J]. Journal of Aircraft,1977,14(8):795-802.

[6] SHINAR J,TABAK R.New results in optimal missile avoidance[J]. Journal of Guidance,Control,and Dynamics,1994,17(5):897-902.

[7] YANUSHEVSKY R.Analysis of optimal weaving frequency of maneuvering targets[J].Journal of Spacecraft and Rockets,2004,41(3):477-479.

[8] IMADO F.Some aspects of a realistic three-dimensional pursuit-evasion game[J].Journal of Guidance，Control,and Dynamics,1993,16(2):125-139.

[9] IMADO F,UEHARA S.High-g barrel roll maneuvers against proportional navigation from optimal control,viewpoint[J]. Journal of Guidance,Control,and Dynamics,1998,21(6):876-881.

[10] IMADO F,KURODA T. Engagement tactics for two missiles against an optimally maneuvering aircraft[J].Journal of Guidance,Control,and Dynamics,2011,34(2):574-582.

[11] BENSON D A,THORVALDSEN G T,RAO A V.Direct trajectory optimization and costate estimation via an orthogonal collocation method[J].Journal of Guidance,Control,and Dynamics,2006,29(6):1435-1440.

[12] HUNTINGTON G T,RAO A V.Optimal reconfiguration of spacecraft formation via a Gauss pseudospectral method[J].Journal of Guidance,Control,and Dynamics,2008,31(3):689-698.

[13] ZHANGK N,CHEN W C. Reentry vehicle constrained trajectory optimization[C]∥Proceedings of AIAA International Space Planes and Hypersonic Systems and Technologies Conference. Reston:AIAA,2011,1:1-16.

[14] DARBY C L,RAO A V.Minimum-fuel low earth orbit aeroassisted orbital transfer of small spacecraft[J].Journal of Spacecraft and Rocket,2011,48(4):618-628.

[15] WEISS M,SHIMA T,CASTANEDA D,et al.Minimum effort intercept and evasion guidance algorithms for active aircraft defense[J]. Journal of Guidance,Control,and Dynamics,2016,39(10):2297-2311.

[16] DIVIA G.Advances in global pseudospectral methods for optimal control[D]. Gainesville：University of Floroda,2011:92-96.

[17] ZARCHAN P.Tactical and strategic missile guidance[M]. 6th ed. Reston:AIAA,2012:117-122.