双碳目标下综合能源系统低碳优化调度：Matlab+Yalmip+Cplex实战解析

简介本资源是一份面向能源系统优化初学者的Matlab实践程序聚焦“双碳”目标下综合能源系统低碳运行调度问题涵盖光伏、风电、热电联产、燃气锅炉、电锅炉、电储能及碳捕集设备等多元组件并嵌入碳交易机制建模以系统总运行成本最小化为优化目标。压缩包为单文件rar格式4KB内含1个核心Matlab脚本.m文件采用YALMIP建模语言调用CPLEX求解器完整实现低碳约束下的多能协同调度建模与求解流程。程序注释详尽、逻辑分层清晰变量命名规范便于理解各设备出力耦合关系、碳流平衡约束及交易成本计算逻辑。目前已有2745人学习下载适合高校能源类专业本科生、研究生开展课程设计、毕业设计或科研入门训练可直接运行复现结果亦可作为拓展碳捕集与碳市场机制建模的学习范例。 双碳目标提出来之后综合能源系统Integrated Energy System, IES的低碳优化调度基本成了电力、能源相关专业硕士博士的“必修课”。我自己做这个方向做了三四年从最开始只会用Cplex解线性规划到后来把碳交易机制、阶梯碳价、P2G电转气、碳捕集这些模块一个个塞进模型里踩过的坑确实不少。这篇就把我之前做的**双碳下的综合能源系统低碳运行优化调度Matlab程序YalmipCplex求解**整个项目拿出来拆一遍从模型设计、约束构建、代码实现到求解器配置尽量把每一个关键决策背后的逻辑讲清楚。适合刚入门做IES优化调度的同学也适合那些模型建好了但求解总出问题的朋友参考。1. 综合能源系统低碳调度的核心逻辑1.1 为什么要做低碳调度低碳体现在哪里传统经济调度只盯着成本最小化目标函数里是购电费用、燃气费用、设备运维费用。但在双碳背景下光看钱不够了你得把碳排放也当成一种“代价”放进优化问题里。最直接的做法是给目标函数加一个碳排放成本项比如碳税机制下每吨CO₂价格固定碳交易机制下碳排放配额不足的部分需要去市场购买超排越多单价越贵。这个项目的核心思路其实就一句话在满足冷热电负荷需求的前提下协调各设备出力使得“运行成本碳交易成本”最小化。听起来简单但真正做起来设备建模、约束线性化、求解器配置每一个环节都有讲究。1.2 典型的IES结构组成我用的系统结构比较典型几乎涵盖了目前论文里最常见的几类设备参考结构如下能源输入侧电网购电、天然气网购气能量转换设备燃气轮机GT、燃气锅炉GB、电制冷机EC、吸收式制冷机AC、余热回收装置WH等能量存储设备蓄电池ES、蓄热罐TST等耦合环节电、热、冷三种能量通过热电联产、热泵等设备耦合在一起这里有一个特别容易忽略的点电力系统的P2G电转气设备。它可以把富余风电转化为天然气既能促进新能源消纳又能在碳交易机制下获得额外收益因为它消耗了CO₂。我后面那个版本就加了P2G结果实验效果比不加要好不少尤其是风电渗透率高的场景。1.3 模型选型的几个关键考虑做这类项目模型不是越复杂越好关键看你的研究目标和可用数据。我的经验是分三步第一确定时间尺度。一般用典型日比如冬季典型日、夏季典型日跑24小时步长1小时。也有用4个典型日分别代表四季的更精细的还有用多场景随机优化的但后两种对数据和求解时间的要求高很多。第二决定要不要考虑不确定性。不考虑的话用确定性优化就行Cplex轻松解考虑风电出力的不确定性就要用鲁棒优化、分布鲁棒或者场景法模型复杂度直接上一个台阶。我这次做的是确定性优化打底后面再往鲁棒方向扩展。第三碳排放的建模方式。这块选择直接决定了目标函数的形态。碳税机制最简单线性项碳交易机制要分超额、余量、阶梯碳价三种情形阶梯碳价是分段线性函数需要引入0-1变量处理Yalmip里用implies约束或者大M法都行。2. 目标函数设计与碳交易建模细节2.1 目标函数的总体构成我的目标函数是运行成本加碳交易成本运行成本里又细分了购电成本、购气成本、设备运维成本、弃风惩罚成本。写成表达式大概是这样的min F F_buy_e F_buy_g F_om F_curt F_co2F_buy_e从电网购电的费用分时电价下按时段单价计算F_buy_g购买天然气的费用注意燃气轮机和燃气锅炉都用同一气源时按总耗气量统一计算F_om各设备的运行维护费用一般用单位出力的运维系数乘以出力F_curt弃风惩罚费用是为了让优化结果尽量消纳风电F_co2碳排放成本这是低碳调度的灵魂2.2 碳交易成本的三级阶梯模型碳交易机制最关键的是怎么给碳排放定价。我这里用的是目前论文主流的三级阶梯碳价模型。概念解释一下首先政府给系统发一个碳排放配额实际排放量低于配额的那部分可以拿到市场上卖钱高于配额的部分要花钱买而且超得越多单价越贵。三级阶梯模型是这样处理的设定一个初始配额E0然后设三个区间——配额内、超排1倍以内、超排1倍以上。不同区间对应不同的碳价。写成数学表达就是一个分段函数我直接把分段线性化的做法整理成表格碳排放区间碳交易价格数学处理方式实际排放 ≤ 配额按奖励价购入线性项回收收益为负成本配额 实际排放 ≤ 配额×2按碳价基础价线性项0-1变量标记实际排放 配额×2按碳价倍数价分段线性需引入辅助变量实现的时候要注意一个坑分段线性函数如果直接用Yalmip的pwf函数某些老版本会生成额外的二进制变量降低求解效率。我后来改用大M法手写分段约束Cplex求解速度快了一倍多代码也更好排查。2.3 碳排放源的核算细节碳排放源通常是三个从电网购电对应的间接排放、燃气轮机燃烧天然气的直接排放、燃气锅炉燃烧天然气的直接排放。前两者的排放系数我有段时间一直用固定值后来发现不够严谨——电网排放因子在不同季节、不同负荷水平下其实有差异如果数据允许建议至少区分峰谷时段的排放因子。电转气P2G环节的碳排放核算是另一个容易出错的地方。P2G消耗CO₂合成天然气这部分CO₂从哪来如果直接从系统碳排放中扣减得明确来源。最严谨的做法是设置一个碳捕集装置CCS把烟气里的CO₂分离出来供给P2G这样系统的净碳排放量才会明显下降。如果模型里只有P2G没有CCS那就默认P2G从外部环境获取CO₂不能扣减系统排放否则模型会出现“凭空减排”的不合理结果。3. 约束条件从设备特性到系统平衡3.1 功率平衡约束约束条件里第一个要写的就是各能源母线的功率平衡。电、热、冷三条母线各一条等式约束这是保证模型物理意义正确的基础。电功率平衡P_gt(t) P_wt(t) P_buy(t) P_dis(t) P_p2g_in(t) P_load(t) P_ec(t) P_ch(t) P_p2g_out(t)热功率平衡Q_gt_recover(t) Q_gb(t) Q_dis_heat(t) Q_load(t) Q_ac(t) Q_ch_heat(t)冷功率平衡Q_ec(t) Q_ac(t) C_load(t)这里想提醒一点P2G设备是耗电的它的输入电功率在电平衡里应该作为负荷项放等号右侧它的输出是天然气会影响气平衡。不少初学的朋友把P2G的输入输出搞反了日志里怎么看怎么不对劲。3.2 设备出力上下限与爬坡约束每台设备都要有出力上下限约束这是最基础的。燃气轮机和燃气锅炉还得加爬坡约束——因为天然气机组响应速度有限相邻时段出力变化不能太大。爬坡约束的表达式是-P_gt_ramp ≤ P_gt(t) - P_gt(t-1) ≤ P_gt_ramp这个约束在Yalmip里直接用循环写就行24时段跑下来速度不受影响。不过要注意t1时需要单独处理要么用初始状态要么特殊设定第一个时段的爬坡量为0。3.3 储能设备约束与SOC连续性蓄电池和蓄热罐的建模也是个性化很强的环节。核心约束是SOC荷电状态/蓄热状态的递推关系SOC(t) SOC(t-1) η_ch * P_ch(t) - P_dis(t) / η_dis这里特别注意充放电效率的方向。充电时实际存入电池的能量等于充电功率乘以充电效率放电时实际放出的能量等于放电功率除以放电效率这两个效率如果搞反了能量就不守恒了。还要注意几个隐含约束同一个时刻不能同时充放电。这个“互斥约束”可以用二元变量处理P_ch(t) ≤ M * u(t) P_dis(t) ≤ M * (1 - u(t))如果嫌麻烦也有人不设互斥约束靠分时电价自然引导——电价低时充电、电价高时放电两者同时发生的优化解几乎没有但这种做法在学术上不够严谨。我建议还是老老实实加二元变量一个u(t)的事求解器也不怕。3.4 燃气轮机热电比约束燃气轮机有热电联产特性它的电出力和热回收出力之间不是完全独立的通常有一个热电比范围。简化版模型里可以用固定热电比即Q_gt_recover(t) α * P_gt(t)但更精细的模型应该考虑热电比的可调区间。可调热电比约束写成α_min * P_gt(t) ≤ Q_gt_recover(t) ≤ α_max * P_gt(t)这个约束是线性的加入后模型求解不受影响但设备运行区域的刻画更真实了。特别是夏季热负荷很低、冬季热负荷很高的场景固定热电比会导致约束过紧或过松调度结果偏离实际。4. MatlabYalmipCplex环境配置与代码实现4.1 Yalmip安装与Cplex配置这套工具链的配置踩过坑的人不少我先把最省事的路线说清楚。第一步安装Cplex。IBM的Cplex有免费学术版去IBM官网注册学术账号后下载。用的版本建议选12.10或20.1太老的版本在新版Matlab上经常出现兼容性问题太新的版本比如22.1又对Matlab版本有要求。我试过Cplex 20.1配Matlab R2020b稳得不行后来换到R2022bCplex 22.1反而出现了一个诡异的报错折腾半天才定位是cplexlink没编译好最后是重装才解决。第二步安装Yalmip。Yalmip就是个纯Matlab工具箱从GitHub或者官网下载解压后把文件夹整个路径加到Matlab搜索路径里就行。官网下载的是压缩包解压后目录里会有一个yalmip文件夹用addpath(genpath(你的路径/yalmip))设置路径。环境弄好后命令行里跑一下诊断命令确认安装成功yalmiptest这个命令会检测所有已安装的求解器看到Cplex对应的状态是successful就说明配置OK了。4.2 用Yalmip建模的核心流程Yalmip建模的套路非常固定用熟了之后基本就是一套模板。核心步骤四个定义决策变量用sdpvar定义连续变量binvar定义0-1变量写目标函数把所有成本项加起来写约束条件用循环或者向量化方式添加约束调用求解器求解设置好solvesdp或optimize的参数下面这段是我项目中目标函数和约束构建的骨架代码实际运行时直接复制改参数就行% 决策变量定义 P_gt sdpvar(1, 24); % 燃气轮机出力 P_buy sdpvar(1, 24); % 购电量 P_wt sdpvar(1, 24); % 风电出力 P_p2g sdpvar(1, 24); % P2G耗电量 SOC sdpvar(1, 24); % 蓄电池SOC u_ch binvar(1, 24); % 充电状态0-1变量 % 目标函数运行成本 碳交易成本 Cost_buy sum(price_e .* P_buy); % 购电成本 Cost_gas sum(price_g .* (P_gt P_gb)); % 购气成本 Cost_om sum(om_gt .* P_gt om_gb .* P_gb); % 运维成本 Cost_curt sum(penalty .* (P_wt_forecast - P_wt)); % 弃风惩罚 Cost_co2 sum(co2_price .* E_co2_total); % 碳交易成本 Objective Cost_buy Cost_gas Cost_om Cost_curt Cost_co2; % 约束条件 Constraints []; % 电功率平衡 Constraints [Constraints, P_gt P_wt P_buy P_dis - P_p2g - P_ec - P_ch P_load]; % 爬坡约束 for t 2:24 Constraints [Constraints, -P_ramp P_gt(t) - P_gt(t-1) P_ramp]; end % 储能SOC递推 for t 2:24 Constraints [Constraints, SOC(t) SOC(t-1) eta_ch * P_ch(t) - P_dis(t)/eta_dis]; end % 充放电互斥 Constraints [Constraints, P_ch M * u_ch, P_dis M * (1 - u_ch)]; % 求解 ops sdpsettings(solver, cplex, verbose, 2, showprogress, 1); optimize(Constraints, Objective, ops);4.3 Cplex求解性能调优经验Cplex的求解速度和几个设置项关系很大。我自己实践下来最影响性能的是MIP间隙容忍度mip.tolerances.mipgap和是否开启预处理preprocessing.reduce。如果你对最优解的精度要求不是特别高可以把MIP间隙设为1%求解速度会有质的提升ops sdpsettings(solver, cplex, ... cplex.mip.tolerances.mipgap, 0.01, ... cplex.preprocessing.reduce, 1, ... verbose, 2);另外一个容易忽略的点是变量范围的设定。如果能在定义变量时直接给出上下限比如sdpvar(1,24)后立刻加约束0 P_gt P_gt_maxCplex的预处理阶段就能把变量边界收紧从而减少搜索空间。不要偷懒把所有约束都堆在后面写。还有一个建议不要过度使用implies这种逻辑约束。Yalmip的implies确实方便但它内部会生成辅助二进制变量用得太多会让模型膨胀得很厉害。能用大M法写出来的逻辑约束尽量用大M法。我自己对比过同样的模型把implies改写成大M约束后求解时间从45秒降到了12秒左右。5. 仿真结果分析与典型场景对比5.1 基础场景运行结果解读跑完模型后第一步要做的是画各设备出力曲线图。我一般用Matlab的plot画电出力和热出力的堆叠图看每个时段各种设备的出力占比从中能读出很多信息。以冬季典型日为例我跑出来的结果有几个非常明显的规律夜间23:00-6:00电价低系统倾向于多购电给蓄电池充电同时风电出力大时P2G启动把多余的电力变成天然气储存起来。白天8:00-18:00电价高燃气轮机和蓄电池放电成为主力购电量明显减少。热负荷高峰时间段余热回收和燃气锅炉两者配合燃气锅炉主要在燃气轮机热回收不足时补足。这里有个很有意思的现象加了碳交易成本后系统的购电策略会发生反转。在高排放因子的时段即使电价低系统也不一定愿意多购电因为购电对应的间接碳排放成本高。这就体现了低碳调度和经济调度的根本区别——低碳约束会重塑系统的运行策略。5.2 碳交易价格灵敏度分析只跑一组参数是不够的论文里通常要做一个碳交易价格的灵敏度分析。具体做法是把碳价从低到高设定多个档位比如50、100、150、200、250元/吨分别跑一次优化然后对比总成本、碳排放量、各设备出力结构的变化。我跑下来发现了一个规律这里以我自己的数据为例当碳价低于100元/吨时系统的调度策略几乎不受影响碳交易成本在总成本里面占比很小。当碳价上升到150元/吨以后系统开始明显减少燃气轮机出力转而增加购电和P2G投入碳排放量随之下降。当碳价超过250元/吨时碳排放量下降的边际效应就很弱了因为此时所有“便宜”的减排手段都已经用尽剩下的减排潜力需要依赖更大规模的设备投资比如增大储能配置、扩建P2G那不是调度层面能解决的。这个规律其实很符合经济学的边际递减原理。做报告或者写论文时画一张“总成本-碳价”曲线再画一张“碳排放-碳价”曲线两张图并排放在一起整个研究的核心结论就一目了然了。5.3 不同方案对比有无碳交易机制为了突出碳交易机制的效果我习惯设两个对照方案方案一基准方案只考虑经济成本不考虑碳排放方案二碳交易方案在目标函数中加入碳交易成本两个方案跑完直接对比购电结构、燃气轮机利用小时数、总碳排放量三个指标。我常用的展示方式是画柱状图对比总成本和总碳排放的双轴图能直观看出减碳的成本代价。结果通常表现为方案二相比方案一总运行成本会小幅上升比如3%-8%但总碳排放量会明显下降比如10%-20%。这就是碳交易机制在起作用——用可控的成本代价换取明显的减排效果。这个“减碳代价”的量化本身也很有研究价值。6. 常见问题排查与调试经验6.1 Yalmip常见报错速查表这三年用下来Yalmip的报错信息我几乎都见过一遍。大部分报错其实根源就那么几个整理成表格给大家参考报错信息根本原因解决办法No suitable solver for problem type模型里有Cplex不支持的约束类型检查是否有非线性项sdpvar变量之间是否有乘法Index exceeds array bounds变量维度不匹配用size()检查每个sdpvar的维度特别是循环里的索引Warning: Solver not applicable求解器未正确安装或未添加路径先跑yalmiptest确认求解器状态Infeasible problem约束过紧或模型物理矛盾先删除碳交易约束跑一遍确定基础模型可行性再逐步添加约束定位冲突NaN or Inf in model变量未初始化或除数为0检查初始条件赋值和效率参数是否有0值这里面最要命的是Infeasible problem因为Yalmip本身的报错信息很笼统你得自己去找是哪组约束导致不可行。我的经验是二分法先把所有约束屏蔽一组一组放回去放到哪组开始报不可行问题就在哪组约束里。这个方法虽然土但效率非常高。6.2 求解时间过长怎么办Cplex求解MIP问题遇到大规模模型时求解时间会飙升。按照我的经验可以从三个方向排查第一检查0-1变量的数量。每多一个二进制变量MIP分支定界的搜索空间就翻一倍。如果模型里0-1变量超过1000个就要考虑是否能把某些决策变量改成连续变量。比如充放电互斥约束在部分场景下可以用“充放电同时发生的概率极低”为由通过设置高惩罚系数替代二进制变量。第二检查约束的稀疏性。Cplex对稀疏矩阵的处理效率远高于稠密矩阵这是个我很深的体会比如我曾把循环逐时段约束改成24行×变量数的矩阵乘约束之后速度提升非常明显。Yalmip的向量化其实不太可控所以我一般在Matlab里预先算好一些系数矩阵然后用矩阵乘法生成约束这样传给Cplex的模型更稀疏。第三检查是否有多余的冗余约束。有些约束从物理上是满足的但你把它显式写进去Cplex在预处理阶段就会花费大量时间去推导。比如电功率平衡约束已经隐含了储能SOC的上下限关系如果再把SOC约束重复加一遍就属于冗余约束。6.3 四个季度典型日如何处理很多IES调度的研究要用春夏秋冬四个典型日分别跑然后汇总年运行数据。要注意的是不同季节的设备参数和负荷曲线不同不能直接复制代码跑四遍就完事有三个地方要改一是负荷曲线。冬季负荷热气电少夏季冷负荷高电负荷也可能高春秋则相对温和。用不同的负荷输入文件区分。二是风电出力曲线。不同季节的风资源差异很大需要根据实际数据设定否则结果没参考意义。三是电网购电价格。如果当地分时电价政策在不同季节有差异电价数组也要相应修改。这些基础数据没弄好的话后面求解结果再漂亮也是空中楼阁。我有个习惯是会先把负荷曲线和风电曲线画出来看一眼再跑模型如果曲线形状明显不合理比如冬季冷负荷反而比夏季高会直接停下来检查数据来源。7. 从单目标到多目标低碳调度的扩展方向7.1 碳排放量作为约束条件的处理方式前面讲的是把碳排放成本放进目标函数碳排放是“代价”。另一种做法是把碳排放量作为硬性约束总碳排放量 ≤ 碳排放上限这种做法适合有明确减碳目标的场景比如上级规定了年度碳排放总量配额。不过这种方法有个问题碳排放上限取值不同时最优解的变化是不连续的需要跑多次优化观察“碳约束的紧松程度”非常费时间。7.2 多目标优化的两种实现思路如果想同时优化“经济成本最小化”和“碳排放最小化”两个目标可以往多目标优化方向做。这里我只提两种最常规的思路第一种是线性加权法妥协规划。给两个目标各分配一个权重把双目标合成单目标min w₁ * F_cost w₂ * F_co2这个方法简单但权重怎么定是个问题而且权重变化时目标函数量纲差异成本是万元碳排放是吨会导致加权失真。第二种是ε约束法。把其中一个目标比如碳排放设为约束上限ε然后最小化经济成本。反复取不同的ε值求解就能得到帕累托前沿。这个方法在写论文时很常见画出的帕累托曲线非常直观——每个点代表一个“经济-低碳”权衡下的最优方案。我建议做研究时用ε约束法论文汇报效果好物理意义也清晰。7.3 不确定性建模的引入进阶一些的方向是考虑风电出力的不确定性。风电出力预测误差导致调度方案不鲁棒是实际运行中很大的痛点。目前主流做法是两阶段鲁棒优化——第一阶段决定机组的启停状态和储能调度计划第二阶段在不确定集内调整设备出力以保证可行性。Cplex也能解两阶段鲁棒问题但需要采用CCG列与约束生成算法迭代求解Yalmip在这个场景下通常配合gurobi或cplex的LP求解接口来写循环。这部分内容和这篇的基础模型是叠加关系在把基础模型吃透之后再做扩展会顺很多。现在我个人的体会是低碳调度的门槛其实不在求解器使用而在建模的完整性和数据的合理性。很多同学拿到代码跑通一遍就当作完成了任务但实际答辩或者审稿时最容易出问题的反而是碳排放核算方法、碳交易参数的取值依据、系统边界的定义这些“细节”。建议做这个方向的读者先把自己系统里的碳排放来源一个个搞清楚、把碳交易机制的参数来源标注清楚再往模型里加复杂度。这样得到的结论才经得住推敲代码也才真正有复用价值。本文还有配套的精品资源点击获取