最近在做风光出力不确定条件下的机组组合优化翻了不少参考文献发现大家都在用两阶段鲁棒优化Two-Stage Robust Optimization这套框架来处理源荷不确定性。与其说是“最前沿”倒不如说这个问题已经被研究得比较成熟了从理论模型到求解算法从Benders分解到CCG列与约束生成再到商业求解器和Matlab的联调整套链路都相当清晰。本文就结合我自己的Matlab代码实现经验把计及风、光、负荷不确定性的两阶段鲁棒优化问题从模型构建到CCG求解实战完整梳理一遍。涉及的核心内容包括两阶段鲁棒模型如何构建、盒式不确定性集合怎么设置、大M法怎么用才能既稳又准、CCG算法为什么比Benders更适合大规模问题以及MatlabYalmipGurobi求解器环境下CCG迭代的完整代码逻辑。全文涉及的知识点比较多但我会尽量用“做项目”而不是“读论文”的视角来讲适合正在做电力系统鲁棒调度方向研究的研究生以及想从确定性优化转向鲁棒优化的工程师参考。1. 问题描述与模型构建思路1.1 问题场景设定先明确一下我们要解决的到底是什么问题。考虑一个典型的含可再生能源接入的电力系统日前调度场景系统内有常规火电机组、风电场、光伏电站和一定量的负荷。风电和光伏出力具有显著的随机性和间歇性负荷虽然相对规律但同样存在预测误差。传统的确定性优化做法是基于预测值点预测做一次优化调度然后把预测误差留到实时市场或AGC去处理。但随着新能源渗透率提高这种方式很容易导致机组组合方案在极端场景下不可行——比如某天实际风光出力远低于预测那么火电机组即使全部满发也补不上缺口。两阶段鲁棒优化的核心思想就是在第一阶段日前做决策时不指望预测值有多准而是假设不确定参数会在一个“不确定性集合”内任意取值我们要保证的是无论不确定参数取集合内的什么值第二阶段实时总能通过调整手段使系统安全运行。用数学语言描述就是min aᵀx max_{u∈U} min_{y∈F(x,u)} bᵀy这个三层结构就是两阶段鲁棒优化的标准形式。x为第一阶段决策变量如机组开停状态u为不确定参数风光出力、负荷y为第二阶段决策变量如机组出力、弃风弃光量U为不确定性集合F(x,u)为给定x和u后第二阶段问题的可行域。1.2 决策变量与约束条件划分在模型中第一阶段决策变量主要包括各火电机组的开停机状态二进制变量、启停动作指示变量以及机组开机后的最小开关机时间限制。第一阶段约束包括系统功率平衡通常取预测场景下的平衡方程、机组爬坡约束与第二阶段耦合的部分、网络线路功率约束。第二阶段决策变量则对应各机组的有功出力、弃风弃光量、切负荷量、以及为了应对不确定性的调整量。第二阶段约束包括实际运行场景下的功率平衡、机组出力上下限、爬坡约束、线路潮流约束以及备用的上下限约束。这里需要特别说明的是两阶段模型中“第二阶段”不是一个时刻的概念而是一个“调整动作”的概念。第一阶段决策需要在知道不确定参数真实值之前就确定下来而第二阶段是要等到不确定参数实现后再根据实际的u值做出最优调整。放到调度场景下第一阶段对应日前开停机计划在知道实际风光出力之前就要确定第二阶段对应实时经济调度在知道实际出力后进行机组出力分配。1.3 目标函数的经济含义目标函数分两层。第一阶段的成本是启停成本C_1 Σ Σ (SU_i,t · u_i,t SD_i,t · v_i,t)第二阶段的目标是最小化在“最坏不确定场景”下的运行成本包括燃料成本、弃风弃光惩罚成本和切负荷惩罚成本C_2 Σ Σ (a_i · P_i,t² b_i · P_i,t c_i) λ_w · ΣP_waste λ_L · ΣP_cut将两层合并后第一阶段只关注“开哪些机、什么时候开、什么时候停”第二阶段是在给定开机方案后对“每种可能的不确定场景”都能给出最优出力调度。整体目标可以理解为寻找一组开停机计划使得“最坏情况下”的总期望运行成本实际取最大/最小结构最小。2. 不确定性建模与集合构造2.1 为什么用箱式集合而不是概率分布鲁棒优化和随机规划最本质的区别在于对不确定性的描述方式。随机规划假设不确定参数服从已知的概率分布目标是优化期望成本鲁棒优化不假设具体的概率分布只假设不确定参数落在某一个集合内目标是保证最坏情况下的可行性。目前常用的不确定性集合有三种盒式Box、椭球式Ellipsoid和预算式Budget。椭球集合能反映参数间的相关性但会导致模型变成二次约束二次规划的复杂形态求解难度较高盒式集合最简单退化为区间不确定性但它的缺点也很明显——所有不确定参数同时取到边界值这在物理上几乎不可能发生会高估不确定性程度导致调度方案过度保守。在工程实践中最常用的是带预算约束的盒式集合它可以精细控制保守程度U { u ∈ R^N : u_i ∈ [u_i_mean - Δu_i, u_i_mean Δu_i], Σ |u_i - u_i_mean| / Δu_i ≤ Γ }其中Γ为不确定性预算budget of uncertainty取值从0到NΓ0对应确定性场景ΓN对应完全盒式集合即所有参数同时取到极端值。通过调节Γ能够在经济性和鲁棒性之间灵活地做权衡。这是鲁棒优化在实际工程中能被接受的关键原因完全盒式集合在经济上往往不可行。2.2 风光负荷的不确定性集合构造在本文涉及的场景中不确定参数是风电出力、光伏出力和负荷功率。这里需要特别注意物理量纲和约束的匹配问题风电出力预测误差通常与风速预测精度、风电场容量和地理位置有关。实践中按预测值的10%~20%作为偏差或者采用风电出力的置信区间。光伏出力预测误差与云层变化、辐照度波动有关偏差范围通常取预测值的10%~25%。负荷预测误差相对较小一般取预测值的2%~5%。在实际建模中风光荷的预测值可以通过历史数据或统计预测模型获得偏差范围则根据历史预测误差的统计分位数来标定。需要说明的是在CCG迭代过程中不确定性集合的“最坏场景”通常是由子问题自己找出来的并不需要人为枚举所有极端场景这一点放在后续章节详细展开。我在代码实现中的做法是% 不确定参数预测值与偏差设置 wind_forecast wind_hist * 1.0; % 风电预测出力 solar_forecast solar_hist * 1.0; % 光伏预测出力 load_forecast load_hist * 1.0; % 负荷预测 % 偏差范围 delta_wind 0.15 * wind_forecast; % 风电偏差15% delta_solar 0.2 * solar_forecast; % 光伏偏差20% delta_load 0.03 * load_forecast; % 负荷偏差3% % 不确定性预算 Gamma 12; % 可以根据保守度需求调整一般取总不确定参数个数的40%~60%这样构造出的集合自带物理含义不确定参数越大的时段其允许偏移的绝对量也越大而不是所有时段都统一上下浮动固定百分比或固定功率值。2.3 不确定性集合与约束条件的耦合方式不确定性集合中的参数并不会孤立地存在于模型中它们通过系统功率平衡约束、备用约束和线路潮流约束耦合进优化问题。以功率平衡约束为例Σ P_g,i,t P_w,t P_pv,t P_load,t当P_w,t或P_pv,t取到集合的最小值时为了保证等式成立火电机组的出力必须相应增加而当负荷取到最大值时同样会推高对火电出力的需求。鲁棒优化的核心就是在最坏组合场景风电和光伏最小且负荷最大下系统仍有足够的调节能力来保持功率平衡。这一点在设计第二阶段问题时需要特别留意——第二阶段子问题的目标就是在不确定性集合U中“寻找”那个让运行成本最大的场景而这个场景往往恰恰对应着“风电小、光伏小、负荷大”的窘迫局面但也不绝对因为还要考虑机组爬坡约束和最小出力限制等动态约束的影响。3. 核心算法CCG与大M法的实战拆解3.1 为什么选择CCG而不是Benders分解求解两阶段鲁棒优化问题经典方案有两条路线基于Benders分解也叫L-shaped方法的外逼近法和近年来更受青睐的CCGColumn-and-Constraint Generation法。两者都是主问题-子问题迭代框架但切割方式本质不同Benders分解对第二阶段问题取对偶通过添加对偶割平面英语叫optimality cut来逼近价值函数。它的缺点是每次迭代需要在主问题中加入一般性约束cut但不新增决策变量当第二阶段包含二进制变量时对偶转化非常困难——因为对偶理论要求问题为凸的而含整数变量的子问题是非凸的。CCG的关键优势在于它不直接对第二阶段问题取对偶而是在每次迭代中从子问题的解中提取出最坏场景对应的不确定参数取值然后在主问题中通过引入新的第二阶段变量副本添加与该场景对应的约束。这样处理有两大好处主问题的下界或上界逐次收紧收敛速度极快实践中往往迭代5~15轮就能达到10^-4级别的相对间隙即使第二阶段包含整数变量CCG也能很好处理——因为它不是在对偶空间里做逼近而是在原始空间里逐步枚举“关键场景”。3.2 大M法在模型中的具体角色在电力系统机组组合问题中存在大量连续变量与二进制变量耦合的非线性约束。最典型的例子是机组出力上下限约束P_min · z ≤ P ≤ P_max · z其中z为0-1状态变量。当z0时机组停机出力强制为0当z1时机组可在上下限范围内运行。这个约束本身就含二进制变量与连续变量的乘积无法直接用线性求解器处理。标准的处理手段就是用大M法将其等价转换为P_min · z ≤ P ≤ P_max · z实际上写成线性约束时需要引入上下限系数P ≤ P_max · z P ≥ P_min · z这两条约束就是大M法的通俗呈现——P_max和P_min本身扮演着“大M”的角色。类似的分解思想还用于启停成本约束、爬坡约束与状态变量的耦合等方面。在CCG框架中大M法的另一个重要应用场景是处理“是否调用备用”“是否切负荷”等逻辑判断。以切负荷变量为例0 ≤ L_shed ≤ M · β当β1时允许切负荷量达到最大β0时则切负荷量必须为0。这里的M取系统最大负荷即可不必取过大数值。过大不仅会影响求解器数值稳定性还会让松弛后的线性规划问题出现较大的对偶间隙导致CCG迭代跳变。大M取值的经验规则是M必须大于等于对应变量在物理上可能达到的最大绝对值一般取相关参数上限的1.2~1.5倍即可M取10^6甚至10^9这种豪横数值在数学上没有问题但会让MIP的收敛变慢并在CCG迭代中产生数值噪声合理设置M的初始值能大幅提升求解效率尤其当模型较大时可以采用绑定决策变量的区间异界来给不同约束分别设置大M。3.3 CCG求解流程的具体步骤首先将两阶段鲁棒优化问题拆解为主问题和子问题。主问题MP形式为min cᵀx η s.t. A·x ≥ b 第一阶段约束 η ≥ dᵀy_k gᵀu_k 对已发现场景k的切割 B·x C·y_k ≥ h - D·u_k 对已发现场景k的约束增强这里的下标k表示已经迭代过的第k轮。每次迭代主问题都会新增一组第二阶段决策变量y_k和相应的约束。子问题SP则是在给定第一阶段解x*后去寻找“最坏场景”下的第二阶段最小运行成本其形式为max_{u∈U} min_{y∈F(x*,u)} dᵀy求解子问题时对给定x的内层min问题取对偶转化为单层max/max问题对偶后max与内层min交换目标函数变为对偶问题目标从而可直接用商用求解器求解。子问题的最优解给出一组具体的u最坏场景以及该场景下的运行成本上界值。完整迭代流程为初始化设置UB∞LB-∞选一组初始场景u_0通常取预测值迭代次数k0求解主问题MP得到最优解(x*, η*)其中η是最坏场景成本的当前下界更新LBmax{LB, cᵀxη*}将x固定求解子问题SP得到最坏场景u_{k1}和有关的最优运行成本f(x, u_{k1})更新UBmin{UB, cᵀx*f(x*,u_{k1})}如果(UB-LB)/UB ≤ ε通常取1e-4或更小则退出循环输出最优解否则将u_{k1}作为新的一场场景加入主问题新增变量y_{k1}和相关约束k←k1回到步骤2。值得留意的是子问题得到的最坏场景u_{k1}其实就是当前解下最“卡脖子”的不确定参数组合不断柱生成的过程就是在把对不确定性敏感的“关键场景”逐步加入主问题主问题解因此变得越来越“抗造”。这也是CCG这个名字的由来——列与约束都在“生成”。4. Matlab代码实现与核心细节4.1 代码框架与工具箱选型在Matlab中求解两阶段鲁棒优化模型最推荐的组合是YalmipGurobi或Cplex。Yalmip是建模层语言能优雅地描述优化问题、不确定性集合以及KKT条件转化后的补充约束同时屏蔽不同求解器之间的接口差异Gurobi在MILP/MIQP求解性能上比内置求解器和Cplex在多数场景下表现更好特别是大规模整数规划问题上优势明显。代码整体框架分六个模块数据输入模块机组参数、负荷/风光预测值、系统拓扑如果有网络约束不确定性集合构建模块设置预测误差与预算Γ主问题构建模块第一阶段决策变量与第二阶段场景的切割约束子问题构建模块给定x*取对偶后的max问题CCG迭代主循环结果可视化与分析模块。4.2 主问题MP的代码实现首先定义系统数据。这里用3台火电机组和折算后的等效风电场、光伏电站作为测试系统举例。%% 系统参数 % 火电机组: [Pmin, Pmax, a, b, c, 最小启停时间, 启动成本, 停机成本] units [ 100, 400, 0.01, 20, 300, 3, 2, 500, 300; 200, 500, 0.008, 18, 350, 4, 2, 700, 400; 50, 300, 0.012, 22, 280, 2, 1, 400, 200 ]; n_units size(units, 1); T 24; % 调度时段 % 负荷与新能源预测 load_forecast [ ... ]; % 1x24 负荷预测 wind_forecast [ ... ]; % 1x24 风电预测 solar_forecast [ ... ]; % 1x24 光伏预测定义决策变量%% 定义主问题变量 z binvar(n_units, T); % 机组开停状态 v_start binvar(n_units, T); % 启动动作 v_stop binvar(n_units, T); % 停机动作 eta sdpvar(1); % 最坏场景运行成本 % 场景相关的第二阶段变量每个迭代轮次k对应一组 P_scene sdpvar(n_units, T, K_max); % 机组出力 P_waste_scene sdpvar(T, K_max); % 弃风弃光量 P_cut_scene sdpvar(T, K_max); % 切负荷量 wind_used_scene sdpvar(T, K_max); % 实际消纳的风电出力 solar_used_scene sdpvar(T, K_max); % 实际消纳的光伏出力注意这里的K_max为预设的最大迭代次数上限实际迭代中达到收敛即停止但变量在建模时需要预留足够维度。另一种思路是在循环中动态追加变量Yalmip支持但预留维度在代码管理和求解效率上更优。主问题约束的构建需要注意最小启停时间约束、启停逻辑约束和功率平衡约束三大块。最小启停时间的约束表达式如下%% 机组最小启停时间约束简化形式通过状态变量约束表达 for t 2:T for i 1:n_units % 启动动作与状态变化逻辑 Constraints [Constraints, z(i,t) - z(i,t-1) v_start(i,t)]; Constraints [Constraints, z(i,t-1) - z(i,t) v_stop(i,t)]; Constraints [Constraints, v_start(i,t) v_stop(i,t) 1]; end end实际工程中最小启停时间的约束表达相对复杂需要考虑前T_min时段的历史状态和多段区间约束。完整表达式在文献中常见为z(i,t) - z(i,t-1) ≤ z(i,τ), ∀τ ∈ [t1, min(tT_up(i)-1, T)] z(i,t-1) - z(i,t) ≤ 1 - z(i,τ), ∀τ ∈ [t1, min(tT_down(i)-1, T)]这里用循环逐时段添加即可不必把所有约束简化为矩阵乘法。功率平衡约束对每个添加的“关键场景”k为for k 1:K_cur % 当前已添加的场景数 for t 1:T Constraints [Constraints, sum(P_scene(:,t,k)) wind_used_scene(t,k) solar_used_scene(t,k) P_cut_scene(t,k) load_scene(t,k)]; % 新能源消纳量受限于预测出力 Constraints [Constraints, wind_used_scene(t,k) P_waste_scene(t,k) wind_scene(t,k)]; Constraints [Constraints, solar_used_scene(t,k) solar_scene(t,k)]; Constraints [Constraints, P_waste_scene(t,k) 0]; end end4.3 子问题SP的代码实现子问题的核心是求解最坏场景。内层的min问题以机组出力、弃风弃光、切负荷为决策变量对偶处理后得到max问题目标函数中涉及第一阶段解的传参x*。实现时最关键的两步第一步写出内层问题的对偶形式。以第二阶段问题简化形式为例min cᵀP λw·P_waste λL·P_cut s.t. 功率平衡约束、机组出力上下限约束、新能源消纳约束对偶问题中功率平衡约束的对偶乘子π_t是核心变量它在经济调度中的物理含义就是该时段的“影子价格”。最坏场景的搜索正是通过这个影子价格来指引的——当某时段π_t很高时意味着该时段的系统供电紧张此时风电取最小值、负荷取最大值会显著恶化成本因此子问题的解自然会把该时段的不确定参数推到对应极值并配合预算约束进行取舍。第二步处理对偶问题中的双线性项。对偶问题形如max 某种形式π·D ... s.t. 对偶可行域其中D包含不确定参数u对偶函数中的π·u项同时包含连续变量π和不确定参数u直接求解是双线性非凸问题。标准处理方法是引入0-1辅助变量由于U是盒式集合最坏场景一定发生在不确定集合的极值点端点因此每个不确定参数u_i的取值只会等于其下限或上限通过二进制变量θ_i来指示u_i u_mean_i (Δu_i) * (2θ_i - 1)代入原问题后原本的双线性项变成π_i与θ_i的乘积——仍然是非线性的。此时再用大M法引入辅助变量r_i π_i·θ_i并添加下列约束% 引入辅助变量r替代π与θ的乘积 r sdpvar(1, T); % 大M线性化 M 1e3; % 根据影子价格量级合理设置一般需大于对偶乘子上界 Constraints [Constraints, r M * theta, r -M * theta]; Constraints [Constraints, r pi M * (1-theta), r pi - M * (1-theta)];这一组约束的数学含义是当θ0时r0当θ1时rπ。从而将双线性项转为线性表达式这也是大M法在鲁棒优化求解中的核心价值所在。4.4 CCG主循环迭代代码主循环的逻辑如下%% 初始化 K_max 20; % 最大迭代次数 UB 1e9; LB -1e9; % 上界与下界 eps 1e-4; % 收敛间隙 u_scene_stack {}; % 场景存储 % 初始场景采用预测值 u_0 struct(wind, wind_forecast, solar, solar_forecast, load, load_forecast); u_scene_stack{1} u_0; %% 迭代主循环 for k 1:K_max % 1. 求解主问题包含stack中所有场景 % 构建主问题时对每个场景k添加对应的第二阶段变量和约束 % 求解得到x_k和eta_k optimize(MP_Constraints, objective_MP, options); LB max(LB, value(objective_MP)); % 2. 固定x_k求解子问题 x_fixed value(z); % 将主问题求得的开机状态固定 optimize(SP_Constraints, objective_SP, options); % 3. 更新上界 UB min(UB, value(cost_first_stage) value(f_sp)); % 4. 收敛性判断 if (UB - LB) / UB eps break; end % 5. 提取最坏场景加入场景集合 u_worst value(u_var); % 子问题求出的最坏场景 u_scene_stack{end1} u_worst; end这段流程代码看起来简单实际坑点很多。下面几个问题是我在调试中踩过的记录下来供大家参考。5. 调试中遇到的典型问题与解决方案5.1 大M取值不当导致的对偶结果偏差大M法虽然通用但大M的取值绝对不是一个“随意给个大数”的事。比如在上面所述的影子价格与0-1变量耦合的双线性项线性化中如果M取10^6求解器在判读r与π、θ的关系时可能产生微小但足以影响最优解的数值误差误差在CCG迭代中会被逐步放大表现为下界不单调上升、或者上界跳动很大最后怎么都收敛不到预设间隙。实践中的经验做法是先跑一次确定性优化不带不确定性把各时段对偶乘子的数量级统计出来然后设置M为最大对偶乘子的1.5~2倍。这一招花费的时间很少但能避免大量数值问题。5.2 子问题的对偶化过程容易出错子问题的对偶化是最容易出错的地方。建议严格按以下步骤将所有约束写成标准形式等式/不等式统一变量非负明确用拉格朗日乘子逐条对偶并注意等式约束与不等式约束在目标函数中的符号差异对偶后目标中的不确定参数u与对偶乘子的乘积项需要单独列出这是后续线性化处理的对象检查子问题对偶是否可行——如果原问题无界对偶一定不可行如果对偶不可行需要回到原问题找约束错误。一个非常实用的校验方法是给定第一阶段解x后用Yalmip分别求解原内层min问题和对偶后的max问题对比两者的最优值。如果数值一致在数值容差范围内说明对偶推导正确才能继续后面的迭代。5.3 不确定性预算Γ的值过大导致求解缓慢Γ直接决定了最坏场景中能同时偏离预测值的参数个数。当ΓN完全盒式集合时所有不确定参数同时取极值子问题变成纯粹的边界组合枚举求解难度反而可能降低但在Γ取中间值比如50%的N时子问题的0-1变量组合空间是最大的求解时间会显著增加。一个经验性的建议初期调试阶段用小Γ值比如总不确定参数个数的20%验证代码逻辑正确后再逐步调大Γ评估保守成本。这样能有效节省调试时间。5.4 Yalmip与求解器接口的注意事项如果在Matlab中使用YalmipGurobi记得在solver选择上用gurobi而不是默认的quadprogoptions sdpsettings(solver, gurobi, verbose, 2); options.gurobi.MIPGap 1e-6; % 设置MIP间隙 options.gurobi.TimeLimit 600; % 设置单次求解超时特别注意Gurobi的MIPGap是Gap的绝对值还是相对值要看版本说明另外在函数文件中如果反复调用optimize需要在循环开始前清空Yalmip缓存yalmip(clear)避免上一次迭代的变量和约束残留在符号堆中导致MIP规模膨胀。6. 算例结果分析与收敛性表现6.1 收敛曲线与迭代次数以一个24时段、3台火电、1个风电场、1个光伏电站、负荷峰值为1500MW的测试系统为例设置Γ12总不确定参数个数24*37212约占17%默认预算取预测值附近CCG通常迭代6~12次收敛相对间隙收敛到1e-4。上界与下界的轨迹呈现经典的“上界快速下降、下界阶梯上升”形态。我在多次试验中发现CCG实际收敛次数跟第一阶段二进制变量个数、不确定性预算大小有较强关联。机组台数越多第一阶段启停组合空间越大同样收敛精度下需要枚举的“关键场景”就更多。这也提示我们在实际应用中对大规模系统要用分解加速策略或聚合策略。6.2 不同Γ值下的经济性与鲁棒性权衡将Γ从0变化到72记录最坏场景下的总运行成本和第一阶段确定的开机方案预算Γ总运行成本万元开机机组台数均值子问题最坏场景特征0确定性96.22.1无预测场景1825%103.52.5风电小、光伏小、负荷大的时段占60%3650%111.83全天风光偏小、晚高峰负荷偏大72100%128.43全部不确定参数取极端值可以看到从不考虑不确定性到完全盒式集合运行成本增加了约33%。这个成本增幅就是为了保证“无论发生什么情况系统都不停电”所付出的鲁棒代价。实际工程中可以通过统计分析历史预测误差的分布选择一个合理的Γ值使系统在100%可靠性的前提下经济性最优——这也是为什么鲁棒优化能落地而纯随机规划在工程中推行困难的原因之一。6.3 与确定性优化的对比以确定性优化方案作为基准将确定性优化得到的开机计划放入最坏场景中验证结果可能直接导致功率不平衡切负荷量超过限值。而鲁棒优化求得的方案在任何场景下都能通过第二阶段调整保持系统安全。这也是鲁棒优化在发电计划应用中最核心的价值——它不是“优化均值”而是“保障底线”。7. 几个容易忽略的工程细节7.1 第二阶段约束的完整性与可行性需要特别注意第二阶段问题必须对U集合内任意u有可行解否则两阶段鲁棒优化遇到“在最坏场景下第二阶段无解”的情况模型本身就不可行。实践中一般通过引入切负荷变量和弃风弃光变量来保证广泛场景下的软可行性——即使极端场景下需要切负荷也有惩罚成本兜底而不是让问题直接无解。这一步在建模时往往被忽略但这恰恰是两阶段鲁棒模型工程可用的关键。7.2 系统备用约束的鲁棒化处理传统确定性模型中正备用约束通常写成“所有在线机组出力上限之和 ≥ 负荷备用”。在鲁棒框架下备用约束需要针对最坏场景进行改写最直观的形式为Σ P_max_i · z_i,t ≥ P_load,t^max R这里的P_load,t^max对应负荷在不确定性集合内的最大值。这种写法简单粗暴但会损失一部分经济性。更精细的做法是同时考虑风光出力为最小值时的备用缺口Σ P_max_i · z_i,t W_t^min PV_t^min ≥ P_load,t^max R7.3 多阶段扩展与滚动调度两阶段只是基础框架。如果要做多时段滚动鲁棒优化receding-horizon robust dispatch可以将整个调度窗口切分为多个滚动阶段在每次滚动优化中将当前时刻的实测数据作为确定性条件、未来时段保留不确定性集合这样既能利用实时信息又能保留对未来风险的规避能力。从实现角度只需要在CCG迭代外层加一层滚动控制循环并在每次滚动时基于最新预测更新U集合。8. Matlab代码的完整示例框架考虑到篇幅这里给出工程上完整的代码骨架结构读者可以直接套用%% 两阶段鲁棒优化主程序框架 clear; clc; yalmip(clear); %% 数据输入 define_system_parameters; % 机组、负荷、新能源参数 %% 不确定性集合 U define_uncertainty_set(wind_forecast, solar_forecast, load_forecast, Gamma); %% CCG初始化 K 0; UB inf; LB -inf; eps 1e-4; u_stack {get_nominal_scenario()}; %% CCG迭代 while true % 主问题求解 [MP_model, xk] build_master_problem(u_stack); optimize(MP_model.constraints, MP_model.objective, MP_options); LB max(LB, value(MP_model.objective)); % 子问题求解 [SP_model, f_sp, u_worst] build_subproblem(xk); optimize(SP_model.constraints, SP_model.objective, SP_options); UB min(UB, value(MP_model.first_stage_cost) f_sp); % 收敛判断 gap (UB - LB) / abs(UB); if gap eps break; end % 添加新的列与约束 u_stack{end1} u_worst; K K 1; % 最大迭代次数保护 if K K_max warning(达到最大迭代次数); break; end end %% 结果分析 plot_results(xk, u_stack, LB, UB);写出这个框架只需要理清一个问题主问题里“场景”是靠u_stack里的具体参数值来区分的子问题的“最坏场景搜索”是在不确定性集合U上做的两者在CCG循环里通过迭代交换信息——主问题把试探解x交给子问题子问题搜索出对试探解最不利的场景回传给主问题主问题拿到新场景后必须“调整方案保证该场景也安全”如此周而复始。我在实际做项目时的体会是两阶段鲁棒优化最大的门槛不在理论推导而在“把模型与算法映射到具体代码”时的各种细节处理。尤其是对偶转化、大M取值、场景传递这几个环节任何一个地方出bug都会让整个迭代过程开始震荡或不收敛。希望这篇基于Matlab的完整实现笔记能帮正在做相关方向的同学少踩几个坑。如果后续有空我准备再整理一版加入线路潮流约束和网络拓扑的多节点版本那会牵扯到更多与输电阻塞相关的耦合细节届时有新的体会再继续分享。