简介：面向电力系统与配电网研究者的MATLAB实现资源，聚焦最小路法和非序贯蒙特卡洛算法在配电网可靠性评估中的应用。资源包内含IEEE RBTS可靠性测试系统的原始参数文档（PDF与Excel格式）、IEEE33节点系统参数，以及RBTS BUS6的MATLAB数据文件；同时提供两套完整程序：基于最小路法的可靠性评估程序，以及采用节点影响分析法判断受影响负荷的蒙特卡洛主程序。压缩包体积仅1.34MB，主要为PDF、Excel和MATLAB脚本，便于离线查阅与直接运行。已有3816人学习下载，适合电气工程专业学生、电网科研人员及从事可靠性分析的工程师。文件组织清晰，PDF参数文档方便核对原始数据，Excel表格可直接修改导入，MATLAB脚本注释完整，既可用于计算可靠性指标，也能作为理解两类算法原理、改造适应自有网络拓扑的参考模板。整体精简而完整，适合作为课程设计或科研起步的实用参考。
1. 配电网可靠性评估的两种算法：最小路法和蒙特卡洛法是不是只能二选一
我去年调一个乡镇10kV配网的可靠性数据，SAIFI怎么算都跟真实统计对不上。后来发现问题不在数据脏，而在方法选错了：一条带联络开关的网格状线路，我硬用最小路法去解析，结果永远解释不了负荷转移那几小时。配电网可靠性评估的核心问题其实很朴素——一年里一个用户平均停几次电、每次停多久。最小路法把拓扑拆成每条负荷到电源的必经之路，用枚举方式算出期望值；蒙特卡洛法则模拟几十年甚至上百年的随机故障，再统计停电事件。两者不是二选一的竞争关系，而是不同网架、不同数据条件下的互补工具。这篇文章写给要做配网规划、供电可靠性管理或新能源接入方案比选的工程师，读完你能判断该用哪个，并且能直接跑通一个最小可用的评估脚本。
2. 最小路法到底怎么算：从找最小路到输出SAIFI的全过程
2.1 最小路法的建模思路：先画出每个负荷点的“命脉”
最小路并不是图论里的最短路径，而是电力实际流通的那条必经之路：从变电站母线出发，经过各段线路、开关、配变，直到某个负荷点。对辐射状配电网来说，这条路径通常唯一，所以最小路法才特别适合。方法的核心是给每个负荷点分类元件：在最小路上的元件故障，负荷点只能等它修好，停电时长等于修复时间；不在最小路上的元件故障，负荷点不一定完全停电，可能只经历一次开关操作期间的短时停电，甚至因为熔断器或分支保护而完全不受影响。
这个分法实际就是故障模式影响分析（FMEA）的变体。它的好处是模型透明，每个元件故障影响谁、影响多久都能摊开来查。坏处也很明显：遇到多电源、有联络开关或者能转供的闭环网架，最小路的“唯一性”被打破，必须额外定义转供时间，否则算出来的指标会偏大。所以在跑最小路法之前，我会先花半小时把网架图和保护配置画清楚，尤其是隔离开关、熔断器、联络开关的位置，这些直接决定非最小路元件的故障影响范围。
下面的表格是影响时间分类，我每次建模都会对照它过一遍：
| 元件位置 | 对负荷点的影响 | 影响时间取值 |
|---|---|---|
| 负荷点最小路上的元件 | 停电至修复完成 | 修复时间 r |
| 分支上的元件，且分支有熔断器/保护 | 不影响其他分支负荷 | 0 |
| 分支上的元件，没有可靠保护 | 上游开关动作后恢复 | 开关操作时间 s |
| 馈线入口与上游主网元件 | 整条馈线失压 | 随场景取 r 或 s |
这个表不是死的，具体取决于保护配置。写代码前把每个元件的影响负荷点和影响时间列出来，是最小路法最重要的准备工作，比写代码本身更容易出错。
2.2 用Python实现FMEA式的最小路法：一个可以直接抄的脚本
最小路法的常规落地做法是枚举所有元件故障，逐个判断它对每个负荷点的影响。下面这段代码我用一个带两个负荷点的简化辐射网演示，数据写在字典里，便于你改成自己的网架。
import networkx as nx # 拓扑：S为电源，LP1、LP2为负荷点 G = nx.Graph() G.add_nodes_from(["S", "A", "B", "C", "LP1", "LP2"]) edges = [ ("S", "A", 0.20, 5.0), # 进线：故障率0.2次/年，修复5小时 ("A", "B", 0.10, 4.0), # 主干线段 ("B", "LP1", 0.10, 3.0), # 分支到LP1 ("A", "C", 0.30, 5.0), # 分支到C ("C", "LP2", 0.15, 3.0), # 分支到LP2 ] for u, v, lam, r in edges: G.add_edge(u, v, lam=lam, r=r) # 每条边故障时，对每个负荷点的影响时间（小时） # 影响时间为0表示不影响该负荷点 impact_time = { ("S", "A"): {"LP1": 5.0, "LP2": 5.0}, ("A", "B"): {"LP1": 4.0, "LP2": 4.0}, ("B", "LP1"): {"LP1": 3.0, "LP2": 0.0}, ("A", "C"): {"LP1": 0.0, "LP2": 5.0}, ("C", "LP2"): {"LP1": 0.0, "LP2": 3.0}, } loads = {"LP1": 120, "LP2": 80} # 每个负荷点的用户数 bd = {lp: 0.0 for lp in loads} # 年停电小时数期望 U ld = {lp: 0.0 for lp in loads} # 年停电次数期望 lambda for (u, v) in impact_time: lam = dict(G[u][v])["lam"] for lp, hours in impact_time[(u, v)].items(): if hours > 0: ld[lp] += lam bd[lp] += lam * hours # 系统指标：SAIFI / SAIDI / CAIDI total_customers = sum(loads.values()) SAIFI = sum(loads[lp] * ld[lp] for lp in loads) / total_customers SAIDI = sum(loads[lp] * bd[lp] for lp in loads) / total_customers CAIDI = SAIDI / SAIFI if SAIFI > 0 else 0 for lp in loads: print(f"{lp}: lambda={ld[lp]:.3f} 次/年, U={bd[lp]:.3f} 小时/年") print(f"SAIFI={SAIFI:.3f} 次/(用户·年), SAIDI={SAIDI:.3f} 小时/(用户·年), CAIDI={CAIDI:.3f} 小时/次")这段代码的关键不是networkx画图，而是impact_time这个字典。它就是你把2.1节的保护配置表翻译成程序的过程。lam是元件故障率，单位是“次/年”；hours是影响时间，单位是“小时”。把每条边的lam * hours累加到负荷点上，就得到该负荷点的年停电小时期望值 U；把所有lam累加，得到年停电次数期望 lambda。最后用负荷点的用户数做加权，就是标准的SAIFI和SAIDI。
跑这组参数，LP1的停电频率约0.40次/年、年停电时长约1.7小时，LP2约0.75次/年、3.35小时/年，系统SAIFI约0.54次/(用户·年)，SAIDI约2.36小时/(用户·年)。读者可以拿这个结果做基准，验证自己写的版本有没有对齐。
实际工程里网架会比这个demo大得多。我一般不会手写impact_time，而是从电网GIS导出馈线拓扑，再用一段脚本自动生成“元件—影响负荷点”矩阵。但自动化的前提还是要有一套明确的开关逻辑，否则矩阵生成出来也是错的。
2.3 最小路法的适用边界：辐射网好用，遇到转供要小心
最小路法最大的优势是快、透明、结果确定。同样的网络跑一百次，结果都一样，不需要担心随机误差。所以它特别适合规划阶段的方案比选：十个接线方案，每个都算一遍SAIFI/SAIDI，几分钟就能排个高低。
但边界条件也很硬。第一，多电源或环网运行时，负荷点有两条及以上供电路径，最小路不再唯一，传统最小路法会低估可靠性。第二，网架中存在联络开关时，非故障段可能通过转供恢复供电，这段恢复时间既不是修复时间，也不是0，必须在故障影响分析里显式建模。第三，分布式电源接入后可能形成孤岛运行，这种状态用最小路法表达非常别扭。
所以我在实际项目里的判断是：纯辐射状、没有分布式电源、转供关系简单的馈线，直接上最小路法；出现环网转供、微网孤岛、电气距离复杂的情况，就换蒙特卡洛法，或者在最小路法里手动加入“转供恢复时间”这个修正项。别把最小路法当作万能工具，它在简单场景下是利器，在复杂场景下就是坑。
3. 蒙特卡洛法模拟配电网：序贯抽样的Python实现与收敛判据
3.1 序贯与非序贯蒙特卡洛：怎么选
蒙特卡洛法的本质是用大量随机样本来逼近真实的故障-修复过程。配电网里最常用的两种是序贯和非序贯。非序贯方法对每个元件按“当前时刻处于故障状态的概率”独立抽样，生成一个系统状态，再统计该状态下的停电情况。它不考虑这些事件在时间上的先后，计算快，但很难处理时变负荷、天气变化和修复顺序等带时间约束的因素。
序贯方法则沿着时间轴推进，先抽样每个元件的正常运行时间，再抽样故障持续时间，形成一条连续的状态演化序列。它保留了故障发生的先后顺序，也天然兼容“多个元件同时故障”“故障后转供”“负荷随时间波动”这些场景。代价是计算量大得多，尤其当元件数量多、故障率低的时候，可能跑几千“系统年”才能让指标稳定下来。
我的选择原则是：只算年期望值、网架规模大、需要快速筛选场景，用非序贯；要做精细化评估、研究某个具体故障序列或验证转供策略，用序贯。实际工作中可靠性数据往往带季节性，故障率和修复时间在不同天气下差别很大，所以我个人更常写序贯版本，宁可多等几分钟。
下面用一个带两负荷点的辐射网演示序贯蒙特卡洛的完整流程，这段代码可以直接改成你自己的网络。
3.2 一个可运行的序贯蒙特卡洛脚本：从元件状态序列到系统指标
这个脚本的思路是：对每条边，用指数分布抽样正常运行持续时间和故障修复时间，生成若干年的故障事件列表；把所有边的事件按时间排序；然后逐个事件把故障边从网络里移除，用图的连通性判断每个负荷点是否失电。为了方便展示，我假设故障事件在时间上不重叠，这对故障率低的小型网络是合理的简化。
import numpy as np import networkx as nx # 与2.2节相同的网络和参数 edges = [ (("S", "A"), 0.20, 5.0), (("A", "B"), 0.10, 4.0), (("B", "LP1"), 0.10, 3.0), (("A", "C"), 0.30, 5.0), (("C", "LP2"), 0.15, 3.0), ] G = nx.Graph() G.add_nodes_from(["S", "A", "B", "C", "LP1", "LP2"]) for (u, v), lam, r in edges: G.add_edge(u, v, lam=lam, r=r) loads = {"LP1": 120, "LP2": 80} source = "S" years = 20 rng = np.random.default_rng(2024) # 固定种子方便复现，正式评估不要固定 def simulate_component(lam, r, horizon_hours, rng): """生成单个元件的故障时段，单位统一为小时""" t = 0.0 events = [] while t < horizon_hours: t += rng.exponential(8760 / lam) # 正常工作时间 if t >= horizon_hours: break down = rng.exponential(r) # 故障修复时间 events.append((t, min(down, horizon_hours - t))) t += down return events horizon = years * 8760 all_events = [] for (u, v), lam, r in edges: for start, down in simulate_component(lam, r, horizon, rng): all_events.append((start, start + down, (u, v), down)) all_events.sort(key=lambda x: x[0]) # 统计负荷点停电次数和年停电小时 outage_count = {lp: 0 for lp in loads} outage_hours = {lp: 0.0 for lp in loads} for start, end, edge, down in all_events: u, v = edge G.remove_edge(u, v) for lp in loads: # 故障期间负荷点是否与电源不再连通 if not nx.has_path(G, source, lp): outage_count[lp] += 1 outage_hours[lp] += down G.add_edge(u, v) # 输出结果，与最小路法结果对拍 for lp in loads: avg_freq = outage_count[lp] / years avg_hours = outage_hours[lp] / years print(f"{lp}: 年均停电 {avg_freq:.3f} 次, 年停电时长 {avg_hours:.3f} 小时") safi = sum(loads[lp] * outage_count[lp] for lp in loads) / (sum(loads.values()) * years) saidi = sum(loads[lp] * outage_hours[lp] for lp in loads) / (sum(loads.values()) * years) print(f"模拟SAIFI={safi:.3f}, 模拟SAIDI={saidi:.3f}")这段代码里最需要注意的单位统一：lam是次/年，r是小时，而抽样时需要把故障率换算成小时尺度，所以正常工作时间用8760 / lam。rng.exponential抽的都是指数分布变量，这是标准做法，默认元件寿命和修复时间服从指数分布。
运行后你会看到模拟值接近最小路法的结果，但不会完全一致。比如我曾用同一组参数跑到SAIFI=0.55左右、SAIDI=2.2~2.5小时，这个波动就是抽样误差。如果你把years改成200，误差会缩小，但计算时间也会线性增加。对辐散状的网络，序贯蒙特卡洛真正的费用在多次连通性判断上，节点多了以后不建议每次都全图搜索。
3.3 收敛判据与仿真年限：抽多少次才能信
很多新手拿到蒙特卡洛结果就问“这个数能不能用”。我一般不会只看一次运行结果，而是看两个东西：年度指标的变异系数，以及置信区间。变异系数 CV 是年指标标准差除以均值，CV 越小说明结果越稳定。常见做法是把 CV 控制在0.05以内，否则就继续增加仿真年限。
一段简单的判据代码可以这样写：
import numpy as np # yearly_values 是一个数组，每个元素是某一年的SAIDI或SAIFI def check_convergence(yearly_values, target_cv=0.05): mean = np.mean(yearly_values) std = np.std(yearly_values, ddof=1) cv = std / mean if mean > 0 else float('inf') return cv < target_cv, cv这里的yearly_values怎么来？在序贯模拟中，可以把整个长模拟切分成一年一段，每一段单独统计一次停电次数和停电时长；在非序贯模拟中，每次抽样生成一个系统状态，相当于一年。把这些年值积累下来，CV就自然算出来了。
要特别提醒：别为了CV好看而用固定随机种子把结果卡死。固定种子只能用来复现脚本，不能证明模型正确。正确的做法是跑不同年数、不同随机种子，看结果是否在合理区间内波动。如果SAIDI每次变化超过20%，说明仿真年限不够，而不是方法有问题。对故障率偏低的配电网，模拟几千年系统年仍然不收敛的情况我都遇到过，这种时候不要盲目加年限，先检查是不是有某个小概率高影响事件没建模，再考虑用重要抽样这类方差缩减技术。
4. 两法对拍：用最小路法校准蒙特卡洛的五个关键参数
4.1 用最小路法做基准：给蒙特卡洛结果套置信区间
最小路法算出的是一个确定性期望值，蒙特卡洛算的是一堆随机样本的均值。两者既然面对同一张网络和同一组参数，结果就不该差太远。所以我的习惯是：每次跑完蒙特卡洛，都拿最小路法结果当靶子，看模拟结果是否落在置信区间内。置信区间可以简单用经验公式算：
import numpy as np def conf_interval(values, alpha=0.05): mean = np.mean(values) std = np.std(values, ddof=1) n = len(values) z = 1.96 if alpha == 0.05 else 2.58 return mean, mean - z*std/np.sqrt(n), mean + z*std/np.sqrt(n)如果你把模拟年份切分成多个观测样本，这里的values就是每年的SAIFI或SAIDI。以2.2节最小路法得到的SAIFI=0.54、SAIDI=2.36小时为基准，如果蒙特卡洛的95%置信区间是0.48~0.60和2.1~2.6小时，那这组模拟是可信的；如果区间完全偏离靶心，就要回头查代码或参数。
对拍还有一个用处：当最小路法因为网架复杂而无解时，它可以作为蒙特卡洛结果的一个近似参考点。反过来，蒙特卡洛也能暴露最小路法没有建模的那些随机波动。所以两种方法在我的项目里总是成对出现，不是互斥。
4.2 五个关键参数：取值量级与对指标的影响
不管用哪种方法，下面这五个参数决定成败。我把常用量级和影响方向整理成表，供参考：
| 参数 | 常见取值量级 | 主要影响 |
|---|---|---|
| 线路故障率 lambda | 架空线0.05~0.5次/km·年，电缆0.01~0.1次/km·年 | 近似正比影响SAIFI和SAIDI |
| 故障修复时间 r | 架空线3~6h，电缆4~8h | 正比影响SAIDI和CAIDI |
| 隔离开关/联络开关操作时间 s | 人工0.5~2h，自动开关2~10min | 影响转供场景的SAIDI |
| 负荷点用户数与负荷分布 | 按配变容量和低压台区折算 | 改变系统指标的加权权重 |
| 仿真年限/抽样数 | 序贯至少1000~2000系统年，非序贯更高 | 决定蒙特卡洛抽样误差 |
故障率数据我一般优先用当地供电公司发布的可靠性统计，没有就按设备型号和运行年限经验估算。修复时间也要看抢修资源，偏远地区普遍偏高。开关操作时间则要区分“人工到现场操作”和“遥控自动操作”，两者能差一个数量级，最容易被拍脑袋带偏。
4.3 故障率和转供时间的影响差异：为什么不能只看平均值
先说故障率敏感性：如果某条馈线的故障率从0.1次/年变成0.2次/年，而其他参数不变，SAIFI近似翻倍，SAIDI也接近翻倍。这是因为两个指标都跟故障次数线性相关。但如果你在北方的冬天把架空线路故障率调高到0.4，修复时间从4h变成8h，SAIDI的变化会比SAIFI更夸张，CAIDI也会明显升高。这个差异说明，SAIFI主要反映设备故障频率，SAIDI还受抢修资源和开关策略影响，做敏感性分析时不能只看一个数。
转供参数的影响方向则相反：它几乎不影响SAIFI，因为停电次数已经发生，但能显著减少SAIDI。一个原本要停电6小时的非故障段，如果通过联络开关在1小时内转供，SAIDI就降了5小时。所以在有联络开关的网架里，对拍时如果发现蒙特卡洛SAIDI比最小路法小很多，先别急着说代码错了，很可能是因为最小路法把转供时间设成了0，而蒙特卡洛的开关模型把它放进了流程里。
我建议每个项目至少做三组敏感性：一组把故障率整体上下调20%，一组把修复时间上下调50%，一组把转供时间在“人工2h/自动10min/不转供”三种模式中切换。这样能快速找到指标波动的最大来源，也能判断手上数据够不够支撑评估精度。
5. 避坑指南：可靠性评估中常见的6个翻车现场
5.1 现象：最小路法算出的SAIFI比统计值小一截
常见做法是只把负荷点最小路上的元件，或者只把主干道的故障率累加，忽略了分支线故障和馈线出口断路器对短时停电的贡献。我接过一个项目，进线故障率明明有0.2次/年，但团队只统计了馈线主干段，结果SAIFI少了将近两成。原因不是公式错，而是故障元件的影响范围没列全。解决办法是回到FMEA表，把每个元件影响哪些负荷点、影响多久重新过一遍，尤其是变电站出线以前的元件和分支线路，必须显式写进impact_time。
5.2 现象：蒙特卡洛每次跑出来的SAIFI差很多
第一次跑0.50，第二次跑0.62，第三次又回到0.53。如果你只跑了500个系统年，这种现象在故障率低的配电网里非常正常，因为年故障事件本来就稀疏，偶发一次大故障就能把均值拉高。解决思路是先算变异系数，低于0.05再定数据；再就是增加仿真年限，或者改用分层抽样让偶发大故障更平均地出现在样本里。不要试图用固定随机种子消除波动，那是自欺欺人。
5.3 现象：把年故障率当成单次概率，结果直接指数爆炸
在非序贯蒙特卡洛里，有人会写if np.random.uniform() < lam: 故障。这行代码看似合理，实际上完全错误，因为lam是“每件每年故障多少回”，不是“某次抽样的概率”。用0.1次/年去做单次概率，相当于默认每年发生十几次故障，得到的SAIFI会比真实值大一个数量级。正确做法是用泊松分布抽样一年内的故障次数np.random.poisson(lam)，或者干脆按时间轴做指数分布抽样。我见过太多翻车都出在这一行，值得反复提醒。
5.4 现象：变压器低压侧故障被计入系统指标
配电网可靠性评估的统计边界一般是变电站10kV出线开关以下到配变低压侧出口。如果把低压线路、用户内部故障也算进来，SAIDI会高得离谱。原因是这些部分不属于中压配网运维单位的管理范围，数据口径不一致。解决办法是在建模一开始就建立边界规则：只统计10kV馈线、分段开关、联络开关、配变高压侧，低压侧一律不参与系统指标汇总。如果用户要求评估到户，就单独做低压可靠性模块，不要混在中压模型里。
5.5 现象：模拟了几百年，收敛判据还是不达标
这种情况往往不是样本量不够，而是模型里有小概率高影响事件，比如一条主干线一年发生一次台风倒杆，一停就是十几个小时。这种事件对SAIDI方差贡献极大，普通随机抽样需要极多年份才能稳定它。解决方法是把天气状态单拎出来：正常天气用一个故障率，恶劣天气用更高的故障率和更长的修复时间，再做状态抽样。也可以引入重要抽样，但工程上先把天气状态建好基本就能解决问题。
5.6 现象：联络开关的转供时间没进模型
有联络开关的网架，非故障段往往可以在隔离后转供复电，停电时间不是修复时间，而是“隔离+合联络开关”的操作时间。我在一个新区网架上吃过亏：最小路法模型里没写转供，结果SAIDI算出来2.8小时，实际运行只有1.9小时。反过来，如果把自动化开关的转供时间设成0，又会把SAIDI算得过于乐观。正确做法是给每条联络开关设置一个转供时间，在故障影响分析里对每个负荷点取min(修复时间, 隔离时间+转供时间)。这个逻辑要写清楚，别一个开关时间套全部线路。
6. 把评估做成参数驱动的工具：批量场景跑法与置信区间报告
前面几张图的代码都是一次性脚本，但实际做规划时往往要比较几十组方案。我现在的做法是把所有参数放到一张CSV里，一行动一个场景，脚本循环读取后输出指标表。核心思路是：拓扑结构写在图构建函数里，所有故障率、修复时间、操作时间、负荷用户数都从参数表读入，这样换网架只需要改一行拓扑，换数据只需要改CSV。
给一个最小的参数驱动骨架：
import pandas as pd import numpy as np # params.csv包含列: scene_name, lam, r, switch_time, customers df = pd.read_csv("params.csv") results = [] for _, row in df.iterrows(): lam = row["lam"] r = row["r"] # 调用你自己的最小路法或蒙特卡洛函数 saifi = run_min_path(...) # 具体函数见第2、3章 saidi = run_monte_carlo(...) results.append({"scene": row["scene_name"], "SAIFI": saifi, "SAIDI": saidi}) result_df = pd.DataFrame(results) result_df.to_csv("reliability_report.csv", index=False)批量场景跑完，还要给每个SAIDI补一个置信区间，用4.1节的conf_interval就行。报告里我会把每次蒙特卡洛的均值、置信区间、样本年限写在一起，不让别人只看到一个干巴巴的点估计。那种“最小路法算出来0.54，蒙特卡洛算出来0.55”的结论，只有带上置信区间才有意义。
我现在的习惯是：拿到任务先不写代码，先画包含联络开关、分段开关和保护配置的网架图。画完再决定用最小路法还是蒙特卡洛法，然后才动手写脚本。这个顺序让我少翻了好几次车。做配电网可靠性评估，最重要的不是算法多高深，而是边界、参数和转供逻辑都清清楚楚。希望这个习惯也能帮到你。
本文还有配套的精品资源，点击获取
