一種近似解析的再入滑翔飛行器可達域快速生成方法

李兆亭,周 祥,張洪波,湯國建

國防科技大學空天科學學院,長沙 410073

0 引言

再入滑翔飛行器的可達域是指在滿足過程約束條件下,飛行器從初始狀態出發,能夠到達著陸或者交班狀態下所能允許的區域范圍。它綜合反映了飛行器的縱向機動能力與橫向機動能力,體現了飛行器的整體機動性能?？蛇_域的求解可為飛行器軌跡規劃與制導、應急著陸點選擇、制導律性能評估提供依據,具有重要意義。

可達域求解方法依賴于軌跡規劃方法的發展。在基于經典的剖面軌跡規劃方法研究條件下,發展得到了若干種可達域求解方法[1-3]。如He等[4]推導得到一個具有不確定性效應的再入走廊,其中攻角剖面可調節,縱向阻力剖面設計為上、下擬合安全邊界的插值結果,橫向升力-阻力剖面則利用準平衡滑翔條件得到,最后通過比例積分微分(PID)跟蹤器來跟蹤剖面,得到可達域。同時,隨著優化算法的成熟與發展,具有剖面形式的軌跡規劃問題亦可通過相關的優化算法進行規劃求解,如梁巨平等[5]選取傾側角為時間的分段常值函數,采用遺傳算法來求解一系列的軌跡優化問題。王濤等[6]提出一種基于Gauss偽譜法的再入可達域計算方法,采用固定的攻角剖面,僅對傾側角進行單變量尋優。

此外,還有一些研究者打破了剖面假設、擬平衡滑翔假設等藩籬,設計了其他類型的可達域求解方法,如基于虛擬目標點[7]、基于直接法[8-11]等。相比基于剖面的求解方法,這些新方法不需要設計相關剖面,其可達域范圍更大且更接近實際值。如藺君等[12]對攻角、傾側角進行參數化,利用帶約束的差分進化算法求解滿足再入過程約束和終端約束的再入軌跡。章吉力等[13]對考慮禁飛區條件下的空天飛機再入可達區域問題進行研究,并基于極限繞飛軌跡提出一種不限禁飛區位置的可達區域求解方法。吳楠等[14]基于平衡滑翔假設和最優化飛行準則,通過數值積分獲得最大縱程和橫程彈道,基于橢圓近似法利用3個末端點構建可達域橢圓邊界。

在基于直接法的可達域計算方法中,偽譜法由于具有良好的適用性和強大的求解能力成為其中一類重要的方法。如岳彩紅等[15]基于改進高斯偽譜法對高超聲速變形飛行器再入可達域進行了求解。Wang等[16]在準平衡滑行條件下,使用高斯偽譜法計算得到了返回航天器的可達域。此類可達域計算方法通常將可達域計算歸結為兩類問題,即最大/最小縱程軌跡規劃問題和給定縱程下的最大/最小橫程軌跡規劃問題。

本文借鑒了上述基于直接法的可達域計算思路,針對描述可達域的極值軌跡規劃問題,利用極大值原理推導出攻角應為當前最大升阻比攻角的相關推論,并仿真驗證了該推論的有效性。其次設計了給定縱程下不同橫程極值軌跡的傾側角變化規律,并通過梯度下降的優化方法計算得到相關參數。最終,提出了一種近似解析的可達域求解方法,并通過仿真驗證了所提方法的有效性。

1 再入動力學模型及其約束

在進行可達域的求解時,將飛行器的運動轉換到換極坐標系[17]中描述。在換極坐標系中,再入飛行器起點的經緯度均為0,終點的經度和緯度分別描述了再入縱程和再入橫程。利用這些特性,可方便地標定經緯度等參數的變化范圍,簡化彈道規劃算法。

在換極坐標系下,可以推導得到飛行器再入運動方程,如下:

(1)

(2)

其中:ρ為大氣密度,M為飛行器質量,S為參考面積,CL和CD分別為升力系數和阻力系數,均與攻角α有關。

同時,再入飛行中還需考慮多種過程約束,如熱流密度約束、過載約束、動壓約束等,如式(3)所示。

(3)

本文中,仿真所用的目標飛行器為CAV-H,其相關參數和氣動數據可以參考文獻[18]。

2 極值軌跡規劃問題最優特性分析

將可達域求解的兩類軌跡統稱為極值軌跡,意在說明該軌跡在縱、橫程方面具有的最優性。同時將兩類軌跡規劃問題統稱為極值軌跡規劃問題。

為降低推導的復雜度,本節中軌跡優化問題暫不考慮熱流、過載、動壓等約束。

2.1 最大縱程軌跡規劃問題

對于該問題而言,其升力幾乎全部用于縱向的調整,無側向機動,因此不妨令傾側角始終為0。同時,假設θ為小角度,忽略地球自轉項。在這些假設條件下,最大縱程軌跡優化問題就等價于最大路程軌跡優化問題。

記飛行器在橫向平面內的路程為S,則有:

(4)

與公式(1)中的有關式子聯立,假設θ為小量,并忽略二階小量,可得:

(5)

可將上述最大縱程軌跡優化問題近似轉化為問題P0,如下所示:

(6)

利用極大值原理求解上述問題P0的最優性條件。其哈密頓函數為:

(7)

其極值條件為:

(8)

關于控制量的部分取為極小,有:

(9)

考慮到控制量u1、u2為飛行器的升阻力加速度,應當滿足一定的關系式,有:

u1=ku2

(10)

式中:升阻比k在不同速度和攻角下為不同的值。此時控制量只有攻角。

將θ=x4/x3代入,則式(9)轉化為:

u1(λ4-λ3θ)-u2(λ4θ+λ3)

(11)

考察CAV-H飛行器的氣動情況,如圖1所示。飛行器的升阻比隨攻角具有相同的變化規律。其升阻比總是從A點的最小攻角開始,逐漸增大至B點的當前最大升阻比攻角,而后逐漸減小,直到端點C。

圖1 不同速度、高度下的升阻比隨攻角變化

(12)

因此不妨對式(11)分情況討論:當攻角位于AB段時,因參數(λ4-λ3θ)和(λ4θ+λ3)值大小的不同,其極小值一般取為A點或B點。當攻角位于BC段時,其極小值一般取為C點或B點。因此,最大縱程軌跡的攻角應為最大攻角、最小攻角和當前最大升阻比攻角三者之一。

2.2 給定縱程下的最大/最小橫程軌跡規劃問題

引入θ小角度假設,忽略地球自轉項,則給定縱程下的最大橫程軌跡規劃問題P1可以描述為式(12),式中λe為給定的縱程。而對于最小橫程問題而言,其性能指標取為J=φ(tf)+(λ(tf)-λe)2,其他完全一致。

求解問題P1的最優性條件,其哈密頓函數為:

(13)

其極值條件為:

(14)

提取關于控制量的部分,并記有,

(15)

則應當取極小值的部分為:

(16)

當ε+σ=-π/2+2kπ,k=0,1,…,N時,式(16)進一步轉化為:

(17)

2.3 結論驗證

綜上所述,對于這兩類極值軌跡規劃問題而言,其攻角應為最大攻角、最小攻角和當前最大升阻比攻角三者之一。

為檢驗上述結論的有效性,采用偽譜法對飛行器的可達域進行求解。采用兩個算例進行換極坐標系下可達域求解,算例1的初始高度為60 km,算例2的初始高度為65 km,兩個算例的初始速度均為6.5 km/s,初始速度傾角均為0°,終點速度要求大于1 km/s,終點速度傾角的大小不做限制要求,終點高度為30 km。

不同算例下各條極值軌跡的升阻比如圖2所示。在圖2中,短劃線代表當前最大升阻比,點虛線代表最小攻角的升阻比,下方的點劃線代表最大攻角的升阻比,上方的實線代表偽譜法規劃結果的升阻比。

圖2 不同算例下各極值軌跡的升阻比變化情況

可以看出,不同算例下飛行器的各極值軌跡幾乎總是沿著當前最大升阻比攻角前進。除了開始時刻和終端時刻略有差別。這種差別可能由于優化方法中一些終端狀態約束設置的不同引起。

從物理意義上,阻力代表了機械能的損耗速率。升力則有助于飛行器的上升,可增大速度傾角,使高度減小速率變緩甚至高度增加,從而增加縱程。從動力學角度來說,高度越高,阻力越小,也會使得機械能損失速率減小。那么,升阻比就代表了損失相同機械能的條件下,飛行器獲得的縱程增益。升阻比小意味著其縱程增益小,當升阻比為0時,其縱程增益為0,該運動退化為一個具有阻力加速度的拋物運動。顯然,對于構成可達域邊界的每一條軌跡,選擇一個最大升阻比攻角,無疑是一個最優的選擇。

因此,在計算可達域的極值軌跡時,不妨直接令攻角為當前最大升阻比攻角,積分到期望狀態條件,可以得到一條近似的極值軌跡。

3 近似解析的再入可達域生成方法

對于最大縱程軌跡規劃問題,可以以當前最大升阻比攻角和0傾側角積分得到其軌跡。但對于給定縱程下的最大/最小橫程軌跡規劃問題而言,還需得到傾側角的變化規律。

3.1 傾側角近似與參數計算

傾側角難以直接計算得到,在2.2節的分析中可知,其與拉格朗日乘子λ2和λ3直接相關,而拉格朗日乘子的估計一直是優化求解的難題。這里,考慮到文獻[19]中可達域的求解結果,其傾側角的變化非常有規律,因此不妨直接設計傾側角的變化規律,以得到近似的可達域。

設計不同極值軌跡的傾側角隨時間的曲線為不同斜率、不同初值的直線。直線的起點為不同的初始傾側角,直線的斜率為常值,直線的終點為0。直線的起點,即初始傾側角一般取為90°i/N,i=1,…,N,N為給定縱程下最大橫程軌跡的數目。第i條軌跡的斜率ki按照如下給定:

ki=(i-1)(kmax-kmin)/(N-1)+kmin

(18)

其中:kmax為最小縱程軌跡傾側角變化的斜率,kmin為次最大縱程軌跡傾側角變化的斜率。

對應的傾側角變化規律為:

σ=90°i/N+kit

(19)

最小縱程軌跡對應的傾側角變化規律為:

(20)

最大縱程軌跡對應的傾側角變化規律為:

σ=0, 0≤t≤tmax

(21)

其中:tmax可通過令傾側角為0并以當前最大升阻比攻角進行動力學積分得到。tmin可以通過求解傾側角符合上述變化規律的最小縱程問題得到。tmin的詳細求解流程如下所示:

1) 給定一個tmin的估計值t0;

2) 令tmin=t0和tmin=t0+10,并積分計算2種情況下的縱程,計算縱程對于tmin的梯度g,尋優方向為負梯度方向-sgn(g);

3) 回溯直線搜索方法確定尋優步長Δt,計算得到參數tmin=tmin-sgn(g)Δt;

4) 如果尋優步長小于某一閾值ε,即Δt≤ε,則停止,并輸出tmin;否則令t0=tmin,進入步驟2)。

如上所述,在求解得到參數tmax和tmin后即可確定kmax和kmin,則通過式(18)、(19)可以得到任意一條最大橫程軌跡的傾側角變化率。其給定縱程下最小橫程軌跡的傾側角變化規律與最大橫程一致,但符號是相反的。

3.2 近似解析的可達域生成方法

最終,形成了一種可達域的近似解析求解方法。在此方法中,最大縱程軌跡通過令傾側角為0,以當前最大升阻比攻角進行動力學積分得到。給定縱程下的最大最小橫程軌跡,通過以設計的傾側角變化規律和當前最大升阻比攻角進行動力學積分得到。其流程如圖3所示。

圖3 一種近似解析的可達域計算流程

如圖3所示,首先以0傾側角和當前最大升阻比攻角積分得到最大縱程軌跡,并得到最大飛行時間tmax。然后以式(20)所示的傾側角變化規律和當前最大升阻比攻角進行積分,并利用3.1節所述的優化方法計算得到最小飛行時間tmin。最后以式(18)、(19)所述的傾側角變化規律和當前最大升阻比攻角進行積分,得到一側的極值軌跡。改變傾側角的符號,積分得到另一側的極值軌跡,即可形成再入飛行器的可達域。

4 仿真校驗

為驗證所提方法的有效性,本節選取3個不同的飛行器起點狀態,分別利用所提方法和偽譜法計算得到各自的可達域。將得到的兩組結果作為對比,驗證本文所提方法的有效性和合理性。

可達域計算的參數設置如表1所示。在表1中,3個起點狀態被分別賦予了不同的速度與高度值,期望的終點速度大于1 km/s,期望的終點高度在10 km～30 km之間,同時不對終點的速度傾角做限制。得到的可達域近似結果如圖4所示。

表1 不同算例下的可達域求解參數設置

圖4 可達域計算結果對比

圖4給出了不同起點狀態下的可達域計算結果。其中實線代表偽譜計算得到的可達域邊界,點劃線代表由本方法計算得到的可達域邊界,虛線為由本方法計算得到的各條極值軌跡?？梢娤鄬τ趥巫V法計算的結果,本方法得到的可達域與之極為接近,僅存在較小的差異。

同時,考慮到如圖3所示的可達域計算流程,主要包含了參數tmin的優化求解、極值軌跡的積分等過程。且參數tmin的優化求解也僅僅借助于積分運算與梯度下降算法,因此對于整個求解過程而言,其不包含復雜的優化計算過程,具備工程試驗的可操作性。

5 結論

針對再入滑翔飛行器的可達域計算問題,提出一種近似解析的可達域計算方法。首先針對2類極值軌跡優化問題進行了相關推導,指出無論是最大縱程軌跡優化問題,還是給定縱程下的最大/最小橫程軌跡優化問題,其攻角應當是當前最大升阻比攻角。通過仿真驗證了上述推論的有效性。然后針對性地設計了各條極值軌跡的傾側角變化規律,給出了最終解析近似的可達域生成方法流程。最后通過與偽譜法計算得到的可達域的對比,驗證了本文所提方法的有效性。