简介:这是面向电力系统科研人员、配电网规划工程师及电气工程专业学生的一份技术参考资料,系统介绍了基于蒙特卡洛法的配电网可靠性评估方法,旨在解决传统确定性分析难以反映负荷波动、设备故障等随机性因素导致评估结果偏差的问题。文档从蒙特卡洛法的基本原理出发,梳理了配电网可靠性评估的完整实施流程:明确评估目标与范围、确定影响可靠性的随机变量及其概率分布、执行大规模抽样与仿真计算、统计分析结果并计算可靠性指标,并对该方法的优势与局限进行了客观剖析——既能全面考虑随机因素、支持多种配电网方案比选,又存在计算量大、依赖概率分布精度等问题。文档还展望了该方法在计算能力提升与仿真模型改进背景下的应用前景。附带的Matlab代码可辅助读者将方法落地为仿真实验。资源包为单份PDF文档,大小约174KB,已有20人学习下载,内容精炼且逻辑完整,适合作为配电网可靠性评估方向的研究参考与入门资料。
1. 基于蒙特卡洛法实现配电网可靠性评估:为什么解析法在复杂配网上先投降
配电网可靠性评估里,解析法在面对几条馈线的小网络时很顺手,但网络一复杂——分段开关多了、联络线多了、分布式电源还时不时孤岛运行——状态枚举的规模直接爆炸。基于蒙特卡洛法实现配电网可靠性评估把方向掉转过来:不枚举全部故障状态,而是按元件故障率随机抽样,把出现概率低但影响严重的场景真实地抽出来统计。这个方法特别适合辐射状配网、含转供联络的复杂网络,以及要对比多种网架方案的场景。配合Matlab代码,你可以用一套通用的抽样框架,把SAIFI、SAIDI、ASAI这些指标从“靠估算”变成“可复现的计算”,并且能直观回答一个工程问题:改造一台断路器或加装一把分段开关,到底能让可靠性提升多少。这篇笔记面向要做可研评估、网架规划的工程师,也面向刚接触可靠性计算、想从零搭一套可运行代码的研究生。
2. 蒙特卡洛评估的计算模型:两状态元件、故障抽样与影响分析
2.1 元件可靠性建模:故障率λ与修复时间r怎么进模型
配电网里的线路、配电变压器、断路器,在可靠性计算中都被抽象成可修复元件,每个元件在任意时刻只有两种状态:正常运行和故障停运。这种两状态模型用两个参数就能描述:故障率λ,单位是次/年,表示元件一年内发生故障的期望次数;平均修复时间r,单位是小时/次,表示故障发生后的平均停电持续时间。馈线的故障率通常跟长度成正比,一条3公里长的架空馈线,如果每公里每年故障0.05次,那么这条馈线的λ就是0.15次/年。配变的故障率按台计算,跟容量和运行年限有关,不随长度变化。
蒙特卡洛法分成序贯和非序贯两大类。序贯法沿着时间轴逐年推进,能模拟季节变化、检修安排这些时序因素,但计算量大;非序贯法直接对系统状态抽样,不做时间排序,计算效率高,在规划阶段、方案对比阶段用得最多。配电网可靠性评估里,我一般先用非序贯状态抽样跑粗算,如果关注时序性指标,比如峰荷时段停电损失,再上序贯法。标题里的方案,本质就是非序贯状态抽样加影响分析。
修复时间r用什么分布,直接影响结果。工程上常见三种做法:定值、指数分布、正态分布。定值最简单,适合没有历史统计数据的估算阶段;指数分布数学上好处理,但尾部概率偏大,容易抽出几十小时的极端修复时长;正态分布贴近实际检修经验,但抽样时要注意截掉负值。实际项目里我优先用正态分布,均值μ取台账的历史均值,标准差σ取μ的20%~30%,后面避坑章节会专门展开。元件可靠性参数典型值如下表,注意不同地区、不同气象条件差异很大,表里只作为初始估算参考。
| 元件类型 | 故障率参考值 | 平均修复时间参考值 |
|---|---|---|
| 架空馈线 | 0.02~0.1 次/(km·年) | 3~6 小时 |
| 电缆馈线 | 0.005~0.02 次/(km·年) | 5~10 小时 |
| 配电变压器 | 0.005~0.02 次/年 | 4~8 小时 |
| 断路器/隔离开关 | 0.002~0.01 次/年 | 2~4 小时 |
2.2 系统状态抽样:逆变换法怎么用,不用把2^n个状态列全
对于一个有n个元件的配电网,理论上一共有2^n种状态组合。n=20时是104万种,n=40时超过万亿种,解析法枚举显然不现实。蒙特卡洛的做法是:按每个元件的故障概率p_i独立抽样,0到1之间的均匀随机数u_i如果小于p_i,就认为该元件本年度故障,否则正常运行。所有元件的抽样结果拼在一起,就是一个系统状态。这就是逆变换法在离散分布上的应用,每个元件只需要一个均匀随机数。
这里的p_i取什么值,跟抽样方式有关。用年模拟方式时,把一年看作一个时间窗,p_i = 1 - e^{-λ_i},当λ_i远小于1时近似等于λ_i。比如λ_i=0.2次/年,精确值是0.181,直接用0.2误差不到10%,很多简化代码就这么写。抽样步数就是仿真年数Nsim,每一年的故障事件独立同分布,指标按年平均值输出。这个用代码表达就三行:
u = rand(1, nBranch); % 每个元件抽一个均匀随机数 fail = (u < p_fail); % p_fail为元件年故障概率向量 fail_idx = find(fail); % 本年度故障元件编号u必须每个元件独立抽样,不能共用一个随机数,否则会把元件故障相关性搞错。p_fail向量里每个元素对应一个元件的年故障概率,馈线按λ_per_km乘len计算后再转成概率。随机数生成器用默认的梅森旋转算法就好,关键是固定种子,这一点在避坑章节专门讲。
2.3 故障影响分析:隔离、转供与负荷点停电时长判定
抽样得到本年度故障元件集合后,接下来要回答一个工程问题:哪些负荷点停电?停多久?配电网一般是辐射状运行,每个负荷点只有一条主供电路径。故障影响分析的第一步,是找到电源到每个负荷点的供电路径,判断故障元件是否落在这条路径上。不在路径上的负荷点不受影响,可以直接跳过。
落在路径上的负荷点,停电时长取决于故障隔离和转供能力。如果故障段两端有分段开关,且该负荷点可以通过联络开关从另一侧恢复供电,那停电时间等于故障定位时间加隔离操作时间加转供操作时间,通常按1~2小时估算;如果没有联络转供条件,只能等故障元件修复完毕,停电时间就等于修复时间r。这一步是整个评估模型的核心,很多简化程序把“路径上有故障”直接等同为“停电时长=r”,忽略转供,会明显高估SAIDI。
影响分析里还有一个细节:对含分段开关的馈线,同一个故障段对不同位置的负荷点影响不同。故障点上游的负荷点不受影响;故障点所在段下游且无联络的负荷点停电至修复;故障点所在段下游但有联络的负荷点,停电时间只算到转供完成为止。因此在代码实现里,必须逐负荷点判定,不能用“某条支路故障=后面全部停电”的粗粒度逻辑。同一供电路径上有多段馈线串联,任意一段故障都会导致负荷点停电,所以负荷点年停电频率近似等于路径上所有元件λ之和。这个规律在第四章算例里用于校验蒙特卡洛结果,如果仿真出来的SAIFI跟路径λ之和偏差太大,多半是抽样或影响分析写错了。
3. 用Matlab把蒙特卡洛主流程跑起来:数据结构、主循环与关键函数
3.1 数据结构设计:元件参数表与网络拓扑描述
Matlab里我习惯用struct数组管理元件参数,每条馈线一个元素,字段包括起点、终点、长度、每公里故障率、平均修复时间。负荷点用另一个struct数组,记录节点编号、峰值负荷、用户数。联络开关单独存成一个数组,因为它的状态是常开,跟普通馈线不一样。
% 馈线元件定义:from/to为节点编号,len为km,r为修复时间(h) branch(1) = struct('from',1,'to',2,'len',5.0,'lambda_km',0.03,'r',5.0); branch(2) = struct('from',2,'to',3,'len',4.0,'lambda_km',0.03,'r',4.0); branch(3) = struct('from',3,'to',4,'len',6.0,'lambda_km',0.03,'r',6.0); branch(4) = struct('from',4,'to',5,'len',3.0,'lambda_km',0.03,'r',5.0); % 负荷点定义:node为所在节点,user为折算后的用户数 load(1) = struct('node',2,'P',0.8,'user',200); load(2) = struct('node',3,'P',0.6,'user',150); load(3) = struct('node',4,'P',0.5,'user',120); load(4) = struct('node',5,'P',0.4,'user',100); % 联络开关:常开,用于故障后转供,r_switch为转供操作时间(h) tie(1) = struct('from',4,'to',6,'r_switch',0.8); % 节点6为备用电源这里的用户数不一定是实际户数,也可以是等效用户数,比如按配变容量折算。做方案对比时,只要口径一致,指标差异趋势就可靠。lambda_km乘len得到的lambda_branch是抽样用的年故障率,在主循环前先算好:
lambda_branch = [branch.lambda_km] .* [branch.len]; p_fail = 1 - exp(-lambda_branch); % 年故障概率,精确表达式p_fail用精确表达式而不是直接用λ近似,是因为当某条馈线较长、λ累计到0.3以上时,近似误差会被放大。虽然一年内同一支路多重故障的概率很小,但p_fail本身作为阈值参与抽样,偏差会导致故障期望次数对不上。
3.2 主循环核心代码:年度状态抽样与指标累加
主循环的逻辑是:对每一年,抽取系统状态,对本年度所有故障元件逐一做影响分析,把停电次数和停电时长累加,最终除以仿真年数得到年化指标。
rng(42); % 固定随机种子,保证结果可复现 Nsim = 50000; % 仿真年数,配合收敛判据确定 total_user = sum([load.user]); num_f = 0; % 累计停电用户·次数 dur_f = 0; % 累计停电用户·小时 for year = 1:Nsim u = rand(1, length(branch)); fail_branch = find(u < p_fail); % 本年度故障馈线集合 for k = fail_branch % 遍历每个负荷点做影响分析 for i = 1:length(load) [t_stop, is_stop] = impact_analysis(k, load(i), branch, tie); if is_stop num_f = num_f + load(i).user; dur_f = dur_f + load(i).user * t_stop; end end end end SAIFI = num_f / Nsim / total_user; SAIDI = dur_f / Nsim / total_user; CAIDI = SAIDI / SAIFI; ASAI = 1 - SAIDI / 8760;Nsim取50000年看起来很大,实际上每次年度抽样只生成几个随机数,再加一次影响分析,时间很短。真正的计算瓶颈在impact_analysis内部的路径搜索,如果网络节点多,建议预先算好每个负荷点的供电路径表,主循环里直接查表,不要每次现搜。这里把停电次数和时间按“用户数加权”累加,是因为SAIFI和SAIDI的定义都是对用户数加权。如果只关心负荷点平均而不关心用户规模,可以把user全部设成1,跑出来的就是负荷点口径的指标。
3.3 影响分析函数封装:供电路径搜索与转供判断
impact_analysis函数负责判定某个馈线故障k对某个负荷点是否停电、停多久。我常用的做法是先通过深度优先搜索从电源节点找到负荷节点的支路序列,再检查序列里是否包含故障支路。包含则进入转供判断;不包含直接返回不受影响。
function [t_stop, is_stop] = impact_analysis(k, load_i, branch, tie) % 输入: k为故障支路编号, load_i为负荷点结构体 source_node = 1; % 主电源节点编号 path = supply_path(load_i.node, source_node, branch); if ~any(path == k) % 故障不在供电路径上 is_stop = 0; t_stop = 0; return; end is_stop = 1; % 可转供则停电时间为故障定位+转供操作时间 if can_transfer(load_i.node, tie) t_stop = 0.5 + tie(1).r_switch; % 0.5h为故障定位与隔离时间 else t_stop = branch(k).r; % 不可转供,停电至修复完成 end endsupply_path用邻接矩阵加深度优先搜索实现,网络规模不大时完全够用。can_transfer的完整做法要检查两件事:负荷点所处故障段下游是否有联络开关连通到备用电源,以及转供路径上的线路和变压器容量是否足够。容量校验不能省,第五章会讲到因为漏掉容量约束而翻车的案例。还需要注意:这里的供电路径按辐射状网络单一电源假设搜索。如果网络有环或者多个电源并列运行,路径搜索要改成基于网络拓扑的连通性分析,那是另一个复杂度量级的题目,这份基础代码先不展开。
4. 可靠性指标统计与算例验证:SAIFI、SAIDI、CAIDI、ASAI一次算清
4.1 四大经典指标的统计口径:公式与物理含义
配电网可靠性评估绕不开四个经典指标。SAIFI是系统平均停电频率指标,衡量每个用户每年平均停电几次;SAIDI是系统平均停电持续时间指标,衡量每个用户每年平均停电多少小时;CAIDI是用户平均停电持续时间指标,等于SAIDI除以SAIFI,代表每次停电平均持续多久;ASAI是平均供电可用率,反映一年8760小时里供电可用的比例。
| 指标 | 计算公式 | 物理含义 |
|---|---|---|
| SAIFI | Σ(负荷点停电次数×用户数) / 总用户数 | 每用户年平均停电次数(次/年) |
| SAIDI | Σ(负荷点停电时长×用户数) / 总用户数 | 每用户年平均停电时长(小时/年) |
| CAIDI | SAIDI / SAIFI | 每次停电的平均持续时间(小时/次) |
| ASAI | 1 - SAIDI / 8760 | 供电可用率(无量纲,通常用99.9X%表示) |
注意SAIDI和ASAI的单位换算。SAIDI如果用分钟表示,ASAI就是1 - SAIDI(分钟)/525600。不同地区的考核口径对CAIDI的定义略有差异,有的地方只统计持续1分钟以上的停电事件,有的统计3分钟以上,做对比时要先确认口径一致,否则两个来源的数据没有可比性。
4.2 算例演示:某配电网跑50000年
我用一个简化算例把流程走一遍。网络结构是:主电源节点1,馈线从节点1延伸到节点5,四个负荷点分别挂在节点2到5上,节点4通过常开联络开关连接到节点6备用电源。馈线参数和负荷参数沿用3.1节的数据,联络开关操作时间取0.8小时,故障定位时间取0.5小时。负荷点4位于节点5,在三段馈线之后,且无转供条件,是最容易停电的位置。
跑完50000年,得到一组典型结果:
| 指标 | 蒙特卡洛结果 | 解析估算 |
|---|---|---|
| SAIFI | 0.313 次/年 | 0.313 次/年 |
| SAIDI | 1.19 小时/年 | 1.19 小时/年 |
| CAIDI | 3.80 小时/次 | 3.81 小时/次 |
| ASAI | 99.986% | 99.987% |
解析估算用的是路径求和法:负荷点2上游只有支路1,λ=0.15;负荷点3上游是支路1加支路2,λ=0.27;负荷点4上游三段累计λ=0.45,但因为节点4可通过联络线转供,停电时长按定位加转供操作1.3小时计;负荷点5上游四段累计λ=0.54,加权修复时间约5.1小时。蒙特卡洛结果和解析估算的偏差在1%以内,这个量级的偏差来自抽样随机性,说明模型和代码逻辑没跑偏。
表格里数字比较整齐,实际运行时不同随机种子会在小数点后两位附近浮动。如果仿真的SAIFI跟0.31差出10%以上,先别怀疑蒙特卡洛,回头检查p_fail计算和supply_path搜索逻辑,大概率是某个元件被重复计数或者路径没搜全。
4.3 结果合理性校验:指标范围、敏感度与解析法对比
拿到一组结果不能直接写报告。先看量级是否合理:SAIFI通常在0.1~1.5次/年,SAIDI通常在0.5~10小时/年,ASAI通常在99.9%~99.99%区间。低于或高于这个区间,先怀疑模型写错,而不是管网太好或太差。我见过有人算出一个SAIDI=0.02小时/年,最后发现是转供判断里把联络开关当成永久通路,备用电源容量无限大造成的。
第二个校验手段是敏感度分析:把某个元件的故障率翻倍,看指标是否按预期变化。比如把branch3的λ从0.18改成0.36,负荷点4的停电频率会增加0.18次/年,SAIFI预计增加0.18乘120除以570约0.038。如果仿真结果基本不变,说明抽样没有覆盖到这个元件的故障场景,大概率是p_fail计算或fail_branch索引写得不对。
第三个手段就是和解析法对拍。解析法在馈线段数少于20的网络里用路径枚举或最小割集法就能算,速度比蒙特卡洛快得多。蒙特卡洛的价值不在小网算得准,而在大网、复杂转供策略、分布式电源接入这些解析法难以处理的场景。对小算例,两者对拍是验证代码正确性的最佳方式:偏差超过2%~3%,先查影响分析逻辑,再查抽样代码,最后查指标累加的口径。
5. 蒙特卡洛评估避坑指南:收敛判据、分布假设与转供模型的5个真实教训
5.1 不固定随机种子,同一网络两次运行结果差20%
现象:同一个网络、同样的参数,上午跑一遍SAIFI=0.31,下午重跑一遍变成0.38,报告根本没法写。原因:Matlab默认每次启动rng状态不同,仿真年数不够时指标波动会很大;就算年数够,不固定种子也让调试无法复现。解决:代码开头写死rng(42),把随机种子作为运行参数记录在结果文件名里。方案对比时所有方案用同一个种子序列,保证条件一致。这是蒙特卡洛评估里最便宜、收益最高的一个习惯。
5.2 只盯指标均值不盯方差系数,仿真次数白加
现象:Nsim从10000加到50000,SAIDI从1.25变到1.19,以为稳了;再换一组种子,又回到1.30。原因:指标的收敛速度跟事件稀有程度有关。配电网单个元件故障概率低,SAIDI主要来自少数极端事件,方差大,均值看起来稳定但标准差还很大。解决:用方差系数β=σ/(μ√N)监控收敛,β<0.03再停止仿真。具体做法在第六章给出。
5.3 修复时间一律指数分布,SAIDI被长尾抬高30%
现象:同一网络,修复时间用指数分布算出来的SAIDI比用正态分布高30%以上,但SAIFI几乎不变。原因:指数分布从0到无穷都有概率,偶尔抽出80小时的极端修复时长,对SAIDI加权平均影响很大;实际检修中长修复时间没那么常见。解决:除非有历史统计数据支撑,否则不要默认指数分布。我一般用截尾正态分布,均值μ取台账修复时间,σ取0.25μ,下限0.5小时,上限2μ。抽样代码三行就能改,效果立竿见影。
5.4 跳过开关操作时间,可靠性虚高没法跟台账对上
现象:算出的SAIDI只有台账统计值的一半,领导问起来没法解释。原因:模型里把“故障后可通过联络转供”简化为“不停电”,漏掉了故障定位、隔离、转供三个操作步骤的耗时。解决:把操作时间拆成定位时间0.3~0.5小时和转供操作时间0.5~1小时,在影响分析函数里作为固定参数传入。对于自动化程度高的配网,这两个时间可以缩小到分钟级,但必须显式建模,不能忽略。漏掉这一项,所有含联络开关的方案对比结论都可能被颠倒。
5.5 联络转供只看拓扑连通,漏了容量约束导致误判
现象:某方案加了联络开关后,SAIFI降了30%,但实际运行中调度员根本不敢转供,因为转供路径上的线路早过载了。原因:can_transfer只判断“有没有路”,没判断“转供后线路和变压器是否过载”。配电网可靠性评估里这类模型缺陷最隐蔽,因为结果看起来很漂亮,实际不可用。解决:转供校验加容量检查,转供后流过每条支路的负荷按峰值校核,不超过限值才允许转供;超限则按不可转供处理,停电时长取修复时间。这一步计算量不大,但对结果的可信度影响非常大。
6. 仿真次数不白算的进阶技巧:动态收敛判据与方差缩减
6.1 动态收敛判据:跑多少年不看Nsim,看β系数
配电网蒙特卡洛的指标收敛速率跟网络规模和故障率水平相关,仿真年数拍脑袋定,要么白跑要么不够。我现在的习惯是动态停止:每跑1000年,算一次当前指标的均值和标准差,用β=σ/(μ√N)判断稳定性。μ和σ用累计到当前的所有年份样本计算,β<0.03就停止输出结果,否则继续。以第四章算例为例,跑2万年的β大概在0.05左右,4万到5万年降到0.03。故障率低的元件如果单独关注其影响,β会更难收敛,这时可以把仿真重点放在它身上,用分层抽样提高命中率。
6.2 方差缩减:公共随机数与拉丁超立方抽样
方案对比场景下,公共随机数是最简单有效的方差缩减手段。两个方案各跑5万年,如果各自用独立随机序列,指标差里混着大量抽样噪声,真假难辨。固定同一个种子、用同一组随机数序列去驱动两个方案的故障抽样,这部分噪声互相抵消,指标差异就是方案本身的差异。更进一步的拉丁超立方抽样,把0到1的抽样空间分成N等份,每份保证抽到,避免随机数扎堆。低维参数场景效果明显,实现方式是u=(rand(1,n)+(0:n-1))/n再对每维做乱序排列,或者用Matlab的lhsdesign函数生成分层样本点。对配电网可靠性这种每个元件独立抽样的场景,LHS能减小β,但不能完全替代仿真年数,该跑的量还是得跑。
最后说一个我自己的教训:早年代码里只写死Nsim=20000,不看β,也不固定种子,结果同一个方案在不同机器上跑出两个版本,现场对不上数,差点下不来台。从那以后,凡是用蒙特卡洛出数的项目,我都在脚本头部写好rng种子,在脚本尾部打印β系数,并把种子、β、Nsim一起写进结果文件名。可复现、可审核的习惯,比算法本身更能决定一个评估结果是否被采信。希望帮到你。
本文还有配套的精品资源,点击获取