做配电网可靠性评估的同行应该都有过这种纠结手头有故障率、修复时间这些基础参数想算一算系统全年平均停电多少小时、频率是多少却总在“用什么方法算”这一步卡住。公式推导不难难的是让模型贴近实际运行状态——比如修复时间有没有碰到负荷高峰开关动作是自动的还是人工去现场合的这些都会明显影响最终指标但常规的统计计算方法很难把它们收进同一个框架里。这个课题当时立项时核心目标很明确用序贯蒙特卡洛模拟法把配电网的全年运行状态逐时推演出来再通过Matlab把评估程序整体落地。项目做下来我对这套方法的体会是它不像解析法那样直接给一个公式解而是用“仿真年份叠加”的思路把时间、负荷、开关动作、故障修复这些因素全部揉进状态转移过程里。写这篇文章我会把整个研究过程的逻辑、程序实现的关键模块、以及我踩过的几个坑完整梳理一遍给想做配电网可靠性仿真、或者正在写蒙特卡洛相关代码的同学提供一个可以直接参考的工程化版本。1. 从解析公式到全年时序状态序贯蒙特卡洛到底弥补了什么1.1 传统故障枚举模型的盲区用一个简单例子说清楚很多人接触配电网可靠性分析最早学的都是故障模式及后果分析法FMEA和基于最小割集的解析法。这类方法思路很直接把每个元件可能发生的故障一一列出来通过网络的连通性判断哪些负荷点会失去供电然后累加得到负荷点的年故障率和年停电时间。这套方法在拓扑简单、元件数量少的系统里很好用但它的核心假设是“故障事件相互独立、修复过程不跟时间挂钩”。举个例子一条馈线在某年的运行过程中只发生一次故障但这一次故障是上午9点发生、维修2小时后恢复还是凌晨2点发生、维修8小时后恢复对用户的感受差异很大。前者可能赶上工厂开工损失几十度电后者多半只影响居民夜间照明负荷本身也处在低谷。可如果用解析法算SAIFI和SAIDI只是把“故障率”和“平均修复时间”相乘完全不体现这两个场景的区别。更现实的问题是现在配电网里大量接入分布式电源、储能和可转供负荷运行方式随时在变。解析法要处理多重故障、时变负荷、保护配合这些因素模型复杂度会指数上升最后落地的公式又长又难维护。1.2 “序贯”二字背后的三个关键能力序贯蒙特卡洛模拟法本质上是用“随机抽样状态持续时间推演”来模拟系统在一段长时间内的运行历史。每一年的仿真不是简单抽一个状态就算完而是把每个元件从正常运行到故障、再从故障到修复的整个过程按时间轴推进。具体到配电网我觉得它的价值主要体现在三点。第一支持时变负荷。仿真里可以把一年切成8760个小时每个小时的负荷水平按典型日负荷曲线取值。故障发生在哪个小时、修复占用哪些小时意味着损失的供能量ENS计算得更准不是简单用平均负荷蒙一个数。第二支持更真实的开关动作逻辑。配电网故障处理不是“抽掉一根线就完了”而是要先断开故障点附近的开关、隔离故障区段甚至合上联络开关把非故障区域从备用电源转供。不同开关有自动、手动之分动作时间从几秒到几十分钟不等。这些逻辑在解析法里很难用闭式公式表达但在序贯仿真里就是一段状态转移判断。第三支持评价“时间相关”的可靠性措施。比如配网自动化改造、储能应急供电、分布式电源孤岛运行这些方案它改变的不是元件的平均参数而是故障后的恢复路径和恢复速度。这类效果的差异只有把时间过程仿真出来才能体现。1.3 它和非序贯模拟、解析法的边界在哪里技术方案选型时经常有人问既然序贯蒙特卡洛这么好是不是所有评估场景都该用我的结论是取决于你要回答什么问题。方法核心原理优点典型适用场景解析法FMEA/最小割集枚举故障模式按逻辑关系叠加指标计算快、结果确定、可解释性强网络规模较小、结构简单、只需年均指标非序贯蒙特卡洛独立抽样系统各元件状态直接统计失效概率实现简单、适合求可靠性概率只关心概率型指标、不考虑时间和顺序序贯蒙特卡洛按时间轴推演状态持续时间统计故障频率与时长信息最全、可反映时间和控制逻辑需要评估指标分布、优化运行方式、分析复杂配电网我个人的经验是如果只是给某个新接入的项目算一个“年停电小时数”交差解析法完全够用没必要上仿真。但如果研究的是开关策略优化、分布式电源对供电可靠性的影响或者要做规划方案对比那非时序模型往往会忽视最关键的时间维度这时候序贯蒙特卡洛就是更稳妥的选择。2. 序贯蒙特卡洛模拟在配电网评估中的完整落地框架2.1 状态持续时间采样三套时钟如何驱动仿真年序贯仿真的第一步是把每个元件看成“运行—故障—修复—运行”循环的独立设备。工程上最常用的寿命模型是指数分布因为此时失效率是常数采样公式非常简便。对于第i个元件假设其年故障率为那么它从“正常运行”到“发生故障”的时间间隔TTFTime To Failure可以通过反变换法采样得到% 单位lambda_i 单位为 次/年结果 TTF 单位为 年 ttf_i -log(rand()) / lambda_i;如果希望仿真以小时为基本时间单位则把故障率换算成“次/小时”。比如一条线路年故障率是0.1次/年换算成小时为0.1/8760那么一次随机抽样得到的TTF就在逐年时间轴上表示该线路下次故障的等待时间。修复时间TTRTime To Repair同理如果修复时间服从均值为MTTR的指数分布% MTTR 单位为 小时 ttr_i -MTTR_i * log(rand());这里要特别注意单位协调TTF通常以年或小时记录TTR通常以小时记录。程序里的时间轴建议统一换算成小时否则在跨年事件处理时非常容易出边界错误。在配电网可靠性仿真里除了“故障-修复”这对时钟还要考虑计划检修。计划检修有一个固定的持续窗口大部分情况下检修期间允许提前转移负荷所以不能简单把计划检修当成一次停电事件来统计。我的做法是把计划检修的影响折算到开关操作判定里如果检修的元件本身已经有联络转供路径且切换动作可以提前完成那么对负荷点的停电贡献就接近零如果检修期间正好叠加了另外的故障那就会形成N-2类重叠事件这部分要单独记录。2.2 负荷曲线、开关自动化和转供路径如何进入状态推演仿真里的“一天”不是简单地存在“正常运行”和“故障停运”两个状态而是每一小时都有一组系统状态。为了不让程序膨胀成不可控的“小时级全状态枚举”我会在仿真循环内部用一个事件调度表来推进时间。基本流程是先给每个元件采样出首次故障时刻找到最早发生故障的元件把当前仿真时间跳到该时刻然后执行故障处理逻辑判断故障元件位置沿供电方向搜索最近的断路器或分段开关下发跳闸指令断开故障隔离区段检查与备用电源之间的联络开关是否具备合环条件如果负荷可以从对侧转供则合上联络开关恢复非故障区段故障元件修复完成后重新合上原分段开关恢复原始运行状态。其中负荷曲线是影响ENS指标的关键。仿真里可以用“典型日曲线x季节权重”的方式描述8760小时负荷也可以用更精细的96点曲线。每次故障事件发生后把故障持续时间切成整数小时分别乘上对应小时的负荷值再累加到年能量不足指标里。这个处理思路有一点像财务管理里的“逐日盯市”把每一笔停电时间精确落到它所影响的高峰还是低谷。2.3 辐射型配网下故障影响范围的判断拓扑关联分析配电网绝大多数是辐射型结构判断一次故障影响哪些负荷点本质上是一个“电源可达性”分析。我的程序里维护了一张基于邻接表的网络拓扑图每个节点记录上游供电电源、下游子节点、开关状态等信息。故障发生后影响范围分三层第一层故障点上游靠最近的保护开关断开从而被切除的区域第二层故障点下游失去主电源供电、但在联络开关合上后可以由备用电源供电的区域第三层故障点下游既无备用电源、也不能转供的末端区域这部分负荷只能等待故障修复。判定的算法很简单就是从变电站电源节点出发做广度优先搜索遇到断开的开关就停止扩展。每断一个开关或合一个联络开关之后重新进行一次连通性分析就能得到当前时刻的失电负荷集合。这个方案的优点是逻辑清晰、不容易漏判缺点是每次开关操作都要重新遍历网络。对于几百个节点的配电网每秒跑几千次遍历完全没问题。3. Matlab程序实现从网络参数到指标统计的完整拆解3.1 输入数据组织节点、支路、开关、负荷四张表写Matlab程序之前建议先设计好数据结构。我把整个配电网的参数拆成四张表全部用struct保存既方便后续扩展也方便做敏感性分析时批量修改参数。第一张是节点表记录节点编号、节点类型、所带用户数、峰值负荷、年用电量。第二张是支路表记录每条线路的起点、终点、长度、单位长度故障率、平均修复时间、是否有分段开关。第三张是开关表记录开关所在支路、开关类型断路器/负荷开关/联络开关、动作模式自动/手动、动作时间。第四张是联络线表记录可供转供的备用电源路径。% 参数结构示意 network.node.id (1:nNode); network.node.peakLoad peakLoad; % kW network.node.customers nCust; % 户数 network.branch.from branchFrom; network.branch.to branchTo; network.branch.lambda branchLambda; % 次/年 network.branch.MTTR branchMTTR; % 小时 network.switch.branchIdx swBranchIdx; network.switch.type swType; % 1-断路器 2-负荷开关 3-联络开关 network.switch.auto swAuto; % 1-自动 0-手动这里有一个容易被忽略的点同一段线路两端的故障率在不同供电方向下对负荷点的贡献是不同的。以某条支路故障为例故障点在线路首端和末端时被隔离开的负荷范围可能完全不同。我是用线路长度算单位长度故障率再按“故障点落在支路上的位置”做概率化处理而不是简单把整条线路当成一个点元件。这样累加出来的负荷点故障率更接近实际运行统计值。3.2 主仿真循环年推演、事件排序和状态更新主循环分成两层。外层是仿真年数内层是单个年里的故障事件推进。每一年的结果独立记录最后再对多年结果做平均和方差统计。这样设计的最大好处是可以利用多个仿真年的样本独立性来计算指标置信区间而不是只给一个均值就完事。伪代码如下% 核心逻辑示意 for simYear 1:maxYears % 初始化该年所有元件处于正常运行状态 for i 1:nBranch ttf(i) -log(rand) * YEAR_HOURS / branchLambda(i); end % 找出第一个发生故障的元件 [minTtf, idx] min(ttf); % 推进到故障时刻执行故障隔离、转供、统计 while minTtf YEAR_HOURS % 1. 时间推进 currentHour minTtf; % 2. 拓扑分析断开保护开关隔离故障区段 % 3. 判断联络开关是否可合恢复部分负荷 % 4. 累加负荷点停电次数和停电时长 % 5. 为故障元件采样修复时间更新 ttf % 6. 重新寻找下一个最小事件时刻 end end这个流程看起来简单真正复杂的是第4步和第6步之间的耦合。修复时间可能跨越到第二年需要把未完成的修复事件“结转”到下一年初处理同一年里还可能发生同一元件再次故障的叠加事件。因此我在代码里维护了一个“事件列表”每一轮只取时间上最早的事件处理完后再把由该事件引发的新事件如修复完成时间、重合闸动作时间插入列表。这种事件驱动的思路比“逐小时扫描”在计算效率上高很多。3.3 可靠性指标口径SAIFI、SAIDI、CAIDI、ASAI、ENS的统计方式指标计算是整个程序最后收口的环节。很多初学者在这里容易把“用电户数”和“负荷点数量”搞混。在配电网可靠性里一个负荷点可能带几十个用户指标里的用户数必须以计费用户数为准。我程序里的统计矩阵是“负荷点x年份”的停电次数矩阵和停电时间矩阵。每个负荷点记录下全年发生过的所有停电事件包括故障停电、计划检修停电、转供成功但经历短时停电的过程。指标计算公式如下SAIFI系统平均停电频率指标所有用户停电次数之和 / 总用户数单位次/(户·年)SAIDI系统平均停电持续时间指标所有用户停电时间之和 / 总用户数单位小时/(户·年)CAIDI用户平均停电持续时间指标SAIDI / SAIFI单位小时/次ASAI平均供电可用率指标(总用户数×8760 - 用户停电小时总数) / (总用户数×8760)×100%ENS系统年电量不足指标各负荷点停电时长×对应时段负荷功率的累加单位MWh/年Matlab里统计时要注意停电次数和停电时间的累加要在“事件发生时”同步完成不能在仿真结束后再事后补算。一个常见的问题是当仅仅发生“短时停电”而没有形成连续停电时SAIFI和SAIDI统计口径不同瞬时性故障引起的重合闸动作对SAIFI的贡献取决于是否计入瞬时停电通常瞬时故障不计入SAIFI只在单独统计“用户短时停电次数”里体现。我在程序里单独设置了一个标志位来控制是否累计瞬时性故障并把它作为输出结果里的附属指标。4. 收敛性控制、效率优化与可信度检验4.1 收敛判据跑多少年才算“算准了”蒙特卡洛仿真的随机波动不会因为样本量无限增加而消失但会随样本量增加而逐渐收敛到真值附近。所以要判断程序是否跑够了不能只看平均值是否变化很小要结合方差信息。常用的收敛判据有两个。一个是指标均值的相对变化率小于某个阈值比如连续几百个仿真年的SAIDI均值变化不超过0.5%另一个是变异系数CV计算公式为CV 标准差 / (均值 × 仿真年数)一般要求CV小于0.02到0.05工程精度已经足够。对中等规模配电网通常需要仿真几千到上万个年。如果系统里故障率特别低比如某条线路几十年才一次故障低频事件对SAIDI均值的影响会拖慢收敛速度这时候要考虑用方差削减技术。我在项目里的做法是每500年输出一次当前累计指标、标准差和CV并让程序在满足收敛条件后自动停止。不要一上来就固定跑10万年那样即耗时间又没必要也不要只看一次仿真的结果就下结论特别是做方案对比时两次仿真的随机波动可能超过方案本身带来的指标差异。4.2 提升收敛效率的三个工程化手段针对蒙特卡洛仿真速度慢的问题我试验过几种方法真正有效的有三个。第一个是拉丁超立方抽样LHS。普通的蒙特卡洛抽样是完全随机撒点可能会在某个区间内过度集中。LHS的思路是把[0,1]区间均等分成若干层在每个子区间里分别各抽一次保证样本在整个分布空间内分布更均匀。我把这个方法用在元件TTF采样上能在相同仿真年数下让SAIDI的方差下降大约20%~40%。第二个是公共随机数法。如果研究的目标是对比两种方案比如“有没有配网自动化”两种方案都应该使用同一组随机数序列这样两者差异主要来自方案本身的差异而不是来自随机噪声。Matlab里可以预先用RandStream保存一组随机流方案A和方案B共用这一组流结果对比会干净很多。第三个是并行计算。Matlab的parfor可以直接把独立仿真年分发到多个Worker上跑。要注意两点一是每个Worker之间必须使用独立的随机流否则并行后结果和串行结果会失真二是要控制进内存的数据量建议每个年只累积必要的矩阵不要存储整整几千年的逐小时状态快照。4.3 用解析结果反查仿真程序的对错程序写完后最重要的自检方法是先用解析公式快速算出几个核心指标再和仿真输出做对比。最简单的一个校验是用所有元件年故障率乘上负荷点平均用户数和平均修复时间粗略估计SAIFI和SAIDI的期望范围。如果仿真输出和这个粗糙的解析估值差了一个数量级那大概率不是收敛问题而是程序里有逻辑错误。常见错误包括单位时间轴没统一、线路故障率重复累加、转供成功时仍然把负荷统计为停电等。我自己的调试流程是先构造一个只有两三条馈线的小型测试系统手动枚举所有可能故障算出理论上的故障率和停电时长再让程序跑同一个系统。这个测试用例虽然简单却能覆盖切断开关、合联络开关、恢复供电的关键路径。只要小系统能对上换成大规模网络后起码可以说程序逻辑没有方向性错误。5. 我在Matlab实现中踩过的坑和最终经验5.1 随机数种子和可复现性看似小事坑起人来很要命蒙特卡洛程序最大的特点是随机性但这也意味着调试时如果每次结果都不一样你很难判断是代码改对了还是运气变好了。我强烈建议在程序入口固定随机数种子比如rng(2024, twister);这样每次跑出来的指标完全一致方便定位逻辑变化。做方案对比时再换不同的种子来检验结论的稳健性。这里还有一个细节如果程序里同时开了parfor固定主线程的随机数种子并不能保证每个Worker的抽样序列是可复现的。正确做法是用RandStream为每个Worker显式创建独立随机流并把这个流作为参数传入每个并行迭代里。我最初没有处理这个问题并行优化前后的结果始终对不上排查了很长时间才发现是并行随机流管理的问题。5.2 修复时间跨越仿真年的“结转”处理序贯仿真里一个元件的故障修复时间有可能跨越到下一个仿真年。例如第9999小时发生故障修复时间有5小时那么有4个小时的停电影响发生在第10000小时之后。如果程序在每个仿真年结束时就粗暴地将所有元件恢复为“正常运行”状态那这4个小时的停电就会被静默丢掉全年SAIDI和ENS都会偏低。我在程序里单独用了一个“跨年事件队列”在每年的最后时刻把所有尚未修复完的元件状态保存下来下一年初始化时直接从这个状态开始继续推演。这个改动看起来细节但对含低频长修复故障的输电网以及修复时间很长的大型变压器等设备影响会非常明显。5.3 瞬时性故障、重合闸和开关操作时间怎么建模配电线路故障并不全是永久性故障。架空线路上有相当比例的故障是瞬时性的比如树枝碰到导线重合闸动作后就能恢复。如果程序把所有故障都按永久故障处理SAIFI会被人为抬高这会严重误导“要不要上配网自动化”的决策。我的处理方式是每一条支路设置“瞬时故障占比”参数。每次采样出故障后先按概率判断是瞬时还是永久。瞬时故障只触发一次重合闸动作累加一次“用户短时停电次数”但不改变SAIDI和ENS永久故障才进入完整的隔离、转供、修复流程。开关操作时间也要细分断路器动作时间按毫秒到秒计通常认为不会引起负荷点停电自动分段开关的动作时间按几秒到几十秒计对整个年停电时间贡献很小手动操作开关需要运行人员到现场这部分时间通常在几分钟到半小时之间是SAIDI的大头。建模时要把这些时间分别累加统计口径才准确。5.4 针对大规模网络的代码优化心得配电网的规模一旦上百个节点逐小时扫描和频繁拓扑遍历就会成为运行瓶颈。我在程序里做了三个优化。第一把仿真时间轴的最小单位设为小时不设更小的分钟级模拟。开关动作时间本身只有分钟级计入小时级统计时向上取整对8760小时总时长的误差可以忽略。第二把“网络拓扑是否连通”的判断改成增量式更新。原本每次开关变位都要对整个网络做一次BFS优化后只保存每个节点当前所属的连通域编号开关变位时只做局部重标号代码速度提升了好几倍。第三把与小时负荷相关的统计计算向量化。比如故障持续时长是10小时与其在循环里逐小时做10次乘法不如先生成一块10x负荷点数的矩阵一次性累加Matlab的矩阵运算比循环快一个量级。做到最后我最大的体会是序贯蒙特卡洛模拟的难点不是数学公式而是工程建模的细致程度。它逼着你把“故障率”“修复时间”“开关动作”“负荷曲线”这些参数真正变成一套逻辑自洽的运行规则。同样的网络参数相同开关策略不同仿真出来的可靠性指标就会明显不同。这也是这个方法在配电网规划、运行方式优化上不可替代的原因。如果后面要继续往下做我觉得最有价值的方向是把序贯模拟和优化算法结合比如用仿真结果作为适应度函数去搜索分段开关的最优配置位置或者评估分布式电源不同渗透率方案下的可靠性提升效果。这个方向跑通之后这套Matlab程序就不再只是一个评估工具而是变成方案决策的实验台。