做电力系统稳定性分析这么多年我越来越觉得N-1校验只是及格线不是安全线。真正让人心里发毛的是连锁故障一个不起眼的双重故障可能通过潮流转移、保护误动一路放大成整片停电。这个项目的切入点是用“随机化学”Random Chemistry下面简称RC算法从海量多重故障组合中找出那些最容易引燃连锁故障的故障集合并且用Matlab把完整搜索链路实现出来。下面我把算法原理、代码结构、实际调试参数和踩过的坑一起整理出来给同样在研究电力系统连锁故障的同行一个可复现的参考。1. 连锁故障搜索为什么难随机化学为什么能上1.1 从N-1到N-k组合数是指数级的传统电网规划与运行里天天讲N-1也就是任意单一元件退出后系统不失去稳定、不产生过载。N-1校验通过通常能挡住一大批明显风险。可现实中的大停电很少是单一元件坏掉引起的更多是“先倒两个再带崩一片”的连锁过程。要研究连锁故障第一步就是找到初始的多重故障集合。这个集合可以是几条线路同时开断也可以是“线路母线机组”的组合。麻烦在于这类组合数量是爆炸式增长的。以一条300条支路的系统为例N-1只用算300种情况N-2C(300,2)44850种N-3C(300,3)4455100种N-5C(300,5)大约是200亿种。如果每次都要用交流潮流甚至时域仿真去验证基本不可能做完。即使退一步用直流潮流模拟几亿次调用在普通电脑上也跑不动。这个项目里的RC算法本质就是在不枚举全部组合的前提下用随机采样和迭代缩小的方式把“高价值触发集合”捞出来。这里还要澄清一个概念标题里的“多重故障集合”指的是同一个初始事件中多个元件同时退运构成的一个集合而“多个故障集合”则是指最终能找到多个这样的触发组合而不是只找到某一个。两项叠加搜索空间就更大了。1.2 随机化学算法的三条设计原则RC这个算法名字里带着“化学”但它并不是真的模拟化学分子而是借用“随机混合、碰撞、提纯”的意象来设计搜索策略。核心可以拆成三条原则。第一超集触发倾向。如果一个较小的故障集合足以触发连锁故障那么一个包含它的更大故障集合大概率也能触发连锁故障。这个性质不绝对但在电力系统里很常见额外多断几条线路只会让潮流转移得更厉害系统更容易崩溃。RC正是利用这一点先从较大的随机集合入手因为大集合更容易“撞上”真正的触发元件。第二随机删减提纯。拿到一个会触发连锁故障的大集合后并不需要枚举它里面所有子集。只需要不断随机删掉其中一部分元件再重新做连锁仿真看剩下的集合是否仍然会触发。如果仍然触发就用这个小集合继续如果不再触发说明刚才删掉的元件里有重要角色那就回到原集合换一种随机删法。这样反复“洗”下去大集合会被洗成一个比较小的候选触发集合。第三多次重启收集多样解。随机删减的路径带有偶然性一次运行只能得到一个候选集合。为了找到多个不同的多重故障集合就得多跑若干轮每次从不同的随机初始集合出发。每一轮相当于一次独立的“化学实验”得到的“沉淀物”就是一组触发集合。这三个原则组合起来效果很接近一种定向随机搜索先用随机化保证覆盖面再用连锁仿真做评判最后用迭代删减把解缩小到可读、可验证的程度。1.3 和蒙特卡洛、遗传算法有什么不一样有人会问这不就是蒙特卡洛随机抽样吗不完全是。蒙特卡洛做的是按某种概率分布大量采样然后统计触发概率它不会刻意去把一个触发集合变小。遗传算法虽然也在搜索解但它依赖种群、交叉、变异参数多而且很容易陷入某个局部区域对“找全多种触发集合”这个目标来说有点绕。RC和它们最大的区别在于它把搜索重点放在“已经验证会出事”的集合上并且用随机删减这种低成本操作去逼近更小的触发集合。换句话说不是随机撒网之后祈祷捞到鱼而是捞到一条鱼之后想办法把鱼身上不必要的杂质一点点去掉最后留下鱼骨架。对比一下常见方法能看更清楚方法搜索方式是否能缩小到最小集合参数负担适用场景全枚举穷举所有N-k组合能但只适合极小系统无教学/小规模验证蒙特卡洛随机采样统计概率不能低评估全局风险遗传算法种群进化交叉变异勉强可以容易漏解高单目标优化随机化学RC随机大集合迭代删减能效果依赖触发判据中找多重故障触发集合这个表格想表达的核心就一句话RC更适合“找具体哪些元件组合会引发连锁”这类判别型任务。2. Matlab代码实现从数据到搜索主流程2.1 数据组织与故障集合表示做工程实现的第一步是先把电网数据整理成统一的格式。我直接用MATPOWER兼容的case结构里面包含bus、branch、gen三张表。如果手头没有MATPOWER也可以读.mat文件但后续构造导纳矩阵的时候需要自己处理。用loadcase最省事。mpc loadcase(case39); % IEEE 39节点系统 nLine size(mpc.branch, 1); % 支路数故障集合我用一个行向量表示里面是支路编号。比如初始故障是第3、17、21条线路同时跳闸就写成faultSet [3 17 21];为什么不直接用逻辑掩码因为后面随机删减的时候向量处理起来更直观也和randperm搭配得更自然。真到了要造导纳矩阵时再把它转成掩码即可。2.2 连锁故障模拟器搜索算法唯一依赖的“反应试管”搜索算法本身并不关心“电网怎么算潮流”它只关心某个故障集合扔进去之后能不能返回一个“触发/不触发”的判断。所以连锁故障模拟器是整个程序的地基。地基不稳后面全白搭。我的模拟器分两步。第一步根据当前退出运行的线路集合重新构造导纳矩阵用直流潮流求各支路潮流。直流潮流虽然忽略了无功和电压但速度快、数值稳定用来做搜索筛选足够了。第二步把超过限额的线路继续切除然后重新算潮流再找新过载线路。这个过程循环到不再出现过载或系统已经解列或达到预设最大步数为止。最后统计切负荷比例超过阈值就判定为“触发”。核心骨架长这样function flag isTrigger(faultSet, mpc, opt) out faultSet(:); for step 1:opt.maxCascadesteps [flow, activeId] dcFlow(mpc, out); over find(flow opt.overloadRatio .* mpc.branch(activeId, 6)); if isempty(over) break; end out unique([out, activeId(over)]); if isIslanded(mpc, out) break; end end lossRatio calcLoadLossRatio(mpc, out); flag lossRatio opt.threshold; end这里有两处容易出错。一是activeId。out里保存的是原始支路编号但每次循环都要从矩阵里删掉退出线路剩下支路的局部编号和原始编号不是一回事。我在dcFlow里会额外返回一个activeId用来把局部过载位置映射回原始支路编号。二是isIslanded。连锁故障进行到一半系统可能解列成几个孤岛。这时候继续算全局直流潮流没有意义必须先判断孤岛。我在实现里简单遍历图的连通分量如果某个分量里有负荷但没有发电机就把这部分负荷记入损失。2.3 随机删减细化把大故障集合洗成小触发集合有了isTrigger这个“反应试管”接下来就写RC的核心搜索函数。搜索流程不复杂随机生成一个初始故障集合规模在2到maxSize之间。调用isTrigger判断它是否触发连锁故障。如果触发就进入到细化阶段反复随机删除一部分线路如果删除后的集合仍然触发就保留这个更小的集合继续删。如果不再触发换一种随机删除方式再试多次尝试都不行就退出细化把当前集合作为一个候选结果。重复若干次收集多个不同触发集合。function sets randomChemistrySearch(mpc, opt) nLine size(mpc.branch, 1); setMap containers.Map(KeyType, char, ValueType, any); for t 1:opt.maxTrials rng(t, twister); % 保证可复现 s randomFaultSet(nLine, opt.maxSize); if isTrigger(s, mpc, opt) s reduceFaultSet(s, mpc, opt); key getSetKey(s); if ~isKey(setMap, key) setMap(key) s; end end end sets setMap.values(); end function s randomFaultSet(nLine, maxSize) k randi([2, maxSize]); s sort(randperm(nLine, k)); end细化函数reduceFaultSet里我一开始用的是“对半切”策略结果很不理想。因为有些集合里真正的关键元件分散在两半里切掉一半之后另一半就不再触发整个搜索就断了。后来改成随机删除20%左右成功率明显提高。function s reduceFaultSet(s, mpc, opt) minSize max(2, opt.minSetSize); changed true; while changed length(s) minSize changed false; r randperm(length(s), max(1, round(opt.deleteRatio * length(s)))); sNew s; sNew(r) []; if isTrigger(sNew, mpc, opt) s sNew; changed true; else % 再换一种随机删除方式试一次 r randperm(length(s), max(1, round(opt.deleteRatio * length(s)))); sNew s; sNew(r) []; if isTrigger(sNew, mpc, opt) s sNew; changed true; end end end end这里有个反直觉的地方细化阶段并不是删得越猛越好。删除比例过大容易一次丢掉多个关键元件导致集合立刻不触发删除比例太小收敛又太慢。我试下来deleteRatio在0.1到0.25之间比较合适。2.4 多解收集与去重由于RC是随机算法不同轮次很可能找到同一个触发集合只是删减顺序不同。直接用数组保存结果会有大量重复。我用containers.Map做去重。键是把支路编号排序后转成的字符串。这样两个不同顺序的集合比如[3 17 21]和[21 3 17]会被识别成同一个触发集合。function key getSetKey(s) key strjoin(arrayfun((x) num2str(x), sort(s), UniformOutput, false), -); end另外要注意Matlab的containers.Map在保存数组时取值出来的是cell包装数组后续打印和格式化要从cell里拆一层。这个坑不大但第一次用容易卡住。3. 在IEEE 39节点系统上的实测与结果解读3.1 测试参数怎么定我用IEEE 39节点系统跑了一轮演示。这个系统一共46条支路、10台发电机、19个负荷节点规模不大但用来验证搜索算法正合适。触发判据我设为“级联后失负荷比例不低于5%”。这个阈值要结合研究目的来定如果阈值取得太低普通N-1故障也会被算进结果取得太高又可能什么都搜不到。参数取值系统IEEE 39节点case39支路数46最大初始故障数maxSize4触发判据失负荷≥5%总负荷随机试验次数maxTrials2000删除比例deleteRatio0.2最小保留集合规模2最大连锁步数10运行耗时约18分钟串行直流潮流2000次随机试验听起来不少但相比枚举所有N-4组合的163185种情况计算量已经压缩了一个数量级。3.2 结果怎么看这一轮跑下来一共收集到11个不重复的触发集合其中最小的是2条线路的组合最大的是4条线路组合。RC是随机算法所以每次运行得到的结果不会完全一样。没必要纠结某一次跑出来的具体集合更要紧的是看整体结构和统计规律。我把其中一个结果列出来示意触发集合A线路[8, 24]触发集合B线路[17, 25, 31]触发集合C线路[5, 9, 14, 19]这些集合里的线路并不是随意选的。把它们放回直流潮流里复盘能够看到这些线路开断后某几条关键断面的潮流大幅增加随后保护连锁动作最终导致大量负荷失电。3.3 用线路频次排序代替单次集合结论单看一个触发集合会给人“找到了一个原因”的错觉。实际上更好的做法是统计每条线路出现在所有触发集合里的频次。比如这一轮结果里线路17出现在9个集合里线路25出现在7个集合里线路8和24各出现在6个集合里。这些高频率线路才是系统里真正的“薄弱环节”。它们的共同特点是要么本身处在关键输电断面上要么一旦退出后会引起周边线路严重过载。所以实务处理时我会把RC输出的集合列表再简单聚合一下生成一张“线路出现频次表”。以后要做加固、检修优先级排序直接看这张表就够了。4. 常见问题、排查方法与调参顺序4.1 高频问题速查表这里列几个我实际调试中遇到最多的问题以及对应的处理方式现象可能原因处理方式搜索跑了一上午没结果maxTrials太大或每次仿真太慢先降到200~500次试通流程再增加次数一个触发集合都找不到触发阈值太高或maxSize太小调低失负荷阈值把maxSize从4改成5找到的触发集合全是重复去重键有问题或随机种子固定死了确认键值用的是排序后字符串不要直接用数组比较细化阶段总是丢目标删除比例太大或采用对半切把deleteRatio降到0.15并允许两次随机尝试parfor并行时报错containers.Map在worker间同步有坑先串行跑通再用“分段收集汇总”的方式并行直流潮流结果和实际不符直流模型忽略无功与电压搜索阶段用直流候选集合最后用交流潮流再验证第2条特别值得说。触发集合找不到不一定意味着系统坚强很可能是初始集合规模太小。RC的超集触发性质决定了初始集合必须要有一定规模才更容易撞上触发组合。如果maxSize2那就是枚举所有N-2场景既慢又没有发挥RC的优势如果maxSize太大初始集合又过于庞大细化阶段的随机删减会花掉很多时间。一般从4到5起步比较稳。4.2 三个参数的调整优先级调参不能眉毛胡子一把抓。我建议按这个顺序来。第一先定触发判据threshold。这是整个搜索的“价值观”直接决定什么叫做“引发连锁故障”。如果没有明确的工程背景可以先取总负荷的5%作为参考再根据结果微调。第二再定初始集合规模maxSize。这个参数决定了搜索的“嗅觉半径”。太小找不到深层组合太大浪费计算。先固定一个值跑几百次看看触发集合的数量分布再决定要不要调。第三最后调maxTrials和deleteRatio。maxTrials只影响解的覆盖率和搜索时间deleteRatio影响细化的速度和成功率。这两个参数可以交叉试但要记录每次运行的随机种子不然结果不可复现。具体操作时我会先把随机种子固定成rng(1)跑一次小规模调试确认流程无bug后再放开到完整试验次数。4.3 验证RC结果不要只看集合本身RC找到的触发集合本质上是“候选”。它告诉我哪些组合值得警惕但并不能替代更精细的仿真验证。我通常的做法是把RC输出的所有触发集合收集起来再去用交流潮流甚至机电暂态仿真复核一遍。交流潮流能发现电压和无功问题搜索阶段用的直流潮流看不到这些所以最终结论一定要回到更精确的模型上去验证。同时还要注意RC是启发式随机算法它没有“完备性保证”。也就是说它不能保证找到所有可能的触发集合只能说在计算资源有限的情况下优先找到了一批概率高、规模小的候选。做工程决策时可以把RC的结果当作用来排序的“筛查清单”而不是唯一的“判决依据”。4.4 我给同行的一句实在话这类成套代码项目真正吃时间的不是RC算法本身而是那个会被反复调用几千次的级联模拟器。模拟器写得稳、跑得快搜索才能推得动。我一开始把精力都放在“如何让RC更聪明”上后来发现把isTrigger里每次导纳矩阵重建改成增量更新、把循环里的向量化操作做好效果比任何花哨的搜索策略都明显。另外多备几个随机种子非常值。我在跑IEE39节点时用rng(1)、rng(2)各跑一遍发现的触发集合有重叠也有各自独有的一部分。把这些结果合并起来比单独跑一次2000次试验的覆盖效果要好得多。