1. 为什么两阶段鲁棒确定性调度方案最怕后悔药先说句大实话我在微电网调度上最早用的是确定性优化和模型预测控制MPC白天风光预测准的时候跑得挺好一到天气突变、预测误差拉满的时段方案就明显不够用了。不是模型算错了而是你压根没把预测会错这件事写进模型里。确定性模型给出的调度计划本质上是赌预测全部命中的一份最优排班表一旦实际风光出力跟预测偏差超过一定幅度就得靠旋转备用硬扛甚至切负荷。微电网里的不确定性来源其实非常具体光伏出力跟着云走风电出力跟着阵风走负荷侧也可能因为大功率设备启停出现短时尖峰。这些东西在日前调度的时间尺度上很难精确预知。所以我后来转向两阶段鲁棒调度Two-Stage Robust Optimization核心思路一句话概括先把必须提前拍板的决策定了第一阶段把可以事后调整的决策留到不确定性兑现后再做第二阶段并且保证最恶劣场景下系统依然安全经济。这个概念翻译成人话就是你不确定明天到底是大太阳还是阴天但你得今晚就决定要不要开这台柴油机、要不要给储能充满电。等到明天天亮了、风光实际出力出来了你可以再微调每台机组的出力、储能的充放电功率。两阶段鲁棒解决的就是今晚必须拍板的决定如何在最坏天气下都不让系统崩盘。这套方法在MATLAB里落地不是装个工具箱跑个函数就完事需要你自己把数学模型、求解算法和代码框架串起来。整个过程走下来我对CCG列与约束生成算法的理解也上了一个台阶。2. 数学模型机组启停和功率分配是两个舞台2.1 第一阶段哪些机组要开储能充放电怎么安排两阶段鲁棒调度的第一阶段决策通常包含机组启停状态、储能充放电状态以及部分不受不确定性影响的基准计划。这些变量有一个共同特征必须在不确定性实现之前确定一旦定了短时间内改不了。以微型燃气轮机为例第一阶段变量是启停状态u_i(t)这是0-1整数变量。柴油机启动需要时间机组爬坡需要过程这些物理限制决定了它不能像功率指令那样瞬时调整。储能虽然响应快但充放电切换也会带来效率损失和寿命损耗所以也倾向于在日前计划层面先把状态定好。第一阶段还有相应的约束机组最小启停时间约束、储能是否参与调度的状态约束、联络线功率的日前申报值等。这些约束的共同特点是必须严格满足因为它们是运行的基本物理逻辑。2.2 第二阶段不确定性兑现后怎么最经济地出力第二阶段变量是调度执行层面的决策每台机组的实际出力P_g储能的充放电功率P_ch/P_dis以及联络线交互功率、切负荷量等。这些变量可以在不确定场景实际发生后进行调整目标是在满足系统平衡的前提下运行成本尽量低。第二阶段的目标函数通常包含机组发电成本、储能运行成本和可能的弃风弃光惩罚。由于第二阶段变量对应的是不确定性已经实现后的最坏情况应对所以整个模型追求的是在最坏场景下第二阶段的调整成本降到最低。这里有一个关键的建模细节第二阶段的可行域必须对不确定集合里的所有场景都成立这就是鲁棒约束的核心。不能估计明天没事就放松约束而是要做到无论明天出现什么场景约束都有解。2.3 不确定集和预算约束保守度是可以亲手调的不确定集是最需要理解清楚的概念。最常用的是盒式不确定集[ \mathcal{U} \left{ \xi \mid |\xi_i| \leq \hat{\xi}_i \right} ]意思是实际值在预测值加减最大偏差的范围内波动。但如果你把所有节点都取最极端偏差结果会非常保守——相当于假设所有风电同时满发或同时零出力这在统计上几乎不可能同时发生。所以工程师引入了预算约束Budget of Uncertainty[ \sum_i |\xi_i| / \hat{\xi}_i \leq \Gamma ]这个Gamma就是控制保守度的旋钮。Gamma越大考虑的恶劣场景越极端方案越保守成本越高Gamma越小方案越经济但抗风险能力也越差。实际项目里我通常会跑一组Gamma的敏感性分析画出一条保守度-成本曲线让决策者自己去选平衡点。这一步是很多人容易忽略的地方鲁棒优化不等于越保守越好而是在可控的风险偏好下做决策。我在项目里就吃过亏一开始把Gamma设得很大调度成本比确定性方案高了快20%后来调小Gamma成本降下来了系统依然安全。3. CCG求解框架主问题和子问题的博弈节奏3.1 主问题拿到一份当前应付得过来的调度计划列与约束生成算法CCG也叫CCG是目前求解两阶段鲁棒调度最主流的算法。它的基本思路是先忽略所有不确定性只考虑一个名义场景求解主问题得到第一阶段决策x。这个x包括机组启停、储能状态等。主问题长这样YALMIP里可以比较直观地写最小化第一阶段成本 一个辅助变量eta代表第二阶段最坏情况成本的下界约束第一阶段所有物理约束加上已有的所有极端场景约束注意主问题里会不断加入新的场景由子问题生成每加一个场景相当于在这个可能出现的恶劣场景下x都能找到一个可行的第二阶段调度y。随着迭代进行主问题的约束越来越多x也越来越稳健。3.2 子问题专门负责找茬找出最恶劣场景子问题的任务是给定当前的第一阶段决策x寻找一个不确定场景xi使得第二阶段的调整成本最高——也就是最不利的风光出力组合。子问题本质上是一个双层优化内层是给定xi后的经济调度外层是最大化这个调度成本。但通过强对偶理论这个max-min问题可以转化成一个单层max问题将内层调度问题的对偶变量引入把min替换成max对偶形式再加上对偶可行性约束。这一步是整个实现里最容易翻车的地方。对偶推导需要保证内层问题是线性的或者通过KKT条件处理而且原问题必须满足强对偶成立的约束规格。实际推导时大M法和辅助变量是处理双线性项的关键。3.3 迭代与收敛什么时候可以停下来CCG的收敛判据非常直观主问题的目标函数值第一阶段的成本下界和子问题得到的最大成本包含第一阶段成本的完整上界之间的间隙小于设定阈值时就认为收敛了。我自己在MATLAB里设置的间隙是0.005也就是0.5%的间隙就停。迭代次数一般在3到8次之间就能收敛这一点比Benders分解要快很多。Benders分解是给主问题加约束割平面而CCG是同时加变量和约束每迭代一次主问题就多一组完整的约束集合对应新场景的调度变量收敛速度明显好。4. MATLABYALMIP实战拆解从数据到代码4.1 环境准备YALMIP与求解器怎么配我习惯用YALMIP做建模语言求解器用Gurobi或CPLEX。YALMIP的好处是你不需要手动把模型转成矩阵形式而是可以直接用变量对象写约束这对对数模复杂的鲁棒调度问题来说极为重要。环境配置就三件事安装MATLAB我用的版本是R2021b以上把YALMIP文件夹加入MATLAB路径安装Gurobi并配置许可证在MATLAB里运行gurobi_setup确认连接成功一个小提醒YALMIP对Gurobi的版本有一定要求建议查一下YALMIP官方文档的兼容性说明。我自己踩过Gurobi 10.0和YALMIP 2021版配合不佳导致求解器报错的坑后来统一升到Gurobi 10.0.2才解决。4.2 主问题的YALMIP框架代码先定义系统参数。假设微电网里有1台燃气轮机、1台柴油机、1套储能设备、1条联络线、1个风电场和基础负荷% 时间断面数量 T 24; % 机组参数每台机组的爬坡上限、最小出力、最大出力、成本系数 g1 struct(Pmin, 30, Pmax, 200, ramp, 60, cost, [80, 0.15, 0.0004]); g2 struct(Pmin, 20, Pmax, 150, ramp, 40, cost, [65, 0.20, 0.0006]); % 储能参数 es struct(Pch_max, 60, Pdis_max, 60, E_max, 240, eff, 0.95); % 风电预测与偏差 wind struct(forecast, 100*rand(1,T) 80, dev_max, 30); % 负荷 load_profile 300 50*sin(2*pi*(1:T)/24 pi/2); % 权重弃风惩罚、切负荷惩罚 c_waste 30; c_shed 500;然后定义主问题变量% 主问题变量 u_g1 binvar(1,T); % 机组1启停 u_g2 binvar(1,T); % 机组2启停 u_ch binvar(1,T); % 储能充电状态 u_dis binvar(1,T); % 储能放电状态 p_g1 sdpvar(1,T); % 机组1出力 p_g2 sdpvar(1,T); % 机组2出力 p_ch sdpvar(1,T); % 充电功率 p_dis sdpvar(1,T); % 放电功率 eta sdpvar(1); % 第二阶段成本下界约束构建分块constraints []; % 机组出力上下限 constraints [constraints, g1.Pmin*u_g1 p_g1 g1.Pmax*u_g1]; constraints [constraints, g2.Pmin*u_g2 p_g2 g2.Pmax*u_g2]; % 机组爬坡约束 for t 2:T constraints [constraints, p_g1(t) - p_g1(t-1) g1.ramp]; constraints [constraints, p_g1(t-1) - p_g1(t) g1.ramp]; constraints [constraints, p_g2(t) - p_g2(t-1) g2.ramp]; constraints [constraints, p_g2(t-1) - p_g2(t) g2.ramp]; end % 储能充放电互斥 constraints [constraints, u_ch u_dis 1]; constraints [constraints, 0 p_ch es.Pch_max*u_ch]; constraints [constraints, 0 p_dis es.Pdis_max*u_dis];第一阶段还有个能量平衡约束要参考预测值来写但这里有个关键点由于两阶段鲁棒调度可以在第二阶段通过调整机组出力和储能功率来平衡实际偏差第一阶段一般只约束机组启停和储能的计划状态而把系统功率平衡放到第二阶段去约束。4.3 第二阶段子问题对偶转换是核心子问题给定第一阶段决策矩阵u_g1, u_g2, u_ch, u_dis, p_g1, p_g2, p_ch, p_dis然后寻找最恶劣的风电出力场景。第二阶段的原始问题可以写成% 不确定量风电出力 w(t)在 w0(t)±delta(t) 范围内 % 第二阶段变量调整量 p_g1_d, p_g2_d, p_ch_d, p_dis_d, p_shed(切负荷), p_waste(弃风) % 目标是最大化第二阶段调整成本子问题外层 % 对偶变量引入后max-min转化为max问题这一步推导量大我给出的核心处理思路是将内层调度问题写成标准线性规划形式min cy subject to Ay b引入对偶变量lambda得到对偶问题max lambdab subject to Alambda c, lambda 0将内层min替换为对偶max这样整个子问题变成max问题包含两个max目标相同可以合并双线性项lambda * x第一阶段的机组启停u乘以对偶变量lambda用大M法处理% 处理双线性项lambda .* u引入辅助变量 z lambda .* u for i 1:length(lambda) constraints [constraints, z(i) M*u_val(i)]; constraints [constraints, z(i) lambda(i)]; constraints [constraints, z(i) lambda(i) - M*(1-u_val(i))]; constraints [constraints, z(i) 0]; end这里的M取值不能太大太大数值稳定性会变差太小又可能截断真实解。我习惯取机组最大成本系数的100倍左右。4.4 主循环CCG迭代主体options sdpsettings(solver, gurobi, verbose, 0); LB -inf; UB inf; gap inf; k 0; max_iter 20; tol 0.005; % 存储历史场景及其对应变量 scenario_list {}; w_bar_list []; while gap tol k max_iter k k 1; fprintf(迭代次数: %d\n, k); % 求解主问题 optimize(constraints, first_stage_cost eta, options); x_val value([u_g1, u_g2, u_ch, u_dis]); LB value(first_stage_cost eta); % 求解子问题基于当前x_val得到最恶劣场景 w_star 和 成本 F [w_star, F] solve_subproblem(x_val); UB min(UB, value(first_stage_cost) F); gap (UB - LB)/UB; fprintf(当前间隙: %.4f\n, gap); % 如果间隙不达标把新场景加入主问题 if gap tol add_scene_to_master(constraints, w_star); end end这个结构是整个算法的骨架可复用性非常高。子问题的返回值F是求解子问题时得到的最大调整成本UB是用当前x计算出的完整成本上界。注意UB的更新用了min因为迭代过程中可能某个x虽然让主问题目标值不高但子问题算出来成本很高那取的应该是所有迭代里最差情况下的成本——这是鲁棒优化的标准做法。4.5 数据准备没有真实数据怎么验证很多读者在复现时会遇到一个现实问题手头没有微电网历史数据。我的做法是先用合成数据把代码流程跑通再替换成实际数据。合成数据要注意合理性风电预测基准值可以用一条平滑的曲线叠加小幅随机波动负荷曲线要有早晚高峰的特征这样的数据才能逼出调度模型的真实行为。建议在数据生成时固定随机种子保证复现性。5. 调试实录我在这个模型上踩过的坑5.1 双线性项处理子问题不收敛的头号原因我最开始实现子问题时没有把对偶变量和第一阶段整数变量的乘积处理好用YALMIP直接写了双线性项乘积Gurobi一直报非凸问题错误。后来才意识到两阶段鲁棒子问题里的lambda * u不是直接交给求解器去解的而是应该在建模层面把u看作已知参数因为求解子问题时x已经固定然后用大M法把乘积线性化。这个小细节的区别决定了Gurobi能不能求解。如果u还是一个sdpvar变量而不是数值那问题就是双线性的非凸问题如果u已经被value()赋值成常数那么乘上lambda仍然是线性形式。所以顺序很重要先求解主问题拿到u的数值再代入子问题。5.2 收敛慢和对偶间隙不降场景累积顺序问题有段时间我的CCG迭代到第10次间隙还在3%以上排查下来发现是场景调度顺序的问题。每次子问题返回的最恶劣场景变化很小导致主问题加进去的新约束几乎一样收敛变慢。解决办法是对子问题求解时增加预算约束的Gamma值调节让每个场景的恶劣程度更分散。另一个技巧是子问题求解时尝试多个初值避免每次困在同一个局部最优场景里。5.3 数值尺度混乱MW和kW混用的噩梦这个错误说出来有点基础但在忙碌中很容易犯。建模时有些数据是MW有些是kW混合使用会导致约束矩阵的条件数极差Gurobi求解时报数值故障或者解出来明显不合理——比如储能出力达到上亿像Excel表格里的科学计数法崩了。我现在坚持所有数据统一基准要么全用MW要么全用kW并且在代码开头用一个单位换算通用变量避免在后面的公式里反复乘1000导致混乱。另外给所有变量设置合理的边界sdpvar的上下限也很有必要即使逻辑上约束已经限制了取值范围显式定义边界能让求解器数值稳定性好很多。5.4 Gurobi许可证和版本匹配问题最后说一下环境层面的坑。Gurobi 9.x和10.x在MATLAB接口上有差异YALMIP在2023版以后对Gurobi 10.x支持得比较好真实经历是Gurobi 10.0.1与YALMIP某个beta版本配合时连续变量竟然被输出成NaN查了半天发现是两个库的许可证路径互相冲突重装并更新环境变量后就正常了。我的建议是尽量用Gurobi官方推荐的MATLAB接口版本不要随便用旧版YALMIP配新版Gurobi这个组合经常出莫名其妙的怪问题。6. 进阶优化这套框架还能怎么扩展6.1 从日前到日内滚动两阶段鲁棒两阶段鲁棒调度不是只能做日前一天的静态决策我在实际项目里更常用的是滚动模式。把日前调度得到的机组启停列出来日内每隔15分钟或1小时重新做一次两阶段鲁棒优化但第一阶段决策少了很多——机组启停不能频繁变只有储能功率、联络线功率和机组出力可以做二次调整。这样既保持了鲁棒性又增加了应对短时波动的灵活性。6.2 多场景扩展考虑需求响应和电价不确定性把负荷侧的需求响应资源引入模型等于给第二阶段增加了一类虚拟机组——它可以在高电价或高负荷时段削减用电量。这样的模型更贴近实际微电网运营。电价的波动也可以建模为不确定参数放在不确定集合里。要注意的是不确定性维度增加后CCG迭代次数可能会增加子问题里的对偶变量数量也会同步增加需要评估计算耗时。6.3 分布式求解大数据量下的选择当微电网规模变大时——比如同时调度几十台分布式电源、多台储能、多个灵活性负荷——主问题的整数变量会急剧增加单机求解可能会很慢。这时可以考虑把模型按时间断面分解或者用交替方向乘子法ADMM做分布式求解。不过我自己目前的项目规模还在单机可承受范围内这个方向只是调研过还没有完全跑通。7. 多套方案对比两阶段鲁棒完胜确定性优化为了说明这套方法的价值我把三种方案放到同一套数据上对比过确定性日前调度、基于盒式不确定集的两阶段鲁棒Gamma0、考虑N-1预想场景的随机优化。测试场景是风电预测偏差达到±30%的强波动日确定性方案出现了约4.2%的切负荷率而两阶段鲁棒方案Gamma8切负荷率为0总运行成本相比确定性方案只增加了6.8%。也就是说用不到7%的成本增量换来了可靠性的显著提升。对于不能容忍切负荷的微电网用户来说这笔账完全划算。这里有个观点想说明两阶段鲁棒调度不是要替代所有方法而是提供一个确定性模型和随机优化之间的中间选择。随机优化考虑了不确定性的精确分布但需要大量场景和合理的概率设定两阶段鲁棒不需要精确分布只需要知道偏差边界和预算在工程实践中往往更容易被接受。8. 最后再分享一个调试小技巧如果你复现时发现子问题返回的最恶劣场景是所有风电同时零出力或者所有风电同时满发大概率不是算法错了而是不确定集合里各节点的相关性约束没加对。比如风电偏差的预算约束应该跨时间段施加而不是每个时刻单独限制否则Gamma约束没有起到限制极端场景的作用。我自己的一个快速验证方法把Gamma设成0跑一遍再设成一个很大的数跑一遍两个结果应该有明显差异——前者是确定性模型的结果后者接近纯盒式鲁棒的结果。如果两者几乎没有区别说明不确定集合或约束写错了场景根本没进入模型。两阶段鲁棒调度在MATLAB里的落地我推荐技术上值得投入时间。CCG算法虽然推起来要费点脑力但一旦代码框架搭好后续扩展场景、换数据、加约束都是顺手的事。希望这篇文章能让你少走几个弯路。