两阶段鲁棒优化Matlab实现：基于CCG算法的微电网调度模型

做微电网和电力系统调度的朋友十有八九都遇到过这个场景风光出力曲线看着挺平滑但一到实际运行风电一阵风过来可能直接掉一半光伏一片云飘过来瞬间降功率负荷更是神仙难测。确定性优化算出来的调度方案放在理想场景里很好看但拿到现实中根本不敢直接执行——这就是我为什么把目光投向计及风、光、负荷不确定性的两阶段鲁棒优化。最近我用Matlab完整实现了一套两阶段鲁棒优化调度模型求解算法选的是大M法加CCG列与约束生成也就是在未知气象和负荷偏差最恶劣的情况下先做日前决策、再做日内调整把“最坏情况下的最小成本”算出来。这篇就把我的建模思路、求解细节、代码框架和踩过的坑一次性说清楚。适合正在做鲁棒优化方向的研究生以及做微电网、综合能源系统调度的工程师参考。如果说随机优化是“猜一个概率分布然后求期望”那鲁棒优化就是“假设最坏情况会发生然后保证系统扛得住”。两阶段鲁棒优化更进一步把决策拆成“现在必须定的”和“看到不确定量之后可以调整的”两部分既保留鲁棒性又不至于像传统单阶段鲁棒那样保守到成本失控。下面我从问题定位开始把整个项目拆开讲。1. 问题定位不确定性建模与两阶段决策逻辑1.1 为什么确定性优化在新能源场景下不够用传统经济调度和机组组合模型一般是确定性优化给定风电、光伏、负荷的预测曲线在满足功率平衡、机组出力上下限、爬坡等约束的前提下最小化发电成本。模型本身很成熟求解也快。但问题是预测值和实际值之间有偏差而且偏差在新能源高渗透率场景下非常大。风电的短期预测误差可能达到额定容量的20%到30%光伏更受云层影响负荷预测虽然相对准但峰谷时段偏差仍然不可忽略。如果调度方案完全按照预测值来做实际运行中功率平衡必然被打破。要么机组来不及调整导致频率越限要么被迫弃风弃光要么切负荷。所以工程上通常的做法是预留备用容量但备用容量取多少很讲究取小了不够用取大了经济性差。确定性模型没办法系统地回答“备用到底该留多少”而两阶段鲁棒优化天然把这个问题转化为“第二阶段实时调整”的范畴让模型自己算出最恶劣场景下需要多少调整能力。1.2 两阶段“先决策-后调整”架构的物理含义两阶段鲁棒优化对应到电力系统调度里物理含义非常清晰。第一阶段对应日前调度机组启停计划、机组开机状态、储能是否参与调峰这样的离散决策必须提前一天定下来因为机组启停需要时间不可能等风突然停了再开机。这些决策在不确定性实现之前做出用数学语言说就是“here-and-now”决策。第二阶段对应日内实时调度在风光负荷的实际出力偏差暴露之后在已经确定的机组启停状态基础上调整各机组实际出力、储能充放电功率、弃风弃光量和切负荷量。这些变量可以在看到不确定量之后调整叫“wait-and-see”决策。这个两阶段结构的好处在于它不会像单阶段鲁棒优化那样要求“一个方案在所有场景下都可行”而是允许第二阶段有调整动作只要调整代价可控即可。这就把保守度降下来了成本和鲁棒性之间有了一个可调节的旋钮——不确定集的预算Γ。1.3 不确定集构造从盒式集合到预算约束不确定集是鲁棒优化的灵魂。它刻画了不确定参数可能出现的范围。最常见的构造方式是盒式不确定集也就是每个不确定参数都在预测值加减一个偏差范围之内w(t) ∈ [w_pred(t) - Δw(t), w_pred(t) Δw(t)] s(t) ∈ [s_pred(t) - Δs(t), s_pred(t) Δs(t)] l(t) ∈ [l_pred(t) - Δl(t), l_pred(t) Δl(t)]如果只用一个盒式集模型会在所有不确定参数同时取最极端值时寻找最优解。但实际中风和光不会同时拉满偏差负荷也不一定正好在最坏方向这种“所有坏事一起发生”的情形概率极低。直接按盒式集做结果会非常保守成本高到没法用。所以通常要在盒式集基础上加一个预算约束限制总偏离程度∑ |w(t) - w_pred(t)| / Δw(t) ∑ |s(t) - s_pred(t)| / Δs(t) ∑ |l(t) - l_pred(t)| / Δl(t) ≤ ΓΓ就是预算控制整个调度周期内允许出现偏差的“总量”。Γ等于0时模型退化成确定性优化Γ等于所有时段数之和时模型允许所有参数在所有时段都取极端值退化成纯盒式鲁棒。实际使用中Γ一般取值的经验区间是先跑一个确定性模型再逐步增大Γ看成本曲线的拐点在哪里选取“成本增加可接受、鲁棒性明显提升”的那个值。这个不确定集的构造方式直接决定了后面CCG子问题的复杂度和求解效率所以在建模阶段就要想清楚。2. 数学模型两阶段鲁棒优化的完整约束体系2.1 第一阶段模型日前机组组合与储能计划第一阶段决策变量主要包含机组启停状态、启停动作以及预测场景下的基准出力、储能基准充放电功率。目标函数里除了常规的运行成本还要加上第二阶段最恶劣场景下的调整成本。写成数学形式就是min ∑(启动成本 空载成本 基准燃料成本) max min ∑(调整成本 弃风弃光惩罚 切负荷惩罚)第一阶段约束包括系统功率平衡约束基准场景机组出力上下限P_min * x(t) ≤ P(t) ≤ P_max * x(t)这里x(t)是启停状态0-1变量这个约束就是典型的大M逻辑约束——如果机组停机出力必须为0机组爬坡约束最小启停时间约束可选如果先期验证可以简化掉储能SOC递推约束、充放电功率上下限约束第一阶段模型中基准场景下的功率平衡可以看作“计划平衡”但真实平衡要等第二阶段拿到不确定参数后才会最终体现。前一阶段的解决定了哪些机组在线、哪些机组停机这决定了第二阶段“能够调整的空间”。2.2 第二阶段模型实时经济调整与惩罚成本第二阶段是在给定第一阶段决策和不确定参数实现值之后寻找最小调整成本的运行策略。这里的关键是引入虚拟变量来保证可行性——切负荷变量和弃风弃光变量。如果不加这两个变量在某些极端场景下第二阶段可能无解整个模型就崩了。加进去之后模型允许在极端场景下弃掉一点风光或者切掉少量负荷但惩罚系数要设得足够大让模型只在“实在没办法”的时候才用。第二阶段约束的核心是实际功率平衡∑ P_actual(t) (w(t) - curtail_w(t)) (s(t) - curtail_s(t)) P_dis(t) l(t) - LS(t) P_ch(t)其中P_actual(t)是在基准出力基础上的调整量受机组出力上下限和爬坡约束限制。第二阶段变量的可行域本质上由第一阶段决策变量x决定如果某台机组停机它的第二阶段调整量必须为0。这也是为什么第一阶段决策那么重要——它决定了第二阶段可行域的大小也就是系统的灵活性。2.3 整体模型的min-max-min表述把两个阶段合起来两阶段鲁棒优化的标准形式就是min { c^T x max min d^T y } x∈X w∈U y∈Ω(x,w)这里x表示第一阶段变量w表示不确定参数风光负荷偏差y表示第二阶段变量Ω(x,w)表示给定x和w后第二阶段的可行域。对这个结构直观理解就是先做日前决策x然后自然界挑一个最“坏”的风光负荷场景w我们再在这个场景下做最经济的日内调整y。三层结构对应的就是日前计划、最恶劣场景识别、日内再调度三个层次。整个算法要做的事情就是同时求出最优的x、最恶劣的w、以及对应的最优y。3. 求解算法CCG与Big-M法的配合逻辑3.1 主问题-子问题分解思想对比Benders分解直接求解min-max-min三层模型非常困难尤其是当x包含0-1变量、w是连续区间、y又有大量连续变量时几乎不可能用一个求解器直接端到端算出来。主流做法是将其分解为主问题和子问题迭代求解。CCG列与约束生成和Benders分解的核心思路类似主问题只决策第一阶段变量x和一个辅助变量ηη表示第二阶段成本的估计子问题在给定x后求解最恶劣场景下的第二阶段成本然后把结果反馈给主问题。区别在于Benders方法只向主问题添加一条割平面约束不添加新的决策变量而CCG每次迭代会把最恶劣场景对应的第二阶段变量和约束完整地添加到主问题里。后者添加的信息量更大所以收敛速度通常快很多尤其适合第二阶段可行域约束比较复杂的电力系统调度问题。我实测下来CCG在多数24时段算例中6到12次迭代就能收敛而Benders往往需要几十次。3.2 子问题对偶转化与双线性项处理CCG的关键难点在子问题。给定第一阶段解x*后子问题是max min d^T y w∈U y∈Ω(x*,w)内层是一个线性规划LP外层是对不确定参数w求max。如果直接用数值方法内层LP要随着w的变化反复求解没法做。标准解法是对内层LP取对偶。假设内层约束是G y ≥ h - E x* - M w对偶变量为λ那么内层LP的对偶是max (h - E x* - M w)^T λ s.t. G^T λ ≤ d λ ≥ 0把对偶问题代回子问题内层的min就变成了max于是子问题变成max (h - E x* - M w)^T λ s.t. G^T λ ≤ d λ ≥ 0 w ∈ U这里出现了一个麻烦目标函数中w和λ相乘产生双线性项w·λ。这不是线性规划普通MILP求解器不能直接处理。解决双线性项正是大M法在这个项目里最核心的用途。3.3 Big-M法线性化双线性项的关键约束与M值选择处理双线性项w·λ工程上比较实用的做法是借助大M法加二进制变量。对盒式不确定集因为目标函数关于w是线性的子问题的最优解会落在U的顶点上。也就是说最恶劣场景下每个不确定参数要么取下界、要么取上界在预算约束限制下一部分取极端值其余取预测值。于是我们可以把每个w_i表示成w_i w_i_pred δ_i * (w_i_ub - w_i_pred) δ_i ∈ {0, 1}δ_i为1表示该参数取上界为0表示取下界或保持预测值。如果考虑负荷、风电、光伏多个参数只是变量个数多一些形式不变。然后处理w_i λ_j的乘积项。令r_ij δ_i · λ_j这是0-1变量乘连续变量可以用大M法线性化0 ≤ r_ij ≤ λ_j_max · δ_i λ_j - λ_j_max · (1 - δ_i) ≤ r_ij ≤ λ_j第一条约束保证当δ_i等于0时r_ij为0第二条约束保证当δ_i等于1时r_ij等于λ_j。这样一来ritz项就被替换成了线性约束子问题转化为一个标准的MILP可以直接丢给Gurobi或者CPLEX求解。这里有一个非常实际的细节λ_j_max怎么取。理论上λ_j由对偶约束G^T λ ≤ d和λ ≥ 0界定但实际中很难直接得到一个紧的上界。我的经验是先不引入大M单独求解一个以λ_j为目标的最大化LP得到一个有限上界再放大约10%到20%作为λ_j_max。如果图省事也可以根据成本系数数量级估算比如燃料成本和惩罚成本的量级是10^3到10^5那λ_j_max取10^6左右基本不会出问题。大M值的选择是另一个坑。M太大比如直接设1e9求解器会出现数值病态表现为UB和LB来回震荡、最优解出现非物理小数M太小又可能把真正可行的解砍掉。我的调试方法是以模型中的最大成本系数或最大出力值的1到2个数量级为起点从1e4、1e5、1e6逐个试观察目标值和迭代次数的变化选一个“再调大结果也不再变化”的最小值。3.4 CCG迭代主流程整个CCG算法的主循环可以概括为以下步骤初始化设置LB为负无穷、UB为正无穷迭代次数k1设定最大迭代次数和收敛容差ε。求解主问题得到第一阶段解x_k和辅助变量η_k令LB等于主问题目标值。固定x_k求解子问题得到最恶劣场景w_k和第二阶段最优成本Q(x_k)。计算当前上界UB min(UB, c^T x_k Q(x_k))。如果UB - LB ≤ ε终止迭代输出x_k和调度方案。否则将w_k对应的第二阶段变量y_{k1}及约束添加到主问题中新增割平面约束η ≥ d^T y_{k1}k k 1回到第2步。这里需要注意主问题每迭代一次都会多出一组完整的第二阶段变量和约束所以主问题的规模会逐渐增大。但因为CCG收敛很快增加的规模在24时段微电网算例中完全可接受。4. Matlab代码实现从建模到跑通4.1 工具链选择YALMIP GurobiMatlab本身提供了intlinprog可以求解MILP但CCG迭代要反复求解主问题和子问题intlinprog在大M展开后的MILP上性能不够理想。我用的组合是YALMIP做建模语言外接Gurobi做求解器。YALMIP的语法和Matlab原生优化工具箱很接近但建模灵活度高很多动态添加约束非常方便特别适合CCG这种迭代框架。如果你没有Gurobi的license也可以用CPLEX或者SCIP。再退一步小规模算例用YALMIP加默认的intlinprog也能跑只是速度慢一些。OOMAO之类的仿真工具箱在这个场景中用不上不需要额外装。4.2 不确定集与基础参数的初始化代码模型的数据结构建议用结构体组织方便后续扩展。先定义基础参数机组台数、时段数、负荷基准值、风光预测值和偏差范围。这一部分简单但重要数据写错后面全白搭。% 基础算例参数6节点、24时段 T 24; n_G 4; % 机组数量 n_W 1; % 风电场数量 n_S 1; % 光伏电站数量 % 预测曲线示意数据实际按算例填 load_pred 100 20 * sin((1:T)/T * pi); % 负荷预测MW wind_pred 30 10 * cos((1:T)/T * 2 * pi); solar_pred max(0, 20 * sin((1:T)/T * pi)); % 偏差范围取预测值的15%~20% delta_load 0.15 * load_pred; delta_wind 0.20 * wind_pred; delta_solar 0.20 * solar_pred; % 预算约束参数 Gamma 6;不确定集合的定义用YALMIP的sdpvar来做w sdpvar(1, T, full); s sdpvar(1, T, full); l sdpvar(1, T, full); Uncertainty_Con []; for t 1:T Uncertainty_Con [Uncertainty_Con, ... w(t) wind_pred(t) - delta_wind(t), w(t) wind_pred(t) delta_wind(t), ... s(t) solar_pred(t) - delta_solar(t), s(t) solar_pred(t) delta_solar(t), ... l(t) load_pred(t) - delta_load(t), l(t) load_pred(t) delta_load(t)]; end如果需要预算约束注意它的绝对值形式。YALMIP里可以用norm函数但在大M线性化的子问题里绝对值会引入额外复杂性。更稳妥的做法是把预算约束直接写成前后两部分的偏差分解形式或者干脆先不加预算约束、只做盒式集等整套框架跑通后再加预算约束验证保守度调节。4.3 主问题代码框架主问题用YALMIP建模。第一阶段变量包括机组启停状态xbinvar和基准出力Pg0sdpvar还有一个辅助变量eta表示第二阶段成本估计。x binvar(n_G, T, full); % 机组启停 Pg0 sdpvar(n_G, T, full); % 基准出力 eta sdpvar(1, 1); Constraints_MP []; % 第一阶段约束出力上下限、爬坡、功率平衡等 for t 1:T Constraints_MP [Constraints_MP, ... sum(Pg0(:, t)) wind_pred(t) solar_pred(t) load_pred(t)]; for i 1:n_G Constraints_MP [Constraints_MP, ... Pg_min(i) * x(i, t) Pg0(i, t) Pg_max(i) * x(i, t)]; if t 1 Constraints_MP [Constraints_MP, ... Pg0(i, t) - Pg0(i, t-1) R_up(i), ... Pg0(i, t-1) - Pg0(i, t) R_down(i)]; end end end Objective_MP sum(C_fixed * x(:)) sum(C_fuel * Pg0(:)) eta;在CCG迭代中每轮会把新场景对应的第二阶段变量和约束追加进Constraints_MP。这里有一个小技巧第二阶段变量用元胞数组Y存储每组Y{k}对应第k个最恶劣场景这样主问题每次迭代只是多添一组变量和约束不会干扰之前的变量。4.4 子问题与线性化CCG中最容易写错的环节子问题的输入是第一阶段解x_fixed输出是当前最恶劣场景下的第二阶段成本Q_k和最恶劣参数w_k。这里最容易出错的地方有两个。一个是不要在子问题里保留x为优化变量应该把x的数值直接代入约束右边项。否则子问题又和第一阶段耦合在一起CCG的分解就失去了意义求解结果也会不对。第二个就是双线性项的处理。下面是核心代码思路以简单box集合为例function [Q_k, w_k, s_k, l_k] solve_subproblem(x_fixed) % 定义第二阶段变量 y机组出力调整量、弃风弃光量、切负荷量 deltaP sdpvar(n_G, T, full); curt_w sdpvar(1, T, full); curt_s sdpvar(1, T, full); LS sdpvar(1, T, full); % 不确定参数变量box集合 w_var sdpvar(1, T, full); s_var sdpvar(1, T, full); l_var sdpvar(1, T, full); % 对偶变量 lambda以及二进制变量 delta % 这里简化示意只对风电的极端取值做线性化光伏和负荷同理 delta_w binvar(1, T, full); r_w sdpvar(1, T, full); % r_w(t) delta_w(t) * lambda相关项 Constraints_SP []; % 不确定集约束省略... % 大M线性化约束 M 1e6; for t 1:T Constraints_SP [Constraints_SP, ... 0 r_w(t) M * delta_w(t), ... lambda_w(t) - M * (1 - delta_w(t)) r_w(t) lambda_w(t)]; end % 功率平衡约束、机组调整上下限约束等省略... optimize(Constraints_SP, Objective_SP, sdpsettings(solver, gurobi)); Q_k value(Objective_SP); w_k value(w_var); s_k value(s_var); l_k value(l_var); end上面是一个骨架完整实现时还需要把对偶变量的推导展开。我在实际项目里是把整段对偶约束写在子函数里再把线性化约束批量生成。这段代码写清楚之后CCG的整个循环就稳了一大半。4.5 迭代循环与结果分析主循环代码比较固定但有几个细节值得注意。首先是上下界的更新顺序。每次迭代先解主问题得到LB再解子问题得到UB顺序不要颠倒。其次是终止条件我一般设容差为1e-3或1e-4太严的意义不大因为模型本身是保守估计没必要追求小数点后三位的高精度。LB -1e9; UB 1e9; k 0; tol 1e-3; max_iter 30; while (UB - LB) / abs(UB) tol k max_iter k k 1; % 求解主问题 optimize(Constraints_MP, Objective_MP, options); LB value(Objective_MP); x_k value(x); % 求解子问题 [Q_k, w_k, s_k, l_k] solve_subproblem(x_k); UB min(UB, sum(C_fixed * x_k(:)) sum(C_fuel * Pg0_k(:)) Q_k); % 添加割平面 y_new sdpvar(...); Constraints_MP [Constraints_MP, 第二阶段约束 using w_k, s_k, l_k, ... eta 第二阶段成本表达式 using y_new]; end跑完之后重点看两个东西一是收敛曲线也就是每代的UB和LB差距是否单调收窄二是子问题返回的最恶劣场景也就是w_k、s_k、l_k的具体取值。这个最恶劣场景非常值得分析它往往不是所有参数同时取极端的“完美风暴”而是在预算约束下让系统最难受的组合。理解了这一点才算真正理解了鲁棒优化的价值。另外第二阶段结果里的弃风弃光量和切负荷量也要单独查看。如果最优解里频繁出现切负荷说明机组组合给第二阶段的调整空间不够要么增加在线机组要么增大备用容量要么调大不确定集的预算Γ让它更保守。这个反馈信号比只看总成本要直观得多。5. 调试心得与常见问题排查5.1 子问题无界或不可行的常见原因子问题求解出现unbounded是CCG项目里最让人头疼的问题。我排查时一般按三个顺序检查。第一步检查对偶约束的方向。内层是min问题约束是G y ≥ h - E x - M w那么对偶变量λ应该满足G^T λ ≤ d且λ ≥ 0。很多同学把不等号方向搞反或者把λ的范围写成无约束结果子问题目标值变成负无穷或者正无穷UB直接崩掉。第二步检查对偶转换前的原LP是否确实有界。可以先固定一个任意的w单独用linprog求解内层LP如果它本身无界说明第二阶段约束少了边界条件。常见原因是机组调整量deltaP没有加上下限或者切负荷变量LS没有限制不超过负荷值。第三步检查大M法的辅助变量有没有漏约束。r_ij δ_i λ_j如果用大M法线性化两条约束缺一不可。少了一条求解器就可能把r_ij放到任意值子问题结果完全失真。5.2 大M值过大导致的数值病态大M值不能拍脑袋乱设。我一开始图省事统一用1e9结果Gurobi频繁报数值警告UB和LB锯齿状震荡怎么调都不收敛。后来把M值降到了1e5、1e6问题立刻消失。判断M值是否合适的直观方法看主问题的最优解里有没有“微小的非物理量”。比如机组出力出现1e-7这种值或者某个0-1变量的整数解变成0.999999。这些都是数值病态的典型症状。解决方式除了调M值还可以在sdpsettings里开启数值增强选项比如gurobi.NumericFocus, 3能缓解但不治本。5.3 CCG迭代不收敛或振荡如果迭代超过30次还不收敛大概率不是算法问题而是模型写错了。我遇到过一种情况主问题添加割平面时把约束η ≥ d^T y_{k1}写成了等式η d^T y_{k1}结果主问题被绑死在某一个场景上LB永远不会上升。还有一种情况是添加第二阶段变量时y_{k1}没有和对应的最恶劣场景w_k建立约束关联。看起来割平面加上了但新变量实际上是个“自由人”不反映任何场景主问题永远得不到有效反馈。调试这类问题我的方法是在迭代循环里打印每一代的关键信息当前代k、LB、UB、Gap、最恶劣场景中哪些参数取了极端值。只要打印出来几代之内的规律一眼就能看出来。如果发现某代UB比上一代还差十有八九是前面说的“自由人”问题。5.4 典型问题速查表现象可能原因排查方向子问题unbounded对偶约束方向错误、λ范围错误检查G^T λ ≤ d和λ ≥ 0子问题infeasible第二阶段约束过紧、缺少切负荷/弃风弃光松弛变量加入虚拟松弛变量并设大惩罚UB和LB锯齿状震荡大M值过大数值病态缩小M值开启NumericFocus迭代次数异常多主问题割平面约束写成等式检查η ≥ d^T y_{k1}最恶劣场景不符合直觉预算Γ太小或不确定集约束写松了检查预算约束和区间上下界模型运行极慢大M线性化引入大量二进制变量先降时段数到6验证逻辑后再扩规模5.5 小规模测试用例的构造技巧整套代码刚开始跑的时候别直接上24时段、几十台机组的大算例。我推荐的调试路径是先用6个时段、2到3台机组、1个风电场、1个光伏站的小系统把整个CCG循环跑通确认UB/LB收敛、最恶劣场景的输出合理再逐步扩大规模。小系统的另一个好处是可以手工核对结果——比如把预算Γ设成0结果应该和确定性优化完全一致这个对拍测试能快速帮你发现建模错误。我在调试时还常用一个“暴力对照法”固定小系统后把不确定集的顶点全部枚举出来每个顶点单独跑确定性优化取最大值和CCG子问题返回的最恶劣场景成本对比。如果两者一致说明大M线性化和对偶转换是对的如果不一致一定是某个环节写错了。这个方法虽然笨很多时候比看半天公式管用。跑完这套模型我最大的体会是两阶段鲁棒优化的效果很大程度上取决于你对不确定世界的刻画。Γ取0它就是一个普通的经济调度Γ取到顶它就退回传统最坏情况鲁棒优化。中间那段才是它真正值钱的地方。另外一个小技巧CCG调试时先用小规模系统把逻辑跑通再加复杂度不然一旦数值出问题根本分不清是模型错还是求解器精度问题。希望这篇能帮到正在调试两阶段鲁棒优化的你少走点弯路。