简介面向配电网可靠性评估的MATLAB实现基于最小路法与蒙特卡洛算法适合电气工程专业学生、研究人员及相关工程技术人员使用可用于配电网规划、运行方式分析等场景。压缩包约1.34MB以PDF参数文档、Excel原始数据表和MATLAB程序文件为主要类型涵盖IEEE RBTS系统原始参数、IEEE33节点系统参数以及IEEE RBTS BUS6的matlab数据。资源提供两套完整的评估程序最小路算法程序基于网络拓扑搜索完成可靠性指标计算非序贯蒙特卡洛主程序通过节点影响分析法判断受影响负荷再结合随机模拟完成评估。资源按参数准备、算法实现与主程序三部分组织便于按需查阅两种方法互为补充便于读者对比算法差异、学习可靠性指标计算流程与程序实现思路。已有3816人学习下载适合Matlab编程基础较好、需要完成课程设计或课题验证的读者。1. 配电网可靠性评估为什么绕不开最小路法和蒙特卡洛法先分清“算得准”和“算得动”配电网可靠性评估这件事最让人难受的不是不会算而是同一个辐射网换一种算法结果差出两倍两边还都能自圆其说。最小路法和蒙特卡洛法是这个领域最常被摆在一起的两条技术路线一个把拓扑拆成串联元件用故障率直接叠加一个把元件状态当随机变量靠大量抽样逼近长期平均表现。做规划改造要结果做运行分析要过程下面按两条路线各自怎么落地、参数怎么设、坑在哪展开适合正在做配网网架改造、供电可靠性指标测算的工程师参考。2. 最小路法先求路径再叠加指标用 DFS 把拓扑变成可靠性表达式2.1 最小路的本质负荷点只关心从电源到它的必经元件辐射状配电网里某个负荷点的电是从变电站母线出发沿着馈线一路穿过多级线路段、分段开关、配变才进到用户表计。这一路上任何元件故障在没有任何备用电源的场景下负荷点就直接失电。所以这个负荷点的可靠性本质上由这条“从电源到负荷点”的通路上所有元件的串联可靠性决定。所谓最小路就是从电源到负荷点的不含重复节点的简单路径。在纯辐射网中它往往只有一条叫“最小”是为了和闭环运行、多电源转供时存在的多条候选路径区分开。有了最小路计算就退化成串联系统的可靠性折算。设通路上第 i 个元件的年故障率为 λ_i次/年单次故障修复时间为 r_i小时/次那么负荷点年故障频率 λ_LP Σλ_i年停电小时数 U_LP Σ λ_i r_i可用度 A_LP Π A_i工程上常近似写成 A 1 - U/8760。这套公式看起来简单真正容易出错的是“最小路之外的那批元件”。它们不在负荷点的必经通路上但支线故障可能通过保护配合影响上游处理不好就会把指标算虚高。这也是为什么最小路法的实现重点不是求路径本身而是路径之外的影响折算。2.2 用 DFS 求最小路径一份可直接套用的邻接表实现求最小路最常见的是深度优先搜索DFS把馈线拓扑建成长邻接表从电源节点出发一路搜到负荷点。下面是按邻接表实现的最小路搜索辐射网和开环运行的网架都能直接套。def find_minimal_paths(graph, source, target): graph: {节点: [相邻节点, ...]}返回所有不含重复节点的路径。 paths [] def dfs(node, path): if node target: paths.append(path[:]) return for nxt in graph.get(node, []): if nxt not in path: # 防止绕圈不重复经过节点 path.append(nxt) dfs(nxt, path) path.pop() dfs(source, [source]) return paths这段代码的逻辑很直白从 source 出发每到一个新节点就继续向下扩展直到到达 target 才把当前路径存档。nxt not in path是唯一的循环保护保证路径里不会出现重复节点也就是我们说的“最小路”的数学定义。参数上只有三个对象graph 是节点邻接表source 是变电站低压母线节点target 是负荷点节点。辐射网通常每个负荷点只返回一条路径如果网架里有联络开关处在合闸位置、形成闭环这段代码会返回多条路径这时要以正常运行方式为准把合闸的联络开关在建模阶段先断开。路径拿到手后要折算成串联可靠性指标下面的函数把一条路径上所有支路的 λ 和 r 累加起来def path_reliability(path, comps): comps: {(u,v): {lam: 次/年, r: 小时/次}}按串联模型累加。 lam_sum 0.0 u_sum 0.0 for i in range(len(path) - 1): key (path[i], path[i 1]) lam comps[key][lam] r comps[key][r] lam_sum lam u_sum lam * r # 单位小时/年 return lam_sum, u_sum逻辑说明串联系统里故障频率相加年停电小时数也相加这是两状态模型下最基本的折算。参数说明里最需要注意的是 r 用“小时/次”而不是“天/次”否则 U 会整体差 24 倍。另外如果馈线装有重合闸并且重合成功率为 p应该把 λ 先乘上 (1-p) 再做累加因为一次永久性故障才会真正造成用户停电。2.3 最小路之外的元件怎么折算分支贡献与 FMEA 累加辐射网里大部分元件其实不在某个负荷点的最小路上但它们故障时可能通过“越级”影响上游负荷点。常见处理是以故障模式影响分析FMEA为核心把每条支路可能造成的停电影响列成一张表再逐负荷点聚合。下面是一张三支路馈线的典型影响表。故障元件所在位置影响的负荷点停电时间(h)L2主干LP1、LP2r_L2 3L4分支出口有熔断器LP2r_L4 2L4熔断器拒动分支LP1、LP2r_L4 上级开关整定 0.5逻辑说明L2 在主干上故障只能靠上游开关切除所以两个负荷点都停到修复完L4 在分支且出口熔断器正常动作时只有 LP2 停电LP1 不受影响一旦熔断器拒动故障就要靠上级分段开关切除LP1 也会被拖进去。这张表每个元件只出现一行但“影响范围”一列随保护配合变化聚合时对每个负荷点把相关行的 λ 和停电时间相乘求和就得到该负荷点的年停电小时数。注意最小路法对故障影响表的依赖远大于对求路径的依赖。程序写得再快这张表漏了一行“拒动”结果照样失真先核对表再谈性能。3. 蒙特卡洛法当解析规则失效时用随机抽样把缺供电时间“数”出来3.1 序贯与非序贯蒙特卡洛按你的指标粒度选型解析法的前提是各元件独立、故障率和修复时间是常数这对规划阶段的年均指标够用但一碰到时序负荷曲线、天气相关的故障率波动、检修计划和转供策略解析表达式就会变成黑匣子里的玄学。蒙特卡洛法换了个思路按照元件的故障率和修复时间抽样把一年的运行状态一遍遍模拟出来最后对上千个模拟年的结果求平均。蒙特卡洛内部还分两派。非序贯抽样状态抽样每个模拟年只抽一次元件的“完好/故障”状态适合算 LOLP、EENS 这类概率型指标序贯抽样则按时间轴排故障和修复事件能自然给出停电持续时间、停电频次分布SAIDI 和 CAIDI 的估计就要靠它。选型的标准很简单只要年平均指标非序贯就够了要看“停电多久、一年几次”这些过程性指标老老实实写序贯。3.2 非序贯抽样最小代码状态抽样加连通性检查非序贯抽样的核心是两件事按每个元件的稳态可用度抽状态然后判断负荷点是否与电源保持连通。下面这份代码把这两件事都做了。import random def reachable(states, graph, source, target): states: {(u,v): True/False}沿带电支路做深度优先搜索。 stack [source] seen set() while stack: node stack.pop() if node target: return True for nxt in graph.get(node, []): if states.get((node, nxt), True) and nxt not in seen: seen.add(nxt) stack.append(nxt) return False def mc_non_sequential(branches, graph, load_points, source, samples20000, seed2024): random.seed(seed) lp_keys list(load_points.keys()) prob_out {k: 0.0 for k in lp_keys} for _ in range(samples): states {} for bid, br in branches.items(): lam br[lam] # 次/年 r br[r] # 小时/次 p_up 8760.0 / (8760.0 lam * r) # 稳态可用度 st random.random() p_up states[(br[u], br[v])] st states[(br[v], br[u])] st # 无向边两个方向同一状态 for k in lp_keys: if not reachable(states, graph, source, k): prob_out[k] 1.0 for k in lp_keys: prob_out[k] / samples print(fLOLP {k}: {prob_out[k]:.5f}, fEENS {prob_out[k] * 8760.0 * load_points[k]:.4f} MWh/yr)逻辑说明每个模拟年对每条支路独立抽样p_up 是两状态马尔可夫模型的稳态可用度推导过程是 年可用小时数(8760-λr) 除以 8760。抽样后把支路状态写进 states且正反两个方向都写保证从任意一端都能遍历到。reachable从电源节点做 DFS只要有一条带电通路能到负荷点就算供电正常。参数说明branches 是支路字典lam单位是次/年r单位是小时/次load_points 给的是每个负荷点的年平均负荷MWEENS 用失电概率乘 8760 小时再乘负荷得到单位 MWh/年。samples默认两万次实际够不够要用下一章的方差系数判据说话不要拍脑袋定。3.3 序贯抽样的核心循环把 TTF/TTR 排成事件时间线序贯蒙特卡洛的骨架是把每条支路的下一次故障时刻排成一张事件表取最早事件处理处理完再排下一次。下面的代码演示这个循环产出的是按时间排序的故障-修复事件时间线。import math import random def sequential_timeline(branches, horizon_years, seed7): random.seed(seed) end horizon_years * 8760.0 # 模拟总时长单位小时 next_t {} for bid, br in branches.items(): lam br[lam] # 次/年 ttf_h -math.log(1 - random.random()) * 8760.0 / lam next_t[bid] ttf_h timeline [] while True: bid min(next_t, keynext_t.get) t next_t[bid] if t end: break repair_over t branches[bid][r] # 固定修复时间 timeline.append((bid, t, repair_over)) ttf_h -math.log(1 - random.random()) * 8760.0 / branches[bid][lam] next_t[bid] repair_over ttf_h # 修好后再开始下一次寿命 return timeline逻辑说明指数分布抽样公式是 TTF -ln(1-u)/λ注意这里 λ 以“次/年”为单位算出来是年要乘 8760 换成小时。每次取事件表中最小时刻的支路记录故障开始和修复完成两个时间点然后把该支路下一次故障时刻更新为修复完成之后的再失效时间。参数说明branches 结构和非序贯版本一致r在这里按“小时/次”直接参与时间轴。这段代码是骨架真正的工程版本还需要把 timeline 接到拓扑分析上——每个故障时间窗内判断负荷点能否通过联络开关转供、转供需要多长操作时间然后把停电时长累积到 SAIDI。修复时间也可以换成对数正态抽样只要把repair_over t branches[bid][r]改成从对数正态分布随机采样即可方差大会让收敛变慢这点后面专门讲。4. 参数与收敛判据抽样次数、故障率和修复时间到底怎么定4.1 SAIFI/SAIDI/EENS 指标与算法的对应关系选算法前先明确要输出什么指标两者不是互相替换的关系而是各管一段。下表是我常用的对应关系。指标含义最小路法蒙特卡洛法SAIFI系统平均停电频率次/用户·年直接累加 λ 再按用户加权序贯法可以非序贯只能算概率SAIDI系统平均停电持续时间h/用户·年累加 λr 再按用户加权序贯法最准CAIDISAIDI / SAIFI平均每次停电时长派生计算序贯法天然输出EENS期望缺供电量MWh/年需负荷-故障矩阵非序贯乘 8760 小时即可停电频次分布每户每年停多少次很难给出分布序贯法直接统计说明几句规划阶段要快速比选网架方案最小路法几分钟出一版 SAIFI/SAIDI够用运行风险分析要看到“哪些用户停得久”必须上序贯蒙特卡洛。特别提醒别拿非序贯抽样去算 SAIFI它给的是失电概率不是停电频率概率除以时间不是频率这个转换在数学上不成立。4.2 收敛判据与抽样次数方差系数 β 的实用阈值蒙特卡洛的结果是随机估计量不能跑完固定次数就交差得看估计值的相对波动。工程上通用判据是方差系数 β σ/μ建议 β 小于 0.05 才停止抽样。对 LOLP 这类伯努利估计量β 有解析表达式。import math def variance_coefficient(p_lolp, n): LOLP 估计的方差系数建议 β 0.05 再停止抽样。 if p_lolp 0.0: return 0.0 return math.sqrt((1.0 - p_lolp) / (n * p_lolp))逻辑说明n 是已完成的抽样次数p_lolp 是当前算出的失电概率估计值。β 反比于根号下 n 与 p 的乘积所以 p 越小需要的 n 越大。参数经验当 p 在 1e-3 量级时β0.05 对应的 n 大约 40 万次p 到 1e-4 量级就要 400 万次。现实中不少项目跑一万次就收工LOLP 还在抖动这不是代码问题是样本量数学上就不够。实际做法是每两万次算一次 β满足阈值就停顺便把上一次的均值留作对照。提示β 只反映随机抽样误差不反映模型误差。断路器拒动率没建模β 再小也算不对。4.3 可靠性参数从哪里来故障率、修复时间、开关操作时间参数是可靠性评估的粮食。配电网常见元件的典型参数范围下面这张表基本够起步用。元件 / 场景λ 典型值修复或操作时间架空裸导线0.1 ~ 0.3 次/(km·年)3 ~ 6 h/次交联电缆0.03 ~ 0.08 次/(km·年)6 ~ 12 h/次定位困难配电变压器0.005 ~ 0.02 次/台·年4 ~ 8 h/次断路器拒动率0.5% ~ 2%每次故障不需要自动化分段开关操作不计故障率0.1 ~ 0.3 h/次人工操作联络/分段开关不计故障率1 ~ 2 h/次参数说明有两个关键点。第一λ 的单位必须跟元件粒度匹配线路按“每公里每年”算之前把长度乘进去变压器按“每台每年”开关一般只计自身故障率约 0.01~0.03 次/台·年不要把它的人工操作时间混进 λ。第二开关操作时间直接决定转供策略的价值——自动化开关能把主干停电时间从“修复时间”压缩到“操作时间”这是网架改造里 SAIDI 下降的主要来源算方案时必须把“可转供”和“不可转供”两类负荷点分开处理。5. 配电网可靠性评估避坑五个最容易被算错的环节按最小路法和蒙特卡洛法做过几轮真项目之后会发现出错点其实很集中。下面五条是我踩过或者见别人踩过的坑每条按现象、原因、解决三步拆开。5.1 分支元件的故障影响被重复叠加进上游负荷点现象一条带多个分支的馈线最小路法算出来的 SAIDI 比实际运行统计高 40% 以上第一反应是故障率给大了但对账下来发现不是。原因把不在负荷点最小路上的分支元件也按“串联”直接加进上游负荷点的 λ等于假设每次分支故障都会拖累主干完全忽略了熔断器和分段开关的隔离作用。解决先做故障影响分类——分支出口有熔断器且正常动作时分支故障只影响支线自身负荷点上游只计熔断器动作的短暂过程没有熔断器时才按上游开关动作时间计入。代码实现上就是把分支元件的 r 换成 min(r, 开关操作时间)并只在对应的负荷点上累加。5.2 蒙特卡洛收敛慢一味加样本数没有用现象抽样从 1 万次加到 10 万次LOLP 还在 1e-3 量级抖动算力烧了不少结果还是不稳。原因配电网元件故障率很低失效事件本身就是稀有事件普通抽样下相对方差与 1/p 成正比p 越小收敛越慢加样本是线性收益治标不治本。解决先用 4.2 的 β 判据确认是不是真没收敛再用重要抽样把失效概率抬高、统计时再按似然比还原如果是两个方案做对比用公共随机数同一个种子消除两边的共同方差收敛速度能快一个量级。5.3 故障率单位不统一路径归并后指标整体偏大或偏小现象同一份台账里架空线有的按“条”写 λ0.2有的按“公里”写 λ0.2电缆按“环”计变压器按“台”计合并到系统指标后 SAIFI 明显失真但逐项看又挑不出错。原因λ 必须随元件粒度配套混用是最隐蔽的错误特别是线路长度参与折算时漏乘长度会让主干贡献被严重低估。解决建台账时强制统一口径——线路一律折算到“次/(km·年)”长度单独一列变压器按台开关按台但只计自身故障。跑数之前写一段校验把每条支路的 λ 乘上长度得到“支路年故障率”再看全网的加权平均故障率是否落在合理区间内。5.4 联络开关闭合后最小路出现环导致递归出不来现象DFS 求最小路时程序死循环或者同一个负荷点返回两条“等价”路径后续串联折算不知道取哪条。原因配网存在联络线合环运行或检修转供时简单路径不再唯一最小路法赖以成立的“辐射状唯一路径”前提被打破。解决按正常运行方式建图把联络开关一律置为断开再跑最小路需要校核转供场景时改用手动枚举合闸开关组合的方式并且给 DFS 加路径深度上限别让它在环里空转。我的习惯是拓扑建模阶段就把开关状态作为图的边属性存好断开状态直接不进邻接表从源头避免环路。5.5 保护拒动和重合闸成功率被忽略指标虚高现象算出的 SAIFI 比电网实际统计值低 30% 左右一开始怀疑是故障率数据偏乐观后来发现是保护环节压根没建模。原因断路器存在拒动率故障时会越级跳闸把故障影响面从一段分支扩大到上级母线用户停电范围成倍增加重合闸成功率也会决定“一次故障”到底折算成一次还是两次用户停电。解决在最小路法里把越级影响的概率折算进 λ故障率乘上拒动率作为“扩大影响”的附加项重合闸按成功率把永久性故障次数乘 (1-成功率) 累加成功率高的线路段指标改善非常明显。6. 进阶打法用最小路法给蒙特卡洛降维留一手手算校核6.1 手算小网作为常驻验算基准不管项目规模多大我手里永远留一个三条支路的手算小网当“标准答案”。比如一条 S—L1—N1—L2—N2—L3—LP 的链式馈线三条支路的 λ 分别为 0.2、0.15、0.1 次/年r 分别为 4、3、2 小时/次。手算结果λ_LP 0.45 次/年U_LP 0.2×4 0.15×3 0.1×2 1.45 小时/年可用度 A 1 - 1.45/8760 0.999834。每次上新代码先拿这个小网跑一遍非序贯蒙特卡洛的 LOLP 必须落在 1.45/8760 附近对不上就先查参数单位再查拓扑建模最后才查抽样逻辑。6.2 混合法主干路径解析、分支元件抽时序纯最小路法处理不了时序负荷纯蒙特卡洛在支路多时又跑得慢常见的折中是混合法用最小路法把馈线拆成主干路径和分支集合主干上元件少、故障影响规律直接用解析式算分支元件数量大、修复时间离散度高交给序贯蒙特卡洛抽样。主干贡献的 U 是确定性数值分支贡献是随机估计量两者相加就是负荷点的年停电小时数。这样抽样维度从全网点数降到分支点数β 收敛到 0.05 需要的样本量能小一个数量级。实现上只要把前面path_reliability和sequential_timeline两个函数接到同一个负荷点聚合逻辑里主干走解析、分支走时序输出端统一按用户数加权。我现在接到配电网可靠性评估的活儿第一件事不是开写脚本而是把馈线拓扑画出来手算一条三节点路径摆在那当“标准答案”。不管后面用最小路还是蒙特卡洛只要算出来跟这个对不上先查参数单位再查建模最后才碰抽样逻辑——这个顺序帮我挡掉过不少返工。希望帮到你。本文还有配套的精品资源点击获取