做电力系统可靠性分析这些年我越来越确认一个观点连锁故障的真正源头往往不是一个元件而是多个元件同时退出才把系统逼到崩溃边缘。这阵子我一直在用 Matlab 实现一种叫“随机化学”Random Chemistry的算法用来从海量 N-k 组合里快速识别能够引发连锁故障的多重故障集合。几个测试系统跑下来踩了不少坑也攒了些能直接复用的经验。这篇文章就把这个算法的原理、代码框架、实操流程以及我踩过的雷一次讲清楚给正在做电力系统脆弱性、连锁故障或防灾减灾相关课题的朋友做个参考。1. 这个课题到底在解决什么问题1.1 连锁故障的“点火开关”不在单个元件上先聊一个经常被忽略的事实电力系统正常运行的时候单个元件退出通常不会导致系统崩溃因为调度规程里早就留足了 N-1 安全冗余。也就是说任意一条线路、一台变压器或一台发电机跳掉之后剩下的网络仍然能保持稳定运行。真正麻烦的是多重故障。我复盘过不少国内外的大停电事件几乎每次都能在故障初始序列里找到不止一块“多米诺骨牌”可能是同一走廊上的几条线路同时被极端天气打断也可能是保护装置误动叠加一次检修隔离错误甚至是对侧变电站的母线故障连带切掉了多条馈线。这些初始的“多重故障集合”一旦形成就相当于直接绕过了系统的 N-1 防线让潮流传导路径瞬间收窄剩余元件承受的转移潮流立刻越限于是保护一个接一个动作级联跳闸就开始了。所以课题的核心问题从来不是“哪个元件最重要”而是“哪一组元件的组合退出会触发后续一连串的灾难”。单个元件的重要度可以用介数、潮流灵敏度之类的指标算但多重故障集合之间的耦合效应没法用简单的单要素排序来等价描述。这是整个研究存在的前提也是为什么我们需要专门的搜索算法。1.2 多重故障集合的筛选难在哪里从数学上看这个问题是一个典型的组合爆炸问题。假设系统里有 L 条线路那么 N-2 故障组合数是 C(L, 2)N-3 是 C(L, 3)随着 k 增大候选集合规模增长得非常快。一个中等规模的实际系统线路数量通常在 1000 到 5000 条量级C(3000, 3) 已经是 45 亿量级根本不可能逐组做连锁故障仿真。常规的做法是启发式搜索或者智能优化比如遗传算法、粒子群算法去找“最坏故障组合”。这类方法的问题在于它们倾向于收敛到局部最优解而实际电网里能够引发连锁故障的组合往往不是唯一的它们分布在故障空间的不同角落。如果只找到少数几个危险组合漏掉其他同样危险的候选对一个防御决策来说就是隐患。这里还有一个更隐蔽的难点连锁故障本身是一个复杂动态过程。同一个初始故障集合在不同的运行方式、不同的负荷水平、不同的保护定值下可能表现出完全不同的后果。这意味着我们不仅要在故障空间里搜还要在运行场景空间里验证。这两个维度叠加之后计算负担会非常重。所以任何实用的方法都必须采取“快速筛选 精确验证”的思路先用便宜的手段把大量无害组合过滤掉只留下少量可疑组合交给精细仿真。1.3 “随机化学”算法名字的由来我第一次听到“随机化学”Random Chemistry这个名字的时候还以为是搞化学动力学的人来跨界。实际上这个名字很形象在化学反应里大量分子在容器里随机碰撞绝大多数碰撞都只是擦肩而过只有极少数碰撞具有足够能量和正确取向才能打破化学键、引发新物质生成。把电力系统的故障组合想象成分子把连锁故障想象成“化学反应”事情就顺了。系统里有海量故障组合大部分组合像无效碰撞一样不会引发任何后果少部分组合是真正的“点火”组合。随机化学算法要做的事情就是通过随机碰撞的方式去快速发现这些稀少的危险组合并且一旦发现就尝试“浓缩”它——从组合里逐个抽掉元件看是不是仍然有触发能力从而把一个比较大的初始故障集合缩小到“准最小触发集合”。这个策略的价值在于它把“在巨大空间里找少数危险点”变成了“随机采样 快速验证 逐步降维”不需要穷举不需要梯度信息对离散的、非线性程度很高的连锁故障过程也适用。用 Matlab 实现起来思路非常直接剩下的问题就是工程细节怎么处理才能算得又快又稳。2. 随机化学算法的核心原理与设计取舍2.1 算法整体框架碰撞、点火、浓缩我实现的随机化学算法大致分成四个阶段每个阶段对应一个关键函数逻辑很清晰。第一阶段是“产生碰撞”从所有线路中随机抽取 k 条线路构成一个候选故障集合k 通常取 2 到 4对应 N-2、N-3、N-4 故障。这个阶段不需要任何电网知识纯粹是均匀随机采样保证故障空间各区域都有机会被覆盖到。第二阶段是“点火检测”把这个候选故障集合注入潮流计算让故障线路退出运行然后跑一轮连锁故障快速筛选器。筛选器输出一个布尔量这个初始故障集合是否引发了至少一次后续的级联越限或切负荷。如果是我们就认为这次“碰撞”是一次成功的“点火”。第三阶段是“浓缩”对每一个成功点火的故障集合尝试从中随机摘除一个元件然后重新跑点火检测。如果摘除之后仍然能引发连锁故障说明这个元件不是必须的如果摘除之后连锁反应消失了说明这个元件是关键成员。不断重复摘除过程直到集合里每个元件都不可或缺就得到了一个“准最小触发集合”。这一步是随机化学算法的精髓它把随机搜索的粗结果提炼成了有物理意义的答案。第四阶段是“统计与输出”把多轮采样得到的所有准最小触发集合汇总统计每个元件出现在高危集合里的频率。频率越高说明该元件越容易成为连锁故障扳机的一部分。最后按频率输出高风险元件排名同时保留完整的触发组合列表方便后续精确仿真验证。整个算法的计算开销主要集中在第二阶段的点火检测上其他阶段都是轻量级操作。所以“筛选器”的效率和精度直接决定了整个项目能不能用。2.2 为什么随机采样能顶住组合爆炸很多人看到“随机”两个字会本能怀疑海量组合空间里随便抽几万个样本真能找到那些稀少的危险组合吗这里有个简单的概率数学可以解释。假设危险组合占所有组合的比例是 p每次采样独立且均匀那么采样 M 次之后一次坏组合都没找到的概率是 (1-p)^M。换句话说至少找到一个危险组合的概率是[ P 1 - (1-p)^M ]举个具体数字。假设一个系统里真正能引发连锁故障的组合占比是 0.1%也就是 p 0.001那么采样 5000 次时理论上的命中概率是[ 1 - (0.999)^{5000} \approx 0.993 ]也就是说99.3% 的概率能至少命中一个危险组合。如果采样 10000 次概率更是逼近 99.99%。这个计算告诉我们一个关键结论随机化学算法的效率高度依赖于危险组合在空间中的密度。密度高少量样本就能命中密度极低比如低于百万分之一那么再增加样本数也不合算。在实际电网里经过电力系统运行方式多年优化危险的多重故障组合确实不算特别罕见尤其在重负荷、极端天气、检修叠加等特殊运行方式下危险组合密度会显著上升。所以随机化学算法在这种场景下非常适用。相比枚举法随机采样的优势在于计算量不随系统规模爆炸式增长而是由样本数 M 和单次筛选成本决定。你可以把 99% 的计算资源花在最可疑的少数样本上而不是把所有候选组合都过一遍。2.3 关键参数怎么定样本数、初始故障阶数、判定阈值实际操作中最影响结果质量的参数有三个。第一个是初始故障阶数 k。我建议根据系统本身的冗余水平来定。对大部分高压输电系统N-1 是安全底线N-2 已经算严重事件N-3 及以上通常对应极端外部灾害条件。所以我一般把 k 的下限定为 2上限根据系统规模和计算预算从 3 起步有余力再试 4。k 取得太大候选空间指数膨胀采样命中率反而会下降k 取太小又容易漏掉必须三个元件同时故障才能触发的情形。建议做一组 k2、3、4 的对比实验观察输出组合的重叠度就能找到适合你系统的平衡点。第二个是样本数 M。我推荐用上面的概率公式反推先粗略估计 p比如根据历史事故或筛试结果假设 p 在 0.1% 到 1% 之间然后设定目标命中概率不低于 95%算出 M 的下限。通常 M 在几千到几万量级比较合理样本太多计算时间不可接受太少又缺乏统计意义。第三个是连锁故障判定阈值这是最需要小心的地方。筛选器判断一条线路是否过载跳闸时需要比较实际潮流和线路容量。我建议把阈值设在额定容量的 0.8 到 1.0 之间不要直接用 1.0。原因有两个一是实际运行中线路保护定值一般留有余度过载到 100% 但不一定立即跳闸二是直流潮流模型本身忽略无功和电压计算结果到底准不准确得靠这一段余量补偿。阈值设得太保守会把大量无害组合误判成危险组合导致后续“浓缩”阶段效率变低设得太宽松又会漏掉真实危险组合。务必要结合你的系统模型做灵敏度实验。这些参数之间是耦合的不要一个一个孤立地调。我的习惯是先把 p 粗估出来定 M再把阈值从 0.8 到 1.2 扫一遍看同一批样本下命中数量的变化曲线选择曲线斜率开始变缓的那个点作为阈值。3. Matlab 代码实现的整体架构3.1 模块划分与数据流这一节说一下我用 Matlab 实现时如何组织代码。整个项目我拆成了四个模块结构清楚以后扩展也方便数据层负责读取电网数据生成节点导纳矩阵、线路容量向量、发电机出力与负荷数据。评估层实现直流潮流求解器和连锁故障快速筛选器是整个算法的心脏。搜索层实现随机采样、点火检测、故障集合浓缩、多轮统计。输出层生成风险元件排名、危险故障集合清单、可视化图表。数据流是单向的数据层把电网模型打包成一个结构体 mpc传给评估层评估层接收一个故障集合数组快速返回连锁故障判定结果搜索层通过循环调用评估层完成随机化学迭代输出层接收搜索层的统计结果生成报告。我强烈建议用 Matlab 的 struct 把电网参数组织成一个对象避免到处传散装变量。比如mpc.branch ...; % 支路表状态、from、to、r、x、b、容量 mpc.bus ...; % 节点表类型、注入功率、电压等 mpc.gen ...; % 发电机表节点位置、有功、无功、上下限 mpc.branch_status ones(size(mpc.branch, 1), 1);这样写的好处是函数签名简洁后续不管换算例还是增加模型细节改动都集中在数据层。3.2 电网模型与快速判定器DC 潮流与 PTDF连锁故障筛选器必须快但也不能太粗糙。综合考虑之后我选择以直流潮流DC Power Flow为核心搭建判定器配合功率传输分布因子PTDF做增量计算。这套组合的精度对连锁故障的宏观过程识别已经足够。直流潮流的本质是忽略无功、电压幅值和网损把有功潮流近似成线性方程。虽然精度不如交流潮流但它的优势非常突出计算稳定、速度快、不会出现交流潮流在严重越限时的不收敛问题。在筛选阶段我们需要的只是判断“故障后哪些线路会过载、是否可能引发下一轮跳闸”这个粒度上直流潮流的误差是可以接受的。PTDF 是更进一步的加速手段。它描述的是“某一节点注入功率变化时各线路潮流的变化量”是一个线性灵敏度矩阵。利用 PTDF当某条线路开断后不需要重新求解全系统潮流而是可以快速估算每条线路的潮流变化。实际编码时我建议在初始化阶段就调用 Matpower 或者自己写一个小函数把 PTDF 矩阵算好保存起来后续每次连锁故障迭代都复用这个矩阵能省非常多时间。我在实际测试中发现如果系统规模一两百节点直接用直流潮流重解也没问题但到了上千节点的系统反复调用线性求解器会越来越慢PTDF 增量法的收益就很明显了。3.3 连锁故障演化逻辑的实现要点连锁故障快速判定器的核心逻辑就是一个循环迭代过程。我把它总结成下面几条规则初始故障发生时将对应线路状态置为 0表示断开。重新计算潮流算出所有仍运行线路的有功潮流。把潮流超过容量阈值建议留 10% 到 20% 的裕度的线路标记为“待切除”。如果待切除集合为空说明系统达到了新的稳态连锁故障停止如果非空则一次性全部切除再回到第 2 步。为了防止无限循环设置最大迭代次数同时统计累计切负荷量或失去的连通性。当切负荷量超过某个阈值就认为连锁故障是灾难性的。实现时有一个小细节值得注意不要在一次迭代里只切一条过载线路因为实际电网中保护动作通常有先后顺序但时间尺度很短如果一条一条切计算次数会爆炸。批量切除在工程上更高效而且对识别高危故障集合的结论影响很小。另外还有一个必须处理的情况是网络解列。某些线路断开后系统会分裂成若干孤岛有些孤岛可能失去发电机或负荷。我的筛选器里特意加了一个孤岛检测函数如果某个孤岛里发电机总出力小于负荷则按比例削减负荷并把削减量计入损失指标。这个机制虽然简单但能让连锁故障的严重性评估更贴近实际避免“只算线路过载不算供电损失”的片面性。4. 从零跑通一遍完整流程4.1 准备数据从 Matpower 算例到邻接关系推荐直接用 Matpower 自带的算例来测试算法比如 case3939 节点系统、case118118 节点系统。这些算例是公开标准算例自带支路参数和负荷数据方便我们对标验证。第一步是把系统数据读进来然后提取我们需要的字段。Matpower 的 case 文件本质是一个结构体支路表里每一行代表一条线路或变压器字段包括 from 节点、to 节点、电阻、电抗、电纳和长期载流量限值。注意变压器支路的电抗可能为负且不一定有明确的容量限值需要单独处理。第二步是构建快速筛选所需的索引结构。比如br mpc.branch; L size(br, 1); from_idx br(:, 1); to_idx br(:, 2); rate br(:, 6); % 长期容量如果为 0 需要人为设定有些算例里支路容量字段是 0这在实际中不可能出现。我一般根据该支路在其他运行条件下的典型潮流来估算一个合理值或者直接参考潮流计算结果取基准运行潮流的 1.2 到 1.5 倍作为容量。这一步如果偷懒后面的连锁故障判定会失真得很厉害。第三步是算好 PTDF 矩阵。可以用 Matpower 内部函数也可以自己用 N 节点、L 线路的增广矩阵求伪逆核心是得到 L×N 的灵敏度矩阵。算一次后面全程复用。4.2 核心代码骨架随机采样与故障缩减下面给出一套可以直接改的代码骨架分三个函数。第一个函数是随机生成候选故障集合function cand generate_candidate(L, k, M) % 随机生成 M 组 k 阶故障组合 cand zeros(M, k); for i 1:M cand(i, :) randperm(L, k); end end第二个函数是“单次隔离检测”也就是上一节说的连锁故障筛选器function hit cascade_filter(branch_status, ptdf, flow_base, rate, threshold, max_iter) % branch_status: 1xL 的 0/1 向量1 表示线路投入 % ptdf: 功率传输分布因子矩阵 % flow_base: 初始运行方式下各线路潮流 % rate: 线路容量 % threshold: 过载判定阈值比如 0.9 或 1.0 status branch_status; flow flow_base; lost_load_ratio 0; for iter 1:max_iter % 根据 PTDF 计算当前开断下的潮流近似值 % 工程简化在基础潮流上叠加开断线路影响的增量 overload find(status 1 abs(flow) threshold * rate); if isempty(overload) break; end status(overload) 0; % 切除过载线路 % 更新潮流近似值这里省略了孤岛处理和潮流修正细节 end hit (lost_load_ratio 0.01) || any(status 0 flow_base 0); % 说明hit1 表示该故障集合引发了连锁后果 end实际使用中纯 PTDF 的线性叠加在连锁故障后期误差会累积所以我通常在这个函数里加一个判断如果开断线路数小于 5用 PTDF 快速估算一旦开断线路数增多就重新跑一次完整直流潮流来校正。这种“混合模式”比我一开始全部用 PTDF 要稳得多。第三个函数是故障集合浓缩function core shrink(fault_set, ...) % 输入一组能引发连锁故障的故障集合 % 尝试逐个剔除元件保留必须保留的最小集合 core fault_set; for tries 1:length(fault_set) r randi(length(core)); trial core; trial(r) []; if cascade_filter(set_to_status(trial), ...) core trial; % 删掉后仍危险保留精简集合 else % 删掉后不危险说明该元件必须保留 end end end注意这里的随机剔除顺序会影响最终得到的“最小集合”形态因为可能存在多个等价的最小集合。所以多跑几轮采样会对覆盖不同形态有帮助。4.3 结果输出关键故障集合与高风险元件排序算法跑完之后我最关心的输出有两个一个是具体的危险故障集合列表另一个是元件风险频率排序。危险故障集合列表建议保存成 cell 数组每一行是一组“准最小触发集合”旁边记录它的连锁后果指标比如切负荷比例、切除线路条数、迭代轮数。这样后续精确仿真可以直接按表复算不用再搜一次。元件风险频率排序用 Matlab 的直方图命令即可完成。把所有采样轮次里出现过的高危集合拆开统计每条线路出现在高危集合里的次数除以高危集合总数得到频率值。这个频率值虽然不是严格意义上的概率但作为一个相对风险指标非常直观。我还会额外输出一张热力图横轴是故障集合索引纵轴是线路编号热力值表示“该线路是否出现在该高危集合中”。这张图能很快看出高危集合是不是集中在某个区域或某几个关联走廊上对规划部门做风险防护很有用。5. 实战中躲不开的问题5.1 解出来的组合“不危险”或“不完整”这是我最开始跑随机化学算法时最容易遇到的困惑明明程序报告某个故障集合触发了连锁故障但换用更精确的交流潮流仿真验证时却发现后果并不严重甚至不连锁。反过来也有一些被筛选器放过的组合在精细仿真里却是真正的大麻烦。原因要从筛选器的模型精度找。直流潮流忽略电压和无功当系统面临严重有功缺额时可能低估电压失稳型连锁反应而当线路重载到接近热稳定极限时交流潮流中电压下降会导致线路潮流“虚高”直流模型又可能低估。所以我后来形成了固定动作任何被随机化学筛选器判定为高危的组合都要用交流潮流或者准稳态仿真复核一遍。另一个问题是“不完整”。随机采样本身决定了它只能保证以一定概率找到一部分危险组合不可能保证穷尽。如果某个危险组合在空间里覆盖得非常稀疏或者被参数设置卡在阈值边缘就可能始终没被采样到。解决思路是多换几个随机种子、提高样本数、扩大初始 k 的取值区间把几轮结果合并起来看。5.2 计算量还是太大怎么优化筛选器的循环迭代是最大的计算瓶颈。我踩过的坑和优化方案可以整理成几条不要每次迭代都重新组装导纳矩阵。把节点导纳矩阵拆成“固定部分 线路开断修正项”开断时只更新修正项能节省大量时间。能用 PTDF 增量计算的地方绝不要重新解方程组。只有开断数多了之后才切回完全潮流求解。样本循环尽量向量化。随机生成 M 组候选集合之后可以批次检查一部分组合而不是一组一组串行跑。Matlab 的 parfor 并行工具箱也值得用上M 很大时收益显著。如果初始故障是同一走廊的几条线路它们的连锁故障行为往往高度相似。我增加了一个“去重”机制在进入筛选前先判断故障集合之间是否过于相似相似的就按一个代表样本处理能少算不少。5.3 随机数影响结果稳定性你第一次跑和第二次跑输出结果可能不一样因为随机抽样本质如此。这个问题不算 bug但如果不处理对工程汇报很不友好。我的做法是固定随机种子。在脚本最前面加一句rng(2025);这样同一份数据、同一批参数任何时候重跑都是一模一样的结果。但这个做法也有代价如果种子选得不巧恰好覆盖不到某个角落结果反而有偏差。所以我建议综合处理出成果、做对比时固定种子保证可复现探索性研究时多跑几个种子把结果并集作为最终候选集合。如果并集数量太大再按风险频率排序截取前 N 组。5.4 使用经验与避坑清单最后列一份我实际使用过程中总结的避坑清单每一条都是真金白银换来的容量字段千万不能省。case 文件里很多支路容量是 0直接算会把所有路线都判为过载连锁故障一跑到底结果完全失真。阈值不要定得太低。0.7 以下会把系统里本来正常运行的重载线路全部当作过载切除连锁故障很快就蔓延全系统输出高危集合数量巨大没有区分度。缩减阶段要保留“组合语义”。有些线路单独断开时没有风险三个同时断开才有风险这种联合效应在浓缩步骤里非常脆弱。我建议在剔除元件时先随机试如果剔掉任何一个都不行就保留整个原始集合不要强行剥离只剩单根线路。考虑到保护隐性故障会更有意义。实际连锁故障里经常出现保护误动、拒动、隐性故障等问题但标准算例里没有这些信息。如果需要贴近现实可以在筛选器里按概率给每条线路加一个隐性故障标志虽然会让结果更加复杂但可信度会明显提升。算例规模从小起步。先用 case9、case39 验证算法逻辑再上 case118 和大系统不要一开始就挑战几千节点的算例否则定位 bug 会很难受。6. 把这个算法往后延伸的几条路6.1 两阶段筛选随机化学粗筛 精确仿真精算随机化学算法最大的价值在于“压缩候选空间”它本身不必输出最终结论。我目前最推荐的使用方式是把随机化学算法作为第一道筛子粗筛出几千组高危故障集合第二阶段再用交流潮流甚至时域仿真对这几千组逐组精算得到更准确的后果指标。这样既保住了采样效率又不牺牲模拟精度。 两阶段流程里第一阶段的阈值可以适当放宽宁可误报多一些也不要漏报第二阶段再严格核算把真正的高危集合定下来。这样虽然第一阶段耗时多一点但第二阶段的计算量大幅缩减整体效率仍然远高于直接枚举。6.2 与其他脆弱性指标结合随机采样是“盲搜”但如果能利用网络结构信息给采样概率加一点偏置往往能更快命中高危故障。比如可以先计算电气介数、线路关键度指标将高风险线路的采样权重调高同时保留一定比例的均匀采样覆盖全局形成“偏好随机采样 均匀随机采样”的混合机制。这样既保留随机化学算法的全局探索能力又能利用已有知识加速收敛。我试过在混合采样比例为 70% 偏好 30% 均匀时命中高危组合的效率比纯均匀采样提高了约 40% 左右代价是需要多算一次指标但整体看很划算。6.3 对我个人来说最有价值的用法对我来说这个算法最大价值不在于寻找“唯一最坏故障”而是输出一张“风险热区地图”告诉运行人员哪些线路组合需要重点关注、哪些地区容易成为连锁故障的火种源。调度部门可以据此优化运行方式、调整保护定值配合规划部门可以据此规划加强线路走廊、增加断面冗余。换句话说随机化学算法是一个很好的“风险认知工具”它帮你从噪声里看清事故模式而不是代替你做出最终决策。在做电网安全分析项目时我越来越明白真正难的往往不是找到答案而是知道该问什么问题。随机化学算法恰好帮我更快地提出那些关键问题。