☰
基于蒙特卡洛的IEEE33配电网概率潮流计算与风光出力不确定性分析
2026/10/3 10:30:50 网站建设 项目流程

做配电网规划、调度或者新能源接入评估的朋友,应该都有过这种体验:潮流计算跑出来永远是那组"标准答案",负荷多大就对应多大的潮流,光伏发多少就是多少。但真实系统不是这样运转的,光照会飘过一片云,风速会突然降下来,电动汽车充电高峰可能恰好撞上晚峰负荷。最近我在做电力系统不确定性分析,重点就是用蒙特卡洛法处理概率潮流计算,拿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]区间,正好完美匹配这个物理约束。规范化之后,光伏出力标幺值p(p = 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倍左右。很多风电资源评估资料里,典型k=2(此时Weibull退化为瑞利分布)、c=6到8 m/s,我用k=2、c=6.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就用前面公式从均值标准差反推。例如多云场景取alpha=2.5、beta=5.83,装机0.6MW;风电参数取k=2、c=6.5、v_ci=3、v_r=12、v_co=25、装机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, et="b") # 保持原有平衡节点不变即可 pp.create_sgen(net, 17, p_mw=0.0, name="pv_node17") pp.create_sgen(net, 32, p_mw=0.0, name="wind_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, axis=0) voltage_std = np.nanstd(voltage_valid, axis=0) # 越限概率 p_over = np.mean(voltage_valid > 1.07, axis=0) p_under = np.mean(voltage_valid < 0.93, axis=0) p_voltage_violation = p_over + p_under # 支路重载概率:负载率超过80% p_load_heavy = np.mean(loading_valid > 80, axis=0)

这些统计量做出来后,概率潮流的核心产出就齐了。均值告诉你"最可能的运行状态",标准差告诉你"波动有多大",越限概率告诉你"风险有多大"。这三样正是确定性潮流给不了的。

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_jobs=8)( 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这个平台让我用很小的成本把整套方法练了一遍,后面再换真算例网架,流程完全一致,只是数据和规模变了。如果你也在做新能源接入评估,建议先跑通这个小系统,把风光采样、潮流循环、统计输出的链路打通,再去碰复杂的真实网架,能少走很多弯路。

需要专业的网站建设服务?

联系我们获取免费的网站建设咨询和方案报价,让我们帮助您实现业务目标

立即咨询