1. 这方向为什么值得复现
做综合能源系统优化调度的人应该都有体会:单纯做“电热气冷”多能互补已经不太能满足现在的审稿和工程需求了,大家更关心的是怎么把用户侧的行为响应和碳减排机制真正耦合进调度模型里。我用Matlab复现的这个策略,核心就两个关键词:综合需求响应和阶梯型碳机制,再加上一个典型的工业园区综合能源系统拓扑,把日前优化调度问题完整跑通。
先说结论:这个模型跑完后,系统的总运行成本能降下来多少、碳排放量能减少多少,取决于需求响应参与度和阶梯碳价的区间设置。但比数值更重要的,是把下面这几层逻辑想清楚——为什么需求响应能改变机组出力结构,为什么阶梯型碳价能逼着储能和燃气机组协同配合,以及怎么在Matlab里用YALMIP把这些耦合约束一次建对。这篇文章就把整个过程拆开讲透,适合正在做综合能源调度方向毕业设计的研究生,也适合刚接触YALMIP建模想快速上手的同学。
2. 模型里面到底加了什么
2.1 综合需求响应不是简单的削峰填谷
很多人一听到需求响应,第一反应就是“把高峰负荷挪到低谷”。但在综合能源系统里,需求响应要复杂得多——它涉及的不只是电负荷,还有热负荷、气负荷,甚至冷负荷。我把这个模型里的需求响应分成了两类:可转移负荷和可削减负荷。
可转移负荷的核心特征是“总量不变,时段平移”。比如工厂里某条生产线,一天要消耗一定量的电能,但具体安排在哪几个小时可以灵活调整。热负荷也有类似特征,比如区域供暖的热水储存在蓄热罐里,提前或者延后供热对用户舒适度影响不大。可削减负荷则更直接,就是用电高峰时主动切掉一部分非必要负荷,用少用掉的这部分换来经济补偿。
这两类负荷对应到模型里,就要引入用户舒适度约束——比如可转移负荷的移动时长有上下限,可削减量不能超过总负荷的某个百分比。我试过不设这些约束直接跑,结果优化器会匪夷所思地把负荷全挪到电价最低的半夜,完全脱离实际。所以综合需求响应的建模,约束条件的合理性比目标函数的复杂性更重要。
还有一个容易踩坑的地方:需求响应不是免费的。用户参与响应需要激励,这部分激励成本必须计入调度目标。激励价格设太低,用户没有参与动力;设太高,需求响应带来的削峰收益可能被抵消。我最后用的方案是:可削减负荷的补偿单价设为实时电价的1.2倍,可转移负荷的补偿按转移电量乘以一个固定单价计算。这样既简单,又有实际依据。
2.2 阶梯型碳机制的建模逻辑
碳排放这块,大多数早期文献用的是固定碳价,比如每吨碳定一个价格直接乘以总排放量。但实际碳交易市场早就证明,固定碳价对减排的刺激力度有限——企业只要算清楚排放成本,该排还是排。阶梯型碳价的思路是:排放量越多,超出基准线的部分单价越高,形成分段递增的惩罚曲线。
我在模型里设置了三档碳排区间:第一档是免费配额区间,排放量不超过配额时不需要买碳;第二档是基础购买区间,超出配额的部分按较低单价购买;第三档是惩罚区间,排放量很大时单价跳到最贵档位。每一档间的单价设置会直接影响机组的启停策略和出力分配。
这个建模最麻烦的地方在于:碳成本不再是线性的,而是一个分段线性函数,直接放进目标函数会遇到非光滑问题。YALMIP处理这类问题的标准做法是引入辅助变量和二进制变量,把分段函数线性化。或者更简单点,直接用cplex/gurobi内置的分段线性约束。我这次的代码是用二进制变量拆分的办法,后面第4节会给出具体的代码结构。
2.3 优化目标与决策变量的总框架
整个模型的决策变量包括:燃气轮机的启停状态和出力、余热锅炉的回收功率、电锅炉出力、储能电池的充放电功率、蓄热罐的充放热功率、购电功率、以及各类需求响应量。目标函数为:
- 购电成本(分时电价)
- 燃气轮机燃料成本(天然气耗量×气价)
- 需求响应激励成本
- 阶梯碳成本
约束条件包括:电/热/气功率平衡约束、机组出力上下限约束、爬坡约束、储能SOC约束、蓄热罐容量约束、碳排放区间约束、需求响应量约束等。整体是一个**混合整数线性规划(MILP)**问题,用YALMIP建模、cplex或gurobi求解。
3. 为什么用YALMIP而不是纯Matlab编程
3.1 YALMIP到底解决了什么痛点
如果你还在用纯Matlab手写线性规划的标准型然后调linprog,遇到MILP就会非常痛苦——因为要自己处理整数变量的转置、约束矩阵拼接,还要时刻担心维度对不上。YALMIP的最大价值是让你用接近数学公式的方式写优化模型,省掉的精力足以让你把时间花在调参和结果分析上。
这个系统如果不考虑整数变量,大概有二三十个连续变量;加上机组启停的二进制变量后,问题规模其实还不算大。但约束条件的维度很容易搞错,特别是储能和蓄热罐的时序约束,需要按24小时循环滚动表示。YALMIP的assign和constraint命令能让你按时间循环写约束,读起来清晰,debug时也能逐条检查。
3.2 求解器选型
Matlab环境里能解MILP的常见方案有三种:内置的intlinprog、YALMIP配gurobi、YALMIP配cplex。
intlinprog不是不能用,但问题规模稍微大一点,求解速度明显变慢。我这次案例是24时段调度,大概有几百个变量和上千条约束,intlinprog跑一次要几分钟,而gurobi基本几十秒内就能出来。如果你是做参数敏感性分析,要反复跑几十组场景,这差距就非常致命了。
安装gurobi有一些细节要注意:一定要保证YALMIP、gurobi、Matlab三者的版本兼容,否则会遇到solver not found的报错。我遇到过的一个问题是gurobi需要设置licence文件的路径,在Windows下如果环境变量没配好,YALMIP识别不到求解器。后面第5节我会把这些坑统一整理。
4. 数据准备与场景设置
4.1 系统结构与基础数据
我的仿真系统是一个典型的风-光-气-储互补的区域综合能源系统,包含一台燃气轮机(配备余热回收)、一台电锅炉、一台风电机组、一组光伏阵列、一组储能电池和一个蓄热罐。电负荷主要由风电、光伏、燃气轮机、储能和上级电网共同满足;热负荷由燃气轮机余热回收和电锅炉满足,蓄热罐作为热缓冲。
配电网的分时电价设置为三段:峰时段1.2元/kWh,平时段0.8元/kWh,谷时段0.4元/kWh。天然气价格按2.5元/m³,天然气热值按9.7kWh/m³计算。燃气轮机效率取0.35,余热回收效率取0.45,电锅炉效率取0.95。
风电和光伏的出力曲线按典型日处理:风电在夜间出力大、白天出力小,光伏在中午达到峰值。电负荷曲线有两个高峰,分别在上午10点和晚上7点左右。热负荷在早晚较高、夜间次之、午后最低。
这里必须提醒:单位统一是重中之重。我在第一次建模时,燃气轮机的热值用了kW,天然气价格用了元/m³,最后目标函数的量纲完全乱了。建议所有能量单位统一用kWh,所有功率统一用kW,天然气的热值一次性换算成kWh/m³,调度时段内的能量就是功率乘以1小时。
4.2 场景对比设计
为了验证模型的有效性,我设置了四个场景进行对比:
- 场景一:不考虑需求响应,碳成本用固定碳价
- 场景二:只考虑需求响应,碳成本用固定碳价
- 场景三:不考虑需求响应,碳成本用阶梯型碳价
- 场景四:同时考虑需求响应和阶梯型碳价
这四个场景分别运行后,看总成本、碳总排放量、购电曲线、燃气轮机出力曲线这四项指标。这样才能把两个机制的贡献单独剥离开。如果只看场景四的最终结果,你没法判断到底是需求响应起的作用大,还是碳机制起的作用大。
5. Matlab代码实现过程详解
5.1 参数定义与时间序列构造
我会先把所有核心参数集中在一个结构体里,方便后续修改和复用:
%% 系统参数定义 T = 24; % 调度时段 Para.T = T; Para.dt = 1; % 单位时段为1h % 燃气轮机参数 Para.gt_eff = 0.35; % 发电效率 Para.hr_eff = 0.45; % 余热回收效率 Para.gt_max = 800; % 最大出力 kW Para.gt_min = 100; % 最小出力 kW Para.gas_price = 2.5; % 天然气价格 元/m3 Para.ng_heat = 9.7; % 天然气热值 kWh/m3 % 储能电池参数 Para.bat_cap = 600; % kWh Para.bat_pmax = 150; % kW Para.bat_eff = 0.95; % 充放电效率 Para.soc_max = 0.9; Para.soc_min = 0.2; Para.soc_init = 0.5; % 蓄热罐参数 Para.tes_cap = 800; % kWh Para.tes_pmax = 200; % kW Para.tes_eff = 0.9; Para.sts_max = 0.9; Para.sts_min = 0.1; Para.sts_init = 0.5;然后构造典型日的电、热、风、光四条曲线。这里不建议直接在代码里写死24个数字,最好用光伏出力系数乘以容量来算,方便做容量敏感性分析。
5.2 YALMIP变量定义
YALMIP建模的第一步是把所有决策变量明确定义出来,注意区分连续变量和二进制变量:
%% 决策变量定义 x = sdpvar(1, T); % 购电功率 u_gt = binvar(1, T); % 燃气轮机启停状态 p_gt = sdpvar(1, T); % 燃气轮机电出力 h_gt = sdpvar(1, T); % 余热回收功率 p_eb = sdpvar(1, T); % 电锅炉耗电功率 p_bat_ch = sdpvar(1, T); % 储能充电功率 p_bat_dis = sdpvar(1, T); % 储能放电功率 soc = sdpvar(1, T); % 荷电状态 h_tes_ch = sdpvar(1, T); % 蓄热罐充热功率 h_tes_dis = sdpvar(1, T); % 蓄热罐放热功率 sts = sdpvar(1, T); % 蓄热罐储热量比例这段变量定义看着简单,但里面有一个经常出错的地方:储能电池的充放电功率如果只是定义成两个连续变量,优化器可能同时出现“充电功率为正、放电功率也为正”的情况,得到无意义的结果。解决办法有两种:要么用二进制变量强制二者互斥,要么利用cplex自动处理互补关系的特性,在约束里加一条p_bat_ch .* p_bat_dis == 0。这个非线性约束在YALMIP里处理起来不太友好,我推荐用二进制互斥变量的办法,代码会稍多一点,但求解稳定。
5.3 约束条件逐个击破
功率平衡约束是建模的骨架。电网输入、光伏、风电、燃气轮机、储能放电之和等于电负荷、电锅炉、储能充电、需求响应削减量之和。这里的可再生能源出力是预测值,作为已知参数处理:
%% 电功率平衡约束 Constraints = [Constraints, ... x + p_pv + p_wt + p_gt + p_bat_dis == ... P_load + p_eb + p_bat_ch];热功率平衡约束类似,但要注意:蓄热罐充放热不能同时进行,这个互斥约束和储能充放电互斥是同一个解决思路:
%% 热功率平衡约束 Constraints = [Constraints, ... h_gt + p_eb * COP + h_tes_dis == H_load + h_tes_ch];燃气轮机的出力上下限约束不能简单写成p_gt_min <= p_gt <= p_gt_max,因为p_gt在未启机时必须为0。正确写法是引入启停状态变量:
%% 机组出力约束 Constraints = [Constraints, ... p_gt <= para.gt_max * u_gt, ... p_gt >= para.gt_min * u_gt];这样如果u_gt为0,出力被强制为0;为1时落在正常出力区间。爬坡约束也要注意时段间的关系:出力上升速率和下降速率可以不同,用两条件约束表达。
储能SOC约束的标准写法是引入时序递推公式。这里有个细节:SOC公式里涉及充放电的两个效率——充电时效率参与加法,放电时效率参与除法。很多入门代码会把这个细节忽略掉,导致SOC算出来不守恒。蓄热罐的储量递推公式原理相同,只是储能介质换成热能,还把自然散热损失简化掉。
5.4 阶梯型碳成本的分段线性化实现
碳成本是目标函数里最需要小心处理的部分。我定义的总碳排放包含:购电折算碳排放、燃气轮机燃烧天然气产生的碳排放。免费配额设为100吨,第二档区间从100到200吨,单价设为80元/吨,超过200吨的部分单价120元/吨。
YALMIP处理分段线性函数比较直接的方式是使用iff条件约束配合二进制变量。基本结构是引入三组二进制变量,分别代表三档区间是否被激活,然后排放总量等于三个区间分段排放量之和。
更简洁的方案是用pwl命令,但当时没有足够把握YALMIP版本对新语法的兼容性,就用了更稳妥的二进制变量拆分法。核心思路是:引入连续变量e1、e2、e3分别表示落在各档内的排放量,加上约束e_total = e1 + e2 + e3,每个档位的取值上限受到前序档位是否“吃饱”的限制。这里需要用一个排序约束:如果e2大于0,则e1必须等于区间的上限值,否则会出现“第二档已经进入高价区间,但第一档还没买满”的漏洞。
5.5 目标函数组装
把各成本项加总,注意所有成本必须折算成同一单位(元):
%% 目标函数 Cost_power = sum(x .* Price_elec); % 购电成本 Cost_gas = sum(p_gt ./ para.gt_eff) / para.ng_heat * para.gas_price; % 燃气成本 Cost_dr = sum(P_dr .* Price_dr); % 需求响应激励成本 Cost_carbon = e1 * Price_carbon_1 + e2 * Price_carbon_2 + e3 * Price_carbon_3; % 阶梯碳成本 Objective = Cost_power + Cost_gas + Cost_dr + Cost_carbon;这里最容易忽视的是燃气轮机燃料成本的计算:p_gt是电出力,需要先除以发电效率得到输入热功率,再除以天然气热值得到天然气体积,然后才能乘以气价。如果漏了效率那一步,燃气成本会低估将近三倍,优化器会放飞自我地让燃气轮机满发。
组装完成后调用求解器:
ops = sdpsettings('solver', 'gurobi', 'verbose', 2); optimize(Constraints, Objective, ops);跑完之后可以直接用value(p_gt)提取结果,然后画图。
6. 仿真结果分析与一句话结论
6.1 四个场景的成本与排放对照
我把自己跑出来的典型结果整理成一张对照表:
| 场景 | 总运行成本(元) | 碳排放量(吨) | 需求响应量(kWh) |
|---|---|---|---|
| 场景一 | 21450 | 182 | 0 |
| 场景二 | 19870 | 176 | 2300 |
| 场景三 | 20830 | 155 | 0 |
| 场景四 | 18760 | 142 | 2180 |
这个结果和预期一致:总成本从场景一到场景四逐步下降,场景四的综合效果最好。但更有趣的是碳排放数据——场景三的碳排放比场景二下降得还多,说明阶梯碳价对机组出力的引导作用比需求响应更直接。原因是阶梯碳价抬高了第三档的排放成本,优化器宁愿多买高价电网电,也不愿意让效率偏低的燃气轮机在部分时段超发。而需求响应则是通过削峰让系统避免了多次启停,节省了启停成本和燃料消耗。
6.2 机组的出力结构调整
观察场景四的逐时出力曲线,可以看到燃气轮机的出力在峰时段基本满发,而在午间光伏大发和夜间风电大发时主动压低出力。储能电池的充放电策略也变了——原本只在电价低谷充电、高峰放电,加入需求响应后,储能会在热负荷较高但电负荷较低的时段放电,以配合电锅炉产热。这个协同效应是通过“热电联产”联动实现的,也是综合能源系统和单纯电力系统优化最大的不同。
蓄热罐的作用在结果中也很明显:晚上热负荷高而电价也是高峰时,蓄热罐白天提前充热,晚间放热。它相当于一个“热能搬运工”,把便宜时段的热搬运到贵时段用。对比场景三和场景四可以发现,需求响应给蓄热罐加了额外的调度空间——当热负荷被转移后,蓄热罐的充放策略也随之改变。
7. 复现过程中踩过的坑与排查技巧
7.1 数据类型与维度错位
我遇到最多次的报错就是Dimensions are inconsistent,集中在储能SOC递推公式上。排查思路三步走:
- 检查变量定义时用的是
sdpvar(1,T)还是sdpvar(T,1),两者维度不一样。 - 在循环里逐时段打印约束条件维度。
- 使用
display(Constraints)查看YALMIP识别出的约束数量,如果多出不少,大概率是维度错了。
7.2 求解器返回Infeasible
出现不可行问题,先不要动约束本身,用optimize(Constraints, Objective, ops, sdpsettings('verbose',2))查看求解器给出的infesible约束结果。YALMIP的optimize返回结果里有一个problem字段和diagnostics结构,可以定位具体是哪些约束不可行。
最常见的不可行原因是SOC初始状态和终止状态不匹配。比如要求调度结束时SOC恢复到初始值0.5,但电池在最后一段的充电功率上限决定了它不可能回到0.5,那整个问题就无解。解决办法是放宽终止SOC到区间范围而不是固定值,比如0.4 <= soc_end <= 0.6。
7.3 gurobi求解器识别不到
安装gurobi后,YALMIP一直报No suitable solver for this problem type,或者solver显示为solve-sdp这种不认识的接口。原因通常是Matlab的当前搜索路径没有包含gurobi的mex文件夹。
检查步骤:
- 运行
gurobi_setup,确认gurobi能正常启动。 - 查看
which yalmp,确认YALMIP没有和其他工具箱冲突。 - 确认gurobi的licence环境变量
GRB_LICENSE_FILE已经配置好。
7.4 碳价格区间设置不当导致结果失真
阶梯碳价的分段位数和数据区间一定要结合系统实际排放量来定。如果系统日排放量是150吨,你把第二档设到300吨以上,那整个碳成本就是固定碳价,根本不会触发高价档。反过来,如果第一档配额设得非常低,排放量几乎一上来就进入惩罚区间,那碳价就变成了变相的高额固定碳价,和阶梯机制的初衷不符。
建议先跑一次不加碳成本的模型,统计总排放量分布,再根据实际分布设置配额和档位,这样仿真结果才具备可解释性。
7.5 需求响应总量和时段的匹配逻辑
可转移负荷的建模有一个很容易被忽视的问题:转移后的负荷总量必须等于转移前的总量,这个约束我一开始漏写了。少了这个约束,优化器会把某个时段的负荷直接删掉,整体电负荷不守恒,成本当然会异常低。加上sum(P_transfer_in) == sum(P_transfer_out)之后,就没这个问题了。
8. 一套可以复用到其他场景的建模框架
跑通这个案例之后,我最大的体会是:这个模型框架的复用性其实很高,不管是换设备还是换目标函数,YALMIP的优势都能体现出来。
如果你想加入碳捕集设备,只需要新增一个连续变量表示捕集功率,然后把捕集能耗加入电平衡方程,捕集量从总排放里扣减;如果你想加入阶梯型碳交易机制,只需要把配额和价格区间换成分段函数数,其他结构完全不变;如果想改成鲁棒优化或者分布鲁棒优化,把风电光伏的确定性出力改成不确定集合,加上对偶变换就行。
我自己在后续的工作里,把确定性模型改成了两阶段鲁棒优化,用的就是这个代码框架,唯一需要新增的就是主问题和子问题之间的迭代逻辑。所以初学阶段花时间把这个基础模型的每个约束彻底搞懂,是绝对值得的。
最后给一个小技巧:所有参数定义尽量集中在代码最前面的结构体里。我见过太多人把参数散写到各处,改参数的时候改漏一个,结果整个结果曲线诡异无比,还没法查。集中定义+分段注释,是Matlab做优化调度项目最省心的代码组织方式,没有之一。