鄭奕揚，倪 何，金家善

(1.海軍92118部隊，浙江 舟山 316000；2.海軍工程大學 動力工程學院，武漢 430033；3.海軍工程大學 船舶與海洋學院，武漢 430033)

蒸汽動力系統(tǒng)因其組成設(shè)備多、系統(tǒng)耦合關(guān)系復(fù)雜的特點，在不同的運行工況下都保持穩(wěn)定運行是相當困難的，因此蒸汽動力系統(tǒng)的運行穩(wěn)定性評估相當重要.目前，穩(wěn)定性評估在多系統(tǒng)耦合的復(fù)雜大系統(tǒng)中的應(yīng)用停留在理論研究的層面上[1]，不能全面評估系統(tǒng)的真實運行穩(wěn)定性.而系統(tǒng)在一段運行時間內(nèi)參數(shù)的時間序列則真實反應(yīng)了系統(tǒng)的運行狀態(tài)，分解參數(shù)的時間序列，對重構(gòu)后得到的趨勢項和擾動項進行預(yù)測分析可以從中提取反應(yīng)系統(tǒng)運行穩(wěn)定性的評估指標.

在非平穩(wěn)時間序列趨勢項提取方面，基于經(jīng)驗?zāi)B(tài)分解(EMD)及基本模式分量(IMF)重構(gòu)是一種常用的算法[2-4].由于傳統(tǒng)的經(jīng)驗?zāi)B(tài)分解不能完全提取所有IMF分量中包含的趨勢信息，出現(xiàn)了大量的改進算法[5-7].文獻[8]提出了集合經(jīng)驗?zāi)B(tài)分解(EEMD)、奇異值分解(SVD)和排列熵(PE)的趨勢提取方法，并驗證了該方法相對于傳統(tǒng)分解算法的優(yōu)越性，為本文研究提供了思路.但是，EEMD算法在原信號的基礎(chǔ)上添加了不同的白噪聲并通過重復(fù)集合平均進行抵消，分解效果取決于添加白噪聲的去除，重構(gòu)誤差較大.同時，在排列熵的計算上較為粗糙，無法保證以此為依據(jù)進行選取分量的合理性.在時間序列分析中，非平穩(wěn)時間序列的常用預(yù)測模型為整合滑動平均自回歸模型(ARIMA)[9]；由于該模型比較成熟且具有良好的通用性，所以本文在計算失穩(wěn)概率(PI)時，采用ARIMA模型對運行參數(shù)時間序列的趨勢項和擾動項進行預(yù)測.

綜上所述，本文提出了一種名為MSOP的單參數(shù)運行穩(wěn)定性評估方法，該方法是一種按照特定流程進行的復(fù)合評估方法，由中值回歸經(jīng)驗?zāi)B(tài)分解(MREMD)、奇異值分解、基于最佳算法參數(shù)排列熵(OAPPE)重組和失穩(wěn)概率(PI)計算等4個要素組成，其名稱MSOP即為這4個要素的英文首字母組合.該方法經(jīng)實際案例驗證，可以根據(jù)運行參數(shù)的時間序列來定量評估蒸汽動力系統(tǒng)單參數(shù)的運行穩(wěn)定性，可為蒸汽動力系統(tǒng)的整體運行穩(wěn)定性評估提供底層評估指標輸入，對蒸汽動力系統(tǒng)的在線運行穩(wěn)定性管理和決策研究提供技術(shù)基礎(chǔ).

1 基于MSOP的運行穩(wěn)定性評估方法

1.1 MREMD分解

MREMD[10]是經(jīng)驗?zāi)B(tài)分解的改進算法，通過優(yōu)化均值包絡(luò)的生成方式來抑制端點效應(yīng).對于輸入長度為N的時間序列s(t)，MREMD分解的具體步驟如下.

步驟1在s(t)的左右端點處采用自回歸(AR)模型進行延拓處理，找出時間序列所有極值點構(gòu)成的序列X={x1,x2,…,xk}，計算極值點序列的統(tǒng)計有義值.

P{|xi||(|xi|≤xσ)；i=1,2,…,k}≥0.68

(1)

式中：k為極值點的個數(shù)；xσ為極值點序列的統(tǒng)計有義值；P{·}為概率.式(1)的數(shù)學意義為“全部極值點的絕對值小于統(tǒng)計有義值的概率不小于0.68”.

(2)

(3)

步驟3輸出時序s1(t)為1階IMF分量，并判斷下列終止條件是否成立.

(4)

式中：σ0和σ1分別為時間序列s(t)和s1(t)的標準差；θ1為時間序列s1(t)的局部限制參數(shù)序列.P{θ1|θ1≤θ0}的數(shù)學意義為“θ1中各參數(shù)不大于θ0中對應(yīng)參數(shù)的統(tǒng)計概率不小于0.95”.

若終止條件式(4)成立，則直接進入步驟4，否則將s1(t)作為新的輸入時序重復(fù)步驟1～2，計算得到2階IMF分量；如此反復(fù)，直至滿足終止條件式(4)，記錄此時IMF分量的最大階數(shù)為n，并輸出各階IMF分量sg(t)(g=1,2,…,n).

步驟4計算第n階殘余分量rn(t)=sn-1(t)-sn(t)，若rn(t)為單調(diào)函數(shù)或常數(shù)，則結(jié)束分解過程；否則將rn(t)作為新的輸入時序重復(fù)步驟1～3，最終得到的分解結(jié)果為

(5)

將每一個IMF分量看作一個行向量，得到IMF分量矩陣如下：

SN×(n+1)=

(6)

1.2 SVD分解

由于趨勢項具有低復(fù)雜度和低頻等特點，舍去最為復(fù)雜的高頻1階IMF分量，得到矩陣SN×n=[s2(t)Ts3(t)T…sn(t)Trn(t)T]，計算SN×n的協(xié)方差矩陣Cn×n：

Cn×n=

(7)

式中：E(·)為取均值計算.E(SN×n)為矩陣SN×n的均值矩陣，形式如下：

E(SN×n)=

對矩陣Cn×n進行SVD分解[11]，重構(gòu)后的奇異值分量矩陣QN×K為

(8)

式中：q1,q2,…,qK為分解得到的奇異值分量；U為奇異矩陣；K為非零奇異值個數(shù).

1.3 基于OAPPE的趨勢項重構(gòu)

排列熵反映了時序的復(fù)雜程度，根據(jù)排列熵的大小可以對時間序列的所有極值點進行排序，從分解得到的奇異值分量矩陣QN×K中得到時間序列的預(yù)測趨勢項，具體步驟如下.

步驟1對所計算的極值點時間序列X={x1,x2,…,xk}進行相空間重構(gòu).

(9)

式中：m為嵌入維數(shù)；τ為延遲時間.

步驟2將每一行升序排列獲得元素初始位置的索引值序列，計算m!種排列方式在索引值序列中出現(xiàn)的概率Pz，并由此計算排列熵Hpe:

(10)

步驟3根據(jù)文獻[12]提出的互信息算法確定延遲時間τ.首先，確定分段數(shù)d=1.87(n-τ-1)0.4，將一個二維坐標軸空間均勻劃分成d2個網(wǎng)格，并將每個網(wǎng)格記為Rf(f=1,2,…,d2)；然后，取時序X的前k-τ項{x1,x2,…,xk-τ}和后k-τ項{xτ+1,xτ+2,…,xk}，構(gòu)成A、B兩個與延遲時間τ相關(guān)的時序空間，計算兩個空間的互信息值I(τ)，則有：

(11)

式中：PA(Rf)和PB(Rf)為時序空間A和B中的點全部落在網(wǎng)格Rf中的概率；PAB(Rf)為時序空間A和B中的點同時落在網(wǎng)格Rf中的概率.

互信息值I(τ)實際上是聯(lián)合分布PAB(Rf)與邊緣分布PA(Rf)、PB(Rf)的相對熵，反映了兩個空間系統(tǒng)的統(tǒng)計相關(guān)性，互信息值越小,兩個系統(tǒng)間的相關(guān)性就越小，本文取I(τ)的第1個極小值點對應(yīng)的時間τ為最佳延遲時間.

步驟4在獲得最佳延遲時間的基礎(chǔ)上，采用偽近鄰法確定最小嵌入維數(shù).對于時序X={x1,x2,…,xk}，當嵌入維數(shù)為m時，重構(gòu)時序X的集合Ω可以表示為

Ω={Xi|Xi={xi,xi+τ,…,xi+(m-1)τ}；

i=1,2,…,n-(m-1)τ}

(12)

(13)

(14)

式中：ath和bth為根據(jù)計算精度要求自定的閾值.

將集合Ω中所有元素的近鄰點計算出來，組成集合Ψ可以表示為

i=1,2,…,n-(m-1)τ}

(15)

根據(jù)式(14)篩選集合Ψ中所有偽近鄰點，并計算偽近鄰率(FNNP)FΨ，可以表示為

FΨ=NΨ/k

(16)

式中：NΨ近鄰點集合Ψ中偽近鄰點的個數(shù).

當偽近鄰率小于5%或者偽近鄰率不再隨著嵌入維數(shù)m的增大而減小時，則該點對應(yīng)的嵌入維數(shù)即為最小嵌入維數(shù)[13].

步驟5對于奇異值分量矩陣QN×K中所有的奇異值分量，按照步驟1～4計算其排列熵，然后采用K-means++算法[14]對計算得到的排列熵進行聚類分析，選擇排列熵值小的奇異值分量進行疊加重構(gòu)，獲得時間序列的預(yù)測趨勢項.

1.4 基于ARIMA模型的失穩(wěn)概率計算

根據(jù)所獲得的預(yù)測趨勢項，預(yù)測運行參數(shù)時間序列在未來一小段時間內(nèi)的趨勢走向，計算這段時間內(nèi)參數(shù)超出穩(wěn)定范圍的失穩(wěn)概率，并以此作為其運行穩(wěn)定性的評估指標.考慮到復(fù)雜熱力系統(tǒng)的運行參數(shù)都具有趨勢項非平穩(wěn)的特點，所以選擇在非平穩(wěn)時間序列分析中較為常用的ARIMA模型進行失穩(wěn)概率計算.ARIMA模型通過多次差分將非平穩(wěn)時間序列變?yōu)槠椒€(wěn)時間序列，在計算時有3個關(guān)鍵參數(shù)：自回歸項數(shù)、滑動平均項數(shù)和差分次數(shù)，其計算流程和計算方法在文獻[9]中已有詳細描述，本節(jié)不再贅述.

對于預(yù)測趨勢線上任一時間點的參數(shù)y，采用正態(tài)分布模型計算其失穩(wěn)概率Py：

Py=P{y≥ymax∪y≤ymin}

(17)

1.5 穩(wěn)定性評估流程圖

綜上所述，基于MOSP的蒸汽動力系統(tǒng)單參數(shù)運行穩(wěn)定性評估流程，如圖1所示.

圖1 蒸汽動力系統(tǒng)單參數(shù)運行穩(wěn)定性評估流程圖Fig.1 Flow chart of single parameter operation stability assessment of steam power system

對于輸入的時間序列依次進行MREMD分解與SVD分解；通過最佳算法參數(shù)獨立確定方法計算奇異值分量矩陣的排列熵；以排列熵為依據(jù)選取最低排列熵值的奇異值分量進行疊加重構(gòu)，獲得原時間序列的趨勢項與擾動項；通過ARIMA模型預(yù)測未來一段時間內(nèi)參數(shù)的變化趨勢，并計算預(yù)測趨勢上各點的失穩(wěn)概率.

2 蒸汽動力系統(tǒng)單參數(shù)運行穩(wěn)定性評估

為了驗證所提方法的有效性，取某型蒸汽動力系統(tǒng)在一段時間內(nèi)的微過熱蒸汽壓力為驗證數(shù)據(jù)，在MATLAB 2016環(huán)境下編寫代碼完成評估全過程，各步驟的運行結(jié)果如下.

(1)對輸入的運行參數(shù)(微過熱蒸汽壓力)進行MREMD分解，獲得的7個IMF分量(s1～s7)和殘余分量(re)如圖2所示，其中：t為時間.

圖2 由MREMD分解得到的IMF分量Fig.2 IMF components obtained by MREMD decomposition

(2)在去除1階IMF分量后繼續(xù)進行SVD分解，獲得的7個奇異值分量(q1～q7)如圖3所示.

圖3 由SVD分解得到奇異值分量Fig.3 Singular value components obtained by SVD decomposition

(3)為計算各奇異值分量的排列熵，通過互信息法確定最佳延遲時間，計算當τ∈[0,100]s時，各階奇異值分量的互信息值(I1～I7)隨τ的變化情況,如圖4所示.

圖4 互信息值隨延遲時間的變化Fig.4 Mutual information versus delay time

(4)以最佳延遲時間為輸入，計算當m∈[0,100]時，各階奇異值分量的偽近鄰率(F1～F7)隨m的變化情況，如圖5所示.

圖5 偽近鄰率隨嵌入維數(shù)的變化Fig.5 False nearest neighbor percents versus embedding dimension

(5)根據(jù)1.3節(jié)給出的依據(jù)選取最小嵌入維數(shù)，并計算各奇異值分量的排列熵，計算結(jié)果如表1所示.

表1 各奇異值分量的最優(yōu)排列熵算法參數(shù)及排列熵Tab.1 Optimal permutation entropy algorithm parameters and permutation entropy of each singular value component

(6)根據(jù)排列熵對各奇異值分量進行K-means聚類，并選擇排列熵值小的分量進行重組.根據(jù)分類結(jié)果，選取q4、q5、q6和q7進行疊加重構(gòu)，得到的微過熱蒸汽壓力的預(yù)測趨勢項與擾動項如圖6所示.其中：pst為微過熱蒸汽的壓力值；psr為微過熱蒸汽壓力的擾動值.

圖6 微過熱蒸汽壓力數(shù)據(jù)的預(yù)測趨勢項與擾動項Fig.6 Predicted trends and disturbance items of slightly superheated steam pressure data

(7)選擇自回歸項數(shù)為3、滑動平均項數(shù)為2、差分次數(shù)為1，建立ARIMA(3,1,2)模型對微過熱蒸汽壓力進行預(yù)測.為檢驗?zāi)Ｐ偷念A(yù)測效果，取圖6(a)中原信號數(shù)據(jù)的前80%作為訓練集，后20%作為驗證集，對模型預(yù)測能力進行檢驗，對比結(jié)果如圖7所示.

圖7 ARIMA模型預(yù)測能力檢驗Fig.7 Predictive ability test of ARIMA model

由圖7可見，模型輸出的預(yù)測趨勢項均落在驗證數(shù)據(jù)的95%置信區(qū)間內(nèi)，說明模型參數(shù)選擇合理、預(yù)測精度較高.

(8)該型蒸汽動力系統(tǒng)的微過熱蒸汽壓力整定值為2.45 MPa，取其±3%的波動區(qū)間為穩(wěn)定區(qū)間，應(yīng)用ARIMA(3,1,2)模型進行預(yù)測，計算得到的微過熱蒸汽壓力在未來60 s的變化趨勢及失穩(wěn)概率如圖8所示，其中：P′為失穩(wěn)概率.

由圖8可知，微過熱蒸汽壓力在未來60 s內(nèi)的失穩(wěn)概率較小，均在0.2以下，說明其運行穩(wěn)定性較好，失穩(wěn)風險較低.

3 結(jié)語

針對目前蒸汽動力系統(tǒng)缺少定量運行穩(wěn)定性評估方法的問題，提出一種基于MSOP的單參數(shù)穩(wěn)定性評估方法，并通過對某型蒸汽動力系統(tǒng)微過熱蒸汽壓力的評估分析，驗證了該方法的科學性.本文的研究成果具有一定的理論創(chuàng)新性和工程應(yīng)用價值，通過增加數(shù)據(jù)接口讀取實際系統(tǒng)的實時運行數(shù)據(jù)，實現(xiàn)蒸汽動力系統(tǒng)的在線運行穩(wěn)定性評估，為操作人員的使用管理提供參考.同時,本研究也可為蒸汽動力系統(tǒng)的全系統(tǒng)運行穩(wěn)定性提供底層評估指標，通過建立蒸汽動力系統(tǒng)運行穩(wěn)定性評估指標體系，實現(xiàn)整個系統(tǒng)的運行穩(wěn)定性評估.