蔡霽霖，郝麗麗，張柯琪

（1. 南京工業(yè)大學(xué)電氣工程與控制科學(xué)學(xué)院，江蘇省 南京市 211816；2. 國網(wǎng)常州供電公司，江蘇省 常州市 210024）

0 引言

隨著新能源發(fā)電滲透率增加，傳統(tǒng)發(fā)電廠年平均利用小時(shí)數(shù)逐年下降，但棄風(fēng)、棄光、棄水現(xiàn)象仍時(shí)有發(fā)生，引發(fā)了電力系統(tǒng)中發(fā)電資源是否冗余的討論［1-2］。因此，作為發(fā)電充裕度量化手段的電力系統(tǒng)可靠性評(píng)估得到關(guān)注。

早期的電力系統(tǒng)可靠性評(píng)估主要依賴解析法，其計(jì)算復(fù)雜度取決于電力系統(tǒng)解空間的復(fù)雜度［3］。為了緩解狀態(tài)空間過大給解析法帶來的效率下降問題，文獻(xiàn)［4］提出一種基于影響增量和高階故障狀態(tài)削減的可靠性評(píng)估解耦方法。文獻(xiàn)［5］則提出一種基于卷積法和傳輸點(diǎn)輸電能力-負(fù)荷解耦模型的可靠性評(píng)估方法。然而，隨著大量可再生能源并網(wǎng)，系統(tǒng)內(nèi)隨機(jī)變量大幅增加，其相關(guān)性結(jié)構(gòu)也愈發(fā)復(fù)雜［6-7］，傳統(tǒng)解析法的瓶頸無法徹底解決。因此，基于蒙特卡洛采樣（Monte Carlo sampling，MCS）思想的可靠性評(píng)估模擬法受到了關(guān)注［8-11］。

雖然理論上只要樣本容量足夠大，MCS 法得到的可靠性指標(biāo)將無限逼近真實(shí)值，但對(duì)大規(guī)模高可靠性電力系統(tǒng)進(jìn)行可靠性評(píng)估時(shí)，對(duì)每一個(gè)樣本都運(yùn)行一次計(jì)及網(wǎng)絡(luò)安全約束和隨機(jī)變量相關(guān)性的最優(yōu)潮流算法的代價(jià)較高。因此，需盡可能縮減MCS法的樣本容量。為此，許多研究將MCS 法與各類方差削減方法結(jié)合起來，其中效果最為突出的是重要性采樣（importance sampling，IS）法［12-14］。文獻(xiàn)［15］利用IS 法修正前后停運(yùn)率差異大的元件，構(gòu)成重要狀態(tài)空間來提升采樣效率。文獻(xiàn)［16］將元件狀態(tài)空間分為低階和高階子空間，并采用IS 法計(jì)算高階空間可靠性指標(biāo)。文獻(xiàn)［17］通過在IS 法修正后的系統(tǒng)狀態(tài)函數(shù)中增加額外控制項(xiàng)，進(jìn)一步減小IS 法估計(jì)可靠性指標(biāo)時(shí)的方差。為同時(shí)處理連續(xù)和離散變量，文獻(xiàn)［18］提出一種基于交叉熵思想的IS 法，求取各類隨機(jī)變量的重要性分布來替代各自的原始分布，提高故障樣本的出現(xiàn)概率。文獻(xiàn)［19］在文獻(xiàn)［18］的基礎(chǔ)上，針對(duì)系統(tǒng)中同一節(jié)點(diǎn)常接有多臺(tái)參數(shù)型號(hào)相同機(jī)組的情況，利用多項(xiàng)分布實(shí)現(xiàn)隨機(jī)變量降維。

然而，IS 法對(duì)于偏移后的分布函數(shù)仍采用簡單隨機(jī)采樣獲取樣本，因此依舊存在樣本重復(fù)度高的可能，故無法保證用少量隨機(jī)樣本表征概率分布全貌的效果。實(shí)際上，提高樣本對(duì)概率分布的表征能力同樣能降低估計(jì)量方差，而分層抽樣法就是按照這一思路構(gòu)建的一類方法［20］，如拉丁超立方采樣（Latin hypercube sampling，LHS）法［21-23］。

從理論上來說，IS 法改變概率分布，不改變樣本抽取方法，LHS 法恰好與之相反，故將IS 法與LHS 法結(jié)合，應(yīng)該能提升指標(biāo)計(jì)算效率。文獻(xiàn)［24］對(duì)此進(jìn)行了嘗試，但只考慮了離散變量，未涉及LHS 法在連續(xù)變量相關(guān)性控制時(shí)如何與IS 法配合。事實(shí)上，LHS 法與IS 法結(jié)合存在以下問題：LHS 法在對(duì)樣本進(jìn)行相關(guān)性控制時(shí)，需使用相關(guān)系數(shù)矩陣；而IS 法則需要預(yù)先指定偏移分布的帶參數(shù)函數(shù)形式［25］，但只有少量函數(shù)形式能夠方便地提取出相關(guān)系數(shù)矩陣，可這些函數(shù)形式不一定能完全表征各類隨機(jī)變量的復(fù)雜邊際分布和相關(guān)性結(jié)構(gòu)特征。

針對(duì)上述問題，本文提出了一種非參數(shù)重要性分層抽樣法，該方法主要由IS 法中偏移分布的非參數(shù)構(gòu)建方法和LHS 法相關(guān)性控制部分與IS 法偏移分布兼容的實(shí)現(xiàn)方法2 個(gè)部分組成。將該方法應(yīng)用于含多新能源廠站的電力系統(tǒng)可靠性評(píng)估，并利用實(shí)際新能源廠站的歷史出力數(shù)據(jù)和IEEE RTS-79測試系統(tǒng)設(shè)計(jì)了算例，驗(yàn)證了該方法的有效性。

1 小概率事件分析中的IS 和LHS

1.1 小概率事件分析的數(shù)學(xué)模型

在小概率事件分析中，最受關(guān)注的是對(duì)事件發(fā)生的概率進(jìn)行估計(jì)。本文關(guān)注的電力系統(tǒng)可靠性評(píng)估中，離散變量（發(fā)電機(jī)組、線路運(yùn)行狀態(tài)等）與連續(xù)變量（負(fù)荷、風(fēng)電出力等）相互獨(dú)立，在此前提下事件發(fā)生概率的數(shù)學(xué)表達(dá)式如下。

式中：S(x)為系統(tǒng)狀態(tài)評(píng)價(jià)函數(shù)；Sth為人為設(shè)定的閾值；γ 表示S（x）小于Sth的概率；E(·)表示求期望值；I(·)表示一個(gè)指示函數(shù)，括號(hào)內(nèi)表達(dá)式成立時(shí)取1，反之取0；x 為系統(tǒng)內(nèi)所有隨機(jī)變量組成的向量，簡稱系統(tǒng)狀態(tài)向量；xdv為所有離散變量組成的向量；xcv為所有連續(xù)變量組成的向量，三者間的關(guān)系為x=[ xdv，xcv]；Ωdv和Ωcv分別為xdv和xcv的取值空間；p(xdv)和f (xcv)分別為離散變量和連續(xù)變量的概率密度函數(shù)（probability density function，PDF）。

要估計(jì)式（1）中γ 的值，一個(gè)最直接的方法就是采用MCS 法，由此構(gòu)建出的估計(jì)量如下。

式中：βMCS為γMCS的方差系數(shù)；var(·)表示求方差。

1.2 基于交叉熵思想的IS 法

IS 法的基本思想是對(duì)式（1）做如下等價(jià)變換：

由式（4）可知，IS 法本質(zhì)上是一個(gè)兩階段過程：第1 階段對(duì)應(yīng)式（4）前2 行，是對(duì)隨機(jī)向量的概率分布進(jìn)行修正并得到IS-PDF；第2 階段對(duì)應(yīng)式（4）第3 行，用MCS 法修正分布采樣，并估計(jì)關(guān)心的指標(biāo)數(shù)值。

文獻(xiàn)［18］指出，IS 法理論上最優(yōu)的IS-PDF 是無法得到的，故許多學(xué)者提出了逼近最優(yōu)分布的方法［26-27］。其中，基于交叉熵思想的修正方法是最為廣泛使用的方法之一［28-29］，其基本思想是預(yù)先為h(xdv)和g(xcv)選取含待定參數(shù)的函數(shù)形式，將待定參數(shù)作為優(yōu)化變量，以該函數(shù)形式與理論最優(yōu)分布的交叉熵最小為目標(biāo)函數(shù)，構(gòu)建并求解優(yōu)化模型，從而得到待定參數(shù)最優(yōu)值，并據(jù)此確定h(xdv)和g(xcv)的表達(dá)式［18］。

1.3 LHS 法

LHS 法也分為2 個(gè)階段：獨(dú)立分層采樣和相關(guān)性控制［23］。第1 階段進(jìn)行獨(dú)立分層采樣時(shí)，暫不考慮隨機(jī)向量x 中各變量間的相關(guān)性，只用各變量的邊際分布進(jìn)行采樣。記x 中的變量個(gè)數(shù)為K，所需樣本容量為N。則第1 階段結(jié)束時(shí)將會(huì)得到一個(gè)各行之間沒有相關(guān)性的樣本矩陣X。在第2 階段的相關(guān)性控制過程中，利用連續(xù)變量間相關(guān)系數(shù)矩陣，在X 中加入指定的相關(guān)性，具體步驟可參見文獻(xiàn)［23］。

2 非參數(shù)重要性分層抽樣法的算法原理設(shè)計(jì)

傳統(tǒng)IS 法的重點(diǎn)在于找到一個(gè)效果較好的IS-PDF，將樣本集中到對(duì)指標(biāo)貢獻(xiàn)更大的“重要性區(qū)域”，而在該區(qū)域內(nèi)的采樣仍然采用簡單隨機(jī)采樣，存在出現(xiàn)大量重復(fù)或相似樣本的可能性，故其加速計(jì)算的潛力仍有挖掘空間。因此，本文將IS 法與LHS 法結(jié)合，提出一種非參數(shù)重要性分層抽樣（nonparametric stratified importance sampling，NSIS）法來進(jìn)一步挖掘IS 法的潛力。

2.1 NSIS 法設(shè)計(jì)思路

1）采用混合高斯模型（GMM）這一非參數(shù)化模型作為IS-PDF，GMM 的表達(dá)式如下［25］。

2）設(shè)計(jì)一種與GMM 兼容的改進(jìn)LHS 方法，選擇GMM 的理由如下。

①從保證IS 法性能不受損失的角度：理論最優(yōu)IS-PDF 一般是一個(gè)非常規(guī)分布，例如存在多個(gè)不相連的、偏度不為0 的尖峰。如果預(yù)先選擇一個(gè)確定性帶參函數(shù)，則其在空間的基本形狀是確定的。以正態(tài)分布密度函數(shù)為例，其有且只有一個(gè)尖峰，調(diào)節(jié)其參數(shù)都只能改變尖峰的位置和寬度，而無法改變尖峰的個(gè)數(shù)和偏度。因此，選擇確定性帶參函數(shù)從理論上就失去了無限逼近最優(yōu)IS-PDF 的可能性。而GMM 作為一個(gè)非參數(shù)模型，只要給予足夠分量，理論上能擬合各種邊際分布和相關(guān)性結(jié)構(gòu)，故其存在無限逼近最優(yōu)IS-PDF 的可能。

②從與LHS 法結(jié)合的難易度角度：如1.3 節(jié)所述，LHS 法需要進(jìn)行相關(guān)性控制，且控制過程中，變量間相關(guān)系數(shù)矩陣是一個(gè)重要參數(shù)，但大多數(shù)確定性帶參的高維分布函數(shù)無法顯式地提取相關(guān)系數(shù)矩陣。而GMM 的每個(gè)分量都顯式包含相關(guān)系數(shù)矩陣。因此，相較其他分布，選擇GMM 為IS-PDF 的基本形式，更容易讓IS 法與LHS 法結(jié)合。

然而，即便GMM 擁有以上優(yōu)點(diǎn)，其與LHS 法的結(jié)合仍不是無縫銜接。這是由于GMM 的每個(gè)分量都包含一個(gè)相關(guān)系數(shù)矩陣，而LHS 法進(jìn)行相關(guān)性控制時(shí)，只需要一個(gè)相關(guān)系數(shù)矩陣。因此，兩者的結(jié)合還需經(jīng)過特定的設(shè)計(jì)才能達(dá)成，這也是本文所提NSIS 法的算法核心所在。

2.2 NSIS 法流程

NSIS 法也按照IS 法慣例分為2 個(gè)階段，其一般性算法流程如下。其中，步驟3 中利用LHS 法采樣GMM 的過程是NSIS 法的核心步驟，實(shí)現(xiàn)了IS 法與LHS 法的結(jié)合。

步驟1：開始NSIS 法的第1 階段。記小概率事件中，xcv包含的變量個(gè)數(shù)為Kcv。以連續(xù)變量歷史數(shù)據(jù)為原始樣本，利用赤池信息準(zhǔn)則和最大期望（EM）算法［25］，估計(jì)GMM 分量個(gè)數(shù)nGM以及各分量參數(shù)，即ωk、μk和Σk。

步驟2：從每個(gè)Σk中分離出標(biāo)準(zhǔn)差參數(shù)向量σk和線性相關(guān)系數(shù)矩陣ρk，即

先用LHS 法對(duì)所有離散變量xdv采樣Nfir次，得到樣本矩陣Xdv，共Nfir列，每列即為xdv的一個(gè)樣本，其行數(shù)即為離散變量個(gè)數(shù)。

再用LHS 法對(duì)GMM 采樣Nfir次，具體方法如下。

1）按照GMM 各分量的權(quán)重分配每個(gè)分量采樣的樣本數(shù)，即第k(k=1，2，…，nGM)個(gè)分量需采樣ωkNfir次。初始化k=1。

2）對(duì)于第k 個(gè)分量，利用μk和σk中包含的信息，構(gòu)建Kcv個(gè)連續(xù)變量的邊際IS-PDF（均為一維正態(tài)分布函數(shù)），即

式中：Xk，e為第e 個(gè)變量利用GMM 第k 個(gè)分量采樣出的ωkNfir個(gè)樣本組成的向量；P 為一個(gè)行向量，由連續(xù)正整數(shù)序列{1，2，…，ωkNfir}打亂順序所得；ε 為與P 維數(shù)相同的隨機(jī)數(shù)向量，其每個(gè)元素都服從區(qū)間（0，1）上的均勻分布且互相獨(dú)立；Φ-1k，e(·)表示φk，e(·)對(duì)應(yīng)的累積分布函數(shù)的反函數(shù)。

Kcv個(gè)變量全部采樣完畢后，將Xk，e（e=1，2，…，Kcv）組合成樣本矩陣Xk，即

步驟4：根據(jù)交叉熵IS 法思想，采用所得樣本更新IS-PDF 參數(shù)集{ωk，μk，Σk|k=1，2，…，nGM}。更新公式將在第3 章電力系統(tǒng)可靠性評(píng)估的具體應(yīng)用場景下與其他離散變量的參數(shù)更新公式一并給出。

步驟5：判斷樣本矩陣X 中故障樣本占比是否滿足要求，若是，轉(zhuǎn)步驟6，否則轉(zhuǎn)步驟2。

步驟6：進(jìn)入NSIS 法的第2 階段。將第2 階段采樣樣本總數(shù)記為Nsec。同樣按照式（7）—式（10）所述方法對(duì)當(dāng)前IS-PDF 分層采樣Nsec個(gè)樣本，并采用該樣本按式（4）最后一行估計(jì)待求指標(biāo)并輸出。

3 NSIS 法在含可再生能源電力系統(tǒng)可靠性評(píng)估中的應(yīng)用

在含可再生能源的電力系統(tǒng)中，系統(tǒng)狀態(tài)x=[ xdv，xcv]可實(shí)例化為：

式中：xtg為系統(tǒng)中各火電機(jī)組狀態(tài)組成的向量；xhg為系統(tǒng)中各水電機(jī)組狀態(tài)組成的向量；xtl為系統(tǒng)中各輸電線路狀態(tài)組成的向量；Pwf為系統(tǒng)中各風(fēng)電場最大可用功率組成的向量；Lnet為系統(tǒng)各負(fù)荷節(jié)點(diǎn)的凈負(fù)荷組成的向量，凈負(fù)荷的具體含義將在介紹Lnet的概率分布模型時(shí)給出。

式（1）在系統(tǒng)缺電概率（loss of load probability，LOLP）指標(biāo)計(jì)算這一場合下實(shí)例化為：

3.1 含可再生能源電力系統(tǒng)可靠性評(píng)估基本模型

在計(jì)算可靠性指標(biāo)之前，需要明確式（18）中各函數(shù)的具體表達(dá)式，包括x 中各分量的原始PDF 和電力系統(tǒng)狀態(tài)評(píng)價(jià)函數(shù)S（x）。

3.1.1 電力系統(tǒng)變量的原始PDF 模型

下面針對(duì)xtg、xhg、xtl的PDF 以及Pwf與Lnet的聯(lián)合PDF 分別說明。

1）火電機(jī)組狀態(tài)變量的PDF

xtg包含多個(gè)元素，可表示為：

由于水電機(jī)組任意時(shí)刻的最大可能出力受全年不同時(shí)期來水量和水頭影響，因此需采用多狀態(tài)模型建模。xhg，j的不同取值代表水電機(jī)組的最大可能出力處于不同水平，各狀態(tài)具體數(shù)值可通過對(duì)水電機(jī)組的歷史實(shí)際出力數(shù)據(jù)進(jìn)行K-means 聚類得到。

式中：Lb為第b 個(gè)節(jié)點(diǎn)的負(fù)荷需求；Ppv，b為第b 個(gè)節(jié)點(diǎn)所接光伏電站的出力。

根據(jù)第b 個(gè)節(jié)點(diǎn)的光伏出力和負(fù)荷的歷史數(shù)據(jù)，利用式（26）得到Lnet，b的歷史數(shù)據(jù)。在構(gòu)建聯(lián)合PDF 時(shí)，將Lnet，b視作隨機(jī)變量，而不再單獨(dú)考慮Lb和Ppv，b。

對(duì)光伏采用凈負(fù)荷處理方式的理由為：與風(fēng)電不同，光伏與負(fù)荷的功率隨時(shí)間的變化均具有較為固定的模式，如光伏夜間出力為0，晴天時(shí)最大出力出現(xiàn)在12：00—14：00 之間，而夜間負(fù)荷也處于低谷期，峰荷則一般出現(xiàn)在傍晚。為了保留形狀相對(duì)固定的負(fù)荷和光伏出力日內(nèi)曲線在時(shí)間上的對(duì)應(yīng)關(guān)系，本文選擇文獻(xiàn)［30］中的凈負(fù)荷法。在此基礎(chǔ)上采用Copula 思想建立各風(fēng)電場出力與各節(jié)點(diǎn)凈負(fù)荷的聯(lián)合PDF。本文采用文獻(xiàn)［6］中提出的R-vine Copula 函數(shù)法擬合相關(guān)性結(jié)構(gòu)，采用文獻(xiàn)［31］中介紹的GMM 法擬合各變量的邊際分布。由此得到聯(lián)合PDF 表達(dá)式如下。

式中：crv(·)表示R-vine Copula 密度函數(shù)；F(Pwf)和F(Lnet)分別為Pwf和Lnet的累積分布函數(shù)；f (Lnet，b)表示第b 個(gè)節(jié)點(diǎn)的凈負(fù)荷的原始邊際PDF；f (Pwf，w)表示第w 座風(fēng)電場出力的原始邊際PDF。R-vine Copula 密度函數(shù)和GMM 的參數(shù)估計(jì)方法可分別參考文獻(xiàn)［6］和文獻(xiàn)［31］，此處不再重復(fù)說明。

3.1.2 電力系統(tǒng)狀態(tài)評(píng)價(jià)函數(shù)定義

S（x）的定義需依托最小切負(fù)荷模型的目標(biāo)函數(shù)值，其約束條件包括潮流約束、節(jié)點(diǎn)電壓約束、輸電線路傳輸功率約束等。

在電力系統(tǒng)規(guī)劃和可靠性分析中，常采用基于直流潮流的最小切負(fù)荷模型，其數(shù)學(xué)表達(dá)式已在多篇文獻(xiàn)中給出［18-19，25］，此處不再贅述，只給出該模型的目標(biāo)函數(shù)，其表達(dá)式為：

式中：Gi為第i 臺(tái)火電機(jī)組的裝機(jī)容量。

需要說明的是，在NSIS 法第2 階段正式估計(jì)可靠性指標(biāo)時(shí)，Sth=0。但在第1 階段求IS-PDF 參數(shù)時(shí)，Sth會(huì)動(dòng)態(tài)調(diào)整，以加快參數(shù)尋優(yōu)的速度，調(diào)整方式將在下節(jié)介紹。

3.2 基于NSIS 法的電力系統(tǒng)可靠性評(píng)估

3.2.1 NSIS 法計(jì)算電力系統(tǒng)可靠性指標(biāo)過程

1）選取IS-PDF 的函數(shù)形式

對(duì)于離散變量xtg、xhg和xtl，其IS-PDF 選取與原分布相同的函數(shù)形式，即

式中：ηtg，i、ηhg，j，m（m=1，2，…，o( j)）和ηtl，l分別為第i臺(tái)火電機(jī)組、第j 臺(tái)水電機(jī)組和第l 條輸電線路需要優(yōu)化的參數(shù)。

對(duì)于連續(xù)變量Pwf和Lnet，按照式（6）選擇GMM作為聯(lián)合IS-PDF 的基本函數(shù)形式，即

2）優(yōu)化IS-PDF 參數(shù)

在確定各變量的基本函數(shù)形式后，即可對(duì)ISPDF 中的相關(guān)參數(shù)進(jìn)行優(yōu)化。為了提高參數(shù)優(yōu)化效率，可以引入迭代過程，在迭代中動(dòng)態(tài)調(diào)整Sth值進(jìn)行參數(shù)的逐步更新，其具體步驟如下：

步驟1：按式（7）—式（16）所述方法對(duì)當(dāng)前ISPDF 采樣Nfir次。初次迭代時(shí)，采用原始PDF 進(jìn)行采樣。

步驟2：對(duì)每個(gè)樣本計(jì)算其S（x）數(shù)值，從小到大排序得到序列{S(1)，S(2)，…，S(Nfir)}，提取排在第αNfir位的數(shù)值S(αNfir)，其中α 為處于缺電狀態(tài)的樣本比例。

步驟3：如果S(αNfir)≤0，說明以當(dāng)前IS-PDF采樣得到的樣本中，能使系統(tǒng)處于缺電狀態(tài)的樣本比例已達(dá)到α，可停止迭代并進(jìn)入第2 階段；反之，令Sth=S(αNfir)，據(jù)此更新IS-PDF 的參數(shù)，并返回步驟1。此處，α 的取值越大，IS-PDF 的參數(shù)越接近最優(yōu)值，但對(duì)應(yīng)的迭代次數(shù)、所消耗的樣本數(shù)量和計(jì)算時(shí)間也會(huì)增加。因此，為了兼顧IS-PDF 參數(shù)尋優(yōu)的代價(jià)和后續(xù)可靠性指標(biāo)計(jì)算的加速效果，α 的典型值可取0.05 或0.1［14］。

參數(shù)更新公式如下。

式中：下標(biāo)s 表示第s 個(gè)樣本對(duì)應(yīng)的值；rk，s為更新連續(xù)變量IS-PDF 參數(shù)時(shí)使用的中間變量。

另外，需要說明的是，IS-PDF 的參數(shù)更新必須按照式（37）—式（40）的順序進(jìn)行。式（37）中用到的ωk是來自上一輪迭代的結(jié)果，而式（38）—式（40）中用到的參數(shù)都是本輪迭代更新過后的數(shù)值。如果第k 個(gè)分量的權(quán)重ωk小于1/Nfir，則刪除該分量。

3）利用IS-PDF 估計(jì)可靠性指標(biāo)

記第2 階段正式估計(jì)可靠性指標(biāo)時(shí)的樣本容量為Nsec，將Sth固定為0。LOLP 的樣本估計(jì)式為：

此外，電力系統(tǒng)各類非時(shí)序可靠性指標(biāo)都可由本文所提算法進(jìn)行估計(jì)，包括電力不足期望值（expected power not supplied，EPNS）γEPNS、切負(fù)荷持續(xù)時(shí)間期望值（expected duration of load curtailments，EDLC）γEDLC、電量不足期望值（expected energy not supplied，EENS）γEENS等，其樣本估計(jì)的表達(dá)式與LOLP 的表達(dá)式十分相似，具體可表示為［32］：

式中：T 為統(tǒng)計(jì)可靠性指標(biāo)的時(shí)段長度，常取1 年，即8 760 h。

若需要用方差系數(shù)判斷上述可靠性指標(biāo)估計(jì)的精度，則需計(jì)算數(shù)次可靠性指標(biāo)，并將每一次結(jié)果視作一個(gè)樣本，從而直接計(jì)算可靠性指標(biāo)的樣本方差和方差系數(shù)，對(duì)應(yīng)公式如下。

式中：β 為可靠性指標(biāo)的方差系數(shù)；γz表示重復(fù)計(jì)算過程中，第z 次計(jì)算得到的可靠性指標(biāo)的數(shù)值；nγ為可靠性指標(biāo)重復(fù)計(jì)算的次數(shù)。

NSIS 法計(jì)算含可再生能源電力系統(tǒng)的靠性指標(biāo)的算法流程如附錄A 所示。

3.2.2 計(jì)及水文特性的電力系統(tǒng)可靠性指標(biāo)算法設(shè)計(jì)

如果水電機(jī)組來水量和水頭在豐水期和枯水期差距較大，可按如下方法計(jì)算可靠性指標(biāo)：

步驟1：將豐水期和枯水期占全年總時(shí)長的比例記為pwet和pdry。

步驟2：利用3.2.1 節(jié)所述NSIS 法分別計(jì)算在豐水期和枯水期電力系統(tǒng)的可靠性指標(biāo)，并記為γwet和γdry。

步驟3：對(duì)豐水期和枯水期的可靠性指標(biāo)進(jìn)行加權(quán)求和，得到全年可靠性指標(biāo)γall，即

其中，在豐水期和枯水期時(shí)段內(nèi)利用NSIS 法計(jì)算可靠性指標(biāo)的流程區(qū)別只在于：隨機(jī)變量xhg、Pwf和Lnet的原始PDF 參數(shù)是用各自時(shí)段內(nèi)的歷史數(shù)據(jù)擬合得到的。

4 算例分析

4.1 算例環(huán)境和基本參數(shù)設(shè)置

硬件環(huán)境：CPU 為Intel core i7-5500U、內(nèi)存為16 GB 的微機(jī)。軟件環(huán)境：MATLAB 2016a。

采用IEEE RTS-79 標(biāo)準(zhǔn)測試系統(tǒng)的基本網(wǎng)架結(jié)構(gòu)和數(shù)據(jù)作為基礎(chǔ)［33］，在該系統(tǒng)的節(jié)點(diǎn)1、2 和13處分別接入一座風(fēng)電場，在節(jié)點(diǎn)18 處接一座光伏電站，構(gòu)成改進(jìn)IEEE RTS-79 系統(tǒng)。

根據(jù)測試系統(tǒng)的數(shù)據(jù)可知，所有離散元件的狀態(tài)中，發(fā)生概率的最小值為pdv，min=3.423 5×10-4，因此建議NSIS 法中第1 階段樣本數(shù)Nfir≥29 210。故可取Nfir=30 000。此外，令反映第1 階段缺電樣本比例要求的參數(shù)α=0.05。

在標(biāo)準(zhǔn)IEEE RTS-79 系統(tǒng)中，所有節(jié)點(diǎn)負(fù)荷完全正相關(guān)，但節(jié)點(diǎn)18 接入光伏電站后，其凈負(fù)荷與其他節(jié)點(diǎn)負(fù)荷不會(huì)繼續(xù)完全正相關(guān)。因此，本文將除節(jié)點(diǎn)18 外所有節(jié)點(diǎn)負(fù)荷之和視作一個(gè)隨機(jī)變量，將節(jié)點(diǎn)18 的凈負(fù)荷單獨(dú)視作一個(gè)隨機(jī)變量，再加上3 座風(fēng)電場的隨機(jī)出力，故本文算例共設(shè)計(jì)有5 個(gè)連續(xù)隨機(jī)變量。

風(fēng)電場、光伏電站的某年度歷史出力數(shù)據(jù)見文獻(xiàn)［34］。其中，3 座風(fēng)電場的原始邊際PDF 見圖1。

圖1 3 座風(fēng)電場原始邊際PDF 的GMM 擬合結(jié)果Fig.1 GMM fitting results of original marginal PDF for three wind farms

4.2 歷史數(shù)據(jù)相關(guān)性控制

后續(xù)算例中將會(huì)討論風(fēng)電場間相關(guān)性、風(fēng)電出力與負(fù)荷的相關(guān)性對(duì)電力系統(tǒng)可靠性的影響，因此需要人為給定和改變連續(xù)變量的相關(guān)系數(shù)矩陣。具體操作方式如下。

首先，假設(shè)各風(fēng)電場出力變量間的相關(guān)系數(shù)相同，均為ρW；各風(fēng)電場出力與各負(fù)荷變量的相關(guān)系數(shù)也相同，均為ρLW；2 個(gè)負(fù)荷變量的相關(guān)系數(shù)記為ρL。按照風(fēng)電場1、風(fēng)電場2、風(fēng)電場3、節(jié)點(diǎn)18 凈負(fù)荷、系統(tǒng)內(nèi)剩余負(fù)荷的順序排列5 個(gè)連續(xù)隨機(jī)變量，則其線性相關(guān)系數(shù)矩陣ρ 可表示為：

然后，參考式（13）和式（14）所述的相關(guān)性控制方法，調(diào)整歷史數(shù)據(jù)的排序，即可在不改變各連續(xù)變量的邊際PDF 的前提下，讓其擁有指定的相關(guān)系數(shù)矩陣。隨后，這些經(jīng)過相關(guān)性控制的歷史數(shù)據(jù)就作為式（27）的輸入數(shù)據(jù)，用于建立電力系統(tǒng)連續(xù)變量的聯(lián)合PDF。

4.3 NSIS 法的正確性驗(yàn)證

為驗(yàn)證NSIS 法計(jì)算結(jié)果的正確性，需將其與其他方法的計(jì)算結(jié)果作對(duì)比。本節(jié)將基于簡單隨機(jī)采樣的MCS 法記為方法1，LHS 法記為方法2，文獻(xiàn)［13］中提出的基于交叉熵的IS 法記為方法3。

采用NSIS 法與上述幾種方法來計(jì)算改進(jìn)IEEE RTS-79 系統(tǒng)的LOLP、EDLC、EPNS、EENS等可靠性指標(biāo)。對(duì)于方法3，將其第1 階段（預(yù)抽樣）的樣本容量設(shè)置為30 000，與NSIS 法保持一致。將所有方法的收斂條件均設(shè)置為所計(jì)算的可靠性指標(biāo)的方差系數(shù)小于等于0.02。各方法重復(fù)50 次后的計(jì)算結(jié)果平均值如表1 所示。

由表1 可知，在相同收斂條件下，NSIS 法與其他3 種方法計(jì)算出的改進(jìn)IEEE RTS-79 系統(tǒng)的4 種可靠性指標(biāo)基本一致，進(jìn)一步驗(yàn)證了NSIS 法的正確性。

表1 4 種方法評(píng)估改進(jìn)IEEE RTS-79 系統(tǒng)可靠性的結(jié)果對(duì)比Table 1 Result comparison of four methods for evaluating reliability of improved IEEE RTS-79 system

4.4 NSIS 法的性能驗(yàn)證

4.3 節(jié)提到的4 種方法按照方差系數(shù)小于等于0.02 的要求重復(fù)計(jì)算50 次可靠性指標(biāo)，其消耗樣本和計(jì)算時(shí)間的平均值如表2 所示。

表2 4 種方法評(píng)估改進(jìn)IEEE RTS-79 系統(tǒng)可靠性的計(jì)算成本對(duì)比Table 2 Comparison of calculation costs of four methods for evaluating reliability of improved IEEE RTS-79 system

由表2 可知，在同樣的精度要求下，NSIS 法消耗的樣本數(shù)和計(jì)算時(shí)間遠(yuǎn)小于MCS 法和LHS 法。由于包含了分層抽樣思想，在計(jì)算4 種可靠性指標(biāo)時(shí)，NSIS 法相較于單純的交叉熵IS 法，計(jì)算用時(shí)縮短了20.3%～30.4% 不等，計(jì)算效率得到了顯著提升。

此外，為驗(yàn)證NSIS 法在不同可靠性水平的電力系統(tǒng)中的適應(yīng)性，將上述4 種方法的樣本容量均設(shè)置為80 000 次，對(duì)LOLP 指標(biāo)的方差系數(shù)進(jìn)行關(guān)于負(fù)荷水平（即系統(tǒng)可靠性水平）的靈敏度分析。其中，根據(jù)4.3 節(jié)計(jì)算過程可知，方法3 和NSIS 法在計(jì)算該測試系統(tǒng)時(shí)，需要在預(yù)抽樣階段迭代2 次，共消耗2Nfir=60 000 次樣本。此外，方法1 和方法3 采用式（3）計(jì)算方差系數(shù)，方法2 和NSIS 法采用式（45）計(jì)算方差系數(shù)。結(jié)果如圖2 和圖3 所示。

由圖2 可知，隨著負(fù)荷水平的降低，系統(tǒng)內(nèi)發(fā)電資源的充裕度將逐漸提高，在可靠性水平上表現(xiàn)為LOLP 絕對(duì)值下降，即系統(tǒng)缺電這一小概率事件發(fā)生的概率進(jìn)一步降低。

圖2 LOLP 關(guān)于負(fù)荷水平的關(guān)系圖Fig.2 Relation diagram of LOLP with respect to load level

由圖3 可知，在LOLP 數(shù)值下降的情況下，各方法采樣20 000 次后的計(jì)算結(jié)果的方差系數(shù)均逐漸變大，但NSIS 法相較于其他方法而言，方差系數(shù)總是最小的（即精度最高），且受可靠性水平的影響也最不明顯。這說明NSIS 法比其他3 種方法更適合用來評(píng)估現(xiàn)代高可靠性電力系統(tǒng)的可靠性。

圖3 方差系數(shù)關(guān)于負(fù)荷水平的關(guān)系曲線Fig.3 Relation curves of coefficient of variance with respect to load level

4.5 變量相關(guān)性對(duì)系統(tǒng)可靠性水平的影響

儲(chǔ)能設(shè)備和需求響應(yīng)技術(shù)的大規(guī)模應(yīng)用，可實(shí)現(xiàn)源荷匹配，促進(jìn)新能源消納。這一背景下，含儲(chǔ)能的新能源場站出力與計(jì)及需求響應(yīng)的負(fù)荷之間的相關(guān)性水平可能會(huì)有明顯提升，因此，需分析源荷相關(guān)性水平對(duì)系統(tǒng)可靠性水平的影響規(guī)律。

將式（47）中的ρL固定為0.6，讓?duì)袻W和ρW在區(qū)間［0，0.6］內(nèi)變化，并求出對(duì)應(yīng)的LOLP 值。繪制三者的函數(shù)關(guān)系如圖4 所示。需要說明的是，圖4 中留白部分的任意一點(diǎn)的坐標(biāo)代入式（47）后得到的矩陣都是非正定的，而相關(guān)系數(shù)矩陣必須是正定的，故留白區(qū)域無法進(jìn)行可靠性計(jì)算，也就沒有對(duì)應(yīng)的LOLP 值。

圖4 各變量相關(guān)性水平對(duì)系統(tǒng)可靠性的影響Fig.4 Influence of correlation level of each variable on system reliability

由圖4 可知，風(fēng)電場出力與負(fù)荷的相關(guān)性越高，系統(tǒng)的可靠性越高。而各風(fēng)電場出力的相關(guān)性越高，系統(tǒng)的可靠性越低。

4.6 計(jì)及水文特性的可靠性指標(biāo)評(píng)估

本節(jié)中將改進(jìn)IEEE RTS-79 測試系統(tǒng)中節(jié)點(diǎn)22 所接的6 臺(tái)裝機(jī)容量為50 MW 的火電機(jī)組替換為相同容量的水電機(jī)組，利用3.2.2 節(jié)所述方法計(jì)算計(jì)及水文特性的可靠性評(píng)估指標(biāo)。

首先，根據(jù)中國青海地區(qū)實(shí)際水電出力數(shù)據(jù)，利用K-means 方法對(duì)該地區(qū)豐水期、枯水期水電機(jī)組出力水平進(jìn)行聚類。其中，青海地區(qū)水電豐水期為每年4 月1 日至10 月31 日，剩余時(shí)間均屬于枯水期，故式（46）中的pwet和pdry分別為0.586 3 和0.413 7。

在豐水期、枯水期時(shí)段內(nèi)將水電機(jī)組出力數(shù)據(jù)分別聚類為3 類，即對(duì)應(yīng)3 個(gè)出力水平。此外，還應(yīng)考慮機(jī)組故障狀態(tài)下的出力為0。因此，水電機(jī)組的出力水平在豐、枯水期均可用4 種狀態(tài)模型表征，其數(shù)值和對(duì)應(yīng)的概率大小如表3 所示。

表3 豐水期、枯水期水電機(jī)組出力水平和對(duì)應(yīng)概率Table 3 Power output levels and corresponding probabilities of hydro power units during wet and dry periods

表3 中出力水平以標(biāo)幺值形式給出，對(duì)應(yīng)的基準(zhǔn)值為水電機(jī)組的裝機(jī)容量。按照3.2.2 節(jié)所述方法，分別采用NSIS 法計(jì)算豐、枯水期的LOLP 指標(biāo)，其數(shù)值分別為8.587 0×10-5和1.268 3×10-3。利用式（46）求取電力系統(tǒng)全年LOLP 指標(biāo)為5.750 4×10-4。與表1 對(duì)比可知，由于火電機(jī)組被含有多個(gè)降額狀態(tài)的水電機(jī)組代替，系統(tǒng)發(fā)電資源可用容量必然減少，可靠性下降，故LOLP 值增大。

5 結(jié)語

本文以含可再生能源電力系統(tǒng)的可靠性指標(biāo)計(jì)算方法為研究對(duì)象，提出了一種非參數(shù)重要性分層抽樣法。該方法考慮到IS 法與LHS 法優(yōu)缺點(diǎn)互補(bǔ)的特性，利用兩者的結(jié)合實(shí)現(xiàn)了電力系統(tǒng)可靠性指標(biāo)計(jì)算的進(jìn)一步加速。

本文首先分析了IS 法與LHS 法結(jié)合時(shí)面臨的矛盾和困難。然后，選擇GMM 作為非參數(shù)形式的IS-PDF，利用GMM 各分量均為易于提取相關(guān)系數(shù)矩陣的多元正態(tài)分布的特點(diǎn)，為GMM 定制了計(jì)及相關(guān)性控制的LHS 法采樣方案，從而在不犧牲IS 法與LHS 法各自效率的同時(shí)，實(shí)現(xiàn)了兩者的結(jié)合，即本文所提出的NSIS 法。設(shè)計(jì)算例驗(yàn)證了所提NSIS 法在含可再生能源電力系統(tǒng)可靠性評(píng)估中的正確性和有效性，相關(guān)結(jié)論如下：

1）將NSIS 法和MCS 法、LHS 法以及基于交叉熵的IS 法同時(shí)用于計(jì)算測試系統(tǒng)的可靠性指標(biāo)。結(jié)果顯示NSIS 法計(jì)算結(jié)果與其他3 種算法基本一致，但效率高于其他3 種算法，說明了NSIS 法的正確性和有效性。

2）NSIS 法的精度受可靠性指標(biāo)絕對(duì)值大小的影響較小，因此適用于高可靠性電力系統(tǒng)。

3）新能源出力與負(fù)荷的正相關(guān)性有助于系統(tǒng)可靠性的提升，而新能源出力之間的相關(guān)性過高則會(huì)損害系統(tǒng)可靠性。因此，可從這一結(jié)論出發(fā)，利用大規(guī)模儲(chǔ)能技術(shù)和需求響應(yīng)技術(shù)控制新能源出力間、新能源出力與負(fù)荷的相關(guān)性，從而提升系統(tǒng)的可靠性。

本文在撰寫過程中得到江蘇省配電網(wǎng)智能技術(shù)與裝備協(xié)同創(chuàng)新中心開放基金資助項(xiàng)目(XTCX202001)和南京工業(yè)大學(xué)引進(jìn)人才科研啟動(dòng)專項(xiàng)經(jīng)費(fèi)(39810127)資助，特此感謝！

附錄見本刊網(wǎng)絡(luò)版（http：//www.aeps-info.com/aeps/ch/index.aspx），掃英文摘要后二維碼可以閱讀網(wǎng)絡(luò)全文。