馮子鑫 薛文超 張 冉 齊洪勝
近些年來,由于可以顯著降低太空探索的成本,關于可回收火箭的研究受到廣泛的關注.20 世紀50~60 年代,美國在航天飛機上首次實現航天運載器的部分重復使用,隨后進行大量的可重復使用運載器的研究[1].21 世紀以來,SpaceX、Blue Origin、Maste Space System 等商業公司同樣開展大量關于可回收火箭的研究[2-3].動力下降作為可回收火箭中的一個重要問題,要求利用稠密大氣的空氣阻力和火箭發動機的推力來實現火箭的定點軟著陸,是火箭著陸的最后一環.由于環境存在許多干擾,在動力下降段中常常存在大量的不確定性,尤其是未知的風場,這給動力下降段的制導控制帶來挑戰.
目前針對火箭動力下降問題的制導方法主要可以分為三類: 解析制導方法,基于軌跡優化的制導方法,基于學習的制導方法.解析制導方法是指制導指令由火箭狀態量解析表示的一系列方法,包括多項式制導[4]、重力轉彎制導[5]和近似最優解析制導[6](零控位置/速度誤差制導)等.這些制導方法具有形式簡單、易于實現等優點,已經在阿波羅登月、嫦娥登月等實際場景中得到應用.但是為了獲得解析的制導指令,往往使用較為簡單的問題模型,難以處理大氣層內的動力下降問題.基于學習的制導方法[7-8]則是近年來隨著深度學習興起而發展起來的方法,通過訓練可以顯著提高方法對初始狀態和模型不確定性的適應性.但是基于學習的方法具有難以設計獎勵函數、可解釋性差等缺點,因此目前還難以在實際場景中進行應用.
基于軌跡優化的制導方法[9-10]將控制和狀態作為待優化變量,在一定的性能指標下通過求解優化問題來獲得制導指令,其主要分為間接法和直接法.間接法是將動力下降問題建模為最優控制問題,通過極大值原理將最優控制問題轉化為兩點邊值問題進行求解[11-12].但是對于較為復雜的問題模型來說,求解轉化后的兩點邊值問題同樣較為困難,因此間接法更多只能用在一些較為特殊的、化簡后的模型.直接法則是指通過對原始最優控制問題進行直接描述和尋優的方法.凸優化是一種在動力下降問題中得到廣泛應用的直接法,基于凸優化的軌跡規劃將原始的非凸問題轉化或近似轉化為凸問題,再通過離散化將無窮維的連續系統轉化為有限維的離散系統,最后通過成熟的凸優化求解算法進行快速求解.文獻[13-14]首次在動力下降問題中提出無損凸化的概念,通過引入松弛變量將控制幅值約束轉化為凸約束并推廣到更一般的情形,從理論上嚴格證明轉化前后最優解的等價性.文獻[15]則將上述結果推廣到同時具有控制和狀態約束的情況,同樣證明轉化前后最優解的等價性.但是上述無損凸化的方法只對特定的約束形式成立,難以推廣到動力著陸段的一般情形.文獻[16]提出序列線性化的方法,針對非凸的、一般的非線性不等式,每次迭代使用上一次迭代的軌跡并在其附近進行一階近似展開,求解近似的凸優化問題,使用信賴域技術,通過多步迭代使近似最優解收斂.隨后同樣有一系列工作使用凸優化作為軌跡優化的工具[17-19].目前,基于凸優化的軌跡規劃方法已經被廣泛應用于可回收火箭的動力著陸問題中.但是上述方法針對的均是確定性問題,未考慮不確定性的影響.
對于軌跡規劃中不確定性的研究則相對較少.一種考慮不確定性的方法是基于敏感度矩陣(系統狀態關于未知參數的雅可比矩陣)的方法[20-23].文獻[20]將敏感度矩陣作為增廣狀態并在原始的最優控制問題中增加關于增廣狀態的初始和終端約束,通過求解增廣的最優控制問題來進行魯棒軌跡規劃.文獻[21]則設計一個基于敏感度矩陣和風場觀測器的臨近最優制導律,并且為處理真實風場的不確定性,加入擾動觀測補償器以提升魯棒性.文獻[22]考慮存在參數不確定性的軌跡跟蹤問題,通過最小化狀態敏感度和輸入敏感度進行魯棒軌跡規劃.文獻[23]將敏感度矩陣作為增廣狀態,對不等式約束進行一階近似并使用敏感度矩陣計算一階項,將原問題轉化為一個易求解的優化問題.此外,基于多項式混沌展開(Polynomial chaos expansion,PCE)的方法[24]則是軌跡規劃中處理不確定性的一種典型方法,其主要思想為,將存在不確定性的隨機過程展開為一系列正交的隨機多項式基底函數的加權和,從而建立不確定輸入與不確定輸出的映射關系.文獻[25]則將PCE 與凸優化的方法相結合,在僅考慮初值不確定性的情況下,通過序列線性化的方法降低求解問題的計算代價.文獻[26]通過引入譜分解技術來自適應更新基函數,以處理長時間的優化場景可能導致PCE 方法發散的問題.文獻[27]通過特征選擇和序列采樣提出序列增強PCE 方法,相比于原始的PCE 方法具有更高的計算效率和精度.
可以看出,現有對存在不確定性的動力下降問題或者軌跡規劃的研究要么對不確定性(或者風場)的建模較為簡單[23-24],要么使用隨機的PCE 方法,通過期望與方差等統計特性來衡量性能指標.這種方法需要把原始的隨機問題轉化為更高維的確定性問題,大大增加求解的復雜度.此外,PCE 方法的最優指標是在期望和方差等統計學意義下的平均最優,難以應對極端場景.本文則使用一種確定性的魯棒優化(Robust optimization,RO)方法[23],該方法轉化后的問題維度遠低于PCE 方法,并且考慮在參數最差時的最優性,可以處理更加極端的情況.本文的主要工作包括: 首先,建立更加貼近實際場景的風場模型并在確定性框架下給出動力下降的RO 問題;其次,使用一種對不等式約束采取一階近似的RO 方法,得到一個可以易于求解的單階段RO 算法;然后,通過理論分析單階段RO 算法在終端約束范圍較小時可能無解的缺陷,進一步提出一種改進的多階段RO 算法;最后,通過仿真說明多階段RO 算法對未知不確定風場同時具有較高的最終狀態精度和較強的魯棒性.
大氣層內的可回收火箭的動力下降問題是指火箭在發動機推力、重力和空氣阻力的影響下從高空狀態被引導到指定的著陸點.在本文中,x軸、y軸和z軸分別指向東、北和天空.為了簡便,文中的下標 (·)x、(·)y和 (·)z分別代表對應量在x軸、y軸和z軸上的分量.考慮如下大氣層內火箭的動力學模型
其中,r=[rx,ry,rz]T∈R3,v=[vx,vy,vz]T∈R3分別是位置向量和速度向量,g代表重力加速度,m∈R是火箭質量,κ∈R代表火箭恒定的燃料消耗系數,u=[ux,uy,uz]T∈R3是系統的控制輸入,代表火箭引擎提供的推力,其滿足如下幅值約束
其中,D(rz,v,vw)代表空氣阻力,vw(rz)∈R3代表風場向量,ρ(rz)∈R是大氣密度,Sref∈R,CD∈R分別代表火箭的參考面積和空氣阻力系數.
設初始時刻為t0=0,該問題的初始條件為
終端時刻tf自由,在終端時刻需要滿足終端條件
rf,vf代表期望的終端位置和速度.性能指標為終端時刻和推力二次型積分之間的加權和
其中,γ代表加權權重.
記X=[rT,vT,m]T,X0=[,,m0]T,那么動力學模型和初值可以寫為如下緊湊形式
標稱優化(Nominal optimization,NO)問題1.給定初始狀態X0以及目標終端狀態Xf,選擇合適的控制律u,在vw=(rz)的情況下,滿足式(3)、式(5)和式(6)的同時最小化性能指標(4).
NO 問題使用標稱風場信息進行計算,因此可以直接使用標準優化算法進行求解.基于NO 問題可以得出如下NO 算法.
NO 算法.輸入X0,rf,vf,求解NO 問題,輸出最優控制和相應的狀態軌跡與終端時刻.
第1.1 節中的NO 問題使用標稱風場,但在實際場景中,風場具有很強的不確定性.為了在模型中考慮不確定性的影響,在風場模型中引入不確定參數.
假設風場在y軸方向存在不確定性,不確定性風場為是未知參數,d代表未知參數的維度.具體而言
未知參數s的取值范圍滿足
此外,由于實際場景中z軸方向的風速可以忽略不計,因此假設z軸方向不存在風.對于x軸方向的風場,為了簡單起見,本文假設x軸方向的風場是確定已知的,滿足
其中,k,b均為已知參數.最后,記vw(rz,s)代表上述定義的、包含不確定參數s的風場.
在第1.2 節給出的不確定風場vw(rz,s)下,考慮如下魯棒優化問題:
魯棒優化(RO)問題2.
其中,Lri,Lvi代表事先確定的終端約束范圍.
注1.RO 問題2 在滿足對任意參數都不違背約束條件的可行解中尋求最小化優化指標的解,這使得RO 問題2 的最優解即使在最差的參數條件下也能擁有較好的效果.
注2.注意到RO 問題2 的終端約束和NO 問題1 的終端約束不同,這是因為在風場存在不確定性的情況下,使火箭嚴格滿足給定的終端條件式(3)是難以實現的,而要求火箭盡可能落在給定終端狀態的一定范圍之內更加具有可行性.
注3.在式(9)中選擇不同的約束范圍Lri,Lvi,是因為第1.2 節定義的風場不確定性出現在y軸方向,所以相比于y軸方向而言,火箭由風場不確定性導致的、在x軸和z軸方向的偏差較小,自然也應該取不同的約束范圍.
一般而言,直接獲得RO 問題2 的最優解是十分困難的,特別是當系統的動態方程較為復雜時.接下來將簡要介紹一種對不等式約束采用一階近似的RO 方法,該方法可以將RO 問題2 近似轉化為一個易求解的最優控制問題,并且轉化后的問題維度低于PCE 方法.該方法更詳細的推導過程見文獻[23].考慮如下一般形式的非線性優化問題
在上述優化問題中,通過等式約束確定一個X與u,s之間的隱式關系式,不妨設為
將F(X,u,s)=0兩端同時關于s求微分,可得Xs(u,s)滿足
其中,變量下標代表關于該變量求偏導數.
另一方面,在隱式關系X=X(u,s)下,令
其中,Sτ由式(8)定義.將(u,s)一階泰勒展開可得
其中,?s代表關于s的梯度算子,〈·,·〉 代表內積.因此
其中,等號成立是因為赫爾德不等式(Holder inequality).為了求取?s(u,),將式(12)關于s求偏導可得
其中,ei∈Rm是m階單位矩陣的第i列,m是G的維數.將式(14)代入式(13)中,再結合式(11)和原始問題式(10),可以得到原始問題在一階近似意義下的魯棒優化問題
其中,gi代表G的第i個元素,FX,Fs,GX和Gs均使用參數標稱值計算,Xs是新的待優化變量,代表通過等式約束F(X,u,s)=0和式(11)所導出的X關于s的偏導數.
注4.式(15)中不等式約束的第二項可以看作是式(10)的不等式約束對應于參數變化的“安全裕量”,它代表gi對于參數改變的敏感度.雖然在計算中使用的是參數標稱值,但是通過在式(15)中添加安全裕量來盡可能保證實際的參數不會導致優化結果違背不等式約束.
文獻[23]指出,雖然上述魯棒優化方法不能保證不等式約束G(X(u,s),u,s)≤0對?s∈Sτ都一定成立,但是在一定的連續性假設下,它能保證超出約束范圍的值具有一個關于τ2的上界.
為了敘述簡便,令
?v代表風和火箭的相對速度,分別代表 ?v在x,y,z軸上的分量.由于風場在z軸方向風速為零,因此 ?vz=vz.另外,令
那么式(2)中的空氣阻力可以表示為
此外,對于一個矩陣P,引入符號P(i:j,k:l)代表由P的第i行到第j行、第k列到第l列組成的子矩陣,P(i:j,:)代表由P的第i行到第j行構成的子矩陣.將上述魯棒優化方法應用到RO 問題2 中,可以得到:
RO 問題3.
RO 問題3 使用標稱風場信息進行計算,因此可以直接使用傳統優化算法求解.基于RO 問題3可以得出如下單階段RO 算法.
單階段RO 算法.輸入X0,rf,vf,Lri,Lvi,i=x,y,z,求解RO 問題3,輸出最優控制、相應的狀態軌跡與終端時刻.
第1.3 節得到了參數的變化導致的終端約束函數變化,即注4 所說的安全裕量,也即式(17)中的‖rsi(tf)A‖1,‖vsi(tf)A‖1.由式(13)的推導可知,這種方法在每個時刻總會選取最差的參數(盡管實際場景中很難出現這種情況),并且隨著時間推移,這種最差的參數對約束函數帶來的變化也會逐漸累積,最終會導致安全裕量比較大.這就導致RO 問題3 出現矛盾: 一方面,終端約束范圍Lri和Lvi,i=x,y,z越小越好,因為這代表著能使火箭更加精準地以給定速度落到給定地點;另一方面,如果Lri和Lvi比安全裕量還小,那么式(17)中的不等式約束將無法滿足,導致優化問題無解,因此獲取更多安全裕量的信息十分重要.事實上,本節將會定量地給出安全裕量的一個上界.為此,首先給出一個假設和一個引理.
在假設1 下,可以證明有如下引理成立.
引理1.若假設1 成立,在RO 問題3 中,Xs等式約束內的fX(4:6,:)和fs(4:6,:)滿足
由式(26)和假設1,有
其中,第二個不等式是因為式(7).因此
對于fX(4:6,4:6),可以得到
對于fs(4:6,:),在式(23)中,如果記
那么可以得到
引理1 說明式(17)中的fX,fs都是有界的,因此可以通過=fXXs+fs求出rsi,vsi,i=x,y,z的上界,進而求出式(17)中安全裕量的界.具體而言,給出如下定理.
定理1.若Xs滿足
那么對?tf ≥0,式(17)中 的‖rsi(tf)A‖1和‖vsi(tf)A‖1有界,滿足
其中,i=x,y,z,C=1+(tf+1)M1+2M2.
證明.定理中的結果需要證明rs和vs所有行向量的1 范數有界,因此只需證明 ‖rs‖∞和‖vs‖∞有界.為了證明 ‖rs‖∞和‖vs‖∞有界,首先證明它們的Frobenius 范數(F 范數)有界,再通過矩陣范數之間的關系(對任意n階方陣A,有‖A‖∞≤得到‖rs‖∞和‖vs‖∞的界.為了簡便,在本證明中,令 |·|代表向量或者矩陣的2范數,對于標量則代表取絕對值,并記rsj=rs(:,j),vsj=vs(:,j),j=1,2,3 代表rs和vs的各個列向量,注意和rsi,vsi區分,i=x,y,z.它們代表rs和vs的各個行向量.將式(19)拆為三個列向量對應的微分方程可得
兩端同時取2 范數并由引理1 可得
或者等價地
〈·,·〉 代表內積.對上式兩端在 [0,t] 上積分,注意到sj(0)=0,可得
對于上式中的 |rsj|2,有
將上式代入式(30)中,得到
由Gronwall 不等式可得
回顧vsj代表vs的列向量,由F 范數的定義可知
又因為矩陣的2 范數總是小于矩陣的F 范數,所以
由式(33)即可得知定理1 的式(29)成立.對于rs,由式(31)可得
類似于式(32)和式(33)的推導,有
由式(34)即證定理1 的式(28)成立.
從定理1 可以看出,式(28)~(29)中給出的界是隨時間指數增長的,對于一些需要較長時間才能完成的優化場景,顯然會導致必須選取較大的終端約束范圍優化問題才能有解,但是大的終端約束范圍又會降低落點的精度.對于優化時間較長的場景,為了通過縮短優化時長使問題具有可行性,本文提出一種多階段RO 算法.其主要思想是通過將整個優化進程拆分為多個串聯的子階段,在每個子階段中分別求解RO 子問題,縮短優化時長,進而避免優化時間較長致使RO 問題3 無解.
具體而言,在如圖1 所示的K階段RO 算法流程圖中,選取每個階段的預估終端時刻t1,t2,···,tK,使得每個階段的預估時長為標稱軌跡時長的1/K.將標稱軌跡在t1時刻的狀態rNO(t1)、vNO(t1)作為第一階段的目標終值,將初始狀態X0作為第一階段的初值,使用單階段RO 算法進行優化,將最優解輸入到火箭中,得到火箭在第一階段末的位置在第二階段中,則以為第二階段的初值,以rNO(t2),vNO(t2)作為第二階段的目標終值,同樣使用單階段RO 算法求解,如此進行下去,直到第K階段結束.需要注意的是,最后一個階段的目標終值滿足

圖1 多階段RO 算法流程圖Fig.1 Flow chart of multi-stage RO algorithm
算法1.多階段RO 算法
注5.關于階段數的選擇,在實際場景中應綜合考慮多方面因素.一方面,階段數量越多,每個階段的優化時長就越短,根據本節的分析,就可以求解更小終端約束范圍的優化問題,進而使最終狀態更加精確.另一方面,階段數量進一步增加也會帶來一些問題.首先,階段數增多會使每個階段的優化時長縮短,縮短優化時長是為了縮小約束范圍Lri,Lvi,i=x,y,z.但是從定理1 可以看出,當優化時長tf減小并趨于零,式(17)中的安全裕量同樣會減小并趨于零,若Lri,Lvi同樣減小并趨于零,此時RO 問題3 將近似等價于NO 問題,這表明當優化時間過短時,RO 算法相比于NO 算法便不再具有優勢.其次,當階段數增加時,整個優化進程需要求解的優化問題也會增加,這對計算性能有著更大的需求.因此,在實際場景中,設計算法時應綜合考慮火箭初始狀態、能夠容許的終端約束范圍大小、計算能力等條件來選擇階段數.本節算法將每個階段的預估優化時長取為相等,正是考慮到每個階段的優化時長不能過長、也不能過短的因素.
首先說明仿真所使用的風場模型.本文通過實際測量出的6 個時刻風場數據擬合出6 個風速與高度之間的函數關系,將這些擬合函數各個參數對應的最大、最小值的平均值作為參數標稱值.x軸方向的風場參數選為擬合出的標稱參數.對于y軸方向的不確定風場,將各個擬合參數對應的最大、最小值分別作為該參數取值范圍的上、下界,擬合出的y軸風場與高度之間的函數關系如圖2 所示.如果令風場序號0 代表標稱風場,風場序號1 到6 代表實際風場,那么各個實際風場的參數與標稱參數在無窮范數意義下的距離如圖3 所示.可以看出在6 個實際風場中,風場1、2、6 與標稱參數差距較大,風場4、5 與標稱參數差距較小.因此可以認為,相比于標稱風場而言,風場1、2、6 是較差的風場,風場4、5 是較好的風場.

圖2 y 軸的六個實際風場和一個標稱風場Fig.2 Six actual wind fields and one nominal wind field on the y-axis

圖3 實際風場參數與標稱風場參數的距離Fig.3 Distances between parameters of actual wind fields and parameters of nominal wind field
本節的仿真參數如表1 所示,圖4 則給出在表1的仿真參數條件下隨機生成的一組風場,可以看出本節的參數條件較好地覆蓋了6 個實際風場.

表1 仿真參數值Table 1 Simulation parameter values

圖4 在表1 的仿真參數條件下隨機生成的一組風場Fig.4 A set of wind fields randomly generated under the simulation parameter conditions of Table 1
此外,本節使用Matlab 軟件包GPOPS 來求解最優控制問題.GPOPS 是一個基于偽譜法和序列二次規劃方法來求解最優控制問題的算法,具有求解速度快、精度高等特點.
注6.在表1 中,單階段RO 算法的終端約束范圍大于多階段RO 算法,這是因為當單階段RO 算法和多階段RO 算法取相同的終端約束范圍參數時,單階段RO 算法將由于約束范圍太小而無法求解.
在本節,將通過對比單階段RO 算法和多階段RO 算法來顯示多階段RO 算法在最終落點狀態精度上的優越性.單階段RO 算法和多階段RO 算法在不同風場下最終位置和期望位置的誤差如圖5 所示,最終速度和期望速度的誤差如圖6 所示,它們在各個實際風場下的軌跡如圖7 所示,控制–時間關系則如圖8、圖9 所示.從圖5 和圖6 可以看出,多階段RO 算法無論是位置還是速度,最終狀態的精度均高于單階段RO 算法.并且隨著風場的變化,多階段RO 算法的最終位置誤差波動也明顯小于單價段RO 算法.這顯示多階段RO 算法相比于單階段RO 算法在最終落點狀態精度上的優越性.

圖5 單階段RO 算法和多階段RO 算法在不同風場下的最終位置誤差Fig.5 The terminal position errors of single-stage RO algorithm and multi-stage RO algorithm under different wind fields

圖6 單階段RO 算法和多階段RO 算法在不同風場下的最終速度誤差Fig.6 The terminal velocity errors of single-stage RO algorithm and multi-stage RO algorithm under different wind fields

圖7 單階段RO 算法和多階段RO 算法在各個實際風場下的軌跡Fig.7 The trajectories of single-stage RO algorithm and multi-stage RO algorithm under each actual wind field

圖8 單階段RO 算法的控制–時間圖Fig.8 Control-time diagram of single-stage RO algorithm

圖9 多階段RO 算法的控制–時間圖Fig.9 Control-time diagram of multi-stage RO algorithm
此外,單階段RO 算法和多階段RO 算法的運行時間如表2 所示.從表2 可以看出,由于多階段RO 算法需要求解更多優化問題,所以具有更高的計算代價,這也表明在設計算法時階段數不宜過多.

表2 單階段RO 算法和多階段RO 算法的運行時間Table 2 Running time of single-stage RO algorithm and multi-stage RO algorithm
在本節,將通過單階段RO 算法與多階段RO算法在不同風場下的終端約束滿足情況以及單階段RO 算法與NO 算法的對比,來顯示RO 算法的有效性.單階段RO 算法在各個風場下最終位置、最終速度和相應的約束范圍如圖10、圖11 所示,其中,虛線代表該方向對應的約束范圍,圓圈代表對應的位置或速度.多階段RO 算法在各個風場下最終位置、最終速度和相應的約束范圍如圖12、圖13所示,其中,虛線代表該方向對應的約束范圍,叉號代表對應的位置或速度.從圖10~圖13 可以看出,無論是單階段RO 算法還是多階段RO 算法,它們在不同的風場下均未違背對應的終端狀態約束.

圖10 單階段RO 算法在不同風場下的最終位置與約束范圍Fig.10 The terminal position and its constraint range of single-stage RO algorithm under different wind fields

圖11 單階段RO 算法在不同風場下的最終速度與約束范圍Fig.11 The terminal velocity and its constraint range of single-stage RO algorithm under different wind fields

圖12 多階段RO 算法在不同風場下的最終位置與約束范圍Fig.12 The terminal position and its constraint range of multi-stage RO algorithm under different wind fields

圖13 多階段RO 算法在不同風場下的最終速度與約束范圍Fig.13 The terminal velocity and its constraint range of multi-stage RO algorithm under different wind fields
圖14 則對比在不同風場下單階段RO 算法與NO 算法的最終位置誤差.從圖14 可以看出,除了標稱風場(風場0),單階段RO 算法在其他風場下的最終位置精度都高于NO 算法,并且在較差的風場條件(風場1、2、6)下,兩算法精度差距較大.這表明NO 算法在面對未知的風場時,受風場變化的影響較大,而RO 算法則對變化的風場具有更好的魯棒性.上述討論說明了RO 算法的有效性.

圖14 單階段RO 算法與NO 算法在不同風場下的最終位置誤差Fig.14 The terminal position errors of single-stage RO algorithm and NO algorithm under different wind fields
本文考慮處于不確定風場中可回收火箭動力下降的魯棒制導問題,主要工作包括:
1)建立一個包含不確定性的風場模型用來刻畫未知風場,在該模型的基礎上給出火箭動力下降的魯棒優化問題.進一步,使用一種對不等式約束采取一階近似并將一階項作為安全裕量加入約束中的魯棒優化方法,得到一個易于求解的單階段魯棒優化算法.
2)從理論上定量地給出安全裕量的上界,借助該上界提出一種多階段魯棒優化方法,該方法相比于單階段魯棒優化算法,可以求解具有更小的終端約束范圍的優化問題.
3)通過仿真對比分析單階段、多階段魯棒優化算法和標稱風場下的優化算法,結果表明提出的多階段魯棒優化方法在落點精度較高的同時,對不同的風場具有較強的魯棒性.