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