做配电网规划、生产运维或者新能源接入评估的同行,应该对“可靠性评估”这四个字不陌生。不管是10kV馈线改造、分段开关位置优化,还是分布式光伏大规模接入后的影响分析,最终都得回答一个很朴素的问题:用户一年平均停几次电、每次停多久、一年里有多少时间用不上电。这几个问题的答案,就是SAIFI、SAIDI、ASAI这些指标。而算准它们,恰好是很多项目和论文的分水岭——解析法算起来快但拓扑一复杂就发怵,非序贯蒙特卡洛简单但没法还原时间轴上的出力与负荷变化。这段时间我把一套基于序贯蒙特卡洛模拟法的配电网可靠性评估流程用Matlab完整跑通了,从指标建模、元件状态抽样、故障影响分析到收敛判据和结果统计,每个环节都有不少具体实现细节和能提前避开的坑。这篇就把整套思路和关键代码实现逻辑整理出来,给打算用Matlab做可靠性评估研究或工程评估的读者一个可以直接参考的框架。
1. 配电网可靠性评估的完整思路:指标口径与方案选择
1.1 四个绕不开的可靠性指标
配电网可靠性评估,本质是把“供电质量”量化成几个可计算的统计量。最常用的是SAIFI、SAIDI、CAIDI、ASAI,再加一个ENS(缺供电量)。先说清楚计算口径,后面代码才能写对。
| 指标 | 计算式 | 物理含义 |
|---|---|---|
| SAIFI | 用户停电次数总和 / 用户总数 | 平均每户每年停电几次,单位次/(户·年) |
| SAIDI | 用户停电时户数总和 / 用户总数 | 平均每户每年累计停电多少小时,单位h/(户·年) |
| CAIDI | SAIDI / SAIFI | 平均每次停电持续多久,单位h/次 |
| ASAI | 1 - SAIDI / 8760 | 供电可用率,一般用小数表示,如0.9996 |
| ENS | Σ(失负荷功率 × 停电时长) | 系统全年缺供电量,单位kWh |
这几个指标里,SAIFI和SAIDI是核心输出,CAIDI只是两者相除,ASAI本质上是SAIDI的另一种表达。做可靠性研究时,电力公司考核指标一般挂在SAIDI和ASAI上,所以代码里这两个一定不能算错。
工程上统计时还要区分“故障停电”和“预安排停电”。如果只评估设备故障对供电可靠性的影响,预安排停电不纳入;如果是综合评估,两个都要算。我建议一开始先只做故障停电,把框架跑通后再加计划检修模块,否则模型细节一多,连自己都分不清指标偏差到底是模型问题还是逻辑问题。
1.2 为什么解析法在配电网面前经常“算不动”
可靠性评估有两类主流方法:解析法和蒙特卡洛模拟法。解析法通过枚举故障事件、建立故障模式后果分析表来计算指标,经典的最小割集法、故障树法都属于这一类。对输电系统或简单辐射网络,解析法很快,公式一套就能出结果。
但配电网一旦复杂起来,解析法就头疼了:多分段多联络的拓扑结构、分段开关和联络开关的配合逻辑、继电保护的动作范围、分布式电源和储能带来的时序行为……这些因素组合在一起,枚举所有可能的系统状态会呈指数爆炸。更麻烦的是,DG出力和负荷是随时间变化的量,解析法很难把“上午十点的光伏出力”嵌入到一年8760小时的积分里去。所以工程上做配电网评估,更多转向模拟法。
1.3 序贯蒙特卡洛与普通蒙特卡洛的区别在哪
蒙特卡洛也分两类。非序贯抽样版的思路是:直接按元件的故障概率抽一个系统状态,比如抽到某个元件故障,就分析这个状态下的失负荷情况,把海量抽样结果求平均。它速度快、实现简单,但有先天缺陷——它算不出“频率”的时间语义。SAIFI、SAIDI要求知道故障何时发生、持续多久,非序贯法只能估计概率性的时间占比,想获得频率指标必须额外做转换。另一大缺陷是它无法处理时序特性:光伏中午出力最大、夜间为零,负荷昼夜峰谷不同,这些信息非序贯法根本没用上。
序贯蒙特卡洛(Sequential Monte Carlo,也叫时序蒙特卡洛)的思路就直白得多:既然真实系统是“年复一年、每时每刻”运行的,那就把一年8760小时、甚至几千上万个8760小时直接模拟一遍。对每个元件抽样它的故障时刻和修复时刻,形成一条“运行—故障—修复—运行”的状态时间轴;多个元件的时间轴交织在一起,系统就随着时间推进经历各种状态。模拟结束后,统计所有停电事件,SAIFI、SAIDI这些指标自然就出来了。这就是序贯法相对普通蒙特卡洛的核心价值,也是为什么项目标题里明确把“序贯”两个字点出来的原因。
2. 序贯蒙特卡洛模拟法的关键机制:抽样、事件推进与故障分析
2.1 元件状态持续时间抽样:TTF和TTR的数学本质
配电网里的大多数元件(馈线段、变压器、开关)在可靠性评估里都抽象成“运行/故障”两状态模型。元件故障率λ(次/年),平均修复时间r(小时)。若假设故障发生过程为泊松过程,元件正常运行持续时间TTF服从指数分布,抽样公式是:
TTF = -ln(U) / λ
其中U是(0,1)均匀分布随机数。这个公式的原理可以这样理解:指数分布的概率积分变换,把均匀随机数映射到指数分布样本,ln括号前的负号保证结果为正。修复持续时间TTR也常按指数分布抽样,但工程上很多人直接用平均修复时间r来近似,因为修复过程受备件、交通、天气影响,并不严格服从指数分布。实测下来,对于中压配电网的简单评估,TTR取固定均值的误差通常可以接受,但前提是全程都取固定值,不能一会儿指数抽样一会儿取均值,那样统计口径就乱了。
这里有一个新手极容易忽略的坑:单位统一。TTF抽样里λ按“次/年”,TTF单位就是年;TTR按“小时”给,两个量加在一起时如果不换算,时间轴直接错乱。我习惯把所有时间统一折算成小时,λ也换算成“次/h”,这样整个模拟循环里不需要中途再做单位转换,出错概率小很多。
2.2 事件驱动推进与故障影响分析
序贯MC的代码核心是一个“下一次事件”循环。我给每个元件记录它下一次状态变化的时刻,流程是这样的:
- 初始化:每个元件处于运行态,抽样各自首次故障时间TTF;
- 找到所有元件中TTF最小的那个,作为下一个事件时刻;
- 推进系统时间到该时刻,触发该元件的故障或修复事件;
- 如果是故障事件:把元件状态置为故障,做故障影响分析,统计失电用户数和失负荷量,然后为该元件抽样修复完成时间TTR;
- 如果是修复事件:把元件状态置为运行,重新抽样下一次故障时间TTF;
- 重复执行,直到模拟时间到达设定年限。
故障影响分析是仿真的“重心”,它决定哪些用户这次停电、停多久。对一个辐射状配电网络,故障发生后首先要区分三块区域:
- 故障点上游:由变电站母线供电的区域,一般不受影响;
- 故障点下游、可通过分段开关隔离并转供的区域:只停“倒闸操作时间”,比如几十分钟;
- 故障点下游、不可转供的区域:停到故障修复完成为止。
这个逻辑在Matlab里通常用图搜索实现。把配电网抽象成图,故障元件对应的边删除后,从电源节点做可达性搜索,能到达的节点视为有电;再逐个检查联络开关路径和容量约束,放宽后得到失电节点集合。如果忽略了联络开关和转供容量,SAIDI会偏高;但如果默认所有下游都能转供,结果又会偏乐观。实际代码里建议把“是否考虑转供”做成一个开关参数,先跑通纯辐射状最简情形,再打开联络开关。
2.3 收敛判据:模拟多少年才够
模拟不是跑固定年数就行,而是要等指标估计值收敛。常用判据是方差系数β:
β = σ / (μ · √N)
其中σ是样本标准差、μ为指标均值、N为模拟年数。工程上一般要求SAIDI和SAIFI的β小于0.02~0.05才可以停止。很多人第一次跑序贯MC会惊讶:“怎么要模拟几千上万年?”但仔细想想就明白了,元件故障率是每年0.1次量级,很多年份压根不停电,样本方差天然很大,收敛确实慢。这不是Matlab实现的问题,而是方法本身的统计代价。
| 仿真年限 | β估计量级 | 结论 |
|---|---|---|
| 500年 | 0.3~0.5 | 波动很大,不可用 |
| 2000年 | 0.1左右 | 趋势可见,精度不足 |
| 10000年 | 0.02~0.05 | 工程可用 |
| 50000年 | 0.01以下 | 研究级精度,计算量巨大 |
另有一个容易被忽略的细节:模拟开始前要设置预热(burn-in)阶段。初始状态全部元件可用,如果从零时刻开始统计,前一段的可靠性会被高估。我一般在总仿真时间的前5%~10%不参与统计,代码里用一个逻辑标志控制统计起点,而不是跑完再截断,这样内存更友好。
3. Matlab代码实现:从零搭建序贯MC评估程序
3.1 输入数据与拓扑建模
Matlab代码的第一步是把配电网拓扑和参数组织成数据结构。我推荐用矩阵和稀疏邻接表组合的方式:
branch矩阵:每行一个馈线区段,包含首节点、末节点、区段长度、单位长度故障率、平均修复时间;loadData数组:每个节点的用户数和平均负荷(功率);- 联络开关、分段开关位置,以及转供容量上限;
- 电源节点编号。
示例定义如下:
% 节点数N,区段数M % [首节点, 末节点, 长度km, 故障率次/(km*年), 修复时间h, 是否可控开关] branch = zeros(M, 6); % [用户数, 平均负荷kW] loadData = zeros(N, 2); sourceNode = 1; % 电源点节点编号把所有时间量统一为小时,故障率统一换算成“次/h”:
lambda_h = branch(:,4) / 8760 .* branch(:,3);这个换算看似简单,实际犯错的概率很高。曾见过第一版程序里TTF按天、TTR按小时,算出来的SAIDI直接飙到几百小时,后来逐项排查才发现是单位混用。
3.2 主循环与关键函数设计
代码主体分三块:抽样模块、事件推进模块、故障分析模块。核心主循环用Matlab风格写成下面这样:
% 参数准备 NumElements = size(branch,1); State = zeros(NumElements, 1); % 0运行, 1故障 NextTime = -log(rand(NumElements,1)) ./ lambda_h; % 首次故障时间 Tsim = SimYears * 8760; % 总仿真时间(h) burnStart = Tsim * 0.05; % 预热期 t = 0; % 统计量 FailCount = zeros(N, 1); % 节点停电次数 FailDur = zeros(N, 1); % 节点累计停电时间(h) FailEnergy = 0; % 累计缺供电量(kWh) while 1 [t, k] = min(NextTime); % 最早事件时刻及元件编号 if t > Tsim break; end if State(k) == 0 % ----------------- 发生故障 ----------------- State(k) = 1; if t > burnStart [lostNodes, lostPower] = FaultImpactAnalysis(k, branch, loadData); FailCount(lostNodes) = FailCount(lostNodes) + 1; FailDur(lostNodes) = FailDur(lostNodes) + repair_time(k); FailEnergy = FailEnergy + lostPower * repair_time(k); end % 修复完成后,该元件下次事件为修复完成 NextTime(k) = t + branch(k,5); % 固定修复时间 else % ----------------- 修复完成 ----------------- State(k) = 0; % 重新抽样下一次故障时间 NextTime(k) = t + (-log(rand) / lambda_h(k)); end end这段代码的关键在于min(NextTime)找最早事件,通过State(k)区分故障和修复两种分支。真正项目里,我还会把“当前系统时间是否要跨年”单独处理,因为如果需要统计逐年指标,跨年时的节点需要拆开累计。
修复时间这里用了固定值,如果要指数抽样,把branch(k,5)替换成-log(rand) * branch(k,5)即可。两种写法都行,但必须在同一套代码里统一。
3.3 故障影响分析的实现思路
FaultImpactAnalysis函数是影响范围判断的核心。最简单可靠的方式是删除故障边后做BFS可达性搜索:
function [lostNodes, lostPower] = FaultImpactAnalysis(faultBranch, branch, loadData) % 构建邻接矩阵,把故障边临时移除 % 从电源节点出发做BFS,标记能到达的节点 % 不可达节点即为失电节点 % 再考虑联络开关转供:如果某失电区域能通过联络路径可靠转供, % 且转供容量满足,则从失电集合中移除,停电时长改为倒闸时间 end如果网络规模不大,几十个节点,Matlab自带的graph对象配合bfsearch函数非常方便。手写邻接矩阵+队列也很快,而且不容易被对象封装的细节坑到。我做大规模测试时,会优先写BFS的可达性版本,把转供逻辑作为一个可选开关,先不做任何转供,验证纯辐射网结果的合理性,再打开转供逻辑逐步加复杂度。
3.4 指标统计与一个33节点算例的仿真效果
统计完成后,把累加量换算成网级指标,这一步最容易出错的是SAIFI的分子。注意FailCount是每个节点的停电次数,最终统计要乘以该节点用户数再求和:
totalCustomer = sum(loadData(:,1)); SAIFI = sum(FailCount .* loadData(:,1)) / totalCustomer; SAIDI = sum(FailDur .* loadData(:,1)) / totalCustomer; CAIDI = SAIDI / SAIFI; ASAI = 1 - SAIDI / 8760;用一个经典的33节点辐射状配电网做测试(这类算例在很多教科书和论文中都有,参数大同小异),馈线故障率按0.1~0.2次/(km·年)、修复时间4~6小时设置,每个负荷节点用户数按节点功率折算。在SimYears=5000~10000年、预热期5%的设置下,常见结果区间是:SAIFI约0.5~1.5次/(户·年),SAIDI约2~6小时/(户·年),ASAI大于99.9%。如果你的结果跑出来差数量级,大概率不是算法复杂性问题,而是参数、单位或统计公式的问题,先回查这三处。
4. 踩坑实录:收敛慢、结果不稳、逻辑出错的排查技巧
4.1 收敛慢的五种诱因与加速方案
序贯MC最让人烦躁的就是“跑不完”。这里列我实际踩过的问题和对应处理方案:
| 常见问题 | 表现 | 解决方法 |
|---|---|---|
| 固定仿真年数不够 | 结果每次跑出来差很大 | 改用方差系数β作为停止判据 |
| 随机数序列质量差 | 指标波动大、不平稳 | 用rng(seed)固定种子,或改用Sobol低差异序列 |
| 重复做故障影响分析 | 相同故障区段反复计算 | 对故障区段的影响集合做缓存,映射复用 |
| 元件数量多、模拟年数大 | 整体耗时不可接受 | 按母线和馈线类型做元件聚合 |
| 单线程跑 | CPU利用率低 | 用parfor按年拆分并行,注意随机流隔离 |
收敛慢的背后是方差大。如果项目允许,优先使用低差异序列替代普通rand,实测下来同样精度下收敛速度能提升20%到50%。低差异序列在Matlab里的实现也不复杂,几十行内就能嵌入现有抽样逻辑。
4.2 故障影响分析结果偏小的隐藏原因
SAIDI和SAIFI偏低时,先别急着调参数,最容易出问题的地方藏在影响分析逻辑里:
- 默认所有故障都能通过联络开关转供,忽略了转供容量约束。实际配电网中联络线容量有限,重负荷区域根本转带不了;
- 忽略了保护配合关系。支路故障本该由支路开关切除,如果模型里一跳闸就动主断路器,影响范围被明显扩大;
- 分布式电源孤岛成功概率被估计过高,没有考虑逆变器故障穿越能力和孤岛功率平衡;
- 只看设备故障,没把计划检修叠加进去。真实系统中很多停电是预安排停电和故障停电共同作用的结果。
每种情况都要对照实际网络参数逐一确认,不要为了“指标好看”而偷偷调系数,这会让评估结果完全失去工程参考意义。
4.3 随机数与统计逻辑的低级错误清单
这部分都是“不报错但结果全错”的隐性坑,整理成清单供自查:
rand()返回0导致-log(0)产生Inf,时间轴直接崩溃。建议对抽样加一个下限保护,比如U = max(rand, 1e-10);- 没有设置随机种子,两次仿真结果不可复现,论文里的图没法稳定复现;
- 统计起始点没有跳过预热阶段,前几百年的“完全正常供电”把可靠性指标拉高;
- SAIFI分子忘记乘以用户数,直接用停电次数除以总用户数,算出来的值小到离谱;
- 修复时间一会儿用指数抽样、一会儿用固定均值,统计口径完全混乱;
- 转供后的负荷点停电时长仍然用修复时间,导致SAIDI偏大。转供场景下停电时长应该是倒闸操作时间,不是修复时间。
调试时有个笨办法但非常好用:先搭一个三节点的极简网络,手算出SAIFI和SAIDI的预期值,再跑程序对比。对得上再上大算例。这一步能过滤掉绝大多数低级错误,省下的排查时间远大于搭建小样例的时间。
5. 这套代码还能往哪儿扩展
5.1 分布式电源与时序负荷建模
序贯蒙特卡洛最有价值的扩展方向就是接入时序模型。光伏出力按8760小时曲线生成,晴雨天用云量系数修正;负荷用典型日负荷曲线描述;储能按充放电策略建模。具体做法是在主循环里增加一个“当前时刻”查询函数,根据时刻查表得到该小时DG出力和负荷倍率,再叠加到故障影响分析里。
DG接入后的一个重要场景是孤岛:主网故障时段内,某区域由DG独立供电,该区域用户不停电或只有短时停电。这个场景在非序贯蒙特卡洛里很难表达,但序贯MC只需在故障影响分析里加一个“是否能形成稳定孤岛且功率平衡”的判断。判断逻辑不复杂,却能让评估结果明显更接近实际——尤其当馈线上光伏渗透率超过30%时,忽略孤岛会让SAIDI严重高估。
5.2 与规划优化算法结合
可靠性评估很少是最终目标,更多时候是作为评估模块嵌入网架规划、分段开关优化、检修策略优化。比如用遗传算法迭代时,每一代都要调用序贯MC评估候选方案,这会带来巨大的计算量。我处理这类问题有三个基本策略:
- 方案差异不大时只评估变化区域的元件,做增量评估,而不是每次全网络重跑;
- 把MC评估函数写成
parfor并行版本,或者对热点区域单独编译成mex函数; - 对相同拓扑的故障影响集合做缓存,命中后直接复用,避免反复BFS。
整体思路是“评估结果准确”优先,性能优化在准确性验证通过后再做。很多人一上来就优化性能,结果指标口径错了还浑然不知,最后优化得越快,错得越离谱。
最后再说两句
这套序贯蒙特卡洛配电网可靠性评估程序完整跑通之后,我最大的体会是:可靠性指标看着只是几个数字,真正让它靠谱的是每一步的物理含义都讲清楚。抽样公式里U为什么不能取0,SAIFI分子为什么必须按时户次统计,故障分析里转供容量到底怎么约束——每一个细节都会直接反映在最终指标上。用Matlab做实现的读者,我建议先在一个10节点以内的小网络上手动算一遍预期指标,再上大算例,调试效率会高很多。等框架稳定之后,再往里面加DG、储能、自动开关这些扩展模块,改动量其实不大。也正是因为这个框架足够模块化,后面好几个衍生研究都是从这一版代码上长出来的。