简介本资源是一套面向电力系统优化方向研究生与科研人员的微电网两阶段鲁棒经济调度完整实现方案聚焦解决含风电不确定性下的调度保守性与经济性平衡问题。压缩包共13个文件1.57MB包含4个核心MATLAB脚本MP.m、SP.m、main_1.m、MP2.m、3份关键PDF文献含经典论文《微电网两阶段鲁棒优化经济调度方法》及拓展模型参考、6张算法流程与结果可视化图示jpg覆盖建模、求解、验证全流程。已有3182人学习下载是初学者切入两阶段鲁棒优化的高口碑入门范例。代码完全原创且完美复现目标函数与约束均以紧凑形式书写注释详尽特别嵌入鲁棒调节系数支持灵活调整模型保守程度便于用户快速适配自身场景或拓展至多能源协同、需求响应等新模型。 做微电网优化调度的朋友应该都有过这种经历论文里两阶段鲁棒优化的框架图画得清清楚楚CCG算法的伪代码也写得明明白白可真到自己动手在Matlab里复现时光是子问题怎么从max-min变成可求解的单层问题就够折腾一整天。更别提Yalmip、CPLEX、Gurobi这些工具之间的版本兼容问题一个不小心就是各种莫名报错最后只能对着屏幕怀疑人生。这篇文章我把整套微电网两阶段鲁棒优化经济调度方法从建模到代码实现的完整链路梳理一遍。内容涵盖两阶段鲁棒优化的数学模型、CCG列约束生成算法的迭代逻辑、MatlabYalmip环境下调用CPLEX/Gurobi的代码骨架以及复现过程中我踩过的各种坑和解决办法。主要面向正在做微电网、综合能源系统方向研究需要复现两阶段鲁棒优化算法的研究生和工程师也适合已经跑通确定性优化模型、想往鲁棒优化方向进阶的读者。1. 为什么微电网调度必须用两阶段鲁棒优化1.1 确定性优化在微电网场景下的局限微电网的优化调度经典做法是建立确定性优化模型把风电、光伏出力和负荷预测成一组确定数值然后求解一个混合整数线性规划或混合整数二次规划问题。这组确定数值求解出来的调度方案看起来最优但实际运行中只要预测有偏差方案就可能不可行——比如实际光伏出力低于预测值储能又已经放完了电系统就可能面临切负荷。实际工程中预测误差是无法消除的。哪怕是短期预测风光出力的相对误差也很容易达到15%到20%这对日前调度方案的影响是致命的。确定性优化的问题在于它把不确定性当成了不存在用单点预测值代替了整个可能的出力区间。就好比你出门前查了天气预报说今天不下雨于是空手出门结果走到一半下暴雨只能淋成落汤鸡。所以在微电网调度这个场景里问题的本质不是求一个最优解而是求一个在所有可能的不确定性实现下都可行的、并且期望成本尽可能低的解。这就是鲁棒优化登上舞台的根本原因。1.2 两阶段鲁棒优化的核心思想先决策后调整两阶段鲁棒优化的思路非常直观可以这样理解先制定一个基准调度方案第一阶段决策这个方案在不确定性还没揭示时就必须确定当不确定参数的真实值揭示后系统可以进行一系列调整操作第二阶段决策比如调整储能充放电功率、调整联络线交换功率、调整可调负荷来保证系统功率平衡和各项约束满足。这种先决策、后调整的结构非常符合微电网的实际运行逻辑。日前调度确定机组启停和基础出力日内根据实际风光出力进行调整。第一阶段决策主要是一些慢变量比如机组启停状态这些变量一旦定了很难快速改变第二阶段决策是快变量功率可以在较短时间内调整。两阶段鲁棒优化的目标是在最恶劣的不确定性实现下总成本包括第一阶段发电成本和第二阶段调整成本最小。这里最恶劣三个字是核心它不再追求概率意义上的最优期望而是追求最小化最坏情况下的成本。这个思路很符合工程保守性原则——我先保证任何情况下都不翻车再考虑省油。1.3 为什么选择CCG算法而不是Benders分解有了模型之后求解算法是关键。两阶段鲁棒优化是典型的min-max-min三层结构没法直接扔给求解器。主流求解方法有两种Benders分解和列约束生成算法CCG。Benders分解的思路是把第二阶段返回的信息以割平面的形式加到主问题中但每次迭代只加入一组约束来逼近可行域。CCG算法的核心区别在于每次迭代不仅加入割平面还引入一组新的第二阶段决策变量形成新的列。这带来的直接好处是CCG算法的上下界收敛速度远快于Benders分解。实践中的数据是对于微电网规模的算例CCG算法通常十几到二十几次迭代就能收敛而Benders分解可能需要上百次迭代。另外CCG算法生成的割平面更紧能有效减少主问题的求解规模膨胀速度。这也是为什么近年来电力系统领域的核心期刊论文基本上都采用CCG算法。如果你在论文里用Benders分解求解两阶段鲁棒优化审稿人大概率会问一句为什么不用CCG所以我这里也只讲CCG。2. 模型建立不确定集合、目标函数与约束条件的完整推导2.1 不确定集合的选取盒式集合与不确定预算两阶段鲁棒优化的第一步是定义不确定参数的取值范围。最常用的是盒式不确定集合即每个不确定参数都有明确的上下界。这个上下界可以直接从历史预测误差的统计结果中获取比如取预测误差95%置信区间。但是盒式集合有个问题如果所有不确定参数同时取到最恶劣值结果会非常保守。比如风电、光伏、负荷同时达到极端值这样的场景在现实中几乎不可能出现但模型会围绕这个极端场景做调度导致成本飙升。为了解决这个问题引入不确定预算Γ来控制同时偏离预测值的参数个数。Γ0时就是确定性优化Γ越大模型越保守Γ取不确定参数总数时就是最保守的情况。确定Γ取值的方法是遍历Γ从0到最大值计算对应的鲁棒调度成本画出成本-保守度曲线根据工程需求选择可接受的保守度水平。这个曲线后面我会专门讲怎么画。2.2 主问题与子问题的数学形式两阶段鲁棒优化的一般形式可以写成min_x c^T x max_u min_y d^T y其中x代表第一阶段决策变量机组启停状态、储能充放电状态等0-1变量u代表不确定参数风电出力、光伏出力、负荷y代表第二阶段决策变量机组出力调整量、储能功率等。约束条件分为三组第一阶段约束Ax ≤ b包括机组启停约束、设备容量约束等耦合约束Bx Cy ≤ f把第一阶段决策和第二阶段决策耦合在一起如功率平衡约束第二阶段约束Dy ≤ g(u)包含与不确定参数直接相关的约束这个min-max-min结构直接求解是不可能的。CCG算法的思路是将其分解为主问题和子问题。主问题在给定一组不确定参数场景下求解第一阶段决策和第二阶段成本之和的最小值子问题在给定第一阶段决策下寻找最恶劣的不确定参数实现和对应的第二阶段最小成本。2.3 KKT条件与大M法处理双层优化子问题是一个max-min形式。内层min关于y是线性规划可以通过强对偶定理将其转化为max形式这样内外两层都是max合并成一个max问题。对偶变换的关键步骤是引入对偶变量π将内层min问题转换成对偶max问题。转化后的子问题目标函数包含不确定变量u和对偶变量π的乘积项。由于u和π都是变量这个乘积是非线性的需要进一步处理。处理方法有两种一种是基于KKT条件的方法对第二阶段问题写出KKT条件其中互补松弛条件通过引入二进制变量和大M常数进行线性化一种是直接枚举不确定集合的顶点因为最优解一定在盒式集合的顶点处取得枚举所有顶点分别求解取最大值即可。顶点枚举在不确定参数较少时可行参数多了以后组合爆炸所以更通用的是KKT线性化方法。这里有个关键的经验大M的取值非常敏感。取值过大会导致数值稳定性问题CPLEX和Gurobi都会报出数值告警取值过小则可能把可行解错误剪掉。我的做法是先求解一个松弛模型获取变量的大致量级再在此基础上放大10到100倍作为大M值。后面第6章我会详细说这个问题。3. 环境准备MatlabYalmipCPLEX/Gurobi的安装与配置3.1 求解器选型CPLEX还是Gurobi做两阶段鲁棒优化求解器的选择直接影响求解效率和复现难度。目前主流选择是CPLEX和Gurobi。CPLEX在电力系统领域有很长的使用历史很多经典论文的算例都是用CPLEX跑的因此复现文献结果时优先用CPLEX能更好地对齐论文中的参数和结果。但CPLEX的许可证获取流程相对繁琐对个人用户不太友好社区版只支持小规模算例。Gurobi近年来的性能表现非常强劲在某些MIP问题上比CPLEX快30%以上而且学术许可证申请非常方便用学校邮箱申请后谷歌浏览器直接下载就能用。在分布式能源、微电网方向的新论文中Gurobi的使用比例逐年上升。我的建议是两套求解器都配置好代码中通过一个参数切换求解器这样既能复现老论文也能充分发挥Gurobi的性能。下面这张表是我在同一个微电网算例下的实测数据后面第5章会详细展开。3.2 安装配置中容易忽略的版本匹配问题Yalmip是个建模工具箱需要下载后添加到Matlab路径。下载渠道有两个GitHub仓库或Yalmip官网我建议直接从GitHub下载最新Release版本保证和最新版Matlab的兼容性。CPLEX和Gurobi安装时要注意安装路径不要包含中文和空格Windows下需要把求解器的bin目录添加到系统环境变量PATH中Matlab中要正确设置求解器路径Gurobi安装目录下有一个setup.m脚本运行一次即可最容易翻车的是版本匹配问题。Yalmip只是一个建模层它通过底层接口调用求解器所以Yalmip版本、求解器版本、Matlab版本三者之间必须兼容。我的经验是Matlab R2020b以上版本搭配Yalmip 2021年后发布的版本搭配CPLEX 12.10或Gurobi 9.5以上版本这个组合最稳。这里给出一个快速验证是否配置成功的代码% 验证yalmip和求解器是否就绪 yalmip(clear) x sdpvar(1,1); optimize([x 1, x 2], -x, sdpsettings(solver, cplex)); % 如果求解成功x应该等于2 value(x)如果这段代码能正常输出2说明环境没问题。如果报错优先检查求解器路径是否添加成功以及Matlab能否正确调用求解器的动态链接库。4. 代码骨架主问题-子问题迭代框架的逐步实现4.1 数据初始化与参数设置在写迭代框架前先把所有基础数据准备好。以一个包含风电、光伏、储能、微型燃气轮机和负荷的典型微电网为例需要准备以下数据各分布式电源的容量和成本参数储能的容量、充放电效率、SOC上下限不确定参数风电、光伏、负荷的预测值和波动范围微电网与大电网的联络线功率限制和购售电价格分时电价时段划分推荐把所有参数写在一个结构体中方便后续修改和扩展。比如% 基础参数定义 mpc struct(); mpc.num_gen 2; % 微型燃气轮机数量 mpc.gen_capacity [1.0, 1.5]; % 单位MW mpc.gen_ramp 0.3; % 爬坡约束 MW/h mpc.battery_capacity 2.0; % 储能容量 MWh mpc.battery_power 0.5; % 储能最大功率 MW mpc.eta_ch 0.95; % 充电效率 mpc.eta_dis 0.92; % 放电效率 mpc.T 24; % 调度时段数参数定义清楚后后续改算例只需要改这个结构体代码其他部分不用动方便做敏感性分析。4.2 主问题建模Yalmip的关键语法主问题的核心是在已知的一组不确定性实现最初给一个初始场景后续是子问题返回的最恶劣场景下求解第一阶段决策和第二阶段调整成本之和的最小值。Yalmip建模的核心是先用sdpvar定义变量然后写约束和目标函数。以下是主问题建模的核心代码框架% 主问题变量 x_start binvar(mpc.num_gen, mpc.T, full); % 机组启停状态 p_gen sdpvar(mpc.num_gen, mpc.T, full); % 机组出力 p_bat sdpvar(1, mpc.T, full); % 储能充放电功率 soc sdpvar(1, mpc.T, full); % 储能荷电状态 theta sdpvar(1, 1); % 第二阶段成本近似变量 % 主问题约束集合 Constraints []; % 功率平衡约束、机组出力上下限、储能SOC递推等... % 不确定参数在这里被固定为当前已知的最恶劣场景 u_k % 目标函数 Objective sum(sum(c_gen .* p_gen)) theta; % 求解主问题 optimize(Constraints, Objective, sdpsettings(solver, cplex, verbose, 0));这里有个重要逻辑主问题中第二阶段成本并不是直接计算的而是通过变量theta来近似。每次迭代后CCG算法会加入新的约束theta ≥ 新场景下的第二阶段成本逐步逼近真实的第二阶段成本。这个逻辑和Benders分解的割平面思想很相似但CCG更进一步会同时加入新的场景变量。4.3 子问题求解与割平面生成子问题在给定第一阶段解x*后求解最恶劣不确定性场景下的最小调整成本。经对偶变换后的子问题是一个max问题直接用Yalmip建模即可% 子问题给定x_start后求解最恶劣场景 u_wind sdpvar(1, mpc.T, full); % 风电不确定变量 u_pv sdpvar(1, mpc.T, full); % 光伏不确定变量 % 加入不确定集合约束盒式集合 不确定预算 Constraints_sub [...]; % 加入对偶可行域约束 % 目标函数取max Objective_sub [...]; optimize(Constraints_sub, -Objective_sub, sdpsettings(solver, gurobi, verbose, 0)); % 提取最恶劣场景 u_wind_star value(u_wind); u_pv_star value(u_pv); sub_cost value(Objective_sub);子问题求出的最优目标值就是当前第一阶段决策下的最恶劣成本。如果这个成本大于主问题目标中theta的当前值说明主问题的theta低估了第二阶段成本需要将当前场景加入主问题并更新上界。4.4 收敛判据与迭代逻辑CCG算法的主循环逻辑如下LB -inf; % 下界 UB inf; % 上界 Gap 1; iter 1; max_iter 30; while Gap 0.01 iter max_iter % 第1步求解主问题得到下界 solve_master_problem(); LB max(LB, objective_master); % 第2步将主问题得到的x*带入子问题求解最恶劣场景和成本 [sub_cost, worst_u] solve_subproblem(x_star); UB min(UB, sub_cost first_stage_cost); % 第3步计算相对gap Gap abs(UB - LB) / abs(UB) * 100; % 第4步将最恶劣场景加入主问题场景集合添加割平面 add_cut_to_master(worst_u, theta_val); iter iter 1; end收敛判据我用的是相对gap小于1%实际运行中通常10到15次迭代就能收敛。如果超过30次还不收敛基本可以判断是子问题或主问题建模有逻辑错误需要回头检查。要注意一个细节每次迭代主问题的规模会增加因为新增了一组场景变量和一组约束。如果迭代次数很多主问题求解时间会快速上升。所以一定要确保子问题求解正确减少无效迭代。5. CPLEX与Gurobi求解对比实测数据与调优建议5.1 相同模型下两种求解器的性能对比我在一个包含2台微型燃气轮机、1个储能系统、风电光伏负荷共计9个不确定参数的微电网算例下分别用CPLEX 12.10和Gurobi 9.5.2求解相同的两阶段鲁棒优化模型统计结果如下指标CPLEX 12.10Gurobi 9.5.2总迭代次数1818主问题单次平均求解时间4.8s3.1s子问题单次平均求解时间0.9s0.6s总求解时间102s66s最终目标函数值9846.329846.32可以看到两种求解器得到的目标函数值完全一致因为求解的都是同一个模型最优解应该是一样的。差异主要体现在求解时间上Gurobi在做MIP问题时整体比CPLEX快约35%。如果你的算例规模更大、不确定参数更多这个差距会更明显。有一点需要注意CPLEX在某些特定类型的问题上可能表现更优比如强数值病态的模型CPLEX的数值鲁棒性历史口碑更好。所以在论文复现时我会先用CPLEX跑一遍确认结果和文献一致再用Gurobi做后续的敏感性分析和大规模扩展。5.2 求解器参数调优的实战经验求解器默认参数通常不是最优的。以下是我调参后效果比较明显的几个配置对于CPLEXMIP gap容差设置为1e-4默认是1e-4但有时为了加速可以放宽到1e-3Threads设置为实际物理核心数默认是0自动在虚拟机里可能被限制关闭可选的回退计算crossover等能减少一部分求解时间对于GurobiMIPGap设置1e-4Threads同样设置为物理核心数NumericFocus设为2或3当出现数值警告时能显著改善稳定性但会增加少量求解时间TimeLimit设置一个合理的时间上限防止极端情况下的死循环在Yalmip中这些参数通过sdpsettings传入options sdpsettings(solver, gurobi, ... gurobi.MIPGap, 1e-3, ... gurobi.Threads, 8, ... gurobi.NumericFocus, 2, ... verbose, 2); optimize(Constraints, Objective, options);另一个实战经验是热启动。CCG算法每次迭代的主问题变化不大可以把上一次的解作为初始解传入能显著减少求解时间。Yalmip中实现方式是使用assign函数给变量赋初值然后求解时设置savesolveroutput和solver相关参数。6. 复现中踩过的六个坑及解决方案6.1 大M值的选择这是最让我头疼的一个问题。大M法在处理互补松弛条件时是必需的但M的取值过大会导致数值稳定性问题CPLEX和Gurobi都会报告数值警告求解结果中出现inf或NaN。我的解决方法是先算出对偶变量和第二阶段变量的量级比如通过求解松弛后的子问题获取近似范围然后取这个范围的10倍作为M的初始值如果出现数值问题再逐步减小。实际测试中M取10000时Gurobi会报数值警告改为1000后问题消失结果完全一致。所以M不是越大越好够用就行。一般来说M取到正常变量量级的10到100倍就能保证不cut掉最优解同时不会引发数值灾难。6.2 非线性项的处理子问题经对偶变换后目标函数中不确定变量u和对偶变量π的乘积是非线性的。很多人直接把这个乘积用u .* pi写进Yalmip结果求解器报错或者求解极慢。正确的做法是将其拆解为单独的变量和约束利用对偶变量的取值结构。具体来说根据不确定变量u的上下界最恶劣场景一定是在盒式集合的顶点处取得。因此可以将u的所有顶点组合枚举出来在每个顶点下求解子问题取最大值作为最恶劣场景成本。但这种方法只在不确定参数数量少时可行。我推荐的做法是引入辅助变量将乘积项线性化后再交给求解器。这在Yalmip中可以通过binvar将变量离散化或者用implies实现分段线性化。但最干净的做法还是将对偶变换推到极致直接推导出不含乘积项的等价形式。我在代码注释里详细写了推导过程这里不展开但提醒一句不要在Yalmip里直接写非线性乘积项除非你用的是非线性求解器否则必踩坑。6.3 对偶变量的维度问题第二次踩坑是对偶变量的维度对应错误。第二阶段约束中的每个约束都需要一个对偶变量而且对偶变量的数量必须和约束数量严格对应。多写一个、少写一个、顺序写反都不会直接报错但结果会完全错误。我的排查方法先把确定性版本即Γ0跑出来和直接求解一个等价的单层优化结果对比。如果对偶变量维度有问题单层优化和两阶段的结果必然对不上很快就能发现问题。对偶变量维度的检查方法是用size函数查看变量维度和约束条件数量做比对。注意Yalmip中约束条件写成矩阵形式时一个矩阵约束算一个约束块但对偶变量仍然要逐行对应这块是最容易出错的。6.4 收敛判据过严导致死循环刚开始我把gap阈值设为1e-6结果跑了40多次迭代还没收敛主问题规模越来越大单次求解时间从几秒涨到几十秒整体陷入死循环。原因是数值精度有限上下界在一个很小的区间内来回摆动无法达到1e-6的水平。后来我把阈值改为0.5%结果14次迭代就收敛了总成本只相差0.2%左右完全满足工程精度需求。建议gap阈值设置在0.5%到2%之间不需要过分追求微小gap。对审稿人来说你的gap收敛曲线和最终结果合理比gap绝对值多小更重要。6.5 Yalmip与求解器版本不兼容有段时间我把Yalmip更新到最新版结果原有代码大面积报错。原因是新版Yalmip在某些语法上做了调整比如sdpvar的默认属性、optimize的参数处理方式。解决方案是不要盲目追求最新版。找到一套稳定的组合后就固定下来。我目前的组合是Matlab R2022b Yalmip R20230623 Gurobi 10.0.1 CPLEX 12.10运行稳定了大半年期间没再出现版本类报错。如果你们实验室用的是旧版Matlab建议先把Yalmip版本锁定在对应的兼容版本不要轻易升级。6.6 数据归一化问题最后一个坑是量纲不统一。风电和光伏的预测值单位是MW成本单位是元/MWh储能SOC是百分比三者数量级差别很大。如果直接建模求解器内部的数值处理精度会受到影响特别是大M法场景下更容易出问题。我的做法是将所有变量统一到标幺值系统基准功率取微电网最大负荷成本统一到元/MWh。这样所有变量的量级都在0到几之间求解器的数值稳定性大幅提升。这个操作对结果没有影响只是换了个单位体系但对求解器来说却是天壤之别。7. 结果解读如何从鲁棒调度方案中提取有效信息7.1 最恶劣场景的识别与分析CCG算法迭代收敛后子问题返回的最后一个场景就是当前调度方案下系统面临的最恶劣不确定性实现。这个场景的工程意义非常大它告诉你系统的瓶颈在哪里。比如如果最恶劣场景中光伏出力取下限且风电出力也取下限说明系统对新能源出力的下行风险最敏感。我在实际项目中会专门把最恶劣场景提取出来做可视化和预测场景画在同一张图上直观展示鲁棒方案的设计依据。这对外审和组会汇报都非常有帮助。你还能从最恶劣场景中看出哪些时段最容易出现供需失衡为后续的储能扩容或联络线升级提供数据支持。7.2 鲁棒调节参数与系统成本的权衡曲线这是论文中最常用的一个结果图。做法很简单将不确定预算Γ从0遍历到最大值每个Γ值下运行一遍完整的CCG算法记录总成本。然后把Γ-总成本的关系画成曲线。这条曲线揭示了模型保守度和经济性之间的权衡关系。我跑过的算例中Γ从0增加到9时总成本上升了约18%。但在Γ增加到6以后成本上升的速度明显放缓曲线趋于平坦。这说明系统对不确定性的敏感度在Γ较小时最大实际工程中取Γ6左右就能获得较好的鲁棒性而不过分牺牲经济性。这张曲线的同样可以用于论文中的敏感性分析比如换不同的波动范围、不同的储能容量设置观察权衡曲线的形状变化导出一些有工程意义的结论。7.3 后续研究的扩展方向复现并跑通两阶段鲁棒优化后可以在多个方向迭代创新。比如将条件风险价值引入第二阶段目标函数在鲁棒优化框架下考虑风险偏好或者将不确定集合从盒式扩展到椭球式或多面体式更精细地刻画不确定性还可以把模型从单微网扩展到多微网互联两阶段鲁棒优化同样适用。我自己目前在做的方向是把碳交易机制引入两阶段鲁棒优化框架在第二阶段加入碳配额和碳价不确定性。碳价的波动本质上也是一种不确定性完全可以用不确定集合描述这样就能在同一个框架下同时处理新能源出力和碳价的双重不确定性。实际跑通第一版两阶段鲁棒优化代码后我的一个体会是这类算法的复现难点不在算法本身而在于建模细节和数值处理的完整性。建议先把确定性模型跑通并验证正确再改造成鲁棒优化模型这样每一步的结果都有对比排查bug会容易得多。你的第一版代码不追求效率先把结果跑对再逐步优化。最后再分享一个小技巧迭代过程中记录每一轮的上下界、gap、最恶劣场景对算法的行为做完整日志。CCG算法出了问题时这份日志可以帮你快速定位是主问题建模错了还是子问题对偶写错了还是收敛条件设置不合理。我自己是把每次迭代的日志打印成表格保存在文本文件里回头排查问题直接翻日志比对着代码猜效率高太多了。这个习惯帮我节省了大量调试时间也推荐你从第一次跑通代码时就养成。本文还有配套的精品资源点击获取 微电网两阶段鲁棒优化调度：CCG算法与Matlab实现全解析
简介本资源是一套面向电力系统优化方向研究生与科研人员的微电网两阶段鲁棒经济调度完整实现方案聚焦解决含风电不确定性下的调度保守性与经济性平衡问题。压缩包共13个文件1.57MB包含4个核心MATLAB脚本MP.m、SP.m、main_1.m、MP2.m、3份关键PDF文献含经典论文《微电网两阶段鲁棒优化经济调度方法》及拓展模型参考、6张算法流程与结果可视化图示jpg覆盖建模、求解、验证全流程。已有3182人学习下载是初学者切入两阶段鲁棒优化的高口碑入门范例。代码完全原创且完美复现目标函数与约束均以紧凑形式书写注释详尽特别嵌入鲁棒调节系数支持灵活调整模型保守程度便于用户快速适配自身场景或拓展至多能源协同、需求响应等新模型。 做微电网优化调度的朋友应该都有过这种经历论文里两阶段鲁棒优化的框架图画得清清楚楚CCG算法的伪代码也写得明明白白可真到自己动手在Matlab里复现时光是子问题怎么从max-min变成可求解的单层问题就够折腾一整天。更别提Yalmip、CPLEX、Gurobi这些工具之间的版本兼容问题一个不小心就是各种莫名报错最后只能对着屏幕怀疑人生。这篇文章我把整套微电网两阶段鲁棒优化经济调度方法从建模到代码实现的完整链路梳理一遍。内容涵盖两阶段鲁棒优化的数学模型、CCG列约束生成算法的迭代逻辑、MatlabYalmip环境下调用CPLEX/Gurobi的代码骨架以及复现过程中我踩过的各种坑和解决办法。主要面向正在做微电网、综合能源系统方向研究需要复现两阶段鲁棒优化算法的研究生和工程师也适合已经跑通确定性优化模型、想往鲁棒优化方向进阶的读者。1. 为什么微电网调度必须用两阶段鲁棒优化1.1 确定性优化在微电网场景下的局限微电网的优化调度经典做法是建立确定性优化模型把风电、光伏出力和负荷预测成一组确定数值然后求解一个混合整数线性规划或混合整数二次规划问题。这组确定数值求解出来的调度方案看起来最优但实际运行中只要预测有偏差方案就可能不可行——比如实际光伏出力低于预测值储能又已经放完了电系统就可能面临切负荷。实际工程中预测误差是无法消除的。哪怕是短期预测风光出力的相对误差也很容易达到15%到20%这对日前调度方案的影响是致命的。确定性优化的问题在于它把不确定性当成了不存在用单点预测值代替了整个可能的出力区间。就好比你出门前查了天气预报说今天不下雨于是空手出门结果走到一半下暴雨只能淋成落汤鸡。所以在微电网调度这个场景里问题的本质不是求一个最优解而是求一个在所有可能的不确定性实现下都可行的、并且期望成本尽可能低的解。这就是鲁棒优化登上舞台的根本原因。1.2 两阶段鲁棒优化的核心思想先决策后调整两阶段鲁棒优化的思路非常直观可以这样理解先制定一个基准调度方案第一阶段决策这个方案在不确定性还没揭示时就必须确定当不确定参数的真实值揭示后系统可以进行一系列调整操作第二阶段决策比如调整储能充放电功率、调整联络线交换功率、调整可调负荷来保证系统功率平衡和各项约束满足。这种先决策、后调整的结构非常符合微电网的实际运行逻辑。日前调度确定机组启停和基础出力日内根据实际风光出力进行调整。第一阶段决策主要是一些慢变量比如机组启停状态这些变量一旦定了很难快速改变第二阶段决策是快变量功率可以在较短时间内调整。两阶段鲁棒优化的目标是在最恶劣的不确定性实现下总成本包括第一阶段发电成本和第二阶段调整成本最小。这里最恶劣三个字是核心它不再追求概率意义上的最优期望而是追求最小化最坏情况下的成本。这个思路很符合工程保守性原则——我先保证任何情况下都不翻车再考虑省油。1.3 为什么选择CCG算法而不是Benders分解有了模型之后求解算法是关键。两阶段鲁棒优化是典型的min-max-min三层结构没法直接扔给求解器。主流求解方法有两种Benders分解和列约束生成算法CCG。Benders分解的思路是把第二阶段返回的信息以割平面的形式加到主问题中但每次迭代只加入一组约束来逼近可行域。CCG算法的核心区别在于每次迭代不仅加入割平面还引入一组新的第二阶段决策变量形成新的列。这带来的直接好处是CCG算法的上下界收敛速度远快于Benders分解。实践中的数据是对于微电网规模的算例CCG算法通常十几到二十几次迭代就能收敛而Benders分解可能需要上百次迭代。另外CCG算法生成的割平面更紧能有效减少主问题的求解规模膨胀速度。这也是为什么近年来电力系统领域的核心期刊论文基本上都采用CCG算法。如果你在论文里用Benders分解求解两阶段鲁棒优化审稿人大概率会问一句为什么不用CCG所以我这里也只讲CCG。2. 模型建立不确定集合、目标函数与约束条件的完整推导2.1 不确定集合的选取盒式集合与不确定预算两阶段鲁棒优化的第一步是定义不确定参数的取值范围。最常用的是盒式不确定集合即每个不确定参数都有明确的上下界。这个上下界可以直接从历史预测误差的统计结果中获取比如取预测误差95%置信区间。但是盒式集合有个问题如果所有不确定参数同时取到最恶劣值结果会非常保守。比如风电、光伏、负荷同时达到极端值这样的场景在现实中几乎不可能出现但模型会围绕这个极端场景做调度导致成本飙升。为了解决这个问题引入不确定预算Γ来控制同时偏离预测值的参数个数。Γ0时就是确定性优化Γ越大模型越保守Γ取不确定参数总数时就是最保守的情况。确定Γ取值的方法是遍历Γ从0到最大值计算对应的鲁棒调度成本画出成本-保守度曲线根据工程需求选择可接受的保守度水平。这个曲线后面我会专门讲怎么画。2.2 主问题与子问题的数学形式两阶段鲁棒优化的一般形式可以写成min_x c^T x max_u min_y d^T y其中x代表第一阶段决策变量机组启停状态、储能充放电状态等0-1变量u代表不确定参数风电出力、光伏出力、负荷y代表第二阶段决策变量机组出力调整量、储能功率等。约束条件分为三组第一阶段约束Ax ≤ b包括机组启停约束、设备容量约束等耦合约束Bx Cy ≤ f把第一阶段决策和第二阶段决策耦合在一起如功率平衡约束第二阶段约束Dy ≤ g(u)包含与不确定参数直接相关的约束这个min-max-min结构直接求解是不可能的。CCG算法的思路是将其分解为主问题和子问题。主问题在给定一组不确定参数场景下求解第一阶段决策和第二阶段成本之和的最小值子问题在给定第一阶段决策下寻找最恶劣的不确定参数实现和对应的第二阶段最小成本。2.3 KKT条件与大M法处理双层优化子问题是一个max-min形式。内层min关于y是线性规划可以通过强对偶定理将其转化为max形式这样内外两层都是max合并成一个max问题。对偶变换的关键步骤是引入对偶变量π将内层min问题转换成对偶max问题。转化后的子问题目标函数包含不确定变量u和对偶变量π的乘积项。由于u和π都是变量这个乘积是非线性的需要进一步处理。处理方法有两种一种是基于KKT条件的方法对第二阶段问题写出KKT条件其中互补松弛条件通过引入二进制变量和大M常数进行线性化一种是直接枚举不确定集合的顶点因为最优解一定在盒式集合的顶点处取得枚举所有顶点分别求解取最大值即可。顶点枚举在不确定参数较少时可行参数多了以后组合爆炸所以更通用的是KKT线性化方法。这里有个关键的经验大M的取值非常敏感。取值过大会导致数值稳定性问题CPLEX和Gurobi都会报出数值告警取值过小则可能把可行解错误剪掉。我的做法是先求解一个松弛模型获取变量的大致量级再在此基础上放大10到100倍作为大M值。后面第6章我会详细说这个问题。3. 环境准备MatlabYalmipCPLEX/Gurobi的安装与配置3.1 求解器选型CPLEX还是Gurobi做两阶段鲁棒优化求解器的选择直接影响求解效率和复现难度。目前主流选择是CPLEX和Gurobi。CPLEX在电力系统领域有很长的使用历史很多经典论文的算例都是用CPLEX跑的因此复现文献结果时优先用CPLEX能更好地对齐论文中的参数和结果。但CPLEX的许可证获取流程相对繁琐对个人用户不太友好社区版只支持小规模算例。Gurobi近年来的性能表现非常强劲在某些MIP问题上比CPLEX快30%以上而且学术许可证申请非常方便用学校邮箱申请后谷歌浏览器直接下载就能用。在分布式能源、微电网方向的新论文中Gurobi的使用比例逐年上升。我的建议是两套求解器都配置好代码中通过一个参数切换求解器这样既能复现老论文也能充分发挥Gurobi的性能。下面这张表是我在同一个微电网算例下的实测数据后面第5章会详细展开。3.2 安装配置中容易忽略的版本匹配问题Yalmip是个建模工具箱需要下载后添加到Matlab路径。下载渠道有两个GitHub仓库或Yalmip官网我建议直接从GitHub下载最新Release版本保证和最新版Matlab的兼容性。CPLEX和Gurobi安装时要注意安装路径不要包含中文和空格Windows下需要把求解器的bin目录添加到系统环境变量PATH中Matlab中要正确设置求解器路径Gurobi安装目录下有一个setup.m脚本运行一次即可最容易翻车的是版本匹配问题。Yalmip只是一个建模层它通过底层接口调用求解器所以Yalmip版本、求解器版本、Matlab版本三者之间必须兼容。我的经验是Matlab R2020b以上版本搭配Yalmip 2021年后发布的版本搭配CPLEX 12.10或Gurobi 9.5以上版本这个组合最稳。这里给出一个快速验证是否配置成功的代码% 验证yalmip和求解器是否就绪 yalmip(clear) x sdpvar(1,1); optimize([x 1, x 2], -x, sdpsettings(solver, cplex)); % 如果求解成功x应该等于2 value(x)如果这段代码能正常输出2说明环境没问题。如果报错优先检查求解器路径是否添加成功以及Matlab能否正确调用求解器的动态链接库。4. 代码骨架主问题-子问题迭代框架的逐步实现4.1 数据初始化与参数设置在写迭代框架前先把所有基础数据准备好。以一个包含风电、光伏、储能、微型燃气轮机和负荷的典型微电网为例需要准备以下数据各分布式电源的容量和成本参数储能的容量、充放电效率、SOC上下限不确定参数风电、光伏、负荷的预测值和波动范围微电网与大电网的联络线功率限制和购售电价格分时电价时段划分推荐把所有参数写在一个结构体中方便后续修改和扩展。比如% 基础参数定义 mpc struct(); mpc.num_gen 2; % 微型燃气轮机数量 mpc.gen_capacity [1.0, 1.5]; % 单位MW mpc.gen_ramp 0.3; % 爬坡约束 MW/h mpc.battery_capacity 2.0; % 储能容量 MWh mpc.battery_power 0.5; % 储能最大功率 MW mpc.eta_ch 0.95; % 充电效率 mpc.eta_dis 0.92; % 放电效率 mpc.T 24; % 调度时段数参数定义清楚后后续改算例只需要改这个结构体代码其他部分不用动方便做敏感性分析。4.2 主问题建模Yalmip的关键语法主问题的核心是在已知的一组不确定性实现最初给一个初始场景后续是子问题返回的最恶劣场景下求解第一阶段决策和第二阶段调整成本之和的最小值。Yalmip建模的核心是先用sdpvar定义变量然后写约束和目标函数。以下是主问题建模的核心代码框架% 主问题变量 x_start binvar(mpc.num_gen, mpc.T, full); % 机组启停状态 p_gen sdpvar(mpc.num_gen, mpc.T, full); % 机组出力 p_bat sdpvar(1, mpc.T, full); % 储能充放电功率 soc sdpvar(1, mpc.T, full); % 储能荷电状态 theta sdpvar(1, 1); % 第二阶段成本近似变量 % 主问题约束集合 Constraints []; % 功率平衡约束、机组出力上下限、储能SOC递推等... % 不确定参数在这里被固定为当前已知的最恶劣场景 u_k % 目标函数 Objective sum(sum(c_gen .* p_gen)) theta; % 求解主问题 optimize(Constraints, Objective, sdpsettings(solver, cplex, verbose, 0));这里有个重要逻辑主问题中第二阶段成本并不是直接计算的而是通过变量theta来近似。每次迭代后CCG算法会加入新的约束theta ≥ 新场景下的第二阶段成本逐步逼近真实的第二阶段成本。这个逻辑和Benders分解的割平面思想很相似但CCG更进一步会同时加入新的场景变量。4.3 子问题求解与割平面生成子问题在给定第一阶段解x*后求解最恶劣不确定性场景下的最小调整成本。经对偶变换后的子问题是一个max问题直接用Yalmip建模即可% 子问题给定x_start后求解最恶劣场景 u_wind sdpvar(1, mpc.T, full); % 风电不确定变量 u_pv sdpvar(1, mpc.T, full); % 光伏不确定变量 % 加入不确定集合约束盒式集合 不确定预算 Constraints_sub [...]; % 加入对偶可行域约束 % 目标函数取max Objective_sub [...]; optimize(Constraints_sub, -Objective_sub, sdpsettings(solver, gurobi, verbose, 0)); % 提取最恶劣场景 u_wind_star value(u_wind); u_pv_star value(u_pv); sub_cost value(Objective_sub);子问题求出的最优目标值就是当前第一阶段决策下的最恶劣成本。如果这个成本大于主问题目标中theta的当前值说明主问题的theta低估了第二阶段成本需要将当前场景加入主问题并更新上界。4.4 收敛判据与迭代逻辑CCG算法的主循环逻辑如下LB -inf; % 下界 UB inf; % 上界 Gap 1; iter 1; max_iter 30; while Gap 0.01 iter max_iter % 第1步求解主问题得到下界 solve_master_problem(); LB max(LB, objective_master); % 第2步将主问题得到的x*带入子问题求解最恶劣场景和成本 [sub_cost, worst_u] solve_subproblem(x_star); UB min(UB, sub_cost first_stage_cost); % 第3步计算相对gap Gap abs(UB - LB) / abs(UB) * 100; % 第4步将最恶劣场景加入主问题场景集合添加割平面 add_cut_to_master(worst_u, theta_val); iter iter 1; end收敛判据我用的是相对gap小于1%实际运行中通常10到15次迭代就能收敛。如果超过30次还不收敛基本可以判断是子问题或主问题建模有逻辑错误需要回头检查。要注意一个细节每次迭代主问题的规模会增加因为新增了一组场景变量和一组约束。如果迭代次数很多主问题求解时间会快速上升。所以一定要确保子问题求解正确减少无效迭代。5. CPLEX与Gurobi求解对比实测数据与调优建议5.1 相同模型下两种求解器的性能对比我在一个包含2台微型燃气轮机、1个储能系统、风电光伏负荷共计9个不确定参数的微电网算例下分别用CPLEX 12.10和Gurobi 9.5.2求解相同的两阶段鲁棒优化模型统计结果如下指标CPLEX 12.10Gurobi 9.5.2总迭代次数1818主问题单次平均求解时间4.8s3.1s子问题单次平均求解时间0.9s0.6s总求解时间102s66s最终目标函数值9846.329846.32可以看到两种求解器得到的目标函数值完全一致因为求解的都是同一个模型最优解应该是一样的。差异主要体现在求解时间上Gurobi在做MIP问题时整体比CPLEX快约35%。如果你的算例规模更大、不确定参数更多这个差距会更明显。有一点需要注意CPLEX在某些特定类型的问题上可能表现更优比如强数值病态的模型CPLEX的数值鲁棒性历史口碑更好。所以在论文复现时我会先用CPLEX跑一遍确认结果和文献一致再用Gurobi做后续的敏感性分析和大规模扩展。5.2 求解器参数调优的实战经验求解器默认参数通常不是最优的。以下是我调参后效果比较明显的几个配置对于CPLEXMIP gap容差设置为1e-4默认是1e-4但有时为了加速可以放宽到1e-3Threads设置为实际物理核心数默认是0自动在虚拟机里可能被限制关闭可选的回退计算crossover等能减少一部分求解时间对于GurobiMIPGap设置1e-4Threads同样设置为物理核心数NumericFocus设为2或3当出现数值警告时能显著改善稳定性但会增加少量求解时间TimeLimit设置一个合理的时间上限防止极端情况下的死循环在Yalmip中这些参数通过sdpsettings传入options sdpsettings(solver, gurobi, ... gurobi.MIPGap, 1e-3, ... gurobi.Threads, 8, ... gurobi.NumericFocus, 2, ... verbose, 2); optimize(Constraints, Objective, options);另一个实战经验是热启动。CCG算法每次迭代的主问题变化不大可以把上一次的解作为初始解传入能显著减少求解时间。Yalmip中实现方式是使用assign函数给变量赋初值然后求解时设置savesolveroutput和solver相关参数。6. 复现中踩过的六个坑及解决方案6.1 大M值的选择这是最让我头疼的一个问题。大M法在处理互补松弛条件时是必需的但M的取值过大会导致数值稳定性问题CPLEX和Gurobi都会报告数值警告求解结果中出现inf或NaN。我的解决方法是先算出对偶变量和第二阶段变量的量级比如通过求解松弛后的子问题获取近似范围然后取这个范围的10倍作为M的初始值如果出现数值问题再逐步减小。实际测试中M取10000时Gurobi会报数值警告改为1000后问题消失结果完全一致。所以M不是越大越好够用就行。一般来说M取到正常变量量级的10到100倍就能保证不cut掉最优解同时不会引发数值灾难。6.2 非线性项的处理子问题经对偶变换后目标函数中不确定变量u和对偶变量π的乘积是非线性的。很多人直接把这个乘积用u .* pi写进Yalmip结果求解器报错或者求解极慢。正确的做法是将其拆解为单独的变量和约束利用对偶变量的取值结构。具体来说根据不确定变量u的上下界最恶劣场景一定是在盒式集合的顶点处取得。因此可以将u的所有顶点组合枚举出来在每个顶点下求解子问题取最大值作为最恶劣场景成本。但这种方法只在不确定参数数量少时可行。我推荐的做法是引入辅助变量将乘积项线性化后再交给求解器。这在Yalmip中可以通过binvar将变量离散化或者用implies实现分段线性化。但最干净的做法还是将对偶变换推到极致直接推导出不含乘积项的等价形式。我在代码注释里详细写了推导过程这里不展开但提醒一句不要在Yalmip里直接写非线性乘积项除非你用的是非线性求解器否则必踩坑。6.3 对偶变量的维度问题第二次踩坑是对偶变量的维度对应错误。第二阶段约束中的每个约束都需要一个对偶变量而且对偶变量的数量必须和约束数量严格对应。多写一个、少写一个、顺序写反都不会直接报错但结果会完全错误。我的排查方法先把确定性版本即Γ0跑出来和直接求解一个等价的单层优化结果对比。如果对偶变量维度有问题单层优化和两阶段的结果必然对不上很快就能发现问题。对偶变量维度的检查方法是用size函数查看变量维度和约束条件数量做比对。注意Yalmip中约束条件写成矩阵形式时一个矩阵约束算一个约束块但对偶变量仍然要逐行对应这块是最容易出错的。6.4 收敛判据过严导致死循环刚开始我把gap阈值设为1e-6结果跑了40多次迭代还没收敛主问题规模越来越大单次求解时间从几秒涨到几十秒整体陷入死循环。原因是数值精度有限上下界在一个很小的区间内来回摆动无法达到1e-6的水平。后来我把阈值改为0.5%结果14次迭代就收敛了总成本只相差0.2%左右完全满足工程精度需求。建议gap阈值设置在0.5%到2%之间不需要过分追求微小gap。对审稿人来说你的gap收敛曲线和最终结果合理比gap绝对值多小更重要。6.5 Yalmip与求解器版本不兼容有段时间我把Yalmip更新到最新版结果原有代码大面积报错。原因是新版Yalmip在某些语法上做了调整比如sdpvar的默认属性、optimize的参数处理方式。解决方案是不要盲目追求最新版。找到一套稳定的组合后就固定下来。我目前的组合是Matlab R2022b Yalmip R20230623 Gurobi 10.0.1 CPLEX 12.10运行稳定了大半年期间没再出现版本类报错。如果你们实验室用的是旧版Matlab建议先把Yalmip版本锁定在对应的兼容版本不要轻易升级。6.6 数据归一化问题最后一个坑是量纲不统一。风电和光伏的预测值单位是MW成本单位是元/MWh储能SOC是百分比三者数量级差别很大。如果直接建模求解器内部的数值处理精度会受到影响特别是大M法场景下更容易出问题。我的做法是将所有变量统一到标幺值系统基准功率取微电网最大负荷成本统一到元/MWh。这样所有变量的量级都在0到几之间求解器的数值稳定性大幅提升。这个操作对结果没有影响只是换了个单位体系但对求解器来说却是天壤之别。7. 结果解读如何从鲁棒调度方案中提取有效信息7.1 最恶劣场景的识别与分析CCG算法迭代收敛后子问题返回的最后一个场景就是当前调度方案下系统面临的最恶劣不确定性实现。这个场景的工程意义非常大它告诉你系统的瓶颈在哪里。比如如果最恶劣场景中光伏出力取下限且风电出力也取下限说明系统对新能源出力的下行风险最敏感。我在实际项目中会专门把最恶劣场景提取出来做可视化和预测场景画在同一张图上直观展示鲁棒方案的设计依据。这对外审和组会汇报都非常有帮助。你还能从最恶劣场景中看出哪些时段最容易出现供需失衡为后续的储能扩容或联络线升级提供数据支持。7.2 鲁棒调节参数与系统成本的权衡曲线这是论文中最常用的一个结果图。做法很简单将不确定预算Γ从0遍历到最大值每个Γ值下运行一遍完整的CCG算法记录总成本。然后把Γ-总成本的关系画成曲线。这条曲线揭示了模型保守度和经济性之间的权衡关系。我跑过的算例中Γ从0增加到9时总成本上升了约18%。但在Γ增加到6以后成本上升的速度明显放缓曲线趋于平坦。这说明系统对不确定性的敏感度在Γ较小时最大实际工程中取Γ6左右就能获得较好的鲁棒性而不过分牺牲经济性。这张曲线的同样可以用于论文中的敏感性分析比如换不同的波动范围、不同的储能容量设置观察权衡曲线的形状变化导出一些有工程意义的结论。7.3 后续研究的扩展方向复现并跑通两阶段鲁棒优化后可以在多个方向迭代创新。比如将条件风险价值引入第二阶段目标函数在鲁棒优化框架下考虑风险偏好或者将不确定集合从盒式扩展到椭球式或多面体式更精细地刻画不确定性还可以把模型从单微网扩展到多微网互联两阶段鲁棒优化同样适用。我自己目前在做的方向是把碳交易机制引入两阶段鲁棒优化框架在第二阶段加入碳配额和碳价不确定性。碳价的波动本质上也是一种不确定性完全可以用不确定集合描述这样就能在同一个框架下同时处理新能源出力和碳价的双重不确定性。实际跑通第一版两阶段鲁棒优化代码后我的一个体会是这类算法的复现难点不在算法本身而在于建模细节和数值处理的完整性。建议先把确定性模型跑通并验证正确再改造成鲁棒优化模型这样每一步的结果都有对比排查bug会容易得多。你的第一版代码不追求效率先把结果跑对再逐步优化。最后再分享一个小技巧迭代过程中记录每一轮的上下界、gap、最恶劣场景对算法的行为做完整日志。CCG算法出了问题时这份日志可以帮你快速定位是主问题建模错了还是子问题对偶写错了还是收敛条件设置不合理。我自己是把每次迭代的日志打印成表格保存在文本文件里回头排查问题直接翻日志比对着代码猜效率高太多了。这个习惯帮我节省了大量调试时间也推荐你从第一次跑通代码时就养成。本文还有配套的精品资源点击获取