做过配电网可靠性分析的朋友应该都遇到过这个问题常规可靠性评估里故障率是一个长期平均值算出来的结果是“年化期望”但台风一来这个结果根本没法用来做应急准备。台风天气下配电网故障建模与场景生成解决的就是这个矛盾——把台风过程的动态影响映射成设备故障概率再用场景生成的手段输出一批“哪几条线可能一起跳、影响多少负荷”的具体情形。这套方法在某个33节点配电网算例上落地用Matlab实现了从风场构建、故障率时变计算、蒙特卡洛抽样到应急响应指标评估的全流程尤其重点处理了“如何通过故障特征引入配电网应急响应内容”这一环。这篇文章我把项目的整体设计、核心参数、代码结构和踩坑过程都拆开讲适合电力系统方向的研究生、配电网规划工程师以及做韧性评估的研发人员参考。1. 项目到底在做什么从台风到配电网故障场景1.1 台风凭什么能让配电网“成片倒”先说说台风对配电网的破坏机理这是后面所有建模的物理基础。台风导致配电网大量故障的路径其实很直接强风作用下杆塔可能倾斜甚至倒塔导线会发生大幅舞动甚至断线周边树木倒伏压线暴雨伴随的低气压环境会让空气湿度接近饱和绝缘子表面闪络概率显著上升分支线上的跌落保险、避雷器也更容易损坏。配电网的电压等级低、绝缘裕度小很多架空线路的通道净空不足加上多数馈线是单电源辐射状结构一旦某一段故障下游负荷全部跟着失电。所以台风场景下的故障不是孤立事件而是成组出现、沿着台风路径方向和强度衰减方向高度聚集。我见过一次沿海区县的台风过程最大风速达到过25 m/s一条10 kV馈线上同时跳闸的支路有7条平时一年的故障次数都没这一晚多。这就是常规可靠性评估用单一平均故障率根本覆盖不了的现象。台风过程是一个短时、高破坏强度的动态过程风速在空间上随时间推移故障概率也在时空上快速变化。因此建模时不能把故障率当成常数必须引入“时变故障率”和“空间关联性”。1.2 33节点配电网为什么大家爱拿它当试验田这个项目选择某个33节点标准测试系统不是随手选的而是因为它太适合做这类故障场景分析。这个33节点配电网的结构特点很典型1号节点作为变电站出口平衡节点32条支路构成辐射状馈线基准电压12.66 kV总负荷接近3700 kW和2300 kvar左右。它规模不大但足够反映线路分段、分支负荷、末端电压偏低这些配电网的真实问题。与动辄几百上千节点的实际馈线相比33节点系统数据公开透明潮流计算方便故障后负荷损失也能手工核对。我经常跟人说做算法验证最怕发现问题之后不知道是你代码写错了还是数据错了33节点系统的好处就在这里是“标准答案”代码跑完能自查。台风故障建模要落地最终要回答“哪些负荷没电”这类问题而33节点系统每条支路挂接的负荷都明明白白。把台风场景叠加到这个网络上故障集合对应明确的失电区域后续的应急响应分析才可复现、可对比。当然方法框架并不局限于这个算例把节点坐标、支路参数和负荷数据替换成更大的系统整体流程完全照用。1.3 故障场景生成与“应急响应”的关系很多人拿到这个项目标题时第一反应是“场景生成”和“应急响应”是两件独立的事。实际上二者是由“故障特征”串联起来的。故障场景生成解决的是“台风来了之后系统可能处于什么样的故障状态”。通过蒙特卡洛抽样我们可以得到成千上万个可能的故障集合比如“线路2和线路9同时跳闸”“线路15单独跳闸”等。每个故障集合对应一个场景带有概率权重。但应急响应部门真正需要的不是几十万个抽样场景而是“这次可能有多少用户失电、失电负荷有多大、哪些重要负荷受影响、大概多久能恢复”。这些信息必须从故障场景中提炼成故障特征才能输入到应急决策模型里。这里有一个很关键的思路我在项目里专门做了设计每个故障场景不但记录故障线路编号还同步计算出该场景下的失电负荷、失电用户数、关键负荷占比、可转供容量等特征字段。应急响应模块拿到这些字段后再做抢修优先级排序和恢复时间评估。也就是说故障特征是从场景生成到应急响应的“接口协议”接口定义得好不好直接决定整套流程的可用性。2. 故障建模的核心逻辑和参数设计2.1 风速场模型与坐标映射台风场景生成的起点是算出每个节点在特定时刻所处的风速。为了简化处理项目采用圆对称风场模型并在此基础上修正移动速度的影响。圆对称模型的基本思想是台风中心有一个最大风速半径在半径之内风速随半径增大而增大在半径之外则逐渐衰减。常用的一种形式是当距离 r 小于等于最大风速半径 Rmax 时风速 V(r) Vmax * (r / Rmax)当距离 r 大于最大风速半径 Rmax 时风速 V(r) Vmax * (Rmax / r)^x其中 x 通常取0.4~0.6。这个模型虽然简单但胜在参数少、计算方便适合作为第一版实现。如果你拿到了更精细的气象再分析数据可以换成更复杂的经验风场模型但在Matlab代码框架上底层函数只需要修改输入输出就行。坐标映射是第一个容易出问题的地方。项目采用自定义平面坐标把33个节点放在一张 xy 平面图上另外定义台风中心路径为一段随时间变化的点序列。在每一时刻依次计算每个节点和台风中心的距离 r代入上面的公式就得到该节点的瞬时风速。如果你用的是经纬度坐标必须先做投影转换否则距离算出来的单位是“度”风速公式完全没法用。我习惯在程序开头统一转成以km为单位的平面坐标这样输入的长度、风速、时间维度都对得上。台风中心路径不是静止的。实际仿真中把整个台风过境过程按15分钟到30分钟的时间步长离散每个步长内更新一次台风中心位置然后重新计算全部节点的风速。这里建议时间步长不要超过30分钟否则会漏掉阵风峰值后面计算故障率时误差会被指数放大。2.2 线路时变故障率函数怎么选风速和故障率之间的关系是整个建模的核心。我把它分成三种常用形式第一种是分段线性模型认为风速低于某个阈值时故障率基本不变超过阈值后故障率随风速线性上升到一定值后饱和。这种模型直观易标定适合快速验证。第二种是指数模型形式为 λ(v) λ0 * exp(α * max(v - v0, 0))。其中 λ0 是正常天气下的故障率v0 是开始明显影响的风速阈值α 是拟合系数。指数模型能反映强风下故障率急剧上升的特点台风场景下更贴近实际项目里我采用了这一种。第三种是多项式模型一般用于拟合历史故障统计数据如果有足够多的历史数据支撑会更准确但多数研究场景下数据不足反而不如前两种实用。在Matlab里实现时可以把线路的时变故障率做成一个矩阵行是线路编号列是时间步。核心代码风格可以这样写% 线路故障率时变计算 lambda0 0.1; % 常态故障率单位: 次/(km·年) alpha 0.15; % 指数拟合系数 v0 10; % 风速阈值单位: m/s % windSpeedMat 是 [n_branch x n_time] 的风速矩阵 deltaV max(windSpeedMat - v0, 0); lambdaMat lambda0 / 8760 * exp(alpha * deltaV); % 换算为 次/(km·小时)这里特别提醒一下单位换算。很多人直接用“次/(km·年)”的故障率去乘以步长结果小得离谱。正确做法是先除以8760变成小时级的故障率再与时间步长相乘。否则概率算出来几乎为零后面蒙特卡洛抽样全废。另外不同设备元件的故障率参数是不一样的。架空导线的故障率、杆塔的故障率、配电变压器的故障率在台风下的响应系数也不同。项目里至少要把线路分成三类参数主干线、分支线和配变支路分别设置 λ0 和 α这样算出来的场景才不至于把所有线路都当成同一种脆弱性。2.3 故障场景生成的蒙特卡洛流程故障场景生成的核心是蒙特卡洛随机抽样基本逻辑不复杂。在每一个时间步长 dt 内一条故障率为 λ 的线路发生故障的概率为p 1 - exp(-λ * dt)当 λ 很小时p 近似等于 λ * dt但台风场景下 λ 可能比常态高出几十倍直接用近似式会低估概率所以还是保留指数形式更稳。抽样流程是对每条线路生成一个[0,1]区间均匀分布的随机数如果小于 p那就认为该线路在这个时间步发生故障。重复 N 次得到 N 个场景。每个场景的本质是一个“故障线路编号集合”例如场景1{线路2, 线路9}场景2{线路15}场景3{线路4, 线路16, 线路27}抽样代码可以写成N_scen 5000; dt 0.25; % 单位: 小时 (15分钟) rng(2024); % 固定随机种子保证可复现 scenCell cell(1, N_scen); for m 1:N_scen faultLines []; for k 1:n_branch lambda_k mean(lambdaMat(k, :)); % 取整个台风过程的平均故障率做简化 pf_k 1 - exp(-lambda_k * dt); if rand pf_k faultLines [faultLines; k]; end end scenCell{m} faultLines; end严格一点的实现会按时间步推进在每一时刻判断故障发生再记录故障时刻和恢复时刻。但初期跑框架时可以先按全过程平均故障率抽样逻辑更清楚。抽样出来的 N 个场景不能全部拿去做应急分析否则计算量太大而且很多场景非常相似。通常用同步回代削减法或K-means聚类把相似场景合并得到 S 个典型场景及对应的概率权重。场景削减不是单纯“减少数量”它的目标是让削减后的典型场景集合在期望缺供电量、失电负荷分布等关键指标上与原场景集合的误差控制在可接受范围内。2.4 如何把故障特征引到应急响应内容里这一节是标题里特别点出来的内容我的处理方式是直接在场景数据结构里嵌入故障特征字段。首先每个故障场景生成后立即计算以下特征故障线路集合失电负荷总量MW需要结合配电网拓扑和潮流计算结果失电用户数估算值关键负荷标记比如线路是否带医院、水厂、通信基站等重要用户可转供容量即相邻联络开关是否存在、转供后是否过载。有了这些特征应急响应模块就能做三件事。第一抢修优先级排序。按“重要用户、失电负荷、影响用户数”加权排序先处理影响最严重的故障线路。第二恢复时间估算。给每条故障线路配一个修复时间分布按分布抽样得到恢复时间再考虑抢修队数量有限这一约束输出“失电负荷随时间变化曲线”。第三应急指标统计。包括最大同时失电负荷、系统平均停电持续时间SAIDI、缺供电量ENS等。这些指标就是台风场景下应急响应内容的具体落点。这里有个实操心得故障特征计算要和潮流分析联动。有些故障线路跳闸后下游负荷并不会全部失电因为可以通过联络开关从相邻馈线转供。如果忽略可转供容量算出来的失电负荷会严重偏高应急响应方案也会偏保守。所以在场景结构体中除了故障线路集合我还保存了“实际失电节点集合”和“失电负荷总量”这两个字段是在潮流计算之后写进去的。3. Matlab实现路径与关键代码拆解3.1 整体程序架构一个完整的Matlab实现程序结构大致分为七块主脚本 main_typhoon_fault_analysis.m负责串联整个流程get_typhoon_parameters.m设置台风路径、最大风速、最大风速半径等参数wind_field_at_nodes.m计算所有节点的风速时变矩阵failure_rate_with_time.m计算线路时变故障率矩阵scenario_sampling.m蒙特卡洛场景抽样scenario_reduction.m场景削减calculate_emergency_metrics.m计算应急响应指标并输出结果。数据组织上我习惯用矩阵和结构体混搭% 某33节点配电网部分支路数据 % bus_i bus_j r(pu) x(pu) length(km) branch [ 1 2 0.0922 0.0470 0.30; 2 3 0.4930 0.2511 0.50; 3 4 0.3660 0.1864 0.40; ... ]; % 节点坐标矩阵第二列是x坐标第三列是y坐标单位km nodeCoord [ 1 0.0000 0.0000; 2 0.3000 0.1000; ... ]; % 台风路径结构体 typhoon.path [0.5 0.2; 0.7 0.3; ...]; % 每个时间步的台风中心坐标 typhoon.Vmax 40; % 最大风速 m/s typhoon.Rmax 30; % 最大风速半径 km结构体存台风参数矩阵存节点和支路数据cell数组存场景这样可以很自然地把不同数据类型组织起来。3.2 故障概率采样与场景削减的代码逻辑场景抽样代码前面已经给过一版简化实现。在实际项目中我会把“故障时刻”也记下来方便后续恢复时间计算scenCell{m} struct(... faultLines, faultLines, ... faultTimes, faultTimes, ... outageLoadMW, 0, ... affectedUsers, 0, ... criticalLoadFlag, 0);抽样完成后做潮流计算用实际失电节点集合更新outageLoadMW和affectedUsers。场景削减我经常用同步回代削减法思路是反复合并距离最小的两个场景并把被合并场景的概率权重累加到保留场景上。在Matlab里可以直接调用聚类函数做初步压缩但严谨研究建议自己写同步回代削减。特征之间的距离可以用Jaccard距离即两个场景故障集合的交并比。如果两个集合的故障线路几乎一致它们距离就小可以合并。削减后的典型场景数量控制在10到20个之间。数量太少会丢掉小概率高风险场景数量太多又会拖慢应急决策模块。我一般先用削减前后ENS误差做校验误差不超过5%就认为削减结果可以接受。3.3 结果可视化和指标计算台风场景生成类项目结果可视化很重要。至少需要画三张图。第一张是台风路径与33节点位置的叠加图figure; hold on; plot(typhoon.path(:,1), typhoon.path(:,2), r-, LineWidth, 2); scatter(nodeCoord(:,1), nodeCoord(:,2), 40, filled); xlabel(x (km)); ylabel(y (km)); title(某台风路径与33节点配电网位置示意);第二张是线路故障概率的空间分布图可以用线路线宽或颜色深浅来表示故障概率高低。这样能直观看出台风路径附近的线路故障风险更高。第三张是典型场景的失电负荷柱状图每个柱子对应一个典型场景柱高表示该场景下的失电负荷旁边标注场景概率。指标计算方面重点算缺供电量ENSENS Σ_s p_s * Σ_i P_loss(i) * T_recovery(i)其中 p_s 是场景概率P_loss(i) 是第 i 个失电负荷节点的有功T_recovery(i) 是恢复时间。应急响应模块输出最大失电负荷、平均恢复时间、关键用户失电时长等指标这些都可以汇总到表格里输出。3.4 参数标定和仿真设置的心得参数标定是整个项目最花时间的部分也是最影响结果可信度的部分。我的建议是分三层台风参数、故障率参数、场景削减参数。台风参数方面最大风速和最大风速半径要基于历史台风数据标定最好找接近当地典型强度的台风过程做参考。如果没有实测数据至少要在仿真里做敏感性分析看结果对参数的敏感程度。故障率参数方面λ0 可以根据历史跳闸率统计α 比较难标定通常在两倍量级范围内调参并结合敏感性分析说明结果稳健性。仿真设置上几个经验值可以参考初始场景数取5000左右统计稳定性够用削减后保留场景数10到20个时间步长15到30分钟全过程时长一般覆盖台风过境的24到48小时每组仿真前用不同的随机种子运行3到5次检查指标波动幅度。4. 实操中踩过的坑与排查技巧4.1 风速模型选错导致故障率失真我最早用最简单的圆对称公式风速偏大故障率直接高了一个数量级生成的场景几乎每条线路都在跳一看就不合理。后来发现问题是忽略了最大风速半径之外的快速衰减以及台风移速叠加。排查方法很简单把风场计算得到节点风速画成时间曲线和气象站实测数据放在一起对比。实测可用的就对比不可用的就对比文献里典型的台风风速曲线形态。只要风速场形状不对后面的故障率必然失真这个环节值得第一个检查。4.2 场景削减时场景数怎么定场景削减的典型错误是贪图计算速度把场景砍得太狠。某个测试里我把5000个场景直接削减到5个结果ENS和未削减前差了将近20%尾部风险场景全部被抹掉。反过来削减到50个场景应急响应模块每轮计算要跑很久。实操建议先保留10个、15个、20个、30个场景做对比实验计算削减前后ENS误差曲线找到误差小于5%的最少场景数。大多数情况下15到20个场景就足够稳定但这个数值要在你自己的算例上验证不能直接抄。4.3 应急响应指标和故障场景的衔接问题最容易掉进去的坑是场景抽样后只记录故障线路编号没有及时计算失电负荷和用户数导致应急响应模块无从下手。解决这个问题一定要在场景数据结构里预留字段scenCell{m}.outageLoadMW outageLoadMW; scenCell{m}.affectedUsers affectedUsers; scenCell{m}.criticalLoadFlag criticalLoadFlag;另一个衔接问题是忽略转供。33节点配电网系统中一些联络开关可以让部分负荷从其他馈线恢复供电。应急响应里如果不考虑转供能力恢复时间会被高估。处理方法是把故障场景传给潮流计算模块用实际潮流结果决定哪些节点真正失电再更新故障特征。4.4 Matlab内存与随机数管理场景数一多Matlab经常报内存不足。原因往往是某个中间变量不小心把 N_scen 和 n_branch 都用满矩阵存瞬间爆炸。解决办法是场景用cell数组存稀疏结构风速场不要存全过程三维矩阵只保存节点风速矩阵33×时间步即可。这个规模的内存占用非常小。随机数管理同样是个细节。我最开始为了“复现”在每次抽样前都调用一次rng(seed)结果每个场景都抽到了相同的随机数序列生成的场景几乎全部一样。正确做法是在主脚本开头设置一次随机种子循环里直接用rand这样才能得到独立抽样。如果要做并行计算建议给每个并行worker单独配置随机种子流避免重复场景。最后再分享一点个人体会这个项目做完之后我最大的感受是台风故障建模的难点不在于把故障率公式写得多漂亮而在于让场景生成结果真正能被应急决策用起来。故障特征字段的设计一定要前置在写蒙特卡洛抽样之前就想清楚“后面应急响应要什么数据”。另外整个流程跑通之后后续扩展其实很方便可以把修复时间分布改成更精细的抢修队调度模型也可以把静态场景生成升级为时序恢复过程分析甚至接入移动储能应急接入策略。33节点算例只是一个起点方法框架本身是可以复用的。 
长河君 · 2026/10/21 14:04:48 · 阅读 1.2k 台风场景下配电网故障建模与应急响应：从时变故障率到场景生成
做过配电网可靠性分析的朋友应该都遇到过这个问题常规可靠性评估里故障率是一个长期平均值算出来的结果是“年化期望”但台风一来这个结果根本没法用来做应急准备。台风天气下配电网故障建模与场景生成解决的就是这个矛盾——把台风过程的动态影响映射成设备故障概率再用场景生成的手段输出一批“哪几条线可能一起跳、影响多少负荷”的具体情形。这套方法在某个33节点配电网算例上落地用Matlab实现了从风场构建、故障率时变计算、蒙特卡洛抽样到应急响应指标评估的全流程尤其重点处理了“如何通过故障特征引入配电网应急响应内容”这一环。这篇文章我把项目的整体设计、核心参数、代码结构和踩坑过程都拆开讲适合电力系统方向的研究生、配电网规划工程师以及做韧性评估的研发人员参考。1. 项目到底在做什么从台风到配电网故障场景1.1 台风凭什么能让配电网“成片倒”先说说台风对配电网的破坏机理这是后面所有建模的物理基础。台风导致配电网大量故障的路径其实很直接强风作用下杆塔可能倾斜甚至倒塔导线会发生大幅舞动甚至断线周边树木倒伏压线暴雨伴随的低气压环境会让空气湿度接近饱和绝缘子表面闪络概率显著上升分支线上的跌落保险、避雷器也更容易损坏。配电网的电压等级低、绝缘裕度小很多架空线路的通道净空不足加上多数馈线是单电源辐射状结构一旦某一段故障下游负荷全部跟着失电。所以台风场景下的故障不是孤立事件而是成组出现、沿着台风路径方向和强度衰减方向高度聚集。我见过一次沿海区县的台风过程最大风速达到过25 m/s一条10 kV馈线上同时跳闸的支路有7条平时一年的故障次数都没这一晚多。这就是常规可靠性评估用单一平均故障率根本覆盖不了的现象。台风过程是一个短时、高破坏强度的动态过程风速在空间上随时间推移故障概率也在时空上快速变化。因此建模时不能把故障率当成常数必须引入“时变故障率”和“空间关联性”。1.2 33节点配电网为什么大家爱拿它当试验田这个项目选择某个33节点标准测试系统不是随手选的而是因为它太适合做这类故障场景分析。这个33节点配电网的结构特点很典型1号节点作为变电站出口平衡节点32条支路构成辐射状馈线基准电压12.66 kV总负荷接近3700 kW和2300 kvar左右。它规模不大但足够反映线路分段、分支负荷、末端电压偏低这些配电网的真实问题。与动辄几百上千节点的实际馈线相比33节点系统数据公开透明潮流计算方便故障后负荷损失也能手工核对。我经常跟人说做算法验证最怕发现问题之后不知道是你代码写错了还是数据错了33节点系统的好处就在这里是“标准答案”代码跑完能自查。台风故障建模要落地最终要回答“哪些负荷没电”这类问题而33节点系统每条支路挂接的负荷都明明白白。把台风场景叠加到这个网络上故障集合对应明确的失电区域后续的应急响应分析才可复现、可对比。当然方法框架并不局限于这个算例把节点坐标、支路参数和负荷数据替换成更大的系统整体流程完全照用。1.3 故障场景生成与“应急响应”的关系很多人拿到这个项目标题时第一反应是“场景生成”和“应急响应”是两件独立的事。实际上二者是由“故障特征”串联起来的。故障场景生成解决的是“台风来了之后系统可能处于什么样的故障状态”。通过蒙特卡洛抽样我们可以得到成千上万个可能的故障集合比如“线路2和线路9同时跳闸”“线路15单独跳闸”等。每个故障集合对应一个场景带有概率权重。但应急响应部门真正需要的不是几十万个抽样场景而是“这次可能有多少用户失电、失电负荷有多大、哪些重要负荷受影响、大概多久能恢复”。这些信息必须从故障场景中提炼成故障特征才能输入到应急决策模型里。这里有一个很关键的思路我在项目里专门做了设计每个故障场景不但记录故障线路编号还同步计算出该场景下的失电负荷、失电用户数、关键负荷占比、可转供容量等特征字段。应急响应模块拿到这些字段后再做抢修优先级排序和恢复时间评估。也就是说故障特征是从场景生成到应急响应的“接口协议”接口定义得好不好直接决定整套流程的可用性。2. 故障建模的核心逻辑和参数设计2.1 风速场模型与坐标映射台风场景生成的起点是算出每个节点在特定时刻所处的风速。为了简化处理项目采用圆对称风场模型并在此基础上修正移动速度的影响。圆对称模型的基本思想是台风中心有一个最大风速半径在半径之内风速随半径增大而增大在半径之外则逐渐衰减。常用的一种形式是当距离 r 小于等于最大风速半径 Rmax 时风速 V(r) Vmax * (r / Rmax)当距离 r 大于最大风速半径 Rmax 时风速 V(r) Vmax * (Rmax / r)^x其中 x 通常取0.4~0.6。这个模型虽然简单但胜在参数少、计算方便适合作为第一版实现。如果你拿到了更精细的气象再分析数据可以换成更复杂的经验风场模型但在Matlab代码框架上底层函数只需要修改输入输出就行。坐标映射是第一个容易出问题的地方。项目采用自定义平面坐标把33个节点放在一张 xy 平面图上另外定义台风中心路径为一段随时间变化的点序列。在每一时刻依次计算每个节点和台风中心的距离 r代入上面的公式就得到该节点的瞬时风速。如果你用的是经纬度坐标必须先做投影转换否则距离算出来的单位是“度”风速公式完全没法用。我习惯在程序开头统一转成以km为单位的平面坐标这样输入的长度、风速、时间维度都对得上。台风中心路径不是静止的。实际仿真中把整个台风过境过程按15分钟到30分钟的时间步长离散每个步长内更新一次台风中心位置然后重新计算全部节点的风速。这里建议时间步长不要超过30分钟否则会漏掉阵风峰值后面计算故障率时误差会被指数放大。2.2 线路时变故障率函数怎么选风速和故障率之间的关系是整个建模的核心。我把它分成三种常用形式第一种是分段线性模型认为风速低于某个阈值时故障率基本不变超过阈值后故障率随风速线性上升到一定值后饱和。这种模型直观易标定适合快速验证。第二种是指数模型形式为 λ(v) λ0 * exp(α * max(v - v0, 0))。其中 λ0 是正常天气下的故障率v0 是开始明显影响的风速阈值α 是拟合系数。指数模型能反映强风下故障率急剧上升的特点台风场景下更贴近实际项目里我采用了这一种。第三种是多项式模型一般用于拟合历史故障统计数据如果有足够多的历史数据支撑会更准确但多数研究场景下数据不足反而不如前两种实用。在Matlab里实现时可以把线路的时变故障率做成一个矩阵行是线路编号列是时间步。核心代码风格可以这样写% 线路故障率时变计算 lambda0 0.1; % 常态故障率单位: 次/(km·年) alpha 0.15; % 指数拟合系数 v0 10; % 风速阈值单位: m/s % windSpeedMat 是 [n_branch x n_time] 的风速矩阵 deltaV max(windSpeedMat - v0, 0); lambdaMat lambda0 / 8760 * exp(alpha * deltaV); % 换算为 次/(km·小时)这里特别提醒一下单位换算。很多人直接用“次/(km·年)”的故障率去乘以步长结果小得离谱。正确做法是先除以8760变成小时级的故障率再与时间步长相乘。否则概率算出来几乎为零后面蒙特卡洛抽样全废。另外不同设备元件的故障率参数是不一样的。架空导线的故障率、杆塔的故障率、配电变压器的故障率在台风下的响应系数也不同。项目里至少要把线路分成三类参数主干线、分支线和配变支路分别设置 λ0 和 α这样算出来的场景才不至于把所有线路都当成同一种脆弱性。2.3 故障场景生成的蒙特卡洛流程故障场景生成的核心是蒙特卡洛随机抽样基本逻辑不复杂。在每一个时间步长 dt 内一条故障率为 λ 的线路发生故障的概率为p 1 - exp(-λ * dt)当 λ 很小时p 近似等于 λ * dt但台风场景下 λ 可能比常态高出几十倍直接用近似式会低估概率所以还是保留指数形式更稳。抽样流程是对每条线路生成一个[0,1]区间均匀分布的随机数如果小于 p那就认为该线路在这个时间步发生故障。重复 N 次得到 N 个场景。每个场景的本质是一个“故障线路编号集合”例如场景1{线路2, 线路9}场景2{线路15}场景3{线路4, 线路16, 线路27}抽样代码可以写成N_scen 5000; dt 0.25; % 单位: 小时 (15分钟) rng(2024); % 固定随机种子保证可复现 scenCell cell(1, N_scen); for m 1:N_scen faultLines []; for k 1:n_branch lambda_k mean(lambdaMat(k, :)); % 取整个台风过程的平均故障率做简化 pf_k 1 - exp(-lambda_k * dt); if rand pf_k faultLines [faultLines; k]; end end scenCell{m} faultLines; end严格一点的实现会按时间步推进在每一时刻判断故障发生再记录故障时刻和恢复时刻。但初期跑框架时可以先按全过程平均故障率抽样逻辑更清楚。抽样出来的 N 个场景不能全部拿去做应急分析否则计算量太大而且很多场景非常相似。通常用同步回代削减法或K-means聚类把相似场景合并得到 S 个典型场景及对应的概率权重。场景削减不是单纯“减少数量”它的目标是让削减后的典型场景集合在期望缺供电量、失电负荷分布等关键指标上与原场景集合的误差控制在可接受范围内。2.4 如何把故障特征引到应急响应内容里这一节是标题里特别点出来的内容我的处理方式是直接在场景数据结构里嵌入故障特征字段。首先每个故障场景生成后立即计算以下特征故障线路集合失电负荷总量MW需要结合配电网拓扑和潮流计算结果失电用户数估算值关键负荷标记比如线路是否带医院、水厂、通信基站等重要用户可转供容量即相邻联络开关是否存在、转供后是否过载。有了这些特征应急响应模块就能做三件事。第一抢修优先级排序。按“重要用户、失电负荷、影响用户数”加权排序先处理影响最严重的故障线路。第二恢复时间估算。给每条故障线路配一个修复时间分布按分布抽样得到恢复时间再考虑抢修队数量有限这一约束输出“失电负荷随时间变化曲线”。第三应急指标统计。包括最大同时失电负荷、系统平均停电持续时间SAIDI、缺供电量ENS等。这些指标就是台风场景下应急响应内容的具体落点。这里有个实操心得故障特征计算要和潮流分析联动。有些故障线路跳闸后下游负荷并不会全部失电因为可以通过联络开关从相邻馈线转供。如果忽略可转供容量算出来的失电负荷会严重偏高应急响应方案也会偏保守。所以在场景结构体中除了故障线路集合我还保存了“实际失电节点集合”和“失电负荷总量”这两个字段是在潮流计算之后写进去的。3. Matlab实现路径与关键代码拆解3.1 整体程序架构一个完整的Matlab实现程序结构大致分为七块主脚本 main_typhoon_fault_analysis.m负责串联整个流程get_typhoon_parameters.m设置台风路径、最大风速、最大风速半径等参数wind_field_at_nodes.m计算所有节点的风速时变矩阵failure_rate_with_time.m计算线路时变故障率矩阵scenario_sampling.m蒙特卡洛场景抽样scenario_reduction.m场景削减calculate_emergency_metrics.m计算应急响应指标并输出结果。数据组织上我习惯用矩阵和结构体混搭% 某33节点配电网部分支路数据 % bus_i bus_j r(pu) x(pu) length(km) branch [ 1 2 0.0922 0.0470 0.30; 2 3 0.4930 0.2511 0.50; 3 4 0.3660 0.1864 0.40; ... ]; % 节点坐标矩阵第二列是x坐标第三列是y坐标单位km nodeCoord [ 1 0.0000 0.0000; 2 0.3000 0.1000; ... ]; % 台风路径结构体 typhoon.path [0.5 0.2; 0.7 0.3; ...]; % 每个时间步的台风中心坐标 typhoon.Vmax 40; % 最大风速 m/s typhoon.Rmax 30; % 最大风速半径 km结构体存台风参数矩阵存节点和支路数据cell数组存场景这样可以很自然地把不同数据类型组织起来。3.2 故障概率采样与场景削减的代码逻辑场景抽样代码前面已经给过一版简化实现。在实际项目中我会把“故障时刻”也记下来方便后续恢复时间计算scenCell{m} struct(... faultLines, faultLines, ... faultTimes, faultTimes, ... outageLoadMW, 0, ... affectedUsers, 0, ... criticalLoadFlag, 0);抽样完成后做潮流计算用实际失电节点集合更新outageLoadMW和affectedUsers。场景削减我经常用同步回代削减法思路是反复合并距离最小的两个场景并把被合并场景的概率权重累加到保留场景上。在Matlab里可以直接调用聚类函数做初步压缩但严谨研究建议自己写同步回代削减。特征之间的距离可以用Jaccard距离即两个场景故障集合的交并比。如果两个集合的故障线路几乎一致它们距离就小可以合并。削减后的典型场景数量控制在10到20个之间。数量太少会丢掉小概率高风险场景数量太多又会拖慢应急决策模块。我一般先用削减前后ENS误差做校验误差不超过5%就认为削减结果可以接受。3.3 结果可视化和指标计算台风场景生成类项目结果可视化很重要。至少需要画三张图。第一张是台风路径与33节点位置的叠加图figure; hold on; plot(typhoon.path(:,1), typhoon.path(:,2), r-, LineWidth, 2); scatter(nodeCoord(:,1), nodeCoord(:,2), 40, filled); xlabel(x (km)); ylabel(y (km)); title(某台风路径与33节点配电网位置示意);第二张是线路故障概率的空间分布图可以用线路线宽或颜色深浅来表示故障概率高低。这样能直观看出台风路径附近的线路故障风险更高。第三张是典型场景的失电负荷柱状图每个柱子对应一个典型场景柱高表示该场景下的失电负荷旁边标注场景概率。指标计算方面重点算缺供电量ENSENS Σ_s p_s * Σ_i P_loss(i) * T_recovery(i)其中 p_s 是场景概率P_loss(i) 是第 i 个失电负荷节点的有功T_recovery(i) 是恢复时间。应急响应模块输出最大失电负荷、平均恢复时间、关键用户失电时长等指标这些都可以汇总到表格里输出。3.4 参数标定和仿真设置的心得参数标定是整个项目最花时间的部分也是最影响结果可信度的部分。我的建议是分三层台风参数、故障率参数、场景削减参数。台风参数方面最大风速和最大风速半径要基于历史台风数据标定最好找接近当地典型强度的台风过程做参考。如果没有实测数据至少要在仿真里做敏感性分析看结果对参数的敏感程度。故障率参数方面λ0 可以根据历史跳闸率统计α 比较难标定通常在两倍量级范围内调参并结合敏感性分析说明结果稳健性。仿真设置上几个经验值可以参考初始场景数取5000左右统计稳定性够用削减后保留场景数10到20个时间步长15到30分钟全过程时长一般覆盖台风过境的24到48小时每组仿真前用不同的随机种子运行3到5次检查指标波动幅度。4. 实操中踩过的坑与排查技巧4.1 风速模型选错导致故障率失真我最早用最简单的圆对称公式风速偏大故障率直接高了一个数量级生成的场景几乎每条线路都在跳一看就不合理。后来发现问题是忽略了最大风速半径之外的快速衰减以及台风移速叠加。排查方法很简单把风场计算得到节点风速画成时间曲线和气象站实测数据放在一起对比。实测可用的就对比不可用的就对比文献里典型的台风风速曲线形态。只要风速场形状不对后面的故障率必然失真这个环节值得第一个检查。4.2 场景削减时场景数怎么定场景削减的典型错误是贪图计算速度把场景砍得太狠。某个测试里我把5000个场景直接削减到5个结果ENS和未削减前差了将近20%尾部风险场景全部被抹掉。反过来削减到50个场景应急响应模块每轮计算要跑很久。实操建议先保留10个、15个、20个、30个场景做对比实验计算削减前后ENS误差曲线找到误差小于5%的最少场景数。大多数情况下15到20个场景就足够稳定但这个数值要在你自己的算例上验证不能直接抄。4.3 应急响应指标和故障场景的衔接问题最容易掉进去的坑是场景抽样后只记录故障线路编号没有及时计算失电负荷和用户数导致应急响应模块无从下手。解决这个问题一定要在场景数据结构里预留字段scenCell{m}.outageLoadMW outageLoadMW; scenCell{m}.affectedUsers affectedUsers; scenCell{m}.criticalLoadFlag criticalLoadFlag;另一个衔接问题是忽略转供。33节点配电网系统中一些联络开关可以让部分负荷从其他馈线恢复供电。应急响应里如果不考虑转供能力恢复时间会被高估。处理方法是把故障场景传给潮流计算模块用实际潮流结果决定哪些节点真正失电再更新故障特征。4.4 Matlab内存与随机数管理场景数一多Matlab经常报内存不足。原因往往是某个中间变量不小心把 N_scen 和 n_branch 都用满矩阵存瞬间爆炸。解决办法是场景用cell数组存稀疏结构风速场不要存全过程三维矩阵只保存节点风速矩阵33×时间步即可。这个规模的内存占用非常小。随机数管理同样是个细节。我最开始为了“复现”在每次抽样前都调用一次rng(seed)结果每个场景都抽到了相同的随机数序列生成的场景几乎全部一样。正确做法是在主脚本开头设置一次随机种子循环里直接用rand这样才能得到独立抽样。如果要做并行计算建议给每个并行worker单独配置随机种子流避免重复场景。最后再分享一点个人体会这个项目做完之后我最大的感受是台风故障建模的难点不在于把故障率公式写得多漂亮而在于让场景生成结果真正能被应急决策用起来。故障特征字段的设计一定要前置在写蒙特卡洛抽样之前就想清楚“后面应急响应要什么数据”。另外整个流程跑通之后后续扩展其实很方便可以把修复时间分布改成更精细的抢修队调度模型也可以把静态场景生成升级为时序恢复过程分析甚至接入移动储能应急接入策略。33节点算例只是一个起点方法框架本身是可以复用的。