做配电网仿真的朋友应该都有同感：只要系统里接了光伏和风电，传统的一次潮流计算就有点"不够用了"。晴天、大风天、用电低谷、负荷尖峰，每个时刻的注入功率都在变，节点电压也跟着来回波动。你想知道某个节点电压超限的风险到底有多大，光算一次静态潮流给不出答案。这个项目就是围绕这个问题展开的——在IEEE33节点系统上，用蒙特卡洛概率潮流来量化风光出力不确定性对电网的影响，把"会不会越限"从猜测变成概率数字。
我先说结论：概率潮流不是要替代传统潮流，而是给运行分析多加一个维度。传统潮流给你一个确定的工作点，概率潮流给你一条完整的概率分布曲线。对于配电网规划、分布式电源接入评估、无功补偿配置这类场景，这个分布信息比单点结果有价值得多。本文会从数学建模、系统搭建、代码实现到结果解读，完整走一遍基于蒙特卡洛的概率潮流计算流程，适合刚接触概率潮流的电力专业学生、做分布式电源接入评估的工程师，以及想用Python快速验证配电网方案的开发者参考。
1. 概率潮流到底解决什么问题
1.1 为什么单次潮流计算不够用
传统潮流计算的基本思路是：给定网络拓扑、节点注入功率和平衡节点电压，通过求解功率平衡方程得到各节点电压幅值和相角。这是确定性计算，输入是什么，输出就是什么。问题在于，当系统中含有风电和光伏时，输入本身就是随机变量——风速每小时都可能从3 m/s跳到12 m/s，云层飘过光伏出力瞬间掉一半，这种情况下你算一次潮流，只能代表某个特定时刻的运行状态。
举个例子：某节点接了200 kW光伏，你按额定出力算潮流，电压合格；但你按阴天实际出力算，电压可能偏低。那到底该按哪个工况设计？如果按最恶劣工况设计，投资浪费；如果按理想工况设计，风险兜不住。工程上需要的是一个区间、一个分布，而不是一个点。
概率潮流（Probabilistic Load Flow, PLF）就是干这个事的。它的核心思想是：把输入随机变量（风、光出力、负荷波动）的概率分布，通过潮流方程传播到输出端，得到节点电压、支路功率、网损等指标的概率分布或统计特征值。这样就能回答"节点18的电压低于0.93 p.u.的概率是3.2%"这类问题，决策就有了量化依据。
1.2 蒙特卡洛方法在这个场景的优势
求解概率潮流有两类主流方法：解析法和模拟法。解析法基于线性化潮流方程或半不变量法，速度快，但要处理随机变量间的相关性比较麻烦，而且当系统重载、方程非线性程度高时误差会变大。模拟法的代表就是蒙特卡洛，思路更直接：按输入变量的概率分布海量抽样，对每个样本做一次确定性潮流计算，最后对结果做统计。
蒙特卡洛的优势主要体现在三个方面。第一是原理简单，几乎不需要对潮流方程做任何化简，牛顿-拉夫逊、PQ分解、前推回代都直接能用，精确度只取决于抽样规模。第二是灵活性极强，不管输入变量是正态分布、威布尔分布还是贝塔分布，不管变量之间是否相关，都能处理。第三是结果信息完整，可以得到完整的经验分布函数，而不只是均值和方差。缺点也明显——计算量大。一次抽样算一次潮流，抽样一万次就是一万次潮流，在33节点这种小型系统上压力不大，但在几百上千节点的大规模系统上就需要并行计算、拉丁超立方抽样等提速手段了。
对IEEE33节点这种经典测试系统来说，蒙特卡洛完全是够用的。系统规模小，潮流计算快，哪怕抽样两万次，在普通PC上也就是几分钟的事。这也是我推荐新手从蒙特卡洛入手的原因——先搞清楚概率潮流的逻辑，再去研究解析法，会顺很多。
2. 风光不确定性的数学建模
2.1 风速、光照与出力的随机模型
要做概率潮流，第一步是给风光出力建立概率模型。风电出力通常由风速决定，工程上习惯用威布尔分布描述风速的统计规律，其概率密度函数为：
f(v) = (k/c) * (v/c)^(k-1) * exp(-(v/c)^k)
其中 v 是风速，k 是形状参数，c 是尺度参数。实际项目里，k 常取2左右（也就是瑞利分布的特殊情形），c 可以根据当地年平均风速估算。得到风速分布之后，再通过风机的功率特性曲线把风速映射为出力：切入风速以下出力为0，额定风速以上出力为额定值，中间段近似线性。
光伏出力建模类似，但逻辑稍微不同。光照强度的随机性常用贝塔分布描述，归一化的光照强度 r（0到1之间）的概率密度为：
f(r) = Γ(α+β) / (Γ(α)Γ(β)) * r^(α-1) * (1-r)^(β-1)
其中 α 和 β 是形状参数，可以根据光照均值与方差反推。光伏出力 P 近似等于光照强度乘以额定容量再乘一个效率系数，用 P = r * P_N * η 即可。这个模型忽略温度影响，工程精度够用。
2.2 蒙特卡洛抽样如何生成出力场景
模型建好之后，就是蒙特卡洛的核心环节——抽样。用计算机生成服从上述分布的随机数。以Python为例，scipy.stats 库提供了现成的方法：weibull_min.rvs(c, scale=c_param, size=n) 可以直接抽风速样本，betabinom 或者 beta.rvs(alpha, beta, size=n) 可以抽光照样本。把风速样本代入风机功率曲线，把光照样本代入光伏出力公式，就得到一组风光出力场景。
这里有一个容易忽略的点：抽样不是抽完就完事了。如果只是各自独立抽样，得到的场景集会高估极端组合出现的概率——实际中往往大风天和晴天存在一定互补性。比较严谨的做法是引入相关系数，用拉丁超立方采样 + Cholesky分解处理相关性，或者用Copula函数联合建模。如果实在不想做复杂相关性处理，至少也要先把独立抽样的结果算一遍，作为上下限参考。
2.3 随机变量的数量级选择
抽样规模直接影响结果精度，也直接影响计算耗时。理论上样本数越多，经验分布越接近真实分布。但工程上不需要追求极致精度——出力概率模型本身也有误差，算到小数点后第三位意义不大。根据我的经验，节点电压的均值和标准差在抽样5000次左右开始趋于稳定，分位数（比如5%和95%分位）在1万到2万次时比较可靠。IEEE33节点系统单次潮流计算在毫秒级，抽样1.5万次总耗时大概几十秒到几分钟，完全可控。
有个检验抽样规模是否足够的经验方法：把样本量从5000增加到10000，再看结果指标变化。如果节点电压均值和标准差的波动小于0.1%，说明样本量够了；如果还在明显跳动，就继续加量。
3. IEEE33节点系统的搭建与仿真实现
3.1 标准算例的电气参数
IEEE33节点系统是配电网研究里最常见的测试系统，原始模型来源是单馈线辐射状配电网，基准电压12.66 kV，基准功率10 MVA，总负荷约3715 kW加2300 kvar。系统包含33个节点、32条支路，另有5条联络开关支路，正常运行时联络开关打开，呈辐射状结构。
在仿真前要把这些参数整理成程序能读的结构。常见的数据格式是三个表：节点表（节点编号、节点类型、有功负荷、无功负荷）、支路表（首端节点、末端节点、电阻、电抗）、发电机表（接入节点、出力类型、额定容量）。注意IEEE33的标准参数里节点编号从0到32还是从1到33，不同版本不一样，自己用的时候要统一，别混。
我建议直接把数据写成CSV或Python字典形式，方便后续多次调用。这个环节看着琐碎，但其实特别重要——后面所有潮流计算的结果，都是建立在这一步数据准确的基础上。我见过不少人在这一步抄错一个阻抗值，导致后面结果怎么算怎么不对劲，排查半天。
3.2 风光接入位置与容量怎么选
风光接入节点的选择，直接影响概率潮流的结果。不同的接入位置，对电压分布的影响差异很大。常规做法是选在馈线末端附近——那些地方本来就是电压薄弱区，接上波动性电源后电压波动更明显，做概率潮流分析更有代表性。我在这个项目里把风电（额定容量约400 kW）接在节点17附近，光伏（额定容量约500 kW）接在节点32附近，这两个节点离电源点远、电压偏低，是典型的风险监测点。
接入容量也要合理。总光伏加风电容量不能超过系统总负荷太多，否则可能出现大量潮流倒送，反而掩盖了我们要观察的"波动影响"。按传统配电网经验，分布式电源总渗透率控制在20%到30%比较合适，既能看出概率潮流的效果，也不至于让系统状态过于极端。
3.3 核心流程与代码实现
蒙特卡洛概率潮流的完整流程可以归纳为以下几步：
- 初始化：读入网络数据，设定抽样总数N、接入风光节点及参数。
- 生成场景：按风速威布尔分布抽样，通过风机功率曲线得风出力；按光照贝塔分布抽样，得光伏出力。
- 更新负荷：把风光出力作为负的有功注入叠加到对应节点上。
- 潮流计算：对每个场景执行确定性潮流算法，记录关键节点电压、支路功率。
- 统计分析：对全部结果做概率统计，输出均值、标准差、概率分布、越限概率等。
下面给出核心Python实现思路，不是完整工程代码，但结构可以直接用：
import numpy as np from scipy import stats N = 15000 # 抽样次数 n_nodes = 33 # 风电机组接入节点、额定参数 wind_node = 17 wind_rated = 0.4 # MW v_in, v_rated, v_out = 3.0, 12.0, 25.0 # 切入/额定/切出风速 # 光伏接入节点、额定参数 pv_node = 32 pv_rated = 0.5 # MW pv_eta = 0.95 # 逆变效率 # 1. 抽取风速和辐照度样本 k_shape, c_scale = 2.1, 7.5 # 威布尔参数 v_samples = stats.weibull_min.rvs(k_shape, scale=c_scale, size=N) r_samples = stats.beta.rvs(5.0, 4.0, size=N) # 光照归一化 # 2. 风速 -> 风电出力 P_wind = np.zeros(N) P_wind[(v_samples >= v_in) & (v_samples < v_rated)] = \ wind_rated * (v_samples[(v_samples >= v_in) & (v_samples < v_rated)] - v_in) / (v_rated - v_in) P_wind[v_samples >= v_rated] = wind_rated # 3. 辐照度 -> 光伏出力 P_pv = r_samples * pv_rated * pv_eta # 4. 对每个样本，把注入功率拼入节点功率列向量，调用潮流函数 for i in range(N): P_inj[wind_node] = P_load[wind_node] - P_wind[i] P_inj[pv_node] = P_load[pv_node] - P_pv[i] V_mag, V_ang, P_line, Q_line = run_power_flow(..., P_inj, Q_inj) V_result[:, i] = V_mag P_line_result[:, i] = P_line这里的run_power_flow就是普通的牛顿-拉夫逊潮流函数，参数收敛精度设为1e-8即可。如果你的项目需要更高的抽样效率，可以把for i in range(N)改成向量化批处理，把潮流算子针对矩阵输入重写，速度能提升一个量级，但编程复杂度也上来了。新手做课题或方案验证，逐样本循环完全够用。
4. 结果分析：怎么读概率潮流算出来的数据
4.1 节点电压的概率分布怎么看
跑完蒙特卡洛后，每个节点都会得到 N 个电压幅值样本。最常见的做法是画电压分布直方图和累计概率曲线，然后提取几个关键统计量：均值、标准差、5%分位数、95%分位数。对于接入风电或光伏的节点，电压分布通常是偏态的——光伏大发时电压抬升，光伏不出力时电压跌落，两头的样本都会累积，分布可能呈现双峰趋势。
这里有一个分析技巧：别只盯着电压期望值。期望值正常不代表风险低。举个例子，某节点电压期望0.97 p.u.，看似合格，但它的5%分位可能低到0.88 p.u.，已经超过低压限值。这类信息只有在概率视角下才能暴露出来。所以在报告里，我建议同时给出"期望值加误差棒"和"越限概率"两个维度，比如：
- 节点17电压期望 0.962 p.u.，标准差 0.018 p.u.，低于 0.93 p.u. 的概率为 4.7%
- 节点32电压期望 0.971 p.u.，标准差 0.021 p.u.，高于 1.05 p.u. 的概率为 1.2%
这样才算把概率潮流的价值发挥出来了。
4.2 支路潮流与网络损耗的越限风险
支路功率的分析和电压类似，重点看支路负载率。IEEE33系统的支路容量在数据表里不总是显式给出，通常需要自己设一个基准容量，比如按导线载流量折算。算完所有样本后，可以得到每条支路的有功功率概率分布，挑负载率超过80%或100%的支路作为瓶颈支路，对其进行标记。在概率潮流框架下，"支路过载风险"是一个概率数字，比如"支路6-7过载概率0.8%"，这就比单纯说"这条线重载"更精确。
网络损耗方面也值得统计。由于风光出力的波动性，同一系统不同场景下的网损差异可能很大。你可以计算网损的期望值和概率分布，用来评估分布式电源接入后对整体经济性的影响。我在这个项目中看到的一个典型现象是：光伏接入后，日间网损期望值下降，但重负荷加上弱出力场景下网损可能比不接分布式电源还高。这就是"波动性电源带来的双重影响"，做规划时必须权衡。
4.3 影响范围的界定方法
所谓"影响范围"，指的是风光出力不确定性对系统运行指标的波及程度。一个实用定义是：如果某节点电压的标准差因为接入风光而显著增大（比如比未接入时增大超过0.005 p.u.），就认为该节点落在影响范围内。按照这个标准，电压波动影响通常集中在接入点附近2到3个节点之内，离接入点越远影响越弱。但这个规律在馈线末端可能会被打破——末端节点自身电压基数低，对注入功率变化更敏感，即使离接入点不太近，波动也可能不小。
所以做影响范围分析时，不要只看电气距离远近，要结合节点本身的电压灵敏度来判断。实用做法是跑两组仿真：一组是无风光接入的确定性潮流，一组是含风光接入的概率潮流，对比每个节点电压分布的差异。差异显著的节点就是重点监测对象，后续可以在这些节点配备动态无功补偿或者储能，改善电压质量。
5. 常见问题与排查技巧实录
5.1 蒙特卡洛抽样结果不够稳定
很多人跑完发现方差大、分布曲线毛糙，第一反应是增加抽样次数。但如果抽样次数已经到两三万还不见改善，就要检查随机数种子和抽样方法了。这里一个关键点是：直接用基础的np.random做独立抽样，在高维场景下分布覆盖能力差，尾部概率估计不准。改用拉丁超立方采样能明显改善稳定性——它把每个输入变量的分布区间等分成 N 个小区间，在每个区间内各抽一个点，确保抽样点均匀覆盖整个概率空间。
另外要注意随机数种子固定。固定种子可以让实验结果可复现，这一步在写论文和做方案评审时特别重要。不固定种子的话，前后两次运行结果稍有差异，会给对方造成算法不稳的错觉。
5.2 潮流计算不收敛怎么办
蒙特卡洛里最头疼的就是某个场景潮流不收敛。根源通常是抽样抽到了极端场景——比如风电出力很高而负荷很低，节点电压被抬得过高，牛顿-拉夫逊法迭代发散。排查方法有三个方向：
第一，检查潮流算法的收敛门槛，把迭代次数上限从默认的10次加到50次，收敛精度设为1e-8，很多发散其实是迭代次数不够。第二，对极端出力场景做限幅，把节点注入功率限制在系统允许的最大/最小范围内。第三，如果某个场景反复不收敛，直接跳过该场景并计数，在统计时扣除这些无效样本。但要注意，如果跳过比例超过1%，说明抽样分布或系统参数设置有问题，不能简单无视。
5.3 风光出力数据的相关性处理
做独立抽样时，风电和光伏的出力是互不相关的，但实际情况未必如此。有些地区白天风小、晚上风大，光伏和风电天然负相关；另一些地区二者同向波动。忽略相关性会导致概率潮流结果偏差，尤其在评估系统总出力波动时。
处理相关性的一个可行方案是先构造相关系数矩阵，用Cholesky分解对独立正态样本做线性变换，再把变换后的样本映射回目标分布。这在工程上有一个成熟的组合方案：先用拉丁超立方生成独立均匀分布样本，变换成标准正态样本，施加相关性后，再通过等概率变换映射回威布尔和贝塔分布。这套流程做下来，能在保留各变量边缘分布特性的同时，还原它们之间的相关性。
5.4 结果呈现与报告输出的建议
仿真计算只是项目的一半，怎么把结果讲清楚同样关键。我的习惯是输出三类材料。第一是概率密度图，每个关键节点画一张，直观展示电压波动形态；第二是概率表，列出各节点的均值、标准差、5%/95%分位数和越限概率；第三是越限风险热力图，在IEEE33节点拓扑图上用颜色标注各节点的电压越限概率，一眼就能看出风险集中在哪个区域。
如果要把结果写进技术报告，核心结论不要只给数字，要结合工程判断。比如"节点32电压高于1.05 p.u.的概率为1.2%"之后，建议再补一句"建议在该节点配置储能或无功补偿，将高电压风险控制在0.5%以内"，这样就把概率潮流的结果落地到实际工程措施上了。
我在实际做这个项目的过程中，最大的体会是：蒙特卡洛概率潮流的编程门槛其实不高，难的是对结果的解读和建模选择的合理性判断。先花时间把风光出力模型的参数选准，再去调抽样规模和潮流算法，路就走顺了。后续如果想做更深入的扩展，可以考虑在同一个IEEE33平台上加入储能优化调度，对比储能投运前后的电压越限概率变化——这正好能承接本文的概率分析框架，让仿真平台继续发挥价值。
