前些日子有同行问我现在做综合能源系统优化运行是不是不加碳交易机制和需求响应就感觉少了点什么。这问题半对半错。错的是“跟风”对的是“这两个机制确实改变了对象系统的运行逻辑”。尤其是当你想用Matlab把一套调度优化模型完整落地会发现碳成本一旦进入目标函数机组出力的经济信号就完全变了需求响应变成可调变量后电负荷也不再是一根刚性曲线。这篇博文就针对“碳交易机制下考虑需求响应的综合能源系统优化运行Matlab代码”这条主线把模型怎么建、约束怎么列、代码怎么组织、结果怎么验一条线拆开来讲梳理出可以直接迁移到你自己项目里的思路。1. 综合能源系统优化为什么要把碳交易和需求响应拉进来1.1 传统成本调度模式的边界在哪传统IES优化调度目标函数基本就是运行成本最小上级电网购电费、天然气购气费、设备运维费零星再加点弃风弃光惩罚。这个框架本身没问题但它隐含了两个前提第一碳排放没有价格系统自然偏向选择便宜但未必清洁的供能路径第二用户负荷是刚性的不管电价怎么波动需求侧都不响应。这两个前提在现实中越来越站不住。碳排放一旦进入结算体系就变成了一项真实的运营成本你不能再让燃气轮机为了填负荷就无脑满发。需求侧也一样用户手里的电锅炉、空调、充电桩、可中断产线只要给足补偿或电价信号是有调节空间的。如果你把这两块排除在模型之外算出来的“最优”调度方案很可能在真实系统里根本执行不下去或者经济账并不成立。我调试过不少套IES模型最直观的感受是只做成本最小化时储能在低谷充电、高峰放电就结束了加入需求响应后储能开始和负荷侧协调动作再由碳排放成本参进来燃气轮机的运行区间都需要重新平衡。三个机制会相互影响所以标题把“碳交易机制”和“需求响应”并列放在“综合能源系统优化运行”前面不是修饰是问题本身的扩展。1.2 一个24小时调度场景的直观感受拿一个典型园区IES举例。系统里有光伏、燃气轮机、蓄电池、电锅炉和热负荷园区和上级电网之间有购售电通道。传统方法中某个时刻系统内部发电小于电负荷就买电有富余就卖电。这个逻辑很直接但不会回答下面这类问题高峰时段电价很高储能到底是放电满足内部负荷还是索性卖给电网如果碳价爬升燃气轮机发一度电的隐性成本变高要不要降低出力、转购绿电或者让需求侧多削一部分负荷用户侧可削减负荷的补偿成本是0.8元/kWh而电网高峰电价是1.2元/kWh调度该不该把这段负荷压下来这些问题你没法凭经验拍脑袋因为每个时段都在变光伏出力、负荷、电价、碳价都在动。你需要把它写成一个优化问题让求解器在约束空间里自动做权衡。加入碳交易和需求响应之后系统从一个“负荷跟踪型”的模型变成一个“成本信号响应型”的模型这正是这套代码和传统机组组合模型的核心差异。1.3 优化问题的数学形态变化把这三件事写进模型后问题的一般形态是这样目标函数由购电成本、天然气成本、设备运维成本、需求响应补偿成本、碳交易成本或收益共同构成约束包含电/热/冷能量平衡、设备出力上下限、爬坡率、储能SOC递推、购售电互斥、碳排放配额、负荷削减上限等。连续变量负责功率和能量0-1变量负责设备启停、购售电状态、储能充放互斥所以整体是一个混合整数线性规划也就是MILP。问题规模一般不会太大典型的调度周期是24小时或一周逐小时建模变量数在几百到几千之间。Matlab下用YALMIP作为建模层调用Gurobi或CPLEX求解是学术界和工程界都很成熟的一条路线。接下来我们就从建模颗粒度开始一步步把这套系统变成可计算的数学对象。2. 建模型的颗粒度把握设备模型、能量平衡与状态变量怎么定2.1 设备建模清单与取舍综合能源系统的设备五花八门但运行优化不是做设备级仿真不需要管电磁暂态或者热力学内部的微分方程。调度优化的核心是回答“每个小时该出力多少、存储多少、买多少、卖多少”。所以我一般按三类来收拢供能设备、储能设备、转换设备。我在常见算例里会用这样一份设备清单设备类型决策变量核心约束建模要点光伏P_pv0 ≤ P_pv ≤ 预测值按不可控电源处理相当于负负荷燃气轮机CHPP_chp、H_chp、启停0-1出力上下限、爬坡约束、热电比电热联产是IES耦合的关键燃气锅炉H_gb0 ≤ H_gb ≤ 额定容量补足热负荷缺口蓄电池P_ch、P_dis、SOC充放互斥、SOC递推、容量限制24小时周期初末电量一致蓄热罐H_st_ch、H_st_dis类似储能SOC热力惯性强调度意义重要电锅炉/热泵H_hp0 ≤ H_hp ≤ 额定容量、电热转换用耗电换热是电热耦合的另一种形态这套清单已经能覆盖园区级、社区级、微能源网级的大多数场景。你要做的不是把设备机理写得多细而是把设备之间的耦合关系抓住。CHP机组同时产电产热蓄电池和蓄热罐分别平移电、热负荷电锅炉把多余的电转换成热——这些耦合关系才是“综合能源系统”区别于单一电网调度的核心。2.2 能量平衡等式与不等式约束电、热两个能量载体在每个调度时段都要满足平衡。电平衡是光伏风电出力加上CHP发电、电网购电、储能放电等于电负荷减去需求响应削减量、储能充电、电锅炉耗电、电网售电。热平衡是CHP余热、燃气锅炉产热、蓄热罐放热等于热负荷加蓄热罐充热。写成等式约束大概长这样P_pv(t) P_chp(t) P_buy(t) P_dis(t) P_load(t) - P_dr(t) P_ch(t) P_hp(t) P_sell(t) H_chp(t) H_gb(t) H_hs_dis(t) H_load(t) H_hs_ch(t)这里有个易错点所有不等式约束比如设备出力上限、储能SOC区间、爬坡率都必须和等式中的单位保持一致。一个非常常见的错误是温度或者热量单位混用kW和kWh搞混后面整个结果都会变形。我习惯全部统一成“功率单位kW能量单位kWh时间步长1小时”这样功率乘时间就是能量简单直接。2.3 状态变量怎么定偏不把储能SOC写成标量新手最容易掉坑的地方是把储能SOC当成一个标量去做。正确做法是SOC是时间序列变量维度为T124小时的调度周期就要定义E(1)到E(25)其中E(1)和E(T1)都必须赋值代表调度周期开始和结束时的电量为同一值。否则求解器为了省钱会把最后一小时的电量全部放光得到一组物理上不合理但数学上“最优”的解。类似的充放电互斥需要用到0-1变量。蓄电池充电P_ch和放电P_dis之间必须加一个互斥约束u_ch u_dis ≤ 1同时P_ch有上限u_ch乘额定充电功率P_dis有上限u_dis乘额定放电功率。这也是典型的MILP表达方式。如果不加互斥模型会出现同时充电又放电的荒唐场景目标函数必然被这种“偷电”行为拉低。3. 碳交易与需求响应的数学化表达配额结算、阶梯碳价和可调度负荷3.1 碳排放流计算与配额结算逻辑碳交易机制落到IES里本质上是在目标函数里增加一项“碳排放资产结算”。计算碳排放和配额的方式在不同文献里有细节出入但通行做法分三步第一步核算系统总碳排放量。系统外购电量对应的上游发电碳排放按购电量和电网排放因子折算天然气消耗量乘燃料碳排放系数加上燃气锅炉和燃气轮机的天然气消耗如果系统内有绿电或者生物质排放系数按0处理。总体排放量像是E_total E_grid E_gas E_grid sum( P_buy(t) * factor_grid * dt ) E_gas sum( gas_use(t) * factor_gas )第二步计算系统允许的免费配额。很多场景是按“供能量乘配额系数”来给比如按购电量、CHP发电量和供热量的合理当量给一个免费配额E_quota。这样做的意思是系统的正常供能活动会获得一定数额的免费排放空间并不是每排放一吨都要花钱。第三步看总排放和配额的差。若总排放量低于免费配额富余配额可以在碳市场出售给系统带来收益若高于配额就要购买不足部分。这一正一负就构成了目标函数中的碳成本项if E_total E_quota: carbon_income price * (E_quota - E_total) else: carbon_cost price * (E_total - E_quota)3.2 阶梯碳价的分段线性建模现实中的碳价并不总是单一常数很多算例会做阶梯递增排放超出配额越多单位碳价越高。这样做的好处是能让模型在优化时主动控制高碳排放路径。阶梯碳价写进MILP是一个典型的分段线性成本建模把超过配额的部分拆成几段独立变量。假设超额部分分成三段第一段0到500吨单价60元/吨第二段500到1200吨单价75元/吨第三段超过1200吨的部分单价90元/吨。那么把超额排放变量E_exc拆成三个非负连续变量e1、e2、e3约束E_exc e1 e2 e3 e1 500 e2 700 碳成本 60*e1 75*e2 90*e3这个结构很干净。即使应用场景不同改区间边界和单价就行。注意不要在这里引入e1和e2之间的“优先级”逻辑分段线性建模本身不要求小段必须填满如果优化结果出现e20而e1不满额说明分段的目标斜率不是递增的那就需要检查价格大小关系。3.3 需求响应的三类负荷模型需求响应也不是一个笼统的“负荷可以变”落到可计算层面我习惯按三种可调度能力分开建模。第一类是可削减负荷。给定每个时段的可削减上限和补偿单价削减变量P_cut(t)从0到上限之间连续可调削减成本为补偿单价乘削减量。这类最常用于商业空调、照明、部分工业负荷。第二类是可转移负荷。比如某台生产设备总用电量必须满足但允许在可接受时间窗内前后移动。我一般用一个平移变量P_shift(t)允许正负然后加一个全时段总和为0的约束保证转移前后总用电量守恒同时限制每个时段转移量的上下限。这样电平衡等式里就多出P_shift项需求响应成本里则按转移量乘一个单位转移成本。第三类是可中断负荷。类似可削减但通常带有最大中断次数或最长连续中断时间的附加约束补偿单价更高。加了事件类约束后模型会多几个0-1变量不属于纯LP但仍属于MILP框架内。三类负荷的共同点是它们都以成本进入目标函数以可调变量进入平衡约束。这样需求响应就不是一个强制目标而是让求解器在“削减负荷省下的购电费”和“支付给用户的补偿成本”之间自动做权衡。4. 从参数到求解器的Matlab代码骨架变量声明、约束装配、数据接口4.1 代码目录规划和参数命名习惯一套能反复用的Matlab优化项目不建议把几千行堆在一个文件里。我的组织习惯是IES_CET_DR/ ├── data/ # 负荷、光伏、电价等输入数据 │ ├── load_profile.mat │ ├── pv_profile.mat │ └── system_params.xlsx ├── scripts/ │ ├── params_define.m # 设备参数、碳价、需求响应参数 │ ├── build_variables.m # YALMIP变量声明 │ ├── build_constraints.m # 装配等式不等式约束 │ ├── build_objective.m # 目标函数 │ └── solve_and_export.m # 求解与结果导出 └── results/ ├── schedule_result.mat └── figures/参数命名我有一个私心推荐所有功率相关参数统一带单位后缀或注释比如Pb_max表示与电网交互功率上限E_rated表示储电额定容量所有价格参数统一以“元/kWh”或“元/吨”为基准。项目做久了你会发现跨天调试时最耗时间的就是单位不一致带来的“玄学错误”。4.2 YALMIP变量与约束装配代码片段用YALMIP做建模层的核心语法很简单就是sdpvar定义连续变量、binvar定义0-1变量然后用约束表达式和for循环装配。24小时的典型变量定义如下T 24; dt 1; % 时间步长小时 Pbuy sdpvar(1, T); % 向电网购电 Psell sdpvar(1, T); % 向电网售电 Pchp sdpvar(1, T); % 燃气轮机发电 Hchp sdpvar(1, T); % 燃气轮机余热 Hgb sdpvar(1, T); % 燃气锅炉产热 Pch sdpvar(1, T); % 蓄电池充电 Pdis sdpvar(1, T); % 蓄电池放电 Ebat sdpvar(1, T1); % 蓄电池SOC时序 Pcut sdpvar(1, T); % 可削减负荷 Pshift sdpvar(1, T); % 可转移负荷净转移量 uChp binvar(1, T); % 燃气轮机启停 uBuy binvar(1, T); % 购电状态 uSell binvar(1, T); % 售电状态 uCh binvar(1, T); % 储电充电状态 uDis binvar(1, T); % 储电放电状态约束装配阶段我习惯用一个C变量累积拼接。买电卖电互斥约束的写法是C []; C [C, Pbuy 0, Pbuy Pbmax*uBuy]; C [C, Psell 0, Psell Psmax*uSell]; C [C, uBuy uSell 1];储能约束的写法要注意索引错位。Ebat的维度是T1所以SOC递推要针对t1到T写C [C, Ebat(1) Ebat0]; C [C, Ebat(T1) Ebat0]; % 周期始末一致 for t 1:T C [C, Ebat(t1) Ebat(t) Pch(t)*eta_ch*dt - Pdis(t)/eta_dis*dt]; C [C, El_min Ebat(t1) El_max]; end这里的eta_ch和eta_dis分别是充放电效率。很多文献会把充放电效率直接做一个合成效率但拆开建模更贴近电池实际特性。注意单位是功率kW乘时间h得到能量kWhSOC才对应得上。4.3 目标函数构建与求解器调用目标函数就是把各项成本加起来。购电成本是分时电价乘购电量售电收入是上网电价乘售电量天然气成本是单位气价乘耗气量运维成本按设备出力乘单位维护成本需求响应补偿按削减量、转移量乘对应单价碳交易成本用前面分段线性碳价处理。一个典型的目标函数表达式Cost_buy sum(buy_price .* Pbuy) * dt; Income_sell sum(sell_price .* Psell) * dt; Cost_gas sum(gas_price .* gas_use_total) * dt; Cost_om sum(om_chp .* Pchp) * dt sum(om_gb .* Hgb) * dt; Cost_dr sum(c_cut .* Pcut) * dt sum(c_shift .* abs(Pshift)) * dt; Cost_carbon lambda1*e1 lambda2*e2 lambda3*e3; objective Cost_buy Cost_gas Cost_om Cost_dr Cost_carbon ... - Income_sell;求解器调用是YALMIP的标准三步配置求解器参数、调用optimize、取变量值。ops sdpsettings(solver, gurobi, verbose, 2); ops.gurobi.MIPGap 1e-4; sol optimize(C, objective, ops); if sol.problem 0 Pbuy_opt value(Pbuy); Pchp_opt value(Pchp); Ebat_opt value(Ebat); else disp([求解失败错误码, num2str(sol.problem)]); end实际工作中我常用Gurobi。Gurobi在MILP场景的并行性能和数值稳定性都很靠得住如果机器没有Gurobi授权也可以用开源的SCIP或者Matlab内置的intlinprogYALMIP会自动检测可用求解器。不过做阶梯碳价和0-1变量互斥这种场景有个专业的MIP求解器会让调试效率高很多。4.4 结果导出先表格后图片求解完成后不要只盯着一个总成本数字看。我通常导出三类结果文件一是各个时段的设备出力表二是成本构成表三是关键曲线图。保存可以用writematrix或writetable到results目录曲线直接Matlab绘图。多方案场景下我会在不同case之间统一图表样式这样对比成本构成时视觉上已经能看出哪个方案碳成本降了多少、需求响应补偿占了多少。代码最后阶段再跑一个简单的合理性检查比如再算一遍电平衡余量确认每条母线每个小时都满足等式数值精度在1e-6量级以内就可以认为这一次求解是可信的。5. 跑完算例后怎么判断结果对错目标构成、削峰填谷效果和几个隐蔽错误5.1 算例方案怎么设计才有说服力不管是做研究还是做工程项目一套模型都要回答“加入碳交易和需求响应到底带来了什么改变”。我常用的做法是设置三个递进方案方案是否含碳交易是否含需求响应说明Case 1否否传统成本最小化基准Case 2是否单独观察碳价信号的作用Case 3是是完整机制下协同优化做一个20节点以内的园区案例设备参数采用典型数据光伏额定600kW燃气轮机300kW蓄电池200kWh燃气锅炉400kW电锅炉200kW分时电价峰谷比为3比1配额系数设为0.45吨/MWh初始碳价60元/吨。三个case跑完后对比总成本、碳排放量、购电曲线和需求响应补偿量。你会发现Case 2比Case 1的碳排放量明显下降代价是总成本略升Case 3则可能把总成本再拉回来一部分因为需求响应缓解了高峰购电压力。这一升一降的过程就是模型有效的直观证据。5.2 读结果的三张图和三个判断我拿到一套优化结果不会先看总成本而是按顺序看三张图。第一张是电平衡图。把光伏出力、CHP出力、购电、储能充放电、负荷响应后的曲线叠在一起逐小时验证电平衡等式成立。只要这张图里出现“某时刻负荷加上用电设备不等于供给”基本就是约束漏写或者数据口径错了。第二张是储能SOC曲线。它应当是一天一个周期起点电量等于终点电量期间有充有放波形平滑。如果SOC前十几个小时一直满着最后几小时持续放光说明你没有加周期始末一致约束如果SOC频繁跳变可能是充放互斥0-1约束出了问题。第三张是成本构成堆叠柱状图。重点看碳成本和需求响应补偿在总成本里的占比。正常情况下碳成本应当是正值如果算出来是负的大额收入先别高兴检查是否把碳排放配额设置为过大或者碳价分段边界有误。5.3 我踩过的几个隐蔽“雷”第一个坑是购售电互斥缺失。没有uBuy加uSell小于等于1之前我有一版结果里系统某个时段同时高价卖电、低价买电目标函数居然变小了因为模型在“倒买倒卖”套利。加了互斥约束后结果立刻恢复正常。第二个坑是蓄热罐初始状态没设。热惯性系统如果初始蓄热量留了很大的自由度结果会利用没有成本的初始热量来满足前几个时段的热负荷造成热平衡图完美但实际工程不可能执行。后来我会在参数里固定蓄热罐初始蓄热量或者对初始时段做热平衡校验。第三个坑是阶梯碳价的分段区间上限写反了。第一段上限500第二段上限700看起来是对的但第二段的意思是“从500到1200”所以表达式里要用第二段上限减去第一段上限写作700。这个数字其实取的是区间宽度。我初期在这类宽度定义上吃过亏现在的习惯是把每段区间宽度单独定义成一个参数比如seg2_width 700改起来方便得多。6. 把代码改到自己的场景维度、非凸、求解器与版本这些坑6.1 维度一致性是第一个拦路虎我自己接手不少别人写的Matlab代码第一头疼的就是维度错位。YALMIP对维度敏感sdpvar(1,T)和sdpvar(T,1)虽然数值内容一样但和矩阵约束拼接时可能直接报维度不匹配。推荐统一用行向量也就是sdpvar(1,T)所有参数也用1×T的行向量配合矩阵点乘和sum函数最不容易出错。另外储能SOC维度是T1而其他功率变量是T做约束时经常出现索引差一。建议所有递推类约束都用for循环写虽然慢一点但可读性强也方便排查。6.2 0-1变量带来的非凸性和求解器选择含有0-1变量的模型天然非凸LP求解器是直接拒绝的。模型跑之前我通常会先看一眼变量构成只要出现binvar就必须配置MIP求解器。Gurobi对MIPGap的控制参数很细我会设置一个合理的MIPGap比如1e-4避免默认值太松导致结果成本虚低。如果模型规模增长到几千个0-1变量求解时间会明显变长。此时不要急着加约束先看看能不能把一些互斥约束用Big-M弹性处理或者把部分0-1变量松弛成[0,1]之间的连续变量前提是目标函数会推动它取到边界。这个技巧需要反复试验但效果立竿见影。6.3 Matlab版本和YALMIP版本的兼容关系YALMIP是一个持续迭代的开源工具Matlab官方版本也在不断更新。旧版YALMIP在Matlab 2023b之后的某些环境中会报出奇怪的mex函数错误。遇到这种情况我的处理顺序是先升级YALMIP到最新发布版再检查求解器接口是否匹配。Gurobi每一代版本也会更换Matlab接口文件你需要在Gurobi安装目录里运行matlab设置脚本把路径加到matlabpath里。装好之后一定用一个小例子自测比如求解一个简单的线性规划确认YALMIP能正确识别出Gurobi并返回最优解。别等到大规模算例跑挂才发现接口没配好。6.4 参数敏感性扫描的快速实现模型稳定后做碳价、配额系数、需求响应补偿单价的敏感性分析非常常见。这个不需要改主代码把solve部分包成一个函数输入是碳价和配额系数输出是目标成本、碳排放量和DR削减量。然后用for循环在参数网格上跑一遍把结果存成表格。注意每次循环都要用yalmip(clear)清理旧的变量和约束空间否则上一次的约束会叠加到下一次求解里。我比较喜欢的一个扩展是把“配额系数从0.3到0.6”和“碳价从40到100”做成二维扫描画一个成本等高线图。那一瞬间你能非常清楚地看到哪个参数对系统运行成本影响最大也方便确认模型没有隐藏的数值病态。最后再分享一个个人体会这套代码真正难的不是把公式敲进去而是你要对每一项成本、每一条约束的物理意义有画面感。碳交易成本不是抽象的数字它对应着你实际买进卖出的配额需求响应补偿不是可有可无的项它代表用户侧真实的调节意愿。建模时多问自己一句“这个约束在真实系统里如果不加会发生什么”很多调试困难都能提前避开。把这个思路带走你手头那套IES模型会比你现在以为的更经得起推敲。 碳交易与需求响应下综合能源系统优化运行Matlab实现
前些日子有同行问我现在做综合能源系统优化运行是不是不加碳交易机制和需求响应就感觉少了点什么。这问题半对半错。错的是“跟风”对的是“这两个机制确实改变了对象系统的运行逻辑”。尤其是当你想用Matlab把一套调度优化模型完整落地会发现碳成本一旦进入目标函数机组出力的经济信号就完全变了需求响应变成可调变量后电负荷也不再是一根刚性曲线。这篇博文就针对“碳交易机制下考虑需求响应的综合能源系统优化运行Matlab代码”这条主线把模型怎么建、约束怎么列、代码怎么组织、结果怎么验一条线拆开来讲梳理出可以直接迁移到你自己项目里的思路。1. 综合能源系统优化为什么要把碳交易和需求响应拉进来1.1 传统成本调度模式的边界在哪传统IES优化调度目标函数基本就是运行成本最小上级电网购电费、天然气购气费、设备运维费零星再加点弃风弃光惩罚。这个框架本身没问题但它隐含了两个前提第一碳排放没有价格系统自然偏向选择便宜但未必清洁的供能路径第二用户负荷是刚性的不管电价怎么波动需求侧都不响应。这两个前提在现实中越来越站不住。碳排放一旦进入结算体系就变成了一项真实的运营成本你不能再让燃气轮机为了填负荷就无脑满发。需求侧也一样用户手里的电锅炉、空调、充电桩、可中断产线只要给足补偿或电价信号是有调节空间的。如果你把这两块排除在模型之外算出来的“最优”调度方案很可能在真实系统里根本执行不下去或者经济账并不成立。我调试过不少套IES模型最直观的感受是只做成本最小化时储能在低谷充电、高峰放电就结束了加入需求响应后储能开始和负荷侧协调动作再由碳排放成本参进来燃气轮机的运行区间都需要重新平衡。三个机制会相互影响所以标题把“碳交易机制”和“需求响应”并列放在“综合能源系统优化运行”前面不是修饰是问题本身的扩展。1.2 一个24小时调度场景的直观感受拿一个典型园区IES举例。系统里有光伏、燃气轮机、蓄电池、电锅炉和热负荷园区和上级电网之间有购售电通道。传统方法中某个时刻系统内部发电小于电负荷就买电有富余就卖电。这个逻辑很直接但不会回答下面这类问题高峰时段电价很高储能到底是放电满足内部负荷还是索性卖给电网如果碳价爬升燃气轮机发一度电的隐性成本变高要不要降低出力、转购绿电或者让需求侧多削一部分负荷用户侧可削减负荷的补偿成本是0.8元/kWh而电网高峰电价是1.2元/kWh调度该不该把这段负荷压下来这些问题你没法凭经验拍脑袋因为每个时段都在变光伏出力、负荷、电价、碳价都在动。你需要把它写成一个优化问题让求解器在约束空间里自动做权衡。加入碳交易和需求响应之后系统从一个“负荷跟踪型”的模型变成一个“成本信号响应型”的模型这正是这套代码和传统机组组合模型的核心差异。1.3 优化问题的数学形态变化把这三件事写进模型后问题的一般形态是这样目标函数由购电成本、天然气成本、设备运维成本、需求响应补偿成本、碳交易成本或收益共同构成约束包含电/热/冷能量平衡、设备出力上下限、爬坡率、储能SOC递推、购售电互斥、碳排放配额、负荷削减上限等。连续变量负责功率和能量0-1变量负责设备启停、购售电状态、储能充放互斥所以整体是一个混合整数线性规划也就是MILP。问题规模一般不会太大典型的调度周期是24小时或一周逐小时建模变量数在几百到几千之间。Matlab下用YALMIP作为建模层调用Gurobi或CPLEX求解是学术界和工程界都很成熟的一条路线。接下来我们就从建模颗粒度开始一步步把这套系统变成可计算的数学对象。2. 建模型的颗粒度把握设备模型、能量平衡与状态变量怎么定2.1 设备建模清单与取舍综合能源系统的设备五花八门但运行优化不是做设备级仿真不需要管电磁暂态或者热力学内部的微分方程。调度优化的核心是回答“每个小时该出力多少、存储多少、买多少、卖多少”。所以我一般按三类来收拢供能设备、储能设备、转换设备。我在常见算例里会用这样一份设备清单设备类型决策变量核心约束建模要点光伏P_pv0 ≤ P_pv ≤ 预测值按不可控电源处理相当于负负荷燃气轮机CHPP_chp、H_chp、启停0-1出力上下限、爬坡约束、热电比电热联产是IES耦合的关键燃气锅炉H_gb0 ≤ H_gb ≤ 额定容量补足热负荷缺口蓄电池P_ch、P_dis、SOC充放互斥、SOC递推、容量限制24小时周期初末电量一致蓄热罐H_st_ch、H_st_dis类似储能SOC热力惯性强调度意义重要电锅炉/热泵H_hp0 ≤ H_hp ≤ 额定容量、电热转换用耗电换热是电热耦合的另一种形态这套清单已经能覆盖园区级、社区级、微能源网级的大多数场景。你要做的不是把设备机理写得多细而是把设备之间的耦合关系抓住。CHP机组同时产电产热蓄电池和蓄热罐分别平移电、热负荷电锅炉把多余的电转换成热——这些耦合关系才是“综合能源系统”区别于单一电网调度的核心。2.2 能量平衡等式与不等式约束电、热两个能量载体在每个调度时段都要满足平衡。电平衡是光伏风电出力加上CHP发电、电网购电、储能放电等于电负荷减去需求响应削减量、储能充电、电锅炉耗电、电网售电。热平衡是CHP余热、燃气锅炉产热、蓄热罐放热等于热负荷加蓄热罐充热。写成等式约束大概长这样P_pv(t) P_chp(t) P_buy(t) P_dis(t) P_load(t) - P_dr(t) P_ch(t) P_hp(t) P_sell(t) H_chp(t) H_gb(t) H_hs_dis(t) H_load(t) H_hs_ch(t)这里有个易错点所有不等式约束比如设备出力上限、储能SOC区间、爬坡率都必须和等式中的单位保持一致。一个非常常见的错误是温度或者热量单位混用kW和kWh搞混后面整个结果都会变形。我习惯全部统一成“功率单位kW能量单位kWh时间步长1小时”这样功率乘时间就是能量简单直接。2.3 状态变量怎么定偏不把储能SOC写成标量新手最容易掉坑的地方是把储能SOC当成一个标量去做。正确做法是SOC是时间序列变量维度为T124小时的调度周期就要定义E(1)到E(25)其中E(1)和E(T1)都必须赋值代表调度周期开始和结束时的电量为同一值。否则求解器为了省钱会把最后一小时的电量全部放光得到一组物理上不合理但数学上“最优”的解。类似的充放电互斥需要用到0-1变量。蓄电池充电P_ch和放电P_dis之间必须加一个互斥约束u_ch u_dis ≤ 1同时P_ch有上限u_ch乘额定充电功率P_dis有上限u_dis乘额定放电功率。这也是典型的MILP表达方式。如果不加互斥模型会出现同时充电又放电的荒唐场景目标函数必然被这种“偷电”行为拉低。3. 碳交易与需求响应的数学化表达配额结算、阶梯碳价和可调度负荷3.1 碳排放流计算与配额结算逻辑碳交易机制落到IES里本质上是在目标函数里增加一项“碳排放资产结算”。计算碳排放和配额的方式在不同文献里有细节出入但通行做法分三步第一步核算系统总碳排放量。系统外购电量对应的上游发电碳排放按购电量和电网排放因子折算天然气消耗量乘燃料碳排放系数加上燃气锅炉和燃气轮机的天然气消耗如果系统内有绿电或者生物质排放系数按0处理。总体排放量像是E_total E_grid E_gas E_grid sum( P_buy(t) * factor_grid * dt ) E_gas sum( gas_use(t) * factor_gas )第二步计算系统允许的免费配额。很多场景是按“供能量乘配额系数”来给比如按购电量、CHP发电量和供热量的合理当量给一个免费配额E_quota。这样做的意思是系统的正常供能活动会获得一定数额的免费排放空间并不是每排放一吨都要花钱。第三步看总排放和配额的差。若总排放量低于免费配额富余配额可以在碳市场出售给系统带来收益若高于配额就要购买不足部分。这一正一负就构成了目标函数中的碳成本项if E_total E_quota: carbon_income price * (E_quota - E_total) else: carbon_cost price * (E_total - E_quota)3.2 阶梯碳价的分段线性建模现实中的碳价并不总是单一常数很多算例会做阶梯递增排放超出配额越多单位碳价越高。这样做的好处是能让模型在优化时主动控制高碳排放路径。阶梯碳价写进MILP是一个典型的分段线性成本建模把超过配额的部分拆成几段独立变量。假设超额部分分成三段第一段0到500吨单价60元/吨第二段500到1200吨单价75元/吨第三段超过1200吨的部分单价90元/吨。那么把超额排放变量E_exc拆成三个非负连续变量e1、e2、e3约束E_exc e1 e2 e3 e1 500 e2 700 碳成本 60*e1 75*e2 90*e3这个结构很干净。即使应用场景不同改区间边界和单价就行。注意不要在这里引入e1和e2之间的“优先级”逻辑分段线性建模本身不要求小段必须填满如果优化结果出现e20而e1不满额说明分段的目标斜率不是递增的那就需要检查价格大小关系。3.3 需求响应的三类负荷模型需求响应也不是一个笼统的“负荷可以变”落到可计算层面我习惯按三种可调度能力分开建模。第一类是可削减负荷。给定每个时段的可削减上限和补偿单价削减变量P_cut(t)从0到上限之间连续可调削减成本为补偿单价乘削减量。这类最常用于商业空调、照明、部分工业负荷。第二类是可转移负荷。比如某台生产设备总用电量必须满足但允许在可接受时间窗内前后移动。我一般用一个平移变量P_shift(t)允许正负然后加一个全时段总和为0的约束保证转移前后总用电量守恒同时限制每个时段转移量的上下限。这样电平衡等式里就多出P_shift项需求响应成本里则按转移量乘一个单位转移成本。第三类是可中断负荷。类似可削减但通常带有最大中断次数或最长连续中断时间的附加约束补偿单价更高。加了事件类约束后模型会多几个0-1变量不属于纯LP但仍属于MILP框架内。三类负荷的共同点是它们都以成本进入目标函数以可调变量进入平衡约束。这样需求响应就不是一个强制目标而是让求解器在“削减负荷省下的购电费”和“支付给用户的补偿成本”之间自动做权衡。4. 从参数到求解器的Matlab代码骨架变量声明、约束装配、数据接口4.1 代码目录规划和参数命名习惯一套能反复用的Matlab优化项目不建议把几千行堆在一个文件里。我的组织习惯是IES_CET_DR/ ├── data/ # 负荷、光伏、电价等输入数据 │ ├── load_profile.mat │ ├── pv_profile.mat │ └── system_params.xlsx ├── scripts/ │ ├── params_define.m # 设备参数、碳价、需求响应参数 │ ├── build_variables.m # YALMIP变量声明 │ ├── build_constraints.m # 装配等式不等式约束 │ ├── build_objective.m # 目标函数 │ └── solve_and_export.m # 求解与结果导出 └── results/ ├── schedule_result.mat └── figures/参数命名我有一个私心推荐所有功率相关参数统一带单位后缀或注释比如Pb_max表示与电网交互功率上限E_rated表示储电额定容量所有价格参数统一以“元/kWh”或“元/吨”为基准。项目做久了你会发现跨天调试时最耗时间的就是单位不一致带来的“玄学错误”。4.2 YALMIP变量与约束装配代码片段用YALMIP做建模层的核心语法很简单就是sdpvar定义连续变量、binvar定义0-1变量然后用约束表达式和for循环装配。24小时的典型变量定义如下T 24; dt 1; % 时间步长小时 Pbuy sdpvar(1, T); % 向电网购电 Psell sdpvar(1, T); % 向电网售电 Pchp sdpvar(1, T); % 燃气轮机发电 Hchp sdpvar(1, T); % 燃气轮机余热 Hgb sdpvar(1, T); % 燃气锅炉产热 Pch sdpvar(1, T); % 蓄电池充电 Pdis sdpvar(1, T); % 蓄电池放电 Ebat sdpvar(1, T1); % 蓄电池SOC时序 Pcut sdpvar(1, T); % 可削减负荷 Pshift sdpvar(1, T); % 可转移负荷净转移量 uChp binvar(1, T); % 燃气轮机启停 uBuy binvar(1, T); % 购电状态 uSell binvar(1, T); % 售电状态 uCh binvar(1, T); % 储电充电状态 uDis binvar(1, T); % 储电放电状态约束装配阶段我习惯用一个C变量累积拼接。买电卖电互斥约束的写法是C []; C [C, Pbuy 0, Pbuy Pbmax*uBuy]; C [C, Psell 0, Psell Psmax*uSell]; C [C, uBuy uSell 1];储能约束的写法要注意索引错位。Ebat的维度是T1所以SOC递推要针对t1到T写C [C, Ebat(1) Ebat0]; C [C, Ebat(T1) Ebat0]; % 周期始末一致 for t 1:T C [C, Ebat(t1) Ebat(t) Pch(t)*eta_ch*dt - Pdis(t)/eta_dis*dt]; C [C, El_min Ebat(t1) El_max]; end这里的eta_ch和eta_dis分别是充放电效率。很多文献会把充放电效率直接做一个合成效率但拆开建模更贴近电池实际特性。注意单位是功率kW乘时间h得到能量kWhSOC才对应得上。4.3 目标函数构建与求解器调用目标函数就是把各项成本加起来。购电成本是分时电价乘购电量售电收入是上网电价乘售电量天然气成本是单位气价乘耗气量运维成本按设备出力乘单位维护成本需求响应补偿按削减量、转移量乘对应单价碳交易成本用前面分段线性碳价处理。一个典型的目标函数表达式Cost_buy sum(buy_price .* Pbuy) * dt; Income_sell sum(sell_price .* Psell) * dt; Cost_gas sum(gas_price .* gas_use_total) * dt; Cost_om sum(om_chp .* Pchp) * dt sum(om_gb .* Hgb) * dt; Cost_dr sum(c_cut .* Pcut) * dt sum(c_shift .* abs(Pshift)) * dt; Cost_carbon lambda1*e1 lambda2*e2 lambda3*e3; objective Cost_buy Cost_gas Cost_om Cost_dr Cost_carbon ... - Income_sell;求解器调用是YALMIP的标准三步配置求解器参数、调用optimize、取变量值。ops sdpsettings(solver, gurobi, verbose, 2); ops.gurobi.MIPGap 1e-4; sol optimize(C, objective, ops); if sol.problem 0 Pbuy_opt value(Pbuy); Pchp_opt value(Pchp); Ebat_opt value(Ebat); else disp([求解失败错误码, num2str(sol.problem)]); end实际工作中我常用Gurobi。Gurobi在MILP场景的并行性能和数值稳定性都很靠得住如果机器没有Gurobi授权也可以用开源的SCIP或者Matlab内置的intlinprogYALMIP会自动检测可用求解器。不过做阶梯碳价和0-1变量互斥这种场景有个专业的MIP求解器会让调试效率高很多。4.4 结果导出先表格后图片求解完成后不要只盯着一个总成本数字看。我通常导出三类结果文件一是各个时段的设备出力表二是成本构成表三是关键曲线图。保存可以用writematrix或writetable到results目录曲线直接Matlab绘图。多方案场景下我会在不同case之间统一图表样式这样对比成本构成时视觉上已经能看出哪个方案碳成本降了多少、需求响应补偿占了多少。代码最后阶段再跑一个简单的合理性检查比如再算一遍电平衡余量确认每条母线每个小时都满足等式数值精度在1e-6量级以内就可以认为这一次求解是可信的。5. 跑完算例后怎么判断结果对错目标构成、削峰填谷效果和几个隐蔽错误5.1 算例方案怎么设计才有说服力不管是做研究还是做工程项目一套模型都要回答“加入碳交易和需求响应到底带来了什么改变”。我常用的做法是设置三个递进方案方案是否含碳交易是否含需求响应说明Case 1否否传统成本最小化基准Case 2是否单独观察碳价信号的作用Case 3是是完整机制下协同优化做一个20节点以内的园区案例设备参数采用典型数据光伏额定600kW燃气轮机300kW蓄电池200kWh燃气锅炉400kW电锅炉200kW分时电价峰谷比为3比1配额系数设为0.45吨/MWh初始碳价60元/吨。三个case跑完后对比总成本、碳排放量、购电曲线和需求响应补偿量。你会发现Case 2比Case 1的碳排放量明显下降代价是总成本略升Case 3则可能把总成本再拉回来一部分因为需求响应缓解了高峰购电压力。这一升一降的过程就是模型有效的直观证据。5.2 读结果的三张图和三个判断我拿到一套优化结果不会先看总成本而是按顺序看三张图。第一张是电平衡图。把光伏出力、CHP出力、购电、储能充放电、负荷响应后的曲线叠在一起逐小时验证电平衡等式成立。只要这张图里出现“某时刻负荷加上用电设备不等于供给”基本就是约束漏写或者数据口径错了。第二张是储能SOC曲线。它应当是一天一个周期起点电量等于终点电量期间有充有放波形平滑。如果SOC前十几个小时一直满着最后几小时持续放光说明你没有加周期始末一致约束如果SOC频繁跳变可能是充放互斥0-1约束出了问题。第三张是成本构成堆叠柱状图。重点看碳成本和需求响应补偿在总成本里的占比。正常情况下碳成本应当是正值如果算出来是负的大额收入先别高兴检查是否把碳排放配额设置为过大或者碳价分段边界有误。5.3 我踩过的几个隐蔽“雷”第一个坑是购售电互斥缺失。没有uBuy加uSell小于等于1之前我有一版结果里系统某个时段同时高价卖电、低价买电目标函数居然变小了因为模型在“倒买倒卖”套利。加了互斥约束后结果立刻恢复正常。第二个坑是蓄热罐初始状态没设。热惯性系统如果初始蓄热量留了很大的自由度结果会利用没有成本的初始热量来满足前几个时段的热负荷造成热平衡图完美但实际工程不可能执行。后来我会在参数里固定蓄热罐初始蓄热量或者对初始时段做热平衡校验。第三个坑是阶梯碳价的分段区间上限写反了。第一段上限500第二段上限700看起来是对的但第二段的意思是“从500到1200”所以表达式里要用第二段上限减去第一段上限写作700。这个数字其实取的是区间宽度。我初期在这类宽度定义上吃过亏现在的习惯是把每段区间宽度单独定义成一个参数比如seg2_width 700改起来方便得多。6. 把代码改到自己的场景维度、非凸、求解器与版本这些坑6.1 维度一致性是第一个拦路虎我自己接手不少别人写的Matlab代码第一头疼的就是维度错位。YALMIP对维度敏感sdpvar(1,T)和sdpvar(T,1)虽然数值内容一样但和矩阵约束拼接时可能直接报维度不匹配。推荐统一用行向量也就是sdpvar(1,T)所有参数也用1×T的行向量配合矩阵点乘和sum函数最不容易出错。另外储能SOC维度是T1而其他功率变量是T做约束时经常出现索引差一。建议所有递推类约束都用for循环写虽然慢一点但可读性强也方便排查。6.2 0-1变量带来的非凸性和求解器选择含有0-1变量的模型天然非凸LP求解器是直接拒绝的。模型跑之前我通常会先看一眼变量构成只要出现binvar就必须配置MIP求解器。Gurobi对MIPGap的控制参数很细我会设置一个合理的MIPGap比如1e-4避免默认值太松导致结果成本虚低。如果模型规模增长到几千个0-1变量求解时间会明显变长。此时不要急着加约束先看看能不能把一些互斥约束用Big-M弹性处理或者把部分0-1变量松弛成[0,1]之间的连续变量前提是目标函数会推动它取到边界。这个技巧需要反复试验但效果立竿见影。6.3 Matlab版本和YALMIP版本的兼容关系YALMIP是一个持续迭代的开源工具Matlab官方版本也在不断更新。旧版YALMIP在Matlab 2023b之后的某些环境中会报出奇怪的mex函数错误。遇到这种情况我的处理顺序是先升级YALMIP到最新发布版再检查求解器接口是否匹配。Gurobi每一代版本也会更换Matlab接口文件你需要在Gurobi安装目录里运行matlab设置脚本把路径加到matlabpath里。装好之后一定用一个小例子自测比如求解一个简单的线性规划确认YALMIP能正确识别出Gurobi并返回最优解。别等到大规模算例跑挂才发现接口没配好。6.4 参数敏感性扫描的快速实现模型稳定后做碳价、配额系数、需求响应补偿单价的敏感性分析非常常见。这个不需要改主代码把solve部分包成一个函数输入是碳价和配额系数输出是目标成本、碳排放量和DR削减量。然后用for循环在参数网格上跑一遍把结果存成表格。注意每次循环都要用yalmip(clear)清理旧的变量和约束空间否则上一次的约束会叠加到下一次求解里。我比较喜欢的一个扩展是把“配额系数从0.3到0.6”和“碳价从40到100”做成二维扫描画一个成本等高线图。那一瞬间你能非常清楚地看到哪个参数对系统运行成本影响最大也方便确认模型没有隐藏的数值病态。最后再分享一个个人体会这套代码真正难的不是把公式敲进去而是你要对每一项成本、每一条约束的物理意义有画面感。碳交易成本不是抽象的数字它对应着你实际买进卖出的配额需求响应补偿不是可有可无的项它代表用户侧真实的调节意愿。建模时多问自己一句“这个约束在真实系统里如果不加会发生什么”很多调试困难都能提前避开。把这个思路带走你手头那套IES模型会比你现在以为的更经得起推敲。