閆 瑞，泮斌峰

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

1 引言

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

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

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

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

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

2 滑翔段縱向運動方程

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

式中：

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

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

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

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

3.1 縱向軌跡跟蹤控制律設計

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

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

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

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

橫截條件為式（6）：

極值條件為式（7）：

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

定義式（9）

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

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

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

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

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

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

表1 FDGMRES算法Table 1 FDGMRES algorithm

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

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

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

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

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

定義航向角誤差為導彈當前航向角與導彈相對目標的視線角之差，如式（15）所示：

式中，導彈當前位置相對目標點 ( θT，φT) 的視線角 ψLOS為式（16）：

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

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

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

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

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

4 仿真驗證與分析

4.1 無擾動的軌跡跟蹤控制

以滑翔式再入飛行器X-33為例進行仿真驗證，該飛行器質量為35 828 kg，參考面積為149.388 1m2，參考軌跡來自基于凸優(yōu)化方法的三維再入軌跡規(guī)劃［16］?；瓒纬跏紶顟B(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)量進行無量綱化，相應的標準值如第2節(jié)所述。設置仿真所需參數(shù)：采樣周期Δl為50 m，預測區(qū)間長度 Δs為250Δl，h為1E-7，kmax為1000?？紤]以跟蹤高度為主，依據(jù)Bryson準則選取Q、R分別為diag（［3E10，0，0］）、10，選取制導指令更新周期為1 s，仿真步長為0.1 s。

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

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

圖5 傾側角變化曲線Fig.5 Bank angle history

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

4.2 有擾動的軌跡跟蹤控制

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

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

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

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

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

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

5 結論

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

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

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