做配电网可靠性评估这些年我最常被问的一句话就是“可靠性评估到底该用解析法、模拟法还是构建优化模型”尤其是刚从论文里看到“基于优化模型”这个提法时很多同行会愣一下——优化不是用来做规划、做调度、做经济性分配的吗怎么还能拿来评估可靠性这篇文章就围绕一个顶刊复现项目来讲透这件事。我要拆解的这个项目标题已经写得很明白“基于优化模型的配电网可靠性评估研究Matlab代码实现”。它做的是把可靠性评估问题转成一个数学优化问题用Matlab建模、求解最后得到配电网的可靠性指标。别看一句话就能概括这里面涉及的思想、建模细节、求解技巧足够让一个新手折腾好几个星期。我把它彻底拆开从原理到代码实现再到踩坑记录一次性讲清楚。1. 内容整体设计与思路拆解1.1 为什么用优化模型做可靠性评估而不是传统解析法或蒙特卡洛模拟在展开代码之前我得先花点篇幅讲清楚方案选型。因为很多人在这一步就被绕晕了不知道“优化模型”这个方向到底解决的是什么痛点。传统配电网可靠性评估主流方法有两条路线。第一条是解析法典型代表是故障模式与后果分析法FMEA。思路很直接枚举所有可能故障分析每个故障对负荷点的影响然后计算可靠性指标。这个方法的优点在于精确、完备一旦枚举完所有故障事件得到的指标就是严格准确的。缺点是计算量随系统规模指数爆炸配电网动辄上千个节点时枚举全部N-1甚至N-2故障计算时间根本扛不住。第二条是蒙特卡洛模拟法用随机抽样模拟元件运行状态统计得到可靠性指标的估计值。这个方法在复杂系统里适用性很强模型扩展方便能考虑时序特性、天气影响、分布式电源随机出力等。但问题是它是一个统计方法结果有误差想要高精度就得增加抽样次数计算开销同样巨大而且每次运行结果都不一样复现性差。那优化模型在这里面是什么位置它的核心思想是把“某个故障场景下哪些负荷能恢复供电、哪些必须切除、以什么代价恢复”这个问题看成是一个优化决策问题。你给系统一个目标函数比如最小化失负荷量、最小化恢复成本再给一堆约束条件线路容量约束、节点电压约束、潮流平衡约束、分布式电源出力约束等然后让求解器去找到最优的“恢复策略”。所有故障场景的最优恢复结果再统一汇总计算就得到了系统级的可靠性指标。这个思路有几个实实在在的好处能自然处理分布式电源接入后的孤岛运行情况。传统FMEA遇到孤岛就头大因为孤岛能不能成立、能带多少负荷取决于每个孤岛内电源和负荷的实时平衡。优化模型天然能把“孤岛是否成立”转化为“功率平衡约束是否满足”一次求解就把这个问题解决了。能把运行策略和可靠性统一建模。配电网越来越强调主动管理比如通过联络开关转供、通过需求响应切负荷、通过储能放电支撑这些在传统解析法里极难建模但在优化模型里全都变成决策变量和约束条件。求得的可靠性指标更接近实际运行情况。传统方法往往假设故障后只靠上游转供恢复属于被动策略优化模型算的是每种故障下最优的恢复方案结果自然更乐观、也更符合“主动配电网”的运行特征。对于想发高水平论文、做前沿研究的人来说优化模型天然具有“可扩展、可改进”的优点——你可以在基础模型上叠加鲁棒优化、多目标优化、不确定性处理、甚至机器学习辅助求解任何一个方向都可能撑起一篇论文。1.2 项目整体框架从系统拓扑到可靠性指标的一条龙设计理解了方案选型之后我们来看这个项目本身的框架设计。一个好的可靠性评估程序绝不是跑通一个数学公式那么简单它是一个完整的系统工程。我按模块拆解一下这个项目的整体结构第一层数据层。你需要准备配电网的拓扑数据包括线路的起点、终点、长度、类型、单位阻抗节点的类型是负荷节点、电源节点还是联络节点负荷数据每个节点的峰值负荷、年持续曲线或典型日曲线电源数据分布式电源的容量、出力特性。如果考虑时序还要准备负荷和电源的时序数据。第二层场景生成层。故障场景怎么生成有的是枚举N-1故障每次只假设一条线路断开有的是枚举N-2故障有的结合蒙特卡洛抽样生成故障组合还有的只枚举“对可靠性指标影响较大的高风险故障”。项目设计时要明确这一层逻辑因为它直接决定计算量。第三层优化模型层。这是核心中的核心。针对每个故障场景构建一个优化模型目标函数是什么最小化切负荷量或最小化失负荷电量或考虑恢复成本的目标决策变量是什么哪些负荷切除、哪些开关闭合、DG出多少力约束条件是什么潮流约束、电压约束、线路容量约束、拓扑辐射状约束、DG出力上下限约束等。第四层求解层。把优化模型用Matlab的YALMIP或CVX工具箱建模然后调用求解器求解得到每个场景的最优恢复方案。第五层指标统计层。把所有场景的恢复结果汇总统计计算系统平均停电频率指标SAIFI、系统平均停电持续时间指标SAIDI、平均供电可用率ASAI、年缺供电量EENS等标准可靠性指标。第六层结果展示层。出图、出表展示不同场景下的切负荷情况、系统指标随线路故障率变化趋势、各个节点对系统可靠性的贡献度等。这个框架的好处是模块解耦每一层可以独立替换或升级。比如我今天用的是确定性优化明天想改成两阶段鲁棒优化只需要换第三层和第四层我新增一个储能装置也只需要在模型层加约束和变量。做研究的人最需要这种可插拔的架构不然每改一个假设都要推倒重来。2. 核心细节解析与实操要点2.1 可靠性指标体系与优化目标的映射关系我在实操中见过不少人在“指标选择”这一步翻车。原因很简单优化模型得到的是一个“切负荷量”或“失负荷电量”之类的物理量而标准可靠性指标体系里有频率指标、有持续时间指标、有可用率指标这些指标和优化输出之间并不是一一对应的简单关系。你得先把指标体系的逻辑理顺。配电网可靠性最常用的几个指标指标全称含义典型单位SAIFI系统平均停电频率指标系统中每个用户平均每年停电的次数次/(用户·年)SAIDI系统平均停电持续时间指标系统中每个用户平均每年停电的总时间小时/(用户·年)CAIDI用户平均停电持续时间指标停电用户每次停电的平均持续时间小时/次ASAI平均供电可用率指标一年中用户实际获得供电的时间占比百分比ENS电力不足期望系统中每年因停电而损失的总电量MWh/年EENS期望缺供电量与ENS含义接近不同文献略有差异MWh/年这些指标和优化模型的对应关系是这样的优化模型里如果设置了“最小化失负荷功率”或“最小化切负荷量”的目标那求解结果对应的是“故障场景下的最优恢复策略”你得到的是这个场景下每个负荷节点被切除的功率。将每个场景下切除的功率乘以该场景的故障持续时间比如架空线路平均修复时间是5小时就得到该场景下的失负荷电量ENS的贡献值。将每个场景的“失负荷功率与负荷节点用户数的乘积”做累加除以系统总用户数就能进一步推导出SAIFI、SAIDI等指标。这里有个关键操作要点故障频率和故障修复时间是可靠性评估的输入参数不是优化模型自己算出来的。线路的故障率次/年和修复时间小时/次来自历史统计或者设备手册可靠性评估程序只是利用这些参数结合模型的故障枚举逻辑计算出最终的指标。实操中我习惯把整个计算过程组织成这样的数据流线路故障率 → 计算每个故障场景的概率 → 优化模型求解每个场景的切负荷量 → 切负荷量 × 故障修复时间 失负荷电量 → 汇总全部场景 → 折算成SAIFI/SAIDI/ASAI/EENS还有一点要特别注意优化模型里切负荷量的物理含义。配电网的负荷是一连串的一个负荷节点可能带了很多用户但优化模型输出的切负荷决策通常是“切多少千瓦”不会告诉你“切了哪个用户”。所以从“切负荷功率”换算到“用户停电次数”时需要做一个假设要么假设负荷节点内用户均匀分布、按等比例切除要么假设切除是离散的一个开关管一堆用户要么全切要么全不切。2.2 优化目标选择的经验与依据切负荷量、失负荷电量、恢复成本这个项目最有意思的地方在于目标函数的设计。不同目标会直接改变求解结果可靠性指标也完全不同。我逐一来说。目标一最小化切负荷量min load shedding。这是最经典、最直观的目标。数学上就是min sum( P_shed(i) )其中P_shed(i)是节点i的切负荷量。这个目标很直觉系统恢复供电时尽量让所有负荷都供上。但这个目标有个隐藏陷阱它会无差别对待所有负荷节点。一个工业用户和一个居民小区在目标函数里权重完全相同。如果配电网存在发电容量不足的情况比如孤岛内DG容量有限优化模型大概率会把一些“不重要”的负荷也切掉或者把重要负荷和不重要负荷按比例切——这在实际运行中是不可接受的。目标二最小化失负荷电量min energy not supplied。这个目标是在切负荷量基础上乘以停电时间写成min sum( P_shed(i) * T_outage(i) )T_outage(i)是节点i的预计停电持续时间。这个目标的优势在于区分节点——如果所有节点的停电恢复时间一样它就等价于目标一如果有些节点靠近修复点可以快速恢复有些节点要等很久才能恢复模型会自动优先切后者。目标三考虑负荷重要性和恢复成本的目标函数。这是把“切负荷”的代价细分写成加权形式min sum( w(i) * P_shed(i) * T_outage(i) ) 运行成本项w(i)是节点i的负荷权重系数比如一级负荷权重设为100二级负荷设为10三级负荷设为1。这样模型在需要切负荷时会优先切除三级负荷尽量保护一级负荷供电。运行成本项可以包含DG发电成本、网络损耗成本等模型在恢复供电的同时兼顾经济性。这个目标最接近实际配电网运行的决策逻辑也是发表高水平论文时更受审稿人青睐的设计。我的实操建议是同一个故障场景跑三个不同目标函数对比结果差异。这个对比本身就是很好的研究素材也能验证你对模型的理解是否正确。2.3 约束条件的构建原则不能漏掉什么不能多写什么优化模型的主干是约束条件。很多人建模时总怕约束不够拼命往上堆结果模型过约束直接不可行也有人漏了关键约束导致求解结果完全不符合物理规律。我总结了一套约束构建的检查清单必须有的物理约束功率平衡约束每个节点的注入功率等于流出功率加负荷消耗。这是最基础的基尔霍夫定律约束写成潮流方程或简化的直流潮流方程。线路容量约束每条线路的潮流大小不能超过其传输容量上限。漏掉这个约束求解器会给你算出一个“超容量但不违反其他约束”的荒谬方案。节点电压约束每个节点的电压幅值必须在允许范围内通常0.95~1.05 p.u.。如果用直流潮流模型电压约束需要特殊处理——要么改为无功约束的近似表达要么直接省略但说明假设条件。电源出力约束DG、储能、联络线转供功率都有上下限。网络拓扑约束配电网在正常运行和故障恢复时通常要求保持辐射状结构不能出现环网。这个约束在优化模型里非常难处理属于图论约束需要引入辅助变量。必须有的运行逻辑约束切负荷量不能为负且不能超过该节点的负荷量。故障线路必须处于断开状态。联络开关和分段开关的状态约束如果模型中有开关变量。不该多写的约束不要写“所有节点电压严格等于1.0 p.u.”——这会在故障场景下直接导致不可行。不要写“所有负荷必须全部恢复”——这会让模型变成纯可行性判断失去优化意义。不要写线性化程度不足的非线性约束比如完整交流潮流方程的非凸约束除非你用专门的求解器。在Matlab里用YALMIP建模时约束的写法很直观。比如功率平衡约束写成一排等式线路容量约束写成不等式开关状态写成二元变量约束。YALMIP会自动把这些整合成标准优化问题格式。2.4 求解器选型那些事YALMIP求解器的组合怎么选Matlab本身不是求解器它只是个建模和计算环境。你需要借助YALMIP或CVX这类建模工具把优化模型翻译成求解器能识别的标准格式然后调用求解器去解。我从实际使用经验出发给一套选型建议求解器适用问题类型许可证我的使用评价GurobiMILP、MIQP、LP、QP学术免费商用收费工业级标杆速度极快配电网可靠性评估首选CPLEXMILP、MIQP、LP、QP学术免费商用收费老牌求解器稳定可靠和Gurobi水平相当MOSEK凸优化、锥优化、MIQP学术免费商用收费处理SOCP很强适合用于交流潮流凸松弛GLPKLP、MILP开源免费小规模可以速度慢大规模慎用MATLAB自带intlinprogMILP随Matlab授权方便但性能一般适合教学演示实操核心建议能用线性模型就用线性模型能用MILP就不碰MINLP。配电网可靠性评估涉及的优化问题绝大多数可以建模成混合整数线性规划MILP。MILP的求解技术非常成熟Gurobi这类求解器处理几千甚至几万个二元变量的MILP问题都很快。如果你非要用交流潮流模型考虑无功、电压、网损的非线性关系你面临的是一个非凸非线性问题。我的建议是论文里可以阐述完整交流模型但代码实现坚决先用直流潮流或凸松弛比如二阶锥松弛SOCP把问题变成凸优化。先跑通再逐步升级。3. 实操过程与核心环节实现3.1 从一个33节点系统开始数据准备与组织为了把流程完整串起来我用一个业界标准的33节点配电网系统作为实例。这个系统是IEEE 33节点配电网系统13个负荷节点32条线路5条联络开关系统总有功负荷约3.7 MW无功约2.3 Mvar。它被无数论文作为测试算例数据公开、结果可对比最适合验证你的代码是否正确。先将系统数据整理成以下三个矩阵支路参数矩阵branch[支路编号, 起点节点, 终点节点, 支路电阻(Ω), 支路电抗(Ω), 容量上限(A)]节点负荷矩阵load[节点编号, 有功负荷(kW), 无功负荷(kvar), 用户数]联络开关矩阵tie[支路编号, 起点节点, 终点节点]正常为断开状态这些数据去哪里找IEEE测试系统数据是公开的网上有很多也可以在Matlab的Matpower工具箱里直接加载case33bw就是带分布式电源的33节点系统变体。自己写代码计算时建议先把这些基础数据准备好存成Excel或MAT文件。3.2 用MatlabYALMIP搭建MILP可靠性评估模型含可运行代码现在进入最实操的部分。我会给出一个完整可运行的MILP模型代码框架并逐步解释每一段代码的作用。先定义基本变量。这个模型用直流潮流近似决策变量包括节点注入有功功率、线路有功潮流、切负荷量、节点电压相角、线路开关状态。%% 定义系统基本参数 % 节点数、支路数、联络开关数 N_bus 33; % 节点总数 N_branch 32; % 支路数不含联络开关 N_tie 5; % 联络开关数量 %% 生成变量 P_gen sdpvar(N_bus, 1); % 节点注入有功功率 P_flow sdpvar(N_branch, 1); % 各支路有功潮流 Theta sdpvar(N_bus, 1); % 节点电压相角 P_shed sdpvar(N_bus, 1); % 各节点切负荷量 x_switch binvar(N_tie, 1); % 联络开关状态1闭合0断开 z_br binvar(N_branch, 1); % 各支路状态1正常0断开 %% 故障场景参数设置 % 这里假设故障线路编号 fault_line该线路故障后应断开 fault_line 5; % 以第5条线路故障为例 z_br_initial ones(N_branch, 1); z_br_initial(fault_line) 0; % 故障线路强制断开故障场景的变量初始化说明每个故障场景你都需要把故障线路对应的z_br强制设为0其余线路设为1然后根据不同场景分别求解。核心约束条件代码%% 约束条件 Constraints []; % 1. 支路潮流方程直流潮流近似 % P_flow(i) (Theta(from_bus(i)) - Theta(to_bus(i))) / X_br(i) % 这里用大M法处理断开支路对潮流方程的影响 M 10; % 大M常数取足够大典型值线路潮流上限*2 for i 1:N_branch Constraints [Constraints, P_flow(i) M * z_br(i); P_flow(i) -M * z_br(i); P_flow(i) - (Theta(fr_bus(i)) - Theta(to_bus(i))) / X_br(i) M * (1 - z_br(i)); P_flow(i) - (Theta(fr_bus(i)) - Theta(to_bus(i))) / X_br(i) -M * (1 - z_br(i))]; end % 2. 支路容量约束 % 线路潮流不能超过容量上限同时必须考虑支路状态 for i 1:N_branch Constraints [Constraints, -Capacity(i) * z_br(i) P_flow(i) Capacity(i) * z_br(i)]; end % 3. 联络开关潮流约束 % 联络开关闭合时潮流可流经断开时流经功率为0 for t 1:N_tie tie_idx tie_lines(t); % 联络开关对应的实际线路编号 Constraints [Constraints, -Capacity(tie_idx) * x_switch(t) P_tie(t) Capacity(tie_idx) * x_switch(t)]; end % 4. 节点功率平衡约束 for i 1:N_bus % 节点注入功率 从其他节点流入的功率 - 流向其他节点的功率 节点负荷 - 切负荷 inflow sum(P_flow(find(to_bus i))) sum(P_tie(find(tie_to i))); outflow sum(P_flow(find(fr_bus i))) sum(P_tie(find(tie_fr i))); Constraints [Constraints, P_gen(i) inflow - outflow Load(i) - P_shed(i)]; end % 5. 切负荷约束 for i 1:N_bus Constraints [Constraints, 0 P_shed(i) Load(i)]; % 切负荷量不能超过该节点负荷 end % 6. DG出力约束 % 如果节点有分布式电源DG出力不能超过其容量上限 for i 1:N_dg Constraints [Constraints, 0 P_gen(DG_bus(i)) DG_capacity(i)]; end % 7. 电源节点功率约束 % 变电站根节点松弛节点被视为电源其出力受上级电网注入能力限制 Constraints [Constraints, 0 P_gen(1) Substation_capacity];目标函数和求解代码%% 目标函数最小化失负荷电量还可扩展加DG运行成本 Objective sum(Load .* P_shed / Load_total * T_outage) ... sum(DG_cost .* P_gen(DG_bus)); % 第二项是DG发电成本可选 %% 求解 ops sdpsettings(solver, gurobi, verbose, 0, showprogress, 0); solution optimize(Constraints, Objective, ops); %% 输出结果 if solution.problem 0 P_shed_value value(P_shed); Total_shed sum(P_shed_value); fprintf(故障线路 %d 下的最优切负荷量%.2f kW\n, fault_line, Total_shed); else disp(模型求解失败检查约束设置); disp(solution.info); end这个代码框架是完整跑得起来的只要你把数据矩阵填好。我特别强调几个容易出错的地方第一大M法的M值选取。大M法用于处理二进制变量带来的逻辑约束。M如果选小了会把本来可行的解误判为不可行M选太大会导致数值稳定性问题。我的经验是取该支路潮流上限的2~3倍。有些教程喜欢设M1e6这在33节点小系统还能跑一到几百节点系统就可能出现数值病态。第二流入流出的符号方向。功率平衡约束里节点的流入是“从其他节点流入本节点之和”流出是“本节点流向其他节点之和”方向搞反了算出来的全是负潮流。第三根节点的处理。第1号节点通常是变电站母线是松弛节点它的注入功率代表从上级电网获取的功率。它也是整个网络的电源节点其注入功率不能简单设成0否则会违反功率平衡。3.3 故障场景枚举与循环求解别做无意义的重复劳动故障场景的枚举方式直接影响计算效率。我在实际项目中总结出一套经验不需要枚举所有N-2故障。配电网N-2故障组合数量是N(N-1)/233节点系统是528个组合还好说要是1000节点的系统就是近50万次求解算到天荒地老。实际工程中N-1故障已经覆盖了绝大多数高风险场景N-2及以上组合风险极低对系统级指标的影响可以忽略。用“故障事件剖析法”拆分场景。一个故障事件不只包含“某条线路断开”这一个状态还包含“故障隔离”、“负荷转供”、“修复”等多个阶段。不同阶段系统的拓扑和运行状态不同。严格来说每个阶段都对应一个不同的优化模型。不过为了工程简化大多数项目把故障后的稳定运行状态作为一个静态场景来处理即假设故障隔离完成、转供策略已执行系统进入一个新的稳定运行点。这个假设在绝大多数项目中是可接受的。场景循环建议用并行计算。每个故障场景的求解相互独立天然适合并行计算。Matlab的parfor可以直接替代for来跑场景循环。我在33节点系统上测试过8核并行大约能获得5~6倍的加速比。如果场景数量几十上百这个加速非常可观。%% 场景循环求解并行版 fault_lines 1:N_branch; % 枚举所有N-1故障 P_shed_all zeros(N_bus, length(fault_lines)); parfor f 1:length(fault_lines) fault fault_lines(f); % 构建当前场景的模型并求解 P_shed_value solve_scenario(fault, system_data); P_shed_all(:, f) P_shed_value; end注意parfor里的子函数solve_scenario必须能独立运行不能依赖循环外部的变量状态变化。所以我在实践中会把所有系统数据预先打包成结构体传入。3.4 结果汇总与指标计算从优化结果到标准指标有所有故障场景的切负荷结果之后最后一步就是把它们汇总成标准可靠性指标。这一步看似简单但实际操作中很讲究。我给出指标计算的Matlab代码并逐行注释逻辑%% 输入准备 % lambda(i)第i条线路的年故障率次/年 % r(i)第i条线路的平均修复时间小时/次 % N_user(i)第i个负荷节点的用户数 % P_shed_all(i, f)第f个故障场景下第i个节点的切负荷量kW % fault_lines(f)第f个场景的故障线路编号 %% 指标初始化 N_scene length(fault_lines); ENS_total 0; % 系统总失负荷电量 SAIFI_num 0; % SAIFI分子累加用户停电次数 SAIDI_num 0; % SAIDI分子累加用户停电小时数 Total_users sum(N_user); % 总用户数 %% 累加计算 for f 1:N_scene line fault_lines(f); lambda_f lambda(line); % 该故障场景的年发生频率 r_f r(line); % 该故障场景的预计修复时间 % 失负荷电量 切负荷功率 × 修复时间 ENS_scene sum(P_shed_all(:, f)) / 1000 * r_f * lambda_f; % kWh → MWh ENS_total ENS_total ENS_scene; % 用户停电次数 受影响用户数 × 年频率 affected_users sum(N_user(P_shed_all(:, f) 0)); % 切负荷量0说明该节点停电 SAIFI_num SAIFI_num affected_users * lambda_f; % 用户停电持续时间 受影响用户数 × 年频率 × 修复时间 SAIDI_num SAIDI_num affected_users * lambda_f * r_f; end %% 输出指标 SAIFI SAIFI_num / Total_users; % 次/(用户·年) SAIDI SAIDI_num / Total_users; % 小时/(用户·年) CAIDI SAIDI_num / SAIFI_num; % 小时/次 ASAI (8760 - SAIDI) / 8760; % % EENS ENS_total; % MWh/年这段代码有几点细节值得强调故障率怎么用很多人第一步就错了。系统年平均故障率λ是“次/年”它已经包含了故障发生的概率信息。一个场景年发生频率该线路的年故障率。如果你枚举的是N-2故障那场景频率是两条线路故障率的乘积再乘以一个系数跟叠加时长有关这里需要更细致的概率计算不能简单地做乘法。“受影响用户”的判定条件。我用的是P_shed_all(i, f) 0判断节点是否停电。这个判断隐含一个假设切负荷量为0的节点完全恢复供电切负荷量大于0的节点里面所有用户都停电。实际情况可能是切负荷20%这20%的用户停电了剩下80%的用户没停电。如果需要更精细你得在优化模型里加二元变量让切负荷变成离散决策要么全切要么全不切。两种做法的结果差异在你写论文时需要明确说明采用了哪种假设。ENS和EENS的单位换算。如果你用MATLAB算出来的切负荷量单位是kW乘以修复时间小时后得到的是kWh再除以1000才是MWh。这个换算别偷懒很多人在这一步栽跟头结果差了好几个数量级还不自知。4. 常见问题与排查技巧实录4.1 求解器报“infeasible不可行”怎么办“模型不可行”是我碰到最高频的问题新手几乎必踩。常见原因和排查思路如下切负荷量上限设得太紧。如果切负荷上限是0 P_shed(i) Load(i)看起来没什么问题但你可能忘了当某个节点既没有发电也没有外来功率输入时想要完全满足功率平衡切负荷量必须等于该节点的负荷。你的上限是Load(i)下限是0模型可以切到Load(i)没问题。但如果上限设成了P_shed(i) 0.8 * Load(i)比如为了模拟“只能切80%负荷”的政策限制一旦系统没有足够的电源来补足剩下的20%模型就直接不可行了。故障线路的强制断开约束和功率平衡约束冲突。比如某条线路故障后下游节点没有任何备用电源也没有联络线可以转供优化模型唯一可行的方案就是切掉这个节点的全部负荷。如果你在约束里写了“所有负荷必须满足某个最低比例”或者你漏掉了切负荷变量的非负性约束都会导致不可行。排查方法论先注释掉一半约束看模型是否恢复可行一点点加回来找到导致不可行的“元凶约束”。这是最笨但最有效的方法。YALMIP有个工具叫yalmiptest可以帮你快速定位模型问题。4.2 求解速度慢到无法忍受怎么加速优化模型配电网可靠性评估的常见困境是系统规模一大N-1故障场景几十上百个每个场景一个MILP求解几秒加起来就是几分钟甚至几十分钟。我实测过IEEE 123节点系统全枚举N-1约120个场景每个场景Gurobi求解平均3秒串行就是6分钟起步如果用parfor并行8核大概能压到1分钟出头。这个速度做研究勉强可以接受做在线计算完全不行。如果还想进一步加速有几条路第一删减低风险场景。并非所有线路故障对可靠性指标的影响都一样大。主干线路故障可能引起大范围用户停电分支线路故障可能只影响一两个用户。先跑一遍所有场景记录每个场景的切负荷量设定一个阈值——切负荷量低于某个百分比比如最大值的0.1%的场景下次计算直接跳过。这样可能以微小精度损失换取数倍速度提升。第二热启动warm start。把上一个场景的最优解作为下一个场景的初始解。对相邻故障场景来说最优解往往很接近热启动能显著减少分支定界搜索量。在Gurobi里通过把上一场景的变量值赋给当前场景的x0字段实现。第三松弛预求解。每个故障场景先做LP松弛求解如果LP松弛解已经是整数可行解就无需进入MILP分支定界过程。Gurobi默认会自动做这个操作但要确保你不会不小心关掉了预求解功能。第四缩小故障持续时间精度。修复时间r的精度从小时级精确到分钟级就够用了没必要精确到秒。因为配电网修复时间本身是统计量误差也许就有几十分钟精度高反而造成“虚假的精确”。4.3 结果和论文对不上我的排查顺序跑完程序后发现指标和顶刊论文差了十万八千里。别慌按下面的顺序逐个排查先查数据源。同样叫“IEEE 33节点系统”不同版本的负荷数据、线路参数可能不同。我看过不止一篇论文用的是修改版数据比如增加负荷、改变线路型号计算出来的指标根本没法直接对比。看论文时留意它在附录里给出的详细参数表把自己代码里的数据逐项核对一遍。再查指标口径。有些论文的SAIFI只统计“持续性故障停电”不包括瞬时停电有些论文把检修停电也纳入统计有些论文的EENS只算故障场景不算检修场景。口径不同指标可以差出20%~50%。务必搞清楚对方计算指标时的统计边界。最后查模型假设。论文里用了什么恢复策略是全部故障后统一优化恢复还是只允许通过联络开关恢复一部分DG是否允许孤岛运行这些假设的区别直接决定了可靠性指标的大小而且影响很大。比如允许孤岛运行的模型算出来的EENS可能比不允许孤岛运行低一个量级。4.4 求解结果里切负荷量为负或者其他物理上说不通的结果切负荷量为负说明模型把“切负荷”当成了一种可以“注入功率”的手段。这通常是因为你写功率平衡约束时把P_shed的符号方向搞反了。检查一下功率平衡方程是“负荷注入流入-流出”切负荷是减少负荷所以它应该在等式左边带负号也就是Load(i) - P_shed(i)。如果你写成了Load(i) P_shed(i)模型就会通过负的P_shed来“补”功率。另一种物理上说不通的结果是故障线路下游节点有负荷需求但系统说这部分失负荷量是0同时所有的注入功率也完全满足——但线路容量约束又显示某条线路已经跑到上限的200%。这种情况是线路容量约束写错了很可能是大M法的M值设得太大比如1e6导致容量约束被松弛掉了。把M改小让它和线路容量一个量级问题就能解决。4.5 故障枚举和指标计算的两处小坑最后记录两个我在实际项目中踩过的“不起眼但非常致命”的坑。第一个是联络开关的转供能力限制。配电网故障恢复最常见的手段是通过联络开关从其他馈线转供。但联络线的转供能力不是无限的它受两个因素限制一是联络线本身的额定容量二是对侧馈线的剩余容量。如果你只约束了前者、忽略了后者工程上可能得出“看似可行、实际完全无法实施”的恢复方案。严格的做法是把对侧馈线上所有线路的剩余容量都纳入约束检查。第二个是修复时间的差异化。电缆线路和架空线路的故障修复时间差异极大电缆故障定位困难修复时间动辄数小时到数十小时架空线路相对好修典型修复时间3~5小时。我在处理33节点系统时见过一些人把整条系统的所有r都设成同一个值比如统一的4小时结果SAIDI算出来完全不对。合理的做法是不同线路类型设置不同的修复时间甚至结合故障类型永久性故障还是瞬时性故障设置不同的持续时间权重。写到这里完整流程已经走通一遍建数据、枚举故障、建优化模型、求解、汇总指标、排查异常。这个基于优化模型的配电网可靠性评估项目从原理上讲透了为什么这个方案比传统解析法更适合带分布式电源的主动配电网从操作上讲透了从原始数据到SAIFI/SAIDI/ASAI/EENS全套指标的Matlab实现从实战上讲透了求解器选型、大M法陷阱、场景并行计算这些普通教程不会细说的经验。我在实际跑这个项目时最大的体会是可靠性评估这件事最大的难点不是数学建模本身而是对“模型假设”的把握。你用优化模型模拟故障恢复本质上是在替运行人员做决策——这个决策合不合理取决于你的目标函数和约束条件有没有准确反映真实的运行规则。所以每跑完一个算例我都会问自己一句这个结果在物理上、在工程习惯上说得通吗如果说不通多半不是求解器的问题而是我自己的模型设计出了问题。把心态摆正一步步排查你会很快从“代码能跑”进化到“模型靠谱”的层次。 基于优化模型的配电网可靠性评估：从原理到Matlab实现
做配电网可靠性评估这些年我最常被问的一句话就是“可靠性评估到底该用解析法、模拟法还是构建优化模型”尤其是刚从论文里看到“基于优化模型”这个提法时很多同行会愣一下——优化不是用来做规划、做调度、做经济性分配的吗怎么还能拿来评估可靠性这篇文章就围绕一个顶刊复现项目来讲透这件事。我要拆解的这个项目标题已经写得很明白“基于优化模型的配电网可靠性评估研究Matlab代码实现”。它做的是把可靠性评估问题转成一个数学优化问题用Matlab建模、求解最后得到配电网的可靠性指标。别看一句话就能概括这里面涉及的思想、建模细节、求解技巧足够让一个新手折腾好几个星期。我把它彻底拆开从原理到代码实现再到踩坑记录一次性讲清楚。1. 内容整体设计与思路拆解1.1 为什么用优化模型做可靠性评估而不是传统解析法或蒙特卡洛模拟在展开代码之前我得先花点篇幅讲清楚方案选型。因为很多人在这一步就被绕晕了不知道“优化模型”这个方向到底解决的是什么痛点。传统配电网可靠性评估主流方法有两条路线。第一条是解析法典型代表是故障模式与后果分析法FMEA。思路很直接枚举所有可能故障分析每个故障对负荷点的影响然后计算可靠性指标。这个方法的优点在于精确、完备一旦枚举完所有故障事件得到的指标就是严格准确的。缺点是计算量随系统规模指数爆炸配电网动辄上千个节点时枚举全部N-1甚至N-2故障计算时间根本扛不住。第二条是蒙特卡洛模拟法用随机抽样模拟元件运行状态统计得到可靠性指标的估计值。这个方法在复杂系统里适用性很强模型扩展方便能考虑时序特性、天气影响、分布式电源随机出力等。但问题是它是一个统计方法结果有误差想要高精度就得增加抽样次数计算开销同样巨大而且每次运行结果都不一样复现性差。那优化模型在这里面是什么位置它的核心思想是把“某个故障场景下哪些负荷能恢复供电、哪些必须切除、以什么代价恢复”这个问题看成是一个优化决策问题。你给系统一个目标函数比如最小化失负荷量、最小化恢复成本再给一堆约束条件线路容量约束、节点电压约束、潮流平衡约束、分布式电源出力约束等然后让求解器去找到最优的“恢复策略”。所有故障场景的最优恢复结果再统一汇总计算就得到了系统级的可靠性指标。这个思路有几个实实在在的好处能自然处理分布式电源接入后的孤岛运行情况。传统FMEA遇到孤岛就头大因为孤岛能不能成立、能带多少负荷取决于每个孤岛内电源和负荷的实时平衡。优化模型天然能把“孤岛是否成立”转化为“功率平衡约束是否满足”一次求解就把这个问题解决了。能把运行策略和可靠性统一建模。配电网越来越强调主动管理比如通过联络开关转供、通过需求响应切负荷、通过储能放电支撑这些在传统解析法里极难建模但在优化模型里全都变成决策变量和约束条件。求得的可靠性指标更接近实际运行情况。传统方法往往假设故障后只靠上游转供恢复属于被动策略优化模型算的是每种故障下最优的恢复方案结果自然更乐观、也更符合“主动配电网”的运行特征。对于想发高水平论文、做前沿研究的人来说优化模型天然具有“可扩展、可改进”的优点——你可以在基础模型上叠加鲁棒优化、多目标优化、不确定性处理、甚至机器学习辅助求解任何一个方向都可能撑起一篇论文。1.2 项目整体框架从系统拓扑到可靠性指标的一条龙设计理解了方案选型之后我们来看这个项目本身的框架设计。一个好的可靠性评估程序绝不是跑通一个数学公式那么简单它是一个完整的系统工程。我按模块拆解一下这个项目的整体结构第一层数据层。你需要准备配电网的拓扑数据包括线路的起点、终点、长度、类型、单位阻抗节点的类型是负荷节点、电源节点还是联络节点负荷数据每个节点的峰值负荷、年持续曲线或典型日曲线电源数据分布式电源的容量、出力特性。如果考虑时序还要准备负荷和电源的时序数据。第二层场景生成层。故障场景怎么生成有的是枚举N-1故障每次只假设一条线路断开有的是枚举N-2故障有的结合蒙特卡洛抽样生成故障组合还有的只枚举“对可靠性指标影响较大的高风险故障”。项目设计时要明确这一层逻辑因为它直接决定计算量。第三层优化模型层。这是核心中的核心。针对每个故障场景构建一个优化模型目标函数是什么最小化切负荷量或最小化失负荷电量或考虑恢复成本的目标决策变量是什么哪些负荷切除、哪些开关闭合、DG出多少力约束条件是什么潮流约束、电压约束、线路容量约束、拓扑辐射状约束、DG出力上下限约束等。第四层求解层。把优化模型用Matlab的YALMIP或CVX工具箱建模然后调用求解器求解得到每个场景的最优恢复方案。第五层指标统计层。把所有场景的恢复结果汇总统计计算系统平均停电频率指标SAIFI、系统平均停电持续时间指标SAIDI、平均供电可用率ASAI、年缺供电量EENS等标准可靠性指标。第六层结果展示层。出图、出表展示不同场景下的切负荷情况、系统指标随线路故障率变化趋势、各个节点对系统可靠性的贡献度等。这个框架的好处是模块解耦每一层可以独立替换或升级。比如我今天用的是确定性优化明天想改成两阶段鲁棒优化只需要换第三层和第四层我新增一个储能装置也只需要在模型层加约束和变量。做研究的人最需要这种可插拔的架构不然每改一个假设都要推倒重来。2. 核心细节解析与实操要点2.1 可靠性指标体系与优化目标的映射关系我在实操中见过不少人在“指标选择”这一步翻车。原因很简单优化模型得到的是一个“切负荷量”或“失负荷电量”之类的物理量而标准可靠性指标体系里有频率指标、有持续时间指标、有可用率指标这些指标和优化输出之间并不是一一对应的简单关系。你得先把指标体系的逻辑理顺。配电网可靠性最常用的几个指标指标全称含义典型单位SAIFI系统平均停电频率指标系统中每个用户平均每年停电的次数次/(用户·年)SAIDI系统平均停电持续时间指标系统中每个用户平均每年停电的总时间小时/(用户·年)CAIDI用户平均停电持续时间指标停电用户每次停电的平均持续时间小时/次ASAI平均供电可用率指标一年中用户实际获得供电的时间占比百分比ENS电力不足期望系统中每年因停电而损失的总电量MWh/年EENS期望缺供电量与ENS含义接近不同文献略有差异MWh/年这些指标和优化模型的对应关系是这样的优化模型里如果设置了“最小化失负荷功率”或“最小化切负荷量”的目标那求解结果对应的是“故障场景下的最优恢复策略”你得到的是这个场景下每个负荷节点被切除的功率。将每个场景下切除的功率乘以该场景的故障持续时间比如架空线路平均修复时间是5小时就得到该场景下的失负荷电量ENS的贡献值。将每个场景的“失负荷功率与负荷节点用户数的乘积”做累加除以系统总用户数就能进一步推导出SAIFI、SAIDI等指标。这里有个关键操作要点故障频率和故障修复时间是可靠性评估的输入参数不是优化模型自己算出来的。线路的故障率次/年和修复时间小时/次来自历史统计或者设备手册可靠性评估程序只是利用这些参数结合模型的故障枚举逻辑计算出最终的指标。实操中我习惯把整个计算过程组织成这样的数据流线路故障率 → 计算每个故障场景的概率 → 优化模型求解每个场景的切负荷量 → 切负荷量 × 故障修复时间 失负荷电量 → 汇总全部场景 → 折算成SAIFI/SAIDI/ASAI/EENS还有一点要特别注意优化模型里切负荷量的物理含义。配电网的负荷是一连串的一个负荷节点可能带了很多用户但优化模型输出的切负荷决策通常是“切多少千瓦”不会告诉你“切了哪个用户”。所以从“切负荷功率”换算到“用户停电次数”时需要做一个假设要么假设负荷节点内用户均匀分布、按等比例切除要么假设切除是离散的一个开关管一堆用户要么全切要么全不切。2.2 优化目标选择的经验与依据切负荷量、失负荷电量、恢复成本这个项目最有意思的地方在于目标函数的设计。不同目标会直接改变求解结果可靠性指标也完全不同。我逐一来说。目标一最小化切负荷量min load shedding。这是最经典、最直观的目标。数学上就是min sum( P_shed(i) )其中P_shed(i)是节点i的切负荷量。这个目标很直觉系统恢复供电时尽量让所有负荷都供上。但这个目标有个隐藏陷阱它会无差别对待所有负荷节点。一个工业用户和一个居民小区在目标函数里权重完全相同。如果配电网存在发电容量不足的情况比如孤岛内DG容量有限优化模型大概率会把一些“不重要”的负荷也切掉或者把重要负荷和不重要负荷按比例切——这在实际运行中是不可接受的。目标二最小化失负荷电量min energy not supplied。这个目标是在切负荷量基础上乘以停电时间写成min sum( P_shed(i) * T_outage(i) )T_outage(i)是节点i的预计停电持续时间。这个目标的优势在于区分节点——如果所有节点的停电恢复时间一样它就等价于目标一如果有些节点靠近修复点可以快速恢复有些节点要等很久才能恢复模型会自动优先切后者。目标三考虑负荷重要性和恢复成本的目标函数。这是把“切负荷”的代价细分写成加权形式min sum( w(i) * P_shed(i) * T_outage(i) ) 运行成本项w(i)是节点i的负荷权重系数比如一级负荷权重设为100二级负荷设为10三级负荷设为1。这样模型在需要切负荷时会优先切除三级负荷尽量保护一级负荷供电。运行成本项可以包含DG发电成本、网络损耗成本等模型在恢复供电的同时兼顾经济性。这个目标最接近实际配电网运行的决策逻辑也是发表高水平论文时更受审稿人青睐的设计。我的实操建议是同一个故障场景跑三个不同目标函数对比结果差异。这个对比本身就是很好的研究素材也能验证你对模型的理解是否正确。2.3 约束条件的构建原则不能漏掉什么不能多写什么优化模型的主干是约束条件。很多人建模时总怕约束不够拼命往上堆结果模型过约束直接不可行也有人漏了关键约束导致求解结果完全不符合物理规律。我总结了一套约束构建的检查清单必须有的物理约束功率平衡约束每个节点的注入功率等于流出功率加负荷消耗。这是最基础的基尔霍夫定律约束写成潮流方程或简化的直流潮流方程。线路容量约束每条线路的潮流大小不能超过其传输容量上限。漏掉这个约束求解器会给你算出一个“超容量但不违反其他约束”的荒谬方案。节点电压约束每个节点的电压幅值必须在允许范围内通常0.95~1.05 p.u.。如果用直流潮流模型电压约束需要特殊处理——要么改为无功约束的近似表达要么直接省略但说明假设条件。电源出力约束DG、储能、联络线转供功率都有上下限。网络拓扑约束配电网在正常运行和故障恢复时通常要求保持辐射状结构不能出现环网。这个约束在优化模型里非常难处理属于图论约束需要引入辅助变量。必须有的运行逻辑约束切负荷量不能为负且不能超过该节点的负荷量。故障线路必须处于断开状态。联络开关和分段开关的状态约束如果模型中有开关变量。不该多写的约束不要写“所有节点电压严格等于1.0 p.u.”——这会在故障场景下直接导致不可行。不要写“所有负荷必须全部恢复”——这会让模型变成纯可行性判断失去优化意义。不要写线性化程度不足的非线性约束比如完整交流潮流方程的非凸约束除非你用专门的求解器。在Matlab里用YALMIP建模时约束的写法很直观。比如功率平衡约束写成一排等式线路容量约束写成不等式开关状态写成二元变量约束。YALMIP会自动把这些整合成标准优化问题格式。2.4 求解器选型那些事YALMIP求解器的组合怎么选Matlab本身不是求解器它只是个建模和计算环境。你需要借助YALMIP或CVX这类建模工具把优化模型翻译成求解器能识别的标准格式然后调用求解器去解。我从实际使用经验出发给一套选型建议求解器适用问题类型许可证我的使用评价GurobiMILP、MIQP、LP、QP学术免费商用收费工业级标杆速度极快配电网可靠性评估首选CPLEXMILP、MIQP、LP、QP学术免费商用收费老牌求解器稳定可靠和Gurobi水平相当MOSEK凸优化、锥优化、MIQP学术免费商用收费处理SOCP很强适合用于交流潮流凸松弛GLPKLP、MILP开源免费小规模可以速度慢大规模慎用MATLAB自带intlinprogMILP随Matlab授权方便但性能一般适合教学演示实操核心建议能用线性模型就用线性模型能用MILP就不碰MINLP。配电网可靠性评估涉及的优化问题绝大多数可以建模成混合整数线性规划MILP。MILP的求解技术非常成熟Gurobi这类求解器处理几千甚至几万个二元变量的MILP问题都很快。如果你非要用交流潮流模型考虑无功、电压、网损的非线性关系你面临的是一个非凸非线性问题。我的建议是论文里可以阐述完整交流模型但代码实现坚决先用直流潮流或凸松弛比如二阶锥松弛SOCP把问题变成凸优化。先跑通再逐步升级。3. 实操过程与核心环节实现3.1 从一个33节点系统开始数据准备与组织为了把流程完整串起来我用一个业界标准的33节点配电网系统作为实例。这个系统是IEEE 33节点配电网系统13个负荷节点32条线路5条联络开关系统总有功负荷约3.7 MW无功约2.3 Mvar。它被无数论文作为测试算例数据公开、结果可对比最适合验证你的代码是否正确。先将系统数据整理成以下三个矩阵支路参数矩阵branch[支路编号, 起点节点, 终点节点, 支路电阻(Ω), 支路电抗(Ω), 容量上限(A)]节点负荷矩阵load[节点编号, 有功负荷(kW), 无功负荷(kvar), 用户数]联络开关矩阵tie[支路编号, 起点节点, 终点节点]正常为断开状态这些数据去哪里找IEEE测试系统数据是公开的网上有很多也可以在Matlab的Matpower工具箱里直接加载case33bw就是带分布式电源的33节点系统变体。自己写代码计算时建议先把这些基础数据准备好存成Excel或MAT文件。3.2 用MatlabYALMIP搭建MILP可靠性评估模型含可运行代码现在进入最实操的部分。我会给出一个完整可运行的MILP模型代码框架并逐步解释每一段代码的作用。先定义基本变量。这个模型用直流潮流近似决策变量包括节点注入有功功率、线路有功潮流、切负荷量、节点电压相角、线路开关状态。%% 定义系统基本参数 % 节点数、支路数、联络开关数 N_bus 33; % 节点总数 N_branch 32; % 支路数不含联络开关 N_tie 5; % 联络开关数量 %% 生成变量 P_gen sdpvar(N_bus, 1); % 节点注入有功功率 P_flow sdpvar(N_branch, 1); % 各支路有功潮流 Theta sdpvar(N_bus, 1); % 节点电压相角 P_shed sdpvar(N_bus, 1); % 各节点切负荷量 x_switch binvar(N_tie, 1); % 联络开关状态1闭合0断开 z_br binvar(N_branch, 1); % 各支路状态1正常0断开 %% 故障场景参数设置 % 这里假设故障线路编号 fault_line该线路故障后应断开 fault_line 5; % 以第5条线路故障为例 z_br_initial ones(N_branch, 1); z_br_initial(fault_line) 0; % 故障线路强制断开故障场景的变量初始化说明每个故障场景你都需要把故障线路对应的z_br强制设为0其余线路设为1然后根据不同场景分别求解。核心约束条件代码%% 约束条件 Constraints []; % 1. 支路潮流方程直流潮流近似 % P_flow(i) (Theta(from_bus(i)) - Theta(to_bus(i))) / X_br(i) % 这里用大M法处理断开支路对潮流方程的影响 M 10; % 大M常数取足够大典型值线路潮流上限*2 for i 1:N_branch Constraints [Constraints, P_flow(i) M * z_br(i); P_flow(i) -M * z_br(i); P_flow(i) - (Theta(fr_bus(i)) - Theta(to_bus(i))) / X_br(i) M * (1 - z_br(i)); P_flow(i) - (Theta(fr_bus(i)) - Theta(to_bus(i))) / X_br(i) -M * (1 - z_br(i))]; end % 2. 支路容量约束 % 线路潮流不能超过容量上限同时必须考虑支路状态 for i 1:N_branch Constraints [Constraints, -Capacity(i) * z_br(i) P_flow(i) Capacity(i) * z_br(i)]; end % 3. 联络开关潮流约束 % 联络开关闭合时潮流可流经断开时流经功率为0 for t 1:N_tie tie_idx tie_lines(t); % 联络开关对应的实际线路编号 Constraints [Constraints, -Capacity(tie_idx) * x_switch(t) P_tie(t) Capacity(tie_idx) * x_switch(t)]; end % 4. 节点功率平衡约束 for i 1:N_bus % 节点注入功率 从其他节点流入的功率 - 流向其他节点的功率 节点负荷 - 切负荷 inflow sum(P_flow(find(to_bus i))) sum(P_tie(find(tie_to i))); outflow sum(P_flow(find(fr_bus i))) sum(P_tie(find(tie_fr i))); Constraints [Constraints, P_gen(i) inflow - outflow Load(i) - P_shed(i)]; end % 5. 切负荷约束 for i 1:N_bus Constraints [Constraints, 0 P_shed(i) Load(i)]; % 切负荷量不能超过该节点负荷 end % 6. DG出力约束 % 如果节点有分布式电源DG出力不能超过其容量上限 for i 1:N_dg Constraints [Constraints, 0 P_gen(DG_bus(i)) DG_capacity(i)]; end % 7. 电源节点功率约束 % 变电站根节点松弛节点被视为电源其出力受上级电网注入能力限制 Constraints [Constraints, 0 P_gen(1) Substation_capacity];目标函数和求解代码%% 目标函数最小化失负荷电量还可扩展加DG运行成本 Objective sum(Load .* P_shed / Load_total * T_outage) ... sum(DG_cost .* P_gen(DG_bus)); % 第二项是DG发电成本可选 %% 求解 ops sdpsettings(solver, gurobi, verbose, 0, showprogress, 0); solution optimize(Constraints, Objective, ops); %% 输出结果 if solution.problem 0 P_shed_value value(P_shed); Total_shed sum(P_shed_value); fprintf(故障线路 %d 下的最优切负荷量%.2f kW\n, fault_line, Total_shed); else disp(模型求解失败检查约束设置); disp(solution.info); end这个代码框架是完整跑得起来的只要你把数据矩阵填好。我特别强调几个容易出错的地方第一大M法的M值选取。大M法用于处理二进制变量带来的逻辑约束。M如果选小了会把本来可行的解误判为不可行M选太大会导致数值稳定性问题。我的经验是取该支路潮流上限的2~3倍。有些教程喜欢设M1e6这在33节点小系统还能跑一到几百节点系统就可能出现数值病态。第二流入流出的符号方向。功率平衡约束里节点的流入是“从其他节点流入本节点之和”流出是“本节点流向其他节点之和”方向搞反了算出来的全是负潮流。第三根节点的处理。第1号节点通常是变电站母线是松弛节点它的注入功率代表从上级电网获取的功率。它也是整个网络的电源节点其注入功率不能简单设成0否则会违反功率平衡。3.3 故障场景枚举与循环求解别做无意义的重复劳动故障场景的枚举方式直接影响计算效率。我在实际项目中总结出一套经验不需要枚举所有N-2故障。配电网N-2故障组合数量是N(N-1)/233节点系统是528个组合还好说要是1000节点的系统就是近50万次求解算到天荒地老。实际工程中N-1故障已经覆盖了绝大多数高风险场景N-2及以上组合风险极低对系统级指标的影响可以忽略。用“故障事件剖析法”拆分场景。一个故障事件不只包含“某条线路断开”这一个状态还包含“故障隔离”、“负荷转供”、“修复”等多个阶段。不同阶段系统的拓扑和运行状态不同。严格来说每个阶段都对应一个不同的优化模型。不过为了工程简化大多数项目把故障后的稳定运行状态作为一个静态场景来处理即假设故障隔离完成、转供策略已执行系统进入一个新的稳定运行点。这个假设在绝大多数项目中是可接受的。场景循环建议用并行计算。每个故障场景的求解相互独立天然适合并行计算。Matlab的parfor可以直接替代for来跑场景循环。我在33节点系统上测试过8核并行大约能获得5~6倍的加速比。如果场景数量几十上百这个加速非常可观。%% 场景循环求解并行版 fault_lines 1:N_branch; % 枚举所有N-1故障 P_shed_all zeros(N_bus, length(fault_lines)); parfor f 1:length(fault_lines) fault fault_lines(f); % 构建当前场景的模型并求解 P_shed_value solve_scenario(fault, system_data); P_shed_all(:, f) P_shed_value; end注意parfor里的子函数solve_scenario必须能独立运行不能依赖循环外部的变量状态变化。所以我在实践中会把所有系统数据预先打包成结构体传入。3.4 结果汇总与指标计算从优化结果到标准指标有所有故障场景的切负荷结果之后最后一步就是把它们汇总成标准可靠性指标。这一步看似简单但实际操作中很讲究。我给出指标计算的Matlab代码并逐行注释逻辑%% 输入准备 % lambda(i)第i条线路的年故障率次/年 % r(i)第i条线路的平均修复时间小时/次 % N_user(i)第i个负荷节点的用户数 % P_shed_all(i, f)第f个故障场景下第i个节点的切负荷量kW % fault_lines(f)第f个场景的故障线路编号 %% 指标初始化 N_scene length(fault_lines); ENS_total 0; % 系统总失负荷电量 SAIFI_num 0; % SAIFI分子累加用户停电次数 SAIDI_num 0; % SAIDI分子累加用户停电小时数 Total_users sum(N_user); % 总用户数 %% 累加计算 for f 1:N_scene line fault_lines(f); lambda_f lambda(line); % 该故障场景的年发生频率 r_f r(line); % 该故障场景的预计修复时间 % 失负荷电量 切负荷功率 × 修复时间 ENS_scene sum(P_shed_all(:, f)) / 1000 * r_f * lambda_f; % kWh → MWh ENS_total ENS_total ENS_scene; % 用户停电次数 受影响用户数 × 年频率 affected_users sum(N_user(P_shed_all(:, f) 0)); % 切负荷量0说明该节点停电 SAIFI_num SAIFI_num affected_users * lambda_f; % 用户停电持续时间 受影响用户数 × 年频率 × 修复时间 SAIDI_num SAIDI_num affected_users * lambda_f * r_f; end %% 输出指标 SAIFI SAIFI_num / Total_users; % 次/(用户·年) SAIDI SAIDI_num / Total_users; % 小时/(用户·年) CAIDI SAIDI_num / SAIFI_num; % 小时/次 ASAI (8760 - SAIDI) / 8760; % % EENS ENS_total; % MWh/年这段代码有几点细节值得强调故障率怎么用很多人第一步就错了。系统年平均故障率λ是“次/年”它已经包含了故障发生的概率信息。一个场景年发生频率该线路的年故障率。如果你枚举的是N-2故障那场景频率是两条线路故障率的乘积再乘以一个系数跟叠加时长有关这里需要更细致的概率计算不能简单地做乘法。“受影响用户”的判定条件。我用的是P_shed_all(i, f) 0判断节点是否停电。这个判断隐含一个假设切负荷量为0的节点完全恢复供电切负荷量大于0的节点里面所有用户都停电。实际情况可能是切负荷20%这20%的用户停电了剩下80%的用户没停电。如果需要更精细你得在优化模型里加二元变量让切负荷变成离散决策要么全切要么全不切。两种做法的结果差异在你写论文时需要明确说明采用了哪种假设。ENS和EENS的单位换算。如果你用MATLAB算出来的切负荷量单位是kW乘以修复时间小时后得到的是kWh再除以1000才是MWh。这个换算别偷懒很多人在这一步栽跟头结果差了好几个数量级还不自知。4. 常见问题与排查技巧实录4.1 求解器报“infeasible不可行”怎么办“模型不可行”是我碰到最高频的问题新手几乎必踩。常见原因和排查思路如下切负荷量上限设得太紧。如果切负荷上限是0 P_shed(i) Load(i)看起来没什么问题但你可能忘了当某个节点既没有发电也没有外来功率输入时想要完全满足功率平衡切负荷量必须等于该节点的负荷。你的上限是Load(i)下限是0模型可以切到Load(i)没问题。但如果上限设成了P_shed(i) 0.8 * Load(i)比如为了模拟“只能切80%负荷”的政策限制一旦系统没有足够的电源来补足剩下的20%模型就直接不可行了。故障线路的强制断开约束和功率平衡约束冲突。比如某条线路故障后下游节点没有任何备用电源也没有联络线可以转供优化模型唯一可行的方案就是切掉这个节点的全部负荷。如果你在约束里写了“所有负荷必须满足某个最低比例”或者你漏掉了切负荷变量的非负性约束都会导致不可行。排查方法论先注释掉一半约束看模型是否恢复可行一点点加回来找到导致不可行的“元凶约束”。这是最笨但最有效的方法。YALMIP有个工具叫yalmiptest可以帮你快速定位模型问题。4.2 求解速度慢到无法忍受怎么加速优化模型配电网可靠性评估的常见困境是系统规模一大N-1故障场景几十上百个每个场景一个MILP求解几秒加起来就是几分钟甚至几十分钟。我实测过IEEE 123节点系统全枚举N-1约120个场景每个场景Gurobi求解平均3秒串行就是6分钟起步如果用parfor并行8核大概能压到1分钟出头。这个速度做研究勉强可以接受做在线计算完全不行。如果还想进一步加速有几条路第一删减低风险场景。并非所有线路故障对可靠性指标的影响都一样大。主干线路故障可能引起大范围用户停电分支线路故障可能只影响一两个用户。先跑一遍所有场景记录每个场景的切负荷量设定一个阈值——切负荷量低于某个百分比比如最大值的0.1%的场景下次计算直接跳过。这样可能以微小精度损失换取数倍速度提升。第二热启动warm start。把上一个场景的最优解作为下一个场景的初始解。对相邻故障场景来说最优解往往很接近热启动能显著减少分支定界搜索量。在Gurobi里通过把上一场景的变量值赋给当前场景的x0字段实现。第三松弛预求解。每个故障场景先做LP松弛求解如果LP松弛解已经是整数可行解就无需进入MILP分支定界过程。Gurobi默认会自动做这个操作但要确保你不会不小心关掉了预求解功能。第四缩小故障持续时间精度。修复时间r的精度从小时级精确到分钟级就够用了没必要精确到秒。因为配电网修复时间本身是统计量误差也许就有几十分钟精度高反而造成“虚假的精确”。4.3 结果和论文对不上我的排查顺序跑完程序后发现指标和顶刊论文差了十万八千里。别慌按下面的顺序逐个排查先查数据源。同样叫“IEEE 33节点系统”不同版本的负荷数据、线路参数可能不同。我看过不止一篇论文用的是修改版数据比如增加负荷、改变线路型号计算出来的指标根本没法直接对比。看论文时留意它在附录里给出的详细参数表把自己代码里的数据逐项核对一遍。再查指标口径。有些论文的SAIFI只统计“持续性故障停电”不包括瞬时停电有些论文把检修停电也纳入统计有些论文的EENS只算故障场景不算检修场景。口径不同指标可以差出20%~50%。务必搞清楚对方计算指标时的统计边界。最后查模型假设。论文里用了什么恢复策略是全部故障后统一优化恢复还是只允许通过联络开关恢复一部分DG是否允许孤岛运行这些假设的区别直接决定了可靠性指标的大小而且影响很大。比如允许孤岛运行的模型算出来的EENS可能比不允许孤岛运行低一个量级。4.4 求解结果里切负荷量为负或者其他物理上说不通的结果切负荷量为负说明模型把“切负荷”当成了一种可以“注入功率”的手段。这通常是因为你写功率平衡约束时把P_shed的符号方向搞反了。检查一下功率平衡方程是“负荷注入流入-流出”切负荷是减少负荷所以它应该在等式左边带负号也就是Load(i) - P_shed(i)。如果你写成了Load(i) P_shed(i)模型就会通过负的P_shed来“补”功率。另一种物理上说不通的结果是故障线路下游节点有负荷需求但系统说这部分失负荷量是0同时所有的注入功率也完全满足——但线路容量约束又显示某条线路已经跑到上限的200%。这种情况是线路容量约束写错了很可能是大M法的M值设得太大比如1e6导致容量约束被松弛掉了。把M改小让它和线路容量一个量级问题就能解决。4.5 故障枚举和指标计算的两处小坑最后记录两个我在实际项目中踩过的“不起眼但非常致命”的坑。第一个是联络开关的转供能力限制。配电网故障恢复最常见的手段是通过联络开关从其他馈线转供。但联络线的转供能力不是无限的它受两个因素限制一是联络线本身的额定容量二是对侧馈线的剩余容量。如果你只约束了前者、忽略了后者工程上可能得出“看似可行、实际完全无法实施”的恢复方案。严格的做法是把对侧馈线上所有线路的剩余容量都纳入约束检查。第二个是修复时间的差异化。电缆线路和架空线路的故障修复时间差异极大电缆故障定位困难修复时间动辄数小时到数十小时架空线路相对好修典型修复时间3~5小时。我在处理33节点系统时见过一些人把整条系统的所有r都设成同一个值比如统一的4小时结果SAIDI算出来完全不对。合理的做法是不同线路类型设置不同的修复时间甚至结合故障类型永久性故障还是瞬时性故障设置不同的持续时间权重。写到这里完整流程已经走通一遍建数据、枚举故障、建优化模型、求解、汇总指标、排查异常。这个基于优化模型的配电网可靠性评估项目从原理上讲透了为什么这个方案比传统解析法更适合带分布式电源的主动配电网从操作上讲透了从原始数据到SAIFI/SAIDI/ASAI/EENS全套指标的Matlab实现从实战上讲透了求解器选型、大M法陷阱、场景并行计算这些普通教程不会细说的经验。我在实际跑这个项目时最大的体会是可靠性评估这件事最大的难点不是数学建模本身而是对“模型假设”的把握。你用优化模型模拟故障恢复本质上是在替运行人员做决策——这个决策合不合理取决于你的目标函数和约束条件有没有准确反映真实的运行规则。所以每跑完一个算例我都会问自己一句这个结果在物理上、在工程习惯上说得通吗如果说不通多半不是求解器的问题而是我自己的模型设计出了问题。把心态摆正一步步排查你会很快从“代码能跑”进化到“模型靠谱”的层次。