基于Nataf變換和準(zhǔn)蒙特卡洛模擬的隨機(jī)潮流方法

方斯頓，程浩忠，徐國(guó)棟，曾平良，姚良忠

（1.上海交通大學(xué) 電力傳輸與功率變換控制教育部重點(diǎn)實(shí)驗(yàn)室，上海 200240；2.中國(guó)電力科學(xué)研究院，北京 100192）

0 引言

近年來(lái)，電力市場(chǎng)化改革[1]和以風(fēng)電為代表的可再生能源接入[2]給電力系統(tǒng)運(yùn)行和規(guī)劃帶來(lái)更多不確定性，隨機(jī)潮流PLF（Probabilistic Load Flow）計(jì)及這些不確定因素的影響，可綜合分析系統(tǒng)的運(yùn)行風(fēng)險(xiǎn)和薄弱節(jié)點(diǎn)，因而受到廣泛關(guān)注[3-5]。

PLF 首先由 Borkowaska 于1974 年提出[6]，40 年來(lái)，國(guó)內(nèi)外學(xué)者對(duì)其進(jìn)行了大量?jī)?yōu)化與改進(jìn)，目前主要可分為解析法[7-10]和模擬法[11-13]兩大類(lèi)。解析法可分為卷積法[7]、半不變量[8-9]、點(diǎn)估計(jì)[10]。 卷積法利用輸入變量分布函數(shù)的卷積運(yùn)算得到輸出變量的概率分布，計(jì)算量大。半不變量通過(guò)線(xiàn)性化潮流方程建立的線(xiàn)性關(guān)系將卷積運(yùn)算轉(zhuǎn)化為半不變量的代數(shù)運(yùn)算，計(jì)算效率高但準(zhǔn)確性差，得到的變量概率分布在首尾可能出現(xiàn)小于零的部分。點(diǎn)估計(jì)是已知輸入變量的若干階矩信息求取輸出變量矩特征的數(shù)學(xué)方法，計(jì)算效率高，但高階矩誤差較大。

除卷積法外，解析法計(jì)算效率高，但誤差均較大。而模擬法可準(zhǔn)確地對(duì)實(shí)際物理過(guò)程建模，且建模質(zhì)量不受變量分布性質(zhì)、系統(tǒng)非線(xiàn)性的影響，因而比其他PLF方法更準(zhǔn)確，其中蒙特卡洛模擬MCS（Monte Carlo Simulation）應(yīng)用最為廣泛?；诤?jiǎn)單隨機(jī)抽樣SRS（Simple Random Sampling）的 MCS 收斂慢，獲得準(zhǔn)確解的計(jì)算代價(jià)大，因此，如何在保證精度的同時(shí)提高計(jì)算效率一直是MCS的研究熱點(diǎn)。研究表明，樣本的差異性[14]是影響MCS計(jì)算精度的主要因素。文獻(xiàn)[11-13]采用各種基于拉丁超立方采樣LHS（Latin Hypercube Sampling）的方法改進(jìn) MCS，但由于無(wú)法保證序列的低差異性，因此無(wú)法克服MCS收斂性的瓶頸。而以Sobol序列為代表的低差異序列LDSs（Low-Discrepancy Sequences）MCS 在計(jì)算精度和效率上均優(yōu)于LHS，已在電子電路設(shè)計(jì)中取得了應(yīng)用[14]，但在隨機(jī)潮流中的應(yīng)用還沒(méi)有見(jiàn)諸報(bào)道。

另一方面，為得到含相關(guān)性的樣本序列，隨機(jī)排序[15]、Cholesky 分解[16]、遺傳算法[17]等改進(jìn)排序方法均在MCS中得到應(yīng)用。隨機(jī)排序在樣本較小時(shí)效果差；Cholesky分解僅針對(duì)相關(guān)系數(shù)矩陣正定的情況；遺傳算法有迭代操作，結(jié)果較準(zhǔn)確但計(jì)算耗時(shí)長(zhǎng)。因此，本文提出采用奇異值分解的方法對(duì)樣本進(jìn)行排序，計(jì)算量與Cholesky分解相當(dāng)，并可處理相關(guān)系數(shù)矩陣非正定的情況。

為提高M(jìn)CS的計(jì)算效率并處理相關(guān)系數(shù)矩陣非正定的情況，本文提出一種基于Nataf變換和準(zhǔn)蒙特卡洛模擬NQM（Nataf transformation based Quasi Monte Carlo simulation method）的PLF方法。首先利用Sobol方法產(chǎn)生服從均勻分布的LDSs，而后利用基于奇異值分解的Nataf變換生成服從輸入變量聯(lián)合分布的序列，再經(jīng)過(guò)若干次確定型潮流計(jì)算得到輸出變量的概率分布特性。對(duì)IEEE 30節(jié)點(diǎn)系統(tǒng)和我國(guó)某實(shí)際大區(qū)域電網(wǎng)的仿真驗(yàn)證了本文方法的有效性。

1 基于奇異值分解的Nataf變換

Nataf變換是在已知輸入變量邊緣分布的情況下重構(gòu)聯(lián)合分布的數(shù)學(xué)方法[18]。現(xiàn)將其過(guò)程簡(jiǎn)述如下：對(duì)任意 n 維輸入變量 X= [X1，…，Xi，…，Xn]T，設(shè)其相關(guān)系數(shù)矩陣為[ρXij]n×n，其中 ρXij定義如式（1）所示。

其中，σi和 σj分別為 Xi、Xj的標(biāo)準(zhǔn)差；cov（Xi，Xj）為變量Xi、Xj的協(xié)方差；ρXij為相關(guān)系數(shù)。隨機(jī)向量X可用Nataf變換轉(zhuǎn)化為標(biāo)準(zhǔn)正態(tài)分布向量Z=[Z1，…，Zi，…，Zn]T：

其中，F(xiàn)i為Xi的累積分布函數(shù)；Φ為標(biāo)準(zhǔn)正態(tài)分布的累積分布函數(shù)。設(shè)變換以后Z的相關(guān)系數(shù)矩陣為[ρZij]n×n，由 Nataf變換性質(zhì)，[ρXij]n×n與[ρZij]n×n可根據(jù)式（3）相互轉(zhuǎn)化：

其中，fXiXj（Xi，Xj）為 Xi、Xj的聯(lián)合概率密度函數(shù)；μ、σ分別為相應(yīng)變量的期望值與標(biāo)準(zhǔn)差；φZ(yǔ)iZj（Zi，Zj，ρZij）為相關(guān)系數(shù)為ρZij的二維標(biāo)準(zhǔn)正態(tài)分布概率密度函數(shù)。本文采用文獻(xiàn)[19]中方法對(duì)式（3）進(jìn)行求解以求得[ρZij]n×n。 而由于[ρZij]n×n可能為非正定或非滿(mǎn)秩陣，此時(shí)其Cholesky分解不存在，無(wú)法采用文獻(xiàn)[16]的方法對(duì)樣本進(jìn)行排序。但一般相關(guān)系數(shù)矩陣都為對(duì)稱(chēng)陣，存在奇異值分解，本文通過(guò)證明定理1說(shuō)明奇異值分解排序方法的過(guò)程。

定理1 設(shè)K為n維獨(dú)立標(biāo)準(zhǔn)正態(tài)分布向量，矩陣UρZ∑1/2ρZ由向量Z的相關(guān)系數(shù)矩陣ρZ的奇異值分解生成：

其中，UρZ為酉矩陣；∑ρZ為由矩陣奇異值構(gòu)成的對(duì)角矩陣，并按從大到小排列。則由式（5）定義的Z*相關(guān)系數(shù)矩陣等于ρZ：

證明如下。

對(duì)?z*i，其期望值滿(mǎn)足式（6）：

相關(guān)系數(shù)矩陣為：

因此，可根據(jù)定理1的步驟生成相關(guān)系數(shù)為ρZ的正態(tài)分布向量Z*，再通過(guò)Nataf反變換，可得到服從任意分布的隨機(jī)向量X。

2 基于Sobol序列的準(zhǔn)蒙特卡洛方法

2.1 MCS及其誤差

在確定型潮流計(jì)算中，節(jié)點(diǎn)注入功率S與節(jié)點(diǎn)電壓U、支路潮流B關(guān)系如式（8）所示。

其中，S∈Ω為節(jié)點(diǎn)注入功率向量，Ω為節(jié)點(diǎn)注入功率空間；U為節(jié)點(diǎn)電壓向量；B為支路潮流向量；f、g分別為由潮流方程導(dǎo)出的節(jié)點(diǎn)電壓函數(shù)和支路潮流函數(shù)。而在PLF計(jì)算中，S、U和B均為隨機(jī)變量，Ω為隨機(jī)空間。在S的分布特性已知時(shí)，U和B的數(shù)字特征可由如式（9）所示的Stieltjes積分求取。

其中，EUk、EBk分別為 U 和 B 的 k階矩；H（S）為 S的累積分布函數(shù)。這類(lèi)積分當(dāng)f、g形式復(fù)雜時(shí)很難求解，而MCS提供了求取這一積分的數(shù)值方法，求取方法如式（10）所示。

其中，Si為根據(jù)節(jié)點(diǎn)注入功率分布產(chǎn)生的序列；n為MCS次數(shù)。 為便于分析，將式（10）簡(jiǎn)化為式（11）：

其中，Qn為n次模擬得到的數(shù)值結(jié)果；p代表式（9）中的被積函數(shù)，pi則代表式（10）中的 fk（Si）ΔH（Si）或 gk（f（Si））ΔH（Si）。 根據(jù)大數(shù)定律，MCS 估計(jì)得到的Qn會(huì)以概率1收斂到其真實(shí)值Q。

其中，P（·）為事件的概率。

積分的估計(jì)誤差則可由定理2給出。

定理2[20]如果p的方差有限，即：

則MCS誤差為：

由定理 2，當(dāng) p確定時(shí)，MCS以 O（n-1/2）的速度收斂于真實(shí)值。欲提高M(jìn)CS的精度，則可以通過(guò)減小σ（p）或增大n來(lái)實(shí)現(xiàn)。以L(fǎng)HS為代表的方差縮減技術(shù)可以減小σ（p），因此相比于簡(jiǎn)單隨機(jī)抽樣，LHS可以提高M(jìn)CS的精度，但從式（14）可以看出，其并不能改變收斂速度（仍然為O（n-1/2）），因而無(wú)法從本質(zhì)上提高M(jìn)CS的效率。下文將引入低差異序列改進(jìn)MCS的樣本生成方法，提高其收斂速度。

2.2 低差異性序列與準(zhǔn)蒙特卡洛

基于低差異性序列的MCS稱(chēng)為準(zhǔn)蒙特卡洛模擬 QMCS（Quasi Monte Carlo Simulation），目前已在電子電路設(shè)計(jì)中取得應(yīng)用[14]，精度和效率均優(yōu)于基于拉丁超立方的方法。本文采用該方法進(jìn)行PLF計(jì)算。首先，本文通過(guò)引入文獻(xiàn)[20]的結(jié)論說(shuō)明序列差異性與MCS精度的關(guān)系。

定理3[20]MCS誤差可定義為：

結(jié)合式（15），基于 LHS的 MCS收斂速度為O（n-1/2），與式（14）一致。 文獻(xiàn)[23]指出，任意維度均勻分布序列最小Dn*如式（18）所示。滿(mǎn)足該式的序列稱(chēng)為L(zhǎng)DSs。

對(duì)比式（17）和（18），序列的 L∞-差異性由 O（n-1/2）變?yōu)?O（n-1），再結(jié)合式（15），引入低差異序列后 MCS的收斂速度可由 O（n-1/2）變?yōu)?O（n-1），從而大幅加快MCS的收斂。

2.3 NQM誤差定義

為比較結(jié)合Nataf變換的SRS、LHS和NQM的誤差，本文以SRS計(jì)算50000次的結(jié)果作為標(biāo)準(zhǔn)，定義估計(jì)誤差為：

其中，εγs為相對(duì)誤差指標(biāo)；n為向量的維度；γ為輸出變量類(lèi)型，包括節(jié)點(diǎn)電壓幅值、電壓相角、支路有功潮流、支路無(wú)功潮流；s為數(shù)值特征，包括期望值μ與標(biāo)準(zhǔn)差 σ；αγsc為某采樣規(guī)模下得到的輸出結(jié)果；αγsm為SRS模擬50000次的結(jié)果。

因MCS收斂具有隨機(jī)性，因此本文對(duì)在每一采樣規(guī)模下均計(jì)算100次，將誤差的平均值εˉγs作為誤差分析指標(biāo)。

2.4 NQM計(jì)算流程

綜合前文，NQM的計(jì)算流程如圖1所示。簡(jiǎn)述如下：

a.輸入基礎(chǔ)數(shù)據(jù)，包括系統(tǒng)參數(shù)、節(jié)點(diǎn)注入功率的分布及其相關(guān)系數(shù)矩陣ρX、LDSs序列長(zhǎng)度；

b.利用式（3）將ρX轉(zhuǎn)化為標(biāo)準(zhǔn)正態(tài)向量的相關(guān)系數(shù) ρZ；

c.利用文獻(xiàn)[21]中方法生成0-1分布的獨(dú)立LDSs序列矩陣Un×M，并利用標(biāo)準(zhǔn)正態(tài)分布的逆函數(shù)將其轉(zhuǎn)化為正態(tài)分布的LDSs序列矩陣Dn×M；

d.利用步驟a中的方法將Dn×M轉(zhuǎn)化為相關(guān)系數(shù)矩陣等于ρZ的正態(tài)分布序列矩陣

f.根據(jù)Xn×M進(jìn)行M次確定型潮流計(jì)算，并統(tǒng)計(jì)得到各輸出變量的數(shù)字特征及概率分布。

圖1 計(jì)算流程圖Fig.1 Flowchart of calculation

3 算例分析

3.1 算例介紹

本文所提算法分別采用IEEE 30節(jié)點(diǎn)系統(tǒng)和我國(guó)某大區(qū)域電網(wǎng)500 kV等值網(wǎng)絡(luò)進(jìn)行驗(yàn)證，該網(wǎng)絡(luò)含節(jié)點(diǎn)1594個(gè)、發(fā)電機(jī)535臺(tái)、線(xiàn)路3359條。2個(gè)算例分別采用SRS、LHS和NQM進(jìn)行計(jì)算并分析誤差。仿真平臺(tái)基于MATLAB2013a，Intel Core dual i7-3820 3.6 GHz，RAM 8 GB。算例考慮的輸入變量包括風(fēng)電機(jī)組出力和負(fù)荷。假設(shè)風(fēng)速滿(mǎn)足威布爾分布 W（c，ξ）=W（10.7，3.97），風(fēng)電場(chǎng)功率因數(shù)為 0.9，對(duì)于 IEEE 30節(jié)點(diǎn)系統(tǒng)，節(jié)點(diǎn) 3、4、6、28分別接入容量為15 MW的風(fēng)電場(chǎng)；實(shí)際電網(wǎng)算例中的節(jié)點(diǎn)101—150均裝設(shè)300 MW的風(fēng)電場(chǎng)。

風(fēng)電機(jī)出力的互相關(guān)系數(shù)均為0.9。負(fù)荷均為恒功率因數(shù)0.9，有功負(fù)荷服從期望為額定有功功率、標(biāo)準(zhǔn)差為期望5%的正態(tài)分布，互相關(guān)系數(shù)均為0.5。有功出力與風(fēng)速關(guān)系滿(mǎn)足：

其中，PWR為風(fēng)電場(chǎng)額定功率；vci、vr和 vco分別為切入風(fēng)速、額定風(fēng)速和切出風(fēng)速，大小分別為2.5 m/s、13 m/s和 25 m/s。

3.2 NQM誤差分析

分別利用SRS、LHS和NQM這3種方法對(duì)2個(gè)測(cè)試系統(tǒng)進(jìn)行計(jì)算。按式（19）計(jì)算輸出變量誤差，分別表示為：節(jié)點(diǎn)電壓幅值期望誤差，標(biāo)準(zhǔn)差誤差；節(jié)點(diǎn)電壓相角期望誤差，標(biāo)準(zhǔn)差誤差；支路有功期望誤差，標(biāo)準(zhǔn)差誤差；支路無(wú)功潮流期望誤差，標(biāo)準(zhǔn)差誤差。 IEEE 30 節(jié)點(diǎn)系統(tǒng)輸出結(jié)果見(jiàn)圖2、3。實(shí)際網(wǎng)絡(luò)計(jì)算結(jié)果見(jiàn)表1。

由圖2、3及表1中數(shù)據(jù)可得：

a.當(dāng)采樣規(guī)模相近時(shí)，NQM和LHS的精度遠(yuǎn)高于SRS；

b.在輸出變量的誤差上，當(dāng)系統(tǒng)規(guī)模較小時(shí)，NQM和LHS求得的期望值結(jié)果相近，但NQM收斂速度快于LHS，而在輸出變量標(biāo)準(zhǔn)差的誤差上，NQM精度遠(yuǎn)高于同樣規(guī)模的LHS；

c.當(dāng)系統(tǒng)規(guī)模擴(kuò)大時(shí)，由NQM得到的變量期望值及標(biāo)準(zhǔn)差均優(yōu)于LHS，且收斂更快。這說(shuō)明NQM方法同樣適用于實(shí)際大系統(tǒng)。

圖2 IEEE 30節(jié)點(diǎn)系統(tǒng)電壓幅值相角誤差曲線(xiàn)Fig.2 Error curves of bus voltage magnitude and angle of IEEE 30-bus system

圖3 IEEE 30節(jié)點(diǎn)系統(tǒng)支路潮流誤差曲線(xiàn)Fig.3 Error curves of branch power flows of IEEE 30-bus system

表1 實(shí)際網(wǎng)絡(luò)平均誤差比較（N=1000）Table 1 Comparison of average errors of a practical power grid（N=1000）

用NQM還可以方便地得到輸出變量的概率分布曲線(xiàn)。選取IEEE 30節(jié)點(diǎn)系統(tǒng)支路4-6有功潮流及實(shí)際網(wǎng)絡(luò)中某線(xiàn)路無(wú)功潮流為研究對(duì)象，樣本規(guī)模為50000的SRS作為標(biāo)準(zhǔn)結(jié)果，LHS計(jì)算800次和NQM計(jì)算500次的隨機(jī)變量概率密度函數(shù)和累積分布函數(shù)比較分別如圖4、5所示。

從圖4、5中可得，無(wú)論系統(tǒng)大小，采樣規(guī)模為500的NQM得到了與SRS計(jì)算50000次相近的結(jié)果，明顯優(yōu)于LHS計(jì)算800次的結(jié)果。這再次說(shuō)明LDSs可克服蒙特卡洛收斂性的瓶頸，從而提高M(jìn)CS的精度，節(jié)省PLF的計(jì)算時(shí)間。

3.3 NQM收斂趨勢(shì)分析

為獲得3種方法對(duì)特定變量的收斂趨勢(shì)，本文以IEEE 30節(jié)點(diǎn)系統(tǒng)節(jié)點(diǎn)14的電壓和支路6-8有功功率、實(shí)際網(wǎng)絡(luò)某線(xiàn)路的無(wú)功功率為例；期望值與標(biāo)準(zhǔn)差隨采樣規(guī)模的變化趨勢(shì)如圖6所示，圖中電壓幅值均值、電壓幅值標(biāo)準(zhǔn)差為標(biāo)幺值。

圖6中數(shù)據(jù)說(shuō)明，3種方法的收斂趨勢(shì)是相同的，但NQM收斂明顯快于LHS和SRS，因此NQM可以在較短時(shí)間內(nèi)獲得較高的計(jì)算精度。

圖4 IEEE 30節(jié)點(diǎn)系統(tǒng)支路4-6有功潮流概率分布曲線(xiàn)Fig.4 Active power probability distribution curves of branch 4-6 of IEEE 30-bus system

圖5 實(shí)際系統(tǒng)某支路無(wú)功潮流概率分布曲線(xiàn)Fig.5 Reactive power probability distribution curves of a branch of a practical power grid

3.4 NQM計(jì)算時(shí)間分析

3種方法在不同采樣規(guī)模下對(duì)IEEE 30節(jié)點(diǎn)系統(tǒng)和實(shí)際網(wǎng)絡(luò)進(jìn)行計(jì)算所花費(fèi)的時(shí)間如表2所示。

圖6 3種方法的結(jié)果比較Fig.6 Comparison of simulative results among three methods

從表2中數(shù)據(jù)可得，3種方法的計(jì)算時(shí)間相當(dāng)，且均與采樣規(guī)模成正比，證明計(jì)算時(shí)間大部分消耗在潮流計(jì)算上，而樣本生成不會(huì)對(duì)計(jì)算時(shí)間產(chǎn)生大的影響。但由于LHS需要采樣和排序的過(guò)程，相對(duì)計(jì)算時(shí)間略長(zhǎng)。

表2 不同采樣規(guī)模下3種方法計(jì)算時(shí)間Table 2 Computational time of three methods under different sample sizes

4 結(jié)論

以拉丁超立方為代表的偽隨機(jī)采樣無(wú)法保證序列的低差異性，這限制了MCS精度的提高。本文針對(duì)該問(wèn)題，提出一種基于NQM的PLF方法。對(duì)IEEE 30節(jié)點(diǎn)系統(tǒng)和某實(shí)際電網(wǎng)的仿真證明了本文方法的有效性，并得到以下結(jié)論：

a.用LDSs序列可以克服MCS收斂性的瓶頸，因此NQM的收斂速度快于SRS和LHS；

b.使用奇異值分解對(duì)隨機(jī)變量進(jìn)行排序可在不增加計(jì)算量的情況下獲得含相關(guān)性的隨機(jī)序列；

c.SRS、LHS和NQM在期望值的計(jì)算結(jié)果上比較接近，而在標(biāo)準(zhǔn)差上NQM遠(yuǎn)遠(yuǎn)優(yōu)于另外2種方法，在相同采樣規(guī)模下可得到較精確的結(jié)果，因而生成的概率分布也較精確；

d.SRS、LHS和NQM 3種方法計(jì)算時(shí)間接近，而LHS由于采樣方式的原因消耗時(shí)間略長(zhǎng)于其余2種方法；

e.對(duì)我國(guó)某實(shí)際大區(qū)域電網(wǎng)的測(cè)試表明本文提出的算法對(duì)于大規(guī)模系統(tǒng)仍然適用。