做配电网可靠性评估最头疼的往往不是指标定义本身而是面对一个实际的中压馈线系统时用解析法算到怀疑人生或者用非时序抽样算完心里完全没底。我自己在做一个带分布式电源和联络开关的10kV馈线可靠性项目时被故障模式影响分析法FMEA的枚举量折磨过之后彻底转向了序贯蒙特卡洛模拟法。这个方法最吸引我的地方是它能保留完整的时间序列信息天然支持时变负荷、储能充放电策略、分布式电源出力波动这些复杂场景而且只要代码框架搭好换网络拓扑只是改数据表的事。这篇内容适合三类人一是刚接触配电网可靠性评估、想搞懂序贯蒙特卡洛到底怎么落地成Matlab代码的研究生二是做配电网规划的工程师想把手头的可靠性计算从静态校核升级到时序仿真三是想在论文里复现可靠性指标的科研人员。我会把元件时序状态建模、故障影响分析、指标统计、收敛判据、代码骨架和实际踩过的坑全部都拆开讲。1. 为什么配电网可靠性评估最终选了序贯蒙特卡洛1.1 解析法算不动的那几类场景配电网可靠性评估的传统路径是解析法典型代表是FMEA和最小割集法。对于简单的辐射状配电网FMEA的思路很直接逐条枚举馈线段的故障分析每个故障对下游负荷点的影响再叠加得到指标。它的优点是计算快、结果确定一套网络参数算完不会变。但麻烦也在这里——一旦系统规模变大或者出现多电源转供、分段开关联动、分布式电源孤岛运行、储能参与黑启动这些时序相关策略枚举的故障场景会呈爆炸式增长。我见过一份40个节点的网络光考虑两两组合故障就列了上千行FMEA表还没算分布式电源出力变化的时段。另一个问题是解析法很难处理时变负荷。可靠性指标比如ENS期望缺供电量本质上是功率对时间的积分如果负荷曲线是季节性的、峰谷分时的解析法的解析表达式会变得非常繁琐。你可能需要把一年切成几十个时段每个时段重新做一次系统状态分析。实际测试下来发现这样做精度上去了但代码复杂度也上去了维护起来很痛苦。1.2 序贯蒙特卡洛与非序贯、解析法的边界蒙特卡洛模拟法分为非序贯和序贯两种。非序贯蒙特卡洛的做法是从各元件的概率分布中直接抽样一组状态比如正常运行或者故障形成系统当前时刻的“快照”然后对这个快照做静态影响分析统计指标。这种方法计算效率高适合大规模系统的年度指标估算但丢失了状态之间的时间关联无法回答“停电发生在几点”“两次停电之间隔了多久”这类时序问题。序贯蒙特卡洛则是按照时间顺序模拟每个元件从正常运行到故障、再到修复、再故障的完整生命周期然后把所有元件的时间轴对齐形成系统层面的状态转移序列。它的核心优势恰恰是保留了“状态持续多长时间”的信息。这意味着你可以在仿真过程中叠加负荷曲线、分布式电源出力曲线、储能SOC变化曲线甚至评估不同控制策略对可靠性指标的影响。代价也很明显计算量大收敛速度通常慢于解析法单次仿真结果有波动必须依赖统计收敛判据。1.3 这个方法的代价与适用前提不要以为用了序贯蒙特卡洛就万事大吉。它的代价有两个第一是计算耗时尤其当系统规模大、元件多的时候一个10万年的仿真可能需要几十分钟甚至更久第二是结果不是精确值而是带置信区间的估计值。所以适用前提是你确实需要时间序列信息或者解析法已经复杂到不可维护。如果只是算一个纯放射状、无联络、负荷恒定的网络用解析法反而更合适。我在项目中采用序贯蒙特卡洛是因为需要对比“分段开关加联络开关”策略在不同负荷水平下的可靠性改善效果这种场景下时序仿真几乎是唯一的选择。2. 时序状态模型元件怎么“坏”系统怎么“转”2.1 可靠性参数与指数分布抽样原理序贯蒙特卡洛的第一步是给每个元件建立状态持续时间模型。工程上最常用的是两状态模型正常运行和故障修复。对应两个随机变量元件正常运行持续时间TTFTime To Failure和故障修复时间TTRTime To Repair。元件可靠性参数通常由历史统计得到包括年平均故障率λ次/年和平均修复时间MTTR小时。在指数分布假设下故障率恒定TTF服从参数为λ的指数分布其概率密度为f(t)λe^(-λt)期望E[TTF]1/λ单位是年TTR通常也假设服从指数分布期望E[TTR]MTTR单位是小时。实际统计中TTR可能更接近对数正态分布但为了简化计算多数文献和代码仍采用指数分布误差在可接受范围。指数分布抽样非常简单用的是逆变换法。对于累积分布函数F(x)1-e^(-λx)令U为一个[0,1]区间均匀分布的随机数则X-ln(1-U)/λ就服从此指数分布。由于1-U与U同分布可直接写成X-ln(U)/λ。Matlab里rand函数生成的是(0,1)开区间的均匀随机数正好可以直接用。写成匿名函数就是ttf_sample (lambda_year) -log(rand) / lambda_year * 8760; % 单位小时 ttr_sample (mttr_hour) -log(rand) * mttr_hour; % 单位小时这里有个易错点如果λ的单位是“次/年”那么TTF抽样的自然单位是年要乘以8760换算成小时如果λ直接给的是“次/小时”比如λ0.065/8760那抽样结果自然就是小时。我在代码里统一用小时做仿真时间单位避免后面做时间轴运算时单位混乱。2.2 网络拓扑与故障影响分析上游隔离、下游转供配电网不同于输电网绝大多数工况下是辐射状运行。这意味着潮流方向是单向的故障影响范围可以基于“源-荷”路径关系判断。我在代码里把网络表示成了一张带父子关系的树状结构每个馈线段或者叫支路都有一个首端节点和一个末端节点从电源点开始按拓扑顺序编号。故障影响分析的核心逻辑可以概括成几句话故障发生在某个馈线段k上时从该馈线段到电源点之间的负荷即上游负荷在完成了故障隔离之后可以由主电源恢复供电。这一段停电时间只包括故障定位、隔离和恢复操作的时间一般是0.5到1小时而不是整个修复时间。故障点之后、与该馈线段直接相连的下游负荷如果存在联络开关且转供容量足够则停电时间为“隔离合联络开关”的操作时间也就是通常说的转供时间。如果没有联络开关或者联络开关的容量不够那么这些下游负荷只能等故障元件修复完毕停电时间是修复时间TTR。如果故障发生在某个分支线比如带熔断器的支路且熔断器能正确动作那么主馈线上的其他负荷不受影响只是该分支负荷停电。为了在Matlab里实现这个逻辑我维护了一个负荷点上游元件的列表做法是对每个负荷点反向搜索父节点直到电源节点。这个搜索在预处理阶段做一次之后运行仿真时直接查表开销很小。具体代码框架在后面第3章展开。2.3 负荷时序与故障恢复策略的建模序贯蒙特卡洛的另一个好处是可以把负荷建模成随时间变化的曲线。最简单的做法是给每个负荷点配一条24小时归一化负荷曲线再乘以该负荷点的年峰值负荷。仿真时钟按小时推进或按事件推进时提取对应小时的负荷值用于停电电量计算。这里有个细节按事件推进时如果一次停电跨过了多个小时就需要对停电时间段内的每个小时累计缺供电量。我在程序里是这么处理的记录停电开始时刻t_start停电结束时刻t_end然后在这段时间内按小时采样负荷曲线逐小时累加功率与1小时的乘积。虽然这会多一层循环但计算结果更接近真实。故障恢复策略也就是转供策略跟网络结构强相关。如果要做精细化需要把联络开关的额定容量和下游负荷的峰值都纳入判断转供容量不足时只能转供一部分负荷。这个可以在故障影响分析函数里加一段负荷约束判断。对于一般的教学和研究场景可以先假设联络开关容量足够跑通主流程再加约束。我自己的做法是先不加约束跑一版结果再逐步增加容量约束对比不同转供能力对指标的边际影响。3. Matlab代码框架从数据准备到主循环3.1 输入数据怎么组织元件表、拓扑表、负荷表写Matlab代码之前最关键的是把数据结构设计清楚。我习惯用三种表结构体数组或table类型都可以branch表存馈线段load表存负荷点link表存联络开关。branch表包括字段编号、首端节点、末端节点、长度km、单位长度故障率次/km·年、修复时间小时/次、线路类型主干/分支、末端是否有分段开关。这里单位长度故障率要和长度相乘得到该元件的年故障率λ。举个例子一段2km的架空馈线单位长度故障率0.065次/km·年那么该段λ0.13次/年。load表包括字段负荷点编号、挂接节点、用户数、年平均负荷kW、24小时归一化负荷曲线。归一化负荷曲线其实是个24维数组方便程序里按小时索引。link表包括字段联络开关编号、连接的首端节点和末端节点、额定容量kVA、闭合操作时间小时。这一表是可选的没有联络就留空。拓扑关系的核心是两个数组parent_node和child_list。我通过遍历branch表从电源点出发按广度优先建立父子关系。这个关系是后面故障影响分析的依据。还有一个预处理是建立“每个负荷点的上游branch集合”这个也可以用反向查父节点做。3.2 序贯蒙特卡洛主循环的代码骨架下面给出一个可运行的最小版本主循环骨架。这里采用最常用的“状态持续时间抽样法”即对每个元件抽样一系列“运行—修复”状态持续时间沿时间轴组合成系统状态转移序列。为了方便展示我做了两个简化假设元件修复期间其他元件不发生新故障仿真年数内不考虑设备老化引起的故障率变化。这两个假设在常规计算中是可接受的后面的进阶说明里我会讲怎么去掉第一个假设。%% 序贯蒙特卡洛主循环简化版 rng(2024); % 固定随机种子结果可复现 % 元件参数 n_branch size(branch, 1); lambda branch.lambda_per_km .* branch.length; % 次/年 mttr branch.mttr_hour; % 小时 % 仿真参数 simu_years 50000; % 仿真年限 hours_per_year 8760; n_load length(load); % 统计量初始化 LP_outage_count zeros(n_load, 1); % 停电次数 LP_outage_duration zeros(n_load, 1); % 停电小时数 LP_ens zeros(n_load, 1); % 缺供电量 kWh for yr 1:simu_years % 首次故障时间抽样小时 ttf -log(rand(n_branch, 1)) ./ lambda * hours_per_year; sim_time 0; while sim_time hours_per_year % 找到最早发生故障的元件 [t_min, k] min(ttf); if sim_time t_min hours_per_year break; end sim_time sim_time t_min; % 抽样该元件的修复时间 ttr -log(rand) * mttr(k); % 故障影响分析累加停电统计 [count_k, dur_k, ens_k] outage_analysis(k, sim_time, ttr); LP_outage_count LP_outage_count count_k; LP_outage_duration LP_outage_duration dur_k; LP_ens LP_ens ens_k; % 修复完成后重新抽样该元件的下次故障时间 ttf(k) -log(rand) / lambda(k) * hours_per_year; % 注意其他元件的故障时间保持不变 end end % 可靠性指标计算 SAIFI sum(LP_outage_count .* load.user) / sum(load.user); SAIDI sum(LP_outage_duration .* load.user) / sum(load.user); CAIDI SAIDI / SAIFI; ASAI 1 - SAIDI / (simu_years * hours_per_year); ENS sum(LP_ens);这里最小事件时间的选择用min(ttf)就可以每次只推进到最早故障事件。每次修复完成后只重新抽样故障元件的TTF其他元件维持原来的剩余TTF这正是事件导向仿真的标准做法。需要指出的是第一次进while循环时其他元件的TTF是完整抽样值随着每次事件推进这些TTF会自然老化因为它们是基于绝对时间抽样的。每次事件后的ttf数组不用全部减去t_min因为只用最小值决定下个事件绝对时间和剩余时间是等效的。如果你习惯用剩余时间建模也可以在每个循环里让ttfttf-t_min结果等价但少了随机数调用速度略快。3.3 故障影响评估函数怎么写故障影响分析函数是整个代码里最核心的部分它决定指标算得对不对。我的outage_analysis函数输入是故障元件编号k、故障发生时刻sim_time和修复时长ttr输出是本轮故障对每个负荷点累计的停电次数向量、停电时长向量和缺供电量向量。实现逻辑分为几步第一步根据branch的父子关系找到该故障元件对应的子树范围。具体做法是从故障元件末端节点往下遍历找出所有受影响的节点集合。第二步对所有负荷点判断该负荷点是否在受影响的子树上。这里我用一个n_load×n_branch的逻辑矩阵load_upstream_precompute提前算好每个负荷点的上游branch集合这样在仿真循环里只需要做逻辑索引不用重复遍历。第三步对每个受影响负荷点根据它的位置决定停电时间如果负荷点在故障元件的上游停电时间为ISOLATION_TIME比如0.5小时如果负荷点在故障元件下游但通过联络开关转供成功停电时间为ISOLATION_TIMESWITCH_TIME比如0.51小时如果负荷点在故障元件下游且无转供或转供容量不足停电时间为ttr。如果负荷点在故障元件所在的区段内比如故障元件就是负荷点直接挂接的支路停电时间为ttr。第四步计算缺供电量。缺供电量等于停电时间内每个小时的负荷功率之和。如果停电时间段跨越了峰值时段ENS会显著增大这就是时序仿真与恒定负荷解析法的差别所在。为了方便阅读我把这一步写成一个独立函数函数的头部大致是这样function [count_vec, dur_vec, ens_vec] ... outage_analysis(k, sim_time, ttr, branch, load, upstream_map)由于upstream_map在仿真开始前已经算好这一步就是个查表和累加的过程计算速度很快。我实际测试过一个60节点的馈线系统5万年的仿真时长主循环部分跑下来大约40秒瓶颈主要在Matlab的循环本身而非算法。如果你对性能有更高要求可以改成MEX或者vectorized聚合不过对一般研究用途这个速度已经足够。4. 可靠性指标怎么统计从负荷点到系统4.1 负荷点指标与系统指标定义配电网可靠性指标分负荷点指标和系统指标两层。负荷点指标是基础系统指标是所有负荷点指标的汇总。负荷点指标常用的有三个平均故障率λ_LP次/年、平均停电持续时间r_LP小时/次、年平均停电时间U_LP小时/年。三者满足U_LP λ_LP × r_LP。比如某负荷点一年平均故障0.3次每次平均停电2.5小时那么年平均停电时间就是0.75小时。系统指标则是加权汇总。其中最重要的四个SAIFI系统平均停电频率指标每个用户每年平均停电次数单位次/用户·年。计算方式是所有负荷点的停电次数乘用户数之和除以总用户数。SAIDI系统平均停电持续时间指标每个用户每年平均停电小时数单位小时/用户·年。计算方式是所有负荷点的停电时长乘用户数之和除以总用户数。CAIDI用户平均停电持续时间指标单位次停电的平均持续时间等于SAIDI除以SAIFI单位小时/次。ASAI供电可用率一年内实际供电小时数占全年小时数的比例等于1减去SAIDI除以8760。ENS期望缺供电量和AENS平均系统缺供电量则从能量维度反映可靠性计算方式是把停电期间的用户负荷逐小时累加。这些指标各有侧重SAIFI和SAIDI关注频次和时长ASAI关注整体可用率ENS关注经济损失。工程做规划时往往四个都要看比如某个改造方案可能让SAIFI没变但SAIDI显著下降这代表故障次数没少但每次停电时间因为转供能力增强而缩短了这也是有价值的改善。4.2 仿真钟推进中的统计逻辑统计逻辑的要点是“停电事件归因”。序贯蒙特卡洛里一次系统状态转移可能同时影响多个负荷点但每个负荷点的停电时长可能不同。我实现的方法是在outage_analysis函数里直接对每个负荷点做独立累加返回的是向量而不是标量。这样主循环里就是简单的向量加法。这里有一个细节需要仔细处理如果同一负荷点一年内被两个不同的故障元件影响出现两次停电那它的停电次数就该加2停电时长累加两次。由于每个故障元件是独立事件这个累加逻辑是自然的。麻烦的是如果两个故障元件在同一时间段内同时故障重叠故障负荷点的停电时段可能重叠这时需要避免重复计数。在简化假设修复期间其他元件不故障下不会出现重叠所以我在主循环里的统计逻辑是干净的。如果做严格模型需要把重叠故障根据实际恢复时间重新划分停电区间这个我在后面的进阶说明里展开。4.3 收敛性判断与仿真年限设置蒙特卡洛最核心的问题不是“跑多少年”而是“结果收敛了没有”。我见过一些论文正文写着仿真10万次但没有任何收敛判据这个审稿人一问就得露馅。常用的收敛判据是相对方差系数β定义为估计值的标准差除以均值再除以样本量的平方根β σ_hat / (μ_hat × √N)其中σ_hat和μ_hat是已仿真样本的方差估计和均值估计N是仿真年数。工程上一般要求β小于0.05科研论文里更严格的要求是小于0.03甚至0.01。需要每跑一段时间就检查一次SAIFI或SAIDI的β达到要求就提前终止可以节省大量计算时间。我在Matlab里每1000年检查一次记录所有关键指标的均值序列观察是否进入平稳波动。实际经验是SAIFI通常收敛比较快因为故障次数本身就是个明显的计数SAIDI和ENS收敛会慢一些因为少数长时停电会在后期把均值拉高导致方差波动。如果你的系统里有低频率但长修复时间的元件一定要多跑几年再下结论。我碰过硬要省时间只跑1万年结果SAIDI的β还在0.2左右根本没有统计意义。当时就吸取了教训——仿真年限不是拍脑袋定的是收敛判据说了算。5. 算例验证一个可复现的小型配电网5.1 算例网络结构与参数为了验证代码正确性我给出一套简单但完整的算例读者可以直接照着搭。网络结构是“两段式馈线加联络”10kV母线引出主干馈线L1和L2两条馈线之间在末端通过联络开关连接正常运行时联络开关断开形成两个独立的辐射状馈线。馈线L1上挂三个负荷点LP1、LP2、LP3用两个分段开关分成三段馈线L2上挂两个负荷点LP4、LP5分为两段。每个负荷点通过一段长度为0.1km的支线接入主干馈线支线装有熔断器。具体参数如下元件长度(km)故障率(次/km·年)修复时间(h)备注L1段12.00.0655主干L1段21.50.0655主干L1段31.00.0655主干L2段12.50.0655主干L2段21.80.0655主干各支线0.10.0653熔断器保护负荷参数每个负荷点用户数200户峰值负荷300kW归一化负荷曲线采用典型居民日负荷曲线。分段开关的隔离操作时间0.5h联络开关的闭合操作时间1h。5.2 仿真结果与分析用上述参数跑5万年仿真固定随机数种子为2024得到典型结果如下同一个网络参数不同随机种子会有微小波动但应落在统计误差范围内指标无联络开关有联络开关SAIFI (次/用户·年)0.2850.281SAIDI (小时/用户·年)1.430.97CAIDI (小时/次)5.023.45ASAI0.9998370.999889ENS (MWh/年)0.820.55从结果可以看到联络开关对SAIFI的影响很小因为故障频次本身不因联络而改变但对SAIDI的改善非常明显从1.43小时降到0.97小时降幅约32%。原因很直观接入联络后L1下游区段和L2下游区段的负荷在故障时可以通过对侧电源转供停电时间从修复时间5小时缩短到转供操作时间1.5小时左右。这正是序贯蒙特卡洛方法最擅长体现的恢复策略价值——解析法算这类问题也不是不行但要处理转供路径和容量约束公式推导会变得很繁琐。5.3 和解析法结果的对比为了验证仿真代码的正确性我把同一网络用FMEA解析法重新算了一遍。对于这种规模很小的辐射状网络解析法可以直接枚举每个元件的故障按区域划分停电范围得到精确值。两种方法的结果对比如下SAIFI仿真值0.285解析值0.283误差约0.7%。SAIDI仿真值1.43解析值1.41误差约1.4%。误差来源主要是仿真抽样波动收敛后两者应趋于一致。这个对比的意义在于它证明代码里的故障影响分析逻辑没有系统性错误随机抽样流程也没有偏差。对于刚接触蒙特卡洛的同学我强烈建议先在小网络上用解析法交叉验证一次自己的代码确认无误后再放到大网络里跑否则出了错很难定位是指标统计问题还是抽样逻辑问题。6. 工程中踩过的坑和参数选择的经验6.1 随机数设置与结果复现蒙特卡洛仿真必须固定随机数种子否则每次跑完结果都不一样排查问题的时候会非常痛苦。Matlab里用rng(2024)设置种子可以保证同一环境下结果完全可复现。如果你要在论文里报告仿真结果一定要写明随机数种子和仿真年数并且建议用多组种子各跑几遍取均值和标准差这样可以同时报告平均值和误差带。还有一个容易被忽略的问题Matlab的rand函数默认是梅森旋转算法这个算法生成的伪随机数序列质量很高但对并行计算不友好。如果你用Parallel Computing Toolbox做多核并行仿真每个worker会以不同的随机数流运行这时要用RandStream显式创建独立的随机数流否则并行结果可能因为随机数冲突而出现莫名的相关性。实际测试中我遇到过并行跑出来的指标序列明显比串行更平滑就是因为随机数流没有正确隔离。6.2 开关动作时间和修复时间别混用这是很多刚上手的人最容易犯的错误也是最隐蔽的问题故障影响分析里对不同的负荷区域要用不同的停电时间不能一律用修复时间。我见过一个同学的代码故障影响分析函数写得比较粗糙所有下游负荷一律按修复时间计算停电时长完全没考虑分段开关隔离和联络开关转供。结果算出来的SAIDI比实际值高出两三倍整个可靠性评估失去了意义。正确做法是区分三类停电时间上游负荷的停电时间是“故障定位隔离”时间0.5h转供成功负荷的停电时间是“故障定位隔离合联络开关”时间约1.5h只有无法转供的故障区段负荷才用修复时间5h。这里还要注意一个顺序问题实际抢修中隔离和转供通常先于修复完成。所以上游和转供负荷在故障元件还没修好时就恢复供电了。如果把修复时间当作所有负荷的停电时间就完全抹掉了分段开关和联络开关的价值这样的评估没法指导规划。6.3 收敛判据不能只看仿真年数仿真年数不等于收敛。常规的5000年仿真可能在某些指标上已经收敛但如果你关心的ENS或者极端事件的方差很大可能5万年还在缓慢漂移。我建议每次都监控指标的方差系数β并且把指标随仿真年数的变化曲线画出来看它是否进入随机震荡区间。这里有个经验技巧如果仿真时间宝贵可以优先优化样本效率。一种常见的做法是“控制变量法”就是把负荷点和元件的抽样结果存下来用公共随机数比较不同方案这样方案间的差异估计会更灵敏而不是被公共随机噪声淹没。我在对比“有无联络开关”两种方案时就用同一组随机故障序列驱动两种方案结果差异的可信度显著提升。这一点在学术研究中非常实用因为审稿人最关心的往往是“差异是否显著”公共随机数法能让差异更干净地暴露出来。6.4 进阶如何处理重叠故障和元件状态耦合我前面提到了简化假设元件修复期间其他元件不发生新故障。严格来说这个假设是有误差的尤其当系统元件数量多、总故障率高的时候。要精确处理可以在每个元件进入修复状态后检查其他未故障元件的最小TTF是否小于当前修复剩余时间。如果是则说明在修复期间发生了新的故障需要按时间先后分解事件。实现方式是维护一个“修复完成时刻”列表每次推进到最早事件时不仅要考虑新故障还要考虑修复完成。这个事件驱动逻辑比简单找min(ttf)复杂一些但原理清楚。我在项目里试过严格版和简化版的对比在一个60节点网络里严格版的SAIFI和简化版相差不超过2%因为大多数元件的MTTR只有几小时而MTTF是几十天甚至几年修复期间发生新故障的概率确实很低。如果不需要抠到极致精度简化版完全够用但如果系统中存在大量修复时间特别长的电缆故障或者做的是高风险场景评估建议还是上严格版。6.5 代码性能优化与批处理技巧Matlab跑蒙特卡洛循环性能瓶颈主要在每个仿真年的while循环和函数调用开销。我做了两个优化第一把所有负荷点的上游关系预先计算成逻辑索引矩阵outage_analysis里直接按列索引省掉字符比较第二主循环里避免频繁创建临时数组所有统计量都预先分配好。如果系统特别大单次仿真要跑很久可以考虑两个方向一是把每个仿真年拆成独立batch用parfor并行跑最后合并统计量二是使用向量化的随机数生成比如一次性生成1000年的元件状态序列避免每步生成一个随机数。但向量化方案需要小心内存占用一万个元件乘以一万年的矩阵就是上亿个元素很可能把内存打爆。我实际使用中还是倾向于事件驱动循环加并行batch稳定、可控、不爆内存。