做配电网可靠性评估的人,应该都经历过那个阶段:理论公式背得滚瓜烂熟,一到实际算例就不知道从哪下手。尤其是最小路算法,原理讲起来三句话,但真正要把一个IEEE标准测试系统的所有参数整理好、把拓扑和算法对到代码里去跑出可信的结果,起码要折腾好几天。这也是我最初整理这套Matlab代码的动机——把配电网可靠性评估里最常用的IEEE RBTS(Roy 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 | 系统平均停电持续时间,小时/户·年 |
| CAIDI | SAIDI / SAIFI | 用户每次停电平均持续时间,小时/次 |
| ASAI | (8760 - Σ(U_i × N_i)) / 8760 | 供电可用率,无量纲 |
| ENS | Σ(U_i × L_i) | 期望缺供电量,千瓦时/年 |
| AENS | ENS / Σ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.2 |
| LP2 | 约0.35 | 约1.8 | 约5.1 |
| LP3 | 约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.m:
T = 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函数里增加一个分支:如果负荷点附近有分布式电源,且孤岛运行策略允许,那么某些原本不能转供的故障,就可以通过孤岛运行恢复供电。
再比如多联络开关的优化问题,也是在这套代码基础上做的。最小路法算出来的每个负荷点的停电时间,本质上就是评估不同联络开关配置下的可靠性差异,你把开关配置变成一个可调参数,跑一个循环,就能做简单的可靠性比选。这些扩展的工程价值比单纯跑通一个算法大得多。
我在实际使用这套代码时最深的体会是:可靠性评估的核心不是“会跑代码”,而是能正确判断每个元件故障后系统的真实响应过程,最小路法只是帮我们把这种判断过程结构化、自动化了。参数可以换、拓扑可以改、指标可以加,但“元件分类—时间赋值—指标叠加”这个逻辑骨架是通用的,这也是这套代码真正能沉淀下来的原因。