基于鄰域最優控制律的再入軌跡在線修正

李世杰，李 炯，劉 滔，陳文鈺

(1.空軍工程大學 防空反導學院，西安 710051；2.中國空氣動力研究與發展中心 計算空氣動力研究所，四川 綿陽 621000；3.中國空氣動力研究與發展中心 空氣動力學國家重點實驗室，四川 綿陽 621000)

0 引 言

高超聲速飛行器再入大氣層，初始能量巨大，飛行速度較高，過程非常復雜，必須尋找一種快速可行軌跡優化的方法。再入軌跡優化問題實質上是一個帶有狀態約束和控制約束條件的非線性最優控制問題，當前求解最優控制問題的方法主要有間接法和直接法兩類[1-4]。其中，間接法利用Pontryagin極小值原理和一階最優必要條件，將最優問題轉換為兩點或多點邊值問題。這種方法控制參數較少，方程積分性好，不易發散，但對初值敏感。

間接法滿足最優性一階必要條件，且其解的精度較高。鄰域優化算法就是在其基礎上發展起來的[5-6]，該方法相比于一階變分的最優控制方法，具有二階精確度，且最為重要的特點是，在一階最優彈道小鄰域范圍內，能夠一次計算滿足一族初始或終端約束的最優彈道族，非常適用于形成優化彈道數據。

相比間接法，直接法對初值不敏感，其利用參數化方法將連續空間的最優控制問題轉化為NLP問題，然后再利用NLP問題求解工具[7]進行求解。直接法在解決實際復雜問題上具有優勢[8-9]，其典型代表就是偽譜法。近年來，由于偽譜法的快速求解效率，已逐漸成為求解軌跡優化問題方法的研究熱點[10-15]。

文獻[12]基于Radau多段偽譜法對最優機動突防彈道進行了求解，將突防脫靶量提高到百米量級，充分發揮了高超聲速飛行器的縱向機動性能，同時也證明了偽譜法處理彈道優化問題的優越性。文獻[13]基于偽譜法研究了高超聲速助推-滑翔飛行器的最大化射程以及最大化橫程的彈道優化問題。文獻[14]研究了基于間接Legendre偽譜法的最優制導律，通過Legendre偽譜法離線獲得參考軌跡，通過LQR方法在LG配點上求得控制量的修正值。文獻[15]在文獻[14]的基礎上，研究基于間接Radau偽譜法的軌跡跟蹤控制律，兩者的不同之處在于配點的選擇不同。文獻[14-15]中制導律的計算不需迭代及求解Riccati方程等復雜數值計算，能快速獲得控制量的修正值，但未解決增益矩陣的選擇問題，對計算結果影響較大。本文在此基礎上，提出基于鄰域最優控制的再入軌跡制導律，首先基于分段Radau偽譜法獲得最優參考軌跡，基于參考軌跡將最優控制問題轉化為狀態量偏差與協態量偏差的兩點邊值問題，然后求解LGR(Legendre-Gauss-Radau)配點上的最優反饋控制律，解決了LQR法的增益矩陣選擇問題，對初始狀態量的較大范圍偏差具有良好的魯棒性，并且能夠滿足實時性的需求。

1 再入軌跡制導問題

1.1 再入滑翔運動模型

忽略地球自轉的影響，高超聲速飛行器再入飛行的三自由度運動模型建立如下：

(1)

式中：飛行器的狀態量x=(r,φ,θ,v,γ,χ)，分別為地心矢徑、經度、緯度、速度、航跡傾角、航跡偏角；控制量u=[α,β]，分別為攻角和傾側角；m為飛行器質量；g=μ/r2為重力加速度，μ=398 603.2 km3/s2為地球引力常數)。

升力L和阻力D分別為

(2)

CL(α)=CL0+CLαα

CD(α)=CD0+CDαα+CDα2α2

(3)

式中：CL0=-0.207 0；CLα=1.676；CD0=0.078 54；CDα=-0.352 9；CDα2=2.040；α為攻角。

1.2 再入約束條件

再入約束包括過程約束、終端約束及控制量約束，其中過程約束主要有動壓約束、過載約束、狀態量約束及熱流密度約束等。

熱流密度約束：Q=QaQr≤Qmax。

式中：Qa=h0+h1α+h2α2+h3α3，h0=1.067，h1=-1.101，h2=0.698 8，h3=-0.190 3；Qr=CρNvM，N=0.5，M=3.07，C=9.289×10-9(kJ/m3.57)/kg。

其他變量約束-89°≤χ≤89°,-89°≤γ≤89°。

初始條件：r0=(Re+79.248) km,θ0=0°,φ0=0°,v0=6.492 km/s,γ0=-1.5°,χ0=0°。

終端約束條件包括位置、速度和角度等，即rf=(Re+24.384) km，vf=0.762 km/s，γf=-10°。其中：Re是平均地球半徑。

優化軌跡以最大化經度為目標函數，同時滿足上述一系列約束條件：

J=min{-θf}

(4)

1.3 最優參考軌跡的生成

最優軌跡生成即求解控制量u=[α,β]，使目標函數最小(或最大)，并且狀態滿足運動方程及各種約束條件。軌跡優化問題可當成是一般的最優控制問題。首先需將最優控制問題的時域[t0,tf]轉化為[-1,1]區間，對時間變量做變換，即

(5)

為不失一般性，考慮Bolza形式的一般最優控制問題，性能泛函指標為

(6)

式中：狀態變量x(t)∈Rn；初始時間和終端時間t0,tf(自由或固定)滿足動力學微分方程約束。

(7)

初始狀態x(t0)=x0，終端邊界條件約束為

N[x(tf),tf]=0

(8)

針對上述典型的兩點邊值最優控制問題，采用文獻[16]中的分段Radau偽譜法進行求解。該方法在求解中采用Legendre多項式進行狀態量、控制量和協態量估計與離散化，將微分約束轉化為代數約束，將最優控制問題轉化為非線性規劃問題，再采用非線性規劃求解器SNOPT[7]等進行求解。針對上述飛行器再入過程，采用文獻[16]中的分段Radau偽譜法計算最優參考軌跡。

2 再入軌跡在線修正問題

2.1 鄰域最優控制問題

由于隨機擾動及系統誤差，飛行器的實際飛行軌跡與最優參考軌跡并不一致。為解決此問題，一般是重新規劃最優軌跡，但這種方法計算量較大，花費時間較多，難以實現；另一種方法是采用鄰域最優控制理論，通過計算獲得控制量的修正值。鄰域最優控制問題描述如下：

(9)

設u*(t)為最優控制，x*(t)為由此產生的最優曲線，存在與之相對應的協態λ*(t)，根據極小值原理，滿足正則方程：

(10)

哈密爾頓函數為

H(x,u,λ,t)=L(x,u,t)+λTf(x,u,t)

(11)

最優控制條件為

(12)

若終端狀態受約束，則邊界約束條件為

(13)

若終端狀態自由時，則邊界約束條件為

(14)

不論終端狀態，終端邊界條件可統一為

ψ[x(tf),λ(tf)]=0

(15)

若終端時刻自由時，確定tf的終端橫截條件為

(16)

考慮飛行器實際飛行軌跡與最優參考軌跡的細小初始偏差δx(t0)，初始偏差會產生狀態量、協態量及控制量與最優參考軌跡的偏差，分別為δx，δλ和δu，對其進行局部線性化得

(17)

當Huu為非奇異矩陣時，可得

(18)

式(17)中，δx(t)和δλ(t)的狀態方程及初始條件δx(t0)=δx0與終端約束條件0=[ψxδx+ψλδλ]t=tf形成了δx(t)和δλ(t)的線性兩點邊值問題。

2.2 鄰域最優控制律設計

最優參考軌跡是在域[-1,1]中求解的，所以鄰域最優控制設計在此域中，狀態和協態偏差量微分及邊界條件轉換為

(19)

δx(τ0)=δx0

(20)

ψxδx(τf)+ψλδλ(τf)=0

(21)

偽譜法是將狀態變量與控制變量在LGR點上進行離散[16]，在區間Sk,(k=1，2，…，K)中，利用Lagrange 插值方法對狀態量偏差量δx(τ)和協態變量偏差δλ(τ)在LGR點進行近似：

(22)

(23)

利用LGR點處近似值的微分仍為代數形式，對狀態量偏差量和協態變量偏差微分得

(24)

(25)

(26)

(1+τi)[PNk(τi)-PNk-1(τi)]，PNk(τi)為Nk階Legendre多項式。從而有

(27)

式中：i=1，2，…，N

(28)

式中:

為計算簡便，令

(29)

(30)

(31)

(32)

式中:In與0n分別為n×n階單位矩陣和n×n階零矩陣，則有

(33)

(34)

邊界條件為

δx(t0)=δx0

(35)

ψxδxN+ψλδλN=0

(36)

為解方程式(33～34)，將其與邊界條件合并，則有

(37)

式中:P1與P2為n×n的N階矩陣，P1=[0n,…,0n,ψx]，P2=[0n,…,0n,ψλ]。令ZT=[XT,ΛT]，V=[V1,…,Ve]，可得

V1δX1+VeZe=0

(38)

Ze=-VeV1δX1=WδX1

(39)

式中：為MATLAB中的左除算子。令Z=[δX1Ze]T，則可得

(40)

式中：Wx和Wλ為[InW]T的分塊矩陣，且均為nN×n階。因而有

δXi=WxiδX1

δΛi=WλiδΛ1

(41)

式中：Wxi和Wλi分別為矩陣Wx和矩陣Wλ的部分。則有

(42)

從而可得

(43)

基于最優軌跡的在線修正求解過程為：首先，采用文獻[16]中的偽譜法計算最優參考軌跡X*(τ)和參考控制量U*(τ)，然后根據實際軌跡與最優參考軌跡的狀態偏差量，通過鄰域最優控制律計算獲得LGR點處的修正控制量δu(τi)，節點處修正后的控制量為u(τi)=u*(τi)+δu(τi)，最后通過插值可獲得整個軌跡修正后的控制量。

這種鄰域最優控制算法在執行過程中不需要迭代運算和積分運算，只是進行矩陣運算，減小了計算量，具有較強的實時性。

2.3 初始偏差下的鄰域最優控制律驗證

針對上述的高超聲速飛行器再入過程，首先生成最優參考軌跡，采用文獻[16]中的分段Radau偽譜法進行計算，在聯想CPU 3.4 GHz Intel Core i7計算機上運行，仿真結果如圖1所示。

圖1 基于分段Radau偽譜法的仿真結果Fig.1 Simulation results based on multi-stage Radau pseudospectral method

圖1為飛行器飛行狀態地心矢徑、緯度、經度、速度、航跡傾角、航跡方位角以及控制變量攻角與傾側角隨時間變化的曲線。整個飛行過程為1 656 s，最終經度為60.56°，仿真計算時間為6.87 s。在給定的初始條件下，整個飛行中狀態量及控制量滿足各自的上下邊界約束條件，飛行器最終狀態滿足給定的終端約束條件。

為了驗證鄰域最優控制律的有效性，令狀態量初始擾動為|Δr(0)|≤1 000 m；|Δφ(0)|≤0.5°;|Δθ(0)|≤0.5°；|Δv(0)|≤100 m/s;|Δγ(0)|≤0.5°；|Δχ(0)|≤0.5°。

為了便于比較，令case 1和case 2分別為狀態量初始擾動取最大正值和最大負值。close 1和close 2分別為擾動case 1和case 2時采用上述鄰域最優控制律求得的修正控制量。與其對應，open 1和open 2分別為未采用鄰域最優控制律求得的偏差量,仿真結果如圖2所示。

圖2 飛行器狀態與標準參考軌跡的比較Fig.2 Comparison of aircraft state and standard reference trajectory

close 1和close 2的偏差量最終都趨近于零。而open 1和open 2中的偏差量均有較大值，且有的狀態量的偏差絕對值逐漸增大。仿真結果表明，鄰域最優控制律對飛行器初始狀態量的較大范圍偏差具有良好的魯棒性。

考慮到整個再入過程中高超聲速飛行器的熱流密度、動壓、過載等過程約束，分別設定Qmax=79.5 kW/m2，qmax=0.013 4 MN/m2，nmax=2.5g，驗證如圖3所示。圖3中normal為初始狀態量無擾動時采用文獻[16]分段Radau偽譜法仿真的結果。close 1與close 2分別為初始擾動取case 1和case 2時采用鄰域最優控制律仿真的結果。與其對應，open 1與open 2分別為未采用鄰域最優控制律仿真的結果。

圖3 過程約束情況Fig.3 The case of the process constraint

從圖3可以看出，文獻[16]方法對再入最優參考軌跡的求解是有效的。close 1與close 2曲線表明，存在初始狀態量大范圍擾動時的實際飛行軌跡滿足動壓、過載、熱流等路徑約束條件。open 1與open 2曲線表明，當不采用本文鄰域最優控制律時，熱流約束超過了限制值。仿真結果表明,本文設計的鄰域最優控制律在最優軌跡跟蹤中的有效性。

3 結 論

本文在多段Radau偽譜法生成最優參考軌跡的基礎上，針對軌跡制導中存在初始擾動的問題，提出一種基于間接Radau偽譜法的鄰域最優控制律，得出以下結論：

(1) 多段Radau偽譜法能夠快速求解獲得一條精確的再入最優軌跡，并滿足路徑約束條件及終端約束條件。

(2) 本文提出了基于Radau偽譜法的鄰域最優控制律，對該制導律的性能進行了仿真驗證。首先過程采用文獻[16]中的分段Radau偽譜法獲得了最優參考軌跡，然后驗證了在狀態量具有較大初始擾動的情況下，本文基于Radau偽譜法的鄰域最優控制律可以在線修正最優參考軌跡，減少了初始擾動引起的偏差，證明了其具有良好的魯棒性。通過進行矩陣運算，能夠快速獲得控制量的修正值，證明了其具有良好的實時性。

(3) 當終端約束條件發生變化，重新規劃新的軌跡不能滿足實時性要求時，如何通過鄰域最優控制律使實際軌跡能夠適應終端約束條件的變化，是下一步的研究內容。