简介：资源包内含423个文件，整体约8.07MB，主体为262个m格式脚本，辅以mat数据文件、c/h源文件以及zbak备份等，构成一套完整的MATPOWER交流级联故障模型实践环境。模型面向电力系统弹性分析场景，适合电气工程专业学生、电网运行研究人员以及MATLAB电力仿真初学者，用于模拟电网遭受扰动后连锁故障的传播、保护动作与系统恢复过程。通过调整仿真参数、设置不同故障场景，读者可观察电压、频率变化及潮流分布调整，进而评估电网在极端事件下的稳定性与恢复力。压缩包保留了MATPOWER相关源文件、示例脚本、说明文档与工具文件，目录结构清晰，便于按模块调用和二次开发。已有50人学习浏览，对于希望借助开源工具开展电网弹性定量分析的读者，这套文件提供了可直接运行的建模基础与实验参考。
1. 交流级联故障模型与电网弹性分析：为什么静态N-1校验不够用
一次极端天气同时断掉多条输电线路，系统往往不是立刻全黑，而是在几分钟内因为潮流转移、线路过载、保护动作而连锁跳闸。传统N-1静态校验只回答“单条线断开后系统是否越限”，回答不了“多条线断开后会崩成什么样”。基于MATLAB的MATPOWER交流级联故障模型，做的就是这件事：用交流潮流反复迭代“断线→重新计算→切除过载线路”的过程，量化负荷损失规模和系统能扛住多大的初始扰动，直接服务于电网弹性分析。适合电力系统研究生、规划工程师和调度分析人员，用来评估极端场景下的系统脆弱点，而不是做常规可靠性计算。
2. AC级联故障模型的两层设计：潮流内核与级联循环，以及它和DC模型的差距
2.1 级联失效的闭环逻辑：扰动、过载识别、切除策略、终止条件
级联故障模型在软件层面要拆成两层：底层是一次次交流潮流计算，顶层是“切除决策”循环。前者负责给出当前拓扑下各支路的有功、无功和电压，后者负责判断哪些支路已经超出热稳定限额、切除它们会怎样影响下一轮潮流。MATPOWER在这两层里各有对应：runpf()负责潮流，算例里的RATE_A、RATE_B、RATE_C三列是过载判断的参照基准；循环本身则需要自己写，MATPOWER没有把所有场景都替你包好。
一次完整的级联循环按下面的顺序推进。初始扰动一般是手工指定一组断线支路，或者从随机抽样中获得；模型把对应支路的BR_STATUS置0，调用runpf()求交流解。潮流收敛后，把所有支路的潮流绝对值与其RATE_A比较，超过阈值（比如1.0表示满载，1.2表示允许短时过载20%）的支路进入候选切除集合。然后按“过载最严重优先”的顺序切除全部或部分候选支路，进入下一轮潮流计算。终止条件有四类：没有新支路过载，系统达到稳定态；潮流不收敛，判定系统崩溃或已分裂成不可计算的孤岛；达到预设的最大迭代轮数；或者全部负荷已丢失。稳定态的判据要写得比潮流收敛严格得多，因为潮流“收敛了”不代表系统安全，电压低于0.8 p.u.的母线和无功严重越限的发电机组，都说明系统其实已经处于不健康状态。
切除策略对结果影响非常大。每轮把过载支路全部切掉，模拟的是“保护完全正确动作”的乐观情形；每轮只切过载最严重的1到2条，模拟的是保护时序有先后、系统有一定缓冲的中间情形。实际电网里保护动作时间并不一致，所以我一般把“每轮最多切除数”做成可调参数，默认值取3，观察不同切除策略下的最终损失差异。
2.2 为什么用AC潮流而不是DC潮流：模型精度与计算代价的取舍
直流潮流的级联模型在很多论文里出现过，求解快、数值稳定，网络参数变化后不需要重新迭代求解。但直流模型把电压、无功、网损和电压相角之间的耦合全忽略了，这在高负载率场景下会给出偏乐观的结论。典型例子是：某条线路断开后，潮流转移到相邻回路，相邻回路的电压可能被拉低，无功损耗剧增，进而让本不越限的有功潮流因为电压下降而进一步上升——这种“电压-无功-有功”的连锁关系，直流模型完全捕捉不到。
交流模型把这些耦合接进来了，代价是收敛性问题。级联过程中系统不断被推向重载甚至接近电压崩溃点，牛顿法很容易不收敛。这时候要区分“数学上无解”和“算法上没找到解”：后者常常可以通过改初值、换快速分解法pf.alg选项或者给某些母线加少量无功补偿来缓解。工程上我的经验是，先让交流模型跑通case9和case39这类标准算例，再往自己的网络上迁移，不要一上来就上一个几百节点的规划模型，否则你会花大量时间在被潮流不收敛本身干扰的排错上。
AC和DC的对比可以从下面的角度把握：
| 对比维度 | AC潮流（牛顿法/快速分解法） | DC潮流（线性化） |
|---|---|---|
| 变量范围 | 电压幅值、相角、有功、无功 | 仅相角和有功 |
| 对电压崩溃的反映 | 能反映电压失稳过程 | 无法反映 |
| 对无功越限和补偿的敏感度 | 敏感，可作为切除判据 | 不敏感 |
| 数值稳定性 | 重载下易不收敛 | 稳定 |
| 单次求解耗时 | 慢，量级为毫秒到秒 | 很快 |
| 适合场景 | 弹性分析、事故后状态评估 | 大规模静态安全分析早期筛选 |
如果只是做几万次蒙特卡洛的粗筛，用DC先砍掉大量无关场景，再对少数高风险场景跑AC，是性价比最高的组合。但你要清楚，DC筛掉的场景里可能包含“电压崩溃”导致的级联，这些恰恰是弹性分析最关心的尾部风险。
2.3 把交流级联结果翻译成弹性指标
级联仿真跑完后，输出不能只是一张“哪些线路被切了”的清单，需要量化为弹性指标，才能横向比较不同方案。常用指标有三类：
| 指标 | 定义 | 计算方式 |
|---|---|---|
| EDNS | 事件造成的期望缺供电量 | 各级联轮次损失负荷的累计值（MWh） |
| 幸存负荷比例 | 最终仍带电负荷与初始总负荷之比 | 1 - 总丢失负荷 / 初始负荷 |
| 级联轮数 | 从初始扰动到稳定态经历的切除次数 | 循环迭代次数 |
| 韧性三角形面积 | 供给损失在时间轴上的积分，衡量事件影响深度 | trapz(时间序列, 失去的负荷) |
| 关键线路集 | 初始断线导致损失最大的线路组合 | 对每条线路做N-1或N-2扫描后的EDNS排序 |
弹性分析和传统可靠性分析的最大区别在于：它不关注随机故障概率乘以后果的期望值，而关注“极端但低概率”的事件发生后，系统能不能把损失控制在一定范围。所以用这个模型做出来的EDNS不应该被平均化处理，而要看尾部分位数，比如99%分位数下的负荷损失。这也是为什么“只跑一个小规模算例给一个期望值”的用法，会让这个模型的价值大打折扣。
3. 在MATLAB里把AC级联主循环跑起来：从case9到case39
3.1 最小实现：带过载切除的仿真主循环
下面这段代码是完整的AC级联故障主循环，输入一个MATPOWER算例名和初始断线支路集合，输出每一轮被切除的支路、最终负荷损失和系统状态。这是我这套方案的最小骨架，直接复制到MATLAB里，配好MATPOWER路径就能跑。
function out = ac_cascade_sim(caseName, initOutages, opt) % AC级联故障仿真主循环 % caseName: 算例名，如 'case9' 或 'case39' % initOutages: 初始断开的支路行号向量，如 [10 23] % opt.maxIter: 最大迭代轮数, 默认50 % opt.tripLimit: 每轮最多切除系数, 0表示不限制 % opt.overloadTh: 过载阈值, 默认1.0（RATE_A的倍数） mpc0 = loadcase(caseName); mpopt = mpoption('out.all', 0, 'verbose', 0); baseMVA = mpc0.baseMVA; if nargin < 3 || isempty(opt) opt.maxIter = 50; opt.tripLimit = 3; opt.overloadTh = 1.0; end % RATE_A为0时视为容量无穷大，避免除零 rateA = mpc0.branch(:, 6); rateA(rateA <= 0) = Inf; out.tripHistory = {}; % 每轮切除的支路行号 out.lostLoad = []; % 每轮丢失的有功负荷 out.status = 'stable'; % 最终状态: stable/collapse/maxIter baseLoad = sum(mpc0.bus(:, 3)); % 初始总有功负荷 % 初始断线 mpc = mpc0; mpc.branch(initOutages, 11) = 0; for k = 1:opt.maxIter % 交流潮流计算 mpc2 = runpf(mpc, mpopt); if mpc2.success == 0 out.status = 'collapse'; % 潮流不收敛，按全系统损失处理 out.lostLoad(end+1) = baseLoad; break; end % 还原成与输入case一致的支路顺序再取潮流结果 mpc2 = int2ext(mpc2); branch = mpc2.branch; % 过载支路识别：abs(PF) / RATE_A loading = abs(branch(:, 12)) ./ rateA; overloaded = find(loading > opt.overloadTh & branch(:, 11) ~= 0); overloaded = overloaded(:)'; if isempty(overloaded) out.status = 'stable'; break; end % 按过载程度排序，优先切除最严重的 [~, idx] = sort(loading(overloaded), 'descend'); if opt.tripLimit > 0 trip = overloaded(idx(1:min(opt.tripLimit, numel(idx)))); else trip = overloaded(idx); end out.tripHistory{end+1} = trip; % 计算本轮到下轮之间的负荷损失： % 被切除支路直接导致母线孤立时，该母线负荷计入损失 % 这里按“潮流不收敛后统一计损”简化，逐母线计损在第4章展开 mpc.branch(trip, 11) = 0; % 记录当前仍带电的总负荷（视为未损失） % 简化处理：只统计所有发电机出口有功之和作为运行负荷 out.lostLoad(end+1) = max(0, baseLoad - sum(mpc2.gen(:, 2))); end if k == opt.maxIter out.status = 'maxIter'; end out.finalServed = max(0, baseLoad - sum(out.lostLoad)); out.edns = sum(out.lostLoad) * 0.5; % 每轮按半个评估周期折算 end这段代码的逻辑要点是：每轮先用runpf求交流解，再对支路潮流做热稳定判断。过载判断用的是第12列PF与RATE_A的比值，注意不是用电流，因为MATPOWER的交流潮流结果直接给出支路有功，按有功与限额比是工程上最快的近似。int2ext这一步不能省：runpf返回的是MATPOWER内部节点排序后的结果，不还原的话，支路行号和输入算例对不上，后面所有统计都是错的。
参数上最值得调的是tripLimit。设成1模拟“保护按严重程度一台一台动作”，级联轮数会变长，但能观察到多轮转移过程；设成很大比如999，就是一次性切除所有过载支路，过程快但可能跳过中间状态。初次跑建议先设3，看tripHistory的变化再调整。
3.2 关键参数怎么调：电压门槛、迭代上限与RATE_A处理
除了上面的基础参数，实际建模中有三个参数最容易影响结果，值得单独说。
第一是电压考核门槛。MATPOWER的潮流计算不把电压幅值当作约束，它只负责求出方程的解。所以你会看到某些算例里系统“收敛了”，但大量母线电压只剩0.75 p.u.——这已经不是稳态运行点，而是濒临电压崩溃的亚稳态。我一般会在每轮潮流后追加一个电压检查：任何母线电压低于0.85 p.u.（具体取值看你的电压标准，国内一般按0.9 p.u.考虑）时，把该母线的部分恒阻抗负荷切掉，再重新算潮流，模拟低压减载装置动作。这个机制对弹性分析的EDNS结果影响很大，不加它的话，电压崩溃场景会被误报为“系统还能运行”。
第二是迭代上限。级联过程在真实的保护动作中通常持续几十秒到几分钟，映射到仿真循环里一般不超过二十轮。设得太大（比如100）会让程序在重载系统上反复震荡，最后在“过载→切除→潮流不收敛”之间卡死；设得太小（比如5）又可能漏掉缓慢的无功转移过程。建议从30起步，观察你的算例通常在几轮内进入稳定态或崩溃，再收紧。
第三是RATE_A的预处理。标准算例里RATE_A、RATE_B、RATE_C不是全部填了实数，有些是0，代表该方向容量未限制或未填写。如果直接拿0做分母，loading是无穷大，所有连接该支路的线路都会被误切除。所以代码里我先把rateA(rateA <= 0)置为Inf，再参与比较。这个坑在实际工程数据里极其常见，因为很多电网导出数据里热稳限额没有填全。
3.3 批量N-k扫描和蒙特卡洛抽样
单次仿真只能回答“这一个初始断线组合会怎样”，弹性分析需要的是“哪些组合最危险”。批量做法是写一个外层脚本，对初始断线集合做枚举或抽样，把每次仿真结果收进结构体数组。
function res = batch_cascade(caseName, k, Nsample, opt) % 批量抽样N-k初始断线，做级联仿真 % k: 初始断线数量 % Nsample: 抽样次数 % 返回res结构体数组，字段包括edns、tripHistory、status mpc0 = loadcase(caseName); nl = size(mpc0.branch, 1); % 支路总数 baseLoad = sum(mpc0.bus(:, 3)); % 初始总有功负荷 rng(2025); % 固定随机种子，保证可复现 res = struct('init', {}, 'edns', {}, 'status', {}, 'nTrips', {}); for s = 1:Nsample init = randperm(nl, k); % 随机抽取k条支路 out = ac_cascade_sim(caseName, init, opt); res(s).init = init; res(s).edns = out.edns / baseLoad; % 归一化到标幺 res(s).status = out.status; res(s).nTrips = numel([out.tripHistory{:}]); end % 按EDNS降序排列，头部就是高风险场景 [~, order] = sort([res.edns], 'descend'); res = res(order); end这个脚本有两个工程细节。第一是固定随机种子rng(2025)，蒙特卡洛结果必须可复现，否则两次跑出不同结论，很难跟人解释。第二是randperm(nl, k)是无放回抽样，不会出现重复断线同一支路的情况，符合实际。输出按EDNS降序排，排在前面的就是需要进一步人工分析的关键场景。Nsample的取值可按算例规模掌握：case39上取2000次足够看出趋势，case118上建议5000到10000次，因为支路数量增加后组合空间增长很快。
4. 从仿真曲线到电网弹性结论：如何用级联结果排序关键线路
4.1 定义“关键线路”：初始断线频率与损失程度双重排序
级联仿真的原始结果是一堆tripHistory和edns，它们还不能直接指导加固决策。工程上用得最多的是“关键线路排序”：先定义每条支路在随机N-2场景中的出现频率，再乘上该场景的EDNS，得到“脆弱度评分”。这样排序能把“每次都导致大停电的线路”和“频繁出现但后果不严重的线路”区分开。
function [score, cnt] = line_vulnerability(res, nl) % 从批量仿真结果中统计每条支路的脆弱度 % res: batch_cascade的输出 % nl: 支路总数 % score: 每条支路的脆弱度评分 = 参与度 x 平均EDNS cnt = zeros(nl, 1); ednsSum = zeros(nl, 1); for s = 1:numel(res) % 初始断线里的支路 init = res(s).init; cnt(init) = cnt(init) + 1; ednsSum(init) = ednsSum(init) + res(s).edns; end score = zeros(nl, 1); avgEdns = ednsSum ./ max(cnt, 1); appearRate = cnt / numel(res); % 脆弱度 = 出现频率 x 平均损失 score = appearRate .* avgEdns; score(cnt == 0) = 0; end这段统计做了两件事：一是累加每条支路被当作初始断线的次数和对应场景的EDNS；二是把“出现频率”和“平均损失”相乘，得到脆弱度。纯按损失排序会把一些很少断但一断就全黑的支路排太靠前，纯按频率排序会漏掉小概率大后果事件，两者相乘是工程上的折中。拿到score后，我一般把排序前10%的支路挑出来，手工查看它们的物理位置和连接关系，再决定加固优先级。
4.2 韧性曲线与韧性三角形面积的计算
级联仿真的另一个输出是“时间轴上的系统状态”。把每轮的负荷损失累加绘成曲线，就是韧性曲线。这里要分清两个概念：模型输出的“级联轮次”本身不是真实时间，它只是一个事件序列。要画时间轴，需要给每轮分配一个持续时间。最简单的做法是假设每轮保护动作间隔为常数（比如30秒），更精细的做法是给不同切除动作分配不同时间，但这需要保护定值数据，通常拿不到。
function area = resilience_area(lostLoadMW, dtPerRound) % 计算韧性三角形面积：负荷损失对时间的积分 % lostLoadMW: 每轮损失的负荷向量 % dtPerRound: 每轮持续时间, 秒 t = (0:numel(lostLoadMW)-1) * dtPerRound; area = trapz(t, lostLoadMW); % 面积单位: MW*秒 endtrapz是MATLAB里最容易被低估的函数，它对非均匀时间步长同样有效，适合用来处理后续你可能会遇到的“不同级联阶段时间步长不同”的情况。这个面积值相当于事件造成的“能量损失”，比单纯看最大损失更能反映事件对用户的总体影响。如果你要跟同行对比，记得统一dtPerRound的取值并写进报告，否则不同格式之间比面积没有意义。
4.3 机组越限和多孤岛场景的处理边界
第3章的简化实现里，潮流不收敛就按“全部负荷损失”处理，这在网络被切成多个孤岛时过于粗糙。实际电网断线后常常形成多个孤岛，其中一部分孤岛内有发电机组、且能自行维持电压和频率，另一部分孤岛则全部失去电源直接黑掉。把“部分存活”的场景计成“全黑”，会高估EDNS，进而让关键线路排序失真。
处理办法是：在runpf不收敛时，不要急着判全黑，先用图连通性分析找出所有连通分量；对每个分量检查是否包含发电机；无源分量直接判黑，有源分量单独提取出来做一次潮流。MATPOWER的ext2int机制在内部会处理孤岛节点重排，但对外暴露的接口在不同版本有差异，所以我的做法是自己在主循环外做一个预处理，先把网络抽象成图结构再逐分量求解。这一步虽然代码量不大，但逻辑分支多，建议单独写一个函数而不是堆在主循环里。第一次做多孤岛分析时一定要先画图检查孤岛划分是否正确——这一步也是整个模型里最让我觉得“玄学”的地方：看起来正确的矩阵操作，结果却经常因为在孤岛边界上少处理了一条支路而翻车。
5. 避坑：MATPOWER级联仿真里最容易翻车的5个地方
5.1 现象：runpf返回success=0且报“矩阵奇异”
原因：断线后系统被切成了多个孤岛，其中某个孤岛没有发电机（或没有松弛母线），潮流方程组的雅可比矩阵奇异。这不是潮流算法的问题，是模型本身在该拓扑下无解。
解决：先做连通分量分析。凡是不含发电机的孤岛，直接按全黑处理并从拓扑中摘除；含发电机的孤岛单独提取成独立算例，分别求解后汇总结果。最简单的一步是：断开线路后，先执行conncomp检查每个连通分量里的GEN_BUS分布，再把无源分量甩掉。
5.2 现象：某条支路的loading变成无穷大或被误判为严重过载
原因：RATE_A填的是0，被当成“限额为零”，导致abs(PF)/0=Inf。大量自定义算例和从其他格式转换来的数据里，热稳限额并没有全部填入，RATE_A为0的含义是“未限制”，不是“不能带任何功率”。
解决：在潮流计算前做一次数据清洗：rateA(rateA <= 0) = Inf;。注意同时把RATE_B和RATE_C一并检查，因为后续如果做N-1越限统计，RATE_B通常是短时限额，比RATE_A高约20%，你可能会用到它做短时过载评估。
5.3 现象：输出结果的支路顺序和case文件对不上，切错了线
原因：runpf返回的是MATPOWER“内部节点重排后”的结果（包括孤岛重排和节点编号压缩），与loadcase读入时的顺序不一致。直接拿内部顺序的branch(:,12)和外部顺序的rateA做除法，结果张冠李戴。
解决：每次runpf后先做int2ext(mpc2)还原成外部顺序，再取潮流结果。这条属于“不踩永远不知道”的坑，我早期在case39上吃过亏——因为39节点系统支路较多，顺序错乱一两行很难发现，但结果就是关键线路排序表整个不可信。血泪经验是：级联仿真里对任何“支路编号”都要有一个明确的约定，从头到尾只用外部顺序，中间转换一律在取数环节完成。
5.4 现象：重载条件下AC潮流反复不收敛，但系统并非真的崩溃
原因：牛顿法在接近电压崩溃点时迭代发散，这更多是数值问题而非物理问题。尤其是在级联仿真的中间轮次，系统短暂处于高负载状态，实际还可能通过电压调节稳定下来，但牛顿法已经无法收敛。
解决：新版本MATPOWER可以在mpoption里切换潮流算法，比如用pf.alg从NR改成FD（快速分解法），或者给发电机母线设更合理的电压初值VG。另一个更直接的办法是给所有负荷母线加一个很小的并联电纳（比如1e-5），改善迭代初值附近的雅可比矩阵条件数。但这会让潮流结果有微小偏差，建议只把它当作中间轮次的计算手段，最终稳定态仍用严格不修改参数的runpf复算一遍。
5.5 现象：蒙特卡洛跑几千上万次，MATLAB越来越慢甚至卡死
原因：每轮runpf都在新开求解器、打印日志、做大量结构体复制；批量循环里又没有专门关掉输出和检查选项。
解决：批量仿真专用的mpoption要设成mpoption('out.all', 0, 'verbose', 0, 'pf.tol', 1e-6)。另外把算例结构体一次性载入内存，所有轮次复用，不要每轮loadcase。如果装了Parallel Computing Toolbox，把外层for改成parfor，注意结构体数组的写入要做成索引式，别让MATLAB因为广播变量把内存撑爆。这一步能让批量耗时从一小时压到几分钟，属于投入产出比最高的优化。
6. 验证与进阶：用case9校准模型，再上case39做关键线路筛选
6.1 验证步骤：先复现经典断线场景，再做收敛性检查
任何级联模型在用于自己的算例之前，都值得先做一次“已知答案”验证。做法是拿case9这个只有9节点、3台发电机、6条负荷支路的小算例，人为指定一组断线，手工推演一遍潮流转移方向，再和模型输出的每轮过载支路对比。case9的总有功负荷大约在300MW量级，几条主要输电线的传输功率都能从潮流结果里直接读出来，手推误差不会太大。这个验证方法虽然土，但能一次性暴露两三类错误：支路顺序不对、过载阈值方向反了、切除后又恢复运行（BR_STATUS没保持）。跑通case9之后再上case39，这时你已经有把握确认主循环逻辑本身没有低级错误。
6.2 一个具体技巧：用双阈值跑敏感性，输出“断线级联表”
级联仿真最怕的不是跑不出来，而是跑出来以后没法解释给调度人员听。我习惯在报告阶段加一个“双阈值敏感性”分析：把过载阈值分别设成1.0和1.2跑两组批量仿真，对比同一初始断线集合下的最终EDNS差异。若两者差异很大，说明系统里存在大量“刚好站在过载边缘”的支路，保护定值的小调整就能显著改变事故后果；若差异很小，说明系统是硬性结构脆弱，保护定值帮不上忙。这种对比比单独一组仿真结果有说服力得多，它把“偶然性”从级联过程里分离出来了。
检查模型可信度的另一个习惯是打印“断线级联表”——把每次仿真里tripHistory的支路按轮次展开成表格，列是轮次，行是断线支路编号。这个表能让懂电网的人一眼看出级联传播路径是否合理，比如是否从重载区域向轻载区域蔓延，是否存在莫名其妙跨越大区域的跳闸。级联传播路径不合理，往往意味着算例数据有错，而不是模型问题。
我自己每次拿到一个新算例，一定先跑一遍case9的对照场景，再跑case39全批量扫描，最后才敢把模型挂到实际网络上。级联故障模型的黑匣子成分比普通潮流大得多——因为你永远不知道哪一轮的切除动作放大了误差。上面这套流程能把这个黑匣子打开一半，剩下一半靠你一遍遍和实际故障案例对照日志。这个模型在电网弹性分析上的价值不在计算本身，而在于它帮助你建立一个共识：哪些线路值得加固，哪些区域值得做孤岛运行预案。希望帮到你。
本文还有配套的精品资源，点击获取
