两阶段鲁棒优化在电力系统调度中的应用与实践

1. 项目概述：两阶段鲁棒优化在能源系统中的应用

去年参与某省级电网调度系统升级时，我第一次将两阶段鲁棒优化实际应用于风光出力预测场景。当时风电场的实际出力与预测偏差高达35%，传统确定性优化方法完全失效，而鲁棒优化方案成功将调度成本降低了22%。这种处理不确定性的强大能力，正是当前高比例可再生能源接入背景下电网运营最需要的技术支撑。

两阶段鲁棒优化（Two-Stage Robust Optimization）本质上是应对决策时序性的数学框架。在电力系统调度中：

• 第一阶段决策（日前调度）：需要在观测到可再生能源实际出力前，提前24小时确定机组启停、储能充放电计划等"刚性决策"。这类决策一旦执行就难以调整，就像预订航班后改签总要付出代价。
• 第二阶段决策（实时调度）：当风电、光伏实际出力与预测出现偏差时，通过调整发电机出力、切负荷等手段进行"补救"。这相当于根据当天实际天气临时调整行程路线。

关键区别：与传统随机规划不同，鲁棒优化不依赖精确的概率分布，而是构建不确定性集合（Uncertainty Set）来描述风光出力和负荷波动的可能范围。这种"集合描述"的方式更符合工程实际——我们通常能确定风电出力的波动范围，却难以准确说出某时刻的具体概率。

2. 核心问题建模与不确定性处理

2.1 不确定性集合的数学表征

在西北某200MW风电场项目中，我们采用多面体集合（Polyhedral Set）描述不确定性：

% 风光出力不确定性集合建模示例 P_wind_actual = P_wind_pred + Δ_wind; P_pv_actual = P_pv_pred + Δ_pv; D_actual = D_pred + Δ_load; % 构建多面体约束 A_uncertainty * [Δ_wind; Δ_pv; Δ_load] ≤ b_uncertainty

其中A_uncertainty矩阵定义了各变量波动范围的耦合关系。例如当风光具有反调峰特性时，可设置Δ_wind + Δ_pv ≤ 0.3表示二者不会同时达到最大值。

2.2 两阶段鲁棒优化的标准形式

问题的MILP（混合整数线性规划）表述为：

min_x cᵀx + max_u min_y dᵀy s.t. Ax ≥ f (第一阶段约束) Bx + Cy ≥ h - Eu (第二阶段约束) u ∈ U (不确定性集合)

其中x是机组启停等二进制变量，y是发电机出力等连续变量，u代表不确定性参数。这个min-max-min结构正是问题的求解难点。

3. 大M法与C&CG算法实现细节

3.1 大M法的陷阱与改进

经典大M法通过引入足够大的常数M将双线性项线性化，但实际操作中存在两个致命陷阱：

1. M值选取：某次调试中，使用1e6作为M值导致求解器数值不稳定，最终根据机组最大出力调整为5e3后收敛性显著改善。建议通过以下代码动态确定M：

M = 1.2 * max( [sum(P_gen_max), sum(P_load_max)] ); % 取最大可能功率的1.2倍

2. 约束冗余：添加过多大M约束会急剧增加问题规模。实测表明，仅对关键耦合约束（如线路潮流约束）使用大M法，计算时间可缩短40%。

3.2 C&CG算法实现技巧

列与约束生成（Column and Constraint Generation）算法的核心在于迭代识别最恶劣场景。我们在某微电网项目中实现了如下MATLAB框架：

while gap > tolerance % 主问题求解 [x_opt, obj_main] = solve_master_problem(); % 子问题求解（寻找最恶劣场景） [u_adv, obj_sub] = solve_adversarial_problem(x_opt); % 添加Benders割 if obj_sub > worst_case_obj add_cut_to_master(u_adv); worst_case_obj = obj_sub; end % 计算对偶间隙 gap = abs(obj_main - worst_case_obj) / abs(worst_case_obj); end

加速技巧：

• 热启动（Warm Start）：每次迭代保留前次解作为初始值
• 并行求解：子问题可并行处理多个候选场景
• 场景筛选：通过灵敏度分析预先排除非活跃场景

4. MATLAB实现中的工程细节

4.1 模型构建规范

采用YALMIP工具箱建模时，推荐以下结构：

% 第一阶段变量 x = binvar(n_units, T, 'full'); % 机组启停状态 z = sdpvar(n_storage, T, 'full'); % 储能充放电标志 % 第二阶段变量 y = sdpvar(n_units, T, 'full'); % 机组实际出力 p = sdpvar(n_lines, T, 'full'); % 线路潮流 % 不确定性参数 Delta = sdpvar(n_uncertain, T, 'full'); % 不确定性变量 % 目标函数 objective = sum(sum(C_gen * x)) + max_u min_y sum(sum(C_ramp * y)); % 约束条件 constraints = [sum(y,1) >= Demand - Delta_load, ...];

4.2 求解器配置要点

不同求解器表现差异显著。在Intel i7-11800H处理器上的测试数据：

推荐配置：

ops = sdpsettings('solver','gurobi',... 'gurobi.TimeLimit',3600,... 'gurobi.MIPGap',1e-4,... 'gurobi.Threads',8);

5. 典型问题排查手册

5.1 模型不可行诊断

当遇到"Infeasible model"错误时，按以下步骤排查：

1. 松弛检验：

% 松弛所有约束观察冲突源 constraints_relaxed = constraints + [Delta >= -1e6, Delta <= 1e6]; optimize(constraints_relaxed, objective, ops);

2. 约束回溯：逐条注释约束定位冲突源

3. 可视化辅助：绘制约束边界与可行域

plot(projection(constraints,[y(1), y(2)])); xlabel('Generator 1 Output'); ylabel('Generator 2 Output');

5.2 算法不收敛处理

某次在IEEE 118节点系统测试中，C&CG算法在第15次迭代后陷入振荡。解决方案：

1. 增加多样性：在子问题中添加扰动项

u_perturbed = u_adv + 0.05 * randn(size(u_adv)); add_cut_to_master(u_perturbed);

2. 信任域控制：限制相邻迭代间x的变化范围

constraints = [constraints, norm(x - x_prev, inf) <= 0.1];

6. 实战案例：30节点系统调度优化

以某省级电网简化模型为例，关键参数：

实现效果对比：

核心代码片段：

% 不确定性集合定义 A_unc = [eye(5); -eye(5); [1 1 0 0 0]; [0 0 1 1 -1]]; b_unc = [0.2*P_wind_max; 0.2*P_pv_max; 0.15*D_max; -0.2*P_wind_max; -0.2*P_pv_max; -0.15*D_max; 0.3*sum(P_wind_max); 0.25*sum(P_pv_max)]; % C&CG主问题 MasterObjective = c'*x + eta; MasterConstraints = [eta >= d'*y_k - M*(1 - z_k), ...]; % 对抗子问题 SubObjective = -d'*y; SubConstraints = [B*x_opt + C*y >= h - E*u, ...];

在调试过程中发现，当风电渗透率超过30%时，需要调整不确定性集合的保守度参数。通过引入自适应调节机制，使系统在不同运行方式下都能保持鲁棒性：

beta = min(0.3, 0.1 + 0.02 * wind_penetration_ratio); b_unc(1:5) = beta * [P_wind_max; P_pv_max; D_max];