做配电网规划或供电可靠性分析的人大概率都绕不开一个问题某条馈线上接了哪些用户、线路多长、故障率多高一年下来用户平均停几次电、停多久、少供了多少电。工程上这个问题有一整套标准考核口径就是SAIFI、SAIDI、ENS这些可靠性指标。而当你面对的配电网越来越复杂时变负荷、联络转供、分布式电源这些因素全都掺进来的时候解析法会变得非常吃力这时候序贯蒙特卡洛模拟法就是最趁手的一把螺丝刀。这篇文章就围绕“用Matlab实现基于序贯蒙特卡洛模拟法的配电网可靠性评估”这件事把从建模、抽样、事件驱动模拟到指标统计的完整链路拆开讲一遍代码结构可以直接拿去改自己的算例适合正在做配电系统可靠性课题的研究生也适合电网公司或设计院里需要做供电能力评估和网架规划的工程师参考。1. 为什么选序贯蒙特卡洛方案选型的底层逻辑1.1 可靠性评估到底在算什么配电网可靠性评估的核心说白了就是回答三个问题每个负荷点一年大概停几次电、平均每次停多久、一年总共损失多少电量。围绕这三个答案行业里定义了一整套指标常用的有系统平均停电频率SAIFI、系统平均停电持续时间SAIDI、用户平均停电持续时间CAIDI、平均供电可用率ASAI以及缺供电量ENS。这些指标是电网公司内部考核、规划投资决策、供电合同约定里的硬数据。比如规划一条新的10kV馈线评审时一定会问这条馈线投运后系统SAIDI能改善多少下游新增用户对供电可靠性有多大的边际影响这些问题不靠拍脑袋必须拿出量化计算结果。计算这些指标本质上是一个概率分析问题。每个元件馈线段、配电变压器、断路器、隔离开关都有随机的故障和修复行为这些行为的随机组合决定了用户是否停电、停多久。可靠性评估的任务就是把这种随机性通过模型和算法翻译成期望值指标。1.2 解析法、非序贯与序贯蒙特卡洛的取舍行业内处理这个问题有三条技术路线我把它们的本质区别做个对比方法核心思路优势劣势解析法FMEA、最小割集、网络等值枚举所有故障模式按概率加权求和计算快、结果精确、适合简单辐射网系统规模大时枚举组合爆炸难处理时序因素非序贯蒙特卡洛按元件状态概率随机抽样一个系统快照能处理大规模网络、任意拓扑只“拍快照”不“拍电影”无法直接算停电频率和持续时间序贯蒙特卡洛按时间轴模拟每个元件的“正常-故障-修复-正常”全过程天然含时序能处理时变负荷、控制策略、复杂维修逻辑计算量大结果带统计波动解析法在单馈线、固定负荷、无控制策略的场景下非常好用手算公式都能出结果。但一旦系统里有时变负荷曲线、联络开关倒闸操作、分布式电源孤岛运行这些时间相关因素解析法的建模复杂度会急剧上升。非序贯蒙特卡洛虽然在网络规模上不设限但它抽样的是“瞬时状态”丢失了状态转移的时间信息而SAIFI和SAIDI恰恰需要停运频率和持续时间这就必须靠序贯方法才能直接算出来。1.3 什么情况下才值得用序贯蒙特卡洛我不建议一上来就无脑上序贯蒙特卡洛。如果只是评估一个简单的单电源辐射网负荷固定不变也没有任何转供策略解析法几分钟就算完了没必要开一个跑十万年仿真的重型武器。但如果系统满足下面任一条件序贯蒙特卡洛就是更合理的选择负荷随时间变化典型日曲线、季节波动存在修复顺序、计划检修、天气相关故障等时间逻辑有联络开关转供、储能调度、分布式电源孤岛等运行策略系统规模大、拓扑复杂解析枚举变得不现实。工程上写代码前先想清楚这一点能省大量时间。2. 可靠性评估的建模要点从元件到系统2.1 元件可靠性模型两状态与指数分布抽样在配电网可靠性评估里最常用的元件模型是两状态马尔可夫模型元件要么在正常运行状态要么在故障修复状态。从正常到故障由故障率λ驱动从故障到正常由修复率μ驱动μ 1/MTTRMTTR是平均修复时间。指数分布是这里最常见的假设故障前时间TTF服从参数为λ的指数分布修复时间TTR服从参数为μ的指数分布。指数分布的“无记忆性”意味着元件不会因为已经运行了很久就更容易老化故障这在长周期平均可靠性评估中是一种合理且工程上被广泛接受的近似数据也最容易获取——很多标准手册直接给出每公里线路的年故障次数和平均修复时间。抽样时用逆变换法。对于故障率λ单位次/年先换成小时单位再通过均匀随机数U转换function ttf sample_ttf(rate_per_year) % 根据年故障率抽样故障前时间返回单位小时 rate_per_hour rate_per_year / 8760; u rand(); u max(u, eps); % 防止 rand()0 导致 log(0) ttf -log(u) / rate_per_hour; end function ttr sample_ttr(mttr_hours) % 根据平均修复时间小时抽样修复时长 repair_rate 1 / mttr_hours; u rand(); u max(u, eps); ttr -log(u) / repair_rate; end这里可以直接用Matlab内置的exprnd但我更推荐显式写逆变换抽样一是少一次函数调用开销二是代码里能清楚看到指数分布抽样原理后面要扩展成威布尔分布或其他分布时也方便改。2.2 负荷与时变因素的建模思路如果只算固定负荷下的可靠性指标负荷模型很简单就是一个常数功率乘以停电时长就是缺供电量。但实际系统的日负荷曲线、季节性变化、用户类型差异都会影响停电损失。这就是序贯蒙特卡洛的一个突出优势它天然有一条时间轴可以在每个小时点上应用对应的负荷值。典型做法是构造一条8760点的年小时负荷曲线可以是实测历史数据也可以由“典型日曲线×季节系数×随机波动”合成。模拟过程中某个负荷点停电时长的统计是精确到小时的把停电时间段与负荷曲线对齐累加得到的就是真正的缺供电量。非序贯方法很难把这种时间匹配做好这也是很多论文和工程报告选择序贯模拟的直接原因。2.3 网络拓扑、保护与转供逻辑建模配电网不同于输电网通常是辐射状运行故障处理有一套固定逻辑故障发生后断路器跳闸然后通过分段开关和隔离开关逐步隔离故障段再通过联络开关对故障段下游可转供区域恢复供电最后对故障元件进行修复。在序贯蒙特卡洛模拟中每个时刻要判断“如果某个元件正处于故障状态那么哪些负荷点受影响是否已经被转供恢复”。这个过程本质上是网络连通性和转供能力分析。工程实现上通常有两种层次第一种是简化层次预计算每个元件故障导致的负荷点影响集合模拟时直接查表不考虑转供时序。这种适合代码原型。第二种是精细层次在故障发生和修复的时间点触发一次网络拓扑搜索考虑开关状态、转供路径容量、操作时间动态确定影响范围。代码复杂度高很多但能分析转供策略对可靠性的改善效果。我实际做项目时会先把简化层次跑通确认指标量级正确再加转供和操作逻辑避免一开始就被复杂的开关时序淹没。2.4 可靠性指标体系算哪些、怎么算设备和负荷建模完成之后需要明确最终要输出什么指标。下面这张表是配电网可靠性评估的标准输出指标名称计算公式含义SAIFI系统平均停电频率指标用户停电总次数 / 总用户数每个用户年平均停电次数SAIDI系统平均停电持续时间指标用户停电总时长 / 总用户数每个用户年平均停电小时数CAIDI用户平均停电持续时间SAIDI / SAIFI每次停电平均持续小时数ASAI平均供电可用率1 - SAIDI / 8760一年中供电可用的比例ENS缺供电量各负荷点缺供电量之和系统年缺供电量kWhAENS平均系统缺供电量ENS / 总用户数平均每户年缺电量所有指标都是“用户数加权”的不是简单对负荷点取平均。这一点特别容易写错比如一个工业用户和一个居民用户停电相同时间但工业用户停电次数和电量损失的口径完全不同计算系统指标时必须乘以各自的用户数。3. Matlab代码实现模块拆解与关键代码3.1 代码总体结构设计与数据组织整套代码我习惯拆成四个模块数据定义、元件生命周期生成、事件驱动模拟器、指标统计与收敛判断。系统规模还比较小时用结构体数组组织元件数据足够如果要做上百节点的配电网建议用表格table或者直接定义成class方便扩展。元件数据至少要包含这些字段元件编号、类型馈线段/变压器/断路器、年故障率λ、平均修复时间MTTR、上下游节点编号以及预计算的影响负荷点编号集合。一个典型的定义方式% 定义一个简化配电网3个负荷点、4条馈线段的算例 element struct(); element(1).id 1; element(1).type line; element(1).lambda 0.12; % 次/年 element(1).mttr 5; % 小时 element(1).from_node 0; % 电源侧节点 element(1).to_node 1; element(1).affected_loads [1]; % 拓扑用 graph 对象表示方便做连通性分析 % 节点0为电源 G graph(); G addedge(G, 0, 1); G addedge(G, 1, 2); G addedge(G, 1, 3);拓扑用Matlab自带的graph对象管理非常方便后面的连通性检查直接用conncomp或reachable就能实现。数据量更大时建议把元件参数统一存成CSV或Excel用readtable读进来写代码时不用反复改结构体。3.2 元件生命周期生成把随机过程变成事件列表每个元件都沿着“运行-故障-修复-运行”循环运动。生成生命周期时先给每个元件抽样第一个故障时间然后模拟过程就是不断回答两个问题下一次故障发生在什么时候、修复需要多久。这里有一个关键选择事件驱动还是时间步进。我强烈推荐事件驱动。配电网元件的平均无故障时间通常是几十小时到几千小时而平均修复时间只有几小时如果用小时步长模拟大量时间步里系统状态根本没变化纯属浪费计算量。事件驱动只在状态变化的时刻推进时间系统有一万个元件也就处理一万次左右的故障事件效率高几个数量级。3.3 事件驱动模拟主循环实现主循环是整套程序的心脏。逻辑是维护每个元件的next_fault时间和next_repair时间每一轮找出全局最早事件然后推进到该时刻处理状态变化和指标记录。故障发生时受影响负荷点的停电次数加一并记录停电开始时刻修复发生时累加停电时长和缺供电量。% 初始化 rng(2024); % 固定随机种子保证可复现 n_years 2000; sim_hours n_years * 8760; n_ele length(element); n_lp 3; % 负荷点数量 lp_load [200, 150, 180]; % 负荷点平均负荷kW lp_cust [250, 180, 220]; % 负荷点用户数 next_fault zeros(n_ele, 1); next_repair inf(n_ele, 1); for i 1:n_ele next_fault(i) sample_ttf(element(i).lambda); end outage_count zeros(n_lp, 1); % 停电次数 outage_duration zeros(n_lp, 1); % 停电时长小时 outage_ens zeros(n_lp, 1); % 缺供电量kWh outage_start inf(n_lp, 1); % 当前停电开始时刻 current_t 0; while current_t sim_hours [t_fault, idx_fault] min(next_fault); [t_repair, idx_repair] min(next_repair); if t_fault t_repair if t_fault sim_hours break; end current_t t_fault; % 故障元件进入故障状态 affected element(idx_fault).affected_loads; for j affected if outage_start(j) inf outage_count(j) outage_count(j) 1; outage_start(j) current_t; end end % 抽样修复时间 next_repair(idx_fault) current_t sample_ttr(element(idx_fault).mttr); next_fault(idx_fault) inf; else if t_repair sim_hours break; end current_t t_repair; % 修复完成恢复受影响负荷点 affected element(idx_repair).affected_loads; for j affected if outage_start(j) ~ inf outage_duration(j) outage_duration(j) current_t - outage_start(j); outage_ens(j) outage_ens(j) (current_t - outage_start(j)) * lp_load(j); outage_start(j) inf; end end % 重新抽样该元件下一次故障时间 next_fault(idx_repair) current_t sample_ttf(element(idx_repair).lambda); next_repair(idx_repair) inf; end end这段代码有个细节值得注意outage_start(j) inf 的判断。它维护了每个负荷点当前的停电状态解决了多个元件同时故障叠加影响同一负荷点时重复计数的问题。当A元件故障导致负荷停电后B元件又故障理论上该负荷点的停电时长不应从B故障时刻重新累计停电次数也不应再加一次。只有等所有故障元件都修复、负荷点彻底恢复供电后下一次停电才会重新计数。这是初写序贯蒙特卡洛最容易踩的坑之一很多代码算出的SAIFI虚高都是这里逻辑出了问题。3.4 故障影响分析通过图连通性快速判断主循环里用到的affected_loads可以提前预计算也可以在故障发生时动态分析。预计算的实现方式是对每个元件将其从图中移除然后寻找与电源节点失去连通的负荷节点function affected find_affected_loads(G, source_node, load_nodes, from, to) % 移除故障支路找出受影响负荷节点 G2 rmedge(G, from, to); bins conncomp(G2); source_bin bins(source_node); load_nodes load_nodes(:); affected load_nodes(bins(load_nodes) ~ source_bin); end这个函数的核心是Matlab的conncomp它返回图中每个节点所属的连通分量编号。电源节点所在的分量之外的负荷节点就是故障发生后失电的节点。这个逻辑默认了不会通过联络开关转供属于基础版。要加转供逻辑只需要在判断后用额外的图搜索检查哪些失电节点能通过联络线回到电源再把可转供节点从影响集合里剔除。很多实际算例里故障影响集合并不只是简单的“下游全部停电”因为配网有熔断器和分段开关分支线故障可能被熔断器隔离只影响很小范围。做精细评估时需要在元件数据里为每个元件指定其保护装置动作边界或者通过二次拓扑搜索判断隔离后的供电范围。这个话题比较大这里先点到为止。3.5 可靠性指标计算与收敛性判断模拟结束后把累计量按仿真年数归一化按用户数加权计算系统指标total_cust sum(lp_cust); SAIFI sum(outage_count .* lp_cust) / total_cust / n_years; SAIDI sum(outage_duration .* lp_cust) / total_cust / n_years; CAIDI SAIDI / SAIFI; ASAI 1 - SAIDI / 8760; ENS sum(outage_ens) / n_years; AENS ENS / total_cust; fprintf(SAIFI %.4f 次/户年\n, SAIFI); fprintf(SAIDI %.4f 小时/户年\n, SAIDI); fprintf(CAIDI %.4f 小时/次\n, CAIDI); fprintf(ASAI %.6f\n, ASAI); fprintf(ENS %.2f kWh/年\n, ENS);这里SAIFI和SAIDI分别是停电次数指标准确率、停电时长估计准确率的关键指标。仿真年数直接决定结果稳定性。不能随手填一个100年就交给领导科学性上站不住脚。规范做法是计算方差系数CV来判断收敛一般要求关键指标的CV小于5%function cv calc_cv(annual_values) n length(annual_values); cv std(annual_values) / (abs(mean(annual_values)) * sqrt(n)); end做法是在主循环里记录每年的负荷点停电次数和停电时长逐年统计出SAIFI和SAIDI的年度序列每跑完200年检查一次CV。CV达标就停止不达标就继续累加仿真年数。有人会担心仿真时间太长实际上对于中等规模配电网几千年的仿真在Matlab里通常几分钟就能完成瓶颈往往在于故障影响分析的频繁图搜索而不是事件循环本身。4. 实操踩坑记录与调试技巧4.1 先用能手算的小系统验证代码正确性这是我最想强调的一点。写完代码第一件事不要直接跑IEEE 33节点或者RBTS配电网先手工构造一个简单到能口算的系统比如一条单馈线带两个负荷点中间一段线路故障率固定、修复时间固定手算SAIFI和SAIDI的理论值。然后跑代码对比结果。数值对上了再逐步加复杂度。举一个具体例子馈线全长2公里故障率0.1次/公里年所以整条馈线等效故障率λ0.2次/年MTTR5小时。末端有100个用户平均负荷50kW。手算结果SAIFI0.2次/户年SAIDI0.2×51.0小时/户年ENS1.0×5050 kWh/年。如果代码跑2000年仿真的结果和这个值在统计误差范围内一致那核心逻辑基本没问题。这套校验方法能在一小时之内筛掉80%的隐性bug。4.2 常见问题速查表症状可能原因解决办法结果数量级完全不对故障率、修复时间的单位混用年/小时统一换算所有抽样函数内部除以8760SAIFI虚高多个元件同时故障时负荷点停电次数被重复累加用outage_start状态维护防止重复计数每次运行结果差异很大没固定随机种子主程序开头调用rng(seed)仿真很久不收敛仿真年数太少或方差系数目标过严用CV判断动态决定仿真时长rand()取到0导致infrand()理论上能返回0导致log(0)加max(u, eps)保护系统指标与手算不符指标计算时用户数权重没乘检查公式中lp_cust加权项修复完成后负荷点没恢复事件顺序处理有误修复时刻和故障时刻同值用严格不等号或加微小时间扰动这些坑我基本都在真实项目里踩过。最隐蔽的是第一个单位混用。故障率0.2次/年如果直接当成次/小时抽样出的TTF就会比真实值大8760倍结果就是跑一万年仿真也几乎看不到一次故障指标全部趋近于零。4.3 性能优化仿真跑得快的几个实用手段序贯蒙特卡洛模拟最大的软肋是计算量大但Matlab里有很多办法能把计算时间压下来。首先记住一个原则主循环内部尽量不要动态增长数组不要在循环里反复调用eval、subsref这类开销大的操作预分配所有数组。其次预计算元件的影响负荷点集合而不是每次故障都在循环里查图能显著提速。再次大规模系统可以按年拆分成独立任务用parfor并行跑最后把各worker的统计量汇总。蒙特卡洛天然适合这种并行拆分因为不同年份的模拟彼此独立。如果你做的算例特别大比如几百条馈线、几千个元件建议把影响负荷点集合的计算和主循环分开优化。前者涉及图搜索算法后者是高强度循环计算两个部分瓶颈不同分别剖析才有意义。用Matlab的profiler跑一遍通常能立刻定位到耗时最大的函数。4.4 从算例到论文/报告结果如何落地跑出一堆指标之后工作并没有结束。实际工程报告和论文里还需要做敏感性分析、方案对比和可视化。我会额外做两个图一个是系统SAIDI随仿真年份增长的收敛曲线证明结果稳定另一个是各负荷点停电次数和停电时长的柱状图直观暴露网络中的薄弱环节。这两个图说服力很强评审专家一眼就能看出你的代码是收敛的、结果是可信的。这个方法框架的可扩展性也值得一提。同一套序贯蒙特卡洛骨架改一下元件抽样分布就能考虑老化元件在故障影响分析里加点逻辑就能评估分布式电源孤岛、储能应急供电把负荷曲线改成随机场景就能做考虑不确定性的风险评估。我一直觉得写这套代码最值钱的不是那几行指标公式而是事件驱动的骨架和故障影响分析的接口设计它们决定了你能在多大程度上扩展这个程序去回答更多工程问题。5. 几点个人心得与调试手记代码写到今天我最大的体会是序贯蒙特卡洛模拟的数学原理并不高深真正拉开距离的是工程细节。指标定义、单位换算、状态维护、收敛判断每一个细节都可能导致结果差出数量级。别嫌麻烦一定要先搭一个能手算的小系统把所有逻辑验一遍再上规模。另外一个建议是把代码模块化做好。模拟器、抽样器、拓扑分析器、指标计算器各归各的文件以后换系统、换参数、加逻辑都方便。很多人图省事把几百行塞在一个脚本里调试三天就开始怀疑人生这个亏我吃过。最后分享一个调试技巧在事件驱动主循环里加一个“抽帧输出”比如每处理1000个事件就打印一次当前时间和累计指标能帮你快速发现“卡死”还是“正常推进”。如果模拟时间推进非常缓慢说明事件生成逻辑有问题。这个技巧在排查问题上比任何断点都好用。这套方法后续还能往很多方向扩展比如做配电网规划方案的可靠性比选、评估自动化开关带来的可靠性提升、分析分布式电源接入对供电可用率的影响。框架搭好了这些研究就只是往里填模块的事。
基于序贯蒙特卡洛的配电网可靠性评估：Matlab实现全解析
做配电网规划或供电可靠性分析的人大概率都绕不开一个问题某条馈线上接了哪些用户、线路多长、故障率多高一年下来用户平均停几次电、停多久、少供了多少电。工程上这个问题有一整套标准考核口径就是SAIFI、SAIDI、ENS这些可靠性指标。而当你面对的配电网越来越复杂时变负荷、联络转供、分布式电源这些因素全都掺进来的时候解析法会变得非常吃力这时候序贯蒙特卡洛模拟法就是最趁手的一把螺丝刀。这篇文章就围绕“用Matlab实现基于序贯蒙特卡洛模拟法的配电网可靠性评估”这件事把从建模、抽样、事件驱动模拟到指标统计的完整链路拆开讲一遍代码结构可以直接拿去改自己的算例适合正在做配电系统可靠性课题的研究生也适合电网公司或设计院里需要做供电能力评估和网架规划的工程师参考。1. 为什么选序贯蒙特卡洛方案选型的底层逻辑1.1 可靠性评估到底在算什么配电网可靠性评估的核心说白了就是回答三个问题每个负荷点一年大概停几次电、平均每次停多久、一年总共损失多少电量。围绕这三个答案行业里定义了一整套指标常用的有系统平均停电频率SAIFI、系统平均停电持续时间SAIDI、用户平均停电持续时间CAIDI、平均供电可用率ASAI以及缺供电量ENS。这些指标是电网公司内部考核、规划投资决策、供电合同约定里的硬数据。比如规划一条新的10kV馈线评审时一定会问这条馈线投运后系统SAIDI能改善多少下游新增用户对供电可靠性有多大的边际影响这些问题不靠拍脑袋必须拿出量化计算结果。计算这些指标本质上是一个概率分析问题。每个元件馈线段、配电变压器、断路器、隔离开关都有随机的故障和修复行为这些行为的随机组合决定了用户是否停电、停多久。可靠性评估的任务就是把这种随机性通过模型和算法翻译成期望值指标。1.2 解析法、非序贯与序贯蒙特卡洛的取舍行业内处理这个问题有三条技术路线我把它们的本质区别做个对比方法核心思路优势劣势解析法FMEA、最小割集、网络等值枚举所有故障模式按概率加权求和计算快、结果精确、适合简单辐射网系统规模大时枚举组合爆炸难处理时序因素非序贯蒙特卡洛按元件状态概率随机抽样一个系统快照能处理大规模网络、任意拓扑只“拍快照”不“拍电影”无法直接算停电频率和持续时间序贯蒙特卡洛按时间轴模拟每个元件的“正常-故障-修复-正常”全过程天然含时序能处理时变负荷、控制策略、复杂维修逻辑计算量大结果带统计波动解析法在单馈线、固定负荷、无控制策略的场景下非常好用手算公式都能出结果。但一旦系统里有时变负荷曲线、联络开关倒闸操作、分布式电源孤岛运行这些时间相关因素解析法的建模复杂度会急剧上升。非序贯蒙特卡洛虽然在网络规模上不设限但它抽样的是“瞬时状态”丢失了状态转移的时间信息而SAIFI和SAIDI恰恰需要停运频率和持续时间这就必须靠序贯方法才能直接算出来。1.3 什么情况下才值得用序贯蒙特卡洛我不建议一上来就无脑上序贯蒙特卡洛。如果只是评估一个简单的单电源辐射网负荷固定不变也没有任何转供策略解析法几分钟就算完了没必要开一个跑十万年仿真的重型武器。但如果系统满足下面任一条件序贯蒙特卡洛就是更合理的选择负荷随时间变化典型日曲线、季节波动存在修复顺序、计划检修、天气相关故障等时间逻辑有联络开关转供、储能调度、分布式电源孤岛等运行策略系统规模大、拓扑复杂解析枚举变得不现实。工程上写代码前先想清楚这一点能省大量时间。2. 可靠性评估的建模要点从元件到系统2.1 元件可靠性模型两状态与指数分布抽样在配电网可靠性评估里最常用的元件模型是两状态马尔可夫模型元件要么在正常运行状态要么在故障修复状态。从正常到故障由故障率λ驱动从故障到正常由修复率μ驱动μ 1/MTTRMTTR是平均修复时间。指数分布是这里最常见的假设故障前时间TTF服从参数为λ的指数分布修复时间TTR服从参数为μ的指数分布。指数分布的“无记忆性”意味着元件不会因为已经运行了很久就更容易老化故障这在长周期平均可靠性评估中是一种合理且工程上被广泛接受的近似数据也最容易获取——很多标准手册直接给出每公里线路的年故障次数和平均修复时间。抽样时用逆变换法。对于故障率λ单位次/年先换成小时单位再通过均匀随机数U转换function ttf sample_ttf(rate_per_year) % 根据年故障率抽样故障前时间返回单位小时 rate_per_hour rate_per_year / 8760; u rand(); u max(u, eps); % 防止 rand()0 导致 log(0) ttf -log(u) / rate_per_hour; end function ttr sample_ttr(mttr_hours) % 根据平均修复时间小时抽样修复时长 repair_rate 1 / mttr_hours; u rand(); u max(u, eps); ttr -log(u) / repair_rate; end这里可以直接用Matlab内置的exprnd但我更推荐显式写逆变换抽样一是少一次函数调用开销二是代码里能清楚看到指数分布抽样原理后面要扩展成威布尔分布或其他分布时也方便改。2.2 负荷与时变因素的建模思路如果只算固定负荷下的可靠性指标负荷模型很简单就是一个常数功率乘以停电时长就是缺供电量。但实际系统的日负荷曲线、季节性变化、用户类型差异都会影响停电损失。这就是序贯蒙特卡洛的一个突出优势它天然有一条时间轴可以在每个小时点上应用对应的负荷值。典型做法是构造一条8760点的年小时负荷曲线可以是实测历史数据也可以由“典型日曲线×季节系数×随机波动”合成。模拟过程中某个负荷点停电时长的统计是精确到小时的把停电时间段与负荷曲线对齐累加得到的就是真正的缺供电量。非序贯方法很难把这种时间匹配做好这也是很多论文和工程报告选择序贯模拟的直接原因。2.3 网络拓扑、保护与转供逻辑建模配电网不同于输电网通常是辐射状运行故障处理有一套固定逻辑故障发生后断路器跳闸然后通过分段开关和隔离开关逐步隔离故障段再通过联络开关对故障段下游可转供区域恢复供电最后对故障元件进行修复。在序贯蒙特卡洛模拟中每个时刻要判断“如果某个元件正处于故障状态那么哪些负荷点受影响是否已经被转供恢复”。这个过程本质上是网络连通性和转供能力分析。工程实现上通常有两种层次第一种是简化层次预计算每个元件故障导致的负荷点影响集合模拟时直接查表不考虑转供时序。这种适合代码原型。第二种是精细层次在故障发生和修复的时间点触发一次网络拓扑搜索考虑开关状态、转供路径容量、操作时间动态确定影响范围。代码复杂度高很多但能分析转供策略对可靠性的改善效果。我实际做项目时会先把简化层次跑通确认指标量级正确再加转供和操作逻辑避免一开始就被复杂的开关时序淹没。2.4 可靠性指标体系算哪些、怎么算设备和负荷建模完成之后需要明确最终要输出什么指标。下面这张表是配电网可靠性评估的标准输出指标名称计算公式含义SAIFI系统平均停电频率指标用户停电总次数 / 总用户数每个用户年平均停电次数SAIDI系统平均停电持续时间指标用户停电总时长 / 总用户数每个用户年平均停电小时数CAIDI用户平均停电持续时间SAIDI / SAIFI每次停电平均持续小时数ASAI平均供电可用率1 - SAIDI / 8760一年中供电可用的比例ENS缺供电量各负荷点缺供电量之和系统年缺供电量kWhAENS平均系统缺供电量ENS / 总用户数平均每户年缺电量所有指标都是“用户数加权”的不是简单对负荷点取平均。这一点特别容易写错比如一个工业用户和一个居民用户停电相同时间但工业用户停电次数和电量损失的口径完全不同计算系统指标时必须乘以各自的用户数。3. Matlab代码实现模块拆解与关键代码3.1 代码总体结构设计与数据组织整套代码我习惯拆成四个模块数据定义、元件生命周期生成、事件驱动模拟器、指标统计与收敛判断。系统规模还比较小时用结构体数组组织元件数据足够如果要做上百节点的配电网建议用表格table或者直接定义成class方便扩展。元件数据至少要包含这些字段元件编号、类型馈线段/变压器/断路器、年故障率λ、平均修复时间MTTR、上下游节点编号以及预计算的影响负荷点编号集合。一个典型的定义方式% 定义一个简化配电网3个负荷点、4条馈线段的算例 element struct(); element(1).id 1; element(1).type line; element(1).lambda 0.12; % 次/年 element(1).mttr 5; % 小时 element(1).from_node 0; % 电源侧节点 element(1).to_node 1; element(1).affected_loads [1]; % 拓扑用 graph 对象表示方便做连通性分析 % 节点0为电源 G graph(); G addedge(G, 0, 1); G addedge(G, 1, 2); G addedge(G, 1, 3);拓扑用Matlab自带的graph对象管理非常方便后面的连通性检查直接用conncomp或reachable就能实现。数据量更大时建议把元件参数统一存成CSV或Excel用readtable读进来写代码时不用反复改结构体。3.2 元件生命周期生成把随机过程变成事件列表每个元件都沿着“运行-故障-修复-运行”循环运动。生成生命周期时先给每个元件抽样第一个故障时间然后模拟过程就是不断回答两个问题下一次故障发生在什么时候、修复需要多久。这里有一个关键选择事件驱动还是时间步进。我强烈推荐事件驱动。配电网元件的平均无故障时间通常是几十小时到几千小时而平均修复时间只有几小时如果用小时步长模拟大量时间步里系统状态根本没变化纯属浪费计算量。事件驱动只在状态变化的时刻推进时间系统有一万个元件也就处理一万次左右的故障事件效率高几个数量级。3.3 事件驱动模拟主循环实现主循环是整套程序的心脏。逻辑是维护每个元件的next_fault时间和next_repair时间每一轮找出全局最早事件然后推进到该时刻处理状态变化和指标记录。故障发生时受影响负荷点的停电次数加一并记录停电开始时刻修复发生时累加停电时长和缺供电量。% 初始化 rng(2024); % 固定随机种子保证可复现 n_years 2000; sim_hours n_years * 8760; n_ele length(element); n_lp 3; % 负荷点数量 lp_load [200, 150, 180]; % 负荷点平均负荷kW lp_cust [250, 180, 220]; % 负荷点用户数 next_fault zeros(n_ele, 1); next_repair inf(n_ele, 1); for i 1:n_ele next_fault(i) sample_ttf(element(i).lambda); end outage_count zeros(n_lp, 1); % 停电次数 outage_duration zeros(n_lp, 1); % 停电时长小时 outage_ens zeros(n_lp, 1); % 缺供电量kWh outage_start inf(n_lp, 1); % 当前停电开始时刻 current_t 0; while current_t sim_hours [t_fault, idx_fault] min(next_fault); [t_repair, idx_repair] min(next_repair); if t_fault t_repair if t_fault sim_hours break; end current_t t_fault; % 故障元件进入故障状态 affected element(idx_fault).affected_loads; for j affected if outage_start(j) inf outage_count(j) outage_count(j) 1; outage_start(j) current_t; end end % 抽样修复时间 next_repair(idx_fault) current_t sample_ttr(element(idx_fault).mttr); next_fault(idx_fault) inf; else if t_repair sim_hours break; end current_t t_repair; % 修复完成恢复受影响负荷点 affected element(idx_repair).affected_loads; for j affected if outage_start(j) ~ inf outage_duration(j) outage_duration(j) current_t - outage_start(j); outage_ens(j) outage_ens(j) (current_t - outage_start(j)) * lp_load(j); outage_start(j) inf; end end % 重新抽样该元件下一次故障时间 next_fault(idx_repair) current_t sample_ttf(element(idx_repair).lambda); next_repair(idx_repair) inf; end end这段代码有个细节值得注意outage_start(j) inf 的判断。它维护了每个负荷点当前的停电状态解决了多个元件同时故障叠加影响同一负荷点时重复计数的问题。当A元件故障导致负荷停电后B元件又故障理论上该负荷点的停电时长不应从B故障时刻重新累计停电次数也不应再加一次。只有等所有故障元件都修复、负荷点彻底恢复供电后下一次停电才会重新计数。这是初写序贯蒙特卡洛最容易踩的坑之一很多代码算出的SAIFI虚高都是这里逻辑出了问题。3.4 故障影响分析通过图连通性快速判断主循环里用到的affected_loads可以提前预计算也可以在故障发生时动态分析。预计算的实现方式是对每个元件将其从图中移除然后寻找与电源节点失去连通的负荷节点function affected find_affected_loads(G, source_node, load_nodes, from, to) % 移除故障支路找出受影响负荷节点 G2 rmedge(G, from, to); bins conncomp(G2); source_bin bins(source_node); load_nodes load_nodes(:); affected load_nodes(bins(load_nodes) ~ source_bin); end这个函数的核心是Matlab的conncomp它返回图中每个节点所属的连通分量编号。电源节点所在的分量之外的负荷节点就是故障发生后失电的节点。这个逻辑默认了不会通过联络开关转供属于基础版。要加转供逻辑只需要在判断后用额外的图搜索检查哪些失电节点能通过联络线回到电源再把可转供节点从影响集合里剔除。很多实际算例里故障影响集合并不只是简单的“下游全部停电”因为配网有熔断器和分段开关分支线故障可能被熔断器隔离只影响很小范围。做精细评估时需要在元件数据里为每个元件指定其保护装置动作边界或者通过二次拓扑搜索判断隔离后的供电范围。这个话题比较大这里先点到为止。3.5 可靠性指标计算与收敛性判断模拟结束后把累计量按仿真年数归一化按用户数加权计算系统指标total_cust sum(lp_cust); SAIFI sum(outage_count .* lp_cust) / total_cust / n_years; SAIDI sum(outage_duration .* lp_cust) / total_cust / n_years; CAIDI SAIDI / SAIFI; ASAI 1 - SAIDI / 8760; ENS sum(outage_ens) / n_years; AENS ENS / total_cust; fprintf(SAIFI %.4f 次/户年\n, SAIFI); fprintf(SAIDI %.4f 小时/户年\n, SAIDI); fprintf(CAIDI %.4f 小时/次\n, CAIDI); fprintf(ASAI %.6f\n, ASAI); fprintf(ENS %.2f kWh/年\n, ENS);这里SAIFI和SAIDI分别是停电次数指标准确率、停电时长估计准确率的关键指标。仿真年数直接决定结果稳定性。不能随手填一个100年就交给领导科学性上站不住脚。规范做法是计算方差系数CV来判断收敛一般要求关键指标的CV小于5%function cv calc_cv(annual_values) n length(annual_values); cv std(annual_values) / (abs(mean(annual_values)) * sqrt(n)); end做法是在主循环里记录每年的负荷点停电次数和停电时长逐年统计出SAIFI和SAIDI的年度序列每跑完200年检查一次CV。CV达标就停止不达标就继续累加仿真年数。有人会担心仿真时间太长实际上对于中等规模配电网几千年的仿真在Matlab里通常几分钟就能完成瓶颈往往在于故障影响分析的频繁图搜索而不是事件循环本身。4. 实操踩坑记录与调试技巧4.1 先用能手算的小系统验证代码正确性这是我最想强调的一点。写完代码第一件事不要直接跑IEEE 33节点或者RBTS配电网先手工构造一个简单到能口算的系统比如一条单馈线带两个负荷点中间一段线路故障率固定、修复时间固定手算SAIFI和SAIDI的理论值。然后跑代码对比结果。数值对上了再逐步加复杂度。举一个具体例子馈线全长2公里故障率0.1次/公里年所以整条馈线等效故障率λ0.2次/年MTTR5小时。末端有100个用户平均负荷50kW。手算结果SAIFI0.2次/户年SAIDI0.2×51.0小时/户年ENS1.0×5050 kWh/年。如果代码跑2000年仿真的结果和这个值在统计误差范围内一致那核心逻辑基本没问题。这套校验方法能在一小时之内筛掉80%的隐性bug。4.2 常见问题速查表症状可能原因解决办法结果数量级完全不对故障率、修复时间的单位混用年/小时统一换算所有抽样函数内部除以8760SAIFI虚高多个元件同时故障时负荷点停电次数被重复累加用outage_start状态维护防止重复计数每次运行结果差异很大没固定随机种子主程序开头调用rng(seed)仿真很久不收敛仿真年数太少或方差系数目标过严用CV判断动态决定仿真时长rand()取到0导致infrand()理论上能返回0导致log(0)加max(u, eps)保护系统指标与手算不符指标计算时用户数权重没乘检查公式中lp_cust加权项修复完成后负荷点没恢复事件顺序处理有误修复时刻和故障时刻同值用严格不等号或加微小时间扰动这些坑我基本都在真实项目里踩过。最隐蔽的是第一个单位混用。故障率0.2次/年如果直接当成次/小时抽样出的TTF就会比真实值大8760倍结果就是跑一万年仿真也几乎看不到一次故障指标全部趋近于零。4.3 性能优化仿真跑得快的几个实用手段序贯蒙特卡洛模拟最大的软肋是计算量大但Matlab里有很多办法能把计算时间压下来。首先记住一个原则主循环内部尽量不要动态增长数组不要在循环里反复调用eval、subsref这类开销大的操作预分配所有数组。其次预计算元件的影响负荷点集合而不是每次故障都在循环里查图能显著提速。再次大规模系统可以按年拆分成独立任务用parfor并行跑最后把各worker的统计量汇总。蒙特卡洛天然适合这种并行拆分因为不同年份的模拟彼此独立。如果你做的算例特别大比如几百条馈线、几千个元件建议把影响负荷点集合的计算和主循环分开优化。前者涉及图搜索算法后者是高强度循环计算两个部分瓶颈不同分别剖析才有意义。用Matlab的profiler跑一遍通常能立刻定位到耗时最大的函数。4.4 从算例到论文/报告结果如何落地跑出一堆指标之后工作并没有结束。实际工程报告和论文里还需要做敏感性分析、方案对比和可视化。我会额外做两个图一个是系统SAIDI随仿真年份增长的收敛曲线证明结果稳定另一个是各负荷点停电次数和停电时长的柱状图直观暴露网络中的薄弱环节。这两个图说服力很强评审专家一眼就能看出你的代码是收敛的、结果是可信的。这个方法框架的可扩展性也值得一提。同一套序贯蒙特卡洛骨架改一下元件抽样分布就能考虑老化元件在故障影响分析里加点逻辑就能评估分布式电源孤岛、储能应急供电把负荷曲线改成随机场景就能做考虑不确定性的风险评估。我一直觉得写这套代码最值钱的不是那几行指标公式而是事件驱动的骨架和故障影响分析的接口设计它们决定了你能在多大程度上扩展这个程序去回答更多工程问题。5. 几点个人心得与调试手记代码写到今天我最大的体会是序贯蒙特卡洛模拟的数学原理并不高深真正拉开距离的是工程细节。指标定义、单位换算、状态维护、收敛判断每一个细节都可能导致结果差出数量级。别嫌麻烦一定要先搭一个能手算的小系统把所有逻辑验一遍再上规模。另外一个建议是把代码模块化做好。模拟器、抽样器、拓扑分析器、指标计算器各归各的文件以后换系统、换参数、加逻辑都方便。很多人图省事把几百行塞在一个脚本里调试三天就开始怀疑人生这个亏我吃过。最后分享一个调试技巧在事件驱动主循环里加一个“抽帧输出”比如每处理1000个事件就打印一次当前时间和累计指标能帮你快速发现“卡死”还是“正常推进”。如果模拟时间推进非常缓慢说明事件生成逻辑有问题。这个技巧在排查问题上比任何断点都好用。这套方法后续还能往很多方向扩展比如做配电网规划方案的可靠性比选、评估自动化开关带来的可靠性提升、分析分布式电源接入对供电可用率的影响。框架搭好了这些研究就只是往里填模块的事。