基于高斯過程回歸的連續式風洞馬赫數控制

高 赫,劉學軍,郭 晉,呂宏強

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

0 引 言

風洞是能人工產生和控制氣流，以模擬飛行器或物體周圍氣體流動，并可量度氣流對物體的作用以及觀察物理現象的一種管道實驗設備，它是進行空氣動力實驗最常用、最有效的工具之一[1]。在進行風洞實驗時，保持風洞實驗段馬赫數的穩定，以測量飛行器在不同迎角下的氣動性能，對于研制高性能飛行器具有重要意義。為了獲得準確的實驗數據，階梯變迎角實驗模式被廣泛應用在風洞實驗中。階梯變迎角模式基于迎角的變化范圍間歇地改變迎角以保持馬赫數的穩定，但這種方式只能獲得少量特定迎角處的實驗數據，實驗效率不高，而且在每一個迎角處都要等待馬赫數穩定后再測量相應的氣動性能，這樣就不可避免地延長了實驗時間，造成較高的能源消耗，對于風洞這種高耗能設備來說是不經濟的。連續變迎角實驗模式具有更高的實驗效率和較低的運行成本，并且可以獲得更密集的實驗數據，這些優點促使風洞實驗方式開始從階梯變迎角向連續變迎角方式轉變。但與此同時,連續變迎角實驗模式的快速時變特性也給風洞馬赫數的控制帶來了新的挑戰，需要研究新的控制算法與策略。

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

針對PID控制算法存在的缺陷，近年來出現了很多基于預測的控制方法，這些方法統稱為模型預測控制(Model Predictive Control，MPC)算法[4]。許多研究者已經通過仿真實驗證明，使用模型預測控制可以很好地對連續變迎角實驗模式下的馬赫數進行控制。文獻[5]使用統一模型預測控制對阿姆斯特丹航空航天實驗室擁有的高速風洞進行馬赫數控制，實驗結果表明該控制方法與傳統PID控制算法相比，抗迎角擾動能力提高4倍并且馬赫數的控制性能提高30%～60%。文獻[6]出于馬赫數高精度控制的考慮，選取穩定段總壓和實驗段靜壓作為被控量，設計了基于迎角補償的模型預測控制器，使馬赫數的控制精度達到了大飛機研制要求。但是模型預測控制方法需要在線求解約束優化問題，計算量大，只適用于具有高性能計算機的環境[7]。

中國航空工業空氣動力研究院的0.6 m×0.6 m連續式跨聲速風洞將壓縮機轉速作為控制馬赫數穩定的唯一驅動方式。為了解決PID控制方法在連續變迎角實驗模式下存在的缺陷，同時盡可能節約能源，其目前主要采用前饋控制與PID反饋控制結合的策略。前饋部分的設計需要先通過固定轉速，連續變化迎角，得到迎角與實際馬赫數的關系曲線；然后再通過固定迎角，增加轉速，得到轉速與實際馬赫數的關系曲線；最后通過這兩條曲線得到前饋系數，即轉速變化量與迎角變化量的比值，這樣就可以根據迎角變化量提前對轉速進行補償。但是這種方式存在兩個明顯缺點：1)為了得到前饋系數需要進行兩次實驗，增加了實驗次數；2)不同迎角下，轉速與實際馬赫數的關系曲線是存在一定差別的，導致前饋系數計算不準確。

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

1 高斯過程回歸(GPR)

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

1.1 預測

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

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)為均值函數，k(x,x′)為協方差函數。

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

yi=f(xi)+ε

(4)

(5)

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

對于給定的測試樣本x*，其對應的預測值f*與觀測值y的聯合先驗分布為：

(6)

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

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

(7)

其中，

(8)

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

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

(9)

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

1.2 訓練

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

(10)

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

(11)

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

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

(12)

θ*=argminθL(θ)

(13)

2 實驗設計

2.1 數據預處理

使用J7標模在連續變迎角實驗模式下得到的數據進行建模，其中迎角變化范圍為-5°～11°。在實際實驗過程中，由于迎角回零操作或者其他一些復雜因素的影響可能會產生一些異常數據，這些異常數據會嚴重影響建模的準確性，所以在建模前需要將這些數據進行刪除。刪除前與刪除后的數據分別如圖1、圖2所示，圖中的虛線表示馬赫數±0.001誤差帶。

圖1 包含異常點的數據Fig.1 Data with outliers

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

圖2 清洗后的數據Fig.2 Cleaned data

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

(14)

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

2.2 初始超參數選擇

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

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

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

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

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

本文最終選擇標準高斯分布作為超參數的先驗分布，從中隨機采樣20次，得到20組不同的初始超參數。對于每組超參數使用共軛梯度下降算法進行500次迭代優化，從中選取使負對數邊際似然最小的取值作為最終的初始超參數。

2.3 協方差函數的選擇

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

S-GPR假設目標函數具有如下形式：

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

(15)

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

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

(16)

(17)

(18)

(19)

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

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

3 實驗結果及分析

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

圖4 輸入向量散點圖Fig.4 Scatter plot of input vector

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

圖5 實際馬赫數散點圖Fig.5 Scatter plot of actual Mach numbers

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

3.1 隨機劃分數據集

使用機器學習算法執行預測任務有一個基本的假設，即訓練數據與測試數據是獨立同分布的。基于這個假設我們將2578個樣本隨機打亂，使用常用的機器學習模型分別進行5折交叉驗證，最終得到預測結果的RMSE(Root Mean Square Error)如表1所示。其中，RMSE的計算公式為：

(20)

表1 預測結果對比Table 1 Comparison of prediction results

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

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

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

3.2 分組劃分數據集

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

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

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

(a)馬赫數0.8下的預測結果

(b)馬赫數0.9下的預測結果

(c)馬赫數1.0下的預測結果

(d)馬赫數1.1下的預測結果圖7 轉速預測曲線Fig.7 Predictive curve of rotational speed

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

3.3 控制策略

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

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+τ))為該迎角下保持馬赫數穩定所需要的壓縮機轉速；Ufb為PID反饋控制的輸出信號，τ為PID控制算法的延遲時間，其中前饋部分先于PID反饋τ時間進行控制。

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

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

圖8 提高預測精度Fig.8 Improving predictive accuracy

圖8中，紅色曲線表示真實的轉速曲線，藍色曲線表示預測結果，陰影部分為預測值95%的置信區間。左圖為只有14個樣本點的實驗結果，右圖是在左圖預測結果不確定性最大的區域額外加入14個訓練樣本后的實驗結果。可以發現右圖轉速曲線的預測結果有一定改善，并且預測結果的置信度得到了顯著提高，因而PID反饋控制的精度得到了保證。

4 結 論

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