基于MSOP的蒸汽動力系統單參數運行穩定性評估方法

鄭奕揚，倪 何，金家善

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

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

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

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

1 基于MSOP的運行穩定性評估方法

1.1 MREMD分解

MREMD[10]是經驗模態分解的改進算法，通過優化均值包絡的生成方式來抑制端點效應.對于輸入長度為N的時間序列s(t)，MREMD分解的具體步驟如下.

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

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

(1)

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

(2)

(3)

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

(4)

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

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

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

(5)

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

SN×(n+1)=

(6)

1.2 SVD分解

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

Cn×n=

(7)

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

E(SN×n)=

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

(8)

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

1.3 基于OAPPE的趨勢項重構

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

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

(9)

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

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

(10)

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

(11)

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

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

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

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

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

(12)

(13)

(14)

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

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

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

(15)

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

FΨ=NΨ/k

(16)

式中：NΨ近鄰點集合Ψ中偽近鄰點的個數.

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

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

1.4 基于ARIMA模型的失穩概率計算

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

對于預測趨勢線上任一時間點的參數y，采用正態分布模型計算其失穩概率Py：

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

(17)

1.5 穩定性評估流程圖

綜上所述，基于MOSP的蒸汽動力系統單參數運行穩定性評估流程，如圖1所示.

圖1 蒸汽動力系統單參數運行穩定性評估流程圖Fig.1 Flow chart of single parameter operation stability assessment of steam power system

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

2 蒸汽動力系統單參數運行穩定性評估

為了驗證所提方法的有效性，取某型蒸汽動力系統在一段時間內的微過熱蒸汽壓力為驗證數據，在MATLAB 2016環境下編寫代碼完成評估全過程，各步驟的運行結果如下.

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

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

(2)在去除1階IMF分量后繼續進行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 偽近鄰率隨嵌入維數的變化Fig.5 False nearest neighbor percents versus embedding dimension

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

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

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

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

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

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

由圖7可見，模型輸出的預測趨勢項均落在驗證數據的95%置信區間內，說明模型參數選擇合理、預測精度較高.

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

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

3 結語

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