本立言，謝祥華，張 銳

(上海微小衛(wèi)星工程中心，上海 201203)

0 引 言

London[1]首次提出氣動輔助變軌方式，這對于航天器軌道轉(zhuǎn)移任務(wù)而言，可以節(jié)約大量的燃料，吸引了眾多學(xué)者的研究。氣動輔助變軌是指具有一定形狀的航天器經(jīng)過大氣層，通過調(diào)整升力系數(shù)以及阻力系數(shù)，來降低航天器能量或者改變軌道平面。在氣動輔助變軌問題中，主要熱點集中于航天器在大氣層內(nèi)的飛行控制，Miele等[2]研究了航天器在大氣層中的飛行規(guī)律，給出了9種軌道轉(zhuǎn)移最優(yōu)控制指標(biāo)，其中比較重要的有燃料最省指標(biāo)和時間最短指標(biāo)。

在氣動輔助變軌問題中，傳統(tǒng)數(shù)值方法可分為直接法和間接法。直接法通過不同的離散化方法對控制量進(jìn)行參數(shù)化，把原優(yōu)化問題轉(zhuǎn)化為帶有約束的參數(shù)優(yōu)化問題，然后采用非線性規(guī)劃算法進(jìn)行求解。這類方法目前的研究熱點是偽譜法[3-6]，該方法雖然具有求解精度較高、尋優(yōu)參數(shù)少、收斂性好的優(yōu)點，但仍然需要猜測初值，嚴(yán)重依賴非線性規(guī)劃算法的求解能力。對于氣動輔助變軌問題，控制律較為震蕩[7]，數(shù)值經(jīng)驗表明，偽譜法在求解此類問題時效率以及計算精度都出現(xiàn)明顯的降低[8]，由于偽譜法很難捕獲到控制跳躍點，會出現(xiàn)虛擬數(shù)值振蕩[9-10]。為得到精確的控制律，常用的改進(jìn)算法有hp自適應(yīng)偽譜法[11-13]，通過劃分網(wǎng)格和增加配點個數(shù)，來控制參數(shù)優(yōu)化問題規(guī)模和提高逼近精度，但失去了全局快速收斂優(yōu)勢。

間接法首先根據(jù)極大值原理將最優(yōu)控制問題轉(zhuǎn)化為兩點邊值問題[14]，然后利用各種數(shù)值方法求解[15-17]，常用的有打靶法。盡管采用間接法可以得到精確的最優(yōu)控制律，但這種方法求解過程比較復(fù)雜，需要獲取目標(biāo)函數(shù)或者終值誤差函數(shù)對于狀態(tài)變量以及協(xié)態(tài)變量的梯度矩陣，而且對初值非常敏感，由于協(xié)態(tài)變量沒有明確的物理含義，因此初值難以猜測，收斂范圍小，給求解問題帶來了很大的困難。

針對上述研究現(xiàn)狀，本文提出一種基于無損卡爾曼濾波(UKF) 參數(shù)估計的氣動輔助變軌問題求解方法，首先將氣動輔助變軌所對應(yīng)的兩點邊值問題轉(zhuǎn)化為參數(shù)估計問題，然后用UKF濾波器進(jìn)行參數(shù)估計，從而得到精確解。由于UKF濾波算法基于估計理論，不依賴梯度信息，避免梯度矩陣的推導(dǎo)，而且克服了猜測協(xié)態(tài)變量初值的困難，降低了求解氣動輔助變軌問題的難度。數(shù)值仿真結(jié)果表明，該方法求解效率高，具有良好的魯棒性，可以很好的解決氣動輔助變軌問題。

1 過程描述

氣動輔助軌道變軌過程如圖1所示。在初始軌道施加一個反向的速度增量ΔV1后，航天器進(jìn)入轉(zhuǎn)移軌道，該軌道的遠(yuǎn)地點為變軌點1，航天器于特征點2進(jìn)入大氣層后，在大氣層內(nèi)通過調(diào)整升力系數(shù)，在特征點3飛出大氣層，進(jìn)入轉(zhuǎn)移軌道，該軌道的遠(yuǎn)地點為特征點4，在特征點4施加速度增量ΔV2，使航天器進(jìn)入目標(biāo)軌道，從而完成氣動輔助軌道轉(zhuǎn)移。記μ為地球引力常數(shù)，r1和r2分別為初始軌道和目標(biāo)軌道的地心距，ra為大氣層邊界地心距。

為了獲取足夠氣動力，航天器必須進(jìn)入大氣層，即轉(zhuǎn)移軌道近地點必須確保在大氣層內(nèi)，令近地點地心距為rpe，則rpe

(1)

由于rpe值對優(yōu)化結(jié)果并不敏感，在本文算例中，若rpe相差50 km，則ΔV1僅改變9 m/s，因此可認(rèn)為固定不變。記大氣入口點的地心距r0，速度V0以及航跡角γ0分別滿足

r0=ra

(2)

(3)

(4)

根據(jù)切入點約束條件，保證轉(zhuǎn)移軌道遠(yuǎn)地點高度為r2，根據(jù)動量守恒以及能量守恒定律，則大氣出口點的地心距rf、速度Vf和航跡角γf滿足

rf=ra

(5)

(6)

特征點4變軌前后速度分別為

(7)

(8)

則第二次施加速度增量

(9)

速度增量ΔV2的大小與大氣出口點的航跡角極為敏感，當(dāng)γf=0°，ΔV2最優(yōu)，為了確保航天器能順利飛出大氣層，則需要稍大于0°，因此認(rèn)為γf為一稍大于0°的固定值，所以，在航天器氣動輔助軌道轉(zhuǎn)移的過程中，主要優(yōu)化過程集中于大氣飛行段，在滿足入口r0，V0，r0出口rf，Vf，γf條件下，通過調(diào)整升力系數(shù)，使航天器盡快從大氣層中逃逸，避免氣動加熱時間過長，因此將轉(zhuǎn)移時間最優(yōu)作為性能指標(biāo)。

2 優(yōu)化模型

假設(shè)航天器在大氣內(nèi)飛行僅受氣動力和地心引力作用，則平面內(nèi)運動方程為

(10)

式中：r,V,γ以及m分別表示航天器的地心距、速度、航跡角以及質(zhì)量，D和L分別表示氣動阻力以及氣動升力，其表達(dá)式分別為

(11)

式中：S為氣動參考面積，ρ為大氣密度，定義ρ0為海平面大氣密度，β為指數(shù)因子，h為航天器地心高度，則ρ表達(dá)式為

ρ=ρ0exp(-βh)

(12)

(13)

式中：CD0和K分別表示零升阻力系數(shù)和誘導(dǎo)阻力因子。

則轉(zhuǎn)移時間最優(yōu)問題性能指標(biāo)可以記為

(14)

記狀態(tài)變量x=[r,V,γ]T，協(xié)態(tài)變量λ=[λr,λV,λγ]T，則哈密頓函數(shù)H為：

(15)

狀態(tài)與協(xié)態(tài)方程：

(16)

初始條件：

(17)

末端條件：

(18)

哈密頓函數(shù)H的末值約束：

H(tf)=0

(19)

(20)

則最優(yōu)控制律如下所示：

若λV<0，則

(21)

式中：CLS具體表達(dá)式為

(22)

若λγ≥0，則

(23)

3 利用UKF參數(shù)估計求解氣動輔助平面變軌問題

本節(jié)將求解最優(yōu)控制的兩點邊值問題改寫為參數(shù)估計問題，如果選擇估計參數(shù)為初始協(xié)態(tài)變量λ(t0)=[λr(t0),λV(t0),λγ(t0)]T以及轉(zhuǎn)移時間tf，則輸出末端約束條件有：r(tf)，V(tf)，γ(tf)以及H(tf)，若初始估計參數(shù)不當(dāng)，則出現(xiàn)碰撞地球的情況，速度接近為 0 m·s-1，導(dǎo)致積分無法進(jìn)行，為了避免上述情況，將積分停止條件設(shè)置為當(dāng)前速度到達(dá)目標(biāo)速度Vf，則上述問題可以降低一維，避免對轉(zhuǎn)移時間的猜測，減少求解問題的難度，則估計參數(shù)wk與期望輸出dk分別為

wk=[λrk(t0),λVk(t0),λγk(t0)]T

(24)

dk=[rk(tf),γk(tf),Hk(tf)]T

(25)

記xk為初始狀態(tài)變量，則dk是xk和wk的函數(shù)，可通過以xk和wk為初始條件，積分方程(16)求解，記這個函數(shù)為

dk=G(xk,wk)

(26)

將上述參數(shù)估計問題寫成狀態(tài)空間表達(dá)式：

(27)

式(27)表示了一個狀態(tài)轉(zhuǎn)移矩陣為單位陣的靜態(tài)過程，rk為過程噪聲，期望輸出dk則與wk滿足確定的非線性關(guān)系，ek為觀測噪聲。

則該參數(shù)估計問題就可以用各種濾波算法求解，這里直接給出基于UKF濾波器的參數(shù)估計算法流程，如圖2所示，其中：

(28)

式中:N是w的維數(shù)，η是尺度參數(shù)，ε決定了無損變換的σ點相對于w當(dāng)前均值的分布范圍，一般設(shè)為小量，取值范圍為[10-4，1]。常量κ一般取為0或3-N,β是與w先驗分布的相關(guān)常量，對于高斯分布，β=2是最優(yōu)的，ρRLS為遺忘因子，用于防止因模型誤差較大造成的濾波發(fā)散，其取值范圍為(0，1]。

濾波參數(shù)對濾波收斂的過程影響很大，如選取不合適，會導(dǎo)致濾波收斂過程的振蕩幅度很大，不能較快收斂，甚至發(fā)散。本文通過大量數(shù)值仿真，給出了可調(diào)參數(shù)的推薦值，如表1所示，其中尺度參數(shù)ε隨著濾波進(jìn)行，該值隨觀測誤差的減小而減小，有利于算法的快速收斂。

表1 濾波器參數(shù)Table 1 The value of UKF parameter estimation

圖2 UKF參數(shù)估計流程圖Fig.2 The flow chart of UKF parameter estimation

4 算例與分析

本節(jié)給出了一個利用UKF參數(shù)估計算法求解氣動輔助平面變軌的算例，算例中參數(shù)的具體數(shù)值見表2。為了增強(qiáng)算法求解效率，將各物理量進(jìn)行歸一化處理，其中，長度量綱選取106m，時間量綱選取103s。積分方法選擇取四階龍格庫塔方法，積分步長為0.5 s。

圖4 控制量變化曲線Fig.4 Optimal control vs time

圖5 哈密頓函數(shù)變化曲線Fig.5 Hamilton function vs time

對比文獻(xiàn)[3,7]的結(jié)果，兩者狀態(tài)變化趨勢是一致的，則驗證了該算法的正確性。圖5中哈密頓函數(shù)值在整個飛行過程中近似保持不變，則驗證了所得結(jié)果的最優(yōu)性。航天器飛出大氣層外在r2處速度為7495.7 m·s-1，航跡角為0.4°，則整個轉(zhuǎn)移過程中需要特征速度增量為霍曼轉(zhuǎn)移所需特征速度增量的49.9%，接近理想的49.3%。

為了研究UKF參數(shù)估計法的求解性能，將計算結(jié)果與hp自適應(yīng)Radau偽譜法(簡稱偽譜法)的計算結(jié)果作比較，具體結(jié)果見表3。

表3 優(yōu)化結(jié)果的對比Table 3 Comparation of the optimal result

圖6給出了兩種算法得到的控制律，其中，實線為UKF參數(shù)估計法的計算結(jié)果，虛線為偽譜法的計算結(jié)果。比較升力系數(shù)變化，兩種方法得到的控制律有一定差別，可以解釋為兩種方法得到不同的局部最優(yōu)解，具體表現(xiàn)為性能指標(biāo)的數(shù)值差別，偽譜法通過延長轉(zhuǎn)移時間，抬高了航天器的最低飛行高度，降低了峰值熱流。在實際工程應(yīng)用中，為了使航天器在極端條件不被燒蝕，必須考慮熱流的限制，保證航天器可靠性，可以作為下一步研究工作。從圖6可以看出，UKF參數(shù)估計法得到的控制律過度光滑，精度高，沒有出現(xiàn)虛擬數(shù)值震蕩。

在相同計算條件下(Intel G870 3.10 GHz)，UKF參數(shù)估計法求解效率優(yōu)于偽譜法，而且利用UKF參數(shù)估計法，算法結(jié)構(gòu)簡單，不需要強(qiáng)大的非線性規(guī)劃軟件包。

圖6 最優(yōu)控制律的對比Fig.6 Comparation of the optimal control

5 結(jié) 論

本文提出了一種基于UKF參數(shù)估計的氣動輔助變軌求解算法。該方法克服猜測協(xié)態(tài)變量初值的困難，協(xié)態(tài)變量初值可以隨機(jī)選取，同時避免了打靶法對梯度矩陣的推導(dǎo)，因此降低了求解氣動輔助變軌問題的難度。數(shù)值仿真表明：該方法結(jié)構(gòu)簡單，求解效率高，具有良好的魯棒性，可以作為求解氣動輔助變軌的一個有力工具。