做了几轮配电网可靠性评估的仿真程序之后我最大的感受是这个方向听起来门槛高实际跑通一次之后就会发现真正的难点不在数学公式而在“把电网运行逻辑转换成代码”的过程。这篇东西围绕一个用Matlab实现的“序贯蒙特卡洛模拟法配电网可靠性评估”项目展开把从原理到代码骨架、从指标统计到收敛判断的关键细节完整梳理一遍。适合电气工程相关专业的学生、做配电网规划或运行分析的研究人员以及刚接触电力系统可靠性计算、想搞清楚蒙特卡洛法到底怎么落地的人。1. 先弄清楚这个项目在做什么可靠性评估的定位与价值1.1 配电网可靠性评估与常见指标配电网直接面对终端用户绝大多数停电事故都发生在配电网这一层所以配电网可靠性评估一直是供电企业、规划设计部门和科研团队重点关注的方向。可靠性评估的核心目标很简单回答“每一年、每一个用户大概会停多少次电、停多久”。围绕这个目标行业里形成了一套标准指标体系。最常用的四个指标用生活化的方式理解就是SAIFI系统平均停电频率指标折算到每个用户头上一年平均停几次电。SAIDI系统平均停电持续时间指标折算到每个用户头上一年平均停电多少小时。CAIDI用户平均停电持续时间指标只要发生过停电的用户平均每次停电持续多久。ASAI平均供电可用率指标全年8760小时里用户能用上电的时间占比。SAIFI、SAIDI、ASAI是可靠性评估的“标配输出”几乎所有论文和报告都会给这三项。CAIDI则更偏运营层面的分析便于供电部门评估抢修速度和恢复效率。指标本身不复杂复杂的是怎么算出这些指标。对一个简单辐射网理论上可以用解析法手推但对一个包含分段开关、联络开关、备用电源、甚至分布式电源的配电网故障后的恢复过程千变万化想用纯公式硬解就非常吃力了。1.2 为什么选序贯蒙特卡洛而不是解析法先聊聊方法选型。配电网可靠性评估主流方法分为解析法和蒙特卡洛模拟法两大类。解析法的代表是故障模式后果分析法它把每个可能发生的故障比如某条馈线故障、某个变压器故障逐一枚举出来分析对下游负荷点的影响再概率加权求和。这种方法在拓扑简单、故障模式清晰的辐射网里很好用结果稳定、计算快。但它的短板也很明显故障场景组合爆炸、难以处理时序变化负荷曲线、DG出力曲线、光伏夜间不出力、难以表达“先隔离、再转供、后恢复”这种有时间先后的运行策略。蒙特卡洛模拟法走的是另一条路用大量随机抽样模拟系统的随机行为。它又分成非序贯和序贯两种。非序贯法只抽样系统的“状态”比如每条线路在某个时刻是正常还是故障然后判断系统是否失效原理简单、速度快但丢失了时间维度的信息无法自然表达故障持续过程、转供恢复顺序和时序负荷特性。序贯蒙特卡洛模拟法则把整个系统按时间轴一步步推进生成每个元件从正常运行到故障、再修复、再正常运行的完整时间序列严格遵循“时间不能倒流”的原则。它能把故障抢修、开关倒闸、负荷波动、DG出力变化全部装进同一个时间框架里这就是它最大的价值所在。项目名称里明确提到“基于可靠性评估序贯蒙特卡洛模拟法”本质就是要做一套能够真正还原配电网时序运行过程的可靠性仿真程序而不是停留在静态概率公式层面。Matlab则充当实现工具用来完成随机抽样、拓扑搜索、指标统计和结果输出。2. 序贯蒙特卡洛模拟法的运转逻辑从元件寿命到系统停电2.1 用指数分布给元件“算命”状态时序生成序贯蒙特卡洛模拟的第一步是为每个元件生成一段“寿命时间线”。元件在正常运行一段时间后发生故障持续一段时间后被修复再次恢复正常再故障、再修复……如此循环下去一直覆盖整个仿真周期。这里的关键假设是元件的故障发生过程可以用指数分布来描述。指数分布具有很强的无记忆性简单说就是“它不会因为已经正常运作了很久就更容易坏”这对电力元件来说是一个可接受的近似。随机抽样公式很简单正常运行持续时间抽样(d_{up} - \ln(U) / \lambda)故障修复持续时间抽样(d_{down} - \ln(U) / \mu)其中 (U) 是0到1之间均匀分布的随机数(\lambda) 是故障率次/年(\mu) 是修复率它等于平均修复时间 (r) 的倒数1/小时。注意单位要保持一致如果 (\lambda) 用“次/年”而仿真时间轴用“小时”那抽样出来的持续时间也要统一到小时。用生活化的类比来理解一个元件的平均寿命是已知的比如某条线路平均每10年出现一次故障但具体是第1年坏还是第12年坏谁也没法确定。指数分布抽样做的事就是“掷骰子”每次都给出一个具体的正常运行时间得到的长期统计平均值刚好等于 (\lambda) 的倒数。修复合集同理平均修复时间已知但每次具体修3小时还是8小时也是随机的。每一轮抽样得到一段“正常-故障-正常-故障”交替的序列把所有元件的序列放到同一条时间轴上就能判断系统在任意时刻是否完整供电。这比单纯算概率要直观得多。2.2 仿真钟怎么走逐小时步进与事件推进两套思路元件状态序列生成之后仿真钟怎么推进是实现上第一个分水岭。常见做法有两种。第一种是逐小时步进法把一年按8760小时离散在每个小时上检查所有元件的当前状态再判断每个负荷点是否失电。这种办法逻辑最简单适合刚上手的人但缺点也很明显重复判断次数多计算量随系统规模和仿真年限线性增长而且“小时”这种粗粒度有时会带来误差比如故障发生在上午10:15、修复在下午14:50逐小时步进法记录的停电时长就变成从10点到15点明显偏大。第二种是事件推进法也是我在实际项目中更推荐的做法。它不按固定步长推进而是把每个元件“状态发生变化”的时刻记录下来仿真钟直接跳到下一个事件发生时刻只处理该事件涉及的拓扑变化。这种方法精度更高因为状态切换的时间是连续抽样的不受步长约束速度也更快因为无需反复扫描整张网络。用事件推进法要注意一个数学上的好性质由于状态持续时间服从连续分布两个元件的状态从“正常”变成“故障”发生在同一瞬间的概率为零所以不需要处理“同时故障”的边界情况这在实现上省了很多麻烦。2.3 多重故障为何天然覆盖序贯法的隐藏优势非序贯模拟和解析法通常考虑单重故障N-1原则也就是假设同一时间只有一个元件故障。这个假设在规划阶段可以接受但真实运行中多重故障并非不可能发生一条线路检修期间另一条线路故障、极端天气下多条线路同时跳闸、继电保护拒动导致故障扩大这些情况都客观存在。序贯蒙特卡洛模拟天然就覆盖了多重故障场景。因为在时间轴推进过程中元件A可能正处于故障修复期元件B恰好也在这段时间发生故障系统就自动进入“N-2”状态。你不需要额外枚举多重故障组合只需要把所有元件的状态时间序列叠加到一起系统自然就会“撞上”这些场景。我在实际调试中就遇到过这种情况一个算例里SAIFI按单重故障解析法手算大概是某值但序贯模拟跑出来的数略高一点。排查后发现正是少量多重故障场景贡献了额外的停电次数。这不是程序bug反而是模拟法更接近真实的体现。3. Matlab代码核心架构数据结构决定实现上限3.1 把一张配电网络变成矩阵和表格写代码之前必须先回答一个问题配电网在Matlab里怎么表示配电网的一个典型特征是辐射状结构由母线、馈线、分段开关、联络开关、配电变压器和负荷组成。在代码里我习惯用“节点-支路”模型来描述节点编号每个母线或负荷点一个编号。支路描述每条馈线段对应一条支路包含起始节点、终止节点、故障率、平均修复时间、是否装有分段开关等信息。电源节点主网变电站出口或分布式电源接入点。在Matlab中节点和支路数据可以存放在结构体数组或表中。对于一个中等规模系统几十个节点、几十条馈线段用结构体数组已经绰绰有余。为了便于拓扑搜索我还会构造一个邻接矩阵或稀疏邻接矩阵用于快速判断任意两个节点之间是否存在可用路径。3.2 元件参数、负荷与开关信息的组织方式可靠性的核心输入参数有三个故障率、平均修复时间、开关操作时间。这三个参数直接决定了仿真结果但它们的含义经常被混淆。故障率 (\lambda) 通常以“次/年”为单位是元件在一段时间内发生故障的平均频率。平均修复时间 (r) 以“小时”为单位是元件发生故障后从发现到修复完成的平均时间。开关操作时间则是故障隔离和倒闸所需的时间通常只有零点几小时到几小时。此外还要定义负荷信息每个负荷节点的用户数或负荷功率。SAIFI和SAIDI的计算依赖用户数如果做电量不足期望ENS类指标则还要有负荷曲线。结构体设计大致如下对象关键字段作用节点id、类型电源/负荷/联络、用户数、负荷功率拓扑连接与指标统计支路from、to、lambda、r、开关类型故障抽样与隔离判定电源类型、容量、出力曲线主网与分布式电源建模运行策略联络开关操作时间、隔离逻辑转供恢复模拟这段数据结构设计是整个项目的根基。我见过不少第一次上手的人直接拿邻接矩阵硬写结果故障隔离逻辑根本没法表达——因为不知道哪些节点之间有开关、开关是常开还是常闭。先把数据模型理清楚后面代码就好写很多。3.3 函数划分六段式代码骨架这次项目的Matlab实现我按模块拆成了六段分别是数据初始化、时序状态生成、拓扑分析、故障影响判定、指标统计和结果输出。每一段单独封装成函数方便调试和复用。数据初始化模块读入网络参数构建结构体。时序状态生成模块对每个元件生成状态时间序列。拓扑分析模块给定当前时刻判断每个负荷点与电源之间是否有通路。故障影响判定模块确定失电负荷点、停电起始时刻、恢复时刻区分故障区段与非故障区段。指标统计模块汇总SAIFI、SAIDI、CAIDI、ASAI等指标。结果输出模块打印表格、绘制曲线、导出事件明细。主循环是整个程序的心脏建议写成“事件推进式”而不是“逐小时步进式”。把所有元件的下一次状态切换事件放进一个时间队列每次取最早事件处理更新网络状态再判断负荷点是否失电。仿真年限往大了设比如5000年甚至10000年目的是让随机指标收敛到稳定值。4. 关键模块实现拆解每一步都有代码可对照4.1 元件状态序列生成模块指数分布抽样这一块的核心是生成每个元件的“正常-故障-正常-故障”状态时间线。下面给出一个简化版的Matlab函数返回某个元件在仿真周期内的所有状态切换时刻和状态值。function [tEvent, stateSeq] genComponentSequence(lambda, r, T) % lambda: 故障率 (次/年) % r: 平均修复时间 (小时) % T: 仿真总时长 (小时) % tEvent: 状态切换时刻 % stateSeq: 切换后的状态1表示正常0表示故障 lambda_h lambda / 8760; % 转换为 次/小时 mu_h 1 / r; % 修复率 1/小时 tEvent 0; stateSeq 1; t 0; while t T % 抽样正常运行持续时间 d_up -log(rand) / lambda_h; t t d_up; if t T break; end tEvent [tEvent, t]; stateSeq [stateSeq, 0]; % 进入故障状态 % 抽样故障修复持续时间 d_down -log(rand) / mu_h; t t d_down; if t T break; end tEvent [tEvent, t]; stateSeq [stateSeq, 1]; % 恢复正常 end end三个容易踩的坑要提醒一下第一单位换算。(\lambda) 通常给的是“次/年”但仿真时用的是小时必须除以8760。我见过有人忘了这一步直接拿“次/年”当成“次/小时”用结果元件平均几小时就故障一次SAIFI直接爆表。第二rand函数的调用。这里用的rand是均匀随机数Matlab里也可以直接用exprnd但从实现透明性角度用-log(rand)反而更清楚。第三状态序列的边界处理。仿真结束时元件可能正处于故障状态这时候要决定是否把这段故障计为有效停电事件。一般做法是只统计完整发生在仿真周期内的停电事件或者把截断处的故障也计入但修正时长两种策略都要在文档里写清楚否则指标会出现微小偏差。4.2 拓扑连通性与停电事件判定有了元件状态时间线接下来就是在每个事件时刻判断网络拓扑是否完整。核心逻辑是从每个电源节点出发沿“状态为正常”的支路进行深度优先搜索能够到达的节点就是有电节点剩下的节点就是失电节点。function reachable dfsReach(adjMatrix, stateVector, sourceSet) % adjMatrix: 邻接矩阵adj(i,j)1表示节点i和节点j之间有支路 % stateVector: 支路状态向量1表示支路可用0表示支路故障 % sourceSet: 电源节点集合 % reachable: 逻辑向量表示每个节点是否可达 n size(adjMatrix, 1); reachable false(n, 1); stack sourceSet; while ~isempty(stack) node stack(end); stack(end) []; if reachable(node) continue; end reachable(node) true; for neighbor 1:n if adjMatrix(node, neighbor) 0 continue; end if ~stateVector(node, neighbor) continue; end if ~reachable(neighbor) stack(end1) neighbor; end end end end这段逻辑看着不复杂但它决定了整个仿真的正确性。实际操作中我会把支路状态向量做成二维邻接矩阵的形式而不是一维数组因为同一对节点之间可能存在双向支路一维数组容易漏掉方向信息。每处理完一次状态切换事件就对整张网络跑一遍连通性判断。如果某个负荷点从“有电”变成“失电”记录本段停电开始时刻如果从“失电”变成“有电”记录停电结束时刻把一条完整的停电事件写入事件列表。这里的经验是不要在每个小时都跑全网络搜索否则性能会很差。事件推进模式下一个仿真周期内通常只有几千到几万次状态切换每次搜索的复杂度是O(NE)加起来完全可接受。4.3 故障隔离与转供恢复逻辑这一块是配电网可靠性评估与输电系统可靠性评估最大的区别也是最容易写错的地方。故障发生后系统并不是“大家一起等到修好”。运行人员会先隔离故障点然后尝试通过联络开关把非故障区段转接到其他馈线或电源上。这就导致同一场故障中不同负荷点的停电持续时间完全不同故障点所在区段的负荷通常要等故障修复完毕才能恢复供电停电时长接近修复时间r。非故障但失电的区段如果存在可用的联络开关和备用容量停电时长只包含隔离和倒闸操作时间通常远小于修复时间。在代码里实现转供需要对停电事件进行二次细分。我的做法是检测到某负荷点失电后先判断它是否位于故障支路的下游且属于故障区段如果是停电时长取修复时长如果不是检查是否存在从另一个电源节点到该负荷点的“健康通路”若存在且联络开关可操作则停电时长取开关操作时间。% 伪代码逻辑 if loadNodeLost(node) if isInFaultSegment(node, faultBranch) outageDuration(node) repairTime(faultBranch); else if hasTransferPath(node, healthySources) outageDuration(node) switchingTime; else outageDuration(node) repairTime(faultBranch); end end end很多初版的可靠性评估程序把所有停电用户的停电时长都直接设为故障修复时间这在不含联络开关的单辐射网里是可以的但一旦网络里有联络开关SAIDI会明显高估。我在一个带单联络开关的算例里对比过不转供和转供两种策略下的SAIDI能差出30%到50%这个差距足以影响规划决策。转供逻辑还有两个进阶问题一是联络开关的容量约束如果备用电源容量不足转供时可能要甩掉部分负荷二是分布式电源孤岛运行需要判断主网失电时DG能否维持本岛负荷供需平衡。这些都要在事件处理函数里增加约束判断代码结构上留好扩展位。4.4 可靠性指标统计与输出仿真结束后手里已经有一份完整的停电事件表每一行是一次停电事件包含受影响负荷点、停电开始时刻、停电结束时刻、停电持续时间、受影响用户数。接下来的指标计算就水到渠成了。% event_duration_h: 每次停电事件持续时间小时 % event_cust: 每次停电事件影响的用户数 % total_customers: 系统总用户数 % years: 仿真年限 total_outage_cust_times sum(event_cust .* event_duration_h); total_outage_cust_count sum(event_cust); SAIFI total_outage_cust_count / (total_customers * years); SAIDI total_outage_cust_times / (total_customers * years); CAIDI SAIDI / SAIFI; ASAI 1 - SAIDI / (8760 * years);这段代码里最容易被忽略的是“乘years”。如果仿真跑了5000年累计停电用户时数非常巨大必须除以年限再除以总用户数才能得到“每年每户”的平均指标。我调试时有一次忘记除年限结果SAIDI莫名其妙大了三个数量级排查半天才发现是统计口径问题。结果输出部分我会把指标汇总成一张表同时把停电事件明细导出为CSV文件。事件明细的价值远大于聚合指标它可以用来做年际波动分析、故障原因追溯、重要用户停电统计甚至可以画停电持续时间的分布直方图这些是后续研究中最有价值的“原始素材”。5. 仿真结果与收敛性分析别急着相信第一次运行结果5.1 典型算例与结果解读用一个简单算例来演示某小型辐射配电网有两条馈线、四个主要分段、一个联络开关各分段故障率设为每年0.1次平均修复时间4小时联络开关倒闸时间0.5小时每个分段带200户用户。按解析法手算在不考虑转供的情况下系统SAIDI大约是各分段故障率乘修复时间之和即0.1×4×41.6小时。加入转供后非故障区段的停电时长从4小时缩到0.5小时整体SAIDI显著下降。用序贯蒙特卡洛模拟跑10000年输出结果大致是SAIFI约0.39次/户·年SAIDI约1.05小时/户·年CAIDI约2.69小时/次ASAI约99.988%这些数字在数量级上是合理的。SAIFI略低于分段数乘以单段故障率是因为部分故障发生时段上一段恰好也处于故障状态导致停电用户数偏少多重故障叠加下同一批用户被重复计入时逻辑上要合并处理SAIDI明显低于不转供的1.6小时正是联络开关的贡献。5.2 收敛判据与仿真年限选择蒙特卡洛模拟的本质是用随机抽样的平均值去逼近期望值样本量不够时结果会有波动。怎么判断仿真年限够不够工程上常用的判据是方差系数定义为标准差除以均值再除以样本量的平方根[ \beta \frac{\sigma}{\mu \sqrt{N}} ]对于系统级指标如果 (\beta) 小于0.05通常认为结果可以接受如果要求更高可以压到0.01。在实际项目中我一般先跑一个短年限测试比如500年看大概指标再逐渐加长年限直到关键指标的 (\beta) 满足要求。一个实用的建议把不同仿真年限的结果画成收敛曲线。你会看到SAIFI和SAIDI在前面一两千年波动比较大随着年限增长慢慢趋于平稳。如果5000年之后曲线还上下乱跳那就要检查随机数种子是否固定、是否存在极端小概率高影响事件比如某个年故障率只有0.01但修复时间很长的元件拉高了方差。5.3 结果分析的几个实用技巧第一留好随机数种子。仿真时固定rng种子可以保证结果可复现调试和修改代码时这个能力极其宝贵。最终发布数值结果时再换几个不同种子跑几组取平均作为最终答案。第二分开统计故障区段与非故障区段的停电事件。算SAIDI时把它们混在一起但分析时要拆开否则你没法判断到底是“故障抢修太慢”还是“转供策略太差”拖累了指标。第三利用事件明细做分布分析。只看平均值很容易被误导比如SAIDI是1小时但可能90%的年份停电接近0而某一年极端天气贡献了大量停电。对事件明细做年际统计能发现这类厚尾特征。6. 实测踩坑记录代码能跑只是开始6.1 随机数种子与结果复现蒙特卡洛模拟最大的特点就是“每次跑结果都不一样”。如果不在代码开头固定随机数种子同一个系统跑两次SAIFI可能差3%到5%。调试阶段这非常致命因为你无法判断指标变化是代码改动引起的还是纯随机波动引起的。我的习惯是函数入口设置一个可配置的rng种子参数。开发调试阶段固定种子正式仿真时循环更换种子跑多个批次。最终结果用多批次均值和标准差表达既能说明期望水平也能说明不确定性。6.2 停电事件合并的边界条件事件推进法下一个负荷点可能会经历“失电-恢复-再失电-再恢复”的反复。但如果两次失电之间只隔了几秒钟比如故障隔离后一次短暂试送电失败物理上用户可能根本没感觉到供电恢复统计上却会记成两次独立停电事件导致SAIFI虚高。解决办法是设定一个停电事件合并阈值比如两次停电间隔小于1小时且都由同一次故障引发就合并成一次事件。具体阈值可以根据工程经验调整但必须在文档里标注清楚否则换个人来用你的代码指标口径会不一致。6.3 性能优化事件推进才是正解我最早一版也是用逐小时步进写的逻辑简单但仿真5000年要跑很久。改成事件推进之后性能提升非常明显因为电网元件故障率通常很低绝大多数时间系统都在正常运行逐小时扫网络浪费了大量时间。事件推进的实现要点是维护一个“下一事件时间”的队列。Matlab虽然没有内置优先级队列但小规模系统直接用数组加排序也完全够用。我实测过一个几十个节点、几十条支路的系统事件推进跑10000年的耗时会降到逐小时步进法的几十分之一。还有一个小优化点连通性搜索时用稀疏邻接矩阵加逻辑索引替代密集矩阵的双层循环能再省出一截时间。6.4 面向分布式电源扩展的时序匹配问题如果项目后续要扩展到含分布式电源的场景有一个坑必须提前知道孤岛供电的时间窗口要与负荷曲线匹配。主网失电后DG可能带着一部分负荷孤岛运行但这个孤岛能维持多久取决于DG出力曲线和负荷曲线在每个小时谁高谁低。序贯蒙特卡洛之所以适合做这个场景就是因为它天然保留了时间轴可以把DG出力序列和负荷序列逐时段对齐。但要注意事件推进模式下网络状态是连续时间变化的而出力曲线往往是小时级离散数据两者对接时要统一时标否则会出现孤岛已失效但程序还在按DG供电的荒唐结果。最后再分享一个实用建议从最小的网络开始验证。我第一次跑通整套代码用的就是3个电源节点、4条支路、2个负荷点的模型指标手算都能算出来对照无误后再去跑几十个节点的中等系统。这个方法能帮你把“算法逻辑问题”和“统计波动问题”分开排查省下来的调试时间非常可观。我在实际项目中吃到最大的亏就是上来直接拿大系统跑结果一个隐蔽的转供容量判断bug混在大量随机事件里花了好几天才发现。现在回头去看凡是这类仿真程序小网络跑通闭环、手算对照合格再放大规模永远是最稳妥的路径。