方斯頓，程浩忠，徐國棟，曾平良，姚良忠

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

0 引言

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

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

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

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

為提高MCS的計算效率并處理相關(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)過若干次確定型潮流計算得到輸出變量的概率分布特性。對IEEE 30節(jié)點系統(tǒng)和我國某實際大區(qū)域電網(wǎng)的仿真驗證了本文方法的有效性。

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

Nataf變換是在已知輸入變量邊緣分布的情況下重構(gòu)聯(lián)合分布的數(shù)學(xué)方法[18]?，F(xiàn)將其過程簡述如下：對任意 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ù)。隨機向量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)差；φZiZj（Zi，Zj，ρZij）為相關(guān)系數(shù)為ρZij的二維標(biāo)準(zhǔn)正態(tài)分布概率密度函數(shù)。本文采用文獻[19]中方法對式（3）進行求解以求得[ρZij]n×n。 而由于[ρZij]n×n可能為非正定或非滿秩陣，此時其Cholesky分解不存在，無法采用文獻[16]的方法對樣本進行排序。但一般相關(guān)系數(shù)矩陣都為對稱陣，存在奇異值分解，本文通過證明定理1說明奇異值分解排序方法的過程。

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

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

證明如下。

對?z*i，其期望值滿足式（6）：

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

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

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

2.1 MCS及其誤差

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

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

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

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

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

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

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

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

則MCS誤差為：

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

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

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

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

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

對比式（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計算50000次的結(jié)果作為標(biāo)準(zhǔn)，定義估計誤差為：

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

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

2.4 NQM計算流程

綜合前文，NQM的計算流程如圖1所示。簡述如下：

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

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

c.利用文獻[21]中方法生成0-1分布的獨立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進行M次確定型潮流計算，并統(tǒng)計得到各輸出變量的數(shù)字特征及概率分布。

圖1 計算流程圖Fig.1 Flowchart of calculation

3 算例分析

3.1 算例介紹

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

風(fēng)電機出力的互相關(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)系滿足：

其中，PWR為風(fē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種方法對2個測試系統(tǒng)進行計算。按式（19）計算輸出變量誤差，分別表示為：節(jié)點電壓幅值期望誤差，標(biāo)準(zhǔn)差誤差；節(jié)點電壓相角期望誤差，標(biāo)準(zhǔn)差誤差；支路有功期望誤差，標(biāo)準(zhǔn)差誤差；支路無功潮流期望誤差，標(biāo)準(zhǔn)差誤差。 IEEE 30 節(jié)點系統(tǒng)輸出結(jié)果見圖2、3。實際網(wǎng)絡(luò)計算結(jié)果見表1。

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

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

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

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

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

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

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

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

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

3.3 NQM收斂趨勢分析

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

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

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

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

3.4 NQM計算時間分析

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

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

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

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

4 結(jié)論

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

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

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

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

d.SRS、LHS和NQM 3種方法計算時間接近，而LHS由于采樣方式的原因消耗時間略長于其余2種方法；

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