鄭奕揚(yáng)，倪 何，金家善

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

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

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

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

1 基于MSOP的運(yùn)行穩(wěn)定性評(píng)估方法

1.1 MREMD分解

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

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

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

(1)

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

(2)

(3)

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

(4)

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

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

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

(5)

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

SN×(n+1)=

(6)

1.2 SVD分解

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

Cn×n=

(7)

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

E(SN×n)=

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

(8)

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

1.3 基于OAPPE的趨勢(shì)項(xiàng)重構(gòu)

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

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

(9)

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

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

(10)

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

(11)

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

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

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

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

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

(12)

(13)

(14)

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

將集合Ω中所有元素的近鄰點(diǎn)計(jì)算出來(lái)，組成集合Ψ可以表示為

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

(15)

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

FΨ=NΨ/k

(16)

式中：NΨ近鄰點(diǎn)集合Ψ中偽近鄰點(diǎn)的個(gè)數(shù).

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

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

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

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

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

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

(17)

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

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

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

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

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

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

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

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

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

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

(3)為計(jì)算各奇異值分量的排列熵，通過(guò)互信息法確定最佳延遲時(shí)間，計(jì)算當(dāng)τ∈[0,100]s時(shí)，各階奇異值分量的互信息值(I1～I(xiàn)7)隨τ的變化情況,如圖4所示.

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

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

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

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

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

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

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

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

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

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

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

由圖8可知，微過(guò)熱蒸汽壓力在未來(lái)60 s內(nèi)的失穩(wěn)概率較小，均在0.2以下，說(shuō)明其運(yùn)行穩(wěn)定性較好，失穩(wěn)風(fēng)險(xiǎn)較低.

3 結(jié)語(yǔ)

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