徐文螢, 江駒, 甄子洋, 萬(wàn)月豐
(南京航空航天大學(xué) 自動(dòng)化學(xué)院,江蘇 南京 210016)
近空間可變翼飛行器(near-space vehicle,NMV)是一種翼形可變的工作于近空間的飛行器,其外形結(jié)構(gòu)呈三角形,采用翼身融合構(gòu)型。為了適應(yīng)不同飛行環(huán)境,飛行器通過(guò)伸縮翼來(lái)改變氣動(dòng)外形,完成不同的飛行任務(wù),具有重要的研究意義與軍事應(yīng)用價(jià)值,受到了國(guó)內(nèi)外眾多學(xué)者和研究機(jī)構(gòu)的重視[1-3]。
近空間飛行器軌跡優(yōu)化一般是非線性、帶有控制約束及狀態(tài)約束的最優(yōu)控制問(wèn)題[4]。由于本文研究對(duì)象的特殊性,在控制變量中加入小翼伸縮量,所以該問(wèn)題可以轉(zhuǎn)變?yōu)槭芏嗉s束條件下的高度非線性最優(yōu)控制問(wèn)題。該類最優(yōu)控制問(wèn)題根據(jù)求解方式分成兩大類。1)間接法:將問(wèn)題轉(zhuǎn)化為Hamilton邊值問(wèn)題而后求解,此類方法在滿足最優(yōu)性的一階必要條件同時(shí)求解的精度高,但存在收斂域小x(t)、難以估算共軛變量初值等缺點(diǎn);2)直接法:將時(shí)序連續(xù)的問(wèn)題轉(zhuǎn)化為離散的NLP問(wèn)題來(lái)獲取最優(yōu)解[5],存在求解精度較低等缺點(diǎn)。
近年來(lái),在解決軌跡優(yōu)化問(wèn)題上,與其他方法相比,直接偽譜法表現(xiàn)出全局收斂及高效率、高精度,因此在航空航天領(lǐng)域得到了廣泛的應(yīng)用[6-7]。文獻(xiàn)[8-9]采用Gauss偽譜法解決在多約束下的滑翔軌跡優(yōu)化問(wèn)題;文獻(xiàn)[10]應(yīng)用Gauss偽譜法解決了臨近空間飛行器上升段的軌跡優(yōu)化問(wèn)題。文獻(xiàn)[11]應(yīng)用一種改進(jìn)的點(diǎn)位經(jīng)緯度變換條件GPM的方法,提升了航跡規(guī)劃。文獻(xiàn)[12] 采用Gauss偽譜法針對(duì)可變翼飛行器小翼伸縮與最省燃油軌跡進(jìn)行優(yōu)化。
上述優(yōu)化方法存在一定的局限性:1)當(dāng)目標(biāo)函數(shù)復(fù)雜多變時(shí),須增加優(yōu)化函數(shù)的配點(diǎn)數(shù),才可在一定程度上提升優(yōu)化精度;2)隨著配點(diǎn)數(shù)的增多,狀態(tài)微分矩陣的維數(shù)也隨即增大,導(dǎo)致計(jì)算量大、效率低;3)收斂速度會(huì)因控制變量的不連續(xù)而降低。針對(duì)上述局限性,本文提出一種改進(jìn)的多段整合優(yōu)化的高斯偽譜法優(yōu)化策略,首先在時(shí)間區(qū)間上進(jìn)行逐一分段,而后采用低階多項(xiàng)式插值逼近狀態(tài)變量,通過(guò)逐步增加分段的個(gè)數(shù),提高逼近精度。
考慮多約束條件的非線性最優(yōu)控制問(wèn)題,即Bolza問(wèn)題[13]:
(1)
式中:x(t)、u(t)分別為狀態(tài)變量和控制變量(優(yōu)化變量);t0、tf為起始時(shí)間和末端時(shí)間;E為等式約束;C表示不等式約束。
針對(duì)Bolza問(wèn)題通過(guò)引入時(shí)間變量τ,將時(shí)間區(qū)間[t0,tf]投射到[-1,1][14],采用非等間距的 Gauss 偽譜法求解:
(2)
將時(shí)序連續(xù)的無(wú)限維最優(yōu)控制問(wèn)題轉(zhuǎn)換得到有限維的NLP問(wèn)題:
(3)
采用全局Lagrange 多項(xiàng)式插值,動(dòng)力學(xué)微分方程和控制變量離散化,路徑約束、邊界約束、終端狀態(tài)約束與目標(biāo)函數(shù)的精確近似??梢越LP問(wèn)題:
(4)
由于采用一般的高斯偽譜法進(jìn)行優(yōu)化時(shí),是將整個(gè)時(shí)間區(qū)間內(nèi)的狀態(tài)變量用一個(gè)全局Lagrange函數(shù)逼近,一般需要較多的配點(diǎn)才能達(dá)到所需的精度要求。為了減少配點(diǎn)數(shù),降低優(yōu)化時(shí)間,提高優(yōu)化精度,提出采用分段優(yōu)化,即將時(shí)間段劃分p個(gè)區(qū)間進(jìn)行優(yōu)化。記為:
[t0,t1],[t1,t2],…,[tp-1,tp]
(5)
式中:tp=tf,[tp-1,tp](p=1,2,…,P)表示第p個(gè)階段的時(shí)間段,記為階段p,其中tp(p=1,2,…,P-1)稱為分割點(diǎn)。由于對(duì)時(shí)間段劃分了p個(gè)區(qū)間后,其存在p-1個(gè)間斷點(diǎn),在每個(gè)間斷點(diǎn)上需要考慮其系統(tǒng)狀態(tài)的連續(xù)性問(wèn)題,因此分段優(yōu)化需要滿足以下連接約束[15-16]:
x(p+1)(tp)-x(p)(tp)=0,p=1,2,…,p-1
(6)
分段示例如圖1所示。圖1 中可以看出存在3個(gè)階段,其中包含2個(gè)連接約束點(diǎn)。每個(gè)連接點(diǎn)即為前一段的結(jié)束時(shí)刻和下一段的初始時(shí)刻??梢栽诿恳粋€(gè)優(yōu)化階段選擇合適的LG配點(diǎn)數(shù)np,(p=1,2,…,P-1),在每段上,根據(jù)LG配點(diǎn)擬合狀態(tài)量與控制量,給出每個(gè)子段間約束條件并建立目標(biāo)函數(shù)等,最后然后采用SQP算法求解。分段后的Bolza型問(wèn)題表示為:
圖1 分段連接點(diǎn)示意Fig.1 Segmental connection point diagram
性能指標(biāo)函數(shù):
min:Jp=Φ(x(1)(t0),t0,x(p)(tp),tp)+
(7)
狀態(tài)方程:
p=1,2,…,P
(8)
邊界約束:
E(x(1)(t0),t0,x(p)(tp),tp)=0
(9)
路徑約束:
C(x(p)(t),u(p)(t),t)≤0,t∈[tp-1,tp],
p=1,2,…,P
(10)
連接約束:
x(p+1)(tp)-x(p)(tp)=0
(11)
式中:x(p)∈Rn,u(p)∈Rm分別是階段p=1,2,…,P上的狀態(tài)變量和待優(yōu)化的控制變量。
針對(duì)爬升過(guò)程中小翼伸縮最省燃料軌跡優(yōu)化,綜合考慮爬升過(guò)程中聲速、密度、引力及推力變化等因素的影響,根據(jù)發(fā)動(dòng)機(jī)工作狀態(tài)及飛行狀態(tài)特性對(duì)爬升段合理劃分,建立多子段的優(yōu)化模型。
本文的研究對(duì)象相關(guān)外形參數(shù)如下:飛行器總質(zhì)量1.002×105kg,機(jī)身總長(zhǎng)60.69 m,機(jī)翼可伸縮面積30 m2,翼弦長(zhǎng)30 m,俯仰力矩慣性積8.47×106kg·m2。與一般的飛行器相比,機(jī)翼翼梢裝有可伸縮小翼。根據(jù)飛行需求,通過(guò)伸縮小翼改變翼型,引起不同的氣動(dòng)效應(yīng),適應(yīng)不同的飛行環(huán)境,完成各種飛行任務(wù)。
1)動(dòng)力學(xué)模型:可變翼飛行器爬升段數(shù)學(xué)模型[17]為:
(12)
(13)
式中:V、h、γ、m、α分別為飛行器速度、高度、航跡角、質(zhì)量、迎角;L、D、T分別為升力、阻力、發(fā)動(dòng)機(jī)推力;Isp、g0、r、μ分別為燃料比沖、海平面重力加速度、地球半徑、引力參數(shù);Ma為馬赫數(shù)。
2) 發(fā)動(dòng)機(jī)模型:
T=0.5ρV2(S+Sxiaoyi)CT
(14)
式中:S為飛行器有效參考面積;Sxiaoyi小翼的伸縮面積。CT為推力系數(shù),且:
(15)
式中β為發(fā)動(dòng)機(jī)節(jié)流閥調(diào)定值狀態(tài)。
3) 大氣環(huán)境模型[18]:
(16)
ρ=ρ0e-h/7 315.2
(17)
式中:r0為地球半徑;g0即標(biāo)準(zhǔn)重力加速度;ρ0=1.226 6 kg/m3,即地球表面的大氣密度。本文數(shù)值仿真驗(yàn)證中采用的標(biāo)準(zhǔn)重力加速度g0=9.806 65 m/s2,地球半徑r0=6 356 766 m。
4) 氣動(dòng)力模型:
飛行器飛行過(guò)程中受到阻力D和升力L的計(jì)算公式為:
(18)
式中:CL和CD分別為升力系數(shù)和阻力系數(shù),其中小翼伸縮時(shí)的升阻力系數(shù)通過(guò)插值法得到。圖2給出小翼伸、縮狀態(tài)下氣動(dòng)參數(shù)變化與α、Ma的三維關(guān)系曲面。
圖2 小翼伸縮前后升阻力系數(shù)Fig.2 Wing lift drag coefficient before and after expansion
結(jié)合近空間飛行環(huán)境模型,通過(guò)高精度插值法獲得小翼伸縮時(shí)的氣動(dòng)參數(shù),建立氣動(dòng)參數(shù)擬合數(shù)據(jù)庫(kù),根據(jù)建立的氣動(dòng)參數(shù)擬合數(shù)據(jù)庫(kù)查詢各項(xiàng)氣動(dòng)數(shù)據(jù),進(jìn)行飛行器非線性數(shù)學(xué)模型的后續(xù)計(jì)算,具體獲取過(guò)程如圖3所示。
圖3 氣動(dòng)數(shù)據(jù)獲取流程Fig.3 Flow chart of pneumatic parameter acquisition
綜上,通過(guò)氣動(dòng)數(shù)據(jù)庫(kù)查詢獲得各項(xiàng)氣動(dòng)數(shù)據(jù)之后,代入到各氣動(dòng)力和氣動(dòng)力矩中,即可進(jìn)行近空間飛行器非線性數(shù)學(xué)模型的后續(xù)計(jì)算。
近空間可變翼飛行器在爬升模態(tài)下需要經(jīng)歷平流層、中間層和部分熱層。隨著自身速度和高度的增加,飛行狀態(tài)發(fā)生顯著變化,飛行約束條件愈加苛刻。針對(duì)爬升模態(tài)下可變翼飛行器的參數(shù)狀態(tài)變化、發(fā)動(dòng)機(jī)推力變化等特性,將爬升模態(tài)劃分3個(gè)爬升子段進(jìn)行研究。
近空間可變翼飛行器從起飛段加速達(dá)到一定速度后開(kāi)始拉升進(jìn)入爬升段階段。在飛行爬升過(guò)程中,飛行器受到的發(fā)動(dòng)機(jī)推力較大,燃料消耗嚴(yán)重,飛行器的質(zhì)量、重心、慣性力矩等都會(huì)產(chǎn)生明顯的變化。根據(jù)爬升段發(fā)動(dòng)機(jī)工作狀態(tài)變化及飛行器高度、速度變化將其劃分成3個(gè)子段,表2給出了飛行器各分段的狀態(tài)及對(duì)飛行器的飛行要求。
表2 各分段的飛行器狀態(tài)及飛行要求Table 2 Aircraft status and flight requirements for each segment
依據(jù)飛行器在爬升時(shí)的狀態(tài)變化及發(fā)動(dòng)機(jī)的推力特性,將上升軌跡合理分段:將原最優(yōu)控制問(wèn)題轉(zhuǎn)化為多段最優(yōu)控制問(wèn)題后,采用改進(jìn)的多段高斯偽譜法進(jìn)行整合并進(jìn)行優(yōu)化計(jì)算。
2.4.1 邊界約束
近空間飛行器一般是通過(guò)助推器助推到預(yù)定的高度和速度后,啟動(dòng)發(fā)動(dòng)機(jī),給定飛行器初始與末端的狀態(tài):
(19)
2.4.2 路徑約束
飛行器在加速爬升過(guò)程中,受到大氣壓強(qiáng)作用,會(huì)破壞機(jī)體產(chǎn)生形變并摩擦生熱。因此,飛行過(guò)程中動(dòng)壓、熱流率及過(guò)載須滿足一定的約束[18]。
1)熱流率約束:
(20)
式中:C=7.968 6×10-5;p=0.5;q=3.15。
2)動(dòng)壓約束:
(21)
3)過(guò)載約束:
(22)
2.4.3 控制量與狀態(tài)量約束
飛行器執(zhí)行機(jī)構(gòu)在爬升過(guò)程中受限,因而控制量也受到一定的范圍內(nèi)約束,且爬升中滿足機(jī)動(dòng)要求狀態(tài)量相應(yīng)受到約束。對(duì)于可變翼飛行器,存在迎角、發(fā)動(dòng)機(jī)節(jié)流閥、小翼伸縮面積,u=[αβSxiaoyi],3個(gè)控制量和速度、高度、航跡角、質(zhì)量,x=[Vhγm],4個(gè)狀態(tài)量。
(23)
2.4.4 間斷點(diǎn)約束
針對(duì)劃分的3個(gè)子階段,出現(xiàn)的2個(gè)間斷點(diǎn),需要滿足如下間斷點(diǎn)約束條件:
(24)
飛行器在爬升飛行過(guò)程中,動(dòng)能和勢(shì)能增長(zhǎng)迅速,需要為其提供大量的能量。基于本文所研究的近空間可變翼飛行器,其在各個(gè)子段爬升過(guò)程中質(zhì)量的減少近似為燃料的消耗。由于飛行器執(zhí)行任務(wù)過(guò)程中攜帶的燃油有限,為了能節(jié)省更多燃料用于后續(xù)任務(wù),因此建立以最省燃料為優(yōu)化指標(biāo)進(jìn)行爬升,即飛行器在爬升末端質(zhì)量保持最大:
(25)
首先針對(duì)近空間可變翼飛行器劃分的P個(gè)子段,確定每個(gè)子段的初始配點(diǎn)數(shù)np,然后在每個(gè)子段上進(jìn)行離散化處理,將時(shí)間段[tp-1,tp]變換到[-1,1]上:
(26)
式中:τ表示時(shí)間t在[tp-1,tp]上轉(zhuǎn)化后的時(shí)間變量,p=1,2,…,P。
分別將x=[Vhγm]在每子段LG配點(diǎn)上進(jìn)行離散化:
(27)
(28)
對(duì)式子(27)狀態(tài)變量x兩邊分別求導(dǎo)可得:
(29)
(30)
其次同樣將控制變量u=[αβcSxiaoyi]在LG配點(diǎn)上進(jìn)行離散化:
(31)
(32)
最后同理將爬升段各段的約束條件、邊界點(diǎn)連接約束及性能指標(biāo)表示為[20-21]:
(33)
(34)
(35)
根據(jù)離散結(jié)果最后得到改進(jìn)的多段高斯偽譜法優(yōu)化模型:
(36)
式中:由于選取的優(yōu)化指標(biāo)為最省燃油爬升飛行,所以Jp可由式(25)Jm變換得到,C為不等式形式下的路徑約束包括熱流、動(dòng)壓、過(guò)載等,E1、E2為等式形式下的初始和終端的邊界約束。
改進(jìn)多段高斯維譜法優(yōu)化策略流程如圖4。
圖4 優(yōu)化策略算法流程Fig.4 Flow chart of the optimization strategy algorithm
利用上述改進(jìn)的多段高斯維譜法優(yōu)化策略進(jìn)行飛行器爬升段優(yōu)化。根據(jù)明確的物理特性對(duì)時(shí)間分段后,通過(guò)選取少量的LG配點(diǎn)進(jìn)行優(yōu)化,一方面有利于提高計(jì)算效率,另一方面保證優(yōu)化的精度與收斂性。
針對(duì)上述改進(jìn)的優(yōu)化策略,對(duì)小翼伸縮量及最省燃油軌跡進(jìn)行仿真實(shí)驗(yàn),每段參數(shù)設(shè)計(jì)如下:
1)各爬升子段初值與終值設(shè)置如下:
爬升速度:V0=590 m/s,Vtf=3 100 m/s;
爬升高度:h0=4 500 m,hf=36 000 m;
航跡角:γ0=0.12 rad,γtf=0 rad;
飛行器總質(zhì)量:m0=100 200 kg,max(mtf);
爬升時(shí)間:t0=0 s,tf=310 s。
2)控制變量u=[α,βc,Sxiaoyi]約束如下:
迎角:αmin=-3°,αmax=8°;
油門開(kāi)合狀態(tài):βmin=0,βmax=1;
小翼伸縮面積:Sxiaoyimin=0 m2,Sxiaoyimax=30 m2。
3)狀態(tài)變量限制范圍如下:
速度:Vmin=590 m/s,Vmax=3 150 m/s;
高度:hmin=4 500 m,hmax=36 500 m;
航跡角:γ0=0 rad,γtf=0.785 3 rad;
總質(zhì)量:mmax=100 200 kg。
4)飛行器在各個(gè)子段爬升過(guò)程中,每個(gè)子段的初始與前一子段的終端值必須保持平滑連接,因此各子段連接點(diǎn)的值域約束范圍如下:
動(dòng)壓約束:qmin=4×104Pa,qmax=1.6×104kPa;
過(guò)載約束:nmax=3.5。
6)其他參數(shù):飛行器參考面積S=360 m2,地表標(biāo)準(zhǔn)重力加速度g0=9.806 65 m/s2,地表大氣密度值ρ0=1.226 6 kg/m3,地球半徑r0=6 356 766 m,每個(gè)子段初始LG配點(diǎn)數(shù)均為5。
為了獲得爬升高精度優(yōu)化結(jié)果,當(dāng)進(jìn)行單段優(yōu)化時(shí),為了實(shí)現(xiàn)誤差低于10%,配點(diǎn)數(shù)N必須大于80,增大優(yōu)化時(shí)間,影響爬升優(yōu)化效率。采用多段優(yōu)化,將爬升時(shí)間劃分為3個(gè)子段[0,60]、[60,260]、[260,310]、LG最終收斂配點(diǎn)個(gè)數(shù)分別為n1=20,n2=35,n3=20,總優(yōu)化時(shí)長(zhǎng)為178.9 s,最大相對(duì)誤差為6.73%,在同等配點(diǎn)數(shù)情況下,多段優(yōu)化逼近比未分段的逼近精度提高了3.27%。仿真結(jié)果如圖5。
由圖5 (a)、(b)仿真速度、高度可知,飛行器在每個(gè)子段末都能穩(wěn)定達(dá)到最終收斂值或處于約束條件內(nèi)。隨著速度的不斷增加,初始段高度值先穩(wěn)定而后迅速增加。這種先加速后爬升的飛行軌跡有2大優(yōu)勢(shì):1)先加速可使飛行器快速達(dá)到安全速度;2)先加速后爬升可以滿足最省燃油優(yōu)化指標(biāo)。 圖5(c)為航跡角曲線,隨著爬升進(jìn)行,每子段內(nèi)航跡角基本都是先增大后減小,這也說(shuō)明速度矢量先是向垂直方向增加,然后往水平方向移動(dòng),先加速向上爬升而后爬升趨于緩慢,與每個(gè)階段速度、高度曲線相對(duì)應(yīng)。最終航跡角收斂為零,飛行器達(dá)到巡航速度高度,進(jìn)入平穩(wěn)巡航飛行階段。
圖5 分段爬升優(yōu)化狀態(tài)變量曲線Fig.5 Segmental climbing optimization state variable curves
圖6是3個(gè)優(yōu)化控制變量變化曲線。隨著飛行高度、速度的不斷增加,迎角在不斷減少,最后穩(wěn)定在一定的幅值范圍內(nèi)。每個(gè)階段中發(fā)動(dòng)機(jī)節(jié)流閥狀態(tài)都相對(duì)穩(wěn)定在一定幅度范圍內(nèi),但在階段連接點(diǎn)處存在跳變情況,跳變?cè)颍?)由于在階段的連接處,為了飛行狀態(tài)的平滑連接,油門存在突變;2)小翼在階段連接處伸縮變化引起油門大幅度的跳變。小翼在初始爬升段的時(shí)候,由于初始爬升速度較小,此時(shí)小翼伸縮對(duì)飛行器的升阻比影響較大,小翼縮回會(huì)更大程度上減少阻力的產(chǎn)生,優(yōu)化表明,爬升初始段,在速度較小的情況下,小翼選擇收回,更易于飛行器加速爬升減少飛行燃油消耗。切換過(guò)后,油門狀態(tài)保持較小的恒定值,對(duì)于加速爬升油量控制具有可觀的影響。飛行器在爬升至150 s左右時(shí),加速到高超聲速階段,由于發(fā)動(dòng)機(jī)推力在做轉(zhuǎn)換,且小翼伸縮量對(duì)升阻比影響較大,伸出小翼的同時(shí)氣動(dòng)焦點(diǎn)后移,使得穩(wěn)定性與機(jī)動(dòng)特性得到增強(qiáng)。因此小翼在150 s處從伸出到緩慢收回,此過(guò)程對(duì)飛行器油量節(jié)省也具有一定的影響。在最后一個(gè)子段,由于小翼的伸出可以增大阻力特性,為了降低油量消耗,采用伸出的小翼,此策略意在配合阻力使得加速度逐漸減少,最后結(jié)合油門使得飛行器達(dá)到最終收斂巡航速度和高度。
圖6 分段爬升優(yōu)化控制變量曲線Fig.6 Segmented climbing optimal control variable curves
圖5 (d)為質(zhì)量曲線,與固定翼狀態(tài)和伸縮翼對(duì)比可知,小翼伸縮可節(jié)省燃油爬升,在高超聲速段,小翼變化對(duì)阻力變化明顯,找到合適的小翼狀態(tài)會(huì)縮減燃料消耗。整個(gè)爬升過(guò)程中耗油約12 t,與小翼收回相比節(jié)省燃料約800~1 000 kg,相比完全展開(kāi)小翼節(jié)省燃料約1 600~2 000 kg,相當(dāng)于爬升段耗油量的10%,為飛行完成下個(gè)任務(wù)節(jié)省了足夠的油量。
1) 為了能模擬真實(shí)飛行環(huán)境與氣動(dòng)數(shù)據(jù),本文綜合考慮近空間區(qū)域聲速、密度、引力及推力等眾多變化因素對(duì)飛行器的影響,通過(guò)氣動(dòng)參數(shù)擬合數(shù)據(jù)庫(kù),使得建模更為準(zhǔn)確。
2) 設(shè)計(jì)并改進(jìn)多段高斯偽譜法優(yōu)化策略,有效解決爬升模態(tài)下的多段小翼伸縮軌跡優(yōu)化,在減少計(jì)算耗時(shí)的同時(shí)提高了優(yōu)化精度,進(jìn)一步提升多配點(diǎn)問(wèn)題的優(yōu)化效率。
3) 針對(duì)可伸縮小翼這一特性,對(duì)飛行器爬升模態(tài)下進(jìn)行最省燃油優(yōu)化,應(yīng)用伸縮翼達(dá)到節(jié)省大量的燃油,得到高效飛行收益,為今后可變翼飛行研究提供了參考,具有一定的工程應(yīng)用價(jià)值。