1. 项目概述与核心痛点拆解1.1 为什么需要两阶段鲁棒优化风电、光伏大规模并网之后电网调度面临的最直接问题就是不确定性。天气一变风功率曲线跟着变光照被云层遮挡光伏出力直接跳水负荷侧又有用户行为的随机波动。传统的确定性调度模型把所有参数都当成已知量来处理一旦实际风光出力偏离预测值制定的发电计划就可能不满足安全约束严重时甚至要拉闸限电。鲁棒优化Robust Optimization的思路和随机规划不一样它不假设不确定参数的概率分布而是用一个不确定集合来刻画参数的波动范围寻求的是最恶劣场景下依然可行的调度方案。这种思路在电力系统里特别实用因为风电、光伏、负荷的预测误差虽然有一定统计规律但准确建模概率分布往往很困难而且随机规划求解规模庞大。两阶段鲁棒优化的“两阶段”对应的是调度决策的两个层次第一阶段是机组启停、日前预调度等“现在就要定下来”的决策第二阶段是风光负荷不确定参数显现之后机组出力调整、切负荷、弃风弃光等“看情况再做”的决策。这种先后递进的决策结构非常契合电力系统“日前计划实时调整”的运行机制也是这个项目采用两阶段鲁棒优化模型的根本原因。从适用范围来说这个项目适合三类读者一是刚接触鲁棒优化、想看懂两阶段模型怎么建和怎么解的电力系统方向研究生二是从事调度算法开发的工程师需要一个可落地的参考框架三是想了解大M法和CCG算法在Matlab/YALMIP环境下怎么实现的控制类、运筹优化类学习者。1.2 建模语言与求解环境选择算例代码基于Matlab平台建模使用YALMIP工具箱求解器一般采用Cplex或Gurobi。YALMIP是Matlab环境下非常成熟的优化建模语言能用贴近数学表达的方式描述优化问题对新手友好对熟悉数学建模的人而言效率也高。这里有一个容易被忽略的细节同一套模型不同求解器性能差异很大。两阶段鲁棒优化在CCG迭代过程中会反复求解混合整数线性规划MILP如果模型规模大、迭代次数多求解器的选择很关键。Cplex和Gurobi是当前两个主流选择Gurobi在纯MILP求解上通常稍占优势但Cplex在部分电力系统模型的特殊结构上表现也不错。我建议学生阶段两个求解器都装上遇到求解瓶颈时可以对比。如果不想装商业求解器YALMIP也支持开源的SCIP等求解器但遇到大规模MILP时性能差距明显项目复现阶段建议优先考虑商业求解器。2. 两阶段鲁棒优化模型构建2.1 模型总体框架与决策变量划分两阶段鲁棒优化的核心思想是“在这里决策在那里应对”。第一阶段决策变量通常在不确定参数实现之前就必须确定包括机组启停状态0-1变量机组日前预调度出力系统的旋转备用容量安排第二阶段决策变量则是在不确定参数已知后做出调整包括机组实际出力调整量切负荷量弃风弃光量用数学语言描述两阶段鲁棒优化模型的标准形式如下$$ \min_{x} \left( c^T x \max_{u \in U} \min_{y \in \Omega(x,u)} b^T y \right) $$其中外层 min 是第一阶段的调度决策内层 max-min 是第二阶段问题的双层结构不确定性 u 试图最大化运行成本而运行层面 y 会在给定 u 的条件下争取最小化调整成本。每一轮 CCG 迭代就是在主问题中加入由第二阶段返回的“最恶劣场景”对应的约束不断逼近真实最优解。这种结构的妙处在于第一阶段求得的调度方案不只在预测场景下可行而且对不确定集合内所有可能的出力场景都能保证通过第二阶段的调整来恢复可行性。这比传统的确定性“预留备用”方法要精确得多因为它明确地考虑了所有极端场景的可行性约束。2.2 盒式不确定集合与预算约束风电、光伏、负荷的波动范围用不确定集合来描述。最常用的是盒式不确定集合Box Uncertainty Set$$ u_i \in [\bar{u}_i - \hat{u}_i, \bar{u}_i \hat{u}_i] $$其中 \bar{u}_i 是预测值\hat{u}_i 是最大偏差。方法很简单把所有可能的值都限定在一个区间内但这么做的代价是适应最坏情况的同时会过于保守。比如如果所有风光场站同时都取最大正偏差会导致系统预留大量备用经济性很差。所以实际工程中一般不直接用纯盒式集合而是引入预算约束Budget Constraint来控制参数的“同时波动”程度$$ \sum_i \frac{|u_i - \bar{u}_i|}{\hat{u}_i} \leq \Gamma $$这个约束的含义很直观所有不确定参数同时偏离预测值的总程度是有限的现实中不会所有风电场、光伏电站同时出现最大预测误差。\Gamma 越小模型越乐观调度方案经济性好但抗风险能力弱\Gamma 越大模型越保守安全性高但运行成本上升。具体取值需要结合历史预测误差数据的分布来标定一般经验值是取总场站数量的50%-70%。考虑到Matlab代码实现的便捷性实际建模中可以简化为到每个节点或每个场站单独设置波动偏差只需要在代码的输入数据中定义每台风机、光伏电站的预测上下限即可不需要额外的概率分布假设。2.3 目标函数与约束条件的数学表达在工程项目中第一阶段目标通常是最小化机组运行成本和启停成本第二阶段是调整成本与弃风弃光、切负荷惩罚费用之和。数学表达为$$ \min \sum_{t} \sum_{g} \left( C_{g}^{op} P_{g,t} C_{g}^{su} u_{g,t}^{su} C_{g}^{sd} u_{g,t}^{sd} \right) \max_{u \in U} \min \sum_{t} \left( \sum_{g} C_{g}^{adj} \Delta P_{g,t} C^{curtail} P^{curtail}{t} C^{load_shed} P^{shed}{t} \right) $$约束条件包括功率平衡约束各时刻总发电功率等于负荷功率加网损网损通常简化处理机组出力上下限约束机组出力必须在技术出力范围内机组爬坡约束相邻时刻出力变化量不能超过爬坡速率限制最小启停时间约束机组不能频繁启停有最小的开机/停机持续时间支路潮流约束传统方法用直流潮流近似转化为节点相角的线性约束如果考虑网络安全还要加入线路传输容量约束备用容量约束系统必须有一定比例的旋转备用以应对不确定性第二阶段约束还包括机组实际出力可调范围与第一阶段出力之间的关系以及切负荷量和弃风弃光量不能超过实际负荷和风光出力值。值得注意的是这个问题的约束规模会随着机组数量、节点数量、时段数量的增加而急剧增长典型的算例规模是24时段乘以数个节点乘以数台机组最后形成的MILP问题规模可能非常大。因此建模时要注意避免不必要的冗余约束能用等式约束的地方就不要扩大成不等式。3. 大M法与KKT条件转换3.1 为什么需要大M法做线性化两阶段鲁棒优化模型中第二阶段内部其实是一个双层优化问题max-min 结构。直接求解这种双层问题极其困难而且决策变量之间的乘积会引入非线性项例如0-1变量与连续变量的乘积互补松弛条件中的乘积项分段线性函数中的选择逻辑大M法Big-M Method是处理这些混合整数规划和非线性项的经典工具。核心思想是引入一个足够大的常数 M通过适当的约束构造将一个非线性的、逻辑性的条件转化为等价的线性不等式组。在电力系统鲁棒优化的具体场景中最经典的大M法应用是处理“互补松弛条件”的线性化。当我们将第二阶段 min 问题通过KKT条件转化为 max 问题时KKT条件中包含了形如“拉格朗日乘子 × 约束不等式松弛量 0”的互补松弛条件这是非线性且非凸的。引入大M和二进制变量就可以用以下方式线性化这类条件若存在互补条件 0 ≤ λ ⊥ μ ≥ 0则可引入辅助二进制变量 z ∈ {0,1}构造$$ \lambda \leq M z $$ $$ \mu \leq M(1-z) $$这个做法的物理含义很直接两个非负量不可能同时为正至少要有一个为零。大M就像是一个“开关”边界决定了某个量是否被允许非零。很多同学在这里会踩一个大坑大M的取值过大或过小都会严重破坏求解性能。M太小会错误地切掉可行解M太大会导致数值病态求解器精度下降甚至得到违反约束的错误解。实际调试过程中需要根据具体问题的数据量级反复尝试一般建议从量级上比目标函数中系数大100到1000倍左右开始试探具体操作我在第5节会更详细展开。3.2 第二阶段问题的对偶转换与M取值技巧另一种常用方法是直接将第二阶段问题取对偶把 max-min 问题转化为 max-max 问题即极大化一个对偶问题。这需要对原问题的拉格朗日函数求关于内层变量的极值得到对偶约束条件。对偶转换的数学操作比较繁琐尤其当约束条件较多的时候拉的乘积项会非常多极容易出错。我的经验是先手写推导一遍对偶问题再用Matlab的符号工具对部分表达式做验证。需要注意的是对偶转换的前提条件是内层问题必须是凸的且满足强对偶条件。在电力系统运行约束中大多数线性化约束满足这个要求。在项目实现中M的取值直接决定求解质量和速度。经验法则是首先统计所有决策变量和参数的数值范围找出最大量级以这个量级为基础取 M 1000 * max(|变量|) 作为初值测试求解后检查互补松弛条件是否严格成立即每个乘积项都接近0如果发现乘积项明显不为零则逐步增大M如果求解时间过长或不收敛则适当减小M还有一个更精细的调法不同约束可以用不同数量级的M。例如机组出力相关的约束用 M1e4而线路潮流的约束用 M1e5这比全局统一用一个大M要好得多能明显改善求解器的数值稳定性。4. CCG算法设计与收敛性分析4.1 CCG算法的核心思想与迭代框架CCGColumn-and-Constraint Generation算法全称列与约束生成算法是目前求解两阶段鲁棒优化问题的主流算法之一。与经典的Benders分解法相比CCG在处理带整数变量的鲁棒问题时有明显优势。CCG算法的思想可以类比为一个动态调整预算的过程先按预测场景做一个初始计划然后让不确定系统找出当前计划最危险的那个场景把应对这个场景所需的调整能力补进计划中如此反复直到计划能应对所有可能场景。具体框架如下第0步设定一个初始的预测场景通常取预测值求解主问题第一阶段问题得到调度计划 x第1步固定 x求解子问题第二阶段问题找出最恶劣场景 u*并得到对应的最优调整量第2步将新找到的 u* 对应的第二阶段决策变量和约束添加到主问题中第3步重新求解主问题更新生产计划回到第1步反复迭代直到目标函数上界和下界的间隙小于设定阈值用数学式表达主问题为$$ \min_{x, y, \eta} c^T x \eta $$ $$ s.t. \quad Ax \leq b, \quad \eta \geq b^T y_k, \quad Fx Gy_k \leq h - Eu_k^{*}, \quad \forall k \leq K $$其中 K 是已迭代次数u_k* 是第 k 次迭代得到的极端场景。主问题规模会随着迭代次数K逐步增大因为每轮都要加入一组新的变量和约束。而子问题是给定 x 后求解$$ \max_{u \in U} \min_{y \in \Omega(x,u)} b^T y $$这一步需要结合大M法或对偶大M把内层min问题转化为单层MILP/MILP问题来处理。4.2 CCG与Benders分解的对比分析我在实际项目对比中发现CCG相比Benders分解在很多场景下都有明显优势具体对比如下对比维度CCG算法Benders分解收敛速度快通常10次以内即可收敛较慢需要生成大量切平面对整数变量的支持良好可直接处理第二阶段整数变量较差需要额外的处理技巧主问题规模随迭代增加新增变量约束增长相对缓慢新增切平面实现复杂度中等需要维护场景集合较简单但推导复杂鲁棒问题适用性专门针对两阶段鲁棒设计传统分解方法非鲁棒专用CCG之所以收敛快是因为它每轮都在主问题中加入完整的新变量和新约束相当于在当前最恶劣场景下重新优化了所有变量而Benders分解只添加一个割平面信息量相对有限。在我的多个测试案例中CCG通常在3-8次迭代就达到10^-4的收敛精度而Benders可能需要数十次。不过CCG也有代价就是每轮迭代后主问题规模增长明显如果迭代次数多后期的MILP求解难度也会很大。因此对大规模算例可以先用Benders做一个粗糙的下界估计再切换到CCG追求精度这种混合策略在实际项目中很实用。4.3 收敛判据设置收敛判据是CCG算法的关键细节。标准做法是设定相对间隙Relative Gap$$ Gap \frac{|UB - LB|}{|LB|} \leq \epsilon $$其中上界UB来自子问题求得的最恶劣场景下的总成本下界LB来自主问题的最优目标函数值。工程上我建议将 \epsilon 设置为1e-3或1e-4。如果太宽松比如1e-2得到的调度方案可能在实际最恶劣场景下不够安全如果太严格比如1e-6迭代次数会显著增加而最后的调度方案差异其实很小。我在实际调参中发现1e-3和1e-4的结果差异通常不足0.1%但迭代次数会多出1-2轮。另外要留意一种特殊情况如果子问题是无界的目标函数趋于无穷说明主问题给出的调度方案在某些场景下根本不可行需要先检查主问题是否考虑了所有必要约束。这种问题通常和M取值不当或第二阶段对偶推导错误有关需要回到模型本身去排查。5. Matlab代码实现与关键细节5.1 算例数据准备与参数定义以6节点系统作为标准算例包含3台火电机组、1个风电场、1个光伏电站和若干负荷这是目前教材和论文中最常用的验证系统规模适中既能完整展示算法流程又不会因为复杂度过高导致调试困难。时间尺度选24小时步长1小时。典型参数设置如下参数类型具体数值火电机组台数3台风电场数量1个装机容量150MW光伏电站数量1个装机容量100MW峰荷约300MW风功率最大偏差预测值的±20%光伏功率最大偏差预测值的±25%负荷偏差预测值的±5%爬坡速率10-30MW/h有了这些参数后还需要对风光出力预测值做归一化处理。在Matlab代码中通常会使用以下结构来定义不确定集% 定义不确定集参数 u_wind_pred wind_forecast; % 风电预测值(每个时段) u_wind_dev 0.2 * u_wind_pred; % 风电波动范围 u_solar_pred solar_forecast; % 光伏预测值 u_solar_dev 0.25 * u_solar_pred; % 光伏波动范围 u_load_pred load_forecast; % 负荷预测值 u_load_dev 0.05 * u_load_pred; % 负荷波动偏差 Gamma 20; % 预算约束参数表示24小时内最多有20个时段会同时偏离预测值这里预算参数 \Gamma 的物理含义值得细说如果把所有时段都纳入不确定集合最恶劣场景可能是所有时段风光出力全部取最小、负荷全部取最大但这显然过于悲观。\Gamma 限定了这种“同时发生偏差”的场景数量工程上一般取总时段数的50%-80%。5.2 主问题与子问题的YALMIP建模要点在YALMIP中主问题的建模相对直接。关键是定义sdpvar对象和binvar对象并设置约束条件。下面是一段主问题建模的核心框架% 主问题变量 x_start binvar(n_gen, T, full); % 机组启停状态 x_pre sdpvar(n_gen, T, full); % 日前预调度出力 eta sdpvar(1, 1); % 第二阶段成本上界 % 主问题约束 Constraints []; % 机组出力上下限约束 for t 1:T for g 1:n_gen Constraints [Constraints, ... P_min(g) * x_start(g,t) x_pre(g,t) P_max(g) * x_start(g,t)]; end end % 功率平衡约束 for t 1:T Constraints [Constraints, ... sum(x_pre(:,t)) u_wind_pred(t) u_solar_pred(t) sum(load_base(:,t))]; end % 第二阶段成本约束 Constraints [Constraints, ... eta cost_second(n, y_vars)];子问题方面我们需要对第二阶段问题进行转换。假设第二阶段问题不含整数变量或已经通过大M法将互补约束线性化那么可以将原问题转化为对偶问题进行求解。这里用KKT条件大M法实现% 第二阶段子问题给定 x_pre求解最恶劣场景 % 引入不确定变量 u_wind, u_solar, u_load u_wind sdpvar(T, 1); u_solar sdpvar(T, 1); u_load sdpvar(T, 1); % 不确定集约束 Constraints_uncertain []; for t 1:T Constraints_uncertain [Constraints_uncertain, ... u_wind_pred(t) - u_wind_dev(t) u_wind(t) u_wind_pred(t) u_wind_dev(t)]; Constraints_uncertain [Constraints_uncertain, ... u_solar_pred(t) - u_solar_dev(t) u_solar(t) u_solar_pred(t) u_solar_dev(t)]; Constraints_uncertain [Constraints_uncertain, ... u_load_pred(t) - u_load_dev(t) u_load(t) u_load_pred(t) u_load_dev(t)]; end % 预算约束 Constraints_uncertain [Constraints_uncertain, ... sum(abs(u_wind - u_wind_pred) ./ u_wind_dev) ... sum(abs(u_solar - u_solar_pred) ./ u_solar_dev) ... sum(abs(u_load - u_load_pred) ./ u_load_dev) Gamma];这里最容易被忽视的地方是预算约束中的绝对值项不是线性约束需要引入辅助变量进行线性化处理。具体做法是引入非负变量表示正负偏差构造线性等式化形式% 引入正负偏差辅助变量 d_wind_pos sdpvar(T, 1); d_wind_neg sdpvar(T, 1); d_solar_pos sdpvar(T, 1); d_solar_neg sdpvar(T, 1); d_load_pos sdpvar(T, 1); d_load_neg sdpvar(T, 1); % 偏差定义 Constraints_uncertain [Constraints_uncertain, ... u_wind - u_wind_pred d_wind_pos - d_wind_neg]; Constraints_uncertain [Constraints_uncertain, ... d_wind_pos 0, d_wind_neg 0]; % 预算约束线性形式 Constraints_uncertain [Constraints_uncertain, ... sum((d_wind_pos d_wind_neg) ./ u_wind_dev) ... sum((d_solar_pos d_solar_neg) ./ u_solar_dev) ... sum((d_load_pos d_load_neg) ./ u_load_dev) Gamma];5.3 CCG主循环代码实现详解CCG的迭代入口如下% 初始化 UB inf; LB -inf; iter 0; max_iter 20; tol 1e-3; scenarios {}; % 保存所有极端场景 while (abs(UB - LB) / max(1, abs(LB))) tol iter max_iter iter iter 1; % 求解主问题考虑已有场景 [x_result, LB] solve_master_problem(scenarios); % 给定第一阶段的解求解子问题 [u_worst, second_cost, feasibility] solve_subproblem(x_result); if feasibility 1 % 计算上界 total_cost first_stage_cost(x_result) second_cost; UB min(UB, total_cost); % 将最恶劣场景添加到主问题场景集中 scenarios{end1} u_worst; else % 主问题不可行处理添加可行性割 scenarios{end1} generate_feasibility_scene(x_result); end fprintf(迭代: %d, LB: %.4f, UB: %.4f, Gap: %.6f\n, ... iter, LB, UB, abs(UB - LB) / max(1, abs(LB))); end这里的LB来自主问题最优值包含已经迭代的历史场景集合UB来自固定第一阶段方案后最恶劣场景下的完整成本。当我打印迭代记录时会特别关注每个阶段的LB和UB的变化轨迹如果LB经过多次迭代还在原地踏步说明主问题没有捕获到新的有效约束如果UB持续下降说明子问题在找的更恶劣场景有效。有一个实践经验是在第一轮迭代前可以把初始场景设为预测场景值这样主问题先求解出正常的确定性调度方案再让子问题去“攻击”该方案。这种方式能让前几轮迭代更有方向性避免初始场景过于极端导致主问题无解。5.4 结果可视化与调试技巧CCG算法收敛后需要对结果进行可视化检查这是论文写作和工程交流中不可或缺的一环。核心可视化内容包括各机组24小时的调度计划柱状图或堆叠图最恶劣场景下风电、光伏、负荷与实际出力的对照曲线各阶段迭代间隙的收敛曲线切负荷量与弃风弃光量的时间分布以机组出力计划为例代码实现% 绘制机组出力堆叠图 figure; bar(T, x_result, stacked); xlabel(时段/h); ylabel(出力/MW); legend(机组1, 机组2, 机组3); title(两阶段鲁棒优化机组调度计划); grid on;在调试中我发现YALMIP和求解器之间偶尔会出现变量名冲突或MATLAB工作区变量被意外覆盖的问题。每次迭代循环前建议对工作区关键变量做备份同时使用optimize函数的返回值检查求解状态。例如diagnostic optimize(Constraints, Objective, options); if diagnostic.problem ~ 0 warning(求解失败: %s, yalmiperror(diagnostic.problem)); end这样不仅代码健壮性更高排查起问题来也更快。6. 算例测试与结果分析6.1 CCG迭代收敛过程案例分析以某6节点系统为例设置风功率最大偏差20%、光伏25%、负荷5%预算因子 \Gamma0.7即24小时约17个时段可同时波动运行CCG算法后收敛记录如下迭代次数下界LB上界UB间隙Gap112853.6215542.3517.30%214132.0815387.228.90%314765.4315196.732.84%414955.7715110.461.02%515038.2115082.530.29%615061.3815070.120.06%可以看到前两次迭代Gap下降非常快这是因为前几次加入的最恶劣场景对调度方案的约束效果显著。到第4轮以后Gap下降趋势放缓到第6轮已经达到低于1e-3的收敛精度。总求解时间在2分钟以内使用I5处理器Gurobi求解器对于研究算例来说完全可接受。实际调试中如果前两轮Gap没有明显下降优先排查子问题是否真的找到了极端场景。一个有效技巧是人工设置一个明显偏差极大的场景手动加入主问题观察LB是否显著上升。如果LB没有反应说明主问题可能没有正确耦合第二阶段变量。6.2 不同预算因子对调度方案的影响保持其他参数不变只改变 \Gamma 值观察总成本的变动趋势。\Gamma 反映了决策者的风险偏好程度\Gamma0 相当于完全忽略不确定性退化为确定性优化\Gamma24考虑所有时段全波动或按总偏差上限调整相当于最保守策略。\Gamma 值总成本元相对\Gamma0增幅计算耗时0确定性10842.98-20s0.2511830.399.1%50s0.513424.6123.8%95s0.715070.1239.0%120s1.0完全保守17543.7861.8%180s这个表非常直观地说明了一个工程权衡鲁棒性不是免费的提高抗风险能力意味着成本显著上升。在实际电力系统运行中调度部门需要结合对预测精度的判断选择一个合理的 \Gamma 值既不盲目乐观也不过度保守。6.3 风电、光伏、负荷不确定性对结果的对比分析为了分析不同类型不确定性的影响幅度可以分别只让一种参数波动其他固定为预测值观察各场景下的成本差异不确定类型确定性成本鲁棒优化成本成本增加比例仅风电波动10842.9812654.3116.7%仅光伏波动10842.9811892.079.7%仅负荷波动10842.9811347.664.7%三者同时波动10842.9815070.1239.0%可以看到风电波动的影响最大因为风电场装机容量大且波动区间大负荷波动由于数据本身预测精度较高波动区间设置为±5%所以影响最小。这个结果也间接验证了模型对不同类型不确定源的处理能力对项目工作量的分配提供了参考如果把精力集中在提高风电预测精度上比提高负荷预测精度能获得更大的成本收益。7. 常见问题与排查技巧实录7.1 求解不收敛的常见原因与对策CCG算法不收敛这个问题我在调试中碰到不下十次。总结下来主要原因集中在以下几类第一类是子问题无界或不可行。常见原因是主问题给出的调度方案在某些场景下不能满足负荷平衡约束尤其当遇到极端场景时系统所有机组的调整能力加起来都不够。解决方法是检查第二阶段模型中是否加入了“切负荷”和“弃风弃光”这两个松弛变量。真实系统中极端天气下为了保安全切负荷是允许的但会有高额惩罚成本。这两个变量相当于给了模型一个“兜底选项”能有效避免子问题因找不到可行解而崩溃。第二类是主问题规模增长过快导致求解器内存不足。CCG每轮迭代都会向主问题添加一组新变量和新约束迭代到10轮以后主问题可能有数千个变量求解时间显著增加甚至出现求解器卡死。应对策略是设置最大迭代次数一般20次就够了和求解时间限制以及定期对已有场景做精简例如删除对结果影响极小的历史场景。第三类是大M取值不当导致的数值问题。这个问题最隐蔽。表现为求解器报告“infeasible”但不给任何参考原因或者结果违反明显约束。我遇到过一次某个约束的M值设了1e6结果求解器在精度上出了问题实际上得到的解并不满足互补松弛条件。后来把M降到1e4问题立即恢复稳定。经验是M值够用即可不要贪大。7.2 YALMIP编程中常见语法与数值陷阱YALMIP虽然上手快但有一些细节必须注意。一个非常常见的错误是变量维度的隐式扩展。例如两个向量相加时如果一个是列向量、一个是行向量YALMIP会报维度错误但有时它会自动进行隐式广播这会导致约束实际构造得和预期不符。所以建模前一定要检查变量声明使用size函数确认维度避免在循环中出现矩阵维度意外扩增。另一个常见坑是optimize的默认参数设置。求解MILP时如果不设置求解器参数YALMIP会使用求解器的默认配置这些默认配置通常不是最优的。建议显式设置几个关键参数options sdpsettings(solver, gurobi, ... gurobi.mipgap, 0.001, ... gurobi.timelimit, 600, ... verbose, 1);高斯提示就是prove here的mipgap设置要与外层CCG的收敛精度匹配。如果外层要让Gap小于1e-3而内层MILP的mipgap也设为1e-3那么两个误差叠加后外层Gap可能无法达到预期。建议内层mipgap设置为外层精度的五分之一到十分之一。7.3 求解效率优化建议两阶段鲁棒优化的计算瓶颈基本都在MILP求解上。提高效率的经验法则是优先采用CCG而不是Benders因为前者迭代次数更少在CCG中延用上次求解得到的MILP解作为热启动点可以显著加速使用支持多线程的求解器Gurobi默认可利用多核如果模型规模非常大考虑将全时段的耦合约束爬坡约束等进行松弛或者减少精确建模的时段数先用粗粒度模型测试算法流程再逐步精细化还有一个实践中很管用的做法先在确定性场景下用一个较小的MILP模型验证求解器安装和YALMIP的建模流程确认无误后再扩展到完整的两阶段鲁棒优化。这种“从小到大、从简到全”的自底向上调试策略极大节省了我调试复杂模型的时间。8. 进一步扩展与个人实操体会8.1 模型扩展方向探讨这个两阶段鲁棒优化框架有很好的扩展性。方向上可以考虑把盒式不确定集合换成多面体不确定集合或数据驱动的凸包不确定集合后者能利用历史数据的分布信息降低解的保守性将交流潮流约束纳入第二阶段的可行性校验此时模型变为混合整数非线性规划MINLP需要结合凸松弛或线性化技术处理加入储能系统作为新的灵活性资源储能的作用在鲁棒优化框架下尤其值得研究因为它可以在时间维度上转移能量能显著缓解不确定性的影响改为分布鲁棒优化Distributionally Robust Optimization将不确定参数的分布信息引入模型兼顾鲁棒性和经济性8.2 我在代码调试过程中的经验总结第一次完整跑通这个项目的代码时前后花了接近两周时间中间踩过不少坑其中最有价值的一条是一定要在写完整代码之前先把数学模型中每个约束的维度、上下界和物理含义理清楚并用小规模算例手工验证一遍。很多看上去莫名其妙的bug比如约束没生效、变量没关联、求解器报无界等根源都在于建模时某些约束遗漏或者下标对应错误。还有一个小技巧使用YALMIP时可以通过assign函数给变量赋初值然后用value函数查看优化后的变量结果。这个调试工具非常有用能帮你逐步追踪每一类变量是否符合物理直觉而不是只看最终的目标函数值就草草收手。另外在打印迭代日志时我会在每一轮输出最恶劣场景的取值向量例如“第3轮最恶劣场景风电12时段取最小值、光伏17时段取最小值、负荷5时段取最大值”。这种可视化能直观地验证子问题是否在寻找逻辑上合理的极端场景。如果发现最恶劣场景的逻辑明显不合理比如风电和光伏同时取极小值但负荷却没有变化那就需要回头检查不确定集约束是否建对预算约束是否生效。8.3 基于这个框架可以继续学习的方向如果你是这个领域的初学者跑通这个例子之后建议按以下顺序继续深入动手修改不确定集合的类型比如把对称盒式改成不对称的多面体集合感受不同不确定集合对解的保守性和计算复杂度的影响尝试用Benders分解同样实现一遍两阶段鲁棒优化和CCG做对比这样能直观理解两类算法的本质差别把模型耦合到实际的IEEE 30节点或118节点系统上体验大规模系统下求解压力的骤增以及相应的加速手段尝试加入需求响应、储能、备用联动等运行措施观察它们对系统鲁棒性和经济性的影响鲁棒优化这个方向入门时公式推导看起来有些繁琐但一旦用代码实现了完整的求解框架后续很多工作都可以基于这个框架做扩展。我在实际教学中发现能独立复现这个项目的学生再去读高水平的鲁棒优化论文基本不会再卡壳。希望这篇记录能帮你避开我当年踩过的那些坑用更短的时间跑通全部流程。
