尧图网站建设 
简介面向电力系统可靠性研究领域的学者与工程人员这份基于序贯蒙特卡洛模拟法的配电网可靠性评估Matlab实现可直接用于分析配电网在随机故障下的供电可靠性水平。压缩包共包含3个文件均为Matlab脚本.m整体仅3KB轻量但逻辑完整。资源内置IEEE RBTS可靠性测试系统原始参数、IEEE33节点系统参数及RBTS BUS6参数的Matlab数据文件便于研究者快速替换拓扑与负荷数据。评估主程序采用节点影响分析法判断故障影响的负荷范围再通过序贯蒙特卡洛算法模拟时序运行状态从而输出可靠性指标。目前已有2802人学习下载适合需要开展配电网可靠性定量评估、熟悉序贯蒙特卡洛实现细节的电气工程学生与研究人员使用。 做配电网可靠性评估很多人一开始会扑向解析法想着把故障枚举、最小割集、故障树全用上结果拓扑稍微复杂一点就陷入保护动作逻辑的泥潭。我自己就干过这种事写一个带分支线和多级隔离开关的辐射网故障影响分析绕来绕去改了三版逻辑最后差点把项目做完。后来换成序贯蒙特卡洛模拟法整体思路瞬间通透不用穷举状态让元件按照自己的寿命规律“演”一遍又一遍把停电事件统计出来指标自然而然就有了。这篇文章用MATLAB把序贯蒙特卡洛模拟法的配电网可靠性评估完整拆解一遍覆盖原理、数据结构、核心抽样式、FMEA表构建、主循环代码、收敛判断和提速技巧。适合电力系统方向的研究生、做配电网规划或运行分析的工程师以及需要快速验证“某个网架方案可靠性到底行不行”的设计人员。不需要额外工具箱MATLAB基础功能就能跑通。1. 序贯蒙特卡洛法解决什么问题1.1 解析法在配电网里的尴尬解析法的核心是把所有可能的状态列出来算概率再加权。发输电系统用这套问题不大因为元件数量少、事件相对清晰但配电网是辐射结构加上断路器、隔离开关、熔断器、分段开关这一堆设备故障后的隔离和恢复路径非常依赖网络拓扑。举个简单例子主馈线中段短路保护动作后故障点上游有一部分负荷可以通过隔离开关快速恢复故障点下游则要等修复这个“哪些负荷停多久”的判断本身就带时序逻辑。解析法做这类分析不是不行但每次调整保护配置或网架结构都要重新梳理所有故障场景代码量和出错概率直线上升。更麻烦的是一旦考虑分布式电源孤岛、储能调度这类时间相关的恢复策略解析法的状态空间会爆炸因为它难以描述“时间先后”带来的影响。1.2 序贯蒙特卡洛到底在做什么序贯蒙特卡洛的思路用大白话说就是给每个元件编一本“寿命剧本”记录它什么时候坏、坏了多久、修好之后又什么时候再坏然后把所有元件剧本叠在同一个时间轴上看每一时刻系统处于什么状态。重复模拟几万甚至几十万年把每次停电事件的频率和时长统计出来就能得到可靠性指标。它和非序贯蒙特卡洛最大的区别在于是否关心时间顺序。非序贯法只抽样“某个元件此刻是坏还是好”相当于给系统拍一张瞬时快照适合发输电系统的充裕度评估序贯法则关注元件状态的演变过程能精确表达“故障发生—保护动作—隔离—修复—恢复供电”这个完整链条正好命中配电网可靠性评估的需求。1.3 为什么配电网特别适合时序模拟配电网里的大多数可靠性指标本质上都和时间强相关。SAIFI关心停电次数SAIDI关心停电时长CAIDI是两者比值这些指标离开了时间轴就说不清楚。序贯蒙特卡洛天然输出事件发生时刻和持续时间统计指标只是水到渠成的事。另一个优势是它不依赖“故障事件彼此独立”这类强假设。多重故障、计划检修叠加、天气导致的故障率时变、多状态元件这些在解析法里异常棘手的情况在序贯模拟中只是往抽样模型里加规则而已。所以我一直觉得配电网可靠性评估用序贯蒙特卡洛不是“退而求其次”而是更贴近工程实际的选择。2. 评估指标与模型假设2.1 负荷点和系统两级指标配电网可靠性评估通常分两级看。负荷点指标服务于单个用户系统指标服务于整个供电区域的横向对比。负荷点指标包括年平均故障次数 λ次/年年平均停电时间 U小时/年平均故障持续时间 r小时/次r U / λ系统指标则是在所有负荷点指标基础上加权汇总最常用的几个SAIFI系统平均停电频率指标单位次/年SAIDI系统平均停电持续时间指标单位小时/年CAIDI用户平均停电持续时间指标CAIDI SAIDI / SAIFIASAI供电可用率通常用百分比表示ENS年期望缺供电量单位kWh/年AENS平均系统缺供电量单位kWh/户/年这些指标的物理意义很简单把一次停电当成一个事件累加次数和时长按用户数或负荷量加权平均。真正考验代码实现的是怎么准确统计每次停电事件的次数和时长而不是指标公式本身。2.2 两状态模型和抽样式描述元件状态最常用的模型是两状态马尔可夫模型正常运行和故障停运。元件在正常状态停留的时间TTF和故障状态停留的时间TTR都服从指数分布。TTF -ln(U1) / λTTR -ln(U2) / μ其中 λ 是年故障率μ 是修复率U1、U2 是0到1之间的均匀分布随机数。注意这里的细节如果 λ 单位是“次/年”那么 TTF 的单位就是“年”如果 μ 单位是“次/小时”TTR 单位就是“小时”。混用单位是新手最容易踩的坑后面排查部分会再强调。指数分布有一个无记忆性所以每次抽出的 TTF 和 TTR 是相互独立的不需要记录上一次修好多久了。这也是整个序贯模拟能够大幅简化的基础。2.3 故障影响的FMEA思维FMEA失效模式与影响分析在配电网可靠性里的用法是把“某个元件故障后哪些负荷点停多久”事先梳理成表。比如某个馈线区段故障保护断路器跳开下游负荷全部失电接着运行人员操作隔离开关隔离故障段这时故障点上游且能切换的负荷短时停电后恢复故障点下游的负荷则要等故障元件修复才能恢复供电。这张FMEA表是连接“元件状态抽样”和“可靠性指标统计”的桥梁。模拟中一旦抽样出某个元件发生故障直接查表就能知道每个负荷点的停电时长。这样把时序模拟和网络拓扑解耦代码结构清晰计算速度也快得多。2.4 工程上的常用简化实际工程计算中通常先做几个合理简化元件只考虑正常运行和故障停运两态忽略计划检修或者作为修改后的参数叠加进来保护装置默认可靠动作不考虑变电站内部母线故障的细节假设修复完成后立刻恢复全部负荷。这些简化让模型门槛降低又不至于影响核心指标的准确性。如果你的场景需要精细建模比如配电网自愈或微网孤岛那就在FMEA表里补充控制策略的时延参数原理不变。3. MATLAB实现步骤拆解3.1 数据结构设计写序贯蒙特卡洛的程序数据结构设计比代码本身更影响维护成本。我习惯分三块组织元件参数的向量包括年故障率 lambda 和平均修复时间 MTTR小时用普通数组存避免循环时反复访问结构体字段。负荷点的基础数据包括用户数、平均负荷。和元件参数一样用数组存。FMEA表用一个二维矩阵或 cell 数组描述元件故障对负荷点的影响。行对应元件列对应负荷点元素是该负荷点的停电时长。MATLAB里的结构体读起来直观但循环体内部频繁访问结构体字段会拖慢速度。跑几十万年仿真时建议把字段拆成独立向量性能差距在数据量大时非常明显。3.2 指数分布抽样与事件推进核心抽样式很简单用 MATLAB 内置的 rand 函数生成均匀分布随机数再取负对数换算成指数分布抽样% 生成指数分布抽样值 % lambda: 故障率(次/年)返回 TTF(年) ttf -log(rand) / lambda(k); % MTTR_h: 平均修复时间(小时)返回 TTR(年) ttr -log(rand) * MTTR_h(k) / 8760;这里的单位换算很重要ttf 以年为单位ttr 从小时换算成年。如果想全部统一成小时两者都换算到小时维度也可以关键是统计阶段保持单位一致。事件推进有两种方案。最简单的是按年循环每一年内对每个元件独立推进时钟更严谨的是全程时间推进总模拟时长拉通彻底避免年份边界截断问题。建议直接用全程时间推进代码反而更简单。3.3 搭建配电网FMEA表FMEA表是整个程序里最需要动脑子的部分。以图1所示的简化辐射网为例电源 ── C1 ── L1 ── S1 ── L2 ── S2 ── L3 ── T1 ── LP1 │ │ L4 L5 │ │ LP2 LP3C1是出口断路器S1、S2是隔离开关L1到L5是馈线区段T1是配电变压器LP1到LP3是负荷点。假设故障后运行人员操作隔离开关隔离故障隔离操作时间为1小时元件修复时间按各自MTTR计算。这张网络的FMEA表可以写成如下形式故障元件LP1LP2LP3C1修复时间修复时间修复时间L1修复时间修复时间修复时间L2修复时间隔离时间隔离时间L3修复时间隔离时间修复时间L4隔离时间修复时间隔离时间L5隔离时间隔离时间修复时间T1修复时间修复时间修复时间这张表的逻辑是主馈线L1故障时故障段前没有可操作隔离开关全部负荷停到修复完成L2故障时S1、S2配合隔离故障段LP2和LP3不经过故障段可以通过切换恢复LP1在L2下游必须等修复L3故障时故障段位于LP1上游所以LP1要等修复LP3位于L3下游分支也要等修复LP2则可以短时切换恢复。每个网络的FMEA表都要根据实际保护配置重新推导这是不能偷懒的步骤。3.4 主循环与指标统计主循环推荐采用每个元件独立推进、时间轴拉通的方式%% 参数设置 rng(2024); % 固定随机种子结果可复现 total_year 100000; % 模拟总年限 total_hours total_year * 8760; % 转换为小时 K length(lambda); % 元件数量 N_LP length(LP_user); % 负荷点数量 % 指标累积变量 N_evt zeros(1, N_LP); % 负荷点停电次数 U_hour zeros(1, N_LP); % 负荷点停电小时数 ENS_kwh zeros(1, N_LP); % 负荷点缺供电量 for k 1:K t 0; % 当前元件的时间轴 while t total_hours % 抽样正常运行时间 ttf_h -log(rand) / lambda(k) * 8760; t t ttf_h; if t total_hours break; end % 抽样修复时间 ttr_h -log(rand) * MTTR_h(k); % 根据FMEA表统计负荷点影响 for i 1:N_LP outage_h FMEA(k, i); U_hour(i) U_hour(i) outage_h; N_evt(i) N_evt(i) 1; ENS_kwh(i) ENS_kwh(i) LP_P(i) * outage_h / 8760; end % 修复完成后继续推进时间 t t ttr_h; end end这个循环的逻辑很直白每个元件在自己的时间轴上不断经历“运行—故障—修复”的循环每发生一次故障就查FMEA表把停电次数、停电小时数、缺供电量累加到对应负荷点。因为元件之间独立推进所以即便某个修复过程跨越了模拟总时长前面已经发生并统计完毕的故障事件也不受影响。指标计算直接套公式%% 负荷点指标 lambda_lp N_evt / total_year; % 次/年 U_lp U_hour / total_year; % 小时/年 r_lp U_lp ./ max(lambda_lp, 1e-12); % 小时/次 %% 系统指标 total_user sum(LP_user); SAIFI sum(lambda_lp .* LP_user) / total_user; SAIDI sum(U_lp .* LP_user) / total_user; CAIDI SAIDI / SAIFI; ASAI (total_year * 8760 * total_user - sum(U_hour .* LP_user)) / ... (total_year * 8760 * total_user); ENS sum(ENS_kwh) / total_year;3.5 收敛性判断不能省略蒙特卡洛模拟最大的陷阱是“看起来跑完了实际上没收敛”。判断收敛最常用的指标是方差系数也就是变异系数。做法是把总模拟分成若干块每块独立计算某个指标得到一串样本然后算样本的标准差和均值看比值是否小于一个阈值。block_num 50; block_year 2000; SAIDI_block zeros(1, block_num); for b 1:block_num % 跑一个block记录该块的SAIDI SAIDI_block(b) simulate_block(block_year); end beta std(SAIDI_block) / mean(SAIDI_block) / sqrt(block_num); if beta 0.05 disp(结果已收敛); end0.05是很多人用的经验阈值要求高就压到0.01。我自己的习惯是先跑小规模测试确定量级再整批跑正式结果既省时间又不至于盲目加大模拟量。4. 典型算例与结果校验4.1 算例系统参数以3.3节那个简化辐射网为例元件参数参考常见配电网可靠性测试系统中的量级设定元件故障率次/年平均修复时间小时断路器 C10.0024馈线区段 L1~L30.105馈线分支 L4~L50.085配电变压器 T10.015200负荷点数据按表设定隔离操作时间取1小时负荷点用户数平均负荷kWLP15030LP23020LP340254.2 运行结果参考量级按上述参数模拟10万年典型结果大致落在以下量级指标参考值SAIFI0.2~0.5 次/年SAIDI1~3 小时/年CAIDI4~8 小时/次ASAI99.95%~99.99%ENS几十到一百多 kWh/年不同网络拓扑、不同隔离策略会让数值明显变化但量级基本一致。如果你的结果偏离这个范围很远优先检查单位换算和FMEA表。4.3 验证实现正确性的三个套路第一个套路是手算简单场景对比。挑一个单负荷点、双元件的小系统用解析法手算理论指标再跑蒙特卡洛对比误差应该在模拟误差范围内。第二个套路是故意放大故障率。把所有元件的故障率同时放大100倍此时故障事件非常频繁蒙特卡洛结果和解析法结果的相对误差会缩小因为事件样本数足够多。如果放大后结果差异还很大说明逻辑有bug而不是抽样误差。第三个套路是检查事件次数。如果元件k的故障率是λ模拟Y年那么预期的故障次数大约是λ×Y次。程序跑完后统计每个元件的故障次数和这个预期值对比偏离太远说明抽样式或循环推进写错了。5. 踩坑记录与提速建议5.1 常见问题速查表写序贯蒙特卡洛的过程中我见过的问题集中在下面几个地方问题现象可能原因解决办法结果每次跑都不一样且差异很大没有固定随机种子或模拟年限不足用rng固定种子检查方差系数停电次数偏少TTF单位错误用了次/小时却没换算统一时间单位建议全用小时SAIDI偏小隔离操作时间没有计入FMEA表FMEA表补充短时停电事件仿真耗时过长循环内重复访问结构体字段用向量存储参数减少字段访问修复时间出现异常大值指数分布长尾本身就存在确认抽样公式是-log(rand)不是log(rand)5.2 经验提速技巧指数分布抽样直接用-log(rand)不要用exprnd。虽然数学上等价但exprnd有额外的函数调用和参数检查开销在几十万次抽样的规模下累积差距很可观。FMEA表尽量在模拟前一次性构建不要在每次故障事件发生时实时搜索网络拓扑。我见过有人用邻接矩阵DFS搜索故障影响范围思路没错但每次故障都搜一遍图模拟速度会慢到让人怀疑人生。稳妥做法是离线构建FMEA表模拟期间只查表。如果网络特别大可以用“批处理随机数”技巧。先一次性生成大矩阵rand(M, N)再按索引取用减少 rand 函数的调用次数。不过这个优化对中小规模网络收益有限代码可读性会下降按需取舍。5.3 重要抽样加速收敛序贯蒙特卡洛的一个明显痛点是配电网元件故障率通常很低0.01到0.1次/年模拟中大量时间消耗在“什么也没发生”的时段。为了加速收敛可以用重要抽样人为把故障率放大k倍让故障事件更频繁地出现统计指标时把贡献除以k作为修正。scale 20; % 放大倍数 lambda_scaled lambda * scale; % 用放大后的故障率抽样 % 统计时停电次数贡献除以scale N_evt(i) N_evt(i) 1 / scale;这个技巧能显著缩短达到相同方差系数的模拟时间代价是引入偏差修正逻辑需要你清楚自己的抽样分布变化。我建议先用无偏的普通模拟跑通验证再引入重要抽样优化大规模算例。6. 一些个人体会把序贯蒙特卡洛在配电网可靠性评估里跑通之后我最大的感受是它其实没有想象中那么高不可攀核心逻辑就是“元件故障—查FMEA表—累计指标”难点反而在参数单位、FMEA构建和收敛判断这些看起来不起眼的细节上。跑通一次之后后续加分布式电源、储能孤岛、保护策略调整都是在现有框架上扩展事件规则整个代码结构不用推倒重来。最后分享一个小技巧做项目汇报时除了给出最终的SAIFI、SAIDI数值建议顺手画出“指标随模拟年限变化的收敛曲线”。一张收敛曲线比任何口头说明都有说服力能直接告诉对方你的结果可信度有多高也能帮你提前发现“跑得不够久”的问题。在做完这些之后你会发现所谓可靠性评估其实等同于一场足够耐心的统计实验。本文还有配套的精品资源点击获取 基于序贯蒙特卡洛的配电网可靠性评估及MATLAB实现
简介面向电力系统可靠性研究领域的学者与工程人员这份基于序贯蒙特卡洛模拟法的配电网可靠性评估Matlab实现可直接用于分析配电网在随机故障下的供电可靠性水平。压缩包共包含3个文件均为Matlab脚本.m整体仅3KB轻量但逻辑完整。资源内置IEEE RBTS可靠性测试系统原始参数、IEEE33节点系统参数及RBTS BUS6参数的Matlab数据文件便于研究者快速替换拓扑与负荷数据。评估主程序采用节点影响分析法判断故障影响的负荷范围再通过序贯蒙特卡洛算法模拟时序运行状态从而输出可靠性指标。目前已有2802人学习下载适合需要开展配电网可靠性定量评估、熟悉序贯蒙特卡洛实现细节的电气工程学生与研究人员使用。 做配电网可靠性评估很多人一开始会扑向解析法想着把故障枚举、最小割集、故障树全用上结果拓扑稍微复杂一点就陷入保护动作逻辑的泥潭。我自己就干过这种事写一个带分支线和多级隔离开关的辐射网故障影响分析绕来绕去改了三版逻辑最后差点把项目做完。后来换成序贯蒙特卡洛模拟法整体思路瞬间通透不用穷举状态让元件按照自己的寿命规律“演”一遍又一遍把停电事件统计出来指标自然而然就有了。这篇文章用MATLAB把序贯蒙特卡洛模拟法的配电网可靠性评估完整拆解一遍覆盖原理、数据结构、核心抽样式、FMEA表构建、主循环代码、收敛判断和提速技巧。适合电力系统方向的研究生、做配电网规划或运行分析的工程师以及需要快速验证“某个网架方案可靠性到底行不行”的设计人员。不需要额外工具箱MATLAB基础功能就能跑通。1. 序贯蒙特卡洛法解决什么问题1.1 解析法在配电网里的尴尬解析法的核心是把所有可能的状态列出来算概率再加权。发输电系统用这套问题不大因为元件数量少、事件相对清晰但配电网是辐射结构加上断路器、隔离开关、熔断器、分段开关这一堆设备故障后的隔离和恢复路径非常依赖网络拓扑。举个简单例子主馈线中段短路保护动作后故障点上游有一部分负荷可以通过隔离开关快速恢复故障点下游则要等修复这个“哪些负荷停多久”的判断本身就带时序逻辑。解析法做这类分析不是不行但每次调整保护配置或网架结构都要重新梳理所有故障场景代码量和出错概率直线上升。更麻烦的是一旦考虑分布式电源孤岛、储能调度这类时间相关的恢复策略解析法的状态空间会爆炸因为它难以描述“时间先后”带来的影响。1.2 序贯蒙特卡洛到底在做什么序贯蒙特卡洛的思路用大白话说就是给每个元件编一本“寿命剧本”记录它什么时候坏、坏了多久、修好之后又什么时候再坏然后把所有元件剧本叠在同一个时间轴上看每一时刻系统处于什么状态。重复模拟几万甚至几十万年把每次停电事件的频率和时长统计出来就能得到可靠性指标。它和非序贯蒙特卡洛最大的区别在于是否关心时间顺序。非序贯法只抽样“某个元件此刻是坏还是好”相当于给系统拍一张瞬时快照适合发输电系统的充裕度评估序贯法则关注元件状态的演变过程能精确表达“故障发生—保护动作—隔离—修复—恢复供电”这个完整链条正好命中配电网可靠性评估的需求。1.3 为什么配电网特别适合时序模拟配电网里的大多数可靠性指标本质上都和时间强相关。SAIFI关心停电次数SAIDI关心停电时长CAIDI是两者比值这些指标离开了时间轴就说不清楚。序贯蒙特卡洛天然输出事件发生时刻和持续时间统计指标只是水到渠成的事。另一个优势是它不依赖“故障事件彼此独立”这类强假设。多重故障、计划检修叠加、天气导致的故障率时变、多状态元件这些在解析法里异常棘手的情况在序贯模拟中只是往抽样模型里加规则而已。所以我一直觉得配电网可靠性评估用序贯蒙特卡洛不是“退而求其次”而是更贴近工程实际的选择。2. 评估指标与模型假设2.1 负荷点和系统两级指标配电网可靠性评估通常分两级看。负荷点指标服务于单个用户系统指标服务于整个供电区域的横向对比。负荷点指标包括年平均故障次数 λ次/年年平均停电时间 U小时/年平均故障持续时间 r小时/次r U / λ系统指标则是在所有负荷点指标基础上加权汇总最常用的几个SAIFI系统平均停电频率指标单位次/年SAIDI系统平均停电持续时间指标单位小时/年CAIDI用户平均停电持续时间指标CAIDI SAIDI / SAIFIASAI供电可用率通常用百分比表示ENS年期望缺供电量单位kWh/年AENS平均系统缺供电量单位kWh/户/年这些指标的物理意义很简单把一次停电当成一个事件累加次数和时长按用户数或负荷量加权平均。真正考验代码实现的是怎么准确统计每次停电事件的次数和时长而不是指标公式本身。2.2 两状态模型和抽样式描述元件状态最常用的模型是两状态马尔可夫模型正常运行和故障停运。元件在正常状态停留的时间TTF和故障状态停留的时间TTR都服从指数分布。TTF -ln(U1) / λTTR -ln(U2) / μ其中 λ 是年故障率μ 是修复率U1、U2 是0到1之间的均匀分布随机数。注意这里的细节如果 λ 单位是“次/年”那么 TTF 的单位就是“年”如果 μ 单位是“次/小时”TTR 单位就是“小时”。混用单位是新手最容易踩的坑后面排查部分会再强调。指数分布有一个无记忆性所以每次抽出的 TTF 和 TTR 是相互独立的不需要记录上一次修好多久了。这也是整个序贯模拟能够大幅简化的基础。2.3 故障影响的FMEA思维FMEA失效模式与影响分析在配电网可靠性里的用法是把“某个元件故障后哪些负荷点停多久”事先梳理成表。比如某个馈线区段故障保护断路器跳开下游负荷全部失电接着运行人员操作隔离开关隔离故障段这时故障点上游且能切换的负荷短时停电后恢复故障点下游的负荷则要等故障元件修复才能恢复供电。这张FMEA表是连接“元件状态抽样”和“可靠性指标统计”的桥梁。模拟中一旦抽样出某个元件发生故障直接查表就能知道每个负荷点的停电时长。这样把时序模拟和网络拓扑解耦代码结构清晰计算速度也快得多。2.4 工程上的常用简化实际工程计算中通常先做几个合理简化元件只考虑正常运行和故障停运两态忽略计划检修或者作为修改后的参数叠加进来保护装置默认可靠动作不考虑变电站内部母线故障的细节假设修复完成后立刻恢复全部负荷。这些简化让模型门槛降低又不至于影响核心指标的准确性。如果你的场景需要精细建模比如配电网自愈或微网孤岛那就在FMEA表里补充控制策略的时延参数原理不变。3. MATLAB实现步骤拆解3.1 数据结构设计写序贯蒙特卡洛的程序数据结构设计比代码本身更影响维护成本。我习惯分三块组织元件参数的向量包括年故障率 lambda 和平均修复时间 MTTR小时用普通数组存避免循环时反复访问结构体字段。负荷点的基础数据包括用户数、平均负荷。和元件参数一样用数组存。FMEA表用一个二维矩阵或 cell 数组描述元件故障对负荷点的影响。行对应元件列对应负荷点元素是该负荷点的停电时长。MATLAB里的结构体读起来直观但循环体内部频繁访问结构体字段会拖慢速度。跑几十万年仿真时建议把字段拆成独立向量性能差距在数据量大时非常明显。3.2 指数分布抽样与事件推进核心抽样式很简单用 MATLAB 内置的 rand 函数生成均匀分布随机数再取负对数换算成指数分布抽样% 生成指数分布抽样值 % lambda: 故障率(次/年)返回 TTF(年) ttf -log(rand) / lambda(k); % MTTR_h: 平均修复时间(小时)返回 TTR(年) ttr -log(rand) * MTTR_h(k) / 8760;这里的单位换算很重要ttf 以年为单位ttr 从小时换算成年。如果想全部统一成小时两者都换算到小时维度也可以关键是统计阶段保持单位一致。事件推进有两种方案。最简单的是按年循环每一年内对每个元件独立推进时钟更严谨的是全程时间推进总模拟时长拉通彻底避免年份边界截断问题。建议直接用全程时间推进代码反而更简单。3.3 搭建配电网FMEA表FMEA表是整个程序里最需要动脑子的部分。以图1所示的简化辐射网为例电源 ── C1 ── L1 ── S1 ── L2 ── S2 ── L3 ── T1 ── LP1 │ │ L4 L5 │ │ LP2 LP3C1是出口断路器S1、S2是隔离开关L1到L5是馈线区段T1是配电变压器LP1到LP3是负荷点。假设故障后运行人员操作隔离开关隔离故障隔离操作时间为1小时元件修复时间按各自MTTR计算。这张网络的FMEA表可以写成如下形式故障元件LP1LP2LP3C1修复时间修复时间修复时间L1修复时间修复时间修复时间L2修复时间隔离时间隔离时间L3修复时间隔离时间修复时间L4隔离时间修复时间隔离时间L5隔离时间隔离时间修复时间T1修复时间修复时间修复时间这张表的逻辑是主馈线L1故障时故障段前没有可操作隔离开关全部负荷停到修复完成L2故障时S1、S2配合隔离故障段LP2和LP3不经过故障段可以通过切换恢复LP1在L2下游必须等修复L3故障时故障段位于LP1上游所以LP1要等修复LP3位于L3下游分支也要等修复LP2则可以短时切换恢复。每个网络的FMEA表都要根据实际保护配置重新推导这是不能偷懒的步骤。3.4 主循环与指标统计主循环推荐采用每个元件独立推进、时间轴拉通的方式%% 参数设置 rng(2024); % 固定随机种子结果可复现 total_year 100000; % 模拟总年限 total_hours total_year * 8760; % 转换为小时 K length(lambda); % 元件数量 N_LP length(LP_user); % 负荷点数量 % 指标累积变量 N_evt zeros(1, N_LP); % 负荷点停电次数 U_hour zeros(1, N_LP); % 负荷点停电小时数 ENS_kwh zeros(1, N_LP); % 负荷点缺供电量 for k 1:K t 0; % 当前元件的时间轴 while t total_hours % 抽样正常运行时间 ttf_h -log(rand) / lambda(k) * 8760; t t ttf_h; if t total_hours break; end % 抽样修复时间 ttr_h -log(rand) * MTTR_h(k); % 根据FMEA表统计负荷点影响 for i 1:N_LP outage_h FMEA(k, i); U_hour(i) U_hour(i) outage_h; N_evt(i) N_evt(i) 1; ENS_kwh(i) ENS_kwh(i) LP_P(i) * outage_h / 8760; end % 修复完成后继续推进时间 t t ttr_h; end end这个循环的逻辑很直白每个元件在自己的时间轴上不断经历“运行—故障—修复”的循环每发生一次故障就查FMEA表把停电次数、停电小时数、缺供电量累加到对应负荷点。因为元件之间独立推进所以即便某个修复过程跨越了模拟总时长前面已经发生并统计完毕的故障事件也不受影响。指标计算直接套公式%% 负荷点指标 lambda_lp N_evt / total_year; % 次/年 U_lp U_hour / total_year; % 小时/年 r_lp U_lp ./ max(lambda_lp, 1e-12); % 小时/次 %% 系统指标 total_user sum(LP_user); SAIFI sum(lambda_lp .* LP_user) / total_user; SAIDI sum(U_lp .* LP_user) / total_user; CAIDI SAIDI / SAIFI; ASAI (total_year * 8760 * total_user - sum(U_hour .* LP_user)) / ... (total_year * 8760 * total_user); ENS sum(ENS_kwh) / total_year;3.5 收敛性判断不能省略蒙特卡洛模拟最大的陷阱是“看起来跑完了实际上没收敛”。判断收敛最常用的指标是方差系数也就是变异系数。做法是把总模拟分成若干块每块独立计算某个指标得到一串样本然后算样本的标准差和均值看比值是否小于一个阈值。block_num 50; block_year 2000; SAIDI_block zeros(1, block_num); for b 1:block_num % 跑一个block记录该块的SAIDI SAIDI_block(b) simulate_block(block_year); end beta std(SAIDI_block) / mean(SAIDI_block) / sqrt(block_num); if beta 0.05 disp(结果已收敛); end0.05是很多人用的经验阈值要求高就压到0.01。我自己的习惯是先跑小规模测试确定量级再整批跑正式结果既省时间又不至于盲目加大模拟量。4. 典型算例与结果校验4.1 算例系统参数以3.3节那个简化辐射网为例元件参数参考常见配电网可靠性测试系统中的量级设定元件故障率次/年平均修复时间小时断路器 C10.0024馈线区段 L1~L30.105馈线分支 L4~L50.085配电变压器 T10.015200负荷点数据按表设定隔离操作时间取1小时负荷点用户数平均负荷kWLP15030LP23020LP340254.2 运行结果参考量级按上述参数模拟10万年典型结果大致落在以下量级指标参考值SAIFI0.2~0.5 次/年SAIDI1~3 小时/年CAIDI4~8 小时/次ASAI99.95%~99.99%ENS几十到一百多 kWh/年不同网络拓扑、不同隔离策略会让数值明显变化但量级基本一致。如果你的结果偏离这个范围很远优先检查单位换算和FMEA表。4.3 验证实现正确性的三个套路第一个套路是手算简单场景对比。挑一个单负荷点、双元件的小系统用解析法手算理论指标再跑蒙特卡洛对比误差应该在模拟误差范围内。第二个套路是故意放大故障率。把所有元件的故障率同时放大100倍此时故障事件非常频繁蒙特卡洛结果和解析法结果的相对误差会缩小因为事件样本数足够多。如果放大后结果差异还很大说明逻辑有bug而不是抽样误差。第三个套路是检查事件次数。如果元件k的故障率是λ模拟Y年那么预期的故障次数大约是λ×Y次。程序跑完后统计每个元件的故障次数和这个预期值对比偏离太远说明抽样式或循环推进写错了。5. 踩坑记录与提速建议5.1 常见问题速查表写序贯蒙特卡洛的过程中我见过的问题集中在下面几个地方问题现象可能原因解决办法结果每次跑都不一样且差异很大没有固定随机种子或模拟年限不足用rng固定种子检查方差系数停电次数偏少TTF单位错误用了次/小时却没换算统一时间单位建议全用小时SAIDI偏小隔离操作时间没有计入FMEA表FMEA表补充短时停电事件仿真耗时过长循环内重复访问结构体字段用向量存储参数减少字段访问修复时间出现异常大值指数分布长尾本身就存在确认抽样公式是-log(rand)不是log(rand)5.2 经验提速技巧指数分布抽样直接用-log(rand)不要用exprnd。虽然数学上等价但exprnd有额外的函数调用和参数检查开销在几十万次抽样的规模下累积差距很可观。FMEA表尽量在模拟前一次性构建不要在每次故障事件发生时实时搜索网络拓扑。我见过有人用邻接矩阵DFS搜索故障影响范围思路没错但每次故障都搜一遍图模拟速度会慢到让人怀疑人生。稳妥做法是离线构建FMEA表模拟期间只查表。如果网络特别大可以用“批处理随机数”技巧。先一次性生成大矩阵rand(M, N)再按索引取用减少 rand 函数的调用次数。不过这个优化对中小规模网络收益有限代码可读性会下降按需取舍。5.3 重要抽样加速收敛序贯蒙特卡洛的一个明显痛点是配电网元件故障率通常很低0.01到0.1次/年模拟中大量时间消耗在“什么也没发生”的时段。为了加速收敛可以用重要抽样人为把故障率放大k倍让故障事件更频繁地出现统计指标时把贡献除以k作为修正。scale 20; % 放大倍数 lambda_scaled lambda * scale; % 用放大后的故障率抽样 % 统计时停电次数贡献除以scale N_evt(i) N_evt(i) 1 / scale;这个技巧能显著缩短达到相同方差系数的模拟时间代价是引入偏差修正逻辑需要你清楚自己的抽样分布变化。我建议先用无偏的普通模拟跑通验证再引入重要抽样优化大规模算例。6. 一些个人体会把序贯蒙特卡洛在配电网可靠性评估里跑通之后我最大的感受是它其实没有想象中那么高不可攀核心逻辑就是“元件故障—查FMEA表—累计指标”难点反而在参数单位、FMEA构建和收敛判断这些看起来不起眼的细节上。跑通一次之后后续加分布式电源、储能孤岛、保护策略调整都是在现有框架上扩展事件规则整个代码结构不用推倒重来。最后分享一个小技巧做项目汇报时除了给出最终的SAIFI、SAIDI数值建议顺手画出“指标随模拟年限变化的收敛曲线”。一张收敛曲线比任何口头说明都有说服力能直接告诉对方你的结果可信度有多高也能帮你提前发现“跑得不够久”的问题。在做完这些之后你会发现所谓可靠性评估其实等同于一场足够耐心的统计实验。本文还有配套的精品资源点击获取