劉 旭,李 響,張后軍,郭宇恒,王曉鵬
(1.北京理工大學宇航學院,北京 100081;2.北京機電工程研究所,北京 100083;3.上海機電工程研究所,上海 200233)
傳統的再入軌跡優化不考慮不確定性因素,即默認優化模型中所有參數都是準確的[1-3],而實際飛行過程中不確定性因素客觀存在,如飛行器氣動參數誤差、大氣參數誤差等,這些誤差的影響一般通過后續的控制系統設計消除。近年來,有學者在飛行軌跡優化階段就將這些不確定性因素考慮進去,提出了飛行軌跡的魯棒優化(Robust design optimization,RDO)[4]與可靠性優化(Reliability-based design optimization,RBDO)[5]。RDO在最小化目標函數期望的同時,盡量使目標函數對不確定性因素不敏感;RBDO著重于在不確定性因素影響下將失效概率限制在某個值以內。通過RDO或RBDO優化后得到的飛行軌跡對不確定性因素的敏感度降低,自身就有一定抗干擾的能力。
文獻[4]指出在考慮不確定性的情況下,相較于飛行器靜態魯棒優化(如氣動優化、結構優化等),飛行器動態魯棒優化(如飛行軌跡優化)的研究比較有限,因此提出一種不確定條件下的飛行軌跡魯棒優化方法;文獻[6]指出飛機軌跡規劃和預測中的不確定性給未來空中交通管理系統帶來了重大挑戰,因此提出一種考慮氣象不確定性的基于最優控制的飛機軌跡規劃方法;文獻[7]注意到不確定性(尤其是風場不確定性)對于無人機動態滑翔軌跡的影響,在量化不確定性后提出一種動態滑翔軌跡魯棒優化方法;火星大氣密度具有的強不確定性是火星進入段飛行任務的主要干擾源,對此文獻[8]提出一種火星大氣進入段的軌跡魯棒優化,得到的最優軌跡在不確定條件下的魯棒性明顯優于確定性優化,而文獻[5]則是提出用序列可靠性優化(Reliability-based sequential optimization,RBSO)來降低參數不確定性對火星進入段軌跡的影響,同樣得到優于確定性優化的結果;文獻[9]應用RBDO優化飛行器在不規則小行星軟著陸的軌跡,將不確定動態模型的軌跡優化轉化為帶有可靠性約束的參數優化,并應用序列優化和可靠性評估(Sequential optimization and reliability assessment,SORA)流程求解;文獻[10]利用脫敏最優控制(Desensitized optimal control,DOC),將標稱軌跡設計和標稱軌跡跟蹤結合在一起,得到的火星下降段脫敏軌跡對不確定性和擾動的敏感性大幅降低。
在飛行軌跡魯棒優化或可靠性優化的實際計算中,不確定性傳播(Uncertainty propagation,UP),即如何根據輸入不確定量計算得到輸出不確定量的統計特征(一般指均值和方差)是關鍵步驟。傳統蒙特卡洛(Monte Carlo,MC)仿真法易于操作,對各種復雜問題均有很強的適應性,得到了較廣泛的應用,也常用于檢驗其他方法的準確性。但MC仿真法的精度依賴于大量樣本點上的仿真次數,為得到可靠結果需要的高昂計算代價限制了它在軌跡魯棒優化中的應用。混沌多項式展開(Polynomial chaos expansion,PCE)是一種基于Askey正交多項式隨機展開的不確定性傳播分析法,能高效地計算輸出不確定性變量的統計特征[11],近年來在飛行軌跡魯棒優化中得到了較廣泛的應用[4,6-8]。基于可靠性的飛行軌跡優化[5,9]多采用最可能失效點(Most probable point,MPP)法進行不確定性分析,包括一次可靠度法(First-order reliability method,FORM)[12]和二次可靠度法(Second-order reliability method,SORM)[13],其基本思路是在MPP點進行局部展開來近似極限狀態函數,并以此計算失效概率。文獻[14]引入人工神經網絡模型,用訓練后的神經網絡作為復雜系統的替代模型,改善單獨使用MC法的不足。文獻[15] 通過Kriging代理模型來構建目標變量和不確定性輸入參數的響應關系,以此來分析目標變量的不確定性信息。
除了前述PCE方法與基于MPP的方法,協方差分析描述函數技術(Covariance analysis describing function technique,CADET)是另外一種高效的隨機系統統計特征計算方法。該方法最早由Gelb提出[16],先采用統計線性化方法將非線性系統轉化為近似線性系統,再用協方差分析法導出系統隨機狀態量均值和協方差的傳播方程,只需一次求解就能確定特定隨機系統的統計特征,在飛行器制導系統性能統計分析中得到了廣泛應用。文獻[17]應用CADET進行單發武器的射擊效能分析;文獻[18] 基于CADET建立了制導炸彈半實物仿真系統的誤差傳播理論,分析了仿真姿態角和姿態角速度的誤差對仿真結果的影響;文獻[19]基于閉環控制系統下線性化的相對運動動力學模型,采用CADET對定義航天器誤差橢球的協方差矩陣進行分析;文獻[20] 針對火星探測捕獲的誤差傳播,利用CADET分析單一誤差和綜合誤差對不同捕獲策略入軌精度的影響。現有的文獻中,CADET主要用于制導精度分析,本研究將CADET引入再入軌跡魯棒優化中,實現了再入軌跡魯棒優化的高效計算。
不失一般性,先以傳統確定性條件下Bolza型飛行軌跡優化問題構建飛行軌跡魯棒優化模型,可以描述為:將u(t)作為控制量,使得

(1a)

(1b)
C[x(t),u(t),t]≤0
(1c)
H[x(t0),u(t0),x(tf),u(tf),tf]=0
(1d)
式中:x(t)是狀態向量;u(t)是控制向量;f是動態約束向量函數;C是過程約束向量函數;H是端點約束向量函數。
現在考慮不確定性因素的影響,將隨機量s引入式(1),得到:將u(t)作為控制量,使得

u(t),t]dt
(2a)

(2b)
C[x(t,s),u(t),t]≤0
(2c)
H[x(t0,s),u(t0),x(tf,s),u(tf),tf]=0(2d)
由于引入了隨機量s,式(1)中的狀態量x(t)由確定量轉變為隨機量x(t,s)。相應地,式(1)中確定性的目標函數式(1a)、動力學方程式(1b)、過程約束式(1c)和端點約束式(1d)分別轉變為隨機條件下的式(2a)、(2b)、(2c)和(2d)。
由于含有隨機量s,式(2)中的目標函數式(2a)和約束條件式(2c-2d)與傳統軌跡優化中的目標函數式(1a)和約束條件式(1c-1d)有本質的不同。現有的軌跡優化方法都是針對如式(1)所示的確定性問題發展起來的,不能直接處理式(2)所示的含有隨機量的優化問題。為此,引入飛行器氣動和結構優化中魯棒優化的理念[21-22],將傳統確定條件下飛行軌跡優化模型轉化為如式(3)所示不確定條件下的飛行軌跡魯棒優化模型:將u(t)作為控制量,使得
min μ(J)+kJ·σ(J)
(3a)

(3b)
μ(C[x(t,s),u(t),t])+kC·σ(C[x(t,s),
u(t),t])≤0
(3c)
μ(H[x(t0,s),u(t0),x(tf,s),u(tf),tf])=0
(3d)
σ(H[x(t0,s),u(t0),x(tf,s),u(tf),tf])≤εH
(3e)
式中:μ(·)和σ(·)分別表示隨機量的均值和標準差;kJ和kC是權重系數;εH是隨機端點約束標準差的容許值。
式(3)所表達的具體含義如下:
式(3a)與目標函數式(1a)對應,是隨機目標函數的均值μ(J)和標準差σ(J)的加權和,通過最小化隨機目標函數的均值來表達確定性目標函數式(1a)的設計意圖,同時通過最小化隨機目標函數的標準差來限制它的離散程度。目標函數的最優性和魯棒性可以通過選取權重系數達到設計要求,理論上kJ值越大,目標函數的離散程度越小。
式(3c)與過程約束式(1c)對應,式(1c)是不等式約束,在魯棒優化中,不等式約束通過隨機約束的均值μ(C)和標準差σ(C)的加權和來表達,kC越大,所得最優軌跡的魯棒性越高,滿足過程約束的概率越大,但求解難度也會相應增加。
式(3d)、(3e)與端點約束式(1d)對應,式(1d)是等式約束,在魯棒優化中,等式約束通過隨機約束的均值式(3d)和標準差式(3e)兩組約束來表達,這與不等式約束是不同的。均值約束要求等于原問題中的值,而標準差約束要求小于εH,εH越小,魯棒性越高。
總體來看,魯棒優化模型式(3)通過要求隨機優化模型式(2)中隨機量的均值滿足確定性優化模型式(1)中相應的約束,來表達原有的設計意圖;同時通過約束隨機優化式(2)中隨機量的標準差,限制其離散程度,提高最優軌跡的魯棒性。這樣從統計角度來看,式(3)得到的最優軌跡在隨機擾動的影響下偏離設計值的程度要小于式(1),從而達到魯棒優化設計的目標。
值得注意的是,魯棒優化模型中式(3b)目前仍是含有隨機量s的動力學方程,現有軌跡優化方法(如偽譜法等)難以直接處理這類問題,在下節中將討論如何利用協方差分析描述函數技術(CADET)將含有隨機量的魯棒優化問題式(3)轉化為等效的確定性優化問題,從而可以采用偽譜法等現有軌跡優化方法求解。
對于式(3)所示的魯棒優化模型,如何高效計算式(3a,3c-3e)中的統計特征(均值μ(·)和標準差σ(·))是關鍵步驟。蒙特卡洛法程序結構簡單,易于實現,得到了廣泛應用。根據大數定律,樣本數N→∞時,蒙特卡洛法的結果以概率1收斂于真實值。一般認為,樣本點數量達到104時,蒙特卡洛法的結果才是穩定收斂和可靠的[23],工程上通常要成千上萬樣本點才能得到滿足精度要求的結果。對單次隨機過程的統計特征分析而言,蒙特卡洛法是可以接受的;而軌跡優化計算是一個迭代過程,在每一迭代步都需要計算均值μ(·)和標準差σ(·)的值,在這種情況下,蒙特卡洛法將導致巨大的計算量,在實際計算中是不可行的。根據上述情況,本文采用CADET計算隨機系統的統計特征,并將含有隨機量的式(3)轉化為等效的確定性優化模型。
一般連續時域非線性隨機系統可以表示為:

(4)
式中:x(t,s)是n維隨機狀態向量;w(t)是m維強迫函數向量,包含作用于系統的隨機擾動s和控制輸入u;f(x,t)是x(t,s)的非線性矢量函數;G(t)是n×m的連續函數矩陣。
協方差分析描述函數法,就是先利用描述函數理論對非線性系統進行統計線性化,然后對線性化后的系統進行協方差分析。所謂統計線性化,是根據狀態量x(t,s)的概率密度函數形式,在較大的x(t,s)變化范圍內,用一個線性函數來逼近非線性函數f(x,t)。其本質是對系統的每一個非線性用一個或幾個近似的、對輸入幅度靈敏的線性增益來代替[24]。這種方法的優點是不要求f(x,t)連續可微,使許多不能用真線性化進行線性化處理的實際函數有了線性化的可能[25]。

(5)

(6)

(7)
式中:b(t)是w(t)的確定性分量,Q(t)是w(t)隨機分量的譜密度矩陣。只需確定初始條件m(0)和p(0),利用式(7)進行一次計算就能獲得任意時刻系統狀態量的統計特征,而且精度與幾千次的蒙特卡洛仿真相當,極大程度上節省了計算時間與成本。
由式(6)可知,描述函數的計算依賴于狀態量的概率密度函數,一般來說十分復雜,甚至不能求得解析式。工程上常假定x(t,s)服從聯合正態分布,從而描述函數的計算可簡化為
(8)
這里假設隨機狀態量x(t,s)服從正態分布,只有對高斯輸入的線性系統才是嚴格成立的,而對非高斯輸入的線性和非線性系統往往是近似的。對于非線性系統,盡管輸入是高斯的,輸出一般來說是非高斯的。由中心極限定理可知:當經過低通濾波器后,隨機過程趨于高斯過程。因此,當信號通過系統傳遞時,可以借助系統的線性部分來保證非高斯的非線性輸出結果近似是正態的[24]。前述高斯假設的有效性已經得到了廣泛的研究和檢驗。
在附錄A中,給出了采用CADET方法得到的再入軌跡非線性動力學模型的統計線性化模型。
式(3)所示的軌跡魯棒優化模型盡管表達了均值、標準差等引入隨機量后的性能信息,但式(3b)仍是關于x(t,s)的隨機微分方程,含有隨機量s,且其余各式中關于均值和標準差的求解方法也沒有給出,式(3)是無法用于實際計算的。現在應用CADET,將式(3)中含有的隨機變量s消去,得到等效的確定性優化模型,具體步驟為將式(3b)轉化為關于x(t,s)均值和協方差的常微分方程,同時將式(3)中目標函數、過程約束和端點約束表示為關于x(t,s)均值和協方差的函數,從而得到等效的確定性軌跡魯棒優化模型:將u(t)作為控制量,使得
minJμ(m(t),p(t),u(t),t)+kJ·Jσ(m(t),
p(t),u(t),t)
(9a)

(9b)
(9c)
Cμ(m(t),p(t),u(t),t)+kC·Cσ(m(t),
p(t),u(t),t)≤0
(9d)
Hμ(m(t),p(t),u(t),t)=0
(9e)
Hσ(m(t),p(t),u(t),t)≤εH
(9f)
式中:m(t)和p(t)分別是隨機狀態量x(t,s)的均值和協方差,根據CADET中的式(7)計算;(·)μ和(·)σ分別表示用于計算(·)的均值和標準差的函數。
對比式(9)和式(3),可以發現,采用CADET后,含有隨機量s的魯棒優化模型式(3)轉化為不含隨機變量的等效確定性優化模型式(9),應用現有的軌跡優化算法完全可以求解。換言之,通過求解確定性優化問題式(9),就能得到不確定條件下,魯棒優化意義下的最優控制量和最優軌跡。
文獻[26]中的再入軌跡優化問題是一例經典的軌跡優化問題,希望在滿足終端和過程約束的條件下,最大化終點的橫向距離(即緯度),在確定性條件下可以描述為:將α(t)和σ(t)作為控制量,使得
maxJ=φ(tf)
(10a)
(10b)
(10c)
(10d)
式中:r是到地球球心的距離;θ是經度;φ是緯度;v是速度;γ是航跡角;ψ是航向角;α是攻角;σ是傾側角;m是質量;μ是地心引力常數;Re是地球半徑;D是阻力加速度;L是升力加速度;q是熱流密度。式(10c)表示端點約束,式(10d)表示過程約束。關于該模型中各參數的詳細說明可參見文獻[26],這里不再詳述。
進一步,阻力加速度和升力加速度表示為:
(11)
式中:CD和CL分別是阻力系數和升力系數;ρ是大氣密度S是特征面積;h是高度;H是密度尺度高;ρ0是海平面大氣密度。
在傳統飛行軌跡優化中,氣動參數CD·和CL·視為確定值,根據實驗得到的氣動數據求得。
由于制造工藝誤差、測量誤差等各方面因素,將氣動參數CD·和CL·視為不確定性量更為合理。本文中假設:
(12)
式中:(·)N表示(·)的參考值;s1,s2表示隨機變量。本文仿真實驗中,s1~N(0.01,0.012),s2~N(-0.01,0.012)且s1,s2相互獨立。需要說明的是,這里s1,s2的均值不為0,目的是使干擾條件更苛刻,以顯示方法的有效性。關于s1,s2的概率分布規律等因素本身也是值得研究的問題,本研究著眼于再入軌跡魯棒優化方法,因此在這里對此不作過多討論。
引入隨機變量s1,s2后,確定性軌跡優化問題式(10)轉化為如式(2)所示的隨機軌跡優化問題,在此基礎上進一步轉化為式(3)所示的軌跡魯棒優化問題,具體如下:
式(10)中的目標函數式(10a)變為:
maxJ=μ(φ(tf,s))-kJ·σ(φ(tf,s))
(13)
動力學方程式(10b)轉化為含隨機變量的動力學方程:

(14)
式中:x(t,s)表示隨機狀態向量;f表示式(10b)的右端項函數。
相應地,式(10)中的確定性端點約束式(10c)和過程約束式(10d)轉化為用均值和標準差描述的統計意義上的約束。其中,等式約束式(10c)轉化為由均值和標準差表示的兩組約束式(15a)和式(15b),不等式約束(10d)轉化為式(15c):
(15a)
(15b)
(15c)
需要說明:式(15b)僅對原確定性問題中有終端約束的狀態變量施加了標準差約束;所考慮的氣動不確定性不會影響到初始狀態,因此式(15b)中沒有設置對初值標準差的約束;此外,用來限制離散程度的標準差容許值和權重系數由設計者根據實際情況設定,本研究中設置如下:εhf=4500 m,εvf=250 m/s,εγf=3°,kJ=3,kCr=3,kCθ=3,kCv=3,kCγ=3,kCq=3。
綜上,再入軌跡魯棒優化模型(RO)表示為:將α(t)和σ(t)作為控制量,使得
(16)
至此,已構建起2個再入軌跡優化問題,即傳統確定性優化問題DO(式(10))以及同時考慮均值和標準差約束的魯棒優化問題RO(式(16))。在下節中,將分別求解這兩個優化問題,并對結果進行對比分析。
魯棒優化模型式(16)含有式(14),而式(14)是含有隨機變量s的隨機飛行動力學方程,在實際計算求解式(16)時,需要采用CADET方法將式(14)轉化為形如式(7)的關于均值和協方差的確定性動力學方程,具體的轉化結果可參考附錄A。
3.3.1最優控制量的對比分析
采用偽譜法求解這兩個軌跡優化問題(關于偽譜法的相關論述這里不再展開,可參看相關文獻[27-29]),得到傳統優化DO的最優控制量:αDO(t)和σDO(t),魯棒優化RO的最優控制量:αRO(t)和σRO(t),如圖1所示,為了更清楚地顯示,兩個小方框給出了相應時間段的放大圖。

圖1 最優控制量(方框內為放大圖)Fig.1 Optimal control
觀察圖1中控制量的時間歷程曲線,可得如下信息:
魯棒優化最優控制量αRO(t)、σRO(t),分別與對應的確定性優化的最優控制量αDO(t)、σDO(t)的時間歷程趨勢基本一致,但αDO(t)和σDO(t)相對平緩,而αRO(t)和σRO(t)的波動較大,且越接近終端時刻波動越明顯。
αRO(t)和σRO(t)所需的時間比αDO(t)和σDO(t)的更長(RO為2289.4 s,DO為2200.1 s)。
上述現象表明,魯棒優化為了降低不確定性,比傳統確定性優化付出了更多的控制成本,這是符合客觀實際的。此外,需要說明的是,圖1顯示攻角在末端變化較快,這是圖形比例原因導致的(總的飛行時間較長),攻角實際最大變化率絕對值為3.43°/s,工程上是完全可以實現的。
3.3.2目標和約束的魯棒性對比分析
為驗證兩種優化模型得到的最優控制量在不確定因素影響下的效果,將αD0(t)和σD0(t)、αRO(t)和σRO(t)分別代入式(14)所示的含隨機變量的動力學方程中,進行10000次的MC仿真,并對仿真結果進行對比分析,主要從三個方面進行:目標函數的魯棒性對比;終端約束的魯棒性對比;過程約束的魯棒性對比。
首先給出相應的統計結果:
目標函數終端時刻緯度的統計特性及散布情況如表1和圖2所示;終端約束終端時刻高度、速度和航跡角的統計特性及散布情況如表2、表3和圖3所示;過程約束熱流密度約束的情況如圖4和表4所示。圖2和圖3中橫坐標表示MC仿真次數,對圖表的解釋及分析見后。

表1 終端時刻緯度的統計特性Table 1 Statistical property of latitude at terminal time

圖2 終端時刻緯度的散布情況Fig.2 Dispersion of latitude at terminal time

表2 終端時刻高度、速度和航跡角的均值和相對誤差Table 2 Mean and relative error of height,speed and flight path angle at terminal time

表3 終端時刻高度、速度和航跡角的標準差Table 3 Standard deviation of height,speed and flight path angle at terminal time

圖3 終端時刻高度、速度和航跡角的散布情況Fig.3 Dispersion of height,speed and flight path angle at terminal time

圖4 熱流密度的散布情況Fig.4 Dispersion of heat flux density
通過分析上述圖2到圖4中的結果及表1到表4中的數據,可以得到以下結論:

表4 熱流密度約束的滿足情況Table 4 Satisfaction of heat flux density constraint
1)目標函數結果對比
表1中第一、二行數據顯示,對于目標函數φ(tf)的均值,DO大于RO,而對于φ(tf)的標準差,RO小于DO。
φ(tf)的均值越大而標準差越小越好,因為極大化φ(tf)是設計目標,而較小的標準差則意味著抗干擾能力較強。因此上述現象表明,RO為提高系統對不確定性因素的魯棒性,犧牲掉一部分φ(tf)的均值,這是符合設計預期的。
為了顯示出均值和標準差的綜合效果,表1中第三行數據顯示了均值和標準差的加權和,kJ取3時,RO表現優于DO。
圖2顯示了氣動不確定條件下,終端時刻緯度的散布情況,可以直觀地看出RO中的樣本點分布更為緊密,表明RO的目標函數對不確定性因素的敏感程度更低。
2)終端約束結果對比
表2中數據顯示,在氣動不確定因素影響下,DO的終端高度、速度和航跡角的均值與約束要求的值有較大偏差,最大達到51.40%。而RO的終端狀態均值與終端約束要求的值基本吻合。這說明,在氣動不確定因素影響下RO能使終端狀態以約束要求值為中心散布,而DO的終端狀態散布偏離了約束要求值。
表3中數據顯示,RO中的高度和航跡角的標準差小于DO中對應量的標準差,這說明,在氣動不確定影響下,RO中的高度和航跡角的分布更緊密。此外,盡管DO中速度的標準差小于RO中速度的標準差,但由于DO中速度的均值偏離了要求的均值,較小的標準差并未帶來實際的收益,這在圖3中得到了體現,具體如下。
圖3表明了氣動不確定條件下,高度、速度和航跡角終端狀態的散布情況,紅實線是設計要求,表示均值,上、下紅虛線表示均值加減標準差,則落入上、下虛線之間部分的樣本點越多,說明抗干擾能力越強,魯棒性越好。DO的10000個樣本點中分別有5937(高度)、6599(速度)和6310(航跡角)個在這個范圍內,占比為59.37%、65.99%和63.10%;而RO的10000個樣本點中,分別有6883(高度)、6766(速度)和8267(航跡角)個在這個范圍內,占比為68.83%、67.66%和82.67%。顯然,RO所得到的結果比DO的魯棒性更好。
3)過程約束結果對比
本算例中共有7個過程約束,其中熱流密度的過程約束是主動約束,因此這里以熱流密度過程約束來討論不確定因素的影響。
圖4a和圖4b分別顯示了DO和RO中q(t)時間歷程樣本的分布情況,紅線代表飛行過程中熱流密度的約束上限。對比圖4a和圖4b,直觀地顯示出DO中大量熱流密度時間歷程樣本超出了約束限制,而在RO中,絕大多數樣本都是滿足約束的。表4中給出了具體值,在DO的10000條q(t)時間歷程樣本中,有17條完全滿足約束,占比為0.17%;而在RO的10000條q(t)時間歷程樣本中,有9933條完全滿足約束,占比為99.33%。上述結果表明,在氣動不確定條件下,RO中的熱流密度過程約束有極強的抗干擾能力,而DO中的熱流密度過程約束的抗干擾能力很弱。
總的來說,在不確定性因素的影響下,再入軌跡魯棒優化的最優軌跡相較于確定性優化的最優軌跡表現出更強的抗干擾能力,而對于不同的物理性能指標,抗干擾的能力有所不同,在本算例中,對熱流密度的抗干擾能力最強。
傳統的再入軌跡優化是在確定性條件下展開的,本研究將不確定性參數引入再入軌跡優化中,構建了再入軌跡魯棒優化模型,使得到的飛行軌跡具有一定的抗干擾能力。將CADET作為不確定性分析手段,引入再入軌跡魯棒優化的實際計算中,實現了再入軌跡魯棒優化的高效計算。仿真結果表明,相比于確定性軌跡優化,本文方法得到的最優軌跡在不確定性因素的干擾下具有更強的魯棒性。
本研究的著眼點是再入軌跡魯棒優化方法,而軌跡優化中的不確定性來源眾多,將來還應對輸入不確定性的合理建模展開深入研究。此外,協方差分析法的精度會隨著系統非線性程度的增加而降低,對于非高斯輸入的有效性也有待考察,因此準確高效的不確定性傳播方法依然是未來關注的重點。