基于增益-調(diào)度控制和高保真飛機(jī)模型的實(shí)時(shí)四維航跡生成

Olusyo Ojemu , Mhdi Mhfouf , Lohithksh M. Miyr , Arr Al-Hindi ,Mihl Weiszer , Jun Chen ,*

? 2021 THE AUTHORS. Published by Elsevier LTD on behalf of Chinese Academy of Engineering and Higher Education Press Limited Company This is an open access article under the CC BY-NC-ND license(http://creativecommons.org/licenses/by-nc-nd/4.0/).

1. 引言

近來(lái)的研究表明，為了實(shí)現(xiàn)高效運(yùn)作，滿足機(jī)場(chǎng)不斷增加的飛機(jī)和客流量的需求，有必要對(duì)每架航班實(shí)施門到門的全程引導(dǎo)[1,2]。這種全程引導(dǎo)被稱為四維航跡（4DT）引導(dǎo)。四維航跡是指飛機(jī)從登機(jī)口到跑道的時(shí)空導(dǎo)航（包括滑行和回推路徑）[3]。四維航跡引導(dǎo)系統(tǒng)不僅能優(yōu)化地面運(yùn)動(dòng)的各個(gè)階段，而且還能協(xié)調(diào)機(jī)隊(duì)的行動(dòng)。四維航跡引導(dǎo)系統(tǒng)能減少高達(dá)55%的滑行延誤[4]。因此，許多空中交通管理人員目前正在研究如何提高機(jī)場(chǎng)的數(shù)字化水平，包括新的地面決策支持和導(dǎo)航技術(shù)，以應(yīng)對(duì)擁堵并提高機(jī)場(chǎng)運(yùn)行的魯棒性[1]。到2035年，乘客數(shù)量將增至2012年公布數(shù)量的1.5倍[1,5]，為了完成管理各個(gè)機(jī)場(chǎng)的飛機(jī)地面運(yùn)動(dòng)這一日益艱巨的任務(wù)，實(shí)現(xiàn)四維航跡算法變得更加緊迫。此外，飛機(jī)排放占全球溫室氣體排放的3% [5,6]，這為地面運(yùn)動(dòng)優(yōu)化提供了重要?jiǎng)恿Α?/p>

人們針對(duì)機(jī)場(chǎng)地面決策支持和導(dǎo)航系統(tǒng)的可行性進(jìn)行了大量的研究[7-10]。這些研究大致可以分為兩類。第一類包括優(yōu)化機(jī)場(chǎng)地面作業(yè)的路徑和調(diào)度的算法。這類研究更關(guān)注整個(gè)機(jī)場(chǎng)的運(yùn)行模式，而不是特定的飛機(jī)。這類算法主要解決的是滑行時(shí)間的最小化問題[4,11]，并不主動(dòng)或直接優(yōu)化燃油效率及其排放。少數(shù)例外情況，如主動(dòng)路徑和調(diào)度算法，直接考慮燃油消耗和排放[5,12]。然而，這些算法以簡(jiǎn)化的飛機(jī)運(yùn)動(dòng)方程和燃料排放模型為基礎(chǔ)，生成的四維航跡可能不符合操作限制，而且由于其與預(yù)設(shè)的四維航跡不完全一致，遵循這些軌跡可能會(huì)增加燃料的消耗和排放。另外，考慮計(jì)算的復(fù)雜性，用來(lái)生成四維航跡的算法通常是離線執(zhí)行的[5,9]，因此很難在出現(xiàn)全新事件時(shí)自動(dòng)生成新的軌跡。

第二類研究與控制算法有關(guān)，這些算法針對(duì)特定的飛機(jī)，可以實(shí)現(xiàn)預(yù)定的輸出以及研究人為因素的影響。Haus等[7]研究了自動(dòng)化系統(tǒng)如何幫助飛行員在滑行過程中進(jìn)行決策。為了成功生成四維航跡，必須將這兩種方法緊密結(jié)合起來(lái)進(jìn)行使用。具體來(lái)說，調(diào)度程序會(huì)根據(jù)某種最優(yōu)調(diào)度方法來(lái)確定最優(yōu)時(shí)限和路徑點(diǎn)[11]，而特定的飛機(jī)則以最優(yōu)方式遵循該最優(yōu)調(diào)度。

通常情況下，飛機(jī)能夠如實(shí)地遵循既定調(diào)度，這需要由飛行員來(lái)控制環(huán)路。有時(shí)，調(diào)度程序會(huì)給出一個(gè)最優(yōu)速度軌跡（如參考文獻(xiàn)[5,13]），但要求飛行員嚴(yán)格遵循速度軌跡是不現(xiàn)實(shí)的，因?yàn)轱w行員在滑行過程中的工作量非常大，而且在滑行操作中經(jīng)常要執(zhí)行一系列復(fù)雜的項(xiàng)目檢查工作。為了解決上述問題，研究人員開展了在滑行中使用全自動(dòng)系統(tǒng)的研究[8]。然而，由于在地面作業(yè)中采用自動(dòng)滑行系統(tǒng)仍具有很多問題，全自動(dòng)滑行模式?jīng)]有得到廣泛應(yīng)用[8]。但是，正如下一代航空運(yùn)輸系統(tǒng)（NextGen）[2]和歐洲單一天空空中交通管理研究（Single European Sky Air Traffic Management Research, SESAR）報(bào)告[1]中討論的那樣，為了滿足嚴(yán)格的四維航跡要求，就必須引入全自動(dòng)滑行系統(tǒng)。采用半自動(dòng)滑行模式可以規(guī)避這一問題。例如，參考文獻(xiàn)[8]主要研究了以人為中心的四維（4D）地面導(dǎo)航系統(tǒng)，其中飛行員仍負(fù)責(zé)控制飛行，信息通過滑道照明元件被傳遞給飛行員，以幫助其決策。這種方法被稱為followthe-greens方法[14]，目前倫敦希思羅機(jī)場(chǎng)和新加坡樟宜機(jī)場(chǎng)等全球各大機(jī)場(chǎng)都在使用該方法[7]。然而，這種半自動(dòng)滑行模式存在一些缺陷。第一，由于沒有采用全自動(dòng)滑行，這種模式并不能完全跟隨所要求的四維航跡。此外，此類系統(tǒng)在設(shè)計(jì)最優(yōu)控制系統(tǒng)時(shí)，通常會(huì)做出許多與實(shí)際情況不符的假設(shè)，如飛行員以怠速滑行。參考文獻(xiàn)[15]指出，滑行過程中的功率設(shè)置并不是固定的，通常取決于飛行員的操作行為。如果飛行員偏離這些假設(shè)，就會(huì)出現(xiàn)次優(yōu)滑行。另外，許多引導(dǎo)飛行員滑行的算法并沒有主動(dòng)降低燃料燃燒或發(fā)動(dòng)機(jī)的排放。例如，參考文獻(xiàn)[7]提出的引導(dǎo)系統(tǒng)確保了軌跡能夠被嚴(yán)格遵循，但并未考慮燃料和排放問題。

本文提出了一種新的自動(dòng)化系統(tǒng)。該系統(tǒng)能夠完全實(shí)現(xiàn)四維航跡導(dǎo)航，并且主要用于兩種情況：①作為反應(yīng)型決策支持工具，生成能夠解決全新事件的新軌跡；②充當(dāng)半自動(dòng)和全自動(dòng)滑行的自動(dòng)駕駛系統(tǒng)。本文所采用的方法是求解一個(gè)離線優(yōu)化問題，即找到基于比例-積分-微分（PID）的控制系統(tǒng)的增益，使特定目標(biāo)函數(shù)最小化。目標(biāo)函數(shù)一經(jīng)確定，就能和調(diào)度算法緊密結(jié)合，并且該方法考慮了高保真飛機(jī)模型的燃油消耗和排放。具體而言，該調(diào)度算法考慮了基于飛機(jī)模型的約束，確定了最優(yōu)的滑行路線，然后控制飛機(jī)的滑行，此過程利用了最優(yōu)PID增益，使其以最優(yōu)方式遵循調(diào)度。這種方法的優(yōu)點(diǎn)在于使用了一種直觀、簡(jiǎn)單而高效的比例-積分-微分控制系統(tǒng)。由于控制器的優(yōu)化在在線部署之前是離線運(yùn)行的，因此易于部署，且不需要在線優(yōu)化，減少了在線計(jì)算的繁重的工作量。因此，每當(dāng)有干擾使調(diào)度算法的原始計(jì)劃不可行時(shí)，通過該算法可以在線生成一個(gè)新的、高效的四維航跡。

就像飛機(jī)在飛行模式中的自動(dòng)駕駛一樣，該控制系統(tǒng)能被直接應(yīng)用在飛行控制裝置中。這種方法允許全自動(dòng)滑行，而不需要飛行員進(jìn)行干預(yù)。此外，該控制系統(tǒng)可以為飛行員提供決策支持，其方式類似于前文的follow-the-greens方法。特別地，該方法可以幫助飛行員確定適當(dāng)?shù)目刂戚斎虢M合，以確保目標(biāo)函數(shù)（燃料消耗）在半自動(dòng)滑行時(shí)最小化。參考文獻(xiàn)[16]全面回顧了最優(yōu)控制算法的使用。

本文其余部分的結(jié)構(gòu)如下：第2節(jié)提供了利用大型噴氣式飛機(jī)模型驗(yàn)證該算法的詳細(xì)過程；第3節(jié)對(duì)該算法進(jìn)行了描述并解釋了如何將其用于飛機(jī)的自動(dòng)滑行；第4節(jié)概述了將該方法應(yīng)用在復(fù)雜飛機(jī)模型上的結(jié)果；第5節(jié)總結(jié)全文，概述了所取得的成果，并對(duì)未來(lái)的研究方向提出建議。

2. 飛機(jī)模型建立

本節(jié)介紹了研究中使用的非線性波音747-100飛機(jī)模型的數(shù)學(xué)方程。該數(shù)學(xué)模型是由剛體運(yùn)動(dòng)的一般方程導(dǎo)出的狀態(tài)空間模型公式。根據(jù)牛頓力學(xué)，對(duì)于運(yùn)動(dòng)中的剛體，線速度和角速度（分別用V和Ω表示）的變化率與作用在該剛體上的合力和力矩有關(guān)，如下式所示：

式中，F(xiàn)= [Fx Fy Fz]T和M= [L M N]T分別是沿著機(jī)體坐標(biāo)系中x、y、z軸的力和力矩向量；V=[u v w]T是線速度向量，其中u、v、w分別是x、y、z軸上的速度；Ω= [p q r]T是角速度向量，其中p是滾轉(zhuǎn)角速度，q是俯仰角速度，r是橫擺角速度；m是飛機(jī)的質(zhì)量；t是時(shí)間；I是恒定慣性張量（飛機(jī)的力矩和慣性的乘積）[17]。式（1）可以改寫為如下狀態(tài)空間公式：

式（2）為標(biāo)準(zhǔn)狀態(tài)空間形式，一般可以寫成：

式中，x∈RK是相應(yīng)時(shí)間導(dǎo)數(shù)，用x·表示。其中R代表實(shí)數(shù)；K代表狀態(tài)的總數(shù)。u∈RW代表維度W的時(shí)變輸入；υ∈RZ是維度Z的擾動(dòng)向量（如大氣擾動(dòng)）。式（3）是非線性時(shí)變狀態(tài)方程，該方程對(duì)飛行器所受的力和力矩進(jìn)行了隱式建模。

式（3）的非線性時(shí)不變等效模型一般能模擬包括飛機(jī)在內(nèi)的任何實(shí)際的物理動(dòng)力系統(tǒng)。因此，本研究假設(shè)為非線性時(shí)不變模型。盡管慣性矩取決于許多其他因素，如飛機(jī)的質(zhì)量等，但它與時(shí)間（即飛機(jī)的機(jī)齡）無(wú)關(guān)。在本研究使用的波音747-100飛機(jī)模型中，狀態(tài)變量的個(gè)數(shù)為13。具體來(lái)說，這些狀態(tài)包括標(biāo)準(zhǔn)的12維變量（控制機(jī)身角速率和加速度）以及一個(gè)被稱為“預(yù)推力”（pre-thrust）的附加變量。這些狀態(tài)用于描述發(fā)動(dòng)機(jī)的動(dòng)態(tài)性能。2.2節(jié)對(duì)狀態(tài)方程作了簡(jiǎn)要描述。

2.1. 波音747-100飛機(jī)模型

波音747-100是一種四扇噴氣式大型運(yùn)輸機(jī)，于1967年投入使用。波音747-100配備一個(gè)舵柄和一個(gè)方向舵，在滑行時(shí)便于控制方向。假定推力由兩個(gè)內(nèi)部發(fā)動(dòng)機(jī)產(chǎn)生。波音747-100還配備了內(nèi)外克魯格襟翼和帶有四個(gè)副翼的可移動(dòng)穩(wěn)定器，由于在滑行過程中沒有使用它們，因此它們沒有被包含在模型中。研究中用于波音747-100建模的數(shù)據(jù)來(lái)源于參考文獻(xiàn)[18,19]。

作用在機(jī)體坐標(biāo)系上的力和力矩的設(shè)定如圖1所示。圖2顯示了在本研究中用于建模的飛機(jī)部件。發(fā)動(dòng)機(jī)和燃油系統(tǒng)是根據(jù)國(guó)際民用航空組織（ICAO）的數(shù)據(jù)庫(kù)進(jìn)行建模的。

2.2. 狀態(tài)方程

如前文所述，非線性狀態(tài)方程由13個(gè)非線性函數(shù)組成。這些狀態(tài)又細(xì)分為以下幾個(gè)子狀態(tài)：

式中，x1、x2、x3、x4、x5分別代表旋轉(zhuǎn)速度、平移速度、姿態(tài)、位置和預(yù)推力子向量。x1= [p q r]T、x2= [V α β]T、x3= [φ θ ψ]T、x4= [h x y]T、x5=預(yù)推力。x1取決于滾轉(zhuǎn)角速度p、俯仰角速度q和橫擺角速度r；x2是線速度向量V、迎角α和側(cè)偏角β的函數(shù)；姿態(tài)向量x3是滾轉(zhuǎn)角φ、俯仰角θ和偏航角ψ的函數(shù)；位置向量x4由高度h和位置x、y組成。旋轉(zhuǎn)速度和平移速度分別是在飛機(jī)的機(jī)體坐標(biāo)系和氣流坐標(biāo)系內(nèi)定義的，而姿態(tài)和位置是在地面坐標(biāo)系中定義的。最后一個(gè)變量表示推力，在飛機(jī)的機(jī)體坐標(biāo)系中也有給出。將向量狀態(tài)導(dǎo)數(shù)與向量狀態(tài)和輸入聯(lián)系起來(lái)的向量函數(shù)f如下式所示：

圖1. 飛機(jī)的力和力矩的設(shè)定。

圖2. 飛機(jī)滑行部件及其設(shè)定。ICAO：國(guó)際民用航空組織。

式中，f1、f2、f3、f4和f5分別代表旋轉(zhuǎn)、平移、姿態(tài)、位置和預(yù)推力向量與狀態(tài)變量之間的向量子函數(shù)。除了預(yù)推力變量外，每個(gè)子函數(shù)都源自參考文獻(xiàn)[20]，子函數(shù)詳見附錄A。

在波音747-100飛機(jī)模型中，發(fā)動(dòng)機(jī)增壓比（EPR）或預(yù)推力決定了飛機(jī)產(chǎn)生的推力。這一推力由時(shí)間常數(shù)為5 s的一階動(dòng)態(tài)系統(tǒng)表示。為了使其與其他狀態(tài)方程兼容，并且為了便于模型的建立，需將其轉(zhuǎn)換為狀態(tài)空間形式。對(duì)于給定的靜態(tài)EPR（即穩(wěn)態(tài)值）[這一數(shù)值取決于節(jié)流閥調(diào)定（thr）]，預(yù)推力可由以下公式得到：

式中，g代表靜態(tài)EPR，是與節(jié)流閥調(diào)定和給定狀態(tài)有關(guān)的函數(shù)。如2.3節(jié)所示，已知預(yù)推力，可以計(jì)算出發(fā)動(dòng)機(jī)產(chǎn)生的推力。狀態(tài)方程適用于任何飛機(jī)，但是每一架飛機(jī)有不同的力和力矩值。對(duì)力和力矩分量以及所使用的數(shù)值推導(dǎo)見2.3節(jié)。

2.3. 力和力矩的推導(dǎo)

本節(jié)給出了波音747-100飛機(jī)模型的力和力矩方程，包括氣動(dòng)元件、發(fā)動(dòng)機(jī)元件和起落架元件，其數(shù)據(jù)來(lái)自參考文獻(xiàn)[18,19]。

2.3.1. 氣動(dòng)力和力矩

當(dāng)飛機(jī)滑行時(shí)，特別是在相對(duì)高的滑行速度下，氣動(dòng)力和力矩是不可忽略的。力和力矩的氣動(dòng)力分量都是在氣流坐標(biāo)系上計(jì)算的，然后被轉(zhuǎn)換到機(jī)體坐標(biāo)系上。氣動(dòng)力數(shù)據(jù)集以無(wú)單位特征系數(shù)C*的形式得到。如附錄A所示，下標(biāo)*是指定的6種氣動(dòng)力分量的6個(gè)下標(biāo)之一。

2.3.2. 發(fā)動(dòng)機(jī)力和力矩

波音747-100飛機(jī)模型有4臺(tái)Pratt和Whitney JT9D-3發(fā)動(dòng)機(jī)，起飛推力為43 500 lbf (193 500 N)。在整個(gè)模擬過程中假定只有兩臺(tái)發(fā)動(dòng)機(jī)運(yùn)行。因?yàn)榛羞^程中通常不使用反向推力[它只在著陸到大約100 kn（51 m·s-1）時(shí)使用]，所以飛機(jī)發(fā)動(dòng)機(jī)的建模過程不考慮反向推力。按以下步驟計(jì)算發(fā)動(dòng)機(jī)產(chǎn)生的最終推力。

（1）首先得到節(jié)流閥角度，然后利用這個(gè)角度計(jì)算發(fā)動(dòng)機(jī)油門桿角度（PLA）。需要注意的是，本文已經(jīng)對(duì)開發(fā)的波音747-100飛機(jī)模型做了一些修改，以規(guī)范節(jié)流閥調(diào)定系統(tǒng)，使輸入系統(tǒng)中的是一個(gè)在0和1之間的數(shù)字。在波音747-100飛機(jī)模型中，介于60°～130°的節(jié)流閥角度被轉(zhuǎn)換為0～1的數(shù)字。使用以下公式對(duì)其進(jìn)行標(biāo)準(zhǔn)化。

（2）給定發(fā)動(dòng)機(jī)油門桿角度（或使thr規(guī)范化為0～1，如模型所示），以確定靜態(tài)EPR。與thr和EPR相關(guān)的多項(xiàng)式模型如下式所示：式中，k0、k1、k2、k3是多項(xiàng)式系數(shù)，如表1所示。

表1 計(jì)算EPR的多項(xiàng)式函數(shù)的系數(shù)

（3）如前文所述，發(fā)動(dòng)機(jī)包括一個(gè)瞬態(tài)系統(tǒng)，該系統(tǒng)由一個(gè)時(shí)間常數(shù)為5 s的一階動(dòng)態(tài)系統(tǒng)表示。最后一步計(jì)算出的靜態(tài)EPR表示EPR的參考水平或穩(wěn)態(tài)。值得強(qiáng)調(diào)的是，瞬態(tài)系統(tǒng)可由狀態(tài)空間形式表示，因此如前文所述，需要增加一個(gè)狀態(tài)預(yù)推力。從本質(zhì)上講，預(yù)推力變量代表了動(dòng)態(tài)EPR。

（4）根據(jù)參考文獻(xiàn)[18,19,21]，美國(guó)國(guó)家航空航天局（NASA）提供了隨著飛機(jī)馬赫數(shù)增加，在動(dòng)態(tài)EPR中增加增量EPR的圖。然而，本次模型并未考慮EPR的增量[因?yàn)榛兴俣韧ǔＰ∮?0 kn（21 m·s-1），所以速度增加帶來(lái)的EPR的增量可以忽略不計(jì)]。

（5）當(dāng)EPR已知（EPR本質(zhì)上是動(dòng)態(tài)的）時(shí)，最后一步是計(jì)算每個(gè)發(fā)動(dòng)機(jī)推力的凈推力。凈推力是與發(fā)動(dòng)機(jī)EPR和環(huán)境條件有關(guān)的函數(shù)。凈推力的推導(dǎo)公式如下：

式中，斜率為90°，截距為-90°。

由于在所開發(fā)的模型中使用國(guó)際單位制（SI）進(jìn)行計(jì)算，因此以磅（lb, 1 lb = 0.453592 kg）為單位的推力乘以因子4.44822后被轉(zhuǎn)換為牛頓。上述步驟中推導(dǎo)出的推力僅適用于一臺(tái)發(fā)動(dòng)機(jī)。綜上所述，計(jì)算每臺(tái)發(fā)動(dòng)機(jī)推力的過程如圖3所示。

怠速推力：在節(jié)流閥調(diào)定為零時(shí)，EPR的近似值為1.02，這使得發(fā)動(dòng)機(jī)產(chǎn)生正的凈推力。同樣值得注意的是，發(fā)動(dòng)機(jī)會(huì)在飛機(jī)的重心附近產(chǎn)生力矩。在雙引擎滑行模式（使用內(nèi)引擎，關(guān)閉外引擎）中，發(fā)動(dòng)機(jī)推力產(chǎn)生的凈力矩為零，而在單引擎滑行模式中，發(fā)動(dòng)機(jī)將產(chǎn)生一個(gè)力矩使飛機(jī)偏航。這個(gè)力矩由力和發(fā)動(dòng)機(jī)離飛機(jī)重心的距離的乘積進(jìn)行計(jì)算。根據(jù)參考文獻(xiàn)[18,19]，這個(gè)距離等于39.167 ft（11.938 m）。

圖3. 計(jì)算凈推力的步驟。

燃油流量和排放：由于本研究的目的之一是尋找更有效的方法來(lái)幫助飛機(jī)滑行，因此準(zhǔn)確測(cè)量發(fā)動(dòng)機(jī)燃油消耗和排放至關(guān)重要?；贗CAO燃油和排放數(shù)據(jù)庫(kù)[22]，我們建立了一個(gè)發(fā)動(dòng)機(jī)油耗模型。ICAO數(shù)據(jù)庫(kù)已被廣泛應(yīng)用于研究在不同推力設(shè)置下的燃油消耗和發(fā)動(dòng)機(jī)排放[15]。該數(shù)據(jù)庫(kù)提供了1980年以后生產(chǎn)的大約500臺(tái)發(fā)動(dòng)機(jī)的列表，以及它們?cè)?種不同推力設(shè)置下的燃料消耗和排放情況。數(shù)據(jù)庫(kù)中所有發(fā)動(dòng)機(jī)的額定推力分布如圖4所示。

圖4. ICAO數(shù)據(jù)庫(kù)中發(fā)動(dòng)機(jī)額定推力的分布。

除參考文獻(xiàn)[22]所示的排放和燃料統(tǒng)計(jì)數(shù)據(jù)外，ICAO數(shù)據(jù)庫(kù)還包含與發(fā)動(dòng)機(jī)有關(guān)的其他重要信息，如停產(chǎn)狀態(tài)、燃料類型、發(fā)動(dòng)機(jī)制造商以及發(fā)動(dòng)機(jī)是否停用。數(shù)據(jù)庫(kù)列表中不包括JT9D-3發(fā)動(dòng)機(jī)，可能是由于該發(fā)動(dòng)機(jī)是在1980年之前制造的。為了便于計(jì)算JT9D-3發(fā)動(dòng)機(jī)的燃料燃燒和排放，我們對(duì)數(shù)據(jù)庫(kù)進(jìn)行了全面的統(tǒng)計(jì)分析。先前，已經(jīng)有人使用過模糊邏輯模型來(lái)分析這類數(shù)據(jù)[23-25]。研究發(fā)現(xiàn)，發(fā)動(dòng)機(jī)的額定推力和停產(chǎn)狀態(tài)對(duì)于預(yù)測(cè)燃料流量很重要。因此，這兩個(gè)變量可以被用來(lái)預(yù)測(cè)任意額定推力和停產(chǎn)狀態(tài)已知的發(fā)動(dòng)機(jī)燃油流量，如JTD9-3發(fā)動(dòng)機(jī)的情況。不同推力值（每臺(tái)發(fā)動(dòng)機(jī)4個(gè)推力值）下的燃料流量以及發(fā)動(dòng)機(jī)停產(chǎn)狀態(tài)構(gòu)成一個(gè)矢量，然后利用多項(xiàng)式回歸找出對(duì)應(yīng)函數(shù)，對(duì)燃料流量和排放進(jìn)行預(yù)測(cè)。將數(shù)據(jù)分為訓(xùn)練數(shù)據(jù)集和測(cè)試數(shù)據(jù)集（分別為60%和40%）。多項(xiàng)式回歸分析的結(jié)果如圖5（a）、（b）所示。JT9D-3發(fā)動(dòng)機(jī)的額定推力約為193.5 kN，發(fā)動(dòng)機(jī)目前已停產(chǎn)。將這兩個(gè)變量輸入多項(xiàng)式回歸模型中以預(yù)測(cè)燃料流量和排放。這些預(yù)測(cè)的結(jié)果如圖6（a）、（b）所示。圖6（a）中的燃料流量ffuel采用多項(xiàng)式方程計(jì)算，其系數(shù)如式（10）所示。圖6（b）中的一氧化碳（CO）排放指數(shù)fCO由式（11）計(jì)算得出。

2.3.3. 起落架力和力矩

參考文獻(xiàn)[18,19]中給出了用于計(jì)算起落架力和力矩的數(shù)據(jù)和公式。飛機(jī)起落架已經(jīng)簡(jiǎn)化為一個(gè)三輪車模型。在三輪車模型中，假定飛機(jī)具有三個(gè)起落架，將每個(gè)起落架視為一個(gè)非線性油壓桿模型。飛機(jī)有兩個(gè)主起落架（左和右）和一個(gè)前起落架。前起落架用于轉(zhuǎn)向，但僅能通過兩個(gè)主起落架進(jìn)行制動(dòng)。

圖5. 給定額定推力和非工作狀態(tài)下燃料流量的預(yù)測(cè)結(jié)果。（a）訓(xùn)練數(shù)據(jù)結(jié)果；（b）測(cè)試數(shù)據(jù)結(jié)果。RMS：均方根。

圖6. 波音747-100飛機(jī)模型的JT9D-3發(fā)動(dòng)機(jī)的預(yù)估。（a）燃料流量；（b）一氧化碳（CO）排放。

當(dāng)飛機(jī)在地面上運(yùn)動(dòng)時(shí)，地面、輪胎和油壓桿之間復(fù)雜的相互作用會(huì)產(chǎn)生影響飛機(jī)運(yùn)動(dòng)的力和力矩。由起落架產(chǎn)生的力和力矩分別在起落架的局部坐標(biāo)系上計(jì)算，并將其轉(zhuǎn)換到機(jī)體坐標(biāo)系上。同其他部分的力和力矩一樣，這些力和力矩取決于飛機(jī)的狀態(tài)。計(jì)算起落架力和力矩的步驟如下。

·計(jì)算油壓桿的壓縮量和壓縮率。對(duì)于每個(gè)起落架（前起落架、左起落架和右起落架）、壓縮量和壓縮率的推導(dǎo)如圖7所示。

·計(jì)算垂向力。確定油壓桿壓縮量和壓縮率后，就可以計(jì)算出力的三個(gè)分量（每個(gè)分量在起落架的局部機(jī)架中）。第一個(gè)分量用于計(jì)算垂向力Fz。垂向力承載飛機(jī)的重量，并在滑行期間為乘客提供舒適感。計(jì)算每個(gè)起落架的垂向力需先確定每個(gè)起落架的阻尼力和彈簧力，原因是起落架由非線性質(zhì)量阻尼彈簧系統(tǒng)表示。彈簧力可從參考文獻(xiàn)[18,19]中獲得。與用于計(jì)算靜態(tài)EPR的發(fā)動(dòng)機(jī)數(shù)據(jù)一樣，需對(duì)起落架數(shù)據(jù)（主起落架和前起落架的彈簧力和阻尼常數(shù)）進(jìn)行多項(xiàng)式擬合。

圖7. 油的壓縮量和壓縮率的推導(dǎo)。ΔSi和ΔS·i分別表示第i個(gè)起落架油的壓縮量和壓縮率（假定為三輪車模型，i = 1, 2, 3）。給定ΔSi和ΔS·i，圖中所示為求出垂向力的必要步驟。垂向力沿著支柱從輪胎作用到飛機(jī)的機(jī)身。θB和φB分別是機(jī)身的俯仰角和傾斜角。xi、yi、zi是第i個(gè)起落架到機(jī)身軸重心的距離。

·計(jì)算側(cè)向力。飛機(jī)滑行運(yùn)動(dòng)需要包括轉(zhuǎn)彎的一系列操作。轉(zhuǎn)彎運(yùn)動(dòng)是飛機(jī)不同部件的側(cè)向力作用的結(jié)果。計(jì)算起落架側(cè)向力的步驟如下：第一步需根據(jù)作用在輪胎上的垂向力確定輪胎的撓度。作用在輪胎上的垂向力與飛機(jī)的姿態(tài)有關(guān)。下一步需計(jì)算總側(cè)向力，總側(cè)向力是輪胎和運(yùn)動(dòng)方向之間的角度與輪胎撓度的乘積。最大總側(cè)向力是輪胎垂向力的60%。

·計(jì)算阻力。阻力產(chǎn)生制動(dòng)力，該制動(dòng)力由主起落架施加，并在滑行期間使飛機(jī)停下或減速。在起落架中，阻力與制動(dòng)力和摩擦力有關(guān)。制動(dòng)力與常數(shù)KB= 0.263、飛行器質(zhì)量和制動(dòng)踏板偏轉(zhuǎn)量有關(guān)。制動(dòng)力的最大值取決于飛機(jī)的質(zhì)量、滾動(dòng)摩擦（對(duì)于干滑行道，μB= 0.4）和最大制動(dòng)力常數(shù)（KBM= 0.834 + 4.167 ×μB）。摩擦力包含兩個(gè)分量[爆發(fā)力（bof）和恒定的滾動(dòng)摩擦力]。爆發(fā)力取決于飛機(jī)的地面速度（VG），其確定方法如下：

制動(dòng)踏板的偏轉(zhuǎn)已歸一化在0到1之間。

·將力轉(zhuǎn)換到機(jī)體坐標(biāo)系。以上在每個(gè)起落架中計(jì)算出的力都要轉(zhuǎn)換到機(jī)體坐標(biāo)系中，以便計(jì)算機(jī)身的力和力矩。參考文獻(xiàn)[19]所示的轉(zhuǎn)換方程式如下（僅對(duì)小角度有效）：

式中，θB和φB分別是機(jī)身的俯仰角和傾斜角；下角標(biāo)1、2和3分別代表三個(gè)輪胎。Fμi、FGi和FSi分別是第i個(gè)輪胎的輪胎阻力、垂直油壓桿支撐力和輪胎側(cè)向力；δ是前輪轉(zhuǎn)向角。通過三個(gè)起落架的力可以得到起落架的總力，如下所示：

力矩方程式如下：

式中，XL和YL是飛機(jī)重心到完全伸展開的起落架末端的距離。hB= 17 + ΔSi，是飛機(jī)重心到輪胎與跑道接觸所產(chǎn)生的法向力、側(cè)向力和阻力的垂直距離。

2.3.4. 重力模型

有兩種方法可以將重力產(chǎn)生的力加入飛機(jī)模型：第一種方法是利用平移加速函數(shù)[請(qǐng)參見式（5）]；第二種方法明確說明重力，然后將它們計(jì)入作用于飛機(jī)機(jī)身中的總力和力矩。我們使用后一種方法對(duì)這架飛機(jī)進(jìn)行建模。假定重力直接作用在飛機(jī)的重心上，表明沒有產(chǎn)生力矩。機(jī)體坐標(biāo)系中沿x軸（Fgx）、y軸（Fgy）和z軸（Fgz）方向的力由下式給出：

式中，GSA是地面傾斜角度。一個(gè)重要的問題是如何將滑行道的地面坡度加入重力模型。如果飛機(jī)正在下坡行駛（正斜角），則氣體將沿飛機(jī)的x軸產(chǎn)生正向力。此時(shí)地面傾斜角已忽略不計(jì)。

3. 飛機(jī)導(dǎo)航優(yōu)化

導(dǎo)航系統(tǒng)由圖8所示的方框組成，該系統(tǒng)允許飛機(jī)以最佳的方式在地面上移動(dòng)。在每個(gè)機(jī)場(chǎng)中，調(diào)度程序通常都會(huì)提供航路點(diǎn)列表以及到達(dá)每個(gè)航路點(diǎn)的時(shí)間期限。這些航路點(diǎn)通常由最優(yōu)調(diào)度算法[如具有時(shí)間窗的k條最快路徑（k-QPPTW）算法[5] ]或機(jī)場(chǎng)多目標(biāo)A*（AMOA*）算法[26]決定，且飛機(jī)需要以最佳方式遵循該調(diào)度。導(dǎo)航系統(tǒng)由一個(gè)外環(huán)控制器（由速度和航向算法表示）組成，用來(lái)確定參考點(diǎn)。與現(xiàn)有的開環(huán)方法[27-29]相比，這些參考點(diǎn)根據(jù)到下一個(gè)航路點(diǎn)的距離和時(shí)間限制進(jìn)行在線計(jì)算。內(nèi)環(huán)控制系統(tǒng)改變控制輸入端，使飛機(jī)遵循參考點(diǎn)。當(dāng)沒有自動(dòng)駕駛時(shí)，內(nèi)環(huán)系統(tǒng)控制器通常由飛行員替代。內(nèi)環(huán)系統(tǒng)控制器可以引導(dǎo)飛行員改變控制器輸入，從而完成航路點(diǎn)之間的最佳滑行（使用最小數(shù)量的燃料）。因此，本部分的主要目標(biāo)是開發(fā)一個(gè)最佳的控制器，以便在滿足計(jì)劃的同時(shí)使飛機(jī)（無(wú)論是通過完全自動(dòng)化實(shí)現(xiàn)還是作為決策支持）通過滑行道。這里使用的方法是通過調(diào)整PID控制器的參數(shù)，使所設(shè)定的目標(biāo)最小化，比如在遵循計(jì)劃的同時(shí)盡可能減少所消耗的燃料。

3.1. 外環(huán)控制系統(tǒng)

外環(huán)控制系統(tǒng)決定了遵循指定計(jì)劃時(shí)的航向和速度。本研究中未對(duì)外環(huán)控制系統(tǒng)進(jìn)行優(yōu)化。以下例子給出了確定參考航向和速度的過程。

示例1：假設(shè)調(diào)度程序向外環(huán)控制器提供以下信息，即在50 s內(nèi)沿著500 m的距離直線移動(dòng)。

外環(huán)控制系統(tǒng)將使用以下步驟確定速度和航向參考。

圖8. 導(dǎo)航系統(tǒng)的方框圖。

（1）確定航行段的類型。第一步是確定到達(dá)指定航路點(diǎn)的移動(dòng)是轉(zhuǎn)彎還是僅為直線段。在步驟2中計(jì)算出直線段的距離。

（2）計(jì)算到下一個(gè)航路點(diǎn)的距離（D）。

式中，（xp,yp）是航路點(diǎn)的坐標(biāo)；（xc,yc）是當(dāng)前飛機(jī)位置的坐標(biāo)。

（3）假設(shè)飛機(jī)繼續(xù)以當(dāng)前的速度移動(dòng)，計(jì)算預(yù)測(cè)的距離。預(yù)測(cè)距離（Dpred）如下：

式中，tr是到達(dá)航路點(diǎn)的剩余時(shí)間；Sc是當(dāng)前的速度。將從到航路點(diǎn)的距離減去預(yù)測(cè)距離可以得到距離誤差變量（E）。

（4）計(jì)算參考值。

式中，Sref是內(nèi)部控制器的參考值；γ是確定參考速度的變化速率的常數(shù)。γ可以進(jìn)一步被優(yōu)化，但當(dāng)它未被優(yōu)化時(shí)，經(jīng)反復(fù)試驗(yàn)可得γ= 10時(shí)的工作效果良好。

（5）在某些情況下，在末端指定速度約束（例如，當(dāng)從直線段移動(dòng)到轉(zhuǎn)彎段時(shí)，將在該轉(zhuǎn)彎段處使用恒定速度）。在這種情況下，利用以下算法可以計(jì)算參考速度。

（i）根據(jù)當(dāng)前速度（Sc）得到最大減速度（1 m·s-2）和最大加速度（1 m·s-2），計(jì)算飛機(jī)從當(dāng)前速度達(dá)到要求速度所需的時(shí)間（td）。這一時(shí)間可以通過方程td=ΔS/1求得，其中ΔS是Sc和航路點(diǎn)指定速度（Sf）之間的差的絕對(duì)值。

（ii）給定td并假定加速度或減速度恒定，計(jì)算該加速度或減速度滿足航路點(diǎn)速度約束所需的距離（Dd）。計(jì)算方法如下：

（iii）計(jì)算緩沖區(qū)距離（Db），以確定在滿足航路點(diǎn)的速度約束時(shí)行駛的總距離。計(jì)算方法如下：

上述方程確定了通過加速或減速滿足航路點(diǎn)的速度約束所需的距離。

（iv）如果抵達(dá)航路點(diǎn)所需的剩余時(shí)間大于td，則遵循正?？刂瞥绦騕如步驟（i）～（iii）中的定義]，用D-Db代替D，用t-td代替t。

（v）如果剩余時(shí)間小于或等于td，則設(shè)置Sref=Sf。

上述提出的導(dǎo)航算法已經(jīng)在以下幾種情況下進(jìn)行了測(cè)試，其中包括在50 s內(nèi)從原點(diǎn)移動(dòng)到(500 m, 0 m)處，且初始速度為5 m·s-1，最終速度為5 m·s-1。在上述情況下，該算法的測(cè)試結(jié)果如圖9所示。

由圖9分析可得，導(dǎo)航系統(tǒng)能夠在規(guī)定的截止時(shí)間前到達(dá)指定的航路點(diǎn)。對(duì)于轉(zhuǎn)彎段，同樣使用上述導(dǎo)航算法，但以不同的方式計(jì)算距離。具體地說，假設(shè)轉(zhuǎn)彎速度恒定（每秒4°），然后計(jì)算完成轉(zhuǎn)彎所需的時(shí)間和開始轉(zhuǎn)彎的時(shí)間。

3.2. 內(nèi)環(huán)控制系統(tǒng)

到目前為止，討論只涉及了如何確定參考（即航向和速度）。內(nèi)部控制回路負(fù)責(zé)改變飛機(jī)輸入（機(jī)舵、節(jié)流閥和制動(dòng)器），以使飛機(jī)遵循參考點(diǎn)。內(nèi)環(huán)控制由三套PID控制器組成，每個(gè)PID控制器分別操縱節(jié)流閥、制動(dòng)器和機(jī)舵。對(duì)于制動(dòng)器，只使用了比例控制器。為了避免同時(shí)使用制動(dòng)器和節(jié)流閥的情況，控制塊設(shè)定飛機(jī)僅在速度大于參考速度時(shí)使用制動(dòng)器，而在當(dāng)前速度小于參考速度時(shí)使用節(jié)流閥。PID控制系統(tǒng)的配置如圖10所示。

圖9. （a）在規(guī)定最終速度的情況下，飛機(jī)的參考速度和實(shí)際速度（5 m·s-1）；（b）隨時(shí)間變化的距離。

圖10. 在內(nèi)部控制回路中使用的PID控制系統(tǒng)的配置。P：比例，PI：比例-積分。

PID的參數(shù)需根據(jù)機(jī)場(chǎng)滑行的飛機(jī)情況進(jìn)行優(yōu)化。前文討論的示例1是一個(gè)典型的情況。在MATLAB（MathWorks，美國(guó)）中使用遺傳算法（GA）進(jìn)行優(yōu)化。性能標(biāo)準(zhǔn)取決于執(zhí)行指定任務(wù)時(shí)所消耗的燃料。綜上所述，優(yōu)化算法所使用的目標(biāo)函數(shù)采用如下方程式表示：

式中，J是目標(biāo)；θPID表示PID增益；T表示完成由調(diào)度程序確定的計(jì)劃所需的時(shí)間；而f(throttle)是根據(jù)圖5得出的燃料消耗模型。如式（22）所示，節(jié)流閥大小取決于PID增益。本研究中的GA將種群大小設(shè)置為1000。平均而言，這些算法在第90萬(wàn)代時(shí)近似收斂。上述方法的一個(gè)缺點(diǎn)是，PID增益對(duì)于另一種方案而言可能不是最佳的。為了找到新情況的最優(yōu)增益，優(yōu)化算法需要重新運(yùn)行，因此當(dāng)存在大量不同情況時(shí)，可能會(huì)產(chǎn)生顯著的計(jì)算成本。以英國(guó)曼徹斯特機(jī)場(chǎng)作為案例研究，調(diào)度程序包括大約9萬(wàn)個(gè)場(chǎng)景[圖11（a）]。顯然，獲得9萬(wàn)套最優(yōu)的PID增益是不可行的。為了解決這個(gè)問題，本文提出對(duì)情況列表進(jìn)行聚類分析。假設(shè)集群的中心代表集群成員。優(yōu)化只在表示集群中心的情況上進(jìn)行，并且在整個(gè)集群成員中利用集群中心的PID增益。曼徹斯特機(jī)場(chǎng)案例的集群分析結(jié)果如圖11（b）所示。

圖11.（a）場(chǎng)景的散射圖；（b）集群分析場(chǎng)景，其中用紅點(diǎn)表示的6個(gè)集群是各集群的中心。

表2包括了在10個(gè)隨機(jī)選擇的情況中與其他兩種方法相比的結(jié)果。第一個(gè)是從6個(gè)控制器中隨機(jī)選擇的一個(gè)控制器，第二個(gè)是為所選擇的所有情況選擇一個(gè)固定的控制器。

從表2中可以看到，使用聚類策略比對(duì)所有場(chǎng)景使用一個(gè)控制器具有更好的性能。對(duì)于控制策略違反調(diào)度程序或飛機(jī)限制的情況，應(yīng)包括懲罰項(xiàng)（消耗2000 kg燃料）。1000個(gè)隨機(jī)選擇方案的三種不同控制策略所消耗的總?cè)剂先绫?所示?？梢钥闯觯骄?，基于聚類的方法始終比使用固定控制器的策略性能更好。

表2 三種控制策略的結(jié)果的比較

表3 隨機(jī)選擇的1000個(gè)情況的總?cè)剂舷牧俊？刂破?在固定策略方法中提供了所有控制器中的最佳性能

4. 結(jié)果

本文在兩種情況下對(duì)提出的方法進(jìn)行了測(cè)試，并在以下小節(jié)中進(jìn)行討論。

4.1. 沿著矩形路徑滑行

第一種情況是沿著矩形方向進(jìn)行滑行，如圖12所示。

選擇飛機(jī)要遵循的路線圖是為了考慮飛機(jī)在滑行期間的限制，如最大加速度、減速度和轉(zhuǎn)彎速度的限制。該路線圖包括兩種滑行道。第一種是沿著矩形路徑滑行，并假設(shè)沒有地面坡度；第二種是沿著矩形路徑滑行，并假設(shè)地面坡度為2°。坡度因素很重要，因?yàn)轱w機(jī)在滑行過程中會(huì)經(jīng)過地面坡度不同的滑行道。飛機(jī)的起點(diǎn)為矩形底部直線的左端。從該點(diǎn)開始，路徑被分為8個(gè)不同的部分。路徑的各段測(cè)量位置軌跡如圖13所示。模擬數(shù)據(jù)段的截止時(shí)間見表4。

表4 模擬數(shù)據(jù)段的截止時(shí)間

滑行的結(jié)果如圖14所示，由圖可知飛機(jī)能夠在遵守時(shí)間計(jì)劃的同時(shí)遵循指定的路線。

4.2. 曼徹斯特國(guó)際機(jī)場(chǎng)滑行

在英國(guó)的一項(xiàng)真實(shí)滑行研究中對(duì)本文提出的算法進(jìn)行測(cè)試，該研究中的路線為從英國(guó)曼徹斯特國(guó)際機(jī)場(chǎng)登機(jī)口滑行到跑道等待點(diǎn)（準(zhǔn)備起飛）。機(jī)場(chǎng)滑行道布局和調(diào)度程序的優(yōu)化路線如圖15所示。航路點(diǎn)的坐標(biāo)見表 5。

表5 從調(diào)度程序算法輸出的航路點(diǎn)的地理坐標(biāo)

如圖16所示，在曼徹斯特國(guó)際機(jī)場(chǎng)滑行時(shí)，與傳統(tǒng)的Ziegler-Nichols PID方法相比，文章提出的方法可以節(jié)省約11%的燃料。

圖12. 用于測(cè)試的矩形圖。

圖13. 矩形路線試驗(yàn)中飛機(jī)的地理坐標(biāo)。（a）無(wú)坡度滑行；（b）以2°坡度滑行。

圖14. 通過坡度為0°的矩形段滑行的輸入。（a）節(jié)流閥；（b）制動(dòng)器；（c）機(jī)舵。

圖15. 曼徹斯特國(guó)際機(jī)場(chǎng)的滑行道和跑道布局。優(yōu)化的滑行路線以紅色線突出顯示。

圖16. 在曼徹斯特國(guó)際機(jī)場(chǎng)的航路點(diǎn)導(dǎo)航中使用優(yōu)化的PID的燃料消耗與使用Ziegler-Nichols方法的PID的燃料消耗的對(duì)比。（a）通過曼徹斯特國(guó)際機(jī)場(chǎng)滑行時(shí)所消耗的總耗油量；（b）通過曼徹斯特國(guó)際機(jī)場(chǎng)滑行時(shí)的燃料消耗率。

5. 結(jié)論

本文提出了一種針對(duì)機(jī)場(chǎng)飛機(jī)在實(shí)際滑行過程中的最佳導(dǎo)航方法。雖然該算法是基于波音747-100飛機(jī)模型進(jìn)行開發(fā)和測(cè)試的，但所提出的方法是通用的。該方法包括使用GA優(yōu)化PID增益，從而使飛機(jī)以最優(yōu)的方式進(jìn)行滑行。這種方法可以用于全自動(dòng)滑行的飛機(jī)，也可以作為飛行員的決策支持，即在飛行員以最佳的方式執(zhí)行操作時(shí)，向其提供能夠滿足調(diào)度程序的時(shí)間限制的最佳控制輸入。在將其用于解決實(shí)際滑行問題時(shí)，與傳統(tǒng)的Ziegler-Nichols PID控制器相比，該方法可以節(jié)省高達(dá)11%的燃料。本文的研究結(jié)果是在世界范圍內(nèi)實(shí)現(xiàn)飛機(jī)從門到門運(yùn)行中的一個(gè)重要飛躍。這種方法不僅有效地提高了飛機(jī)的操作效率，而且還有助于減少飛機(jī)的燃料消耗（燃料消耗是造成溫室氣體排放的重要因素）。未來(lái)的研究需量化飛機(jī)地面運(yùn)動(dòng)中許多不確定性的影響，如飛行員的操作和飛機(jī)的不同配置。

Acknowledgements

This work was funded by the UK Engineering and Physical Sciences Research Council (EP/N029496/1, EP/N029496/2, EP/N029356/1, EP/N029577/1, and EP/N029577/2).

Compliance with ethics guidelines

Olusayo Obajemu, Mahdi Mahfouf, Lohithaksha M.Maiyar, Abrar Al-Hindi, Michal Weiszer, and Jun Chen declare that they have no conflict of interest or financial conflicts to disclose.

Appendix A. Supplementary data

Supplementary data to this article can be found online at https://doi.org/10.1016/j.eng.2021.01.009.