劉艷麗，余貽鑫

(天津大學(xué)智能電網(wǎng)教育部重點(diǎn)實(shí)驗(yàn)室，天津 300072)

馬爾可夫模型是可修復(fù)系統(tǒng)可靠性和可用性分析中常用的工具[1-2]．對(duì)于電網(wǎng)這一典型大系統(tǒng)，馬爾可夫模型還被應(yīng)用到連鎖故障預(yù)測(cè)、風(fēng)電預(yù)測(cè)、概率安全性評(píng)估和風(fēng)險(xiǎn)評(píng)估中[3-6]．概率安全性評(píng)估、風(fēng)險(xiǎn)評(píng)估等諸多研究中應(yīng)用馬爾可夫模型時(shí)，主要用其求解與時(shí)間相關(guān)的狀態(tài)概率，而求解的關(guān)鍵在于狀態(tài)轉(zhuǎn)移概率矩陣的計(jì)算．對(duì)于由 N 個(gè) m 狀態(tài)元件組成的系統(tǒng)，系統(tǒng)狀態(tài)數(shù)為 mN．對(duì)于小系統(tǒng)而言，系統(tǒng)結(jié)構(gòu)狀態(tài)數(shù)目少，可直觀確定狀態(tài)轉(zhuǎn)移概率矩陣的構(gòu)成和非零元素的位置，非零元素的數(shù)值即狀態(tài)轉(zhuǎn)移率也容易獲得．然而，大系統(tǒng)有巨大數(shù)量的結(jié)構(gòu)狀態(tài)，不能直觀確定狀態(tài)轉(zhuǎn)移概率矩陣的構(gòu)成和非零元素的位置，且若元件狀態(tài)數(shù)m 較大時(shí)，很難快速準(zhǔn)確地獲得狀態(tài)轉(zhuǎn)移率．大系統(tǒng)狀態(tài)轉(zhuǎn)移概率矩陣的形成是十分復(fù)雜、繁重且耗時(shí)的工作，這大大限制了基于馬爾可夫模型的概率安全性評(píng)估[5]等方法在實(shí)際工程中的應(yīng)用．可見，研究大系統(tǒng)狀態(tài)轉(zhuǎn)移概率矩陣的快速形成方法具有重要的工程應(yīng)用價(jià)值．

狀態(tài)轉(zhuǎn)移概率矩陣形成方法的相關(guān)研究很少．通常情況下，均是基于狀態(tài)空間圖獲得狀態(tài)轉(zhuǎn)移概率矩陣[7]，其優(yōu)點(diǎn)是狀態(tài)及其相互轉(zhuǎn)移有清晰的圖形表示，但建立大系統(tǒng)的狀態(tài)空間圖是不可能的，故該方法對(duì)大系統(tǒng)的應(yīng)用十分困難．Amoia 等[8-9]提出了面向計(jì)算機(jī)的基于Kronecker 代數(shù)形成大系統(tǒng)狀態(tài)轉(zhuǎn)移概率矩陣的方法，但其計(jì)算時(shí)間隨著系統(tǒng)元件的增多而指數(shù)增長(zhǎng)，無(wú)法滿足實(shí)際工程應(yīng)用的要求．

筆者提出一個(gè)快速形成馬爾可夫模型狀態(tài)轉(zhuǎn)移概率矩陣的方法．定義元件狀態(tài)轉(zhuǎn)移率矩陣和系統(tǒng)狀態(tài)數(shù)組，將系統(tǒng)狀態(tài)轉(zhuǎn)換為便于計(jì)算機(jī)存儲(chǔ)與處理的數(shù)組，有效地描述了系統(tǒng)狀態(tài)之間的轉(zhuǎn)移；基于元件狀態(tài)轉(zhuǎn)移率矩陣和系統(tǒng)狀態(tài)數(shù)組提出不受系統(tǒng)狀態(tài)和元件狀態(tài)數(shù)目限制快速計(jì)算狀態(tài)轉(zhuǎn)移率的方法，挖掘矩陣中非零元素的分布規(guī)律提出非零元素的快速定位方法，進(jìn)而快速形成狀態(tài)轉(zhuǎn)移概率矩陣的稀疏存儲(chǔ)；針對(duì)兩狀態(tài)元件組成的系統(tǒng)，提出基于給定的系統(tǒng)狀態(tài)排序和服務(wù)狀態(tài)集數(shù)組快速定位矩陣中非零元素的方法．對(duì)電力系統(tǒng)概率安全性評(píng)估的應(yīng)用效果顯著，證實(shí)了方法的實(shí)用性和有效性．

1 元件狀態(tài)轉(zhuǎn)移率矩陣和系統(tǒng)狀態(tài)數(shù)組

1.1 元件狀態(tài)轉(zhuǎn)移率矩陣

定義元件狀態(tài)轉(zhuǎn)移率矩陣T，即若某元件在系統(tǒng)中有m 個(gè)狀態(tài)，該元件的狀態(tài)轉(zhuǎn)移率矩陣為m × m維矩陣，第i 行、第j 列對(duì)應(yīng)的元素為元件從元件狀態(tài)i向元件狀態(tài)j 的轉(zhuǎn)移率．

為方便理解，均以兩狀態(tài)元件(即 m = 2)為例說明相關(guān)的概念和方法．元件有運(yùn)行和停運(yùn)兩個(gè)狀態(tài)，0 表示運(yùn)行，1 表示停運(yùn)．設(shè)元件事故率為λ，修復(fù)率為μ，則元件狀態(tài)轉(zhuǎn)移率矩陣T 為

1.2 系統(tǒng)狀態(tài)數(shù)組

若第i 個(gè)元件有mi個(gè)狀態(tài)，用0，1，…，mi-1 表示，各元件狀態(tài)可重復(fù)組合[10]數(shù)為可知系統(tǒng)共有維狀態(tài)數(shù)組，以此可快速生成所有系統(tǒng)狀態(tài)．

以N 個(gè)兩狀態(tài)元件組成的系統(tǒng)為例，如圖1 所示，將各元件狀態(tài)可重復(fù)組合可得2N個(gè)系統(tǒng)狀態(tài)數(shù)組，這些系統(tǒng)狀態(tài)數(shù)組對(duì)應(yīng)2N個(gè)系統(tǒng)狀態(tài)．

圖1 系統(tǒng)狀態(tài)數(shù)組Fig.1 System configuration array

可見，通過系統(tǒng)狀態(tài)數(shù)組可將系統(tǒng)狀態(tài)轉(zhuǎn)換為便于計(jì)算機(jī)存儲(chǔ)與處理的數(shù)組并快速生成所有系統(tǒng)狀態(tài)．

1.3 系統(tǒng)狀態(tài)轉(zhuǎn)移的描述

基于元件狀態(tài)轉(zhuǎn)移率矩陣和系統(tǒng)狀態(tài)數(shù)組可有效描述系統(tǒng)狀態(tài)之間的轉(zhuǎn)移．以可修復(fù)兩元件系統(tǒng)為例進(jìn)行說明，元件有運(yùn)行和停運(yùn)兩個(gè)狀態(tài)，元件1的故障率和修復(fù)率分別為λ1和μ1；元件2 的故障率和修復(fù)率分別為λ2和μ2．

由第1.1 節(jié)得元件1 和元件2 的狀態(tài)轉(zhuǎn)移率矩陣T1和T2為

由第1.2 節(jié)可知，可用一個(gè)1× 2 維狀態(tài)數(shù)組表示系統(tǒng)狀態(tài)，數(shù)組中第1 個(gè)元素代表元件1 的狀態(tài)，第2 個(gè)元素代表元件2 的狀態(tài)，將各元件狀態(tài)進(jìn)行可重復(fù)組合生成4 個(gè)狀態(tài)數(shù)組[0 0]、[1 0]、[0 1]和[1 1]，對(duì)應(yīng)元件全部運(yùn)行、元件1 停運(yùn)元件2 運(yùn)行、元件1運(yùn)行元件2 停運(yùn)和元件全部停運(yùn)等4 個(gè)系統(tǒng)狀態(tài)．

基于元件狀態(tài)轉(zhuǎn)移率矩陣和系統(tǒng)狀態(tài)數(shù)組可快速獲得狀態(tài)之間的轉(zhuǎn)換關(guān)系，有：[0 0]→[1 0]元件1從狀態(tài)0 轉(zhuǎn)移到狀態(tài)1，即元件1 停運(yùn)，由T1可知狀態(tài)轉(zhuǎn)移率為λ1；[1 0]→[0 0]元件1 從狀態(tài)1 轉(zhuǎn)移到狀態(tài)0，即元件1 恢復(fù)運(yùn)行，由T1可知狀態(tài)轉(zhuǎn)移率為μ1；[0 0]→[0 1]元件2 從狀態(tài)0 轉(zhuǎn)移到狀態(tài)1，即元件2 停運(yùn)，由T2可知狀態(tài)轉(zhuǎn)移率為λ2；[0 1]→[0 0]元件2 從狀態(tài)1 轉(zhuǎn)移到狀態(tài)0，即元件2 恢復(fù)運(yùn)行，由T2可知狀態(tài)轉(zhuǎn)移率為μ2；[1 0]→[1 1]元件2 從狀態(tài)0 轉(zhuǎn)移到狀態(tài)1，即元件2 停運(yùn)，由T2可知狀態(tài)轉(zhuǎn)移率為λ2；[0 1]→[1 1]元件1 從狀態(tài)0 轉(zhuǎn)移到狀態(tài)1，即元件1 停運(yùn)，由T1可知狀態(tài)轉(zhuǎn)移率為λ1；[1 1]→[0 1]元件1 從狀態(tài)1 轉(zhuǎn)移到狀態(tài)0，即元件1 恢復(fù)運(yùn)行，由T1可知狀態(tài)轉(zhuǎn)移率為μ1；[1 1]→[1 0]元件2 從狀態(tài)1轉(zhuǎn)移到狀態(tài)0，即元件2 恢復(fù)運(yùn)行，由T2可知狀態(tài)轉(zhuǎn)移率為μ2；[0 0]和[1 1]之間無(wú)轉(zhuǎn)移，狀態(tài)轉(zhuǎn)移率為0．

與狀態(tài)空間圖相比，元件狀態(tài)轉(zhuǎn)移率矩陣和系統(tǒng)狀態(tài)數(shù)組有效地描述了狀態(tài)間的轉(zhuǎn)換，便于計(jì)算機(jī)處理，更適用于大系統(tǒng)．

2 狀態(tài)轉(zhuǎn)移概率矩陣的快速計(jì)算方法

2.1 基本假設(shè)

取與文獻(xiàn)[9，11]相同的假設(shè)如下：

(1) 系統(tǒng)的狀態(tài)取決于各元件的狀態(tài)，這與系統(tǒng)狀態(tài)數(shù)組的定義剛好吻合；

(2) 鑒于實(shí)際工程應(yīng)用中，常采用離散時(shí)間的馬爾可夫模型分析電網(wǎng)等大規(guī)模系統(tǒng)，所以本文研究的對(duì)象是離散時(shí)間的馬爾可夫模型；

(3) 當(dāng)i≠j 且狀態(tài)i 和j 中有2 個(gè)或2 個(gè)以上元件的狀態(tài)相異時(shí)，狀態(tài)i 到狀態(tài)j 的狀態(tài)轉(zhuǎn)移率為0．

2.2 狀態(tài)轉(zhuǎn)移率的快速準(zhǔn)確計(jì)算

系統(tǒng)狀態(tài)i 和j 對(duì)應(yīng)的狀態(tài)數(shù)組為 Ai和 Aj，可快速計(jì)算狀態(tài)轉(zhuǎn)移率．

(1) Ai和Aj有2 個(gè)或2 個(gè)以上元素相異時(shí)，pij＝0；

(2) Ai和Aj有1 個(gè)元素相異時(shí)，計(jì)算狀態(tài)轉(zhuǎn)移率pij的公式為

式中：Tk表示元件k 的狀態(tài)轉(zhuǎn)移率矩陣；Ai(k)、Aj(k)表示Ai和Aj中的第k 個(gè)元素，即系統(tǒng)狀態(tài)i 和j 中元件k 的狀態(tài)；TkAi(k)Aj(k)表示元件k 從狀態(tài)Ai(k)到Aj(k)的狀態(tài)轉(zhuǎn)移率；N 為系統(tǒng)中元件的個(gè)數(shù)．

(3) Ai和Aj有零個(gè)元素相異時(shí)，即狀態(tài)i 向自身發(fā)生轉(zhuǎn)移時(shí)，則計(jì)算狀態(tài)轉(zhuǎn)移率pij的公式為

該方法不受系統(tǒng)狀態(tài)和元件狀態(tài)數(shù)目的限制，無(wú)需狀態(tài)空間圖便可快速準(zhǔn)確地計(jì)算狀態(tài)轉(zhuǎn)移率．

2.3 狀態(tài)轉(zhuǎn)移概率矩陣中非零元素的分布規(guī)律

系統(tǒng)狀態(tài)轉(zhuǎn)移概率矩陣中有大量的零元素，若能找到矩陣中非零元素的分布規(guī)律，將大大有助于實(shí)現(xiàn)非零元素的快速定位，進(jìn)而實(shí)現(xiàn)矩陣的稀疏存儲(chǔ)．

基于第2.2 節(jié)給出的系統(tǒng)狀態(tài)轉(zhuǎn)移率的計(jì)算方法以及諸多實(shí)例系統(tǒng)驗(yàn)證，筆者發(fā)現(xiàn)系統(tǒng)狀態(tài)轉(zhuǎn)移概率矩陣中的非零元素存在如下分布規(guī)律．

若系統(tǒng)中有N 個(gè)元件，元件有m 個(gè)狀態(tài)，系統(tǒng)狀態(tài)i 對(duì)應(yīng)的狀態(tài)數(shù)組為Ai，更改Ai某個(gè)元素形成的狀態(tài)數(shù)組為Ai'，則

(1) Ai'對(duì)應(yīng)的系統(tǒng)狀態(tài)是i 可轉(zhuǎn)移的系統(tǒng)狀態(tài)；

(2) i 可向N×(m ? 1)個(gè)不同于自身的系統(tǒng)狀態(tài)發(fā)生轉(zhuǎn)移；

(3) 系統(tǒng)狀態(tài)轉(zhuǎn)移概率矩陣中，每行均有N ×(m ? 1) + 1個(gè)非零元素；

(4) 特別地，若系統(tǒng)中元件僅有運(yùn)行和停運(yùn)2 個(gè)狀態(tài)，則任一系統(tǒng)狀態(tài)可向N 個(gè)不同于自身的系統(tǒng)狀態(tài)發(fā)生轉(zhuǎn)移；系統(tǒng)狀態(tài)轉(zhuǎn)移概率矩陣中，每行有N＋1 個(gè)非零元素．

2.4 狀態(tài)轉(zhuǎn)移概率矩陣中非零元素的快速定位

基于第2.3 節(jié)所述狀態(tài)轉(zhuǎn)移概率矩陣中非零元素的分布規(guī)律得出非零元素的快速定位方法如下．

(1) 系統(tǒng)中有N 個(gè)元件，每個(gè)元件有m 個(gè)狀態(tài)，系統(tǒng)狀態(tài)集M 中系統(tǒng)狀態(tài)的個(gè)數(shù)為NS，從1 至NS對(duì)M 中的狀態(tài)進(jìn)行編號(hào)；

(2) ?i ∈ M，依次改變Ai的N 個(gè)元素，確定可向i 演變的 N ×(m ? 1)個(gè)不同于i 的系統(tǒng)狀態(tài)，這些狀態(tài)構(gòu)成的集合為L(zhǎng)i；

(3) M 中可向i 演變的系統(tǒng)狀態(tài)構(gòu)成的集合為Mi＝{M∩Li，i}；

(4) 統(tǒng)計(jì)Mi中系統(tǒng)狀態(tài)的個(gè)數(shù)為Ni，Ni即為系統(tǒng)狀態(tài)轉(zhuǎn)移概率矩陣第i 行中非零元素的個(gè)數(shù)；

(5) ?j ∈ Mi，系統(tǒng)狀態(tài)轉(zhuǎn)移概率矩陣第i 行、第j 列即為非零元素的位置．

2.5 狀態(tài)轉(zhuǎn)移概率矩陣的計(jì)算流程

狀態(tài)轉(zhuǎn)移概率矩陣的計(jì)算流程如圖2 所示．

圖2 狀態(tài)轉(zhuǎn)移概率矩陣的計(jì)算流程Fig.2 Flow chart of state transition probability matrix

以9 節(jié)點(diǎn)電力系統(tǒng)(如圖3 所示)為例說明上述狀態(tài)轉(zhuǎn)移概率矩陣的計(jì)算流程．

圖3 9節(jié)點(diǎn)系統(tǒng)Fig.3 9-bus system

步驟1確定系統(tǒng)狀態(tài)集和元件狀態(tài)轉(zhuǎn)移率矩陣．

考慮表1 所示的6 條線路，每條線路有運(yùn)行和停運(yùn)2 個(gè)狀態(tài)，分別用0 和1 表示．將各元件狀態(tài)可重復(fù)組合生成26＝64 個(gè)系統(tǒng)狀態(tài)數(shù)組，與之對(duì)應(yīng)的系統(tǒng)狀態(tài)構(gòu)成系統(tǒng)狀態(tài)集．這些狀態(tài)及其對(duì)應(yīng)的狀態(tài)數(shù)組中存儲(chǔ)的元件狀態(tài)如表2 所示．元件的事故率和修復(fù)率見表1，元件i 的狀態(tài)轉(zhuǎn)移率矩陣為

步驟2對(duì)狀態(tài)集中的狀態(tài)進(jìn)行編號(hào)．對(duì)狀態(tài)集中的64 個(gè)狀態(tài)進(jìn)行編號(hào)如表2 所示．

步驟3確定每個(gè)狀態(tài)可轉(zhuǎn)移的狀態(tài)．以狀態(tài)10 為例，狀態(tài)10 對(duì)應(yīng)的狀態(tài)數(shù)組為[1 0 0 1 0 0]，依次改變狀態(tài)數(shù)組中的元素獲得6 個(gè)狀態(tài)數(shù)組：[0 0 0 1 0 0]，[1 1 0 1 0 0]，[1 0 1 1 0 0]，[1 0 0 0 0 0]，[1 0 0 1 1 0]，[1 0 0 1 0 1]．對(duì)照表2 可得，以上6 個(gè)狀態(tài)數(shù)組對(duì)應(yīng)的狀態(tài)為5、24、27、2、30、31．則狀態(tài)10 可轉(zhuǎn)移的狀態(tài)有 7 個(gè)，它們是(按狀態(tài)編號(hào)從小到大排序)2、5、10、24、27、30、31．

步驟4確定每行非零元素的位置．仍以狀態(tài)10 為例．狀態(tài)10 對(duì)應(yīng)狀態(tài)轉(zhuǎn)移概率矩陣的第10行，則第10 行中2、5、10、24、27、30、31 列存在非零元素．

步驟5根據(jù)式(2)和式(3)計(jì)算非零元素的數(shù)值．以狀態(tài)轉(zhuǎn)移概率矩陣第10 行為例，有p10,2＝μ4，p10,5＝μ1，p10,10＝1 ?(μ1＋μ4＋λ2＋λ3＋λ5＋λ6)，p10,24＝λ2，p10,27＝λ3，p10,30＝λ5，p10,31＝λ6．

步驟6形成系統(tǒng)狀態(tài)轉(zhuǎn)移概率矩陣．

表1 線路及其編號(hào)Tab.1 Line and its order

循環(huán)執(zhí)行步驟3～步驟5 形成示例系統(tǒng)狀態(tài)轉(zhuǎn)移概率矩陣的稀疏存儲(chǔ)如表3～表6 所示．采用Matlab 編程，在聯(lián)想Dual Core CPU、1.66,GHz、2,G內(nèi)存計(jì)算機(jī)上計(jì)算示例系統(tǒng)的狀態(tài)轉(zhuǎn)移概率矩陣耗時(shí)0.000,33,s，與文獻(xiàn)[9]耗時(shí)0.001,4,s 相比，本文方法計(jì)算速度更快．

需要說明的是，文獻(xiàn)[9]是現(xiàn)有大系統(tǒng)馬爾可夫模型狀態(tài)轉(zhuǎn)移概率矩陣形成方法中較好的方法，故在本節(jié)和第3 節(jié)中將其作為比較對(duì)象對(duì)本文方法進(jìn)行分析．

2.6 兩狀態(tài)元件組成的系統(tǒng)的狀態(tài)轉(zhuǎn)移概率矩陣的快速計(jì)算

鑒于兩狀態(tài)元件組成的系統(tǒng)常被作為研究對(duì)象，筆者對(duì)其做了進(jìn)一步研究，提出基于給定系統(tǒng)狀態(tài)排序和服務(wù)狀態(tài)集數(shù)組快速定位狀態(tài)轉(zhuǎn)移概率矩陣中非零元素的方法．仍以圖4 所示系統(tǒng)為例說明如下．

1) 給定系統(tǒng)狀態(tài)排序

系統(tǒng)有N 個(gè)元件，k 個(gè)元件故障對(duì)應(yīng)的狀態(tài)組成集合Bk．按下列規(guī)則排序系統(tǒng)狀態(tài)集中的狀態(tài)：

(1) 按B0至BN排序；

(2) 按如下規(guī)則排序Bk中狀態(tài)：

① k＝0 時(shí)，B0＝{1}；

② k＝1 時(shí)，將A1中第1 個(gè)零元素至最后一個(gè)零元素依次改為1，形成狀態(tài)2 至狀態(tài)N＋1；

③ k≥2 時(shí)，Bk-1中編號(hào)最小的狀態(tài)為s，狀態(tài)數(shù)為N'，Ai中非零元素后有Ni個(gè)零元素．

將As中非零元素后第一個(gè)零元素至最后一個(gè)零元素依次改為1，形成狀態(tài)s＋N'至狀態(tài)s＋N'＋Ns-1，記為狀態(tài)Ps至狀態(tài)Qs；

將As＋i中非零元素后第一個(gè)零元素至最后一個(gè)零元素依次改為1，形成狀態(tài)Ps＋i至狀態(tài)Qs＋i，Ps＋i＝Qs＋i-1＋1，Qs＋i＝Qs＋i-1＋Ns＋i，i＝1，2，…，N'-1．

示例系統(tǒng)的系統(tǒng)狀態(tài)排序見表2．

表2 狀態(tài)及其編號(hào)以及服務(wù)狀態(tài)集及其編號(hào)Tab.2 State and its order and service set and its order

表3 9節(jié)點(diǎn)系統(tǒng)的狀態(tài)轉(zhuǎn)移概率矩陣的第1行到第17行Tab.3 1st row to 17th row of state transition probability matrix of 9-bus system

表4 9節(jié)點(diǎn)系統(tǒng)的狀態(tài)轉(zhuǎn)移概率矩陣的第18行到第32行Tab.4 18th row to 32nd row of state transition probability matrix of 9-bus system

表5 9節(jié)點(diǎn)系統(tǒng)的狀態(tài)轉(zhuǎn)移概率矩陣的第33行到第48行Tab.5 33rd row to 48th row of state transition probability matrix of 9-bus system

2) 生成服務(wù)狀態(tài)集數(shù)組

服務(wù)狀態(tài)集(service set)用S 表示，第n 個(gè)服務(wù)狀態(tài)集用Sn表示．除S1外，Sn中僅有一個(gè)狀態(tài)xn，中最后一個(gè)元素為1．S 的定義為

示例系統(tǒng)的服務(wù)狀態(tài)集見表2．

定義服務(wù)狀態(tài)集數(shù)組C：C 是一個(gè)行數(shù)組，其列數(shù)為n，存儲(chǔ)的元素為xn．

3) 基于C 快速定位狀態(tài)轉(zhuǎn)移概率矩陣的非零元素

設(shè)i 可轉(zhuǎn)移的狀態(tài)構(gòu)成的狀態(tài)集為Mi＝{Sup，i，Sdown}．若Ai中有Nup個(gè)1，Ndown個(gè)0，則Sup中有Nup個(gè)狀態(tài)，由將Ai中某個(gè)1 改為0 后對(duì)應(yīng)的狀態(tài)構(gòu)成；Sdown中有Ndown個(gè)狀態(tài)，由將Ai中某個(gè)0 改為1 后對(duì)應(yīng)的狀態(tài)構(gòu)成．可基于C 直接確定Mi如下所示．

(1) i＝1 時(shí)，Nup＝0，Sup＝?；Ndown＝N，Sdown＝S2＝{2，3，…，N＋1}，Mi＝{1，2，…，N＋1}．

(2) ?i∈B1，Nup＝1，Sup＝{1}；Ndown＝N-1，i＝2 時(shí)，Sdown＝S3；i＝N＋1 時(shí)，Sdown＝{C(1，3)，C(1，4)，…，C(1，i)}；i 為2 和N＋1 以外的狀態(tài)時(shí)，Sdown＝Si＋1∪{C(1，j)-(C（1，2)-i），j＝3，4，…，i}．其中，Si＋1＝{C(1，i) ＋1，C(1，i)＋2，…，C(1，i＋1)}．

(3) ?i∈B2，Nup＝2，Ndown＝N-2．判定i 所屬服務(wù)狀態(tài)集的編號(hào)，用w 表示，即i∈Sw．則有Sup＝{w-1，C(1，2)-(C（1，w)-i）}；按w＝N＋1，i＝C(1，w)，i＝C(1，w)-1 和一般情況等4 種情況計(jì)算Sdown.

(4) ?i∈B3∪B4∪…∪BN，仍可基于C 直接確定Mi，鑒于篇幅限制，本文不再贅述．

表6 9節(jié)點(diǎn)系統(tǒng)的狀態(tài)轉(zhuǎn)移概率矩陣的第49行到第64行Tab.6 49th row to 64th row of state transition probability matrix of 9-bus system

3 方法分析和應(yīng)用實(shí)例

3.1 方法分析

用 Matlab 編程，基于聯(lián)想 Dual Core CPU、1.66,GHz、2,G 內(nèi)存計(jì)算機(jī)比較了本文方法和文獻(xiàn)[9]方法，如表7 所示．結(jié)果表明了本文方法的優(yōu)越性．

進(jìn)一步地，通過Matlab 中cftool 曲線擬合工具，基于表7 數(shù)據(jù)按 ebxy a= (y 表示計(jì)算時(shí)間，x 表示系統(tǒng)規(guī)模，a 和b 表示擬合系數(shù))進(jìn)行指數(shù)曲線擬合，擬合結(jié)果如表8 所示．?dāng)M合誤差eSSE取為擬合數(shù)據(jù)和原始數(shù)據(jù)對(duì)應(yīng)點(diǎn)的誤差的平方和．結(jié)果進(jìn)一步證實(shí)了本文方法的有效性．對(duì)應(yīng)的狀態(tài)轉(zhuǎn)移概率矩陣，便于應(yīng)用的同時(shí)，可通過并行計(jì)算進(jìn)一步加速狀態(tài)轉(zhuǎn)移概率矩陣的形成．

表7 本文方法與現(xiàn)有方法計(jì)算效率對(duì)比Tab.7 Comparison of calculating efficency of method presented with existing method

表8 本文方法與現(xiàn)有方法的擬合曲線參數(shù)Tab.8 Parameters of fitting curve of method presented with existing method

3.2 應(yīng)用實(shí)例

此外，本文方法便于編程，且無(wú)需形成全狀態(tài)空間對(duì)應(yīng)的狀態(tài)轉(zhuǎn)移概率矩陣便可形成任意狀態(tài)子集

第2 節(jié)所述方法已被成功應(yīng)用于輸電系統(tǒng)概率靜態(tài)和動(dòng)態(tài)安全性綜合評(píng)估[12]．該評(píng)估模型計(jì)及了風(fēng)電和負(fù)荷引起的節(jié)點(diǎn)注入的不確定性、系統(tǒng)靜態(tài)和動(dòng)態(tài)安全約束以及網(wǎng)絡(luò)拓?fù)錉顟B(tài)之間的演變，即事故的不確定性和元件的可修復(fù)性，通過求解一個(gè)線性向量微分方程獲得 “到不安全時(shí)間” 的概率分布以評(píng)估系統(tǒng)的安全性，模型詳見文獻(xiàn)[12]．然而，由于模型求解方面存在無(wú)法快速形成微分方程系數(shù)矩陣的稀疏存儲(chǔ)等問題，致使模型無(wú)法應(yīng)用于復(fù)雜大系統(tǒng)．基于本文方法快速形成了微分方程系數(shù)矩陣的稀疏存儲(chǔ)，推進(jìn)了該模型的實(shí)用化．

以新英格蘭10 機(jī)39 節(jié)點(diǎn)系統(tǒng)為例(如圖4 所示)，考慮變壓器線路之外的34 條線路，線路有運(yùn)行和停運(yùn)兩個(gè)狀態(tài)，用于安全性評(píng)估的狀態(tài)集由B0、B1、B2、B3和B4構(gòu)成，總計(jì)52,956 個(gè)狀態(tài)．基于本文方法快速形成微分方程系數(shù)矩陣的稀疏存儲(chǔ)有：微分方程系數(shù)矩陣每行有35 個(gè)非零元素，基于服務(wù)狀態(tài)集數(shù)組直接確定非零元素的位置，非零元素的計(jì)算詳見文獻(xiàn)[12]．鑒于篇幅限制，未列出微分方程系數(shù)矩陣的計(jì)算結(jié)果．采用Matlab 編程，在聯(lián)想Dual Core CPU、1.66,GHz、2,G 內(nèi)存計(jì)算機(jī)上進(jìn)行概率靜態(tài)和動(dòng)態(tài)安全性綜合評(píng)估耗時(shí)1.1,s．

圖4 新英格蘭10機(jī)39節(jié)點(diǎn)系統(tǒng)Fig.4 New England 10 generators and 39-bus system

4 結(jié) 語(yǔ)

為了推進(jìn)基于馬爾可夫模型的概率安全性評(píng)估等方法在實(shí)際工程中的應(yīng)用，對(duì)大系統(tǒng)馬爾可夫模型狀態(tài)轉(zhuǎn)移概率矩陣的快速形成方法進(jìn)行了研究．通過狀態(tài)轉(zhuǎn)移概率矩陣和系統(tǒng)狀態(tài)數(shù)組對(duì)系統(tǒng)狀態(tài)及其之間的轉(zhuǎn)移進(jìn)行了有效描述，將系統(tǒng)狀態(tài)轉(zhuǎn)換為便于計(jì)算機(jī)存儲(chǔ)與處理的數(shù)組；提出快速準(zhǔn)確計(jì)算狀態(tài)轉(zhuǎn)移率以及快速定位狀態(tài)轉(zhuǎn)移概率矩陣中非零元素的方法，從而快速形成了狀態(tài)轉(zhuǎn)移概率矩陣的稀疏存儲(chǔ)．該方法計(jì)算準(zhǔn)確、耗時(shí)短、便于編程，使自動(dòng)快速形成大系統(tǒng)狀態(tài)轉(zhuǎn)移概率矩陣的稀疏存儲(chǔ)成為可能，為運(yùn)用馬爾可夫模型分析電網(wǎng)等大規(guī)模系統(tǒng)提供了有利工具．

［1］郭永基. 可靠性工程原理[M]. 北京：清華大學(xué)出版社，2002.Guo Yongji.Principle of Reliability Engineering[M].Beijing：Tsinghua University Press，2002(in Chinese).

［2］Shunji O，Toshio N. Bibliography for reliability and availability of stochastic systems[J].IEEE Transactions on Reliability，1976，25(4)：284-287.

［3］吳文可，文福拴，薛禹勝，等. 基于馬爾科夫鏈的電力系統(tǒng)連鎖故障預(yù)測(cè)[J]. 電力系統(tǒng)自動(dòng)化，2013，37(5)：29-37.Wu Wenke，Wen Fushuan，Xue Yusheng，et al. A Markov chain based model for forecasting power system cascading failures[J].Automation of Electric Power Systems，2013，37(5)：29-37(in Chinese).

［4］Carpinane A，Langella R，Giorgio M，et al. Very short-term probabilistic wind power forecasting based on Markov chain models[C] //IEEE International Conference on Probabilistic Methods Applied to Power Systems.Singapore，2010：107-112.

［5］Wu Felix，Tsai Y K. Probabilistic dynamic security assessment of power systems (Part I)：Basic model[J].IEEE Transactions on Circuits and Systems，1983，30(3)：148-159.

［6］李文沅. 電力系統(tǒng)風(fēng)險(xiǎn)評(píng)估-模型、方法和應(yīng)用[M].北京：科學(xué)技術(shù)出版社，2006.Li Wenyuan.Risk Assessment of Power Systems-Models，Methods，and Applications[M]. Beijing ：Science and Technology Press，2006(in Chinese).

［7］張雪松，王 超，程曉東. 基于馬爾可夫狀態(tài)空間法的超高壓電網(wǎng)繼電保護(hù)系統(tǒng)可靠性分析模型[J]. 電網(wǎng)技術(shù)，2008，32(13)：94-99.Zhang Xuesong，Wang Chao，Cheng Xiaodong. Reliability analysis model for protective relaying system of UHV power network based on Markov state space method [J].Power System Technology， 2008 ， 32(13)：94-99(in Chinese).

［8］Amoia V，Santomauro M. Computer oriented reliability analysis of large electrical system[C] //Summer Symposium on Circuit Theory. Praga，Czech，1977：317-325.

［9］Amoia V，Micheli G，Santomauro M. Computer-oriented formulation of transition probability matrices via Kronecker algebra [J].IEEE Transactions on Reliability，1981，30(2)：123-132.

［10］焜劉嘉 ，王家生. 應(yīng)用概率統(tǒng)計(jì)[M]. 北京：科學(xué)出版社，2004.Liu Jiakun，Wang Jiasheng.Applied Probability and Statistics[M]. Beijing：Science Press，2004(in Chinese).

［11］Ioannis A，Elias P. Markov processes for reliability analyses of large systems[J].IEEE Transactions on Reliability，1977，26(3)：232-237.

［12］Liu Yanli，Yu Yixin. Probabilistic steady-state and dynamic security assessment of power transmission system[J].Sci China Tech Sci，2013，56(5)：1-10.