简介：面向电力系统可靠性分析初学者以及配电网运维、规划人员的蒙特卡洛仿真程序包，围绕含10个负荷点的配电网展开，通过随机抽样模拟设备故障与修复过程，统计停电时间、停电次数等数据，进而估算SAIDI、SAIFI、ASAI等可靠性指标。程序内置故障模型与恢复策略设定，可结合随附的网络结构与参数表复现算例，帮助理解在复杂网络存在不确定性的情况下如何开展可靠性评估。

压缩包共44个文件、约1.07MB，以Visual Studio工程文件（含sln解决方案、vcxproj项目、cpp源码以及debug目录下的pdb、obj、tlog等编译中间件）和MATLAB脚本为主体，另附xlsx参数表、pdf网络结构说明与程序运行结果截图，源码既可直接编译运行，也便于对照修改。

已有1651人学习下载。读者可据此搭建一套小型配电网可靠性评估流程，理解抽样次数与结果收敛性的关系，并对比不同故障率、修复时间下的指标变化，为规划与风险分析提供参照。

1. 10 个负荷点、8760 小时，这台机器能算到多细

10 个负荷点、8760 小时、几千次状态事件——这套配电网蒙特卡洛可靠性计算程序跑完一遍，给的不只是几个指标数字，还有每个数字背后的方差系数 β。手算解析法的结果和现场台账对不上，这事在可靠性评估里太常见——差的地方不在公式，而在假设：开关操作时间被抹成零，转供被当成百分之百成功，计划检修压根没算进去。解析法给的是期望值，蒙特卡洛给的是带方差的抽样分布，后者能告诉你“这个 SAIDI 是 3.2，还是 3.2±0.3”。

程序面向 10 负荷点的中压辐射状馈线，把线路区段、配变、熔断器、分段开关按两状态马尔可夫模型建模，输出负荷点级的平均停电频率 λ 和年均停电时长 U，再汇总成系统级的 SAIFI、SAIDI、CAIDI、ASAI、ENS。要出可靠性台账的规划人员、做故障影响分析的运维、课程设计里需要一条收敛曲线的同学，都能直接上手。

2. 序贯抽样怎么建模：把 10 个负荷点拆成区段、节点、负荷

2.1 拓扑抽象：三级结构决定故障影响范围

10 个负荷点的馈线看起来简单，真正写代码时会发现麻烦全在“故障发生在哪一段、这段到负荷点之间有没有可隔离的开关”上。常见的做法是把网络拆成三级：节点、区段、负荷点。节点只标记位置（母线、分支点、负荷点接入点），区段是两端节点之间的一段线路，负荷点挂在某个节点上、通过一台配变取电。区段上再挂一个属性，说明这一段的故障靠什么手段隔离。

这里的选型理由很直接：负荷点的“最小路”不手工填，而是从拓扑里用广度优先搜索自动算出来。手工填路径表在 10 个负荷点时看着省事，一旦分支增加到二三十条，改一次拓扑就要重填一遍表，不出三次就会填错。自动算最小路多写了十来行代码，换来的是拓扑可改、可校验。

2.2 两状态模型与指数抽样：为什么不能图省事用正态分布

配电网可靠性评估里最主流的元件模型是两状态（运行、修复）马尔可夫模型。元件在运行状态停留的时间服从参数为 λ 的指数分布，在修复状态停留的时间服从参数为 μ=1/r 的指数分布，两个状态交替，互不相关。指数分布的关键性质是无记忆性：元件已经无故障运行了多久，不影响它下一秒的故障概率，这正好对应“恒定故障率”这个工程假设。

序贯蒙特卡洛的做法是把所有元件的“当前状态剩余时间”放在一张表里，每次取最小值推进，谁的时间先到谁就切换状态，其余元件的剩余时间同步扣减。这样时间轴是一条真实推进的线，而不是一堆彼此独立的抽样拼起来。

import numpy as np def sample_ttf(lam_per_yr: float, rng: np.random.Generator) -> float: """抽样下次故障前的运行时间，返回单位：年""" u = rng.random() # U(0,1) return -np.log(u) / lam_per_yr # 指数分布逆变换，均值 1/λ def sample_ttr(repair_hr: float, rng: np.random.Generator) -> float: """抽样修复时间，返回单位：年""" u = rng.random() return -np.log(u) * repair_hr / 8760.0 # 均值 = repair_hr 小时

sample_ttf的 λ 单位是次/年，传入 0.096 表示这条区段平均十年故障一次，抽样结果直接是年，和主循环的时间轴单位一致，省掉到处换算的麻烦。sample_ttr接收小时、返回年，是为了让两个函数返回值能直接比较大小；如果修复时间也用年作输入，参数表里填 0.00057 这种数，录入时手指一抖就是数量级错误。

提示：逆变换法用的是-ln(U)，U 取不到 0，但rng.random()理论上可能返回恰好 0。稳妥的做法是写u = max(rng.random(), 1e-12)，否则np.log(0)会给一个inf，这个inf会一路传染到min()里，最后表现为“程序跑完了但全是 0”，且不报任何错。

说到抽样，另一条血泪经验是：修复时间别用正态分布。正态分布左侧无界，抽样出负数修复时间时程序不会报错，停电时长会变成负的，SAIDI 被这个负值拉低，指标看起来漂亮得不像话。指数分布天然非负，均值可控，和恒定修复率的假设也一致，没有理由换。

2.3 指标口径先统一：五个指标各自的分子分母

指标算错，八成不是算法错，是口径没对齐。同样是 SAIDI，有的算含计划检修，有的不含；有的按用户数加权，有的按配变容量加权；有的把瞬时停电（小于 3 分钟）也算进去。程序里把这几个口径全部显式化，写在配置里而不是散落在代码里。

CAIDI 的分母是 SAIFI，如果某个负荷点全年一次都没停过，它的 SAIFI 是 0，做除法就炸。程序里统一用if saifi > 0 else 0.0兜底，但这个 0 本身有歧义——“没停电所以没有平均停电时长”和“停电时长为 0”是两回事，输出表里要单独标一个n_outage=0的标记，否则后期做统计的人会把这批样本当成有效数据算进去。

还有一条容易被忽略：ENS 里的负荷取值。用年平均负荷算出来的是能量不足期望，用峰值负荷算出来的是容量不足期望，两者能差出三成。程序默认读avg_load_kw，如果输入表里填的是峰值，指标名就得改，不能还叫 ENS。

3. 跑通第一个算例：从三个 CSV 到收敛输出

3.1 环境与目录：三分钟把骨架搭起来

依赖很轻，numpy 做抽样和统计，pandas 读写配置，tqdm 只是看进度条用的，不想装可以删掉对应的 import。

python -m venv .venv source .venv/bin/activate # Windows 下是 .venv\Scripts\activate pip install numpy pandas tqdm

目录按“配置、代码、输出”三分，输出目录不要和代码目录混在一起，否则清空结果时容易手滑。

mc_reliability/ ├── data/ │ ├── components.csv # 区段与配变：故障率、修复时间、隔离方式 │ └── loadpoints.csv # 10 个负荷点：用户数、平均负荷、有无转供 ├── src/ │ ├── model.py # 拓扑解析、最小路搜索 │ ├── sample.py # 抽样与一年内的状态推进 │ ├── metrics.py # 指标累加、方差系数 │ └── run.py # 入口，处理命令行参数 └── out/ # 运行结果，每次运行前清空

3.2 输入文件写法：字段、量纲、以及最容易填错的那一列

components.csv一行一个元件，type区分线路区段和配变。配变的 from 和 to 填同一个节点，length_km填 0，这样解析逻辑不用写两套。

id,type,from,to,length_km,lambda_per_yr,repair_hr,switch_hr,isolation S1,line,BUS,SW1,1.20,0.096,5.0,0.5,sectionalizer S2,line,SW1,SW2,0.85,0.068,5.0,0.5,sectionalizer S3,line,SW2,SW3,0.60,0.048,5.0,0.5,sectionalizer B1,line,SW1,N1,0.30,0.024,4.0,0.5,fuse B2,line,SW2,N4,0.25,0.020,4.0,0.5,fuse T1,transformer,N1,N1,0.00,0.015,12.0,0.0,none

lambda_per_yr是整段线路的年故障率，不是每公里的故障率。这是最容易翻车的一列：有人按行业手册填 0.08 次/(km·年)，忘了乘length_km，整条馈线的故障率被低估一个数量级，SAIDI 直接掉到 0.3 以下，看起来像“这个网架特别可靠”。程序里加了一条断言：如果所有区段的 λ 都小于 0.01 且length_km全部大于 0，打印一条警告，让填表的人自己确认单位。

loadpoints.csv里has_backup是 0/1 标志，表示这个负荷点有没有可转供的联络路径。它的作用在下一章的停电判定里体现：有转供的负荷点，主干区段故障后停电时长按开关操作时间计，没转供的只能等修复。

lp_id,node,users,avg_load_kw,has_backup LP1,N1,120,320,0 LP2,N2,85,210,1 LP3,N3,60,150,1

users直接决定 SAIFI/SAIDI 的权重。如果一个负荷点挂的是整台配变供几十户，另一个挂的是商业用户按容量折算，权重口径要在表里统一，别一半填户数一半填 kVA。

3.3 运行命令与参数：max-years 和 beta-target 要一起设

python -m src.run \ --data data \ --out out \ --max-years 20000 \ --beta-target 0.05 \ --seed 20240601 \ --beta-metric saidi

两个停止条件的关系要说清楚：max-years是硬上限，beta-target是质量目标，谁先满足谁停。最常见的误用是只设max-years，跑到上限就出结果，β 还在 0.2 上下，这种结果自己跟自己第二次跑都对不上。

3.4 输出三张表：怎么读、先看哪一张

运行结束后out/下有三个文件。lp_index.csv是 10 个负荷点各一行，含 λ、U、SAIFI 贡献、ENS 贡献；system_index.csv只有一行，是系统级汇总；convergence.csv是每年的 β 变化轨迹。

import pandas as pd sys_idx = pd.read_csv("out/system_index.csv").iloc[0] print(f"SAIFI={sys_idx.saifi:.3f} 次/户·年 " f"SAIDI={sys_idx.saidi:.3f} 小时/户·年 " f"ASAI={sys_idx.asai:.5f} beta={sys_idx.beta:.4f} years={int(sys_idx.years)}")

先看system_index.csv里的years和beta，再看指标值。如果years等于max-years而beta大于目标，说明样本不够，指标只能当参考值用。接着看lp_index.csv里各负荷点的 U 分布——如果某个负荷点的 U 比其他高出一个量级，多半是它的最小路上挂了一个isolation=none的区段，故障后只能等修复，这是拓扑本身的性质，不是 bug。

4. 把代码拆开：抽样、判定、累加三个动作

4.1 下次事件推进：一年之内的状态循环

主循环是整份代码的核心，逻辑是：每轮取所有元件里“剩余时间”最小的那个，把时钟推进到它到期，然后切换这个元件的状态。故障状态下，顺带统计它对每个负荷点的影响。

import numpy as np def simulate_one_year(net, rng): """仿真一年，返回 (每个负荷点停电次数, 每个负荷点停电时长/小时)""" # remaining: 元件在当前状态下的剩余时间，单位年 remaining = {cid: sample_ttf(c.lam_per_yr, rng) for cid, c in net.components.items()} state = {cid: "up" for cid in net.components} t = 0.0 cnt = dict.fromkeys(net.loadpoints, 0) dur = dict.fromkeys(net.loadpoints, 0.0) while True: cid = min(remaining, key=remaining.get) dt = remaining[cid] if t + dt >= 1.0: # 事件落在年末之后，本年结束 break t += dt for k in remaining: # 所有元件时钟同步推进 remaining[k] -= dt if state[cid] == "up": # 故障事件 ttr = sample_ttr(net.components[cid].repair_hr, rng) state[cid] = "repair" remaining[cid] = ttr for lp in net.loadpoints: # 逐个负荷点判定是否停电 h = lp_outage_duration(lp, cid, net, ttr * 8760.0) if h > 0: cnt[lp] += 1 dur[lp] += h else: # 修复完成事件 state[cid] = "up" remaining[cid] = sample_ttf(net.components[cid].lam_per_yr, rng) return cnt, dur

t和remaining都以年为单位，t + dt >= 1.0这一句是年度边界，越过就退出，不跨年。这样处理的好处是每个样本年独立，方便做逐年方差统计；代价是年度交界处正在修复的元件状态被截断，对年均指标的影响在千分之一量级，10 个负荷点的算例上可以忽略，但写论文时要在假设里说清楚。

lp_outage_duration返回的是小时，而ttr是年，所以乘了 8760。这一处单位换算是整个代码里最容易写错的地方，改代码时如果只改了一个方向，出来的 SAIDI 会离谱到让人怀疑人生。

4.2 停电判定：最小路加上隔离手段，一共四类情况

判定函数要回答的问题只有一个：元件cid故障后，负荷点lp停多久电。答案取决于两件事——故障是不是在lp的最小路上，以及这段路上有没有能把它隔离掉的开关。

def lp_outage_duration(lp, faulted_cid, net, repair_hr): """返回负荷点因该元件故障造成的停电时长（小时），0 表示不停电""" if faulted_cid not in net.lp_path[lp]: return 0.0 # 不在最小路上，供电不受影响 comp = net.components[faulted_cid] if comp.isolation == "fuse": return repair_hr # 熔断器下游故障，只影响本分支，按修复时间 if comp.isolation == "sectionalizer" and net.lp_backup[lp]: return comp.switch_hr # 分段开关隔离 + 可转供，按开关操作时间 return repair_hr # 无可隔离手段，等修复

net.lp_path在模型初始化时用 BFS 算好，是每个负荷点到电源的区段 ID 集合。net.lp_backup直接来自输入表的has_backup。

这张表就是整套程序的“物理模型”，代码只是它的翻译。改可靠性水平时，改的其实是这张表——比如给某条分支加个熔断器，就是把对应区段的isolation从none改成fuse，重跑一遍看 SAIDI 降多少。

有一处简化要交代清楚：熔断器动作时间和分段开关的倒闸时间在这里是固定值，不参与抽样。实际倒闸时间受操作队到达时间影响，波动不小，要模拟这个波动，把switch_hr也做成指数抽样或者对数正态抽样即可，改动量只有一行。

4.3 指标累加与方差系数：用在线更新省掉一个大数组

最常见的写法是把每年的 SAIDI 追加到一个 list 里，最后统一算均值和方差。10 个负荷点、两万年的算例下这点内存不算什么，但如果负荷点加到几百个、年数加到十万，list 的内存增长就很可观了。用 Welford 在线更新，常数级内存就能同时维护均值和方差。

class Welford: """在线更新均值与方差，避免保存全部样本""" def __init__(self): self.n = 0 self.mean = 0.0 self.m2 = 0.0 def update(self, x): self.n += 1 d = x - self.mean self.mean += d / self.n self.m2 += d * (x - self.mean) @property def var(self): return self.m2 / (self.n - 1) if self.n > 1 else 0.0 @property def beta(self): """方差系数 β = σ / (μ·√n)""" if self.n < 2 or self.mean == 0: return float("inf") return (self.var ** 0.5) / (self.mean * self.n ** 0.5)

beta这个定义要记牢：它衡量的是“均值的相对精度”，不是样本的相对离散度。β=0.05 意味着均值的 95% 置信区间大约在 ±10% 的均值范围内，样本量翻四倍才能让 β 减半，这也是为什么低故障率网络的仿真年数动辄上万。

收敛判断放在每一年的末尾，用beta-target和max-years两个条件做或运算。判断用的指标建议固定成一个，比如 SAIDI，别今天看 SAIFI 明天看 ENS——不同指标的方差差好几倍，混着判会让停止点变得不可解释。

5. 避坑与排查：五个把结果带偏的真实记录

5.1 收敛曲线反复横跳，β 卡在 0.08 下不去

现象：跑到五千年时 β 还在 0.08 附近来回跳，加年数也不见明显下降，指标值在 3.0 到 3.3 之间晃。

原因：判收敛用的是未加权的系统平均停电时长，而影响最大的那几个负荷点（用户数多、又没有转供）本身的停电事件就稀疏，单年样本的方差极大。样本量 n 与 β 的关系是 β ∝ 1/√n，从 0.08 降到 0.05 需要样本量翻 2.5 倍，不是加个一两千年能解决的。

解决：先把max-years提到 50000 跑一次，确认 β 确实随 n 下降，排除代码里 β 用了固定分母的问题；再确认判收敛指标是否是加权后的 SAIDI。如果是低故障率网络（全线年故障次数小于 0.5），改用分层抽样，把主干区段单独放大抽样频率，能把这个算例所需的年数压到五分之一。

5.2 SAIDI 比手算小一半，且找不到哪里错

现象：解析法手算的 SAIDI 是 5.8 小时/户·年，程序跑出来 2.9，正好一半，看起来像哪里漏算了一半负荷点。

原因：lp_outage_duration对主干区段故障一律按switch_hr（0.5 小时）返回，没有判断负荷点有没有转供路径。10 个负荷点里有 5 个has_backup=0，它们本该按修复时间（5 小时）计，被错误地按 0.5 小时计，恰好把总停电时长拉低了一半。

解决：在判定函数里加and net.lp_backup[lp]条件，并在输出里加一列backup_flag，把有转供和无转供的负荷点分开统计。改完之后再用手算的串联系统上界对一遍：把所有开关时间置 0、所有负荷点设为无转供，程序结果应该逼近解析法的纯串联结果，误差在 β 量级以内。

5.3 所有负荷点的停电次数都是 0，但程序正常退出

现象：三个输出文件都生成了，lp_index.csv里 λ 列全 0，system_index.csv里 SAIFI=0，程序没报任何错。

原因：lambda_per_yr那一列填的是每公里故障率（0.08），但代码当成整段年故障率用；如果区段还特别短，或者填表时多除了个 100，抽样出的 TTF 会大到超过一万年，一万年的仿真里一次都抽不中。另一种可能是np.log(0)产生的inf进了min()，被当成“永远不会到期”。

解决：在模型初始化时加一行断言，if max(c.lam_per_yr for c in net.components.values()) < 1e-4: raise ValueError("检查 lambda 单位")；同时在抽样函数里对u做下限截断（max(u, 1e-12)）。排查时最快的办法是打印首年的事件表，看min(remaining)的量级——正常应该在 0.01 到 1 年之间，如果是inf或者几百，问题就定位到了。

5.4 跑一次要几十分钟，加到两万年后内存还在涨

现象：任务管理器里内存曲线一路往上，跑一万年时占用好几个 G，速度也肉眼可见地变慢。

原因：每年的结果都append到 list 里，两万年乘以 10 个负荷点就是 20 万个浮点数，本身不大，但如果同时保存了每年的完整状态快照（某些实现为了“可回溯”会这么干），内存就是元件数乘以年数。另一种情况是内层循环里反复构造字典和 DataFrame，Python 的对象开销远大于数值计算本身。

解决：指标累加换成 4.3 的 Welford 在线更新，只保留标量；每年的中间结果只在--debug模式下落盘；内层循环里用数组而不是字典做索引。改完之后，同样的算例耗时能从分钟级降到秒级，这两万年的样本本身不是瓶颈，瓶颈在解释器开销。

5.5 换个随机种子，SAIDI 差 20%

现象：同一个输入文件，--seed换个值，SAIDI 从 3.12 变成 3.74，两次都“收敛”了。

原因：两次运行都在 β 刚好压线时停止，样本量本身处在临界状态，均值的置信区间还没收窄；如果程序里还用time.time()或者不设种子，问题会被掩盖成“蒙特卡洛本来就有随机性”。

解决：固定种子并把种子值写进system_index.csv；把beta-target从 0.05 收紧到 0.02 重跑一次，观察两次结果的差异是否落到 5% 以内。工程上判断结果稳不稳，最快的办法是两个不同种子各跑一遍，看两次的差值是否小于 β 对应的置信区间半宽——如果大于，说明 β 判据设得太松，或者收敛判据用错了指标。

6. 方差缩减与自检：先估样本量，再动手跑

跑蒙特卡洛最容易犯的错是“先跑再说”。正确的顺序是用一小批样本估出 σ 和 μ，反推需要的样本量，再决定max-years设多少。正态近似下，所需样本年数 n ≈ (σ / (β·μ))²，用两百年的预跑就能把 σ 估出来。

import numpy as np def required_years(sigma, mu, beta_target): """由目标方差系数反推所需样本年数（正态近似）""" return int(np.ceil((sigma / (beta_target * mu)) ** 2)) # 预跑 200 年，收集每年的 SAIDI pre = np.array(pre_saidi_list) # 长度 200 n_need = required_years(pre.std(ddof=1), pre.mean(), 0.02) print(f"目标 beta=0.02，预计需要 {n_need} 个样本年")

这段代码的价值在于把“跑多久”从拍脑袋变成可计算。比如预跑估出 σ=2.1、μ=3.1，目标 β=0.02，算出来需要约 4600 个样本年；而如果目标放宽到 0.05，只需要 730 年，差了六倍多。先算这笔账，再决定要不要上分层抽样——很多时候不是算法不够快，是精度目标定得不合理。

另一个自检手段是把模型退化到解析法能算的极端情形。把switch_hr全部设为 0、所有负荷点has_backup=0、只保留串联关系，此时系统的 SAIDI 应该等于各负荷点最小路上元件 λ×r 之和的加权平均。程序结果和手算结果对不上，问题一定在拓扑解析或判定逻辑里，而不是在抽样。这个自检我一般放在改完判定函数的第一次运行，比逐行 debug 快得多。

进阶一点的方差缩减是分层抽样：把主干区段按故障率高低分成两层，低

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