基于序贯蒙特卡洛的配电网可靠性评估:Matlab实现与避坑指南 做配电网可靠性评估最常被人问的一句话是这套网架一年到底会停几次电、平均停多久、折算到用户侧又是什么水平这些问题靠拍脑袋估不准尤其当网架里出现分段开关、联络线路、分布式电源以后故障影响关系变得特别绕。我前段时间刚好把一套基于序贯蒙特卡洛模拟法的配电网可靠性评估流程用Matlab完整跑通了从元件故障抽样、时序状态推演、负荷点停电记录到指标统计与收敛判断都做了落地。这篇文章把整体思路、代码结构和调试过程中踩过的坑整理出来给正在做可靠性方向课程设计、毕业论文或者规划方案比选的读者做个参照。1. 配电网可靠性评估到底在算什么从停电到指标1.1 停电是怎么发生和传导的配电网可靠性评估不是简单统计某台变压器坏了几次就结束了它关心的是“故障从元件出发经过开关动作与网络重构之后最终影响到了哪些负荷点”。这里面的因果关系链条往往比想象中长。举个最简单的中压辐射网例子一条馈线带着若干个负荷点沿线装了分段开关馈线出口有断路器。当某一段线路发生永久性故障时变电站出口断路器会先跳闸导致整条馈线失压——这是“故障隔离前”的停电影响。随后保护动作完毕分段开关通过人工或遥控方式打开将故障段隔离出来这时故障段上游的负荷点可以由变电站恢复供电故障段下游的负荷点如果存在联络线则可能倒换到相邻馈线继续供电如果没有联络就只能等待故障修复。这就产生了一个关键问题同一个故障事件不同位置的负荷点停电时间和停电次数可能完全不一样。有的用户只停“断路器跳闸到分段开关操作完成”的几分钟有的用户要停到故障修复完成有的用户还要算上联络开关倒闸的时间。可靠性评估要做的就是把每一次故障事件以及它引发的所有负荷点停录事件一一记录下来最后汇总成系统级指标。这也是为什么可靠性评估不能只看元件故障率而要和网络拓扑、开关配置、运行方式绑在一起分析。1.2 核心可靠性指标体系与公式行业里衡量配电网可靠性绕不开下面这张表里的几个经典指标。很多时候一个项目的最终交付就是给业主一张这样的指标汇总表。指标全称公式说明SAIFI系统平均停电频率指标用户停电总次数 / 总用户数平均每年每个用户停几次电SAIDI系统平均停电持续时间指标用户停电总时长 / 总用户数平均每年每个用户停多久CAIDI用户平均停电持续时间指标SAIDI / SAIFI每次停电平均持续多久ASAI平均供电可用率指标用户总供电小时数 / (用户数 × 8760)一年里用户能用上电的时间比例ENS期望缺供电量每次停电负荷 × 停电时长 的累加折算成损失的电量常用于经济性分析SAIFI和SAIDI是绝大多数报告里最显眼的两个数。ASAI则更直观——99.9%的可用率听起来很美好但折算一下就是每年约8.76小时的停电时间对用户来说体感一点也不轻。ENS经常被用来做停电损失评估把电量指标和电价、产损系数相乘就能算出一个网架方案对应的停电经济损失这也是可靠性评估经常被规划项目拉来做方案比选的原因。理解这些指标是理解序贯蒙特卡洛的第一步因为无论用什么算法最终输出都在朝这套指标体系收敛。2. 为什么是序贯蒙特卡洛三种评估路线的横评2.1 解析法靠谱但活儿只能干一半配电网可靠性评估的传统主力是解析法尤其以故障模式与后果分析法FMEA为代表。思路很直接枚举所有可能的元件故障事件对每个事件做网络影响分析再乘上故障率和修复时间叠加得到系统指标。解析法的优势是精确、可解释对简单辐射网和规模很小的配电网非常友好。手工就能列出几条馈线的所有故障模式再根据分段开关位置划出停电范围。但一旦网络规模变大或者引入多电源、分布式电源、储能和复杂的运行策略枚举组合会指数膨胀。比如一条馈线上出现两个元件同时故障或者一个故障叠加一次计划检修类似场景的排列组合数量直接起飞解析法要么算不动要么得做大量简化假设。另外还有更现实的一点解析法处理时变负荷非常别扭。传统FMEA通常假设负荷恒定或者只按峰值负荷计算这样ENS的精度会打折扣。当你需要评估“不同季节、不同时段的负荷水平对可靠性指标的影响”时解析法的公式会越写越复杂最后维护成本高到让人想放弃。2.2 非序贯蒙特卡洛与序贯蒙特卡洛差在一个“时间”蒙特卡洛模拟法给了另一条路不枚举而是抽样。由随机数驱动元件状态通过大量重复模拟来逼近真实可靠性水平。蒙特卡洛内部又分两派非序贯蒙特卡洛和序贯蒙特卡洛。非序贯的思路是“拍快照”随机抽样一年中某一天甚至某个小时判定该时刻系统状态抽样几万次后统计平均。计算效率高但它天然缺乏时间先后顺序难以精确处理需要时间逻辑的细节比如元件故障后修复持续多久、什么时候恢复供电、负荷随时间的变化曲线怎么走。很多情况下非序贯法只能给出期望值给不出时间相关的分布信息。序贯蒙特卡洛则是按照时间轴一步步推演。它给每个元件建立“运行-故障-修复-再运行”的状态时间序列所有元件并行推进从系统全年的时序响应里统计停电事件。这样负荷曲线能真实参与计算分布式电源的出力波动也能按小时建模检修计划、季节差异这类时间性很强的影响都可以塞进去。我所配电网项目里经常被问到“某个方案的年停电时间是多少”序贯蒙特卡洛给出的结果可以直接落到SAIDI、SAIDI这些指标上和实际运行统计口径更接近。2.3 什么时候该果断放弃序贯蒙特卡洛不过序贯蒙特卡洛不是万能药它有两个绕不开的痛点计算量大、收敛慢。因为要逐小时甚至逐分钟推进时序状态一个规模稍大的配电网仿真一万年模拟时间往往意味着上亿步运算每步还要做网络影响分析跑起来确实慢。不少人在项目里一上来就选序贯蒙特卡洛跑到一半发现收敛判据迟迟不达标才意识到前期估算不足。我的判断标准是这样如果网络比较简单、负荷模型可以取恒定值、又需要快速得到比较确定的指标直接用解析法。如果网络中等复杂度、只关心期望指标而不强调时序特性非序贯蒙特卡洛配合方差缩减技巧性价比很高。只有当场景里有明显的时序逻辑——比如时变负荷、储能充放电策略、分布式电源季节性出力、复杂的转供策略——序贯蒙特卡洛才有不可替代的价值。另外正式开跑前最好用非序贯或解析法粗算一轮给序贯仿真提供一个初始解和参考量级能提前发现参数错误避免白跑几十万步才发现模型根本不对。3. Matlab仿真主流程拆解从拓扑参数到指标输出3.1 元件可靠性建模与参数设定无论是哪种蒙特卡洛第一步都是把网络和参数变成程序能处理的数据结构。配电网可靠性仿真至少需要四类输入。第一类是拓扑信息。常见的建模方式是把网络抽象成节点和支路支路类型包括线路、变压器、分段开关、联络开关等。每个负荷点挂在某个节点上并关联用户数、平均负荷和时序负荷曲线。Matlab里我习惯用结构体数组存节点用邻接表表达连接关系这样后续做连通性判断和影响分析都方便。第二类是元件可靠性参数。线路和变压器需要给出故障率与平均修复时间。故障率单位一般用“次/(km·年)”注意乘以线路长度才是单条线路的年故障率。架空线路典型故障率在0.05~0.2次/(km·年)之间电缆会低不少修复时间通常在2~8小时之间。这些参数可以来自电网运行统计也可以参考经典文献的典型值。仿真结果的绝对数值完全由参数决定所以参数来源必须记录清楚。第三类是运行参数。比如分段开关操作时间、联络开关倒闸时间、计划检修周期、保护动作时间等。这些时间通常只有十几分钟到几小时但对SAIDI的影响不小不能忽略。第四类是负荷模型。如果只有全年平均负荷数据可以把负荷曲线设成单位比例曲线再乘以年峰值。更精细的做法是按季节、工作日/休息日、24小时分段给曲线。在Matlab里可以准备一个维数为[小时数, 负荷点数量]的矩阵模拟推到某一小时时直接查表得到该时刻负荷。3.2 序贯抽样的核心代码与状态推进序贯蒙特卡洛最核心的机制是对每个元件抽样它的无故障工作时间TTF和故障修复时间TTR然后让所有元件的状态沿着时间轴交替变化。抽样公式是可靠性工程里的老朋友。% lambda: 元件的年故障率 % MTTR: 平均修复时间单位小时 % rand: [0,1)均匀分布随机数 % 抽样无故障工作时间单位: 年 u1 rand; TTF -log(u1) / lambda; % 抽样修复时间单位: 小时 u2 rand; TTR MTTR * (-log(u2));这里用到的原理是逆变换法当随机变量服从指数分布时通过-ln(U)/λ就能从均匀随机数得到对应的采样值。元件“运行-故障”交替过程接近指数分布假设所以这套抽样公式在工程上是主流选择。主循环的写法有很多种核心思路都是“下一次事件驱动”。每轮循环找到所有元件里最先发生的状态转移把系统时钟推进到那个时刻然后更新发生转移的那个元件状态再做故障影响分析。Clock 0; % 系统仿真时钟单位年 SimYears 5000; % 仿真总年数 N length(lambda); % 元件数量 State ones(1, N); % 1运行, 0故障 NextChangeTime zeros(1, N); % 初始化: 给每个元件抽样首次状态转移时间 for k 1:N NextChangeTime(k) Clock (-log(rand) / lambda(k)); end while Clock SimYears % 找到最早发生转移的元件 [dt, idx] min(NextChangeTime - Clock); Clock NextChangeTime(idx); % 状态翻转: 运行-故障 或 故障-修复 State(idx) 1 - State(idx); if State(idx) 0 % 刚发生故障做故障影响分析 DoOutageAnalysis(idx, Clock); % 安排下一次修复完成时间 NextChangeTime(idx) Clock MTTR(idx) * (-log(rand)) / 8760; else % 刚修复完成安排下一次故障时间 NextChangeTime(idx) Clock (-log(rand) / lambda(idx)); end % 若有计划检修这里还需要插入检修状态及其持续时间 end这段代码里有个小坑时间单位很容易混。故障率λ单位如果是“次/年”-log(rand)/lambda得到的TTF单位就是“年”而修复时间一般按小时给如果MTTR的单位是小时就得除以8760换算成年否则系统时钟会乱掉。我见过不少人在这个单位换算上翻车改起来倒简单但定位起来特别费劲。3.3 故障影响分析谁在停电、停多久抽样完成之后最关键的环节是故障影响分析。这一步要做的事情是给定某个元件在当前网络状态下故障找出哪些负荷点失去供电以及各自要停多久。在辐射状配电网里这个分析可以按“上游、故障段、下游”三个区域来划分。首先找到故障元件所在的馈线以及变电站母线。故障发生后馈线出口断路器跳闸保护动作瞬间整条馈线失压。接下来隔离故障段分段开关打开后故障点上游且与母线连通的部分恢复供电恢复时间是开关操作时间。故障点下游的部分要看有没有联络开关能够倒供。如果存在联络开关且对侧电源容量足够带电合联络开关后下游恢复供电恢复时间为倒闸时间如果下游没有联络通路这部分负荷就只能等故障修复完成再恢复。注意实际项目中“容量足够”这四个字背后也要判断。联络馈线的剩余容量不够合闸后可能过载这时候要么削减部分负荷要么放弃转供。严谨的仿真应该在倒供判断里加入线路容量约束简单实现时可以查一回馈线额定容量与当前负荷之差。故障影响分析的结果要记录到按负荷点维度的数组里本次停电事件持续时间、停电时刻、停电期间缺供电量。由于序贯蒙特卡洛按时间轴推进每个负荷点在仿真年份里的所有停电区间都能完整记录天然支持后续指标统计和时间分布分析。3.4 指标统计与收敛判据仿真结束后指标统计就是纯算术。把仿真期间每个负荷点的停电次数、停电时长、缺供电量累加起来再按用户数加权就能得到表1里的各个指标。totalCustomers sum(CustomerNum); % 总用户数 SAIFI sum(InterruptionCount) / totalCustomers; SAIDI sum(InterruptionDuration) / totalCustomers; CAIDI SAIDI / SAIFI; ASAI 1 - SAIDI / 8760; ENS sum(EnergyNotSupplied); % 单位: kWh/年除了算指标蒙特卡洛仿真一定要做收敛判断否则仿真结束得很随意。这里引入方差系数β一般定义为目标指标估计值的标准差与均值之比。工程上常用的收敛标准是β小于0.02~0.05对应的仿真年数通常在几千到几万不等视网络规模和指标类型而定。% 以SAIFI为例计算逐年样本的方差系数 betaSAIFI std(SAIFI_year) / (mean(SAIFI_year) * sqrt(nYears));更稳妥的做法是边跑边检查实现“自适应停”每仿真200年检查一次β达到阈值就停止达不到就继续。这样比固定仿真年数更科学因为不同网络参数下达到相同精度所需的仿真年数可能差一个数量级。如果每次跑完发现不收敛优先考虑增加仿真年数或降低指标方差而不是盲目修改随机数种子。4. 实操中的坑与排查技巧4.1 仿真年数、方差系数与随机数种子的三角关系第一次跑序贯蒙特卡洛的人最容易上来就设一个“看起来很大”的仿真年数比如1000年然后默认数据够用了。实际上1000年对于结构比较复杂、故障率较低的网络SAIDI的方差系数可能还在0.1以上指标抖动非常明显。我的经验是先在Matlab里做一个“年数-指标”的敏感性测试。从1000年开始每次翻倍仿真年数观察SAIFI和SAIDI估计值的变化幅度。如果2000年的结果和4000年的结果差得不多几千年的量级基本够用如果差异还很大继续翻倍。这个测试本身也花不了多少时间但能避免两种尴尬情况一种是不够收敛就出报告另一种是不必要地堆仿真年数白白让程序跑一下午。随机数种子也需要刻意管理。MATLAB里每次启动rand的状态不同如果不设置全局种子同样是5000年仿真两次结果也会有可观察到的偏差。正规做法是在主程序入口固定种子rng(2024); % 固定随机数种子保证结果可复现这在调试阶段特别重要。当开发人员修改了某段影响分析逻辑后固定种子能保证两次运行只因为代码变化而出现差异排除随机噪声的干扰比对结果才有说服力。4.2 时变负荷处理恒定功率会算错什么很多初版程序为了省事把所有负荷点都设成全年恒定功率结果ENS算出来往往和真实情况偏差很大。原因很直观停电事件发生在凌晨低谷负荷时段损失的电量比发生在晚高峰时段小得多。恒定功率模型相当于把每个时段都按同一个负荷水平计必然扭曲缺供电量。要修这个问题我建议按三条曲线来做负荷模型工作日负荷曲线、周末负荷曲线、节假日负荷曲线。Matlab里可以预先生成一个全年8760小时的负荷比例矩阵每行对应一个小时的负荷倍数然后与负荷点的年峰值功率相乘。系统时钟推进到某小时时直接查表得到该小时每个负荷点的功率值。这个做法的代码量不大但对ENS和SAIDI的结果影响很明显。如果项目实在没有负荷曲线数据至少要分“峰、平、谷”三段负荷水平而不是一个常数走到底。我自己接手过一个评估项目最初的恒定负荷模型ENS比修正后高出约18%工程上这个差距足以改变方案排序了。4.3 联络开关与孤岛逻辑的边界判断另一类高频bug出现在联络转供逻辑里。最常见的错误是“无条件合联络”只要故障段下游有联络开关就认为负荷能全部转供不管联络馈线是否过载、对侧是否还有剩余容量。改进逻辑至少应该加三个判断。第一联络开关对侧的馈线必须处于运行状态不能那边也在检修。第二转供后联络馈线的载流量不能超过上限这个要结合当前时刻负荷算。第三转供操作时间要合理计入一般为几分钟到几十分钟不是瞬时完成。盲目合联络的结果是SAIDI被低估尤其重负荷时段故障时实际调度员根本不敢直接合闸程序却默认合上了。还有边界条件故障元件如果恰好在联络开关所在支路上或者两个联络开关同时可用影响分析逻辑要分清楚谁先合、合哪个不同顺序可能影响后续容量分配。遇到这种逻辑我建议先把情况拆成单故障事件表手算一遍期望结果再用代码跑同一事件对比快速定位逻辑错误。4.4 常见问题速查表现象可能原因排查思路指标明显偏小负荷恒定、忽略开关操作时间检查负荷模型和恢复时间设置SAIDI特别大联络转供逻辑错误或未启用检查故障影响分析中的下游转供判断仿真死活不收方差系数阈值过严、故障率过低增大仿真年数或改用方差缩减技巧两次运行结果不一致随机数种子未固定在入口处固定rng种子ENS异常偏大时间单位混用检查TTF/TTR单位是否统一为年孤立事件频繁出现联络线被认为不可用检查对侧馈线的运行状态判断5. 一个算例看指标怎么变成决策依据5.1 典型三馈线系统的仿真结果长什么样拿一个我实际做过的模拟项目举例。场景是某园区规划三回10kV馈线每回馈线带8个负荷点线路总长约18公里全部为架空线故障率取0.08次/(km·年)平均修复时间4小时。分段开关配置在每回馈线中段两回馈线之间有一处联络开关故障后可倒供部分负荷。所有负荷点年峰值负荷之和约5200kW全年负荷曲线按时段变化。用上面的Matlab流程固定种子跑5000年输出的核心指标大概是指标数值SAIFI1.34 次/(用户·年)SAIDI3.12 小时/(用户·年)CAIDI2.33 小时/次ASAI0.99964ENS约 8720 kWh/年这类绝对值本身并没有标准答案完全依赖参数设定和网架结构但它能反映出方案之间的相对差异。同一模型把联络开关取消掉以后SAIDI直接升到5.8小时以上ENS翻了一倍多这就是可靠性评估最有决策价值的时刻——多花多少钱加一个联络开关能换回多少停电损失一目了然。5.2 怎么从负荷点指标反查网架薄弱环节只看系统级指标不够规划里通常还要看每个负荷点的个体指标。我每次跑完仿真一定会输出一张负荷点指标清单然后按SAIDI从大到小排序。排在前面的负荷点要么处于馈线末端要么处于转供盲区要么紧邻高故障率线段。这类分析对规划很有帮助。比如某负荷点SAIDI明显高于系统平均值一看拓扑它处于故障段下游且没有联络通路所有上游故障都要等修复完成才能送电。对策可以是加装联络开关也可以调整线路分段位置来缩小故障隔离范围。可靠性评估从“算个数”变成“定位短板”靠的就是这些负荷点级别的中间结果。我实际做项目时还会把负荷点ENS按年汇总后再除以该负荷点的年用电量得到“损失电量占比”这个相对量对于比较不同性质的用户比如工业用户与商业用户更公平也方便进一步做停电经济损失评估。最后说几句体会做了这么多轮可靠性仿真我自己最大的感受是这套方法的难点其实不在蒙特卡洛抽样而在故障影响分析逻辑。抽样只是把随机性送进来判断哪些用户停、停多久、什么时候恢复才是决定指标准不准的核心。实际项目里我习惯把故障影响分析单独写成函数用一组手工构造的事件逐项测试确认无误后再接入蒙特卡洛主循环。另外再分享一个小技巧正式跑大批次仿真之前先固定随机数种子跑一轮记录下指标输出以后每次改代码都用同一批随机数做回归对比避免算法改进半天最后发现是随机波动在干扰判断。这套Matlab流程目前还在持续扩展后面计划把储能充放电策略和分布式电源出力的时序模型加进去进一步贴近真实的运行场景。