符 俊，蔡 洪，張士峰，丁洪波

(國(guó)防科技大學(xué)航天與材料工程學(xué)院，長(zhǎng)沙 410073)

0 引言

迄今有很多學(xué)者對(duì)空間攔截問(wèn)題進(jìn)行了研究，一般來(lái)說(shuō)，它們可分為2類:脈沖推力，如Herrick方法、Godal方法、Gauss方法、普適變量法和 Lambert方法等［1］;有限推力，如文獻(xiàn)［2－4］。對(duì)于有限推力攔截情況，往往是在一定的性能指標(biāo)下尋找最優(yōu)的攔截軌道，這是一個(gè)最優(yōu)控制問(wèn)題。

隨著計(jì)算機(jī)水平的高速發(fā)展，求解最優(yōu)控制問(wèn)題的數(shù)值方法得到了廣泛應(yīng)用，其中魯棒性和通用性更好的直接法往往能更好地求解動(dòng)態(tài)優(yōu)化問(wèn)題。它將函數(shù)空間中的最優(yōu)控制問(wèn)題轉(zhuǎn)化為歐式空間中的非線性規(guī)劃問(wèn)題，然后再利用各種非線性算法完成問(wèn)題的求解［5］。

本文基于非線性化C-W方程，利用Legendre偽譜法(Legendre Pseudospectral Method，LPM)對(duì)航天器遠(yuǎn)程攔截軌道進(jìn)行優(yōu)化，同時(shí)基于固定時(shí)間飛行定理討論了一種脈沖推力假設(shè)下的優(yōu)化方案，并將這2種方法進(jìn)行了對(duì)比分析。

1 數(shù)學(xué)模型

1.1 動(dòng)力學(xué)方程

C-W方程常用于描述圓參考軌道下的航天器近距離相對(duì)運(yùn)動(dòng)，這主要是因?yàn)镃-W方程在建立的過(guò)程中作了以下2個(gè)假設(shè):(1)目標(biāo)航天器的軌道為圓軌道;(2)攔截器與目標(biāo)航天器之間的距離遠(yuǎn)小于它們的地心距。但在攔截問(wèn)題中，假設(shè)(2)極大地限制了C-W方程的應(yīng)用，因?yàn)閿r截器和目標(biāo)航天器之間的距離往往很大，故假設(shè)2不再成立。在這種情況下，下面給出適于描述遠(yuǎn)程攔截問(wèn)題的航天器相對(duì)運(yùn)動(dòng)方程(目標(biāo)航天器不機(jī)動(dòng)):

式中 ［x y z vxvyvz］T表示攔截器在目標(biāo)航天器LVLH(Local Vertical，Local Horizontal)坐標(biāo)系中的相對(duì)運(yùn)動(dòng)狀態(tài)矢量;Tx、Ty、Tz為攔截器的推力在LVLH坐標(biāo)系三軸上的分量;m為攔截器的質(zhì)量;Isp為攔截器推進(jìn)劑比沖;a為目標(biāo)航天器的地心距;g為重力加速度;μ為地球引力常數(shù);rc為攔截器的地心距，為目標(biāo)航天器的平均軌道角速度為攔截器推力

在優(yōu)化問(wèn)題中，控制量為 U=［Tx，Ty，Tz］T，性能指標(biāo)應(yīng)為攔截器消耗的燃料最少，即終端時(shí)刻質(zhì)量最大:

式中 tf表示終端時(shí)刻。

1.2 無(wú)量綱化參數(shù)

采用數(shù)值方法進(jìn)行優(yōu)化時(shí)，必須對(duì)物理量進(jìn)行無(wú)量綱化，選取的無(wú)量綱化參數(shù)如下:

式中 Rs為目標(biāo)航天器的地心距;m0為攔截器的初始質(zhì)量。

2 有限推力最優(yōu)攔截

2.1 Legendre偽譜法簡(jiǎn)介

Legendre偽譜法由美國(guó)海軍研究院Fahroo Fariba和Ross I Micheal提出［6］，它屬于直接法中的一種。該方法將狀態(tài)變量和控制變量在一系列LGL(Legendre-Gauss-Lobatto)點(diǎn)上離散，并以離散點(diǎn)為節(jié)點(diǎn)構(gòu)造Lagrange插值多項(xiàng)式來(lái)逼近狀態(tài)變量和控制變量，從而將動(dòng)態(tài)最優(yōu)控制問(wèn)題轉(zhuǎn)化為靜態(tài)參數(shù)優(yōu)化問(wèn)題(NLP問(wèn)題)。Legendre偽譜法通過(guò)對(duì)全局插值多項(xiàng)式求導(dǎo)來(lái)近似狀態(tài)變量對(duì)時(shí)間的導(dǎo)數(shù)，從而將微分方程約束轉(zhuǎn)換為一組代數(shù)約束，這些約束條件加上問(wèn)題本身的約束條件，如邊界約束、路徑約束等，共同構(gòu)成NLP問(wèn)題的約束條件。Legendre偽譜法的求解步驟可參考文獻(xiàn)［3 －4，7］。

2.2 優(yōu)化結(jié)果及分析

以一具體的軌道攔截任務(wù)為例，研究Legendre偽譜法的應(yīng)用。仿真條件:攔截器總質(zhì)量5 00 kg，推進(jìn)劑質(zhì)量4 00 kg，比沖300 s，最大推力2 000 N。初始時(shí)刻 t0，攔截器的軌道根數(shù):a=16 678.137 km，e=0.2，i=28.5°，ω =0，Ω =0，?=1.54°。目標(biāo)航天器的軌道根數(shù):a=6 878.137 km，e=0，i=30°，Ω =0，u=2.05°。

經(jīng)過(guò)計(jì)算，可得到t0時(shí)刻攔截器在目標(biāo)航天器LVLH坐標(biāo)系中的相對(duì)運(yùn)動(dòng)狀態(tài):

其中，前3項(xiàng)表示相對(duì)位置，后3項(xiàng)表示相對(duì)速度，距離單位為km，速度單位為km/s。

仿真計(jì)算中采用60個(gè)LGL點(diǎn)，優(yōu)化結(jié)果見(jiàn)圖1～圖4。

圖1 控制變量的優(yōu)化結(jié)果Fig.1 Optimal results of control variables

圖1(a)為推力隨時(shí)間的變化曲線，可看出在攔截過(guò)程中，攔截器受到的推力主要沿著跡向，在另外2個(gè)方向上的分量很小。發(fā)動(dòng)機(jī)工作時(shí)產(chǎn)生的推力一般沿體系x方向，為了使推力在LVLH坐標(biāo)系三軸上的分量大小滿足要求，可通過(guò)控制攔截器的姿態(tài)來(lái)實(shí)現(xiàn)。在本次仿真中，推力作用的時(shí)間為506.9 s，最大值為2 000 N，符合最大推力約束條件。之后的過(guò)程發(fā)動(dòng)機(jī)停止工作，攔截器在地心引力的作用下飛向目標(biāo)航天器進(jìn)行攔截。圖1(b)為攔截器的過(guò)載變化曲線圖，從圖1(b)可看出，最大過(guò)載出現(xiàn)在y方向，但均不超過(guò)1，滿足容許過(guò)載較小的航天器的變軌要求。

顏曉晨正在試衣服，一個(gè)二十五六歲的長(zhǎng)發(fā)女子走了進(jìn)來(lái)，看了她幾眼，拿了一套顏曉晨試穿的衣服，翻看價(jià)格牌。一個(gè)營(yíng)業(yè)員在接電話，另一個(gè)營(yíng)業(yè)員正低著頭幫顏曉晨整理褲腳，都沒(méi)顧上招呼她，顏曉晨笑著說(shuō)：“全場(chǎng)五折?！?/p>

圖2 狀態(tài)變量?jī)?yōu)化結(jié)果Fig.2 Optimal results of state variables

圖2為攔截器的狀態(tài)變量?jī)?yōu)化結(jié)果圖。其中圖2(a)為相對(duì)距離變化曲線，可看到，經(jīng)歷8 080.5 s后，攔截器與目標(biāo)航天器的相對(duì)距離趨于0，成功實(shí)現(xiàn)攔截。分析圖2(a)可發(fā)現(xiàn)，攔截器與目標(biāo)航天器的相對(duì)距離經(jīng)歷了一個(gè)先增大后減小的過(guò)程，而不是單調(diào)減小，這是因?yàn)樾阅苤笜?biāo)為能量最省，在前面給定的初始相對(duì)運(yùn)動(dòng)狀態(tài)條件下，如果直接控制攔截器朝目標(biāo)航天器飛去，勢(shì)必會(huì)消耗很大的能量，故攔截器與目標(biāo)航天器之間的距離不是單調(diào)減小。圖2(b)為相對(duì)速度變化曲線，圖2(c)為攔截器的質(zhì)量變化圖。攔截器最終質(zhì)量為 233.5 kg，燃料消耗為 266.5 kg。

圖3 相對(duì)運(yùn)動(dòng)軌跡Fig.3 Trajectories of relative motion

圖4 哈密頓函數(shù)隨時(shí)間的變化曲線Fig.4 Time histories of Hamiltonian

圖3為攔截器相對(duì)目標(biāo)航天器的運(yùn)動(dòng)軌跡，從圖3看到，它們之間的相對(duì)運(yùn)動(dòng)主要發(fā)生在x-y平面，是由于初始時(shí)刻兩者的軌道面空間位置差異不大造成的。

在應(yīng)用Legendre偽譜法處理最優(yōu)控制問(wèn)題時(shí)，往往面臨著一個(gè)問(wèn)題:解的最優(yōu)性能否得到保證?這可通過(guò)求解哈密頓函數(shù)的值是否滿足一階最優(yōu)性條件來(lái)進(jìn)行判斷:作為區(qū)別于一般直接數(shù)值解法的重要特性之一，Legendre偽譜法將最優(yōu)控制問(wèn)題轉(zhuǎn)化為L(zhǎng)GL點(diǎn)上的非線性規(guī)劃問(wèn)題后，可將對(duì)該問(wèn)題的求解轉(zhuǎn)化為對(duì)一個(gè)增廣性能指標(biāo)的優(yōu)化求解，并能在求解該問(wèn)題的同時(shí)得到Lagrange乘子，即LGL點(diǎn)上的協(xié)態(tài)變量值，從而計(jì)算哈密頓函數(shù)值以檢驗(yàn)是否滿足一階最優(yōu)性必要條件。如果滿足則可認(rèn)為結(jié)果最優(yōu)，否則就不是最優(yōu)的［8－9］。

本次仿真中的性能指標(biāo)為J=－m(tf)，這是一個(gè)末值型的性能指標(biāo)，且終端時(shí)刻自由。根據(jù)龐特李亞金極大值原理，在優(yōu)化過(guò)程中，哈密頓函數(shù)值應(yīng)恒為0。圖4為哈密頓函數(shù)隨時(shí)間的變化曲線，從圖4可看到，該值一直處于0附近，最大相差為0.7%，滿足最優(yōu)性必要條件，因此解的最優(yōu)性能夠得到保證。

3 基于脈沖推力的最優(yōu)攔截方案

前面利用Legendre偽譜法得到了最優(yōu)攔截問(wèn)題的解。下面討論在脈沖推力假設(shè)下航天器遠(yuǎn)程最優(yōu)攔截方案的研究方法，理論基礎(chǔ)是固定時(shí)間攔截定理。

3.1 Battin-Vaughn 方法

求解固定時(shí)間攔截問(wèn)題的經(jīng)典算法是高斯方法和Lambert飛行時(shí)間定理。通過(guò)對(duì)上述2種算法的改進(jìn)，派生出了很多算法，如海里克方法、歌德方法、普適變量法、p迭代法和Battin-Vaughan方法。其中，普適變量法算法簡(jiǎn)單，具有普適性，適用于所有的圓錐曲線，但在某些情況下，收斂速度太慢;p迭代法用牛頓迭代法來(lái)修正p的試探值，收斂速度很快，其不足在于r1、r2共線情況下不收斂;而B(niǎo)attin-Vaughan方法通過(guò)引入超幾何函數(shù)和連分?jǐn)?shù)進(jìn)行數(shù)值迭代求解，迭代速度與精度均比較理想，適合計(jì)算機(jī)快速求解對(duì)于任意情況都收斂很快，具有幾乎一致收斂的性質(zhì)。據(jù)Klumpp在1986年以來(lái)對(duì)Battin-Vaughan算法就所有合理情況所做的廣泛實(shí)驗(yàn)研究表明，Battin的算法對(duì)于載人和不載人自主制導(dǎo)都能提供所需要的可靠性、快速性和列緊性，并且對(duì)于行星軌道確定以及整個(gè)天體力學(xué)也都提供了所需要的精度和應(yīng)用范圍［10］。

3.2 優(yōu)化方案

在空間攔截任務(wù)中，關(guān)鍵技術(shù)是在一個(gè)時(shí)間窗口(也稱任務(wù)窗口)［t0，t1］中攔截器能命中目標(biāo)并使消耗的能量最小。根據(jù)固定時(shí)間攔截定理，單圈Lambert問(wèn)題存在唯一解。因此，攔截器不同的機(jī)動(dòng)時(shí)刻tman對(duì)應(yīng)著攔截目標(biāo)航天器所需的不同速度沖量Δv，這稱為一種攔截方案，以(tman，Δv)表征。根據(jù)這種思路，在實(shí)際計(jì)算過(guò)程中，首先計(jì)算出一次攔截任務(wù)中任務(wù)窗口內(nèi)的所有機(jī)動(dòng)攔截方案，然后在這些方案中選擇一種最優(yōu)的方案(比如所需特征速度最小或攔截時(shí)間最短)或者是滿足特定攔截任務(wù)要求(在指定時(shí)刻攔截目標(biāo))的方案，這種靈活性是該方法的一大優(yōu)點(diǎn)。

尋求最優(yōu)攔截方案的關(guān)鍵是求解機(jī)動(dòng)時(shí)刻tman，可采用遺傳算法進(jìn)行求解，也可采用枚舉法，這里選擇枚舉法。采用枚舉法尋找最優(yōu)攔截方案，具體要求是在時(shí)間窗口內(nèi)，選擇機(jī)動(dòng)時(shí)刻tman作為迭代變量，求解特征速度最小的機(jī)動(dòng)方案，迭代步長(zhǎng)可視具體情況靈活控制。為克服枚舉法計(jì)算量大的缺點(diǎn)，可通過(guò)變步長(zhǎng)控制計(jì)算量。求解流程見(jiàn)圖5。

圖5 仿真流程Fig.5 Flow of simulation

3.3 仿真分析

仿真條件與2.2節(jié)相同，為與偽譜法的優(yōu)化結(jié)果進(jìn)行對(duì)比，選擇攔截器的任務(wù)時(shí)間窗口為［t0，t1］，其中t1=8 080.5 s。

圖6為機(jī)動(dòng)時(shí)刻與所需速度沖量的關(guān)系，橫坐標(biāo)表示機(jī)動(dòng)時(shí)刻tman，縱坐標(biāo)表示所需的速度增量Δv。t0為零時(shí)刻，可看出，當(dāng)機(jī)動(dòng)時(shí)刻tman=4 328 s時(shí)，所需的速度增量最小，為1 501.73m/s。

圖6 機(jī)動(dòng)時(shí)刻與速度沖量的關(guān)系Fig.6 Relationship of maneuver time and velocity impulse

攔截器最優(yōu)變軌過(guò)程如表1所示，變軌所需的速度沖量為［1 091.58，－930.36，－445.03］(J2000 坐標(biāo)系)，單位為 m/s;變軌點(diǎn)的軌道位置為 ?=241.91°，攔截點(diǎn)的軌道位置為?=86.25°。攔截器在原軌道上運(yùn)行4 328 s后開(kāi)始變軌?？梢钥闯?，變軌前后攔截器的軌道面并沒(méi)有改變，只是軌道的形狀和大小發(fā)生了變化，因此所需的速度增量不是很大。根據(jù)公式Δv=Ispg ln(m+/m－)計(jì)算得到消耗的燃料為199.99 kg。

表1 攔截器變軌前后軌道根數(shù)Table 1 Orbital elements of interceptor

4 2種方法的比較

為便于描述，Legendre偽譜法用方法Ⅰ表示，基于脈沖推力的優(yōu)化方法用方法Ⅱ表示。根據(jù)前面的分析，發(fā)現(xiàn)方法Ⅱ的燃料消耗要小于方法Ⅰ，其中一個(gè)原因?yàn)榉椒á蛲ㄟ^(guò)選擇攔截器與目標(biāo)航天器相對(duì)相位關(guān)系最優(yōu)時(shí)進(jìn)行變軌，如圖7所示。方法I是直接在初始相對(duì)相位關(guān)系下實(shí)施攔截機(jī)動(dòng)變軌，因此使得攔截消耗的能量較大。2種方法的綜合對(duì)比詳情見(jiàn)表2。

圖7 攔截器變軌圖Fig.7 Orbit change of interceptor

表2 2種方法的比較Table 2 Comparison of two kinds of methods

5 結(jié)論

(1)仿真過(guò)程中發(fā)現(xiàn)，Legendre偽譜法精度高，對(duì)變量初值不敏感，收斂半徑大，適于處理航天器遠(yuǎn)程攔截問(wèn)題。

(2)在有限推力攔截中，Legendre偽譜法用來(lái)優(yōu)化推力的大小和方向，從而得到最優(yōu)攔截軌道;而在脈沖推力最優(yōu)攔截方法中，主要是通過(guò)調(diào)整攔截器和目標(biāo)航天器的相位關(guān)系來(lái)搜索最優(yōu)攔截軌道。

(3)與Legendre偽譜法相比，采用基于脈沖的最優(yōu)機(jī)動(dòng)攔截方案具有更大的靈活性，既能夠滿足性能指標(biāo)要求，也能滿足具體攔截任務(wù)要求。

［1］任萱.人造地球衛(wèi)星軌道力學(xué)［M］.長(zhǎng)沙:國(guó)防科技大學(xué)出版社.1988:141-157.

［2］涂良輝，袁建平，羅建軍.基于偽光譜方法的有限推力軌道轉(zhuǎn)移優(yōu)化設(shè)計(jì)［J］.宇航學(xué)報(bào)，2009，29(4):1189-1193.

［3］胡正東，丁洪波，曹淵，等.偽譜法在SGKW軌道快速優(yōu)化中的應(yīng)用［J］.航天控制，2009，27(4):3-7.

［4］宣穎，張為華，張育林.基于Legendre偽譜法的固體運(yùn)載火箭軌跡優(yōu)化研究［J］.固體火箭技術(shù)，2008，31(5):425-429.

［5］羅亞中，唐國(guó)金，田蕾.基于模擬退火算法的最優(yōu)控制問(wèn)題全局優(yōu)化［J］.南京理工大學(xué)學(xué)報(bào)，2005，29(2):144-148.

［6］Fahroo F，Ross I M.Costate estimation by a legendre pseudospectral method［J］.Journal of Guidance，Control and Dynamics，2002，24(2):270-277.

［7］Huntington Geoffrey Todd.Advancement and analysis of a gauss pseudospectral transcription for optimal control problems［D］.Cambridge，Massachusetts Institute of Technology，2007:115-143.

［8］Qi G，Ross I M，Wei K，et al.Connections between the convector mapping theorem and convergence of pseudospectral methods for optimal control［J］.Compute Optimal Application，DOI.10.1007/s10589-007-910-2-4，2007(1):1-29.

［9］丁洪波，蔡洪，張士峰，等.高超聲速滑翔式再入飛行器最大航程飛行軌跡分析［J］.國(guó)防科技大學(xué)學(xué)報(bào)，2009，31(6):67-72.

［10］趙瑞安.空間武器軌道設(shè)計(jì)［M］.北京:中國(guó)宇航出版社，2008:176-178.

［11］Richard H Battin.An introduction to the mathematics and methods of astrodynamics［M］.AIAA Education Series，AIAA，Reston，VA，1999.