简介本资源是面向电力系统优化方向研究生与科研人员的微电网两阶段鲁棒经济调度完整实现方案聚焦解决含不确定性风电出力下的调度保守性与经济性平衡问题。压缩包共13个文件1.57MB含4个核心MATLAB脚本MP.m、SP.m、main_1.m、MP2.m、3份PDF文献含经典论文《微电网两阶段鲁棒优化经济调度方法》及拓展研究、6张关键结果图如目标函数收敛曲线、场景调度对比图等覆盖建模、求解、可视化全流程。已有3182人学习下载印证其作为入门两阶段鲁棒优化的高口碑参考价值。代码完全原创且完美复现注释详尽目标函数与约束均以紧凑矩阵形式表达结构工整特别嵌入鲁棒调节系数接口支持灵活调整保守程度便于适配不同不确定性集或迁移至综合能源系统等拓展场景。 “微电网两阶段鲁棒优化经济调度”这个方向这几年几乎是电力方向硕博论文里的标配。我最近把整套东西从数学推导到MATLABYALMIPCPLEX/Gurobi代码完整复现了一遍跑了标准微网算例也踩了不少坑。这篇文章就把整个复现过程捋一遍两阶段鲁棒模型怎么建、CCG算法怎么推、代码怎么一步步写出来、哪些地方最容易出错。如果你正在写论文、正在复现类似代码或者只是想把“两阶段鲁棒优化”这几个字真正落到能跑的代码上这篇应该能省你不少时间。1. 微电网经济调度为什么要用“两阶段鲁棒”1.1 确定性调度的局限与不确定性来源传统微电网经济调度本质上是在满足功率平衡、设备出力上下限、储能SOC约束等条件下最小化总运行成本。目标函数一般包括燃气轮机燃料成本、与大电网交互的电费、储能老化成本等。只要光伏、风机、负荷都按预测值给定这类问题就是一个标准的混合整数线性规划用YALMIP调CPLEX或者Gurobi几分钟就能解出来。但问题恰恰出在“预测值”这三个字上。光伏出力和风速具有很强的随机性夏天的云层移动可能让光伏在十分钟内掉一半出力负荷曲线也总被各种偶发因素扰动。传统确定性调度只对单一预测场景优化实际运行中一旦风光偏差较大就只能靠备用容量硬扛极端情况下可能违反功率平衡约束甚至导致切负荷。随机优化把不确定性建模为概率分布理论上更精细但实际中很难拿到准确的分布信息概率场景再多也怕分布假设本身不靠谱。鲁棒优化走的是另一条路不去猜概率分布只给定一个不确定集合认为实际值一定落在这个集合内然后求“最坏情况下的最优解”。这种思路在工程上更稳也更好解释所以近年在微电网调度里越来越流行。1.2 两阶段结构的实际工程含义两阶段鲁棒优化里的“两阶段”对应的其实是电力调度的两层时间结构。第一阶段是“日前决策”发生在实际风光出力已知之前需要提前决定燃气轮机的启停状态、与大电网的购售电计划等慢动作变量第二阶段是“实时调整”发生在前一日调度执行过程中看到实际风光出力之后再根据一阶段确定的开机组合调整各机组的实际出力、储能的充放电功率以最小的调整成本保证功率平衡。打个比方一家餐厅前一天晚上就把厨师排班和采购菜单定好了这是第一阶段第二天客人实际来了多少、哪些菜卖完了后厨再临时调整备菜顺序和炒菜节奏这是第二阶段。第一阶段决策要留有足够弹性才能让第二阶段的实时调整不失控。用数学语言说两阶段鲁棒优化的决策顺序就是“在这里优化然后在那里优化最后看最坏情况”。1.3 模型假设、不确定集合与完整数学表达我把复现用的模型定义为一个含光伏、风电、蓄电池、燃气轮机和上级电网互联的交流微电网调度周期为24小时单位调度间隔1小时。不确定性只考虑光伏出力、风电出力和负荷预测偏差。不确定集合采用最常见的盒式预算约束形式实际值等于预测值加偏差偏差绝对值有上限同时整个调度周期内所有不确定量的偏差绝对值之和不超过预算值。预算值的作用很关键它表达了“同一时刻最多只有几个变量同时达到最坏偏差”的保守程度取值越大结果越保守。两阶段模型可以写成这样一个min-max-min结构外层第一阶段求解日前决策和总成本中间层max负责寻找最坏情况下的不确定场景内层min对应第二阶段实时调整成本。目标函数由开机成本、燃料成本、购售电成本、储能成本组成。约束包括功率平衡、光伏/风电出力界限、燃气轮机爬坡和出力上下限、储能SOC递推与容量约束、与大电网交互功率限制等。这个min-max-min结构就是整套代码的核心骨架。很多初学者一看到三层嵌套就懵了实际上后面的CCG算法就是为了把它拆成一个可求解的主问题加一个可求解的子问题。2. 核心求解算法CCG列与约束生成推导2.1 主问题MP与子问题SP的分工两阶段鲁棒优化没办法直接交给求解器因为max-min这个内层结构破坏了MILP的标准形式。工程上最常用的方法是CCG算法也叫列与约束生成。它的思路很朴素把问题拆成主问题MP和子问题SP主问题负责决策第一阶段变量并给出下界子问题负责在给定第一阶段变量下寻找最坏场景并给出上界然后不断把子问题产生的割加到主问题里循环迭代直到上下界收敛。主问题MP的形式是把第二阶段的目标和约束用一个辅助变量替换同时把已经发现的“最坏场景”对应的第二阶段约束显式加入。每轮迭代都会增加一组变量和一组约束这就是“列与约束生成”名字的由来。2.2 子问题的对偶推导与大M处理子问题SP是一个双层结构外层max找最坏的不确定场景内层min求该场景下的实时调整成本。这个结构不能直接求解但只要把内层min用强对偶定理转成max就把“max-min”变成了“max-max”合并成一个max问题也就是一个普通的线性规划。这里有一个常见的坑对偶推导必须在连续变量层面做第二阶段如果有二进制变量比如储能充放电状态就不能直接对偶。我复现时对储能模型做了简化处理用连续变量表示充放电功率再用充放电效率模型规避二进制变量或者把充放电状态放到第一阶段决定。如果你的模型里第二阶段确实有0-1变量那就没法用强对偶得引入KKT条件或线性化技巧复杂度会明显增加。对偶之后原目标里会出现“对偶变量乘以不确定变量”的双线性项。解决方法是用大M法把这个乘积项线性化。大M的取值在数值上很敏感取得太大容易导致求解器数值崩溃取得太小又会错误地限制可行域。实操中我一般先算一下相关约束的量级再放大10到100倍并且不同约束可以设不同的M值。2.3 CCG迭代流程整个迭代过程可以归纳为以下步骤初始化下界LB为负无穷上界UB为正无穷迭代次数k1。求解主问题MP得到第一阶段最优解和最优目标值将目标值更新为当前下界LB。把第一阶段解固定代入子问题SP求解得到最坏场景和对应的第二阶段最优成本。计算上界UB 第一阶段成本 第二阶段最坏情况成本。检查UB - LB是否小于设定的收敛阈值满足则停止并输出当前解。不满足则把当前识别出的最坏场景作为新一列生成一组新约束加入主问题kk1回到第2步。这个框架本身很简单但每一步落实到代码里都有不少细节。接下来我重点讲yalmip代码怎么写。3. 基于MATLABYALMIP的代码实现与关键细节3.1 环境配置YALMIP、CPLEX与Gurobi我用的是MATLAB R2022aYALMIP最新版求解器同时装了CPLEX和Gurobi。先说配置YALMIP本身是一个建模层不负责求解它把模型翻译成求解器能识别的格式。安装YALMIP就是把文件夹放到本地路径然后在MATLAB里执行addpath(genpath(yalmip路径))再用savepath保存路径。CPLEX和Gurobi都是商业求解器但都提供学术免费许可。安装完之后最关键的一步是让MATLAB能识别它们。对CPLEX在MATLAB里执行addpath(cplex安装目录/cplex/matlab)对Gurobi需要先运行gurobi_setup脚本。完成后可以用yalmiptest命令测试哪些求解器可用。我自己的经验是同一个模型在Gurobi上求解速度通常比CPLEX快一些但CPLEX在某些整数问题上也很稳具体选哪个看个人习惯。建议把上述配置写在代码开头% 求解器路径配置按实际安装位置修改 addpath(genpath(D:/Program Files/YALMIP-master)); addpath(C:/Program Files/IBM/ILOG/CPLEX_Studio221/cplex/matlab); % 或 Gurobi % run(C:/gurobi1100/win64/matlab/gurobi_setup.m); yalmip(clear);3.2 不确定集合与主问题的YALMIP建模主问题建模前先把不确定集合相关的预测值和偏差范围定义好。这里我用一个长度为24的向量表示各时刻光伏预测出力偏差最大值为预测值的20%。决策变量方面第一阶段变量包括燃气轮机启停状态、购售电状态等二进制变量第二阶段变量包括各时刻机组出力、储能功率、SOC等连续变量。在YALMIP里二进制变量用binvar声明连续变量用sdpvar声明。主问题的辅助变量eta代表第二阶段最坏情况成本的下界。第一轮迭代不加任何场景约束只加基础约束和第一阶段约束。每轮迭代后把子问题识别出的最坏场景u_k作为新参数生成一组第二阶段约束并加入主问题同时用新的第二阶段变量y_k实例化这些约束。需要注意的是YALMIP的变量可以在循环里动态拼接但每次迭代实例化新的第二阶段变量时要小心命名空间冲突。我习惯用eval或cell数组来管理不同迭代下的变量最省心的方式是直接建立元胞数组比如y{k}表示第k轮新增的第二阶段变量。主问题核心代码片段如下% 主问题变量 z binvar(1, 24); % 燃气轮机启停状态示意 x sdpvar(1, 24); % 一阶段其他决策示意 eta sdpvar(1, 1); % 二阶段成本下界 % 约束容器 Constraints []; % 第一阶段基础约束 % ... Constraints [Constraints, sum(z) 4]; % 示例约束 % CCG迭代中加入的场景约束示意 % for k 1:K % y_k sdpvar(1, 24); % Constraints [Constraints, eta cost_k * 0 f(y_k, u_k)]; % Constraints [Constraints, g(x, y_k, u_k) 0]; % end Objective a*z b*x eta; ops sdpsettings(solver, gurobi, verbose, 0); optimize(Constraints, Objective, ops);3.3 子问题的对偶实现与双线性项线性化子问题在YALMIP里有两种写法。一种是用dual函数直接提取对偶变量手动构造对偶问题另一种是直接调用YALMIP的dualize命令自动生成对偶。自动对偶省事但可读性差出了问题不好排查。我推荐手动对偶虽然推导过程麻烦但你能完全掌控每一条约束的来龙去脉。手动对偶的步骤是先把内层min问题写清楚包括它的约束和变量然后根据标准线性规划对偶规则构造对偶目标和约束最后把外层max与内层max合并得到一个同时包含对偶变量和不确定变量的max问题。求解这个合并问题的最优解就是当前最坏场景下的最坏成本以及对应的最坏场景值。合并后的目标函数中会出现对偶变量lambda与不确定变量u的乘积项。线性化方法是用大M引入辅助变量和辅助约束。这部分是代码里最容易出错的地方我最初就是因为漏了一个M的约束导致迭代不收敛。子问题求解核心代码示意% 子问题变量给定第一阶段的解 x_fixed lambda sdpvar(size_dual, 1); % 对偶变量 u sdpvar(1, 24); % 不确定变量 % 对偶目标包含线性化后的项 Objective_dual d*lambda lambda*A_u*u; % 线性化处理乘积项 % 引入辅助变量theta约束 theta M*u, theta M*lambda 等 for t 1:24 theta{t} sdpvar(1, size_dual); Constraints_sub [Constraints_sub, theta{t} M * u(t)]; Constraints_sub [Constraints_sub, theta{t} M * lambda]; % 注意实际线性化需要根据乘积项符号保留正/负部分 end optimize(Constraints_sub, -Objective_dual, ops); % 因为YALMIP默认求min取负号 UB -value(Objective_dual) first_stage_cost;这里有一个容易被忽略的细节子问题求解完成后要把当前最坏场景u的值保存下来带到主问题里生成新的割。我见过不少实现只存了成本没有存场景导致主问题加的割全是同一个场景迭代永远不收敛。3.4 CCG主循环的收敛判断与性能优化CCG主循环的结束条件一般用相对间隙或绝对间隙。我用的收敛判据是(subopt - obj)/abs(obj) 1e-4同时设置最大迭代次数作为保险。但从实际执行结果看大多数算例在3到5次迭代内就能收敛如果超过15次还没收敛大概率是代码bug或者M取值有问题。性能优化方面有几个小技巧很管用。第一子问题每次求解完可以先把当前最坏场景跑一遍完整确定性优化验证这个场景是否真的会让系统运行成本偏高做一层逻辑自洽性检查。第二在循环开头加入一个“热启动”机制把上一轮的求解结果作为当前轮的初始值可以显著减少重优化时间。第三如果算例规模很大可以考虑把第二阶段连续模型改成线性化后的凸问题避免求解器反复处理非线性。主循环核心代码LB -1e6; UB 1e6; k 1; tol 1e-4; while (UB - LB) / abs(UB) tol k 15 % 求解主问题 optimize(Constraints, Objective, ops); LB value(Objective); x_fixed value(x); % 求解子问题 optimize(Constraints_sub, -Objective_dual, ops); UB first_stage_cost value(Objective_dual); % 保存最坏场景 u_star u_star{k} value(u); % 生成新割并加入主问题 y_k sdpvar(1, 24); Constraints [Constraints, eta ...]; % 新增成本割 Constraints [Constraints, ...]; % 新增可行性约束 k k 1; end4. 典型算例与仿真结果分析4.1 算例参数设置复现用的微网结构是一台200kW燃气轮机、100kW光伏、80kW风机、200kWh储能电池最大充放电功率50kW与大电网联络线功率上限100kW。负荷曲线采用典型工业园区日负荷光伏和风力采用标准日曲线并叠加预测误差。不确定预算值取总时段的一半即12表示“最恶劣情况下一天内最多有12个时段出现极端偏差”。成本参数方面燃气轮机单位发电成本0.6元/kWh向上级电网购电价格分时计价峰时1.2元/kWh谷时0.4元/kWh售电价格固定0.35元/kWh。储能充放电效率均取0.95。4.2 迭代收敛过程以Gurobi作为求解器模型规模约包含1500个变量和2800条约束。运行环境是i7-12700H处理器、16GB内存。实测收敛过程如下表所示迭代次数下界LB元上界UB元间隙Gap15246.35872.810.67%25542.15811.44.63%35659.75772.51.95%45720.35758.20.66%55739.65748.90.16%第5轮迭代之后Gap降到0.16%满足1e-3的收敛阈值程序在第5轮停止。整体运行时间约40秒。我在多台不同配置的电脑上跑过同样代码结果基本一致这也是“完美复现”的一个体现。4.3 调度结果解读与鲁棒性对比从最优调度结果看系统在白天的光伏大发时段倾向于减少燃气轮机出力同时给储能充电在夜间电价低谷时段从电网购电补充负荷并在傍晚高峰时段放电。储能SOC曲线呈现清晰的“谷充峰放”特征。为了说明鲁棒优化的价值我对比了确定性调度不考虑不确定性和两阶段鲁棒调度在最坏场景下的表现。确定性调度在最坏场景下总成本比预测场景高出约11%而鲁棒调度只高出约5%。差距的来源在于鲁棒调度会主动保留更多向上调节空间比如燃气轮机不会在预测场景下满发储能也会预留一部分容量应对突发偏差。此外我测试了不同预算值的效果。预算值从0增加到24总成本单调上升从约5600元上升到约6000元。这说明预算值是一个“保守度旋钮”实际使用时应根据对风险的接受程度动态调整。论文里常用的做法是画成本-预算曲线来展示不确定性对运行经济性的影响这条曲线也是很多审稿人喜欢看的图。5. 常见问题与调试经验实录5.1 求解器报错与YALMIP状态检查YALMIP优化结束后一定要检查求解器返回状态我习惯用以下代码判断optimize(Constraints, Objective, ops); if ~strcmp(info.problem, Successfully solved) yalmiperror(info.problem); error(求解失败); end如果显示数值问题Numerical issues优先检查是不是大M取值过大或约束量级不一致。YALMIP本身不打印详细对偶信息时可以在sdpsettings里开debug参数sdpsettings(debug, 1)它会帮你定位是哪条约束或变量引起的问题。这个开关在模型不收敛时非常有用。5.2 常见问题排查速查表现象可能原因解决方案主问题求解报Infeasible第一阶段约束过紧或割约束写错先去掉所有割约束检查基本模型是否可解子问题最优值比主问题下界还低对偶推导遗漏约束或符号写反把对偶问题与原问题在固定场景下对比验证迭代不收敛且UB持续震荡割对应的最坏场景没正确传入主问题打印每轮u_star确认场景确实变化总不收敛但场景从未变化子问题求解失败一直返回同一个可行解检查子问题求解状态可能M取值太小Gurobi提示数值警告大M值过大或变量量级跨越太大缩小M值或对变量做归一化处理结果与论文对不上计算基准或成本参数不同逐条核对模型表达式和单位统一5.3 我踩过的几个坑和个人心得第一个坑是第二阶段互斥约束的处理。最初我在第二阶段用了储能的充放电0-1状态变量导致子问题无法直接对偶。后来我把储能模型改成“充放电效率连续变量”把充放电状态作为第一阶段决策彻底绕开了这个问题。如果你的模型必须保留第二阶段0-1变量那就得用嵌套KKT条件或者大M法把互补约束线性化代码复杂度会成倍上升我建议能简化就先简化。第二个坑是主问题里的辅助变量eta和最终成本的关系。很多初学者直接把eta当作最终成本参与计算但实际上下界LB应该是主问题的目标值而不是eta本身。这两个概念在迭代初期差别很大只有收敛后才会趋于一致。第三个坑是“可行性割”的处理。当子问题在某个第一阶段解下没有可行解时不能直接加最优割需要加可行性割。我复现时给子问题加入了松弛变量并令松弛变量在目标函数中有足够大的惩罚系数。这样即使某个一阶段解让子问题不可行求解器也能返回一个可行解并生成正确的割约束而不是直接报错终止。最后分享一个调试技巧在CCG主循环跑通之前先关闭不确定集合把预算值设为0让模型退化成普通确定性两阶段问题。这样能先验证模型基本框架的正确性再逐步引入不确定性。等确定性版本完全跑通、结果合理了再打开不确定集合和CCG迭代定位问题会容易得多。我每次新写CCG代码都按这个流程走基本能在半天内把所有低级错误清理干净。本文还有配套的精品资源点击获取