接手微网容量配置项目时我一直被同一个问题反复拷问业主问你这套优化方案的底气从哪来。我用两阶段鲁棒优化算法做微网多电源容量配置在Matlab里搭了完整的代码框架前后折腾了三周才把模型、主子问题迭代和数值稳定性全部跑顺。这篇文章把当时的完整思路和代码骨架整理出来——不只是贴代码更重要的是讲清楚每一步为什么这样做。适合正在做微电网规划、储能配置或者对鲁棒优化算法感兴趣的工程师和研究生参考。1. 为什么微网容量配置必须拥抱两阶段鲁棒1.1 确定性优化的阿喀琉斯之踵先还原一下场景。常规的容量配置思路是拿历史负荷和风光出力曲线塞进一个确定性优化模型里一次求解出光伏装多少、风电装多少、柴油机和储能配多大。这类模型本身没有错问题出在输入数据的单点假设上——负荷取预测均值光伏出力取典型曲线风电出力也取期望值然后求解器给出一个在平均天气剧本下成本最低的方案。这个方案看起来完美但只要你把时间尺度拉长到项目全生命周期就会发现问题微网容量一旦落地基本20年不变而期间的天气模式、负荷增长、电价波动都是不确定的。用一个单点预测值去对抗长期不确定性本质上是在赌未来二十年的剧本不会偏离预测。我在实际项目中见过太多确定性最优解在连续阴雨天里被购电成本压垮的案例——那些看似省了百分之十几投资额的方案在最不利场景下运行一年就能把节省的钱全部亏回去。1.2 两阶段鲁棒的直观思想两阶段鲁棒优化解决的就是这个决策不可逆 环境不可知的矛盾。它的核心结构是min-max-min外层min是投资决策第一阶段内层是两个min和max的组合——在给定投资方案下考虑所有可能的不确定场景max对每个场景都做出最优运行调度min最后找到最坏场景下运行成本最低的投资方案。用开餐厅来类比最直观你决定后厨设备买多少第一阶段花一大笔钱之后每一天的备菜量、进货量则根据当天客流和菜价调整第二阶段属于日常运行。你没法控制客流和菜价但希望在最差的一天——客流爆满且菜价飙涨——餐厅依然不会亏损。鲁棒优化帮你把设备容量定在最坏一天也扛得住的水平而不是平均一天刚好够用的水平。对应到微网第一阶段变量是各电源的安装容量第二阶段变量是每个时段里柴油机出力、储能充放电功率、向主网购电功率不确定性则来自光伏出力、风电出力和负荷预测的偏差。1.3 适用场景与前置条件当然不是所有项目都得用两阶段鲁棒。如果项目负荷曲线平滑、可再生能源占比低、弃光弃风风险小那么确定性优化加灵敏度分析就够了。但以下情况我会强烈建议上两阶段鲁棒可再生能源渗透率超过40%、电网购电容量受限或电价峰谷差极大、业主对停电风险容忍度很低、或者你手里有多套历史气象数据却无法确定未来真实分布。前置条件有两个必须提前评估。第一求解器支撑。虽然模型本身规模不大但加上时间耦合储能SOC和迭代生成的场景后需要Gurobi或CPLEX这种级别的商业求解器开源的CBC和SCS在大型算例下会很吃力。第二建模功底。两阶段鲁棒的难点从来不在算法本身而在把物理系统准确翻译成数学约束尤其是储能充放电的二进制变量、不确定性集合的构造、以及子问题对偶时双线性项的处理。这三块做扎实了代码实现反而是水到渠成的事。2. 数学模型投资变量、运行变量与不确定性集合的三层设计2.1 决策变量的归属数学建模第一步是明确哪些变量属于第一阶段、哪些属于第二阶段这个归属不清后面所有代码都会乱。第一阶段变量投资层是容量决策整个项目周期内不发生改变x_pv光伏安装容量kWx_wt风电安装容量kWx_dg柴油发电机安装容量kWx_bat储能安装容量kWh第二阶段变量运行层随调度时段变化在本文采用24小时典型日的滚动调度框架P_dg(t)柴油机在时段t的出力P_pv(t)、P_wt(t)光伏、风电在时段t的实际出力受不确定性和安装容量共同约束P_bat_ch(t)、P_bat_dch(t)储能充放电功率SOC(t)储能荷电状态连接相邻时段P_grid(t)时段t向主网的购电功率这里有一个容易踩的坑储能容量和功率是两个维度一个决定能量上限一个决定功率上限。简化处理时我通常把储能容量x_bat当作能量上限而充放电功率上限按容量乘以一个倍率系数来换算例如0.5C也就是功率上限 0.5 × x_bat。这个系数要让业主确认因为不同储能厂商的倍率特性差别很大。2.2 不确定性集合到底怎么选两阶段鲁棒的性能高度依赖不确定性集合的形状集合太小则鲁棒效果不足集合太大会让解过度保守。实际工程中我推荐盒式集合加预算约束的组合 U { u | u_nom - u_dev ≤ u ≤ u_nom u_dev, Σ |u_i - u_nom,i| / u_dev,i ≤ Γ }这个设计的妙处在于预算参数Γ。它限制的是所有时段同时偏离预测值的总量而不是逐点都取极端值。否则如果允许24个时段的光伏同时比预测低20%模型求出来的解只会是光伏全装最大没有任何参考价值。Γ取值的经验规律是Γ越小决策越接近确定性优化Γ越大决策越保守。具体怎么定量选我在第6章专门讲。除了盒式集合还有椭球集合适合不确定量服从正态分布的场合和场景集合适合有典型历史场景的场合。椭球集合在求解时需要二阶锥约束对求解器要求更高我在微网规划里用得少主要因为业主更喜欢用偏差百分比这种直白的语言沟通而盒式集合的参数含义一目了然。2.3 完整模型下面给出一个能够直接移植到代码的紧凑模型。目标函数写作min_x c_inv^T x max_{u∈U} min_{y∈Ω(x,u)} c_run^T y其中c_inv是单位容量投资成本向量c_run是运行成本向量包括燃料成本、购电成本、运维成本。第一阶段的约束包括容量上下限x_min ≤ x ≤ x_max投资预算约束c_inv^T x ≤ B_inv第二阶段的运行约束Ω(x,u)包括功率平衡P_pv(t) P_wt(t) P_dg(t) P_bat_dch(t) P_grid(t) P_load(t) P_bat_ch(t)柴油机出力范围与爬坡约束P_dg_min ≤ P_dg(t) ≤ P_dg_max以及相邻时段爬坡限制光伏/风电出力受容量和不确定性共同约束P_pv(t) ≤ u_pv(t) × x_pv储能SOC递推SOC(t1) SOC(t) η_ch × P_bat_ch(t) - P_bat_dch(t) / η_dch储能容量约束SOC_min × x_bat ≤ SOC(t) ≤ SOC_max × x_bat联络线功率上限0 ≤ P_grid(t) ≤ P_grid_max这个模型的阶段耦合点在于第一阶段的容量x通过约束P_dg(t) ≤ x_dg、P_pv(t) ≤ u_pv(t)×x_pv等进入第二阶段从而把投资和运行绑定在一起。鲁棒性就体现在不确定参数u不是固定的而是在集合U内由max层主动寻找最坏取值。3. CCG迭代主问题、子问题与对偶处理的完整逻辑3.1 从min-max-min到主子问题迭代两阶段鲁棒问题不能直接用常规优化器求解因为max-min嵌套结构让模型非凸。工程中最常用的求解框架是列与约束生成算法CCG也叫Zeng-Zhao算法。它的核心思路是不一次性枚举所有不确定性场景而是让子问题主动找出最坏场景把这个场景作为新的约束补充到主问题中反复迭代直到收敛。主问题是最小化投资成本加一个辅助变量ηη代表当前已发现的所有场景中最大的运行成本。主问题形式如下 MP: min c_inv^T x η s.t. (第一阶段约束) 对每个已发现的场景u_i (第二阶段约束针对场景i)且 η ≥ c_run^T y_i子问题则是在给定主问题解x后寻找使运行成本最大的不确定性场景 SP: Q(x) max_{u∈U} min_{y∈Ω(x*,u)} c_run^T yCCG的收敛性已经在理论上被证明当不确定性集合U是多面体时算法有限步内收敛到全局最优。这对我们来说是个好消息意味着不需要担心迭代过程陷入局部最优。3.2 子问题的对偶变换与双线性项处理实现子问题的关键难点在于嵌套结构。对于固定的x*内层是一个以y为决策变量的线性规划LP。LP满足强对偶条件因此可以把内层min问题替换为其对偶问题从而把max-min结构变成一个等价的max问题内层min的对偶形式为 max_{λ,μ} (b C×u D×x*)^T λ e^T μ s.t. A^T λ E^T μ d λ ≥ 0这里λ和μ是对偶乘子分别对应第二阶段的不等式约束和等式约束。合并到外层max后子问题变成 max_{u∈U, λ, μ} (b C×u D×x*)^T λ e^T μ s.t. A^T λ E^T μ dλ ≥ 0u ∈ U问题来了目标函数中出现了λ^T C×u这种双线性项两个变量相乘子问题不再是一个线性规划。处理双线性项的常规手段是McCormick包络或大M法。具体做法是先估计λ的上界λMax根据对偶可行域有限性λ确实存在上界然后引入辅助变量Z_ij λ_i × u_j用一组线性约束包络住这个乘积Z_ij ≥ 0 Z_ij ≤ u_j^max × λ_i Z_ij ≤ λ_i^max × u_j Z_ij ≥ u_j^max × λ_i λ_i^max × u_j - λ_i^max × u_j^max这组约束的逻辑很简单当u_j和λ_i都在各自的上下界内时这四个不等式恰好把乘积的取值框在可行范围内。需要注意的是这个线性化只有在λ和u都有界时才成立因此需要对λ_i的上界做预估。我在实际项目中通常先用一个较小的上界试算再检查求解结果中是否有对偶变量接近边界若有则放大上界重新求解。3.3 CCG与Benders分解的取舍很多教材会拿CCG和Benders分解做对比。直观说Benders分解是把子问题的对偶信息割平面逐步加入主问题而CCG是直接把最坏场景对应的完整变量和约束加入主问题。CCG的优势在于收敛快很多尤其是子问题含有整数变量时更是如此因为CCG不需要引入复杂的组合割。我在微网容量配置中优先用CCG的另一个原因是代码结构更清晰主问题就是一个参数化的场景集合子问题就是一个独立求解的max问题两个部分各自好调试出了问题也容易定位。如果模型规模很大比如源网荷储一体化规划我才会考虑Benders因为这时主问题每轮迭代求解负担更重Benders的轻量割平面可能整体更快。4. Matlab代码落地YALMIP建模与主循环实现4.1 代码架构与文件规划代码我建议按功能拆成五个文件不要把所有逻辑堆在一个脚本里否则调试时你会疯掉main_ccg.mCCG主循环负责迭代控制与结果输出load_system_data.m读取系统参数输出结构体parsolve_MP.m求解主问题solve_SP.m求解子问题返回最坏场景和运行成本build_run_constraints.m构建第二阶段的运行约束被MP和SP共用建模工具我用YALMIP因为它的语法最接近数学表达式省去大量矩阵拼接工作。求解器接口通过sdpsettings统一设置方便在Gurobi和CPLEX之间切换对比。4.2 主问题MP的建模实现主问题MP的完整YALMIP代码如下。核心逻辑是对每一个已发现的场景u_cell{i}复制一份运行变量构建对应的运行约束并通过η约束保证主问题的目标值不小于所有场景的运行成本。function [x_opt, eta, LB] solve_MP(u_cell, par) K length(u_cell); % 第一阶段容量变量 x sdpvar(par.nGen, 1, full); % 运行成本辅助变量 eta sdpvar(1, 1); % 各场景运行变量cell数组 y cell(1, K); CostRun zeros(1, K); Constraints []; % 第一阶段约束容量上下限与投资预算 Constraints [Constraints, par.xMin x par.xMax]; Constraints [Constraints, par.invCost * x par.budget]; % 对每个已发现场景构建运行约束 for i 1:K y{i} sdpvar(par.nRun, 1, full); [cons_i, cost_i] build_run_constraints(x, y{i}, u_cell{i}, par); Constraints [Constraints, cons_i]; Constraints [Constraints, eta cost_i]; end % 目标投资成本 最大运行成本由eta代理 Objective par.invCost * x eta; ops sdpsettings(solver, par.mpSolver, verbose, 0); sol optimize(Constraints, Objective, ops); x_opt value(x); eta value(eta); LB value(Objective); if sol.problem ~ 0 warning(MP求解异常: %s, sol.info); end end这里有一个细节值得强调每轮迭代主问题都会新增一组变量和约束对应新的最坏场景导致主问题规模逐渐增大。如果迭代到七八轮之后求解开始变慢可以考虑上一轮主问题的老场景约束其实已经不再起作用通过解耦或热启动来加速但这属于进阶优化新手先确保逻辑正确再考虑提速。4.3 子问题SP的对偶实现子问题的完整实现比主问题复杂因为需要处理对偶乘子和双线性项。我直接给一个经过验证的YALMIP代码框架其中Z矩阵就是McCormick线性化引入的辅助变量。注意代码中的par.lambdaMax需要根据对偶可行域估计这是整个实现中最容易出错的环节。function [u_worst, obj_sp] solve_SP(x_opt, par) % 不确定性变量 u sdpvar(par.nUnc, 1, full); % 对偶乘子不等式约束对应lambda非负等式约束对应mu自由 lambda sdpvar(par.nIneq, 1, full); mu sdpvar(par.nEq, 1, full); % 双线性项辅助变量 Z sdpvar(par.nIneq, par.nUnc, full); Constraints []; % 对偶可行性约束 Constraints [Constraints, par.A * lambda par.E * mu par.d]; Constraints [Constraints, lambda 0]; % 不确定性集合盒式预算约束 Constraints [Constraints, par.uMin u par.uMax]; Constraints [Constraints, sum(abs(u - par.uNom) ./ par.uDev) par.Gamma]; % McCormick线性化Z(i,j) lambda(i) * u(j) M par.lambdaMax; for i 1:par.nIneq for j 1:par.nUnc Constraints [Constraints, ... Z(i,j) 0, ... Z(i,j) par.uMax(j) * lambda(i), ... Z(i,j) M(i) * u(j), ... Z(i,j) par.uMax(j) * lambda(i) M(i) * u(j) - M(i) * par.uMax(j)]; end end % 对偶目标线性部分 双线性部分 ObjLinear par.b * lambda x_opt * par.D * lambda par.e * mu; ObjBilinear sum(sum(par.C .* Z)); Obj ObjLinear ObjBilinear; % YALMIP默认最小化取相反数实现最大化 ops sdpsettings(solver, par.spSolver, verbose, 0); sol optimize(Constraints, -Obj, ops); u_worst value(u); obj_sp value(Obj); if sol.problem ~ 0 warning(SP求解异常: %s, sol.info); end end代码里的par.C矩阵严格来说应该是一个与约束系数对应的张量结构具体维度取决于你的模型定义。这里为了展示核心逻辑把它简化为一个nIneq × nUnc的矩阵实际工程中需要根据对偶约束的系数结构仔细拼装。我建议你在实现时把小规模算例对偶目标函数的数值和手算结果做对比确认符号和转置方向没有错。4.4 CCG主循环与容差判断主循环代码如下。初始化时我把u_cell的第一个元素设为预测值场景这样第一次主问题求解就等价于一个确定性容量配置问题给后续迭代提供一个可行的起点。%% 两阶段鲁棒优化CCG主循环 clear; clc; close all; addpath(genpath(pwd)); par load_system_data(); par.Gamma 2; par.tol 1e-3; par.maxIter 20; % 初始化以预测场景作为第一个场景 u_cell{1} par.uNom; LB -inf; UB inf; gap inf; k 1; while gap par.tol k par.maxIter fprintf( 第 %d 次迭代 \n, k); % 求解主问题 [x_opt, eta, LB] solve_MP(u_cell, par); % 求解子问题 [u_worst, obj_sp] solve_SP(x_opt, par); % 更新上界当前投资成本 最坏运行成本 UB_candidate par.invCost * x_opt obj_sp; UB min(UB, UB_candidate); % 计算相对间隙 gap (UB - LB) / abs(UB); fprintf(LB %.2f, UB %.2f, gap %.4f\n, LB, UB, gap); % 未收敛则把最坏场景加入主问题 if gap par.tol k k 1; u_cell{k} u_worst; end end %% 输出结果 disp( 鲁棒优化结果 ); disp(最优容量配置 ); disp(x_opt); disp(总成本下界LB ); disp(LB); disp(总成本上界UB ); disp(UB);注意一个数值细节LB和UB的更新方向容易搞混。主问题是原问题的松弛问题因为η只约束了已发现场景所以主问题目标值给出的永远是一个下界而投资成本子问题最优值对应一个可行方案给出的是上界。程序里我让LB直接取主问题目标值UB取最小候选值这个逻辑和CCG理论完全一致。5. 3机微网算例收敛过程与鲁棒解对比分析5.1 测试系统与参数设置为了验证代码正确性我用一个简化的3机微网系统做测试。系统包含光伏、风电、柴油发电机、储能电池以及与主网的联络线。典型日按24个时段建模负荷曲线和风光出力预测曲线采用某园区夏季典型日数据但不透露出处以免卷入数据合规问题。核心参数如下设备容量下限容量上限单位投资成本备注光伏0 kW300 kW4200 元/kW出力上限受u_pv×x_pv约束风电0 kW200 kW6800 元/kW出力上限受u_wt×x_wt约束柴油机0 kW150 kW2800 元/kW燃料成本按实时出力计算储能0 kWh200 kWh1600 元/kWh0.5C倍率SOC范围0.1-0.9购电价格为分时电价高峰时段1.2元/kWh低谷时段0.4元/kWh。不确定性偏差设置为光伏预测上下浮动20%风电预测上下浮动30%负荷预测上下浮动10%。不确定性预算参数Γ取2表示所有时段偏差加权后总量被限制在两倍单时段偏差以内。5.2 迭代收敛过程详解运行主循环后CCG经历了5次迭代收敛相对间隙降到0.3%以下。每次迭代的上下界变化如下表迭代次数LB万元UB万元相对间隙最坏场景特征1468.2533.712.3%夜间高负荷2498.5527.15.4%光伏低出力负荷偏高3512.0522.62.0%风电低出力晚峰高负荷4518.9521.40.48%连续阴天负荷尖峰5520.2520.80.12%与第4轮接近趋于稳定这里的成本口径是全生命周期折算值即投资成本等年值加典型日运行成本折算到年。第一次迭代LB较低是因为主问题里只有一个预测场景约束很少松弛程度大。随着最坏场景不断加入主问题可行域被逐步收紧LB向上攀升UB则保持在上界候选的最小值两条曲线最终靠拢。这个单调收敛的模式是CCG的正常表现看到LB单调上升、UB基本持平说明代码逻辑正确。5.3 鲁棒决策与确定性决策的对比为了体现鲁棒优化的价值我把同样数据跑了一遍确定性优化即固定不确定性为预测值两个方案对比如下项目确定性方案鲁棒方案Γ2变化光伏容量180 kW240 kW33%风电容量90 kW130 kW44%柴油机容量60 kW80 kW33%储能容量100 kWh160 kWh60%总投资折算435万元520万元19.5%最坏日运行成本1.82万元/日1.25万元/日-31.3%平均日运行成本0.92万元/日0.98万元/日6.5%这个对比清晰地展示了两阶段鲁棒的代价与收益投资多花约19.5%但最坏情况下的运行成本下降31.3%代价是平均运行成本小幅上升6.5%。业主可以根据风险偏好选择中庸的Γ值——例如Γ1时投资约增加10%最坏日运行成本下降18%左右。这套模型跑出来的是一整条投资-保守度帕累托曲线而不是只有一个孤零零的最优点这恰恰是两阶段鲁棒相对确定性优化最有说服力的地方。6. 实战中踩过的坑Big-M、求解器与预算参数调优6.1 Big-M参数并不是越大越好子问题中McCormick线性化需要预先估计λ的上界也就是代码里的par.lambdaMax。我最初偷懒给了一个全局统一的M1e5结果子问题求解时间飙升而且几次迭代后数值明显不稳定。排查下来发现过大的M会严重破坏约束矩阵的条件数Gurobi在双重里耗费大量迭代。我的经验做法分三步第一步先解一个去掉双线性项、只保留线性部分的子问题记录λ的数值范围通常就能获得一个合理的上界第二步将这个范围的10倍设为lambdaMax第三步求解完整子问题后检查λ值是否接近lambdaMax如果触界则逐步放大到50倍再试直到没有变量触界为止。这套流程虽然多花一点时间但能避免用拍脑袋的M值去赌数值稳定性。提示McCormick线性化只有在λ和u都有界时才严格成立所以lambdaMax不能省。如果你发现即使把lambdaMax放大到很大某些约束仍然被激活那多半不是M的问题而是对偶模型本身写错了。6.2 求解时间爆炸的常见原因与优化三机微网规模不大单次MP和SP都只需秒级求解所以迭代时间完全可控。但当你把系统扩展到十几台机组、采用多典型日联合规划时主问题每轮会复制成千上万个变量和约束求解时间会指数级恶化。我总结出三个实用的优化手段第一个是场景压缩。粗糙的历史数据往往包含上千个时段我会先做场景聚类把相似时段合并为不超过5到7个典型日每个典型日再降采样为24个时段这样能在保精度的前提下大幅压缩模型规模。第二个是热启动。上一轮主问题的最优解可以作为下一轮主问题的初始可行解传入求解器显著减少分支定界的探索量。第三个是子问题并行——多典型日时各日子的最坏场景互相独立可以并行求解把total wall time压到原来的十分之一。6.3 不确定性预算参数的经验选择Γ是这套方法里最需要跟业主沟通的参数因为它直接决定方案的保守程度。我给项目做参数标定时一般流程是先从历史数据中统计每个时段预测误差的分布取85%分位数作为u_dev然后让Γ从0到TT是时段数以0.5为步长扫描画出总投资-最坏运行成本的帕累托曲线最后和业主确认风险承受能力选择曲线上拐点位置对应的Γ值。集群算例的经验表明Γ取T的15%到25%通常是一个合理区间。太小起不到鲁棒作用太大则方案成本激增且几乎不再有意义地降低最坏运行成本。还有一个容易忽视的点Γ的取值要结合实际情况审核。比如光伏在夜间出力本来就应该为0你不能在夜间时段也给它设置20%的波动偏差否则子问题会生成一个物理上不成立的夜间光伏满发场景导致容量配置虚高。解决办法是在不确定性集合里给每个时段乘一个出力上限系数光伏夜间强制为0。我在实际运行中还发现两阶段鲁棒优化的结果对不确定性集合的形状很敏感但对求解器选择的敏感性相对较低。Gurobi和CPLEX在同一个模型上给出的最优解完全一致差异主要体现在求解时间上。真正影响结果质量的永远是模型的物理准确性和不确定性参数的取值依据。在这两点上多花时间比反复调求解器参数要有价值得多。最后再分享一个小技巧所有CCG代码里我强烈建议在第一轮迭代时手动输入一个已知的极端场景作为初始u_cell元素而不只是用预测场景。比如历史记录里真实发生过的连续阴雨负荷尖峰那天的出力曲线。这样做能让第一次主问题的解更接近最终最优解减少两三轮迭代更重要的是能帮你快速验证子问题找出的最坏场景是否合理——如果这么明显的极端场景子问题都没找到说明你的对偶实现大概率有bug。