简介面向电力系统可靠性分析学习者与配电网运维、规划人员这份程序围绕10负荷点配电网的蒙特卡洛可靠性评估展开可用于模拟设备故障与修复过程统计停电时间、停电次数及SAIDI、SAIFI、ASAI等指标帮助理解小型配电网的薄弱环节。压缩包共44个文件约1.07MB包含C工程源码与解决方案、MATLAB计算脚本、可执行程序、Excel网络参数表、PDF结构说明及运行结果截图并带有编译缓存与日志文件便于直接运行、复现实验和二次修改。已有1650人学习下载。读者可据此搭建从故障建模、随机抽样到指标统计的完整流程对照网络结构与参数核对结果也可作为课程设计、科研仿真和运行策略比选的参考。1. 10 个负荷点的辐射状馈线要回答“分段开关挪一格能省多少 SAIDI”一条 10kV 辐射状馈线10 个负荷点、30 个元件评审会上被问把第 5 段的分段开关挪到第 6 段SAIDI 能降多少解析法在这种场景下容易出问题——分段开关、熔断器、联络开关一组合故障后果表就得手写一遍拓扑改一次重算一次改到第三版就没人敢保证表没错。这份「10 负荷点配电网蒙特卡洛法可靠性计算程序」处理的正是这类重复计算把元件台账年故障率、平均修复时间和网络拓扑输入进去用序贯蒙特卡洛抽样跑几千个仿真年输出 10 个负荷点的 λ、U、r再聚合出 SAIFI、SAIDI、CAIDI、ASAI、ENS。适合配网规划、可靠性专责也适合拿它做课程设计或论文复现——代码三百行左右参数全在文件开头看得懂也改得动。2. 元件三参数与序贯抽样λ、r、U 怎么变成一行行时间2.1 从元件台账到负荷点停电序列配电网可靠性里一个元件只需要三个数就能进模型年故障率 λ次/年、平均修复时间 r小时/次、年停电时间 U λ·r小时/年。前两个是台账能查到的第三个是算出来的很多人写程序时把 U 也当输入结果三个数自相矛盾跑出来的 SAIDI 和手算对不上查半天才发现是台账里 U 用的不是同一年的统计口径。蒙特卡洛抽的是时间不是概率。序贯抽样假设元件无故障工作时间服从指数分布即 TTF ~ Exp(λ)反变换抽样就是 TTF -ln(U)/λU 为 (0,1) 均匀随机数。修复时间同样可以抽指数分布也可以直接取常数 r——这一点常被争论实务上对 10 负荷点这种规模修复时间用常数和用指数分布对 SAIFI 几乎没有影响对 ENS 的影响也在 1% 量级因为多重故障重叠的概率本身就低。我一般直接用常数 r少一层随机性结果反而好解释。真正决定结果对不对的不是分布假设而是“一次故障之后哪些负荷点停、停多久”。这部分不归蒙特卡洛管归拓扑和开关配置管第 3 章单独拆。2.2 序贯抽样和非序贯抽样10 负荷点该选哪个非序贯状态抽样是每个元件独立按 p λr/8760 的概率抽“故障/正常”拼出一个全网状态再用后果表算一次序贯状态持续时间抽样是推进时钟一次只处理最近发生的那次故障。两者不是精度高低的问题是输出能力的问题维度非序贯状态抽样序贯状态持续时间抽样直接输出状态概率、期望值频率、持续时间、时序SAIFI / CAIDI需要额外频率公式推导直接统计次数和时长收敛速度快几千次抽样即可慢通常要几千仿真年代码复杂度低一个循环中要维护事件时钟适合场景大系统、只关心 ENS配网指标考核、含开关动作10 个负荷点的规模非序贯一次抽样只要几百微秒但你要的恰恰是 SAIFI 和 CAIDI 这类带“次数”的指标用非序贯得再补一套频率公式补错的地方比省下来的时间多。所以我这边一律用序贯反正 5000 个仿真年在普通笔记本上几十秒就跑完了。提示如果只想验证 ENS或者系统规模上到几百个负荷点非序贯反而是更划算的选择不必死守序贯。2.3 把元件台账写成程序能吃的表先建元件表。主干线段按长度折算故障率分支线同理配变按台数算。10 个负荷点对应 10 段主干线、10 条分支线、10 台配变一共 30 个元件索引顺序固定下来后面后果矩阵直接用这个顺序。import numpy as np # ---- 元件台账全部按你手里的运行数据替换这里是一组典型中压架空网取值 ---- LAMBDA_LINE 0.10 # 主干线故障率, 次/(km·年) LAMBDA_LAT 0.10 # 分支线故障率, 次/(km·年) LAMBDA_TR 0.02 # 配电变压器故障率, 次/(台·年) R_LINE, R_LAT, R_TR 5.0, 3.0, 15.0 # 平均修复时间, h T_SWITCH 0.5 # 分段开关/联络开关倒闸恢复时间, h SEC_LEN [1.0] * 10 # 10 段主干线长度, km LAT_LEN [0.2] * 10 # 10 条分支线长度, km def build_components(): 返回 (名称, 年故障率数组, 修复时间数组)顺序: S1..S10, B1..B10, T1..T10 names, lam, rep [], [], [] for i in range(10): # 主干线段 索引 0-9 names.append(fS{i1}) lam.append(LAMBDA_LINE * SEC_LEN[i]) rep.append(R_LINE) for i in range(10): # 分支线 索引 10-19 names.append(fB{i1}) lam.append(LAMBDA_LAT * LAT_LEN[i]) rep.append(R_LAT) for i in range(10): # 配变 索引 20-29 names.append(fT{i1}) lam.append(LAMBDA_TR) rep.append(R_TR) return names, np.array(lam), np.array(rep)这段的作用是把“物理台账”翻译成“按索引对齐的数组”。三个要点第一线路故障率必须乘长度直接写 0.1 而不乘公里数是新手最容易犯的错10 段线会让全网故障率凭空放大 10 倍第二三类元件的索引区间要写死在注释里后果矩阵按下标取值索引错位是后面最难查的一类 bug第三T_SWITCH单独拎出来它是开关动作时间不是元件修复时间混在一起后面改不动。抽样函数只做一件事把年故障率换成小时制的无故障工作时间def draw_ttf(lam, rng): 按指数分布抽无故障工作时间(h)。lam 单位为次/年内部换算到小时。 ttf np.full(lam.shape, np.inf) ok lam 0 ttf[ok] rng.exponential(1.0 / lam[ok]) * 8760.0 return ttfrng.exponential(scale)的scale是均值单位跟 λ 走——λ 是次/年抽出来就是年乘 8760 换成小时。这一步漏掉乘 8760程序照样跑只是仿真年里几乎不出故障SAIDI 会低到 0.01 小时看着像“电网特别可靠”实际上是单位错了。lam 0的元件返回inf是为了留出“计划停运单独建模”的口子避免除零告警。3. 拓扑离散化与后果矩阵一次故障折成 10 个停电时长3.1 用区域划分代替逐次拓扑遍历辐射状馈线的结构是变电站母线出来S1 到节点 1节点 1 接 LP1 的分支线 B1 和配变 T1再 S2 到节点 2依此类推末端节点 10 通过联络开关与相邻馈线相连。每个节点装一台分段开关每条分支线首端装熔断器。建模的关键动作是把“拓扑”压成一张 30×10 的后果矩阵 RR[i][j] 表示元件 i 故障时负荷点 j 的停电时长小时0 表示不停。这张表一次算好蒙特卡洛主循环里只做查表和累加复杂度从“每次故障遍历一次图”降到 O(1)。10 个负荷点手写都行但一旦要做“加装分段开关”的灵敏度分析拓扑变了表就得重算所以还是写成函数。有人会问这样是不是把蒙特卡洛做成了确定性方法不是。随机性在“哪个元件、什么时候故障、一年故障几次”后果矩阵只是把一次已发生的故障映射成停电时长这部分本来就是确定性的物理过程。分开之后模型可解释性反而强评审时能指着表说“第 5 段故障LP6 到 LP10 靠联络转供0.5 小时恢复”比一句“程序算的”有说服力。3.2 三条恢复判据熔断器、分段开关、联络开关后果矩阵的每一个非零元素都由下面三条判据之一决定写代码前先把判据说清楚能省掉后面一半的调试时间。第一条是熔断器判据。分支线或配变故障熔断器动作只隔离本分支主干线不跳闸所以只有本负荷点停电时长等于该元件的修复时间分支线 3 小时、配变 15 小时。这条判据在矩阵上的表现是“只有对角元非零”。第二条是分段开关判据。主干线第 k 段故障变电站出线断路器跳闸全网瞬时失电随后拉开故障段两侧的分段开关非故障段恢复供电。所以 LP k 自己要等修复5 小时其余 9 个负荷点等倒闸0.5 小时。这就是为什么辐射状馈线的 SAIFI 里每一段线路的故障率都会计到每一个负荷点头上——短暂停电也算一次持续中断。第三条是联络开关判据。故障段下游的负荷点靠联络开关转供恢复前提是联络容量够、且相邻馈线有裕度。容量不足时下游只能等修复停电时长从 0.5 小时跳到 5 小时SAIDI 能差出一倍。程序里用TIE_OK这个布尔量控制做方案比较时把它翻一次结论会很不一样。注意三条判据都建立在“保护正确动作”的前提下。熔断器拒动、分段开关拒动属于另一类问题得用保护配合的可靠性模型单独算别指望这个程序兜住。3.3 后果矩阵的生成代码TIE_OK True # 联络开关能否转供全部下游负荷容量不足时置 False def build_consequence(rep, n_lp10): R[i][j]: 元件 i 故障时负荷点 j 的停电时长(h)0 表示不停电 R np.zeros((len(rep), n_lp)) # 主干线段 k本段等修复其余段等倒闸无联络时下游等修复 for k in range(n_lp): for j in range(n_lp): if j k: R[k, j] rep[k] elif j k and not TIE_OK: R[k, j] rep[k] # 下游无转供只能等修复 else: R[k, j] T_SWITCH # 上游重合/下游转供均为倒闸时间 # 分支线熔断器隔离只影响本负荷点 for k in range(n_lp): R[n_lp k, k] rep[n_lp k] # 配变同上只影响本负荷点 for k in range(n_lp): R[2 * n_lp k, k] rep[2 * n_lp k] return R参数说明rep是第 2 章生成的修复时间数组长度 30n_lp固定 10但写成参数是为了以后扩到 20 个负荷点时不用改循环体TIE_OK是全局开关控制联络开关是否可用。逻辑上有两处值得多说。一是主干线那段用了j k判断上下游这依赖“负荷点编号沿馈线递增”的假设如果你的实际网络有分支再分支得先把节点按距电源的电气距离排序否则上下游判断会错。二是分支线和配变两段故意只写对角元把“熔断器把故障限制在本负荷点”这个工程假设硬编码进去——这是简化现实里熔断器有误动率但把它加进来会让模型复杂度上一个台阶做规划评估时通常不值。4. 序贯抽样主循环与指标聚合一年几千次掷骰子怎么写才不慢4.1 事件推进法每次只推进到最近一次故障主循环的核心是“事件推进”不是“时间步进”。时间步进要按小时或分钟切片8760 小时一年、5000 年就是四千万次判断纯浪费事件推进每次直接跳到下一个故障时刻30 个元件一年平均出 1.4 次故障一次仿真年只有个位数循环快两个数量级。def simulate_one_year(lam, rep, R, rng, horizon8760.0): 仿真一年返回 (各负荷点停电时长, 各负荷点停电次数) n len(lam) ttf draw_ttf(lam, rng) # 每个元件距下次故障的时间 t 0.0 dur np.zeros(R.shape[1]) cnt np.zeros(R.shape[1], dtypeint) while True: i int(np.argmin(ttf)) # 最近发生的故障 t_fail t ttf[i] if t_fail horizon: # 超出仿真年收工 break t t_fail # 查后果矩阵累加停电时长截断到年末避免跨年溢出 d np.minimum(R[i], horizon - t) dur d cnt (R[i] 0).astype(int) # 该元件修复后重新计时修复时间 新的无故障工作时间 ttf[i] rep[i] rng.exponential(1.0 / lam[i]) * 8760.0 return dur, cnt几个参数值得掰开说。horizon固定 8760闰年不用管误差在千分之一量级。np.minimum(R[i], horizon - t)这行是防跨年污染如果第 8759 小时发生故障、修复要 5 小时只能记 1 小时剩下的记到下一年——严格做法是让状态跨年延续但那要维护跨年状态机。我试过两种写法对 5000 年的均值影响小于 0.3%不值得为此把代码复杂度翻倍。cnt (R[i] 0)统计的是“停电次数”只要后果矩阵里该负荷点停电时长大于 0 就计一次跟时长无关这正是 SAIFI 的定义。这个循环有一个已知简化元件 i 处于修复期时其他元件也可能故障这里没有处理重叠。对 10 负荷点、年故障率 1 次量级的系统重叠概率在千分之几对 SAIFI 和 SAIDI 的影响可以忽略对 ENS 略大一点。要严格处理就得维护“当前故障集合”每次故障后重新计算全网状态代码量大约翻三倍——做规划评估没必要做事故后分析就得补上。4.2 指标聚合从负荷点指标到系统指标USERS np.array([120, 90, 150, 80, 200, 110, 60, 140, 95, 130]) # 各负荷点用户数 LOAD_KW np.array([400, 320, 500, 280, 650, 380, 200, 460, 330, 420]) # 各负荷点负荷功率 def year_indices(dur, cnt, usersUSERS, kwLOAD_KW): N, P users.sum(), kw.sum() saifi (cnt * users).sum() / N # 次/(户·年) saidi (dur * users).sum() / N # h/(户·年) caidi saidi / saifi if saifi 0 else 0.0 asai 1.0 - saidi / 8760.0 ens (dur * kw).sum() # kWh/年 return np.array([saifi, saidi, caidi, asai, ens])聚合的关键是“按用户数加权”不是按负荷点平均。10 个负荷点用户数从 60 到 200 不等直接取算术平均会让小负荷点的权重虚高SAIDI 偏差能到 10% 以上。ENS 则是按负荷功率加权量纲是 kWh/年跟售电量口径一致评审时可以直接换算成停电损失。caidi的分母做了零保护虽然正常运行不会为零但参数调成极端值时比如把所有修复时间设成 0会出现加一行判断比事后调试便宜。4.3 收敛判据、仿真年限与随机种子蒙特卡洛的收敛判据用方差系数 β不是“跑够多少年”。β σ / (√n · μ)σ 是各仿真年 SAIDI 的样本标准差μ 是均值n 是仿真年数。β 小于 0.05 意味着均值估计的相对误差在 5% 以内这个口径在配网可靠性评估里够用了。def run(max_year5000, beta_target0.05, seed20240607): rng np.random.default_rng(seed) names, lam, rep build_components() R build_consequence(rep) hist, beta [], np.inf for y in range(1, max_year 1): dur, cnt simulate_one_year(lam, rep, R, rng) hist.append(year_indices(dur, cnt)) if y 50 and y % 20 0: # 前 50 年不判收敛避免早期方差虚低 arr np.array(hist)[:, 1] # 用 SAIDI 判收敛 beta arr.std(ddof1) / np.sqrt(y) / arr.mean() if beta beta_target: break return np.array(hist), y, betaseed固定住是为了可复现——同一份参数两次跑出完全一样的结果评审时才敢让人复核。y 50这个门槛是血泪经验前 20 年样本太少标准差偏小β 会假性达标程序 200 年就停了换一个种子结果差 15%非常被动。y % 20是判定频率每 20 年判一次省掉重复计算对结果无影响。运行完之后hist的每一行是一年的五项指标取均值就是最终结果取标准差配合 β 就是置信区间。5. 避坑与排查五条血泪记录5.1 现象SAIDI 只有 0.01 小时看起来“电网特别可靠”原因几乎一定是单位错位。draw_ttf里 λ 是次/年抽出来的 TTF 单位是年忘了乘 8760 换成小时于是 30 个元件在 8760 小时里几乎不出故障全年零停电。另一种变体是rng.exponential的scale传成了 λ 而不是 1/λ故障频率被彻底搞反。解决在simulate_one_year开头加一句断言assert np.median(ttf[lam 0]) 8760 * 5中位无故障时间超过 5 年就说明单位大概率错了。再加一条结果侧的检查SAIFI 落在 0.5 到 5 次/年之外就要人工看一眼。5.2 现象把倒闸时间设成 0 之后SAIDI 掉了四成但没人发现原因是对“分段开关动作要时间”缺乏体感。非故障段的负荷点确实能恢复供电但恢复不是瞬时的——断路器重合、就地拉开分段开关、合联络开关一套流程下来 0.5 小时是乐观值老旧开关站 1 小时也正常。设成 0 等于假设全网站域自动化且零动作时间SAIDI 偏低 40% 以上。解决T_SWITCH给一个下限 0.5并且明确它同时进 SAIFI 的统计——0.5 小时远大于 3 分钟的中断判定门槛这些短暂停电是要计入 SAIFI 的。如果故意不计得在cnt累加时加时长过滤这是另一个模型假设要写进报告。5.3 现象分支线故障导致全网停电SAIDI 比手算翻倍原因是没建熔断器模型把分支线故障按主干线故障处理了后果矩阵里那一行非对角元被填上了T_SWITCH。分支线故障率虽然低0.02 次/年但 10 条加起来 0.2 次/年每次全网停 0.5 小时摊到每个负荷点就是 0.1 小时/年叠加起来正好把 SAIDI 抬高一大截。解决检查build_consequence里分支线和配变两段循环确认只写了对角元。一个快速自检办法是打印后果矩阵的行稀疏度主干线行 10 个非零分支线和配变行只有 1 个非零对不上就是这里错了。5.4 现象结果和解析法差 8%怎么调参数都对不上原因通常是重叠故障处理方式不同。解析法故障模式后果分析法在线性叠加假设下忽略了“故障叠加”的场景而序贯蒙特卡洛默认允许元件在修复期内被其他故障“叠加”影响——虽然代码里没显式处理重叠但事件推进本身会产生“上一个故障还没修完下一个故障已经发生”的时间窗口后果矩阵照样累加这在物理上是不对的会让 ENS 偏大。反过来如果你的代码严格禁止重叠比如故障期间时钟暂停ENS 又会偏小。两条路子选一条然后在报告里写清楚别指望不同方法能对到小数点后两位。10 负荷点这个规模5% 到 8% 的差异属于正常范围差 30% 才是代码有 bug。5.5 现象换一个随机种子SAIDI 差 20%β 却显示 0.03原因是收敛判据提前触发了。前 50 年样本量小个别年份出现“故障扎堆”或“全年无故障”标准差被低估β 假性达标或者判收敛用的是 SAIFI 而报告里写的是 SAIDI两个指标的方差不同用错判据等于没判。解决三条一起上——前 50 年不判收敛、y % 20而不是每年判、用你最终要报告的指标通常是 SAIDI做判据。另外固定种子跑两遍两次结果必须完全一致不一致说明程序里有没被种子控制住的随机源比如用了全局np.random。6. 结果校核与提速从“能跑”到“能拿去评审”6.1 用解析法做交叉校核误差控制在 1% 以内程序跑通不等于结果可信得先跟手算对一遍。10 负荷点单电源辐射网在“忽略重叠故障”的假设下解析解是闭式的负荷点 j 的年故障率等于所有会影响它的元件故障率之和年停电时间等于每个元件故障率乘对应停电时长再求和。用第 2、3 章那组参数每个负荷点的 λ 和 U 都能手算出来λ 10 × 0.1 0.02 0.02 1.04 次/年U 9 × 0.1 × 0.5 0.1 × 5 0.02 × 3 0.02 × 15 1.31 小时/年。指标解析值蒙特卡洛5000 年相对误差SAIFI次/户·年1.041.04 0.5%SAIDIh/户·年1.311.31 1%CAIDIh/次1.261.26 1%ASAI0.9998500.999850—ENSkWh/年52405200 左右1%2%ENS 的误差比 SAIDI 大一点是正常的它由少数几次长时停电主导方差天然更大。校核顺序建议是先对单负荷点的 λ 和 U再对系统指标——单点对上了系统还对不上问题一定在加权那一步。6.2 提速把后果矩阵和抽样分开10 负荷点的程序其实不缺速度5000 年几十秒。真要做灵敏度扫描比如扫 20 组开关配置每组跑 5000 年瓶颈在 Python 循环而不是算法。两个不改变结果的提速手段一是把build_consequence移出循环配置变了才重算二是把每年 30 个元件的 TTF 预抽成 (max_event, 30) 的矩阵用累积和加argmin一次性定位所有事件时刻避免每次调用rng.exponential。第二招能把单年耗时压到原来的三分之一代价是代码可读性下降我一般只在需要跑上百组参数时才用。再往上就是换语言或者用numba编译主循环但对 10 负荷点这个规模属于过度设计。真正值得优化的是“把一次仿真拆成多个进程跑不同种子”四核并行直接四倍改三行代码。6.3 灵敏度分段开关和联络开关到底值不值得装程序最大的价值不是算出一个 SAIDI 数字是比较方案。把build_consequence里的分段开关数量做成参数——比如原来 10 段 10 台开关改成每两段一台共 5 台重跑一遍SAIDI 会明显上升反过来加密到 15 台SAIDI 下降但边际收益递减。我做过一组对比开关从 5 台加到 10 台SAIDI 降 32%10 台加到 15 台只再降 9%。这条曲线拿去做投资决策比单点指标有用得多。联络开关的灵敏度更值得跑把TIE_OK从True翻成FalseSAIDI 从 1.31 跳到 3.9 左右差了三倍。这就是为什么规划阶段一定要确认联络容量——很多网架图上画了联络开关实际转供能力只够带三成负荷程序里按“全额转供”算出来的结果评审时被问一句就露馅。从那以后我每次交可靠性评估报告都会强制走一遍解析法校核加一组灵敏度对比哪怕委托方只要求一个数字。单点结果谁都能跑出来能说清“这个数在什么条件下会变、变多少”才是这份程序的用法。希望帮到你。本文还有配套的精品资源点击获取
