航天器遠(yuǎn)程最優(yōu)攔截方法研究①

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

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

0 引言

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

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

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

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

1.1 動力學(xué)方程

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

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

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

式中 tf表示終端時刻。

1.2 無量綱化參數(shù)

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

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

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

2.1 Legendre偽譜法簡介

Legendre偽譜法由美國海軍研究院Fahroo Fariba和Ross I Micheal提出［6］，它屬于直接法中的一種。該方法將狀態(tài)變量和控制變量在一系列LGL(Legendre-Gauss-Lobatto)點上離散，并以離散點為節(jié)點構(gòu)造Lagrange插值多項式來逼近狀態(tài)變量和控制變量，從而將動態(tài)最優(yōu)控制問題轉(zhuǎn)化為靜態(tài)參數(shù)優(yōu)化問題(NLP問題)。Legendre偽譜法通過對全局插值多項式求導(dǎo)來近似狀態(tài)變量對時間的導(dǎo)數(shù)，從而將微分方程約束轉(zhuǎn)換為一組代數(shù)約束，這些約束條件加上問題本身的約束條件，如邊界約束、路徑約束等，共同構(gòu)成NLP問題的約束條件。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。初始時刻 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)過計算，可得到t0時刻攔截器在目標(biāo)航天器LVLH坐標(biāo)系中的相對運動狀態(tài):

其中，前3項表示相對位置，后3項表示相對速度，距離單位為km，速度單位為km/s。

仿真計算中采用60個LGL點，優(yōu)化結(jié)果見圖1～圖4。

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

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

顏曉晨正在試衣服，一個二十五六歲的長發(fā)女子走了進(jìn)來，看了她幾眼，拿了一套顏曉晨試穿的衣服，翻看價格牌。一個營業(yè)員在接電話，另一個營業(yè)員正低著頭幫顏曉晨整理褲腳，都沒顧上招呼她，顏曉晨笑著說：“全場五折。”

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

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

圖3 相對運動軌跡Fig.3 Trajectories of relative motion

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

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

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

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

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

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

3.1 Battin-Vaughn 方法

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

3.2 優(yōu)化方案

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

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

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

3.3 仿真分析

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

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

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

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

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

4 2種方法的比較

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

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

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

5 結(jié)論

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

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

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

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

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

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

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

［5］羅亞中，唐國金，田蕾.基于模擬退火算法的最優(yōu)控制問題全局優(yōu)化［J］.南京理工大學(xué)學(xué)報，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］.國防科技大學(xué)學(xué)報，2009，31(6):67-72.

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

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