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

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

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

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

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

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

1 最優(yōu)IS-PDF的最小割集等式方程組

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

(1)

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

(2)

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

為加快電網(wǎng)可靠性評(píng)估的速度，IS法采用重要抽樣PDF，即使用g(x;v)來進(jìn)行系統(tǒng)狀態(tài)的隨機(jī)抽取，則pF的估計(jì)公式如下：

(3)

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

(4)

當(dāng)g(x;v)滿足式(5)時(shí)，式(3)的統(tǒng)計(jì)量具有零方差，此時(shí)的g(x;v)即為理論最優(yōu)IS-PDF，v為理論上的最優(yōu)IS-PDF參數(shù)，即元件最優(yōu)無效度向量。

(5)

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

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

(6)

對(duì)式(6)最小化以求取參數(shù)v，簡(jiǎn)化后如式(7)：

(7)

上式的等效形式如下：

(8)

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

(9)

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

(10)

基于式(10)可推導(dǎo)出如下的vk計(jì)算公式，其中xi,j為系統(tǒng)狀態(tài)xi中第j個(gè)元件所處的狀態(tài)(1為故障，0為正常):

(11)

由上述可見，CEM每次迭代都需要隨機(jī)抽取N0個(gè)系統(tǒng)狀態(tài)，即xi(1≤i≤N0)，并對(duì)每個(gè)xi進(jìn)行潮流或最優(yōu)負(fù)荷削減計(jì)算以判斷H(xi)的取值，導(dǎo)致計(jì)算成本較高。為避免迭代求解帶來的較大計(jì)算成本，提出如下參數(shù)解析尋優(yōu)方法。

設(shè)Ωf和Ωs分別表示故障和正常系統(tǒng)狀態(tài)子空間，則理論最優(yōu)IS-PDF可表示為：

(12)

可知在理論最優(yōu)IS-PDF下，正常系統(tǒng)狀態(tài)的概率變?yōu)榱悖收舷到y(tǒng)狀態(tài)的概率則按1/pF的比例線性增長(zhǎng)。假設(shè)式中的pF可用一個(gè)預(yù)估值來近似，則任意一個(gè)故障系統(tǒng)狀態(tài)x∈ΩF都滿足下式：

(13)

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

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

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

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

(14)

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

(15)

考慮到下式成立：

(16)

式(15)可簡(jiǎn)化如下:

(17)

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

(18)

2 最小割集辨識(shí)及元件v的最小二乘估計(jì)

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

①初始化：設(shè)置最小割集數(shù)n1=0，r=1。

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

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

④對(duì)每個(gè)r階元件故障事件分別進(jìn)行潮流或最優(yōu)潮流計(jì)算，由故障元件的編號(hào)構(gòu)成集合Ψ。若該r階事件引起電網(wǎng)失負(fù)荷，則判斷是否存在某個(gè)Θk(1≤k≤n1)使得Θk?Ψ，如果存在，則Ψ不是最小割集事件，反之令n1=n1+1,Θn1=Ψ。當(dāng)所有r階事件分析完畢后，令r=r+1，轉(zhuǎn)步驟③。

在求取系統(tǒng)的最小割集之后，還需預(yù)估系統(tǒng)失負(fù)荷概率pF的數(shù)值，否則式(18)無法求解。由可靠性評(píng)估理論可知，失負(fù)荷概率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出現(xiàn)的概率，符號(hào)∩表示最小割集事件同時(shí)(交疊)出現(xiàn)。上式為pF的理論計(jì)算公式，在實(shí)際計(jì)算時(shí)，由于最小割集事件交疊的組合數(shù)太多，采用該式計(jì)算pF雖然理論可行，但非常煩瑣和費(fèi)時(shí)。為克服這一缺陷，應(yīng)用如下近似計(jì)算方法，即分別按(20)和(21)估計(jì)失負(fù)荷概率pF的上界和下界。

(20)

(21)

如果一個(gè)系統(tǒng)中元件的可靠性都較好，則多個(gè)Θk交疊出現(xiàn)的可能性很小，即使忽略交疊也不會(huì)導(dǎo)致較大誤差，此時(shí)可直接采用式(20)的上界估計(jì)。如果考慮多個(gè)(例如2個(gè))Θk同時(shí)發(fā)生的可能，可采用式(20)和式(21)的平均值來近似估計(jì)。

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

(22)

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

(23)

式中S為誤差平方和。

定義H為“最小割集-元件關(guān)聯(lián)矩陣”，且矩陣階數(shù)為n×w。H中的行對(duì)應(yīng)最小割集，而列對(duì)應(yīng)元件，例如第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階單位對(duì)角陣。

最小二乘估計(jì)問題的矩陣表達(dá)形式為：

(27)

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

HTHy-HTb=0。

(28)

對(duì)上式的線性方程組進(jìn)行求解，可解得y，并進(jìn)而得到v：

y=(HTH)-1HTb；

(29)

v=ey。

(30)

利用本節(jié)方法可估計(jì)屬于最小割集中的元件的最優(yōu)無效度，但某些對(duì)系統(tǒng)可靠性指標(biāo)貢獻(xiàn)很小的元件，可能未曾出現(xiàn)在任何一個(gè)搜索到的最小割集中，此類元件其最優(yōu)無效度可取為原始無效度。

3 基于參數(shù)解析尋優(yōu)的重要抽樣算法流程

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

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

圖2 基于參數(shù)解析尋優(yōu)的電網(wǎng)可靠性重要抽樣算法Fig. 2 Algorithm description for the proposed IS method

①初始化：輸入電網(wǎng)元件的電氣參數(shù)和原始無效度u以及網(wǎng)絡(luò)拓?fù)涞葦?shù)據(jù)，設(shè)置允許階數(shù)R，并令N=i=1。

②依據(jù)上一節(jié)的算法搜尋n個(gè)Θk(1≤k≤n)。

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

④按式(29)(30)計(jì)算元件最優(yōu)無效度v。

⑤進(jìn)入主抽樣環(huán)節(jié)，設(shè)置允許的方差系數(shù)β*(通常0.01～0.05之間)。

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

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

(31)

(32)

4 算例分析

4.1 MRBTS可靠性測(cè)試系統(tǒng)的評(píng)估分析

RBTS[19]可靠性測(cè)試系統(tǒng)由于所含元件較少，有利于對(duì)各種不同方法的應(yīng)用效果進(jìn)行直觀比較，但該系統(tǒng)不符合N-1確定性可靠性準(zhǔn)則，為此在該系統(tǒng)的節(jié)點(diǎn)5和6之間增加一條輸電線路L10，并將改變后的系統(tǒng)稱為MRBTS可靠性測(cè)試系統(tǒng)，見圖3。負(fù)荷采用年峰荷模型，基于以下4種方法進(jìn)行電網(wǎng)可靠性評(píng)估，其中方法1～3的收斂條件為EENS指標(biāo)的方差系數(shù)β*=0.01。

圖3 MRBTS可靠性測(cè)試系統(tǒng)的電氣接線圖Fig. 3 Single line diagram of MRBTS

方法1(ISA)：采用IS法，但基于解析參數(shù)優(yōu)化，且參數(shù)優(yōu)化環(huán)節(jié)的枚舉階數(shù)R=3；

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

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

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

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

表1 交叉熵參數(shù)和解析參數(shù)優(yōu)化的結(jié)果比較

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

表2 MRBTS在不同方法下的準(zhǔn)確性和速度比較

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

4.2 IEEE-RTS79可靠性測(cè)試系統(tǒng)

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

方法1(ISA)：采用IS法，但基于解析參數(shù)優(yōu)化，且參數(shù)優(yōu)化環(huán)節(jié)的枚舉階數(shù)R=2；

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

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

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

對(duì)IEEE-RTS79進(jìn)行可靠性評(píng)估后的結(jié)果如表3所示，可見本文ISA法的可靠性指標(biāo)計(jì)算準(zhǔn)確性和其他方法大致相當(dāng)，但可靠性評(píng)估效率卻明顯高于其他方法。雖然ISCE法也采用了重要抽樣，但迭代方式的參數(shù)尋優(yōu)效率顯著低于解析方式參數(shù)尋優(yōu)，例如ISCE法需要3次迭代共隨機(jī)抽取75 000個(gè)樣本用于參數(shù)更新，但I(xiàn)SA法只需枚舉2485個(gè)系統(tǒng)狀態(tài)，可見參數(shù)優(yōu)化階段ISCE法所需樣本遠(yuǎn)大于ISA法，導(dǎo)致ISCE法的總體效率不如ISA法。

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

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

5 結(jié) 語(yǔ)

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