最近在搞电力系统不确定性分析把蒙特卡洛法用在概率潮流计算里折腾了一阵子。说实话一开始我心里是有点抗拒的——这年头机器学习、深度学习满天飞还用这种几十年前的老派“暴力采样”方法总觉得不够高级。但真拿IEEE33节点电网当小白鼠跑完整个流程之后我得说蒙特卡洛在概率潮流计算里的实用价值还真不是花架子。它那种“以量取胜”的朴素逻辑反而是应对风光出力随机性最皮实、最好用的武器尤其是处理那些解析法搞不定的强非线性和复杂分布时优势特别明显。这篇文章我直接摊开讲把蒙特卡洛概率潮流计算的完整流程拆一遍从为什么确定性潮流不够用到风光出力模型怎么搭再到IEEE33节点电网上的代码实现和结果解读全部干货放在台面上。代码我用的Python pandapower都是电力系统分析里最常见的工具你照着跑就能复现。适合做配电网规划、分布式电源接入评估或者正在研究电力系统不确定性方向的朋友看完应该能少走不少弯路。1. 从确定性到概率为什么潮流计算必须考虑不确定性1.1 确定性潮流的盲区传统潮流计算有一个根深蒂固的假设输入是确定的。给定负荷大小、给定发电机出力解一个非线性方程组出来一组节点电压和支路潮流的定值。电网规划环评也好、日常调度也罢长期就是这么用的。但现实里哪来那么多“给定”负荷一直在波动居民用电高峰和低谷能差出两三倍这谁都知道。更麻烦的是风光接入之后不确定性是成倍放大的。一片云飘过来光伏出力可能在几分钟内从满发掉到零一阵妖风过去风机直接停机。你拿一个固定的出力值去算潮流相当于闭着眼睛开车——结果可能在某一刻是准的但绝大多数时间你根本不知道系统到底处于什么状态。确定性潮流给出的是一个点而工程上真正需要知道的是这个点的邻域有多大、边界在哪里。这就是为什么要把确定性潮流扩展成概率潮流。1.2 不确定性分析方法对比为什么选蒙特卡洛概率潮流的本质是把输入随机变量的概率分布通过潮流计算这层非线性变换映射到输出变量电压、支路功率等的概率分布上。学术界折腾出很多方法大体分三类解析法、近似法、模拟法。方法类型代表方法优点缺点解析法Gram-Charlier级数、Cornish-Fisher展开计算快能得到解析表达式对强非线性、多随机变量组合时误差大推导繁琐近似法点估计法、一次二阶矩法计算量小适合快速评估需要假设变量矩存在遇到风光出力的分段函数容易失真模拟法蒙特卡洛法及其变体精度可控实现简单对非线性免疫收敛慢需要的样本量大蒙特卡洛的打法特别直白既然解析上搞不定概率分布是怎么通过非线性方程传播的那就多抽点样本一个一个算算完了做统计。大数定律保证样本均值收敛于期望值频率收敛于概率。这就类比成你想知道一个池塘里鱼的平均重量不用把池塘抽干随机捞三千条上来称一称平均值就和全局均值差不太远了。捞得越多越接近真值。在IEEE33这种规模的配电网里节点数不多单次潮流计算耗时通常只有几毫秒。蒙特卡洛最让人诟病的“计算量大”问题在这里其实被彻底稀释了。三千次潮流算下来也就几十秒完全在可接受范围内。对于要搞强非线性、要处理威布尔分布和Beta分布这些非高斯输入的场景蒙特卡洛的鲁棒性碾压那些花哨的解析法。1.3 概率潮流到底能回答什么问题概率潮流给的不是一堆求不出来的公式而是几个工程上特别想要答案的问题。第一越限概率。某个节点电压低于0.95 pu的概率是多少某条支路过载的概率是多少这在确定性潮流里你只能答“这个工况没越限”但换成概率语言就是“根据我当地的风光资源和负荷特性一年里大概有百分之几的时间会越限”。第二风险识别。在所有样本里哪些节点最容易出事哪些支路承受的压力波动最大这能直接用来指导无功补偿配置和线路扩容优先级比单纯看均值的说服力强得多。第三方案评估。如果我在节点32接入一个光伏电站对末端电压抬升的效果有多强接入点和接入容量的不同组合会带来完全不同的概率分布形态。蒙特卡洛可以帮你把每个候选方案的概率分布都画出来对比。我个人踩过最多的坑就是把概率潮流的结果简单等同于“多算几次潮流取平均”。真不是这样的均值只是初始信息最值、分位数、方差、偏度和越限概率每个维度都有自己的工程含义后面到了结果分析部分我再展开说。2. 搭建实验平台IEEE33节点系统与风光出力模型2.1 IEEE33节点测试系统快速上手IEEE33节点配电网是业内研究分布式电源接入、配电网重构时最经典的小规模测试系统结构是一条馈线带多个分支的辐射状网络总共33条母线、32条支路基准电压12.66 kV总负荷大约3.7 MW加2.3 Mvar。它的经典之处在于保留了真实配电网的主要特点——末端电压低、分支多、线路阻抗R/X比值高而且规模小、算得快非常适合当蒙特卡洛的“小白鼠”。用pandapower加载这个系统只需要两行代码import pandapower as pp import pandapower.networks as pn # 载入IEEE33节点配电网测试系统 net pn.case33bw()跑一个基础潮流看看系统“出厂状态”pp.runpp(net, algorithmnr) print(net.res_bus.vm_pu.min())不接任何分布式电源的原始算例下末端节点约对应节点18和节点32的电压很低通常会在0.93 pu左右已经逼近0.95 pu的低电压警示线了。这就是IEEE33这个系统的“出厂设定”——底子差所以很适合用来验证分布式电源接入的改善效果。注意pandapower不同版本对case33bw的调用方式可能有细微差别。如果你使用的版本里pn.case33bw()报错检查一下网络加载函数名有些版本是create_case33bw或者直接看pandapower.networks模块里有哪些IEEE33相关函数按实际调整即可。2.2 风电出力模型从风速到功率风机出力建模第一步是先描述风速这个随机变量。风速有很多统计模型最常用的是两参数威布尔分布。它的概率密度函数长这样f(v) (k/c) * (v/c)^(k-1) * exp(-(v/c)^k)其中k是形状参数c是尺度参数。k2时退化为瑞利分布实际拟合某个风电场的历史风速数据时k一般在1.5到3之间c大概对应平均风速的量级。我代码里默认取k2、c8这组参数模拟的风速期望值大约在7 m/s左右属于中规中矩的内陆风场水平。风速到功率的转化用的是分段线性模型里面藏着风机的三个关键风速点切入风速一般3 m/s。风速太低时叶片转不起来发电机无法并网出力为零。额定风速一般12到14 m/s。达到额定风速后风机满发。切出风速一般25 m/s。风速太高时必须停机保护出力直接跳回零。功率曲线本质上是一根分段折线低于切入风速出力为0切入与额定之间近似线性爬坡额定到切出之间平顶超过切出又回0。这个分段函数的非线性特征非常明显也正是解析法概率潮流容易在处理这根折线时翻车的原因。2.3 光伏出力模型从辐照度到功率光伏的随机性来源主要是辐照度业内常用Beta分布描述这个分布定义在[0, 1]区间乘以当地最大辐照度就得到实际辐照度值。为什么用Beta分布因为它形状灵活可以拟合从“大晴天”到“多云天”的不同辐照度分布形态。比如alpha2、beta5这组参数分布偏向右侧低辐照度概率高对阴天较多的地区比较贴切alpha5、beta2则反映出日照好的地方那种偏向左侧、高辐照度概率高的特征。辐照度到光伏出力的转化相对直接P_pv P_rated * (G / G_ref)P_rated是光伏额定容量G是实际辐照度G_ref是标准测试条件辐照度取1000 W/m²。这个公式其实忽略了温度对光伏效率的影响属于工程上常用的简化处理。真要精细建模还得叠加一个温度系数修正项温度每升高一度转换效率下降千分之几。但蒙特卡洛本来就是要抽样几千次模型每复杂一点计算量就蹭蹭涨所以第一版先跑简化模型完全没问题后面有必要再升级。实操心得光伏出力模型里还有个容易忽略的点——夜间辐照度是零Beta分布采样永远不可能采出“零”这个值但真实光伏夜里就是零出力。处理方式很简单给辐照度加一个概率质量在零处或者像我这样以固定比例让样本落在“夜间/极端阴天”场景用if判断直接置出力为0。蒙特卡洛的优势就在于这种非连续场景可以轻松处理。2.4 负荷随机性建模风光之外负荷本身也是大随机源。常用简化假设是每个负荷节点的有功功率服从正态分布以原始负荷为均值标准差取均值的5%左右。正态分布的好处是采样方便物理上也说得过去——大量用户独立用电行为的叠加在大数定律下确实趋近正态。更精细的建模还要考虑各节点负荷之间的相关性。同一个小区的两个负荷点用电模式肯定比相距很远的两个负荷点更像。但第一版蒙特卡洛不必急着加相关性先把独立波动的baseline跑出来后面需要再引入相关性矩阵做Cholesky分解。工程迭代要一步一步来一上来就搞高精度模型出了bug都分不清是采样问题还是模型问题。我把这三个分布式电源设置在节点17、22、32其中节点17接风电额定容量0.5 MW节点22和节点32各接光伏额定容量分别0.3 MW和0.2 MW。三个DG合计最大出力1 MW对系统峰值负荷3.7 MW来说渗透率接近三成在配电网研究里属于比较典型的接入场景既能看出DG的支撑效果又不至于把潮流算崩。3. 蒙特卡洛概率潮流计算完整代码实现3.1 代码架构总览整个蒙特卡洛概率潮流程序的架构可以用一句话概括采样、折算、注入、算潮流、存结果、做统计。代码主体是一个for循环循环次数就是采样规模。每轮循环做完以下事情采样一组风速和辐照度折算成风力和光伏出力把出力写入网络中的发电机模型给负荷叠加上随机波动然后解一次潮流把电压、支路潮流、网损结果存起来。循环结束之后对所有结果做统计分析输出概率指标。在写代码之前先把参数和随机数种子设置好import numpy as np import pandapower as pp import pandapower.networks as pn import matplotlib.pyplot as plt from scipy.stats import gaussian_kde np.random.seed(42) # 固定种子保证结果可复现 # 蒙特卡洛参数 n_samples 3000 # 威布尔分布参数风速 k_shape, c_scale 2.0, 8.0 # Beta分布参数辐照度 alpha_irr, beta_irr 2.0, 5.0 g_max 1000 # 最大辐照度 W/m²3.2 风光出力折算函数风速转功率和辐照度转功率的逻辑建议写成独立函数。这不仅是代码洁癖更是工程习惯——后面的优化、替换模型都只需要动函数内部不影响主循环。def wind_power(v, p_rated0.5, v_ci3.0, v_r12.0, v_co25.0): 风速 m/s - 风机出力 MW if v v_ci or v v_co: return 0.0 elif v v_r: return p_rated * (v - v_ci) / (v_r - v_ci) else: return p_rated def pv_power(g, p_rated0.3, g_ref1000.0): 辐照度 W/m² - 光伏出力 MW忽略温度影响的简化模型 return p_rated * max(g, 0.0) / g_ref注意一个细节风速介于切入风速和额定风速之间时我用了线性内插。更精确的风机功率曲线往往是一根S形曲线某些厂商还提供三次方曲线模型。但从蒙特卡洛总体的角度分段线性和S形曲线在输出分布上的差异并不显著影响远小于风速分布参数选取带来的差异。所以第一版用线性就好省事且结果不至于失真。3.3 蒙特卡洛主循环实现主循环是整个程序的核心每一轮都包含五个步骤。我把完整代码贴出来注意里面几个关键写法# 先保存基础负荷后面每轮在基础值上叠波动 base_load_p net.load.p_mw.values.copy() base_load_q net.load.q_mvar.values.copy() # 接入3个分布式电源先全部设0出力 pp.create_sgen(net, 17, p_mw0.0, nameWind) pp.create_sgen(net, 22, p_mw0.0, namePV1) pp.create_sgen(net, 32, p_mw0.0, namePV2) # 结果存储 V_samples np.zeros((n_samples, len(net.bus))) loading_samples np.zeros((n_samples, len(net.line))) loss_samples np.zeros(n_samples) n_failed 0 for i in range(n_samples): # 1. 采样风速和辐照度 v_wind np.random.weibull(k_shape) * c_scale g_irrad np.random.beta(alpha_irr, beta_irr) * g_max # 2. 折算风光出力 p_w wind_power(v_wind, p_rated0.5) p_pv1 pv_power(g_irrad, p_rated0.3) p_pv2 pv_power(g_irrad, p_rated0.2) # 3. 更新DG出力 net.sgen.loc[0, p_mw] p_w net.sgen.loc[1, p_mw] p_pv1 net.sgen.loc[2, p_mw] p_pv2 # 4. 叠加负荷波动5%标准差的正态分布扰动 load_scale np.random.normal(1.0, 0.05, sizelen(net.load)) net.load.loc[:, p_mw] base_load_p * load_scale net.load.loc[:, q_mvar] base_load_q * load_scale # 5. 求解潮流并保存结果 try: pp.runpp(net, algorithmbfw) V_samples[i, :] net.res_bus.vm_pu.values loading_samples[i, :] net.res_line.loading_percent.values loss_samples[i] net.res_line.pl_mw.sum() except: n_failed 1这段代码里有个容易踩的坑pandapower的sgen和load数据框直接用loc按索引赋值是安全的。但如果你在一个仿真里动态增加或删除元素索引就会错位。我这个程序里DG数量和顺序在整个循环里保持不变所以直接loc赋值没问题。另外我在每轮循环都重新赋值net.load.loc[:, p_mw]很多新手会在这儿犯迷糊以为这是在“覆盖”原始网络结构。其实不是这是在给负荷变量赋新值没有破坏网络拓扑定义。运行完毕之后如果还想回到原始网络重新做其他仿真别忘了恢复基础负荷或者干脆把网络重新加载一遍。3.4 为什么算法选bfw而不是nr这里要展开说一下算法选择。配电网和输电网的电气特性差异很大典型配电网线路R/X比值高、节点多、分支密传统的牛顿-拉夫逊法在初始化不好或负荷过重时容易震荡甚至不收敛。pandapower默认的nr算法在IEEE33这种小规模系统上通常能收敛但加上高比例的分布式电源后不收敛的风险会上升。pandapower提供了面向配电网的前推回代法参数写algorithmbfw。这种算法本质上走的是配电网天然的辐射状结构从根节点前推电流、回代电压迭代几次就能收敛速度快、稳定性强在IEEE33这种树上格外合适。我的建议是默认优先尝试bfw如果你是要和别的算法结果互验再用nr跑同一个算例对比。概率潮流一个循环要算几千次潮流算法收敛性的微小差异都会被放大选一个稳的算法能省很多头疼事。实操心得主循环里加了try/except之后我建议顺手打印不收敛的次数。如果n_failed占比超过0.1%大概率是参数设置出了问题比如DG容量过大、负荷波动标准差设得太极端而不是算法的偶发失败。这时候要先回头查输入参数别急着继续加样本量。3.5 样本量怎么定蒙特卡洛的经典问题是样本量取多少才够理论上样本量越大越逼近真分布但计算代价也线性增长。实际做法是“先算后验”跑一组小样本比如500次统计输出变量的均值和方差看它们随样本量的变化是否稳定。我这里的经验是对IEEE33这种小系统3000个样本足够让节点电压均值的变异系数标准差/均值降到0.1%以内。如果你想输出更极端的尾部概率比如电压低于0.85 pu这种小概率事件样本量还要再加大因为概率越小的尾巴需要越多样本才能“采到”。一个粗略的经验公式是要估计p量级的概率样本量最好不低于1/p的10倍。想看到1%概率的越限事件至少抽1万次才稳。3.6 收敛性统计判据如果你不想拍脑袋定样本量可以用一个循环动态判断每跑完一批样本计算当前电压均值与上一批均值之差如果相对误差小于0.5%就停。这属于“在线式”蒙特卡洛工程里更常用但实现上稍微复杂一点。我先保证逻辑简单用固定样本量的版本跑通你后面要精进再改在线式判断。4. 算例结果解读从概率视角重新审视电网状态4.1 电压分布一条分布带而非一根线跑完3000次潮流后最直观的输出就是33个节点电压的统计结果。我先把均值、标准差、0.05分位数和0.95分位数统计出来V_mean V_samples.mean(axis0) V_std V_samples.std(axis0) V_p5 np.percentile(V_samples, 5, axis0) V_p95 np.percentile(V_samples, 95, axis0)画图的话横轴是节点编号纵轴是电压幅值pu把V_p5和V_p95画成两条包络线中间就是90%概率区间。你会发现这条包络带的宽度在末端节点明显变宽——为什么因为DG接入在节点17、22、32附近对近端电压有支撑作用但末端节点离电源远、线路阻抗压降大风光出力波动通过线路传播到末端时不确定性被放大了。这就像平静河流的上游扔了块石头下游的水面波纹不会消失只会扩散铺开。具体到我这组参数不接DG时节点18的电压是固定的约0.93 pu。接入风光后节点18电压变成了一个均值为0.95、标准差约0.011的分布。看起来均值改善了不少但真正要命的还是那个5%分位数——它可能还在0.93以下。这意味着虽然“平均情况”达标了但总有约5%的场景电压还是偏低的。这就是概率潮流在说人话别只盯着均值尾部风险才是决策关键。4.2 电压越限概率风险量化我习惯把越限概率单独画一张柱状图。对IEEE33网络电压运行范围通常按0.95~1.05 pu考核。用代码算一下prob_low (V_samples 0.95).mean(axis0) prob_high (V_samples 1.05).mean(axis0)结果往往会让你吓一跳。某个末端节点的低压越限概率可能高达15%而某个DG接入点附近由于光伏大发时局部电压抬升高压越限概率也会冒出来几个百分点。低压风险来自重负荷风光不出力高压风险来自轻负荷风光满发两种场景可能概率都不高但都真实存在确定性潮流却只能告诉你“当前这个场景没越限”。工程上还有个更细的问题极限场景的持续时间。比如电压低于0.9 pu的概率虽然是1%但每次持续时间是几分钟还是几小时影响完全不同。这个要靠时序蒙特卡洛解决在本文的静态蒙特卡洛里暂时不考虑但你心里要清楚这个边界。4.3 支路潮流与网损分布支路潮流同理。loading_samples里存的是每条支路负载率的百分比超过100%就是过载。统计后你会看到IEEE33网络上靠近电源的那几条主干支路负载率波动幅度最大因为所有DG出力和负荷波动的净功率都要从这些支路流过。它们的过载概率可能并不高但95%分位数负载率可能已经到了85%以上线路扩容需求其实已经提前暴露了。网损分布更有意思。loss_samples记录了每轮潮流下全网有功损耗的数值。我这次跑出来的网损均值大约在0.12 MW标准差0.02 MW概率分布呈右偏态——少数极端场景下网损能翻一倍以上。这种分布形态给规划部门的启示是按“平均损耗”做经济评估会低估极端场景下的运行成本用蒙特卡洛下网损的整个分布做期望成本计算才更贴近实际。实操心得结果解读阶段我强烈建议多画概率密度图而不是只画CDF。CDF虽然在学术论文里常见但给人看的时候PDF能直观展现双峰、偏态这类形态特征。比如风光接入后某些节点电压有可能出现“双峰”分布——一个峰来自光伏大发的高电压场景另一个峰来自光伏停摆的低电压场景。这种双峰特征在均值里完全看不出在PDF图里一眼就能看到。5. 常见问题与工程避坑指南5.1 潮流不收敛、NAN满天飞怎么办最常遇到的报错就是潮流计算不收敛。从我自己的排查经验看80%的原因是DG出力和负荷匹配出现极端组合比如某个极端场景下DG满发而负荷特别低馈线末端电压被抬高到1.1 pu以上潮流发散。排查步骤很固定第一先把algorithm换成bfw配电网场景下它的稳健性远好于nr第二检查DG接入容量是不是超出网络承载能力如果DG总容量超过系统峰值负荷的一半以上建议先把DG容量降下来再跑第三检查负荷波动标准差0.05的正态扰动对IEEE33来说比较温和改成0.1就可能在重负荷场景出现负的有功功率这在实际电网里是有问题的。5.2 样本量不够导致结果方差大一个很典型的症状是前后两次跑同样的代码节点电压的5%分位数差一大截。原因很简单样本量太小。IEEE33节点数少3000次潮流几十秒就跑完我建议直接上5000到10000次反正计算代价不大。如果真嫌慢有几个加速手段。首推拉丁超立方采样Latin Hypercube Sampling它比简单随机采样更聪明把每个输入变量的分布区间等分成N份然后在每一份里采样一次确保样本均匀覆盖整个概率空间。样本数相同的情况下LHS下输出变量方差的收敛速度远超朴素随机采样。用scipy实现LHS只需要几行from scipy.stats import qmc sampler qmc.LatinHypercube(d2, seed42) samples sampler.random(nn_samples) # 将均匀分布的样本转换为威布尔和Beta分布样本 v_wind distributions_weibull.ppf(samples[:, 0]) g_irrad distributions_beta.ppf(samples[:, 1]) * g_max5.3 随机变量相关性被忽略的后果我这版代码里风力和光伏是独立采样的但实际上同一地区的风资源和光照往往存在互补性——白天光照强的时候风速往往偏小夜间光照为零但风速可能抬升。忽略这种负相关性会低估系统电压波动的中间状态、高估极端星座的出现概率。处理方式比较主流的是用Copula函数或者先构造相关正态样本再做变换。第一版独立采样可以做baseline但如果你要评估真实的DG接入场景建议把风光的负相关设成-0.3到-0.5之间结果会更接近实际。5.4 随机种子与结果可复现蒙特卡洛的随机性既是优点也是坑。同一个程序不设置随机种子两次跑出来的结果会在统计波动范围内略有不同。发论文、做报告时被问“为什么我这个墩子电压越限概率是8%你的是15%”就很难受。解决办法很简单程序开头设一个固定的np.random.seed(42)。但还有一个更隐蔽的问题是pandapower在内部可能也会调用随机函数巧了就会影响全局随机状态。所以严谨的做法是每一轮循环结束后不依赖全局随机状态而是显式传入随机数生成器。代码里体现为rng np.random.default_rng(42) v_wind rng.weibull(k_shape) * c_scale g_irrad rng.beta(alpha_irr, beta_irr) * g_max用default_rng替代全局随机函数隔离性更好不同版本的numpy下结果也更稳定。5.5 从静态蒙特卡洛到时序仿真的扩展思路这篇博文里跑的是静态概率潮流也就是把风速、光照、负荷都当独立同分布的随机变量一轮轮独立采样。但在实际工程里风速是有惯性、光伏是白天才有、负荷有日峰谷特性的一个时段的场景和下一个时段强相关。要模拟这种时间相关性需要把蒙特卡洛升级成时序版本用历史风速曲线或ARMA模型生成时间序列的风速样本在每个时刻点跑一次潮流最后统计整个时间窗口内的概率指标。代码架构和静态版几乎一样只是采样部分改用时间序列模型主循环不变。这也是我现在在往前探索的方向——静态方法先把算法和流程跑通时序版本解决工程落地难题。最后分享一点个人体会跑了这几个月蒙特卡洛概率潮流最大的感悟是在这个问题里真正的难点从来不是算法本身而是你对自己的随机模型有多信任。蒙特卡洛就是个“你说输入什么分布我就按什么分布给你采样”的工具如果风速的威布尔参数取错了、光伏Beta分布的形状参数拍脑袋定大了后面算出的所有越限概率都是精致垃圾。所以做概率分析之前先把基础数据的统计特性摸清楚比研究更多采样方法更重要。另外就是结果解释的功力——概率潮流输出的不是一个数而是一堆分布学会用这些分布去讲清楚工程风险才是这个工具最大的价值所在。你上手跑通之后建议也先拿IEEE33把流程走顺再去碰你手头真实电网会从容很多。