1. 从“组合爆炸”说起为什么连锁故障需要找多重故障集合1.1 一个n-2事故为什么会变成全网失稳在电力系统安稳分析里我越来越觉得最难回答的问题不是“某个元件坏了怎么办”而是“哪一组同时发生的故障会直接把系统推进连锁崩溃”。单重故障n-1扫一遍结论清晰日常运行方式基本都覆盖住了。可现实中的大停电绝大多数都不是一个元件自己坏掉而是两个、三个甚至四个故障同时出现或者在一个极短的窗口里接连触发。比如我在仿真演练里复现过类似场景两条220kV线路因为同塔双回故障同时跳闸潮流瞬间压到相邻一条330kV断面上。表面看断面负载率只从40%升到85%还在稳定限额以内但实际系统里对侧水电机组出力跟踪迟缓无功支撑不足母线电压一路掉到0.88p.u.低压减载切掉大量负荷高频机组又跳了两台最后触发低频减载形成大面积失稳。这种事故的本质是“多重故障集合”共同作用的结果。单独断开其中任何一条线路系统都能扛住但把它们放在一起连锁反应就像被点燃的引线一步一步把故障放大。所以电力系统连锁故障研究里识别“哪些多重故障集合是危险集合”比单纯计算某个断面的N-1通过率要关键得多。这也是我这篇笔记想聊清楚的事情怎么用随机化学算法把这种藏在海量组合里的多重故障集合找出来并且用Matlab完整实现、验证、调参。1.2 暴力枚举在真实电网里走不通很多人第一次接触连锁故障时第一反应都是“那就全枚举呗”。对一个只有46条支路的IEEE 39节点系统n-2的组合数是C(46,2)1035n-3是C(46,3)15180n-4是C(46,4)163185。如果是省级电网220kV及以上线路动辄2000条光n-2就有大约200万个组合。每个组合跑一次带保护动作、切负荷、低频减载的连锁仿真就算优化到5毫秒一次单核跑完n-2也要好几个小时n-3直接没法看。更麻烦的是连锁故障的后果不是线性的。A线路和B线路单独断开都没事不代表它们同时断开没事同理三回线里有两回是构成危险集合的必要元素第三回只是“陪跑”你也不知道该从哪个方向去枚举。早年我们常用启发式筛选比如找初始潮流最重、介数中心性最高的线路集合。但我的经验是真正危险的集合往往长在系统最弱的位置而不是最重的位置。一条轻载长线断掉后潮流转移到一条短线上再叠加一台机组出力受限电压越限、保护连锁动作这些组合用“负载率排序”根本捞不到。于是就需要一类不那么依赖人工洞察力的搜索算法随机化学算法就是其中很典型的一种。1.3 随机化学算法提供的是一条概率性的快速通道第一次看到“随机化学”这个名字我以为是模拟分子动力学的新材料算法。后来读了相关文献才发现它把每个候选故障集合当成一个“分子”集合里的每个故障元件当成“原子”。分子被投入反应器做一次连锁故障仿真就相当于一次“化学反应”发生严重连锁故障说明这个分子内含有一个“危险原子团”不发生反应说明这个组合本身不敏感。搜索的目标就是反复让分子分裂、合并、引入新的元素最终把那些必须同时存在才会触发大停电的“最小原子团”分离出来。这个思路最大的优点是你不需要预先知道危险集合的阶数。很多方法要求你先假设故障是n-2还是n-3再在固定阶数里搜随机化学算法则可以从一个较大的随机集合出发通过分裂和重组自动把集合缩到最小。它不追求像穷举那样“一个不漏”而是用可控的仿真次数换取一个很高的发现概率。对于真实的大电网这是唯一现实的做法。这篇笔记就是我完整复现这套算法的记录覆盖数学模型、Matlab实现、参数调节以及各种翻车现场适合正在做电力系统安稳分析、可靠性评估或者刚接触复杂系统搜索算法的同学参考。2. 随机化学算法的原理拆解与关键参数2.1 先建立两个定义故障集合与关键集合为了把问题说严谨先定义符号。假设电力系统里有n个可开断元件包括线路、变压器和部分机组。任意一个多重故障用一个故障集合S表示|S|就是故障阶数。对每个S我们把S内的元件预先退出运行然后跑一次连锁故障仿真得到失负荷比例L(S)。设定一个阈值Lthr如果L(S) ≥ Lthr就称S为“关键集合”否则叫“非关键集合”。工程上更关心的是“最小关键集合”一个集合M是关键集合但M的任何一个真子集都不是关键集合。这是最典型的“多重故障集合”形式。注意它和数学上的“最小集合”不完全一样因为连锁仿真器本身是非线性、非单调的。实际中可能出现集合A关键B也关键但A∪B反而不关键。这种反直觉现象在潮流转移和自动装置动作后并不罕见所以不能用解析方法来求只能靠仿真验证。有了这两个定义随机化学算法的工作就非常聚焦了从随机初始集合里找关键集合再把关键集合尽量缩小到“最小”的状态。缩得越小对调度运行人员的参考价值越大。一个包含5条线的关键集合和一条“缺了它就不行”的3条线最小集合后者才能算是真正的导火索。2.2 反应、分裂、合并三轮搜索循环我在实现时把算法拆成四个模块初始化反应池、反应测试、分裂与合并、收敛输出。初始化阶段随机生成若干初始候选故障集合作为“分子”。集合大小在2到Kmax之间随机取Kmax一般设4到6。因为从工程角度看超过n-6的初始故障集合概率极低再往上搜即便找到也没有太多操作意义。反应测试就是对候选集合调用一次连锁故障仿真。如果没达到阈值说明这个“分子”不活泼放回反应池继续放着如果达到了阈值就进入分裂流程。分裂时最常用的是对半分把集合S随机分成两份T和U分别做仿真。如果某一半仍然关键说明危险原子团被完整保留在这一半里那就抛弃另一半继续对关键的那一半做下一次分裂。如果两半都不关键说明危险原子团被切开了这时候就要做“合并”从集合外部随机补入一些元件或者把T和U按一定比例重新组合成新分子再放进反应器试一次。这个过程听起来很朴实但它能把“指数级”的搜索空间压下来。对每个关键集合典型的缩小过程大约需要O(klogk)次仿真而不是对所有真子集做O(2^k)次检查。而且因为每次分裂用的随机划分都独立重复多轮之后捕捉到某个关键集合的概率会迅速升高。这也解释了为什么随机化学算法很依赖“反复试”它用概率换取了速度。2.3 参数如何取舍从初始集合大小到分裂比例用随机化学算法最忌讳的是拿到代码就乱跑。我整理了一张参数表基本都是我在实验里验证过的范围参数含义常用范围说明Kmax初始故障集合最大阶数4~6太小会漏高阶层集合太大会拉长单次缩小过程splitRatio分裂比例0.4~0.60.5最均衡偏小容易保留核心偏大缩得慢maxShrink单个关键集合最多分裂轮数10~30防止某个奇怪集合无限循环recombineFrac合并时从外部补入元素的比例0.3~0.5太大集合会膨胀太小则跳不出局部outerIters外层随机初始集合个数100~1000越多发现概率越高耗时线性增长lossThr关键集合判定的失负荷阈值5%~20%需要根据仿真严重度分布提前观察seed随机数种子固定值跑实验必须固定否则无法复现这里重点说splitRatio。假设当前集合大小是10危险核心是3个元件。随机选5个作为子集T这5个恰好包含全部3个核心元件的概率是C(7,2)/C(10,5)≈8.3%。因为算法会同时检查T和它的补集U所以单轮找到核心的完整概率大约是16.7%。如果你把splitRatio调到0.2核心被保留在单侧子集里的概率更低连续多轮失败的次数增加搜索效率会明显下降。反过来如果splitRatio调到0.8虽然每次测试的子集包含核心的概率变高了但集合收缩速度太慢需要很多轮才能从10缩小到3。实际跑下来0.5是最省心的起点。2.4 别忽视“假极小”概率性带来的边界随机化学算法本质上是概率性的这一点必须接受。它不能保证找到所有关键集合也不能保证输出的集合一定是最小关键集合。我遇到过不少次这样的情况算法停在一个4元件组合上任何随机分裂都不关键了于是程序把它当成“最小集合”输出但后来我手动做了全子集校验发现去掉其中一个元件后仍然关键。这就是典型的“假极小”不是算法坏只是它没有完备地验证所有子集。解决办法也很直接在所有输出结果上补一步穷举子集校验。对每个输出集合S如果|S|不超过5就把S的所有真子集都跑一遍仿真。因为能走到最后输出的集合数量通常只有几十个这一步计算量完全可控。校验之后把不满足极小性的集合替换成它的关键真子集或者直接丢弃。这个“收尾动作”在论文里往往写得很少工程里却特别重要否则你给调度人员的结果里会混进一堆冗余故障。3. Matlab实现把算法变成可以跑的数据结构3.1 代码工程结构实话说我第一次写这套代码时把所有函数堆在一个main.m里结果改一个参数要滚动三屏仿真器出问题也没法单独调试。后来重构成这样data/存放IEEE 39节点系统数据包括母线、支路、发电、负荷sim/放连锁故障仿真器包括dcFlowWithFault.m、simpleCascadeSim.mrc/放随机化学搜索主体包括randomChemistrySearch.m、shrinkCriticalSet.mvalidate/放暴力枚举、指标计算、结果比对脚本main.m只负责组装数据、设置参数、调用搜索和输出报告。这样做的最大好处是你可以随时替换仿真器。搜索算法只依赖一个“输入故障集合、返回失负荷比例”的函数句柄不关心内部是直流潮流还是交流潮流。先跑通一套简化仿真器再慢慢升级成更精细的AC模型搜索代码一行都不用改。3.2 连锁故障仿真器的取舍先DC后AC做大规模高并发搜索我强烈建议先上直流潮流模型。直流潮流把潮流计算线性化速度极快而且能正确反映支路开断后的主要潮流转移趋势。对“找关键集合”这个目标来说DC模型已经能筛掉绝大多数不敏感组合。等候选集合缩小到几十个以后再用AC潮流、或者含电压稳定指标的模型做二次确认。这种“两级仿真”策略在时间上非常划算。下面是我用来做一次支路开断直流潮流计算的函数骨架function [theta, Pflow] dcFlowWithFault(lines, Bbus0, Pbus, faultLines) % lines: Nx3矩阵, [母线i, 母线j, 电抗x] % Bbus0: 原始节点电纳矩阵, 这里只用于获取节点数 N size(Bbus0, 1); x lines(:, 3); avail ones(size(lines, 1), 1); avail(faultLines) 0; B sparse(N, N); for k find(avail). i lines(k, 1); j lines(k, 2); B(i, i) B(i, i) 1 / x(k); B(j, j) B(j, j) 1 / x(k); B(i, j) B(i, j) - 1 / x(k); B(j, i) B(j, i) - 1 / x(k); end ref 1; % 平衡节点 keep setdiff(1:N, ref); Bred B(keep, keep); Pred Pbus(keep); thetaRed Bred \ Pred; theta zeros(N, 1); theta(keep) thetaRed; Pflow zeros(size(lines, 1), 1); for k find(avail). i lines(k, 1); j lines(k, 2); Pflow(k) (theta(i) - theta(j)) / x(k); end end这段代码的实现思路是每次故障都从原始支路表重新组装节点电纳矩阵。对39节点这种规模速度快到可以忽略不计如果你要上几千节点的系统就得改成稀疏LU分解加局部更新避免每步重建全矩阵。连锁仿真的主体则是一个循环先算潮流找越限比例最高的支路把它按保护动作逻辑退出再重新算潮流直到没有新过载或者系统已经解列。这个模型足够粗暴但对随机化学搜索来讲最重要的是“同一个组合能否稳定触发连锁反应”而不是精确复现某个保护定值。3.3 用“仿真函数句柄”隔离搜索算法搜索算法的核心是一个“缩小关键集合”的函数。我把它设计成接收simFun句柄这样将来换成AC仿真器不需要改搜索代码。下面是我常用的收缩函数function Smin shrinkCriticalSet(S, simFun, nLines, params) S sort(S(:).); for r 1:params.maxShrink m numel(S); if m 1, break; end sp randperm(m); k max(1, round(params.splitRatio * m)); T S(sp(1:k)); U S(sp(k1:end)); if simFun(T) params.lossThr S sort(T); continue; end if ~isempty(U) simFun(U) params.lossThr S sort(U); continue; end % 两侧都不关键随机从外部补入一些元件再试 outside setdiff(1:nLines, S); if isempty(outside), break; end nAdd min(numel(outside), max(1, round(params.recombineFrac * m))); extra outside(randperm(numel(outside), nAdd)); S sort(union(S, extra)); end Smin S; end外层循环则负责生成足够多的随机初始集合重复调用shrinkCriticalSet并做一个去重。这里最需要注意的是去重逻辑。同样的线路组合可能从不同的初始集合出发都会收敛到它如果不做去重后续人工校核的列表会被重复结果塞满。我用的是排序后的线路编号做key再放到containers.Map里查重。3.4 性能优化随机数、缓存与并行这个算法很吃随机数但随机数也是最容易出bug的地方。我踩过一个坑用parfor并行跑外层循环时每个worker都是从同一个全局随机流里取数结果不同次实验之间结果对不上后来排查了半天才意识到是随机流冲突。现在我的习惯是在main.m入口固定一个种子如果是并行则给每个worker单独生成一个子流% 串行实验固定种子 rng(2024); % 并行时给每个worker一个独立子流 spmd rng(2024 labindex * 1000); end仿真缓存是另一个提速关键。同一个子集在分裂和合并过程中会被反复测试很多次尤其当集合缩小到三四个元件时重复率非常高。我在仿真外层包了一层缓存function isCritical cachedSim(cache, simFun, S) key sprintf(%d_, sort(S)); if cache.isKey(key) isCritical cache(key); else isCritical simFun(S); if cache.Count 20000 cache(key) isCritical; end end end用这个简单的记忆化策略我在39节点系统上把仿真次数减少了大约40%。缓存容量要注意控制否则几万个组合存下来内存也会涨得让人头疼。固定容量、只保留最近最关键的组合是在大系统上更稳妥的做法。4. 实验验证在IEEE 39节点系统上跑出来的结果4.1 实验设置与基准我用IEEE 39节点系统做验证系统共46条支路总负荷约6000MW。连锁仿真器用上面说的DC潮流模型过载倍数设为1.05每次最多依次退出8条支路失负荷阈值lossThr设为10%也就是600MW。搜索参数设置为Kmax4splitRatio0.5maxShrink20recombineFrac0.4outerIters200随机种子固定为2024。为了知道算法找得好不好我先把n-2和n-3全部暴力枚举了一遍故障范围组合数单次仿真耗时总耗时n-21035约5 ms约5 sn-315180约5 ms约76 sn-4163185约5 ms约13.6 minn-4全枚举13分钟做一次可以接受但我不想把它作为日常调参基准所以主要用n-2n-3的金标准来评估召回率。暴力枚举在这个规模上的价值不是“实际工程里也这么干”而是让我们对算法有了一个可以对比的尺子。没有这把尺子后面所有参数调优都是盲调。4.2 召回率、虚警率和“最小集合”占比这一轮实验的统计结果如下指标数值RC仿真次数4620输出候选关键集合88最终AC模型校验后仍然关键83最终校验后最小关键集合36暴力枚举金标准关键集合数78找回金标准集合数61召回率78.2%串行运行耗时约24秒和暴力枚举的16215个组合相比随机化学算法只用了不到5000次仿真就找回了大约78%的已知关键集合。剩下的22%主要是随机采样没有覆盖到的犄角旮旯组合或者在高阶n-4里才表现为关键的组合。最终输出的88个集合里有83个在更严格的AC模型下仍然关键虚警率约为5.7%。这个精度在工程上是能接受的因为人工校核只需要看83个集合而不是去翻一万多个组合。值得注意的是最终最小关键集合只有36个不到输出的一半。其余集合虽然关键但都能再删掉一两个元件。这正是我之前说的“假极小”现象。经过子集校验后我发现很多看起来需要4条线路同时断开的组合其实核心只有3条线路额外那条只是“顺带陪跑”。把这些冗余剪掉之后给运行方式人员看的报告清晰了很多。4.3 参数扫描记录哪些参数值得优先调做了一轮参数扫描后我最大的体会是splitRatio和outerIters最值得优先调其他参数只要在合理范围里影响都没有那么剧烈。参数设置召回率观察splitRatio0.361%单轮分裂保留核心的概率太低搜索效率明显退化splitRatio0.578%最优分裂均衡且收缩速度令人满意splitRatio0.770%收缩慢很多结果停在大集合上outerIters5058%初始候选太少漏检严重outerIters20078%性价比最高outerIters50086%召回率提升但耗时也升到约60秒recombineFrac0.266%外部补入元件太少容易困在当前集合里出不去recombineFrac0.478%能有效跳出分裂失败的局部状态recombineFrac0.674%集合经常膨胀需要更多轮数缩小从这张表能看出outerIters直接控制“采样密度”是召回率最敏感的参数splitRatio则决定收缩效率偏大或者偏小都会让结果退化recombineFrac则更像一个调节“跳出能力”的旋钮。实际工程里我会先固定splitRatio0.5然后跑一个小扫描确定合适的outerIters再根据结果微调recombineFrac。5. 常见问题与排坑实录5.1 一个关键集合都找不到先别怀疑算法我最早用随机化学算法跑一个简化14节点系统时整整一晚上什么都没找到一度以为算法实现有问题。后来逐步排查发现是我的lossThr0.3设得太高这个系统不管怎么切最多损失20%负荷算法永远判定不出关键集合。解决办法是从运行数据里先画一个“失负荷严重度分布”随机仿真几百个故障组合看失负荷比例主要集中在什么范围再据此选阈值。如果大量组合都在0%~10%那阈值设在10%可能就会漏掉很多潜在危险集合设在5%左右往往更合适。另外Kmax太小也会导致“一个都找不到”。如果真实危险核心是4条线路的组合而你最多只随机生成3条线路的初始集合那不管跑多少轮都不可能发现它。所以“一个都找不到”时先检查阈值再检查Kmax最后再看outerIters是不是小到几乎没有覆盖度。5.2 结果每次不一样不是程序写错了随机化学算法本身是随机的出现结果差异是正常现象。但如果差异大得离谱比如上一次找到30个关键集合下一次只找到5个那就要检查随机数管理。串行代码里入口处没有固定rng种子就会导致每次运行后台随机流不同。我现在的习惯是每次实验都固定一个种子并且把种子记录到结果文件名里。这样同一个实验重跑三次结果完全一致。如果你用了parfor问题会更隐蔽。parfor每个worker默认拥有独立的随机流但你无法保证它们和主进程的流之间没有重叠。更稳妥的做法是在spmd或parfeval里显式给每个worker赋不同的子流。做完这个改造后并行结果虽然没有完全跨次可复现但每次的统计指标会稳定很多。对于工程报告建议不要只跑一次而是用5个不同种子各跑一遍把结果做并集这样能显著提升召回率。5.3 仿真器成了瓶颈怎么把搜索速度提上去随机化学算法的单次仿真便宜但它会反复调用仿真器。在我的39节点实验里一次DC连锁仿真大约5毫秒4620次也就是24秒左右没有压力。但如果你把仿真器换成AC潮流加动态仿真一次可能变成200毫秒甚至更久这时候搜索速度就完全被仿真器卡死了。我的建议是分两步走。第一搜索阶段坚决用简化模型只要能区分“关键”和“非关键”就够了。第二给simFun做缓存前面提到的cachedSim可以砍掉四成重复计算如果组合规模更大可以考虑用一个固定大小的LRU缓存避免无限增长。第三尽量把连锁仿真的内部循环向量化。比如在计算所有支路过载率时一次性用矩阵运算算出Pflow和limits的比值不要一条条支路循环。39节点上这个优化不明显但节点数上去之后向量化能快一个数量级。5.4 问题速查表把这些经验整理成一张速查表方便以后排查现象可能原因处理方式一个关键集合都找不到阈值设太高 / Kmax太小 / 初始采样太少先画严重度分布调低lossThr增大Kmax和outerIters结果每次差异巨大随机种子未固定 / 并行随机流冲突固定rng种子给worker单独子流多seed并集输出的集合不是最小缺少子集校验对最终输出跑一遍全真子集仿真并替换搜索速度慢仿真器太重 / 重复仿真太多搜索阶段用DC模型加缓存向量化内部循环内存暴涨缓存无限增长给containers.Map设置容量上限或改用LRU策略回头验证时虚警高搜索模型和校核模型差异太大先校准简化模型的过载倍数和切负荷阈值再上搜索最后再分享一个我个人的习惯拿到一份关键集合结果后我不会直接拿去做运行方式调整而是先把每个集合翻译成“调度规程里的语言”。比如“这条500kV线路断开后对侧两个断面的潮流会发生怎样的转移”“哪些电厂的无功支撑会被触发”。这种翻译过程既是在验证算法结果的可操作性也经常能发现仿真器没有建模的二次设备行为。随机化学算法擅长把搜索空间压缩到人工能处理的规模但它不能替你把电网机理想清楚。小系统上先枚举、拿到召回率再上大系统这个流程我建议每一位准备用到这套方法的人都先跑一遍。