做配电网可靠性评估这些年，我最常被问的一句话就是：“可靠性评估到底该用解析法、模拟法，还是构建优化模型？”尤其是刚从论文里看到“基于优化模型”这个提法时，很多同行会愣一下——优化不是用来做规划、做调度、做经济性分配的吗，怎么还能拿来评估可靠性？这篇文章就围绕一个顶刊复现项目来讲透这件事。
我要拆解的这个项目，标题已经写得很明白：“基于优化模型的配电网可靠性评估研究（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这类建模工具把优化模型翻译成求解器能识别的标准格式，然后调用求解器去解。
我从实际使用经验出发，给一套选型建议：
| 求解器 | 适用问题类型 | 许可证 | 我的使用评价 |
|---|---|---|---|
| Gurobi | MILP、MIQP、LP、QP | 学术免费，商用收费 | 工业级标杆，速度极快，配电网可靠性评估首选 |
| CPLEX | MILP、MIQP、LP、QP | 学术免费，商用收费 | 老牌求解器，稳定可靠，和Gurobi水平相当 |
| MOSEK | 凸优化、锥优化、MIQP | 学术免费，商用收费 | 处理SOCP很强，适合用于交流潮流凸松弛 |
| GLPK | LP、MILP | 开源免费 | 小规模可以，速度慢，大规模慎用 |
| MATLAB自带intlinprog | MILP | 随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 用Matlab+YALMIP搭建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倍。有些教程喜欢设M=1e6，这在33节点小系统还能跑，一到几百节点系统就可能出现数值病态。
第二，流入流出的符号方向。功率平衡约束里，节点的流入是“从其他节点流入本节点之和”，流出是“本节点流向其他节点之和”，方向搞反了，算出来的全是负潮流。
第三，根节点的处理。第1号节点通常是变电站母线，是松弛节点，它的注入功率代表从上级电网获取的功率。它也是整个网络的电源节点，其注入功率不能简单设成0，否则会违反功率平衡。
3.3 故障场景枚举与循环求解：别做无意义的重复劳动
故障场景的枚举方式直接影响计算效率。我在实际项目中总结出一套经验：
不需要枚举所有N-2故障。配电网N-2故障组合数量是N(N-1)/2，33节点系统是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法陷阱、场景并行计算这些普通教程不会细说的经验。
我在实际跑这个项目时最大的体会是，可靠性评估这件事，最大的难点不是数学建模本身，而是对“模型假设”的把握。你用优化模型模拟故障恢复，本质上是在替运行人员做决策——这个决策合不合理，取决于你的目标函数和约束条件有没有准确反映真实的运行规则。所以每跑完一个算例，我都会问自己一句：这个结果在物理上、在工程习惯上说得通吗？如果说不通，多半不是求解器的问题，而是我自己的模型设计出了问题。把心态摆正，一步步排查，你会很快从“代码能跑”进化到“模型靠谱”的层次。
