ARTICLE · INTELLIGENCE

战地情报 · 详情页

来自尧图项目组的一线实战观察与深度解析

蒙特卡洛法评估配电网可靠性:指标定义与MATLAB实现

蒙特卡洛法评估配电网可靠性:指标定义与MATLAB实现

简介:面向配电网可靠性评估研究人员与电气工程学习者的一份技术资料,聚焦基于蒙特卡洛法实现配电网可靠性评估。内容从传统确定性评估方法的局限切入,系统讲解蒙特卡洛随机模拟的基本原理,给出确定评估目标与范围、识别随机变量及其概率分布、多次抽样仿真、统计计算可靠性指标的实施步骤,并分析该方法在考虑随机性因素、支持不同配电网方案评估方面的优势,以及计算量大、依赖概率分布精度等不足。资源为PDF格式,共1个文件,压缩包整体约174KB,目前已有20人学习下载。读者通过该文档可掌握蒙特卡洛法应用于配电网可靠性评估的完整理论框架,理解抽样次数与仿真模型精度对结果的影响,为相关课题研究或课程设计提供方法参考与代码实现指引。

1. 蒙特卡洛法评估配电网可靠性:当解析法状态爆炸时,模拟法为什么更实用

配电网可靠性评估,传统做法是解析法,把网络拆成最小割集然后逐条推导。网架一复杂,分段开关、联络开关、分布式电源和负荷转供策略叠加进来,状态空间会膨胀到无法穷举,公式推导很容易变成一场灾难。蒙特卡洛法换了个思路:不追求枚举所有状态,而是按元件的故障规律做大量随机试验,用统计平均逼近真实期望。它特别适合配电网可靠性评估这类“事件稀疏、影响耦合”的问题——哪怕单条馈线故障率只有 0.1 次/年,只要抽样次数够多,故障组合和负荷点停电时长的统计特征就会浮现出来。这篇文章从指标定义、元件建模、MATLAB 代码、抽样收敛到工程坑位,把基于蒙特卡洛法的配电网可靠性评估完整走一遍,适合正在做配电网规划、可靠性指标核算或相关课题仿真的人直接参考。

2. 可靠性指标与元件状态模型:评估前先把尺子和数据备齐

蒙特卡洛只是个计算引擎,装什么指标、用什么元件模型,直接决定输出有没有工程意义。不少初学者一上来就写抽样循环,跑了十万次却说不清 SAIFI 和 SAIDI 的口径差异,结果表根本没法用。这一章先把评估的“尺子”和“原料”定下来。

2.1 七个必算指标:SAIFI、SAIDI、ASAI 的统计口径一张表说清

配电网可靠性评估里,最常用的指标是系统平均停电频率指标(SAIFI)、系统平均停电持续时间指标(SAIDI)、供电可靠率(ASAI),以及负荷侧的缺供电量(ENS)。它们不是并列关系:SAIFI 看“停多少次”,SAIDI 看“停多久”,ASAI 看“全年供电时间的百分比”,ENS 看“损失了多少电量”。四个指标各有侧重,核价、考核、规划选型时看的维度完全不同,建议评估模型一次性全部输出。

指标计算公式物理含义常用口径
SAIFI用户总停电次数 / 总用户数每个用户年均停电几次次/户·年
SAIDI用户总停电时长 / 总用户数每个用户年均停电几小时小时/户·年
CAIFI停电用户总停电次数 / 停电用户数受影响用户平均停几次次/户·年
CAIDI停电用户总停电时长 / 停电用户总次数每次停电平均持续多久小时/次
ASAI用户总供电时间 / 用户总需求时间供电可用率无量纲(常写为%)
ENS负荷点缺供电量之和全年总缺供电量MWh/年
AENSENS / 总用户数每个用户年均缺供电量MWh/户·年

这里有三个容易混的点。第一,CAIFI 和 SAIFI 的分母不同,前者只统计实际被停电的用户,后者统计系统全部用户,所以 CAIFI 通常会大于 SAIFI。第二,CAIDI 在数值上等于 SAIDI / SAIFI,但两个指标来自不同的统计路径,不能在程序里直接拿除法结果冒充 CAIDI,要用原始停电事件逐次累加。第三,ASAI 的口径是“用户总供电小时数”除以“用户总需求小时数”,算年指标时分母是用户数乘 8760 小时,不是简单把 SAIDI 补成 1 减。我建议在代码里分别维护“停电用户数”、“停电用户时数”、“缺供电量”三个累加器,最后统一换算,这样口径不会串。

2.2 二状态与三状态模型:从架空线到配电变压器的状态概率

配电网元件在可靠性评估里通常先建模成二状态:正常运行和故障停运。每个元件有两个关键参数——故障率 λ(次/年)和平均修复时间 r(小时/次)。由此可以算出稳态不可用度 U = λ×r / 8760,这个无量纲概率是后续蒙特卡洛抽样的核心输入,也是一切状态概率的起点。

实际工程里,二状态模型往往不够用。配电变压器有周期性检修,检修期间算不算停运?架空线路在恶劣天气下故障率会上升 10 倍以上,还和季节强相关。这些场景建议扩展成三状态模型:运行、故障停运、计划检修。计划检修状态用检修频率 f 和平均检修时长 s 描述,不可用度变成 (λ×r + f×s) / 8760。注意这里故障率和检修频率不能简单相加,因为两个事件的触发机制不同,前者是随机冲击,后者是计划安排,仿真时需要用不同的随机数流分别抽样。

分布式电源接入后,状态模型还要再加一层。光伏和风电出力不是简单的“可用/不可用”,而是多级出力状态。常见做法是把出力曲线离散成若干区间,每个区间对应一个出力概率,例如把一天 24 小时按光照分成 4 个典型时段,每个时段再按出力水平分成 3 到 5 档。这一层模型不会直接参与元件故障抽样,但会在潮流计算和负荷转供判断里决定分布式电源能否支撑部分负荷。

2.3 数据清单:一份能直接喂给仿真器的可靠性基础数据表

没有干净的数据,再精妙的蒙特卡洛算法也只是在拿噪声做统计。我一般会在写代码之前先整理一张元件参数表,每一行代表一个元件,字段尽量标准化。工程上主馈线的故障率通常在 0.05 到 0.2 次/(km·年)之间,配电变压器的故障率在 0.005 到 0.015 次/年左右,平均修复时间在 2 到 8 小时之间,具体取值以当地供电公司的历史统计为准。

字段名单位典型取值区间说明
id无1~N元件唯一编号,不能重复
type无feeder/transformer/switch/bus元件类型,决定抽样方式
from_node无0~N首端节点编号
to_node无0~N末端节点编号
lengthkm0.5~15仅线路需要,用于折算故障率
lambda次/年0.005~0.2年故障率,需按长度折算
r_repair小时/次2~8平均修复时间
f_maintain次/年0~2计划检修频率,三状态模型用
s_maintain小时/次4~24平均检修时长,三状态模型用
capacityMVA0.1~20容量限值,转供和潮流判断用
active_loadMW0~10负荷节点有功,非负荷节点填 0

2.4 拓扑编码:从节点支路表到供电路径的预备工作

配电网通常是辐射状结构,可靠性评估里最关键的一步是判断“某个元件故障后,哪些负荷点失电”。这依赖网络拓扑信息。最常见的编码方式是用节点-支路表描述网络,每个节点有编号,每条支路记录首末端节点号,然后从电源节点出发做深度优先遍历,生成每个负荷节点的供电路径。供电路径是一条从电源到负荷点的节点序列,路径上任一元件故障,都会导致该负荷点停运(除非有联络开关可以倒供)。

这段预备工作在写蒙特卡洛主循环之前就必须完成。常见做法是预先为每个负荷节点算好“上游支路集合”,存入结构体数组,主循环里只需要查这个集合是否与故障元件集合有交集,就能判断负荷点是否失电。这样做能把每次试验的判断耗时从网络规模 O(N) 降到几乎 O(1),在抽样几万次时差异非常明显。

3. 用 MATLAB 实现最小可复现的蒙特卡洛评估框架:从元件抽样到指标输出

这一章给出一个能直接跑通的最小框架。代码覆盖网络建模、状态抽样、停电分析、指标统计四个环节,读者可以在这个骨架上替换自己的网络参数。

3.1 网络建模:用结构体数组承载节点和支路信息

MATLAB 里最直观的网络建模方式是用结构体数组,每个支路一个元素。下面这段代码定义了一个 3 支路的辐射网,电源节点为 0,负荷节点 2 和 3 分别挂在不同支路末端。

% 支路结构体定义:每个元素是一条支路 % lambda 和 r_repair 是可靠性参数,capacity 是容量限值 clear; clc; br(1) = struct('id', 1, 'from', 0, 'to', 1, 'lambda', 0.10, ... 'r_repair', 5.0, 'capacity', 5.0, 'active_load', 0.0); br(2) = struct('id', 2, 'from', 1, 'to', 2, 'lambda', 0.08, ... 'r_repair', 4.0, 'capacity', 3.0, 'active_load', 1.5); br(3) = struct('id', 3, 'from', 1, 'to', 3, 'lambda', 0.06, ... 'r_repair', 4.5, 'capacity', 2.5, 'active_load', 1.0); % 负荷节点: 节点2有200户, 节点3有150户 cust(2) = 200; cust(3) = 150; cust(1) = 0;

这段代码的结构体里,“from”和“to”是支路两端节点号,“active_load”在这个例子里直接挂在了支路末端节点上。实际工程中负荷数据通常单独存在节点表里,这里为了演示方便直接内嵌在支路中。注意“lambda”和“r_repair”的单位分别是次/年和小时/次,后面计算不可用度时必须除以 8760。

3.2 非时序蒙特卡洛主循环:按不可用度抽样并累计停电统计量

非时序蒙特卡洛的核心是用随机数判定每个元件在“某一时刻”是否处于故障状态。判定依据是稳态不可用度 U = λ×r / 8760。以下代码实现了完整的抽样-判断-统计循环。

% 非时序蒙特卡洛主循环 % 对每个元件生成服从均匀分布的随机数,小于不可用度则判为故障 rng(42); % 固定随机数种子,保证结果可复现 N_sim = 20000; % 抽样次数,工程上建议从 2 万次起步 n_br = numel(br); % 累加器初始化 sum_cust_interrupted = 0; % 用户停电次数累加 sum_cust_hours = 0; % 用户停电时数累加 sum_energy_loss = 0; % 缺供电量累加 for k = 1:N_sim % 步骤1: 元件状态抽样 faulty = []; for i = 1:n_br U = br(i).lambda * br(i).r_repair / 8760; if rand() < U faulty(end+1) = br(i).id; end end % 步骤2: 如果无故障元件,本次试验不停电,直接跳过 if isempty(faulty) continue; end % 步骤3: 对每个负荷节点判断是否受故障影响 % demo_network_int 返回停电用户数、停电时数和缺供电量 [c_off, h_off, e_off] = demo_network_int(br, faulty, cust); sum_cust_interrupted = sum_cust_interrupted + c_off; sum_cust_hours = sum_cust_hours + h_off; sum_energy_loss = sum_energy_loss + e_off; end % 步骤4: 折算成年指标 total_cust = sum(cust(2:end)); % 用户总数(不含电源节点) SAIFI = sum_cust_interrupted / total_cust / N_sim; SAIDI = sum_cust_hours / total_cust / N_sim; ASAI = 1 - sum_cust_hours / (total_cust * 8760 * N_sim); ENS = sum_energy_loss / N_sim; fprintf('SAIFI = %.4f 次/户·年\n', SAIFI); fprintf('SAIDI = %.4f 小时/户·年\n', SAIDI); fprintf('ASAI = %.6f\n', ASAI); fprintf('ENS = %.4f MWh/年\n', ENS);

代码逻辑分四步。第一步对每个元件做一次伯努利试验,故障概率取不可用度,这里最容易犯的错是把 λ 直接当概率用——如果那样,故障率 0.1 的元件会有 10% 的试验都处于故障状态,SAIDI 会放大几万倍。第二步把无故障试验直接跳过,因为无故障不会产生任何停电统计量,这样可以大幅节省计算时间。第三步调用一个网络影响分析函数,传入故障元件集合和用户数据,返回三个累计量。第四步把累计量除以用户数和抽样次数,折算成年指标。

3.3 影响分析函数:如何判断负荷点是否失电并计算停电时长

下面的demo_network_int函数是影响分析的核心:它读取故障元件集合,对每个负荷节点检查供电路径上是否有故障元件,再按是否有联络转供能力计算停电时长。

function [c_off, h_off, e_off] = demo_network_int(br, faulty, cust) % 预定义每个负荷节点的供电路径(按网络结构手工指定或由遍历生成) % 节点2的供电路径: 支路1,2;节点3的供电路径: 支路1,3 path_cell = {[1, 2], [1, 3]}; % 负荷节点2和3的上游支路 load_mw = [1.5, 1.0]; % 节点2和3的有功负荷 c_off = 0; h_off = 0; e_off = 0; for idx = 1:numel(path_cell) path = path_cell{idx}; if isempty(intersect(path, faulty)) continue; % 供电路径上没有故障元件,不影响该负荷点 end % 停电用户数: 该节点用户全部计入 n_cust = cust(idx+1); c_off = c_off + n_cust; % 停电时长: 若无转供能力,按最长修复时间计 % 简化模型: 取故障元件修复时间的最大值 r_max = 0; for f = faulty idx_f = find([br.id] == f); if ~isempty(idx_f) r_max = max(r_max, br(idx_f).r_repair); end end h_off = h_off + n_cust * r_max; % 缺供电量: 停电时长 * 负荷功率 e_off = e_off + load_mw(idx) * r_max; end end

这个函数演示的是不带联络开关的最简场景:供电路径上有故障元件就停电,停电时长按故障修复时间算。真实场景中,分段开关可以在修复前先隔离故障,联络开关可以从相邻馈线倒供非故障段,所以停电时长通常不等于修复时间,而是“隔离操作时间 + 倒供时间 + 不可转供部分的修复时间”。这部分逻辑需要在影响分析函数里根据网络自动化水平细化,后续章节会展开讲相关坑位。

3.4 序贯蒙特卡洛:需要时变负荷和检修计划时的时序推进版本

如果评估对象含时变负荷曲线、分布式电源日出力或计划检修窗口,非时序方法就不够用了,需要序贯蒙特卡洛——也就是逐事件推进时间轴。对每个元件用指数分布抽样故障时刻 ttf,再抽样修复时长 ttr,把所有元件的事件按时间排序,逐年推进并累计停电指标。

% 序贯蒙特卡洛单次试验示例 % TTF 抽样: exprnd(1/lambda), TTR 抽样: exprnd(r_repair) T_year = 8760; t_now = 0; sum_h = 0; while t_now < T_year % 对每个元件抽样下一次故障时刻和修复时长 % 实际代码需要遍历所有元件,这里只示意单元件 ttf = exprnd(1 / br(1).lambda); ttr = exprnd(br(1).r_repair); % 如果故障事件落在年内,累计停电时长 if t_now + ttf < T_year sum_h = sum_h + min(ttr, T_year - (t_now + ttf)); end t_now = t_now + ttf + ttr; end

代码里exprnd(1/lambda)生成平均值为 1/lambda 的指数分布随机数,代表“从当前时刻到下一次故障的时间”。注意 MATLAB 的exprnd参数是均值而不是率,所以必须写成1/lambda,写成exprnd(lambda)会让时间尺度完全错误。序贯法的统计量可以直接映射到日历时间,负荷时变曲线也可以按时间索引插值,这是它能做精细评估的根本原因。代价是计算量比非时序高一个量级,尤其当元件数量超过 100 条支路时,需要配合事件驱动和稀疏化技巧。

4. 抽样次数、收敛判据与参数整定:仿真结果不抖才算真收敛

蒙特卡洛法的结果是一个随机变量,不是确定值。同样是 2 万次抽样,两次独立运行的结果会有差异;差异太大就不能用于决策。这一章回答三个问题:到底要抽多少次、如何判断收敛、哪些参数会直接影响收敛速度。

4.1 故障概率的量级决定了样本量的下限

配电网元件的不可用度通常在 10 的负 4 次方到 10 的负 3 次方量级。以 λ=0.1 次/年、r=5 小时/次为例,U = 0.1×5/8760 ≈ 5.7×10⁻⁵,平均每 17500 次试验才会遇到一次该元件故障。如果一个网络有 20 条支路,系统级故障事件的总概率大约是 1.1×10⁻³,也就是说平均要抽 1000 次试验才能碰到一次有故障的场景。

这就导出一个直观结论:如果抽样次数只有几千,统计结果会由少数的几次故障事件主导,方差极大。工程上我一般把抽样次数设到 2 万次以上,并对关键指标做方差监视。如果某个元件的故障率特别低,比如电缆线路只有 0.01 次/年,那它的贡献可能要 10 万次抽样才能稳定出现,此时可以考虑用重要抽样对其概率进行加权放大。

4.2 方差系数判据:比“抽满 N 次”更可靠的停算条件

固定抽样次数的问题是:同样的次数,在高故障率网络上已经收敛,在低故障率网络上可能还远未收敛。我在实际项目中用方差系数 β 作为停算判据,它的定义是估计值的标准差除以估计值本身,β 越小代表结果越稳定。工程上 β 小于 0.02 就认为结果可信,严谨的项目建议要求 β 小于 0.01。

% 方差系数计算示例: 每 1000 次抽样检查一次收敛情况 % 用停电时长的累计均值作为指标 N_block = 1000; block_means = []; for k = 1:N_sim [~, h_off, ~] = run_one_simulation(); h_values(k) = h_off; % 记录每次试验的停电时数 if mod(k, N_block) == 0 mu_h = mean(h_values(1:k)); sd_h = std(h_values(1:k)); beta = sd_h / (mu_h * sqrt(k)); block_means(end+1) = beta; if beta < 0.02 fprintf('第 %d 次抽样后收敛, beta=%.4f\n', k, beta); break; end end end

计算方差系数时有个容易被忽略的点:分母里的标准差要除以 sqrt(k),因为均值估计的标准差是按 1/sqrt(k) 衰减的。如果直接拿原始样本的标准差除以均值,得到的 beta 会偏大好几倍,导致不必要的过度抽样。另外,停电时长的原始分布高度右偏——大多数试验是 0,偶发一次大停电——所以均值收敛很慢,方差系数也会来回跳动。此时不要只看某一次检查点,建议连续 3 个检查点的 beta 都低于阈值才算真收敛。

4.3 三个必调参数:随机数种子、抽样上限、负荷随机化

第一,随机数种子。MATLAB 里rng(42)固定种子后,每次跑同一段代码结果完全一致,这在调试和验收时是刚需。如果用了rng('shuffle'),两次运行结果必然不同,收敛曲线也会错开,你根本分不清参数的改变是真实效果还是随机波动。建议种子作为函数入参,默认固定值,必要时在批量试验里按种子序号做并行推进。

第二,抽样上限。方差系数判据必须配一个硬性上限,避免低故障率场景下无限跑下去。常见做法是把上限设为 5 万次或 10 万次,达到上限仍未收敛时,程序要输出警告而不是静默结束。这能暴露“网络中有极低概率元件”的真实情况,提示你是否需要改用重要抽样。

第三,负荷随机化。工程评估中负荷不是恒定的,非时序方法里可以在每次试验给负荷乘一个随机系数,模拟峰谷波动。系数建议用正态分布或实际负荷分布采样,均值取年峰值负载率(比如 0.6 到 0.7),标准差取 0.1 左右。这样 ENS 的估计会更接近实际,而不至于系统性偏高或偏低。做完这步,SAIDI 的均值影响不大,但 ENS 的统计意义会明显改善。

5. 配电网蒙特卡洛评估的 5 个常见坑:从数据错位到潮流不收敛

代码跑通只是第一步,结果可信才是最终目标。这一章列出的 5 个坑,每一个我都实际踩过或者帮别人排查过,按“现象→原因→解决”的格式记录。

5.1 现象:SAIDI 算出来比历史报表高一倍,原因:修复时间直接当了停电时长

配电线路故障的处理过程是“故障定位→隔离→修复”,或者“故障定位→隔离→倒供→修复”。关键在隔离和倒供操作通常在修复完成之前就结束了,非故障段的负荷根本不用等修复。如果程序里把负荷点的停电时长直接全部按 r_repair 计,SAIDI 必然虚高。

解决办法是把影响分析函数细化为三段:故障段上游(靠近电源侧)通过重合闸快速恢复,时长远小于修复时间;故障段本身要等修复,时长是 r_repair;故障段下游若无联络开关则也要等修复,若有联络开关则只计倒供操作时间(通常 0.5 到 2 小时)。工程上把这三类时长分别建模,SAIDI 才和报表对得上。我在实际项目里见过把 5 小时的修复时间直接乘以全部下游用户数,最后 SAIDI 做到 2.3 小时,而实际年报只有 0.8 小时,问题就出在这里。

5.2 现象:两次运行结果差到 30% 以上,原因:随机数种子没有固定

如果你的代码里写了rng('shuffle')或者根本没写rng,那么每次运行的随机数流都不一样。改了一个支路的参数,你无法判断结果变化是来自参数调整还是随机波动。更隐蔽的是,有些人固定了rng,但抽样循环里调用了其他依赖随机数的函数(比如randperm或randn),改变了随机数流的位置,导致关键抽样分支的随机序列和上次运行错位。

解决方法是把rng(seed)放在主循环入口,并由输入参数控制;同时把每次试验的随机数状态用rng()输出保存到结果文件里。这样任何一次运行都可以精确重现,排查问题时也能回放当时的随机场景。

5.3 现象:改完数据表结果没变,原因:网络拓扑的派生结构没更新

很多人把支路表直接存在工作区变量里,主函数每次从这份表读数据。问题是供电路径集合、上游支路集合这些派生数据是在初始化阶段算出来的,你改了支路表的lambda或r_repair,这些派生结构一般会更新;但如果你改了网络连接关系(比如新增一段联络线),而初始化阶段只跑了一次,主循环里用的还是旧路径集合——结果当然不变。

解决方法是把网络解析和可靠性参数加载分开成两个函数,每次主程序启动强制重新解析网络生成供电路径,并且不要在同一个脚本里反复手动修改变量。我习惯把网络数据和派生结构的生成封装成一个build_network()函数,任何参数修改后都重新调用一次。这样也方便做多方案对比。

5.4 现象:潮流计算一直不收敛,原因:配电网络不适合用牛顿-拉夫逊法

有些评估框架会把潮流计算嵌进故障影响分析里,用于判断转供后是否越限。配电网的线路电阻较大,X/R 比低,而且部分节点在轻负荷时电压接近额定值附近,牛顿-拉夫逊法在配电网场景下经常出现雅可比矩阵奇异或者迭代发散。这不是算法写错了,而是方法选型问题。

解决方法是改用前推回代法。它对辐射状配电网天然友好,不需要求雅可比矩阵,迭代稳定,速度也快。实现时按拓扑从末端向电源回推电流,再从电源向前推算电压,两轮交替直到电压差小于阈值。前推回代法天然支持分布式电源接入,DG 的注入功率按 PQ 节点处理即可。需要留意的是,如果网络是弱环网,前推回代法需要做环路补偿,工程上要么先解环近似处理,要么直接用商用潮流引擎。

5.5 现象:转供成功率 100%,结果过于乐观,原因:没有容量约束

带联络开关的网络里,故障段下游可以通过联络线倒馈给相邻馈线。如果模型默认联络线容量无限大,下游负荷总能转供成功,那故障停电范围永远只有故障段自身,评估结果会比实际好得多。现实中联络馈线自身也在带负荷,转供后可能过载,电压也可能越限。

解决方法是把转供判断从“是否存在联络开关”升级为“联络馈线是否还有剩余容量,且转供后节点电压不越限”。工程量少时可直接比较负荷与剩余容量;网络规模大时,就要把潮流计算的结果作为判据。这个坑在接入分布式电源后更微妙:DG 出力能减轻转供压力,但可靠性评估和 DG 出力时序耦合,不能简单用一个恒定出力值代表全年。

6. 让评估快一个量级的技巧:预生成故障模式表与向量化改造

蒙特卡洛法最大的痛点就是计算慢,尤其是序贯版本,几万次试验乘上潮流计算,跑一次要几小时。这里给出两个经过验证的提速思路,再加一个结果可信度的收敛方法。

6.1 用故障模式影响分析表把重复计算移出主循环

配电网里的负荷点供电路径相对固定,同一类故障的影响模式其实非常有限。我一般会在蒙特卡洛循环之前,枚举所有单重故障(某些场景加双重故障),逐个计算每个故障事件影响的负荷点集合、停电时长和转供后的恢复情况,形成一张故障模式影响分析表。

% FMEA 表预生成: 对每条支路枚举故障 % 计算受影响负荷点, 不转供时的停电时长, 以及是否有转供能力 fmea_table = struct('faulty_id',{}, 'affected_nodes',{}, ... 'outage_hours',{}, 'can_transfer',{}); for i = 1:numel(br) [affected, hours, transfer] = analyze_single_fault(br, i); fmea_table(i).faulty_id = br(i).id; fmea_table(i).affected_nodes = affected; fmea_table(i).outage_hours = hours; fmea_table(i).can_transfer = transfer; end

蒙特卡洛主循环里就不再需要重复做网络拓扑搜索了,只需在每次试验抽样出故障元件集合后,在 FMEA 表里查对应行,合并影响结果。这个优化对故障率极低的网络特别有效——大量试验是没有故障的,查表操作几乎零成本。对于包含转供逻辑的网络,FMEA 表还可以同时计算转供后的潮流校验结果,把非线性计算从主循环里剥离出去。实际项目中,这一步通常能让单次评估从小时级降到分钟级。

6.2 向量化抽样:把逐元件 for 循环改成一次 rand 调用

MATLAB 的循环性能瓶颈主要来自解释器开销,能用向量化的地方就不要用 for。元件状态抽样是最典型的场景:一次rand()生成全部元件的随机数数组,一次逻辑比较得到全部故障状态。

% 向量化抽样替代逐元件循环 % U_all 是预先算好的所有元件不可用度向量 U_all = [br.lambda] .* [br.r_repair] / 8760; rand_all = rand(1, numel(br)); faulty_mask = rand_all < U_all; faulty_ids = find(faulty_mask);

这段代码把原来需要循环 N 次的抽样压缩成了三行。网络规模 100 条支路时,循环和向量化的耗时差异还不太明显;到 1000 条支路、10 万次试验时,差异就是数量级级别。矩阵运算在 MATLAB 里的底层是高度优化的 BLAS 实现,远优于逐元素解释执行。同样的思路也用在指标累加环节:把每次试验的停电统计量存入数组,试验结束后用cumsum和mean一次算完。

6.3 结果交叉验证:和解析法对拍一遍再上线

提速之后要做的是确认结果没被优化过程改坏。我习惯在完整网络评估前,先用一个极小的网络(比如 3 到 5 个节点)跑一次解析法,用最小割集法手算出 SAIDI 和 SAIFI 的期望值,再和蒙特卡洛结果对比。相对误差在 5% 以内,说明蒙特卡洛程序和指标口径正确;超出这个范围,多半是元件不可用度算法或停电时长分级逻辑有问题。这个方法几乎不花额外时间,却能在早期暴露系统性错误,比对着复杂网络猜原因高效得多。

另一个实用做法是把评估结果和电网年报横向对比。比如供电公司的公开可靠性年报通常会发布 SAIDI 和 SAIFI,如果你的评估区域和它重叠,可以拿同口径做一次校准。要注意的不是数值完全一致——评估用典型参数,年报是实际运行统计——而是量级和方向是否吻合:如果年报 SAIDI 是 1.1 小时,你算出来 15 小时,那一定存在模型层面的错误,不是参数微调能解决的。

我自己的习惯是每次新接一个配电网可靠性评估项目,都先做这两个对拍动作。它们帮我拦下过至少三次严重的数据口径错误:一次是故障率单位写错,一次是修复时间没区分隔离和修复,一次是联络线容量没建模。这些错误如果直接跑大网络,可能跑一周都发现不了。希望这篇笔记能帮你在做配电网可靠性评估时少走这些弯路,把时间花在更实质的网架分析上。

本文还有配套的精品资源,点击获取

RELATED READING

延伸阅读

更多一线实战笔记与深度复盘,助您持续精进