基于間接Radau偽譜法的滑翔段軌跡跟蹤制導律

廖宇新，李惠峰，包為民，2

(1．北京航空航天大學宇航學院，北京100191;2．中國航天科技集團公司科技委，北京100048)

0 引言

在滑翔段飛行過程中，制導系統承擔著引導高超聲速飛行器在不違背各種約束的前提下準確到達目標這一至關重要的作用。因此，研究具有魯棒性和實時性的滑翔段制導律意義重大。

按照制導策略通常將滑翔段制導律分為預測校正制導［1－3］和標稱軌跡跟蹤制導［4－10］。預測校正制導是根據飛行器當前的真實狀態對軌跡進行快速預測，得到終端狀態量，然后根據所得結果與終端約束條件的偏差值對控制量進行校正。預測校正制導不需要存儲標稱軌跡，能適應更大范圍的偏差和擾動。但是，數值預測校正制導對彈載計算機的計算能力要求很高，難以滿足實時性的需求，且不能保證收斂的可行解的存在;而解析預測校正制導所用的模型都經過了簡化處理，解的預測精度和對飛行任務的適應能力均不理想。

標稱軌跡跟蹤制導主要包括標稱軌跡生成和軌跡在線跟蹤兩部分。為了提高標稱軌跡跟蹤制導的自主性、自適應性和魯棒性，對于現有方法的改進主要沿著兩條技術路線展開:一是研究軌跡在線生成算法，二是研究能夠在線實時解算且具有魯棒性和自適應性的標稱軌跡跟蹤方法［4］。文獻［5］以高度、速度、航跡傾角的變化曲線作為標稱軌跡，對降階的運動方程組進行線性化，將二維平面軌跡跟蹤問題轉化為線性時變(Linear Time-Varying，LTV)系統狀態調節器問題來求解，但基于系數凍結法的在線增益調參技術在處理LTV系統狀態調節器問題時存在費時、魯棒性差、穩定性無理論保證等不足。文獻［6］通過求解Riccati矩陣微分方程來求解航天飛機三維平面標稱軌跡跟蹤問題轉化得到的LTV系統狀態調節器問題，但在線求解Riccati矩陣微分方程也存在費時、數值不穩定等不足［8］。文獻［7］將三維平面標稱軌跡跟蹤問題轉化為LTV系統滾動時域控制問題，為避免求解Riccati矩陣微分方程的數值計算負擔，將LTV系統滾動時域控制問題離散化，近似轉化為二次規劃(Quadratic Programming，QP)問題來求解，并證明了閉環系統的穩定性。文獻［7］的方法具有很強的魯棒性，但當離散的階次較高時，推導和計算都比較復雜。文獻［11］利用間接Legendre偽譜法(Indirect Legendre Pseudospectral Method，ILPM)求解LTV系統狀態調節器問題轉化得到的線性兩點邊值問題(Two Point Boundary Value Problem，TPBVP)，不存在求解Riccati矩陣微分方程和在線增益調參的過程。文獻［4，8－10］結合相關策略進一步研究了基于ILPM的高超聲速飛行器滑翔段軌跡跟蹤制導律，驗證了制導律的魯棒性和實時性。

本文提出了間接Radau偽譜法(Indirect Radau Pseudospectral Method，IRPM)來求解LTV系統狀態調節器問題轉化得到的線性TPBVP。文獻［11］中，線性TPBVP經過ILPM離散轉化后，得到的是超定線性方程組，只存在最小二乘解。本文與文獻［11］不同的是，線性TPBVP經過IRPM離散轉化后，得到的是適定線性方程組，能求得唯一的精確解。基于IRPM求解得到的最優反饋控制律，提出了一種滑翔段全狀態標稱軌跡跟蹤制導律。該制導律具有良好的魯棒性和在線執行的能力。

1 標稱軌跡跟蹤制導問題的轉化

假設地球為自轉圓球，高超聲速飛行器滑翔段無動力飛行的三自由度運動方程組參考文獻［3］。飛行器的狀態量x=［r，θ，φ，V，γ，ψ］，分別為地心矢徑、經度、緯度、速度、航跡傾角、航跡偏角;控制量u=［α，σ］，分別為攻角和傾側角。

實際飛行軌跡與標稱軌跡之間存在偏差，其中狀態量偏差量為 δx(t)∈Rn，控制量偏差量為δu(t)∈Rm。基于標稱軌跡(xr，ur)對運動方程組進行線性化，得到LTV系統的狀態方程和初始條件分別為

進一步考慮如下的二次型性能指標

式中:Sf∈Rn×n和Q(t)∈Rn×n是對稱半正定矩陣，R(t)∈Rm×m是對稱正定矩陣。標稱軌跡跟蹤制導問題可以轉化為LTV系統狀態調節器問題:確定最優反饋控制律 δu*(δx*，t)，在滿足式(1)、(2)所示約束的前提下，使式(3)所示性能指標函數取最小值。

令哈密爾頓函數為

式中:λ(t)∈Rn為協態變量。由Pontryagin極大值原理可得線性TPBVP為

一階最優必要條件為

進一步可得最優反饋控制律為

由此可知，通過求解上述線性TPBVP，獲得協態變量 λ(t)，即可得到最優反饋控制律 δu*(δx*，t)。

2 基于間接Radau偽譜法的最優反饋控制律

為了避免后向掃描法中反向積分Riccati矩陣微分方程時的數值敏感性，以及轉換矩陣法中可能產生病態轉換矩陣導致數值計算不穩定的問題［8］，同時滿足實時計算最優反饋控制律的需求，本文采用間接Radau偽譜法來求解上述線性TPBVP，進而獲得最優反饋控制律。

由于Legendre-Gauss-Radau(LGR)點 τk∈［－1，1)，首先需要將時間區間t∈［t0，tf］通過式(9)轉換到［－1，1］

時間區間轉換后的線性TPBVP為

以N個LGR點τk和終端時刻點τf=τN+1=+1為節點，利用Lagrange插值方法對狀態量偏差量δx(τ)和協態變量λ(τ)進行近似

式中:τj∈［－1，1］為插值節點，?j(τ)為Lagrange插值基函數。由此可得，狀態量偏差量和協態變量關于變換后時間的導數在LGR點處的值分別為

將式(15)、(16)代入式(10)，線性TPBVP中的微分方程組轉化為如下的代數方程組

令δx=［(δx(τ1))T，(δx(τ2))T，…，(δx(τN+1))T］T，λ=［λT(τ1)，λT(τ2)，…，λT(τN+1)］T，Ze=［(δx(τ2))T，…，(δx(τN+1))T，λT(τ1)，…，λT(τN+1)］T，Z=［(δx)T，λT］T=，綜合式(11)、(17)、(18)并寫成如下式所示的分塊的形式

式中:E、F、G、H∈RNn×(N+1)n，Z∈R(2N+2)n×1，V∈R(2N+1)n×(2N+2)n，P1=［0n×n，…，0n×n，Sf］ ∈Rn×(N+1)n，P2=［0n×n，…，0n×n，－In×n］∈Rn×(N+1)n，δx0∈Rn×1，V0∈R(2N+1)n×n，Ze∈R(2N+1)n×1，Ve∈R(2N+1)n×(2N+1)n，

式(19)中有(2N+1)n個方程，對應于(2N+1)n個未知數，這表明VeZe=－V0δx0是適定線性方程組。由文獻［10］可知，矩陣V是列滿秩的，且可分解為矩陣V0和方陣Ve，所以Ve是列滿秩的方陣。因此，可求得VeZe=－V0δx0唯一的精確解為

由此可得

式中:W∈R(2N+2)n×n，Wx、Wλ∈R(N+1)n×n，(Wx)j、(Wλ)j∈Rn×n分別為Wx、Wλ中對應于 τj時刻的子矩陣。

將式(26)代入式(8)，當初始時刻的狀態量偏差值δx0已知時，最優反饋控制律在離散節點處的值可以表示為

非離散節點處的標稱控制量和最優反饋控制律求得的控制量偏差值則根據離散節點處的值通過Lagrange插值的方法獲得。

3 閉環制導律

本文研究的滑翔段軌跡跟蹤制導律的執行過程如圖1所示，對應的流程如下:

圖1 閉環制導律的執行過程示意圖Fig．1 Closed-loop guidance law execution diagram

(2)以當前時刻標稱控制量ur與當前制導周期內計算所得的當前時刻最優反饋控制律δu(τ0)之和ur+δu(τ0)作為實際飛行的控制量，轉入步驟(4);

(3)以當前時刻標稱控制量ur與由上一個制導周期的離散節點處的最優反饋控制律(τj)插值所得的當前時刻最優反饋控制律δu(τ)之和ur+δu(τ)作為實際飛行的控制量，轉入步驟(4);

上述滑翔段標稱軌跡跟蹤制導律不需要離線設計誤差反饋增益系數和在線增益調參，在制導指令解算過程中不涉及迭代、在線積分求解Riccati矩陣微分方程等復雜數值運算，計算量小，因此便于在線執行。

4 仿真校驗

4．1標稱軌跡

以CAV-H飛行器［12］為研究對象，質量m為907 kg，參考面積S為0．48387 m2;飛行任務的端點約束如表1所示，過程約束(包括熱流密度Q˙、總過載nz和動壓q)和控制量約束如表2所示，標稱飛行攻角剖面αref為

表1 飛行任務的端點約束Table 1 Endpoint constraints for flight mission

表2 過程約束和控制量約束Table 2 Path constraints and control constraints

軌跡優化時選擇的LGR點數目為40，優化所得軌跡的終端狀態量分別為hf=30000 m、θf=65°、φf=3°、Vf=2500 m/s、γf=－0．8536°、ψf=51．6679°，飛行時間為1481．80 s。狀態量、控制量和過程約束隨飛行時間的變化過程分別如圖2～圖4中的實線所示。可以看出優化所得的飛行軌跡滿足飛行任務的要求且未違背過程約束，這說明優化所得的軌跡是正確和可行的。下文的制導仿真都以本節優化生成的軌跡作為標稱軌跡，

4．2制導律魯棒性驗證

為了驗證本文的滑翔段制導律的魯棒性，在飛行器初始狀態量存在偏差和質量、氣動參數、大氣密

利用Radau偽譜法［13］離線優化獲得滑翔段的標稱軌跡。為滿足飛行軌跡平滑的需求，離線軌跡優化采用的性能指標為度等飛行環境參數有擾動情況下進行制導仿真。

仿真環境為雙核3．00 GHz主頻CPU、3 GB內存及Matlab2012a軟件。制導律解算采用的LGR點數目為10，制導周期和積分步長均取為1 s，權值矩陣Sf、Q、R的取值參考Bryson法則［14］，且在本文所有的仿真中，制導周期、積分步長和權值矩陣的取值都一致。制導律仿真過程中，傾側角σ是主要的軌跡控制量，攻角α只在標稱剖面附近進行小范圍調整，調整范圍如表2中Δα所示。

初始狀態量偏差值、飛行環境參數擾動量和閉環制導仿真所得的終端狀態量相對于標稱軌跡的偏差值如表3所示。其中，組合1和組合2分別對應于 δh0=3000 m、δθ0=0．3°、δφ0=0．3°、δV0=100 m/s、δγ0=0．5°、δψ0=1°和 δh0=－3000 m、δθ0=－0．3°、δφ0=－0．3°、δV0=－100 m/s、δγ0=－0．5°、δψ0=－1°的工況;組合3和組合4分別對應于m(1－5%)、ρ(1－15%)、CL(1－10%)、CD(1－10%)和m(1+5%)、ρ(1+15%)、CL(1+10%)、CD(1+10%)的工況。組合工況下仿真所得的實際狀態量、控制量和過程約束隨飛行時間的變化過程分別如圖2～圖4所示。可以看出，在閉環制導律的作用下，實際飛行軌跡未違背過程約束，終端狀態量的偏差值滿足制導精度的要求。這說明本文的制導律對初始狀態量的較大范圍偏差和飛行環境參數的有限擾動有良好的魯棒性。

從表3可以看出，飛行環境參數擾動對制導精度的影響比初始狀態量偏差的影響大，制導律在飛行器升阻比存在常值偏差情況下的魯棒性有限。這是因為在只存在初始狀態量偏差情況下，制導律解算用的系統模型很準確;而在飛行環境參數存在擾動時，制導律解算所使用的系統模型中的飛行環境參數與真實的飛行環境參數不一致，制導律解算用的系統模型不準確，從而影響制導律的魯棒性和制導精度。同時由于控制量數目的限制和過程約束對控制量調節能力(主要是攻角調節能力)的限制，當飛行環境參數存在擾動時，跟蹤指定的標稱剖面的制導律比跟蹤全狀態標稱軌跡的制導律的魯棒性更強［7］。

4．3制導律實時性驗證

為了驗證本文的滑翔段制導律的實時性，選擇初始狀態量偏差值為δV0=100 m/s的工況，LGR點分別選10、20、30、40和50，對LGR點數目不同的五種工況的仿真結果和制導律解算時間進行對比。從表4可以看出，五種工況的結果都滿足制導精度要求;假定在其它計算誤差相同的情況下，當LGR點的數目增加到一定量級時，終端狀態量偏差值不會再有特別顯著的變化;五種工況下制導律平均解算時間和最大解算時間均小于1 s的制導周期，充分說明制導律滿足實時性要求。隨著LGR點數目增加，制導律平均解算時間也增加，這是因為制導律解算中主要耗時的計算是式(24)中矩陣Ve的求逆，在滿足制導精度要求的前提下，選取少量數目的LGR點能夠減少矩陣Ve的維數，從而減少制導律解算的耗時。

圖2 狀態量隨飛行時間的變化過程Fig．2 Time histories of state variable

表3 初始狀態量存在偏差和飛行環境參數有擾動情況下的終端狀態量偏差Table 3 Terminal states deviations of biased initial states and disturbed flight environment parameters conditions

續表

表4 LGR點數目不同時的終端狀態量偏差和制導律解算時間Table 4 Terminal states deviations and calculation time of guidance law for different numbers of LGR points

圖3 控制量隨飛行時間的變化過程Fig．3 Time history of control variable

圖4 過程約束隨飛行時間的變化過程Fig．4 Time histories of path constraints

5 結論

本文提出了一種新的最優反饋控制求解方法——間接Radau偽譜法，并基于此方法提出了一種高超聲速飛行器滑翔段標稱軌跡跟蹤制導律。制導律解算不存在迭代、在線積分求解Riccati矩陣微分方程等復雜的數值計算，也不包含在線增益調參的過程，便于在線執行;對于不同情況，軌跡跟蹤控制器的結構和參數不需要改變。數值仿真結果表明，制導律能滿足實時性的要求，在初始狀態量存在較大范圍偏差和飛行環境參數存在有限擾動的情況下具有良好的魯棒性;在保證制導精度的前提下，通過選擇合適的離散點數目，能夠減少制導律的解算時間。

［1］ 胡正東，郭才發，蔡洪．天基對地打擊動能武器再入解析預測制導技術［J］．宇航學報，2009，30(3):1039－1044．［Hu Zheng-dong，Guo Cai-fa，Cai Hong．Analytical predictive guidance for space-to-ground kinetic weapon in reentry［J］．Journal of Astronautics，2009，30(3):1039－1044．］

［2］ 水尊師，周軍，葛致磊．基于高斯偽譜方法的再入飛行器預測校正制導方法研究［J］．宇航學報，2011，32(6):1249－1255．［Shui Zun-shi，Zhou Jun，Ge Zhi-lei．On-line predictorcorrector reentry guidance law based on Gauss pseudospectral method［J］．Journal of Astronautics，2011，32(6):1249－1255．］

［3］Xue S B，Lu P．Constrained predictor-corrector entry guidance［J］．Journal of Guidance，Control，and Dynamics，2010，33(4):1273－1281．

［4］ 陳剛，康興無，閆桂榮，等．基于偽譜方法的自適應魯棒實時再入制導律研究［J］．系統仿真學報，2008，20(20):5623－5634．［Chen Gang，Kang Xing-wu，Yan Gui-rong，et al．Real time robust adaptive reentry guidance law based on pseudospectral method［J］．Journal of System Simulation，2008，20(20):5623－5634］．

［5］Roenneke A J，Cornwell P J．Trajectory control for a low-lift re-entry vehicle［J］．Journal of Guidance，Control，and Dynamics，1993，16(5):927－933．

［6］ 南英．航天飛機再入抗干擾軌道優化設計［D］．西安:西北工業 大 學，1990．［Nan Ying．Anti-disturbance reentry trajectory optimization design for the Space Shuttle［D］．Xi’an:Northwestern Polytechnical University，1990．］

［7］Lu P．Regulation about time-varying trajectories:precision entry guidance illustrated［J］．Journal of Guidance，Control，and Dynamics，1999，22(6):784－790．

［8］Tian B L，Zong Q．Optimal guidance for reentry vehicles based on indirect legendre pseudospectral method［J］．Acta Astronautica，2011，68(7):1176－1184．

［9］Zhou H，Tawfiqur R，Wang D W，et al．Onboard pseudospectral guidance for re-entry vehicle［J］．Proceedings of the Institution of Mechanical Engineers，Part G:Journal of Aerospace Engineering，2014，228(11):1925－1936．

［10］ 吳旭忠，唐勝景，郭杰，等．基于滾動時域控制的再入軌跡跟蹤制導律［J］．系統工程與電子技術，2014，36(8):1602－1608．［Wu Xu-zhong，Tang Sheng-jing，Guo Jie，et al．Trajectory tracking guidance law for reentry based on receding horizon control［J］．Systems Engineering and Electronics，2014，36(8):1602－1608．］

［11］Yan H，Fahroo F，Ross I M．Optimal feedback control laws by legendre pseudospectral approximations［C］．The American Control Conference，Arlington，USA，June 25－27，2001．

［12］Phillips T H．A common aero vehicle(CAV)model description and employment guide［R］．Schafer Corporation for AFRL and AFSPC，2003．

［13］Garg D，Patterson M A，Darby C L，et al．Direct trajectory optimization and costate estimation of finite-horizon and infinitehorizon optimal control problems using a Radau Pseudospectral Method［J］．Computational Optimization and Applications，2011，49(2):335－358．

［14］Bryson A E，Ho Y C．Applied optimal control:optimization，estimation and control［M］．Washington DC:Hemisphere Publishing Company，1975．