蒙特卡洛与近似贝叶斯在概率潮流中的MATLAB实现与对比

“蒙特卡洛配合近似贝叶斯做概率潮流”，这标题一眼看去就知道是要在随机性的框架里同时玩两套思路：一套是几乎不用解释的暴力仿真，另一套是这两年在大规模推断里很受捧的“参数反演”思路。真正落到电力系统里，问题就变成一个很实际的事情——负荷波动、新能源出力随风随光照变化，最终节点电压、支路潮流的概率分布到底是什么样？哪些节点最容易越限？如果要给调度留出置信区间，应该怎么算？这篇文章就把我打磨这套MATLAB程序的完整过程拆开讲，从原理、代码骨架到踩坑记录，一次说清楚。

适合什么人看：手里有确定性潮流程序、想升级成概率版本的研究生和工程师，或者对Monte Carlo和贝叶斯推断感兴趣、想在电力系统场景里练习的人。程序本身已经把两种路线放在同一套算例上做了对比，所以读代码的过程基本就是把概率潮流从“做题”变成“上学”的过程。

1. 项目整体设计与思路拆解

1.1 首先说清楚概率潮流到底在算什么东西

常规潮流计算给出的是一个确定的运行点：给定负荷和发电出力，解出各节点电压和支路功率。但实际电网里没有什么是确定的，光伏在云层后面进进出出，充电桩在晚饭时间集体开工，负荷预测永远有误差。这时候问“电压是多少”没有意义，应该问“电压有多大可能落在这个区间里”。

概率潮流就是干这个的。它的输出不是一个数，而是一堆分布：每个节点的电压幅值概率密度、支路潮流的累计分布、越限概率的期望值。有了这些分布，调度员才能回答“在置信度95%下线路不过载的最小备用是多少”这类决策问题。

这套程序的目标很明确：输入随机源的分布参数和潮流网络拓扑，输出关键节点电压、关键支路潮流的概率特征，并提供两种计算路径的结果对比，让使用者自己判断哪种近似方法在本场景下够用。

1.2 为什么选Monte Carlo做基准，又为什么加贝叶斯方法

Monte Carlo模拟（MCS）的原理简单到不像学术界的东西：从输入分布里抽样，把每组样本代入潮流方程求解，攒几千次结果，统计分布。它的优点是无偏、稳定、容易实现，缺点是慢——每次抽样都要解一遍潮流，样本量一大就明显吃时间。

而近似贝叶斯计算（Approximate Bayesian Computation，ABC）解决的是另一个层面的问题。传统概率潮流把输入分布参数当作已知常数，比如认为负荷均值就是100MW、标准差就是10MW。但实际中这些参数本身也是估计出来的，存在不确定性。如果从贝叶斯的视角看，我们应该用一定数量的带有误差的实测数据，去更新这些参数的后验分布，再把这个参数后验传播到潮流结果里。

程序把MCS作为基准解，因为它的结果在样本足够多时可以近似认为是“真实分布”的替代。ABC方法则是那个更聪明但更依赖调参的挑战者——它尝试用近似采样避开复杂的后验解析推导，在不损失太多精度的前提下，给出比纯MCS更灵活的参数不确定性建模能力。

1.3 两种方法在场景上怎么互补

有人会问：既然MCS都能给出结果，为什么还要折腾ABC？关键在于MCS的前提是“输入参数已知且准确”。实际项目中，负荷模型用的均值和标准差通常来自历史数据的点估计，这个点估计本身就有标准误。MCS不会告诉你这个标准误对最终潮流分布的影响，ABC可以。

ABC的基本思路是：对未知参数设定先验分布，从先验中抽取候选参数，用候选参数生成模拟数据，再把模拟数据和观测数据比较，距离足够近就接受这个候选参数，否则丢弃。重复足够多次，接受下来的参数分布就是后验近似。放到概率潮流场景里，候选参数是负荷分布的均值、方差，模拟数据是量测的节点电压或关口功率，距离就是两者差值的加权范数。

这相当于在潮流程序外面套了一层过滤器：MCS告诉你“给定这个参数，潮流结果是什么分布”，ABC告诉你“在观测数据的约束下，参数和潮流结果最可能是哪些组合”。两者的组合是这个程序最大的价值点。

1.4 方案选型背后的取舍

程序最终选择MCS加ABC对比，而不是点估计法或解析法，是因为解析法虽然快，但对非正态分布、非线性潮流方程容易产生较大偏差，尤其当系统接近重载时误差难以控制。MCS慢，但能提供任意精度基准。ABC增加了介于“纯近似”和“纯仿真”之间的一个折中层次，它对应着一大类真实需求：有量测数据但不知道如何反推参数分布的工程场景。

我实际用下来最深的感受是：MCS负责兜底，ABC负责解释“为什么参数这么定”。两者在同一个程序框架里对比，能非常直观地暴露近似方法在哪一段负荷水平下失灵，这比任何理论推导都有说服力。

2. 核心细节解析与实操要点

2.1 输入随机变量怎么建模

概率潮流的精度首先取决于输入建模。不同随机源要用不同分布，这直接影响抽样和ABC先验设计。我在这套程序里默认支持三种场景。

负荷功率一般用正态分布建模，均值取预测值，标准差按“预测误差比例”设定。光伏出力更接近Beta分布，因为光照强度在0到上限之间，Beta分布的两端行为能反映出日照的饱和与零出力区间。风速则用Weibull分布，两参数分别控制尺度与形状。

在MATLAB里建模本身不难，重点是参数别拍脑袋。我的习惯是先用历史数据做一次fitdist，再人工检查拟合是否合理。

下表是程序里常用的建模对照：

注意：Beta分布的两个参数a、b直接决定概率密度形态，小样本拟合很容易过度优化，建议对fit结果做交叉验证。程序里默认开启“尾部分位数校准”，否则ABC跑到后面会发现后验完全偏到数据噪声里去。

2.2 潮流求解核心不能糊弄

概率潮流每一步迭代都要调潮流解算器，这一步的稳定性和速度决定了整个程序的生死。程序里用的是牛顿-拉夫逊法，没有依赖某个封闭工具箱，核心是自己维护的节点导纳矩阵和雅可比矩阵。

牛顿法的关键约束有几个：一是PV节点的无功越限处理，二是PQ节点电压初值要合理，三是收敛判据不能只盯着有功不平衡量，还要看无功。实际工程中遇到最多的问题不是潮流方程本身解不出来，而是某组随机样本把电压初值拉到了非收敛域，比如光伏满发且负荷极低导致的电压抬升。

程序里专门加了一个样本级的容错：每次潮流不收敛时，先补做一次平启动初值重迭代，最多重试三次，仍不收敛就记录样本编号并跳过，最后在结果汇总时给出失败率。这个设计对Monte Carlo阶段尤其重要——你绝不想让一两个坏样本毁掉整批统计。

2.3 Monte Carlo抽样的工程细节

MCS看似简单，但“怎么抽”决定了效率差异。普通平random抽样在最需要精确估计的尾部往往稀疏，所以程序里提供了三种模式：

• 简单随机抽样：对分布最熟悉的人会用randn直接出数，适合快速验证。
• 拉丁超立方抽样：把每个输入维的累积概率分成N个等间隔区间，保证每个区间都至少抽到一个点。对中等规模问题，同样精度下可以减少三四成样本量。
• Sobol序列：低差异序列，样本在整个概率空间里分布更均匀，尾部行为比随机抽样稳定很多。

我强烈建议主算用Sobol，但注意Sobol序列在维度较高时存在相关性结构，需要配合随机化处理。

实际运行时，样本量默认设为2000，这是从程序速度和精度之间折中出来的结果。做方法对比时要响应的不是“结果精度”，而是“达到同样精度所需的样本数或时间”，所以程序输出里会同时记录耗时和两个误差指标。

2.4 近似贝叶斯计算ABC的参数逻辑

ABC在这套程序里的核心是把“潮流结果模拟值”与“观测值”的距离作为接受/拒绝的判据。先用先验分布抽出候选参数θ，再去跑MCS（或简化批次），得到模拟潮流结果Y_sim，然后计算距离D(Y_sim, Y_obs)。如果D小于容差ε，就接受θ。

这里有一个工程上很关键的选择：模拟观测值用什么？我用了两种情况做对比，一是全部节点电压，二是只取关键支路功率，仿真的观测数据是自己用“真实参数”加噪声构造出来的，这在算法验证阶段是唯一可控的做法。

ABC有两个重要的坑。第一个是容差ε怎么定。ε太小，接受率极低，要跑一晚上；ε太大，后验近似像是先验的复制品，没有信息量。程序里采用自适应容差方案：先跑一小批，用接受样本距离的某个分位数作为下一轮ε，迭代收缩。

第二个是距离函数的权重。电压量纲是0.9到1.1的标幺值，支路功率可能是几百兆瓦，如果距离直接相加，后者会完全淹没前者。必须使用标准化权重，让每个观测量在距离里的贡献大致相同。这一点不处理好，ABC后验就会只盯着最大量纲的支路。

2.5 对比指标怎么设置才公平

程序对MCS和ABC的结果做比较时，没有只比均值。均值在概率分布里最容易对齐，但不能反映尾部风险。我设计了三层对比：

• 第一层是各节点电压均值和标准差，用来快速观察两套方法的整体偏差。
• 第二层是越限概率，也就是电压越上/下限或支路过载的概率，这是工程决策真正关心的量。
• 第三层是分位数误差，特别是5%分位数和95%分位数。这一层最苛刻，因为任一近似方法在分布尾部都可能失准。

计算误差时，以够多样本的MCS结果为参考基准，对于每个对比指标都算相对误差百分比。如果某个支路的95%分位误差超过5%，程序会在报告中红字标出。这比单纯给一张概率密度图叠加曲线更实用，因为你能精确定位到哪条支路需要担心。

3. 实操过程与关键环节实现

3.1 程序整体流程步骤

这套程序的外部入口是一个主脚本，逻辑上分成六个阶段。

1. 定义网络：读取节点、线路、变压器参数，构建导纳矩阵。
2. 定义随机源：给每个注入不确定性的节点配置分布类型和参数，同时给出观测数据文件路径。
3. 基准MCS计算：使用Sobol抽样生成输入矩阵，逐一解潮流，得到节点电压矩阵和支路潮流矩阵，计算统计指标。
4. ABC前置：从先验分布抽取候选参数，生成模拟观测，计算距离矩阵，自适应确定容差序列。
5. ABC主循环：迭代多轮，每轮用当前容差过滤候选参数，得到参数后验样本，再把这批参数后验映射回潮流结果分布。
6. 对比与输出：统计误差、生成图表、保存结果结构体。

以下几个函数是整个程序的主干：

function [V, converged] = runPF(net, sample) % net 为网络结构体，sample 为本轮各随机注入的有功无功 % 解潮流，返回节点电压复向量和收敛标志 % 内部使用牛顿-拉夫逊，PV节点越限后自动转PQ % 收敛判据：|dP| < 1e-6, |dQ| < 1e-6 end

function [stats] = calcStats(Vmat, branchP) % Vmat: nsample x nbus % branchP: nsample x nbranch % 返回每个节点的均值、标准差、5%分位、95%分位、越限概率 stats.mean = mean(Vmat, 1); stats.std = std(Vmat, 0, 1); stats.p5 = quantile(Vmat, 0.05, 1); stats.p95 = quantile(Vmat, 0.95, 1); stats.violation = mean(Vmat > vmax | Vmat < vmin, 1); end

function [postParams, postSamples] = runABC(prior, observed, net, options) % 自适应ABC，采用逐轮收缩容差 % prior: 先验分布数组 % observed: 观测数据 % options: 每轮样本数、初始容差、收缩因子、最大轮数 end

3.2 潮流计算函数的关键设计

runPF这个函数是整个程序的性能核心。它不接收外部工具箱依赖，节点数据和雅可比矩阵稀疏结构都在初始化时预计算好了。雅可比矩阵在牛顿法每轮迭代中需要更新，但稀疏模式固定，因此MATLAB可以用sparse预分配索引，避免反复调整内存。

一个容易忽略的细节是节点类型的动态变化。当某个PV节点的无功出力超过上限时，正确的做法是把该节点从PV转成PQ，并固定无功为上限值，下一轮使用新的节点类型继续迭代。很多自编程序在概率潮流场景下崩溃，就是因为这个转换逻辑没写对。

潮流收敛后，结果里还要顺便回传节点参与因子或者灵敏度信息吗？在纯概率框架下不需要，但我在做后处理时发现，用雅可比矩阵的逆可以估计局部灵敏度，对判断“哪个随机输入对电压波动贡献最大”非常有用。

3.3 ABC主循环怎么设计才稳

ABC主循环我写得比较谨慎，因为容差、接受率、后验质量之间存在三角矛盾。第一轮先用很宽的容差，比如距离矩阵的30%分位数，确保能拿到一批可用的参数样本。之后每轮把容差收缩到上一轮的60%，同时增加候选样本数，弥补被过滤掉的部分。

每一轮都会记录接受率和后验参数均值。接受率低于1%时，程序会主动停止收缩并进行一次诊断，因为这说明距离函数设计或先验范围有问题。先验范围太宽，接受率必然低；范围太窄，后验又会被先验边界截断。

for iter = 1:maxIter % 从当前后验（第一轮则是先验）抽样 candidates = drawFromProposal(prior, theta_last, options.propSigma); for k = 1:size(candidates,1) Ysim = quickSimulate(net, candidates(k,:), options.innerMC); dist(k) = weightedDistance(Ysim, observed, options.weights); end epsilon = quantile(dist, options.quantileKeep); accepted = dist <= epsilon; theta_last = candidates(accepted,:); if mean(accepted) < options.minAcceptRate break; end end

注意这里用了“快速模拟”而不是完整潮流解。在ABC前期容差宽时，完全跑2000次潮流太浪费。程序提供innerMC参数，默认值为150，后期再把innerMC逐步提到500，平衡速度和精度。

3.4 参数设置参考与运行时间分配

参数设置的合理范围取决于网络规模和硬件。我用某33节点配电网模型做基准测试时，一组常用参数配置可供参考：

运行时间方面，2000次潮流MCS在普通笔记本上大约需要几十秒到两三分钟，ABC的8轮迭代由于样本量大，通常是MCS时间的三到五倍。如果实在等不起，把支路距离权重里不关心的线路权重设为零，ABC耗时能明显降下来。

3.5 结果怎么量化对比

程序最后输出一个comparison结构体，包含每个节点电压的均值误差、标准差误差、P5/P95分位误差、越限概率绝对值误差。对每条支路，输出以上指标之外还加一个“尾部风险比”：ABC算出的P95减去MCS基准的P95，再除以基准P95。

这个尾部风险比是最直接的工程信号。如果某条支路比例超过10%，就说明ABC在当前参数设置下不适合描述这条支路的重载风险，要么增加innerMC，要么改用更充分的MCS对该支路单独复核。

图表方面不一定画很多花活，我建议至少输出三张：

• 节点电压均值误差散点图，横轴为节点编号，纵轴为ABC相对MCS的偏差百分比。
• 关键支路潮流概率密度叠加图，一条实线是MCS基准，一条虚线是ABC后验结果。
• 越限概率柱状图，对比两套方法的差异。

4. 常见问题与排查技巧实录

4.1 潮流求解在部分样本上不收敛

这是概率潮流最先暴露的问题。原因往往是特殊抽样组合让系统进入重载状态，比如风电全出、负荷低谷导致电压偏高。我用样本级重试加平启动的方式兜底，但如果你发现某批样本的失败率超过1%，一定不能忽视。

排查顺序：先固定不收敛样本，把输入变量打出来，手算一下该场景的无功平衡是否合理；然后看是不是光伏节点在低压穿越逻辑上的无功极限设置问题；最后看网络参数中是否存在零阻抗支路，这类支路会让雅可比矩阵奇异。

经验是：不收敛样本往往不是随机的，而是对应某个物理上不想运行的状态。这一步也叫“概率潮流里的极端事件识别”，值得单独分析。

4.2 样本数到底要多少

很多人一上来就设5000或者10000，其实没必要。确定样本量的方法是看指标稳定度：把MCS分成前1000、2000、3000三组分别统计P95，如果三组间的相对差小于1%，当前样本数就够了。

如果关注的是5%小概率事件，比如过载概率，样本量需要显著提高。简单估计，要稳定测出1%概率事件，至少需要数千个有效样本。程序里的默认2000只适合做方法对比，不是工程定级结论。真实项目中对越限概率做精细评估时，我一般把MCS提高到5000并配合Sobol序列。

4.3 ABC接受率太低或后验崩了

ABC最常见的失败模式是三轮迭代后接受率跌到0.1%以下。不要盲目加样本量，先检查三点：

• 先验范围是否过宽。有个项目把负荷均值先验设为±50%，结果接受率奇低，收窄到±10%后立刻正常。
• 距离权重是否失衡。如果有一条支路的功率量级是其他支路的一百倍，它会在距离函数里完全支配，导致其他约束形同虚设。重新加权是关键。
• 观测数据是否太“干净”。ABC的本质是反演，如果观测数据完全无噪声，任何近似距离都会显得巨大。给观测数据加适量量测噪声反而能让后验更健康。

另外，ABC轮间转移建议使用“以当前后验均值+比例缩放协方差为提议分布”，别用全空间独立抽样。前者在参数维度较高时效率高得多。

4.4 运行时间长到无法接受

这几乎是概率潮流必踩的坑。第一刀切在潮流求解本身，用稀疏矩阵、预分配雅可比、避免循环里反复调用inv。第二刀切在抽样效率，把普通随机抽样换成Sobol后，达到同等P95精度最多能省一半样本。第三刀切在ABC设计，在容差较大的前几轮降低innerMC，后期再加大。

如果还嫌慢，程序提供“灵敏度辅助抽样”：先用少量MCS样本估计各随机源对目标支路的灵敏度，对灵敏度小的输入降低抽样精度甚至可以暂时固定，这样能显著压缩有效维度。实测在某36维随机输入的算例中，这个技巧让总体耗时降到原来的40%左右。

4.5 两套方法结果对不上

对不上很正常，关键要看对不上在什么位置。均值层面对不上，几乎可以断定是ABC后验有偏，检查距离函数或先验。标准差层面有偏差，可能是ABC在收紧容差时把参数多样性压没了——试着增大容差收缩因子到0.7，不要收得太猛。

分位数和尾部对不上，则是ABC这种近似方法最真实的“能力边界”。出现这种情况不要急着怪程序，先看图：如果ABC的P95系统性低于MCS，说明接受样本距离最小化的方向没有覆盖到高风险的参数组合。这时候最有效的不是调ABC，而是改用两阶段抽样，先跑一次MCS识别高风险区，再用ABC在局部区域做后验分析。

5. 工程化落地与个人经验

5.1 从验证程序到生产工具的改造清单

这程序一开始是为了研究对比做的，但后续我发现它可以直接改造成日常分析工具。建议做几件事：

• 网络数据改为从外部文件读取，而不是硬编码在脚本里。这样换算例时不需要改任何逻辑代码。
• 随机源参数单独存成配置表，每个随机源有标签、分布类型、参数、所属节点、量测数据文件名。ABC要先验都会自动匹配。
• 加一份结果摘要导出功能，把越限概率、关键分位数、耗时等信息导出成表格。调度或报告场景非常需要。

有的公司环境里没有完整商业MATLAB授权，只能靠生成的独立程序跑批。这时需要把程序里的可视化部分拆出去，核心计算封装成函数，减少依赖工具箱的数量。代码尽量只依赖基础MATLAB函数，可以大大降低部署成本。

5.2 一些值得记住的经验提醒

第一，任何概率方法的结果都要与物理直觉对得上。曾经有一次算例显示某节点电压P95超过1.1，我第一反应是算法错了，后来查下来是模型里忽略了变压器的抽头调节范围，导致高电压场景被系统性地放大了。概率潮流不是纯数学游戏，网络模型的物理细节最终决定结果可信度。

第二，容差和样本量的组合不要手工瞎调。程序里虽然给了默认值，但换网络、换随机源后，最好先跑一次小规模预扫描，画出“接受率-容差”曲线，找到拐点再批量运行。这比我当初直接跑完整8轮后看到结果不合理而返工要省几个小时。

第三，MCS结果虽然可以当作基准，但它本身也有随机误差。做方法对比时，基准的样本量不能和对比方法“同样水平”，基准应当高出一个量级，否则你对比的其实是两个有噪声估计之间的噪声差。我通常让基准MCS样本数是ABC最终有效样本数的五倍以上。

最后再分享一个我后来才摸索出的技巧：观察ABC每轮参数后验的变化轨迹，而不是只盯着最终结果。如果后验均值随着容差收紧在同一个方向持续漂移，说明先验位置和观测数据之间存在系统性错位。这个漂移方向本身就是有价值的信息，它可以反过来告诉你，某个负荷节点的均值在历史预测里可能被偏置了。程序默认不输出这个轨迹，但我在多个项目里都会手动打开，非常值得看。