做配电网规划、调度或者新能源接入评估的朋友应该都有过这种体验潮流计算跑出来永远是那组标准答案负荷多大就对应多大的潮流光伏发多少就是多少。但真实系统不是这样运转的光照会飘过一片云风速会突然降下来电动汽车充电高峰可能恰好撞上晚峰负荷。最近我在做电力系统不确定性分析重点就是用蒙特卡洛法处理概率潮流计算拿IEEE33节点配电网当测试床把风光出力模型接进去反复采样跑了几千次确定性潮流最后得到的不再是一组数而是一堆概率分布。这篇文章把我的完整思路、代码实现、结果解读和踩过的坑都整理出来给同样在做概率潮流、分布式电源接入评估的朋友做个参考。1. 为什么确定性潮流算不准概率潮流的出发点1.1 确定性潮流的隐含前提与局限普通的潮流计算核心方程就是节点注入功率和电压之间的关系比如P_i jQ_i U_i * (ΣY_ij * U_j)*。给定一组确定的负荷和电源出力牛顿-拉夫逊法或者前推回代法迭代几次就能得到每个节点的电压幅值和相角、每条支路的功率。这个方法本身没有任何问题问题在于输入数据你把负荷取成最大值光伏出力取成典型值算出来的是这个时刻的答案不是这个系统的真实状态。新能源大规模接入之后这种单点计算的局限被放大了。光伏出力随着辐照度剧烈波动一片云的遮挡就能让某台逆变器出力掉去大半风电更不用说风速从切入风速到额定风速之间出力几乎和风速的三次方挂钩换个天气就是两个系统。如果再叠加负荷本身的随机波动你会发现系统的运行状态根本不是一个点而是一团云一个高维的随机变量空间。1.2 概率潮流的思路把输入的不确定性传播给输出概率潮流的想法很直接既然输入是随机的那就用概率分布来描述输入把这个分布通过潮流方程传播出去得到输出的概率分布。这样我们就能回答一类传统潮流完全无法回答的问题节点18的电压越限概率是多少线路5-6的负载率超过80%的可能性有多大系统最薄弱的环节到底在哪概率潮流的求解方法有好多流派比如点估计法、一次二阶矩法、无迹变换、多项式混沌展开还有基于蒙特卡洛模拟的方法。在实际工程里我最终还是选了蒙特卡洛法。理由很朴素第一它几乎不做模型简化潮流方程该是非线性就是非线性不用泰勒展开去近似这对配电网这种经常出现较大电压偏差的场景很重要第二实现极度直观就是采样—算潮流—统计代码量少出问题也好排查第三蒙特卡洛的收敛性和采样次数有关我做几百次大致的趋势就有了做几千次精度完全够用配合并行计算时间完全可以接受。要说缺点也明显就是计算量大。但IEEE33这种规模的配电网单个潮流计算毫秒级完成采5000次也就几十秒完全不是瓶颈。正因为这样蒙特卡洛法在配电网不确定性分析里反而是最实用的起步方案。1.3 为什么选IEEE33节点作为测试对象IEEE33节点系统是配电网领域的标准测试题就像学习机器学习先用MNIST数据集一样。它的规模不大不小33个节点、32条支路、基准电压12.66kV、总负荷大约3715kW加2300kvar有主干、有分支末端电压偏低这种典型配电网问题它全都有。接入分布式电源之后电压抬升、潮流反向、局部过载这些现象都能复现出来。而且数据公开自己手工录入或者从pandapower这类工具里直接加载都行方便复现和对比。用它当小白鼠再合适不过。2. 风光出力模型搭建Beta分布和Weibull分布怎么选参数2.1 光伏出力为什么用Beta分布光伏出力本质上取决于光照辐照度。从物理机制上看辐照度受到太阳高度角、云层厚度、气溶胶等多种因素影响在很多时间尺度上表现出较强的随机性。工程界比较常用的做法是用Beta分布描述辐照度在某个时段内的统计特性再做线性转换得到出力。关键是为什么是Beta而不是正态分布。正态分布虽然好算但它两端是无限延伸的而光伏出力有天然的上下界——不可能小于0也不可能超过装机容量。Beta分布定义在[0,1]区间正好完美匹配这个物理约束。规范化之后光伏出力标幺值pp P / P_nom的概率密度函数可以写成f(p) [p^(α-1) * (1-p)^(β-1)] / B(α, β)其中B(α, β)是Beta函数α和β是两个形状参数。这个结构决定了Beta分布可以呈现左偏、右偏、近似均匀等各种形态对晴天、多云、阴天三种天气都能适应。参数α、β一般不是凭感觉拍的用历史出力数据可以估计出平均值μ和标准差σ然后按下式反推α μ²(1-μ)/σ² - μβ μ(1-μ)²/σ² - (1-μ)我举个例子。假设某光伏电站午间出力标幺值的均值μ0.3标准差σ0.15典型多云天气那么σ²0.0225算下来α2.5β5.83。这个β明显大于α分布曲线偏向低出力侧符合多云场景出力整体偏低的直觉。如果是有太阳的晴天数据μ可能到0.6以上α就会大于β分布右偏。注意α、β必须大于0如果算出来小于零基本就是μ和σ的取值不合理要回头检查数据处理。2.2 风速的Weibull分布与风速-功率转换曲线风力发电的源头是风速风速的长期统计特性在大多数地区都能用两参数Weibull分布描述。概率密度函数为f(v) (k/c) * (v/c)^(k-1) * exp(-(v/c)^k)其中k是形状参数通常取1.5到3之间k越大风速分布越集中c是尺度参数和平均风速正相关可以近似取平均风速的1.13倍左右。很多风电资源评估资料里典型k2此时Weibull退化为瑞利分布、c6到8 m/s我用k2、c6.5 m/s作为基础参数。拿到风速样本之后还要过一道风速-功率转换曲线才能得到风电出力。工程上最常用的是三段式功能关系风速低于切入风速v_ci出力为0风速在v_ci和额定风速v_r之间出力按线性或二次曲线上升风速在v_r和切出风速v_co之间出力保持额定功率超过v_co为了保护风机出力直接切到0。我这里采用线性近似切入风速3m/s、额定风速12m/s、切出风速25m/s装机容量取500kW。这样风速和功率的关系就变成了一个可写进代码的分段函数。值得注意的是用三段式转换之后风电出力的概率分布会出现两个特殊的尖峰一个是0出力附近风速低于切入或高于切出一个是额定出力附近风速在额定与切出之间。这是风电出力的真实特征做概率潮流时不要觉得分布不平滑就是代码错了。2.3 风光相关性先独立假设但要知道它的局限实际中风光出力往往是负相关的因为太阳辐射强的时候通常大气稳定、风速偏低而风大的阴雨天光伏出力又不高。严格的概率潮流应该用Copula理论或者Cholesky分解构造具有相关性的联合采样。但对于工程初筛先用独立采样把框架建立起来完全可行。我自己在研究阶段也先用独立假设跑通了再叠加相关性。后面踩坑清单里我会专门讲相关性带来的结果偏差方向先留个悬念。3. IEEE33节点测试网数据准备与风光接入位置选择3.1 网络基础数据与工具选择IEEE33节点系统的基础数据包括33个节点的负荷功率、32条支路的阻抗参数以及网络拓扑。标准参数在公开文献里都能查到不用自己推导。节点编号一般从0开始0号节点是变电站出口平衡节点1到17构成主干还有从不同位置引出的三个分支比如2-19-20-21-22、3-23-24-25、6-26-27-28-29-30-31-32不同文献编号略有差异但拓扑逻辑一致。工具层面我强烈建议用Python环境选pandapower作为潮流计算引擎。pandapower是专门面向配电网分析的开源库内置了IEEE33等测试系统直接用net pp.networks.case33bw()就能加载省去手工录入数据的大把时间。潮流求解默认用牛顿法对于IEEE33这种规模速度和收敛性都很好。还要配合numpy做采样和统计matplotlib做可视化。3.2 风光接入位置与容量方案分布式电源接入位置直接影响概率潮流的结果这也是分析本身要回答的问题之一。我设计了两个场景场景一光伏和风电集中接在馈线末端节点17接入光伏600kW节点32接入风电500kW用来模拟末端高渗透率最恶劣的情况场景二把相同容量的电源分散接入多个节点比如节点8和22各接300kW光伏、节点25和32各接250kW风电对比分散接入对电压分布的改善效果。接入模型上用sgen静态发电机表示控制方式按最常见的PQ型并网逆变器处理不做电压支撑的PV节点处理。这是有讲究的现在绝大多数分布式光伏逆变器都运行在单位功率因数或指定功率因数下出力只给有功无功按恒功率因数算所以广义上依然是PQ节点跑潮流最容易收敛也最贴近工程现状。3.3 负荷侧要不要做随机化把风光出力设成随机变量之后很多人容易漏掉负荷。负荷本身也有时变性和随机性典型做法是让每个节点的负荷在基准值上叠加一个正态扰动比如p_load_new p_load_base * (1 0.05 * N(0,1))。我一开始只随机化了电源结果电压波动范围比实际情况窄因为负荷的随机响应被完全忽略了。加上负荷扰动之后输出分布的方差明显变大更接近真实运行状态。不过要注意负荷扰动尺度别取太大配电网负荷预报误差一般控制在5%到10%以内比较合理。4. 蒙特卡洛概率潮流代码实现采样—算潮流—做统计的完整链路4.1 采样模块风光出力样本生成函数首先写采样函数。这部分我踩过一个坑就是numpy的weibull函数和手写公式的对应关系。numpy.random.weibull(a, size)生成的是形状参数为a、尺度参数为1的Weibull分布如果要获得尺度参数为c的样本需要在结果上乘以c。很多人把这个c忘了导致风速均值整体偏小。import numpy as np def sample_pv_output(alpha, beta, cap_pv, n_samples): 光伏出力采样Beta分布 p_pu ~ Beta(alpha, beta)乘以装机得到实际出力(MW) p_pu np.random.beta(alpha, beta, n_samples) return p_pu * cap_pv def sample_wind_output(k, c, v_ci, v_r, v_co, cap_wind, n_samples): 风电出力采样Weibull分布风速 三段式风功率转换 返回实际出力(MW) v np.random.weibull(k, n_samples) * c p np.zeros(n_samples) # 线性段切入风速到额定风速 mask_linear (v v_ci) (v v_r) p[mask_linear] cap_wind * (v[mask_linear] - v_ci) / (v_r - v_ci) # 额定段额定风速到切出风速 mask_rated (v v_r) (v v_co) p[mask_rated] cap_wind # 其余情况出力为0 return p这里的光伏参数alpha和beta就用前面公式从均值标准差反推。例如多云场景取alpha2.5、beta5.83装机0.6MW风电参数取k2、c6.5、v_ci3、v_r12、v_co25、装机0.5MW。4.2 主循环赋值—潮流—记录采样函数准备好之后主循环的逻辑很清晰每次采样得到一组光伏出力、风电出力以及可选的负荷扰动把它们写入当前网络调用runpp做确定性潮流计算把节点电压幅值标幺值和支路电流记录下来循环N次。import pandapower as pp from pandapower import networks # 加载IEEE33节点测试系统 net networks.case33bw() # 添加光伏和风电sgen接入节点17和32 pp.create_switch(net, 0, 0, etb) # 保持原有平衡节点不变即可 pp.create_sgen(net, 17, p_mw0.0, namepv_node17) pp.create_sgen(net, 32, p_mw0.0, namewind_node32) # 记录负荷基准值 load_base net.load[p_mw].values.copy() n_samples 3000 n_bus net.bus.shape[0] voltage_records np.zeros((n_samples, n_bus)) loading_records np.zeros((n_samples, net.line.shape[0])) cap_pv 0.6 # 光伏装机MW cap_wind 0.5 # 风电装机MW # 预生成所有采样值 pv_samples sample_pv_output(2.5, 5.83, cap_pv, n_samples) wind_samples sample_wind_output(2, 6.5, 3, 12, 25, cap_wind, n_samples) for i in range(n_samples): # 更新分布式电源出力 net.sgen.loc[0, p_mw] pv_samples[i] # 光伏 net.sgen.loc[1, p_mw] wind_samples[i] # 风电 # 可选给负荷加5%正态扰动 net.load[p_mw] load_base * (1 0.05 * np.random.randn(net.load.shape[0])) try: pp.runpp(net) voltage_records[i, :] net.res_bus.vm_pu.values # 支路电流标幺值也可以记录这里按线路负载率记录 loading_records[i, :] net.res_line.loading_percent.values except pp.LoadflowNotConverged: # 个别样本可能不收敛记录NaN后面统计时跳过 voltage_records[i, :] np.nan loading_records[i, :] np.nan注意net.sgen.loc[0, p_mw]这里我用的是sgen在DataFrame里的行号索引不是节点编号。如果sgen创建顺序是先光伏后风电那就是0和1。实际操作中建议创建之后先打印net.sgen确认行号对应关系避免赋值赋错对象。4.3 统计模块均值、标准差、越限概率蒙特卡洛做完之后核心是把几千组结果压缩成有价值的统计量。我主要关注四类指标节点电压均值、节点电压标准差、节点电压越限概率、支路负载率越限概率。配电网电压合格范围国内工程一般取0.93到1.07 pu分布式电源接入评估时有些地区要求更严格取0.95到1.05。我这里按0.93到1.07做基准。# 剔除不收敛样本 valid ~np.isnan(voltage_records[:, 0]) voltage_valid voltage_records[valid, :] loading_valid loading_records[valid, :] # 均值与标准差 voltage_mean np.nanmean(voltage_valid, axis0) voltage_std np.nanstd(voltage_valid, axis0) # 越限概率 p_over np.mean(voltage_valid 1.07, axis0) p_under np.mean(voltage_valid 0.93, axis0) p_voltage_violation p_over p_under # 支路重载概率负载率超过80% p_load_heavy np.mean(loading_valid 80, axis0)这些统计量做出来后概率潮流的核心产出就齐了。均值告诉你最可能的运行状态标准差告诉你波动有多大越限概率告诉你风险有多大。这三样正是确定性潮流给不了的。4.4 样本量到底取多少用变异系数判断很多人第一次做蒙特卡洛都会问采样次数取多少合适拍脑袋取1000、5000、10000都行但更严谨的做法是用变异系数来控制。对于某个统计量比如节点电压均值其变异系数大约和1/sqrt(N)成正比。工程经验上N取1000到5000就能让均值估计的变异系数降到1%到2%尾部概率越限概率需要更多样本才能稳定特别是越限概率本身很小时。如果你关心的是0.1%的极端事件那就不能只跑5000次因为0.1%概率的事件在5000次采样里平均只出现5次估计非常粗糙。这种场景下最少要几万次采样或者配合重要抽样、子集模拟这类方差削减技术。我自己的习惯是先跑2000次快速看趋势再根据目标精度的变异系数决定是否加量。别一上来就跑十万次等你发现参数写错了再改白白浪费几小时。5. 结果看什么电压分布、越限概率与潮流实况解读5.1 末端电压的分布特征均值塌陷与方差放大跑完3000次采样第一件事就是画电压分布图。把全部节点的电压均值画成曲线你能看到一条从首端到末端逐步下降的曲线0号节点接近1.0主干末端节点17的电压均值可能已经降到0.95左右。这正是IEEE33系统的经典特征重负荷、长主干、末端电压支撑不足。更有意思的是标准差曲线。你会看到从首端到末端电压标准差是逐段放大的末端节点的标准差可能是中段节点的两三倍。这个现象背后的物理原因是末端节点离电源点电气距离远等效阻抗大同样的功率波动在末端引起的电压波动更大。这就是为什么分布式电源接入评估里末端的电压质量问题永远是重点关注对象。5.2 越限概率地图找到系统的危险节点电压均值下降不代表一定越限真正的决策依据是越限概率。我会把所有节点的越限概率画成柱状图或者热力分布曲线发现节点17、32这些末端位置在光伏大发、负荷高峰叠加的时候电压低于0.93的概率能到5%以上。而某些中间分支节点反而可能出现电压偏高越限原因是光伏接入后反向潮流抬高了局部电压。这个信息对规划特别有用。比如末端节点17电压越限概率3.8%决策者可以据此决定是否需要加装调压器、无功补偿或者升级导线截面而不是靠等值年最大负荷最小时的那个潮流结果拍板。同样支路重载概率能暴露哪条线路在风光大发时段容易出现反向过载这对保护定值和设备选型都有直接参考意义。5.3 概率结果和确定性潮流报告的差异有一次我把新能源按额定出力直接接入、跑一次传统确定性潮流结果显示末端电压从0.95抬到了0.99一切正常看起来光伏接入利大于弊。但概率潮流跑完发现电压大于1.07的越限概率在某些节点是4%小于0.93的越限概率在另一些节点是6%。这两个答案看起来矛盾其实各自成立前者是某个极端场景的静态切片后者是整个运行时间尺度内的风险统计。做电网分析如果只看单点确定性潮流很容易被某一个场景骗了。这也是我在文章开头说的确定性潮流给的是一组数概率潮流给的是一张风险地图。6. 实战复盘样本数、相关性、并行加速与踩坑清单6.1 坑一Beta分布参数反推公式用错这是我最早踩的坑。第一次算光伏出力均值μ0.3、标准差σ0.15随手把αμ、βσ填进去结果采样出来的出力均值完全对不上。正确的反推公式我在第2节已经写出来了它的原理是Beta分布的矩估计。如果你不想推导直接用NumPy的scipy.stats.beta.fit去拟合历史数据也行但要知道fit返回的loc、scale参数与标准的[0,1]Beta不太一样需要做归一化处理。6.2 坑二sgen赋值写错行号pandapower里sgen的索引和节点编号不一致这个事真的很坑人。我一开始直接写net.sgen.loc[17, p_mw] pv_samples[i]结果是把光伏出力赋值给了节点17不对行号17对应的是第17个sgen实际可能根本没那个节点。正确做法是创建sgen后用name字段定位比如net.sgen.loc[net.sgen[name] pv_node17, p_mw] pv_samples[i]这样永远不会错。这种小问题在仿真代码里最消耗时间因为程序不报错结果却全是错的。6.3 坑三忽略负荷随机性导致方差被低估前面提到了只随机化电源、不随机化负荷算出来的电压波动范围会比实际窄。这个问题的本质是输入不确定性的维度没给全。概率潮流的意义在于逼近真实世界的随机过程少一个主要随机源输出分布就失真。我建议至少把负荷加上5%到10%的扰动哪怕用最简单的正态分布。如果要有更精细的时变特性那就需要上升到时序概率潮流复杂度高一个量级一般工程用不到。6.4 坑四相关性不做处理导致结果失真风电和光伏独立采样会低估或高估某些风险。举个具体例子如果实际风光出力负相关大风天光伏弱、晴天光伏强但风小独立采样会造成光伏满发同时风电满发出现的概率虚高进而高估电压越上限风险和线路反向过载风险。要严格建模相关性可以用Cholesky分解处理多维正态相关的秩相关系数或者用Copula把边缘分布和相关性结构分开建模。如果只是工程初判也可以取最恶劣的风光同发和风光归零两个边界场景做确定性校核作为概率结果的上下限参考。这不算严谨但在时间有限的工程交付里是实用打法。6.5 加速手段拉丁超立方与并行蒙特卡洛计算量的问题绕不开。在IEEE33上几百几千次还好如果换成几百节点的真实配电网单次潮流计算时间上升到几十毫秒一万次采样就可能要几分钟到十几分钟。两个加速手段非常有效一是拉丁超立方采样LHS。这个方法的本质是把每个输入变量的分布分层然后在每层内均匀采样保证样本更均匀地覆盖整个概率空间。同样3000次采样LHS的方差收敛速度明显优于纯随机蒙特卡洛特别是对尾部概率的估计改善更明显。二是并行计算。蒙特卡洛循环天然适合并行因为每次采样和潮流计算之间没有依赖关系。我一般用joblib或者multiprocessing.Pool把N次循环分配到CPU多核上8核机器能跑出接近6到7倍加速。注意每个子进程里要独立设置随机数种子否则多进程可能拿到相同的采样序列。from joblib import Parallel, delayed def run_single_sample(pv_val, wind_val, net_copy): net net_copy.deepcopy() net.sgen.loc[net.sgen[name] pv_node17, p_mw] pv_val net.sgen.loc[net.sgen[name] wind_node32, p_mw] wind_val # 负荷扰动... try: pp.runpp(net) return net.res_bus.vm_pu.values except: return np.full(net.bus.shape[0], np.nan) results Parallel(n_jobs8)( delayed(run_single_sample)(pv_samples[i], wind_samples[i], net) for i in range(n_samples) )这里我故意用了deepcopy而不是共用一个net对象避免并行写共享内存导致的数据竞争。代价是deepcopy有额外开销但在IEEE33这种规模上完全可接受。如果要追求极致性能可以改成在每个worker里初始化一份net然后派发采样索引减少重复深拷贝。6.6 采样序数与随机数种子最后分享一个工程习惯所有随机过程先固定随机种子。np.random.seed(42)放最前面这样你每次运行代码得到的采样序列完全一致结果可复现。做研究或者写报告的时候可复现性是基本素养。等整体链路验证没问题了再去掉种子做大批量运行获得更稳健的统计结论。我做概率潮流这轮实践的最后发现真正有用的不是那几千张潮流结果表而是系统各处风险到底有多大概率发生这个维度。确定性潮流解决的是某个时刻系统成不成立概率潮流解决的是长期运行下系统安不安全两者是互补关系不是替代关系。IEEE33这个平台让我用很小的成本把整套方法练了一遍后面再换真算例网架流程完全一致只是数据和规模变了。如果你也在做新能源接入评估建议先跑通这个小系统把风光采样、潮流循环、统计输出的链路打通再去碰复杂的真实网架能少走很多弯路。