基于重要抽樣函數參數解析優化的電網可靠性評估

莊汝學，耿 蓮，王 慧，姚浩威，黃 欣，趙凌駿，雙滟杰，趙 淵

(1.蘇州電力設計研究院有限公司，江蘇 蘇州 215000;2.重慶大學 電氣工程學院，重慶 400044)

蒙特卡羅仿真(MCS)因其靈活性和易于實現性在大規模電網可靠性評估中受到廣泛關注[1]。MCS的收斂性主要取決于系統可靠性水平，對于高可靠性系統，為達到給定仿真精度所需樣本容量通常較大，導致仿真時間很長，而方差削減技術則成為加快MCS仿真效率的有效手段。

重要抽樣(IS)法是一種廣為關注的高效方差削減技術，它通過改變隨機變量的概率密度函數(PDF)，使得對電網停電風險起重要作用的系統狀態更易被抽取，從而實現電網可靠性MCS的速度顯著提升[2]。優化求取重要抽樣概率密度函數(IS-PDF)的參數是IS取得良好加速性能的關鍵，如果求取的IS-PDF參數不合適則可能導致相反效果，降低仿真效率。近年來，交叉熵法(CEM)因其在MCS中的高效加速效果而備受關注。CEM首先被用于發電系統充裕性評估[3-5]，然后被進一步應用于旋轉備用充裕性評估[6-7]、組合系統充裕性評估[8-14]和短期運行風險評估[15-17]。

CEM的核心思路為：由于理論上的最優IS-PDF(即零方差PDF)無法獲取，退而求其次，可求取與理論最優IS-PDF之間的交叉熵(衡量兩個PDF相似性的指標)最小的PDF作為實際的IS-PDF。將交叉熵概念引入電網可靠性的MCS后，可有效克服原始MCS收斂慢的缺陷。在現有研究中，CEM法采用迭代算法實現對IS-PDF的參數尋優，每次迭代都隨機抽取固定數量(即預抽樣樣本數)的系統狀態進行系統性能評估，并將其中不良系統狀態用于IS-PDF的參數迭代計算。一方面多次反復迭代導致較大的計算成本；另一方面如果每次迭代所設定的預抽樣樣本數偏小，則在初始幾次迭代中難以保證抽取到數量足夠的不良系統狀態，導致參數迭代更新效率偏低，迭代次數增加，而設定的預抽樣樣本數偏大時，計算成本會顯著增加。

針對這一問題，提出一種全新的IS-PDF參數求取方法，對實際IS-PDF參數進行直接解析求取，避免傳統CEM法中的耗時迭代過程。具體思路包括：首先將每個系統故障狀態的理論最優IS-PDF表達成含元件最優無效度和元件原始無效度的非線性等式方程；由于系統故障狀態數量龐大，建立的等式方程數量太多，導致難以通過求解非線性方程組獲取元件最優無效度，為此進一步引入最小割集概念對系統故障狀態進行合并，大大削減方程數量；最后，由于滿足理論最優IS-PDF的元件最優無效度并不存在，建立的非線性方程組實際上并無準確解，因此對非線性方程組采用最小二乘估計實現元件最優無效度的有效估計。

1 最優IS-PDF的最小割集等式方程組

對于一個包含w個元件的電力系統，以xi表示元件i的狀態，xi=0表示正常，xi=1表示故障，則系統狀態可表示為x=(x1,x2,…,xw)。若各元件的隨機故障相互獨立，則系統狀態x的原始PDF為：

(1)

式中：u=(u1,u2,…,uw)表示原始概率密度函數f(x;u)的參數，ui為元件i的原始無效度。電網的故障概率pF(即失負荷概率)可用下式進行無偏估計：

(2)

式中：H(x)為系統性能測度函數，如果x是失負荷狀態，則H(x)=1，反之，H(x)=0；Ef(·)表示在概率密度函數f(x;u)下求數學期望；Ω為系統狀態空間；N為估計pF所需的樣本容量；xi為通過原始f(x;u)隨機抽取到的第i個系統狀態。

為加快電網可靠性評估的速度，IS法采用重要抽樣PDF，即使用g(x;v)來進行系統狀態的隨機抽取，則pF的估計公式如下：

(3)

式中：xi為采用g(x;v)隨機抽取的第i個系統狀態；W(x;u,v)=f(x;u)/g(x;v)為似然比函數，作用在于保證式(3)的期望值估計具有無偏性，即當N→∞時，pF的估計不會因抽樣函數的改變而出現偏離。g(x;v)通常采用與f(x;u)相同的函數形式：

(4)

當g(x;v)滿足式(5)時，式(3)的統計量具有零方差，此時的g(x;v)即為理論最優IS-PDF，v為理論上的最優IS-PDF參數，即元件最優無效度向量。

(5)

事實上理論最優IS-PDF只是理論存在，首先pF是未知的待估計量；其次即使pF已知，也不存在一個無效度向量v=(v1,v2,…,vw)恰好滿足式(5)。因此尋求一個無效度向量v=(v1,v2,…,vw)，使得g(x;v)和H(x)f(x;u)/pF充分接近，則成為可行的技術思路。

CEM以交叉熵(又名Kullback-Leible距離)D(gop,g)來衡量H(x)f(x;u)/pF和g(x;v)的相似性：

(6)

對式(6)最小化以求取參數v，簡化后如式(7)：

(7)

上式的等效形式如下：

(8)

將式(8)的積分公式離散化后可得：

(9)

式中N0為進行參數優化所需的樣本容量。要從f(x;u)中抽取足夠數量滿足H(xi)=1的xi效率極低，尤其所評估的系統可靠性很好時。為此，CEM采用迭代優化求解思路，即對式(9)再次應用重要抽樣思想，并在首次迭代時以原始f(x;u)作為重要抽樣函數，而在第k次迭代時將第k-1次迭代求得的g(x,vk-1)作為重要抽樣函數，如式(10)所示，其中xi從g(x,vk-1)抽取:

(10)

基于式(10)可推導出如下的vk計算公式，其中xi,j為系統狀態xi中第j個元件所處的狀態(1為故障，0為正常):

(11)

由上述可見，CEM每次迭代都需要隨機抽取N0個系統狀態，即xi(1≤i≤N0)，并對每個xi進行潮流或最優負荷削減計算以判斷H(xi)的取值，導致計算成本較高。為避免迭代求解帶來的較大計算成本，提出如下參數解析尋優方法。

設Ωf和Ωs分別表示故障和正常系統狀態子空間，則理論最優IS-PDF可表示為：

(12)

可知在理論最優IS-PDF下，正常系統狀態的概率變為零，而故障系統狀態的概率則按1/pF的比例線性增長。假設式中的pF可用一個預估值來近似，則任意一個故障系統狀態x∈ΩF都滿足下式：

(13)

對所有x∈ΩF均建立式(13)的非線性等式方程，則可以得到一個含待求參數v的方程組，該方程組所含方程個數等于故障系統狀態的數目。由于屬于ΩF的故障系統狀態數量龐大，一方面不可能也沒必要識別出所有屬于ΩF的x，另一方面即使能完全識別，龐大的非線性方程數目也將導致無法求解。為此引入最小割集概念解決上述難題。割集是一組元件的集合，當它們失效時會導致系統故障，而引起系統故障的失效元件集中的最小子集稱為最小割集。

事實上每個最小割集可代表一系列系統故障狀態。以圖1的5元件系統為例，Θ={1,4,5}為三階最小割集，即當x1=x4=x5=1時，無論元件2、3處于何種狀態組合(注：元件2和3的狀態組合共4種)，系統都會處于故障狀態，可見Θ涵蓋了4個引起系統失效的元件狀態組合。

圖1 用于闡釋最小割集概念的5元件系統Fig. 1 A system with 5 elements for illustrating the concept of the minimum cutset

推而廣之，如果某電網總元件數為w，且元件都為2狀態Markov模型，則任意一個k階最小割集Θ可表征2w-k個引起系統失效的元件狀態組合，而對于其中第d個元件狀態組合，所有屬于Θ的元件其狀態必然為1，即有xd,i=1，?i∈Θ，故參考式(13)可列出如下等式方程：

(14)

上述方程可列寫出2w-k個，所有方程相加可得：

(15)

考慮到下式成立：

(16)

式(15)可簡化如下:

(17)

系統中能夠辨識出的n個最小割集Θk(1≤k≤n)都滿足式(17)，由此得到方程組(18)，即最優IS-PDF的最小割集等式方程組，這是求取元件最優無效度v的重要理論公式。采用此方法，可將方程數目由最初等于故障系統狀態數降低為最小割集數，從而有效減少了方程數目，降低了方程組求解難度。

(18)

2 最小割集辨識及元件v的最小二乘估計

為建立方程組(18)，需要事先知道n個最小割集Θk(1≤k≤n)，因此首要問題是如何以較小的計算成本搜索和辨識出Θk(1≤k≤n)。對于復雜大電網，采用傳統的最小通路法[18]進行最小割集的辨識效率較低?？紤]到高階最小割集事件的識別所需計算成本較高，同時高階最小割集事件的出現概率遠小于低階最小割集，在元件最優無效度的計算中其影響也遠小于低階最小割集，因此基于最小割集的基本概念，采用解析枚舉的方式對最小割集事件進行識別。首先根據計算成本要求事先給定故障元件的最高枚舉階數(即同時故障的最大元件數)R，然后執行以下最小割集辨識算法：

①初始化：設置最小割集數n1=0，r=1。

②對只有第j(1≤j≤w)個元件故障的1階事件進行潮流或最優潮流計算，若該事件引起電網失負荷，則該事件是1階最小割集事件，令n1=n1+1,Θn1={j}。所有1階事件分析完畢后，令r=r+1。

③如果r>R，則最小割集搜索完畢，令n=n1，算法結束，反之轉步驟④。

④對每個r階元件故障事件分別進行潮流或最優潮流計算，由故障元件的編號構成集合Ψ。若該r階事件引起電網失負荷，則判斷是否存在某個Θk(1≤k≤n1)使得Θk?Ψ，如果存在，則Ψ不是最小割集事件，反之令n1=n1+1,Θn1=Ψ。當所有r階事件分析完畢后，令r=r+1，轉步驟③。

在求取系統的最小割集之后，還需預估系統失負荷概率pF的數值，否則式(18)無法求解。由可靠性評估理論可知，失負荷概率pF可由最小割集事件表示為：

pF=P(Θ1∪Θ2…∪Θn)=P(Θ1)+P(Θ2)+…+P(Θn)-(P(Θ1∩Θ2)+…+

P(Θn-1∩Θn))…(-1)n-1(P(Θ1∩Θ2…∩Θn))。

(19)

式中：P(Θk)為最小割集Θk出現的概率，符號∩表示最小割集事件同時(交疊)出現。上式為pF的理論計算公式，在實際計算時，由于最小割集事件交疊的組合數太多，采用該式計算pF雖然理論可行，但非常煩瑣和費時。為克服這一缺陷，應用如下近似計算方法，即分別按(20)和(21)估計失負荷概率pF的上界和下界。

(20)

(21)

如果一個系統中元件的可靠性都較好，則多個Θk交疊出現的可能性很小，即使忽略交疊也不會導致較大誤差，此時可直接采用式(20)的上界估計。如果考慮多個(例如2個)Θk同時發生的可能，可采用式(20)和式(21)的平均值來近似估計。

如第1節所述，式(5)的gop(x;v)只是理論存在，并不存在一個參數v=(v1,v2,…,vw)能恰好滿足式(5)，故由式(5)衍生出的方程組(18)無精確解。為此，采用最小二乘估計，通過式(18)的非線性方程組對v進行近似估計。首先將式(18)按下式轉換為線性方程組：

(22)

將方程組(22)的求解轉化為最小二乘估計問題：

(23)

式中S為誤差平方和。

定義H為“最小割集-元件關聯矩陣”，且矩陣階數為n×w。H中的行對應最小割集，而列對應元件，例如第k行表示最小割集Θk，若某元件i屬于Θk，即i∈Θk，則第k行第i列的元素H(k,i)=1，除此外第k行其余元素都為零。

式(22)可用矩陣形式表示為：

Hy=b;

(24)

y=lnv=[lnv1,lnv2,…,lnvw]T;

(25)

b=H·lnu-E·ln(pF)。

(26)

式中E為w階單位對角陣。

最小二乘估計問題的矩陣表達形式為：

(27)

在極值點處梯度為零，即?S/?y=0，則有：

HTHy-HTb=0。

(28)

對上式的線性方程組進行求解，可解得y，并進而得到v：

y=(HTH)-1HTb；

(29)

v=ey。

(30)

利用本節方法可估計屬于最小割集中的元件的最優無效度，但某些對系統可靠性指標貢獻很小的元件，可能未曾出現在任何一個搜索到的最小割集中，此類元件其最優無效度可取為原始無效度。

3 基于參數解析尋優的重要抽樣算法流程

基于參數解析優化，即基于最小割集辨識和最小二乘估計實現v的優化求取，然后采用重要抽樣進行電網可靠性的仿真計算，整個過程分為2個環節：1)參數優化環節：搜尋n個Θk(1≤k≤n)并建立式(24)～(26)，由式(30)得到IS-PDF的參數v；2)主抽樣環節：以g(x;v)為抽樣函數，隨機抽取N個系統狀態，并判斷每個系統狀態是否失負荷以及負荷削減量，更新可靠性指標統計量。流程圖見圖2，詳細步驟如下。

人工線下單據傳遞模式主要存在兩方面問題：一是，流通過程中可能造成單據受損、遺失的問題；二是，單據遞交人工成本較高，對于部分地理位置較遠的需求單位，每月定期遞交單據將增加人工、差旅和時間成本。

圖2 基于參數解析尋優的電網可靠性重要抽樣算法Fig. 2 Algorithm description for the proposed IS method

①初始化：輸入電網元件的電氣參數和原始無效度u以及網絡拓撲等數據，設置允許階數R，并令N=i=1。

②依據上一節的算法搜尋n個Θk(1≤k≤n)。

③依據Θk(1≤k≤n)建立式(24)～(26)。

④按式(29)(30)計算元件最優無效度v。

⑤進入主抽樣環節，設置允許的方差系數β*(通常0.01～0.05之間)。

⑥采用g(x;v)抽取第i個系統狀態xi。

⑦對xi進行潮流或最優潮流計算，從而得到失負荷標志H(xi)、負荷削減量C(xi)，以及似然比W(xi;u,v)=f(xi;u)/g(xi;v)的數值。

(31)

(32)

4 算例分析

4.1 MRBTS可靠性測試系統的評估分析

RBTS[19]可靠性測試系統由于所含元件較少，有利于對各種不同方法的應用效果進行直觀比較，但該系統不符合N-1確定性可靠性準則，為此在該系統的節點5和6之間增加一條輸電線路L10，并將改變后的系統稱為MRBTS可靠性測試系統，見圖3。負荷采用年峰荷模型，基于以下4種方法進行電網可靠性評估，其中方法1～3的收斂條件為EENS指標的方差系數β*=0.01。

圖3 MRBTS可靠性測試系統的電氣接線圖Fig. 3 Single line diagram of MRBTS

方法1(ISA)：采用IS法，但基于解析參數優化，且參數優化環節的枚舉階數R=3；

方法2(ISCE)：采用IS法，但基于交叉熵參數優化，在參數優化階段，每次迭代所需的樣本容量為20 000；

方法3(CMCS)：原始的非序貫蒙特卡洛仿真，即始終采用原始f(x;u)進行隨機抽樣；

方法4(ENU)：采用狀態枚舉法，且系統狀態枚舉到6階。

將MRBTS系統中部分元件的原始無效度u以及由ISA和ISCE方法在參數優化階段得到的最優無效度v列于表1中。從表1可知，通過ISA和ISCE優化求取的v與原始的u相比都有明顯變化，但二者求取的v總體變化規律大致一樣。ISCE法所得的v存在一個矛盾，即安裝地點、安裝容量以及原始參數u都完全一樣的2個元件，在參數優化后其v卻存在明顯差異，例如1#和2#發電機、6#和9#發電機、1#和6#輸電線路，但采用本研究中提出的ISA法可以有效避免上述問題。其根本原因在于：ISCE法在參數優化階段需要進行多次迭代，且每次迭代都需要對元件狀態進行大量隨機抽樣，由于隨機抽樣的不確定性，即使2個元件完全一模一樣，它們在每次隨機抽樣中也不一定被同時抽取到處于同樣狀態，因此根據式(11)的參數更新公式計算得到的結果就會存在差異。而ISA法并非采用隨機抽樣，而是基于解析思路，通過枚舉系統狀態并尋找其中最小割集，由此建立非線性方程組。由于采用枚舉方式，2個相同的元件能保證以同等的概率被枚舉，因此相同的元件在參數解析優化后其最優無效度也完全相同。

表1 交叉熵參數和解析參數優化的結果比較

對MRBTS進行可靠性評估的結果列于表2，并以ENU法的可靠性指標計算結果作為參照基準進行比較。在可靠性指標的計算準確性上，ISA法得到的LOLP和EENS指標與ENU法相比差異很小，這說明本文方法能保證較好的計算精度。而在所需樣本數和計算耗時上，ISA法由于在可靠性評估階段采用了重要抽樣，所需樣本比CMCS法大大減少，計算速度也得到顯著提升。ISCE法同樣在可靠性評估階段采用了重要抽樣，但由于在參數優化階段采取了迭代尋優方式，并且為了能夠在初始迭代中能抽取到足夠數量的不良系統狀態，每次迭代所需的抽樣樣本數較大，因此ISCE法在可靠性評估階段雖然所需樣本數比ISA法稍小，但總的樣本數和計算耗時卻明顯高得多。結果表明本文ISA法可在保證計算準確性的情況下大幅度提高評估速度，提高了可靠性評估效率。

表2 MRBTS在不同方法下的準確性和速度比較

基于枚舉得到的最小割集，用公式(19)解析計算的LOLP指標值為0.007 9，相比于上表結果明顯偏小，主要在于枚舉階數R=3不能太高，否則參數優化階段的計算成本將顯著增加。另一方面，從上表也可看出，最小割集的階數雖然不高，但并不影響本文方法的計算準確性。其原因如下：參數優化階段的目的在于獲取重要抽樣概率密度函數的參數，當參數計算完畢，才進入真正的大電網可靠性指標計算階段，即主抽樣階段。在主抽樣階段抽取的樣本數有限的情況下，計算準確性常采用可靠性指標統計量的方差系數來衡量。參數優化階段采用不同方式得到重要抽樣函數參數，將使參數優化階段的計算成本出現差異，進而影響整個電網可靠性評估(包含參數優化和主抽樣階段)的計算耗時，但并不影響主抽樣階段可靠性指標的計算準確性，因為計算準確性由主抽樣階段所給定的方差系數決定。最小割集階數越高，得到的重要抽樣函數將使主抽樣階段收斂更快，但也帶來兩個方面的影響：一是參數優化階段需要枚舉和分析的系統狀態數呈指數函數急劇增加，計算成本急劇增長，整個電網可靠性評估效率反而下降；二是最小割集階數過高時，重要抽樣函數的參數變化將進入飽和階段，即參數計算值隨最小割集階數的增加只是微小變動。由此可見，最小割集的階數并非越高越好，同時考慮到對系統可靠性有重要影響的元件在階數不高的割集事件中就能體現出來，因此對階數不高的割集事件中的元件無效度進行優化，突顯它們出現的概率，就能在后續主抽樣階段實現可靠性指標計算的顯著加速。

4.2 IEEE-RTS79可靠性測試系統

將本文方法進一步應用于IEEE-RTS79系統[20]，考慮到年峰荷下該系統的可靠性水平較差，采用原始MCS也能較快滿足收斂條件，因此在90%負荷水平下采用基于同樣的4種方法進行評估分析。

方法1(ISA)：采用IS法，但基于解析參數優化，且參數優化環節的枚舉階數R=2；

方法2(ISCE)：采用IS法，但基于交叉熵參數優化，在參數優化階段，每次迭代所需的樣本容量為25 000；

方法3(CMCS)：原始的非序貫蒙特卡洛仿真，即始終采用原始f(x;u)進行隨機抽樣；

方法4(ENU)：采用狀態枚舉法，且系統狀態枚舉到5階。

對IEEE-RTS79進行可靠性評估后的結果如表3所示，可見本文ISA法的可靠性指標計算準確性和其他方法大致相當，但可靠性評估效率卻明顯高于其他方法。雖然ISCE法也采用了重要抽樣，但迭代方式的參數尋優效率顯著低于解析方式參數尋優，例如ISCE法需要3次迭代共隨機抽取75 000個樣本用于參數更新，但ISA法只需枚舉2485個系統狀態，可見參數優化階段ISCE法所需樣本遠大于ISA法，導致ISCE法的總體效率不如ISA法。

表3 IEEE-RTS79在不同方法下的準確性和速度比較

同樣基于枚舉得到的最小割集直接計算LOLP指標，其數值為0.011 5，與上表結果相比也明顯偏小，原因在于枚舉階數R=2不能太高，以避免預抽樣階段出現較大計算成本。同時可見，本文方法得到的可靠性指標計算準確性并未受割集階數不高的影響。

5 結 語

基于非序貫蒙特卡羅模擬，從解析法的角度將故障系統狀態的最優重要抽樣概率密度表達為含元件最優無效度與原始無效度的非線性方程組；利用最小割集能代表一類故障系統狀態的特性，基于最小割集實現了方程合并，降低了方程組求解的難度；再引入最小二乘估計，以誤差最小為優化目標實現了最優無效度的解析估計。本文方法有效避免了CEM法在迭代參數尋優中隨機抽樣帶來的較大計算負擔，提高了蒙特卡羅仿真的效率。