風(fēng)電、光伏大規(guī)模接入之后電網(wǎng)分析的默認假設(shè)就變了。以前做潮流計算給一組確定的發(fā)電機出力和負荷牛拉法迭代一遍得到一組節(jié)點電壓和支路功率這套流程在傳統(tǒng)火電主導(dǎo)的時代夠用但當(dāng)出力看天吃飯的新能源占比上來輸入側(cè)的確定性假設(shè)本身就站不住腳。概率潮流計算Probabilistic Load FlowPLF正是用來應(yīng)對這種不確定性的工具——它把風(fēng)速、光照、負荷都建模成隨機變量通過大規(guī)模采樣或解析近似得到節(jié)點電壓和支路功率的概率分布而不是一個孤零零的點值。Matlab憑借矩陣計算優(yōu)勢和自帶的統(tǒng)計工具箱一直是做概率潮流最順手的平臺。這篇文章從工程應(yīng)用的角度把含風(fēng)光發(fā)電的概率潮流計算的數(shù)學(xué)模型、三種主流算法、Matlab代碼框架和調(diào)試經(jīng)驗完整過一遍剛接觸這個方向的研究生和做新能源接入評估的工程師都能直接拿去參考。1. 概率潮流到底在算什么從單點結(jié)果到概率分布的思維轉(zhuǎn)換1.1 確定性潮流的“天花板”在哪里傳統(tǒng)的確定性潮流求解的是這樣一組方程節(jié)點注入功率等于電壓與導(dǎo)納矩陣的乘積方程組給定以后迭代求解出唯一一組節(jié)點電壓和相角。問題的關(guān)鍵在于“給定”這兩個字。系統(tǒng)里有風(fēng)電和光伏之后注入功率不再是一個確定數(shù)值。同一座風(fēng)電場年平均風(fēng)速差一點全年發(fā)電量可能差出好幾個百分點同一天內(nèi)云層飄過光伏出力能從額定功率跌到零。此時如果你還用單點出力去做潮流分析得到的結(jié)果只對應(yīng)某一種天氣場景對系統(tǒng)規(guī)劃人員來說參考價值非常有限。更麻煩的是電壓越限這類風(fēng)險恰恰容易出現(xiàn)在極端場景里大風(fēng)天氣下風(fēng)電場滿發(fā)而本地負荷處于低谷多余功率外送導(dǎo)致局部電壓偏高傍晚光伏快速退出負荷卻還在高位電壓又會往下掉。確定性潮流算出來的“正常工況”電壓往往看不出這些邊界風(fēng)險概率潮流正是把這類風(fēng)險顯式地量化出來。1.2 概率潮流的輸入輸出形態(tài)概率潮流做的事情本質(zhì)上是輸入的隨機性到輸出的隨機性的傳遞。輸入側(cè)主要包括三類隨機變量風(fēng)速一般用兩參數(shù)Weibull分布描述概率密度函數(shù)為 f(v) (k/c)(v/c)^(k-1) exp(-(v/c)^k)其中 c 是尺度參數(shù)k 是形狀參數(shù)。光照強度通常用Beta分布描述因為光照強度在0到額定值之間連續(xù)變化Beta分布定義在有界區(qū)間上形態(tài)靈活。負荷功率一般用正態(tài)分布或?qū)?shù)正態(tài)分布波動范圍取均值的3%~8%比較常見。輸出側(cè)就是電網(wǎng)運行人員真正關(guān)心的東西節(jié)點電壓幅值的均值、標(biāo)準(zhǔn)差、概率密度函數(shù)和累計分布函數(shù)支路有功和無功的分布電壓越上/下限概率支路過載概率系統(tǒng)網(wǎng)損的期望值。這些指標(biāo)可以支撐三個層面的決策——規(guī)劃階段評估新能源接入容量是否過度激進運行階段判斷當(dāng)前方式下電壓風(fēng)險水平調(diào)度階段為備用容量和AGC調(diào)節(jié)留出合理區(qū)間。1.3 概率潮流與常規(guī)隨機分析的差別有讀者可能會問這和蒙特卡洛模擬直接撒點有什么區(qū)別蒙特卡洛只是概率潮流的實現(xiàn)手段之一。概率潮流本身是一套完整的方法論框架包含輸入隨機建模、相關(guān)性處理、不確定傳遞、輸出統(tǒng)計四個環(huán)節(jié)。蒙特卡洛是最直觀的傳遞手段點估計法和半不變量法則走解析路線用更少的計算量逼近同樣的統(tǒng)計信息。后面會把這三種方法的Matlab實現(xiàn)逐一展開。2. 先把數(shù)學(xué)模型搭起來風(fēng)光出力的隨機特性描述2.1 風(fēng)速與風(fēng)電出力模型風(fēng)速模型是整個概率潮流里最容易出錯的環(huán)節(jié)因為風(fēng)速分布和風(fēng)機出力之間還隔著一道非線性分段函數(shù)。Matlab里生成Weibull隨機數(shù)直接調(diào)用 wblrnd(scale, shape) 即可但這里有個經(jīng)典坑位Matlab的wblrnd第一個參數(shù)是尺度參數(shù)c單位m/s第二個是形狀參數(shù)k無量綱我見過不止一次有人把參數(shù)填反結(jié)果生成的風(fēng)速樣本整體偏小或偏大風(fēng)電出力分布完全失真。風(fēng)速到電功率轉(zhuǎn)換用最通用的分段線性模型P_w(v) 0, v v_in 或 v v_out P_w(v) P_rated * (v - v_in)/(v_rated - v_in), v_in v v_rated P_w(v) P_rated, v_rated v v_outv_in 一般取3m/sv_rated 取11~13m/sv_out 取25m/s。這個模型雖然簡單但在工程分析里足夠用。更精確的風(fēng)功率曲線可以用廠商實測數(shù)據(jù)插值不過概率潮流關(guān)注的是長期統(tǒng)計分布分段線性模型引入的誤差通常可以接受。Weibull參數(shù)的獲取有兩種途徑一是用當(dāng)?shù)販y風(fēng)塔一整年的小時級平均風(fēng)速數(shù)據(jù)通過極大似然估計反推二是用平均風(fēng)速和標(biāo)準(zhǔn)差近似估算。后者在參考資料匱乏的時候很實用c 約等于平均風(fēng)速的1.12倍k 則通過變異系數(shù)查表或數(shù)值求解。估算出來的參數(shù)用于方案預(yù)研足夠正式工程評估還是建議用實測數(shù)據(jù)。2.2 光照強度與光伏出力模型光伏出力的建模同樣分兩步。第一步用Beta分布描述光照強度的隨機性Beta分布的概率密度函數(shù)為 f(s) (Γ(ab))/(Γ(a)Γ(b)) * s^(a-1) * (1-s)^(b-1)Matlab里用 betarnd(a, b) 直接生成0到1之間的標(biāo)幺值。a和b的取值決定了分布的偏斜程度夏季晴天多曲線右偏a取2~3、b取1~1.5比較貼合陰雨多的地方分布相對左偏參數(shù)相應(yīng)調(diào)整。第二步做光電轉(zhuǎn)換。理想情況下光伏出力與光照強度近似線性P_pv η * S * A其中η是光電轉(zhuǎn)換效率S是實際光照強度A是光伏陣列面積。實際工程里更常用的是容量標(biāo)幺法P_pv P_rated * (S / S_ref)S_ref通常取1000W/m2。溫度對光伏出力的影響也不能完全忽略組件溫度升高導(dǎo)致輸出電壓下降更精細的模型會在上述公式基礎(chǔ)上乘一個溫度修正系數(shù)(1 - β(T - 25))β一般取0.003~0.005 /°C。對于概率潮流溫度修正要不要做得這么細取決于計算目的——如果只是評估年度電壓分布簡化模型足矣如果要做夏季高溫極端場景分析建議把溫度項加進來。2.3 負荷隨機性與相關(guān)性處理負荷波動用正態(tài)分布描述最普遍但需要加截斷處理。我在Matlab里習(xí)慣寫成PD_load PD_mean .* (1 0.05 * randn(n, 1)); PD_load(PD_load 0) 0;不做截斷的話理論上會生成負負荷樣本雖然概率極低但一旦出現(xiàn)就會導(dǎo)致潮流方程出現(xiàn)負注入結(jié)果根本沒法解釋。0.05這個系數(shù)對應(yīng)5%的標(biāo)準(zhǔn)差如果節(jié)點負荷本身波動大可以放寬到8%~10%。另一個容易被忽視的問題是相關(guān)性。風(fēng)電場和光伏電站如果處在同一個區(qū)域電網(wǎng)風(fēng)速和光照之間可能存在相關(guān)性不同節(jié)點的負荷也不會完全獨立。忽略相關(guān)性的后果是低估系統(tǒng)電壓波動的風(fēng)險區(qū)間。Matlab里處理相關(guān)高斯隨機變量最直接的方式是用Cholesky分解給定相關(guān)系數(shù)矩陣R計算 L chol(R)然后把獨立標(biāo)準(zhǔn)正態(tài)樣本矩陣 Z 乘以 L 的轉(zhuǎn)置得到帶相關(guān)性的樣本。對于Weibull和Beta這類非高斯變量更嚴格的做法是引入Copula但工程上如果只是粗略評估Cholesky分解后做概率積分變換也能用。3. 核心方法一蒙特卡洛模擬最穩(wěn)也最燒錢3.1 三步走采樣、潮流計算、統(tǒng)計蒙特卡洛模擬的思路直白到不需要過多解釋既然輸入是隨機變量那就生成大量輸入樣本逐個做確定性潮流計算最后把所有輸出結(jié)果匯集起來做統(tǒng)計分析。三個步驟對應(yīng)三段Matlab代碼每一步都可以調(diào)優(yōu)。第一步是采樣。樣本量太小分布形狀出不來樣本量太大計算耗時線性增長。第二步是潮流計算每次都調(diào)用一次確定性潮流求解器。第三步是統(tǒng)計mean、std、histogram、prctile幾個函數(shù)就夠用。蒙特卡洛最大的優(yōu)點是“無偏”——只要樣本量足夠大輸出分布一定收斂到真實分布不依賴任何線性化假設(shè)。這點在后兩種解析方法里是做不到的。代價就是計算量大。如果單次牛頓法潮流耗時0.01秒5000次就是50秒在多節(jié)點系統(tǒng)里單次潮流可能到0.1秒5000次就是8分鐘級別。3.2 Matlab主程序框架建議用Matpower做潮流計算內(nèi)核自己手寫牛頓法在原理上沒問題但處理PV-PQ節(jié)點轉(zhuǎn)換、無功越限這些細節(jié)時容易踩坑。Matpower自帶IEEE 14節(jié)點等標(biāo)準(zhǔn)算例loadcase一鍵加載runpf閉環(huán)求解省心很多。下面這套代碼是我的常用框架對應(yīng)風(fēng)電和光伏分別接入兩個不同節(jié)點的情況% 蒙特卡洛概率潮流主框架Matpower 統(tǒng)計工具箱 mpc loadcase(case14); base_PD mpc.bus(:, 3); % 保存原始有功負荷 base_QD mpc.bus(:, 4); % 保存原始無功負荷 N 5000; Vr zeros(N, 14); % 電壓幅值記錄矩陣 Pa zeros(N, 41); % 支路有功記錄矩陣case14共41條支路 Sr zeros(N, 14); % 節(jié)點注入視在功率記錄 % 風(fēng)速Weibull參數(shù)wblrnd(尺度c, 形狀k)順序不要寫反 c_w 8.5; k_w 2.2; v_in 3; v_rated 12; v_out 25; Pw_rated 30; % 風(fēng)電場額定功率 MW % 光照Beta分布參數(shù) a_s 2; b_s 1.5; Pp_rated 15; % 光伏電站額定功率 MW for i 1:N % 1. 采樣風(fēng)速 - 風(fēng)電出力 v wblrnd(c_w, k_w); if v v_in || v v_out Pw 0; elseif v v_rated Pw Pw_rated * (v - v_in) / (v_rated - v_in); else Pw Pw_rated; end % 2. 采樣光照強度 - 光伏出力 s_pu betarnd(a_s, b_s); Pp Pp_rated * s_pu; % 3. 把風(fēng)光出力等效為相應(yīng)節(jié)點的負負荷注入 PD_node base_PD; PD_node(9) base_PD(9) - Pw / 100; % 風(fēng)電接入bus9 PD_node(13) base_PD(13) - Pp / 100; % 光伏接入bus13 % 4. 負荷波動對包含負注入的凈負荷施加正態(tài)擾動 PD_load PD_node .* (1 0.05 * randn(14, 1)); QD_load base_QD .* (1 0.05 * randn(14, 1)); PD_load(PD_load 0) 0; QD_load(QD_load 0) 0; mpc.bus(:, 3) PD_load; mpc.bus(:, 4) QD_load; % 5. 確定性潮流計算 res runpf(mpc, mpoption(out.all, 0)); if res.success 0 warning(第 %d 次潮流不收斂, i); continue; end % 6. 記錄輸出量 Vr(i, :) res.bus(:, 8); % 電壓幅值 Pa(i, :) res.branch(:, 14); % 支路有功 end注意代碼第3步和第4步的順序。先把風(fēng)電和光伏出力折算成負的凈負荷再對這種凈負荷施加正態(tài)擾動等價于“風(fēng)光出力隨機 負荷隨機”的疊加方式。如果你先把原負荷擾動完再把風(fēng)光注入單獨減掉那么在注入很大的節(jié)點上凈負荷可能出現(xiàn)負值且分布形態(tài)被扭歪。另外代碼里 Pw/100 是因為case14的基準(zhǔn)容量是100MVA把所有功率統(tǒng)一折算到標(biāo)幺值這一點新手特別容易漏。3.3 收斂性判斷與采樣規(guī)模選擇到底采多少組樣本才算夠經(jīng)驗法則是先看輸出量的均值或標(biāo)準(zhǔn)差隨樣本數(shù)的變化曲線當(dāng)相對波動小于某個閾值時認為收斂。我常用的判據(jù)是電壓均值的無窮范數(shù)誤差連續(xù)兩次采樣之間的變化量小于1e-4就停止。更直接的做法是固定樣本量5000次起步不好再翻倍到10000次。N_max 10000; V_mean_old zeros(1, 14); for i 1:N_max % …… 上述采樣與潮流計算代碼 …… V_mean_new mean(Vr(1:i, :), 1); delta max(abs(V_mean_new - V_mean_old)); if delta 1e-4 i 1000 fprintf(均值收斂于第 %d 次采樣\n, i); break; end V_mean_old V_mean_new; end這個小循環(huán)跑起來有個好處你不用賭樣本量機器自己告訴你夠了。代價是循環(huán)內(nèi)多了mean運算對整體耗時影響很小。蒙特卡洛還天然支持并行化把for改成parfor前提是循環(huán)體內(nèi)不能有依賴全局變量的操作runpf的mpc結(jié)構(gòu)每次都是基于本次采樣數(shù)據(jù)構(gòu)造的滿足parfor要求。我實測過在8核機器上3000次仿真的耗時能壓到原來的三分之一左右。注意parfor里mpc這個變量會被當(dāng)作廣播變量處理數(shù)據(jù)量不大影響有限。4. 省時方案點估計法與半不變量法的Matlab實現(xiàn)4.1 點估計法用少量確定性潮流逼近統(tǒng)計量點估計法的基本思想很取巧輸入隨機變量的分布不參與顯式采樣而是用輸入變量的前幾階矩均值、方差、偏度等構(gòu)造出若干個確定性估計點和對應(yīng)權(quán)重對每個估計點做確定性潮流再對輸出加權(quán)求和得到統(tǒng)計量。最基礎(chǔ)的2點估計法規(guī)則如下對每個輸入隨機變量 x_i取兩個估計點x_{i,1} μ_i σ_i x_{i,2} μ_i - σ_i權(quán)重各取 1/2。然后把第 i 個輸入變量固定在這兩個點上其他輸入變量固定在均值處分別做兩次確定性潮流。如果系統(tǒng)里有 n 個隨機輸入變量總共需要 2n 次潮流計算。相比蒙特卡洛動輒幾千次計算量是天壤之別。輸出變量的期望和方差按下面公式聚合E[Y] ≈ Σ_i Σ_k w_{i,k} * Y(x_{i,k}) E[Y^2] ≈ Σ_i Σ_k w_{i,k} * Y^2(x_{i,k}) Var[Y] E[Y^2] - (E[Y])^22點估計只用到均值和方差對線性系統(tǒng)是精確的對非線性系統(tǒng)會有截斷誤差。想要更高精度可以用3點估計額外引入偏度信息ξ_{i,1} λ3/2 sqrt(λ4 - 3λ3^2/4) ξ_{i,2} λ3/2 - sqrt(λ4 - 3λ3^2/4) ξ_{i,3} 0對應(yīng)的估計點為 x_{i,k} μ_i ξ_{i,k} * σ_i。3點估計的權(quán)重計算鏈條稍長我在Matlab里建議直接用相關(guān)工具箱或者核對Hong在1999年原始文獻的公式避免抄錯。點估計法最大的短板是它只能給出輸出的均值、方差等低階矩不能直接恢復(fù)完整的概率密度分布。如果評估報告里必須畫電壓概率密度曲線點估計法就幫不上忙了。4.2 半不變量法 Gram-Charlier級數(shù)半不變量法走的是“矩-半不變量-級數(shù)展開”的分析路線。先利用輸入隨機變量的概率分布求出各階半不變量然后在期望運行點做一次確定性潮流得到靈敏度矩陣再把輸入半不變量線性映射到輸出最后用Gram-Charlier或Edgeworth級數(shù)擬合輸出分布。Matlab實現(xiàn)的核心步驟是這樣第一步由輸入隨機變量的各階矩計算半不變量。前四階半不變量與矩的關(guān)系為κ1 μ1 κ2 μ2 - μ1^2 κ3 μ3 - 3μ1μ2 2μ1^3 κ4 μ4 - 4μ1μ3 6μ1^2μ2 - 3μ1^4第二步在系統(tǒng)期望運行點做一次牛頓法潮流取得雅可比矩陣求逆得到靈敏度矩陣 S0 J^(-1)。第三步線性映射。對于第i個輸出量Y_i輸入隨機變量第s階半不變量的貢獻為 (S0(i,j))^s 乘以輸入的第s階半不變量再對j求和。公式為 κ_Y,s Σ_j (S0(i,j))^s * κ_W,s。第四步用Gram-Charlier級數(shù)把輸出半不變量轉(zhuǎn)換成概率密度。令 z (Y - μ_Y) / σ_Y則f(z) φ(z) * [1 (κ3/6σ^3) * He3(z) (κ4/24σ^4) * He4(z) ...]其中 φ(z) 是標(biāo)準(zhǔn)正態(tài)密度函數(shù)He3(z)z^3-3zHe4(z)z^4-6z^23。這套方法的計算速度最快適合在線評估場景。代價是靈敏度矩陣來自潮流方程的線性化對非線性強、重尾分布明顯的系統(tǒng)展開到四階筋的精度改善有限偶爾會出現(xiàn)概率密度曲線局部負值的現(xiàn)象這是級數(shù)截斷本身帶來的問題。4.3 三種方法怎么選方法精度計算量實現(xiàn)難度能否重建分布適用場景蒙特卡洛模擬最高無偏極大數(shù)千次潮流低能完整直方圖標(biāo)準(zhǔn)分析、驗證其他方法點估計法中等低階矩精度高小2~3n次潮流中不能只有矩快速評估均值/標(biāo)準(zhǔn)差半不變量法線性化精度極小1次潮流映射高能近似解析分布在線評估、海量場景遍歷從我實際使用的感受來說學(xué)術(shù)論文里最穩(wěn)妥的套路是“蒙特卡洛做基準(zhǔn)、點估計或半不變量做改進方法”。先用蒙特卡洛給出精確結(jié)果再展示改進方法在誤差和耗時上的對比。如果直接上來就做半不變量法而缺失基準(zhǔn)驗證審稿人大概率會追問一句“和蒙特卡洛對比過嗎”。5. 案例實操IEEE 14節(jié)點系統(tǒng)接入風(fēng)光電源5.1 系統(tǒng)改造與參數(shù)設(shè)定這次演示以Matpower自帶的case14為基礎(chǔ)。IEEE 14節(jié)點系統(tǒng)有14個節(jié)點、5臺發(fā)電機系統(tǒng)基準(zhǔn)容量100MVA。我在原始算例基礎(chǔ)上做了三處改動風(fēng)電場接入節(jié)點9額定功率30MW。節(jié)點9原本是純負荷節(jié)點用它接入風(fēng)電后不需要改變發(fā)電機配置直接把注入功率折算成負負荷就行。光伏電站接入節(jié)點13額定功率15MW。同樣處理為負負荷。所有負荷施加5%標(biāo)準(zhǔn)差的正態(tài)擾動截斷到非負。風(fēng)速Weibull參數(shù)取 c8.5m/s、k2.2切入風(fēng)速3m/s、額定風(fēng)速12m/s、切出風(fēng)速25m/s。光照Beta分布參數(shù)取 a2、b1.5。這些參數(shù)偏理想化但演示概率潮流的完整流程足夠了。如果要在實際工程中使用參數(shù)務(wù)必換成現(xiàn)場實測數(shù)據(jù)。5.2 完整代碼實現(xiàn)與運行說明完整代碼在第3章的框架基礎(chǔ)上增加收斂判斷、結(jié)果統(tǒng)計和可視化三個環(huán)節(jié)。我直接貼出循環(huán)結(jié)束后的統(tǒng)計部分% 剔除不收斂樣本假設(shè)存于Vr中不收斂行全為0 Vr_valid Vr(all(Vr 1e-8, 2), :); Pa_valid Pa(all(Vr 1e-8, 2), :); % 節(jié)點電壓統(tǒng)計指標(biāo) V_mean mean(Vr_valid, 1); V_std std(Vr_valid, 1); V_p5 prctile(Vr_valid, 5, 1); V_p95 prctile(Vr_valid, 95, 1); % 示例節(jié)點4的電壓越限概率 prob_low mean(Vr_valid(:, 4) 0.95); prob_high mean(Vr_valid(:, 4) 1.05); fprintf(節(jié)點4電壓均值 %.4f p.u.標(biāo)準(zhǔn)差 %.4f p.u.\n, V_mean(4), V_std(4)); fprintf(電壓低于0.95概率%.4f%%高于1.05概率%.4f%%\n, prob_low*100, prob_high*100); % 支路過載概率有功超過線路容量1.0p.u.基準(zhǔn)100MVA overload_prob mean(max(Pa_valid, [], 1) 1.0); fprintf(支路過載概率%.4f%%\n, overload_prob*100); % 繪制節(jié)點4電壓幅值分布 figure; histogram(Vr_valid(:, 4), 80, Normalization, pdf); xlabel(節(jié)點4電壓幅值 (p.u.)); ylabel(概率密度); title(節(jié)點4電壓幅值概率分布5000次蒙特卡洛);運行這段代碼需要提前確認Matlab環(huán)境具備了統(tǒng)計工具箱wblrnd、betarnd、histogram這些函數(shù)都依賴它和Matpower工具箱。Matlab版本我試過R2021b和R2023a都能跑通新版本沒有遇到兼容性問題。如果你不想裝Matpower也可以自己寫牛頓法潮流函數(shù)但需要注意幾個細節(jié)PV節(jié)點無功越限時要轉(zhuǎn)換成PQ節(jié)點重新迭代平衡節(jié)點的相角要固定雅可比矩陣稀疏化用sparse構(gòu)造不要用滿陣否則系統(tǒng)規(guī)模一大內(nèi)存直接爆掉。5.3 結(jié)果怎么看分布形態(tài)、越限概率與確定性解的差異我這次演示跑出來的典型結(jié)果大致是這樣參數(shù)不同結(jié)果會有波動重點看分布形態(tài)節(jié)點4是系統(tǒng)中比較靠近負荷中心的節(jié)點電壓均值大約在1.01p.u.標(biāo)準(zhǔn)差在0.012p.u.量級。這看起來波動幅度不大但分布尾部確實會越出 [0.95, 1.05] 的常規(guī)運行區(qū)間。支路過載概率非常低在千分位以下這符合case14網(wǎng)架結(jié)構(gòu)相對堅強的特點。一個值得注意的現(xiàn)象是蒙特卡洛采樣得到的電壓均值往往不等于把所有隨機變量固定在期望值時做確定性潮流得到的電壓值。原因是潮流方程關(guān)于注入功率是高度非線性的電壓幅值對注入的響應(yīng)帶有凸性期望值變換到了非線性函數(shù)內(nèi)部就不再等價。這也是概率潮流區(qū)別于“把期望值代入確定性潮流”的根本原因。如果你在報告里寫“風(fēng)光出力取期望潮流算一遍結(jié)果即為系統(tǒng)平均運行狀態(tài)”這在數(shù)學(xué)上是站不住腳的。審稿時這個問題是高頻質(zhì)疑點。6. 常見問題與排查技巧實錄6.1 潮流不收斂怎么辦蒙特卡洛循環(huán)里最煩人的就是跑著跑著某一次潮流不收斂。先用if res.success 0 continue把不收斂樣本剔掉保證主程序不中斷然后回過頭排查不收斂的原因。我從實際調(diào)試經(jīng)驗看排在前面的原因有三個一是風(fēng)光注入功率太大。當(dāng)節(jié)點凈負荷為負且數(shù)值很大時相當(dāng)于一個功率倒送的發(fā)電機節(jié)點潮流方程可能走上一條不收斂的迭代路徑。解決辦法是檢查注入功率是否超過系統(tǒng)承受能力適當(dāng)降低風(fēng)電場額定容量或者給該節(jié)點增加無功補償設(shè)備。二是有功注入過大導(dǎo)致電壓偏高觸發(fā)發(fā)電機無功越限PV-PQ轉(zhuǎn)換反復(fù)震蕩。這種情況可以在潮流計算中打開無功越限處理選項或者調(diào)整該節(jié)點的無功補償容量。三是采樣到了極端惡化的負荷組合。當(dāng)多個節(jié)點負荷同時處于波動上界時系統(tǒng)運行點可能逼近電壓穩(wěn)定邊界。這時需要回溯樣本參數(shù)看看是不是概率分布參數(shù)定得太激進。6.2 計算太慢怎么優(yōu)化蒙特卡洛的耗時大頭在重復(fù)潮流計算。同樣的網(wǎng)絡(luò)導(dǎo)納矩陣結(jié)構(gòu)每次都一樣但Matpower每輪都會重新生成和分解。優(yōu)化手段按收益排序優(yōu)先把mpoption(out.all, 0)設(shè)上關(guān)閉MATPOWER的屏幕輸出5000次仿真能省掉大約20%的IO時間。用parfor替代for這是最直接的提速手段。需要注意parfor里所有變量都必須符合切片規(guī)則我習(xí)慣把每次循環(huán)需要的數(shù)據(jù)預(yù)先構(gòu)造成矩陣循環(huán)內(nèi)只做索引切片。如果自己寫牛頓法可以把雅可比矩陣中與網(wǎng)絡(luò)拓撲相關(guān)的常數(shù)部分離線算好每次迭代只更新與節(jié)點注入相關(guān)的局部元素。這個方法能壓掉不少時間但對代碼能力有一定要求前期不建議折騰。對問題規(guī)模大、采樣次數(shù)要求高的場景考慮用點估計法替代蒙特卡洛做快速預(yù)篩再用蒙特卡洛對高風(fēng)險場景重點驗算。6.3 結(jié)果異常的排查方向碰到概率分布形狀詭異、均值偏移明顯這類問題我通常按下面這個速查表逐項排查現(xiàn)象常見原因排查與解決電壓均值明顯偏低或偏高輸入隨機變量均值參數(shù)與基準(zhǔn)工況不一致先將所有隨機變量固定為期望值跑確定性潮流與Matpower基準(zhǔn)結(jié)果比對風(fēng)電功率樣本出現(xiàn)負值Weibull函數(shù)參數(shù)順序填反檢查wblrnd調(diào)用正確形式為wblrnd(尺度c, 形狀k)分布直方圖出現(xiàn)雙峰負荷截斷過狠導(dǎo)致樣本集中在零附近或風(fēng)光參數(shù)組合形成多模態(tài)檢查輸入樣本直方圖單獨繪制風(fēng)速、光照分布形態(tài)概率密度曲線局部負值半不變量法級數(shù)截斷造成的振蕩增加展開階數(shù)或改用Edgeworth級數(shù)或直接換蒙特卡洛復(fù)核潮流反復(fù)不收斂且集中在特定樣本段該區(qū)間對應(yīng)高滲透率極端場景檢查該樣本的風(fēng)速、光照組合值評估是否超出系統(tǒng)靜態(tài)穩(wěn)定約束蒙特卡洛與點估計法的方差結(jié)果差異大系統(tǒng)非線性強低階矩方法截斷誤差放大以蒙特卡洛為準(zhǔn)增加點估計法的估計點數(shù)核驗還有一個容易忽略的細節(jié)如果風(fēng)光接入節(jié)點原本帶負荷用負負荷等效后負荷波動生成器會對凈負荷做擾動此時同一節(jié)點的注入波動和負荷波動被混在一起。嚴格來說這部分相關(guān)性在概率模型中并未分離。想處理干凈就把“基礎(chǔ)負荷”和“新能源注入”作為兩個獨立的隨機源分開采樣后在節(jié)點注入方程中相加代碼里要預(yù)留對應(yīng)的接口。7. 我踩過的坑和一點個人體會第一次跑通蒙特卡洛概率潮流的時候我用的還是純手寫的牛頓法潮流5000次仿真跑了將近十分鐘Matlab風(fēng)扇嗡嗡轉(zhuǎn)結(jié)果電壓均值比Matpower基準(zhǔn)低了將近2%查了半天才發(fā)現(xiàn)是Weibull尺度參數(shù)c和形狀參數(shù)k填反了。從那以后我養(yǎng)成了一個習(xí)慣任何隨機分布參數(shù)進循環(huán)之前先單獨生成一組樣本畫直方圖目測形態(tài)是否合理。分布參數(shù)錯了后面一切結(jié)果都是空中樓閣這一步省不得。另一個體會是方法論選型的順序。我建議初學(xué)者不要一上來就鉆研半不變量法和Gram-Charlier級數(shù)先用蒙特卡洛把“輸入隨機到輸出隨機”的直覺建立起來看懂電壓分布是怎么來的再去研究怎么用更少的計算量逼近它。順序反了的話公式推了一堆結(jié)果出了偏差你都不知道該懷疑是哪一步。如果后續(xù)想把這套東西擴展到工程應(yīng)用兩個大方向可以考慮一是把風(fēng)光出力之間的空間相關(guān)性特別是同一氣候區(qū)內(nèi)多個風(fēng)電場之間的出力相關(guān)性建進去否則風(fēng)險評估會偏樂觀二是結(jié)合時序運行模擬把風(fēng)光出力的時間相關(guān)性考慮進來這樣得到的電壓越限概率才真正對應(yīng)實際運行中持續(xù)時間的累積風(fēng)險。概率潮流本身解決的是“截面不確定性”問題要和時序信息結(jié)合才完整覆蓋新能源并網(wǎng)評估的整個拼圖。