做配电网可靠性评估的人应该都经历过那个阶段理论公式背得滚瓜烂熟一到实际算例就不知道从哪下手。尤其是最小路算法原理讲起来三句话但真正要把一个IEEE标准测试系统的所有参数整理好、把拓扑和算法对到代码里去跑出可信的结果起码要折腾好几天。这也是我最初整理这套Matlab代码的动机——把配电网可靠性评估里最常用的IEEE RBTSRoy Billinton Test System测试系统参数和最小路算法封装在一起做到直接运行、直接出指标。这套代码后来帮了不少做可靠性规划、写论文、做课程设计的人。不管你是刚接触可靠性评估的研究生还是在电科院、设计院需要快速评估某个配电网方案可靠性的工程师只要跑通这份代码掌握最小路算法的基础逻辑框架后面再改拓扑、加分布式电源或者调整保护策略都是在现有骨架上扩展不用再经历从零开始的痛苦。这篇文章我会把这套代码的几个关键部分重新拆开讲为什么选最小路法和IEEE RBTS系统作为验证平台、系统参数怎么结构化组织、算法每一步在算什么、代码跑通之后怎么核对结果是否正确以及我调这个程序时踩过的几个坑。1. 项目整体设计与方案选型1.1 为什么用解析法而不是蒙特卡洛模拟配电网可靠性评估的方法大致分两派解析法和蒙特卡洛模拟法。解析法的核心思路是把系统的故障事件逐一枚举然后通过概率公式叠加优点是计算速度快、结果确定同一套参数跑多少次结果都一样而且每个指标都能反向追溯到是哪些元件贡献的。Minimal path algorithm就属于这类方法尤其适合配电网这种结构相对简单、保护配合关系明确的辐射状网络。蒙特卡洛模拟法能建模更复杂的时序场景比如分布式电源随机出力、电动汽车负荷波动但代价是收敛慢通常要跑几千甚至上万个样本才能让指标稳定下来对初学者并不友好——你很难判断偏差是算法错了还是样本不够。在这个项目里我选择的是解析法里的最小路法而不是更复杂的故障模式后效分析FMEA或网络等值法。原因很简单FMEA需要把所有故障事件逐条列出来网络一复杂就爆炸网络等值法的思路虽然巧妙但对新手不直观调试起来很头疼。最小路法则刚好卡在“逻辑直观”和“实现简单”的平衡点上对每个负荷点找一条主路径然后分而治之代码量小结果还能用手算验证。1.2 为什么拿IEEE RBTS做验证平台IEEE RBTS是Roy Billinton Test System的缩写由加拿大Saskatchewan大学Billinton团队提出后来被电力系统可靠性领域广泛采用国内很多论文里的“IEEE RTS系统”或者被写错成“RTBS”的指的基本都是这套标准测试系统。它的价值在于所有元件的可靠性参数、系统拓扑、负荷数据都是公开的全世界研究者用同一套数据做研究结果可以互相比较。这套系统有几个版本。最早的是6母线可靠性测试系统RBTS主要面向发输电和配电网两个层面后来扩展到Bus 2、Bus 4等多个配电网馈线系统再往后还有更大规模的RTS-79、RTS-96主要用在发输电可靠性评估。配电网可靠性研究里最常碰到的就是RBTS的配电馈线部分比如Bus 2系统、Bus 4系统——它们包含了好几种典型的配电网络结构有主干线、分支线、隔离刀闸、熔断器、备用电源非常适合演练最小路算法。我在代码里内置的就是RBTS配电馈线系统的参数这样使用者不需要再去翻原始文献找数据拿到代码就能跑。更重要的是因为原始文献里已经给出了标准结果跑出来的指标可以直接对标这是很多自建算例做不到的。1.3 用Matlab承载这套代码的原因选Matlab而不是C或Python主要是三类原因。第一Matlab对矩阵和稀疏矩阵的支持非常成熟配电网的拓扑天然适合用邻接矩阵和稀疏矩阵表达做路径搜索、逻辑索引都很方便。第二Matlab内置了graph对象和shortestpath函数最小路搜索不用自己从头写复杂的图算法两行代码就能搞定。第三Matlab的调试体验好运行过程中能随时在工作区查看中间变量——这对学习和验证算法太重要了你能看到每个负荷点的最小路节点序列、每个元件的分类结果而不是一个黑盒输出最终数值。另外Matlab画图方便跑完之后想输出负荷点故障率柱状图、SAIDI构成饼图都是几行代码的事这对写报告和论文是刚需。2. 系统参数与数据准备2.1 RBTS系统结构速览代码包内置的RBTS配电系统主干结构大致如下一个11kV电源母线向下引出多条主馈线每条馈线上挂着多个负荷点馈线之间在需要的位置通过联络开关连接负荷点附近配置熔断器和隔离开关。这里以代码包中常用的馈线系统为例其拓扑结构和标准文献保持一致。元件类型故障率参数修复/操作时间备注线路每千米0.065次/年修复时间5小时单位长度故障率计算时需要乘线路长度配电变压器0.015次/年修复时间200小时变压器故障后通常直接更换或维修隔离开关操作时间0.5小时隔离故障段故障隔离不影响非故障段转供熔断器可靠动作概率0.95动作后切换时间约1小时保护分支线要特别注意这些是IEEE RBTS标准文献里的典型参数不同版本会有微调。我写的代码包里以data_IEEE_RBTS.m文件里的数值为准文件名本身就提示了参数来源。如果你参考的是其他论文里的参数比较结果前一定要先核对参数是否一致否则指标对不上很正常。2.2 可靠性参数建模的细节整理参数时有一个特别容易踩的坑故障率单位。线路参数0.065次/年这是“每千米每年”的故障率实际计算必须乘上线路长度。比如一段2.5km的线路年故障率是0.065 × 2.5 0.1625次/年。如果你把0.065直接当成整条线路故障率结果会偏差很大而且这个偏差在最小路法里会被串联累加放大最终SAIFI很可能只有标准值的一半不到。另一个细节是修复时间和操作时间必须分开。修复时间对应元件损坏后需要维修或更换的时间通常以小时计比如变压器200小时操作时间对应通过隔离开关、联络开关进行隔离和转供的时间通常只有0.5到1小时。在最小路算法里最小路上元件故障通常按修复时间算而最小路外支路故障如果能够通过开关操作转供按操作时间算。混用这两类时间是新手最容易让结果偏离的地方。此外RBTS原始文献并没有给全所有开关的操作时间比如联络开关的倒闸时间。代码里我按工程惯例取了一个经验值0.5小时并在注释里标明了“经验值可修改”。这种参数在标准文献里缺失很正常关键是代码使用者要知道这类参数是假设值而不是文献原始数据。2.3 代码里如何组织这些参数我见过不少人把系统参数散写在多个脚本里最后自己都找不到哪里改参数。这套代码的做法是建立一个单独的数据文件data_IEEE_RBTS.m所有参数集中管理。%% IEEE RBTS配电系统参数定义部分展示 % line: [from_node, to_node, length_km, fault_rate_per_km_yr, repair_time_h] lines [ 1 2 2.0 0.065 5.0; 2 3 1.5 0.065 5.0; 3 4 2.0 0.065 5.0; 4 5 2.5 0.065 5.0; ]; % transformer: [start_node, end_node, fault_rate_per_yr, repair_time_h] tf [ 5 6 0.015 200.0; ]; % load point: [node, user_count, average_load_kW] loads [ 3 200 350.5; 4 150 280.0; 5 180 400.0; 6 120 180.0; ]; % switch: [switch_node, operate_time_h, can_transfer] switches [ 4 0.5 1; 5 0.5 0; ];矩阵行对应元件列对应参数好处是后续算法里可以直接用矩阵运算批量处理不需要写一堆for循环遍历结构体。如果后面对参数做灵敏度分析改一下矩阵对应位置的数值就行。如果你想扩展系统拓扑比如在节点3和5之间加一条联络线只需要在lines矩阵里加一行同时把switches里相应的转供标志改一下。代码里的最小路搜索是自动识别拓扑的不需要改算法部分。3. 最小路算法核心原理与实现3.1 什么是最小路最小路又叫最小路径指从电源点到某个负荷点的供电路径。在辐射状配电网里每个负荷点的最小路其实就对应从根节点到该负荷点沿着馈线走的那一串元件。比如一个负荷点挂在馈线末端那么从电源母线到它经过的那些线段、变压器、开关就是这个负荷点的最小路上元件。可以打个比方快递从仓库送到一个小区主干道上的每一个路口都可能影响整个小区的收货而小区内部某个楼栋的支路问题只影响那一栋楼。最小路算法就是把这个思路量化主干道上的元件故障负荷点必须等修复支路上的元件故障可能通过熔断、隔离、倒闸快速恢复停电时间就短很多。这里要注意“最小”并不是“最短距离”。配电网里最小路通常就是唯一的供电路径不需要求最短路径。只是在存在多电源、多联络的情况下路径选择需要结合网络拓扑和运行方式判断。代码里我用的是Matlab自带的graph对象和shortestpath函数本质上是在网络图上做可达性搜索找到从电源节点到负荷节点的路径。3.2 元件分类与停电时间计算逻辑对于每个负荷点最小路算法把所有元件分成两类。第一类是最小路径上的元件这些元件与负荷点是串联关系任何一个故障都会直接造成负荷点停电直到修复或更换完毕。设某元件故障率为λ_j、修复时间为r_j则这个元件对负荷点故障率和年停电时间的贡献分别是λ_j和λ_j × r_j。对所有最小路上的元件求和就得到λ_path Σλ_j U_path Σ(λ_j × r_j)第二类是最小路径以外的元件。这些元件故障时负荷点会不会停电取决于保护和开关配置。如果故障支路上有熔断器或隔离开关而且系统可以通过联络开关把负荷转到备用电源那负荷点实际停电时间只是故障隔离和切换的操作时间t_s如果网络根本没有转供通道那负荷点也只能等故障修复这个元件的贡献就按修复时间计算。这就是最小路算法最核心的“分而治之”思想λ_i λ_path λ_off_path U_i U_path U_off_path r_i U_i / λ_i在代码里我通过一个分类标志矩阵来判断每个元件属于哪一类同时给每个元件配置了一个“是否可转供”属性。元件的故障率照常累加但停电时间项根据分类和可转供属性选择乘修复时间还是操作时间。这一层逻辑是整个程序的灵魂后面调bug时也基本都集中在这。3.3 系统指标计算公式负荷点指标算出来之后还要汇总成系统级指标常用的有这几个指标公式含义SAIFIΣ(λ_i × N_i) / ΣN_i系统平均停电频率次/户·年SAIDIΣ(U_i × N_i) / ΣN_i系统平均停电持续时间小时/户·年CAIDISAIDI / SAIFI用户每次停电平均持续时间小时/次ASAI(8760 - Σ(U_i × N_i)) / 8760供电可用率无量纲ENSΣ(U_i × L_i)期望缺供电量千瓦时/年AENSENS / ΣN_i平均每户缺供电量千瓦时/户·年需要提醒的是N_i是负荷点的用户数L_i是该负荷点的平均负荷功率这两个数据的口径在RBTS原文里有明确值。如果用MW还是kW没统一后面ENS会差三个数量级。我的建议是代码里统一用kW和小时结果直接就是kWh。3.4 代码实现关键点最小路搜索和可靠性计算的核心函数大概长这样function [lambda_i, U_i, r_i] min_path_eval(sysData, loadNode) % sysData: 系统参数结构体包含lines, loads等 % loadNode: 负荷点节点编号 % 1. 构建图对象 G graph(sysData.lines(:,1), sysData.lines(:,2)); % 2. 搜索电源节点到负荷节点的最小路 pathNodes shortestpath(G, sysData.sourceNode, loadNode); pathEdgeMask ismember(sysData.lines(:,1), pathNodes(1:end-1)) ... ismember(sysData.lines(:,2), pathNodes(2:end)); pathEdgeIdx find(pathEdgeMask); % 3. 初始化 lambda_i 0; U_i 0; % 4. 最小路径上的元件按故障率累加停电时间按修复时间 for k 1:length(pathEdgeIdx) e pathEdgeIdx(k); lambda_j sysData.lines(e,3) * sysData.lines(e,4); % 长度×单位故障率 r_j sysData.lines(e,5); % 修复时间 lambda_i lambda_i lambda_j; U_i U_i lambda_j * r_j; end % 5. 最小路径外的元件判断影响范围按操作时间或修复时间 offPath setdiff(1:size(sysData.lines,1), pathEdgeIdx); for k 1:length(offPath) e offPath(k); % 判断该元件故障是否影响当前负荷点 if isAffected(sysData, e, loadNode) lambda_j sysData.lines(e,3) * sysData.lines(e,4); if sysData.lines(e,6) 1 % 可以转供 U_i U_i lambda_j * sysData.swOperTime; else % 不能转供 U_i U_i lambda_j * sysData.lines(e,5); end lambda_i lambda_i lambda_j; end end % 6. 计算平均持续时间 r_i U_i / lambda_i; end这段代码是简化逻辑实际里isAffected函数还要判断保护配合和隔离区域但核心骨架就是这样。要注意的是步骤5里那个“是否受影响”的判断非常敏感判断错了结果偏差极大。下一节我会展开讲这个坑。4. 代码运行与结果验证4.1 目录结构与运行方式代码包按功能拆成了几个文件结构很清楚文件作用main_Reliability.m主程序入口依次调用数据读取、算法计算、结果输出data_IEEE_RBTS.m定义系统拓扑和可靠性参数min_path_eval.m单个负荷点的最小路可靠性计算min_path_search.m最小路搜索基于graph对象实现isAffected.m判断某元件故障是否影响指定负荷点plot_results.m生成负荷点指标和系统指标的可视化图表export_results.m将结果导出为Excel或表格数据运行方式很简单在Matlab里打开main_Reliability.m直接点运行按钮或者在命令行输入main_Reliability。程序会按顺序读取参数、遍历所有负荷点计算指标、汇总系统指标最后在工作区输出一个result结构体并自动生成几个图表。需要注意数据文件data_IEEE_RBTS.m是一个脚本而不是函数运行时会把变量直接导入工作区。如果你要修改系统参数只改这个文件就够了。不建议把参数散落到算法文件中否则后面会非常难维护。4.2 输出结果与标准值对标程序运行后你会得到每个负荷点的故障率λ、年平均停电时间U、平均故障持续时间r以及系统级的SAIFI、SAIDI、CAIDI、ASAI和ENS。这里给一个结果结构示意具体数值以本代码包实际运行为准负荷点故障率λ(次/年)停电时间U(小时/年)平均持续时间r(小时/次)LP1约0.23约1.2约5.2LP2约0.35约1.8约5.1LP3约0.29约1.5约5.2系统SAIFI约0.9次/户·年SAIDI约2.6小时/户·年ASAI约0.9997判断结果对不对最直接的办法就是和IEEE RBTS原始文献里的标准值对比。如果偏差超过10%优先检查三个方面线路故障率有没有乘长度、最小路外的元件是否错误地按修复时间计算、用户数N_i有没有和负荷点对应错位。这步验证环节千万别省。我见过不少人在自己的算例上跑出一个数字就觉得万事大吉结果后来发现跟标准算法结果差了两倍最后定位到是参数单位问题。用标准测试系统的好处就是有一个“标准答案”帮你兜底。4.3 结果可视化与成果输出跑完计算后plot_results.m会生成几张图各负荷点故障率柱状图、系统SAIDI构成饼图、负荷点停电时间对比图。这些图对论文和报告很实用。如果你要导出Excel用export_results.mT table(loadPointNames, lambda_i, U_i, r_i, ... VariableNames, {负荷点, 故障率次每年, 年停电时间小时, 平均持续时间小时每次}); writetable(T, reliability_results.xlsx);这个导出功能看起来不起眼但实际做报告时特别好用。不需要手工从命令行复制数据改完参数重跑一遍直接就能拿到一份整理好的结果表。5. 常见问题与调试经验5.1 拓扑编号错位导致路径搜索出错最小路搜索依赖节点编号一旦节点编号在参数矩阵里不一致找出的路径就会莫名多出一截或者少一截最后算出来的指标完全对不上。排查方法是在main_Reliability.m里加一行输出把所有负荷点的最小路节点序列打印出来人工核对一遍。比如负荷点LP2应该在节点1-2-3-4这样的序列上如果打印出来变成了1-2-5-4那一定是节点5和4之间的线路编号写错了。代码包里我预留了一个debug开关打开后会在命令行打印路径节点序列。5.2 非最小路元件的影响范围判断错误这是整个算法里最容易出错、也最难排查的地方我单独拿出来讲。很多初学最小路法的人要么把所有非最小路元件都算成“通过操作时间转供”导致结果偏小要么把所有非最小路元件都算成“停电直到修复”导致结果偏大。正确逻辑必须结合网络结构判断如果故障元件所在的支路与负荷点之间没有隔离开关或熔断器那这个支路故障时负荷点无法被隔离只能等修复这时候就必须按修复时间算如果故障支路可以通过熔断器隔离并且存在联络开关转供通道才可以按操作时间算。举个例子一个纯辐射状馈线没有联络开关那么非最小路上的支路故障下游负荷点根本不可能转供只能停电直到修复。这种情况下把所有非最小路元件都按0.5小时操作时间算SAIDI会明显偏低结果跟实际相差很大。我的建议是在代码里给每段线路、每台变压器加一个“可转供”标志而不是统一用全局开关时间。这样每个元件的停运逻辑都能单独控制结果也更可控。isAffected函数里就是通过这个标志决定用哪个时间的。5.3 参数单位与统计口径不一致这块前面也提到过但值得再强调一次因为它实在太容易踩了。故障率的单位有三种常见写法次/年整条元件、次/km·年每公里、次/百km·年某些文献。如果你从文献里抄参数没注意单位结果很可能差一个数量级。修复时间也有小时、分钟两种写法用混了SAIDI会差60倍。我的经验是所有参数统一在data_IEEE_RBTS.m里转换成“次/年”和“小时”这两个基准单位后续算法里再也不用管单位换算。比如线路参数原始单位是次/km·年进入数据文件前就先乘以长度转成次/年这样算法代码里没有任何单位换算逻辑减少出错机会。5.4 循环慢和向量化最小路算法本身计算量不大普通RBTS几十个元件for循环几秒就跑完。但如果你后面把系统扩展到几百个节点、上千个元件纯for循环就会开始吃力。建议提前做好三件事一是给所有矩阵变量预分配空间不要边循环边扩容二是用Matlab的graph对象做路径搜索而不是自己写Dijkstra三是对“最小路外元件”的批量判定尽量用逻辑索引和元胞数组操作避免每个负荷点都套两层for循环遍历所有元件。Matlab 2022b之后graph和shortestpath的性能提升非常明显如果你用的是新版本不需要自己优化路径搜索代码直接用内置函数就好。5.5 快速核对结果的技巧最后分享一个我调试这类程序时的习惯不要等全部算完再去对总体指标而是先拿一个最简单的负荷点手算一遍再和程序输出对比。比如负荷点LP1如果它的上游就三段线路和一个变压器那手算很简单把三段线路的故障率相加再乘以各自的修复时间加上变压器的故障率乘以200小时得到λ和U然后r U / λ。如果程序输出的比手算的偏差超过0.01说明基础逻辑有问题先不要看系统指标修好单个负荷点再说。这个方法我每次换系统、改参数、重构代码时都会用能省掉后面成倍的排查时间。6. 扩展方向从最小路法到更复杂的场景代码跑通之后很多人的下一步需求是扩展。比如在配电网里加入分布式光伏、储能然后在可靠性评估里体现它们的作用。这时最小路法的基本框架仍然可以用光伏和储能相当于给负荷点增加了一个备用电源会让非最小路元件的“可转供”属性发生改变。你需要做的只是在isAffected函数里增加一个分支如果负荷点附近有分布式电源且孤岛运行策略允许那么某些原本不能转供的故障就可以通过孤岛运行恢复供电。再比如多联络开关的优化问题也是在这套代码基础上做的。最小路法算出来的每个负荷点的停电时间本质上就是评估不同联络开关配置下的可靠性差异你把开关配置变成一个可调参数跑一个循环就能做简单的可靠性比选。这些扩展的工程价值比单纯跑通一个算法大得多。我在实际使用这套代码时最深的体会是可靠性评估的核心不是“会跑代码”而是能正确判断每个元件故障后系统的真实响应过程最小路法只是帮我们把这种判断过程结构化、自动化了。参数可以换、拓扑可以改、指标可以加但“元件分类—时间赋值—指标叠加”这个逻辑骨架是通用的这也是这套代码真正能沉淀下来的原因。