宮梓豐，劉 剛，宋 睿，吳玉彬，蔡光斌

(1. 火箭軍工程大學(xué) 核工程學(xué)院， 陜西 西安 710025； 2. 火箭軍士官學(xué)校 測(cè)試控制系， 山東 青州 262500)

近幾年來，關(guān)于多約束高超聲速滑翔飛行器姿態(tài)控制問題已成為國(guó)內(nèi)外研究的一個(gè)重點(diǎn)[1]。而我國(guó)在高超聲速滑翔飛行器及其再入技術(shù)領(lǐng)域已經(jīng)取得了迅猛發(fā)展，但在試驗(yàn)演習(xí)和實(shí)戰(zhàn)化應(yīng)用方面仍有一定的差距?；栾w行器進(jìn)入再入階段后，不僅需要考慮動(dòng)壓、過載、熱流密度等約束從而對(duì)軌跡進(jìn)行規(guī)劃，還必須綜合實(shí)際打擊環(huán)境進(jìn)行避障與突防[2]。最新的制導(dǎo)方法都可以較為完善地應(yīng)對(duì)再入飛行中發(fā)生的各類情況，但對(duì)于高超聲速滑翔飛行器來說，要發(fā)揮出其重要的戰(zhàn)略突防性能，除卻要考慮相關(guān)制導(dǎo)方法涵蓋的因素，還必須能保證繞過具有地理敏感屬性的空域[3]?，F(xiàn)今對(duì)于禁飛區(qū)與飛行器突防的研究仍然相對(duì)較淺[4-5]。

通過建立高超聲速滑翔飛行器再入段的物理模型，將軌跡優(yōu)化中的過程約束和終端約束考慮在內(nèi)，使得整個(gè)航程滿足禁飛區(qū)約束，并盡可能優(yōu)化控制量是軌跡優(yōu)化與跟蹤控制問題的核心。此問題模型的非線性程度強(qiáng)，且復(fù)雜度高，常規(guī)的變分法、極大值原理等解析法難以快速求解出最優(yōu)解，所以一般利用數(shù)值法求解軌跡優(yōu)化與制導(dǎo)控制問題[6]。文獻(xiàn)[4]完整介紹了偽譜法作為直接法的典型代表方法，采用了Gauss積分規(guī)則的基本原理，而后利用Chebyshev多項(xiàng)式或者Legendre多項(xiàng)式產(chǎn)生配點(diǎn)，最后對(duì)系統(tǒng)進(jìn)行插值處理[7-8]。

采用hp-Radau自適應(yīng)偽譜法將整個(gè)時(shí)間區(qū)間分區(qū)并進(jìn)行網(wǎng)格細(xì)化，同時(shí)避免區(qū)間中產(chǎn)生過于細(xì)密的網(wǎng)格以及難以解算的高次復(fù)雜多項(xiàng)式[9]，綜合評(píng)估網(wǎng)格化密度和插值多項(xiàng)式的冪次，保證系統(tǒng)不僅可以快速收斂，還能夠滿足既定的精度要求[10]。此種方法可以降低初始條件估計(jì)區(qū)間誤差對(duì)結(jié)果的影響，將最優(yōu)控制問題轉(zhuǎn)化為等價(jià)的線性規(guī)劃問題[11-12]。本文經(jīng)過快速軌跡優(yōu)化，設(shè)計(jì)閉環(huán)反饋控制器以實(shí)現(xiàn)對(duì)系統(tǒng)狀態(tài)的跟蹤，通過仿真計(jì)算結(jié)果計(jì)算狀態(tài)量實(shí)時(shí)誤差，利用積分法設(shè)立誤差評(píng)判標(biāo)準(zhǔn)，評(píng)估模型性能得出結(jié)論。

1 再入模型構(gòu)建與控制器設(shè)計(jì)

1.1 滑翔飛行器再入運(yùn)動(dòng)學(xué)模型

高超聲速滑翔飛行器再入過程中會(huì)受到諸多的約束限制，控制系統(tǒng)通道間的耦合現(xiàn)象也較為突出，模型的非線性屬性較強(qiáng)，其動(dòng)態(tài)變化可由一組非線性微分方程表征[6，13]。在模型構(gòu)建中忽略地球的扁率，視其為均勻光滑的球體，可得出滑翔飛行器再入過程中的位置變化方程。

(1)

式中，x、y、z分別為飛行器在再入坐標(biāo)系下的x、y、z軸向所對(duì)應(yīng)的位置分量，v是飛行器的線速度，ζ為航跡傾角，ξ為航跡偏角。

飛行器的速度、航跡偏角和航跡傾角的一次微分滿足：

(2)

式中：Q、S、g、m、φ和σ分別代表飛行器的動(dòng)壓、參考面積、當(dāng)前位置重力加速度、質(zhì)量、緯度和傾側(cè)角；KD、KL、KY分別是擬合而成的飛行器阻力系數(shù)、升力系數(shù)和側(cè)向力系數(shù)，其數(shù)值一般與飛行器運(yùn)行的馬赫數(shù)和攻角相關(guān)；Ω=7.292×10-5rad/s，為地球的自轉(zhuǎn)角速度。

飛行器的運(yùn)動(dòng)學(xué)方程可以表示為：

(3)

式中：α和β是系統(tǒng)的控制量，分別為再入飛行器的攻角與側(cè)滑角；p、q、r分別是再入飛行器的滾轉(zhuǎn)方向角速度、俯仰方向角速度和偏航方向角速度[11]。

飛行器的動(dòng)力學(xué)方程可以表示為：

(4)

式中，b為飛行器的翼展長(zhǎng)度，Ixx、Iyy、Izz分別是再入飛行器相對(duì)于x、y、z軸向的轉(zhuǎn)動(dòng)慣量，Kl、Km、Kn分別是通過風(fēng)洞試驗(yàn)得出的或是擬合而成的滾轉(zhuǎn)、俯仰、偏航乘子系數(shù)。

1.2 限制優(yōu)化問題的約束條件

為使控制問題盡可能切合實(shí)際飛行任務(wù)要求，在軌跡優(yōu)化與制導(dǎo)問題中加入禁飛區(qū)約束。禁飛區(qū)是指飛行器禁止通過的空中區(qū)域，比如防空識(shí)別區(qū)、導(dǎo)彈攔截區(qū)等[14]，在優(yōu)化軌跡的時(shí)候必須要繞過禁飛區(qū)。

通過分析禁飛區(qū)約束的性質(zhì)[15]，設(shè)定禁飛區(qū)即路徑約束為：

C(x(t),u(t),t)≤0,t∈[t0,tf]

(5)

通常情況下，熱流密度、動(dòng)壓和過載被認(rèn)定是路徑約束的一部分，通過不等式組的形式給出[16-17]，其表達(dá)式可以按照通用標(biāo)準(zhǔn)定義，也可根據(jù)飛行器結(jié)構(gòu)特性來單獨(dú)定義。對(duì)于不同的飛行器外形與材料，以及不同的飛行任務(wù)會(huì)有不同的約束指標(biāo)，在此不再贅述。

1.3 優(yōu)化指標(biāo)

航跡參數(shù)的變化對(duì)于規(guī)避禁飛區(qū)有直接的影響，因而嚴(yán)格控制有關(guān)航程參數(shù)相對(duì)重要。為使滑翔飛行器在再入段獲得最大的橫向航程，則確定系統(tǒng)的目標(biāo)函數(shù)為：

(6)

2 問題轉(zhuǎn)化與求解策略

2.1 自適應(yīng)網(wǎng)格細(xì)化方法

數(shù)值法求解最優(yōu)控制問題要先對(duì)系統(tǒng)進(jìn)行時(shí)域轉(zhuǎn)換。其中，通過偽譜法求解需要將控制系統(tǒng)的時(shí)域映射到[-1,1]，所以先對(duì)原時(shí)域做相應(yīng)轉(zhuǎn)換。汲取hp有限元方法的優(yōu)點(diǎn)，hp自適應(yīng)偽譜法將時(shí)域分成了多個(gè)子區(qū)間，再將每個(gè)子區(qū)間上的最優(yōu)控制問題分別轉(zhuǎn)換為非線性規(guī)劃(NonLinear Programming, NLP)問題求解[18]。選取K-1個(gè)配點(diǎn),將最優(yōu)控制問題的整體運(yùn)行時(shí)間即t∈[t0,tf]分成K個(gè)子時(shí)間區(qū)間，將每個(gè)子區(qū)間的時(shí)間定義域[tk,tk+1]轉(zhuǎn)換到區(qū)間[-1,1]上，來滿足Legendre正交多項(xiàng)式的定義區(qū)間[6]。

(7)

通過這種方法可以將最優(yōu)控制問題轉(zhuǎn)換為非線性規(guī)劃問題，可以利用稀疏非線性優(yōu)化器(Sparse Nonlinear OPTimizer, SNOPT)算法求解此類問題。

設(shè)[tk-1,tk]為第k個(gè)子區(qū)間。取相鄰配點(diǎn)之間的中點(diǎn)，以式(8)所示準(zhǔn)則來評(píng)估誤差。

(8)

x(k)表示在第k個(gè)區(qū)間內(nèi)軌跡的狀態(tài)。設(shè)ρ為判斷實(shí)施細(xì)化標(biāo)準(zhǔn)的標(biāo)量。若該狀態(tài)所對(duì)應(yīng)的曲率小于ρ，則認(rèn)為該區(qū)間內(nèi)的軌跡相對(duì)平滑，如果這個(gè)區(qū)間的精度不滿足要求，則可以通過增加區(qū)間中的配點(diǎn)數(shù)來提升精度，即利用“p法”來提升求解精度；若存在大于ρ的元素，則認(rèn)定該區(qū)間的軌跡不平滑，要通過細(xì)化區(qū)間網(wǎng)格的方式提升精度，即采用“h法”來提升求解精度。

2.2 反步法求解控制策略

在得出最優(yōu)軌跡的同時(shí)，飛行器再入的全時(shí)狀態(tài)量也可以同步解算，進(jìn)而得出三通道輸出信號(hào)，控制副翼、方向舵和升降舵的運(yùn)行。

式(3)可以表示為：

(9)

式(4)可以表示為：

(10)

由此可以設(shè)計(jì)控制回路，構(gòu)建姿態(tài)控制器模型，如圖1所示。

圖1 跟蹤控制器結(jié)構(gòu)圖Fig.1 Structure of tracking control system

下面運(yùn)用Backstepping的方式設(shè)計(jì)控制u，實(shí)現(xiàn)對(duì)前一階段得出的參考控制量[σd,αd,βd]進(jìn)行跟蹤。Backstepping法又稱作反演、反推、反步法，其基本思路是根據(jù)系統(tǒng)的階數(shù)，分為若干步驟來實(shí)現(xiàn)控制，最終實(shí)現(xiàn)對(duì)系統(tǒng)目標(biāo)輸出的控制，對(duì)于二階及以上系統(tǒng)的控制優(yōu)勢(shì)明顯。在分步設(shè)計(jì)中，需要設(shè)計(jì)合適的 Lyapunov 函數(shù)，設(shè)計(jì)每一步的虛擬控制函數(shù)需要保證Lyapunov 函數(shù)一階導(dǎo)數(shù)負(fù)定，得到的虛擬控制函數(shù)即為下一階次的控制信號(hào)。從控制系統(tǒng)微分方程組的最低階次開始遞歸，直到最后一步設(shè)計(jì)出最終的控制方案。

考慮以下系統(tǒng)，其中g(shù)1、g2均可逆。

定義誤差變量z1=x1-ξd，ξd為給定的標(biāo)準(zhǔn)狀態(tài)量，對(duì)變量z1求導(dǎo)可得：

(11)

(12)

令x2=z2+x2d，得：

(13)

(14)

(15)

(16)

從式(16)可以看出，選取的Lyapunov函數(shù)一階導(dǎo)數(shù)負(fù)定，中間變量z1和z2都為有界的，使得變量x1和x2對(duì)參考信號(hào)[σd,αd,βd]的跟蹤誤差在可控的范圍內(nèi)，算法實(shí)現(xiàn)了對(duì)參考信號(hào)[σd,αd,βd]的跟蹤。

3 仿真試驗(yàn)

3.1 條件與約束

3.1.1 初始條件

飛行器再入的初速度v0=7.8 km/s，再入高度h0=80 km，初始航跡傾角ζ0=0°，初始航跡偏角ξ0=90°。各個(gè)方向的初始角速度分量均設(shè)為零。

3.1.2 禁飛區(qū)約束

根據(jù)實(shí)際飛行空域識(shí)別區(qū)限制條件，以飛行器再入點(diǎn)對(duì)地面的投影點(diǎn)為原點(diǎn)建立坐標(biāo)系，設(shè)計(jì)如式(17)所示的類圓柱形禁飛區(qū)域。

(17)

式中，lon代表經(jīng)度，lat代表緯度，h為高度。

3.1.3 終端約束

vf=800 m/s，hf=25 km，ζf=0°，ξf=90°。

3.1.4 過程約束

在飛行器再入的全過程中，定義法向過載最大值nmax=2.5g，熱流率最大值Q=1.2×106W/m2。

3.2 仿真結(jié)果與分析

采用上述已知條件，利用hp-Radau自適應(yīng)偽譜法對(duì)再入飛行問題進(jìn)行解算，仿真實(shí)驗(yàn)在MATLAB2018a軟件中運(yùn)行，設(shè)定配點(diǎn)精度為1×10-6。

通過軌跡優(yōu)化，可以得出高超聲速飛行器運(yùn)行的最優(yōu)軌跡，進(jìn)而得出最優(yōu)軌跡高度變化曲線，如圖2所示。

圖2 再入高度曲線Fig.2 Reentry height curve

飛行軌跡呈跳躍式，相較于平滑下降的“錢學(xué)森”彈道，其能耗更小，可以增加飛行航程。不僅如此，飛行器在稠密與稀薄的大氣層交替穿行，可以減少滑行過程中的熱能損耗，使得最優(yōu)控制問題中的熱能約束更容易被滿足，最優(yōu)解的選擇更寬泛。整個(gè)飛行過程中，飛行器速度下降過程偏于平穩(wěn)，除此之外，遵循這種飛行方案可以保證其熱流密度不會(huì)出現(xiàn)間歇尖峰，有利于保護(hù)飛行器的本體結(jié)構(gòu)不被損壞。

最優(yōu)軌跡的控制量(攻角)變化如圖3所示。

圖3 攻角變化曲線Fig.3 Curve of AOA

使用傾斜轉(zhuǎn)彎(Bank To Turn， BTT)飛行控制模式，即在航跡坐標(biāo)系中，設(shè)定攻角與側(cè)滑角為控制量，全程攻角控制在25°以內(nèi)，使得三通道控制信號(hào)耦合程度盡量輕微，從而容易有效實(shí)現(xiàn)系統(tǒng)控制。

加入禁飛區(qū)約束后，其飛行軌跡將要做出必要調(diào)整，原始軌跡與優(yōu)化后的飛行軌跡如圖4、圖5所示。

圖4 優(yōu)化前后軌跡立體圖Fig.4 Trajectory stereogram before and after optimization

圖5 優(yōu)化前后軌跡局部平面圖Fig.5 Trajectory partial layout before and after optimization

由圖4和圖5可知，飛行器沿著優(yōu)化軌跡飛行即可避開禁飛區(qū)，達(dá)成優(yōu)化指標(biāo)的同時(shí)完成突防任務(wù)。

根據(jù)狀態(tài)量與控制量的動(dòng)態(tài)關(guān)系，可以得出滾轉(zhuǎn)、俯仰、偏航三通道的角速度，進(jìn)而導(dǎo)出副翼、方向舵和升降舵的控制信號(hào)，如圖6所示。

圖6 副翼、方向舵和升降舵的控制信號(hào)Fig.6 Control signals of aileron, rudder and elevator

以此控制信號(hào)對(duì)三自由度角速度進(jìn)行跟蹤，跟蹤情況如圖7～9所示。

圖7 滾轉(zhuǎn)角速度的跟蹤情況Fig.7 Tracking of rolling angular velocity

圖8 俯仰角速度的跟蹤情況Fig.8 Tracking of pitch angular velocity

圖9 偏航角速度的跟蹤情況Fig.9 Tracking of yaw angular velocity

由圖7～9可知，系統(tǒng)對(duì)三個(gè)自由度的角速度變量能夠有效實(shí)現(xiàn)實(shí)時(shí)跟蹤，系統(tǒng)設(shè)立的采樣時(shí)間間隔約為40 s，由控制信號(hào)控制得到的角速度誤差因素主要來自目標(biāo)采樣和控制信號(hào)時(shí)間傳導(dǎo)。

對(duì)攻角和傾側(cè)角的跟蹤效果如圖10、圖11所示。

圖10 攻角的跟蹤情況Fig.10 Tracking of angle of attack

圖11 傾側(cè)角的跟蹤情況Fig.11 Tracking of bank angle

由圖10、圖11可以看出，跟蹤控制器已經(jīng)實(shí)現(xiàn)對(duì)于系統(tǒng)控制變量的實(shí)時(shí)跟蹤。在系統(tǒng)存在干擾的情況下，依然可以完成正常輸出，保證高超聲速飛行器再入的姿態(tài)穩(wěn)定。跟蹤信號(hào)的攻角變化較小，因而與之對(duì)應(yīng)的俯仰角速度會(huì)相對(duì)較小，由于采樣時(shí)刻俯仰角速度的突變量較小，所以系統(tǒng)對(duì)攻角的跟蹤效果較好；而傾側(cè)角與滾轉(zhuǎn)角速度和偏航角速度有關(guān)，由于偏航角速度和滾轉(zhuǎn)角速度存在較多的瞬時(shí)振蕩，尤其是在飛行器再入初始階段和飛行軌跡策略重新優(yōu)化的階段，因而跟蹤側(cè)滑角難度相對(duì)較大。由于每個(gè)采樣點(diǎn)間存在一定時(shí)間間隔，系統(tǒng)對(duì)其跟蹤效果會(huì)有一定誤差，但此誤差最終依然可以收斂。系統(tǒng)仿真驗(yàn)證了控制模型能夠達(dá)到預(yù)期輸出，可以有效跟蹤高超聲速飛行器再入時(shí)的控制變量，進(jìn)而實(shí)現(xiàn)機(jī)動(dòng)快速調(diào)姿。

4 結(jié)論

1)由偽譜法解算出的最優(yōu)軌跡狀態(tài)量符合條件限制，三通道信號(hào)雖然存在耦合，但可通過控制攻角范圍使系統(tǒng)能夠?qū)崿F(xiàn)有效跟蹤。

2)通過優(yōu)化原始軌跡，飛行器可以避開設(shè)定的禁飛區(qū)域，始末點(diǎn)的相對(duì)誤差在0.1%以內(nèi)。

3)角速度閉環(huán)跟蹤系統(tǒng)可以基本跟蹤最優(yōu)姿態(tài)變化，但是在突變過于明顯的時(shí)刻誤差較大，存在遲滯現(xiàn)象，這是由于系統(tǒng)采樣時(shí)刻相對(duì)固定，沒有通過感知環(huán)境變化而細(xì)化采樣造成的。

下一步工作將引進(jìn)自適應(yīng)采樣時(shí)間方法，減少因采樣時(shí)間遲滯所造成的誤差，提升跟蹤性能。