半直接配點法在航天器追逃問題求解中的應(yīng)用

郝志偉，孫松濤，張秋華，諶 穎

(1. 哈爾濱工業(yè)大學(xué)航天科學(xué)與力學(xué)系，哈爾濱 150001； 2. 北京控制工程研究所，北京 100190)

0 引 言

兩航天器追逃問題是一個對抗雙方的控制策略選擇問題，由于航天器動力學(xué)模型高維非線性且雙方的對策目標(biāo)相反，因此這樣的問題求解十分困難。在航天器追逃中，追逐航天器希望選擇控制策略使支付函數(shù)(payoff)最小，而逃逸航天器的期望則恰恰相反——選擇控制策略使支付函數(shù)最大，該追逃過程可由微分對策理論描述。微分對策理論是研究兩個或多個決策者的控制作用同時施加于由微分方程描述的運動系統(tǒng)時實現(xiàn)各自最優(yōu)目標(biāo)的對策過程的理論。該理論模型首先由Isaacs[1]提出，隨著研究的深入，Berkovitz[2]和Friedman[3]給出了微分對策鞍點存在定理和最優(yōu)控制策略存在的必要條件。

微分對策中各方策略最優(yōu)的必要條件雖然已知[2-3]，但針對此必要條件的求解存在困難。因為此必要條件求解對應(yīng)一個非線性兩點邊值問題(TPBVP)求解，而非線性TPBVP一般沒有解析解，只能應(yīng)用數(shù)值算法求解。目前，已存在的主要數(shù)值方法包括配點法[4]和多重打靶法[5]，它們各有優(yōu)缺點。配點法將微分方程在配點處近似成代數(shù)方程，從而求解滿足兩點邊值約束和代數(shù)方程約束的變量，進(jìn)而求解兩點邊值問題。這種方法雖然收斂性好，但計算結(jié)果的精度不高，計算量大，且初值需要選擇。多重打靶法通過數(shù)值積分公式在眾多子區(qū)間上將一個微分方程離散成為代數(shù)方程，在任意子區(qū)間上求解一個初值問題，通過不斷優(yōu)化每一個區(qū)間上的初始值，最終滿足兩點邊值條件。與配點法相比，多重打靶法的計算速度快、計算精度較高，但對初值非常敏感，因此收斂性差。航天器追逃問題一般對應(yīng)復(fù)雜的高維非線性動力學(xué)模型，配點法雖對邊值問題收斂性較強，但若隨機給出初值，因模型復(fù)雜追逃問題數(shù)值解仍然不收斂。文獻(xiàn)[6]提出了一種基于多目標(biāo)遺傳算法的初值選擇方法，解決了航天器追逃問題中非線性TPBVP初值的選擇問題，得到了收斂解。這種初值選擇方法，對配點法和多重打靶法均有效。

本文采用半直接配點法求解航天器追逃問題，提出了一種新的數(shù)值求解追逃問題的思路，避免了非線性TPBVP求解，從而提高了數(shù)值求解追逃問題的效率。半直接法思想來源于Conway對導(dǎo)彈攔截問題的研究[7-8]，此方法通過一種變換將微分對策問題轉(zhuǎn)化成一個最優(yōu)控制問題，再將最優(yōu)控制問題構(gòu)造成一個非線性規(guī)劃問題，然后應(yīng)用序列二次規(guī)劃算法(SQP)進(jìn)行求解。文獻(xiàn)[9]和[10]也利用了這種轉(zhuǎn)化思想將微分對策問題轉(zhuǎn)化為一個非線性規(guī)劃問題求解，文獻(xiàn)[10]還通過在目標(biāo)函數(shù)中加入懲罰函數(shù)解決了狀態(tài)量有約束的導(dǎo)彈追逃問題。但是，Conway的方法并不能證明半直接配點法求解微分對策問題得到的解一定是原問題的解。文獻(xiàn)[11]證明了半直接配點法求解微分對策問題的等價性，為半直接配點法求解微分對策問題提供了理論依據(jù)。在求解兩航天器追逃問題的過程中，通過半直接變換將微分對策問題轉(zhuǎn)化為最優(yōu)控制問題，通過Guass-Lobbato配點法數(shù)值求解此問題，提高了數(shù)值方法的收斂性和穩(wěn)定性。

本文針對近地軌道附近時間固定的兩航天器追逃問題，假設(shè)對抗雙方均為連續(xù)小推力、對策時間較短、且瞬時狀態(tài)信息完全已知，以距離為支付建立對策模型，給出了半直接配點法求解此追逃問題的數(shù)值方法。研究最終將給出追逃雙方在對策條件下的最優(yōu)策略，仿真算例表明半直接配點法為航天器追逃問題提供一種有效的求解方法。

1 航天器追逃問題建模

考慮近地軌道附近的兩航天器P與E，O為P與E附近的一顆虛擬參考星。在以下內(nèi)容中，P與E分別表示追逐航天器與逃逸航天器，兩航天器的位置關(guān)系如圖1所示，其中O-XYZ為參考坐標(biāo)系，OX方向為地心連接參考星質(zhì)心的連線方向，OY在參考星軌道平面內(nèi)，延軌道速度方向，OZ軸與OX、OY軸構(gòu)成右手坐標(biāo)系，r(t)為參考星O點瞬時半徑，xi(t)、yi(t)與zi(t)(i=P，E)分別表示兩追逃航天器在參考坐標(biāo)系X，Y和Z方向?qū)?yīng)的分量。建立兩航天器針對參考星的動力學(xué)方程如下：

圖1 兩航天器對策的坐標(biāo)示意圖Fig.1 Coordinate frame sketch of pursuer and evader

(1)

此航天器追逃問題的支付函數(shù)考慮為追逃對策結(jié)束時兩航天器間距離：

2 最優(yōu)策略存在的必要條件

(2)

為了方便，令：

兩航天器初值給定：

(3)

根據(jù)微分對策原理[1]，引入?yún)f(xié)態(tài)量λi(i=P,E)，

λi=[λ1i,λ2i,λ3i,λ4i,λ5i,λ6i]T,i=P,E

Hamilton方程H定義為：

(4)

根據(jù)最優(yōu)策略存在的必要條件[2]，協(xié)態(tài)方程為：

(5)

相應(yīng)的邊值條件為：

(6)

[λ4P(t),λ5P(t),λ6P(t)]T

(7)

[λ4E(t),λ5E(t),λ6E(t)]T

(8)

顯然，式(2)、最優(yōu)策略必要條件(5)、(6)以及兩航天器初值條件(3)組成了一個兩點邊值問題。一般的兩點邊值問題可由多重打靶法求解[8]。第3節(jié)將給出半直接法配點求解此上述兩點邊值問題的思路和過程。

3 半直接配點法

在半直接配點法中，將第2節(jié)最優(yōu)策略必要條件得到的追逐或逃逸航天器的協(xié)態(tài)變量方程加入到狀態(tài)方程中，同時將對應(yīng)航天器的最優(yōu)控制量表達(dá)式代入對應(yīng)的狀態(tài)方程并附上約束條件，將一個微分對策問題轉(zhuǎn)化為一個最優(yōu)控制問題。在求解中，只需求解一個對應(yīng)的支付函數(shù)最大或最小問題即可，相應(yīng)的等價性證明請參考文獻(xiàn)[11]。具體求解步驟如下。

3.1 半直接轉(zhuǎn)換

半直接配點法將微分對策問題轉(zhuǎn)化為最優(yōu)控制問題，這里將逃逸航天器E對應(yīng)的協(xié)態(tài)變量λE加入到狀態(tài)微分方程中。變換的方法簡述如下，令：

(9)

由于已經(jīng)將協(xié)態(tài)變量λE加入到系統(tǒng)方程(9)中，那么相應(yīng)的終端約束方程Ψ定義為：

(10)

式中：

則變換后的支付函數(shù)僅為uP的函數(shù)：

(11)

顯然，微分對策問題已被轉(zhuǎn)化為一個最優(yōu)控制問題(9)～(11)。對應(yīng)追逐航天器P的半直接轉(zhuǎn)化過程同上。接下來，將給出基于配點法求解轉(zhuǎn)化后最優(yōu)控制問題的方法。

3.2 半直接配點法的數(shù)值求解方法

考慮到計算精度，本文采用五階Gauss-Lobbato配點法[12]求解上述最優(yōu)問題。為說明此方法，不妨設(shè)式(9)對應(yīng)的微分方程為

(12)

根據(jù)配點法，將時間區(qū)間[t0,tf]分為N個子區(qū)間(各個子區(qū)間的長度可以不相等)，每個子區(qū)間為[ti,ti+1](i=1,2,…N,t1=t0,tN+1=tf)，這樣就確定了N+1個節(jié)點的位置，并得到初始的離散點。在每個子區(qū)間[ti,ti+1]上，應(yīng)用多項式近似式(12)需先確定配點。由Gauss-Lobbato配點法[13]可知，除端點ti,ti+1外，其余三個配點tc1,tcm和tc2分別為：

(13)

(14)

(15)

將式(14)與式(15)相加可得：

(16)

i=1,2,…,N

(17)

i=1,2,…,N

(18)

式中：符號N表示時間區(qū)間[t0,tf]上的子區(qū)間個數(shù)。

4 數(shù)值仿真結(jié)果

為說明半直接配點法求解航天器追逃問題的特點，將給出2個仿真算例。因為空間站和大量的衛(wèi)星存在于地球低軌道，在仿真中參考星O軌道高度分別選為500 km和1000 km，追逃對策時間均為500 s。

表1給出兩個實例中追逐與逃逸航天器的初值，表2給出了實例中兩航天的單位質(zhì)量加速度，其中g(shù)=9.8×10-3km/s。

為了兼顧計算速度和精度，將時間區(qū)間劃分為10個子區(qū)間進(jìn)行仿真。由子區(qū)間的選擇可知，需要求解的變量為483(2×10×18+(4×10+1)×3)個。

在數(shù)值仿真計算中，應(yīng)用Gauss-Lobatto配點法并采用MATLAB語言編寫的SNOPT求解器進(jìn)行求解[13]。在仿真中，選用Core i7-4790K處理器、內(nèi)存8 G的臺式機進(jìn)行數(shù)值計算。設(shè)約束精度和優(yōu)化允許誤差均為10-9。

在算例1中(參考星軌道高度500 km)，兩航天器的軌跡和最優(yōu)控制量變化曲線如圖2所示。圖2(a)～圖2(f)分別給出了追逐航天器P與逃逸航天器E在固定時間500 s內(nèi)X,Y和Z方向的位置與速度變化軌跡。追逐航天器P的最優(yōu)控制量隨時間變化曲線如圖2(g)所示，逃逸航天器E的最優(yōu)控制量隨時間變化曲線如圖2(h)所示。由圖2可知，通過半直接配點法可獲得航天器追逃問題收斂的解，圖2(a)～圖2(c)顯示追逐航天器P零距離追上逃逸航天器E。在仿真算例1中，數(shù)值計算的約束誤差為1.7532×10-11，計算時間為14.13 s。

在算例2中(參考星軌道高度1000 km)，兩航天器的軌跡和最優(yōu)控制量變化曲線如圖3所示。與算例1相同，圖3(a)～圖3(f)分別給出了追逐航天器P與逃逸航天器E在固定時間500 s內(nèi)X,Y和Z方向的位置與速度變化軌跡。追逐航天器P與逃逸航天器E的最優(yōu)控制量隨時間變化曲線如圖3(g)和圖3(h)所示。由圖3可知，通過半直接配點法可獲得航天器追逃問題收斂的解，圖3(a)～圖3(c)顯示追逐航天器依然可以零距離追上逃逸航天器。在仿真算例2中，數(shù)值計算的約束誤差為2.3584×10-11，計算時間為16.39 s。

由表2可知，在上述兩個算例中，追逐航天器P的最大單位質(zhì)量加速度均比逃逸航天器E的最大單位質(zhì)量加速度大，這是追逐航天器P可以零距離追上逃逸航天器E的主要原因。

5 結(jié) 論

本文采用半直接配點法求解了近地軌道附近的兩航天器追逃問題，得到了收斂的數(shù)值解。半直接配點法將一個微分對策問題轉(zhuǎn)化成了一個最優(yōu)控制問題，應(yīng)用Gauss-Lobbato五階配點法并結(jié)合序列二次規(guī)劃法，最終解決了一個非線性的數(shù)學(xué)規(guī)劃問題。采用半直接配點法求解微分對策問題時，可以避免求解困難的兩點邊值問題。這種算法具有收斂性好，應(yīng)用簡單的特點。數(shù)值仿真實例說明了這種求解方法在解決航天器追逃問題上的有效性。

圖2 算例1：追逐和逃逸航天器的軌跡和最優(yōu)控制變量Fig.2 Test case 1: trajectories and optimal controls of pursuer and evader spacecraft

圖3 算例2：追逐和逃逸航天器的軌跡和最優(yōu)控制變量Fig.3 Test case 2: trajectories and optimal controls of pursuer and evader spacecraft

算例xP/kmyP/kmzp/kmvxP/(km·s-1)vyP/(km·s-1)vzp/(km·s-1)算例1000000算例2000000算例xE/kmyE/kmzE/kmvxE/(km·s-1)vyE/(km·s-1)vzE/(km·s-1)算例123.1430.08-0.011-0.033-0.0130.013算例2-10.1444.970-0.0450.00480.013

表2 兩航天器的單位質(zhì)量最大加速度Table 2 Max acceleration-to-mass of the pursuer and evader spacecraft