去年年底整理实验室研究材料的时候，导师在小组会上问了我一个问题：咱们这个示范区的配电网，可靠性到底能打几分？我当时的反应是愣了几秒——手头有拓扑图、有设备台账、有负荷数据，但真要我说清楚"一年下来平均每户会停几次电、每次停多久、供电可用率是多少"，我居然给不出一个数。也就是从那时候起，我开始认认真真做了一个配电网可靠性评估的Matlab程序，一路踩坑，也一路把以前课本上模糊的概念全部落地成了能算数的代码。这篇文章就把我这段实现之路完整记录下来，包括指标怎么定义、方法怎么选、代码结构怎么搭、算例结果对不对，以及调试过程中踩过的那些坑。如果你是电力系统方向的研究生、配网规划工程师，或者正在被各种可靠性课程设计折磨的本科生，这篇应该能直接当成一份可参考的路线图。
1. 动手前先把指标吃透：可靠性评估到底在算什么
1.1 五个绕不开的指标：从定义到物理意义
很多人一上来就急着写代码，结果算完之后根本不知道结果合不合理。我建议先把指标这关过了，因为后面程序里每一个变量、每一次累加，本质上都是在为这几个指标服务。
配电网可靠性评估用得最多的是这么几个指标：SAIFI（系统平均停电频率指标）、SAIDI（系统平均停电持续时间指标）、CAIDI（用户平均停电持续时间指标）、ASAI（平均供电可用率指标），以及EENS（期望缺供电量）。它们的定义和公式我整理成了表格，方便后面写代码时对照。
| 指标 | 物理含义 | 计算公式 | 典型量级 |
|---|---|---|---|
| SAIFI | 每个用户在一年内平均停电次数 | 用户停电次数总和 / 总用户数 | 0.5~2次/户·年 |
| SAIDI | 每个用户在一年内平均停电小时数 | 用户停电持续时间总和 / 总用户数 | 1~5小时/户·年 |
| CAIDI | 每次停电平均持续多久 | SAIDI / SAIFI | 1~3小时/次 |
| ASAI | 一年内供电可靠时间占比 | 1 - SAIDI / 8760 | 99.9%以上 |
| EENS | 因故障而少送的电量 | 各场景停电功率×停电时间之和 | 视系统规模而定 |
写代码时最容易混淆的是SAIFI和SAIDI。SAIFI只看"次数"，所以某个负荷点只要经历了一次停电，哪怕只停1分钟，也要计入频率；而SAIDI看的是"时长"，停1分钟和停3小时，在SAIDI里的权重完全不同。这两个指标一个管频次、一个管时长，对应到配网运维上，就是两条不同的改进思路：一个靠减少故障次数（比如线路绝缘化改造、防雷改造），另一个靠缩短修复时间（比如加装故障指示器、完善转供能力）。
EENS比较特殊，它需要知道负荷点的平均功率。比如一个节点有200户居民，户均峰值功率按2kW估算，那么该节点负荷就是400kW。如果这个节点在一年里累积停电3小时，那EENS就是1.2MWh。这个指标在做经济性评估时特别有用，因为可以给它乘上电价，算出"少送了多少电费"。
1.2 三种时间参数背后的运行逻辑
程序里反复用到的时间参数有三个：故障修复时间、故障隔离时间、负荷转供时间。很多初学的人把它们混为一谈，这是导致计算结果错误的最常见原因。
我先讲一下配电网故障处置的真实流程。假设一条10kV馈线的某一段线路被外破打断了，变电站出线断路器会先跳闸，整条馈线瞬间失电。这时候调度员并不知道具体故障点在哪，需要根据保护动作信息、故障指示器信号去定位故障区段。找到故障点之后，如果能遥控操作，就把故障点两侧的隔离开关拉开，这叫"故障隔离"。隔离完成之后，如果故障点靠电源侧的区域没有损坏，就可以合上出口断路器先恢复送电；如果故障点后端还有网架结构、存在联络开关，可以通过合联络开关把下游负荷转到另一条馈线上供电，这叫"转供"。最后维修人员到现场把故障元件修好或换掉，整个系统恢复到初始运行方式。
对应到程序里，这几种时间要分开处理：
- 故障修复时间：指故障元件从损坏到修好重新投运的时间，通常是几小时到十几小时，线段越长、交通越不方便，修复时间越长。
- 故障隔离时间：指从故障发生到隔离操作完成的时间，如果是自动化开关，可能是秒级到分钟级；如果是人工操作，可能是0.5~2小时。
- 负荷转供时间：指合上联络开关、把负荷倒到另一条馈线的时间，自动化的按分钟记，手动的按0.5~1小时记。
我在程序里设置了一个底层逻辑：所有非故障区域的停电时间，都不应该超过隔离时间和转供时间的组合；而故障元件所在区域下游、且无法转供的那部分负荷，停电时间才是修复时间。这个逻辑直接决定了SAIDI算得对不对。后面在第3节我会详细展示这段代码怎么实现。
2. 方案选型：为什么我用故障枚举而不是蒙特卡洛
2.1 解析法与模拟法的取舍
配电网可靠性计算方法大体分两类：解析法和蒙特卡洛模拟法。解析法的典型代表是故障模式后果分析法，英文缩写FMEA，基本思路就是"枚举每一个元件故障，分析后果，然后汇总指标"。蒙特卡洛法则完全不同，它是通过生成随机数来抽样模拟元件故障和修复过程，跑成千上万次仿真，统计出可靠性指标的期望值。
这两种方法怎么选？我的经验是，对大多数配电网规划和研究场景，FMEA足够了，而且优势很明显。第一，FMEA结果稳定，同一套参数跑多少遍结果都一样，不会出现蒙特卡洛那种"这次跑和上次跑差几个百分点"的情况；第二，FMEA计算速度快，一个几十条支路的配电网，枚举一遍所有元件故障，在Matlab里几秒钟就能出结果；第三，FMEA天然包含灵敏度信息——你能明确看到"哪条线路故障对SAIFI贡献最大"，这对改造决策非常有价值。
蒙特卡洛的优势在于处理复杂系统，比如考虑时序变化、分布式电源出力波动、元件多重故障等场景。但它的收敛速度比较慢，误差按抽样次数的平方根下降，想提高一位精度就要多100倍计算量。对于一篇课程设计或者一个中等规模的配电网规划项目，我建议优先上FMEA，把模型本身吃透之后再考虑蒙特卡洛扩展。
2.2 网络拓扑怎么存：邻接矩阵、支路表和图对象
程序的数据结构设计，决定了后面代码的复杂度和可维护性。我第一次写的时候用邻接矩阵，后来发现很痛苦——邻接矩阵适合描述两节点之间直接相连的关系，但要表达"支路上串联了开关""这条支路的类型是架空线还是电缆"这种工程属性，它就力不从心了。
我最终采用的是"支路表"结构，这是工程上最常用、也最直观的存储方式。每一行代表一条支路，包含支路编号、首端节点、末端节点、元件类型、故障率、修复时间、是否含开关等字段。这种结构的好处是，你可以把任意复杂的工程属性塞进对应的列里，而且排查数据错误时非常直观。
节点之间怎么连通的？我建议同时用一个Matlab的graph对象来维护纯拓扑关系。它的好处是调用shortestpath、conncomp这些图论函数特别方便，判断"某个节点是否和电源节点连通"只需要一行代码。我在程序里把graph对象当成快速查询工具，把支路表当成权威数据源，两者配合，兼顾了灵活性和准确性。
2.3 分类讨论故障影响范围：FMEA的核心
FMEA的核心是"故障后果分析"。对于一个辐射型配电网，某个元件发生故障后，它的影响范围可以分为三个层次：
第一类，故障元件所在馈线的变电站出口断路器跳闸，导致该馈线所有负荷瞬间失电。这不是最终结果，后续还要看隔离和转供的执行情况。第二类，通过拉开故障段两端的隔离开关，故障段靠电源侧的区域可以快速恢复，这个区域的停电时间只有隔离开关操作时间。第三类，故障段下游的区域，如果存在联络开关且能闭合，则停电时间是转供操作时间；如果不存在联络开关，则只能等故障修复，停电时间就是修复时间。
我在程序里把每个负荷点的停电时间处理成了一个三元组：停电频率+停电时长+缺供电量。例如某负荷点在故障发生后被隔离操作恢复了，那么频率加1，时长加隔离时间；如果该负荷点位于故障下游且无转供路径，频率加1，时长加修复时间。这里有一个细节容易错：即使某个负荷点在故障后很快恢复了供电，只要它经历了停电，SAIFI里就要计入一次。很多新手在写代码时只统计"最终失电时间大于0"的负荷点，导致SAIFI被算小，这是要特别警惕的。
3. Matlab代码怎么写：从数据结构到核心函数
3.1 数据结构设计：元件表、负荷表、系统参数
我习惯先把所有数据定义集中在一个脚本或函数里，方便修改和调试。元件表用Matlab矩阵来表示，每一行是一个元件，列的含义通过注释写清楚。下面是我在测试算例里用的一组结构示例：
% 元件表列定义： % 编号 | 类型(1线路/2变压器) | 首端节点 | 末端节点 | 故障率(次/年) | 修复时间(h) | 是否含开关(0/1) elements = [ 1, 1, 1, 2, 0.100, 5.0, 0; 2, 1, 2, 3, 0.085, 5.0, 0; 3, 1, 3, 4, 0.120, 5.0, 1; 4, 1, 2, 5, 0.070, 5.0, 0; 5, 1, 5, 6, 0.095, 5.0, 1; 6, 2, 6, 7, 0.015, 8.0, 0; ];负荷表单独存储，每一行是一个负荷节点，包含节点编号、用户数量、平均负荷功率，比如loads = [4, 150, 300; 7, 220, 440]表示节点4有150户用户、平均负荷300kW。
我还单独定义了几个全局时间参数：isoTime（故障隔离时间，小时）、tieTime（联络开关转供时间，小时）、annualHours（8760小时）。这些参数独立于元件表，是为了方便做敏感性分析——比如我想看看"如果自动化改造把隔离时间从1小时压缩到0.2小时，SAIDI能改善多少"，只需要改一个变量，程序重跑一遍就行。
3.2 写一个能跑通的故障枚举主循环
主循环的逻辑其实不复杂，就是"遍历元件表，逐个故障，逐个分析"。核心流程是：对每一个元件，把它从拓扑图中去掉，模拟故障状态；然后判断每个负荷点是否失电；再根据隔离、转供逻辑计算失电时长；最后累加到总指标里。下面是主循环的框架代码：
% 初始化累计变量 totalInterruptions = 0; % 用户停电次数总和 totalOutageHours = 0; % 用户停电小时总和 totalMissingEnergy = 0; % 缺供电量总和 totalUsers = sum(loads(:,2)); for k = 1:size(elements, 1) % 获取当前故障元件的两端节点 n1 = elements(k, 3); n2 = elements(k, 4); % 从拓扑图中临时删除故障支路 G_temp = rmedge(G, n1, n2); % 对每个负荷点判断失电情况 for j = 1:size(loads, 1) loadNode = loads(j, 1); users = loads(j, 2); power = loads(j, 3); % 判断电源节点与负荷节点是否连通 if ~isempty(shortestpath(G_temp, sourceNode, loadNode)) continue; % 仍连通，不受影响 end % 到这里说明负荷点失电了，先累计频率 totalInterruptions = totalInterruptions + users; % 计算失电时长（需要调用自定义函数，见3.3节） outageHours = calcOutageDuration(elements, G_temp, k, loadNode); totalOutageHours = totalOutageHours + users * outageHours; totalMissingEnergy = totalMissingEnergy + power * outageHours; end end % 最后计算系统级指标 SAIFI = totalInterruptions / totalUsers; SAIDI = totalOutageHours / totalUsers; CAIDI = SAIDI / max(SAIFI, 1e-9); ASAI = 1 - SAIDI / 8760; EENS = totalMissingEnergy;这里有一个重要的编程坑：rmedge之前必须先把原图备份，因为每个故障场景都要基于原始完整拓扑去删边。我第一版代码就是在循环里直接改了G本身，第二圈循环就出错了。建议在主循环开始前把原始图存一份，比如G_base = G;，每次循环用G_temp = rmedge(G_base, n1, n2)。
3.3 故障影响范围判断的逻辑实现
这一节是最核心、也最容易出错的。上面的主循环里，我只是判断了"负荷点是否失电"，但实际的停电路径还得分三种情况。calcOutageDuration函数要做的事情，就是判断失电负荷点处于哪个区域，进而决定停电时长的取值。这个函数可以拆成三个步骤：
第一步，判断故障元件是否在"电源点到负荷点的唯一供电路径上"。因为辐射型配电网里，电源到每个负荷点的路径是唯一的，故障元件只要不在这条路径上，删掉故障支路根本不会导致负荷失电。这一步其实已经被主循环里的shortestpath查询覆盖了，失电的负荷点，故障元件必然位于其上游路径上。
第二步，判断负荷点在故障点的哪一侧。先搞清楚故障点两侧有没有可操作的隔离开关。方法是：从故障支路的上游端点出发，向电源方向搜索，找到最近的含开关支路；从故障支路的下游端点出发，向负荷方向搜索，找到最近的含开关支路。如果负荷点在"上游隔离开关"和"电源"之间，那它属于第一类（隔离后可恢复）；如果负荷点在故障支路下游，那它属于第二类或第三类。
第三步，对故障支路下游的负荷点，判断是否具备转供条件。方法是从下游隔离开关往下搜索，看有没有联络开关与其他电源点相连。如果存在联络开关且开关闭合后能构成新的供电路径，那它停电时长取tieTime；如果不存在任何联络路径，停电时长取故障元件的修复时间repairTime。
这个搜索逻辑用图论函数实现非常简洁。我写了一个简化版的核心判断代码：
function outageHours = calcOutageDuration(elements, G_temp, faultIdx, loadNode) global isoTime tieTime; % 故障元件的基本信息 n1 = elements(faultIdx, 3); n2 = elements(faultIdx, 4); repairTime = elements(faultIdx, 6); % 如果负荷点可以从电源连通，不需要计算 % （这个判断在主循环里已经做过，这里省略） % 情况1：负荷点在故障点上游（电源侧） % 判断方法：删除故障支路后，如果电源节点连通负荷点，则为上游区域 % 这里由于已经知道失电，说明不连通，所以需进一步判断上游隔离开关位置 % 如果在电源到故障点之间存在隔离开关，则负荷点可在隔离后恢复 if hasUpstreamSwitch(elements, n1) outageHours = isoTime; return; end % 情况2：负荷点在故障点下游，尝试寻找联络开关转供路径 if hasTieSwitchPath(G_temp, n2, loadNode) outageHours = tieTime; return; end % 情况3：无法转供，只能等修复 outageHours = repairTime; end实际工程里，判断hasUpstreamSwitch和hasTieSwitchPath要写不少辅助函数，比如要能够识别"某个节点到电源路径上是否经过开关类型元件""某个故障区域能否通过闭合联络开关连通到备用电源"。这里我建议给自己留一个调试后门——把每个故障元件的元件编号、受影响负荷点、各类时间打印到日志文件里，在开发阶段逐场景核对，远比调试了半天找不到错强得多。
4. 一个完整的小算例：3馈线系统手把手验证
4.1 测试系统参数与期望结果
算法写得再好，不落到具体算例上验证一遍，始终不放心。我自己搭了一个三馈线的测试系统，三条馈线都从同一个110kV变电站的10kV母线引出，每条馈线带2个负荷节点，馈线1和馈线2的末端通过联络开关相连，馈线2和馈线3之间也有一个常开的联络点。我把元件参数列在下面，各位可以照着录入自己的程序：
| 元件 | 类型 | 首端节点 | 末端节点 | 故障率(次/年) | 修复时间(h) | 是否含隔离开关 |
|---|---|---|---|---|---|---|
| E1 | 线路 | 1 | 2 | 0.10 | 5 | 否 |
| E2 | 线路 | 2 | 3 | 0.08 | 5 | 是 |
| E3 | 线路 | 2 | 4 | 0.12 | 5 | 否 |
| E4 | 线路 | 1 | 5 | 0.09 | 5 | 否 |
| E5 | 线路 | 5 | 6 | 0.11 | 5 | 是 |
| E6 | 线路 | 5 | 7 | 0.07 | 5 | 否 |
| E7 | 联络线 | 3 | 8 | 0.01 | 3 | 是（常开） |
| E8 | 联络线 | 7 | 9 | 0.01 | 3 | 是（常开） |
负荷数据：节点4有200户用户、平均负荷400kW；节点3有150户、300kW；节点6有180户、360kW；节点7有120户、240kW。节点编号中8和9是联络线另一侧的备用电源节点，分别由其他馈线供电。隔离时间取0.5小时，转供时间取0.5小时。
4.2 单次故障的推演过程
以元件E5（节点5到节点6的线路）故障为例，手推一遍程序应该给出什么结果。故障发生后，电源节点1到节点6的路径断裂，节点6失电。节点6位于E5下游，检查转供路径：E5下游的节点6，通过E7联络线可以连到节点8，节点8是备用电源，所以节点6的停电时间应该是转供时间0.5小时，而不是修复时间5小时。节点7仍然由原路径供电，因为E5在节点7的上游支路分支点之前，实际经过判断后节点7并没有失去电源连通性，所以不受影响。节点3、节点4也都不受影响。
再看E2（节点2到节点3的线路）故障。节点3失电，且节点3在故障下游。节点3可以通过E6联络线连到节点8，所以停电时间是0.5小时。但如果我把E7的联络线从程序里删掉，节点3就会变成"无转供路径"，停电时间立刻变成5小时修复时间。这就是为什么网架结构对可靠性影响那么大——同样一条线故障，有联络线时SAIDI只增加0.5小时，没联络线时SAIDI就要增加5小时，差别是10倍。
对E1（节点1到节点2的线路）故障，情况比较特殊。节点1是电源，E1故障相当于整条馈线1的电源进线断了，节点2、3、4都失电。节点2在故障点下游，但节点2到节点5之间有没有其他连接？因为我这里没有设置馈线1和馈线2之间的联络，所以节点2、3、4全都无法转供，停电时间均等于E1的修复时间5小时。这说明，辐射型配电网越是靠近电源侧的线路，故障时影响范围越大，这就是所谓"瓶颈元件"——也是规划改造时最值得优先加装自动化开关或联络通道的位置。
4.3 汇总指标与合理性检查
按上面的参数和逻辑，程序跑出来的结果大致是：SAIFI在1.0次/户·年左右，SAIDI在1.5~2.0小时/户·年之间，ASAI在99.98%左右。这个量级和国内许多城市配电网的实际统计值是大体吻合的。
拿到结果之后，我自己会习惯做一个合理性检查：把某一个比较关键元件的故障率调成0，看看SAIFI和SAIDI是不是相应地减少。比如把E1故障率清零，SAIFI应当明显下降，因为E1是馈线主进线，故障影响户数多；如果程序结果显示SAIFI没变化，那一定是拓扑判断逻辑出了bug，或者这个元件根本没被纳入枚举范围。这种"变量扰动检查法"是验证程序正确性的有效手段，比单纯把结果和文献数字对比更能说明问题。
5. 我踩过的坑：常见问题与排查技巧
5.1 结果异常先找这几个原因
调试过程中，我总结出了一个"异常结果排查表"，每次算出来的指标不合理，就按表里的顺序逐项检查：
| 现象 | 可能原因 | 排查方法 |
|---|---|---|
| SAIFI比预期大很多 | 负荷点短暂停电也被重复计次 | 检查是否所有失电负荷点都计了一次频率，是否把"隔离后恢复"的负荷点误加了修复时间 |
| SAIDI比预期小很多 | 把修复时间错误地替换成隔离时间 | 重点检查无联络开关的下游负荷点取的是不是修复时间 |
| EENS为0 | 负荷节点的平均功率数据没有正确读取 | 打印每个负荷节点的功率值，看看是否被初始化为0 |
| 计算时间过长 | 每轮循环都在重复构建最短路径 | 考虑预计算拓扑关系，或用节点邻接表替代graph对象最短路径查询 |
| 结果对某个元件故障率不敏感 | 该元件可能根本没有被枚举进故障集合 | 检查主循环范围，是否漏掉了变压器和联络线等非线路元件 |
还有一个非常隐蔽的问题：开关操作时间与修复时间使用了不同的单位。我曾经在一个工程里把隔离时间定义成分钟，修复时间定义成小时，结果SAIDI莫名其妙小了60倍，查了半天才发现是单位不一致。建议从一开始就统一用"小时"作为时间单位，并在代码注释里大写加粗提醒。
5.2 提高运算速度的几个方法
FMEA本身计算量不大，但当系统规模变大——比如几百条支路、上千个节点时，如果代码写得粗糙，跑起来也会很卡。我有三个提速经验。
第一，不要在主循环里用shortestpath去判断每个负荷点和电源的连通性。更高效的做法是：删除故障支路后，用conncomp函数一次找出所有连通分量，然后查电源节点所在分量，其他分量里的节点就是失电节点。这样一次图遍历能代替对所有负荷点的逐一遍历。
第二，预计算"每个负荷点的供电路径上都有哪些元件"。因为辐射型网络路径唯一，我可以一次性用DFS或shortestpath求出所有负荷点到电源的路径，存成稀疏矩阵。之后判断"某元件故障是否影响某负荷点"，就只需要查这个矩阵的一行，判断故障元件在不在路径集合里，复杂度从O(N)降到O(1)。
第三，多使用向量化操作。比如累计停电时长时，不要在一个循环里一条条累加，而是先记录每个故障场景下各负荷点的停电时长矩阵，最后用sum(A(:))一次性汇总。Matlab的循环开销比较大，能用矩阵运算解决的就别用for。
5.3 避免"看着对，实际错"的校验手段
程序写完之后，最怕的不是报错，而是结果看起来合理、实际内部逻辑有错。我有几种校验手段推荐给大家。
第一种是极小系统手工校验法：构造一个只有一条线路、一个负荷点的最简单系统，手算出SAIFI和SAIDI，再让程序跑一遍。比如线路年故障率0.2次/年，修复时间5小时，负荷点100户。手算SAIFI就是0.2次/户·年，SAIDI就是0.2×5=1.0小时/户·年。如果程序连这个都算不对，那后面复杂系统肯定有问题。
第二种是双实现交叉验证法：我自己写程序时，先用一个最简单的状态搜索方式实现一版"笨方法"，比如完全枚举所有负荷点的失电状态；再用优化后的高效算法实现第二版。两版结果做对比，如果不同，说明优化过程引入了逻辑错误。这个方法虽然笨，但特别稳。
第三种是参数扰动敏感性检查：把某个元件的修复时间从5小时改到10小时，SAIFI应该不变，SAIDI应该变大；把故障率从0.1改到0.01，SAIFI和SAIDI应该近似缩小10倍。如果这些简单的比例关系不成立，那程序内部八成存在错误逻辑。这种检查不会花很多时间，但能帮你提前暴露很多问题。
说到底，配电网可靠性评估程序真正难的地方不在Matlab语法，而在于把电力系统运行逻辑准确地翻译成代码。我刚写的时候也总是想一步到位，结果越急越容易写出一坨逻辑混乱的代码。后来学乖了，先把一条馈线、一个简单开关模型写清楚，再逐步扩展成完整系统，每一步都用上面的校验方法确认没问题后再往前走。这个思路，不管你是做课程设计还是做实际工程，都值得参考。
