薛鵬飛, 劉小雄, 李創(chuàng), 馬青原

(西北工業(yè)大學(xué) 自動(dòng)化學(xué)院, 陜西 西安 710072)

基于高斯偽譜法的飛機(jī)下降段軌跡優(yōu)化

薛鵬飛, 劉小雄, 李創(chuàng), 馬青原

(西北工業(yè)大學(xué) 自動(dòng)化學(xué)院, 陜西 西安 710072)

摘要:以大型運(yùn)輸機(jī)為研究對(duì)象,進(jìn)行風(fēng)干擾下的飛機(jī)下降段軌跡優(yōu)化。通過對(duì)飛機(jī)下降過程的分析,建立了飛機(jī)的運(yùn)動(dòng)學(xué)模型、油耗模型和大氣模型,得出了問題的目標(biāo)函數(shù)和約束。運(yùn)用高斯偽譜法把最優(yōu)控制問題轉(zhuǎn)化為非線性規(guī)劃問題,并用序列二次規(guī)劃求解非線性規(guī)劃問題。仿真計(jì)算得到飛機(jī)最短下降時(shí)間的最優(yōu)軌跡,討論了不同的成本指數(shù)對(duì)下降時(shí)間和燃油消耗的影響。仿真結(jié)果表明,成本指數(shù)越大,飛機(jī)下降時(shí)間越短,燃油消耗越多。

關(guān)鍵詞:最優(yōu)控制; 序列二次規(guī)劃; 最優(yōu)軌跡; 成本指數(shù)

0引言

軌跡優(yōu)化技術(shù)作為未來飛機(jī)空中交通管理和飛行管理的核心技術(shù),已得到越來越多的重視。軌跡優(yōu)化有利于提高飛行器飛行品質(zhì)，以滿足既定任務(wù)要求;軌跡優(yōu)化技術(shù)的應(yīng)用可以緩解機(jī)場交通的擁擠情況,降低飛行過程中的燃油消耗,提高經(jīng)濟(jì)效益，并滿足人們的出行要求。

從早期的龐特里亞金最小值原理到新方法的出現(xiàn),對(duì)飛行軌跡優(yōu)化方法可歸納為兩大類:直接方法和間接方法。間接方法是以變分法和龐特里亞金原理為基礎(chǔ),在飛機(jī)的運(yùn)動(dòng)方程和狀態(tài)變量等約束下,對(duì)性能指標(biāo)函數(shù)求取極值。其優(yōu)點(diǎn)是解的精度較高,最優(yōu)解滿足一階最優(yōu)性必要條件;缺點(diǎn)是推導(dǎo)最優(yōu)解的過程比較復(fù)雜和繁瑣,求解兩點(diǎn)邊值問題時(shí)的收斂域小。近幾年,解決非線性最優(yōu)控制問題的偽譜法發(fā)展迅速,它同時(shí)將控制變量和狀態(tài)變量在時(shí)間區(qū)間上進(jìn)行離散,利用正交多項(xiàng)式對(duì)控制和狀態(tài)變量進(jìn)行逼近,將軌跡優(yōu)化問題轉(zhuǎn)化成了相應(yīng)的參數(shù)優(yōu)化問題。目前,求解參數(shù)優(yōu)化問題的方法有很多,序列二次規(guī)劃(Sequential Quadratic Programming,SQP)法是求解參數(shù)優(yōu)化問題最有效的方法,也是飛行器軌跡優(yōu)化的主流優(yōu)化算法[1]。偽譜法對(duì)數(shù)值迭代的初始值敏感性很低,收斂性非常好,已廣泛應(yīng)用于解決無人機(jī)航跡規(guī)劃、導(dǎo)彈制導(dǎo)和航天器軌道機(jī)動(dòng)等大量的最優(yōu)控制問題,成為目前解決最優(yōu)控制問題和軌跡最優(yōu)化問題最有效的工具之一[2]。

本文以大型運(yùn)輸機(jī)為研究對(duì)象,利用高斯偽譜法對(duì)飛機(jī)下降段進(jìn)行軌跡優(yōu)化,得到了飛機(jī)的最省燃油下降軌跡,并驗(yàn)證了偽譜法的快速收斂性。

1問題描述

1.1飛機(jī)運(yùn)動(dòng)學(xué)模型

考慮飛機(jī)質(zhì)量變化的三自由度運(yùn)動(dòng)學(xué)模型為:

(1)

式中:T,L,D分別為發(fā)動(dòng)機(jī)的推力、升力和阻力;Ffuel為單位燃油消耗量[3];Wx為風(fēng)場沿經(jīng)度方向的分量;Wy為風(fēng)場沿緯度方向的分量。

在一定油門開度下,發(fā)動(dòng)機(jī)的推力隨著飛行高度和飛行馬赫數(shù)的改變而改變,可以表示為:

T=T0(H,Ma)u

式中:u為油門開度。

1.2約束限制

飛機(jī)起飛后,要經(jīng)過特定的航路點(diǎn)才能到達(dá)著陸機(jī)場。在特定的航路點(diǎn),飛機(jī)完成轉(zhuǎn)彎、改變航向。飛機(jī)過航路點(diǎn)的約束一般表示為經(jīng)緯度和高度的約束,即:

式中:ε為約束的容忍誤差。飛機(jī)的航跡角約束為γ∈[-10°,10°];航向角約束為χ∈[-180°,180°];起始點(diǎn)狀態(tài)約束為X(t0)=X0。

1.3目標(biāo)函數(shù)

成本指數(shù)是飛行管理系統(tǒng)實(shí)施垂直導(dǎo)航功能計(jì)算的主要參數(shù)之一[4],是指航班成本中飛行時(shí)間引起的成本與燃油成本的比值,即:

式中:Ci為成本指數(shù);Ctime為飛機(jī)飛行1 h的時(shí)間成本;Cfuel為飛機(jī)消耗1 lb燃油的成本。

飛機(jī)軌跡優(yōu)化的目的是使上述直接成本最小,選擇不同的成本指數(shù)進(jìn)行優(yōu)化,得到在各個(gè)時(shí)刻飛機(jī)的狀態(tài)不同。本文將討論在有風(fēng)的情況下,不同的成本指數(shù)對(duì)飛機(jī)的高度、馬赫數(shù)等狀態(tài)量的影響,以及對(duì)飛機(jī)到達(dá)時(shí)間的影響。

1.4歸一化處理

為了提高算法在求解優(yōu)化問題時(shí)的收斂速度,將飛行器運(yùn)動(dòng)方程歸一化。定義無量綱物理量:

(2)

將式(2)代入式(1)中,可得到歸一化后的運(yùn)動(dòng)方程。

綜上所述,飛機(jī)的軌跡優(yōu)化問題可以歸結(jié)為在飛機(jī)運(yùn)動(dòng)學(xué)模型約束、中間過程約束、邊界條件約束和狀態(tài)變量控制變量約束條件下,以最小時(shí)間成本和到達(dá)時(shí)間為目標(biāo)函數(shù)的最優(yōu)控制問題,即:

(3)

2最優(yōu)控制問題求解

Gauss偽譜法將最優(yōu)控制問題的狀態(tài)變量和控制變量在一系列Legendre Gauss(LG)點(diǎn)上離散,并以離散點(diǎn)為節(jié)點(diǎn)構(gòu)造拉格朗日插值多項(xiàng)式來逼近狀態(tài)變量和控制變量。通過對(duì)全局插值多項(xiàng)式求導(dǎo)來近似狀態(tài)變量對(duì)時(shí)間的導(dǎo)數(shù),從而將微分方程約束轉(zhuǎn)換為一組代數(shù)約束。性能指標(biāo)中的積分項(xiàng)由Gauss積分計(jì)算。經(jīng)過上述變換后,可將最優(yōu)控制問題轉(zhuǎn)化為具有一系列代數(shù)約束的NLP問題[5-6]。

2.1時(shí)域變換

最優(yōu)控制問題的時(shí)間定義為區(qū)間[t0,tf],但偽譜方法的時(shí)間區(qū)域?yàn)閇-1,1],因此要將時(shí)間從[t0,tf]變換到為[-1,1],對(duì)時(shí)間變量t作變換:

(4)

即將t∈[t0,tf]映射到τ∈[-1,1]。

2.2離散化

Gauss偽譜法的Legendre-Gauss配點(diǎn)位于開區(qū)間τ∈(-1,1)中,這些點(diǎn)不包括初始值點(diǎn)和終端值。N個(gè)離散點(diǎn)包括N-2個(gè)內(nèi)部LG配點(diǎn),初始值點(diǎn)τ0≡-1,終端值點(diǎn)τf≡1。利用N-1個(gè)Lagrange插值基函數(shù)組成多項(xiàng)式對(duì)狀態(tài)變量進(jìn)行逼近:

利用N-1個(gè)Lagrange插值基函數(shù)組成的多項(xiàng)式對(duì)控制變量進(jìn)行逼近,表示為:

在配點(diǎn)上的正交分配方程為:

終端狀態(tài)應(yīng)滿足動(dòng)力學(xué)方程約束:

將終端狀態(tài)約束條件離散并用Gauss積分來近似,可得:

式中:wk為Gauss正交權(quán)重;τk為LG配點(diǎn)。

目標(biāo)函數(shù)包括終端性能指標(biāo)和積分型性能指標(biāo),即:

其中,積分型性能指標(biāo)可以用Gauss積分來近似,離散化后的目標(biāo)函數(shù)為:

經(jīng)過上述離散化處理,Gauss偽譜法將連續(xù)型最優(yōu)控制問題轉(zhuǎn)化為離散型NLP問題,轉(zhuǎn)化的NLP問題可以表示為:

(5)

式中：k=1,…,N。

3SQP法求解非線性規(guī)劃問題

SQP的基本思想是:在每一迭代通過求解一個(gè)二次規(guī)劃子問題(QP)來確立一個(gè)下降方向,以減少價(jià)值函數(shù)來取得步長,重復(fù)這些步驟直到求得原問題的解。SQP算法具有整體收斂和局部超一次收斂的特性,被認(rèn)為是目前求解NLP問題的最有效的算法之一[7]?？紤]NLP問題:

式中:f(x),h(x)，g(x),分別為目標(biāo)函數(shù)、等式約束和不等式約束。

在主迭代過程中,搜索方向dk由下述QP子問題確定:

SQP算法序列地求解一系列的二次規(guī)劃子問題,逐步得到原問題的最優(yōu)解。

4算例分析

4.1參數(shù)設(shè)置

某大型運(yùn)輸機(jī)的機(jī)體性能參數(shù)為:最大起飛重量265.5 t;機(jī)翼面積353 m2;最大速度230.5 m/s;最小速度55 m/s;最大航跡角10°;最小航跡角-10°;翼展50.29 m;升限13.715 km;最大迎角10°;最小迎角-10°;最大航向角180°;最小航向角-180°。

飛機(jī)在巡航段結(jié)束后進(jìn)入下降段,主航線下降點(diǎn)的起始航路點(diǎn)為東經(jīng)105.8°、北緯31.57°和高度8 371 m,中間航路點(diǎn)為東經(jīng)104.49°、北緯30.61°和高度3 342 m,終點(diǎn)航路點(diǎn)為東經(jīng)103.55°、北緯30.5°和高度552 m,飛機(jī)速度達(dá)到進(jìn)場速度59.2 m/s。高斯偽譜法具有很高的收斂性,但是選取的高斯配點(diǎn)數(shù)過少,求解問題的誤差較大;選取的高斯配點(diǎn)數(shù)過多,設(shè)計(jì)的變量數(shù)目比較龐大。設(shè)高斯點(diǎn)的個(gè)數(shù)為k,則設(shè)置的變量數(shù)目為10k+22。為了獲得精確解,同時(shí)防止計(jì)算量過大,本文選取高斯節(jié)點(diǎn)數(shù)為14個(gè),應(yīng)用Gauss偽譜法,把最優(yōu)控制問題轉(zhuǎn)換為NLP問題,采用SQP方法求解。

4.2最短時(shí)間軌跡

取Ci=9 999,得到飛機(jī)下降段的最短時(shí)間為30 min,仿真結(jié)果如圖1所示。可以看出,飛機(jī)經(jīng)過3個(gè)航路點(diǎn),高度平緩下降,飛行速度比較大,下降時(shí)間短;航跡角、航向角和迎角滿足約束條件。

圖1　最短時(shí)間軌跡優(yōu)化結(jié)果Fig.1　Optimization results of minimum time

4.3最小直接成本優(yōu)化軌跡

選擇不同的成本指數(shù)Ci,得到不同的到達(dá)時(shí)間和燃油消耗。分別選取Ci為500和1 000,得到對(duì)應(yīng)的最小直接成本軌跡,如圖2所示。從優(yōu)化的軌跡可以看出,飛機(jī)采用階梯下降的策略,航跡角、航向角和迎角滿足約束條件;Ci=500和Ci=1 000得到的時(shí)間均大于最短時(shí)間,飛行速度小于最短時(shí)間軌跡得到的速度。

表1對(duì)比了Ci=500和Ci=1 000時(shí)的燃油消耗和到達(dá)時(shí)間,結(jié)果顯示:飛機(jī)按照Ci=1 000比按照Ci=500優(yōu)化的軌跡進(jìn)行下降,到達(dá)時(shí)間少了0.8 min,燃油消耗多了119.2 kg。根據(jù)1.3節(jié)目標(biāo)函數(shù),假設(shè)燃油消耗的權(quán)重固定,到達(dá)時(shí)間的權(quán)重為Ci,Ci增大,表明到達(dá)時(shí)間的權(quán)重增大。優(yōu)化結(jié)果為到達(dá)時(shí)間減小,仿真結(jié)果與理論分析一致。

表1　不同成本指數(shù)軌跡優(yōu)化結(jié)果

5結(jié)束語

本文利用高斯偽譜法對(duì)大型運(yùn)輸機(jī)下降段的軌跡優(yōu)化問題進(jìn)行研究。研究結(jié)果表明，成本指數(shù)越大,飛機(jī)下降的時(shí)間越短,燃油消耗越多。研究內(nèi)容對(duì)于飛機(jī)控制經(jīng)濟(jì)成本具有一定的實(shí)際意義。

參考文獻(xiàn)：

[1]雍恩米,陳磊,唐國金.飛行器軌跡優(yōu)化數(shù)值方法綜述[J].宇航學(xué)報(bào),2008,29(2):397-406.

[2]胡松啟,陳雨.偽譜法在飛行器軌跡優(yōu)化中應(yīng)用分析[J].火箭推進(jìn),2014,40(5):61-68.

[3]閆國華,魏娜,楊曉軍.基于BADA模型的我國民用飛機(jī)持續(xù)進(jìn)近燃油經(jīng)濟(jì)性研究[C]//第五屆中國智能交通年會(huì)暨第六屆國際節(jié)能與新能源汽車創(chuàng)新發(fā)展論壇.北京:電子工業(yè)出版社,2009:411-414.

[4]楊俊,蘇彬.運(yùn)輸飛行成本指數(shù)計(jì)算[J].飛行力學(xué),2000,18(1):85-88.

[5]Benson D.A Gauss pseudo-spectral transcription for optimal control[D].Massachusetts:Massachusetts Institute of Technology,2005.

[6]雍恩米,唐國金,陳磊.基于Gauss偽譜方法的高超聲速飛行器再入軌跡快速優(yōu)化[J].宇航學(xué)報(bào),2008,29(6):1766-1772.

[7]Dehdari V,Oliver D S,Deutsch C V.Comparison of optimization algorithms for reservoir management with constraints—a case study[J].Journal of Petroleum Science and Engineering,2012,100:41-49.

(編輯：李怡)

Approach trajectory optimization for aircraft using Gauss pseudo-spectral method

XUE Peng-fei, LIU Xiao-xiong, LI Chuang, MA Qing-yuan

(School of Automation, Northwestern Polytechnical University, Xi’an 710072, China)

Abstract:With the large transport aircraft as the research object, this paper focuses on trajectory optimization in approach airspace under disturbance of wind. By analyzing of the approach stage, the dynamitic model, oil consumption model and aerodynamic force model of the aircraft were built. Then, the cost function and constraints of trajectory problem were proposed. The optimal control was translated into a nonlinear programming problem (NLP) by Gauss pseudo-spectral method, and the NLP was solved with the sequential quadratic programming. The optimal trajectory of the shortest descent time was obtained by simulating calculation. Meanwhile, this paper discussed the effect of the different cost indexes on descent time and fuel consumption. Simulation results show that the descent time is less when the cost index is high and fuel consumption is more.

Key words:optimal control; sequential quadratic programming; optimal trajectory; cost index

收稿日期:2015-09-10;

修訂日期:2015-12-11; 網(wǎng)絡(luò)出版時(shí)間:2016-02-29 16:37

基金資助:國家自然科學(xué)基金資助(61374032)；航空科學(xué)基金資助(20150753009)；民機(jī)專項(xiàng)基金資助(MJ-S-2011-16)

作者簡介:薛鵬飛(1991-)，男，河南焦作人，碩士研究生，研究方向?yàn)轱w行控制與軌跡優(yōu)化； 劉小雄(1973-)，男，陜西周至人，副教授，博士，研究方向?yàn)轱w行控制與仿真、容錯(cuò)控制。

中圖分類號(hào)：V249.1

文獻(xiàn)標(biāo)識(shí)碼：A

文章編號(hào)：1002-0853(2016)03-0034-05