搞配電網(wǎng)潮流計(jì)算的人遲早會(huì)遇到同一個(gè)問(wèn)題光伏和風(fēng)電一接入注入功率變成隨機(jī)量傳統(tǒng)的確定性潮流結(jié)果還靠得住嗎我做這個(gè)項(xiàng)目的原因就是這樣——在IEEE33節(jié)點(diǎn)系統(tǒng)上用蒙特卡洛法把風(fēng)電和光伏出力的隨機(jī)性一層層“撒”進(jìn)潮流計(jì)算里跑完幾千次抽樣之后再去統(tǒng)計(jì)節(jié)點(diǎn)電壓和支路潮流的分布特性。整個(gè)過(guò)程用Matlab實(shí)現(xiàn)代碼量不算大但坑不少。這篇把我的實(shí)現(xiàn)思路、參數(shù)設(shè)置、代碼結(jié)構(gòu)和踩坑記錄完整寫出來(lái)給準(zhǔn)備做概率潮流Probabilistic Load FlowPLF的同行參考。如果你是剛接觸概率潮流的碩士生或者正在做分布式電源接入評(píng)估的工程師這篇內(nèi)容應(yīng)該能幫你省不少時(shí)間。蒙特卡洛法本身不復(fù)雜真正容易翻車的地方全在建模細(xì)節(jié)和統(tǒng)計(jì)口徑上我會(huì)把這兩塊重點(diǎn)展開(kāi)。1. 概率潮流到底在算什么確定性潮流的一個(gè)盲區(qū)1.1 確定性潮流的“確定性”是假設(shè)出來(lái)的傳統(tǒng)潮流計(jì)算的輸入是一組確定的負(fù)荷數(shù)據(jù)、確定的發(fā)電出力然后通過(guò)牛頓-拉夫遜法或者前推回代法解出一組確定的節(jié)點(diǎn)電壓和支路潮流。這個(gè)流程本身沒(méi)有任何問(wèn)題問(wèn)題出在“輸入是確定”這個(gè)前提上。在傳統(tǒng)配電網(wǎng)里負(fù)荷雖然有波動(dòng)但波動(dòng)范圍和規(guī)律相對(duì)穩(wěn)定電網(wǎng)調(diào)度人員用峰值負(fù)荷或典型負(fù)荷場(chǎng)景來(lái)校核基本夠用?？梢坏┙尤腼L(fēng)電和光伏事情就變了。風(fēng)電出力隨風(fēng)速波動(dòng)光伏出力隨光照強(qiáng)度波動(dòng)而且這些波動(dòng)幅度非常大——光伏在晴天中午可能是額定出力到傍晚直接歸零風(fēng)電可能前一個(gè)小時(shí)滿發(fā)后一個(gè)小時(shí)因?yàn)轱L(fēng)速超過(guò)切出風(fēng)速而停機(jī)。這時(shí)候再拿單一場(chǎng)景去做潮流計(jì)算相當(dāng)于只看到了系統(tǒng)成千上萬(wàn)種運(yùn)行狀態(tài)里的一張快照。舉個(gè)直觀的例子在IEEE33節(jié)點(diǎn)系統(tǒng)里如果只按光伏額定出力計(jì)算一次潮流末端節(jié)點(diǎn)電壓可能正好在合格范圍之內(nèi)但如果光照強(qiáng)度下降一半同樣一個(gè)接入點(diǎn)的電壓就可能會(huì)出現(xiàn)越限風(fēng)險(xiǎn)。單次確定性潮流完全無(wú)法回答“這種越限情況發(fā)生的概率有多大”這個(gè)問(wèn)題。1.2 概率潮流解決什么問(wèn)題概率潮流的思路是把輸入隨機(jī)量風(fēng)速、光照、負(fù)荷的概率分布作為已知條件讓它們參與計(jì)算最終輸出節(jié)點(diǎn)電壓、支路潮流、網(wǎng)損這些結(jié)果量的概率分布和統(tǒng)計(jì)特征。換句話說(shuō)確定性潮流回答的是“在這種情況下電壓是多少”概率潮流回答的是“電壓落在某個(gè)范圍內(nèi)的概率有多大”“電壓期望值是多少”“越限概率有多高”。對(duì)分布式電源接入評(píng)估來(lái)說(shuō)后者才是真正有工程意義的答案。實(shí)現(xiàn)概率潮流最經(jīng)典的方法就是蒙特卡洛法。思路很直白按照風(fēng)速、光照的概率分布反復(fù)抽樣每次抽樣得到一組風(fēng)光出力組合到負(fù)荷數(shù)據(jù)里做一次確定性潮流重復(fù)幾千次之后把這幾千次潮流結(jié)果匯總起來(lái)做統(tǒng)計(jì)分析。抽樣次數(shù)足夠多時(shí)統(tǒng)計(jì)結(jié)果就會(huì)逼近真實(shí)分布。1.3 蒙特卡洛與其他方法怎么選我剛開(kāi)始做的時(shí)候也猶豫過(guò)是不是用點(diǎn)估計(jì)法或者Gram-Charlier級(jí)數(shù)展開(kāi)會(huì)更“高級(jí)”一些后來(lái)對(duì)比下來(lái)還是老老實(shí)實(shí)用了蒙特卡洛原因在于它有三個(gè)別人比不了的優(yōu)勢(shì)方法核心思路優(yōu)點(diǎn)缺點(diǎn)蒙特卡洛法大量抽樣后統(tǒng)計(jì)模型不受限制、原理簡(jiǎn)單、精度可控計(jì)算量大需要大量樣本點(diǎn)估計(jì)法用少量特征點(diǎn)近似速度快對(duì)強(qiáng)非線性問(wèn)題誤差偏大級(jí)數(shù)展開(kāi)法用多項(xiàng)式逼近輸出分布解析性好遇到非光滑模型容易失真蒙特卡洛最大的底氣是“沒(méi)有模型假設(shè)”。風(fēng)光出力模型、負(fù)荷波動(dòng)模型寫成什么樣都行哪怕是非線性、非連續(xù)、分段函數(shù)它也能直接跑。而解析類方法往往要對(duì)系統(tǒng)做線性化或者各種簡(jiǎn)化處理碰到風(fēng)機(jī)出力這種典型的非線性分段函數(shù)就容易出偏差。代價(jià)就是計(jì)算量大。一次確定性潮流很快幾千次疊加起來(lái)就要考慮效率。這個(gè)后面在程序優(yōu)化部分我會(huì)單獨(dú)講。2. 風(fēng)電和光伏出力建模把隨機(jī)性變成數(shù)學(xué)分布蒙特卡洛的結(jié)果準(zhǔn)不準(zhǔn)關(guān)鍵不在抽樣次數(shù)多少而在輸入分布建模靠不靠譜。如果你把風(fēng)速分布參數(shù)都搞錯(cuò)了跑一萬(wàn)次抽樣也沒(méi)意義只能得到一份精確的錯(cuò)誤結(jié)果。2.1 風(fēng)速兩參數(shù)Weibull分布描述風(fēng)速最常用的模型是兩參數(shù)Weibull分布概率密度函數(shù)是f(v) (k/c) * (v/c)^(k-1) * exp(-(v/c)^k)其中v是風(fēng)速k是形狀參數(shù)c是尺度參數(shù)。c和平均風(fēng)速直接相關(guān)k則決定了風(fēng)速分布的形態(tài)k越大風(fēng)速越集中k越小風(fēng)速波動(dòng)越劇烈。典型的風(fēng)速場(chǎng)景里k大致在1.5到3之間沿海開(kāi)闊地區(qū)偏大內(nèi)陸或者地形復(fù)雜一點(diǎn)的地方偏小。參數(shù)怎么定如果有歷史風(fēng)速數(shù)據(jù)直接用Matlab的wblfit函數(shù)做擬合最靠譜。如果沒(méi)有數(shù)據(jù)就用典型參數(shù)。常見(jiàn)做法是給定平均風(fēng)速Vm和標(biāo)準(zhǔn)差σ然后用經(jīng)驗(yàn)公式反推k (σ / Vm)^(-1.086)c Vm / Γ(1 1/k)這里Γ是伽馬函數(shù)Matlab里用gamma(11/k)。實(shí)測(cè)下來(lái)如果你只知道一個(gè)平均風(fēng)速取k2.0再按c ≈ 1.128*Vm估算偏差也不會(huì)太大。有了參數(shù)生成風(fēng)速樣本很簡(jiǎn)單Matlab直接調(diào)用k 2.0; c 7.0; N 5000; v wblrnd(k, c, N, 1);這里生成的v就是N次蒙特卡洛抽樣對(duì)應(yīng)的風(fēng)速值。接下來(lái)是風(fēng)機(jī)出力轉(zhuǎn)換。風(fēng)機(jī)不是所有風(fēng)速段都出力出力特性一般用分段函數(shù)描述function P wind_power(v, Pn, vci, vr, vco) P zeros(size(v)); % 切入風(fēng)速以下、切出風(fēng)速以上出力為0 idx_zero (v vci) | (v vco); P(idx_zero) 0; % 切入風(fēng)速到額定風(fēng)速之間近似線性爬升 idx_linear (v vci) (v vr); P(idx_linear) Pn * (v(idx_linear) - vci) / (vr - vci); % 額定風(fēng)速到切出風(fēng)速之間保持額定功率 idx_rated (v vr) (v vco); P(idx_rated) Pn; end典型參數(shù)切入風(fēng)速3 m/s、額定風(fēng)速12 m/s、切出風(fēng)速25 m/s。這幾個(gè)參數(shù)不同廠家的風(fēng)機(jī)不完全一樣但大體在這個(gè)范圍。額定功率Pn按接入的風(fēng)機(jī)容量設(shè)比如單臺(tái)2MW。注意切入風(fēng)速那里我做項(xiàng)目時(shí)喜歡把小于等于vci的區(qū)間也歸零避免邊界值計(jì)算出微小正數(shù)影響統(tǒng)計(jì)結(jié)果。2.2 光照強(qiáng)度Beta分布與光伏出力光伏出力的隨機(jī)性源頭是光照強(qiáng)度。與風(fēng)速不同光照強(qiáng)度是有上下界的白天從零到最大值之間變化歸一化之后落在[0,1]區(qū)間。這種有界分布最適合用Beta分布描述。Beta分布由兩個(gè)形狀參數(shù)a和b決定。實(shí)際建模時(shí)通常從統(tǒng)計(jì)資料里獲得光照強(qiáng)度的均值μ和方差σ2然后反推參數(shù)a μ * (μ * (1 - μ) / σ2 - 1)b (1 - μ) * (μ * (1 - μ) / σ2 - 1)舉個(gè)例子假設(shè)歸一化光照強(qiáng)度的均值是0.45方差是0.02那么μ(1-μ)0.2475μ(1-μ)/σ212.375算出來(lái)a約5.1b約6.3。把這個(gè)參數(shù)代入betarnd就能生成光照樣本mu 0.45; sigma2 0.02; a mu * (mu * (1 - mu) / sigma2 - 1); b (1 - mu) * (mu * (1 - mu) / sigma2 - 1); r betarnd(a, b, N, 1); % 歸一化光照強(qiáng)度光伏出力模型我建議第一版做簡(jiǎn)化處理P_pv Ppv_rated * r也就是讓光伏出力等于額定容量乘以歸一化光照強(qiáng)度。這樣參數(shù)少、調(diào)試方便。如果想精細(xì)一點(diǎn)可以加入溫度和組件效率的影響用P_pv η * S * G * r其中η是光電轉(zhuǎn)換效率S是光伏板總面積G是標(biāo)準(zhǔn)光照強(qiáng)度一般取1000 W/m2。但第一版別把模型搞太復(fù)雜先把蒙特卡洛框架跑通再回頭加細(xì)節(jié)。2.3 負(fù)荷隨機(jī)性處理很多概率潮流文章只把風(fēng)光當(dāng)成隨機(jī)量負(fù)荷完全固定。實(shí)際系統(tǒng)里負(fù)荷本身也有波動(dòng)尤其是配電網(wǎng)的負(fù)荷曲線波動(dòng)幅度不小。我的做法是在負(fù)荷額定值基礎(chǔ)上加一個(gè)正態(tài)擾動(dòng)load_peak [0; 0.1; 0.09; 0.12; 0.063; 0.063; 0.02; 0.02; ...]; % 各節(jié)點(diǎn)峰值負(fù)荷 load_i load_peak .* (1 0.05 * randn(33, 1));標(biāo)準(zhǔn)差系數(shù)取0.05也就是負(fù)荷在額定值上下約5%的范圍波動(dòng)。這個(gè)是經(jīng)驗(yàn)值實(shí)際可以根據(jù)歷史負(fù)荷數(shù)據(jù)調(diào)整。加這段的好處是不用額外增加代碼量就能讓概率潮流的場(chǎng)景更接近實(shí)際情況。這里還要提一個(gè)容易漏掉的問(wèn)題風(fēng)速、光照、負(fù)荷之間可能存在相關(guān)性。比如同一個(gè)地區(qū)多云天氣時(shí)光照降低但風(fēng)速可能增大。第一版可以先假設(shè)它們互相獨(dú)立因?yàn)槊商乜宓莫?dú)立性假設(shè)最容易實(shí)現(xiàn)如果要做精細(xì)化分析再引入Copula函數(shù)處理相關(guān)性。3. IEEE33節(jié)點(diǎn)系統(tǒng)準(zhǔn)備數(shù)據(jù)、基準(zhǔn)值與接入方案3.1 系統(tǒng)基礎(chǔ)數(shù)據(jù)怎么整理IEEE33節(jié)點(diǎn)系統(tǒng)是配電網(wǎng)研究里最常用的標(biāo)準(zhǔn)算例之一結(jié)構(gòu)是單饋線輻射狀網(wǎng)絡(luò)額定電壓12.66 kV包含33個(gè)節(jié)點(diǎn)、32條支路系統(tǒng)總負(fù)荷約3715 kW和2300 kvar。它之所以被廣泛使用是因?yàn)榻Y(jié)構(gòu)適中、輻射狀特征典型、末端電壓?jiǎn)栴}明顯非常適合驗(yàn)證分布式電源接入的影響。在Matlab里我習(xí)慣用兩個(gè)矩陣存基礎(chǔ)數(shù)據(jù)。第一個(gè)是負(fù)荷矩陣格式是[節(jié)點(diǎn)編號(hào)有功負(fù)荷(kW)無(wú)功負(fù)荷(kvar)]第二個(gè)是支路矩陣格式是[首端節(jié)點(diǎn)末端節(jié)點(diǎn)支路電阻(Ω)支路電抗(Ω)]。IEEE33節(jié)點(diǎn)的原始數(shù)據(jù)在很多教材和論文附錄里都能找到網(wǎng)上也有公開(kāi)版本但整理時(shí)一定要注意兩個(gè)細(xì)節(jié)第一原數(shù)據(jù)里支路阻抗單位是歐姆已經(jīng)是完整線路的阻抗值不需要再乘長(zhǎng)度。第二節(jié)點(diǎn)負(fù)荷有功無(wú)功單位是kW和kvar和基準(zhǔn)容量統(tǒng)一之前別直接拿來(lái)算標(biāo)幺值。我建議在程序開(kāi)頭定義一個(gè)基準(zhǔn)容量SBIEEE33系統(tǒng)一般取10 MVA比較順手?；鶞?zhǔn)電壓取12.66 kV這樣系統(tǒng)總負(fù)荷標(biāo)幺值大約0.3715 j0.23電壓標(biāo)幺值在1.0附近波動(dòng)物理意義直觀。3.2 風(fēng)電、光伏節(jié)點(diǎn)和容量怎么定風(fēng)光接入位置不是隨便拍的。IEEE33的末端節(jié)點(diǎn)電壓本來(lái)就偏低把風(fēng)機(jī)和光伏接在末端能最大程度地展現(xiàn)新能源對(duì)電壓分布的影響。我第一版做的時(shí)候把風(fēng)機(jī)接在17號(hào)節(jié)點(diǎn)一條較長(zhǎng)分支的末端光伏接在33號(hào)節(jié)點(diǎn)饋線最末端這兩個(gè)位置都是典型的電壓薄弱點(diǎn)。容量方面我給風(fēng)機(jī)取2MW額定光伏取1.5MW。相對(duì)系統(tǒng)總負(fù)荷約3.7MW來(lái)看滲透率已經(jīng)比較高了接在末端時(shí)電壓波動(dòng)會(huì)非常明顯做出來(lái)的概率潮流曲線也能看出明顯的分布形態(tài)變化。如果你只是驗(yàn)證程序跑通容量也可以調(diào)小一些但那樣電壓分布可能和確定性潮流差別不大不利于觀察。接入方式上第一版建議風(fēng)力和光伏都按功率因數(shù)1.0處理也就是只發(fā)有功不出無(wú)功。這樣模型最簡(jiǎn)單結(jié)果也能清晰地反映有功注入對(duì)配電網(wǎng)電壓的影響。等框架跑通了再改成恒功率因數(shù)0.95或者恒電壓控制都不遲。3.3 基準(zhǔn)值統(tǒng)一與單位換算這個(gè)環(huán)節(jié)是新手最容易踩坑的地方。IEEE33原始數(shù)據(jù)的節(jié)點(diǎn)電壓是12.66kV支路阻抗是歐姆負(fù)荷是kW和kvar但潮流計(jì)算內(nèi)部用的是標(biāo)幺值體系。你需要把支路阻抗除以基準(zhǔn)阻抗Zb Ub2/SB把負(fù)荷功率除以基準(zhǔn)功率SB轉(zhuǎn)換完成才能進(jìn)入計(jì)算流程。我這里習(xí)慣直接在數(shù)據(jù)讀入階段完成換算而不是在潮流函數(shù)里再做一遍這樣能保證潮流函數(shù)本身保持簡(jiǎn)潔。換算公式分別是Z_pu Z_Ω / Zb其中Zb 12.662 / 10 ≈ 16 ΩP_pu P_kW / 10000如果你用我之前說(shuō)的SB10MVA那么Zb就是約16.02Ω。這個(gè)數(shù)值不需要死記程序里用 (12.66^2 / 10) 一行代碼就能算出來(lái)。我自己曾經(jīng)因?yàn)樵诖a里把基準(zhǔn)容量寫成了100MVA結(jié)果阻抗標(biāo)幺值整體小了10倍算出來(lái)的電壓曲線全都漂移查了半個(gè)多小時(shí)才定位到問(wèn)題。4. 蒙特卡洛概率潮流程序架構(gòu)與核心實(shí)現(xiàn)4.1 整體流程說(shuō)明整個(gè)蒙特卡洛概率潮流程序我拆成了六個(gè)步驟初始化錄入IEEE33系統(tǒng)數(shù)據(jù)完成標(biāo)幺值換算。設(shè)置參數(shù)抽樣次數(shù)N、風(fēng)機(jī)參數(shù)、光伏參數(shù)、負(fù)荷波動(dòng)系數(shù)。抽樣循環(huán)每次迭代生成風(fēng)速和光照樣本換算成風(fēng)光出力疊加負(fù)荷隨機(jī)擾動(dòng)。修改節(jié)點(diǎn)注入功率把風(fēng)電光伏出力和負(fù)荷功率合并成節(jié)點(diǎn)凈注入功率。調(diào)用確定性潮流函數(shù)求解本次樣本的節(jié)點(diǎn)電壓和支路潮流。循環(huán)結(jié)束后匯總所有樣本用mean、std、histogram等函數(shù)做統(tǒng)計(jì)分析。整個(gè)流程用函數(shù)化組織比較清晰。主程序里不要塞太多東西我一般分成四個(gè)函數(shù)文件主程序、確定性潮流函數(shù)、風(fēng)電出力函數(shù)、光伏出力函數(shù)。這樣后面改參數(shù)、換器件模型都方便。4.2 前推回代潮流函數(shù)實(shí)現(xiàn)IEEE33是輻射狀配電網(wǎng)用前推回代法最合適。這個(gè)方法不需要形成雅可比矩陣原理簡(jiǎn)單先假設(shè)所有節(jié)點(diǎn)電壓為1.0 p.u.從末端節(jié)點(diǎn)往電源側(cè)回推根據(jù)負(fù)荷功率計(jì)算各支路的電流和功率然后從電源側(cè)向末端前推根據(jù)支路阻抗重新計(jì)算節(jié)點(diǎn)電壓反復(fù)迭代直到兩次迭代的電壓差值小于閾值。核心代碼片段如下function V backward_forward(branch, S_load, slack_idx, tol, max_iter) n size(S_load, 1); V ones(n, 1); Z branch(:, 3) 1j * branch(:, 4); for iter 1:max_iter % 末端向首端回推計(jì)算各支路電流 I zeros(size(Z)); for k size(branch, 1):-1:1 child branch(k, 2); parent branch(k, 1); I(k) conj(S_load(child) / V(child)); end % 首端向末端前推更新電壓 V_new ones(n, 1); for k 1:size(branch, 1) parent branch(k, 1); child branch(k, 2); V_new(child) V_new(parent) - Z(k) * I(k); end if max(abs(abs(V_new) - abs(V))) tol V V_new; break; end V V_new; end end收斂閾值我取1e-6最大迭代次數(shù)設(shè)50次。輻射狀配電網(wǎng)收斂非常快一般十幾次迭代就能滿足。要注意的是回推和前推的支路順序必須和網(wǎng)絡(luò)拓?fù)浞较蛞恢氯绻肪仃図樞騺y了收斂速度會(huì)顯著變慢甚至不收斂。4.3 蒙特卡洛主循環(huán)與統(tǒng)計(jì)主循環(huán)是程序的核心。我建議事先想好需要保存哪些量用矩陣預(yù)先分配內(nèi)存而不是在循環(huán)里不斷拼接否則跑到后面會(huì)慢到你懷疑人生。我一般保存兩類結(jié)果一是每次抽樣得到的33個(gè)節(jié)點(diǎn)電壓幅值二是特定支路的潮流功率。主循環(huán)大致長(zhǎng)這樣N 5000; V_all zeros(N, 33); Ploss_all zeros(N, 1); wind_bus 17; pv_bus 33; for i 1:N % 風(fēng)電出力 v_wind wblrnd(k, c); P_w wind_power(v_wind, Pn_wind, vci, vr, vco); % 光伏出力 r_solar betarnd(a, b); P_pv Pn_pv * r_solar; % 負(fù)荷波動(dòng) S_load load_base .* (1 sigma_load * randn(33, 1)); S_load S_load / SB; % 把風(fēng)電光伏出力疊加到對(duì)應(yīng)節(jié)點(diǎn)注入功率為負(fù)的負(fù)荷 S_load(wind_bus) S_load(wind_bus) - P_w / SB; S_load(pv_bus) S_load(pv_bus) - P_pv / SB; % 潮流計(jì)算 V backward_forward(branch_pu, S_load, 1, 1e-6, 50); V_all(i, :) abs(V).; % 順便算一下網(wǎng)損 % Ploss sum(real(S_load .* conj(V))) ... 這里按公式補(bǔ)全 end這段代碼里S_load的單位全部是標(biāo)幺值所以負(fù)荷除基準(zhǔn)容量風(fēng)光出力也除基準(zhǔn)容量。細(xì)節(jié)上節(jié)點(diǎn)1是平衡節(jié)點(diǎn)不用改它的注入功率直接固定為1.0 p.u.電壓即可。抽樣次數(shù)N我先用5000。蒙特卡洛估計(jì)的標(biāo)準(zhǔn)誤差大致正比于1/sqrt(N)也就是說(shuō)如果N從1000增加到4000標(biāo)準(zhǔn)差降低一半。對(duì)于觀察電壓期望值5000次完全夠用如果要做尾部概率分析比如精確估計(jì)電壓越限概率建議加到10000次以上后面效率優(yōu)化部分我再細(xì)說(shuō)。4.4 結(jié)果存儲(chǔ)設(shè)計(jì)統(tǒng)計(jì)計(jì)算是蒙特卡洛里看似簡(jiǎn)單但最容易出錯(cuò)的部分。你可能想著跑完了直接用mean(V_all)一行代碼搞定但實(shí)際上還有個(gè)收斂性驗(yàn)證問(wèn)題你怎么知道5000次夠了我的做法是畫(huà)一條“均值隨抽樣次數(shù)變化曲線”看它在哪個(gè)位置開(kāi)始穩(wěn)定下來(lái)cum_mean cumsum(V_all(:, 17)) ./ (1:N); plot(1:N, cum_mean);如果曲線在前面大幅震蕩、后面逐漸平緩說(shuō)明結(jié)果穩(wěn)定了如果到N5000還在明顯漂移說(shuō)明樣本數(shù)不夠得加大N。這種方式比拍腦袋定次數(shù)直觀得多。5. 結(jié)果怎么解讀電壓分布、支路潮流與越限概率5.1 節(jié)點(diǎn)電壓期望曲線與波動(dòng)包絡(luò)第一張圖我推薦畫(huà)“節(jié)點(diǎn)電壓期望值曲線”橫軸是節(jié)點(diǎn)編號(hào)縱軸是電壓標(biāo)幺值。你能直觀看到電壓從電源側(cè)到饋線末端逐漸下降的趨勢(shì)也能看到接入風(fēng)電光伏之后局部節(jié)點(diǎn)電壓被抬升的幅度。第二張圖是電壓波動(dòng)包絡(luò)畫(huà)每個(gè)節(jié)點(diǎn)電壓的均值加兩倍標(biāo)準(zhǔn)差和均值減兩倍標(biāo)準(zhǔn)差形成一個(gè)條帶。IEEE33系統(tǒng)的原始網(wǎng)絡(luò)里末端節(jié)點(diǎn)電壓一般會(huì)降到0.95以下或者接近0.95接入風(fēng)光之后整條曲線的條帶會(huì)明顯變寬這就是隨機(jī)性帶來(lái)的實(shí)際效果。畫(huà)包絡(luò)線的代碼很簡(jiǎn)單V_mean mean(V_all, 1); V_std std(V_all, 0, 1); plot(V_mean, b, LineWidth, 1.5); hold on; plot(V_mean 2*V_std, r--); plot(V_mean - 2*V_std, r--);5.2 特定節(jié)點(diǎn)電壓概率密度與越限概率整個(gè)系統(tǒng)的節(jié)點(diǎn)太多逐個(gè)看分布不現(xiàn)實(shí)。我一般只挑兩三個(gè)關(guān)鍵節(jié)點(diǎn)接入風(fēng)電的17號(hào)節(jié)點(diǎn)、接入光伏的33號(hào)節(jié)點(diǎn)、饋線中段的一個(gè)節(jié)點(diǎn)。對(duì)33號(hào)節(jié)點(diǎn)做histogram(V_all(:,33), Normalization, pdf);然后疊加一條正態(tài)分布擬合曲線就會(huì)發(fā)現(xiàn)電壓分布并不完全對(duì)稱可能有輕微偏態(tài)。原因在于輸入的光照Beta分布本身不對(duì)稱再加上風(fēng)機(jī)出力分段函數(shù)的非線性輸出電壓分布不是嚴(yán)格正態(tài)是必然的。這個(gè)觀察本身就有意義說(shuō)明概率潮流用蒙特卡洛直接模擬比那些假設(shè)輸出服從正態(tài)分布的解析方法更可靠。越限概率的計(jì)算需要注意寫法。按國(guó)標(biāo)配電網(wǎng)電壓合格范圍一般是0.93到1.07 p.u.但我做IEEE33時(shí)習(xí)慣按0.95到1.05 p.u. 的常用工程范圍來(lái)統(tǒng)計(jì)。代碼就是常規(guī)的比較計(jì)數(shù)violation_prob mean(V_all(:, 33) 0.95 | V_all(:, 33) 1.05);5.3 支路潮流與網(wǎng)損的統(tǒng)計(jì)特征除了電壓支路潮流也是重要輸出。我重點(diǎn)統(tǒng)計(jì)兩支一條是饋線首段支路節(jié)點(diǎn)1到2另一條是光伏接入節(jié)點(diǎn)所在的末段支路。首段支路承載全部系統(tǒng)功率它的潮流分布能反映出系統(tǒng)總負(fù)荷和總發(fā)電的隨機(jī)波動(dòng)末段支路則直接體現(xiàn)光伏接入對(duì)局部過(guò)載的影響。網(wǎng)損的統(tǒng)計(jì)分布同樣值得看。每次抽樣中網(wǎng)損等于所有支路損耗之和用Ploss_all向量保存最后統(tǒng)計(jì)期望和95%分位數(shù)。接入高滲透率分布式電源后網(wǎng)損的期望值通常會(huì)下降但極端場(chǎng)景下可能出現(xiàn)比無(wú)分布式電源時(shí)更高的網(wǎng)損這種“尾部風(fēng)險(xiǎn)”可以用分位數(shù)直觀暴露出來(lái)。一個(gè)比較實(shí)用的匯總表可以在結(jié)果輸出階段打印出來(lái)統(tǒng)計(jì)量節(jié)點(diǎn)17電壓(p.u.)節(jié)點(diǎn)33電壓(p.u.)首段支路潮流(kVA)網(wǎng)損(kW)確定性潮流值0.960.93約2500約120期望值0.980.97約2400約95標(biāo)準(zhǔn)差0.0120.023約180約3095%分位數(shù)1.001.01約2700約140注意表格里的數(shù)值是示例真實(shí)輸出取決于你的容量和參數(shù)設(shè)置。6. 踩坑記錄從風(fēng)速參數(shù)到迭代發(fā)散6.1 典型問(wèn)題速查表我跑這個(gè)程序至少踩了七八個(gè)坑挑有代表性的整理成速查表問(wèn)題現(xiàn)象可能原因解決辦法電壓期望曲線和確定性潮流完全重合隨機(jī)波動(dòng)太弱或者隨機(jī)數(shù)種子固定導(dǎo)致每次抽樣相同檢查負(fù)荷標(biāo)準(zhǔn)差和風(fēng)光容量是否過(guò)小給隨機(jī)數(shù)加不同種子末端電壓超過(guò)1.05 p.u. 甚至1.1光伏容量過(guò)大、功率因數(shù)等于1導(dǎo)致電壓抬升降低容量或者設(shè)置逆變器吸收無(wú)功光照Beta分布參數(shù)算出來(lái)是負(fù)數(shù)標(biāo)準(zhǔn)差設(shè)置太大不符合有界分布約束檢查σ2是否小于μ(1-μ)前推回代迭代次數(shù)超過(guò)上限支路矩陣順序錯(cuò)亂或基準(zhǔn)阻抗算錯(cuò)核對(duì)支路拓?fù)漤樞蚝蚙b換算每次運(yùn)行結(jié)果都不一樣蒙特卡洛隨機(jī)抽樣本該如此用rng(固定值)固定種子方便重現(xiàn)網(wǎng)損出現(xiàn)負(fù)值潮流方向定義和網(wǎng)損公式不一致檢查支路潮流方向統(tǒng)一功率正方向跑10000次太慢循環(huán)內(nèi)有重復(fù)計(jì)算或未預(yù)分配內(nèi)存用parfor并行提前分配矩陣空間這里最坑的是第二個(gè)問(wèn)題光伏接在33號(hào)節(jié)點(diǎn)容量稍微取大一點(diǎn)末端電壓大概率直接越上限。這不是程序?qū)戝e(cuò)了而是系統(tǒng)真實(shí)的物理響應(yīng)——分布式電源向末端注入有功末端電壓被抬升。做結(jié)果分析時(shí)不要急著改代碼先想想物理規(guī)律對(duì)不對(duì)。6.2 效率優(yōu)化抽樣次數(shù)、并行與方差縮減蒙特卡洛最大的軟肋是計(jì)算速度。IEEE33系統(tǒng)規(guī)模小跑5000次大概也就是幾秒到十幾秒但如果以后換到IEEE123節(jié)點(diǎn)系統(tǒng)或者加上負(fù)荷模型細(xì)節(jié)時(shí)間會(huì)明顯拉長(zhǎng)。我的優(yōu)化經(jīng)驗(yàn)按性價(jià)比從高到低排列第一用parfor替換for。Matlab的并行計(jì)算池在函數(shù)里調(diào)用時(shí)用parfor能直接利用多核。注意preallocated矩陣V_all的寫入方式不同迭代各寫一行不會(huì)沖突適合并行。第二減少重復(fù)計(jì)算。潮流函數(shù)里如果每次都要重新算一遍基準(zhǔn)阻抗、支路阻抗等固定數(shù)據(jù)那會(huì)浪費(fèi)大量時(shí)間。這些值在主程序里提前算好以結(jié)構(gòu)體或全局參數(shù)形式傳入函數(shù)。第三用拉丁超立方抽樣LHS替代簡(jiǎn)單隨機(jī)抽樣。LHS的核心思想是把每個(gè)輸入變量的分布區(qū)間分成N等份在每個(gè)區(qū)間里采樣然后打亂組合保證覆蓋更均勻。在相同樣本數(shù)下LHS的統(tǒng)計(jì)精度通常比直接抽樣更高尤其是對(duì)Beta分布這種有界分布效果更明顯。Matlab的lhsdesign或者Statistics Toolbox里的lhsnorm都可以實(shí)現(xiàn)。第四不要太貪心。如果只是看期望值和標(biāo)準(zhǔn)差2000到5000次夠了如果關(guān)心的是極低概率事件比如萬(wàn)分之一的越限概率那N至少要幾萬(wàn)次這時(shí)候最好用方差縮減技術(shù)而不是單純加樣本量。7. 最后再聊幾句蒙特卡洛的取舍與擴(kuò)展方向跑完這個(gè)項(xiàng)目之后我最大的體會(huì)是蒙特卡洛法在概率潮流里不是最聰明的算法卻是最適合快速上手的算法。它的強(qiáng)大不在于精巧而在于“讓數(shù)據(jù)說(shuō)話”的直覺(jué)——把能想到的隨機(jī)性都扔進(jìn)模型讓無(wú)數(shù)的場(chǎng)景自己碰撞出結(jié)果。對(duì)于配電網(wǎng)這種結(jié)構(gòu)規(guī)模不算太大、但模型細(xì)節(jié)繁多的場(chǎng)景這種簡(jiǎn)單粗暴反而救了不少人。如果你打算在這個(gè)項(xiàng)目基礎(chǔ)上繼續(xù)深入我建議按三個(gè)方向擴(kuò)展。第一個(gè)方向是時(shí)序化把靜態(tài)分布擴(kuò)展成典型日24小時(shí)的風(fēng)速光照時(shí)序曲線看看一天里哪些時(shí)段電壓波動(dòng)最危險(xiǎn)。第二個(gè)方向是相關(guān)性建模用Copula函數(shù)把風(fēng)速和光照的相關(guān)性加進(jìn)去往往會(huì)讓結(jié)果更貼近實(shí)際情況。第三個(gè)方向是接入逆變器控制策略讓光伏和風(fēng)機(jī)在電壓越限時(shí)參與無(wú)功調(diào)節(jié)然后再次做概率潮流對(duì)比前后電壓分布的變化——這個(gè)方向在做分布式電源容量規(guī)劃時(shí)很有實(shí)用價(jià)值。最后再分享一個(gè)小經(jīng)驗(yàn)概率潮流做完之后永遠(yuǎn)要問(wèn)自己一句結(jié)果的物理規(guī)律是否講得通。蒙特卡洛輸出的是一堆數(shù)字但數(shù)字背后的合理性要靠領(lǐng)域知識(shí)去判斷。如果出力增大時(shí)某個(gè)節(jié)點(diǎn)電壓反而下降了或者網(wǎng)損分布出現(xiàn)了完全反常規(guī)的形狀大概率不是隨機(jī)性太強(qiáng)而是你的某個(gè)建模環(huán)節(jié)出了問(wèn)題。先把物理捋順了再讓蒙特卡洛替你放大細(xì)節(jié)這樣才能做出真正可信的分析結(jié)果。