序贯蒙特卡洛模拟法做配电网可靠性评估这个项目我前后折腾了快三个月。如果你也是电力系统方向的学生或者刚进入配电网规划岗位八成会遇到类似的需求领导让评估现有馈线的供电可靠性或者做改造方案的横向对比需要一套能算SAIFI、SAIDI、ENS这些指标的工具。Matlab是首选实现语言因为它的矩阵操作、随机数生成、绘图都比其他脚本语言方便而且学术界、工业界都有大量现成代码可参考。这篇博文就围绕“基于可靠性评估序贯蒙特卡洛模拟法的配电网可靠性评估研究Matlab代码实现”来写把原理、代码架构、实操流程和踩坑经验一次说清楚。1. 项目定位配电网可靠性评估到底在“评估”什么1.1 先搞清楚指标SAIFI、SAIDI、CAIDI、ASAI、ENS配电网可靠性评估的核心是回答一个问题在给定网络结构和运行方式下用户一年到底要停几次电、停多久。工程上很少用一个指标概括全部信息而是用一组国际通用的可靠性指标来刻画。SAIFI系统平均停电频率单位是“次/用户·年”计算公式是系统用户停电总次数除以用户总数。这个指标反映停电有多频繁和故障率直接挂钩。SAIDI系统平均停电持续时间单位是“小时/用户·年”计算公式是系统用户停电总时长除以用户总数。它反映停电有多久受修复时间、隔离时间、转供时间共同影响。CAIDI用户平均停电持续时间等于SAIDI除以SAIFI单位是“小时/次”。它衡量一次停电平均持续多久是供电企业服务水平的直观体现。ASAI供电可用率等于1减去SAIDI除以8760是一个接近1的小数比如0.9991。它表示用户一年里可用电时间所占比例。ENS期望缺供电量单位是“MWh/年”计算各类停电场景下缺少的电量总和。这个指标直接和经济损失挂钩规划人员最看重它。这五个指标并不是孤立的它们之间满足“频率×每次时长总时长”的逻辑关系。写报告的时候通常把这五个指标全部列出来因为不同场景下决策视角不同——SAIFI低但SAIDI高说明停电虽然不频繁但一次停很久问题出在检修或转供能力不足反过来则是故障太多但处理很快。1.2 为什么配电网的评估比输电网更麻烦输电网通常是环形网络局部故障可以通过多重联络自动转带评估时用确定性N-1准则往往就够了。配电网不一样绝大多数处于辐射状或弱环状运行馈线之间依靠联络开关连接故障隔离和恢复供电压根结底是一系列复杂的开关动作逻辑。一个10kV馈线上通常有几台分段开关、联络开关还可能有柱上断路器、熔断器。某段线路发生故障后变电站出线断路器先跳闸整条馈线全部失去电源随后调度或自动化系统定位故障通过拉开故障点两侧的隔离开关把故障段隔离再合上联络开关把非故障段转带到其他馈线。这一系列动作的时间和成功率直接决定用户停电时长。解析法处理这种逻辑时公式会变得非常冗长而且一旦线路、开关数量增加状态空间以指数级膨胀手算公式根本不现实。蒙特卡洛模拟法的优势就在这里不追求穷举所有状态而是按照元件的故障规律随机抽取大量“故障—修复”事件用统计平均还原出系统的长期可靠性水平。这也是为什么在配电网可靠性评估领域蒙特卡洛法已经成为事实上的主流工具。1.3 序贯法和非序贯法怎么选蒙特卡洛模拟又分非序贯和序贯两类。非序贯法抽样的是系统状态——比如某时刻若干元件同时处于故障状态它对时间维度处理很粗糙无法表达“先故障后恢复再转供”这类有时间先后逻辑的过程。序贯法全称是序贯蒙特卡洛模拟法它把时间轴真实推进先抽样每个元件的无故障工作时间TTF和故障修复时间TTR再在时间轴上合并所有元件事件逐段判断系统状态。正因为它是按时间走的所以能自然地处理时变负荷、分布式光伏出力的时间相关性问题也能模拟一天内不同时段转供策略差别带来的影响。代价是计算量更大程序也更复杂。如果只是搞个简单辐射网的教学演示非序贯法也能凑合一旦涉及实际馈线、DG、储能序贯法几乎是必选项。这个项目选序贯蒙特卡洛模拟法方向是对的。2. 序贯蒙特卡洛模拟法的核心原理用一个简单例子讲透2.1 元件状态持续时间抽样指数分布与随机数序贯法的底层逻辑是“随机抽样元件状态持续时间”。架空线路、变压器、断路器这些元件在可靠性建模中最常见的是两状态模型正常运行和故障停运。正常运行时间TTF和故障修复时间TTR都是随机变量工程上一般假设它们服从指数分布。指数分布有一个非常好的性质累计分布函数F(t)1-e^(-λt)其中λ是故障率。给定[0,1]区间均匀随机数U令UF(t)反解出t-ln(1-U)/λ。由于U和1-U同分布可以简化为t-ln(U)/λ。这就是蒙特卡洛抽样状态持续时间的核心公式。在Matlab里写出来就两句function ttf sampleTTF(lambdaPerYear) % lambdaPerYear: 元件年故障率单位 次/年 % 返回正常运行时间单位 小时 ttf -log(rand()) / lambdaPerYear * 8760; end function ttr sampleTTR(mttr) % mttr: 平均修复时间单位 小时 % 返回本次修复时间单位 小时 ttr -log(rand()) * mttr; end注意第一段代码里除以“次/年”得到的是“年”乘8760换算成小时。如果直接用小时故障率λhλ/8760那么TTF-ln(U)/λh效果一样只是单位要时刻盯紧这是我最早最容易出错的地方。2.2 把多个元件的状态序列在时间轴上合并单台设备的抽样式子很简单但系统里有几十条线路、多台变压器和开关怎么把它们合成一个全局的过程序贯法是这么做的对每个元件先抽样一次TTF记录下来把它当作该元件的首个事件。找到所有元件事件中时间最早的那一个把系统时钟推进到那个时刻处理该事件比如某条线路故障。被处理的元件进入修复阶段抽样一个TTR修复结束后再继续下一次TTF。其他未发生事件的元件状态保持不变。重复“找最早事件→推进时间→执行故障后果分析→安排下一次事件”这一循环直到模拟的总时间达到预设年限比如10000年。这个方法实现时通常会维护一个事件列表每一行记录“元件编号、事件类型故障/修复、发生时刻”。每次从列表里取最小时刻的事件处理完后用新抽样的事件替换旧事件然后重新排序。Matlab的排序指令sortrows很好用但如果模拟年限长、事件多循环时间会相当可观后面我会说怎么优化。2.3 故障后果分析是精确性的关键系统级模拟只负责生成“何时发生故障”真正决定可靠性指标的是故障后果分析。配电网故障后不同负荷点停电时间来源不一样故障点前端的负荷可能只停几分钟就被重合闸或隔离开关恢复故障点后端的负荷可能要等联络开关转供也可能一直停到故障修复结束。具体分析时要按馈线的拓扑结构做区域划分。以典型辐射状馈线为例某段线路故障后变电站出线断路器先跳闸全线失电。如果重合闸、备自投逻辑存在瞬时性故障可能在几十秒内恢复但在最高一级的长时间可靠性统计里这类短时停电通常也要计入SAIFI。通过隔离开关把故障段隔离后故障点前端的非故障段可以恢复送电停电时间等于开关操作时间。故障点后端的非故障段若存在联络开关且有足够备用容量则可以通过转供恢复停电时间等于转供操作时间如果没有转供路径只能等故障修复完毕停电时间等于TTR。故障区段内的负荷点只能等待修复完成整个过程停电。在实际代码里这种逻辑需要结合开关位置、电源点位置和联络线位置实现搜索算法。最常见的方法是邻接矩阵深度优先搜索先把网络按照开关设备分成若干个区段故障时只需要判断各区段与电源的连通性不必逐个负荷点去做拓扑搜索计算效率高很多。3. Matlab代码实现从数据结构到核心函数的完整拆解3.1 输入数据怎么组织才不容易乱配电网仿真最考验数据结构设计。我的建议是全部结构化存储不要用零散的矩阵变量否则传到后面自己都看不懂。我常用四个结构体% 节点表 node.id % 节点编号 node.type % 类型1-变电站电源, 2-负荷节点, 3-普通连接节点 node.loadMW % 平均负荷单位 MW node.userNum % 该节点折算的用户数 % 支路表线路、变压器统一处理 branch.id % 支路编号 branch.fromNode % 起始节点 branch.toNode % 终止节点 branch.lenKm % 长度单位 km branch.lambdaPerKm % 每公里年故障率 branch.mttr % 平均修复时间单位 小时 branch.isSwitch % 是否含开关 branch.switchType % 开关类型0-无, 1-断路器, 2-隔离开关, 3-联络开关 branch.openTime % 开关操作时间单位 小时 % 联络开关表用于转供路径搜索 tieSwitch.id tieSwitch.branchId tieSwitch.nodeA tieSwitch.nodeB % 模拟参数 simParam.years % 模拟年限 simParam.seed % 随机数种子 simParam.convergeCriterion % 收敛判据阈值把线路、变压器统一到同一张表里对代码实现是最方便的因为它们的可靠性模型本质上相同只是故障率、修复时间参数不同。在实际项目中我需要把GIS系统导出的cad/dwg表格整理成这种结构过程比较繁琐但一次整理清楚后面跑仿真就很顺。3.2 事件表驱动的主循环怎么做用一个N行三列的矩阵events来维护当前的事件表每一行格式是[元件编号, 事件发生时刻(小时), 事件类型]。事件类型用1表示故障开始0表示修复完成。初始化时给每个元件生成一次TTF作为它的第一个故障事件后续当某个元件被处理完就补上它的下一次TTF或TTR事件。主循环大致是这个样子function [metrics] runSequentialSimulation(node, branch, tieSwitch, simParam) rng(simParam.seed); nBranch length(branch); % 初始化事件表 events zeros(nBranch, 3); for k 1:nBranch ttfHours sampleTTF(branch(k).lambdaPerKm * branch(k).lenKm); events(k,:) [k, ttfHours, 1]; % 第一个事件都是故障 end % 累积统计变量 totalStopNum 0; % 总停电用户次数 totalStopHour 0; % 总停电用户小时数 totalLackEnergy 0; % 总缺供电量 MWh totalYears 0; % 已经模拟的年限 while totalYears simParam.years % 取时间最早的事件 [minTime, idx] min(events(:,2)); k events(idx,1); evType events(idx,3); % 更新总模拟时间跨年时统计年度累计值 ... if evType 1 % 故障开始 faultBranch k; [impactedLoads] analyzeFault(node, branch, tieSwitch, faultBranch); % 根据影响的负荷点累加停电次数、停电时间、缺供电量 ... % 下一次事件修复完成 ttrHours sampleTTR(branch(k).mttr); events(idx,:) [k, minTime ttrHours, 0]; else % 修复完成 % 下一次事件下一次故障 ttfHours sampleTTF(branch(k).lambdaPerKm * branch(k).lenKm); events(idx,:) [k, minTime ttfHours, 1]; end end end这个框架简单但实用胜在逻辑清晰。需要注意一个容易忽略的细节模拟到模拟年限的边界时如果跨越了整年分界点需要把前面已经累加的指标先结算一次然后再把总模拟时间清零重新累计。很多新手会漏掉这一步导致每年统计的天平被截断结果偏差很大。3.3 故障后果分析函数怎么搜拓扑故障后果分析是整个程序里最考验逻辑的部分。我采用“区段划分广度优先搜索”的策略具体步骤如下把所有开关包括断路器、隔离开关、联络开关当作天然的分断面网络被它们切分成若干个区段。找到故障区段假设是segment_fault。从变电站电源节点出发做广度优先搜索遇到断开的开关停止扩展。这能得到故障隔离前有电的区段。故障隔离后把故障段两侧最近的开关拉开。此时电源侧的区段都恢复供电非故障侧区段如果通过联络开关能连接到其他电源则标记为转供区段。根据区段类型逐一统计每个负荷点的停电频率和停电时长。这条逻辑里最容易出错的是联络开关转供时还要校验备用容量是否充足。如果另一回馈线自身负载率已经很高转供可能导致过载工程上不能简单认为“一连就通”。在简化模型里可以设置一个可转供容量上限超过则不允许转供此时下游负荷只能等待故障修复。用一个Matlab函数封装function [stopFreq, stopDur, lackEnergy] analyzeFault(...) % 输入故障支路编号、网络拓扑、开关状态 % 输出每个负荷点的停电频率、停电时间、缺供电量 end这个函数的返回值需要设计精确因为后续所有可靠性指标都建立在它的基础上。3.4 指标统计和收敛判据模拟结束后把所有年份累计的停电次数、停电时间、缺供电量分别除以总年份数或者用户数就得到最终的可靠性指标。但这里有个关键问题到底模拟多少年才够蒙特卡洛法天生带有随机误差误差大小用方差系数β来衡量βσ/(μ√N)其中σ是样本标准偏差μ是样本均值N是模拟年限。工程上通常要求β≤0.05对高精度场景要求0.01。我的建议是在主循环里每模拟100年就检查一次β一旦满足收敛条件就提前终止而不是傻乎乎地跑固定的年份。这样能在保证精度的前提下省大量计算时间。beta sigma / (mu * sqrt(N)); if beta 0.05 break; end注意SAIFI、SAIDI、ENS各指标对应的β不同收敛速度也不一样。实际项目中一般只对最关心的那个指标做收敛检查通常是SAIDI或ENS因为它们波动更大。4. 实操复盘一个三馈线测试系统的完整评估4.1 测试系统怎么搭为了验证代码我用一个简化的10kV测试系统三回馈线每回馈线带6个负荷点。三回馈线从同一座110kV/10kV变电站的不同10kV母线引出馈线之间末端通过联络开关两两相连。线路长度约3公里采用架空裸导线故障率取0.1次/(km·年)平均修复时间4小时。变电站出线断路器自动跳闸分段隔离开关操作时间0.5小时联络开关转供时间1小时。每个负荷点用户数从100户到300户不等平均负荷0.2~0.5MW。用上面的参数模拟20000年并设随机数种子为固定值比如2026保证结果可复现。实际运算时间在普通笔记本上大概跑几十秒到几分钟完全在可接受范围。4.2 从参数设置到结果解读参数文件直接写成Matlab脚本最方便我通常单独建一个setup_test_system.m内容和上面3.1节的结构体对应。运行主程序后输出的指标大致如下指标数值单位SAIFI1.32次/用户·年SAIDI7.65小时/用户·年CAIDI5.80小时/次ASAI99.913%—ENS238.4MWh/年这个量级符合10kV架空配电网的实际情况。如果仿真里去掉联络开关转供只靠馈线出线断路器和隔离开关SAIDI会显著上升因为故障点下游负荷只能等修复时间长达4小时增加了联络转供后负荷点平均停电时间降下来不少。这也是为什么许多配电网改造项目的核心就是增加分段开关和联络线。4.3 灵敏度分析让结论更可信项目报告里只放一套结果远远不够通常还要做灵敏度分析。比如修复时间从4小时变成8小时SAIDI和ENS几乎翻倍但SAIFI不变把联络开关操作时间从1小时压到0.5小时CAIDI小幅下降说明这个环节不是瓶颈如果把线路故障率从0.1提到0.2次/(km·年)SAIFI、SAIDI、ENS整体都会上升约一倍。这类分析对规划决策非常有用。因为它能告诉你同样投资是换电缆降低故障率更划算还是增加联络开关提升转供能力更划算。你可以在代码里把这些场景写成批量循环一次跑完自动生成对比表。5. 常见问题与排查技巧实录5.1 结果始终不收敛到底是模拟年数不够还是代码有bug我遇到过最典型的“假不收敛”模拟了5、6万年SAIDI还在波动。排查后发现是随机数种子的影响——换了种子结果差异很大这说明样本量确实不够或者方差太大。但另一个更隐蔽的问题是对数计算时lambda传了0值导致TTF变成无穷大事件表里出现不合理的极大值整个统计被污染。建议从两个方向排查第一检查事件表里是否有明显不合理的时间值比如几百万小时第二削减到单馈线小网络跑一个手工可以算的案例验证主逻辑没问题后再扩大规模。5.2 可靠性指标怎么验证写程序最大的风险是“一本正经地算出错误结果”。我验证程序可靠性的方法是和解析法对比一个小型系统比如一个无限大电源带一条辐射馈线、三个负荷点没有开关转供。这个简单网络的SAIFI和SAIDI理论上等于线路总故障率乘一些修正因子手算小网络能算出来代码跑出来的结果只要跟手算量级和趋势一致就说明整体建模逻辑和FMEA正确。5.3 Matlab性能优化从“一天跑不完”到“半小时跑完”序贯蒙特卡洛模拟天生慢但在Matlab里优化空间很大。第一个方法是事件表用预分配内存不要每循环一次就动态增长。第二个方法是把大量样本一次性向量化抽样批量生成TTF、TTR再把事件排序避免在for循环里反复调用rand()。第三个方法是故障后果分析用稀疏矩阵和图形对象加速不要用嵌套for循环逐节点搜索。我印象最深的是把事件推进逻辑从“逐事件循环”改成“按故障区间批量推进”后计算速度提升了近一个数量级。代价是代码复杂度上升对于教学和一般工程项目先保持简单清晰的版本在确实需要提速时再优化。5.4 负荷点用户数权重容易算错SAIFI和SAIDI都是用户数加权指标。很多初学者直接把每个负荷点的停电次数取平均忘记乘以用户数导致结果整体偏小或偏大。我在代码里会专门维护一个userNum数组统计时严格执行“用户数×停电次数”的汇总再除以总用户数每一步都单独输出检查。6. 扩展方向DG接入、储能和配电网新型态6.1 分布式光伏接入后的时序效应传统可靠性评估假设负荷恒定真实世界的10kV馈线越来越多地接入分布式光伏。光伏出力有强烈的昼夜和季节特征白天可能降低馈线负载从而改善电压、增大转供能力但夜晚又完全无出力。序贯蒙特卡洛模拟法天然适合处理这种时序问题把时钟推进到每个小时光伏出力按实际日照曲线抽样负荷按季节典型曲线变化这样评估出的可靠性指标更贴近现实。6.2 微电网孤岛运行的自愈能力配电网加装储能和微电网控制后外部故障时局部可以离网运行即“孤岛模式”。序贯法的时序模拟可以精确刻画储能SOC在故障发生前的状态进而判断孤岛能维持多久。这在非序贯法里是做不到的。很多做分布式能源规划的朋友来找我咨询我一般建议直接从序贯法入手把光伏、储能按时间步长建模后续扩展很顺滑。6.3 从“评估”到“优化”启发式算法结合有了可靠的单次评估器就可以做配电网规划优化了。常见思路是用遗传算法、粒子群算法搜素网络拓扑或开关配置方案把序贯蒙特卡洛模拟得到的可靠性指标当作适应度函数。虽然计算量大但胜在灵活。我在实际项目里通常会把核心模拟函数封装成black-box接口优化算法每调用一次就完成一次可靠性评估整体循环可控。如果做这块建议先把模拟器性能优化到位否则优化算法跑一组种群就要几天项目根本不具备可行性。最后分享一点个人体会做这个项目踩过的坑很多但最大的体会是配电网可靠性评估的难处不在蒙特卡洛方法本身而在对配电网运行逻辑的理解。你得先清楚断路器怎么跳、隔离开关怎么拉、联络开关什么时候合然后才谈得上写程序。很多论文里的代码库问题不在于抽样公式错而在于故障后果分析过于理想化把转供时间、开关失败概率这些实际因素都省略了。真正能落地的东西恰恰是在这些细节里。如果打算自己写我建议先用一条馈线、两三个负荷点的小例子把主流程跑通再慢慢加开关、加联络、加DG。代码组织上数据结构设计要先想清楚事件表驱动的主循环和故障后果分析函数分开写这样调试起来方便很多。最后再把灵敏度分析和结果可视化补上整个项目就可以拿去应付报告、论文或者实际规划项目了。