曹 瑞，劉燕斌，陸宇平

(1.南京航空航天大學(xué)自動化學(xué)院，南京 211106；2.南京航空航天大學(xué)航天學(xué)院，南京 211106)

0 引 言

高超聲速飛行器因飛行速域和空域大、系統(tǒng)間耦合強烈，所以具有復(fù)雜的氣動模型。此外，由于高超聲速飛行器飛行過程中產(chǎn)生的尾流以及熱膨脹效應(yīng)，改變了飛行器周圍流場從而影響其模型參數(shù)，使得地面與空中實際氣動信息存在差異。這種復(fù)雜模型以及存在的空天差異使高超聲速飛行器的控制更具有挑戰(zhàn)性。因此，考慮采用系統(tǒng)辨識(System identification,SI)方法，根據(jù)飛行測量數(shù)據(jù)建立能夠反映系統(tǒng)實際運動特性且便于控制器設(shè)計的復(fù)雜系統(tǒng)等效模型。

在系統(tǒng)缺乏完整物理描述的情況下，SI可以提供缺失的信息片段，以構(gòu)建具有描述性或預(yù)測性的模型。雖然線性SI問題已經(jīng)得到了很好的解決，但非線性SI問題尚未有很好的解決方案。此外，對高超聲速飛行器系統(tǒng)辨識而言，還面臨著不確定性的問題。引起不確定性的其中一個因素是測量數(shù)據(jù)，其受噪聲以及測量變量本身隨機性的影響,具有可變性，進而導(dǎo)致參數(shù)估計的不確定性。另一個引起不確定性的因素源于模型辨識目的。研究者通過辨識算法獲取復(fù)雜系統(tǒng)的等效簡化模型以便用于控制器設(shè)計、飛行品質(zhì)分析以及故障辨識，而對簡化模型的辨識會不可避免地引入不確定性。

經(jīng)典的SI方法難以解決受上述不確定性影響的系統(tǒng)辨識問題。受Beck的開創(chuàng)性工作以及機器學(xué)習(xí)的指導(dǎo)，研究者認識到可以在貝葉斯原理的基礎(chǔ)上構(gòu)建更具魯棒性的系統(tǒng)辨識方法?？紤]到飛行器模型結(jié)構(gòu)已知，因此其系統(tǒng)辨識問題是典型的白盒辨識問題。而在白盒模型的背景下，貝葉斯方法的應(yīng)用雖未十分廣泛，但在很早之前便有了相關(guān)方面的研究。迄今為止，Beck等學(xué)者對貝葉斯SI進行了最系統(tǒng)的研究。文獻[7]對系統(tǒng)辨識方面的早期工作進行了總結(jié)；文獻[10]介紹了系統(tǒng)辨識向貝葉斯框架的過渡，并使用拉普拉斯近似解決了高維積分的難題。此后，文獻[11]又引入了馬爾可夫蒙特卡洛(Markov chain Monte Carlo,MCMC)方法，作為計算高維積分更通用的方法。文獻[12]采用貝葉斯方法研究了滯后系統(tǒng)模型—Masing模型的辨識問題。此外，貝葉斯方法也在飛行器動力學(xué)辨識中顯示出可觀的前景。在目前的研究中，比較有趣的兩項研究是貝葉斯SI的概率邏輯引入以及選擇信息量最大的訓(xùn)練數(shù)據(jù)以降低MCMC計算負擔(dān)。

對于高超聲速飛行器系統(tǒng)辨識，由于測量數(shù)據(jù)受噪聲干擾以及其等效簡化模型會引入不確定性，所以需要采集大量的實驗數(shù)據(jù)以評估這些不確定性。但是，MCMC算法本身的缺陷導(dǎo)致其難以很好地利用完整的數(shù)據(jù)集。為了解決這一問題，本文將Green提出的數(shù)據(jù)退火概念與MCMC方法相結(jié)合，并引入香農(nóng)熵制定退火策略，提出了基于數(shù)據(jù)退火的馬爾可夫鏈蒙特卡洛(Markov chain Monte Carlo based on data annealing,MCMCDA)算法以更好地分析大數(shù)據(jù)集。并以此為基礎(chǔ)，提出了基于MCMCDA的貝葉斯系統(tǒng)辨識方法，該方法辨識出的是參數(shù)的概率密度信息而非確定值。

最后，通過非線性示例對本文提出的方法進行驗證,并在有/無噪聲干擾情況下利用本文方法對高超聲速飛行器的短周期等效簡化系統(tǒng)進行辨識，然后基于獲得的簡化系統(tǒng)進行動態(tài)逆控制器設(shè)計及仿真，進一步說明本文所提系統(tǒng)辨識方法的有效性。

1 問題描述

1.1 模 型

高超聲速飛行器機身/推進一體化的設(shè)計構(gòu)型，使得飛行迎角和速度會影響發(fā)動機的燃料注入量以及工作狀態(tài)，而燃料注入量作為控制輸入又會影響飛行器的運動狀態(tài)。因此高超聲速飛行器是一個系統(tǒng)狀態(tài)與輸入之間相互影響、相互約束的復(fù)雜系統(tǒng)。這一復(fù)雜的動態(tài)系統(tǒng)模型可以表示為

()=(?())

(1)

式中：?()=[(),…,()),(),…,()]為+維向量，()為系統(tǒng)輸出，(),=1,2,…,為系統(tǒng)輸入，(),=1,2,…,表示系統(tǒng)狀態(tài)，下標(biāo)和分別為系統(tǒng)狀態(tài)階數(shù)和輸入階數(shù)，T是轉(zhuǎn)置符號,(·)表示未知的復(fù)雜非線性函數(shù)。

控制輸入和系統(tǒng)狀態(tài)之間相互耦合，導(dǎo)致其很難同常規(guī)飛行器一樣采用解析的方式進行解耦。因此，高超聲速飛行器這類復(fù)雜的非線性運動學(xué)模型難以應(yīng)用于典型的基于模型的控制器設(shè)計，比如動態(tài)逆控制。為了獲得便于控制器設(shè)計的模型，考慮建立此復(fù)雜系統(tǒng)的等效簡化模型?；陲w行采樣數(shù)據(jù)，采用系統(tǒng)辨識方法獲得能夠反映真實系統(tǒng)動態(tài)響應(yīng)的等效簡化模型，從而將其用于控制器設(shè)計、飛行品質(zhì)評估，以及飛行器故障檢測。

轉(zhuǎn)換式(1)所示的模型，獲得適用于控制器設(shè)計的稀疏模型，如狀態(tài)-空間模型(2)或傳遞函數(shù)模型(3)。

(2)

(3)

為了解決模型簡化和測量噪聲所引入的不確定性，本文采用貝葉斯辨識方法，對飛行器的模型信息進行辨識。通過式(4)量化參數(shù)估計過程中的不確定性，認為未知參數(shù)是服從某種概率分布的隨機變量而非確定值，其結(jié)果是未知參數(shù)的最優(yōu)估計值。

(4)

(5)

由式(4)可知，貝葉斯方法能夠使用觀測數(shù)據(jù)對模型未知參數(shù)的先驗概率分布進行修正，獲得可以正確描述參數(shù)信息的后驗概率分布。

1.2 貝葉斯的高維求解問題

(6)

若考慮上述不確定性，則是以數(shù)據(jù)為條件的隨機向量。在這種情況下，相較于獲得的確定估計值，更希望指定的置信域。從概率的角度來說，的表達式應(yīng)為式(7)，而非式(6)。

(7)

系統(tǒng)辨識的目標(biāo)是獲得一個預(yù)測模型，提供系統(tǒng)輸入，其可以估計或預(yù)測系統(tǒng)輸出。在上述概率環(huán)境中，系統(tǒng)辨識所獲得的預(yù)測模型則需要確定系統(tǒng)輸出的分布情況。假設(shè)將一個新的輸入數(shù)據(jù)應(yīng)用于預(yù)測模型，則考慮所有相關(guān)性的預(yù)測輸出的概率密度為

(8)

此分布的均值將給出預(yù)測的“最佳”估計值，協(xié)方差將建立預(yù)測的置信區(qū)間。式(8)中的預(yù)測輸出密度由下式計算

(9)

(10)

這種積分通常難以處理，而且高維性質(zhì)使其在數(shù)值計算上過于復(fù)雜(>2)。本文采用統(tǒng)計學(xué)中的MCMC算法，通過模擬抽樣近似求解的方式解決貝葉斯估計中的高維積分求解問題。然而，MCMC算法本身存在的一些應(yīng)用問題可能會影響最終結(jié)果的精度。為了改善MCMC算法,獲得更好的辨識結(jié)果，下面將對這些問題進行進一步討論。

1.3 MCMC的局部陷阱和訓(xùn)練數(shù)據(jù)問題

圖1 MCMC陷入“局部陷阱”Fig.1 MCMC becoming stuck in a “l(fā)ocal trap”

另一種是使用文獻[15]中描述的方法，選擇一個信息量大的子數(shù)據(jù)集推斷出參數(shù)分布概率。為判斷數(shù)據(jù)集的信息含量，采用了香農(nóng)熵對數(shù)據(jù)信息含量進行度量，下面對這一過程進行描述。

首先，將后驗參數(shù)分布寫成如下形式：

(11)

(12)

(13)

然后，可以使用香農(nóng)熵的相關(guān)定義來表征該分布的信息含量：

(14)

因此，為了解決MCMC算法易于陷入局部最優(yōu)以及大數(shù)據(jù)下難以很好利用數(shù)據(jù)信息的問題，本文結(jié)合Green提出的數(shù)據(jù)退火概念，提出了一種包含模擬退火和大信息量訓(xùn)練數(shù)據(jù)概念的MCMCDA算法。該算法主要解決局部最優(yōu)問題并克服了上述兩種思路對于大數(shù)據(jù)問題的缺點。之后，將該算法與貝葉斯估計方法相結(jié)合進一步提出基于MCMCDA的系統(tǒng)辨識方法。

2 改進MCMC算法

2.1 MCMCDA算法的基本思路

(15)

(16)

MCMCDA算法采用香農(nóng)熵控制的增加速率，從而控制數(shù)據(jù)信息的引入。下面將討論如何應(yīng)用香農(nóng)熵設(shè)置退火計劃。

..基于香農(nóng)熵的退火計劃

MCMCDA算法的最佳退火計劃是允許香農(nóng)熵所測量的信息量能夠以恒定速率變化，使得對參數(shù)幾乎沒有影響的數(shù)據(jù)被迅速退火，從而允許算法關(guān)注“信息豐富”的數(shù)據(jù)。為便于引入香農(nóng)熵，將目標(biāo)參數(shù)的概率分布(16)改寫為

(17)

(18)

(19)

(20)

(21)

結(jié)合式(14)、(17)和(18)可得目標(biāo)分布的香農(nóng)熵為

(22)

然后，退火計劃(取值)對香農(nóng)熵測量信息量的影響可由下式進行評估

(23)

基于式(19)～(21)，式(23)中的第一項可按下式計算

(24)

(25)

(26)

結(jié)合式(24)、(25)和(26)可得

(27)

因此，如果該算法當(dāng)前使用的是，并且希望在香農(nóng)熵增量Δ不變的情況下“退火”新數(shù)據(jù)，則應(yīng)根據(jù)以下條件更新退火值為+1

(28)

如果考慮算法的初始階段，即僅對第一組數(shù)據(jù)進行退火，則式(28)可簡化為

(29)

選擇一個小的初始可以確保初始參數(shù)分布與先驗相似，使MCMCDA可從目標(biāo)進行有效采樣。經(jīng)大量仿真實驗發(fā)現(xiàn)，設(shè)置初始值為=1×10能夠獲得較好的仿真結(jié)果。

..香農(nóng)熵增量Δ的選擇

(30)

具體而言，+1應(yīng)根據(jù)以下方式選擇

(31)

但受制于以下條件

<+1≤1,-Δ<Δ<Δ

(32)

式中：Δ根據(jù)期望的數(shù)據(jù)信息變化速率進行定義，并且滿足Δ∈(0,1]。

2.2 MCMCDA算法的流程

這里總結(jié)了MCMCDA算法的執(zhí)行流程：

算法1：MCMCDA算法

1)設(shè)置迭代次數(shù)=1，以及參數(shù)初值。

2)設(shè)置=1，=。

3)判斷是否<1。若滿足，進入步驟4；若不滿足，進入步驟7。

4)獲得參數(shù)值。根據(jù)建議的概率分布(·|-1)獲得待選參數(shù)值，即～(·|-1)，由式(33)計算接受概率()：

(33)

6)根據(jù)式(31)和(32)求解+1,+1→，并返回步驟3。

7)判斷目標(biāo)參數(shù)標(biāo)準(zhǔn)差是否滿足()<。若滿足，則終止算法；否則，+→添加更多數(shù)據(jù)，并令+1→，返回步驟2。

本文設(shè)置為10。和根據(jù)實測數(shù)據(jù)設(shè)定。

3 基于MCMCDA的系統(tǒng)辨識方法

本文采用式(34)作為描述輸入數(shù)據(jù)集{0,1,…,,…}與輸出數(shù)據(jù)集{0,1,…,,…}之間的一類普遍結(jié)構(gòu)，下標(biāo)表示采樣時刻。信號集{0,1,…,,…}表示不同時刻下噪聲或者模型簡化對應(yīng)的誤差。本節(jié)以式(34)所示的模型為例，對基于MCMCDA的貝葉斯系統(tǒng)辨識方法的執(zhí)行流程進行詳細敘述。

=(,)+(,)

(34)

式中：(,)表達式如式(3)所示,(,)的表達式與(,)一致,區(qū)別為(,)代表信號至輸出通道的傳遞函數(shù)。信號是一個獨立的隨機過程，∈{,,…,,…}。為降低對系統(tǒng)參數(shù)辨識的影響，將式(34)改寫為

(35)

(36)

與觀測數(shù)據(jù)相關(guān)的似然函數(shù)可以寫為

(37)

式中：?[,,…,],(|)為似然度(|)的初始假設(shè)。然后，假定滿足如下分布：

～(·)

(38)

根據(jù)式(35)，式(36)和式(38)可以得出：

(39)

因此，結(jié)合式(39)并根據(jù)貝葉斯定理，后驗分布(|)的計算式如下

(40)

式中：()為的先驗分布，()是歸一化常數(shù)

(41)

聯(lián)合(36)，(40)和(41)，后驗概率(|)最終為

(,)(,))

(42)

式中：是一個獨立于的常數(shù)，用于對總概率進行歸一化處理。本文中考慮參數(shù)的先驗概率分布()為均勻分布。

然后，將式(42)與式(15)相結(jié)合，進一步獲得算法2中所使用的目標(biāo)參數(shù)分布：

(43)

最后基于式(43)，采用MCMCDA系統(tǒng)辨識方法獲取參數(shù)的概率分布，其執(zhí)行流程如下：

算法2：基于MCMCDA系統(tǒng)辨識方法

1)設(shè)置最大迭代次數(shù)=500，初始引入數(shù)據(jù)長度為，誤差閾值=10。并設(shè)置算法2中每次循環(huán)添加的數(shù)據(jù)量為=03。

2)設(shè)置=1，初始化待辨識參數(shù)為-1。

4)執(zhí)行算法1，獲得待估計參數(shù)的分布。

4 仿真校驗

本節(jié)將通過一個示例模型說明所提出算法的有效性，并將算法應(yīng)用到某高超聲速飛行器短周期等效簡化模型的辨識中。

4.1 隨機系統(tǒng)仿真示例

考慮一個存在外部干擾的系統(tǒng)，采樣數(shù)據(jù)為位移(單位m)信息

(44)

完整的訓(xùn)練數(shù)據(jù)集包含10000個位移測量值，在算法開始時首次引入=1000個數(shù)據(jù)點，設(shè)置Δ=1。圖2展示了算法迭代過程中參數(shù)標(biāo)準(zhǔn)差的變化，橫坐標(biāo)表示迭代次數(shù)。每當(dāng)=1時，都會引入一個新的數(shù)據(jù)段，新數(shù)據(jù)段的長度為=03。當(dāng)某些參數(shù)相對缺乏數(shù)據(jù)信息時，也可通過立即設(shè)置=1實現(xiàn)香農(nóng)熵的期望變化。根據(jù)圖2每次引入新數(shù)據(jù)后參數(shù)標(biāo)準(zhǔn)差的變化可以得知，在很大程度上,引入更多的數(shù)據(jù)始終會降低參數(shù)的不確定性。

圖2 參數(shù)標(biāo)準(zhǔn)偏差變化Fig.2 Variation of parameter standard deviation

圖3展示了基于MCMCDA系統(tǒng)辨識方法的參數(shù)辨識結(jié)果，其輸出的是參數(shù)概率分布。其中,縱軸表示迭代中參數(shù),,在不同值下出現(xiàn)的次數(shù)。

圖3 參數(shù)辨識結(jié)果Fig.3 Parameter identification result

根據(jù)圖3可知，各參數(shù)分布的均值非常接近待辨識參數(shù)的真實值，并且變化區(qū)間較小，驗證了算法的有效性。

4.2 飛行器系統(tǒng)辨識中的應(yīng)用

..短周期簡化模型辨識

本文采用通用高超聲速(Generic hypersonic aerodynamics model example,GHAME)飛行器作為高超聲速飛行器的示例，并且分別用差分進化(Differential evolution,DE)算法、MCMC方法和本文所提算法對GHAME系統(tǒng)進行仿真實驗，驗證本文算法的有效性和可靠性。

以GHAME飛行器在6、15 km狀態(tài)下，迎角、俯仰角速率的測量數(shù)據(jù)和舵面數(shù)據(jù)作為辨識系統(tǒng)的輸入信息,并采用短周期狀態(tài)-空間模型作為GHAME飛行器縱向姿態(tài)運動的等效簡化模型

(45)

其中,,,,,為待辨識參數(shù)。

分別采用本文算法、MCMC方法以及DE算法對模型(45)進行辨識。圖4展示了三種算法下參數(shù)均值的收斂過程，表1描述了最終參數(shù)估計結(jié)果的均值、方差(括號中的數(shù)值表示方差，括號外的為均值)，以及俯仰角速率預(yù)測數(shù)據(jù)與其真實數(shù)據(jù)的誤差。

根據(jù)圖4和表1可知，本文提出的基于MCMCDA的辨識算法與DE算法相比，收斂速度更快，辨識精度更高。與MCMC算法相比，本文算法允許待估計參數(shù)能夠在算法早期階段進行較大變化，以跳出局部最優(yōu)；且當(dāng)參數(shù)變化較為平穩(wěn)時，可以在局部參數(shù)空間內(nèi)進行更精細的搜索，提高辨識精度。為了進一步驗證本文提出的算法在系統(tǒng)存在測量噪聲時的辨識能力，對標(biāo)稱信號引入均值為0，方差為=0.002的高斯白信號作為測量噪聲，并將其作為輸入數(shù)據(jù)。隨后，基于受干擾的測量數(shù)據(jù)，分別采用三種辨識算法進行系統(tǒng)辨識。

圖4 不同算法下辨識迭代過程Fig.4 Identification iterative process under different algorithms

表1 不同算法下得到的參數(shù)平均值與方差(括號內(nèi)為方差)Table 1 The mean and variance of parameters obtained by different algorithms (variance information is in brackets)

圖5展示了500次蒙特卡洛仿真下參數(shù)估計結(jié)果所繪制的箱線圖。很明顯，DE算法的辨識結(jié)果最差，且每次估計結(jié)果相差最大；基于MCMCDA辨識算法獲得的結(jié)果接近實際值，并且在500次估計中獲得的參數(shù)值分布均較為集中?？梢缘贸觯贛CMCDA的系統(tǒng)辨識方法可以對動態(tài)系統(tǒng)進行有效估計，并且辨識結(jié)果具有較好的魯棒性。

圖5 不同算法下500次運行結(jié)果生成的箱線圖Fig.5 The box-line diagram generated by 500 running results under different algorithms

此外，圖6展示了基于MCMCDA的系統(tǒng)辨識方法輸出的等效簡化系統(tǒng)與原復(fù)雜飛行器系統(tǒng)的頻域特性對比圖。圖中實線是真實系統(tǒng)的伯德圖，虛線是辨識參數(shù)均值構(gòu)成的簡化系統(tǒng)伯德圖，陰影區(qū)域是95%置信區(qū)間。

圖6 伯德圖Fig.6 Bode diagram

根據(jù)上述仿真結(jié)果可知，本文算法不僅獲得了目標(biāo)參數(shù)的后驗均值，還可視化了整個模型的不確定性。這驗證了本文算法能夠?qū)?fù)雜系統(tǒng)進行有效簡化。

..控制應(yīng)用

本節(jié)將獲得的簡化模型應(yīng)用于動態(tài)逆控制中。動態(tài)逆控制效果對飛行器氣動導(dǎo)數(shù)信息敏感，因此可以更好地驗證上節(jié)獲得的簡化模型是否能夠有效地應(yīng)用于飛行器控制系統(tǒng)設(shè)計中。

圖7 模型參考動態(tài)逆控制結(jié)構(gòu)Fig.7 Control structure of model reference dynamic inversion

根據(jù)動態(tài)逆理論，得出升降舵控制輸入的表達式如下所示：

(46)

(47)

式中：和為俯仰角速率和迎角的測量信息。

(48)

然后，將獲得的控制輸入和′分別代入原飛行器復(fù)雜模型中進行仿真，獲得的仿真結(jié)果如圖8所示。圖中點劃線為開環(huán)下高階系統(tǒng)包含測量噪聲的短周期信號；實線表示跟蹤指令；虛線為式(47)所示控制輸入下的響應(yīng)曲線;點線為式(48)所示控制輸入′下的響應(yīng)曲線。

圖8 動態(tài)逆控制響應(yīng)曲線Fig.8 Response curve of dynamic inverse control

從圖8可以看出，在開環(huán)情況下飛行器姿態(tài)不穩(wěn)定，加入動態(tài)逆控制輸入后可以跟蹤目標(biāo)指令。此外，將辨識出的噪聲統(tǒng)計信息引入控制設(shè)計可以更好消除噪聲所帶來的影響，進一步改善了控制結(jié)果。因此，本文所提辨識方法獲得的簡化模型可以反映出系統(tǒng)特性，并且適用于控制系統(tǒng)設(shè)計。

綜合上述仿真結(jié)果，驗證了本文提出的基于MCMCDA系統(tǒng)辨識方法的有效性，并且辨識出的簡化模型能夠有效地應(yīng)用于飛行器控制系統(tǒng)設(shè)計之中。

5 結(jié) 論

本文首先提出了MCMCDA算法，該算法設(shè)計了基于香農(nóng)熵的新數(shù)據(jù)引入策略，通過確保目標(biāo)參數(shù)概率分布的香農(nóng)熵變化保持在預(yù)定義的范圍內(nèi)，使得算法能夠以平滑和連續(xù)的方式“吸收”具有較大信息量的訓(xùn)練數(shù)據(jù)。此外，MCMCDA算法融合了模擬退火的思想，能夠在算法早期對待估計參數(shù)進行較大的調(diào)整，解決了MCMC算法中易出現(xiàn)的“局部陷阱”問題。然后，將這一算法應(yīng)用于貝葉斯辨識，進一步提出基于MCMCDA的系統(tǒng)辨識方法。最后，通過仿真驗證了這一系統(tǒng)辨識方法的有效性，并且辨識出的等效簡化模型能夠有效應(yīng)用于控制器設(shè)計之中，取得了較好的控制效果。