基于改進BBO優化的客機空速控制仿真方法

耿 宏，操 正+，郝 磊

(1.中國民航大學 電子信息與自動化學院，天津 300300；2.中國民航大學 基礎實驗中心，天津 300300)

0 引 言

客機空速控制仿真是其飛行仿真的重要組成部分，良好的仿真效果可以準確地描述實際的空速控制特性。目前國內外關于空速控制的仿真多是從設計的層面，按照國家標準或飛行品質標準進行的[1-3]，而對于指定機型的空速控制仿真，則具有較強的針對性，如在某機型的飛行模擬設備中，其空速控制系統就需要逼近該機型實際表現出來的控制特性。一般而言，這些控制特性信息收集在該機型的飛行數據包中，但在實際工程中難以獲取該飛行數據包[4]，因此需要尋求一種易于實現的方法用于仿真指定機型的空速控制律。

在一些工業領域中，系統的歷史運行數據常用于評價或仿真該系統的實際運行過程，或用于系統模型參數的辨識過程。王印松等在其關于火電機組負荷控制的研究中利用系統的歷史數據作為參考，以評價實際系統的輸出特性[5]；代維在其研究中利用系統的歷史數據驅動型腔銑削中主軸功率的仿真和進給速度的優化[6]；Sen.S等則利用金屬塑性行為的歷史數據來估計金屬塑性參數[7]。然而在民航客機上，機載QAR(quick access recorder)能夠完整地記錄飛機在飛行過程中的狀態數據，該數據被廣泛應用于飛行品質分析或模型參數辨識等研究中[8,9]。本文基于上述在研究中融合系統歷史數據的思想，結合具有收斂速度快、精度較高且參數較少等優點的LBBO(Levy flight-biogeography based optimization)智能算法，將QAR數據應用到對指定機型空速控制律的仿真研究中，以改善最終的仿真效果。

1 客機空速控制仿真原理

現代客機通常有自動推力系統和升降舵兩種空速控制方式，根據其原理，建立指定機型的空速控制仿真系統，如圖1所示。Vr為目標空速，Vc為系統輸出的實際空速。當飛機處于升降舵控制空速時，PID(proportion-integration-differentiation)控制器輸出俯仰角指令，G(s)表示角位移控制的內回路；當飛機處于自動推力系統控制空速時，PID控制器輸出油門桿指令，G(s)表示油門伺服器及發動機，此時由自動駕駛儀H(s)控制垂直姿態。然而固定參數的控制器難以在不同的飛行條件下均具有良好的仿真效果，為此結合模糊控制策略[10]，并使模糊控制器的論域隨著控制需求按照一定的準則進行自適應性伸縮變化，以最終實現PID控制參數在線自調整。對于論域的伸縮機制，本文利用仿真系統的輸出與同時態下QAR數據之間的空速偏差和加速度偏差構造目標函數，通過LBBO算法優化論域的伸縮因子α和β，以在有限的迭代次數內獲取使得該偏差指標最小時的論域伸縮機制。最后將優化得到的伸縮因子代入變論域模糊控制器，以得到基于LBBO優化的變論域模糊PID控制系統，即為圖1中虛線框內的部分。

圖1 空速控制仿真原理

2 變論域模糊PID控制

2.1 相對變論域結構

變論域操作過程中，保持模糊論域不變，直接對量化因子和比例因子進行調整，以構成相對變論域模糊控制器。

在模糊控制器輸入端，空速跟蹤誤差值e的基本論域Ve(e)、 模糊論域Ue、實際量化因子ke(e)、 伸縮因子αe(e) 分別為

Ve(e)=[-αe(e)Me,αe(e)Me]

(1)

Ue=[-me,me]

(2)

(3)

αe(e)=1-φ1exp(-e2)

(4)

上式關于誤差e：Me為初始基本論域半區間寬度；me為模糊論域半區間寬度；ke0為初始量化因子；φ1為伸縮因子待定系數。

同理，空速跟蹤誤差變化率ec的基本論域Ve c(ec)、 模糊論域Ue c、實際量化因子ke c(ec)、 伸縮因子αe c(ec) 分別為

Ve c(ec)=[-αe c(ec)Me c,αe c(ec)Me c]

(5)

Ue c=[-me c,me c]

(6)

(7)

(8)

上式關于誤差變化率ec：Me c為初始基本論域半區間寬度；me c為模糊論域半區間寬度；ke c0為初始量化因子；φ2為伸縮因子待定系數。

對模糊控制器輸出端的控制參數修正值，其基本論域Vpara(e,ec)、 模糊論域Upara、實際量化因子kpara(e,ec) (即para分別表示比例環節、積分環節、微分環節)分別有

Vpara(e,ec)=[-βpara(e,ec)Mpara,βpara(e,ec)Mpara]

(9)

Upara=[-mpara,mpara]

(10)

(11)

上式關于各控制參數修正值：Mpara為初始基本論域半區間寬度；mpara為模糊論域半區間寬度；kpara0為初始比例因子；βpara(e,ec) 為對應的伸縮因子，且有

(12)

(13)

(14)

上式φ3、φ4、…、φ8均為伸縮因子待定系數。由于伸縮因子需要滿足單調性和協調性等要求，故將伸縮因子待定系數φ1、φ2、…、φ8的取值范圍均定為(0,1)。

然后，將模糊論域統一劃分為負大、負中、負小、零、正小、正中、正大7個模糊子集，對應地表示為NB、NM、NS、ZE、PS、PM、PB，且均采用三角形隸屬度函數。

為了便于表示，記xi為未經伸縮因子變換時的輸入端模糊量，αi為輸入端的伸縮因子，保持模糊論域不變，則輸入端的隸屬度函數為

(15)

輸出端模糊量是根據模糊規則推理而得到，其對應的伸縮因子將直接影響控制參數修正值的實際量大小，記yi為輸出端模糊量，故輸出端的隸屬度函數為

(16)

式中：an、bn、cn表示三角形隸屬度函數中的常數。

2.2 模糊規則

確定PID控制參數初始值Kp0、Ki0、Kd0之后，利用PID各參數的調節特點，制定出控制參數修正值的模糊量ΔKP、ΔKI、ΔKD分別與空速跟蹤誤差模糊量E及其變化率模糊量Ec之間的模糊規則[11]，見表1。

表1 比例/積分/微分修正值模糊規則

2.3 模糊推理與解模糊化

根據模糊規則建立模糊關系之后，在實際仿真過程中，模糊控制器獲取當前模糊量E及Ec，采用Mamdani推理得到控制參數修正值的模糊量，再利用重心法解模糊得到各參數修正值的實際量ΔKp、ΔKi、ΔKd，進而得到比例系數Kp、積分系數Ki、微分系數Kd的最終值為

(17)

3 LBBO優化伸縮因子

3.1 初始解空間與遷移模型

采用LBBO算法優化變論域過程中的伸縮因子，以相對調整模糊規則，促使仿真系統的輸出向同時態下的QAR數據擬合。本文通過對函數型伸縮因子中的待定系數進行尋優的方式來優化伸縮因子，即φ1、φ2、…、φ8構成了解的8個分向量，并設解的數目為N，則初始解空間Ф可用矩陣表示為

(18)

由于φ均是小于1的正數，故通過(0,1)區間內均勻分布的隨機數來生成初始解空間。

BBO是一種基于自然啟發式的智能算法，為了更貼近自然界棲息地的實際情況以提高尋優性能，選擇非線性三角形遷移模型描述遷移率與物種多樣性之間的關系[12]。設遷入率峰值為I，遷出率峰值為E，物種數目峰值即為解的數目N，則根據遷移模型，得到各解的遷入率λi和遷出率μi分別為

(19)

(20)

其中，i為每次遷移和變異操作之前，對所有解按照適應度值(habitat suitability index，HSI)從低到高重新排列后的解序號值。

3.2 混合遷移和Levy Flight變異

原始BBO算法具有較強的局部搜索能力，為提高其全局搜索能力，在遷移操作上，引入混合遷移機制，在變異操作中，結合具有高頻小步長與低頻大步長特點的Levy Flight游走方式[13]。

遷移操作時，各解根據其遷入率隨機決定接受遷入的分向量，再結合其它各解的遷出率，利用賭輪選擇機制確定遷出量并進行混合遷移，混合遷移算子定義為

φ′i(j)=ρφi(j)+(1-ρ)φh(j)

(21)

ρ=rand()

(22)

φi(j)、φh(j) 分別表示第i、h個解的第j維，rand()表示生成一個[0,1]區間內的隨機數。

變異操作時，首先利用Mantegna算法生成Levy分布[14]，即

(23)

式中：γ屬于常數，v和u均服從正態分布，有

γ∈[0.1,2]

(24)

v～N(0,1)

(25)

u～N(0,σu2)

(26)

(27)

其中，Γ(·)為標準的Gamma函數。

然后在解的各分向量上加上Levy步長以實現變異操作，從而有效減少尋優結果陷入局部最優的現象。在實際仿真過程中，根據Levy分布的特點以及φ的取值要求，需要對Levy步長進行等比例縮放與限幅，即在變異操作過程中有

φ′i(j)=Limited[φi(j)+Limited(·L(u,v))]

(28)

3.3 適應度函數

為提高對指定機型空速控制律的仿真度和有效性，使用經過野值處理、平滑處理及重新采樣等基本操作后的QAR數據，且要求在所選數據對應的飛行過程中，沒有嚴重影響飛行控制性能的故障。

由于在適應度函數的設計過程中需要結合同時態下的真實空速值，而QAR數據中記錄的空速指的是儀表空速，即為飛機儀表所測量的速度。然而在實際飛行過程中，由于受到空氣壓縮性的影響，儀表空速與真實空速之間存在一定的偏差。因此在使用QAR數據中的空速數據之前，需要對其進行修正。

在標準大氣壓下，當飛行過程中全壓和靜壓管所測得的數據精準性較高的時候，儀表空速可較高程度地反映真實空速。因此，可通過由靜壓測量誤差導致的速度指示誤差(ΔVP)得到修正空速(VD)，有

VD=VI+ΔVP

(29)

式中：VI即為儀表空速，ΔVP可以從飛機手冊上的修正曲線中獲取。進而得到在不同高度層下的真實空速VT

(30)

(31)

T和P分別表示當前高度層下的溫度和總壓，T0和P0分別表示海平面標準大氣壓下的溫度和總壓，且有K=1.4，a0=661.475 kn。

設Eerror為仿真系統輸出的空速值Vc與QAR數據修正空速值VT之間的偏差，ECerror為仿真系統輸出的加速度值Ac與QAR數據參考加速度值AT之間的偏差，由于本文中的加速度是為滿足算法需求而定義的物理量，即取單個采樣周期內的空速變化量，其不等同于QAR數據中記錄的加速度。

為了便于計算機計算，令目標函數如式(32)，式中dT即為采樣周期，Eerror(t)、ECerror(t) 分別為第t次采樣時的空速偏差和加速度偏差。采用LBBO算法優化的目的是使得目標函數值減小，然而LBBO算法本質是選擇HSI值最高的棲息地，故構造其適應度函數如式(34)

(32)

η∈(0,1)

(33)

(34)

綜上所述，利用LBBO算法優化伸縮因子的思路流程如圖2所示，其中G表示當前迭代的次數，Gmax為最大迭代次數，φi表示第i個解對應的一組伸縮因子待定系數，φbest用于存儲迭代過程中HSI值最高的一組伸縮因子待定系數，fi表示第i個解所得到的HSI值，fbest用于存儲歷史最高的HSI值。

圖2 伸縮因子優化流程

4 仿真驗證

4.1 JSBSim動力學模型

JSBSim是利用C++程序語言編寫的通用飛行動力學模型，其主要由接口變量、中間變量、計算模型3部分構成。接口變量主要由飛行實時參數集和飛行目標參數集構成，中間變量多指計算模型中的過程量，計算模型主要包含氣動單元(Aerodynamics)、推進單元(Propulsion)、飛行控制單元(flight control system，FCS)。仿真過程中分別通過機型文件(Model.xml)和發動機型號文件(Engine.xml)將Aerodynamics和Propulsion實例化為指定機型，FCS中可以利用JSBSim定義的組件搭建傳統的PID控制器，其能夠基本實現飛行高度控制、垂直速度控制及空速控制等。

4.2 仿真實驗

4.2.1 伸縮因子優化

本文以JSBSim飛行動力學模型為控制對象，聯合Matlab/Simulink工具，以仿真A320飛機在自動駕駛儀控制垂直姿態下降時，自動推力控制空速的模態為例。先后選取6段具有代表性的A320飛機QAR數據作為訓練對象，并分別進行實驗。為了便于分析實驗結果，在每次實驗過程中均由JSBSim本身的FCS實現垂直姿態的控制，并使得JSBSim中的仿真環境盡可能與QAR數據中的相似，主要包括飛行工作狀態即自動駕駛儀(AP)的工作模式、自動推力系統(AT)的工作模式、發動機狀態(ENG STATE)；還包括一些飛行參數即初始氣壓高度(ALT0)、垂直速度(Vs)、剩余機載燃油(FOB)、自然風(WindSpeed/Direction)等。此外，使QAR數據重新采樣周期與仿真系統采樣周期保持一致，即均為dT。由于這些初始設置在6段QAR數據中不盡相同，故此處不作羅列。將上述設置完成之后的JSBSim源碼文件編寫成S函數，并生成Simulink可直接調用的MEX文件。其中輸入參數為油門閥開度值(δ)；輸出參數包括實際空速(Vc)與JSBSim輸出數據采樣計數器(Timer)。此外在仿真過程中，由于推力桿的位置與推力的轉換等環節已直接實例化在JSBSim-Propulsion中，故此處G(s)為油門伺服器，并用時間常數為0.05 s的一階慣性環節表示。最后，再以隨機噪聲信號模擬實際飛行過程中的自然干擾。

模糊控制器參數設置見表2，其中m為各基本變量對應的模糊論域半區間寬度，k0為初始量化因子或比例因子。

表2 模糊控制器初始設置

關于LBBO算法的參數設置：取式(19)、式(20)中的I為1，E為1，N為30；取式(23)中的γ為1.5；另外取Gmax為100。在確保引入的Levy步長不影響原始算法收斂性的情況下，取式(28)中的值為0.0025，且經過測試，取式(32)中的η為0.76時仿真效果較佳。如圖3所示，為6段QAR數據分別作為參考模型輸出時，優化過程中HSI值的變化規律，其中縱軸為歸一化的HSI值(即取實驗過程中最大HSI值為1)，橫軸為迭代次數，最終得到的一組伸縮因子待定系數為：φ1=0.49；φ2=0.58；φ3=0.43；φ4=0.70；φ5=0.53；φ6=0.54；φ7=0.77；φ8=0.54。

圖3 HSI變化規律

4.2.2 驗證示例

第一組實驗：選擇一段區別于伸縮因子優化階段的飛行過程用于驗證，該過程中飛機從3040 ft開始，以-800 ft/min 的垂直速度下降時，自動推力控制空速從180 kts 降至140 kts。其中為了便于結果分析和實際應用，將空速逆轉換為儀表空速。仿真過程中，飛機的垂直姿態控制由JSBSim本身的FCS實現，并參照對應的QAR數據，將JSBSim初始設置為表3狀態，G(s)及干擾信號與伸縮因子優化過程中的相同。在同等仿真環境下使用JSBSim-FCS 中的傳統PID控制器(記作JS-PID)、固定論域的模糊PID控制器(記作FPID)、LBBO優化的變論域模糊PID控制器(記作LBVUFPID)進行對比性實驗。

為描述控制系統的動態特性和靜態特性，選取控制過程中的上升時間tr、最大偏差eσ(即達到控制目標之后的最大超調值)、調節時間ts(2%的誤差)作為仿真系統的控制性能指標，以與QAR數據所表現的控制特性作比較來說明仿真控制器的控制效果。已知QAR數據在該飛行過程中有：tr為73.90 s，eσ為0.92 kts，ts為83.60 s，同理分析對比實驗的對應指標。再分別計算實驗結果與QAR數據之間的均方根誤差(RMSE)和擬合度(R-square)，以作為系統的仿真效果指標，其定義為

表3 JSBSim初始設置

(35)

(36)

表4 第一組實驗指標

圖4 空速變化曲線

第二組實驗：另行選擇兩段不同航班在此模態下的典型飛行過程進行實驗驗證，該部分給出使用LBVUFPID控制器的仿真結果與實際飛行規律之間的比較，其中實驗環境的配置方法和上述一致。實驗結果如圖5(a)(從315 kts減速至275 kts)、圖5(b)(從250 kts減速至220 kts)。

圖5 空速變化曲線

分別計算兩次實驗中的控制性能指標和仿真效果指標，結果見表5。

表5 第二組實驗指標

從本組實驗結果中得知：在該模態下，對多個不同飛行環境下的降速控制，LBVUFPID控制器的輸出規律與QAR數據規律均具有較高的相似程度，從而進一步說明了該方法的有效性。

5 結束語

本文提出一種基于LBBO優化的變論域模糊PID控制方法，用于在無法明確指定機型飛行控制特性的情況下，仿真其空速控制律。結合被仿機型的QAR數據構造目標函數，利用LBBO算法間接優化模糊控制規則，促使仿真系統輸出規律擬合真實飛行數據所呈現的規律。從多角度的實驗結果得知，該方法更為準確地描述了A320自動推力系統在下降階段的空速控制特性，但在確定該方法目標函數中的加權系數η以及Levy步長縮放系數時，存在著一定的主觀性，在下一步研究中需要探索一套科學的方法用于確定這些系數。