这段时间周围好几个做配电网方向的朋友都在聊同一个话题主动配电网故障恢复怎么用两阶段鲁棒优化来做Matlab代码该怎么组织。说实话这个题目听起来确实有些门槛但如果你真的把主动配电网、两阶段鲁棒、恢复策略、Matlab代码这几个关键词拆开看会发现它其实是配电网优化里一条非常清晰的技术路线。我这段时间把整套流程从模型到代码完整跑了一遍把思路、模型、求解框架、代码结构、踩坑记录都整理出来希望能给正在做类似课题的同行一点参考。1. 问题背景为什么配电网恢复需要“两阶段鲁棒”1.1 确定性恢复方案为什么不够用传统配电网故障恢复通常是在故障发生后通过网络重构、分布式电源调度、储能充放电等手段把失电负荷尽量恢复起来。这里面的“传统方案”有一个共同的隐含前提负荷大小、DG出力、储能状态都被假设成已知的固定值然后在确定性模型里求解一个单阶段优化问题。这个假设在早期配电网里勉强能用因为当时配网里分布式电源占比很低负荷预测精度也基本可以接受。但现在情况完全不同光伏出力随天气剧烈波动风电出力更难预测负荷侧还有电动汽车充电这种不确定性很大的用户行为。如果还是用确定性模型去定恢复方案一旦实际场景偏离预测值恢复策略可能直接失效。比如优化结果告诉你某条联络开关需要闭合、某个储能需要放电来支撑电压但实际DG出力远低于预测值开关状态已经固化又没法临时调整结果就是电压越限甚至恢复失败。所以恢复方案不能只针对单一预测场景而要在最不利的不确定性场景下仍然保证安全性同时尽量提高恢复效果。这正是鲁棒优化介入的原因。1.2 两阶段鲁棒优化到底在优化什么两阶段鲁棒优化的核心思想是“先决策、后应对”。放到配电网恢复场景里第一阶段决策是那些一旦确定就很难快速改变的变量比如分段开关和联络开关的状态、储能是否投入运行第二阶段决策是在第一阶段决策确定的开关拓扑下面对不确定性参数的最坏实现运行变量如何调整比如DG出力、储能充放电功率、弃负荷量等。这个过程可以理解为你作为调度员先根据当前掌握的信息把能定的事定下来然后假设老天爷不确定性会尽最大努力给你添麻烦在最坏情况下你仍然能通过运行调整把损失控制在最小。这种“先定预案、再对抗最坏场景”的思路比单纯做确定性优化稳健得多。1.3 这个项目适合谁能学到什么如果你正在研究主动配电网故障恢复、韧性评估或者在做相关课题需要把两阶段鲁棒优化落地成代码那么这个项目会非常对口。即使你不是电力系统出身只要对优化建模和Matlab有一定基础通过这份代码也能把两阶段鲁棒优化的建模思路、CCG求解框架、Yalmip建模方式完整串起来。看完这个项目你能收获三样东西主动配电网恢复问题的标准数学建模方法两阶段鲁棒优化与CCG算法的完整求解框架一套可以直接改参数、换数据、做扩展实验的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. 求解算法CCG列与约束生成3.1 为什么选择CCG而不是Benders分解两阶段鲁棒优化常见的求解方法有两种Benders分解和CCGColumn and Constraint Generation列与约束生成。早期文献多用Benders分解通过主问题与子问题之间传递对偶割平面迭代求解。但Benders在处理包含二阶锥约束的配电网问题时割平面质量往往不稳定收敛速度也不理想。CCG的思路完全不同它不仅往主问题添加割平面还把子问题识别到的最坏场景作为新的列变量加入主问题让主问题的决策空间随着迭代不断扩展。这种做法在离散决策变量多、约束维度高的问题上收敛速度明显更快而且对求解器的友好度更高。我在项目里首选的也是CCG。3.2 主问题与子问题的交互逻辑CCG算法大致分成四步。第一步初始化设定迭代次数 $k1$下界 $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 k1$回到第二步。值得注意的是主问题每迭代一次就会多一组场景约束所以主问题的规模会逐渐增大。好在配电网恢复问题的场景约束数目通常在几十次迭代内就能收敛不会造成不可接受的求解负担。3.3 收敛判据与有限迭代性质CCG的收敛性在理论上有保障当不确定性集合是有限多面体时CCG能在有限次迭代内收敛到全局最优解。实际操作中收敛精度 $\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 主循环与收敛输出把主问题和子问题串起来的是CCG主循环。代码框架如下%% CCG 主循环 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 pu0.951 pu网损相对值1.0001.142开关动作次数79从结果看鲁棒方案牺牲了一定的网损经济性换来了更高的恢复可靠性和电压安全裕度。这个权衡在实际工程中是值得的因为故障恢复阶段的首要目标不是经济性而是供电可靠性。5.3 敏感性分析与计算效率$\Gamma$ 从0增大到10的过程中目标函数值逐步变差因为方案越来越保守但最坏场景下的安全性指标逐步改善。当 $\Gamma$ 超过8之后目标值改善幅度已经非常小说明过度保守并没有带来额外收益实际中可以据此选择一个合适的预算参数。计算效率方面33节点系统规模下CCG算法平均迭代次数在8到14次之间收敛单轮主问题和子问题求解时间在0.5到3秒之间总耗时通常不超过30秒。如果换到更大的系统比如123节点系统迭代次数会增加到20次以上此时需要关注主问题规模的增长速度。6. 常见问题与排查经验实录6.1 收敛慢、上下界震荡怎么处理CCG最常见的现象是前几轮LB和UB快速靠近但到了最后1%的gap时卡住不动。处理方法无非三种把MIP gap从默认值调小一点像我在代码里设置的mipgap1e-4检查子问题是否真的是最优解必要时给求解器设置更严格的容差如果UB反复震荡优先怀疑对偶问题推导有问题尤其是对偶变量符号和约束方向。另外要注意主问题每次新增场景约束后之前保存的所有场景变量都要保持固定否则求解器会调整历史场景对应的运行变量来迎合当前目标导致主问题退化成次优解。6.2 对偶推导容易出错的位置对偶推导的正确性决定整个CCG算法是否可行。最容易出错的是不等式约束的对偶变量非负方向搞反等式约束对偶变量没有设成自由变量目标函数优化方向改变后没有做负号转换内层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设置这些坑我基本都踩过一轮写在这里帮大家节省一些调代码的时间。结尾小经验整个项目跑下来我最大的体会是两阶段鲁棒恢复这个方向数学模型本身并不难理解真正的门槛在于把抽象的三层优化结构干净地落到代码里。CCG框架的代码骨架不到两百行但每一步推导、每个变量索引、每个对偶符号都可能在细节上出问题。我的建议是先从一个很小的系统开始比如把33节点简化为6节点甚至3节点系统手动验证每个中间结果再逐步扩展到完整算例。这样定位问题会快很多对模型的理解也会更扎实。
主动配电网两阶段鲁棒优化故障恢复：Matlab+CCG实现详解
这段时间周围好几个做配电网方向的朋友都在聊同一个话题主动配电网故障恢复怎么用两阶段鲁棒优化来做Matlab代码该怎么组织。说实话这个题目听起来确实有些门槛但如果你真的把主动配电网、两阶段鲁棒、恢复策略、Matlab代码这几个关键词拆开看会发现它其实是配电网优化里一条非常清晰的技术路线。我这段时间把整套流程从模型到代码完整跑了一遍把思路、模型、求解框架、代码结构、踩坑记录都整理出来希望能给正在做类似课题的同行一点参考。1. 问题背景为什么配电网恢复需要“两阶段鲁棒”1.1 确定性恢复方案为什么不够用传统配电网故障恢复通常是在故障发生后通过网络重构、分布式电源调度、储能充放电等手段把失电负荷尽量恢复起来。这里面的“传统方案”有一个共同的隐含前提负荷大小、DG出力、储能状态都被假设成已知的固定值然后在确定性模型里求解一个单阶段优化问题。这个假设在早期配电网里勉强能用因为当时配网里分布式电源占比很低负荷预测精度也基本可以接受。但现在情况完全不同光伏出力随天气剧烈波动风电出力更难预测负荷侧还有电动汽车充电这种不确定性很大的用户行为。如果还是用确定性模型去定恢复方案一旦实际场景偏离预测值恢复策略可能直接失效。比如优化结果告诉你某条联络开关需要闭合、某个储能需要放电来支撑电压但实际DG出力远低于预测值开关状态已经固化又没法临时调整结果就是电压越限甚至恢复失败。所以恢复方案不能只针对单一预测场景而要在最不利的不确定性场景下仍然保证安全性同时尽量提高恢复效果。这正是鲁棒优化介入的原因。1.2 两阶段鲁棒优化到底在优化什么两阶段鲁棒优化的核心思想是“先决策、后应对”。放到配电网恢复场景里第一阶段决策是那些一旦确定就很难快速改变的变量比如分段开关和联络开关的状态、储能是否投入运行第二阶段决策是在第一阶段决策确定的开关拓扑下面对不确定性参数的最坏实现运行变量如何调整比如DG出力、储能充放电功率、弃负荷量等。这个过程可以理解为你作为调度员先根据当前掌握的信息把能定的事定下来然后假设老天爷不确定性会尽最大努力给你添麻烦在最坏情况下你仍然能通过运行调整把损失控制在最小。这种“先定预案、再对抗最坏场景”的思路比单纯做确定性优化稳健得多。1.3 这个项目适合谁能学到什么如果你正在研究主动配电网故障恢复、韧性评估或者在做相关课题需要把两阶段鲁棒优化落地成代码那么这个项目会非常对口。即使你不是电力系统出身只要对优化建模和Matlab有一定基础通过这份代码也能把两阶段鲁棒优化的建模思路、CCG求解框架、Yalmip建模方式完整串起来。看完这个项目你能收获三样东西主动配电网恢复问题的标准数学建模方法两阶段鲁棒优化与CCG算法的完整求解框架一套可以直接改参数、换数据、做扩展实验的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. 求解算法CCG列与约束生成3.1 为什么选择CCG而不是Benders分解两阶段鲁棒优化常见的求解方法有两种Benders分解和CCGColumn and Constraint Generation列与约束生成。早期文献多用Benders分解通过主问题与子问题之间传递对偶割平面迭代求解。但Benders在处理包含二阶锥约束的配电网问题时割平面质量往往不稳定收敛速度也不理想。CCG的思路完全不同它不仅往主问题添加割平面还把子问题识别到的最坏场景作为新的列变量加入主问题让主问题的决策空间随着迭代不断扩展。这种做法在离散决策变量多、约束维度高的问题上收敛速度明显更快而且对求解器的友好度更高。我在项目里首选的也是CCG。3.2 主问题与子问题的交互逻辑CCG算法大致分成四步。第一步初始化设定迭代次数 $k1$下界 $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 k1$回到第二步。值得注意的是主问题每迭代一次就会多一组场景约束所以主问题的规模会逐渐增大。好在配电网恢复问题的场景约束数目通常在几十次迭代内就能收敛不会造成不可接受的求解负担。3.3 收敛判据与有限迭代性质CCG的收敛性在理论上有保障当不确定性集合是有限多面体时CCG能在有限次迭代内收敛到全局最优解。实际操作中收敛精度 $\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 主循环与收敛输出把主问题和子问题串起来的是CCG主循环。代码框架如下%% CCG 主循环 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 pu0.951 pu网损相对值1.0001.142开关动作次数79从结果看鲁棒方案牺牲了一定的网损经济性换来了更高的恢复可靠性和电压安全裕度。这个权衡在实际工程中是值得的因为故障恢复阶段的首要目标不是经济性而是供电可靠性。5.3 敏感性分析与计算效率$\Gamma$ 从0增大到10的过程中目标函数值逐步变差因为方案越来越保守但最坏场景下的安全性指标逐步改善。当 $\Gamma$ 超过8之后目标值改善幅度已经非常小说明过度保守并没有带来额外收益实际中可以据此选择一个合适的预算参数。计算效率方面33节点系统规模下CCG算法平均迭代次数在8到14次之间收敛单轮主问题和子问题求解时间在0.5到3秒之间总耗时通常不超过30秒。如果换到更大的系统比如123节点系统迭代次数会增加到20次以上此时需要关注主问题规模的增长速度。6. 常见问题与排查经验实录6.1 收敛慢、上下界震荡怎么处理CCG最常见的现象是前几轮LB和UB快速靠近但到了最后1%的gap时卡住不动。处理方法无非三种把MIP gap从默认值调小一点像我在代码里设置的mipgap1e-4检查子问题是否真的是最优解必要时给求解器设置更严格的容差如果UB反复震荡优先怀疑对偶问题推导有问题尤其是对偶变量符号和约束方向。另外要注意主问题每次新增场景约束后之前保存的所有场景变量都要保持固定否则求解器会调整历史场景对应的运行变量来迎合当前目标导致主问题退化成次优解。6.2 对偶推导容易出错的位置对偶推导的正确性决定整个CCG算法是否可行。最容易出错的是不等式约束的对偶变量非负方向搞反等式约束对偶变量没有设成自由变量目标函数优化方向改变后没有做负号转换内层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设置这些坑我基本都踩过一轮写在这里帮大家节省一些调代码的时间。结尾小经验整个项目跑下来我最大的体会是两阶段鲁棒恢复这个方向数学模型本身并不难理解真正的门槛在于把抽象的三层优化结构干净地落到代码里。CCG框架的代码骨架不到两百行但每一步推导、每个变量索引、每个对偶符号都可能在细节上出问题。我的建议是先从一个很小的系统开始比如把33节点简化为6节点甚至3节点系统手动验证每个中间结果再逐步扩展到完整算例。这样定位问题会快很多对模型的理解也会更扎实。