微网调度这个活干久了你就知道真正让人头疼的从来不是那些标准参数下的优化计算而是预测不准这四个字。我最早做日前调度方案时用确定性优化把光伏、负荷预测值当成铁板钉钉的输入跑出来的计划漂亮得很——机组启停有序、购电曲线平滑、成本低得让人满意。结果第二天下午一阵云飘过来光伏出力直接从预测的100kW掉到60kW原计划全部作废最后那一小时只能高价买电加切负荷。就是从那时候起我开始认真研究两阶段鲁棒优化并在MATLAB生态里实现了CCG列与约束生成算法来解决微网优化调度问题。这篇文章就把我从建模到MATLAB实现、再到调参踩坑的完整过程写清楚希望能给正在做微网调度、能源系统优化的研究生和工程同行一些参考。1. 为什么微网调度绕不开预测误差——两阶段鲁棒优化的出发点1.1 确定性调度方案的玻璃心先说说传统确定性调度是怎么工作的。微网里通常有光伏、风电、负荷、储能、柴油机可能还连着大电网。日前调度的目标函数一般是总运行成本最小决策变量包括机组启停、储能充放电计划、与大电网的交换功率等。约束条件无非是功率平衡、机组出力上下限、爬坡约束、储能SOC约束。这些在MATLAB里写起来一点都不难YALMIP配上Cplex或者Gurobi一个下午就能搭出完善的模型。问题出在哪呢出在所有输入都是点估计值。光伏出力按预测曲线给负荷按典型日曲线给可实际运行中这些值都会波动。预测值80kW的光伏实际可能只有50kW也可能冲到110kW。负荷同样如此。当你把一组死参数扔进优化模型里得到的调度方案在其他参数环境下可以说是一碰就碎。以功率平衡为例某时段模型假设光伏出力为80kW对应的购电计划是20kW可实际光伏只有50kW那缺口30kW谁来补储能已经按计划充了电柴油机还在最低技术出力以下结果只能紧急从大电网高价买电。这部分额外成本在确定性模型里压根没有体现。所以做微网调度的人迟早要面对一个现实你必须把不确定性放进模型里而不是在求解之后再去做所谓的鲁棒性校验。后者的思路是方案先算出来再嵌套好多场景去测试发现问题了再回头调权重。这种做法在方案简单时勉强能用可一旦机组数量多、储能和多种约束交织你根本不知道该怎么调——因为模型压根没有内生地把最坏情况纳入决策。1.2 随机优化与鲁棒优化一个赌概率一个守底线面对不确定性学界和工程界主要有两条路线随机规划Stochastic Programming和鲁棒优化Robust Optimization。随机规划把不确定参数视为随机变量给每个场景配一个概率。它的优点是结果在期望意义下最优经济性好缺点是——第一你需要准确的概率分布这在实际中很难拿到第二场景树展开之后规模爆炸计算量很大第三哪怕你给了分布决策者往往更关心的是最坏情况下会不会出大事而不是平均情况下平不平滑。鲁棒优化的思路完全相反我不需要概率分布只需要知道不确定参数落在哪个集合内。只要参数在集合里变来变去我的方案都能保证可行而且目标函数被控制在一个可接受的上界内。代价是方案偏保守——因为你是在为最坏情况做防备平均运行成本通常会比随机规划高一点。但这正是微网调度需要的特性。电力系统是实时平衡的如果某个场景下功率不平衡轻则电价惩罚重则影响系统安全。所以工程上宁可多一点保守也不能让方案在常见偏差下崩掉。两阶段鲁棒优化则是将这种思路进一步延伸决策分两阶段做第一阶段先把现在必须定的变量如机组启停状态拍板第二阶段等到不确定性实际发生后再做调整如微调柴油机出力、决定切负荷量。这种先拍板、后补救的结构非常契合日前调度加实时调整的运行机制。2. min-max-min结构拆解第一阶段拍板第二阶段补救2.1 两阶段鲁棒优化的一般数学模型两阶段鲁棒优化在数学上是一个典型的min-max-min三层结构。写成通用形式是这样的$$\min_{x}\left(c^T x \max_{u\in\mathcal{U}} \min_{y\in\Omega(x,u)} d^T y\right)$$其中x是第一阶段决策变量一般包含二进制变量柴油机启停、储能充放电状态和部分连续变量比如事先签订的购电计划。这些决策在不确定性揭晓之前就必须确定。u是不确定参数向量常见的有光伏出力、风电出力、负荷功率。U是不确定集合规定了u的波动范围。y是第二阶段决策变量包括实时调整的机组出力、储能调节量、弃风弃光量、切负荷量等。Ω(x,u)是第二阶段可行域它同时受第一阶段决策x和实际发生的不确定参数u影响。简单理解就是第一阶段以最小化总成本为目标做决策但无法预知u会取什么值。u是自然界或者说是对手选的它专门挑让系统运行成本最大的值来搞你——这就是max那一层。等u真的出现了你再以最小的调整成本去应对——这就是内层min。在微网调度中第一阶段成本c^T x通常包括机组启停成本和计划购电成本第二阶段成本d^T y包括实时购电费用、燃气消耗费用、切负荷惩罚、弃风弃光惩罚等。把这些都放进一个模型里就不是传统MILP能直接求解的了。2.2 不确定性集合的三种建模方式既然鲁棒优化放弃了概率分布那它靠什么描述不确定性靠集合。集合的形状和大小直接决定了方案的保守程度。最基础的是盒式Box集合$$\mathcal{U} {u: u_{\min} \leq u \leq u_{\max}}$$它规定每个不确定参数都在一个区间内变动。盒式集合实现简单但非常保守——它允许所有参数同时取到最坏值。在微网里这意味着光伏出力跌到下限的同时负荷还冲到上限所有倒霉事一块儿发生。现实中这种极端组合发生的概率其实很低按它做方案会过度保守。所以工程上更常用的是带预算约束的盒式集合也叫Budget集合$$\mathcal{U} \left{u: u_i^{\text{nom}} - \hat{u}_i \leq u_i \leq u_i^{\text{nom}} \hat{u}i,; \sum{i}\frac{|u_i - u_i^{\text{nom}}|}{\hat{u}_i} \leq \Gamma\right}$$这里u_i^nom是预测值u_i-hat是最大偏差幅度Γ是预算参数。它表达的意思是系统允许每个参数都在预测值附近波动但所有参数偏离预测值的总幅度不能超过Γ。你可以把Γ理解成同时倒霉的指标——Γ越大允许同时发生的坏偏差越多方案越保守Γ0时退化为确定性优化Γ等于参数个数时退化为盒式集合。我个人的工程经验是Γ取总不确定参数数量的一半左右比较平衡比如有6个不确定时段或不确定节点Γ取3~4。当然具体值要通过灵敏度分析来定后面第6节会细说。预算集合的好处是给了你和乐观/保守之间一个旋钮。这也是两阶段鲁棒在微网调度里比纯盒式集合更受欢迎的原因之一。2.3 为什么第二阶段必须是线性规划这里要强调一个非常重要的点。CCG算法要求第二阶段问题对固定的x和u是一个线性规划LP。为什么因为后面子问题的对偶变换需要用到LP对偶理论如果内层是MILP对偶就无从谈起子问题会变得极其难解。所以建模阶段就要刻意保证第二阶段变量是连续的且约束是线性的。之前有同行问我储能充放电的二进制状态能不能放到第二阶段我一般建议不要。储能状态如果必须用二进制表示就放到第一阶段去定第二阶段只做连续功率调节。万一你的模型确实需要第二阶段有整数变量那CCG的收敛是有风险的可能需要用到更复杂的嵌套分解方法那就超出本文的范围了。3. CCG算法主问题与子问题如何交替逼近3.1 为什么选CCG而不是Benders分解两阶段鲁棒问题直接没法用商业求解器求解因为它本质上有无穷多个场景——u在连续集合里取值。常见的做法是把不确定集合做离散近似但那样精度和计算量很难平衡。CCG的本质是把无穷多场景压缩成迭代中发现的关键场景通过主问题和子问题交替求解来逼近精确解。传统Benders分解也叫L-shaped方法在两阶段随机规划里非常经典。它是通过对偶变量的值生成一条割feasibility cut或optimality cut并把这个割加到主问题里。CCG和Benders最大的区别在于CCG每次识别出一个最坏场景后直接把该场景对应的第二阶段变量和约束当作新的一列加进主问题同时更新约束集。用Zeng和Zhao在2013年那篇经典论文里的说法CCG同时生成列Column和约束Constraint因此得名。实际对比下来CCG的迭代次数通常远少于Benders。我做过的微网测试案例里CCG一般5~8轮收敛Benders往往要30轮以上。这背后的原因简单想也明白Benders是通过割不断逼近第二阶段值函数而CCG是直接把真实场景塞进主问题信息更硬。3.2 主问题MP的构建与下界更新CCG迭代到第k轮时主问题MP长这样$$\min_{x, y_1, \dots, y_k}; c^T x \sum_{i1}^{k} d^T y_i$$$$\text{s.t.}; Ax \geq b$$$$My_i \leq g - Bx - Cu_i^*,\quad \forall i 1,\dots, k$$$$x \in X,; y_i \geq 0,; i1,\dots,k$$其中u_1*, u_2*, ..., u_k*是前k轮子问题返回的最坏场景。注意这里有个关键点主问题里的第二阶段变量y_i有k组每组对应一个已发现的场景。第一阶段变量x在那几组约束里是共享的。目标函数里c^T x只计一次而后面的Σ d^T y_i则是所有已发现场景下的第二阶段成本之和。主问题是原问题的一个松弛因为原问题要求对所有u∈U都可行而主问题只要求对已发现的k个场景可行。所以主问题的最优值一定≤原问题的真正最优值它是全局下界LB。随着迭代次数k增加主问题里的场景越来越多约束越来越紧下界会单调上升。3.3 子问题SP与上界更新子问题要回答两件事给定第一阶段的解x*最坏情况下第二阶段运行成本是多少以及那个最坏场景u*具体长什么样子问题固定x x*写成$$\max_{u\in\mathcal{U}} \min_{y\in\Omega(x^*,u)} d^T y$$$$\Omega(x^,u) {y: My \leq g - Bx^- Cu,; y \geq 0}$$这个max-min问题直接求解很麻烦需要做对偶变换下一节专门讲。求解完子问题后最坏场景u和对应的最坏成本SP(x)就都有了。上界UB的计算必须特别注意子问题算出来的只是第二阶段成本真正的总成本是c^T x* SP(x*)。因为x*是从主问题里取的它本身是可行解所以这个总成本是原问题的可行方案成本构成原问题的上界。每次迭代后比较UB和LB当|UB - LB| / UB ≤ ε一般取1e-3或1e-4时算法收敛返回当前最优解x*。完整流程走一遍就是初始化LB-∞UB∞k1随便给一个初始场景u_1*比如预测值。求解主问题MP得到x*和主问题目标值obj_MP更新LB obj_MP。固定x*求解子问题SP得到最坏场景u_new和第二阶段成本SP_value更新UB min(UB, c^T x* SP_value)。若(UB-LB)/UB ≤ ε终止否则把u_new作为新场景加入主问题kk1回到第2步。这个循环结构我在MATLAB里实现了无数次看似简单但每一步都有隐蔽的坑第6节专门展开讲。4. 子问题求解的硬骨头对偶、双线性项与线性化4.1 从max-min到单层优化子问题是个max-min结构直接枚举u不可行连续集合直接交给求解器也不行。核心解法是拿内层min做对偶。内层是对y的线性规划$$\min_{y}; d^T y \quad \text{s.t.}; My \leq g - Bx^* - Cu,; y \geq 0$$它等价于对偶问题$$\max_{\lambda \geq 0}; \lambda^T (g - Bx^* - Cu) \quad \text{s.t.}; M^T \lambda d$$λ是对偶变量向量。把外层的max和这个对偶问题合在一起子问题变成$$\max_{u\in\mathcal{U},; \lambda \geq 0}; \lambda^T (g - Bx^* - Cu)$$$$\text{s.t.}; M^T \lambda d$$展开目标函数就是$$\lambda^T g - \lambda^T Bx^* - \lambda^T C u$$三项分别是线性的、线性的和双线性的。前两项都好办麻烦的是λ^T C u——λ和u都是变量这俩乘在一起问题就不再是线性规划了。4.2 双线性项的三种线性化路线根据微网问题的规模和u的集合结构有三种常用处理方案第一种是顶点枚举。当不确定参数个数很少比如不超过5个且不确定集合是盒式时可以证明最坏场景一定出现在盒式集合的顶点上也就是每个不确定参数取上限或下限。你只需枚举2^n个顶点对每个顶点解一个LP取最大值。这个方法实现最简单我在建模初期和对照验证时特别喜欢用因为它是精确的而且写代码不容易错。缺点是指数爆炸n20的时候就别指望了。第二种是对偶加大M线性化。这是最通用的做法。对双线性项λ^T C u最常见的线性化手段是引入辅助变量z_{j,i}来表示u_j λ_i的乘积然后用McCormick包络或大M把乘积关系写成线性约束。比如如果u_j在区间[u_min, u_max]λ_i非负且有上界λ_i_max那z_{j,i} u_j λ_i可以近似在凸包意义上精确为z_{j,i} ≥ u_min λ_i 0 \cdot u_j - u_min \cdot 0 u_min λ_iz_{j,i} ≥ u_max λ_i λ_i_max \cdot u_j - λ_i_max \cdot u_maxz_{j,i} ≤ u_min λ_i λ_i_max \cdot u_j - λ_i_max \cdot u_minz_{j,i} ≤ u_max λ_i 0 \cdot u_j - 0 \cdot u_max u_max λ_i这里λ_i_max需要提前预估比如通过求解一个只最大化λ_i的LP给出。实际工程中McCormick包络对盒式-预算集合的处理效果不错但要注意它给出的松弛可能带来一个小gap严格讲不是精确等价。如果对精确性有执念可以引入二进制变量辅助思路是u_j只在有限个候选值里取值或者利用最优解在极端点的性质把双线性项拆成u_j取某个候选值时等于λ_i否则为零再用大M约束。这种方法在不确定集合是box-budget时可以得到精确的MILP子问题。第三种是KKT条件法。把内层min的KKT条件互补松弛写出来替换掉min问题本身于是max-min变成带互补约束的非线性问题再用大M把这些互补条件转成0-1线性约束形成MILP求解。这个方法在原理上很漂亮但实现起来比前两种复杂大M选得不好容易出现数值病态。我一般只在对偶结构特别复杂、McCormick包络撑不住的时候才考虑它。4.3 大M值的工程取值经验这里单独拉一段说大M因为太多人在这一行栽跟头。大M的本质是用一个足够大的常数强制约束成立或失效。大M太小约束该失效的时候没失效子问题会给出错误的最坏场景进而污染整个CCG迭代大M太大数值条件数变差求解器容易给出不稳定的结果。我的经验是不要拍脑袋写M1e6。先分析实际业务量纲。比如微网里功率是kW级别成本是元/kWh级别时间窗是小时级别那么大M取100000可能就够了。更严谨的做法是先对λ的每个分量单独解一个LP求其最大可能值然后乘以u的上下界差得到对应的M。虽然这会让实现代码长一截但换来的是数值稳定。我还习惯在求解完子问题后加一步回验用子问题返回的u*固定住后再解一次内层min LP看得到的成本是否等于外层max给出的值。如果不一致多半就是大M或线性化出了问题。这一步在调试期几乎是必做的。5. MATLAB实现YALMIPGurobi的CCG代码骨架5.1 环境与工具箱准备我在MATLAB里做这套东西的组合是MATLAB R2021b以上 YALMIP Gurobi或Cplex。YALMIP负责建模Gurobi负责解MILP。装好之后可以用yalmiptest命令测试一下求解器是否被正确识别。没有Gurobi的话Cplex也能跑甚至开源的HiGHS在部分纯LP场景下也能凑合但MILP性能明显弱一些。这里多说一句不只一次被问到能不能纯用MATLAB自带的linprog/intlinprog。小规模教学案例可以但CCG子问题一上对偶和大M就变成MILPintlinprog在分支定界性能上比Gurobi差太远。我最早图省事用intlinprog跑一个6节点微网单轮子问题能磨叽几十秒迭代6轮下来接近10分钟换了Gurobi之后整个流程压到几秒钟。工具选型省的时间足够你多调几轮模型。5.2 主循环的工程化封装整个CCG在MATLAB里的代码结构非常清晰。我建议把主问题和子问题各自封装成函数主循环里只负责调用、更新边界和收集场景。教学版的主循环可以按下面这个骨架写function [x_opt, obj, info] ccg_solver() tol 1e-3; LB -inf; UB inf; k 1; % 初始场景用预测值 u_list u_nominal; while true % 1. 求解主问题 [x_k, obj_MP] solve_mp(u_list); LB max(LB, obj_MP); % 2. 固定 x_k求解子问题 [u_worst, sp_cost] solve_sp(x_k); ub_candidate c_first * x_k sp_cost; if ub_candidate UB UB ub_candidate; x_best x_k; end % 3. 收敛判断 if (UB - LB) / abs(UB) tol break; end % 4. 加入新场景 u_list [u_list, u_worst]; k k 1; end x_opt x_best; obj UB; info struct(LB, LB, UB, UB, iter, k); end这个框架需要注意的是UB更新ub_candidate是总成本不是子问题sp_cost本身。这一个细节我当年写错过导致上界一路偏低判断收敛条件的时候一会儿说收敛一会儿又不收敛排查了很久才发现问题。5.3 教学版子问题盒式集合的顶点枚举如果只有3个不确定参数比如光伏、负荷、风电且不确定集合用盒式子问题采用顶点枚举是最稳妥的教学实现。function [u_worst, sp_cost] solve_sp(x_k) % 枚举盒式集合的所有顶点 n length(u_min); vertices u_min; for i 1:n vertices [vertices; ...]; % 生成2^n个顶点的完整清单 end best_cost -inf; for i 1:size(vertices, 1) u vertices(i, :); [cost, ~] solve_inner_lp(x_k, u); % 固定u后解内层LP if cost best_cost best_cost cost; u_worst u; end end sp_cost best_cost; end固定u后solve_inner_lp就是普通的LP可以用YALMIP直接建模function [cost, y_opt] solve_inner_lp(x_k, u) y sdpvar(ng, 1); constraints [M * y g - B * x_k - C * u, y 0]; ops sdpsettings(solver, gurobi, verbose, 0); optimize(constraints, d * y, ops); cost value(d * y); y_opt value(y); end这样做的好处是CCG主循环和子问题逻辑完全透明出了问题一眼就能看出来。缺点是枚举顶点数量2^n只适合不确定参数的小规模情况。5.4 进阶版子问题对偶与大M线性化当不确定参数多了以后顶点枚举撑不住就得用对偶线性化。子问题最终变成一个MILP同样用YALMIP写。关键的YALMIP代码如下function [u_worst, sp_cost] solve_sp_dual(x_k) lambda sdpvar(n_lambda, 1); u sdpvar(n_u, 1); % 双线性项需要线性化引入辅助变量z z sdpvar(n_u, n_lambda, full); % 预算约束 constraints [M_lambda * lambda d, lambda 0, ... u_min u u_max, ...]; if use_budget constraints [constraints, sum(abs(u - u_nominal) ./ u_hat) Gamma]; end % 这里再补上McCormick或大M线性化的约束 objective lambda * g - lambda * B * x_k - sum(sum(z .* C_map)); ops sdpsettings(solver, gurobi, verbose, 0); optimize(constraints, -objective, ops); u_worst value(u); sp_cost value(objective); end注意我在optimize里对目标函数取了负号因为要求max。YALMIP默认是min。对dualize之后的表达式如果你用的是dot product的形式YALMIP会自动识别双线性项并用求解器的内置处理但更可控的做法是手动线性化。McCormick线性化的约束可以按前面4.2节的公式一条条写进去代码长度不小但换来的是你对模型每一步都有掌控权。另外提一下YALMIP有内置的双线性求解器如bonmin、baron但它们不合适大规模MILP场景而且有很多数值陷阱。我做了这么多案例还是觉得手写线性化最稳。6. 实测表现与踩坑清单6.1 收敛速度三到五轮往往就能见分晓拿一个典型微网测试系统说3台柴油机、一套储能、光伏风电负荷共3个不确定参数预算Γ2时间窗24小时。不确定性集合用box-budget第二阶段变量包括机组出力调整、储能实时功率、切负荷量。我跑出来的迭代曲线大致是初始预测场景主问题给出LB大约5.2万UB是6.1万间隙约16%加入第一个最坏场景后LB升到5.6万UB降到5.9万第二轮之后间隙缩到5%以内第四轮基本就到了1%以内。也就是说对这类不太夸张的微网系统CCG在4~6轮内稳定收敛是常态。如果你遇到十几轮还不收敛优先怀疑子问题写错了而不是算法本身有问题。6.2 六个高频坑我一个个踩过坑一UB漏加第一阶段成本。这个前面强调过。子问题返回的SP(x*)只是第二阶段的成本你要算上c^T x*才是完整可行方案的成本。有些人只拿subproblem cost当UB结果UB反而比LB还低收敛判断彻底失去意义。坑二主问题里新场景的约束加重复了。每次迭代把u_worst追加进u_list后如果在写循环时用for i 1:k 建立约束的方式一旦不小心把k写成了迭代次数1就会导致同一个场景被重复添加。重复添加不会改下界但会把主问题规模白白撑大运行时间明显膨胀。建议在每次生成约束时打印一下当前k值用k这个变量写循环不要用size(u_list, 2)这种隐式写法除非你很清楚自己在干什么。坑三大M给得太小。我会在4.3节强调过大M取值。后来我养成一个习惯每次调完参数先把M打印出来验证一下确保子问题的最坏场景在回验环节能对上。坑四Gurobi默认MIPGap太大。Gurobi默认的MIPGap是1e-4看起来很小但在迭代过程中子问题的微小误差会被主问题放大。我在测试中把MIPGap设成1e-6并配合MIPFocus1偏向找可行解发现UB/LB的收敛曲线平滑了不少。代价是单次求解时间变长但迭代次数减少总体反而更快。坑五子问题不可行。微网调度模型如果不做相对完整措施假设第二阶段在某些极端u下可能找不到可行解。解决办法是在第二阶段模型里加上切负荷和弃风弃光变量并给惩罚成本。切负荷惩罚每千瓦时几百元远高于正常购电成本这样既保证可行性也不至于让模型肆无忌惮地切负荷。坑六YALMIP变量重复定义。CCG主循环每轮都会重新创建x、y等sdpvar变量如果上一次循环的变量没有清空YALMIP的符号表达式会不断叠加旧变量定义导致模型越来越大甚至内存爆炸。解决方法是每次循环内重新调用函数把变量定义都放在函数内部不要在主脚本的同一个工作区里反复使用assign。我把主问题和子问题各自封装成独立函数后这个坑再也没出现过。6.3 关于扩展方向从两阶段鲁棒到更前沿的方法两阶段鲁棒CCG并不是终点。我最近在看分布鲁棒优化DRO和自适应鲁棒优化ARO相关的工作它们和CCG的关系其实很密切。DRO是在不确定集合里再套一层概率分布的模糊集合子问题依然是max-min或max-max-min的结构CCG的框架还能继续用。ARO则假设第二阶段决策是不确定参数的一个仿射函数这样的话整个问题可以化成更小的凸优化避免迭代。但这些方法各有取舍DRO对分布信息的需求更高ARO则可能因为仿射策略限制损失最优性。如果是从实际项目的角度我建议先把本文这套两阶段鲁棒CCG吃透跑通主循环再去扩展这些变体。因为无论不确定性集合怎么改、第二阶段策略怎么设计CCG的核心循环——主问题积累场景、子问题发现最坏场景——这个思想是贯穿始终的。我在实际项目里最深的一个体会是CCG本质上是在陪你预演未来。每轮迭代找到的那个最坏场景往往能反推出系统的薄弱环节。比如某个典型的坏场景反复出现说明储能容量在那个时段就是偏紧张或者购电上限设置过低。这套方法不只是输出一个调度方案还能帮你找到系统需要加固的地方这是确定性优化完全做不到的。所以我建议读完这篇的同行别急着上复杂的大系统先拿一个小型微网模型把CCG跑通然后让程序打印每一轮发现的u*你会看到非常直观的坏场景演进过程。等你对主循环的每个边界更新都心里有数了再往里头塞更精细的设备模型——那时候你会发现所谓两阶段鲁棒其实没有想象中那么高不可攀。