刘艳丽，余贻鑫

(天津大学智能电网教育部重点实验室，天津 300072)

马尔可夫模型是可修复系统可靠性和可用性分析中常用的工具[1-2]．对于电网这一典型大系统，马尔可夫模型还被应用到连锁故障预测、风电预测、概率安全性评估和风险评估中[3-6]．概率安全性评估、风险评估等诸多研究中应用马尔可夫模型时，主要用其求解与时间相关的状态概率，而求解的关键在于状态转移概率矩阵的计算．对于由 N 个 m 状态元件组成的系统，系统状态数为 mN．对于小系统而言，系统结构状态数目少，可直观确定状态转移概率矩阵的构成和非零元素的位置，非零元素的数值即状态转移率也容易获得．然而，大系统有巨大数量的结构状态，不能直观确定状态转移概率矩阵的构成和非零元素的位置，且若元件状态数m 较大时，很难快速准确地获得状态转移率．大系统状态转移概率矩阵的形成是十分复杂、繁重且耗时的工作，这大大限制了基于马尔可夫模型的概率安全性评估[5]等方法在实际工程中的应用．可见，研究大系统状态转移概率矩阵的快速形成方法具有重要的工程应用价值．

状态转移概率矩阵形成方法的相关研究很少．通常情况下，均是基于状态空间图获得状态转移概率矩阵[7]，其优点是状态及其相互转移有清晰的图形表示，但建立大系统的状态空间图是不可能的，故该方法对大系统的应用十分困难．Amoia 等[8-9]提出了面向计算机的基于Kronecker 代数形成大系统状态转移概率矩阵的方法，但其计算时间随着系统元件的增多而指数增长，无法满足实际工程应用的要求．

笔者提出一个快速形成马尔可夫模型状态转移概率矩阵的方法．定义元件状态转移率矩阵和系统状态数组，将系统状态转换为便于计算机存储与处理的数组，有效地描述了系统状态之间的转移；基于元件状态转移率矩阵和系统状态数组提出不受系统状态和元件状态数目限制快速计算状态转移率的方法，挖掘矩阵中非零元素的分布规律提出非零元素的快速定位方法，进而快速形成状态转移概率矩阵的稀疏存储；针对两状态元件组成的系统，提出基于给定的系统状态排序和服务状态集数组快速定位矩阵中非零元素的方法．对电力系统概率安全性评估的应用效果显著，证实了方法的实用性和有效性．

1 元件状态转移率矩阵和系统状态数组

1.1 元件状态转移率矩阵

定义元件状态转移率矩阵T，即若某元件在系统中有m 个状态，该元件的状态转移率矩阵为m × m维矩阵，第i 行、第j 列对应的元素为元件从元件状态i向元件状态j 的转移率．

为方便理解，均以两状态元件(即 m = 2)为例说明相关的概念和方法．元件有运行和停运两个状态，0 表示运行，1 表示停运．设元件事故率为λ，修复率为μ，则元件状态转移率矩阵T 为

1.2 系统状态数组

若第i 个元件有mi个状态，用0，1，…，mi-1 表示，各元件状态可重复组合[10]数为可知系统共有维状态数组，以此可快速生成所有系统状态．

以N 个两状态元件组成的系统为例，如图1 所示，将各元件状态可重复组合可得2N个系统状态数组，这些系统状态数组对应2N个系统状态．

图1 系统状态数组Fig.1 System configuration array

可见，通过系统状态数组可将系统状态转换为便于计算机存储与处理的数组并快速生成所有系统状态．

1.3 系统状态转移的描述

基于元件状态转移率矩阵和系统状态数组可有效描述系统状态之间的转移．以可修复两元件系统为例进行说明，元件有运行和停运两个状态，元件1的故障率和修复率分别为λ1和μ1；元件2 的故障率和修复率分别为λ2和μ2．

由第1.1 节得元件1 和元件2 的状态转移率矩阵T1和T2为

由第1.2 节可知，可用一个1× 2 维状态数组表示系统状态，数组中第1 个元素代表元件1 的状态，第2 个元素代表元件2 的状态，将各元件状态进行可重复组合生成4 个状态数组[0 0]、[1 0]、[0 1]和[1 1]，对应元件全部运行、元件1 停运元件2 运行、元件1运行元件2 停运和元件全部停运等4 个系统状态．

基于元件状态转移率矩阵和系统状态数组可快速获得状态之间的转换关系，有：[0 0]→[1 0]元件1从状态0 转移到状态1，即元件1 停运，由T1可知状态转移率为λ1；[1 0]→[0 0]元件1 从状态1 转移到状态0，即元件1 恢复运行，由T1可知状态转移率为μ1；[0 0]→[0 1]元件2 从状态0 转移到状态1，即元件2 停运，由T2可知状态转移率为λ2；[0 1]→[0 0]元件2 从状态1 转移到状态0，即元件2 恢复运行，由T2可知状态转移率为μ2；[1 0]→[1 1]元件2 从状态0 转移到状态1，即元件2 停运，由T2可知状态转移率为λ2；[0 1]→[1 1]元件1 从状态0 转移到状态1，即元件1 停运，由T1可知状态转移率为λ1；[1 1]→[0 1]元件1 从状态1 转移到状态0，即元件1 恢复运行，由T1可知状态转移率为μ1；[1 1]→[1 0]元件2 从状态1转移到状态0，即元件2 恢复运行，由T2可知状态转移率为μ2；[0 0]和[1 1]之间无转移，状态转移率为0．

与状态空间图相比，元件状态转移率矩阵和系统状态数组有效地描述了状态间的转换，便于计算机处理，更适用于大系统．

2 状态转移概率矩阵的快速计算方法

2.1 基本假设

取与文献[9，11]相同的假设如下：

(1) 系统的状态取决于各元件的状态，这与系统状态数组的定义刚好吻合；

(2) 鉴于实际工程应用中，常采用离散时间的马尔可夫模型分析电网等大规模系统，所以本文研究的对象是离散时间的马尔可夫模型；

(3) 当i≠j 且状态i 和j 中有2 个或2 个以上元件的状态相异时，状态i 到状态j 的状态转移率为0．

2.2 状态转移率的快速准确计算

系统状态i 和j 对应的状态数组为 Ai和 Aj，可快速计算状态转移率．

(1) Ai和Aj有2 个或2 个以上元素相异时，pij＝0；

(2) Ai和Aj有1 个元素相异时，计算状态转移率pij的公式为

式中：Tk表示元件k 的状态转移率矩阵；Ai(k)、Aj(k)表示Ai和Aj中的第k 个元素，即系统状态i 和j 中元件k 的状态；TkAi(k)Aj(k)表示元件k 从状态Ai(k)到Aj(k)的状态转移率；N 为系统中元件的个数．

(3) Ai和Aj有零个元素相异时，即状态i 向自身发生转移时，则计算状态转移率pij的公式为

该方法不受系统状态和元件状态数目的限制，无需状态空间图便可快速准确地计算状态转移率．

2.3 状态转移概率矩阵中非零元素的分布规律

系统状态转移概率矩阵中有大量的零元素，若能找到矩阵中非零元素的分布规律，将大大有助于实现非零元素的快速定位，进而实现矩阵的稀疏存储．

基于第2.2 节给出的系统状态转移率的计算方法以及诸多实例系统验证，笔者发现系统状态转移概率矩阵中的非零元素存在如下分布规律．

若系统中有N 个元件，元件有m 个状态，系统状态i 对应的状态数组为Ai，更改Ai某个元素形成的状态数组为Ai'，则

(1) Ai'对应的系统状态是i 可转移的系统状态；

(2) i 可向N×(m − 1)个不同于自身的系统状态发生转移；

(3) 系统状态转移概率矩阵中，每行均有N ×(m − 1) + 1个非零元素；

(4) 特别地，若系统中元件仅有运行和停运2 个状态，则任一系统状态可向N 个不同于自身的系统状态发生转移；系统状态转移概率矩阵中，每行有N＋1 个非零元素．

2.4 状态转移概率矩阵中非零元素的快速定位

基于第2.3 节所述状态转移概率矩阵中非零元素的分布规律得出非零元素的快速定位方法如下．

(1) 系统中有N 个元件，每个元件有m 个状态，系统状态集M 中系统状态的个数为NS，从1 至NS对M 中的状态进行编号；

(2) ∀i ∈ M，依次改变Ai的N 个元素，确定可向i 演变的 N ×(m − 1)个不同于i 的系统状态，这些状态构成的集合为Li；

(3) M 中可向i 演变的系统状态构成的集合为Mi＝{M∩Li，i}；

(4) 统计Mi中系统状态的个数为Ni，Ni即为系统状态转移概率矩阵第i 行中非零元素的个数；

(5) ∀j ∈ Mi，系统状态转移概率矩阵第i 行、第j 列即为非零元素的位置．

2.5 状态转移概率矩阵的计算流程

状态转移概率矩阵的计算流程如图2 所示．

图2 状态转移概率矩阵的计算流程Fig.2 Flow chart of state transition probability matrix

以9 节点电力系统(如图3 所示)为例说明上述状态转移概率矩阵的计算流程．

图3 9节点系统Fig.3 9-bus system

步骤1确定系统状态集和元件状态转移率矩阵．

考虑表1 所示的6 条线路，每条线路有运行和停运2 个状态，分别用0 和1 表示．将各元件状态可重复组合生成26＝64 个系统状态数组，与之对应的系统状态构成系统状态集．这些状态及其对应的状态数组中存储的元件状态如表2 所示．元件的事故率和修复率见表1，元件i 的状态转移率矩阵为

步骤2对状态集中的状态进行编号．对状态集中的64 个状态进行编号如表2 所示．

步骤3确定每个状态可转移的状态．以状态10 为例，状态10 对应的状态数组为[1 0 0 1 0 0]，依次改变状态数组中的元素获得6 个状态数组：[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]．对照表2 可得，以上6 个状态数组对应的状态为5、24、27、2、30、31．则状态10 可转移的状态有 7 个，它们是(按状态编号从小到大排序)2、5、10、24、27、30、31．

步骤4确定每行非零元素的位置．仍以状态10 为例．状态10 对应状态转移概率矩阵的第10行，则第10 行中2、5、10、24、27、30、31 列存在非零元素．

步骤5根据式(2)和式(3)计算非零元素的数值．以状态转移概率矩阵第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形成系统状态转移概率矩阵．

表1 线路及其编号Tab.1 Line and its order

循环执行步骤3～步骤5 形成示例系统状态转移概率矩阵的稀疏存储如表3～表6 所示．采用Matlab 编程，在联想Dual Core CPU、1.66,GHz、2,G内存计算机上计算示例系统的状态转移概率矩阵耗时0.000,33,s，与文献[9]耗时0.001,4,s 相比，本文方法计算速度更快．

需要说明的是，文献[9]是现有大系统马尔可夫模型状态转移概率矩阵形成方法中较好的方法，故在本节和第3 节中将其作为比较对象对本文方法进行分析．

2.6 两状态元件组成的系统的状态转移概率矩阵的快速计算

鉴于两状态元件组成的系统常被作为研究对象，笔者对其做了进一步研究，提出基于给定系统状态排序和服务状态集数组快速定位状态转移概率矩阵中非零元素的方法．仍以图4 所示系统为例说明如下．

1) 给定系统状态排序

系统有N 个元件，k 个元件故障对应的状态组成集合Bk．按下列规则排序系统状态集中的状态：

(1) 按B0至BN排序；

(2) 按如下规则排序Bk中状态：

① k＝0 时，B0＝{1}；

② k＝1 时，将A1中第1 个零元素至最后一个零元素依次改为1，形成状态2 至状态N＋1；

③ k≥2 时，Bk-1中编号最小的状态为s，状态数为N'，Ai中非零元素后有Ni个零元素．

将As中非零元素后第一个零元素至最后一个零元素依次改为1，形成状态s＋N'至状态s＋N'＋Ns-1，记为状态Ps至状态Qs；

将As＋i中非零元素后第一个零元素至最后一个零元素依次改为1，形成状态Ps＋i至状态Qs＋i，Ps＋i＝Qs＋i-1＋1，Qs＋i＝Qs＋i-1＋Ns＋i，i＝1，2，…，N'-1．

示例系统的系统状态排序见表2．

表2 状态及其编号以及服务状态集及其编号Tab.2 State and its order and service set and its order

表3 9节点系统的状态转移概率矩阵的第1行到第17行Tab.3 1st row to 17th row of state transition probability matrix of 9-bus system

表4 9节点系统的状态转移概率矩阵的第18行到第32行Tab.4 18th row to 32nd row of state transition probability matrix of 9-bus system

表5 9节点系统的状态转移概率矩阵的第33行到第48行Tab.5 33rd row to 48th row of state transition probability matrix of 9-bus system

2) 生成服务状态集数组

服务状态集(service set)用S 表示，第n 个服务状态集用Sn表示．除S1外，Sn中仅有一个状态xn，中最后一个元素为1．S 的定义为

示例系统的服务状态集见表2．

定义服务状态集数组C：C 是一个行数组，其列数为n，存储的元素为xn．

3) 基于C 快速定位状态转移概率矩阵的非零元素

设i 可转移的状态构成的状态集为Mi＝{Sup，i，Sdown}．若Ai中有Nup个1，Ndown个0，则Sup中有Nup个状态，由将Ai中某个1 改为0 后对应的状态构成；Sdown中有Ndown个状态，由将Ai中某个0 改为1 后对应的状态构成．可基于C 直接确定Mi如下所示．

(1) i＝1 时，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 时，Sdown＝S3；i＝N＋1 时，Sdown＝{C(1，3)，C(1，4)，…，C(1，i)}；i 为2 和N＋1 以外的状态时，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 表示，即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 种情况计算Sdown.

(4) ∀i∈B3∪B4∪…∪BN，仍可基于C 直接确定Mi，鉴于篇幅限制，本文不再赘述．

表6 9节点系统的状态转移概率矩阵的第49行到第64行Tab.6 49th row to 64th row of state transition probability matrix of 9-bus system

3 方法分析和应用实例

3.1 方法分析

用 Matlab 编程，基于联想 Dual Core CPU、1.66,GHz、2,G 内存计算机比较了本文方法和文献[9]方法，如表7 所示．结果表明了本文方法的优越性．

进一步地，通过Matlab 中cftool 曲线拟合工具，基于表7 数据按 ebxy a= (y 表示计算时间，x 表示系统规模，a 和b 表示拟合系数)进行指数曲线拟合，拟合结果如表8 所示．拟合误差eSSE取为拟合数据和原始数据对应点的误差的平方和．结果进一步证实了本文方法的有效性．对应的状态转移概率矩阵，便于应用的同时，可通过并行计算进一步加速状态转移概率矩阵的形成．

表7 本文方法与现有方法计算效率对比Tab.7 Comparison of calculating efficency of method presented with existing method

表8 本文方法与现有方法的拟合曲线参数Tab.8 Parameters of fitting curve of method presented with existing method

3.2 应用实例

此外，本文方法便于编程，且无需形成全状态空间对应的状态转移概率矩阵便可形成任意状态子集

第2 节所述方法已被成功应用于输电系统概率静态和动态安全性综合评估[12]．该评估模型计及了风电和负荷引起的节点注入的不确定性、系统静态和动态安全约束以及网络拓扑状态之间的演变，即事故的不确定性和元件的可修复性，通过求解一个线性向量微分方程获得 “到不安全时间” 的概率分布以评估系统的安全性，模型详见文献[12]．然而，由于模型求解方面存在无法快速形成微分方程系数矩阵的稀疏存储等问题，致使模型无法应用于复杂大系统．基于本文方法快速形成了微分方程系数矩阵的稀疏存储，推进了该模型的实用化．

以新英格兰10 机39 节点系统为例(如图4 所示)，考虑变压器线路之外的34 条线路，线路有运行和停运两个状态，用于安全性评估的状态集由B0、B1、B2、B3和B4构成，总计52,956 个状态．基于本文方法快速形成微分方程系数矩阵的稀疏存储有：微分方程系数矩阵每行有35 个非零元素，基于服务状态集数组直接确定非零元素的位置，非零元素的计算详见文献[12]．鉴于篇幅限制，未列出微分方程系数矩阵的计算结果．采用Matlab 编程，在联想Dual Core CPU、1.66,GHz、2,G 内存计算机上进行概率静态和动态安全性综合评估耗时1.1,s．

图4 新英格兰10机39节点系统Fig.4 New England 10 generators and 39-bus system

4 结 语

为了推进基于马尔可夫模型的概率安全性评估等方法在实际工程中的应用，对大系统马尔可夫模型状态转移概率矩阵的快速形成方法进行了研究．通过状态转移概率矩阵和系统状态数组对系统状态及其之间的转移进行了有效描述，将系统状态转换为便于计算机存储与处理的数组；提出快速准确计算状态转移率以及快速定位状态转移概率矩阵中非零元素的方法，从而快速形成了状态转移概率矩阵的稀疏存储．该方法计算准确、耗时短、便于编程，使自动快速形成大系统状态转移概率矩阵的稀疏存储成为可能，为运用马尔可夫模型分析电网等大规模系统提供了有利工具．

［1］郭永基. 可靠性工程原理[M]. 北京：清华大学出版社，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］吴文可，文福拴，薛禹胜，等. 基于马尔科夫链的电力系统连锁故障预测[J]. 电力系统自动化，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］李文沅. 电力系统风险评估-模型、方法和应用[M].北京：科学技术出版社，2006.Li Wenyuan.Risk Assessment of Power Systems-Models，Methods，and Applications[M]. Beijing ：Science and Technology Press，2006(in Chinese).

［7］张雪松，王 超，程晓东. 基于马尔可夫状态空间法的超高压电网继电保护系统可靠性分析模型[J]. 电网技术，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］焜刘嘉 ，王家生. 应用概率统计[M]. 北京：科学出版社，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.