楊少瑜，黃國棟，林星宇，樂彥婷，唐俊杰

（1.重慶理工大學電氣與電子工程學院，重慶市巴南區(qū)400054；2.中國電力科學研究院有限公司，北京市海淀區(qū)100192；3.輸配電裝備及系統(tǒng)安全與新技術國家重點實驗室(重慶大學)，重慶市沙坪壩區(qū)400044）

0 引言

截至2018年，全世界范圍內可再生能源發(fā)電的裝機容量占總裝機容量的33%以上，可再生能源全年發(fā)電量約占全球總發(fā)電量的26.2%，其中風力發(fā)電所占全球總發(fā)電的比例超過5%[1]。

隨著大量可再生能源的并網，風速、太陽輻射強度等不確定因素會對電力系統(tǒng)帶來大量的不確定源[2-3]。另外，隨著科技的發(fā)展，負荷類型也在不斷增加，負荷的不確定性也應該引起重視[4]。大量不確定源的引入，導致傳統(tǒng)的潮流分析方法不能夠有效地分析、評估電力系統(tǒng)當前的運行狀態(tài)和所面臨的運行風險。概率潮流（probabilistic power flow，PPF）考慮了輸入變量的概率特性，求取潮流輸出變量的概率信息，能夠有效地分析不確定源給電力系統(tǒng)運行帶來的影響[5]。

概率潮流分析中，對不確定源的概率建模有2個步驟：邊緣分布建模和相關性建模。邊緣分布建模是用特定分布類型和特定分布參數(shù)來描述單個隨機變量概率特性的過程，需要根據(jù)實測樣本，通過數(shù)理統(tǒng)計的方法來實現(xiàn)。常用的不確定源的物理對象和所用分布類型有：用韋布爾分布（Weibull distribution）描述風速的不確定性，正態(tài)分布（Normal distribution）描述負荷的波動規(guī)律等[4,6]。相關性建模是對2個隨機變量的不確定性之間關聯(lián)程度的刻畫，常用皮爾森相關系數(shù)來進行表征[7]。文獻[7]認為在實際應用中，隨機變量之間的相關性必須予以考慮，優(yōu)秀的概率方法必須具備處理隨機變量相關性的能力。

由于實際場景中不確定源種類多樣，所涉及的邊緣分布模型復雜，概率潮流分析時難以直接基于隨機變量概率模型產生計算樣本。對于這個問題，文獻[8]采用Nataf變換，建立了原始分布域與標準正態(tài)分布域之間的轉換關系，使樣本得以在標準正態(tài)分布域進行生成。該方法的關鍵在于求解標準正態(tài)分布域變量的相關系數(shù)。對于標準正態(tài)分布域相關系數(shù)的計算方法，文獻[9]給出了基于數(shù)值積分和二分法的求取方法，其計算過程較為繁瑣，耗時可以占到整個概率潮流分析總耗時的一半以上。文獻[10]給出了經驗公式來快速且精確地實現(xiàn)相關系數(shù)的計算，但是所涉及的分布類型極為有限。由于電力系統(tǒng)規(guī)模龐大，不確定源種類繁多，隨機變量服從的分布類型多樣，這些經驗公式已經不能滿足實際需求。

對于這樣的問題，本文采用拉格朗日插值法對原始域變量累積分布函數(shù)的反函數(shù)進行近似，使用簡化牛頓法對標準正態(tài)分布域相關系數(shù)進行求解，在保證計算精度的情況下，能夠提高標準正態(tài)分布域相關系數(shù)的計算效率。

1 基本模型

1.1 概率潮流模型

對于純交流電力系統(tǒng)，概率潮流的確定性模型的潮流方程可描述為：

式中：Uia和Uia分別為節(jié)點ia和節(jié)點ja的電壓幅值；Pia和Qia為節(jié)點ia上的注入有功功率和注入無功功率；Nb為系統(tǒng)節(jié)點個數(shù)；Giaja和Biaja分別為系統(tǒng)導納矩陣第ia行第ja列元素的實部和虛部；θiaja為節(jié)點ia與節(jié)點ja電壓的相位差。

考慮到新能源的接入，在實際計算中，確定性模型還需要嵌入從一次能源變量到對應電廠發(fā)電功率的轉換。本文考慮了風電場的并網，所使用的風速?風電轉換公式為[11]：

式中：vwind表示風速，m/s；PT表示單臺風機出力，MW；每臺風機消耗無功功率為恒定值0.0002 MV·A。

聯(lián)合公式（1）和公式（2），確定性潮流模型的輸入可以包含：可再生能源電廠的一次能源數(shù)據(jù)如風速，PV節(jié)點發(fā)電機的有功出力和電壓幅值，所有節(jié)點的有功、無功負荷；對于概率潮流分析，本文將電力系統(tǒng)中的隨機變量視為輸入，確定性模型的其余輸入變量均為常數(shù)，整個概率潮流的輸入變量用X表示；對于輸出變量，概率潮流的輸出可以包含PQ節(jié)點電壓幅值與相角，PV節(jié)點電壓相角，發(fā)電機無功出力和支路潮流等；所有輸出用Y表示。對于概率潮流分析，基于公式（1）和公式（2），確定性模型可簡化表示為：

1.2 隨機變量邊緣分布模型

對于電力系統(tǒng)中的隨機變量，風速用韋布爾分布（Weibull distribution）來建模，而負荷的波動用正態(tài)分布來描述（Normal distribution）。本文在只考慮風電廠接入和負荷隨機波動的情況下，所使用的邊緣分布模型如表1所示。另外，為了驗證本文所提出的相關系數(shù)計算法，也會用到表2所示概率分布。

表1 概率潮流中隨機變量的邊緣分布模型Table 1 Marginal distribution model of random variables in PPF calculation

表2 用以驗證本文方法的其他概率分布Table 2 Other probability distributions to verify the effectiveness of the proposed method

1.3 隨機變量線性相關性模型

本文以皮爾森相關系數(shù)作為描述隨機變量之間相關性的指標。對于第i個和第j個隨機變量Xi和Xj，已知其均值與標準差分別為μi，μj，σi和σj；另外，2個變量之積的均值為μij，Xi和Xj的相關系數(shù)ρXij可通過如下方式計算：

對于皮爾森相關系數(shù)，變量進行單調非線性變換（如Nataf變換）前后，變量之間的皮爾森相關系數(shù)數(shù)值會發(fā)生變化。在概率潮流分析中，Nataf變換基于標準正態(tài)分布變量生成用以概率潮流計算的樣本，而標準正態(tài)分布樣本的生成依賴于標準正態(tài)分布變量之間的相關系數(shù)，因此需要對標準正態(tài)分布域相關系數(shù)進行精確計算。

2 概率方法基本理論

2.1 Nataf變換

在概率潮流分析中，具有特定相關性的任意非正態(tài)分布樣本無法直接生成，需要借助于Nataf變換。Nataf變換采用等概率原理，建立原始分布域變量與標準正態(tài)分布域變量的轉換關系。已知m維原始域隨機 變 量X，即X=[X1,···,Xi,···,Xm]，原始域相關系數(shù)矩陣CX，其中第i行第j列元素表示變量Xi和Xj之間的相關系數(shù)，記為ρXi j。

設在標準正態(tài)分布域，有同等維度的標準正態(tài)隨機變量Z。第i個隨機變量Xi和Zi，滿足關系：

對于單峰的非正態(tài)分布，連續(xù)性隨機變量的累積分布函數(shù)的反函數(shù)是單調增函數(shù)，同時Φ是嚴格單調遞增的函數(shù)，那么公式（5）所表述的變換過程一定是非線性的單調變換的過程。因此，在標準正態(tài)分布域和原始域對應的相關系數(shù)數(shù)值不同，若要生成用以概率潮流分析的樣本，需要首先計算出標準正態(tài)分布域的相關系數(shù)。根據(jù)公式（4）對皮爾森相關系數(shù)的定義和概率論的基本知識，原始域隨機變量Xi和Xj之間的相關系數(shù)和標準正態(tài)分布變量Zi和Zj之間的相關系數(shù)ρZi j滿足如下關系：

通過線性變換可以建立獨立標準正態(tài)分布變量U和相關的標準正態(tài)分布變量Z之間的關系：

式中：L為下三角矩陣，可以通過對矩陣CZ進行Cholesky下三角分解得到，即CZ=LLT，同時可以得到：

基于以上理論，生成獨立的標準正態(tài)分布樣本后，使用式(5)—(8)可以生成原始域變量X的樣本，且這些樣本滿足相關系數(shù)矩陣CX。顯然，這個過程涉及到使用式(6)或式(8)求解標準正態(tài)分布域的相關系數(shù)矩陣CZ。而式(6)或式(8)都是含有二重積分的非線性隱式方程。傳統(tǒng)的基于數(shù)值積分和二分法的求解方法能夠較為精確地對ρZi j進行求解，但是耗時及其嚴重，一般不能滿足實時計算的需求。因此有必要對其求解方法加以改進，在保證較高精度的情況下快速完成相關系數(shù)的求解，從而加快概率潮流分析的效率。

2.2 基于拉格朗日插值法和簡化牛頓法的標準正態(tài)分布域相關系數(shù)計算法

拉格朗日插值法通過用多項式函數(shù)，連續(xù)穿過所給的離散點，來近似輸入與輸出之間的函數(shù)關系。對于n+1個給定樣本點(x0,y0),···,(x?,y?),···,(xn,yn)，其中所有的x互不相同?？赏ㄟ^式(9)來近似得到從x到y(tǒng)之間的關系y≈L(x)：

式中l(wèi)?(x)為插值基函數(shù)，其計算方法為

通過式(9)和式(10)可知，拉格朗日函數(shù)L(x)總能寫為

式中a0,a1,···,an均為多項式系數(shù)?，F(xiàn)使用n階拉格朗日插值法（階數(shù)為樣本點數(shù)?1），使用式(11)的形式分別近似式(6)中zi，zj和2個反函數(shù)的關系：

式中A0,A1, ···,An和B0,B1,···,Bn均為多項式系數(shù)。聯(lián)合式(8)和式(12)，式(8)中分子的積分項為

其中C表示組合數(shù)。式(13)中的積分滿足：

其中“!!”表示雙階乘。于是式(13)變?yōu)?/p>

式中：p+k和q?k為偶數(shù)。將式(15)帶入式(8)，式(8)可化為不含積分，僅含求和的表達式，簡化表示為：

式(6)表明對ρZi j的求解是解一元非線性方程的問題。根據(jù)實驗，設定收斂精度后，以ρXi j為迭代初始值，G(ρZi j)的導數(shù)值恒定為1，使用簡化牛頓法，通過若干次迭代便能夠對式（16）精確求解。

在實際應用中，對任意2個變量之間的相關系數(shù)均實施上述過程，最終可得到標準正態(tài)分布域的相關系數(shù)矩陣CZ。另外，相關系數(shù)矩陣為實對稱矩陣，且主對角線元素為1。因此在含有m個隨機變量的情況下，需要計算的相關系數(shù)個數(shù)為0.5m(m?1)，即只需計算矩陣上三角（或下三角）除主對角線以外的元素。最后，此方法是基于拉格朗日插值法和簡化牛頓法的求解法（Lagrange interpolation and simplified Newton method，LISNM），在后文中用LISNM表示。本文所提到的改進Nataf變換，即使用LISNM計算標準正態(tài)分布域相關系數(shù)的Nataf變換方法。

2.3 基于改進Nataf變換的PPF計算方法

在Nataf變換的過程中，LISNM可以對標準正態(tài)分布域的相關系數(shù)進行精確而快速的計算，從而提高整個概率潮流計算的效率。需要注意的是，Nataf變換用于產生具有特定相關性的原始域樣本，而不受到特定概率潮流計算方法的限制，如文獻[8,12-13]在概率潮流中對Nataf變換的使用。

在實際應用中，基于改進Nataf變換的概率潮流計算方法可以參照如下步驟：

步驟1：確定電力系統(tǒng)中的不確定源，即隨機變量X，并通過數(shù)理統(tǒng)計的方法得到其邊緣分布和相關系數(shù)矩陣CX。

步驟2：反復使用LISNM得到標準正態(tài)分布域的相關系數(shù)，最終得到標準正態(tài)分布域相關系數(shù)矩陣CZ。

步驟3：基于特定的概率潮流方法，生成獨立的標準正態(tài)分布樣本，進一步使用公式（7）得到具有相關性的標準正態(tài)分布樣本。

步驟4：使用公式（5），將所有標準正態(tài)分布樣本轉換回原始分布域。

步驟5：將所有原始分布域的樣本，代入確定性模型（3）進行計算，得到相同數(shù)量的輸出樣本。

步驟6：對于輸出樣本，基于所用的概率潮流方法的理論，得到輸出變量的概率信息。

3 算例分析

3.1 LISNM的性能測試

使用表1和表2中的分布類型對LISNM的性能進行測試。設這4類分布的具體信息為：韋布爾分布，形狀參數(shù)和尺度參數(shù)分別為3和10；正態(tài)分布，均值和標準差分別為10和1；均勻分布，下界和上界分別為5和15；對數(shù)正態(tài)分布，對數(shù)均值和對數(shù)標準差分別為3和0.04。為了充分表明所提方法的適用性，設原始域相關系數(shù)數(shù)值R，R=[?0.95,?0.9,?0.8,?0.7,?0.6,?0.5,?0.4,?0.3,?0.2,?0.1,0,0.1,0.2,0.3,0.4,0.5,0.6,0.7,0.8,0.9,0.95]。對于R中所有數(shù)值，任意2個分布之間、每個分布自身之間都將對其進行計算。

使用高精度的基于辛普森數(shù)值積分法和二分法的相關系數(shù)求解法作為精度參考（Simpson’s method and bisection method, SMBM），來測試LISNM的精度。其中SMBM的收斂精度設置為10?4，表示在迭代結束時，SMBM的結果與相關系數(shù)真實值的絕對誤差小于10?4。

定義相關系數(shù)的最大絕對誤差εmax和平均絕對誤差εmean，其計算方法如公式（17）所示：

式中：ρref和ρtested分別表示通過SMBM計算得到的相關系數(shù)和通過LISNM計算得到的相關系數(shù)；Nρ表示相關系數(shù)的個數(shù)。

對于所設的4類變量之間以及每類變量自身之間，分別使用SMBM和LISNM對R中的所有數(shù)值進行計算，得到標準正態(tài)分布域的相關系數(shù)精度和計算總時間的結果如表3所示。

表3 相關系數(shù)計算結果Table 3 Computing results of correlation coefficients

在表3中，LISNM的拉格朗日插值法使用插值點個數(shù)為17，簡化牛頓法迭代初值與原始域相關系數(shù)數(shù)值相同，雅可比矩陣為固定值1。根據(jù)表3的結果，LISNM的計算時間比SMBM小了3個數(shù)量級，計算時間大大減少。而計算精度上，根據(jù)SMBM的精度設定，SMBM結果小數(shù)點后4位是完全精確的，而LISNM計算得到的結果與SMBM相比，最大誤差和平均誤差均在10?4的數(shù)量級。也就表明本文所提方法的精度至少為10?3，從數(shù)值上來看，其結果非常精確。

3.2 基于改進Nataf變換的PPF計算

從3.1節(jié)的結果可以看到，LISNM的計算速度快，精度高。設定以下算例，測試LISNM的概率潮流計算精度的影響。概率潮流方法分別為蒙特卡洛仿真法（Monte Carlo simulation method,MCSM）和Hong氏e method,HPEM）[14]。

基于MCSM的方法在樣本數(shù)足夠大時，其計算結果也是依概率收斂于精確的結果。同時由于計算機只能夠偽隨機生成計算樣本，即使樣本數(shù)量非常大，樣本的相關系數(shù)矩陣與所給矩陣本身便存在一定誤差。通過實驗測試，樣本足夠多時，這個誤差的最大值小于10?3。根據(jù)3.1節(jié)結果可以知道，LISNM計算的相關系數(shù)矩陣誤差和MCSM隨機樣本相關系數(shù)的天然誤差相當。因此，理論上基于MCSM所測得的結果誤差應該更大。

但是對于以HPEM為代表的近似法，在獨立標準正態(tài)分布域采得的樣本是固定的。根據(jù)公式（7），在實施過程中，每個樣本點都會直接受到標準正態(tài)分布域相關系數(shù)矩陣的影響。并且在算例既定的情況下，其原始域樣本數(shù)值僅受到相關系數(shù)矩陣的影響。因此，HPEM的計算結果對相關系數(shù)矩陣的數(shù)值非常敏感。相比于MCSM而言，HPEM的潮流結果誤差更能體現(xiàn)Nataf變換中相關系數(shù)計算精度對概率潮流分析結果精度的影響。

基于標準的IEEE 118節(jié)點系統(tǒng)，接入如文獻[11]完全相同的風電裝置(共14處)。對于所有負荷，其有功負荷均服從正態(tài)分布，均值為算例中對應原始值，標準差為均值的5%。另外，在計算中保證負荷的功率因素與原始算例對應負荷的功率因素相同，從而決定無功負荷。其中有功負荷之間的相關系數(shù)為0.6，負荷與風速之間相關系數(shù)為0.05。另外，本文用相對誤差來描述均值與標準差計算結果的精確度，數(shù)值越小代表該結果越精確，其具體計算方法如公式（18）所示：

式中：εr表示相對誤差；rref和rtested分別代表通過SMBM和LISNM進行相關系數(shù)轉換后的概率潮流計算結果。根據(jù)設定，在使用同一種概率計算方法時，認為基于SMBM的概率分析最終結果是標準結果，考察基于LISNM的概率分析結果與參考結果的偏差（即誤差）。

選取網絡有功損耗、平衡節(jié)點發(fā)電機無功出力作為輸出。分別使用SMBM和LISNM進行相關系數(shù)計算，使用基于100000個簡單隨機樣本的MCSM計算結果的概率密度如圖1所示。對于概率密度精確度的衡量，使用文獻[15]中定義的指標頻率直方圖相似度(frequency histogram similarity index,FHSI)來描述2個概率分布的相似程度。當其中一個概率密度作為標準結果時，F(xiàn)HSI的數(shù)值可以作為另一概率密度的精度指標。對于FHSI，其數(shù)值范圍分布在區(qū)間[0,1]內，數(shù)值越大表明所指結果越精確。根據(jù)文獻[15]所給判據(jù)，當FHSI數(shù)值不小于0.9時，認為所指的概率分布是精確的。

圖1 基于SMBM和LISNM相關系數(shù)計算法的MCSM概率潮流結果Fig.1 PPF results calculated by MCSM based on correlation coefficients obtained by SMBM and LISNM respectively

另外，分別使用MCSM和HPEM計算得到的輸出的均值與標準差誤差如表4所示。在圖1中，2種方法得到的輸出概率密度的相似度超過0.98。根據(jù)表4中的結果，當使用LISNM計算相關系數(shù)，使用MCSM計算概率潮流，最終結果的均值與標準差的誤差均不到0.2%。這表明使用LISNM進行相關系數(shù)計算，且使用MCSM實現(xiàn)概率潮流分析的結果精度非常高，無論是輸出變量的概率分布還是其均值和標準差都能夠滿足實際應用需求。

表4 基于改進Nataf變換的MCSM和HPEM概率潮流計算誤差Table 4 PPF calculation errors of improved Nataf transformation based MCSM and HPEMmethods%

根據(jù)前面的分析，HPEM結果的均值與標準差對標準正態(tài)分布域的相關系數(shù)矩陣元素的數(shù)值更加敏感。但是表4的結果表明，LISNM的誤差對HPEM的概率潮流計算結果的影響非常小。這一方面說明了LISNM對相關系數(shù)計算的精度完全滿足概率潮流計算的要求。另一方面，這表明LISNM的誤差給概率潮流結果帶來的偏差，遠遠小于基于MCSM的概率潮流計算結果本身的不穩(wěn)定性。

3.3 改進Nataf變換的必要性和PPF耗時分析

采用SMBM計算標準正態(tài)分布域相關系數(shù)，且以0.1為收斂精度，記作in-SMBM(inaccurate SMBM)。分別使用100000組簡單隨機樣本的MCSM和HPEM法進行概率潮流計算，其計算結果的誤差見表5。

表5 相關系數(shù)計算不精確時的概率潮流結果誤差Table 5 Calculation error of PPF under inaccurately calculated correlation coefficients%

比較表4與表5中的結果，在相關系數(shù)計算不精確的情況下，概率潮流計算輸出的誤差均有明顯增大。其中均值誤差雖有增加，但是在數(shù)值上可以認為是精確的。但是標準差的誤差均超過了9%，相比于表4中的不到0.2%，其誤差增大明顯，精度已經不能滿足實際應用的需求。因此，不精確的標準正態(tài)分布域的相關系數(shù)將增大概率潮流計算結果的誤差，其中輸出變量的標準差對相關系數(shù)的精度更為敏感。這表明了在概率建模中精確地對標準高斯域相關系數(shù)進行計算的必要性。

概率潮流計算所消耗的時間，主要由相關系數(shù)計算耗時和概率潮流計算耗時組成。表6給出分別使用LISNM和SMBM法計算相關系數(shù)以及分別使用HPEM和MCSM計算概率潮流的耗時情況，并給出了該組合方法完成整個概率潮流分析的總耗時。

根據(jù)表6中不同方法下相關系數(shù)計算耗時的占比可以發(fā)現(xiàn)，相關系數(shù)計算的耗時在整個概率潮流分析中是不可忽略的。對于特定的概率方法，計算相關系數(shù)的耗時甚至占到整個概率分析耗時的80%以上。

表6 相關系數(shù)計算耗時與概率潮流計算耗時的比較Table 6 Comparison of computing time for correlation coefficient calculation and that for PPF calculation

在表6所示的各種PPF方法組合中，相比于SMBM，LISNM對于概率分析效率的提升非?？捎^。

4 結論

1）本文提出的LISNM能夠快速地計算標準正態(tài)分布域的相關系數(shù)，加快了Nataf變換的實施速度，極大地提高了概率潮流分析的效率。同時，LISNM對相關系數(shù)的計算精度高，在提升概率潮流分析效率的同時，保證其結果的精度。

2）基于LISNM的相關系數(shù)計算的誤差，經概率潮流傳遞后，對概率潮流計算結果的影響遠遠小于蒙特卡洛法自身結果的不穩(wěn)定性。而通過Hong氏點估計法的傳遞，LISNM的誤差對概率潮流結果的影響非常小，這也表明了LISNM的精確性。

3）在電力系統(tǒng)概率潮流分析中，相關系數(shù)計算的耗時占比大。提高相關系數(shù)計算的效率，可以極大加快整個概率潮流分析的過程。