☰
微电网两阶段鲁棒优化实战:MATLAB+YALMIP+CPLEX建模与CCG求解全攻略
2026/10/27 9:38:08 网站建设 项目流程

简介:面向微电网经济调度与不确定性优化方向的研究生和工程师,这份代码基于MATLAB+YALMIP+CPLEX完整复现了《微电网两阶段鲁棒优化经济调度方法_刘一欣》一文。针对原文中导致模型难以求解的少量笔误,作者通过邮件与通讯作者确认并做了修正,使代码能够稳定求解任意随机生成的光伏与负荷场景。模型充分考虑分布式电源和负荷的不确定性,经系列推导将两阶段问题转化为主问题与子问题的混合整数线性形式,采用C&CG算法迭代求解,得到最恶劣场景下运行成本最小的调度方案。资源共8个文件,以6个m源码为主,涵盖参数设置、场景生成、主问题、子问题、绘图等模块,另有docx与txt格式的说明文件,压缩包仅272KB。代码注释详细、结构清晰,适合对照文献精读两阶段鲁棒优化的建模、转换与求解流程。目前已有7744人学习浏览,运行中遇到问题还可联系作者获得答疑支持。

1. 微电网两阶段鲁棒优化:为什么用 MATLAB+YALMIP+CPLEX 这一套最稳

微电网两阶段鲁棒优化,本质上是解决“预测不准”的调度问题:风光出力、负荷需求都有波动,你先得在“最坏场景”下做日前决策(第一阶段),等不确定量揭晓后再做实时调整(第二阶段)。用 MATLAB 建模、YALMIP 做抽象层、CPLEX 当底层求解器,是目前复现论文和做工程验证最顺手的组合。YALMIP 把复杂的鲁棒对等转化封装成了几行代码,CPLEX 负责啃下大规模混合整数线性规划这块硬骨头,MATLAB 则帮你把数据处理、结果可视化串起来。这套方案特别适合刚接触鲁棒优化的研究生和做微电网EMS预研的工程师,照着本文的步骤,你可以在自己的机器上从零跑通两阶段鲁棒优化,并清楚每个参数为什么这么设。

2. 两阶段鲁棒优化的建模基础:不确定性集合与CCG框架

2.1 先理清两阶段结构:第一阶段决策什么,第二阶段调整什么

微电网的调度通常分成两个阶段:日前计划(Day-ahead)和实时调整(Real-time)。第一阶段在未知风光出力、负荷确切值之前,决定机组启停、与大电网交互功率的上限、储能充放电的基准计划等“不可变”变量;第二阶段则在这些约束下,针对每个可能出现的不确定场景,给出最经济的调整量,比如机组出力微调、储能SOC修正、弃风弃光量、切负荷量。两阶段鲁棒优化的目标是最小化“第一阶段成本 + 最坏场景下的第二阶段成本”。

为什么采用两阶段而不是单阶段?单阶段随机优化需要假设不确定量的概率分布,计算量也随场景数爆炸;而两阶段鲁棒优化只关心不确定集合内的最坏情况,不需要精确概率,适合工程上保守性要求高的场景。常见的微电网两阶段问题包括:考虑光伏、风电、负荷不确定性的经济调度,考虑N-1安全约束的备用优化,以及储能容量规划与运行联合优化。你手上的标题大概率属于第一类,下面我就按这个方向展开。

2.2 不确定性集合的构造:盒式、椭球式与预算约束的取舍

两阶段鲁棒优化的核心之一是定义不确定量属于哪个集合。最常见的三种:

  • 盒式集合:每个不确定参数独立在上下界内波动,形式为 ( u \in [\hat{u}-\Delta, \hat{u}+\Delta] )。优点是简单,缺点是可能产生极端场景,过度保守。
  • 预算约束集合:在盒式基础上加入总偏差限制,形如 (\sum |u_i-\hat{u}_i| / \Delta_i \le \Gamma)。这个 (\Gamma) 称为不确定性预算,控制同时达到极端值的参数数量,是工程中最常用的折中方案。
  • 椭球式集合:用椭球描述相关性,但会引入二阶锥约束,求解速度明显下降。

我一般推荐从预算约束集合起步。它写起来简单,YALMIP 里甚至可以直接用uncertain和budget模型,或者手动构造。比如风光出力同时偏差,可以用一个总预算限制总的风光偏差绝对值之和。(\Gamma) 的经验取值在 0.5~1.5 之间,具体要靠灵敏度分析确定:把 (\Gamma) 从 0 逐步加到最大,看总成本上升曲线,选取拐点处的值。这一条也直接关系到最终方案的保守性和经济性。

2.3 为什么用 C&CG 而不是 Benders 分解

两阶段鲁棒优化常见的求解框架是 Benders 分解和列与约束生成(C&CG)。Benders 通过不断添加最优割来逼近最坏场景,收敛慢;C&CG 直接在主问题中引入第二阶段变量和约束,每次迭代加入一组新变量,收敛速度快,尤其适合包含二进制变量的问题(比如机组启停)。现在论文里普遍采用 C&CG,YALMIP 虽然不自带 C&CG 实现,但你可以自己写循环:主问题求解一次 MILP,子问题求解一次 max-min,然后提取最坏场景反馈给主问题。

有人问:能不能直接用optimize加'robust'标量直接求解?YALMIP 对带有不确定变量的模型可以自动做对等转化,但只适用于线性或二次约束的一阶段鲁棒问题,两阶段鲁棒需要显式使用C&CG。所以你别指望一行代码搞定,老老实实搭框架。用规范化解法:主问题的目标函数包含第一阶段成本和一个辅助变量 (\eta),子问题是第二阶段在最坏场景下的最小值,把子问题对偶后变成 max 问题,求得最坏场景后把该场景下的约束添加到主问题。

3. 用 YALMIP 声明不确定变量:从约束到可求解模型的三种写法

3.1 写法一:直接调用 robust 模块(适用简单单阶段)

YALMIP 自带鲁棒优化工具箱,在变量声明时用uncertain声明不确定变量,然后用optimize加-robust选项。比如:

sdpvar Pg Cg; sdpvar Pwind_hat; Pwind_hat_val = 50; delta = 5; sdpvar Pwind; uncertain(Pwind); Pwind >= Pwind_hat_val - delta; Pwind <= Pwind_hat_val + delta; % 简单成本约束示例 Constraints = [Pg + Pwind >= 100, Pg >= 0]; Objective = Cg * Pg; options = sdpsettings('solver','cplex','robust.dualize',1); optimize(Constraints, Objective, options);

这段代码里uncertain(Pwind)告诉YALMIP这是一个不确定变量,约束里对它给出上下界,YALMIP 会自动生成对等模型。robust.dualize控制是否对偶化,一般保持默认1即可。这个写法适合验证你的不确定性集合定义是否正确,但不适合两阶段问题。

3.2 写法二:手动 C&CG 主问题建模

两阶段鲁棒的主问题如下所示,第一阶段变量是常规变量,第二阶段变量对每个已找到的场景单独复制一份,约束中包含场景值:

% 第一阶段变量 sdpvar Pg_start Cstart Pgu Pgd Pgrid_buy Pgrid_sell SOC_start; % 第二阶段变量(每个场景 s 复制变量) Pgen_adjust = sdpvar(nscenario,1); SOC_adjust = sdpvar(nscenario,1); % 主问题目标 Objective_MP = Cstart*Pg_start + eta; % 主问题约束(基础约束) Constraints_MP = [...]; % 对每个场景 s 添加第二阶段约束 for s = 1:nscenario Constraints_MP = [Constraints_MP, ... Pgen_adjust(s) + Pwind_scenario(s) + Ppv_scenario(s) ... + Pgrid_buy - Pgrid_sell == load_scenario(s), ... SOC_adjust(s) == SOC_start + 0.5*Pgen_adjust(s)]; end

主问题是一个常规 MILP,直接用 CPLEX 求解。注意eta是辅助变量,它代表最坏场景下的第二阶段成本下界,子问题求解后会更新这个下界。

3.3 写法三:子问题 max-min 的对偶转化

子问题固定第一阶段变量后,变成:

[ Q(x) = \max_{u \in U} \min_{y} ; c^T y ]

这个 max-min 不能直接求解,常见做法是对内层 min 取对偶,变成单层 max 问题。然后在 YALMIP 里写出对偶模型。但手写对偶容易出错,我推荐用 YALMIP 的dualize功能:

% 固定第一阶段变量后的子问题 Constraints_SP = [...]; % 第二阶段约束,含不确定量 u Objective_SP = c'*y; [ Objective_SP_dual, Constraints_SP_dual ] = dualize(Objective_SP, Constraints_SP, y); % 然后对 u 最大化 Objective_SP_dual,且 u 在不确定性集合内 uncertain(u); Constraints_u = [u >= u_min, u <= u_max, ...]; optimize([Constraints_SP_dual, Constraints_u], -Objective_SP_dual, options);

注意dualize要求子问题是线性的,并且已经固定了第一阶段变量。取对偶后如果还有二进制变量,不能直接对偶,需要按场景枚举,或者用强对偶条件处理。多数微电网子问题只包含连续变量(如机组出力、储能、切荷量),所以这个路径可行。

3.4 参数说明:CPLEX 选项与容差

无论主问题还是子问题,CPLEX 的求解参数都很关键:

options = sdpsettings('solver','cplex', 'verbose', 1, ... 'cplex.mip.tolerances.mipgap', 1e-4, ... 'cplex.mip.tolerances.integrality', 1e-6, ... 'cplex.mip.strategy.names', 0);

mipgap是相对 MIP 间隙,工程上取 1e-3~1e-4 就够;取太严会拖慢速度。integrality是整型判定容差,默认即可。names置0可以省内存。另外,如果模型过大,可以试试'cplex.mip.strategy.file'将分支信息写盘,但会明显变慢,不建议一开始就开。

4. 完整求解流程:从数据到 CCG 循环收敛

4.1 算例构造:风电、光伏、负荷的不确定场景

我们用一个小型微电网算例:一个火电机组(成本 220 元/MWh)、一组储能(容量 2MWh,功率 0.5MW)、一个风电场和一个光伏电站,以及一个配电网连接点。单位时间设为 1 小时,调度周期 24 小时。风电出力预测值如下(单位 MW):

Pwind_hat = [15 14 13 12 11 10 12 15 18 20 22 23 22 20 18 16 15 17 19 18 15 12 10 8]; Ppv_hat = [0 0 0 0 0 0 0 1 3 5 7 8 9 8 6 4 2 0 0 0 0 0 0 0]; Load_hat = [12 11 10 9 8 8 9 12 16 18 20 22 21 20 18 17 16 15 16 17 15 13 12 11];

实际运行中,不确定参数可以平移缩放。这里设定每个时刻的风光预测误差为预测值的 ±20%,负荷为 ±10%。你的实际数据可以从微电网历史运行曲线或公开数据集取,但初次跑通用小规模随机数足够。

4.2 第一阶段变量与基础约束初始化

第一阶段变量包括每个时刻的机组启停状态、机组出力基准、与大电网交互功率、储能SOC基准。这些变量绑定在一起,形成主问题的基础约束:

T = 24; Pgub = 10; Pglb = 1; % 机组出力上下限 Ramp = 3; % 爬坡限制 Ecap = 2; Pch_max = 0.5; Pdis_max = 0.5; eta = 0.95; SOC_min = 0.2; SOC_max = 0.9; SOC0 = 0.5; % 第一阶段 sdpvar Pg(T,1) u_on(T,1) Pgrid_buy(T,1) Pgrid_sell(T,1) SOC(T,1); Pg_low = Pgub*u_on; Pg_high = Pglb*u_on; Constraint_basic = []; for t=1:T Constraint_basic = [Constraint_basic, 0 <= Pgrid_buy(t) <= 1, 0 <= Pgrid_sell(t) <= 1, Pg_low(t) <= Pg(t) <= Pg_high(t), SOC_min*Ecap <= SOC(t) <= SOC_max*Ecap, SOC(t) == (t==1 ? SOC0*Ecap : SOC(t-1)) + eta*Pch(t) - (1/eta)*Pdis(t)]; end

注意这里的SOC约束是基础形式,实际需要线性化储能充放电功率的互斥关系。常见做法是引入二进制变量表示充放电状态,或者在SOC变化率中加入松弛。我们先用简单的“SOC变化 = 充电功率 - 放电功率”线性化,其中充放电功率本身属于第二阶段调整变量。

4.3 子问题构建与最坏场景更新

子问题固定第一阶段变量后,决策第二阶段的机组微调、储能功率、弃风切荷,同时要对抗最坏的风光场景。子问题的目标是最小化调整成本,包括微调成本、弃风惩罚、切负荷惩罚。若对偶困难,也可以用枚举法:先采样若干典型极端场景,对每个场景求解一次确定性优化,取成本最大者作为最坏场景。这不是严格鲁棒,但工程上快很多。

在 C&CG 中,每一轮迭代都需要把子问题对偶或枚举得到的场景加入主问题。场景用[Pwind_scene(s,:); Ppv_scene(s,:); Load_scene(s,:)]存储,主问题里增加一组sdpvar第二阶段变量,其约束中包含这些具体数值。迭代停止条件是主目标与子目标之间的间隙足够小,经验值取 0.1%~1%。

4.4 收敛判断与结果输出

C&CG 的收敛性很好,通常 3~5 轮即可收敛到 0.1% 以内。你可以在循环里记录每轮的 lower bound(主问题目标)和 upper bound(子问题目标)。判断条件:

gap = abs(upper - lower) / abs(lower); if gap < 0.001 break; end if iter > 20 warning('CCG not converge within 20 iterations'); break; end

收敛后,主问题解中的第一阶段变量就是你要的日前计划,子问题在最坏场景下的调整量是实时运行策略。我习惯把结果画成三张图:机组出力曲线、储能SOC轨迹、电网交换功率。特别关注储能SOC是否频繁到达边界,经常触界说明不确定性预算设太大或储能容量不足。

5. 避坑专题:YALMIP+CPLEX 跑两阶段鲁棒的 5 个真实翻车点

5.1 翻车一:YALMIP 对不确定变量的约束顺序导致结果完全错误

现象:同样的模型,把不确定性约束写在普通约束之后,目标值突然飙升,甚至无解。

原因:YALMIP 的uncertain变量声明必须在所有约束建立前完成,且不确定约束不能与确定性约束混写。它解析时会把含不确定变量的约束做对等转化,如果顺序不对,部分约束被当作普通非线性约束处理,CPLEX 自然解不了。

解决:代码结构严格按“声明不确定变量 → 定义不确定性集合 → 声明普通变量 → 写普通约束 → 写含不确定量的约束”顺序。我习惯把不确定约束单独放到一个 cell 数组里,最后再拼接。

5.2 翻车二:子问题对偶后出现双线性项,无法求解

现象:dualize之后,目标函数中出现变量 (u) 与对偶变量 (\lambda) 的乘积,CPLEX 报“bilinear objective”。

原因:内层 min 的对偶变量与不确定变量正好相乘,形成对偶目标的二次项,而原问题是线性的,这种双线性来自对偶转化不当。往往是因为子问题中约束右侧包含不确定量,却没有在dualize前显式提取出来。

解决:先做变量替换,把不确定量从约束右侧移到约束左边。例如约束 (a x \le u) 写成 (a x - u \le 0),这样对偶后不确定量进入约束矩阵,不会与对偶变量在目标中直接相乘。另一种更省事的方法是采用“场景枚举法”计算最坏场景,跳过显式对偶,适合变量少的算例。

5.3 翻车三:MP 一直不收敛,lower bound 不上升

现象:主问题目标每轮几乎不变,gap 卡在很大值。

原因:主问题中的第二阶段变量没有和eta建立实际成本联系。eta仅约束为大于等于子问题目标,但主问题内新增场景的变量成本项没有正确计入,导致 MP 的目标永远等于第一阶段成本加一个虚高的eta。

解决:在主问题中,对每个已加入场景写出第二阶段成本项,像这样:

for s = 1:nscenario Constraints_MP = [Constraints_MP, ... eta >= cost_adjust(:,s)' * y(:,s)]; end

同时,目标函数里加入一个极小的正则项 (\varepsilon \sum y^2) 或直接让eta参与最小化。我一般会在每轮打印出eta和新增场景目标值,快速定位。

5.4 翻车四:CPLEX 提示内存不足,或者求解卡在 0 秒无进展

现象:模型规模只有几百个变量,但 CPLEX 内存占用极大,求解时间几十秒甚至几分钟。

原因:YALMIP 默认会生成稠密约束矩阵,尤其是sdpvar标量循环拼接会产生大量附加辅助变量。另外,如果你把cplex.mip.strategy.names=1打开了变量命名,也会拖慢速度。

解决:用矩阵化构建约束,避免for循环逐条拼接;使用repmat和sum表达式减少辅助变量。把cplex.mip.strategy.names设为 0。如果问题实在大,尝试先用确定性等价模型跑一遍,确信逻辑正确后再上鲁棒。

5.5 翻车五:结果对不确定性预算不敏感,成本曲线平得像直线

现象:(\Gamma) 从0扫到1,总成本只变化了不到 0.1%,显然不符合两阶段鲁棒优化的常识。

原因:第二阶段变量往往有足够的自由度来吸收不确定性,导致最坏场景成本等于正常场景成本。比如储能容量过大或机组爬坡能力强,鲁棒优化自然不需要额外的保守性。这不一定是代码错,可能是你的约束太松。

解决:检查第二阶段变量的费用惩罚是否足够大,尤其是弃风和切负荷的惩罚单价,如果惩罚过低,第二阶段不会使用这些“软措施”,也就体现不出鲁棒代价。我常用的惩罚价格为:弃风 400 元/MWh、切负荷 2000 元/MWh、储能寿命折算 80 元/MWh。把这些惩罚设到合理量级后,曲线就会体现出典型的鲁棒经济性。

另外一个隐蔽坑:不确定性集合的单位。如果你的负荷预测误差用了相对标幺值,而调度变量用的是有名值,就会导致误差幅度小到无法影响解。务必在建模最前面统一单位(MW 或 p.u.),我一般先把所有数据转换为以 1MW 为基准的有名值。

6. 进阶验证:参数敏感性分析与工程落地习惯

跑通一次两阶段鲁棒优化不算完,验证的方法和习惯决定最终论文或项目能否站住脚。这里给你两条进阶建议。

第一,做 (\Gamma)-成本敏感性分析,并且画成曲线。(\Gamma) 从 0 到最大值(比如24)每 0.5 算一个点,统计总成本、最坏场景成本、储能充放电次数。正常情况下总成本随 (\Gamma) 单调上升,且上升速度逐渐变缓。如果出现平台期,说明系统对不确定性有足够吸收能力;如果出现跳变,说明某个约束在特定场景下被激活,需要进一步分析是哪个约束。这种曲线是审稿人和专家最认可的一种验证方式。第二,对比两阶段鲁棒与确定性模型在真实历史数据上的表现:将历史数据中的极端日挑出来,分别用鲁棒方案和确定性方案去跑事后模拟,比较日运行成本与违约次数。鲁棒方案应该只多付出 2%~5% 的日均成本,但能避免极端日的大幅切负荷。如果溢价太高,可以调低 (\Gamma);如果仍然高风险,则要考虑升级设备容量而不是继续调参。

还有一个小技巧:用相同的 YALMIP 模型,将不确定集合缩小到原点((\Gamma=0, \Delta=0)),看解是否等于确定性优化。如果不等,说明你的模型中有条件表达式依赖不确定性变量,导致 YALMIP 对等转化走了弯路。这个“退化测试”能在十分钟内发现大部分建模 bug。

控制求解时间也有讲究。对 24 时段的微电网,C&CG 子问题每次只有几百个变量,CPLEX 通常是秒级求解。如果超过 10 秒,检查是不是把整数变量放进了子问题——如果子问题中存在储能充放电状态二进制变量,且没有固定,对偶框架会失效,速度也会爆炸。此时要么把储能二进制变量移到主问题,要么在子问题中先用启发式固定。这个决策要在建模阶段就确定,不要到求解阶段再改。

另一个工程习惯是把所有参数都设置成易于修改的结构体,例如params.delta_wind = 0.2; params.Gamma = 1.2;,这样每次测试不需要改动核心代码。我在实际项目中会把 C&CG 循环封装成一个函数,输入是params和data,输出是解和收敛曲线。这样做的好处是你可以在不同算例之间无缝切换,无论是 5 节点还是 33 节点微电网,只需换data。

最后,我自己的教训是:永远是先跑通确定性模型,再加鲁棒。确定性模型里暴露的所有问题(如负荷不平衡、SOC越界)如果不解决,加鲁棒后会更难 debug。因为鲁棒优化中的错误往往表现为“合理但略微异常”的结果,而不是明显的崩溃错误,非常难查。所以当你看到本文标题,第一反应不要是去找代码,而是先在自己的确定性调频模型里把每个变量的物理含义吃透,然后照上面的四步法手动搭建 C&CG。等你看到第一轮收敛时,那种感觉一定很爽。希望帮到你。

本文还有配套的精品资源,点击获取

需要专业的网站建设服务?

联系我们获取免费的网站建设咨询和方案报价,让我们帮助您实现业务目标

立即咨询