配电网可靠性评估:基于序贯蒙特卡洛模拟法的Matlab实现
2026/9/10 1:23:25 网站建设 项目流程

序贯蒙特卡洛模拟法做配电网可靠性评估,这个项目我前后折腾了快三个月。如果你也是电力系统方向的学生,或者刚进入配电网规划岗位,八成会遇到类似的需求:领导让评估现有馈线的供电可靠性,或者做改造方案的横向对比,需要一套能算SAIFI、SAIDI、ENS这些指标的工具。Matlab是首选实现语言,因为它的矩阵操作、随机数生成、绘图都比其他脚本语言方便,而且学术界、工业界都有大量现成代码可参考。这篇博文就围绕“基于可靠性评估序贯蒙特卡洛模拟法的配电网可靠性评估研究(Matlab代码实现)”来写,把原理、代码架构、实操流程和踩坑经验一次说清楚。

1. 项目定位:配电网可靠性评估到底在“评估”什么

1.1 先搞清楚指标:SAIFI、SAIDI、CAIDI、ASAI、ENS

配电网可靠性评估的核心是回答一个问题:在给定网络结构和运行方式下,用户一年到底要停几次电、停多久。工程上很少用一个指标概括全部信息,而是用一组国际通用的可靠性指标来刻画。

  • SAIFI(系统平均停电频率):单位是“次/用户·年”,计算公式是系统用户停电总次数除以用户总数。这个指标反映停电有多频繁,和故障率直接挂钩。
  • SAIDI(系统平均停电持续时间):单位是“小时/用户·年”,计算公式是系统用户停电总时长除以用户总数。它反映停电有多久,受修复时间、隔离时间、转供时间共同影响。
  • CAIDI(用户平均停电持续时间):等于SAIDI除以SAIFI,单位是“小时/次”。它衡量一次停电平均持续多久,是供电企业服务水平的直观体现。
  • ASAI(供电可用率):等于1减去SAIDI除以8760,是一个接近1的小数,比如0.9991。它表示用户一年里可用电时间所占比例。
  • ENS(期望缺供电量):单位是“MWh/年”,计算各类停电场景下缺少的电量总和。这个指标直接和经济损失挂钩,规划人员最看重它。

这五个指标并不是孤立的,它们之间满足“频率×每次时长=总时长”的逻辑关系。写报告的时候通常把这五个指标全部列出来,因为不同场景下决策视角不同——SAIFI低但SAIDI高,说明停电虽然不频繁但一次停很久,问题出在检修或转供能力不足;反过来则是故障太多但处理很快。

1.2 为什么配电网的评估比输电网更麻烦

输电网通常是环形网络,局部故障可以通过多重联络自动转带,评估时用确定性N-1准则往往就够了。配电网不一样,绝大多数处于辐射状或弱环状运行,馈线之间依靠联络开关连接,故障隔离和恢复供电压根结底是一系列复杂的开关动作逻辑。

一个10kV馈线上通常有几台分段开关、联络开关,还可能有柱上断路器、熔断器。某段线路发生故障后,变电站出线断路器先跳闸,整条馈线全部失去电源;随后调度或自动化系统定位故障,通过拉开故障点两侧的隔离开关把故障段隔离,再合上联络开关把非故障段转带到其他馈线。这一系列动作的时间和成功率直接决定用户停电时长。解析法处理这种逻辑时公式会变得非常冗长,而且一旦线路、开关数量增加,状态空间以指数级膨胀,手算公式根本不现实。

蒙特卡洛模拟法的优势就在这里:不追求穷举所有状态,而是按照元件的故障规律随机抽取大量“故障—修复”事件,用统计平均还原出系统的长期可靠性水平。这也是为什么在配电网可靠性评估领域,蒙特卡洛法已经成为事实上的主流工具。

1.3 序贯法和非序贯法怎么选

蒙特卡洛模拟又分非序贯和序贯两类。非序贯法抽样的是系统状态——比如某时刻若干元件同时处于故障状态,它对时间维度处理很粗糙,无法表达“先故障后恢复再转供”这类有时间先后逻辑的过程。序贯法全称是序贯蒙特卡洛模拟法,它把时间轴真实推进,先抽样每个元件的无故障工作时间(TTF)和故障修复时间(TTR),再在时间轴上合并所有元件事件,逐段判断系统状态。

正因为它是按时间走的,所以能自然地处理时变负荷、分布式光伏出力的时间相关性问题,也能模拟一天内不同时段转供策略差别带来的影响。代价是计算量更大,程序也更复杂。如果只是搞个简单辐射网的教学演示,非序贯法也能凑合;一旦涉及实际馈线、DG、储能,序贯法几乎是必选项。这个项目选序贯蒙特卡洛模拟法,方向是对的。

2. 序贯蒙特卡洛模拟法的核心原理,用一个简单例子讲透

2.1 元件状态持续时间抽样:指数分布与随机数

序贯法的底层逻辑是“随机抽样元件状态持续时间”。架空线路、变压器、断路器这些元件,在可靠性建模中最常见的是两状态模型:正常运行和故障停运。正常运行时间TTF和故障修复时间TTR都是随机变量,工程上一般假设它们服从指数分布。

指数分布有一个非常好的性质:累计分布函数F(t)=1-e^(-λt),其中λ是故障率。给定[0,1]区间均匀随机数U,令U=F(t),反解出t=-ln(1-U)/λ。由于U和1-U同分布,可以简化为t=-ln(U)/λ。这就是蒙特卡洛抽样状态持续时间的核心公式。

在Matlab里写出来就两句:

function ttf = sampleTTF(lambdaPerYear) % lambdaPerYear: 元件年故障率,单位 次/年 % 返回正常运行时间,单位 小时 ttf = -log(rand()) / lambdaPerYear * 8760; end function ttr = sampleTTR(mttr) % mttr: 平均修复时间,单位 小时 % 返回本次修复时间,单位 小时 ttr = -log(rand()) * mttr; end

注意第一段代码里除以“次/年”得到的是“年”,乘8760换算成小时。如果直接用小时故障率λh=λ/8760,那么TTF=-ln(U)/λh,效果一样,只是单位要时刻盯紧,这是我最早最容易出错的地方。

2.2 把多个元件的状态序列在时间轴上合并

单台设备的抽样式子很简单,但系统里有几十条线路、多台变压器和开关,怎么把它们合成一个全局的过程?序贯法是这么做的:

  1. 对每个元件,先抽样一次TTF,记录下来,把它当作该元件的首个事件。
  2. 找到所有元件事件中时间最早的那一个,把系统时钟推进到那个时刻,处理该事件(比如某条线路故障)。
  3. 被处理的元件进入修复阶段,抽样一个TTR;修复结束后再继续下一次TTF。其他未发生事件的元件状态保持不变。
  4. 重复“找最早事件→推进时间→执行故障后果分析→安排下一次事件”这一循环,直到模拟的总时间达到预设年限(比如10000年)。

这个方法实现时通常会维护一个事件列表,每一行记录“元件编号、事件类型(故障/修复)、发生时刻”。每次从列表里取最小时刻的事件,处理完后用新抽样的事件替换旧事件,然后重新排序。Matlab的排序指令sortrows很好用,但如果模拟年限长、事件多,循环时间会相当可观,后面我会说怎么优化。

2.3 故障后果分析是精确性的关键

系统级模拟只负责生成“何时发生故障”,真正决定可靠性指标的是故障后果分析。配电网故障后,不同负荷点停电时间来源不一样:故障点前端的负荷可能只停几分钟就被重合闸或隔离开关恢复;故障点后端的负荷可能要等联络开关转供,也可能一直停到故障修复结束。

具体分析时要按馈线的拓扑结构做区域划分。以典型辐射状馈线为例,某段线路故障后:

  • 变电站出线断路器先跳闸,全线失电。如果重合闸、备自投逻辑存在,瞬时性故障可能在几十秒内恢复,但在最高一级的长时间可靠性统计里,这类短时停电通常也要计入SAIFI。
  • 通过隔离开关把故障段隔离后,故障点前端的非故障段可以恢复送电,停电时间等于开关操作时间。
  • 故障点后端的非故障段,若存在联络开关且有足够备用容量,则可以通过转供恢复,停电时间等于转供操作时间;如果没有转供路径,只能等故障修复完毕,停电时间等于TTR。
  • 故障区段内的负荷点,只能等待修复完成,整个过程停电。

在实际代码里,这种逻辑需要结合开关位置、电源点位置和联络线位置实现搜索算法。最常见的方法是邻接矩阵+深度优先搜索,先把网络按照开关设备分成若干个区段,故障时只需要判断各区段与电源的连通性,不必逐个负荷点去做拓扑搜索,计算效率高很多。

3. Matlab代码实现:从数据结构到核心函数的完整拆解

3.1 输入数据怎么组织才不容易乱

配电网仿真最考验数据结构设计。我的建议是全部结构化存储,不要用零散的矩阵变量,否则传到后面自己都看不懂。

我常用四个结构体:

% 节点表 node.id % 节点编号 node.type % 类型:1-变电站电源, 2-负荷节点, 3-普通连接节点 node.loadMW % 平均负荷,单位 MW node.userNum % 该节点折算的用户数 % 支路表(线路、变压器统一处理) branch.id % 支路编号 branch.fromNode % 起始节点 branch.toNode % 终止节点 branch.lenKm % 长度,单位 km branch.lambdaPerKm % 每公里年故障率 branch.mttr % 平均修复时间,单位 小时 branch.isSwitch % 是否含开关 branch.switchType % 开关类型:0-无, 1-断路器, 2-隔离开关, 3-联络开关 branch.openTime % 开关操作时间,单位 小时 % 联络开关表(用于转供路径搜索) tieSwitch.id tieSwitch.branchId tieSwitch.nodeA tieSwitch.nodeB % 模拟参数 simParam.years % 模拟年限 simParam.seed % 随机数种子 simParam.convergeCriterion % 收敛判据阈值

把线路、变压器统一到同一张表里,对代码实现是最方便的,因为它们的可靠性模型本质上相同,只是故障率、修复时间参数不同。在实际项目中,我需要把GIS系统导出的cad/dwg表格整理成这种结构,过程比较繁琐,但一次整理清楚后面跑仿真就很顺。

3.2 事件表驱动的主循环怎么做

用一个N行三列的矩阵events来维护当前的事件表,每一行格式是[元件编号, 事件发生时刻(小时), 事件类型]。事件类型用1表示故障开始,0表示修复完成。初始化时给每个元件生成一次TTF,作为它的第一个故障事件;后续当某个元件被处理完,就补上它的下一次TTF或TTR事件。

主循环大致是这个样子:

function [metrics] = runSequentialSimulation(node, branch, tieSwitch, simParam) rng(simParam.seed); nBranch = length(branch); % 初始化事件表 events = zeros(nBranch, 3); for k = 1:nBranch ttfHours = sampleTTF(branch(k).lambdaPerKm * branch(k).lenKm); events(k,:) = [k, ttfHours, 1]; % 第一个事件都是故障 end % 累积统计变量 totalStopNum = 0; % 总停电用户次数 totalStopHour = 0; % 总停电用户小时数 totalLackEnergy = 0; % 总缺供电量 MWh totalYears = 0; % 已经模拟的年限 while totalYears < simParam.years % 取时间最早的事件 [minTime, idx] = min(events(:,2)); k = events(idx,1); evType = events(idx,3); % 更新总模拟时间,跨年时统计年度累计值 ... if evType == 1 % 故障开始 faultBranch = k; [impactedLoads] = analyzeFault(node, branch, tieSwitch, faultBranch); % 根据影响的负荷点,累加停电次数、停电时间、缺供电量 ... % 下一次事件:修复完成 ttrHours = sampleTTR(branch(k).mttr); events(idx,:) = [k, minTime + ttrHours, 0]; else % 修复完成 % 下一次事件:下一次故障 ttfHours = sampleTTF(branch(k).lambdaPerKm * branch(k).lenKm); events(idx,:) = [k, minTime + ttfHours, 1]; end end end

这个框架简单但实用,胜在逻辑清晰。需要注意一个容易忽略的细节:模拟到模拟年限的边界时,如果跨越了整年分界点,需要把前面已经累加的指标先结算一次,然后再把总模拟时间清零重新累计。很多新手会漏掉这一步,导致每年统计的天平被截断,结果偏差很大。

3.3 故障后果分析函数怎么搜拓扑

故障后果分析是整个程序里最考验逻辑的部分。我采用“区段划分+广度优先搜索”的策略,具体步骤如下:

  1. 把所有开关(包括断路器、隔离开关、联络开关)当作天然的分断面,网络被它们切分成若干个区段。
  2. 找到故障区段,假设是segment_fault。
  3. 从变电站电源节点出发做广度优先搜索,遇到断开的开关停止扩展。这能得到故障隔离前有电的区段。
  4. 故障隔离后,把故障段两侧最近的开关拉开。此时电源侧的区段都恢复供电;非故障侧区段如果通过联络开关能连接到其他电源,则标记为转供区段。
  5. 根据区段类型逐一统计每个负荷点的停电频率和停电时长。

这条逻辑里最容易出错的是联络开关转供时还要校验备用容量是否充足。如果另一回馈线自身负载率已经很高,转供可能导致过载,工程上不能简单认为“一连就通”。在简化模型里可以设置一个可转供容量上限,超过则不允许转供,此时下游负荷只能等待故障修复。

用一个Matlab函数封装:

function [stopFreq, stopDur, lackEnergy] = analyzeFault(...) % 输入故障支路编号、网络拓扑、开关状态 % 输出每个负荷点的停电频率、停电时间、缺供电量 end

这个函数的返回值需要设计精确,因为后续所有可靠性指标都建立在它的基础上。

3.4 指标统计和收敛判据

模拟结束后,把所有年份累计的停电次数、停电时间、缺供电量分别除以总年份数或者用户数,就得到最终的可靠性指标。但这里有个关键问题:到底模拟多少年才够?

蒙特卡洛法天生带有随机误差,误差大小用方差系数β来衡量,β=σ/(μ√N),其中σ是样本标准偏差,μ是样本均值,N是模拟年限。工程上通常要求β≤0.05,对高精度场景要求0.01。我的建议是在主循环里每模拟100年就检查一次β,一旦满足收敛条件就提前终止,而不是傻乎乎地跑固定的年份。这样能在保证精度的前提下省大量计算时间。

beta = sigma / (mu * sqrt(N)); if beta < 0.05 break; end

注意SAIFI、SAIDI、ENS各指标对应的β不同,收敛速度也不一样。实际项目中一般只对最关心的那个指标做收敛检查,通常是SAIDI或ENS,因为它们波动更大。

4. 实操复盘:一个三馈线测试系统的完整评估

4.1 测试系统怎么搭

为了验证代码,我用一个简化的10kV测试系统,三回馈线,每回馈线带6个负荷点。三回馈线从同一座110kV/10kV变电站的不同10kV母线引出,馈线之间末端通过联络开关两两相连。线路长度约3公里,采用架空裸导线,故障率取0.1次/(km·年),平均修复时间4小时。变电站出线断路器自动跳闸,分段隔离开关操作时间0.5小时,联络开关转供时间1小时。每个负荷点用户数从100户到300户不等,平均负荷0.2~0.5MW。

用上面的参数模拟20000年,并设随机数种子为固定值(比如2026),保证结果可复现。实际运算时间在普通笔记本上大概跑几十秒到几分钟,完全在可接受范围。

4.2 从参数设置到结果解读

参数文件直接写成Matlab脚本最方便,我通常单独建一个setup_test_system.m,内容和上面3.1节的结构体对应。运行主程序后输出的指标大致如下:

指标数值单位
SAIFI1.32次/用户·年
SAIDI7.65小时/用户·年
CAIDI5.80小时/次
ASAI99.913%
ENS238.4MWh/年

这个量级符合10kV架空配电网的实际情况。如果仿真里去掉联络开关转供,只靠馈线出线断路器和隔离开关,SAIDI会显著上升,因为故障点下游负荷只能等修复,时间长达4小时;增加了联络转供后,负荷点平均停电时间降下来不少。这也是为什么许多配电网改造项目的核心就是增加分段开关和联络线。

4.3 灵敏度分析,让结论更可信

项目报告里只放一套结果远远不够,通常还要做灵敏度分析。比如修复时间从4小时变成8小时,SAIDI和ENS几乎翻倍,但SAIFI不变;把联络开关操作时间从1小时压到0.5小时,CAIDI小幅下降,说明这个环节不是瓶颈;如果把线路故障率从0.1提到0.2次/(km·年),SAIFI、SAIDI、ENS整体都会上升约一倍。

这类分析对规划决策非常有用。因为它能告诉你:同样投资,是换电缆降低故障率更划算,还是增加联络开关提升转供能力更划算。你可以在代码里把这些场景写成批量循环,一次跑完自动生成对比表。

5. 常见问题与排查技巧实录

5.1 结果始终不收敛,到底是模拟年数不够还是代码有bug

我遇到过最典型的“假不收敛”:模拟了5、6万年SAIDI还在波动。排查后发现是随机数种子的影响——换了种子结果差异很大,这说明样本量确实不够,或者方差太大。但另一个更隐蔽的问题是对数计算时lambda传了0值,导致TTF变成无穷大,事件表里出现不合理的极大值,整个统计被污染。

建议从两个方向排查:第一,检查事件表里是否有明显不合理的时间值,比如几百万小时;第二,削减到单馈线小网络,跑一个手工可以算的案例,验证主逻辑没问题后再扩大规模。

5.2 可靠性指标怎么验证

写程序最大的风险是“一本正经地算出错误结果”。我验证程序可靠性的方法是和解析法对比一个小型系统:比如一个无限大电源带一条辐射馈线、三个负荷点,没有开关转供。这个简单网络的SAIFI和SAIDI理论上等于线路总故障率乘一些修正因子,手算小网络能算出来;代码跑出来的结果只要跟手算量级和趋势一致,就说明整体建模逻辑和FMEA正确。

5.3 Matlab性能优化:从“一天跑不完”到“半小时跑完”

序贯蒙特卡洛模拟天生慢,但在Matlab里优化空间很大。第一个方法是事件表用预分配内存,不要每循环一次就动态增长。第二个方法是把大量样本一次性向量化抽样,批量生成TTF、TTR,再把事件排序,避免在for循环里反复调用rand()。第三个方法是故障后果分析用稀疏矩阵和图形对象加速,不要用嵌套for循环逐节点搜索。

我印象最深的是把事件推进逻辑从“逐事件循环”改成“按故障区间批量推进”后,计算速度提升了近一个数量级。代价是代码复杂度上升,对于教学和一般工程项目,先保持简单清晰的版本,在确实需要提速时再优化。

5.4 负荷点用户数权重容易算错

SAIFI和SAIDI都是用户数加权指标。很多初学者直接把每个负荷点的停电次数取平均,忘记乘以用户数,导致结果整体偏小或偏大。我在代码里会专门维护一个userNum数组,统计时严格执行“用户数×停电次数”的汇总,再除以总用户数,每一步都单独输出检查。

6. 扩展方向:DG接入、储能和配电网新型态

6.1 分布式光伏接入后的时序效应

传统可靠性评估假设负荷恒定,真实世界的10kV馈线越来越多地接入分布式光伏。光伏出力有强烈的昼夜和季节特征,白天可能降低馈线负载从而改善电压、增大转供能力,但夜晚又完全无出力。序贯蒙特卡洛模拟法天然适合处理这种时序问题:把时钟推进到每个小时,光伏出力按实际日照曲线抽样,负荷按季节典型曲线变化,这样评估出的可靠性指标更贴近现实。

6.2 微电网孤岛运行的自愈能力

配电网加装储能和微电网控制后,外部故障时局部可以离网运行,即“孤岛模式”。序贯法的时序模拟可以精确刻画储能SOC在故障发生前的状态,进而判断孤岛能维持多久。这在非序贯法里是做不到的。很多做分布式能源规划的朋友来找我咨询,我一般建议直接从序贯法入手,把光伏、储能按时间步长建模,后续扩展很顺滑。

6.3 从“评估”到“优化”:启发式算法结合

有了可靠的单次评估器,就可以做配电网规划优化了。常见思路是用遗传算法、粒子群算法搜素网络拓扑或开关配置方案,把序贯蒙特卡洛模拟得到的可靠性指标当作适应度函数。虽然计算量大,但胜在灵活。我在实际项目里通常会把核心模拟函数封装成black-box接口,优化算法每调用一次就完成一次可靠性评估,整体循环可控。

如果做这块,建议先把模拟器性能优化到位,否则优化算法跑一组种群就要几天,项目根本不具备可行性。

最后分享一点个人体会

做这个项目踩过的坑很多,但最大的体会是:配电网可靠性评估的难处不在蒙特卡洛方法本身,而在对配电网运行逻辑的理解。你得先清楚断路器怎么跳、隔离开关怎么拉、联络开关什么时候合,然后才谈得上写程序。很多论文里的代码库,问题不在于抽样公式错,而在于故障后果分析过于理想化,把转供时间、开关失败概率这些实际因素都省略了。真正能落地的东西,恰恰是在这些细节里。

如果打算自己写,我建议先用一条馈线、两三个负荷点的小例子把主流程跑通,再慢慢加开关、加联络、加DG。代码组织上,数据结构设计要先想清楚,事件表驱动的主循环和故障后果分析函数分开写,这样调试起来方便很多。最后再把灵敏度分析和结果可视化补上,整个项目就可以拿去应付报告、论文或者实际规划项目了。

需要专业的网站建设服务?

联系我们获取免费的网站建设咨询和方案报价,让我们帮助您实现业务目标

立即咨询