高超末段機動突防/精確打擊彈道建模與優化

李靜琳， 陳萬春， 閔昌萬

(北京航空航天大學 宇航學院, 北京 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所示?？梢娮顑灆C動的突防效果更好。

圖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.