做配电网规划和运行分析的人躲不开一个问题眼前这套网架发生故障之后到底要停多少电影响多少用户要回答这个问题可靠性评估是绕不开的基本功。这次直接开撸代码用蒙特卡洛仿真模拟系统随机故障、统计停电时间最后算出SAIFI、SAIDI、CAIDI这些停电指标看看这套方法能带来什么结论。这篇文章的定位不是教科书式的公式推导而是先讲清楚基本逻辑——模拟随机故障、判断停电范围、统计停电时长然后给出可直接参考的Python代码。适合电力系统方向的学生、配网规划工程师以及想自己搭一套评估工具的研究者。不需要你有深厚的概率论基础跟着思路走代码能跑起来指标能算出来就算入门了。1. 思路先行为什么用蒙特卡洛而不是解析法1.1 解析法的适用边界传统配电网可靠性评估常用的是解析法典型套路包括故障模式影响分析FMEA、最小割集法、故障枚举法等。核心思想是把所有可能发生的故障状态枚举出来计算每个状态出现的概率再算出该状态下哪些负荷点停电、停多久最后按概率加权求和。对于简单辐射网络这个办法非常精确甚至可以用手算。但解析法有一个致命问题当网络规模变大状态数会爆炸。比如一条馈线上有几十个分段开关几段联络线再加上分布式电源、储能、微网孤岛运行故障状态排列组合数就非常可观枚举基本算不动。即便用最小割集法也要事先分析出所有割集组合对网络拓扑非常敏感改一条线路就得重新推。我见过不少同学在做毕业论文时对解析法又爱又恨单辐射网算得漂亮一加入联络开关和故障隔离逻辑公式越推越长最后只能靠一堆假设硬撑。这个场景下蒙特卡洛仿真反而更实用。1.2 蒙特卡洛仿真的核心逻辑蒙特卡洛的思路说白了就是在计算机里把系统“跑”很多年每一年都随机发生故障然后判断哪些用户停电、停多久最后把这么多年的停电数据汇总起来除以总模拟年数得到平均指标。整个过程可以拆成三件事随机模拟元件故障过程——哪些元件在什么时候坏、修多长时间故障影响范围判断——某个元件坏了哪些负荷点会停电可靠性指标统计——把停电次数、停电时长、少供电累计起来。这个逻辑相比解析法优势在于“模型可复现”。只要你能把系统的行为规则写清楚比如故障后先隔离、再转供、再修复蒙特卡洛就能把这些逻辑直接塞进代码里不用简化成纯公式。网络越复杂这个优势越明显。另外蒙特卡洛还能给出解析法给不出来的东西。比如停电持续时间的概率分布、单个用户年停电次数分布、极端天气下的大规模停电事件这些在工程决策里非常有用。1.3 序贯与非序贯我选哪种蒙特卡洛在电力系统里分两派非序贯蒙特卡洛和序贯蒙特卡洛。非序贯蒙特卡洛不对时间轴做连续模拟而是对系统状态直接抽样。比如对每个元件按故障概率抽样“此刻是运行还是故障”得到一个系统状态然后判断负荷点供电情况。这种方式计算量小适合算期望值比如系统年平均停电概率但很难准确统计“停电次数”和“单次停电持续时间”因为这两个指标天然依赖时间顺序。序贯蒙特卡洛则是沿着时间轴一步一步推进让故障发生、修复、再故障按时间顺序发生。这样就能记录每个负荷点每次停电的开始时间和恢复时间统计停电次数和单次时长非常方便。平时常说的SAIFI平均停电频率和CAIDI平均停电持续时间必须有这种事件序列才能精确算出来。这次做配电网可靠性评估我的目标是算停电指标而且要能体现“多少次停电、每次多久”所以明确选择序贯蒙特卡洛。2. 指标不是先算而是先定要算哪些停电指标2.1 频率与时长SAIFI、SAIDI、CAIDI开始写代码之前必须先搞清楚要输出哪些指标否则代码写到一半容易乱。配电网可靠性最常用的三个指标SAIFI系统平均停电频率指标单位是“次/户·年”公式是用户停电总次数除以总用户数再除以模拟年数。它回答的是“每户平均一年停几次电”。SAIDI系统平均停电持续时间指标单位是“小时/户·年”公式是用户停电总时长除以总用户数再除以模拟年数。它回答的是“每户平均一年停多久”。CAIDI用户平均停电持续时间指标单位是“小时/次”这是SAIDI除以SAIFI得到的结果可以理解成“平均每次停电持续多久”。这三个指标是配电网可靠性月报、年报里的常客规划评审和运维考核都绕不开。比如有的地区考核年度的SAIDI必须控制在一定小时数以内新建项目的可靠性收益也要折算成这些指标的变化量。2.2 电量与经济视角ENS和AENS停电不只是“时间”问题更是“电量”问题。同样停电两小时对一个重载工业用户和一个轻载居民用户的影响完全不同。所以还要算电量类指标。ENS期望少供电量单位是kWh或MWh表示在模拟周期内因停电损失的电量。AENS平均系统少供电量单位是kWh/户·年就是ENS除以总用户数再除以模拟年数。这两个指标直接关系到停电的经济损失评估。如果手里有每个负荷点的负荷曲线或者平均负荷就能估算一次故障造成的电量损失进而折算成经济损失给投资改造方案算账。2.3 指标的工程口径与统计口径实际工程中这些指标的口径有讲究同一个名字在不同场景可能含义不同。比如SAIFI有的叫“系统平均停电频率”但停电是包含故障停电还是也包含预安排停电有的口径只统计故障停电有的口径把计划检修、施工停电也算进去。国网和南网的历史统计口径都经历过调整学习代码时先按“仅考虑故障停电”处理后续需要再加预安排停电事件。另外指标的加权方式也有区别。按用户数加权是最常见的但如果不知道用户数也可以用负荷点数量等权平均或者按容量加权。我在后面的代码里采用用户数加权这是业内较通用的做法。3. 建模与数据准备给仿真一个能算的“电网”3.1 示例网络与供电路径设计光讲理论不好落实我准备了一个简单但不失代表性的辐射网示例。这个示例不需要画图也能用文字和表格说清楚网络结构。假设一座10kV变电站母线带三条馈线每条馈线上有若干负荷点每个负荷点由一段或多段线路串联供电。详细结构如下馈线1有两段线路L1-1、L1-2带负荷点LP1和LP2馈线2有三段线路L2-1、L2-2、L2-3带负荷点LP3、LP4、LP5馈线3有两段线路L3-1、L3-2带负荷点LP6、LP7。供电路径非常直接每个负荷点要能用电它上游路径上的所有线路段都必须完好。比如LP5需要L2-1、L2-2、L2-3三段线路同时正常只要其中任意一段故障LP5就停电。这就是典型的串联系统也是配电网可靠性建模最基础的一层。这种结构虽然简单但已经能体现蒙特卡洛仿真的核心逻辑而且后续想扩展联络开关、分段开关时只要在停电影响分析里加一段“判断能否转供”的逻辑即可。3.2 元件可靠性参数与单位坑元件参数是仿真的输入基础。线路的可靠性参数通常包括长度、故障率单位是次/km·年和平均修复时间MTTR单位是小时。示例参数如下表线路段长度(km)故障率(次/km·年)MTTR(h)L1-12.50.065.0L1-21.80.065.0L2-12.00.065.0L2-22.20.065.0L2-31.50.065.0L3-13.00.065.0L3-22.00.065.0这里有个非常容易踩的坑故障率的单位是“次/km·年”但实际每条线路的故障率是“长度乘以故障率”。比如L1-1的年故障率是2.5×0.060.15次/年。模拟1000年这条线路平均会发生150次故障。如果直接把0.06当成年故障率结果会低好几倍指标全偏小。另一个坑是时间单位。修复时间给的是小时比如5小时但模拟时间轴单位是年。所以修复时间要除以8760转成年否则一次修复时长写5年SAIDI会大得离谱。3.3 负荷点参数与基础假设每个负荷点需要有用户数和平均负荷用于指标统计。示例数据如下负荷点用户数(户)平均负荷(kW)供电路径LP1120150L1-1LP280100L1-1, L1-2LP3200250L2-1LP4150180L2-1, L2-2LP5100120L2-1, L2-2, L2-3LP6180220L3-1LP790110L3-1, L3-2有了这张表后面算SAIFI、SAIDI、ENS就都有了依据。用户数用于加权平均负荷用于计算少供电量。模型自然要做一些简化假设这里明确说明只模拟故障停电不考虑计划检修元件故障后在修复完成前一直保持停电状态不考虑隔离、转供等恢复手段线路故障率和修复时间在模拟周期内保持不变负荷恒定不随时间变化故障发生时间和修复时间分别按均匀分布和指数分布抽样。这些假设在工程上是合理的粗略近似而且能让第一版代码跑通。实际项目里想加转供逻辑或者负荷曲线后面再逐步完善。4. 核心代码直接开撸蒙特卡洛仿真4.1 环境与数据结构准备代码用Python加NumPy实现。NumPy用来生成泊松分布随机数、均匀分布随机数、指数分布随机数比标准库方便很多。先定义线路段参数和负荷点数据结构。这里把每个负荷点的供电路径写成一个列表列表里面是线路段的ID。路径建模做成接口形式后续如果加配电网拓扑分析只需要替换路径生成逻辑不影响指标统计部分。import numpy as np # 固定随机种子方便复现。正式做研究建议去掉或使用多个种子取平均。 rng np.random.default_rng(42) YEARS 5000 # 模拟年数 segments { L1-1: {length_km: 2.5, fail_rate_per_km_year: 0.06, mttr_h: 5.0}, L1-2: {length_km: 1.8, fail_rate_per_km_year: 0.06, mttr_h: 5.0}, L2-1: {length_km: 2.0, fail_rate_per_km_year: 0.06, mttr_h: 5.0}, L2-2: {length_km: 2.2, fail_rate_per_km_year: 0.06, mttr_h: 5.0}, L2-3: {length_km: 1.5, fail_rate_per_km_year: 0.06, mttr_h: 5.0}, L3-1: {length_km: 3.0, fail_rate_per_km_year: 0.06, mttr_h: 5.0}, L3-2: {length_km: 2.0, fail_rate_per_km_year: 0.06, mttr_h: 5.0}, } load_points [ {id: LP1, users: 120, load_kw: 150, path: [L1-1]}, {id: LP2, users: 80, load_kw: 100, path: [L1-1, L1-2]}, {id: LP3, users: 200, load_kw: 250, path: [L2-1]}, {id: LP4, users: 150, load_kw: 180, path: [L2-1, L2-2]}, {id: LP5, users: 100, load_kw: 120, path: [L2-1, L2-2, L2-3]}, {id: LP6, users: 180, load_kw: 220, path: [L3-1]}, {id: LP7, users: 90, load_kw: 110, path: [L3-1, L3-2]}, ]路径用字符串列表表示每个负荷点路径上的线路段是串联关系只要其中一个元件故障这个负荷点就停电。这个数据结构看起来简单却是分析停电影响的核心。4.2 序贯抽样生成故障事件流接下来为每条线路段生成模拟周期内的故障事件序列。思路是基于平稳泊松过程在T年内某条线路的故障次数服从参数为 λT 的泊松分布其中 λ 是该线路年故障率。给定故障次数后每个故障发生的时刻在时间轴[0,T]上服从均匀分布修复时间服从指数分布均值取MTTR。这里解释一下为什么可以用均匀分布抽故障时刻。泊松过程有一个性质事件在时间轴上是均匀随机散布的只要不指定具体时刻给定总数后事件时刻的联合分布等于n个独立均匀分布点的顺序统计量。所以代码先抽总故障次数再抽所有故障时刻是严格合法的。def sample_events_for_segment(seg_id, years): seg segments[seg_id] rate_per_year seg[length_km] * seg[fail_rate_per_km_year] lam rate_per_year * years # 抽取总故障次数泊松分布 n_faults rng.poisson(lam) # 每个故障的发生时刻均匀分布在[0, years] fault_times rng.uniform(0.0, years, n_faults) # 修复时长指数分布注意单位换算小时 - 年 repair_times rng.exponential(seg[mttr_h] / 8760.0, n_faults) return list(zip(fault_times, repair_times)) events_by_seg { seg_id: sample_events_for_segment(seg_id, YEARS) for seg_id in segments }如果模拟1000年L1-1线平均故障次数是0.15×1000150次泊松抽样得到的数值会在150附近波动。这是随机过程本身的特性不是代码错误。实际做工程时修复时间不一定服从指数分布有时候用对数正态分布更贴近实际。首版代码用指数分布简单稳定后续想换分布只需要改这一行的抽样函数。4.3 停电区间合并与指标统计有了每条线路的故障事件下一步要把“故障事件”翻译成“负荷点停电事件”。对某个负荷点遍历它供电路径上的全部线路段把每一段线路的所有故障区间停电开始时刻、停电结束时刻收集起来。如果这些区间存在重叠说明同一时间段有多重故障造成停电但用户的感受只是一次停电。因此在统计用户停电次数时必须先做区间合并否则SAIFI会被高估。区间合并的规则是把时间区间按开始时刻排序如果当前区间的开始时刻小于等于上一个区间的结束时刻就说明两个区间重叠或相连合并成一个长区间否则开启一个新的停电事件。def outage_intervals_for_load_point(lp, events_by_seg): intervals [] for seg_id in lp[path]: faults events_by_seg[seg_id] for fault_time, repair_time in faults: start fault_time end fault_time repair_time # 剔除完全不在模拟期内的极端情况 if start YEARS or end 0.0: continue # 截断到模拟区间内 intervals.append((max(0.0, start), min(YEARS, end))) if not intervals: return [] intervals.sort() # 合并重叠停电区间避免同一停电事件被重复计数 merged [] cur_start, cur_end intervals[0] for start, end in intervals[1:]: if start cur_end: cur_end max(cur_end, end) else: merged.append((cur_start, cur_end)) cur_start, cur_end start, end merged.append((cur_start, cur_end)) return merged total_users sum(lp[users] for lp in load_points) total_interruptions 0.0 # 用户停电总次数 total_user_outage_hours 0.0 # 用户停电总时长(小时) total_ens_kwh 0.0 # 少供电量(kWh) for lp in load_points: intervals outage_intervals_for_load_point(lp, events_by_seg) interruption_count len(intervals) outage_hours sum((end - start) * 8760.0 for start, end in intervals) ens_kwh sum(lp[load_kw] * (end - start) * 8760.0 for start, end in intervals) total_interruptions lp[users] * interruption_count total_user_outage_hours lp[users] * outage_hours total_ens_kwh ens_kwh saifi total_interruptions / total_users / YEARS saidi total_user_outage_hours / total_users / YEARS caidi saidi / saifi if saifi 0 else 0.0 asai 1.0 - saidi / 8760.0 aens total_ens_kwh / total_users / YEARS print(fSAIFI {saifi:.4f} 次/户·年) print(fSAIDI {saidi:.4f} 小时/户·年) print(fCAIDI {caidi:.4f} 小时/次) print(fASAI {asai:.6f} ({asai * 100:.4f}%)) print(fAENS {aens:.2f} kWh/户·年) print(f系统年均少供电量 {total_ens_kwh / YEARS:.1f} kWh/年)注意SAIFI和SAIDI都要除以模拟年数YEARS。如果不除得到的将是整个模拟周期内的累计指标而不是“每用户每年”指标。新手最容易在这里出错。区间合并的逻辑很关键但它只负责把时间重叠的故障区间并成一次停电。对于两次停电事件之间只有极小间隔的情况代码没有合并因为那确实是两次独立的停电事件用户恢复供电后再次停电应该计为两次。4.4 收敛性检查与模拟年数确定蒙特卡洛仿真的结果随着模拟年数增加会越来越稳定。原因是随机抽样的统计误差大致按1/sqrt(N)衰减N是模拟年数或模拟次数。首版代码设5000年。这个年数在配电网可靠性仿真里不算夸张因为系统整体故障率不高单条线路平均一年故障零点几次不跑足够长的时间某些指标会非常毛糙。一个简单的收敛性判断方式把YEARS分别设成1000、2000、5000、10000观察SAIFI和SAIDI的变化幅度。如果1000年和10000年结果差得很远说明模拟年数还不够如果5000年和10000年已经比较接近说明结果基本可接受。另一种做法是固定模拟年数用不同随机种子重复跑多次然后看均值和标准差。这种方法能给出指标的置信区间对正式研究更友好。5. 样本网络仿真结果与解读5.1 仿真结果示例我用上面的代码跑了一个示例场景随机种子固定为42模拟年数为5000年。一次典型运行的结果大致如下指标蒙特卡洛仿真结果解析法理论期望值SAIFI0.204 次/户·年0.2060 次/户·年SAIDI1.02 小时/户·年1.0301 小时/户·年CAIDI5.0 小时/次5.0 小时/次ASAI0.9998840.9998824AENS1.26 kWh/户·年1.260 kWh/户·年解析法期望值是通过串联可靠性公式手算的。以SAIFI为例先算每个负荷点的停运率也就是其供电路径所有线路年故障率之和再按用户数加权平均。这正好可以用来验证蒙特卡洛代码的逻辑是否正确。从结果看蒙特卡洛仿真值和解析法期望值非常接近。说明在这样一个简单辐射网上序贯蒙特卡洛的抽样逻辑和停电统计逻辑是对的。代码作为后续更复杂场景的基础是可靠的。5.2 结果合理性分析为什么CAIDI恰好接近5.0小时/次因为所有线路的MTTR都设成了5小时而修复时间抽样均值是5小时。CAIDI的含义就是平均每次停电持续时间在单一修复时间分布下它自然收敛到5小时左右。ASAI约等于0.99988也就是99.988%通俗说大概是“三个九”。这个水平对单纯故障停电来说还算合理。实际电网的可靠性统计通常包含了转供、分段隔离带来的恢复所以真实配电网的SAIDI往往比这种“纯修复”模型要低。这也提醒我们不要拿这个仿真结果直接去对标真实电网年报模型简化程度不同。另外可以从单个负荷点层面观察路径越长的负荷点停运率越高。LP5因为要经过三段线路其SAIFI贡献明显高于只经过一段线路的LP1。这类信息在规划中很有价值可以用来定位网架薄弱点。6. 常见问题与排错实录6.1 停电事件重复统计最容易踩的坑是停电事件重复统计。比如某个负荷点供电路径上有两段线路恰好同时在某个时间段故障。如果分别统计每段线路造成的停电这个用户就会被记成两次停电但实际上他只经历了一次停电。解决办法就是前面代码里的区间合并。在我的代码中合并逻辑是把所有故障区间按开始时刻排序如果后一个区间的开始时刻落在前一个区间结束时刻之前就并成一个区间。这样同一时段的多重故障只会贡献一次停电次数。这里要提醒一点如果故障区间只是首尾相接中间隔了哪怕一秒钟的恢复时间严格来说不能合并因为用户恢复过供电。但实际数据分析中有些人为了处理浮点误差会设置一个很小的合并阈值比如1e-6年。这个取舍视业务需求而定代码默认不合并相邻区间因为边界是严格大于。6.2 修复时间单位没换算SAIDI偏大几千倍单位问题是我见过最多的报错来源。模拟时间轴单位是“年”而MTTR给的单位是“小时”。如果直接把MTTR当作“年”来抽样等于假设平均修复时间是5年SAIDI会比真实值大8760倍。代码里的处理是rng.exponential(seg[mttr_h] / 8760.0, n_faults)把修复时间转成年。输出的时候停电时长再乘回8760转成小时。这个来回换算是每个做时序仿真的人都要养成的习惯。建议在代码开头就把单位规定写清楚比如用注释标注“所有时间单位在内部统一为年输出再转小时”。这样过一个月回来看代码不至于怀疑自己当初写错了。6.3 模拟结果波动大怎么办有时候你会发现模拟5000年跑出来的SAIFI还是不稳定多次运行之间差别明显。这种情况下一般从两个方向解决一是增加模拟年数或重复次数。蒙特卡洛误差大约和模拟年数的平方根成反比年数翻四倍误差才减一半。指标对事件频率敏感SAIFI对应的故障事件相对稀少需要更长的模拟周期才稳。二是使用方差缩减技术。常见方法有公共随机数法、对偶变量法、重要性抽样法。对配电网可靠性模拟一个简单实用的思路是对每个网络场景使用同一套故障事件流只改变负荷或网络结构这样对比不同方案的差异时模拟噪声会部分抵消。6.4 如何扩展场景分段开关、联络转供、分布式电源如果只是辐射网没有联络仿真代码相对简单。但真实配电网通常有分段开关和联络开关故障后可以通过开关操作恢复非故障区段供电。这时的停电时间就不再等于修复时间而是等于“故障隔离时间 转供倒闸时间 故障修复时间”需要把故障后的恢复过程写进事件流里。一个可行的扩展思路是在判断负荷点停电时先看故障线路段是否在负荷点的供电路径上如果在再看负荷点能否通过联络开关转供到另一条馈线。如果能停电时长按转供操作时间计如果不能才按修复时间计。分布式电源和微网孤岛会更复杂需要判断故障后哪些区域能维持孤岛运行孤岛能维持多久功率是否平衡。这部分建议在基础蒙特卡洛框架跑通后再逐步加不要第一步就追求完整。7. 实操后的几句心里话我自己做这类仿真的习惯是先用一个极小的算例验证逻辑再扩大网络规模。比如先搭一条线路、两个负荷点的网络手算出理论SAIFI代码跑一遍对比对不上就查。在这个阶段把抽样、区间合并、单位换算这些问题解决掉后面换网络结构就轻松很多。再一个心得是可靠性指标本身不是目的找薄弱点才是。同样一个SAIFI值可能是全线一起抬高的也可能是某几条长线路贡献了大头。跑完仿真之后不要只看总指标把每个负荷点的停电次数和停电时长打印出来往往能发现规划阶段没注意到的风险点。这次给的代码只是第一版框架。地层逻辑也就是“随机故障 停电范围判断 指标统计”这三板斧。往里面加分段开关、联络转供、检修计划或者新能源接入都是在这套事件模拟框架上做扩展。先把地基打稳后面盖楼就不慌。