基于高斯過程回歸的連續(xù)式風(fēng)洞馬赫數(shù)控制

高 赫,劉學(xué)軍,郭 晉,呂宏強(qiáng)

(1.南京航空航天大學(xué) 計算機(jī)科學(xué)與技術(shù)學(xué)院,模式分析與機(jī)器智能工業(yè)和信息化部重點(diǎn)實驗室,南京 211106;2.空氣動力學(xué)國家重點(diǎn)實驗室,綿陽 621000;3.軟件新技術(shù)與產(chǎn)業(yè)化協(xié)同創(chuàng)新中心,南京 210023;4.中國航空工業(yè)空氣動力研究院,沈陽 110034;5.南京航空航天大學(xué) 航空學(xué)院,南京 210016)

0 引 言

風(fēng)洞是能人工產(chǎn)生和控制氣流，以模擬飛行器或物體周圍氣體流動，并可量度氣流對物體的作用以及觀察物理現(xiàn)象的一種管道實驗設(shè)備，它是進(jìn)行空氣動力實驗最常用、最有效的工具之一[1]。在進(jìn)行風(fēng)洞實驗時，保持風(fēng)洞實驗段馬赫數(shù)的穩(wěn)定，以測量飛行器在不同迎角下的氣動性能，對于研制高性能飛行器具有重要意義。為了獲得準(zhǔn)確的實驗數(shù)據(jù)，階梯變迎角實驗?zāi)Ｊ奖粡V泛應(yīng)用在風(fēng)洞實驗中。階梯變迎角模式基于迎角的變化范圍間歇地改變迎角以保持馬赫數(shù)的穩(wěn)定，但這種方式只能獲得少量特定迎角處的實驗數(shù)據(jù)，實驗效率不高，而且在每一個迎角處都要等待馬赫數(shù)穩(wěn)定后再測量相應(yīng)的氣動性能，這樣就不可避免地延長了實驗時間，造成較高的能源消耗，對于風(fēng)洞這種高耗能設(shè)備來說是不經(jīng)濟(jì)的。連續(xù)變迎角實驗?zāi)Ｊ骄哂懈叩膶嶒炐屎洼^低的運(yùn)行成本，并且可以獲得更密集的實驗數(shù)據(jù)，這些優(yōu)點(diǎn)促使風(fēng)洞實驗方式開始從階梯變迎角向連續(xù)變迎角方式轉(zhuǎn)變。但與此同時,連續(xù)變迎角實驗?zāi)Ｊ降目焖贂r變特性也給風(fēng)洞馬赫數(shù)的控制帶來了新的挑戰(zhàn)，需要研究新的控制算法與策略。

傳統(tǒng)的PID(Proportion Integral Differential)控制算法及其改進(jìn)形式在階梯變迎角實驗過程中可以很好地滿足風(fēng)洞馬赫數(shù)的控制要求。文獻(xiàn)[2]利用神經(jīng)網(wǎng)絡(luò)對PID控制參數(shù)進(jìn)行在線實時整定，通過仿真實驗證明了該方法可以使馬赫數(shù)控制更加穩(wěn)定、迅速。文獻(xiàn)[3]針對跨超聲速風(fēng)洞馬赫數(shù)范圍較寬的特點(diǎn)，采用了分段控制的思想，在不同的馬赫數(shù)范圍內(nèi)，將模糊控制跟PID控制結(jié)合，利用模糊控制的快速性與PID控制的穩(wěn)定性實現(xiàn)了整個馬赫數(shù)范圍內(nèi)的快速穩(wěn)定控制。由于風(fēng)洞的主體是一個環(huán)形管道，控制量(如壓縮機(jī)轉(zhuǎn)速)與被控量(實驗段馬赫數(shù))在實際物理位置上存在一定的距離，導(dǎo)致控制量的變化不能直接作用于被控量，所以在進(jìn)行連續(xù)變迎角實驗時存在一定的延遲，當(dāng)迎角變化的速度超過0.1°/s時，僅僅使用PID控制或其改進(jìn)算法無法滿足馬赫數(shù)控制精度達(dá)到±0.001的要求。

針對PID控制算法存在的缺陷，近年來出現(xiàn)了很多基于預(yù)測的控制方法，這些方法統(tǒng)稱為模型預(yù)測控制(Model Predictive Control，MPC)算法[4]。許多研究者已經(jīng)通過仿真實驗證明，使用模型預(yù)測控制可以很好地對連續(xù)變迎角實驗?zāi)Ｊ较碌鸟R赫數(shù)進(jìn)行控制。文獻(xiàn)[5]使用統(tǒng)一模型預(yù)測控制對阿姆斯特丹航空航天實驗室擁有的高速風(fēng)洞進(jìn)行馬赫數(shù)控制，實驗結(jié)果表明該控制方法與傳統(tǒng)PID控制算法相比，抗迎角擾動能力提高4倍并且馬赫數(shù)的控制性能提高30%～60%。文獻(xiàn)[6]出于馬赫數(shù)高精度控制的考慮，選取穩(wěn)定段總壓和實驗段靜壓作為被控量，設(shè)計了基于迎角補(bǔ)償?shù)哪Ｐ皖A(yù)測控制器，使馬赫數(shù)的控制精度達(dá)到了大飛機(jī)研制要求。但是模型預(yù)測控制方法需要在線求解約束優(yōu)化問題，計算量大，只適用于具有高性能計算機(jī)的環(huán)境[7]。

中國航空工業(yè)空氣動力研究院的0.6 m×0.6 m連續(xù)式跨聲速風(fēng)洞將壓縮機(jī)轉(zhuǎn)速作為控制馬赫數(shù)穩(wěn)定的唯一驅(qū)動方式。為了解決PID控制方法在連續(xù)變迎角實驗?zāi)Ｊ较麓嬖诘娜毕?，同時盡可能節(jié)約能源，其目前主要采用前饋控制與PID反饋控制結(jié)合的策略。前饋部分的設(shè)計需要先通過固定轉(zhuǎn)速，連續(xù)變化迎角，得到迎角與實際馬赫數(shù)的關(guān)系曲線；然后再通過固定迎角，增加轉(zhuǎn)速，得到轉(zhuǎn)速與實際馬赫數(shù)的關(guān)系曲線；最后通過這兩條曲線得到前饋系數(shù)，即轉(zhuǎn)速變化量與迎角變化量的比值，這樣就可以根據(jù)迎角變化量提前對轉(zhuǎn)速進(jìn)行補(bǔ)償。但是這種方式存在兩個明顯缺點(diǎn)：1)為了得到前饋系數(shù)需要進(jìn)行兩次實驗，增加了實驗次數(shù)；2)不同迎角下，轉(zhuǎn)速與實際馬赫數(shù)的關(guān)系曲線是存在一定差別的，導(dǎo)致前饋系數(shù)計算不準(zhǔn)確。

針對這兩個缺點(diǎn)，本文采用高斯過程回歸(Gaussian Process Regression，GPR)直接對不同馬赫數(shù)下的實驗數(shù)據(jù)進(jìn)行建模，在其他條件一定的情況下，將迎角、實際馬赫數(shù)、實驗段截面積作為模型輸入，壓縮機(jī)轉(zhuǎn)速作為模型輸出，直接預(yù)測給定馬赫數(shù)前提下不同迎角對應(yīng)的壓縮機(jī)轉(zhuǎn)速。根據(jù)預(yù)測得到的轉(zhuǎn)速值，提前對壓縮機(jī)轉(zhuǎn)速進(jìn)行調(diào)節(jié)，可以減輕PID控制的調(diào)節(jié)壓力，從而提高連續(xù)變迎角實驗方式下的馬赫數(shù)控制精度。這種直接從實驗數(shù)據(jù)中學(xué)習(xí)規(guī)律的方法減少了實驗次數(shù)，并且避免了近似計算造成的誤差。另外，GPR輸出的預(yù)測值置信度可以指導(dǎo)新增實驗數(shù)據(jù)的獲取，以進(jìn)一步保證馬赫數(shù)控制精度。

1 高斯過程回歸(GPR)

GPR是基于貝葉斯理論和統(tǒng)計學(xué)習(xí)理論發(fā)展起來的一種全新的機(jī)器學(xué)習(xí)方法，適于處理小樣本、非線性等復(fù)雜問題[8]。GPR不需要顯式的指定目標(biāo)函數(shù)的具體形式，只需假設(shè)其服從某個指定均值函數(shù)和協(xié)方差函數(shù)的高斯過程，目標(biāo)函數(shù)的后驗分布是通過擬合訓(xùn)練數(shù)據(jù)自動學(xué)習(xí)得到的，所以它屬于非參數(shù)模型[9-10]。相較于其他常用的參數(shù)模型，比如神經(jīng)網(wǎng)絡(luò)、支持向量回歸、多項式回歸等，GPR具有容易訓(xùn)練、超參數(shù)自適應(yīng)獲取以及預(yù)測結(jié)果具有概率意義等優(yōu)點(diǎn)，并且已經(jīng)在許多領(lǐng)域得到了成功應(yīng)用[11-13]。

1.1 預(yù)測

高斯過程是一個隨機(jī)變量的集合，并且其中任意有限數(shù)量隨機(jī)變量的聯(lián)合分布為多元高斯分布。從函數(shù)空間的角度出發(fā)，高斯過程可以看作定義在函數(shù)f(x)上的一個分布，它的性質(zhì)完全由均值函數(shù)和協(xié)方差函數(shù)確定:

m(x)=E[f(x)]

(1)

k(x,x′)=E[(f(x)-m(x))(f(x′)-m(x′))]

(2)

f(x)～GP(m(x),k(x,x′))

(3)

式中，x，x′∈Rd,為d維輸入向量，m(x)為均值函數(shù)，k(x,x′)為協(xié)方差函數(shù)。

假設(shè)訓(xùn)練集為{(xi,yi)|i=1,…,n}，針對回歸問題，考慮如下模型：

yi=f(xi)+ε

(4)

(5)

式中，X為特征矩陣，其中每一行表示一個輸入向量；K(X,X)=Kn=(kij)為n×n階對稱正定協(xié)方差矩陣，矩陣元素kij=k(xi,xj)，使用xi與xj之間的距離來度量yi和yj之間的相關(guān)性；In為n階單位矩陣。

對于給定的測試樣本x*，其對應(yīng)的預(yù)測值f*與觀測值y的聯(lián)合先驗分布為：

(6)

其中，K(x*,X)=K(X,x*)T為測試樣本與訓(xùn)練集之間的1×n階的協(xié)方差矩陣；K(x*,x*)為測試樣本之間的協(xié)方差矩陣。

根據(jù)式(6)，并且通過一定的矩陣運(yùn)算，可以得到預(yù)測值f*的條件概率分布為：

(7)

其中，

(8)

cov(f*)=K(x*,x*)-

K(x*,X)×[K(X,X)+In]-1K(X,x*)

(9)

最終，我們將式(8)作為測試樣本的預(yù)測值，式(9)表示預(yù)測的方差，用來衡量該預(yù)測的不確定性。

1.2 訓(xùn)練

GPR的訓(xùn)練其實是選擇合適的協(xié)方差函數(shù)以及確定其最優(yōu)超參數(shù)的過程[8]。不同的協(xié)方差函數(shù)確定了在該高斯過程先驗下目標(biāo)函數(shù)可能具有的性質(zhì)。比如，周期協(xié)方差函數(shù)表示目標(biāo)函數(shù)具有周期性；平方指數(shù)協(xié)方差函數(shù)表示目標(biāo)函數(shù)具有無窮階導(dǎo)數(shù)，即處處光滑；線性協(xié)方差函數(shù)表示目標(biāo)函數(shù)是線性的[14]。在實際應(yīng)用中，使用最廣泛的為平方指數(shù)協(xié)方差函數(shù)，具體形式如下：

(10)

首先根據(jù)多元高斯分布的邊際屬性可得觀測值y的邊緣概率分布：

(11)

式(11)一般稱作邊際似然。通過式(11)可以得到訓(xùn)練集的負(fù)對數(shù)邊際似然函數(shù)為：

L(θ)=-lnp(y|X)

(12)

θ*=argminθL(θ)

(13)

2 實驗設(shè)計

2.1 數(shù)據(jù)預(yù)處理

使用J7標(biāo)模在連續(xù)變迎角實驗?zāi)Ｊ较碌玫降臄?shù)據(jù)進(jìn)行建模，其中迎角變化范圍為-5°～11°。在實際實驗過程中，由于迎角回零操作或者其他一些復(fù)雜因素的影響可能會產(chǎn)生一些異常數(shù)據(jù)，這些異常數(shù)據(jù)會嚴(yán)重影響建模的準(zhǔn)確性，所以在建模前需要將這些數(shù)據(jù)進(jìn)行刪除。刪除前與刪除后的數(shù)據(jù)分別如圖1、圖2所示，圖中的虛線表示馬赫數(shù)±0.001誤差帶。

圖1 包含異常點(diǎn)的數(shù)據(jù)Fig.1 Data with outliers

原始實驗數(shù)據(jù)中，為了檢驗連續(xù)變迎角實驗方式下角速度對馬赫數(shù)控制精度的影響，在某些馬赫數(shù)條件下可能測量了兩組實驗數(shù)據(jù)，其迎角變化速度分別為0.1°/s和0.2°/s，但是由于角速度不影響迎角、實時馬赫數(shù)、實驗段截面積和壓縮機(jī)轉(zhuǎn)速之間的關(guān)系，同時也為了保證不同馬赫數(shù)下的數(shù)據(jù)量盡量均衡，所以對于這種情況只取角速度為0.1°/s的數(shù)據(jù)進(jìn)行建模。在進(jìn)行實際風(fēng)洞實驗時，不同范圍的馬赫數(shù)條件下，實驗段截面積也可能會做相應(yīng)的調(diào)整。本文使用的數(shù)據(jù)中，馬赫數(shù)0.85以上(包括馬赫數(shù)0.85)條件下的截面積是馬赫數(shù)0.85以下的7/6倍，所以我們以馬赫數(shù)0.85為界，其以下的截面積設(shè)為1 m2，以上的截面積設(shè)為7/6 m2(由于最終輸入模型的數(shù)據(jù)是經(jīng)過標(biāo)準(zhǔn)化的，所以這樣設(shè)置是合理的)。經(jīng)過處理后的數(shù)據(jù)集的特征分布直方圖如圖3所示。

圖2 清洗后的數(shù)據(jù)Fig.2 Cleaned data

由圖3可知，不同特征之間的量級是不一樣的，而GPR主要采用基于梯度的優(yōu)化算法進(jìn)行訓(xùn)練。為了加快梯度下降的速度，并且使得算法可以更好地逼近最優(yōu)解，提高預(yù)測精度[15]，所以在將數(shù)據(jù)輸入模型之前，還需要對每一個特征進(jìn)行標(biāo)準(zhǔn)化處理，標(biāo)準(zhǔn)化公式如下所示：

(14)

圖3 特征分布直方圖Fig.3 Feature distribution histogram

2.2 初始超參數(shù)選擇

GPR模型中的超參數(shù)可以通過最大化訓(xùn)練數(shù)據(jù)的邊際似然函數(shù)獲得，但是一般情況下邊際似然函數(shù)都是關(guān)于超參數(shù)的非凸函數(shù)，在使用基于梯度的優(yōu)化算法進(jìn)行求解時，如果超參數(shù)的初值選取不當(dāng)，算法很可能會收斂到局部最優(yōu)解。為了解決該問題，采用最大邊際似然法進(jìn)行初始超參數(shù)的選擇，具體的步驟如下：

① 從特定的超參數(shù)先驗P(θ)中隨機(jī)選擇一個超參數(shù)θ0；

② 將θ0作為初始值，使用共軛梯度下降算法最小化負(fù)對數(shù)邊際似然函數(shù)，經(jīng)過m次迭代，得到優(yōu)化后的超參數(shù)θm；

③ 重復(fù)①、②M次，從M個初始超參數(shù)中選擇使負(fù)對數(shù)邊際似然達(dá)到最小的作為最終的初始超參數(shù)。

文獻(xiàn)[16]通過大量的實驗得出結(jié)論：盡管超參數(shù)先驗分布對協(xié)方差函數(shù)最終優(yōu)化得到的超參數(shù)取值有較大影響，但是對GPR模型最終的預(yù)測精度影響甚微。所以在實際應(yīng)用中選擇比較簡單、常見的概率分布作為超參數(shù)的先驗分布即可，比如標(biāo)準(zhǔn)高斯分布、合適范圍內(nèi)的均勻分布等。

本文最終選擇標(biāo)準(zhǔn)高斯分布作為超參數(shù)的先驗分布，從中隨機(jī)采樣20次，得到20組不同的初始超參數(shù)。對于每組超參數(shù)使用共軛梯度下降算法進(jìn)行500次迭代優(yōu)化，從中選取使負(fù)對數(shù)邊際似然最小的取值作為最終的初始超參數(shù)。

2.3 協(xié)方差函數(shù)的選擇

正如1.2節(jié)所述，協(xié)方差函數(shù)對GPR模型至關(guān)重要，因為它包含了我們對于目標(biāo)函數(shù)特性的假設(shè)。當(dāng)我們對目標(biāo)函數(shù)的先驗知識比較充足時，可以準(zhǔn)確地選擇合適的協(xié)方差函數(shù)來解釋數(shù)據(jù)中的結(jié)構(gòu)。但是通常情況下，我們并不清楚目標(biāo)函數(shù)的一些顯著特征。式(10)所示的平方指數(shù)協(xié)方差函數(shù)的超參數(shù)數(shù)量較少且具有可解釋性，而且已有理論證明，“當(dāng)訓(xùn)練數(shù)據(jù)足夠多時，平方指數(shù)協(xié)方差函數(shù)有能力學(xué)得任意的連續(xù)函數(shù)”[17]。所以當(dāng)我們對目標(biāo)函數(shù)的先驗知識不足時，為了不引入錯誤的先驗知識并且保證預(yù)測的準(zhǔn)確性，可以使用平方指數(shù)協(xié)方差函數(shù)來對數(shù)據(jù)進(jìn)行解釋。下文將使用平方指數(shù)協(xié)方差函數(shù)的GPR模型稱為標(biāo)準(zhǔn)高斯過程回歸(Standard Gaussian Process Regression，S-GPR)。

S-GPR假設(shè)目標(biāo)函數(shù)具有如下形式：

y=f(x1,x2,…,xD)

(15)

式中x1,x2,…,xD表示輸入向量的各個維度。最終的函數(shù)輸出y同時取決于輸入向量的所有維度，這使得S-GPR對于遠(yuǎn)離訓(xùn)練樣本的數(shù)據(jù)預(yù)測效果很差。

為了解決該問題，本文使用加性高斯過程回歸(Additive Gaussian Process Regression，Add-GPR)[18]來對風(fēng)洞實驗數(shù)據(jù)進(jìn)行建模。Add-GPR使用加性協(xié)方差函數(shù)，其可以看作廣義加性模型(Generalized Additive Model，GAM)[19]與S-GPR的推廣。具體的，為輸入向量的每一個維度分配一個基本的協(xié)方差函數(shù)，并且定義不同階的加性協(xié)方差函數(shù)以全面考慮特征之間的交互作用。不同階的協(xié)方差函數(shù)定義如下：

(16)

(17)

(18)

(19)

易知，當(dāng)只考慮特征之間的一階交互時，Add-GPR等價于GAM；如果只考慮特征之間的D階交互，并且將基本協(xié)方差函數(shù)設(shè)為平方指數(shù)協(xié)方差函數(shù)，則Add-GPR就退化為S-GPR。

由于平方指數(shù)協(xié)方差函數(shù)的良好性質(zhì)，所以本文將其作為迎角、實時馬赫數(shù)、實驗段截面積三個輸入維度的基本協(xié)方差函數(shù)，并且使用式(19)所示的完整加性協(xié)方差函數(shù)(此時D=3)作為GPR最終的協(xié)方差函數(shù)對實驗數(shù)據(jù)進(jìn)行建模。

3 實驗結(jié)果及分析

經(jīng)過數(shù)據(jù)清洗后，得到在馬赫數(shù)0.6、0.7、0.75、0.8、0.85、0.9、0.95、1.0、1.05、1.1條件下，迎角連續(xù)變化時的10組控制數(shù)據(jù)，總樣本個數(shù)為2578。其中輸入向量在特征空間的分布情況如圖4所示。

圖4 輸入向量散點(diǎn)圖Fig.4 Scatter plot of input vector

圖4中的每組數(shù)據(jù)的實際馬赫數(shù)都在相應(yīng)的馬赫數(shù)條件附近波動。以馬赫數(shù)0.75條件下的數(shù)據(jù)為例，實際馬赫數(shù)的波動情況如圖5所示。

圖5 實際馬赫數(shù)散點(diǎn)圖Fig.5 Scatter plot of actual Mach numbers

實驗將從兩個角度說明Add-GPR對風(fēng)洞實驗數(shù)據(jù)進(jìn)行建模的有效性,并且提出了基于高斯過程回歸的馬赫數(shù)控制策略。

3.1 隨機(jī)劃分?jǐn)?shù)據(jù)集

使用機(jī)器學(xué)習(xí)算法執(zhí)行預(yù)測任務(wù)有一個基本的假設(shè)，即訓(xùn)練數(shù)據(jù)與測試數(shù)據(jù)是獨(dú)立同分布的?；谶@個假設(shè)我們將2578個樣本隨機(jī)打亂，使用常用的機(jī)器學(xué)習(xí)模型分別進(jìn)行5折交叉驗證，最終得到預(yù)測結(jié)果的RMSE(Root Mean Square Error)如表1所示。其中，RMSE的計算公式為：

(20)

表1 預(yù)測結(jié)果對比Table 1 Comparison of prediction results

由表1可以看出，Add-GPR每一次交叉驗證的預(yù)測精度都要高于其他常用的回歸模型，同時GPR的預(yù)測效果也比其他模型的預(yù)測精度要好。這是因為神經(jīng)網(wǎng)絡(luò)、SVR和多項式回歸屬于參數(shù)化模型[9]，參數(shù)化模型需要對目標(biāo)函數(shù)的形式做有限參數(shù)的假設(shè)，函數(shù)形式確定了，則模型的復(fù)雜度也就相應(yīng)確定了。但是選擇合適的目標(biāo)函數(shù)的形式是一項十分困難的任務(wù)，如果目標(biāo)函數(shù)的形式不符合數(shù)據(jù)的真實分布，則預(yù)測效果會變得很差。相比之下，GPR屬于非參數(shù)模型，其復(fù)雜度隨著樣本數(shù)量的增加而增加，不需要人為做相應(yīng)的假設(shè)，所以GPR模型更加靈活，預(yù)測精度更高[10]。

連續(xù)式跨聲速風(fēng)洞在連續(xù)變迎角模式下，馬赫數(shù)控制要達(dá)到±0.001的精度，則壓縮機(jī)輸出轉(zhuǎn)速與目標(biāo)轉(zhuǎn)速之間的誤差應(yīng)該在±2 r/min之內(nèi)。圖6為Add-GPR 5折交叉驗證實驗中，測試集實際轉(zhuǎn)速與預(yù)測轉(zhuǎn)速之間的絕對誤差分布直方圖。從圖6中可以看出，絕對誤差基本都分布在2 r/min之內(nèi)，可以滿足連續(xù)變迎角方式下的馬赫數(shù)控制要求。

圖6 絕對誤差直方圖Fig.6 Absolute error histogram

3.2 分組劃分?jǐn)?shù)據(jù)集

一個好的預(yù)測模型應(yīng)該對于遠(yuǎn)離訓(xùn)練集的樣本依然有很好的預(yù)測能力。但是隨機(jī)劃分?jǐn)?shù)據(jù)集時，雖然測試樣本對于模型來說是未見的，但其相對于訓(xùn)練集的距離是比較近的。為了進(jìn)一步說明Add-GPR對風(fēng)洞數(shù)據(jù)的建模能力，我們采用分組劃分?jǐn)?shù)據(jù)集的方式產(chǎn)生訓(xùn)練集跟測試集。由于在3.1節(jié)實驗中，S-GPR與Add-GPR的預(yù)測效果最好，所以本節(jié)只對這兩種方法進(jìn)行比較。

由圖4可知，共有10組在不同馬赫數(shù)條件下測得的數(shù)據(jù)，我們以其中9組作為訓(xùn)練集，剩余1組作為測試集共進(jìn)行了4次實驗，每次實驗分別將馬赫數(shù)0.8、0.9、1.0、1.1條件下的數(shù)據(jù)作為測試集，最終得到的轉(zhuǎn)速預(yù)測曲線如圖7所示。

圖7中，左圖為不同馬赫數(shù)下Add-GPR與S-GPR的預(yù)測轉(zhuǎn)速曲線，右圖為兩者預(yù)測的絕對誤差?？梢钥闯?，S-GPR的預(yù)測曲線跟真實的轉(zhuǎn)速曲線相差很大，而Add-GPR的預(yù)測雖然跟真實曲線也存在一定的差距，但是各個攻角處預(yù)測的絕對誤差趨于穩(wěn)定，這說明Add-GPR預(yù)測曲線與真實曲線的變化趨勢是基本一致的。實驗結(jié)果顯示，馬赫數(shù)1.1條件下的預(yù)測相對來說比較差，這主要是因為馬赫數(shù)大于1.1的條件下是沒有實驗數(shù)據(jù)的，所以建模時只能利用馬赫數(shù)小于1.1一側(cè)的數(shù)據(jù)信息，而其他馬赫數(shù)條件的兩側(cè)都存在數(shù)據(jù)，可以利用的信息更多。

(a)馬赫數(shù)0.8下的預(yù)測結(jié)果

(b)馬赫數(shù)0.9下的預(yù)測結(jié)果

(c)馬赫數(shù)1.0下的預(yù)測結(jié)果

(d)馬赫數(shù)1.1下的預(yù)測結(jié)果圖7 轉(zhuǎn)速預(yù)測曲線Fig.7 Predictive curve of rotational speed

接下來將說明基于所提出的控制策略，只需預(yù)測轉(zhuǎn)速曲線與真實轉(zhuǎn)速曲線的變化趨勢一致，就可以滿足馬赫數(shù)控制的要求。

3.3 控制策略

在連續(xù)變迎角實驗過程中，如果我們要控制某一馬赫數(shù)穩(wěn)定，只需將要固定的馬赫數(shù)、連續(xù)變化的迎角以及實驗段截面積輸入訓(xùn)練好的模型，就可以提前得到不同迎角下，為了保持馬赫數(shù)穩(wěn)定，壓縮機(jī)應(yīng)該輸出的轉(zhuǎn)速。此時t時刻產(chǎn)生的總控制信號為：

U(t)=Uff(t)+Ufb(t)

(21)

其中，

Uff(t)=S(α(t+τ))-S(α(t+τ-1))

(22)

KD[e(t)-e(t-1)]

(23)

式(21)～式(23)中，Uff為前饋控制部分的輸出信號，α(t+τ)為t+τ時刻的迎角，S(α(t+τ))為該迎角下保持馬赫數(shù)穩(wěn)定所需要的壓縮機(jī)轉(zhuǎn)速；Ufb為PID反饋控制的輸出信號，τ為PID控制算法的延遲時間，其中前饋部分先于PID反饋τ時間進(jìn)行控制。

由式(22)可知，只要預(yù)測曲線與真實曲線的變化趨勢一致，那么前饋部分輸出的信號就是相同的，即使預(yù)測轉(zhuǎn)速與真實轉(zhuǎn)速之間存在偏差，也不影響最終的控制效果。前饋部分預(yù)測的越準(zhǔn)確，馬赫數(shù)在PID反饋控制的延遲時間內(nèi)產(chǎn)生的偏移就越小，此時通過PID反饋控制就可以達(dá)到馬赫數(shù)的控制精度要求。

訓(xùn)練數(shù)據(jù)很多時，我們可以得到比較準(zhǔn)確的壓縮機(jī)轉(zhuǎn)速預(yù)測曲線。但是很多情況下，我們在某一馬赫數(shù)條件下只能得到有限的實驗數(shù)據(jù)，此時預(yù)測結(jié)果可能無法達(dá)到連續(xù)變迎角模式下的馬赫數(shù)控制精度要求。由于GPR模型是一種概率模型，它可以得到預(yù)測結(jié)果的方差，以評估預(yù)測結(jié)果的置信度，所以當(dāng)控制精度不足時，可以根據(jù)預(yù)測方差來改善預(yù)測精度。在預(yù)測方差比較大的區(qū)域有針對性地增加實驗數(shù)據(jù)來提高預(yù)測的精度，這種有目的的改善可以使用較少的數(shù)據(jù)達(dá)到很好的預(yù)測結(jié)果。我們以馬赫數(shù)0.9條件下的數(shù)據(jù)為例來說明該方法的有效性。實驗結(jié)果如圖8所示。

圖8 提高預(yù)測精度Fig.8 Improving predictive accuracy

圖8中，紅色曲線表示真實的轉(zhuǎn)速曲線，藍(lán)色曲線表示預(yù)測結(jié)果，陰影部分為預(yù)測值95%的置信區(qū)間。左圖為只有14個樣本點(diǎn)的實驗結(jié)果，右圖是在左圖預(yù)測結(jié)果不確定性最大的區(qū)域額外加入14個訓(xùn)練樣本后的實驗結(jié)果?？梢园l(fā)現(xiàn)右圖轉(zhuǎn)速曲線的預(yù)測結(jié)果有一定改善，并且預(yù)測結(jié)果的置信度得到了顯著提高，因而PID反饋控制的精度得到了保證。

4 結(jié) 論

本文提出了一種基于GPR的前饋控制與PID反饋控制結(jié)合的控制策略，該策略可以提高連續(xù)變迎角實驗?zāi)Ｊ较碌鸟R赫數(shù)控制精度,并且具有簡單、節(jié)能等優(yōu)點(diǎn)。前饋部分將迎角、實際馬赫數(shù)和風(fēng)洞實驗段截面積作為Add-GPR的輸入來預(yù)測壓縮機(jī)轉(zhuǎn)速曲線，通過實驗證明了Add-GPR對風(fēng)洞實驗數(shù)據(jù)建模的有效性。同時利用GPR輸出的預(yù)測結(jié)果的置信度,可以指導(dǎo)新增實驗數(shù)據(jù)的獲取，以進(jìn)一步提高馬赫數(shù)控制精度。由于風(fēng)洞實驗條件的限制，本文只考慮了3個因素對壓縮機(jī)轉(zhuǎn)速的影響，而且其中截面積只有兩種不同的取值，這其實限制了模型泛化能力[20]。如果實驗條件允許，可以采集多組不同工況下的實驗數(shù)據(jù)，并且加入其他因素作為模型的輸入，比如溫度、壓強(qiáng)等，這些措施都可以提高轉(zhuǎn)速曲線的預(yù)測精度，進(jìn)而使馬赫數(shù)的控制更加穩(wěn)定。