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

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

0 引 言

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

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

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

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

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

1 高斯過(guò)程回歸(GPR)

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

1.1 預(yù)測(cè)

高斯過(guò)程是一個(gè)隨機(jī)變量的集合，并且其中任意有限數(shù)量隨機(jī)變量的聯(lián)合分布為多元高斯分布。從函數(shù)空間的角度出發(fā)，高斯過(guò)程可以看作定義在函數(shù)f(x)上的一個(gè)分布，它的性質(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}，針對(duì)回歸問(wèn)題，考慮如下模型：

yi=f(xi)+ε

(4)

(5)

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

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

(6)

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

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

(7)

其中，

(8)

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

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

(9)

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

1.2 訓(xùn)練

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

(10)

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

(11)

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

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

(12)

θ*=argminθL(θ)

(13)

2 實(shí)驗(yàn)設(shè)計(jì)

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

使用J7標(biāo)模在連續(xù)變迎角實(shí)驗(yàn)?zāi)Ｊ较碌玫降臄?shù)據(jù)進(jìn)行建模，其中迎角變化范圍為-5°～11°。在實(shí)際實(shí)驗(yàn)過(guò)程中，由于迎角回零操作或者其他一些復(fù)雜因素的影響可能會(huì)產(chǎn)生一些異常數(shù)據(jù)，這些異常數(shù)據(jù)會(huì)嚴(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í)驗(yàn)數(shù)據(jù)中，為了檢驗(yàn)連續(xù)變迎角實(shí)驗(yàn)方式下角速度對(duì)馬赫數(shù)控制精度的影響，在某些馬赫數(shù)條件下可能測(cè)量了兩組實(shí)驗(yàn)數(shù)據(jù)，其迎角變化速度分別為0.1°/s和0.2°/s，但是由于角速度不影響迎角、實(shí)時(shí)馬赫數(shù)、實(shí)驗(yàn)段截面積和壓縮機(jī)轉(zhuǎn)速之間的關(guān)系，同時(shí)也為了保證不同馬赫數(shù)下的數(shù)據(jù)量盡量均衡，所以對(duì)于這種情況只取角速度為0.1°/s的數(shù)據(jù)進(jìn)行建模。在進(jìn)行實(shí)際風(fēng)洞實(shí)驗(yàn)時(shí)，不同范圍的馬赫數(shù)條件下，實(shí)驗(yàn)段截面積也可能會(huì)做相應(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)過(guò)標(biāo)準(zhǔn)化的，所以這樣設(shè)置是合理的)。經(jīng)過(guò)處理后的數(shù)據(jù)集的特征分布直方圖如圖3所示。

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

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

(14)

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

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

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

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

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

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

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

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

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

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

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

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

(15)

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

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

(16)

(17)

(18)

(19)

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

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

3 實(shí)驗(yàn)結(jié)果及分析

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

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

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

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

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

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

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

(20)

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

由表1可以看出，Add-GPR每一次交叉驗(yàn)證的預(yù)測(cè)精度都要高于其他常用的回歸模型，同時(shí)GPR的預(yù)測(cè)效果也比其他模型的預(yù)測(cè)精度要好。這是因?yàn)樯窠?jīng)網(wǎng)絡(luò)、SVR和多項(xiàng)式回歸屬于參數(shù)化模型[9]，參數(shù)化模型需要對(duì)目標(biāo)函數(shù)的形式做有限參數(shù)的假設(shè)，函數(shù)形式確定了，則模型的復(fù)雜度也就相應(yīng)確定了。但是選擇合適的目標(biāo)函數(shù)的形式是一項(xiàng)十分困難的任務(wù)，如果目標(biāo)函數(shù)的形式不符合數(shù)據(jù)的真實(shí)分布，則預(yù)測(cè)效果會(huì)變得很差。相比之下，GPR屬于非參數(shù)模型，其復(fù)雜度隨著樣本數(shù)量的增加而增加，不需要人為做相應(yīng)的假設(shè)，所以GPR模型更加靈活，預(yù)測(cè)精度更高[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折交叉驗(yàn)證實(shí)驗(yàn)中，測(cè)試集實(shí)際轉(zhuǎn)速與預(yù)測(cè)轉(zhuǎn)速之間的絕對(duì)誤差分布直方圖。從圖6中可以看出，絕對(duì)誤差基本都分布在2 r/min之內(nèi)，可以滿足連續(xù)變迎角方式下的馬赫數(shù)控制要求。

圖6 絕對(duì)誤差直方圖Fig.6 Absolute error histogram

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

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

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

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

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

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

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

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

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

3.3 控制策略

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

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為前饋控制部分的輸出信號(hào)，α(t+τ)為t+τ時(shí)刻的迎角，S(α(t+τ))為該迎角下保持馬赫數(shù)穩(wěn)定所需要的壓縮機(jī)轉(zhuǎn)速；Ufb為PID反饋控制的輸出信號(hào)，τ為PID控制算法的延遲時(shí)間，其中前饋部分先于PID反饋τ時(shí)間進(jìn)行控制。

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

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

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

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

4 結(jié) 論

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