基于C/GMRES模型預(yù)測控制的滑翔段軌跡跟蹤研究

閆 瑞，泮斌峰

（西北工業(yè)大學(xué)航天學(xué)院，西安710072）

1 引言

滑翔式高超聲速飛行器作為一種典型的再入飛行器，具有速度快、航程遠(yuǎn)、機(jī)動性強(qiáng)等特點(diǎn)，其再入制導(dǎo)技術(shù)是保證飛行器完成飛行任務(wù)的關(guān)鍵技術(shù)。而滑翔式飛行器的再入飛行過程具有高度和速度變化范圍大、總體約束苛刻、飛行環(huán)境變化劇烈且氣動參數(shù)不確定等特點(diǎn)［1-2］，增大了制導(dǎo)律設(shè)計的難度。

再入制導(dǎo)方法一般可以分為標(biāo)稱軌跡跟蹤制導(dǎo)和預(yù)測校正制導(dǎo)。標(biāo)稱軌跡制導(dǎo)法在制導(dǎo)系統(tǒng)工作前已規(guī)劃好了飛行參考軌跡，計算簡單，對彈載計算機(jī)的存儲容量和計算速度要求較低，更容易實(shí)現(xiàn)；預(yù)測校正制導(dǎo)按落點(diǎn)預(yù)測方法的不同可分為解析預(yù)測制導(dǎo)法和數(shù)值預(yù)測制導(dǎo)法。前者在線計算量小，但是落點(diǎn)的預(yù)測精度不高；后者精度較高，但是數(shù)值積分計算量大，對機(jī)載計算機(jī)要求高［3］。因此研究標(biāo)稱軌跡跟蹤制導(dǎo)具有較大的應(yīng)用價值。

標(biāo)稱軌跡跟蹤制導(dǎo)中的軌跡跟蹤方法有：反饋線性化、線性二次型調(diào)節(jié)器（Linear Quadratic Regulator，LQR）、模型預(yù)測控制［4］等。Bharadwaj等［5］采用近似反饋線性化方法，將三維軌跡跟蹤轉(zhuǎn)化為對橫程和縱程的跟蹤，反饋線性化方法在模型存在較大的不確定性，非線性項(xiàng)無法精確對消，控制性能將變差甚至失穩(wěn)；Dukeman［6］基于LQR方法設(shè)計軌跡跟蹤控制律，根據(jù)二次型性能指標(biāo)離線計算最佳反饋系數(shù)，雖然LQR控制律對初始條件的拉偏具有良好的魯棒性，但是控制效果和精度取決于參考軌跡的特征，當(dāng)軌跡陡峭時，精確度將顯著降低；Lu［7］首先提出采用非線性模型預(yù)測控制完成軌跡跟蹤，與LQR方法相比，模型預(yù)測控制權(quán)重矩陣參數(shù)選擇較不敏感，該方法可以準(zhǔn)確地追蹤三維軌跡并可保證閉環(huán)穩(wěn)定性［8］。

非線性模型預(yù)測控制在每個采樣點(diǎn)通過求解有限時域最優(yōu)控制問題獲得最優(yōu)控制律。當(dāng)采用間接法在線求解最優(yōu)控制律時，需利用極小值原理將有限時域最優(yōu)控制問題轉(zhuǎn)化為兩點(diǎn)邊值問題并求解，但是Riccati矩陣微分方程的傳統(tǒng)求解方法計算量較大，耗時較長［9］。

針對求解控制律時計算量大的問題，本文采用基于C/GMRES算法［10-12］的非線性模型預(yù)測控制方法設(shè)計合理的跟蹤控制律對縱向參考軌跡進(jìn)行跟蹤控制：首先采用延拓法得到用于更新控制序列的微分方程組，使用廣義極小殘差法（Generalized Minimal Residual Methods，GMRES）高效解算近似線性方程，替代求解Riccati微分方程組，從而降低求解的復(fù)雜度，提高計算效率，保證了制導(dǎo)系統(tǒng)的實(shí)時性；然后基于航向角誤差走廊進(jìn)行橫向軌跡控制；最后進(jìn)行了數(shù)值仿真驗(yàn)證。

2 滑翔段縱向運(yùn)動方程

設(shè)再入飛行器滑翔段縱向運(yùn)動方程的狀態(tài)變量X=[ r Vγ]T，控制量u=σ，以剩余航程為自變量的高超聲速滑翔式再入飛行器縱向再入平面無量綱運(yùn)動方程［13］（不考慮地球自轉(zhuǎn)）可表示為式（1）：

式中：

其中，r為無量綱的飛行器到地球中心的距離，無量綱化參數(shù)為地球平均半徑R0=6 378 137m；V為無量綱的飛行器相對地球的速度值，無量綱化參數(shù)為為海平面重力加速度值；γ為當(dāng)?shù)厮俣葍A角；n為有量綱的飛行器質(zhì)量；σ為傾側(cè)角，L和D分別為無量綱的升力和阻力，表示為式（2）：

其中，ρ為有量綱的空氣密度，Sref為有量綱的飛行器參考面積，CL和CD分別為升力系數(shù)和阻力系數(shù)。

3 滑翔段參考軌跡跟蹤控制

滑翔段軌跡跟蹤制導(dǎo)分為縱向制導(dǎo)和橫向制導(dǎo)兩部分，本節(jié)設(shè)計控制律來跟蹤滑翔段縱向參考軌跡，由于縱向運(yùn)動只與傾側(cè)角σ的大小有關(guān)，因此將控制變量改寫為并對其進(jìn)行跟蹤控制；橫向運(yùn)動方面，采用傾側(cè)反轉(zhuǎn)邏輯確定傾側(cè)角的符號，將航向角誤差控制在合理范圍內(nèi)，實(shí)現(xiàn)橫向軌跡控制。

3.1 縱向軌跡跟蹤控制律設(shè)計

采用基于C/GMRES算法的非線性模型預(yù)測控制方法來跟蹤縱向參考軌跡。非線性模型預(yù)測控制問題本質(zhì)上是一類有限時域最優(yōu)控制問題，最小化性能指標(biāo)如式（3）所示：

其中Δs為預(yù)測時域的長度，Xref為參考軌跡狀態(tài)量，uref為參考軌跡傾側(cè)角絕對值，Q和R為權(quán)重系數(shù)。在每一剩余航程點(diǎn)s處，最小化預(yù)測時域τ∈[ s，s+Δs]內(nèi)的性能指標(biāo)J，同時滿足式（1）所示等式約束，以得到最優(yōu)控制輸入u*(τ)，只將初始值u*(s)作為剩余航程為s時的實(shí)際控制輸入u( s)。在下一航程點(diǎn)處，系統(tǒng)重新以新的航程點(diǎn)為初始狀態(tài)，求解最優(yōu)控制輸入，并繼續(xù)將其初始值施加給受控對象，如此循環(huán)，直至完成整個控制過程。

定義哈密爾頓函數(shù)如式（4）所示：

其中，λ為協(xié)態(tài)向量。根據(jù)極小值原理，上述以剩余航程s處的實(shí)際狀態(tài)量X( s) 為初始狀態(tài)的最優(yōu)控制問題可以轉(zhuǎn)化為兩點(diǎn)邊值問題，其中正則方程表示為式（5）：

橫截條件為式（6）：

極值條件為式（7）：

令采樣周期為Δl，將上述連續(xù)兩點(diǎn)邊值問題離散化，得到離散兩點(diǎn)邊值問題，如式（8）所示：

定義式（9）

在每一剩余航程s處求解方程組（10）得到最優(yōu)控制序列U( s)，實(shí)際控制輸入u( s) 取最優(yōu)控制序列的第一個元素( s)。

采用延拓法求解方程組（10），將其重寫為式（11）：

其中As為使F=0的穩(wěn)定矩陣。若FU非奇異，則可得到U( s) 的微分方程，如式（12）所示：

U′(s) 為控制量相對剩余航程的微分。由于式（12）中仍包含計算復(fù)雜的雅可比矩陣FX，F(xiàn)U和與相關(guān)的線性方程，為了減小求解雅可比矩陣和線性方程的計算量，采用雅可比矩陣和向量乘積的前向差分來近似F′，將式（11）近似為線性方程組，如式（13）所示：

其中，h為正實(shí)數(shù)。式（13）可以等價為式（14）：

前向差分近似和GMRES算法組合應(yīng)用于求解關(guān)于U′的近似線性方程（14），稱為FDGMRES算法［14］，表示為FDGMRES（U，X，X′，′，h，kmax），其中表示一個特定的初始猜 測，kmax表示最大迭代次數(shù)。表1中給出了算法的詳細(xì)步驟

表1 FDGMRES算法Table 1 FDGMRES algorithm

通過FDGMRES算法計算得到U′，實(shí)時積分U′得到U( s)。用下綴n表示在剩余航程s0+nΔl處的相應(yīng)變量，s0表示剩余總航程。將FDGMRES算法與上文中基于極小值原理的延拓法相結(jié)合，得到C/GMRES算法，其詳細(xì)過程如表2所述。

應(yīng)當(dāng)注意的是，迭代法僅用于求解關(guān)于U′的線性方程式（14），并且通過積分U′，得到非線性方程組F( U，X) =0的解U，而無需任何線性搜索或牛頓迭代求解。C/GMRES算法在每個采樣點(diǎn)僅求解一次線性方程式（14），因此與牛頓法等迭代方法多次求解線性方程以確定搜索方向相比，C/GMRES算法所需的計算量要少得多。此外，C/GMRES算法不涉及線性搜索，這也與標(biāo)準(zhǔn)優(yōu)化方法有顯著差異。

表2 C/GMRES算法Table 2 C/GMRES algorithm

3.2 基于航向角誤差走廊的橫向軌跡控制

飛行器飛行過程中航向誤差角Δψ應(yīng)一直在一個小范圍內(nèi)變化，使飛行器始終保持向目標(biāo)區(qū)域飛行。本文根據(jù)呂瑞［15］研究所采用的傾側(cè)反轉(zhuǎn)策略和航向角誤差走廊設(shè)計思路，離線規(guī)劃航向角誤差走廊以控制傾側(cè)角反轉(zhuǎn)，實(shí)現(xiàn)對飛行器飛行方向的控制。

定義航向角誤差為導(dǎo)彈當(dāng)前航向角與導(dǎo)彈相對目標(biāo)的視線角之差，如式（15）所示：

式中，導(dǎo)彈當(dāng)前位置相對目標(biāo)點(diǎn) ( θT，φT) 的視線角 ψLOS為式（16）：

采用漏斗形航向角誤差走廊的設(shè)計思想［15］，在再入的初始階段將誤差控制在較小范圍內(nèi)，來減小制導(dǎo)系統(tǒng)的工作負(fù)荷；中間段適當(dāng)放大誤差走廊，將其設(shè)計為速度的分段線性函數(shù)；在終端將誤差縮小到終端約束要求范圍內(nèi)。這種設(shè)計思路既可保證終端航向角的制導(dǎo)精度，也可使傾側(cè)角的反轉(zhuǎn)不過于頻繁。

對于特定飛行器，可通過多次仿真實(shí)驗(yàn)獲得適當(dāng)?shù)暮较蚪钦`差走廊邊界值。最終的航向角誤差走廊設(shè)計為圖1所示分段函數(shù)。

基于航向角誤差走廊的反轉(zhuǎn)邏輯的數(shù)學(xué)表達(dá)式為式（17）：

式中，sign( σi-1)為制導(dǎo)系統(tǒng)前一工作周期傾側(cè)角的符號。ΔψdH為航向角誤差走廊邊界的絕對值。該反轉(zhuǎn)策略表示當(dāng)導(dǎo)彈的航向角誤差的絕對值大于ΔψdH時，改變傾側(cè)角σ的符號；在邊界區(qū)域內(nèi)時，σ保持原來的符號，不進(jìn)行反轉(zhuǎn)；在滑翔段結(jié)束時，航向角誤差將符合終端約束。

圖1 航向角誤差走廊Fig.1 Heading angle error corridor

4 仿真驗(yàn)證與分析

4.1 無擾動的軌跡跟蹤控制

以滑翔式再入飛行器X-33為例進(jìn)行仿真驗(yàn)證，該飛行器質(zhì)量為35 828 kg，參考面積為149.388 1m2，參考軌跡來自基于凸優(yōu)化方法的三維再入軌跡規(guī)劃［16］。滑翔段初始狀態(tài)為：軌道地心距為6 452 735.38 m，速度為7528.77 m/s，速度傾角為-1.328 4°，經(jīng)度為237.58°，緯度為-9.18°，航向角為2°；終端條件：軌道地心距為6 408 761.57 m，速度為1133.37 m/s，經(jīng)度為238.25°，緯度為10.01°，航向角誤差的絕對值小于等于6°。剩余總航程為2138.589 km。對狀態(tài)量進(jìn)行無量綱化，相應(yīng)的標(biāo)準(zhǔn)值如第2節(jié)所述。設(shè)置仿真所需參數(shù)：采樣周期Δl為50 m，預(yù)測區(qū)間長度 Δs為250Δl，h為1E-7，kmax為1000?？紤]以跟蹤高度為主，依據(jù)Bryson準(zhǔn)則選取Q、R分別為diag（［3E10，0，0］）、10，選取制導(dǎo)指令更新周期為1 s，仿真步長為0.1 s。

圖2～圖5分別給出了滑翔段縱向軌跡跟蹤圖像、高度誤差、傾側(cè)角以及航向角誤差隨時間變化過程。

圖2 縱向軌跡跟蹤圖Fig.2 Tracking effect of longitudinal reference trajectory

圖3 高度誤差變化曲線Fig.3 A ltitude error history

圖4 速度傾角變化曲線Fig.4 Velocity path angle history

由圖2可知，基于C/GMRES算法的非線性模型預(yù)測控制律跟蹤效果較好；圖3中高度誤差隨時間在約（-80 m，60 m）范圍內(nèi)變化，并且誤差逐漸趨近于0，高度跟蹤效果較好；圖4給出了速度傾角隨時間變化曲線，速度傾角在約（-8°，1°）范圍內(nèi)變化，再入飛行器滑翔段傾側(cè)角反轉(zhuǎn)次數(shù)為3～6次比較合理［1，4，17］；圖5中傾側(cè)角反轉(zhuǎn)5次，反轉(zhuǎn)次數(shù)較少，保證了制導(dǎo)系統(tǒng)的穩(wěn)定性；圖6中航向角誤差在誤差走廊范圍內(nèi)合理變化，滑翔段結(jié)束時航向角誤差Δψf為-5.75°，符合終端約束要求，橫向制導(dǎo)精度較高。

圖5 傾側(cè)角變化曲線Fig.5 Bank angle history

圖6 航向角誤差變化曲線Fig.6 Heading error history

4.2 有擾動的軌跡跟蹤控制

進(jìn)一步進(jìn)行含有初始狀態(tài)擾動和動態(tài)干擾的蒙特卡洛仿真，檢驗(yàn)算法魯棒性?；枋斤w行器受到的控制力由氣動力提供，所以氣動力的擾動對再入制導(dǎo)的精度有重要的影響。氣動力擾動主要可以分為大氣密度擾動和氣動參數(shù)擾動，且大氣密度擾動可等效為氣動參數(shù)發(fā)生擾動［3］，因此改變動力學(xué)方程中的氣動系數(shù)進(jìn)行仿真驗(yàn)證，采用的偏差量如表3所示。在含有初始狀態(tài)擾動和動態(tài)干擾的條件下進(jìn)行1000次蒙特卡洛打靶仿真，選取其中100次仿真結(jié)果得到高度和航向角偏差隨時間變化圖像如圖7、8所示，1000次打靶的終端高度及航向角誤差如圖9、10所示。

表3 打靶加入的偏差Table 3 Perturbations for shoot

圖7 高度誤差變化曲線Fig.7 Altitude error histories

圖8 航向角誤差變化曲線Fig.8 Heading error histories

圖9 終端高度誤差Fig.9 Term inal altitude error

圖10 終端航向角誤差Fig.10 Term inal heading error

預(yù)先設(shè)定打靶落點(diǎn)高度誤差合理范圍為［-200 m，200 m］。圖7中高度偏差隨著時間逐漸減小，最后穩(wěn)定在約（-200 m，150 m）范圍內(nèi)，符合終端誤差要求；圖8中航向角誤差均在航向角誤差走廊范圍內(nèi)變化，滿足制導(dǎo)要求；圖9中打靶落點(diǎn)在（-200 m，150 m）范圍內(nèi)變化，分布范圍較小且均勻，滿足終端誤差要求；圖10中終端航向角誤差均在（-6°，-6°）范圍內(nèi)變化，符合終端航向角誤差要求。以上仿真分析表明所設(shè)計的制導(dǎo)律對初始狀態(tài)偏差、動態(tài)干擾魯棒性較好，滿足制導(dǎo)要求。

為了驗(yàn)證該控制律的快速性，本文在2 GB內(nèi)存和CPU主頻為2.8 GHz的計算機(jī)上以MATLAB R2018a為仿真平臺，對制導(dǎo)律的計算時間進(jìn)行統(tǒng)計。以1000次蒙特卡洛仿真為例，分別統(tǒng)計每次打靶中的最大制導(dǎo)計算時間和平均制導(dǎo)計算時間，如圖11、圖12所示。

圖11 最大制導(dǎo)計算時間直方圖Fig.11 Histogram of maximum guidance calculation time

圖12 平均制導(dǎo)計算時間直方圖Fig.12 Histogram of average guidance calculation tim e

圖11中最大制導(dǎo)計算時間分布在（0.03 s，0.28 s）范圍內(nèi)，并且在（0.05 s，0.06 s）范圍內(nèi)分布最多，遠(yuǎn)小于制導(dǎo)周期1 s；圖12中平均制導(dǎo)計算時間分布在（0.022 s，0.051 s）范圍內(nèi)，遠(yuǎn)小于制導(dǎo)周期1 s，且在（0.024 s，0.025 s）范圍內(nèi)分布最多。另外，采用信賴域算法求解上述非線性方程組（10）得到最優(yōu)控制律的平均制導(dǎo)計算時間約為3.892 7 s。采用C/GMRES算法求解制導(dǎo)律的平均制導(dǎo)計算時間約為采用信賴域算法求解計算時間的1%左右。以上統(tǒng)計比較充分地證明了基于C/GMRES算法的非線性模型預(yù)測控制律的快速性和實(shí)時性。

5 結(jié)論

本文采用基于C/GMRES算法的非線性模型預(yù)測控制律對再入飛行器滑翔段參考軌跡進(jìn)行跟蹤控制。

1）通過仿真分析驗(yàn)證了該控制方法對初始狀態(tài)擾動和動態(tài)干擾等不確定因素具有較強(qiáng)的魯棒性。

2）1000 次蒙特卡洛仿真的最大制導(dǎo)計算時間在（0.05 s，0.06 s）范圍內(nèi)分布最多，平均制導(dǎo)計算時間都分布在（0.022 s，0.051 s）范圍內(nèi)，遠(yuǎn)小于制導(dǎo)周期1 s，表明求解該控制律所需計算量小，計算效率高且易于實(shí)現(xiàn)，能夠滿足制導(dǎo)實(shí)時性要求。