简介：面向电力工程学习者与从业者，这套MATLAB实现案例聚焦电力系统可靠性评估，覆盖可靠性建模、故障概率计算、负荷分配、状态估计以及SAIFI/SAIDI/ASAI等常用指标计算，可服务于课程设计、科研验证或工程初期分析。压缩包共20个文件，以17个m源码为主，辅以2个asv备份和1个readme说明，整体仅19KB，结构轻量、便于直接加载运行。目前已有3552人学习下载，具备一定参考热度。代码示例涉及元件故障率、修复时间、备用容量等参数建模，并借助MATLAB矩阵运算、蒙特卡洛仿真和潮流计算完成系统在不同故障场景下的行为分析；readme还说明了执行步骤与常见问题，有助于使用者快速上手。对于想掌握可靠性计算思路的读者，既可借此理解贝叶斯公式、寿命分布等概率方法在电力系统中的应用，也能通过修改参数扩展为适配自身研究数据的评估工具。
1. 电力系统可靠性评估的 MATLAB 实现：这套资源解决什么问题
做电网规划或者运行方式分析的人，早晚要面对一个问题：新增一台机组、改一条联络线，系统可靠性到底变了多少？ 这没法靠拍脑袋，得算。网上能查到的可靠性评估资料大多是理论推导，真正能直接跑的 MATLAB 代码反而难找——要么太老，要么只覆盖单一算法。这份“电力系统可靠性评估(matlab实现)”的资源，就是把可靠性评估从“公式能看懂”推到“仿真能出数”的一步：它同时提供了解析法和蒙特卡洛法的 MATLAB 实现，输入发电容量、负荷曲线和元件可靠性参数，就能算出 LOLP、EENS 这类指标。适合正在做电力系统课程设计、毕业论文，或者刚接手规划项目但还没搭过可靠性仿真代码的从业者，你拿到手之后要做的不是重写算法，而是改参数、换数据、读懂每一段计算的边界。
2. 评估模型选型：为什么先定评估算法再写代码
2.1 三种常用算法的适用边界
可靠性评估的算法选型，直接决定你后面写的每行代码。常见做法是把算法分成三类：状态枚举法、非序贯蒙特卡洛法、序贯蒙特卡洛法。这三者不是随便选一个都能用，它们的计算量和能回答的问题差别很大。
| 算法 | 基本原理 | 输出指标 | 优点 | 缺点 |
|---|---|---|---|---|
| 状态枚举法 | 枚举所有元件状态组合，按概率加权 | LOLP、EENS | 理论严谨，结果稳定，无随机误差 | 元件多时组合爆炸，只能截断枚举 |
| 非序贯蒙特卡洛 | 随机抽样生成系统状态，统计切负荷 | LOLP、EENS | 不受元件数量限制，实现简单 | 不能处理时序相关场景，收敛慢 |
| 序贯蒙特卡洛 | 按时间轴模拟元件运行/故障/修复 | LOLP、EENS、LOLE、频率指标 | 能模拟机组检修、负荷时序变化 | 需要时序负荷数据，计算量大 |
选型时先问自己三个问题。第一，系统规模多大？ 节点数少于 30、元件数少的场景，状态枚举法收敛快且结果可复现，适合做教学演示和方案对比。第二，你有没有时序负荷数据？ 如果只有年最大负荷和装机容量，没有 8760 小时的负荷曲线，序贯法根本没数据可喂。第三，你要不要算频率和持续时间类的指标？ 非序贯法算不了期望缺供频率（EFC），序贯法才有这个能力。
这份资源里我印象较深的是它把三种方法都放了进去，但非序贯蒙特卡洛的代码完整度更高，状态枚举给了截断部分，序贯法的数据文件里带了时序负荷样例。你上手时优先跑非序贯模块，出数之后再去看另外两种，这样不容易被算力问题劝退。
2.2 非序贯蒙特卡洛法的实现逻辑拆解
非序贯蒙特卡洛法的核心思路只有一句话：把多状态系统拆成独立抽样事件。接线图里 30 条线路、20 台发电机，每个元件都有“正常运行”和“故障停运”两个状态（发电机可能有三状态，即全额运行、降额运行、停运，这里先简化处理）。
每一次抽样做三件事：用均匀分布随机数发生器的输出决定每个元件状态；把处于故障状态的元件从系统中移除；对剩下的系统做潮流计算，判断有没有切负荷。重复抽样几千次后，统计切负荷次数占比，就是 LOLP 的估计值。
这里最容易误解的是“有序抽样”和“无序抽样”的差别。非序贯法不关心时间顺序，只关心特定时刻系统的静态状态。某线路今天坏还是明天坏，对它来说没有区别——它只判断“此刻有多少元件是坏的”。因此它的随机数发生器只需要产生 [0,1] 均匀分布，然后用阈值映射状态：
% 判断某发电机是否故障，故障概率用 FOR (forced outage rate) 表示 % u 是 rand 函数产生的一个随机数 for i = 1:ngen u = rand; % 产生 [0,1] 均匀分布随机数 if u < gen_FOR(i) % 小于强迫停运率则判定为故障 gen_status(i) = 0; % 0 表示故障 else gen_status(i) = 1; % 1 表示运行 end end这段代码的逻辑看起来简单，但有一个关键点：抽样次数。 每次抽样都是一个独立试验，LOLP 是这些试验结果的均值。抽样太少，结果的置信区间宽得没法用；抽样太多，计算时间又受不了一万次直流潮流的求解。一般我会在代码里加一个统计计数器，每次抽样后检查方差系数是否小于阈值，达到阈值就提前终止。这个收敛判据在后面的章节单独展开。
2.3 元件可靠性参数怎么设才合理
元件参数是可靠性仿真结果的上限决定因素。参数错一个数量级，算出来的 LOLP 就完全是另一套东西。你需要准备的参数表包含这几项：发电机的强迫停运率 FOR（单位是百分比，比如 0.02 表示 2% 时间处于停运），平均修复时间 MTTR（小时），平均无故障运行时间 MTTF（小时），线路的故障率（次/年），以及线路修复时间。
这里有一个参数之间的换算关系，也是初学者最容易翻车的地方。FOR 和 MTTF/MTTR 的关系是 FOR = MTTR / (MTTF + MTTR)。如果你知道某台机组一年平均故障 0.5 次，每次修 10 小时，那么它的可用率是 1 - 10/(8760/0.5) ≈ 0.9994，FOR 约 0.00057。直接用这个 FOR 填到代码里，比从记录的统计数据里硬算更快。
我一般会建议把参数集中放在一个 data.xlsx 里，不要让参数散落在代码各处。 因为后期调整参数跑敏感性分析时，在表格里改一个数字比重编代码要快得多，也减少出错概率。资源里包含的数据文件基本就是这个思路，发电机、线路、负荷三张工作表，改参数后重新跑主脚本即可。
3. 数据准备与核心算法落地：从负荷曲线到切负荷计算
3.1 数据文件的组织方式
把电力系统可靠性评估跑通，第一步不是写潮流函数，而是把数据文件组织好。多数新手拿到资源直接运行主程序，发现报错找不到变量，原因不是代码有 bug，而是数据没按预定的格式放置。这份资源的目录下，数据和代码分开放：一个data文件夹放输入数据，一个src文件夹放可执行脚本，results文件夹放指标输出。
输入数据至少包含三个数组。第一个是发电机数据矩阵gen_data，每行一台机组，列分别是节点编号、额定容量、FOR 值、最小出力、最大出力、发电机类型（火电/水电/可再生能源）。第二个是线路数据矩阵line_data，每行一条线路，列是首端节点、末端节点、电抗值、传输容量上限、故障率、修复时间。第三个是负荷数据矩阵，可以是一个 24 点的日负荷曲线，也可以是 8760 点的年时序曲线，取决于你用哪种算法。
除了这三个核心数据，还有一个容易被忽略的参数向量：系统基准容量。 潮流计算全部采用标幺值，发电机容量和线路容量如果不用同一基准值折算，算出来结果会出现偏差。
常见做法是用一个load_case.m脚本统一读取并处理这些数据，而不是让主程序直接操作原始表格。好处是后续换数据源时，只改load_case.m里的文件路径和行列映射，主程序保持不变。
3.2 切负荷判断与直流潮流求解
非序贯蒙特卡洛每次抽样后，系统的拓扑可能和上一次完全不同。判断“哪个状态会导致切负荷、切多少负荷”，是可靠评估里计算量最大的部分。
工程中最常用的是直流潮流模型：忽略有功损耗和无功分布，只保留 P-θ 线性关系。交流潮流更精确，但迭代计算量大，在蒙特卡洛几千次循环里跑交流潮流，计算时间成倍增加。可靠性评估里 95% 的场景用直流潮流就够了——你评估的是“够不够用”，不是“稳不稳定”。资源里的实现默认用直流潮流，切负荷计算则按“等比例切负荷”方式处理：
% 求解直流潮流：B 为节点电纳矩阵，P 为节点注入功率向量 % theta 为节点相角向量，Pbalance 是每个节点切掉的负荷 theta = B_inv * P_inject; % 节点相角 line_flow = B_f * theta; % 线路有功潮流 % 找到过载线路 overload_idx = find(abs(line_flow) > line_limit); total_shed = 0; while ~isempty(overload_idx) % 等比例切负荷，每次切掉过载线路对应的小量负荷 for k = 1:length(overload_idx) line_id = overload_idx(k); shed_amount = abs(line_flow(line_id)) - line_limit(line_id); total_shed = total_shed + shed_amount; % 将切负荷量分摊到线路两端节点 i_node = line_from(line_id); j_node = line_to(line_id); gen_output(i_node) = gen_output(i_node) + shed_amount/2; gen_output(j_node) = gen_output(j_node) + shed_amount/2; end % 重新求解潮流 theta = B_inv * P_inject; line_flow = B_f * theta; overload_idx = find(abs(line_flow) > line_limit); end这段实现有几个细节需要解释。 第一，切负荷不是一次性切完，而是循环迭代实现——先找到过载最严重的线路，切掉对应分量，再重新算潮流，看是否还有其它线路过载。第二，切负荷量分摊到两端节点，是一种近似处理，实际调度可能更偏向负荷侧切，但如果只做可靠性指标统计，这种近似是业内默认可接受的。第三，循环必须有终止条件，不能无限制切下去，实际代码里会设一个最大迭代次数，超过即终止并报错。
参数说明：B_inv是电纳矩阵的逆矩阵，在每次循环中因为移除了故障线路，这个矩阵会变化，不能只算一次；line_limit是线路容量向量；shed_amount是最小切负荷量，遇到过负荷幅度很大的情况，单次切负荷量不能超过该节点总负荷，否则会出现负负荷的不合理结果。
3.3 抽样循环怎么写才高效
核心循环结构决定了整个仿真运行时间，非序贯蒙特卡洛的基本循环体大概要写成这样：
% 非序贯蒙特卡洛主循环 n_samples_max = 50000; % 最大抽样次数 variance_coef = 0.03; % 方差系数收敛阈值 lolp_accum = 0; % 累计切负荷状态次数 eens_accum = 0; % 累计电量不足期望值 for k = 1:n_samples_max system_state = sample_state(); % 随机抽样各元件状态 [shed, status] = evaluate_state(system_state); % 潮流计算与切负荷判断 if status == 1 % 该状态需要切负荷 lolp_accum = lolp_accum + 1; eens_accum = eens_accum + shed; end % 每50次抽样判断一次收敛 if mod(k, 50) == 0 && k > 500 p_lolp = lolp_accum / k; var_coef = sqrt((1 - p_lolp) / (k * p_lolp)); if var_coef < variance_coef break; end end end LOLP = lolp_accum / k; EENS = eens_accum / k * 8760;这个循环在 MATLAB 里要让“状态抽样”和“潮流计算”尽量函数化，避免在循环里重复读表、重复分配大数组。sample_state()和evaluate_state()做成两个独立函数文件，主循环只做计数和判断。 这么做主要为了两点：一是代码可读性好，二是 MATLAB 的 JIT 编译对结构清晰的代码优化更有效。
variance_coef这个收敛阈值，我常用的经验值是 3%——过高结果波动大，过低计算量爆炸。50000 次是上限保护，一般系统在 3000-10000 次之间就会满足收敛。判断收敛不是只看 LOLP 不再变化，用方差系数是关键：它衡量的是估计值的相对波动，比“连续 1000 次抽样 LOLP 都不变”这种粗糙判断可靠得多。
4. 可靠性指标的计算与解读：LOLP、EENS 不是越大越好
4.1 指标背后的物理含义
算完蒙特卡洛循环后，会得到几个指标，但很多人会忽略它们的经济含义。LOLP（Loss of Load Probability，失负荷概率）衡量的是“在抽样时刻系统无法满足负荷需求的概率”，它是有量纲的概率值，通常以天/年为单位使用。而 LOLE（Loss of Load Expectation）是 LOLP 乘以时间区间得到的期望值，比如年负荷曲线数据下 LOLE 的单位是小时/年。
EENS（期望缺供电量）则直接给出经济损失的计算基础：它表示在评估周期内因供电不足导致用户少用的电量期望值，单位是 MWh/年。规划部门通常会把 EENS 乘以单位电量停电损失费用，估算可靠性改善带来的经济效益。
另一个需要关注的指标是切负荷频率。 如果资源里没有直接给出，你可以用状态枚举法辅助统计，或从序贯法的状态转移记录里提取。注意非序贯法的结果里不应直接给出频率指标，因为它的抽样序列没有时间关联性，切负荷频率是统计意义上的“每个时段”，不是真实的物理事件频率。看到资源里某个函数名直接输出频率指标时，要确认它用的是序贯法还是非序贯法——混用的话结果没有物理意义。
4.2 指标统计的 MATLAB 实现细节
主循环跑完，指标计算本身只有短短几行，但有一个细节影响准确性——EENS 的单位转换。如果不注意每个抽样代表的是“一小时”还是“一个时段”，EENS 的值会差 8760 倍或 24 倍。
% 指标统计输出 LOLP = lolp_accum / k; LOLE = LOLP * 8760; % 单位：小时/年 EENS = eens_accum / k * 8760; % 单位：MWh/年 expected_energy = sum(load_profile); % 系统年总电量需求 % 结果写入结构体，方便后续对比方案 results.LOLP = LOLP; results.LOLE = LOLE; results.EENS = EENS; results.n_samples_used = k; % 实际有效抽样次数 fprintf('LOLP = %.4f\n', LOLP); fprintf('LOLE = %.2f 小时/年\n', LOLE); fprintf('EENS = %.2f MWh/年\n', EENS);这里要注意LOLP是无量纲概率，LOLE是小时数/年，乘的系数取决于你的负荷数据时间颗粒度——如果抽样用的是 8760 点的小时负荷，乘 8760；如果是 365 点日峰值负荷，乘 365 得到的是 “天/年” 的 LOLE 口径。很多论文里的 LOLE 值，如果不注明口径，横向对比时容易产生误读。
4.3 设计人员怎么解读指标数值
算出来的 LOLP=0.001 意味着什么？ 这代表一年 8760 小时里，预计约有 8.76 小时存在失负荷风险，LOLE 约 8.76 小时/年。国内规划常用的可靠性准则大致在“LOLE 小于 24 小时/年”的范围内，低于这个水平可以认为系统供电充裕。
对比方案时，不要只看 LOLP 的绝对值，还要比较相对变化。 增建一台 200MW 机组，LOLP 从 1.2 降到 0.8，绝对变化只有 0.4，但切负荷期望的大幅下降会体现在 EENS 上。EENS 的数值通常波动范围比 LOLP 大，在抽样次数不够多或负荷接近临界水平时尤其敏感。因此我在评估方案时，习惯同时输出两张表：一张是指标绝对值表，一张是相对变化率表。
另外一个容易误用的场景是把指标直接当调度工具。 可靠性指标描述的是长期统计平均特性，异常天气、特殊检修方式等短时场景不在它的适用范围内。如果你要评估“某条线路检修期间系统会不会出问题”，应该用静态安全分析逐时刻校验，而不是重新跑一遍可靠性评估。
5. 常见问题与避坑：从随机数种子到收敛判据的六个排查点
5.1 两次运行结果不一致，无法复现
现象：同一套数据和参数，连续跑两次得到的结果有 5%-10% 的差异，EENS 甚至差 20% 以上。
原因：非序贯蒙特卡洛的随机数发生器每次取的种子不同，抽样序列完全不同。如果抽样次数没有收敛到很低的方差系数，结果波动在所难免。
解决：在脚本开头加一行rng(2024);固定随机数种子，保证结果可复现。 如果固定种子后结果依然波动很大，说明收敛阈值 0.03 太宽松，改为 0.01 或 0.005，同时把最大抽样次数调大。这属于正常代价：可靠性评估要的是稳定比较，不是单纯速度。
5.2 收敛判据永远达不到，循环跑到上限
现象：mod(k, 50) == 0的收敛判断从未触发，程序一直跑到n_samples_max才跳出。
原因：当 LOLP 极低时（比如小于 0.0005），方差系数公式里的分母k * p_lolp很小，需要极大的抽样次数才能收敛。这是一个数学必然，不是代码 bug。
解决：给收敛条件加一个上限保护，同时接受“达到最大抽样次数”也算正常终止。 更合理的方法是同时让方差系数达到 0.05 且 LOLP 值保持稳定时即跳出，而不是只依赖方差系数一个指标。如果你知道该系统的 LOLP 大概在很低量级，可以提前把目标方差系数放宽到 0.05，节约计算时间。
5.3 潮流计算矩阵奇异，直接报错
现象：切负荷计算循环里B_inv矩阵求逆时报 “singular matrix” 警告或直接报错。
原因：当一个区域因线路故障被完全隔离开时，形成了孤岛，该孤岛的子矩阵不可逆。孤岛里电源和负荷都不平衡，直流潮流方程无解——出现这种情况其实是评估里该算“切负荷”的状态，不应该当作崩溃处理。
解决：在做矩阵求逆前判断系统连通性，断开孤岛内所有负荷。 实现上可先对节点做连通性判定：如果有孤岛，孤岛内的全部负荷视为切除，不参与潮流计算。这需要写一个基于图的连通性检测函数，MATLAB 里用graph对象加conncomp函数即可实现。资源里的代码如果没做这一步，建议自行补充，否则遇到 N-2 故障场景容易翻车。
5.4 用了交流潮流，结果反而比直流潮流“好”
现象：把直流潮流换成交流潮流后，LOLP 反而变低了，与直觉相反。
原因：交流潮流考虑无功约束后，部分过载线路可能因为电压支撑不足先触发切负荷，但你的交流潮流模型如果没加无功容量限制和电压上下限，相当于只考虑了有功平衡，把交流当直流用，同时还多了迭代误差。
解决：一致性原则。 评估目的如果是长期可靠性，直流潮流足够；如果必须用交流潮流，要把无功容量、电压控制逻辑都建模进去，否则精度不会提升，只会增加变量。新手不要轻易把评估代码里的潮流模型从直流换成交流，那不是简单替换一个函数的事。
5.5 新增机组后 LOLP 变化不大，方案优化无效
现象：增加一台大容量机组（比如 600MW），LOLP 只降了几个百分点，不确定是不是系统瓶颈在别处。
原因：如果原有系统已经较充裕或输电阻塞是主要限制，增加发电容量对可靠性提升的帮助有限。可靠性评估的输出只有指标，没有告诉你短板在哪里。
解决：对比多个“改造方案”时，同时输出线路负载率和发电出力统计文件。 实际做法是在主循环里额外记录每条线路的最大潮流峰值和过载频次，一并写入结果。看到哪些线路频繁过载，你就知道真正的瓶颈是网络约束而不是装机容量。这份资源的主程序没把线路潮流统计输出，但你在evaluate_state函数里加几行矩阵累加，就能把这个维度补上。
5.6 可再生能源出力模型处理不当
现象：把光伏和风电按固定 FOR 值处理，结果得到极高的 EENS，与实际运行不符。
原因：可再生能源不是简单地“故障/运行”两种状态，出力的随机性不能用 FOR 表示。 如果直接把风电按 FOR=0.5 当 50% 概率停运处理，相当于把间歇性等同于强迫停运，评估结果当然差。
解决：要把可再生能源建模为多状态或直接叠加时序输出曲线。 时序法里直接采 8760 小时的实测出力曲线，非序贯法可以用多状态简化：按出力区间划分，每个区间一个概率值。资源里的发电机数据结构是支持这个扩展的，只需把 FOR 字段替换成多状态表。
6. 把评估从离线统计做到方案对比：扩展出你需要的对比能力
拿这套程序做规划方案对比时，最管用的做法不是改算法，而是在外面包一层“方案循环”。 每个方案用同一个主循环函数，只换gen_data或line_data这两份数据表，然后对比指标。我通常会在compare_plans.m里这样组织：
plans = {'base_case', 'plan_A', 'plan_B'}; for p = 1:length(plans) load_case(plans{p}); % 加载对应方案的数据 [LOLP(p), EENS(p), LOLE(p)] = run_mc_main(); end % 输出对比表 table(plans', LOLP', LOLE', EENS', 'VariableNames', ...);这么写的好处是计算代码完全不动，方案配置全在数据文件里。 实际经验里还有一条很实用：给每个方案固定同一个随机种子，这样不同方案之间的指标差异完全来源于系统变化，抽样噪声被压到最小。 固定种子是方案对比里最重要的技巧，没有之一。
验证这套资源算得准不准，我建议做一个三步交叉验证。 第一步，把所有元件 FOR 设为 0，LOLP 必须为 0。 第二步，把负荷设为极大值，LOLP 必须为 1。 这两个边界测试过了，基本说明潮流和切负荷逻辑没大问题。 第三步，用状态枚举法跑一个小系统对比，解析法结果和蒙特卡洛结果之间的偏差应该在方差范围之内。
回顾我自己做过的一个模拟项目X的可靠性评估，当时就是因为随机数种子没固定，方案 A 比方案 B 多算了几个百分点的提升量，差点做出错误结论。 从那以后，我每次跑对比都强制先写一行rng(固定值)，再跑方案循环。 如果你手头的评估场景没有时序数据需求，这套资源里的非序贯模块足够应付大部分规划和课程设计需求；如果你的课题需要频率指标或考虑时序检修策略，扩展方向也已经有代码骨架可以往上搭。希望这份拆解能帮你把可靠性评估从公式变成能出结果的工具。
本文还有配套的精品资源，点击获取
