陳宇辛，江岳文，2

（1. 福州大學(xué) 電氣工程與自動(dòng)化學(xué)院，福建 福州 350108；2. 福建省電器智能化工程研究中心，福建 福州 350108）

0 引言

為分析可再生能源注入、負(fù)荷變化等不確定因素對(duì)電力系統(tǒng)運(yùn)行狀態(tài)的影響，概率潮流、模糊潮流、區(qū)間潮流等不確定性潮流分析工具被提出［1-2］。其中，區(qū)間潮流建模簡單，只需知道不確定性變量的邊界就可以得到潮流結(jié)果的上、下界，因此受到學(xué)者的廣泛關(guān)注。

文獻(xiàn)［3］將區(qū)間分析方法與潮流計(jì)算原理相結(jié)合，針對(duì)區(qū)間牛頓算子矩陣求逆困難的缺陷，使用Krawcyzk-Moore 迭代法與區(qū)間對(duì)分法收縮初始區(qū)間。由于忽略了變量自身的相關(guān)性，區(qū)間運(yùn)算帶來的區(qū)間擴(kuò)張會(huì)使計(jì)算結(jié)果過于保守，而仿射算術(shù)可以在一定程度上減小保守性［4］。文獻(xiàn)［5］在文獻(xiàn)［3］的基礎(chǔ)上，以仿射算術(shù)替代潮流計(jì)算過程中的區(qū)間加減法運(yùn)算，有效縮小了區(qū)間解的范圍。迭代算法的收斂性依賴于初始區(qū)間的選擇，存在收斂性較差的問題。文獻(xiàn)［6］基于泰勒展開法將區(qū)間潮流問題分解為3 個(gè)確定性問題進(jìn)行求解，避免了迭代運(yùn)算，但計(jì)算二階分量時(shí)仍涉及區(qū)間乘除法運(yùn)算，存在一定的保守性。

區(qū)間潮流問題可以轉(zhuǎn)化為優(yōu)化問題，待求的是潮流狀態(tài)變量的最大值和最小值［7］。文獻(xiàn)［8］采用優(yōu)化算法建立線性規(guī)劃、非線性規(guī)劃和二次規(guī)劃模型，分別求解支路兩端相角差、節(jié)點(diǎn)電壓、支路功率，避免了區(qū)間迭代運(yùn)算。文獻(xiàn)［9］基于仿射算數(shù)重構(gòu)潮流方程，通過優(yōu)化算法壓縮仿射算子中的噪聲元得到潮流區(qū)間。由于在極坐標(biāo)系下仿射算法需要對(duì)三角函數(shù)近似估計(jì)，在直角坐標(biāo)系下電壓方程會(huì)放大電壓幅值的解，仿射優(yōu)化算法不能完全消除保守性。文獻(xiàn)［10］結(jié)合區(qū)間匹配和區(qū)間極值的選取理論，將區(qū)間潮流問題分解為2 個(gè)確定的非線性優(yōu)化問題，但容易忽略少數(shù)極端情況。文獻(xiàn)［11］基于極值理論，以每個(gè)待求狀態(tài)量的上、下界為目標(biāo)函數(shù)，采用非線性優(yōu)化模型進(jìn)行求解，該方法求得的結(jié)果不存在區(qū)間擴(kuò)張問題，是一種精度較高的區(qū)間潮流算法，但對(duì)每個(gè)狀態(tài)變量都分別求解上、下界需要進(jìn)行大量非線性優(yōu)化計(jì)算，計(jì)算效率較低。

上述區(qū)間潮流算法均未考慮變量之間的相關(guān)性，而在實(shí)際電力系統(tǒng)中，各變量之間的相關(guān)性是廣泛存在的［12］，例如同一地區(qū)的負(fù)荷受環(huán)境或社會(huì)因素影響產(chǎn)生相關(guān)性，地理位置相近的風(fēng)電場出力具有較強(qiáng)的相關(guān)性等［13］。若不考慮各變量之間的相關(guān)性，則可能導(dǎo)致區(qū)間潮流計(jì)算結(jié)果不合理。文獻(xiàn)［14］采用橢球模型處理變量間的相關(guān)性，但用橢球模型處理獨(dú)立變量時(shí)會(huì)放大不確定區(qū)域，導(dǎo)致保守性問題。文獻(xiàn)［15］用相關(guān)角表示變量間的相關(guān)性，通過坐標(biāo)系變換實(shí)現(xiàn)去相關(guān)化。文獻(xiàn)［16］以平行四邊形模型表示相關(guān)性，為相關(guān)性區(qū)域提供明確的表達(dá)式。

基于上述分析，為了獲得更精確的潮流區(qū)間邊界，本文首先構(gòu)造平行四邊形模型表示區(qū)間變量間的相關(guān)性；然后，為了提高計(jì)算效率，采用節(jié)點(diǎn)電壓對(duì)節(jié)點(diǎn)功率的靈敏度分析方法，將變化趨勢一致的變量歸類，大幅減少優(yōu)化次數(shù)；最后，建立考慮相關(guān)性和靈敏度分類的改進(jìn)場景優(yōu)化模型求解區(qū)間潮流。算例分析結(jié)果驗(yàn)證了本文方法的有效性。

1 考慮區(qū)間變量相關(guān)性的平行四邊形模型

不同變量的波動(dòng)可能是由相同因素引起的，因此不同的區(qū)間變量間可能存在相關(guān)性。本文采用平行四邊形模型表示區(qū)間變量之間的相關(guān)性［16］。如圖1所示，以2個(gè)具有相關(guān)性的區(qū)間變量為例說明構(gòu)建平行四邊形模型的方法。

圖1 平行四邊形模型Fig.1 Parallelogram model

1）該平行四邊形包含所有的樣本點(diǎn)，且面積最小；

2）該平行四邊形的4 個(gè)頂點(diǎn)均落在矩形的對(duì)角線上；

3）該平行四邊形至少有1 條對(duì)角線上的頂點(diǎn)與矩形的頂點(diǎn)重合。

當(dāng)Xi、Xj間存在相關(guān)性時(shí)，兩者構(gòu)成的不確定域可用由l1、l2、l3、l4這4 條直線圍成的平行四邊形A′BC′D表示，該區(qū)域可用如下形式表示：

式中：ΔXi=Xi-XCi，ΔXj=Xj-XCj，推導(dǎo)過程見附錄A。

當(dāng)存在m個(gè)具有相關(guān)性的區(qū)間變量X1、X2、…、Xm時(shí)，考慮區(qū)間變量兩兩間的相關(guān)性，相關(guān)系數(shù)矩陣為：

則具有相關(guān)性的m個(gè)區(qū)間變量存在的多維平行四邊形域可表示為：

2 考慮靈敏度分類和相關(guān)性的場景優(yōu)化區(qū)間潮流算法

2.1 基于場景優(yōu)化法的區(qū)間潮流計(jì)算

區(qū)間潮流是通過已知的不確定量的上、下界獲取潮流變量上、下界的潮流計(jì)算方法，所得結(jié)果可以囊括電力系統(tǒng)中所有可能的運(yùn)行工況。區(qū)間潮流方程可表示為：

式中：min 子模型對(duì)應(yīng)變量的下界；max 子模型對(duì)應(yīng)變量的上界。

采用式（13）所描述的場景優(yōu)化法分別求解每個(gè)潮流狀態(tài)變量的上、下界，對(duì)于一個(gè)含有nPQ個(gè)PQ節(jié)點(diǎn)的n節(jié)點(diǎn)電力系統(tǒng)而言，求解節(jié)點(diǎn)電壓幅值、相角區(qū)間共需要經(jīng)過2(nPQ+n-1)次非線性優(yōu)化運(yùn)算。而有些狀態(tài)變量在取得最大值或最小值時(shí)的場景是相同的，若分開計(jì)算則會(huì)造成重復(fù)工作，浪費(fèi)過多的計(jì)算時(shí)間成本和計(jì)算資源，因此可將這些在相同場景取得最大值或最小值的變量歸為一類，其最大值或最小值可以在同一次優(yōu)化計(jì)算中得到，為此，本文提出利用靈敏度分類法實(shí)現(xiàn)變量歸類。

2.2 靈敏度分類法

在系統(tǒng)網(wǎng)架結(jié)構(gòu)和元件參數(shù)均確定的情況下，系統(tǒng)狀態(tài)變量的變化是由節(jié)點(diǎn)注入功率的變化引起的，每個(gè)節(jié)點(diǎn)注入功率的變化引起不同狀態(tài)變量變化的程度不同，但變化方向可能相同。場景優(yōu)化法是在節(jié)點(diǎn)注入功率波動(dòng)的區(qū)間范圍內(nèi)尋找使得某一狀態(tài)變量達(dá)到最大值或最小值的場景。若存在一類狀態(tài)變量，每個(gè)節(jié)點(diǎn)注入功率的變化引起這類狀態(tài)變量的變化方向是相同的，即這類狀態(tài)變量隨著各節(jié)點(diǎn)注入功率的變化同時(shí)增大或減小，則這類狀態(tài)變量能在同一場景中取得最大值或最小值，因此，這類狀態(tài)變量的最大值或最小值可在同一次優(yōu)化計(jì)算中獲得。

靈敏度是利用系統(tǒng)中物理量的微分關(guān)系，來獲得因變量對(duì)自變量敏感程度的方法［17］。本文采用節(jié)點(diǎn)電壓對(duì)節(jié)點(diǎn)注入功率的靈敏度來歸類具有相同變化趨勢的狀態(tài)變量。已知用極坐標(biāo)表示的雅可比矩陣的潮流求解方程為：

式中：ΔP、ΔQ分別為節(jié)點(diǎn)有功、無功功率波動(dòng)量的列向量；J11、J12、J21、J22為雅可比矩陣的分塊矩陣；Δθ、ΔU分別為節(jié)點(diǎn)電壓相角、幅值波動(dòng)量的列向量。

對(duì)雅可比矩陣求逆，可得：

記雅可比矩陣的逆矩陣為R，R即為節(jié)點(diǎn)電壓對(duì)節(jié)點(diǎn)注入功率的靈敏度矩陣。設(shè)節(jié)點(diǎn)j注入功率的功率因數(shù)角為φj，可得節(jié)點(diǎn)j有功功率變化對(duì)節(jié)點(diǎn)i電壓影響的方程為：

2.3 考慮靈敏度分類和相關(guān)性的區(qū)間潮流模型

以分類后的每類潮流狀態(tài)變量之和的最大化和最小化為目標(biāo)函數(shù)，考慮輸入變量的相關(guān)性，結(jié)合靈敏度分類法的改進(jìn)場景優(yōu)化法區(qū)間潮流模型為：

式中：Sk為分類后的第k類狀態(tài)變量集合；PG，i、QG，i分別為節(jié)點(diǎn)i處發(fā)電機(jī)注入的有功功率和無功功率；PW，i、QW，i分別為節(jié)點(diǎn)i處可再生能源注入的有功功率和無功功率；PL，i、QL，i分別為節(jié)點(diǎn)i處負(fù)荷的有功功率和無功功率；Ui、Uj分別為節(jié)點(diǎn)i、j處的電壓幅值；Gij、Bij分別為節(jié)點(diǎn)導(dǎo)納矩陣第i行第j列元素的實(shí)部和虛部；θij=θi-θj，θi、θj分別為節(jié)點(diǎn)i、j的電壓相角；SN為系統(tǒng)所有節(jié)點(diǎn)集合；PˉW，i、-PW，i分別為PW，i的上、下界；SW為可再生能源節(jié)點(diǎn)集合；PˉL，i、-PL，i分別為PL，i的上、下界；φi為節(jié)點(diǎn)i處負(fù)荷的功率因數(shù)角；SL為負(fù)荷節(jié)點(diǎn)集合。

3 算例分析

在改進(jìn)的IEEE 39 節(jié)點(diǎn)系統(tǒng)中進(jìn)行區(qū)間潮流計(jì)算，在3.6 GHz 主頻和32 GB RAM 的條件下采用MATLAB 2016b 進(jìn)行編程。改進(jìn)的IEEE 39 節(jié)點(diǎn)系統(tǒng)拓?fù)浣Y(jié)構(gòu)見附錄B，系統(tǒng)中節(jié)點(diǎn)1—29 為PQ 節(jié)點(diǎn)，節(jié)點(diǎn)31 為平衡節(jié)點(diǎn)，其余節(jié)點(diǎn)為PV 節(jié)點(diǎn)。采用標(biāo)幺值進(jìn)行計(jì)算，系統(tǒng)基準(zhǔn)容量為100 MV·A。設(shè)負(fù)荷為恒功率因數(shù)負(fù)荷，風(fēng)電場1 接入節(jié)點(diǎn)4，風(fēng)電場2 接入節(jié)點(diǎn)18，兩風(fēng)電場的有功出力區(qū)間均為［100，300］MW，節(jié)點(diǎn)3 負(fù)荷波動(dòng)范圍為±30%，即［225.4，418.6］MW，其余負(fù)荷波動(dòng)范圍為±10%。設(shè)兩風(fēng)電場出力及節(jié)點(diǎn)3 負(fù)荷三者之間存在相關(guān)性，相關(guān)系數(shù)均為0.6。

3.1 算法有效性分析

采用本文方法與考慮相關(guān)性的蒙特卡羅法進(jìn)行對(duì)比。蒙特卡羅法中，變量在區(qū)間內(nèi)服從均勻分布，抽樣次數(shù)為50000次。以無波動(dòng)時(shí)的潮流結(jié)果計(jì)算靈敏度并進(jìn)行分類，節(jié)點(diǎn)電壓計(jì)算結(jié)果對(duì)比如圖2所示（圖中電壓幅值為標(biāo)幺值，后同）。

圖2 節(jié)點(diǎn)電壓區(qū)間分布對(duì)比Fig.2 Comparison of interval distribution of node voltage

由圖2 可看出，本文方法得到的區(qū)間潮流計(jì)算結(jié)果可以將蒙特卡羅法的計(jì)算結(jié)果包含在內(nèi)。這是由于在有限抽樣次數(shù)下，蒙特卡羅法無法模擬系統(tǒng)所有的運(yùn)行情況，這將導(dǎo)致計(jì)算結(jié)果的區(qū)間偏小，而本文方法在所有可能出現(xiàn)的場景中尋找狀態(tài)變量的最大值與最小值，所得區(qū)間結(jié)果可能在實(shí)際運(yùn)行中出現(xiàn)，不存在保守性問題。

本文方法在計(jì)算過程中將電壓幅值分為7 類，將電壓相角分為1 類，電壓幅值分類情況如表1所示。

表1 PQ節(jié)點(diǎn)電壓幅值分類結(jié)果Table 1 Classification results of voltage magnitude for PQ nodes

電壓相角僅分為1 類，說明所有待求的電壓相角擁有相同的變化趨勢。當(dāng)電壓相角達(dá)到上界時(shí)，節(jié)點(diǎn)3 負(fù)荷以及兩風(fēng)電場有功出力均取區(qū)間上界，其余節(jié)點(diǎn)負(fù)荷均取區(qū)間下界；當(dāng)電壓相角達(dá)到下界時(shí)，情況則相反。電壓相角達(dá)到最大值或最小值的情況一般是在節(jié)點(diǎn)注入有功功率達(dá)到下界或上界時(shí)出現(xiàn)，即負(fù)荷最小、風(fēng)電出力最大或負(fù)荷最大、風(fēng)電出力最小。由于風(fēng)電與負(fù)荷之間正相關(guān)性的約束，節(jié)點(diǎn)3 負(fù)荷與風(fēng)電同時(shí)取最大值或最小值時(shí)節(jié)點(diǎn)相角達(dá)到邊界值。

本文采用靈敏度分類法將有相同變化趨勢的潮流變量進(jìn)行分類，計(jì)算靈敏度需要將潮流方程線性化，且計(jì)算時(shí)采用的是無波動(dòng)情況下的潮流計(jì)算結(jié)果。在節(jié)點(diǎn)注入功率不斷變化的過程中某些狀態(tài)變量對(duì)注入功率靈敏度的正、負(fù)性可能會(huì)改變，存在一定誤差。為分析采用靈敏度分類法所帶來的誤差，將本文方法與文獻(xiàn)［11］中的未分類場景優(yōu)化法進(jìn)行比較。以文獻(xiàn)［11］方法的計(jì)算結(jié)果為基準(zhǔn)，改進(jìn)的IEEE 39節(jié)點(diǎn)系統(tǒng)中各節(jié)點(diǎn)電壓的誤差見附錄C。

由附錄C 可知，采用本文方法求解得到的誤差很小，電壓幅值平均誤差為0.0105%，電壓相角平均誤差為0.001 3%。本文方法對(duì)電壓幅值、相角單次優(yōu)化平均計(jì)算時(shí)間分別為4.442 4、1.142 1 s，文獻(xiàn)［11］方法相應(yīng)用時(shí)分別為4.6376、1.5325 s。采用文獻(xiàn)［11］方法求解節(jié)點(diǎn)電壓幅值和相角區(qū)間分別需要經(jīng)過58次和76次非線性優(yōu)化運(yùn)算，而本文方法對(duì)電壓幅值和相角的優(yōu)化次數(shù)分別為14 次和2 次?？梢姡疚姆椒ㄔ诒ＷC精度的情況下，削減了優(yōu)化次數(shù)，縮短了總體計(jì)算時(shí)間，提升了計(jì)算效率。

為進(jìn)一步驗(yàn)證本文方法的有效性，分別在IEEE 14、IEEE 30、IEEE 57 和IEEE 118 節(jié)點(diǎn)系統(tǒng)中進(jìn)行對(duì)比測試，計(jì)算結(jié)果如表2 所示（表中平均誤差均為標(biāo)幺值）?？梢钥闯觯诓煌瑴y試系統(tǒng)中本文方法均能保證較高精度，本文方法單次優(yōu)化平均計(jì)算時(shí)間與文獻(xiàn)［11］方法相當(dāng)，但減少了優(yōu)化次數(shù)，總時(shí)間大幅縮短。4個(gè)系統(tǒng)中電壓相角均僅分為1類，電壓幅值分類情況見附錄D。總體上，計(jì)算電壓相角的時(shí)間短于計(jì)算電壓幅值的時(shí)間。當(dāng)系統(tǒng)規(guī)模增大至IEEE 118 節(jié)點(diǎn)時(shí)，非線性優(yōu)化模型中的變量及約束條件增多導(dǎo)致計(jì)算時(shí)間大幅增長。

表2 不同測試系統(tǒng)中的計(jì)算結(jié)果對(duì)比Table 2 Comparison of calculation results among different test systems

3.2 相關(guān)性對(duì)潮流的影響分析

為分析輸入變量之間的相關(guān)性對(duì)潮流結(jié)果的影響，將3.1 節(jié)改進(jìn)IEEE 39 節(jié)點(diǎn)系統(tǒng)中的算例與不考慮相關(guān)性時(shí)的計(jì)算結(jié)果進(jìn)行對(duì)比?？紤]相關(guān)性前、后電壓幅值和電壓相角的分布區(qū)間如圖3所示。

由圖3 可以看出，與不考慮相關(guān)性的區(qū)間結(jié)果相比，考慮相關(guān)性的電壓幅值、相角區(qū)間范圍更窄。這是由于考慮相關(guān)性后區(qū)間變量的可行分布區(qū)域變小，使得計(jì)算結(jié)果區(qū)間縮小?？紤]相關(guān)性前、后區(qū)間變量分布區(qū)域?qū)Ρ热鐖D4 所示，設(shè)兩風(fēng)電場出力與節(jié)點(diǎn)3負(fù)荷三者間的相關(guān)系數(shù)均為ρ，平行六面體所包圍的區(qū)域?yàn)閷?duì)應(yīng)相關(guān)系數(shù)下風(fēng)電和負(fù)荷所有可能的取值場景。

圖3 相關(guān)性對(duì)節(jié)點(diǎn)電壓的影響Fig.3 Influence of correlation on bus voltage

圖4 考慮相關(guān)性前、后區(qū)間變量分布區(qū)域?qū)Ρ菷ig.4 Comparison of distribution areas of interval variables between before and after considering correlation

進(jìn)一步分析不同相關(guān)性對(duì)潮流的影響，以節(jié)點(diǎn)3 電壓相角區(qū)間為觀察對(duì)象，設(shè)兩風(fēng)電場出力之間的相關(guān)系數(shù)為0.6，負(fù)荷與兩風(fēng)電場出力之間的相關(guān)系數(shù)均為ρwl，當(dāng)ρwl從-0.6 變化至0.6 時(shí)，節(jié)點(diǎn)3 電壓相角區(qū)間變化情況如圖5所示。

由圖5 可以看出：ρwl從-0.6 增至0 時(shí)，相角區(qū)間并未變化；ρwl從0 增至0.6 時(shí)，隨著正相關(guān)性的增強(qiáng)，區(qū)間逐漸減小。不同正相關(guān)系數(shù)下節(jié)點(diǎn)3 電壓相角達(dá)到最大值和最小值的場景如圖6所示。

圖5 不同相關(guān)系數(shù)下節(jié)點(diǎn)3電壓相角區(qū)間Fig.5 Interval of voltage angle for Bus 3 under different correlation coefficients

如圖6（a）所示，在ρwl=0時(shí)，可行場景的區(qū)域最大，此時(shí)該相角達(dá)到的最大值（最小值）是所有場景中的最大值（最小值）。該相角達(dá)到最大值時(shí)，兩風(fēng)電出力處于最大值，節(jié)點(diǎn)3 負(fù)荷處于最小值；該相角達(dá)到最小值時(shí)，兩風(fēng)電出力處于最小值，節(jié)點(diǎn)3 負(fù)荷處于最大值。在相關(guān)系數(shù)為正，如圖6（b）中ρwl=0.2時(shí)，ρwl=0 時(shí)取得最大值或最小值的場景不在其可行區(qū)域內(nèi)，因此在對(duì)應(yīng)相關(guān)系數(shù)下的可行區(qū)域內(nèi)尋找新的符合要求的場景。隨著正相關(guān)性的增強(qiáng)，可行場景區(qū)域逐漸縮小，且無法取得相關(guān)性較弱時(shí)的最大值或最小值場景，因此相角的區(qū)間范圍也逐漸縮小。而在相關(guān)系數(shù)為負(fù)時(shí)，雖然可行的場景區(qū)域縮小，但在對(duì)應(yīng)相關(guān)系數(shù)下仍能取得ρwl=0 時(shí)使該相角達(dá)到最大值和最小值的場景，如附錄E 所示，因此相角區(qū)間并未變化。

圖6 不同正相關(guān)系數(shù)下節(jié)點(diǎn)3電壓相角取值場景Fig.6 Scenes of voltage angle for Bus 3 under different positive correlation coefficients

4 結(jié)論

本文通過平行四邊形模型構(gòu)造變量間的相關(guān)性區(qū)域，在場景優(yōu)化法區(qū)間潮流計(jì)算的基礎(chǔ)上，采用靈敏度分類法提取具有相同變化趨勢的潮流變量一同計(jì)算，算例分析結(jié)果表明：

1）用靈敏度分類法將變化趨勢相同的狀態(tài)變量進(jìn)行歸類計(jì)算，能在誤差很小的情況下減少非線性優(yōu)化次數(shù)，提升計(jì)算效率；

2）相關(guān)性的強(qiáng)弱決定了區(qū)間變量的分布區(qū)域，相關(guān)性越強(qiáng)則分布區(qū)域越小，考慮相關(guān)性的區(qū)間潮流計(jì)算排除了一些不可能出現(xiàn)的場景，得到的計(jì)算結(jié)果區(qū)間范圍更貼近實(shí)際。

附錄見本刊網(wǎng)絡(luò)版（http：//www.epae.cn）。