这段时间周围好几个做配电网方向的朋友都在聊同一个话题：主动配电网故障恢复怎么用两阶段鲁棒优化来做，Matlab代码该怎么组织。说实话，这个题目听起来确实有些门槛，但如果你真的把主动配电网、两阶段鲁棒、恢复策略、Matlab代码这几个关键词拆开看，会发现它其实是配电网优化里一条非常清晰的技术路线。我这段时间把整套流程从模型到代码完整跑了一遍，把思路、模型、求解框架、代码结构、踩坑记录都整理出来，希望能给正在做类似课题的同行一点参考。
1. 问题背景：为什么配电网恢复需要“两阶段鲁棒”
1.1 确定性恢复方案为什么不够用
传统配电网故障恢复通常是在故障发生后，通过网络重构、分布式电源调度、储能充放电等手段，把失电负荷尽量恢复起来。这里面的“传统方案”有一个共同的隐含前提：负荷大小、DG出力、储能状态都被假设成已知的固定值，然后在确定性模型里求解一个单阶段优化问题。
这个假设在早期配电网里勉强能用，因为当时配网里分布式电源占比很低，负荷预测精度也基本可以接受。但现在情况完全不同：光伏出力随天气剧烈波动，风电出力更难预测，负荷侧还有电动汽车充电这种不确定性很大的用户行为。如果还是用确定性模型去定恢复方案，一旦实际场景偏离预测值，恢复策略可能直接失效。比如优化结果告诉你某条联络开关需要闭合、某个储能需要放电来支撑电压，但实际DG出力远低于预测值，开关状态已经固化又没法临时调整，结果就是电压越限甚至恢复失败。
所以，恢复方案不能只针对单一预测场景，而要在最不利的不确定性场景下仍然保证安全性，同时尽量提高恢复效果。这正是鲁棒优化介入的原因。
1.2 两阶段鲁棒优化到底在优化什么
两阶段鲁棒优化的核心思想是“先决策、后应对”。放到配电网恢复场景里，第一阶段决策是那些一旦确定就很难快速改变的变量，比如分段开关和联络开关的状态、储能是否投入运行；第二阶段决策是在第一阶段决策确定的开关拓扑下，面对不确定性参数的最坏实现，运行变量如何调整，比如DG出力、储能充放电功率、弃负荷量等。
这个过程可以理解为：你作为调度员，先根据当前掌握的信息把能定的事定下来，然后假设老天爷（不确定性）会尽最大努力给你添麻烦，在最坏情况下你仍然能通过运行调整把损失控制在最小。这种“先定预案、再对抗最坏场景”的思路，比单纯做确定性优化稳健得多。
1.3 这个项目适合谁，能学到什么
如果你正在研究主动配电网故障恢复、韧性评估，或者在做相关课题需要把两阶段鲁棒优化落地成代码，那么这个项目会非常对口。即使你不是电力系统出身，只要对优化建模和Matlab有一定基础，通过这份代码也能把两阶段鲁棒优化的建模思路、C&CG求解框架、Yalmip建模方式完整串起来。
看完这个项目，你能收获三样东西：主动配电网恢复问题的标准数学建模方法；两阶段鲁棒优化与C&CG算法的完整求解框架；一套可以直接改参数、换数据、做扩展实验的Matlab代码骨架。
2. 数学建模：从物理问题到标准优化模型
2.1 网络潮流与设备的数学抽象
配电网恢复问题首先要把物理网络抽象成数学模型。潮流方程通常采用DistFlow模型，它通过有功、无功、电压幅值、支路电流之间的递推关系描述网络状态。严格来说DistFlow包含非线性项，直接求解困难，好在工程上普遍采用二阶锥松弛将其转化为凸约束：$P_{ij}^2 + Q_{ij}^2 \leq v_i l_{ij}$ 可写成二阶锥形式 $||[2P_{ij}; 2Q_{ij}; v_i - l_{ij}]||2 \leq v_i + l{ij}$。在Matlab里用Yalmip的cone()函数就能轻松表示。
设备建模方面，分布式电源通常给一个出力约束区间，储能需要建模荷电状态和充放电功率之间的关系，开关状态则是0-1变量。恢复问题最核心的需求是“辐射状拓扑约束”，即配电网在重构后依然保持树状结构。这个约束可以通过生成树计数模型近似实现，即支路数等于节点数减去电源数，同时用方向约束保证连通性。
2.2 目标函数与约束体系
恢复问题的目标函数通常包含三部分：最大化恢复的负荷量（可以按负荷重要性加权重）、最小化网损、最小化弃风弃光量。在实际项目中，最常见的做法是把三者统一成加权失负荷量最小化，再辅以网损惩罚项。
约束体系可以分成四块：潮流约束，即节点有功/无功平衡；设备运行约束，包括DG出力和储能充放电功率约束、储能荷电状态递推关系；运行安全约束，包括节点电压上下限和支路电流上限；网络拓扑约束，即辐射状结构约束。
2.3 不确定性集合：盒式加预算约束
鲁棒优化需要显式建模不确定性参数的可能范围。工程上最常用的是盒式不确定集合，比如光伏出力预测误差在 $[- \Delta, \Delta]$ 范围内，负荷预测误差在 $[-\hat{d}, \hat{d}]$ 范围内。如果所有不确定参数都同时取最坏值，方案会过度保守，所以通常会加入一个预算约束（Budget of Uncertainty），限制同一时刻出现偏差的不确定参数总数不超过某个整数 $\Gamma$。
这个 $\Gamma$ 的设计很有讲究：$\Gamma$ 越小，模型越乐观；$\Gamma$ 越大，模型越保守。具体取多少取决于调度员对风险的态度和实际运行经验。我在项目中一般取 $\Gamma = \lceil \alpha \cdot n_{unc} \rceil$，其中 $\alpha$ 在0.3到0.7之间。
2.4 两阶段模型的紧凑形式
把上面所有元素组合起来，可以写成标准的两阶段鲁棒优化形式：
外层（第一阶段）最小化决策成本加最坏场景下的运行成本；内层引入不确定性变量，在给定第一阶段决策下最大化最坏情况，再内层继续最小化运行成本。这个三层结构就是典型的两阶段鲁棒模型：$\min_x c^T x + \max_{u \in U} \min_{y \in F(x,u)} d^T y$，约束包含第一阶段决策变量的可行域以及第二阶段运行约束。
3. 求解算法：C&CG列与约束生成
3.1 为什么选择C&CG而不是Benders分解
两阶段鲁棒优化常见的求解方法有两种：Benders分解和C&CG（Column and Constraint Generation，列与约束生成）。早期文献多用Benders分解，通过主问题与子问题之间传递对偶割平面迭代求解。但Benders在处理包含二阶锥约束的配电网问题时，割平面质量往往不稳定，收敛速度也不理想。
C&CG的思路完全不同：它不仅往主问题添加割平面，还把子问题识别到的最坏场景作为新的列变量加入主问题，让主问题的决策空间随着迭代不断扩展。这种做法在离散决策变量多、约束维度高的问题上收敛速度明显更快，而且对求解器的友好度更高。我在项目里首选的也是C&CG。
3.2 主问题与子问题的交互逻辑
C&CG算法大致分成四步。第一步，初始化：设定迭代次数 $k=1$，下界 $LB=-\infty$，上界 $UB=+\infty$，收敛精度 $\epsilon$。第二步，求解主问题：主问题包含第一阶段决策变量 $x$、辅助变量 $\eta$ 以及已识别出的有限个最坏场景对应的运行变量和约束，求解得到 $x^{k}$ 和目标值 $\eta^{k}$，更新 $LB = \eta^{k}$。第三步，求解子问题：将 $x^{k}$ 固定代入子问题，求解最坏场景 $u^{k}$ 和对应的运行成本 $Q(x^{k})$，更新 $UB = \min{UB, c^T x^{k} + Q(x^{k})}$。第四步，判断收敛：若 $UB - LB < \epsilon$，停止迭代；否则，把新识别的最坏场景 $u^{k}$ 对应的运行变量和约束加入主问题，更新 $k = k+1$，回到第二步。
值得注意的是，主问题每迭代一次就会多一组场景约束，所以主问题的规模会逐渐增大。好在配电网恢复问题的场景约束数目通常在几十次迭代内就能收敛，不会造成不可接受的求解负担。
3.3 收敛判据与有限迭代性质
C&CG的收敛性在理论上有保障：当不确定性集合是有限多面体时，C&CG能在有限次迭代内收敛到全局最优解。实际操作中收敛精度 $\epsilon$ 一般取 $10^{-3}$ 或者 $10^{-4}$，再结合每轮主问题和子问题的目标值做一个相对误差判断。
我在代码里还额外加了一个调试用的小技巧：打印每一轮的LB和UB变化曲线，如果发现UB在下行过程中突然反弹，说明主问题新增场景约束之后第一阶段决策发生了跳变，这时候要重点检查子问题对偶推导是否正确。
4. Matlab代码实现：从建模到跑通全程拆解
4.1 环境准备与数据组织
代码依赖环境是Matlab外加上Yalmip工具箱，再配一个混合整数求解器，CPLEX或Gurobi都可以。Yalmip的好处在于建模语法简洁，不需要自己处理复杂的矩阵拼接，对于配电网这种约束多、变量维度高的问题尤其省心。
数据组织方面，我建议把节点、线路、DG、储能参数拆成结构体或者表格管理。比如busData保存节点编号、有功/无功负荷；branchData保存支路首末端节点、电阻、电抗；dgData保存DG接入节点、出力上下限；essData保存储能接入节点、容量、充放电效率。代码里用索引访问，后续换算例的时候只改数据文件，不动主逻辑。
4.2 主问题代码实现
主问题的核心逻辑是定义第一阶段决策变量、辅助变量、目标函数以及所有已识别场景对应的约束。下面是一段节选的核心代码：
%% 主问题 MP % 第一阶段决策变量 z_sw = binvar(n_branch, 1); % 支路开关状态，1为闭合 z_ess_on = binvar(n_ess, 1); % 储能启用状态 % 辅助变量 eta = sdpvar(1, 1); % 最坏场景运行成本的替代变量 P_rec = sdpvar(n_bus, n_scen); % 每个场景下的负荷恢复量 % 其他运行变量：P_ij、Q_ij、v_i、i_ij、P_dg、P_ess_ch、P_ess_dis 等 Constraints = []; % 拓扑辐射状约束：闭合支路数 = 节点数 - 根节点数 Constraints = [Constraints, sum(z_sw) == n_bus - 1]; % 支路开关与潮流变量耦合 for k = 1:n_branch Constraints = [Constraints, ... % 潮流变量与 z_sw 的 Big-M 约束 P_ij(k, :) <= M * z_sw(k), ... -M * z_sw(k) <= P_ij(k, :), ... Q_ij(k, :) <= M * z_sw(k), ... -M * z_sw(k) <= Q_ij(k, :)]; end % 目标函数：优先恢复重要负荷，同时惩罚网损 Objective = sum(w_load .* (P_load - sum(P_rec, 2))) + lambda_loss * sum(P_loss) + eta; % 调用求解器 ops = sdpsettings('solver', 'gurobi', 'verbose', 0, 'mipgap', 1e-4); optimize(Constraints, Objective, ops);这段代码里M是大数，注意不要取得太大，否则会引入数值问题，一般取负荷总量的5到10倍即可。目标函数中eta实际上代替了在最坏场景下的运行成本，它会在主问题中通过场景约束被逐步收紧。
4.3 子问题与对偶处理
子问题的结构是给定第一阶段决策后求最坏场景下的运行成本。直接求解max-min问题很困难，常见的做法是用强对偶把内层min转化为max，和外层max合并成一个单层max问题。在配电网恢复中，如果潮流约束是线性化的DistFlow，那么内层是LP，对偶推导并不复杂；如果采用二阶锥松弛，对偶推导会涉及锥对偶，复杂度上了一个台阶。
我在项目里为了兼顾精度和可操作性，采用了“线性DistFlow建模+子问题对偶”的方案。核心代码如下：
%% 子问题 SP：给定 z_sw 后求最坏场景 % 内层对偶变量 lambda_eq = sdpvar(n_eq, 1); % 等式约束对偶 lambda_ineq = sdpvar(n_ineq, 1); % 不等式约束对偶 % 对偶目标：d^T y -> 对偶表达 Objective_dual = lambda_ineq' * (b_ineq + B_u_ineq * u) + lambda_eq' * d_eq; % 对偶可行域约束 Constraints_dual = [A_dual * [lambda_eq; lambda_ineq] >= c_dual, ... lambda_ineq >= 0]; % 外层再对 u 最大化 Constraints_u = [u_min <= u <= u_max, sum(u .* w_u) <= Gamma]; % 合并后的子问题 optimize([Constraints_dual, Constraints_u], -Objective_dual, ops); Q_xk = value(Objective_dual);这里有几个地方特别容易踩坑：一是对偶变量与原始约束方向的匹配关系，写反了会导致对偶目标符号不对；二是不等式约束对偶变量必须非负，等式约束对偶变量自由；三是在配电网模型中，如果含有二阶锥约束，对偶问题里会出现对应的锥约束，直接写成线性约束会出错。
如果实在不想手动推对偶，还有一个思路：在Yalmip里直接把第二阶段原问题对应的KKT条件加上互补松弛约束，用二进制变量和大M法线性化互补条件。这种做法的缺点是引入了大量二进制变量，求解效率会下降，但在验证模型正确性的时候很有用。我用它来交叉验证对偶推导是否正确。
4.4 主循环与收敛输出
把主问题和子问题串起来的是C&CG主循环。代码框架如下：
%% C&CG 主循环 LB = -1e6; UB = 1e6; k = 1; eps = 1e-3; while UB - LB > eps && k <= max_iter % 解主问题 optimize(MP_Constraints, MP_Objective, ops); LB = value(eta); % 固定第一阶段决策，解子问题 x_fixed = value([z_sw; z_ess_on]); optimize(SP_Constraints, SP_Objective, ops); Q_k = value(Objective_dual); UB = min(UB, value(c' * x_fixed) + Q_k); % 添加新的场景约束到主问题 u_k = value(u); add_scenario(MP_Constraints, u_k); fprintf('Iter=%d, LB=%.4f, UB=%.4f, gap=%.4f\n', k, LB, UB, UB-LB); k = k + 1; end主循环里有一个容易忽略的点：每轮子问题求解后要立即保存当前的最坏场景u_k，因为后续求解主问题时，这个场景会作为已知参数被固化到新增约束中。如果值提取写错了位置，会把其他迭代轮次的场景混进来，导致主问题约束混乱。
5. 案例测试与结果分析
5.1 测试场景设计
我选用了标准的IEEE 33节点配电网系统作为测试算例：33个节点、32条分段支路、5条联络开关支路，接入了两个光伏电站和一个储能系统。故障场景设置为某条主干线停运，导致下游区域失电，系统需要通过重构和DG调度恢复负荷。
不确定性参数设置了两类：光伏出力在预测值上下浮动20%到30%，负荷在预测值上下浮动10%。预算约束参数 $\Gamma$ 从0取到10逐一测试，观察恢复方案和计算时间的变化。
5.2 鲁棒方案与确定性方案对比
我在同一故障场景下分别跑确定性恢复模型和两阶段鲁棒恢复模型，然后把两种方案放到最坏场景里进行回代检验。结果差异非常明显：确定性方案在最坏场景下出现了严重的电压越限，多个节点电压低于0.90pu，同时有约15%的已恢复负荷被迫再次切除；鲁棒方案在最坏场景下依然能保持节点电压在0.95pu以上，恢复负荷比例也稳定在90%以上。
| 指标 | 确定性方案（最坏场景回代） | 鲁棒方案（最坏场景回代） |
|---|---|---|
| 恢复负荷比例 | 82.4% | 93.7% |
| 最低节点电压 | 0.883 pu | 0.951 pu |
| 网损（相对值） | 1.000 | 1.142 |
| 开关动作次数 | 7 | 9 |
从结果看，鲁棒方案牺牲了一定的网损经济性，换来了更高的恢复可靠性和电压安全裕度。这个权衡在实际工程中是值得的，因为故障恢复阶段的首要目标不是经济性，而是供电可靠性。
5.3 敏感性分析与计算效率
$\Gamma$ 从0增大到10的过程中，目标函数值逐步变差（因为方案越来越保守），但最坏场景下的安全性指标逐步改善。当 $\Gamma$ 超过8之后，目标值改善幅度已经非常小，说明过度保守并没有带来额外收益，实际中可以据此选择一个合适的预算参数。
计算效率方面，33节点系统规模下，C&CG算法平均迭代次数在8到14次之间收敛，单轮主问题和子问题求解时间在0.5到3秒之间，总耗时通常不超过30秒。如果换到更大的系统，比如123节点系统，迭代次数会增加到20次以上，此时需要关注主问题规模的增长速度。
6. 常见问题与排查经验实录
6.1 收敛慢、上下界震荡怎么处理
C&CG最常见的现象是前几轮LB和UB快速靠近，但到了最后1%的gap时卡住不动。处理方法无非三种：把MIP gap从默认值调小一点，像我在代码里设置的mipgap=1e-4；检查子问题是否真的是最优解，必要时给求解器设置更严格的容差；如果UB反复震荡，优先怀疑对偶问题推导有问题，尤其是对偶变量符号和约束方向。
另外要注意主问题每次新增场景约束后，之前保存的所有场景变量都要保持固定，否则求解器会调整历史场景对应的运行变量来迎合当前目标，导致主问题退化成次优解。
6.2 对偶推导容易出错的位置
对偶推导的正确性决定整个C&CG算法是否可行。最容易出错的是：不等式约束的对偶变量非负方向搞反；等式约束对偶变量没有设成自由变量；目标函数优化方向改变后没有做负号转换；内层min的约束中，不确定性参数在对偶中出现的位置和符号错误。
我的建议是遇到对偶结果可疑时，先把不确定性固定成某个已知场景，和单阶段确定性模型的结果对比。如果两者不一致，基本可以断定是对偶表达式的问题。
6.3 求解器数值问题与参数调整
配电网模型变量数量本身不大，但二阶锥约束和Big-M约束都对数值敏感。Big-M取太大，整数变量的线性松弛会变得非常松散，MIP求解效率直线下降；二阶锥约束的gaptolerance也需要设置合理，默认的1e-6有时候过严，反而导致无解，放松到1e-5或1e-4更容易收敛。
以我的经验，Big-M取值取节点负荷总和的5倍到10倍就够了，不要为了“绝对安全”取到1e6这种量级。
6.4 经验速查表
| 问题 | 原因 | 解决方案 |
|---|---|---|
| 主问题无解 | 辐射状约束太紧或Big-M不匹配 | 检查闭合支路数约束，增大M |
| 子问题无解 | 对偶方向错误或可行域过紧 | 用KKT法交叉验证 |
| 上下界震荡 | 新增场景后历史变量被重算 | 固定所有历史场景变量 |
| 迭代次数过多 | MIP gap过大或预算约束不合理 | 收紧mipgap，调小Gamma |
| 电压越限 | 不确定性集合描述不准确 | 检查Delta和Gamma设置 |
这些坑我基本都踩过一轮，写在这里帮大家节省一些调代码的时间。
结尾小经验
整个项目跑下来，我最大的体会是：两阶段鲁棒恢复这个方向，数学模型本身并不难理解，真正的门槛在于把抽象的三层优化结构干净地落到代码里。C&CG框架的代码骨架不到两百行，但每一步推导、每个变量索引、每个对偶符号都可能在细节上出问题。我的建议是先从一个很小的系统开始，比如把33节点简化为6节点甚至3节点系统，手动验证每个中间结果，再逐步扩展到完整算例。这样定位问题会快很多，对模型的理解也会更扎实。
