简介面向电力系统可靠性研究的一份MATLAB源码资源基于IEEE RTS-24节点可靠性测试系统实现适合电气工程专业学生、科研人员以及进行电网可靠性评估的工程师使用。资源将24个节点的负荷、发电单元等参数和拓扑封装为可直接运行的m文件可快速开展潮流分布、故障恢复、负荷平衡等仿真并支持故障模式与效应分析、可靠性指标计算等扩展研究。整个rar压缩包内共1个m文件包体大小仅2KB轻量紧凑便于阅读和二次修改。目前已有1537人学习/下载借助该文件读者可免去手工搭建RTS-24系统数据的环节直接获得标准测试系统的MATLAB实现进而围绕可用度、故障率、平均修复时间等指标进行实验和结果可视化适用于教学演示、课题预研与论文复现。1. case24_ieee_rts 是什么为什么做了四十年还在用做电力系统可靠性评估的人几乎都会在上手阶段遇到case24_ieee_rts这个名字。这套 24 节点的 IEEE RTSReliability Test System数据最早是 1979 年由 IEEE 可靠性测试工作组公开发布的目的是给学术界和工程界一个统一口径的算例让大家用同一组母线、机组和输电支路去验证可靠性算法。今天你看到的case24_ieee_rts文件通常就是把这份经典数据搬进了 MATPOWER 或自研平台的结构化数据文件包含 24 个母线、32 台发电机组和 38 条支路总装机约 3405MW峰值负荷 2850MW。我刚接触它时最反直觉的体会是文件很小但第一次算出的 LOLE失负荷期望值往往和文献对不上。问题多半不是算法写错了而是负荷时序的粒度、强迫停运率 FOR 的单位、138kV 与 230kV 两电压等级的拓扑处理出了偏差。这篇博客就把这套数据的读取、可靠性指标计算、最小可跑通的评估流程和几个验证技巧一次讲清。适合正在做电网可靠性课题的研究生、调度侧做充裕度评估的工程师以及需要拿标准算例验证新算法的开发人员。2. 拆开 case24_ieee_rts 的文件结构bus / gen / branch 三层2.1 解析 mpc.bus 数据块母线表里有电压等级和节点类型case24_ieee_rts最常见的载体是 MATPOWER 风格的 .m 文件结构是mpc.bus、mpc.gen、mpc.branch三个矩阵。先从 bus 块读起。下面这段 Python 用正则把数据块抠出来按空白切分成数值行import re import numpy as np def load_block(text, name): # 匹配 mpc.name [ ... ]; 之间的全部内容 m re.search(rfmpc\.{name}\s*\s*\[(.*?)\];, text, re.S) if not m: raise ValueError(fblock {name} 未找到) rows [] for line in m.group(1).splitlines(): line line.split(%)[0].strip() # 去掉行尾注释 if line: rows.append([float(x) for x in line.split()]) return np.array(rows) raw open(case24_ieee_rts.m, encodingutf-8).read() bus load_block(raw, bus) gen load_block(raw, gen) branch load_block(raw, branch) print(bus.shape, gen.shape, branch.shape)这段代码的关键在load_block的正则表达式中使用了非贪婪匹配并启用了re.S可以跨行提取整个数组.split(%)是为了去掉 MATPOWER 注释符号后面的内容避免把列说明字符串转成 float 时报错。bus 矩阵每行 13 列标准列顺序是母线编号、节点类型、有功负荷 PD、无功负荷 QD、并联电导、并联电纳、区域编号、电压幅值、电压相角、基准电压 BASE_KV、分区编号、电压上界、电压下界。真正需要盯住的是第 1 列母线编号、第 2 列节点类型和第 3 列有功负荷。RTS-24 中 1~10 号母线属于 138kV 层11~24 号属于 230kV 层两层之间靠若干变压器支路连接。如果你拿到的case24_ieee_rts是别人转换过的格式母线编号和电压等级的对应关系一定要先核对否则后续潮流和削负荷计算会在变压器支路上出现莫名其妙的功率越限。2.2 gen 数据块FOR 并不在 MATPOWER 的原始列里gen 矩阵每行的列数在 MATPOWER 4.x 之后通常是 21 列前几列含义是接入母线号、有功出力 PG、无功出力 QG、无功上限、无功下限、电压幅值设定值、基准容量、开机状态、最大有功 PMAX、最小有功 PMIN。做充裕度计算时PMAX第 9 列就是这台机组能贡献的容量上限。这里有个非常容易踩的坑MATPOWER 自带的case24_ieee_rts潮流案例里没有可靠性计算必须的 FOR强迫停运率和 MTTR平均修复时间。这两个参数存放在 IEEE RTS 原始数据表里需要单独维护成一张表再按机组顺序与 gen 行一一对应。下面的示例是常见的对应方式实际数值务必以你手里的 RTS 可靠性数据表为准机组类型单机容量 / MWFOR 参考量级MTTR / h水电12 / 500.01 左右约 48燃油20 / 1000.10 左右约 50燃煤76 / 1550.02 ~ 0.04约 40燃煤3500.08 左右约 60核电4000.12 左右约 150注意表中 FOR 是长期统计的无故障可用概率的补集不是计划检修停运率。组合这两类停运时要用独立事件的方式计算综合不可用率即FOR_total FOR_forced FOR_maintenance - FOR_forced * FOR_maintenance很多工程失误就出在直接把两个数字相加。2.3 branch 数据块38 条支路里混着线路和变压器branch 矩阵在 MATPOWER 下每行 13 列关键索引是起点母线、终点母线、电阻、电抗、对地电纳、长期载流量 RATE_A、短期载流量 RATE_B、紧急载流量 RATE_C、变比、移相角、投运状态、最小相角差、最大相角差。case24_ieee_rts的 38 条支路由输电线路和变压器支路构成变压器集中在两个电压层的联络位置。计算可靠性时RATE_A 的单位是 MVA而母线负荷 PD 的单位是 MW这在交流潮流里没问题但在做直流潮流削负荷时容易混用。更隐蔽的问题是变压器支路的电抗标幺值基于不同的基准容量case24_ieee_rts的标准 baseMVA 是 100MVA如果你用别的基准重新折算过就要同步修改 RATE_A否则 N-1 扫描时会把变压器容量限制算错。3. 可靠性指标怎么算从状态枚举到蒙特卡洛3.1 状态枚举法在 24 节点上做 N-1 与 N-2 扫描解析法的思路是把系统状态看成一堆元件的组合每个元件只有正常和故障两种状态整个系统的状态概率是各元件概率的乘积。对case24_ieee_rts做状态枚举时先从小到大枚举预想故障集N-1 对应 38 条支路各自退出N-2 对应任意两回线路同时退出。每枚举一个状态就做一次直流潮流检查是否有支路越限或母线失电有则计算需要削减的负荷量。枚举法的概率计算很直接。设第 i 个元件的不可用率为 q_i可用率为 p_i 1 - q_i那么一个第 1、3 条线故障其余正常的系统状态概率就是P(s) q_1 * p_2 * q_3 * p_4 * ... * p_n把所有需要削减负荷的状态概率累加起来就得到失负荷概率 LOLP再乘以对应的小时数就得到 LOLE。这种方法的优点是结果确定、可解释性强缺点是 N-2 及以上的组合数会迅速膨胀32 台机组和 38 条支路组合下来全枚举不现实所以工程上通常只枚举到 N-2高阶风险交给蒙特卡洛补。3.2 非序贯蒙特卡洛用随机抽样逼近系统风险蒙特卡洛思路简单得多每一轮抽样中为每个元件生成一个 [0,1) 均匀随机数如果小于该元件的 FOR就认为它处于故障状态然后判断这轮系统是否切负荷。大量抽样后切负荷的频次和切负荷量的均值就逼近真实风险。下面这段 Python 展示了只考虑发电容量、不建输电网的最简实现适合先验证算法链路import numpy as np rng np.random.default_rng(42) # 固定随机种子便于复现 N 20000 # 抽样轮数每轮代表 1 小时 peak_load 2850.0 # RTS-24 系统峰值负荷单位 MW # 从 gen 表取 32 台机组的 PMAX cap gen[:, 8] # 第 9 列是 PMAX # 外接 FOR 数组长度必须等于机组台数 for_arr np.loadtxt(rts_for.csv) # 每行一台机组的强迫停运率 def sample_loss(): una rng.random(len(cap)) for_arr # True 表示故障 avail_cap cap[~una].sum() # 可用机组容量之和 return max(0.0, peak_load - avail_cap) # 缺额单位 MW loss_samples np.array([sample_loss() for _ in range(N)]) lolp (loss_samples 0).mean() # 失负荷概率 lole lolp * 8760 # 换算为年失负荷小时 eens loss_samples.mean() * 8760 # 期望缺供电量MWh/year print(fLOLP {lolp:.5f}) print(fLOLE {lole:.2f} h/year) print(fEENS {eens:.1f} MWh/year)这段代码的抽样逻辑是rng.random(len(cap))生成一组与机组台数等长的随机数小于 FOR 的机组判为故障cap[~una].sum()把所有未故障机组的容量相加。注意抽样单元是小时所以最终用 8760 把每小时的期望值放大到全年。这里的peak_load固定为 2850MW只反映峰值工况更严谨的做法是把 8760 个时序负荷逐小时代入我一般先用峰值做一次冒烟测试确认代码能跑通再替换成时序负荷。3.3 指标换算LOLP、LOLE、EENS 别搞混很多初学者把 LOLP 和 LOLE 当成一回事其实它们只是同一个风险在不同时间尺度上的投影。下面这张表可以当作换算基准指标全称含义典型单位LOLPLoss of Load Probability某一小时内系统无法满足负荷的概率无量纲LOLELoss of Load Expectation一年内预期失负荷的小时数h/yearEENSExpected Energy Not Supplied一年内预期损失的电量MWh/year三者关系可以写成LOLE LOLP × 8760但这是以小时为抽样单元为前提如果你按天抽样换算系数就是 365。EENS 则进一步考虑切负荷的深度两个系统 LOLE 相同不代表风险相同一个每次只切 50MW、一个每次切 500MWEENS 会差一个数量级。所以工程报告里我通常会同时给 LOLE 和 EENS单给一个容易被质疑。4. 在本地跑通完整可靠性流程从数据文件到风险指标4.1 先做发电充裕度把 32 台机组全口径采进来完整流程的第一步是把 RTS-24 的机组容量、FOR、负荷时序组合成容量 - 负荷比较模型。前面的最小代码只用了峰值负荷这里补上 8760 小时负荷曲线。RTS 原始数据提供按小时、按日、按周区分的负荷系数转换逻辑是对 8760 小时逐小时采样比较当前可用装机与当前负荷判断是否缺电。load_curve np.loadtxt(rts_load_8760.csv) # 长度 8760单位 MW sample_cnt 20000 hour_idx rng.integers(0, 8760, sample_cnt) # 随机抽取小时 loss_list [] for idx in hour_idx: una rng.random(len(cap)) for_arr avail_cap cap[~una].sum() loss_list.append(max(0.0, load_curve[idx] - avail_cap)) loss_arr np.array(loss_list) lo le (loss_arr 0).mean() * 8760 eens loss_arr.mean() * 8760 print(fLOLE {lole:.2f} h/year, EENS {eens:.1f} MWh/year)代码中点在于rng.integers(0, 8760, sample_cnt)先随机选定小时再对该小时做元件抽样这样负荷概率分布和故障概率分布是分开的不会双重随机叠加出偏高的结果。load_curve[idx]取到的是那一小时的实际负荷代入到判断式中。4.2 计入输电网约束为什么必须有直流潮流到此为止只验证了发电容量充裕度case24_ieee_rts的价值更多体现在输电网某条支路退出后即使全网发电容量充足也可能因为断面受限而必须切负荷。判断是否切、切多少的常见做法是跑一次直流潮流再求解一个最小切负荷的线性规划。直流潮流的线性方程是P B * theta加上支路容量约束后切负荷最小化可以写成标准 LP 问题。用 MATLAB 和 MATPOWER 跑某个故障状态下的最小切负荷常见做法是修改mpc.branch中目标支路的BR_STATUS为 0然后调用mpc loadcase(case24_ieee_rts); mpc.branch(k, 11) 0; % BR_STATUS 置 0表示该支路退出 result rundcopf(mpc); % 直流最优潮流 shortfall mpc.total_load - sum(result.bus(:, 3) - result.bus(:, 2));BR_STATUS列是支路的第 11 列置 0 代表停运。rundcopf是 MATPOWER 的直流最优潮流函数它会尝试在满足发电和支路约束的前提下最小化发电成本若系统必须切负荷缺额可以从result.bus中负荷与注入的差看出。注意这里没有启用自动切负荷所以结果只能判断该状态是否安全不能直接得到切负荷量。4.3 MATPOWER 连跑把状态枚举交给脚本循环枚举 N-1 状态时不需要手动改文件用 MATLAB 循环就能批量生成故障状态并记录结果N_line size(mpc.branch, 1); unsafe_states []; for k 1:N_line mpc2 mpc; mpc2.branch(k, 11) 0; % 第 k 回线停运 r rundcpf(mpc2); % 直流潮流不优化 if any(r.branch(:, 6) r.branch(:, 5)) % 实际功率限值 unsafe_states [unsafe_states; k]; end end disp(unsafe_states);这段代码的判断逻辑是r.branch第 5 列是 RATE_A第 6 列是潮流结果对应的功率二者相除超过 1 就是过载。实际可靠性评估不会只做 N-1还会叠加发电机故障所以要把支路停运 机组停运的组合编成双层循环。数据量大时MATLAB 循环性能会成问题我通常会把所有状态改成并行 for 循环parfor或者干脆在 Python 里用mpi4py做分布式扫描。5. 结果对不上文献值时先查这四个细节5.1 用 LOLE 的量级快速定位模型偏差case24_ieee_rts的好处是公开成果多你可以先用文献里的量级验证模型。常见结果是在原始负荷模型下系统的 LOLE 大约在 9 h/year 上下具体数值随负荷曲线处理方式波动。如果你算出来是几十小时甚至上百小时先检查负荷时序是否做了逐小时缩放如果算出来接近 0再检查 FOR 数组是否因为行顺序错位而整体偏小。定位方法很简单把 FOR 全设为 0如果这时 LOLE 还不为 0说明负荷和装机数据本身就不匹配问题在数据不在抽样逻辑。5.2 固定随机种子让每一次改动都可比较蒙特卡洛抽样有方差哪怕算法完全正确两次运行也会差出 3%~5% 的 LOLE。做方案对比时若不固定随机种子你根本分不清改善是来自方案还是抽样噪声。我建议所有采样入口统一使用同一个numpy.random.default_rng(seed)并且把 seed 作为命令行参数传入这样每次跑批只需换 seed代码和结果都能复现。对比两个扩容方案时更讲究一点的做法是用公共随机数两个方案共用同一批随机序列只改变容量数组这样两轮结果的差值主要来自方案差异抽样方差被大幅抵消。实现上就是先一次性生成随机矩阵再分别套用不同的cap数组而不是在循环里各自调用随机数生成器。5.3 别忘了台账里的隐性共用停运RTS-24 最坑人的一点是部分机组存在共用水电头或成组停运关系不能把所有机组当作独立元件。例如同一条水道上的多台水电机组一台检修时其余机组也要降出力成对配置的机组也存在关联停运模式。若你的数据文件里只给了单机 FOR直接采样会把系统可靠性算得偏高。验证方法是做一次全停运事件交叉核对把 N 轮抽样中出现两台及以上相同类型机组同时故障的事件单独打印出来看与实际检修计划是否矛盾。发现关联停运时不能简单把 FOR 相加而应该在抽样时先抽机组群状态再抽单体状态形成两级抽样结构。这一点对 5 年以上经验的人尤其重要因为很多人早期用 case24 时都被它骗过。5.4 用 N-1 结果做二次交叉验证算完蒙特卡洛后再独立跑一遍解析法的 N-1 预想事故分析比较两者给出的高风险支路排序是否一致。如果蒙特卡洛认为某条支路风险最高而 N-1 扫描里它从未引起过载就要检查抽样场景里是否混入了不可能出现的组合。把两个方法的高风险支路输出成两张表逐条核对重叠度是排查数据错位的最快手段。重叠度足够高之后再用 EENS 的分母线结果去和公开的 RTS 可靠性报告比对确认风险在空间分布上也对得上。本文还有配套的精品资源点击获取