在做配电网可靠性评估课题那会儿我最开始用的是故障枚举法。简单说就是把所有可能的元件故障组合列出来逐一算概率和影响听起来很直观可真把系统规模放大到十几条馈线、几十个负荷点之后组合数量直接失控程序的运行时间从几分钟膨胀到几小时而且很多时序因素根本塞不进去。后来换成了序贯蒙特卡洛模拟法配合Matlab做实现才真正把这件事跑通。这篇就聊聊序贯蒙特卡洛模拟法做配电网可靠性评估的核心思路以及Matlab代码落地时那些文档里不会明说但一定会遇到的细节。写这篇文章是想给正在做毕设、写小论文或者工作中需要算配电系统可靠性指标的朋友一条能直接照着走的路径。1. 为什么偏偏是序贯蒙特卡洛解析法的边界在哪里1.1 故障枚举法在配电网面前的尴尬配电网可靠性评估不是新话题经典教材里第一反应都是解析法。最小割集法、故障模式影响分析法这些方法对简单辐射网确实有效枚举出各个元件的故障事件再按开关逻辑把停电范围算出来。问题在于真实配电网的结构远不是一条主馈线带几个分支那么简单。从时间角度看负荷在变、分布式电源出力在变解析法需要把所有可能的状态组合都列为系统状态每个状态还要算稳态概率。元件数量一多状态空间就指数膨胀。我见过一个32节点的小型配网断路器、隔离开关、熔断器都算上完整的故障模式枚举表有几千行手工整理根本不可能。就算强行写代码枚举也会遇到另一个问题很多状态发生的概率极低枚举到它们纯粹是浪费计算资源。更麻烦的是解析法处理时序非常吃力。比如某个负荷点一年被停电三次每次两小时和一年被停电一次、持续六小时两种情况的SAIDI数值可能一样但对用户体验和考核指标的意义完全不同。解析法给出的是期望值它回答不了停电频率和停电时长分别怎么样这一类问题。而序贯蒙特卡洛天然是沿着时间轴走的这些问题都不在话下。1.2 序贯两个字的真正分量很多人把蒙特卡洛等同于随机抽样算期望这个理解没错但序贯两个字才是配电网场景下最值钱的部分。非序贯蒙特卡洛的做法是直接对系统各元件的状态做一次快照式抽样比如线路A在运行、线路B在检修、变压器C停运然后把这一整个状态拿去算影响。这样得到的结果虽然也能出指标但它忽略了状态是怎么演变来的——检修和故障是两回事先后顺序会影响停电范围。序贯蒙特卡洛的意义在于它把时间拆成了连续的事件流。程序从第1小时开始逐条设备用随机数生成它的正常运行多久→故障→修复→再正常运行序列所有设备的事件序列在时间轴上叠加就等于模拟了整个配电网一年、十年、一千年的运行剧本。每一年算一次年度指标最后把多年结果做统计平均。这个思路和拍电影很像解析法是直接给每个演员摆一个姿势拍一张照片序贯蒙特卡洛是把整场戏从头到尾演一遍再剪辑。前者省事但看不到过程后者贵但信息完整。2. 可靠性指标体系动手写代码前先把要算什么搞清楚2.1 三个最常用的持续性指标做配电网可靠性评估绕不开三个经典指标SAIFI、SAIDI、CAIDI。对应中文是系统平均停电频率、系统平均停电持续时间、用户平均停电持续时间。指标公式含义SAIFI年停电用户次数总和 / 总用户数平均每户每年停几次电SAIDI年停电用户时长总和 / 总用户数平均每户每年停多少小时电CAIDI年停电用户时长总和 / 年停电用户次数总和平均每次停电持续多久注意区分SAIDI和CAIDI的分母。SAIDI除以的是系统内全部用户包括这一年根本没有停过电的CAIDI除以的是真正经历过停电的用户次数它衡量的是一旦停电要熬多久。举个例子。假设一个简单配网有两个负荷点LP1接500户LP2接300户。一年内LP1停电3次共5小时LP2停电1次共2小时。那么SAIFI (500×3 300×1) / (500300) 2.25次/户·年SAIDI (500×5 300×2) / 800 3.875小时/户·年CAIDI SAIDI/SAIFI 1.72小时/次。写代码时公式本身不难难的是保证所有累计逻辑口径一致。2.2 电量类指标和负荷点指标SAIFI、SAIDI这类指标只关心用户数不关心损失了多少电。但在工程经济性分析里电量损失才是真金白银。于是有了ENS电量不足期望值和AENS平均系统电量不足值。ENS是把每次停电影响的负荷功率乘以停电时长再求和。它的单位是kWh或MWh直观反映一年因为停电少供了多少电。AENS ENS / 总用户数相当于把电量损失摊到每个用户头上。除了系统级指标负荷点指标也很重要主要三个λ是年停电频率U是年累计停电时间r是单次平均停电时长。它们之间存在 U λ × r 的关系。负荷点指标的价值在于定位薄弱环节——哪个分支的负荷点λ特别高说明这个位置受故障波及面大该考虑加装分段开关或者调整网架结构了。2.3 仿真统计口径逐年算再平均序贯蒙特卡洛输出的是一个年指标序列。假设仿真N年程序应该每年结束时统计当年的SAIFI、SAIDI、ENS存到数组里仿真结束后再求平均。这里有一个容易踩的坑有的实现是先把N年所有停电事件堆在一起最后统一除以N再除以总用户数数学上两者结果一致但后面做方差估计和收敛判断时没有逐年序列你就没法算年指标的标准差收敛性判断就会缺一条腿。3. 元件状态建模8760小时的事件流是怎么造出来的3.1 两状态马尔可夫模型和指数分布抽样配电网里几乎所有元件——架空线路、电缆、变压器、断路器——都可以用两状态模型描述运行状态和故障修复状态。运行一段时间后会故障故障后修复修复完继续运行如此循环。这个模型的关键参数是两个故障率λ次/年和修复率μ次/年或者等价地用MTTF平均无故障工作时间和MTTR平均修复时间来描述MTTF 1/λMTTR 1/μ。工程上最常用的假设是运行时间和修复时间都服从指数分布。指数分布有一个特性叫无记忆性意思是设备已经正常运行了1000小时它接下来继续运行的概率分布和刚投入运行时是完全一样的不会因为年纪大了就更爱出故障。直觉上有点反常识但在电力设备可靠性分析里对于长期平稳运行的元件这个假设有充分的统计依据是行业标准做法。既然是指数分布随机抽样就非常简单。逆变换法直接给出公式TTF -log(rand()) / lambda; % 下次故障前正常工作时间小时 TTR -log(rand()) / mu; % 故障后的修复时间小时rand()产生(0,1)均匀分布随机数-log(rand())就把均匀分布变换成了参数为1的指数分布抽样。除以λ和μ就是对应实际时间尺度的事件间隔。简单得让人怀疑但这就是蒙特卡洛的核心操作一切复杂的系统行为都是从这个最基本的抽样开始的。3.2 多设备事件如何在时间轴上交错单条设备的事件序列好生成从t0开始抽一个TTF时间推进到t1发生故障再抽一个TTR修复到t2再抽一个新的TTF……直到推进到8760小时一年为止。伪代码就是t 0; while t 8760 ttf -log(rand()) / lambda; t t ttf; if t 8760, break; end % 故障还没发生今年就结束了 ttr -log(rand()) / mu; % 记录故障发生在t修复时长为ttr t t ttr; end但一个配电网有很多条线路它们的故障事件是各自独立生成、又在时间轴上互相交错的。最常见的问题就是馈线F1在3月1日故障还没修复完馈线F2在3月5日又故障了。如果F2的下游负荷可以由F1所在的区域转供那F1的故障就会影响F2故障的处理策略。这就是为什么不能把每条线路的故障独立统计必须把所有元件的事件合并按发生时刻排序逐事件推进影响分析。我自己的做法是先把每条线路的事件序列全部生成好展开成一个大表——每行包含故障开始时间、故障设备编号、预计修复时长——然后按故障开始时间排序按顺序一个个处理。排序之前可以先检查有没有重叠这个检查和网络影响分析合并在一起做也行但分开写逻辑更清晰调试时更好定位问题。3.3 开关设备要不要建模线路和变压器是必然要建的两类元件开关类设备比较考验工程判断。断路器、隔离开关、熔断器本身的可靠性参数通常远高于线路故障率小一两个数量级很多研究直接忽略它们的故障只把它们当作影响传播的边界来处理这是合理的简化。但在某些场景下开关必须显式建模当研究重点就是开关拒动、保护误动时或者网架复杂、涉及多级配合时。这个决策没有绝对对错关键是写论文或报告时要把假设写清楚。我建议第一版程序先不建开关故障模型把断路器、隔离开关当成逻辑节点用等基础指标跑通再逐步加复杂度。一开始贪全往往会被调试成本拖垮。4. 故障影响分析整个程序里最耗时的部分4.1 辐射网的连通性判断配电网绝大多数是辐射状结构正常运行时呈树状每个负荷点从唯一路径取电。故障发生后通过开关操作重构网络部分负荷可能转由联络线路供电。影响分析的核心任务就是给定一个元件故障判断哪些负荷点停电停多久。最简单也最可靠的实现是基于图的连通性判断。把配电网抽象成图变压器和母线是节点线路是边负荷挂在节点上。某个线路故障就把它对应的边断开然后从变电站电源节点出发做遍历能到达的节点就是不停电的不能到达的就是受影响的。这样说起来简单实际操作有个经验要分享不要用递归深度优先搜索配电网规模不大还好节点数到了几百个递归可能爆栈。用队列做广度优先遍历或者直接用Matlab的graph对象和conncomp函数几行代码就能找到所有连通分量效率和稳定性都远好于自己手写递归。4.2 按故障位置的三种基本场景找到受影响负荷之后还要按故障位置和开关设置细化停电时间。我总结了三种最典型的情况场景处理方式停电时间构成故障点所在支路下游负荷、无联络转供隔离故障后下游负荷等故障修复修复时长TTR故障点所在支路下游负荷、有联络开关可转供隔离故障合上联络开关转供恢复开关操作时间手动约1~2小时自动约分钟级故障点上游、非本支路负荷保护跳闸-重合闸/备自投动作短时中断或完全不受影响注意同一个故障事件不同位置的负荷点停电时间可能完全不同。所以在程序内部故障影响的输出不能是一个笼统的停电了而应该细化到哪个负荷点停了多久这样才能分别累计SAIDI、ENS等指标。比如一条主馈线中段故障故障点下游的负荷可能要等修复停电时间是5小时故障点上游但和变电站之间有分段开关隔离的负荷可能通过合上联络开关20分钟就恢复而那些没有联络可用的分支末端就只能干等修复。同一个事件算出来的贡献可能是三种不同的时长这种精细度只有序贯蒙特卡洛能自然给出解析法很难表达。4.3 影响分析的性能优化影响分析函数会被调用极其频繁——仿真5000年平均每年100次故障事件就是50万次网络遍历。哪怕一次遍历只花1毫秒也要500秒叠加其他开销整个仿真可能要跑一两个小时。优化方向有几个第一把网络的邻接矩阵、开关位置、负荷点归属等结构信息一次性预计算好不要每次进函数重新解析输入表。第二优先判断故障对哪几个负荷点有影响而不是遍历全图。第三对常见的简单场景做快速通道如果网络不含联络开关且故障在末端支路影响范围就是该支路直接挂的负荷直接查表不用遍历。这个优化能把运行时间缩短一个数量级非常值得做。5. Matlab实现一套可复现的序贯蒙特卡洛仿真器骨架5.1 数据准备分支参数表是程序的宪法不管用什么语言写第一步永远是整理数据。我习惯用一张分支参数表作为输入核心每一行代表一条可故障的线路段。基本字段包括字段含义示例值BranchID支路编号1FromNode首端节点编号2ToNode末端节点编号3Length线路长度km0.8Lambda故障率次/年/整条线路0.08MTTR平均修复时间小时5Users该支路末端直接用户数120PeakLoad该负荷点峰值负荷kW350Lambda可以由每公里故障率 × 长度得到。架空裸导线的故障率一般在0.1~0.3次/年/km之间电缆低一个数量级具体取值参考所在地区的统计数据和相关标准不要自己拍脑袋。这一步做扎实后面所有结果的置信度才有基础。数据我习惯写成Excel或CSVMatlab里用readtable读进来程序主体不硬编码任何拓扑信息。这样换一个算例只需要换输入文件代码一行不动。5.2 主循环框架仿真主逻辑分三层最外层是年循环中间是事件推进最内层是影响分析。核心框架大致长这样rng(2024); % 固定随机种子结果可复现 simYears 2000; % 仿真年数 nBranch height(branchData); indicatorYear zeros(simYears, 3); % 逐年SAIFI SAIDI ENS for year 1:simYears eventList []; for br 1:nBranch lambda branchData.Lambda(br); mu 1 / branchData.MTTR(br); t 0; while t 8760 ttf -log(rand()) / lambda; t t ttf; if t 8760, break; end ttr -log(rand()) / mu; eventList [eventList; t, br, min(ttr, 8760 - t)]; t t ttr; end end eventList sortrows(eventList, 1); % 按故障开始时间排序 % 对eventList逐条调用影响分析函数累计SAIFI/SAIDI/ENS分子 % ... indicatorYear(year, :) ...; end meanValue mean(indicatorYear); % 指标期望值 stdValue std(indicatorYear); % 年指标标准差需要注意如果你的事件列表在生成时没有排除跨年重叠的情况也就是一个故障从去年延续到今年那么更严谨的处理是把这个状态作为今年的初始状态载入。很多教材里的简化做法是假设每年初所有设备都处于正常运行状态即把每一年当成独立同分布的重复这样程序大大简化指标结果也仍然是渐近无偏的但严格说每一年之间失去了连续性。在做论文仿真时我建议先说明采用的是独立年仿真还是连续年仿真两种口径两种口径下误差有细微差别审稿人有时会关心这个。5.3 影响分析函数与指标累计影响分析函数是核心它的输入是某个故障事件输出是每个负荷点是否停电、停多久。我的实现思路是这样的function [affectedLoad, outageHours] impactAnalysis(event, network, branchData) br event(2); % 1. 断开故障支路 % 2. 从电源节点做BFS找到连通区域 % 3. 对每个负荷点判断是否受影响按联络开关位置判断可转供恢复时间 % 4. 返回受影响负荷编号列表和各自的停电时长 end指标累计最容易出错的地方是用户数口径。SAIFI的分子是停电次数 × 该负荷点用户数的累加SAIDI的分子是停电小时数 × 用户数的累加ENS的分子是停电小时数 × 负荷功率的累加。三者累加的是不同的量必须分开维护三个累计变量或者一个结构体数组不要试图用一个变量同时算SAIDI和ENS功率和用户数的关系不是线性的混在一起必错。5.4 代码加速从能跑到跑得快Matlab跑蒙特卡洛最大的敌人是循环慢。三个提速手段是我每次都会用的第一用向量化代替内层循环。元件事件序列的生成可以改成先用rand生成一个大随机数矩阵再一次性做-log变换避免一个元件一个元件地循环。第二用parfor并行年循环。2000年的仿真如果你的机器有六个物理核理论上接近六倍加速。parfor year 1:simYears % 每一年内部的随机数生成互不依赖 end注意parfor要求循环体内不能有依赖前序迭代结果的共享变量所以每一年必须独立生成随机数、独立统计指标最后汇总。随机数种子也要小心处理默认情况下parfor的每个工作进程会使用不同的随机流结果可复现性需要额外设置建议先用普通for把结果调试正确再改parfor提速。第三用计时工具profiler定位瓶颈。有时候你想当然以为影响分析最慢实际上可能是事件列表用append拼出来的动态数组在反复扩张内存。改用预分配、或者用cell数组存储再一次性转换速度差异非常大。6. 收敛性判断跑了多少年才算数6.1 方差系数和相对误差怎么算蒙特卡洛仿真的结果本身是个随机变量仿真年数越多期望值的估计越准。但越准具体是多少需要量化。根据中心极限定理N年仿真得到的某个指标均值其标准误差是σ/√N其中σ是年指标的标准差。更直观的是相对误差η σ / (μ × √N)μ是年指标均值η表示估计值和真实值之间的相对偏差水平。一般要求η不超过5%关键指标最好控制在3%以内。这个公式直接决定了仿真年数怎么选。6.2 仿真年数N的经验取值具体需要多少年取决于系统的波动性。我曾经在一个小型配网上统计过SAIFI年指标的标准差和均值几乎同量级σ/μ接近0.8。要达到5%相对误差N至少是(0.8/0.05)^2 256年。看起来不多但如果系统里有低频重大事件、联络转供逻辑复杂σ/μ可能到2以上N就要上千。我在实际项目里的经验是先跑一个500年的试探性仿真输出逐年指标的均值和方差套公式估算所需N再正式跑。正式仿真N通常取1000到5000年之间具体看系统规模、硬件条件和精度需求。那些动不动说仿真十万年的如果不是系统特别高方差可能就是没用收敛判据单纯图个心理安慰。6.3 方差缩减技巧如果N上去之后运行时间不能接受就得在方差缩减上下功夫。我实测有效的方法有两个一个是对偶变量法。每次抽样生成U序列时同时用1-U生成一个镜像样本两个样本的方差比独立随机抽样小可以显著减少所需N。代价是程序复杂度增加。另一个是条件期望法。对故障率极低的元件不直接做随机抽样——反正它大概率整年不故障——而是把它对指标的期望贡献用解析公式先算出来从年循环中剥离。这样既能保留模型精度又消除了一部分随机波动。数学上这叫控制变量法实现难度适中效果非常明显。7. 实测心得跑通程序只是第一步7.1 先用标准算例验证程序正确性这是我最想强调的一点。配电网可靠性评估领域的公开测试系统比如经典的RBTS母线2和母线4文献里都有公布的可靠性指标基准值。第一次把程序写完之后先别急着往自己的电网数据上套拿公开算例跑一遍把SAIFI、SAIDI、ENS和文献基准值对比。偏差在合理范围内说明程序逻辑大概率是对的如果偏差明显赶紧回去查影响分析函数。这一步能帮你省掉大量无效调试时间。我之前有次程序跑出来的SAIFI偏大20%查了两天才发现是负荷点用户数统计重复计算——某个负荷点同时被两条支路记录。这种问题在实际拓扑短路的自己编的数据里很难发现但在标准算例里指标偏差一眼就能暴露问题。7.2 最容易踩的坑汇总问题表现原因与处理只统计有故障的年份SAIFI、SAIDI虚高把无故障年份的指标按0计入序列而不是直接跳过事件时间取整导致状态跳变停电时长偏大偏小波动用浮点小时精确推进最后展示时再取整动态数组反复拼接程序极慢预分配事件列表或用cell批量处理随机种子未固定两次运行结果不一致入口处rng固定种子parfor场景单独处理联络开关操作时间忽略SAIDI偏大把开关操作时间作为独立参数和修复时间分开7.3 结果的呈现方式仿真结束后除了给均值一定要给标准差和相对误差。一组完整的可靠性评估结果至少包含两部分系统级指标表SAIFI、SAIDI、CAIDI、ASAI、ENS、AENS以及各负荷点或各馈线分支的明细指标表。明细表里的λ、U、r是定位薄弱环节的关键比如发现某条分支的U占整个系统SAIDI的40%那改造优先级就非常明确。参数敏感性分析也很有必要。把线路故障率从基准值上增减30%重新跑仿真观察SAIDI的变化幅度这能说明系统可靠性对哪些参数最敏感。这部分分析对配电网规划的实际价值往往比那几行平均指标大得多。7.4 向分布式电源和时序负荷扩展最后提一个扩展方向。序贯蒙特卡洛真正强大的地方在于它天然适应时序变化——光伏出力白天高晚上零储能削峰填谷电动汽车负荷晚高峰这些特性都是围绕时间轴展开的。解析法要把这些因素塞进去极其困难但序贯蒙特卡洛只需要在每年步进时把该小时的分布式电源出力系数、负荷水平系数乘进去ENS的计算就自然考虑了时序特性。也就是说你花力气搭好的这套仿真器并不会止步于基础可靠性指标它可以平滑过渡到高比例分布式电源接入场景下的可靠性评估这也正是当前配电网规划研究最热的方向之一。最后说一个我自己的习惯。搭这种仿真类的程序我把最容易出错的三个地方——事件序列生成的随机数、影响分析中的连通性判断、指标累计的口径——单独拆成三个函数每个函数都配一个独立的小测试用例。主程序跑之前先跑这三个测试。这三道防线守住之后不管后面怎么改数据、改参数我都能保证核心逻辑是稳的。配电网可靠性评估这项工作最终的可靠程度其实不是取决于累计了多少随机数样本而是取决于你对每个样本背后物理含义的理解有多准确。把那一层想透了仿真年数、代码写法、指标口径都变成自然而然的选择了。