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

(1.空軍工程大學(xué) 防空反導(dǎo)學(xué)院，西安 710051；2.中國(guó)空氣動(dòng)力研究與發(fā)展中心 計(jì)算空氣動(dòng)力研究所，四川 綿陽(yáng) 621000；3.中國(guó)空氣動(dòng)力研究與發(fā)展中心 空氣動(dòng)力學(xué)國(guó)家重點(diǎn)實(shí)驗(yàn)室，四川 綿陽(yáng) 621000)

0 引 言

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

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

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

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

1 再入軌跡制導(dǎo)問(wèn)題

1.1 再入滑翔運(yùn)動(dòng)模型

忽略地球自轉(zhuǎn)的影響，高超聲速飛行器再入飛行的三自由度運(yùn)動(dòng)模型建立如下：

(1)

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

升力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 再入約束條件

再入約束包括過(guò)程約束、終端約束及控制量約束，其中過(guò)程約束主要有動(dòng)壓約束、過(guò)載約束、狀態(tài)量約束及熱流密度約束等。

熱流密度約束：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是平均地球半徑。

優(yōu)化軌跡以最大化經(jīng)度為目標(biāo)函數(shù)，同時(shí)滿足上述一系列約束條件：

J=min{-θf(wàn)}

(4)

1.3 最優(yōu)參考軌跡的生成

最優(yōu)軌跡生成即求解控制量u=[α,β]，使目標(biāo)函數(shù)最小(或最大)，并且狀態(tài)滿足運(yùn)動(dòng)方程及各種約束條件。軌跡優(yōu)化問(wèn)題可當(dāng)成是一般的最優(yōu)控制問(wèn)題。首先需將最優(yōu)控制問(wèn)題的時(shí)域[t0,tf]轉(zhuǎn)化為[-1,1]區(qū)間，對(duì)時(shí)間變量做變換，即

(5)

為不失一般性，考慮Bolza形式的一般最優(yōu)控制問(wèn)題，性能泛函指標(biāo)為

(6)

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

(7)

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

N[x(tf),tf]=0

(8)

針對(duì)上述典型的兩點(diǎn)邊值最優(yōu)控制問(wèn)題，采用文獻(xiàn)[16]中的分段Radau偽譜法進(jìn)行求解。該方法在求解中采用Legendre多項(xiàng)式進(jìn)行狀態(tài)量、控制量和協(xié)態(tài)量估計(jì)與離散化，將微分約束轉(zhuǎn)化為代數(shù)約束，將最優(yōu)控制問(wèn)題轉(zhuǎn)化為非線性規(guī)劃問(wèn)題，再采用非線性規(guī)劃求解器SNOPT[7]等進(jìn)行求解。針對(duì)上述飛行器再入過(guò)程，采用文獻(xiàn)[16]中的分段Radau偽譜法計(jì)算最優(yōu)參考軌跡。

2 再入軌跡在線修正問(wèn)題

2.1 鄰域最優(yōu)控制問(wèn)題

由于隨機(jī)擾動(dòng)及系統(tǒng)誤差，飛行器的實(shí)際飛行軌跡與最優(yōu)參考軌跡并不一致。為解決此問(wèn)題，一般是重新規(guī)劃最優(yōu)軌跡，但這種方法計(jì)算量較大，花費(fèi)時(shí)間較多，難以實(shí)現(xiàn)；另一種方法是采用鄰域最優(yōu)控制理論，通過(guò)計(jì)算獲得控制量的修正值。鄰域最優(yōu)控制問(wèn)題描述如下：

(9)

設(shè)u*(t)為最優(yōu)控制，x*(t)為由此產(chǎn)生的最優(yōu)曲線，存在與之相對(duì)應(yīng)的協(xié)態(tài)λ*(t)，根據(jù)極小值原理，滿足正則方程：

(10)

哈密爾頓函數(shù)為

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

(11)

最優(yōu)控制條件為

(12)

若終端狀態(tài)受約束，則邊界約束條件為

(13)

若終端狀態(tài)自由時(shí)，則邊界約束條件為

(14)

不論終端狀態(tài)，終端邊界條件可統(tǒng)一為

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

(15)

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

(16)

考慮飛行器實(shí)際飛行軌跡與最優(yōu)參考軌跡的細(xì)小初始偏差δx(t0)，初始偏差會(huì)產(chǎn)生狀態(tài)量、協(xié)態(tài)量及控制量與最優(yōu)參考軌跡的偏差，分別為δx，δλ和δu，對(duì)其進(jìn)行局部線性化得

(17)

當(dāng)Huu為非奇異矩陣時(shí)，可得

(18)

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

2.2 鄰域最優(yōu)控制律設(shè)計(jì)

最優(yōu)參考軌跡是在域[-1,1]中求解的，所以鄰域最優(yōu)控制設(shè)計(jì)在此域中，狀態(tài)和協(xié)態(tài)偏差量微分及邊界條件轉(zhuǎn)換為

(19)

δx(τ0)=δx0

(20)

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

(21)

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

(22)

(23)

利用LGR點(diǎn)處近似值的微分仍為代數(shù)形式，對(duì)狀態(tài)量偏差量和協(xié)態(tài)變量偏差微分得

(24)

(25)

(26)

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

(27)

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

(28)

式中:

為計(jì)算簡(jiǎn)便，令

(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)

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

這種鄰域最優(yōu)控制算法在執(zhí)行過(guò)程中不需要迭代運(yùn)算和積分運(yùn)算，只是進(jìn)行矩陣運(yùn)算，減小了計(jì)算量，具有較強(qiáng)的實(shí)時(shí)性。

2.3 初始偏差下的鄰域最優(yōu)控制律驗(yàn)證

針對(duì)上述的高超聲速飛行器再入過(guò)程，首先生成最優(yōu)參考軌跡，采用文獻(xiàn)[16]中的分段Radau偽譜法進(jìn)行計(jì)算，在聯(lián)想CPU 3.4 GHz Intel Core i7計(jì)算機(jī)上運(yùn)行，仿真結(jié)果如圖1所示。

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

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

為了驗(yàn)證鄰域最優(yōu)控制律的有效性，令狀態(tài)量初始擾動(dòng)為|Δ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分別為狀態(tài)量初始擾動(dòng)取最大正值和最大負(fù)值。close 1和close 2分別為擾動(dòng)case 1和case 2時(shí)采用上述鄰域最優(yōu)控制律求得的修正控制量。與其對(duì)應(yīng)，open 1和open 2分別為未采用鄰域最優(yōu)控制律求得的偏差量,仿真結(jié)果如圖2所示。

圖2 飛行器狀態(tài)與標(biāo)準(zhǔn)參考軌跡的比較Fig.2 Comparison of aircraft state and standard reference trajectory

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

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

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

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

3 結(jié) 論

本文在多段Radau偽譜法生成最優(yōu)參考軌跡的基礎(chǔ)上，針對(duì)軌跡制導(dǎo)中存在初始擾動(dòng)的問(wèn)題，提出一種基于間接Radau偽譜法的鄰域最優(yōu)控制律，得出以下結(jié)論：

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

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

(3) 當(dāng)終端約束條件發(fā)生變化，重新規(guī)劃新的軌跡不能滿足實(shí)時(shí)性要求時(shí)，如何通過(guò)鄰域最優(yōu)控制律使實(shí)際軌跡能夠適應(yīng)終端約束條件的變化，是下一步的研究?jī)?nèi)容。