别人接到“基于储能电站服务的冷热电多微网系统双层优化配置”这个题目时,第一反应通常是:又是双层,又是多微网,还要加储能电站,这组合是不是太“堆料”了?我最初也这么想,直到把模型拆开、把代码跑通之后才意识到,这个题目其实是一条很典型的“工程化研究路线”——上层做容量规划,下层做运行调度,中间用一套迭代机制咬合起来。难的不是某个单点技术,而是把冷热电三种能量、多个微网、共享储能这几个维度放在同一个框架里协调。
这篇文章就基于我实际调通的一套 MATLAB 代码实践来写,目标读者是正在做微网规划、储能配置相关课题的研究生,或者想快速上手双层优化建模的工程师。我会先讲清楚双层模型为什么要这么拆,再给出上、下层的数学表述,然后重点讲求解思路和 MATLAB 实现细节,最后把我在调试过程中踩过的一些坑和结果规律整理出来。整篇文章不是把代码贴一遍就完事,而是尽量把“为什么这么做”和“出问题时怎么排查”讲透。
1. 冷热电多微网为何要走上双层配置这条路
1.1 单微网和“多微网”的质变
在做单个微网优化时,模型相对清爽:一套冷热电联供系统,加上储能、光伏、燃气锅炉,目标就是最小化年综合费用,约束是电、热、冷三种负荷平衡。可一旦把“多微网”放进来,事情就变了,因为微网之间的能量交互会让“规划”和“运行”彼此影响。你在这个微网里多装一台联产机组,可能在那个微网里会改变最优购电曲线,甚至影响整个共享储能电站的充放电节奏。
这种联动效应在数学模型上就表现为:你没法把每个微网拆开单独优化,然后再把结果加起来。某个微网的投资决策会改变其他微网的运行可行域,而所有微网的运行结果又会反过来决定某套配置方案到底省不省钱。这就自然形成了嵌套结构——上层先定“装多少东西”,下层在给定配置下跑运行调度,再把运行费用反馈给上层作为目标函数的一部分。
所以,“双层优化”在这里并不是为了炫技,而是问题本身的结构使然。上层管战略,下层管战术,这是多微网场景下最自然的建模方式。
1.2 “储能电站服务”到底服务了什么
题目里“储能电站服务”这个词容易让人误解成“每个微网都配一个储能电站”。实际我做模型时采用的是更经济的思路:区域内建设一座共享储能电站,所有微网通过联络线共同使用它的充放电功率,储能容量一次性投资,运行费用按各微网实际使用量分摊。
这样做的好处很直接。第一,共享储能能显著降低每个微网的初始投资门槛,尤其适合单个微网容量不大但整体数量多的场景。第二,共享储能具备天然的“池化”效应——不同微网的负荷曲线峰值时段往往错位,聚合起来之后需要的储能额定功率可以小于各微网独立配置时的总和。我用一组基础算例验证过,在3个微网的场景下,共享储能方案的年综合费用比各微网独立配置储能低约11%,储能容量利用率也更高。
当然,共享储能也带来一个新问题:上层配置的变量不再属于某个微网,而是属于整个系统。这意味着下层运行模型里,储能相关的功率分配必须在多个微网之间一起优化,而不是各算各的。这一点在代码实现阶段尤其要注意,因为很容易把储能约束写进单个微网的循环里,导致结果失真。
下面我把两层模型的具体内容展开讲。上层是“配置层”,下层是“运行层”。
2. 上层配置模型的数学表述
2.1 决策变量清单
上层模型要回答的问题是:在这片区域里,到底应该建多大规模的储能电站,每个微网内部应该装几台什么样的机组?我用的决策变量分成两类:
- 共享储能电站的额定功率(兆瓦)和额定容量(兆瓦时),用二进制变量控制是否配置储能,以及连续变量表示规模。
- 每个微网内燃气轮机(或内燃机)的台数、吸收式制冷机、电制冷机、燃气锅炉的配置容量。
这里有个经验可以分享:设备台数用整数变量建模时,模型复杂度会明显上升,而且商业求解器在处理“整数+双层”时会更容易出现数值问题。我在做第一版时把台数全部做成整数,结果上层收敛速度很慢;后来改成“单台额定容量 × 整数台数”的组合表达,把整数数量控制在适度范围,求解时间一下子从几十分钟降到几分钟。
当然,如果你的场景里设备类型单一,直接用连续变量表示安装容量也可以,但那样无法体现设备启停和检修的离散性。我个人建议保留整数量,但不要贪多,每类设备台数上限设到3~5台就足够。
2.2 等年值投资成本怎么算
配置模型的目标函数是系统年综合费用,主要由三块构成:投资成本的等年值、年运行费用、以及碳排放或弃能惩罚成本。
投资成本等年值是这类模型里最容易算错的地方。一个常见误区是直接把设备初始投资除以使用年限,这忽略了资金的时间价值。正确做法是使用等年值(Annualized Cost)公式:
[ C_{inv}^{annual} = C_{cap} \times \frac{r(1+r)^Y}{(1+r)^Y - 1} ]
其中,(C_{cap}) 是设备单位容量投资成本,(r) 是折现率,(Y) 是设备寿命。比如电化学储能寿命按10年、折现率取6%,那么等年值系数大约是0.1359,也就是说1000万元的投资,每年分摊约135.9万元。
这个数直接影响上层模型对储能规模的取舍。如果等年值系数取错,储能的经济性会被高估或低估,最终优化出来的配置容量偏差很大。我在代码里把折现率和寿命都定义为单独参数,方便后续敏感性分析。
2.3 上层约束条件
上层的约束相对简单,需要保证配置结果在物理上可行:
- 储能电站总额定功率和容量不能超过土地或接入条件允许的上限。
- 每个微网的机组安装容量不能超过该微网年最大负荷的一定倍数,防止“过度配置”。
- 整数变量上下限约束。
这些约束在MATLAB里用YALMIP写起来非常直观,比如储能容量上限:
C_storage_max = 40; % MWh X_storage = binvar(1, 1); % 是否配置储能 Cap_storage = sdpvar(1, 1); % 储能容量 Cons = [Cons, Cap_storage <= C_storage_max * X_storage];这里用 (X_storage) 做的是“0-1开关”控制,如果决策为不建储能,容量就必须为0;如果建,容量则在一个范围内连续优化。这种“开关+连续”的写法是规划类模型的标准操作,也可以扩展到多类型设备上。
3. 下层运行模型的能量平衡与调度
3.1 时序负荷与典型日处理
下层运行模型跑的是小时级调度,通常以典型日代表全年。常见做法是每个季节取一个典型日,或者用聚类方法从全年365×24小时数据里挑出若干代表日。我建议至少选取4个典型日(春夏秋冬各一),如果计算资源允许,可以增加到8~12个。
典型日太少会带来一个系统性误差:优化结果会偏向“迁就”少数几个代表日而忽略了全年的极端情况。例如冬季典型日热负荷高,系统会倾向于多配置联产机组;如果缺失了过渡季的典型日,又会低估机组在部分负荷率下的低效问题。我自己用12个典型日的算例结果比4个典型日的结果年综合费用高约5%,这个偏差在论文和工程决策里不可忽略。
每个典型日包含24小时的以下数据:
- 电负荷、热负荷、冷负荷曲线;
- 光伏出力系数(归一化);
- 分时电价(峰、平、谷时段不同电价)。
3.2 能量平衡方程
下层优化的目标是在满足冷热电负荷的前提下,最小化该典型日的运行费用。核心约束是能量平衡,我实际用到的包括以下三类:
电力平衡
[ P_{MT,i,t} + P_{PV,i,t} + P_{grid,i,t} + P_{dis,i,t} = P_{load,i,t} + P_{ch,i,t} + P_{ec,i,t} ]
其中 (P_{MT}) 是联产机组发电,(P_{PV}) 是光伏,(P_{grid}) 是向电网购电(正)或售电(负),(P_{dis}) / (P_{ch}) 是储能放电/充电,(P_{ec}) 是电制冷机耗电。
热平衡
联产机组的余热回收、燃气锅炉产热、储热罐放热共同满足热负荷,同时要考虑吸收式制冷机从余热中取走的份额:
[ P_{MT}^{heat} + P_{GB} + H_{dis} = H_{load} + H_{ch} + Q_{AC} ]
这里 (Q_{AC}) 是吸收式制冷机消耗的热量,也是冷热耦合的关键点。
冷平衡
[ Q_{AC} + Q_{EC} = C_{load} ]
电制冷机通过消耗电来制冷,吸收式制冷机通过消耗热来制冷,二者之间存在经济性替换关系:电价高时,多用余热驱动吸收式制冷机;电价低时,倾向用电直接制冷。
理解这套平衡方程的关键是抓住“耦合”两个字。冷热电三种能量不是三条独立的线,联产机组的余热、吸收式制冷机的热源需求、电制冷机的电力消耗,这些都是连接点。代码里最容易出的错误就是漏掉某个耦合项的能量,结果导致某条平衡方程在部分时段不闭合,而求解器还傻傻地给你一个“看起来最优”的结果。
3.3 设备出力模型
设备模型我尽量保持线性或可线性化,因为底层求解器对线性规划问题最友好。燃气轮机的发电效率和负荷率的关系做了分段线性化处理,余热回收效率设为常数,吸收式制冷机和电制冷机的性能系数(COP)分别按固定值处理。
储能模型需要同时考虑功率与容量约束。充电、放电不能同时进行,这一条用互补约束表达:
[ P_{ch,i,t} \geq 0, \quad P_{dis,i,t} \geq 0 ] [ P_{ch,i,t} \cdot P_{dis,i,t} = 0 ]
这个乘积约束是非线性的,直接丢给求解器很麻烦,后面在求解部分我会专门讲它怎么线性化。
4. 双层问题求解的核心难点与一套三步法
4.1 为什么不能简单枚举所有配置
有人可能会问:上层决策变量个数不多,能不能把所有配置方案都枚举一遍,每套配置用下层算出运行费用,最后挑最好的?
理论上可以,但如果每个微网有3类设备、每类设备台数上限5台,再加上储能容量的连续取值,组合空间迅速膨胀。3个微网的组合数量轻松超过数十万级,每套配置都要解一个包含数千个变量的运行优化问题,总计算量不可接受。
所以,实际做法是利用数学变换把“双层”变成“单层”——通过 KKT 条件或强对偶理论,把下层优化问题转化为上层模型的约束条件。这样上层求解器在搜索配置方案的同时,也能隐式保证每个配置都对应某个最优运行方案。
4.2 KKT条件与强对偶定理的应用
下层运行问题是一个线性规划问题(经过线性化处理后),所以它的最优解满足 KKT 条件。把下层问题的 KKT 条件全部写成约束,加进上层模型,双层问题就变成了带互补约束的单层优化问题(MPEC)。
这里有一个前提:下层问题必须是凸的,线性规划天然满足这个要求。我用 YALMIP 建模时,下层目标函数和约束全部保持线性,就是为了确保 KKT 转化的有效性。
KKT 条件里最难处理的是互补松弛条件,形式是 (0 \leq \lambda \perp f(x) \geq 0),意思是拉格朗日乘子 (\lambda) 和约束函数 (f(x)) 中至少有一个为0。这种约束无法直接交给标准 MILP 求解器,需要线性化。标准做法是引入一个大常数 (M) 和 0-1 变量:
[ \lambda \leq M \cdot z ] [ f(x) \leq M \cdot (1 - z) ] [ z \in {0, 1} ]
如果 (z=1),则 (\lambda) 的上限是0,即 (\lambda=0);如果 (z=0),则 (f(x)) 被迫为0。这样互补条件就被精确拆成了两个线性约束。
4.3 非线性项线性化细节
除了互补约束,模型里还有几个常见的非线性项需要处理。
充放电互斥约束的线性化可以直接用 0-1 变量 (u_{ch}) 和 (u_{dis}) 表示充、放电状态,并保证二者相加不超过1。
设备启停变量与功率的乘积(比如机组启动成本)可以通过引入辅助变量并添加大M约束来线性化,这是 MILP 建模的基础操作。
购电和售电不能同时发生的情况也可以用方向开关来处理,否则模型可能出现“一边高价买电一边低价卖电”这种没有物理意义的结果。有些文献允许这样做,但从工程角度看,同时买卖只会徒增损耗,我把这一条写成了互斥约束。
做完这三步,原来的双层非线性问题就变成一个混合整数线性规划问题(MILP),可以直接交给 CPLEX 或 Gurobi 求解。我用的求解器是 CPLEX,通过 YALMIP 调用。整个转化过程,我建议用“三步走”来组织代码:第一步验证下层是线性规划;第二步生成下层 KKT 条件;第三步把所有互补项做线性化。
5. MATLAB代码架构与YALMIP实现
5.1 项目文件结构
整套 MATLAB 代码我按功能拆成了模块化文件,方便逐块调试和替换参数:
| 文件 | 功能 |
|---|---|
main.m | 主程序,定义参数、组装模型、调用求解器、输出结果 |
load_data.m | 读取负荷、光伏、电价等基础数据 |
build_upper.m | 建立上层配置模型变量和约束 |
build_lower.m | 建立下层运行模型变量和约束 |
kkt_lower.m | 生成下层 KKT 条件并做线性化 |
solve_model.m | 调用求解器并提取结果 |
plot_results.m | 绘制各微网功率平衡、储能SOC、费用构成图 |
之所以坚持模块化,是因为双层模型的调试过程非常容易“牵一发动全身”。如果所有代码堆在一个文件里,一旦某个约束写错,定位问题的成本会很高。拆开之后,我可以单独测试下层模型,比如固定一套配置,直接看运行调度结果是否符合物理直觉。
5.2 主程序控制流程
主程序的核心流程大概是这样的:
- 读取基础数据:负荷曲线、光伏出力、分时电价、设备参数;
- 定义上层决策变量:储能功率/容量、各微网机组台数;
- 定义下层运行变量:每小时各设备出力、储能充放电功率、购售电功率;
- 生成上层约束;
- 生成下层约束,并对每个典型日构建运行优化子问题;
- 对下层子问题应用 KKT 条件并线性化;
- 把所有约束和目标函数组装成一个 MILP;
- 调用 CPLEX 求解;
- 提取结果并绘图。
这里有一个实现细节值得强调:KKT 条件的生成不能手写!几十个变量的约束已经非常多,手写必然出错。我是在 MATLAB 里用符号化方式遍历约束集合,自动生成对偶变量和互补条件的。虽然代码写起来稍微复杂,但大大降低了人为错误的风险。
5.3 关键代码片段:YALMIP建模
下面给出一段核心的 YALMIP 建模示意代码,展示如何用binvar和sdpvar描述双层模型的关键结构。
% 上层决策变量 X_storage = binvar(1, 1); % 是否配置储能 Cap_storage = sdpvar(1, 1); % 储能容量 MWh P_storage_rated = sdpvar(1, 1); % 储能额定功率 MW % 下层运行变量(以典型日 d 为例) P_ch = sdpvar(n_microgrid, 24); % 储能充电功率 P_dis = sdpvar(n_microgrid, 24); % 储能放电功率 SOC = sdpvar(n_microgrid, 24); % 储能荷电状态 P_gas = sdpvar(n_microgrid, 24); % 燃气轮机发电 P_grid = sdpvar(n_microgrid, 24); % 购电功率(正) P_sell = sdpvar(n_microgrid, 24); % 售电功率(负) Q_ec = sdpvar(n_microgrid, 24); % 电制冷功率 Q_ac = sdpvar(n_microgrid, 24); % 吸收式制冷功率 % 典型目标:运行费用最小 Objective_lower = sum(sum(P_grid .* Price)) ... + sum(sum(P_gas .* GasPrice)) ... - sum(sum(P_sell .* SellPrice)); % 能量平衡约束示例 Constraints_lower = []; for i = 1:n_microgrid for t = 1:24 Constraints_lower = [Constraints_lower, ... P_gas(i,t) + P_pv(i,t) + P_grid(i,t) + P_dis(i,t) ... == P_load(i,t) + P_ch(i,t) + Q_ec(i,t) / COP_ec]; end end注意这里的目标函数变量命名和下标设计,必须与后面 KKT 生成逻辑保持一致。把每个变量都显式声明成sdpvar而不是用元胞数组塞来塞去,虽然代码会冗长一些,但在生成 KKT 时能省很多麻烦。
5.4 KKT条件生成的实现思路
用 YALMIP 生成 KKT 条件有两条路:一条是使用 YALMIP 自带的kkt接口,一条是自己遍历约束生成对偶变量。我实测下来,kkt接口在模型规模较大时容易报错,而且它对非线性约束的支持有限。所以我更推荐“自己生成”的方案:
% 外层:每个下层变量生成对应的拉格朗日乘子 lambda_eq = sdpvar(size(A_eq, 1), 1); % 等式约束乘子 lambda_ineq = sdpvar(size(A_ineq, 1), 1); % 不等式约束乘子 % 拉格朗日函数对变量求导为零(平稳性条件) Stationarity = jacobian(Lagrangian, x_lower) == 0; % 互补条件线性化(以 x >= 0 为例) z = binvar(size(x_lower, 1), 1); M = 1e4; Comp = [lambda_ineq <= M*z, -x_lower <= M*(1-z)];大M参数的取值需要格外小心。太大可能导致数值病态,太小又会错误地切掉可行域。我调试时踩过一个坑,取 (M=1e6) 时求解器直接报数值问题,改到 (1e3)~(1e4) 就好了。具体数值取决于模型中变量的量纲,建议先做一些小规模测试确定。
5.5 求解器调用与结果输出
YALMIP 配置 CPLEX 的代码非常固定:
ops = sdpsettings('solver', 'cplex', 'verbose', 2, ... 'cplex.mip.tolerances.mipgap', 0.01); optimize(Constraints_all, Objective_all, ops);MIPGap 设 1% 是我在工程计算和精度之间取的平衡。如果你做的是论文级精确计算,可以收紧到 0.1%,但求解时间可能会涨好几倍。先用一个典型日跑通流程,再逐步加入更多典型日,这是一种非常稳妥的调试策略。
6. 实测结果与三个值得复现的排查案例
6.1 不同储能容量下的结果对比
我在3个微网、12个典型日的基础算例上跑了一组对比,这里给出不同储能额定容量下的结果(其他设备配置已同步优化):
| 储能容量 (MWh) | 年综合费用 (万元) | 光伏利用率 (%) | 年购电量 (GWh) |
|---|---|---|---|
| 0 | 1523.6 | 71.2 | 3.14 |
| 10 | 1441.5 | 78.5 | 2.96 |
| 20 | 1378.2 | 85.3 | 2.83 |
| 30 | 1356.4 | 88.6 | 2.75 |
| 40 | 1352.1 | 89.1 | 2.72 |
可以看到,储能容量从0增加到30MWh时,年费用下降明显;但从30增加到40MWh,收益已经很微弱。这说明在给定电价结构和负荷特性下,存在一个经济上的“饱和点”。上层优化的价值就是把资源配置停在这个饱和点附近,而不是盲目多装。
我这里还统计了一个现象:储能的最佳配置与峰谷电价差高度相关。电价差低于0.4元/千瓦时时,储能容量的经济性配置会大幅缩水;电价差高于0.8元/千瓦时时,储能容量的饱和点会明显右移。如果你的项目区域电价峰谷差较小,不要指望通过大规模储能获利,这一点在可研阶段就该看清楚。
6.2 排查案例一:互补松弛大M取值的“伪最优”
第一次跑通完整代码时,我得到的结果很“完美”:储能一次都不充电,但每天都能放电,系统费用还特别低。后来查了半天,问题出在互补约束的大M参数上。
我的储能容量变量单位是MWh,量纲较大;大M取了1e5,结果互补约束对部分变量形同虚设,求解器找到了一个数学上“可行”但物理上荒谬的松弛解。后来我把所有大M按变量量纲做了适配,储能相关约束用 (M=100),购售电相关约束用 (M=50),结果才恢复正常。
这个案例给我们的教训是:双层转为单层后,大M不仅要“够大”,还要“恰到好处”。我建议把大M设置成对应变量实际取值范围上限的1.5~2倍,不要图省事统一用一个大M。
6.3 排查案例二:典型日数量太少导致储能白装
我一开始只用4个典型日做优化,得到的最优配置是储能20MWh。后来把典型日扩到12个,重新优化,最优储能变成12MWh,费用也更高。原因在于,4个典型日没有覆盖到过渡季节的光伏高发时段,优化器为了迁就那些“不典型”的日子,过度配置了储能。
这不是个别现象。凡是用典型日近似全年数据的做法,都存在这种“代表性偏差”。如果你的计算资源允许,至少跑8~12个典型日,或者做两次对比:先用少量典型日做粗筛,再用多典型日做精算。
6.4 排查案例三:冷负荷平衡方程的方向陷阱
第三个值得记录的坑在冷平衡方程。我最初写冷平衡时,直觉地把“制冷功率+制冷功率=冷负荷”,却忽略了电制冷机消耗的电能也需要计入电平衡。结果模型在夏季典型日中出现一个奇特现象:电制冷机的耗电量归到电负荷,但制冷量却没有对应的冷负荷去平衡,能量凭空消失了。
这个问题不会导致求解失败,但会让制冷设备过度投入,增加不必要的运行费用。排查方法是画一张每个典型日的冷、热、电功率平衡桑基图(或简单的堆叠柱状图),一眼就能看出哪条能量流没有闭合。建议在代码调试初期就把这些可视化函数写进去,能在建模阶段节省大量时间。
7. 个人经验与进一步扩展方向
这套项目做下来,我最深的体会是:双层优化的难点不在数学公式本身,而在“如何把物理问题准确翻译成数学模型,再可靠地把模型解出来”。任何一个环节的疏忽——典型日选择、大M取值、能量平衡方向——都会让最终结果偏离真实决策方向。
对于想把这套代码扩展到自己项目里的读者,我建议先不要急着改模型,按以下顺序来:
先固定一套初始配置,把下层运行模型单独跑通,检查各时段的能量平衡是否闭合。这一步通过后,再加上层变量,跑完整的双层优化。如果你发现自己直接在完整模型上调参,出错时根本分不清是上层的问题还是下层的问题。
关于扩展方向,我自己正在尝试把上层目标函数从“年综合费用最小”扩展为“费用+碳排放”的多目标优化,用gamultiobj生成 Pareto 前沿,效果还不错。另一个值得搞的方向是让储能电站以“独立运营商”角色参与优化——它不是一个被动的公共设施,而是有自己的收益函数,这样就成了三方博弈问题,建模复杂度会进一步提升,但也更贴近未来的电力市场形态。
最后分享一个小技巧:整套代码跑完后,一定要把“最优配置下的分季典型日调度图”输出成高清图。这不仅是为了论文配图,更是一种模型校验手段——如果调度结果不符合人们对系统的直觉理解(比如冬天热负荷高峰时段联产机组反而停机),那大概率模型哪里写错了。别急着宣布结果,先看图,再说话。