☰
配电网可靠性评估的序贯蒙特卡洛仿真:从10负荷点模型到工程实践
2026/10/13 11:37:11 网站建设 项目流程

简介:面向电力系统运维与科研人员,这是一套基于蒙特卡洛模拟的10负荷点配电网可靠性计算程序,能够通过大量随机故障-修复序列模拟小型配电网运行状态,计算SAIFI、SAIDI、ASAI等常用可靠性指标,适用于教学演示、方案验证与配电网可靠性评估。压缩包共44个文件,包含cpp/sln/vcxproj等C++工程文件、MATLAB的m脚本、xlsx配电网络参数表、pdf网络结构说明以及运行结果截图,整体仅1.07MB,结构清晰、便于直接运行与二次开发。目前已有1648人学习浏览,适合电力专业学生和配电网规划人员快速入门蒙特卡洛可靠性分析。通过实际10负荷点配电网算例,可掌握设备故障建模、负荷削减统计和可靠性指标计算的完整实现思路,并能修改参数扩展至更大规模网架,为电网规划与运行提供量化参考。

1. 配电网可靠性的蒙特卡洛法:10 个负荷点的小模型,为什么值得跑

蒙特卡洛法在配电网可靠性计算里一直处于“原理简单、落地麻烦”的位置:原理就是反复抽样模拟元件故障,落地却要跟随机数、故障遍历、收敛判断纠缠。这套基于 10 负荷点的蒙特卡洛可靠性计算程序,用序贯模拟把配电网元件的随机故障和修复过程放到时间轴上逐小时推进,直接输出 SAIFI、SAIDI、CAIDI、ASAI 和缺供电量等指标。对刚接触电网可靠性评估的工程师来说,10 负荷点恰好在“能看清逻辑”和“能算完”之间取了个平衡——算得快,还足够暴露配电网模型里最典型的坑。规划方案对比、供电可靠性评估、课程设计验证都用得上。

2. 程序架构与序贯蒙特卡洛原理:故障时间轴怎么叠起来

一套可靠性计算程序,写出来能跑通容易,跑出可信结果很难。按拆过同类程序的经验,这套 10 负荷点程序的核心模块可以分成四块:输入数据准备、元件状态抽样、故障影响分析、指标统计输出。前两块决定仿真逻辑是否严谨,后两块决定指标是否可信。序贯蒙特卡洛的本质是先把单个元件的“正常—故障—修复—正常”过程连成一条时间轴,再把所有元件的时间轴按时间戳合并成系统状态序列。理解这个合并过程,比背住输出指标的公式更重要。

2.1 序贯蒙特卡洛与状态枚举:为什么仿真更贴近实际

状态枚举法理论上更简洁:对每个元件假设只有正常、故障两态,把所有元件状态组合列出来,按概率加权求期望。网络规模小的时候,枚举数量还扛得住;当馈线段超过二十条,组合数量指数上升,枚举法就开始吃力。更麻烦的是,状态枚举法处理不了“今年负荷曲线随时间变化”这类时序问题,而可靠性指标里的缺供电量恰恰高度依赖停电发生时段的负荷水平。

序贯蒙特卡洛的做法是把问题从“枚举所有状态”变成“沿着时间轴往前走”。每年从 1 月 1 日 0 时开始,给每个元件抽样出下一次故障时刻,再抽样出修复持续时长,不断推进,直到年末。负荷点在某次故障期间经历多长时间的停电,哪几个小时停电,都被明确记录。相比状态枚举法直接用期望值累计,序贯蒙特卡洛多出来的计算量换来了两个好处:一是可以接入逐小时负荷曲线,缺供电量算得更准;二是能统计停电次数的分布特征,而不只是一个均值。

两种方法在特定场景下的选择,可以简单看成下面这张表的取舍:

对比项状态枚举法序贯蒙特卡洛
核心输入元件故障概率元件故障率、修复时间、负荷时序
处理时序负荷困难天然支持
输出指标期望值为主期望值、分布、分位数
典型规模组合爆炸风险计算量与仿真年数线性相关
收敛性判断枚举完全则精确需要方差系数或置信区间

这张表只针对辐射网配电网场景。带分布式电源或储能参与转供时,状态枚举法要做时段状态组合,复杂度更高,序贯仿真的优势会更明显。这套程序选择序贯蒙特卡洛,并不是因为状态枚举不好用,而是 10 负荷点恰好能体现时序仿真价值,又不会等太久——一个常规 5000 年仿真在普通办公电脑上也就几分钟量级。

2.2 元件寿命抽样与修复时间建模:指数分布不是唯一选项

序贯蒙特卡洛对每个元件要反复回答两个问题:什么时候坏,坏了修多久。第一个问题的标准答案是按指数分布抽无故障工作时间 TTF,公式是:

TTF = -ln(U) / λ

其中 U 是均匀分布在 (0,1) 上的随机数,λ 是元件失效率,单位必须与 TTF 一致。如果分支表里写的故障率是“次/年”,那么换算是 λ_h = λ_year / 8760。这段代码负责生成单个元件的故障事件序列:

% 生成单条馈线段的一年事件序列 function events = gen_events(lambda_year, r_mean, r_std, seed) rng(seed, 'twister'); lambda_h = lambda_year / 8760; % 年故障率换算为小时 events = []; t = 0; while t < 8760 ttf = -log(rand) / lambda_h; % 指数分布抽样 t = t + ttf; if t >= 8760 break; end % 修复时间按对数正态抽样,避免固定值过于理想 mu_log = log(r_mean^2 / sqrt(r_std^2 + r_mean^2)); sigma_log = sqrt(log(1 + r_std^2 / r_mean^2)); ttr = lognrnd(mu_log, sigma_log); events = [events; t, ttr]; t = t + ttr; end end

这段代码有两个参数需要根据实际数据调整。r_mean 是平均修复时间,工程上通常来自检修记录,单位是小时;r_std 是修复时间的标准差,如果原始数据没给,常见做法是取 r_mean 的 30%~50% 作为初值,做敏感性分析时再上下浮动。修复时间用对数正态分布是为了模拟少数几次长时间停电的肥尾效应,这类长停电对 SAIDI 的影响远大于对 SAIFI 的影响,所以在评估“平均停电时间”时不能忽略。

随机数种子 seed 的处理是这段代码里最容易翻车的地方。rng(seed, 'twister') 放在函数开头,意味着每次传入相同 seed 都会得到完全相同的事件序列,这对调试很有用。若把种子设成年份相关的值,比如 year*100+branch_id,那么每个仿真年、每个元件的序列都不同,但整体可复现。千万别在 while 循环里反复调用 rng,那会破坏随机序列的独立性,结果就是指标振荡得莫名其妙。

2.3 系统事件队列与故障影响分析:辐射网的断电边界怎么划定

每个元件生成独立事件序列后,程序要把它们合并成一条按时间排序的系统事件表。事件表里每条记录至少包含四个字段:发生时刻、元件编号、是故障还是修复完成、连接关系。合并时用队列向前推进,每次取最早事件,更新受影响负荷点。这就是序贯仿真的“系统状态推进”。

时刻(h)元件编号事件类型说明
312.52故障开始支路 2-3 故障,负荷点 3 失电
317.22修复完成支路 2-3 恢复,负荷点 3 复电
2108.74故障开始支路 4-5 故障,负荷点 5、7、9 失电

接下来要对辐射网做故障影响分析。一条馈线故障时,其下游所有负荷点全部失电,上游负荷点不受影响;如果有备用电源转供,失电负荷点在转供容量允许时会部分恢复,这套基础程序假设是纯辐射无联络,所以所有下游负荷点都只能等故障元件修复。

% 辐射网故障影响划分 function affected = find_affected_loads(failed_branch, branch_data, load_data) from = branch_data(failed_branch, 1); to = branch_data(failed_branch, 2); % 从故障末端出发,向下游递归收集负荷点 affected = []; stack = to; while ~isempty(stack) node = stack(end); stack(end) = []; if ismember(node, load_data(:,1)) affected(end+1) = node; % 该节点是负荷点 end next_nodes = branch_data(branch_data(:,1) == node, 2); stack = [stack; next_nodes(:)]; end end

这段代码用栈实现下游递归搜索。核心逻辑是从故障支路末端节点出发,沿父到子的方向把所有子节点收进 affected 数组。要特别注意别把故障支路的起点节点也算进去——起点到电源侧是上游,不受这一条支路故障影响;而终点及其下游才停电。很多初版程序在这写反,导致可靠性指标偏低或偏高。

合并后的系统时间轴按“故障开始”和“修复完成”两类事件驱动。故障开始事件记录断电开始时间,但不立即结算停电时长;修复完成事件到来时才用当前时刻减去上次状态切换时刻,得出停电持续时间。这样处理的好处是避免在多次故障叠加时重复计时。注意事项:如果两个不同元件先后故障且影响同一个负荷点,要在状态表里维护该负荷点是否持续停电;第二个故障发生时停电时长不重置,仍从第一次断电开始计算,直到任一故障修复。这个细节不处理好,缺供电量很容易低估。

3. 输入数据与仿真主循环:从馈线参数到五项可靠性指标

3.1 输入数据结构:节点、馈线、负荷点三张表怎么对齐

任何可靠性程序跑出来的指标,可信度都不会高于输入数据。这套程序的三张核心输入表分别是节点负荷表、馈线参数表和可选负荷曲线表。节点负荷表必须包含节点编号、父节点编号、用户数、平均负荷;馈线参数表包含起点、终点、故障率、平均修复时间。三张表之间的索引关系用节点编号关联,而不允许用行号硬关联,这是后续排查错位问题的关键前提。

下面这组输入是 10 负荷点基础算例的典型组织方式,故障率用“次/年”,修复时间用“小时”:

% 节点负荷表:节点号, 父节点, 用户数, 平均负荷(kW) load_data = [ 1, 0, 120, 350; 2, 1, 80, 220; 3, 2, 95, 280; 4, 1, 60, 180; 5, 4, 75, 210; 6, 4, 90, 260; 7, 5, 50, 140; 8, 6, 70, 200; 9, 7, 45, 130; 10, 8, 85, 240 ]; % 馈线参数表:起点, 终点, 故障率(次/年), 修复时间(h) branch_data = [ 1, 2, 0.10, 5; 2, 3, 0.08, 4; 1, 4, 0.12, 6; 4, 5, 0.09, 5; 4, 6, 0.07, 3; 5, 7, 0.06, 4; 6, 8, 0.08, 4; 7, 9, 0.05, 3; 8, 10, 0.06, 4 ];

负荷点必须同时是节点,节点 1 是电源母线,不挂负荷,所以 load_data 第一行用户数可以设成 0,方便后续统计总用户数时用 sum 过滤。用户数这一列不仅用于 SAIFI 和 SAIDI 的分母,还隐含了“同一节点多用户”的权重。平均负荷列单独用于计算 ENS,如果手头有时序负荷,可以换成 8760 行向量,并让仿真程序按小时索引。

参数设置的另一个容易被忽略点是馈线故障率的量级。主馈线段故障率常取每回线每年 0.1 次左右,分支线略低。实际数据如果是架空线,还会受雷害、外力破坏影响,同一条馈线不同年份波动很大,所以输入数据里不要追求“精确到一个点”,而应该把故障率定义成多年平均值,并在结果输出时标明这是长期均值。

3.2 仿真主循环与可靠性指标统计

主循环的逻辑是:每年年初为所有馈线段生成事件序列,合并后按时间顺序推进,年末统计出该年各负荷点的停电次数和停电时长,然后累加到总样本里。仿真 N 年后,用总累加值除以 N 得到期望指标。这样每一年的样本都独立,方差估计才能成立。

% 主循环:逐年仿真并累计指标 N_year = 5000; total_int = zeros(10,1); total_dur = zeros(10,1); saidi_year = zeros(N_year,1); for year = 1:N_year sys_events = []; for b = 1:size(branch_data,1) ev = gen_events(branch_data(b,3), branch_data(b,4), ... branch_data(b,4)*0.4, year*100+b); sys_events = [sys_events; ev, b*ones(size(ev,1),1)]; end sys_events = sortrows(sys_events, 1); [int_year, dur_year] = simulate_year(sys_events, load_data, branch_data); total_int = total_int + int_year; total_dur = total_dur + dur_year; saidi_year(year) = sum(dur_year) / sum(load_data(:,3)); end % 系统指标 SAIFI = sum(total_int) / sum(load_data(:,3)) / N_year; SAIDI = sum(total_dur) / sum(load_data(:,3)) / N_year; CAIDI = SAIDI / SAIFI; ASAI = (1 - SAIDI / 8760) * 100; ENS = sum(total_dur .* load_data(:,4)) / N_year;

逐年样本之所以能用于收敛判断,是因为每年的事件独立生成、每一年都从零初始化,重复 N 年后可以无偏估计期望值。这里 ENS 用了平均负荷代替逐小时负荷,在负荷曲线平稳的算例里足够;一旦负荷变化剧烈,这种简化会把停电发生在高峰时段的损失低估——有负荷曲线就改成逐年累加时按事件时间索引负荷值。

SAIFI、SAIDI 是系统级指标,计算顺序必须是先各负荷点累加,再除以总用户数,不能先求负荷点平均再算平均,那会丢掉用户数权重。指标口径也要提前统一:

指标计算式含义
SAIFIΣ用户断电次数 / Σ用户数平均每个用户每年停电次数
SAIDIΣ用户停电时间 / Σ用户数平均每个用户每年停电小时数
CAIDISAIDI / SAIFI平均每次停电持续小时数
ASAI(1 - SAIDI/8760) × 100%供电可用率
ENSΣ(负荷点负荷 × 停电时长)年缺供电量

simulate_year 内部维护一个当前各负荷点是否有电的标志数组。每次处理“故障开始”事件时把受影响负荷点标成停电并累加一次次数;处理“修复完成”事件时用当前时间减本段停电起始时间,累加到停电时长。如果程序在故障影响分析时引入了联络转供逻辑,这段恢复判断要拆成“可转供立即恢复”和“等原馈线修复恢复”两种路径。

3.3 收敛判据:仿真年数不是越大越好

仿真年数选多少,不能靠感觉。工程里常用方差系数 β 来判断,定义是样本标准差除样本均值再除根号下样本数,β ≤ 0.05 时认为均值估计稳定。SAIFI、SAIDI 的收敛速度不同,系统级和负荷点级的收敛速度也差很多。一般低压分支线的负荷点,因为故障次数少,收敛明显慢于系统级指标。

% 滑动窗口收敛检查 block = 500; for year = 1:N_year % 已把当年SAIDI存入 saidi_year(year) if mod(year, block) == 0 mu = mean(saidi_year(1:year)); sd = std(saidi_year(1:year)); beta = sd / (mu * sqrt(year)); fprintf('%d年 beta=%.4f SAIDI=%.4f\n', year, beta, mu); if beta < 0.05 break; end end end

这个判据会对长尾分布比较严格:少数年份出现一次长停电,会拉大标准差,让 β 迟迟降不下来。遇到这种情况,别急着加年份,先算一算这一回长停电是由哪个元件、哪次故障引起的,确认它不是输入参数错误后,再判断是否接受更长收敛时间。仿真年数上限要保留——比如设 10000 年,到上限也打印“未收敛”警告,而不是静默结束。默认年限情况在报告里也要写清楚,因为可靠性结果的有效性和样本量直接相关。

4. 避坑与常见问题排查:随机数、编号错位与修复时间的四个典型坑

4.1 随机数种子:换一次种子结果差三成,怎么定位

现象:同一套参数,第一次跑 SAIDI 是 1.82 小时,第二次变成 2.36 小时,差距 30% 以上。

原因:两个典型因素叠加——随机数种子没固定,仿真年数又偏短。序贯蒙特卡洛本身就是随机模拟,种子不同意味着元件的故障年份、修复时长全部重排,短样本下结果波动大是正常的。若是仿真年数已经五千以上还波动这么大,就要怀疑代码里出现了非受控随机源。

解决:先把 rng(固定种子) 放在程序最前面,保证任意两次运行完全一致,这能排除随机性干扰定位逻辑错误。再按 β<0.05 的要求加大仿真年数;通常 10 负荷点网络把种子固定后 SAIDI 的变化仍在百分之几以内才可信。最后检查循环内是否误用 rng 或 randn 重置了随机流——常见的翻车点是某个子函数里写了 rng('shuffle'),这种代码会让随机序列产生重复模式,必须移除。

4.2 负荷点编号错位:指标张冠李戴

现象:节点 7 和节点 9 的停电次数几乎相同,但节点 9 在支路末端,上游任何故障都会影响它,而节点 7 只受一条特定馈线影响,二者统计值不该这么接近。

原因:程序在故障影响分析时用了数组下标而不是节点编号来判断负荷点。例如 load_data 的第几行被直接当作节点号,一旦节点号跳号,所有映射全错。10 负荷点连续编号时问题不明显,改到其他网络立刻暴雷。

解决:强制建 lookup 表,用节点编号查索引。比如 id2idx = containers.Map(load_data(:,1), 1:size(load_data,1)),所有影响判断都先经查询再操作。改完做交叉验证:逐条断开每个馈线段,手工画图推演受影响负荷点,与程序输出对比,偏差为零才算通过。这个校验花不了十分钟,但对后面扩展网络结构非常重要。

4.3 修复时间分布:固定均值导致 SAIDI 偏小

现象:程序输出的停电次数和实际统计接近,但停电时长普遍低两成,报告没法用。

原因:固定修复时间抹掉了修复过程的随机性。配电网的修复时长往往右偏严重,少数几次复杂故障可能断电十几个小时,而固定均值模型最多给到几小时,长尾没模拟出来。

解决:把修复时间改成对数正态分布抽样,均值按历史数据,标准差按均值的 30%~50% 起步。改完对比新旧结果,通常 SAIDI 和 ENS 上升、SAIFI 不变。若差异过大,说明原始均值数据里没有包括极端天气事件,需要重新整理统计数据而不是盲目调分布参数。

% 修复时间抽样:固定值 vs 对数正态 mr = 5; sr = 2.5; mu_log = log(mr^2 / sqrt(sr^2 + mr^2)); sigma_log = sqrt(log(1 + sr^2 / mr^2)); ttr = lognrnd(mu_log, sigma_log);

4.4 仿真年数不足:关键元件故障次数太少

现象:仿真 1000 年结果 SAIFI 0.31,加到 8000 年变 0.37,变幅超过 15%,怎么看都不像收敛。

原因:馈线段年故障率只有 0.05~0.1 次时,1000 年样本里某条支路只故障 50~100 次,抽样波动非常大。这类低频元件恰恰可能位于多个负荷点上游,对系统指标影响不小。

解决:不要只看系统指标收敛,要看关键支路累计故障次数。代码里统计每条馈线年故障总数,小于 200 的在输出结果中标注警告。对年故障率低于 0.05 的支路,可以把该参数适当提高 10%~20% 做敏感性测试,看指标变化幅度,再决定要不要继续增加仿真年份。

4.5 指标与历史统计对不上:先查口径再查算法

现象:仿真出来的 SAIFI 比公开统计值高一倍,第一反应是算法写错了。

原因:统计口径经常不一致。仿真里的“用户”按负荷点用户数累加,统计资料里可能按变压器台区折算;仿真输入故障率如果来自“每公里故障率 × 线路长度”,与统计值默认分母不同,结果自然对不上。这不是算法问题,是口径问题。

解决:先逐项核对故障率单位、用户数定义、修复时间统计范围。建议单独写一段数据口径说明放到报告里,包括故障率统计年份、是否包含极端天气、修复时间是否含故障定位时长。口径对完,再做第 5 章的三重校验,从逻辑层面确认程序没写错。

5. 结果验证与扩展技巧:把 10 负荷点模型推向多馈线联络结构

5.1 结果验证:三重校验法

拿到计算结果后先别急着写报告。第一重校验是单年样本分布:把五千个年 SAIDI 样本画直方图,应呈现明显的长尾,若所有年份几乎一样平,多半是抽样逻辑把故障“漏吃”了。第二重是元器件手工推演:挑一条分支馈线,手算 1000 年期望故障次数,和程序统计对比,偏差超过两个标准差就检查事件合并逻辑。第三重是极简对拍:把所有修复时间设为同值,跑一次状态枚举法求系统期望指标,和蒙特卡洛结果对比,两者差异在 5% 以内即可认为程序忠实于输入参数。

5.2 扩展为带联络的多馈线网络

扩展时需要改两处:故障影响分析不再默认全部下游停电,而是先检查联络开关是否可用、转供容量是否足够;事件恢复逻辑把停电时长拆成“可转供停电时段”和“仍失电时段”。核心新增逻辑是一个转供判断函数:

function ok = can_transfer(lp, tie_status, capacity_remain, load_p) if tie_status(lp) == 0 ok = false; % 联络开关未闭合 elseif capacity_remain < load_p(lp) ok = false; % 余量不够,不可转供 else ok = true; end end

转供容量判断是接入分布式电源或储能时要扩充的重点。若联络馈线本身故障率不低,还要把它作为可失效元件加入事件表;若其故障率比主馈线低两个数量级,通常可以省略。扩展完成后,建议拿原辐射网结果作基准,逐条加联络线对比可靠性提升幅度,既验证程序又便于方案分析。

从那次把随机数固定、编号映射重构、修复时间改对数正态之后,我每次换网络做可靠性评估,都会先把这套三重校验走一遍,再谈优化和扩展。希望帮到你。

本文还有配套的精品资源,点击获取

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

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

立即咨询