基于馬爾可夫蒙特卡洛法的系統(tǒng)辨識方法研究及應用

曹 瑞，劉燕斌，陸宇平

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

0 引 言

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

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

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

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

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

1 問題描述

1.1 模 型

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

()=(?())

(1)

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

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

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

(2)

(3)

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

(4)

(5)

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

1.2 貝葉斯的高維求解問題

(6)

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

(7)

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

(8)

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

(9)

(10)

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

1.3 MCMC的局部陷阱和訓練數據問題

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

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

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

(11)

(12)

(13)

然后，可以使用香農熵的相關定義來表征該分布的信息含量：

(14)

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

2 改進MCMC算法

2.1 MCMCDA算法的基本思路

(15)

(16)

MCMCDA算法采用香農熵控制的增加速率，從而控制數據信息的引入。下面將討論如何應用香農熵設置退火計劃。

..基于香農熵的退火計劃

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

(17)

(18)

(19)

(20)

(21)

結合式(14)、(17)和(18)可得目標分布的香農熵為

(22)

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

(23)

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

(24)

(25)

(26)

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

(27)

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

(28)

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

(29)

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

..香農熵增量Δ的選擇

(30)

具體而言，+1應根據以下方式選擇

(31)

但受制于以下條件

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

(32)

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

2.2 MCMCDA算法的流程

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

算法1：MCMCDA算法

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

2)設置=1，=。

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

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

(33)

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

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

本文設置為10。和根據實測數據設定。

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

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

=(,)+(,)

(34)

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

(35)

(36)

與觀測數據相關的似然函數可以寫為

(37)

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

～(·)

(38)

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

(39)

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

(40)

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

(41)

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

(,)(,))

(42)

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

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

(43)

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

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

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

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

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

4 仿真校驗

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

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

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

(44)

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

圖2 參數標準偏差變化Fig.2 Variation of parameter standard deviation

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

圖3 參數辨識結果Fig.3 Parameter identification result

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

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

..短周期簡化模型辨識

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

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

(45)

其中,,,,,為待辨識參數。

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

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

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

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

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

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

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

圖6 伯德圖Fig.6 Bode diagram

根據上述仿真結果可知，本文算法不僅獲得了目標參數的后驗均值，還可視化了整個模型的不確定性。這驗證了本文算法能夠對復雜系統(tǒng)進行有效簡化。

..控制應用

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

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

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

(46)

(47)

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

(48)

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

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

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

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

5 結 論

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