做配电网规划时最常被问到的三个问题改造这条馈线能降多少户均停电时间新装分段开关值不值接入分布式光伏后可靠性怎么变这三个问题都指向同一个底层工具——配电网可靠性评估。而在评估方法里序贯蒙特卡洛模拟法Sequential Monte Carlo, SMC是能把运行过程讲得最贴近实际的一种它把设备投运、故障、修复的过程按时间顺序一路模拟下去每一处停电都被真实地“演”出来最终统计成SAIFI、SAIDI、ENS这些指标。这篇内容围绕基于可靠性评估的序贯蒙特卡洛模拟法的Matlab实现展开从建模原理、程序架构到算例验证和工程加速完整过一遍。适合正在做配电网可靠性课题的学生也适合实际做规划评估、想知道“怎么算更靠谱”的工程师。1. 为什么配电网可靠性评估绕不开序贯蒙特卡洛1.1 解析法与非序贯抽样卡在哪配电网可靠性评估的传统主力是解析法典型代表是故障模式影响分析FMEA和最小路集法。它们的思路很直接枚举所有可能的元件故障组合再叠加计算对用户的影响。小网络没问题一旦网络稍大故障模式数量会迅速膨胀。辐射状配电网虽然故障影响能被分段开关和断路器局部化但计入隔离开关操作、联络转供、重合闸逻辑之后开关动作顺序的组合数量依然很恐怖。解析法能做到高精度建模成本也高得很快。另一种常用做法是非序贯蒙特卡洛也叫状态抽样法。它的思路是按元件故障概率随机生成系统状态统计大量状态下的期望指标。速度不错但有个硬伤它只回答“一年停电几次、停多久的期望”答不了“故障发生在一天的第几个小时”“修复花了多久”“哪些用户先恢复、哪些用户等最后”这类时序问题。分布式电源、储能、需求响应这些新兴元素偏偏又和时间强耦合非序贯法很难把它们的效益算清楚。1.2 序贯法能多出哪些关键信息序贯蒙特卡洛的核心区别在于它不抽“系统状态”而是抽“每个元件在一段持续时间内的状态变化轨迹”。每个元件会生成一串“正常→故障→修复→正常”的事件序列再把所有元件的事件序列按时间轴归并推进模拟出真实运行过程。只要随机数够好、仿真年份够长统计出来的指标就和真实系统长期运行的期望非常吻合。好处是信息完整。同一场故障中不同负荷点的停电时长可以不同能转供的停1小时不能转供的停4小时序贯模拟天然能区分出来。而且它可以直接输出每个负荷点停电次数和停电时长的分布而不只是期望值。做规划比选时光看SAIDI均值还不够有时还得看“停电超过8小时的用户数占比”“最大停电时长”这类尾部指标这些只有序贯法能给。1.3 什么项目真正需要它我的经验是单纯算几个系统级指标解析法就够了但如果要做下面几类事情序贯法基本是唯一选项。方案比选比较加装分段开关、增加联络、更换导线三种方案的可靠性收益需要稳定的指标体系支撑结论。含分布式电源/储能的评估DG的孤岛运行、储能放电策略都依赖故障时刻和持续时间非序贯法无从下手。用户级指标分析评估特定工业用户或重要负荷的停电风险不能只看全系统均值。蓄意研究故障恢复顺序比如考核调度员的转供策略、不同开关配置下的恢复时间差异。下面用一个常用对比表做个总结方法需要的主要输入输出能力时序信息适用规模主要痛点解析法/FMEA网络结构、故障率、修复时间期望指标部分依赖人工枚举小规模辐射网组合爆炸、建模费时非序贯MCS元件故障概率期望指标无大规模系统丢失时序因果序贯MCS状态持续时间分布、修复、转供策略期望分布完整中大规模系统计算量大、结果波动2. 序贯蒙特卡洛的核心建模逻辑与可靠性指标2.1 元件两状态时序模型与时序抽样原理序贯模拟的第一步是给每个元件建立一个时序状态模型。配电网设备常用二状态马尔可夫模型正常运行UP和故障停运DOWN。两个状态之间有一个转移率λ从正常运行进入故障和一个修复率μ从故障回到正常λ的单位通常是次/年μ的单位换算后是1/MTTR。这两个状态持续时间怎么抽工程上最常见的是指数分布抽样正常持续时间t_run -MTTF × ln(U1)修复持续时间t_repair -MTTR × ln(U2)其中U1、U2是[0,1]区间均匀分布的随机数MTTF 8760 / λ单位是小时。指数分布之所以常用是因为常数故障率假设下元件“无记忆性”抽样公式简单也符合大多数配电网设备在稳定运行期的失效特征。如果设备处于老化期想用威布尔分布也可以但要注意马尔可夫性失效后续状态转移概率的建模会麻烦不少。2.2 故障隔离、转供与负荷恢复的时序处理配电网是辐射状结构故障处理有一套标准时序保护动作、隔离、恢复。断路器/重合闸先断开故障下游全部失电随后通过隔离开关或分段开关隔离故障点如果存在联络开关且容量允许非故障失电区段就能转供停电时间等于开关切换时间如果联络线没有容量或者说转供条件不成立那就只能等故障修复。这段时序逻辑是所有可靠性建模里最容易被写错的地方。序贯模拟中必须严格区分“短时转供停电”和“等待修复停电”。同样是停电一次一个用户停1小时、另一个用户停5小时对SAIDI的贡献完全不同。正确的做法是故障发生后按负荷点在网络中的位置一一定义它的恢复路径和持续时间。对于重要用户可能还涉及备用电源自动投入时间这些都要在模拟前做好配置。2.3 指标定义从SAIFI到ENS的计算口径可靠性指标定义并不复杂但统计口径容易混这里统一列一下指标公式单位含义SAIFI总用户停电次数 / 总用户数次/户·年系统平均停电频率SAIDI总用户停电时间 / 总用户数小时/户·年系统平均停电持续时间CAIDI总用户停电时间 / 总停电次数小时/次每次停电的平均持续时间ASAI(用户总需求小时 - 停电小时) / 用户总需求小时无量纲供电可用率ENS各负荷点平均负荷 × 年停电时间 之和kWh/年年缺供电量注意CAIDI有两条算法路径一个是SAIDI/SAIFI一个是总停电小时数/总停电次数。当某些年没有停电事件时SAIFI为零第一条路径会出现0/0仿真实现里最好用第二种口径更稳。2.4 收敛判据与仿真年限选择序贯蒙特卡洛的结果是随机样本的统计量必须判断收敛。业界常用的判据是方差系数β定义为β σ / (μ × √N)其中σ是样本标准差μ是样本均值N是仿真年数。一般要求β小于5%即相对误差控制在5%以内。实际工程中对SAIDI、SAIFI分别监控β取最差的那个作为停止条件。仿真年份该怎么选我的习惯是至少10000年起步。辐射状配电网事件驱动模拟5万年也就几十秒到几分钟量级Matlab下时间成本完全可接受。太短了指标波动大结论在工程会上站不住脚。还要注意仿真初期的“冷启动”偏差从0时刻开始模拟第一年故障次数会比稳态偏少所以指标最好从第几年之后才开始累计或者干脆把总年限设得足够长来稀释前几百年的偏差。3. Matlab实现的程序架构与关键代码3.1 输入数据组织支路参数与拓扑描述程序第一步是组织数据。我建议用struct数组存放支路和负荷信息字段直接对应模型参数这样代码可读性和扩展性都比较好。支路数组至少包含以下字段% branch(k) 支路k的基本信息 % .from : 起点节点 % .to : 终点节点 % .lambda : 故障率次/年 % .r : 平均修复时间小时 % .ispc : 是否为分段开关 % .hasTie : 下游是否可通过联络转供负荷数组则记录节点位置、用户数、平均负荷% loadBus(m) 负荷点m的信息 % .node : 所在节点 % .numUser : 用户数 % .Pavg : 平均负荷kW如果要做时序负荷曲线可以在loadBus里再加一个8760维的Pcurve字段ENS统计时按故障时刻对应的小时取负荷值会比平均负荷精确不少。工程上很多项目只用平均负荷节约计算量但评估结果会略微低估高峰时段故障的损失。3.2 事件驱动主循环时序模拟的核心骨架仿真主循环用事件驱动方式而不是固定步长扫描。原因很简单配电网元件数量几千个仿真5万年如果按小时扫总共要扫4亿多个时间步而事件驱动只处理故障事件循环次数约等于“元件数 × 故障率 × 仿真年数”差几个数量级。核心代码骨架如下simYears 50000; % 仿真年限 totalHour simYears * 8760; nBranch length(branch); % 生成每个元件第一次故障的时间指数抽样 nextFault zeros(nBranch, 1); for k 1:nBranch nextFault(k) sampleExp(8760 / branch(k).lambda); end % 统计数组 lpNum length(loadBus); accIntNum zeros(lpNum, 1); % 累计停电次数 accIntDur zeros(lpNum, 1); % 累计停电时长 accENS zeros(lpNum, 1); % 累计缺供电量 totalUser sum([loadBus.numUser]); clockHour 0; while clockHour totalHour % 找到最早发生故障的元件 [faultTime, idx] min(nextFault); if faultTime totalHour break; end clockHour faultTime; % 查预计算的影响矩阵统计每个负荷点受影响情况 for lp 1:lpNum cat impactCat(lp, idx); % 0不受影响, 1转供短时, 2等修复 if cat 0 continue; end if cat 1 dur tieHour; % 联络转供操作时间 else dur branch(idx).r switchHour; % 修复时间开关隔离时间 end accIntNum(lp) accIntNum(lp) 1; accIntDur(lp) accIntDur(lp) dur; accENS(lp) accENS(lp) loadBus(lp).Pavg * dur; end % 修复完成安排该元件下一次故障修复时长 新的正常工作时间 nextFault(idx) clockHour branch(idx).r sampleExp(8760 / branch(idx).lambda); end其中sampleExp是自定义的指数分布抽样函数避免依赖统计工具箱function T sampleExp(meanVal) T -meanVal * log(rand()); end这段代码的统计累加用的是“总累计量”适合最后算总体指标。如果要做逐年指标序列和β收敛判断可以在循环里记录每场故障发生时所处的年份按年份分组统计即可。3.3 故障影响分析模块判断停电范围与恢复时间影响矩阵是序贯模拟能否跑快的决定性因素。教学版本可以在每场故障里实时搜索网络从故障支路向上追溯最近断路器向下枚举受影响负荷点。这个思路直观但每步都做拓扑搜索太慢而且要写大量网络遍历代码。工程做法是仿真前生成一张二维影响矩阵impactCat行是负荷点列是支路矩阵值是该支路故障时该负荷点的停电类型。生成逻辑可以概括为% 对每条支路 b % 1) 向上游找到最近的断路器/分段开关位置 % 2) 确定支路 b 下游的所有负荷点集合 % 3) 对每个下游负荷点 % 若存在联络通道且容量允许 - impactCat(lp, b) 1 % 否则 - impactCat(lp, b) 2 % 4) 不在下游集合内的负荷点 - impactCat(lp, b) 0这里最关键的是“下游”关系判定辐射状网络里负荷点是否属于某条支路的下游由根节点到负荷点的路径唯一决定。在Matlab里用有向图从根节点做一次遍历就能得到所有父子关系再结合支路位置判断任意负荷点与支路的上下游关系。这个矩阵只需要建一次仿真中完全变成查表速度能提升几十倍。判断转供能力时还要把联络线容量、对侧馈线负载率一起建模否则容易高估转供效果。3.4 指标统计与结果输出仿真结束后汇总成可靠性指标。用累计量除以总用户数和仿真年数SAIFI sum(accIntNum .* [loadBus.numUser]) / totalUser / simYears; SAIDI sum(accIntDur .* [loadBus.numUser]) / totalUser / simYears; ENS sum(accENS) / simYears; CAIDI sum(accIntDur .* [loadBus.numUser]) / sum(accIntNum .* [loadBus.numUser]);这里要注意accIntDur是和单个负荷点累加的时长乘以该负荷点用户数再求和才等于“用户停电时户数”。很多初学者在这里漏乘用户数导致SAIDI数值差一个数量级。负荷点级指标、系统指标一起输出方便后面定位薄弱环节。4. 算例验证双馈线联络系统算例解读4.1 算例结构与参数设定为了快速验证程序逻辑我搭了一个简化双馈线系统两条馈线L1和L2L1带三个负荷点L2带两个负荷点两馈线末端通过一个联络开关相连。每条支路的故障率统一取0.10次/年修复时间统一取4小时转供操作时间1小时开关隔离时间0.5小时。每个负荷点用户数取50总用户数250。负荷点所在馈线平均负荷(kW)用户数转供通道LP1L1首段下游20050无LP2L1中段下游25050无LP3L1末端30050经联络开关至L2LP4L2首段下游18050无LP5L2末端22050经联络开关至L14.2 收敛特性与指标输出这种小规模算例可以先把手算期望值列出来做程序校验SAIFI 0.18次/户·年SAIDI 0.46小时/户·年ENS约530 kWh/年。跑完50000年仿真结果应该在这个值附近摆动一般偏差在±5%以内。如果偏差明显偏大优先查影响矩阵是否正确、修复时间有没有漏加开关操作时间。仿真收敛方面我监控SAIDI的方差系数β大约在3万多年后降到5%以下。对一个大中型配电网这个收敛速度算正常。指标波动主要来自故障率较低的支路它们在整个仿真期间可能只发生几次故障抽样方差大想压下来只能靠延长仿真总年限或做方差削减。4.3 结果怎么用方案比选与薄弱点定位算例最有价值的应用是把联络开关从系统中摘除重新仿真对比“有联络”和“无联络”两种方案方案SAIFI(次/户·年)SAIDI(h/户·年)ENS(kWh/年)有联络转供0.180.46530无联络摘除联络0.180.81999SAIFI几乎不变因为故障次数和用户范围没变但SAIDI接近翻倍ENS接近翻倍。这说明联络开关的价值主要体现在缩短停电时长而不是减少停电次数。这类结论光看SAIFI指标根本体现不出来正是序贯法的用武之地。再看负荷点级指标LP3和LP5的每次停电都能通过联络转供压到1小时而LP1、LP2、LP4的停电时长都接近修复时间4.5小时。如果你的辖区里有类似LP1这样无转供通道、故障影响面又大的末端用户改造方向就知道该往哪使劲了要么加联络、要么加密分段开关缩小故障范围。这个分析逻辑直接复用到真实10kV馈线就是日常工作。5. 工程实践中的坑位清单与加速思路5.1 最容易翻车的五个细节第一坑是冷启动偏差。从0时刻开始模拟时第一年内设备刚投运故障次数小于稳态期望导致前几年指标偏低。如果总仿真年限只有几百年这个偏差会明显拉低最终结果。应对方法很简单仿真年限设长或者丢弃初始若干年再统计。第二坑是随机数种子问题。Matlab的rand每次运行序列不同导致同参数两次数值差异。做方案比选时这种随机噪声可能掩盖真实差别。我的做法是每个方案用同一个随机数种子让所有候选方案面对同一套故障场景比选结论会干净很多。第三坑是修复时间的分布选择。指数分布抽样会让修复时长方差很大个别长修复事件会把SAIDI尾部拉高。工程上很多项目直接取固定MTTR值结果解释起来更容易代码也更简单。如果研究重点是长修复时间对高可靠性用户的冲击再用指数分布不迟。第四坑是联络容量约束。很多教程代码默认联络开关可以无限转供实际配电网中联络线容量有限对侧馈线也不一定有余量。忽略这个约束SAIDI会被明显低估。建模时至少给每条联络线配一个最大可转供容量故障时比对负荷大小。第五坑是共因故障。台风、外力破坏会让一批元件同时故障简单二状态独立故障模型完全反映不了这种场景。做韧性评估或者极端天气下的可靠性分析需要另外建模共因失效组不能直接用常规序贯模型。5.2 仿真太慢怎么办最小路集、预计算矩阵与并行序贯蒙特卡洛最慢的部分通常是影响分析和随机数生成。影响分析走“预计算影响矩阵”的路线速度提升最明显矩阵一旦建好仿真主循环就是几个矩阵索引操作。其次是随机数指数抽样本质是一次log和一次rand5万仿真年也就几十万次抽样Matlab跑起来没压力。真正的大系统瓶颈在负荷点与支路的规模乘积上这时可以按“同一支路组的负荷点合并统计”做粗粒度简化损失一点用户级精度换来指标计算量指数下降。并行计算方面Matlab的parfor可以用但要注意随机数流每个worker独立生成随机数不控制好RandStream的话并行结果不可复现。我的做法是给每个worker分配独立的随机子流方案比选时仍用同一主种子派生。另外还有一个工程技巧先跑一遍较短年限的预仿真看哪些支路的故障对系统指标贡献最大再对这些支路做重要抽样或者分层抽样。用控制变量法以解析结果作为控制量也能进一步缩小方差。这些方法实现成本略高但对动辄上万元件的大馈线集群收益很可观。5.3 结果如何汇报才可信置信区间与敏感性报告里只给一个SAIDI数值甲方往往会问“准不准”。更专业的做法是给置信区间比如SAIDI 0.46 ± 0.02小时/户·年置信度95%。置信区间的宽度由β和仿真年数决定工程上把β压到5%以下区间宽度就是可接受的。建议在仿真输出里直接把均值、标准差、β、95%置信区间一并打出来。敏感性分析也是重要一环。故障率和修复时间取自历史统计本身就带误差。常规做法是把关键参数上下浮动20%~30%重新仿真后看SAIDI变化范围。如果某个参数对结果影响特别大说明评估结论对这个参数敏感报告中要特别标注“这里的统计数据来自XX厂家的运行统计不同厂家可能有差异”。这些细节做好了评估报告才算完整。我自己的习惯是方案比选时固定随机种子让不同方案共享同一套故障场景SAIDI差0.05小时也能看出趋势平时只报基准指标再用长仿真年限加β判据确保数字站得住脚。序贯蒙特卡洛模拟法本身不难难的是把配电网的操作时序、转供约束、负荷变化这些现场逻辑建模清楚。把这些细节处理到位Matlab代码就只是表达工具了评估结果才能真正支撑工程决策。