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