综合能源系统碳势-价格双响应调度Matlab实现详解

做综合能源系统调度的同学，看到“碳势-能源价格双响应”这个题目，应该马上能意识到这不是普通的单目标优化复现。它把碳排放流里的碳势概念和需求侧价格响应机制揉在一起，在一个Matlab框架里同时解决低碳调度和经济调度的权衡问题。我复现这个项目前后花了三周时间，踩了不少坑，也把整个建模逻辑彻底摸了一遍。这篇就把项目背后的核心思路、Matlab代码实现的具体流程，以及求解过程中遇到的各种问题一次性讲清楚，给准备做同类研究的朋友一个可以直接上手的参考。

一句话说清楚这个项目是干什么的：传统综合能源系统的经济调度只盯着运行成本，低碳调度又只盯着碳排放量，两者经常矛盾——燃气轮机发电便宜但排碳高，电锅炉用电可能干净但成本高。这篇研究通过引入“碳势”作为碳排放强度的观测量，再叠加能源价格的激励，让设备出力和负荷需求在同一个优化模型里同时向“低成本”和“低碳”两个方向调整，最终生成一组兼顾经济和碳排放的调度计划。所有建模和求解都在Matlab里完成，用的是YALMIP加商用求解器这条标准技术路线。

适合谁来读？如果你正在研究综合能源系统、微电网、电力市场或者碳排放流方向，这篇文章对你会比较实用；如果你只是想把一篇EI论文的结果用Matlab复现出来，里面的建模细节和排坑经验也能直接帮你省掉大量试错时间。下面按我自己的复现思路，把整个项目拆开讲。

1. 先把这个项目的研究背景拆透

1.1 “碳势”到底是什么

在做这个项目之前，我先去翻了碳排放流理论的基础。传统碳排放核算都是“总量平均法”，整个系统一年的碳排放除以总的用电量，得出一个平均碳强度——这只能告诉用户平均水平，完全看不出不同节点、不同时段用能的碳排放差异。而碳排放流理论把碳排放当成一种可以随潮流流动的“流体”，通过潮流追踪算法，把每一台机组排放的碳精确分摊到它实际供应的负荷上。这样算出来的结果，才是真正反映“每一度电从哪里来、含多少碳”的精细信息。

这里面的关键概念就是节点碳势，表示某个节点上单位负荷所对应的等效碳排放强度，单位一般是kgCO2/kWh或者tCO2/MWh。电网的节点碳势取决于上游发电组合：风电、光伏占比高的时段，节点碳势自然低；火电、气电占比高的时段，碳势就高。天然气网的碳势相对固定，主要取决于气源的上游排放因子。

在我复现的这套综合能源系统里，电网和气网是耦合的，CHP机组同时消耗天然气、输出电和热，所以它既是气网的负荷，又是电网的电源，碳势的计算需要跨网络追踪。Matlab里面实现的时候，我是分两步做的：第一步先求解基础潮流，得到各支路功率分布和各节点注入功率；第二步根据发电机的碳排放强度数据，用比例分摊法做碳排放流追踪，得到每个节点的碳势值。这里有个容易出错的地方，就是气网和电网的碳势单位不一样，气网按单位体积天然气的隐含碳计算，电网按单位电量的碳强度计算，耦合设备做碳排放分摊时，必须把这两类单位统一，否则后面目标函数里碳势和价格加权时就乱了。我的做法是统一折算成吨标准煤当量下的碳排放因子，再参与计算。

1.2 为什么双响应机制能同时兼顾经济和低碳

还有一个容易误会的点，就是“价格响应”和“碳势响应”到底是不是重复设计。我最早读文献的时候也有这个疑问，一度觉得既然分时电价已经把用户调峰调谷了，再加一个碳势信号很多余。真正把模型建起来之后才理解，这两者服务的维度完全不同。

价格响应是经济利益驱动的。用户在峰时少用电、谷时多用电，本质是为了省电费，它能让系统运行成本下降，但不保证碳排放下降——比如谷时如果正好是火电大发，用户把负荷挪到谷时，成本是低了，碳反而排得更多了。碳势响应则不一样，它给了用户一个实时“用能含碳量”的信号，引导用户在低碳时段多用电、在高碳时段减少用电，目标是降低碳排放，但单独看它，用户没有经济利益驱动，很难自愿响应。

所以原文的思路是把两者做成一个加权需求响应系数。响应系数是归一化后的能源价格和节点碳势的函数，价格高或碳势高的时段，响应系数变大，负荷相应削减或转移；反之负荷增加。我把这个逻辑实现了之后，对比了三种场景：无响应的固定负荷、仅价格响应、价格+碳势双响应。实测下来，双响应场景下的系统总成本和碳排放量都是最低的，碳势信号确实补上了“便宜但高碳”这个盲区。这个对比结果是整个项目最核心的结论之一，也是画图时一定要保留的对比工况。

1.3 复现这个项目能锻炼哪些能力

复现这个项目跟复现普通的“单目标经济调度”不一样，它有明显的递进结构。首先你得掌握综合能源系统的基础建模，包括CHP热电联产、燃气锅炉、电锅炉、储能等设备的数学模型；其次要理解碳排放流理论，知道碳势是怎么算出来的，为什么节点碳势会随时间变化；再次要把双响应机制转化成数学模型，也就是怎么把碳势和价格对负荷的影响写成约束或目标的一部分；最后才是Matlab代码实现，包括YALMIP建模、求解器配置、结果后处理和可视化。

这四层能力不是孤立的。我在复现过程中发现，如果只把代码跑通，不理解碳势的计算原理，后面调试结果异常时完全无从下手；反过来，如果只看公式不动手写代码，对“线性化处理”“0-1变量引入”这类工程细节又很难有深刻体会。所以如果你是给自己定学习目标，建议不要只追求跑通，要把每个模块的公式和代码对应起来，这样收获会大得多。

2. 调度模型的设计思路与数学建模

2.1 系统架构与设备模型

这个项目用的是一个典型的区域级综合能源系统，电、气、热三种能源形态相互耦合。电源侧包括外部电网购电、风电/光伏、CHP机组；热源侧包括CHP余热、燃气锅炉、电锅炉；储能环节包括电储能和蓄热罐。网络侧通过电力和天然气网络连接各类设备和负荷。所有设备的模型最终都要写成可由求解器处理的数学约束。

CHP机组是最核心的耦合设备，它消耗天然气同时产电和产热。这里我用的是定热电比模型，即电输出和热输出之间始终满足一个固定比例关系，表达式是P_heat_chp = k_chp * P_elec_chp。有些论文会用可变热电比，但为了保持模型线性、方便MILP求解，定热电比是更稳妥的选择，而且对论文结论的影响不大。燃气锅炉和电锅炉相对简单，本质就是输入燃料或电能，输出热能，效率系数固定。

储能设备需要重点关注SOC的递推关系。电储能的约束包括充电功率、放电功率、SOC上下限、充放互斥、日初末SOC一致等；蓄热罐同理，需要考虑热损失系数。这里比较容易踩坑的是“充放互斥”约束，如果不加这个约束，求解器可能会让储能同时充电和放电，结果虽然满足功率平衡，但物理上完全不合理。我在代码里用两个0-1变量加一个大M做了互斥约束，实测求解速度也能接受。

风电和光伏的处理相对简单，直接作为负的负荷或者给定功率注入，使用预测值参与调度。如果要考虑不确定性，那就要引入场景法或鲁棒优化，这个项目暂时没有涉足，属于扩展方向。

2.2 目标函数：运行成本与碳交易成本的双目标权衡

调度模型的目标函数是整个项目的灵魂。我复现的版本把两个子目标统一在一个框架里：第一个是系统运行成本，包含从电网购电费用、购买天然气费用、设备运行维护费用；第二个是碳交易成本，采用的是阶梯碳交易机制。

碳交易机制的计算逻辑是这样的：系统会有一个无偿分配的碳排放配额（通常是按历史排放或行业基准确定），实际碳排放低于配额时，可以卖出多余的配额获得收益；高于配额时，超额部分需要购买，而且超额越多，购买价格越高，这就是阶梯碳交易。比如我的算例里设置的是，配额内排放免费，超配额部分在0-1000kg区间时碳价是0.2元/kg，在1000-2000kg区间时跳到0.3元/kg，再往上继续跳。

这样处理之后，目标函数就合成一个单目标：总成本 = 运行成本 + 碳交易成本。实际上模型里只维护了一个目标函数，通过碳交易价格这个杠杆，把“低碳”量化成了“省钱”，这就解决了多目标怎么加权的问题。这里并没有像有些论文那样用权重系数去平衡成本和碳排放，而是直接用碳交易市场机制来协调，我觉得这是这个方案在工程逻辑上最合理的地方。

用Matlab表示目标函数的大致结构是：

% 运行成本：购电 + 购气 + 运维 C_operation = sum(price_e .* P_buy) * delta_t ... + sum(price_g .* G_buy) * delta_t ... + sum(c_om .* P_device) * delta_t; % 碳交易成本：阶梯碳价 C_carbon = carbon_cost(E_total, quota, price_level, price_step); % 总目标 objective = C_operation + C_carbon;

这里carbon_cost我单独写了一个辅助函数，输入总排放、配额、碳价档位和区间宽度，输出碳交易成本。因为它是分段的，所以在YALMIP建模时要用辅助变量加线性约束来实现，后面细说。

2.3 约束条件：从功率平衡到设备运行边界

约束条件方面，首先最基本的当然是电功率平衡和热功率平衡，每个时刻的发电/购电/放电之和等于负荷/充电/电锅炉用电之和，热力侧同理。然后是各设备的出力上下限约束、爬坡约束、储能SOC边界、天然气网络节点压力或流量的简化约束，以及需求侧双响应后的负荷转移约束。

这里要特别说一下双响应负荷约束的实现。基本思路是先计算每个时刻的需求响应系数，该系数在[0,1]区间内，然后让实际参与调度的负荷等于基础负荷乘以（1 - 响应系数*削减比例），同时满足一天总负荷转移量不超过一定比例。我在代码中把这个过程拆成了两步：第一步用上一天的碳势和价格数据计算响应系数，第二步把调整后的负荷作为已知参数代入调度模型。这样做的好处是模型保持线性，求解速度快，坏处是碳势和负荷调整之间没有形成闭环反馈。如果想让模型更精细，可以在迭代中反复更新碳势和负荷直到收敛，这也是后续改进的一个方向。

功率平衡约束的写法：

% 电功率平衡：购电 + CHP电出 + 风电 + 光伏 + 储放电 = 负荷 + 电锅炉 + 储充电 for t = 1:T constraints = [constraints, ... P_buy(t) + P_chp_e(t) + P_wt(t) + P_pv(t) + P_dis(t) ... == P_load(t) + P_eb(t) + P_chg(t)]; end

3. Matlab代码实现全流程

3.1 环境配置与工具箱选型

先说一下我的运行环境：Matlab R2023b，优化建模用的是YALMIP工具箱，求解器用的是IBM CPLEX 12.10。这套组合是学术界做MILP优化的标配，YALMIP负责把数学模型翻译成求解器能识别的标准形式，CPLEX负责真正求解。

如果你是第一次配置环境，有几点需要注意。第一，YALMIP需要下载后手动添加到Matlab路径，用addpath命令加载到工作区即可，不需要额外安装步骤。第二，CPLEX需要用安装包里的MATLAB接口文件，把对应的bin和matlab目录都加到路径里。第三，版本兼容性一定要先验证好，我一开始用的是Matlab R2021b搭配CPLEX 12.8，跑起来没问题；后来换到R2023b，CPLEX的接口没更新，就一直报Unable to load cplexmex，折腾了半天才发现是接口文件版本不匹配。所以在配置环境这一步，建议先跑一个最简单的测试：

x = sdpvar(1,1); optimize(x >= 1, x);

能正常返回结果，再继续往下写。

关于Gurobi，它和CPLEX在这个项目里是可以互换的。我个人觉得Gurobi对MILP的默认参数调优更好，求解速度可能会快10%-20%，但CPLEX在学术界的安装授权更省心，所以如果你的机器空间足够，两个都装上最好，写代码的时候用sdpsettings('solver','cplex')或sdpsettings('solver','gurobi')一行就能切换。

3.2 基础数据与算例构建

数据和模型代码分开是最重要的工程习惯。我的做法是建立一个initData.m脚本，里面把所有算例参数集中定义，包括24小时的负荷曲线、风电光伏预测出力、分时电价、天然气价格、设备容量、效率系数、碳交易参数等。调度主程序只需要调用这个脚本即可。

这里分享一个数据构建的小技巧：24小时的负荷数据和预测出力，不要手写，用正弦函数加噪声生成比较方便。比如电负荷的基准曲线可以设计为早晚两个高峰的形态，热负荷设计为冬季的夜间偏高形态。风电预测出力用夜大昼小的曲线，光伏用正午大、早晚为0的曲线。这样生成的算例逻辑合理，画图也好看。

另外，所有物理量纲一定要统一。我在初版建模时，因为把kW和MW混用，导致求解结果数量级乱七八糟，功率平衡约束“不满足”但求解器还返回了“最优解”。后来我统一把功率、热功率都换算成MW，把碳排放量统一用吨表示，能量用MWh，这样再也没出过类似问题。

3.3 从YALMIP建模到求解输出

建模这块，我的代码结构大致是这样的：先用sdpvar定义决策变量，包括CHP电出力和热出力、燃气锅炉出力、电锅炉功率、储能充放电功率、蓄热罐蓄放热功率、购电功率等；然后用optimize函数传入目标函数和约束集合；最后用value()函数提取结果。

写约束时我比较推荐一个习惯：每个约束都用cell数组累积起来，避免最后传参数时漏掉某一条。核心代码框架：

T = 24; constraints = []; % 决策变量 P_chp_e = sdpvar(1, T); P_chp_h = sdpvar(1, T); P_gb = sdpvar(1, T); P_eb = sdpvar(1, T); P_chg = sdpvar(1, T); P_dis = sdpvar(1, T); SOC = sdpvar(1, T+1); u_chg = binvar(1, T); % 充电状态标志 u_dis = binvar(1, T); % 放电状态标志 % 约束：CHP热电比 constraints = [constraints, P_chp_h == k_chp * P_chp_e]; % 约束：储能SOC递推与互斥 for t = 1:T constraints = [constraints, SOC(t+1) == SOC(t) + P_chg(t)*eta_chg - P_dis(t)/eta_dis]; constraints = [constraints, P_chg(t) <= M * u_chg(t)]; constraints = [constraints, P_dis(t) <= M * u_dis(t)]; constraints = [constraints, u_chg(t) + u_dis(t) <= 1]; end % 目标函数 objective = sum(price_e .* P_buy) + sum(price_g .* G_buy) + ...; % 求解 ops = sdpsettings('solver', 'cplex', 'verbose', 2); optimize(constraints, objective, ops);

算例配置为24小时调度，决策变量上百个，CPLEX通常几秒到十几秒就能求解完毕，作为EI论文复现的算例规模是完全够用的。求解完成后，我一般会保存一个result.mat文件，包含各设备的出力序列、碳势序列、负荷调整前后曲线、成本和碳排放总量。然后单独写一个plotResults.m脚本做画图，这样每次画图不需要重新求解模型，节约大量时间。

4. 代码中的关键细节与调参经验

4.1 碳势迭代计算的一个优化技巧

前面提到，碳势可以直接由潮流结果算出来，逻辑上并不复杂。但实际运行时会面临一个问题：需求侧双响应改变了负荷分布，负荷分布又反过来影响潮流和碳势，碳势变化后又会影响下一轮需求响应系数。这是一个典型的耦合迭代关系。

我采用的务实方案是迭代计算：先用固定负荷算一次碳势，更新负荷，再算一次碳势，重复这个过程。实测下来，迭代3-5次后碳势就基本稳定了，扰动幅度能在1%以内。这个方案的优点是计算量可控，代码实现简单；缺点是它不是严格的全局最优解，跟原论文可能有点差别。我在复现报告里专门说明了这个处理方式，结论部分的对比图表质量没有受影响。

另外，碳势的数值波动很容易出现异常，比如某几个节点突然出现负碳势。这个通常是因为追踪算法里出现了循环扣除或者数值精度问题。我的排查方法是打印每个节点的碳势分布，如果出现负数，检查对应的潮流方向是不是和追踪方向一致。这个问题在直流潮流模型里很少见，但如果用交流潮流就会出现得比较多，所以新手建议先用直流潮流模型。

4.2 YALMIP建模中的几个常见坑

第一个坑是关于变量类型。YALMIP的binvar是二值变量（0-1），intvar是整数变量，sdpvar是连续变量。对于储能充放互斥这类约束，用binvar；对于阶梯碳交易的分档判断，也可以用binvar。很多新手会把binvar误写成sdpvar再加范围0到1的约束，这样求解出来的解是连续的，虽然很多时候也能得到一个结果，但严格来说已经不是MILP了，充放电互斥效果可能不彻底。

第二个坑是大M取值问题。互斥约束里的大M不能取太大，太大容易导致数值稳定性问题，最佳实践是取该变量物理上限的1.1倍，比如充放电功率上限是20MW，M就取22。取值太大会让求解器在预求解阶段就把约束“软化”，出现一些微妙的不满足；太小又会让可行域被错误压缩，所以我一般会按物理量级来定M，而不是随便用一个1e6。

第三个坑是用function句柄封装非线性函数后直接传给YALMIP，YALMIP对这类函数的求导和梯度处理比较有限，所以尽量做线性化。本项目的碳交易阶梯成本和双响应系数都是分段的，我都按线性约束处理了。凡是能用线性约束表达的，绝不引入非线性，这是用YALMIP做调度的铁律。

4.3 结果可视化：怎么把调度结果讲清楚

复现EI论文，图表质量直接决定整套工作的可信度。我最后交付的图有六张：一张系统结构拓扑图、一张24小时电热负荷调整前后对比、一张各设备出力堆叠图、一张节点碳势时变曲线、一张碳排放量和运行成本的场景对比柱状图、一张碳交易价格灵敏度分析曲线。

画图时我比较注重颜色的区分度，采用红绿蓝橙紫六种色系，线宽统一在1.5，字体用Times New Roman 10.5磅。堆叠图用area函数，碳势曲线用plot函数。Matlab默认的配色在多曲线场景下很难区分，建议用paruly色带或者自定义RGB颜色。比如设备出力堆叠图，从上到下按“购电-CLP-风电-光伏-储能放电”的顺序堆叠，图例清晰，颜色由深到浅，读者一眼就能看出每种能源在全天不同时段的贡献。

5. 常见问题与排查经历

5.1 求解器报“Infeasible Problem”

这是所有做优化调度的人都会遇到的头号问题。我的排查顺序是这样的：第一步检查约束是否存在互相矛盾，比如SOC末状态要求等于初状态，但充放电效率太低导致物理上达不到，这种就要调整参数或者放宽末状态约束；第二步检查大M是否太小，互斥约束和阶梯分段约束都可能因此造成不可行；第三步把目标函数临时改为常数0，只求一个可行解，看问题出在哪个约束。如果可行解都找不到，说明约束本身就有问题，可以逐条注释约束来定位。

5.2 碳势数值异常或总碳排放对不上

排查步骤：先检查潮流结果是否正确，再检查发电机碳排放强度是否被正确赋值，最后检查气网和电网的碳势单位是否统一。我调试过程中出现过一次碳势整体偏高，原因是我把天然气的碳排放因子单位写错了，导致CHP购气对应的碳排放被重复计入两次。这个错误直到我把“系统总碳排放=电网侧碳排放+气网侧碳排放-耦合设备重复部分”这个等式列出来作校验时才发现。

5.3 求解时间过长的调参方案

如果模型求解时间超过5分钟，一般有两类原因：一是整数变量过多，MILP的搜索空间爆炸。这时要检查有没有可以用连续变量替代的整数变量，或者把不必要的时间耦合约束放松。二是求解器参数没调好，我通常会把CPLEX的MIP gap容忍度设为0.5%而不是默认的0%，在工程和学术精度要求下，这个设置能明显加速。还可以设置求解时间上限，比如180秒，到时取当前最优解，这样代码即使遇到复杂算例也不会卡死。

5.4 复现过程中最容易被忽视的数据核对表

下面这个表格是我复现时整理的核对表，每次跑完算例我都会过一遍：

6. 复现路线与扩展思路

6.1 从论文到代码的落地顺序

这里给一个我实际走的复现路线，按这个顺序来应该能少走弯路。第一步，先把论文里的系统拓扑画清楚，确认有哪些设备、哪些能源形式、哪些耦合关系；第二步，把设备模型和约束列成数学公式，统一变量符号和单位；第三步，先在Matlab里实现一个最简单的“固定负荷+经济调度”，保证求解器能正常求解；第四步，加入碳势计算模块，固定负荷情况下输出碳势结果并和论文对照；第五步，加入价格响应和碳势响应，实现双响应负荷模型；第六步，加入碳交易成本和整体目标函数，形成完整版本；最后才是做场景对比、灵敏度分析、画图和写复现报告。

每走一步都要及时验证中间结果，而不是等所有代码写完一次性跑。这个习惯帮我节省了大量时间，因为错误在一开始就能被发现，越到后期排查成本越高。这个项目本身对Matlab的依赖非常重，建议你先把YALMIP和求解器环境一次性配好，后面所有调试都在这套环境里进行，避免频繁切换。

6.2 这个模型可以怎么扩展

复现之后如果想做自己的研究或论文，可以从几个方向扩展。第一是考虑不确定性，给风电光伏出力加场景法或者鲁棒优化，这也是目前IES调度的主流方向；第二是增加时间尺度联动，从日前调度扩展到日内滚动，甚至加实时反馈环节；第三是引入更多碳市场机制，比如碳捕集、绿证交易、综合需求响应；第四是把单区域扩展为多区域互联，研究碳势在不同区域之间的转移和碳泄漏问题。在代码层面，这套YALMIP建模框架的可扩展性很好，新增一个设备就是新增几组约束和变量，结构上不会伤筋动骨。

另外还有一个很实用的工程扩展：把Matlab的调度结果输出为Excel表格，再结合Python的绘图库做更漂亮的图表，或者写一个简单的GUI让参数可调，方便做灵敏度分析。这些扩展在工程上都不难，关键是核心优化模型要稳定可靠。

我复现完这个项目最大的感受是，一套好的调度模型，设计逻辑远比代码技巧重要。碳势-价格双响应这个机制，乍看只是把两个因素组合起来，但真正实现的时候你会发现，它逼着你把“碳排放是怎么产生的”“价格是怎么影响行为的”“两者如何在数学上统一”这几个问题想得特别透彻。Matlab和求解器只是工具，真正值钱的是对物理过程的理解和模型结构的把握。

最后分享一个小建议：做这类复现项目，不要只追求跑通出图，一定要自己动手把每个公式和每段代码对应起来，把关键参数改一改，看看结果的变化是否符合物理直觉。比如把碳交易价格调高一倍，再看系统是不是会更倾向于用电锅炉替代燃气锅炉；把分时电价的峰谷差拉大，再看负荷转移幅度是不是明显增加。只有做到这一步，你才真正把这个项目消化成了自己的东西。希望这篇分享能让你在复现时少踩几个坑，有疑问欢迎交流。