两阶段鲁棒优化算法解决微网多电源容量配置问题理论推导看起来干净利落但真要拿Matlab落地跑通你会发现坑全藏在模型转换和迭代框架里。这篇文章我按照自己的复现经验把min-max-min三层结构、CCG列与约束生成迭代、以及最容易出错的几个工程细节全部摊开讲。适合正在做微电网规划、新能源容量配置或者被导师丢了一堆文献让你复现两阶段鲁棒模型的朋友。先说我对这套方法的总评价它解决的核心问题只有一个——“容量拍板的时候未来负荷和风光出力到底按什么取值”。确定性优化默认未来已知随机优化要求你给出概率分布而鲁棒优化只要求一个不确定集合然后直接把你拖进最坏情况里做决策。对微电网这种“晴天多赚、阴天保命”的分布式系统这套思路非常实用。1. 为什么微电网容量配置方案经常“纸上完美、现实翻车”1.1 确定性优化的先天短板我最早做微电网多电源容量配置时用的就是最传统的确定性优化。拿到一份典型日照数据、一条典型负荷曲线然后调一个足够好的求解器算出光伏、风电、储能、柴油发电机分别该装多少。结果方案从报表上看非常漂亮总成本最低、新能源渗透率接近70%、储能充放电次数合理。可只要用另一组历史数据回代验证立刻露馅——多云天光伏出力掉一半傍晚负荷尖峰时段系统就得切负荷柴油机被逼着连续满发好几天。问题出在哪确定性模型本质上是在解这样一个问题$$\min_{x,y} \quad C_{inv}(x) C_{oper}(y)$$其中所有约束里的光伏出力、负荷大小都是固定常数。它隐含的假设是未来所有时刻的出力曲线和负荷曲线跟建模时输入的那条曲线完全一样。这在真实系统里当然不成立。你把“典型”当成了“唯一”优化器自然会把所有鸡蛋放在一个篮子里专挑表观成本最低的组合完全没有给意外留备份。做容量配置不是做设备选型清单它本质上是做不确定性下的投资决策。光伏多装100kW在理想日照下是纯赚在连续阴雨天里就是闲置资产储能多装200kWh在电价波动大时能套利但负荷平稳时就是沉没成本。确定性模型看不到这些它只会在单一场景下找最优。1.2 两阶段鲁棒优化到底在优化什么两阶段鲁棒优化算法把信息结构改成了这样$$\min_{x} \left{ C_{inv}(x) \max_{u \in \mathcal{U}} \min_{y \in \Omega(x,u)} C_{oper}(y) \right}$$拆开读就是三个层次。第一层我现在拍板装多少光伏、风电、储能、柴油机这是“第一阶段决策”对应的是投资成本第二层容量定下来之后大自然会在一个你可控的不确定性集合U里挑最不利的光伏出力、负荷、电价组合这是“max”层代表最坏场景第三层最坏场景确定后调度员在这个场景下做最省钱、最安全的运行调度这是“min”层运行成本在这里产生。总目标就是“投资成本 最坏情况下的运行成本”最小化。请注意这里不是说系统一定朝最坏情况发展而是说优化结果必须保证“哪怕出现最坏情况运行层也有能力兜住”这给容量配置上了一层保险。为什么不直接用随机优化随机优化的数学形式其实更优雅但工程上有个要命的前提——你得给不确定性变量建立准确的概率分布。光伏预测误差、负荷随机波动在不同地区、不同季节的分布特性差异极大你很难获得一个可信的联合分布。鲁棒优化的门槛就低得多我只要知道光伏出力大概在预测值的上下20%范围内波动负荷在预测值的上下10%范围内波动就够了。这个集合的边界比分布容易获取得多这正是鲁棒优化在工程场景里最大的实用性优势。2. 两阶段鲁棒模型的完整数学表达2.1 第一阶段决策变量与投资成本微网多电源容量配置中最常见的电源组合是“光伏风电储能柴油发电机”再加上与配电网的联络线功率。第一阶段决策变量就是各设备的安装容量$$x [P_{PV}^{cap}, ; P_{WT}^{cap}, ; P_{ES}^{cap}, ; P_{DG}^{cap}]^T$$光伏和风电的容量单位是kW储能是kWh或kW如果同时限制功率和容量就两个变量柴油机是kW。投资成本通常按“等年值”折算把一次性建设成本平摊到设备寿命期内$$C_{inv}(x) CRF \cdot \sum_j c_j x_j$$其中$CRF \frac{r(1r)^n}{(1r)^n-1}$r是折现率n是设备寿命年c_j是单位容量投资成本。这里要说清楚一个建模细节第一阶段只需要决策容量不需要管逐时刻出力因为那是第二阶段的事。很多刚上手的朋友会在第一阶段就把储能充放电功率、柴油机启停状态都作为变量建模结果模型规模爆炸而且变量层级全乱。正确的姿势是第一阶段只有4到8个容量变量其余全部放第二阶段。2.2 第二阶段运行调度的约束体系第二阶段是给定容量x和不确定性u之后求解24小时或8760小时的最优运行策略。以典型日24个时段为例核心约束有几块。功率平衡约束是整个模型的骨架$$P_{t}^{PV} P_{t}^{WT} P_{t}^{DG} P_{t}^{grid} P_{t}^{dis} P_{t}^{load} P_{t}^{ch} P_{t}^{curt}$$左边是源右边是荷。储能充电和弃光弃风也算“荷”这个细节容易被忽略导致功率平衡永远对不上。光伏和风电的出力受到容量和不确定性共同约束$$0 \le P_{t}^{PV} \le P_{PV}^{cap} \cdot \mu_{t}^{PV} \cdot (1 \zeta_{t}^{PV})$$其中$\mu_t^{PV}$是归一化预测出力曲线$\zeta_t^{PV}$是不确定扰动。这里的约束形式是“上限约束”而不是“等式约束”允许限电这更符合实际运行逻辑。储能系统是第二个重点$$SOC_t SOC_{t-1} \eta_{ch} P_t^{ch} - \frac{P_t^{dis}}{\eta_{dis}}$$$$0.2 \le SOC_t \le 0.9, \quad 0 \le P_t^{ch}, P_t^{dis} \le P_{ES}^{cap}$$SOC上下限不能取0和1这是工程经验锂电池深度放电会加速老化建模时尽量留出安全余量。柴油机约束相对简单但也最贵$$0 \le P_t^{DG} \le P_{DG}^{cap}$$如果考虑启停状态和爬坡率就需要引入0-1变量和爬坡约束。这里提前提醒第二阶段一旦引入0-1机组组合变量子问题就不再是线性规划CCG中的对偶转换会变得非常棘手。我在第6节会专门说这个坑怎么绕。2.3 不确定性集盒子加预算的工程化处理不确定性集是两阶段鲁棒模型的灵魂它定义了“系统要扛住哪种程度的风险”。最常用的构造是“盒式预算约束”$$\mathcal{U} \left{ \zeta ;|; \zeta_t^{lo} \le \zeta_t \le \zeta_t^{hi}, ; \sum_{t} |\zeta_t| \le \Gamma \right}$$盒式约束保证了每个时刻的扰动不超过上下界预算约束Γ控制了“所有时刻同时达到极端值”的程度。为什么需要Γ因为实际系统中光伏出力不会连续24小时一直处于极端偏低状态。如果所有时刻都允许取到最差值得到的容量配置会保守到浪费钱。Γ的取值相当于在鲁棒性和经济性之间拧旋钮Γ0退化为确定性模型Γ24表示每个时刻都允许极端结果最保守。在Matlab里实现时这个集合直接写成一组线性约束即可% 光伏扰动上下界负荷扰动上下界 U_PV_lo -0.2 * P_PV_forecast; U_PV_hi 0.2 * P_PV_forecast; U_load_lo -0.1 * P_load_forecast; U_load_hi 0.1 * P_load_forecast; % 预算约束 sum(abs(zeta)) Gamma cons [cons, sum(abs(zeta_PV)) Gamma_PV, sum(abs(zeta_load)) Gamma_load];用绝对值函数建模可能会导致求解器报“非凸”建议引入辅助变量做线性化令$z_t |\zeta_t|$然后加约束$z_t \ge \zeta_t$、$z_t \ge -\zeta_t$、$z_t \le \zeta_t M(1-b)$这种或者直接用YALMIP内置的norm(zeta,1)处理预算约束。实际测试中norm函数在YALMIP里会被自动重写为MILP形式效率不差。3. 求解路径为什么用列和约束生成而不是Benders3.1 主问题与子问题的博弈结构两阶段鲁棒模型内层的“max-min”结构无法直接丢给任何商业求解器求解必须用分解算法。目前最主流的是列和约束生成CCG。它的核心思想是不要一开始就把所有不确定性场景列全而是在迭代过程中由子问题逐步“挑”出最恶劣的场景塞回主问题。到了第k轮迭代主问题长这样$$\min_{x, \theta, y_1 \dots y_k} \quad c^T x \theta$$s.t.$$\theta \ge d^T y_j, \quad j1,\dots,k$$$$F y_j \le g H u_j - E x, \quad j1,\dots,k$$$$A x \ge b, \quad y_j \ge 0$$关键就在θ和每个u_j绑定的y_j变量块。每轮子问题返回一个新的极端场景u_{k1}主问题就新增一个对应的运行变量块y_{k1}并强制θ大于等于这个场景下的运行成本。换句话说主问题每次都在问“到目前为止让我最难堪的这几批极端场景我该选什么容量才能让总成本最低”子问题则在反问“你定的这套容量在不确定性集合里到底哪个场景最致命”3.2 子问题对偶变换与双线性项处理给定第一轮主问题的解$x^*$子问题是$$\max_{u \in \mathcal{U}} \min_{y \ge 0} ; d^T y$$s.t.$$F y \le g H u - E x^*$$如果内层y全是连续变量就可以用强对偶把min变成max$$\max_{u \in \mathcal{U}, \lambda \le 0} \quad (g H u - E x^*)^T \lambda$$s.t.$$F^T \lambda \le d$$此时目标里出现了$u^T H^T \lambda$这是不确定变量u和对偶变量λ的乘积是双线性项不是凸的。处理方式有两种。如果u是连续变量用McCormick包络线性化。令$w u \cdot \lambda$由于u有界且λ在上述对偶可行域中也有界可以用四个不等式精确松弛$$w \ge u_{lo} \lambda \lambda_{lo} u - u_{lo} \lambda_{lo}$$$$w \ge u_{hi} \lambda \lambda_{hi} u - u_{hi} \lambda_{hi}$$$$w \le u_{hi} \lambda \lambda_{lo} u - u_{hi} \lambda_{lo}$$$$w \le u_{lo} \lambda \lambda_{hi} u - u_{lo} \lambda_{hi}$$如果你的不确定性变量被建模成0-1变量也就是表示“这个场景是否被选中为极端场景”那就用大M线性化更直接引入辅助变量w加约束$w \le M \lambda$、$w \ge \lambda - M(1-u)$、$w \le M u$。M的取值很有讲究我后面专门讲。经过这一通操作子问题就从一个max-min双层问题变成了一个单层MILP可以直接交给Gurobi或CPLEX求解。这一步是整个复现里最容易写错的地方数学符号稍有疏漏求解器返回的结果就会千奇百怪。3.3 CCG与Benders分解的实际差异Benders分解是另一个选择它也是主问题、子问题迭代但回传的是对偶割约束主问题不新增变量块。理论上一代经典算法实操中收敛曲线却是“锯齿状”的经常看到gap下降一点又回升。CCG每轮回传的是“完整场景”而不是“割平面信息”。主问题新增的不只是一条约束还包括一组和该场景绑定的运行变量。这样主问题的下界提升非常快。我手上的微网案例Benders跑到30多轮还在振荡CCG不到15轮就稳定收敛了。而且CCG对子问题整数变量的容忍度更高更适合工程复现。引用一个经典结论Zeng和Zhao在2013年证明当不确定性集合是有限离散集合时CCG在有限步内精确收敛对连续不确定性集合它也能在较弱的条件下收敛到最优值。工程上用的盒式预算集合虽然连续但最恶劣场景保证出现在不确定性集合的顶点上实际测试收敛非常稳定。4. Matlab代码实现与CCG迭代框架4.1 代码整体结构与建模工具选型我建议用YALMIP做建模层求解器用Gurobi或CPLEX不要直接用MATLAB自带求解器。原因很简单MILP的求解速度差距巨大微网容量配置测试一次需要迭代十几次每次都要解一个不小的MILP自带linprog和intlinprog在这种规模下性能不够看。工程上建议拆成四个文件文件职责data_def.m定义负荷曲线、光伏/风电归一化出力、设备成本、折现率、不确定性上下界等所有参数build_mp.m构建CCG主问题输入为当前所有极端场景u集合输出为YALMIP约束和变量句柄build_sp.m构建子问题输入为容量x输出为对偶化的MILP模型ccg_run.m主循环负责MP-SP迭代、gap判断、结果输出我强烈不建议把所有内容焎在一个大脚本里。因为CCG每轮都要重建MP代码耦合在一起会非常痛苦排查问题的时候你会恨不得重新写一遍。4.2 主问题建模的YALMIP实现要点MP的核心代码骨架如下function [sol_x, sol_theta] solve_mp(u_history, data) % 第一层预定义 x_cap sdpvar(4, 1); % [PV; WT; ES; DG] 容量 theta sdpvar(1, 1); % 最坏运行成本代理变量 cons []; % 容量上界约束 cons [cons, 0 x_cap data.ub_cap]; % 遍历历史极端场景每个场景对应一个 y 变量块 nT 24; y cell(length(u_history), 1); for j 1:length(u_history) uj u_history{j}; % 第二阶段运行变量 y{j} sdpvar(6, nT); % 行1:PV出力, 2:WT出力, 3:DG出力, 4:储放, 5:储充, 6:弃电 cons [cons, y{j} 0]; % 功率平衡: PVWTDGdis load ch curt cons [cons, sum(y{j}([1:3,4],:),1) data.load y{j}(5,:) y{j}(6,:)]; % 新能源能力约束 cons [cons, y{j}(1,:) x_cap(1) * data.mu_pv .* (1 uj.zeta_pv)]; cons [cons, y{j}(2,:) x_cap(2) * data.mu_wt .* (1 uj.zeta_wt)]; % 储能SOC约束在此省略用线性叠加... % theta 要大于等于该场景运行成本 cons [cons, theta data.c_oper * y{j}(:)]; end obj data.c_inv * x_cap theta; ops sdpsettings(solver, gurobi, verbose, 0); sol optimize(cons, obj, ops); sol_x value(x_cap); sol_theta value(theta); end有几个容易踩的点第一每个历史场景u_j都要有自己独立的y_j变量块不能共享y否则等于强制所有场景用同一套调度策略模型就退化成了确定性优化。第二容量变量x_cap在MP里是变量但在SP里要固定。YALMIP里用assign(x_cap, x_k)配合value(x_k)实现。第三theta 运行成本这条约束是CCG的核心它把max层变成一个可求解的θ代理。没有这条约束主问题本质上就是一个普通的多场景优化跟鲁棒优化没什么关系。4.3 子问题对偶化的建模与求解子问题在第k轮固定$x_k$构造对偶化后的MILPfunction [u_new, obj_sp] solve_sp(x_k, data) % 把第一层容量固定 assign(x_cap, x_k); nT 24; u_zeta_pv sdpvar(nT, 1); u_zeta_load sdpvar(nT, 1); u_zeta_wt sdpvar(nT, 1); % 对偶变量 lambda注意维度要与第二阶段约束行数一致 lambda sdpvar(data.n_cons_sp, 1); cons_sp []; % 对偶可行域 F*lambda d cons_sp [cons_sp, data.F * lambda data.d]; cons_sp [cons_sp, lambda 0]; % 不确定集合 cons_sp [cons_sp, data.u_lo u_zeta_pv data.u_hi]; cons_sp [cons_sp, data.u_lo u_zeta_load data.u_hi]; cons_sp [cons_sp, norm(u_zeta_pv, 1) data.Gamma_pv]; cons_sp [cons_sp, norm(u_zeta_load, 1) data.Gamma_load]; % 双线性项用McCormick线性化 % 此项要针对每个 u_i * lambda_j 构建辅助变量 w_ij % 实际代码建议写一个循环自动生成避免手打错误 obj_sp (data.g - data.E * x_k) * lambda sum_linearized_terms; ops_sp sdpsettings(solver, gurobi, verbose, 0); sol optimize(cons_sp, -obj_sp, ops_sp); % YALMIP默认最小化所以取负 u_new.zeta_pv value(u_zeta_pv); u_new.zeta_load value(u_zeta_load); obj_sp value(obj_sp); end这段代码省去了线性化辅助变量的完整展开因为不同模型的F、H矩阵结构差别很大生成辅助变量时建议用一个循环对齐索引。4.4 CCG主循环与收敛判据主循环代码大概是Kmax 30; tol 1e-3; LB -1e9; UB 1e9; k 0; % 用中性场景启动避免MP第一次就不可行 u_history{1} struct(zeta_pv, zeros(24,1), zeta_load, zeros(24,1)); while k Kmax (UB - LB) / max(abs(UB), 1e-6) tol k k 1; % 求解MP [x_k, theta_k] solve_mp(u_history, data); LB data.c_inv * x_k theta_k; % 求解SP [u_new, obj_sp] solve_sp(x_k, data); if isempty(u_new) || isnan(obj_sp) % 子问题不可行的处理一般给SP加松弛变量或可行性割 fprintf(iter %d: SP infeasible, add feasibility cut\n, k); continue; end UB min(UB, data.c_inv * x_k obj_sp); % 新场景加入历史集合 u_history{k1} u_new; fprintf(iter%d LB%.2f UB%.2f gap%.4f\n, ... k, LB, UB, (UB - LB) / abs(UB)); end收敛判据建议用相对gap不要用绝对差。因为大容量微网的投资成本动辄几百万绝对差很容易设置成过严或过松。我实际经验是tol取1e-3在工程上已经足够再往下压只是浪费时间——最后几轮gap下降极慢但容量方案的差别已经在一台小型柴油机的零头以内了。5. 典型仿真结果与参数敏感性分析5.1 鲁棒解到底在“多花多少钱”和“多扛多少风险”之间怎么走我在一组仿真数据上跑过对比假设一个峰值负荷2.5MW的孤立微网候选电源包括光伏、风电、锂电池储能和柴油机成本参数按目前中国市场主流价格估算。结果大致是这样的方案光伏(kW)风电(kW)储能(kWh)柴油机(kW)总等年成本(万元)最坏日运行成本(万元)确定性方案3100350120042064028.5鲁棒 Γ42900500160055070016.2鲁棒 Γ82700620200070076012.8鲁棒 Γ16240070024008808209.6这里有几个反直觉的点。第一鲁棒方案的光伏容量反而比确定性方案更小。原因是确定性方案默认光伏按预测曲线满发所以敢多装鲁棒方案在子问题里会主动挑“光伏连续低谷”场景光伏装太多在极端场景下不仅没用还要被计提高弃光成本。因此在不确定性下光伏的边际价值下降了而柴油机的边际价值上升了。第二储能装机的增长速度比柴油机快。因为储能在极端场景里不仅能当“电量缓冲”还能缓解柴油机爬坡压力在光伏、负荷双扰动的情况下几乎是万能解。第三总成本上升但并不是线性上升。Γ从0增加到4成本上升约10%从4增加到8再上升约9%到16时趋缓。这说明“过度保守”的边际代价在递减也就是说优化器会自动在成本和风险之间找平衡点不会无脑堆冗余。5.2 不确定性预算Γ的“魔鬼细节”Γ取值直接决定系统“敢冒多大风险”。Γ0时方案和确定性一致Γ24时每个时段都允许极端偏差储能和柴油机容量会大到离谱。但真实数据往往显示Γ取到时段数的40%左右就已经能覆盖绝大多数历史糟糕场景了。有一个容易被忽略的细节Γ是对“总偏差”的约束不是对“同时偏差数量”的约束。它给的是绝对值之和的上界所以在极端场景里优化器可能会选择“光伏连续低出力5小时”而不是“光伏、负荷同时极端波动2小时”。这两者对应的容量配置差异很大。如果你想控制同时发生的波动源数量应该把光伏、负荷、风电分成三个独立的Γ参数而不要共用一个总Γ。实操建议先用历史天气和负荷数据统计光伏出力的月均偏差计算实际“最坏三天”的偏差累积量再映射成Γ取值。这样比拍脑袋设参数可靠得多方案也不至于过度保守。6. 复现代码时最容易踩的坑6.1 子问题对偶符号不一致导致结果完全跑偏对偶变换那一步是全书最优雅也最容易出错的地方。常见线性规划原问题写成min c^T y约束F y ≤ gy ≥ 0对偶变量λ应该满足λ ≤ 0对偶约束是F^T λ ≤ d对偶目标才是(g H u - E x)^T λ。很多人随手写成λ ≥ 0结果对偶可行域整个反了目标符号也反了算法收敛到一个完全错误的容量方案上。我的检查手段很简单先固定一个已知x手动算一个最简单的单时段SP解析解然后跟代码输出的SP目标值对一下。对得上再继续对不上就先别跑CCG主要归因基本都在对偶符号上。6.2 主问题中历史场景变量块没有正确隔离CCG每轮迭代都要给主问题增加一个新的y_j块。如果你在循环里直接把新的y追加到原有的约束集合中而没有把上一轮的变量清空那每轮MP的规模会乱套甚至出现变量重复定义。YALMIP的optimize调用不会自动清理旧变量必须手动控制。建议每轮彻底重建y {}然后重新循环所有历史场景。这样代码开销稍大但逻辑清爽不容易在凌晨两点把自己绕进去。6.3 Big-M选值不当导致子问题松弛太大或数值异常如果子问题双线性项采用大M线性化M值的选取直接影响MILP的性能。M1e6这种“爽快值”会让辅助变量w的松弛空间巨大求解器在分支定界时频繁卡壳。更麻烦的是M过大时对偶变量λ的可行域在数值意义上被“撑爆”导致子问题返回一个客观存在但完全不合理的目标值。我的经验是M取值要结合具体物理量纲估算。光伏出力上界是容量×归一化功率量纲最多几百kW柴油机成本也就是每kWh几元那么λ的量纲就是元/kWhu与λ乘积的量纲是元M取10倍到100倍的该量纲数量级就足够了。能不用Big-M的地方尽量用McCormick包络它对连续u的处理更精确数值表现也更稳定。6.4 子问题不可行时到底该加什么割CCG主问题在第k轮给出容量x_k如果这个容量在某个不确定性场景下连功率平衡都满足不了子问题就会不可行。此时不能直接无视否则UB永远不会更新主问题也不知道要怎么改容量。工程上最简单的处理方法是在SP里对所有运行约束加松弛变量并给一个大惩罚系数这样SP永远不会硬性不可行。但代价是会低估最坏运行成本导致UB偏小。我自己的做法是第一遍不加松弛直接Gurobi求解如果返回状态是infeasible再调可松弛版本把罚函数项的目标值作为UB的修正参考值。这套逻辑在CCG迭代里稳定运行了很多次没有出现过发散。6.5 储能SOC跨时段约束在CCG框架里的变量归属储能SOC方程是跨时段的SOC_t关联SOC_{t-1}这导致它在每个场景的y_j块内部约束是关联的没问题。但如果你试图把SOC初始值当成第一阶段变量或者把第一阶段的储能容量同时放在SOC约束里就会让第一阶段和第二阶段变量耦合过深MP的θ割会失真。建议SOC_t全部放在第二阶段y_j块内部初始SOC设为常数0.5不要做变量。储能容量只出现在第一阶段的投资成本和第二阶段充放电功率上限约束里。这样两个阶段的边界才清晰CCG的收敛速度也会快很多。7. 留给后续扩展的思考跑通两阶段鲁棒容量配置之后我最大的体会是这只是一个很基础的起点工程上值得做的扩展还有很多。如果数据充足可以尝试分布鲁棒优化用Wasserstein距离构造一个以经验分布为中心的不确定球既能利用历史分布信息又保留了鲁棒优化“最坏情况”的思想。这类模型在三阶段微网配置里也开始有应用。再者储能设备的循环寿命衰减会影响容量配置的长期经济性。你可以把储能衰减建模成一个随充放电深度变化的退化项加入第二阶段的运行成本这样优化器会主动避免“每日满充满放”的激进策略。我实测过加入衰减项后储能最优容量会下降10%到15%但系统等年总成本反而更优。两阶段鲁棒优化本身是个框架微网容量配置只是其中一个应用。把CCG的主循环、子问题转换、割生成逻辑吃透之后你完全可以迁移到电网规划、天然气网络扩容、蓄能电站设计等更多场景。模型是骨架算法是血脉而真正让结果可信的始终是你对工程约束细节的尊重程度。希望这篇文章能帮你少走弯路早点跑通第一版鲁棒容量方案。
