先说结论这类题目复现起来其实不难真正劝退人的不是算法本身而是建模细节。很多人拿到“考虑综合需求响应和碳交易机制的综合能源系统优化调度策略”这个标题第一反应是打开 Matlab 直接写求解代码结果半天过去连决策变量都没理清楚。我复现过几个类似的综合能源系统优化调度模型之后最大的体会是70% 的工作量在数学建模20% 在数据整理真正的求解代码反而只占 10%。这篇内容就按我实际动手的顺序把这个题目的建模思路、Matlab 实现骨架、以及那些论文里不会写但特别容易踩的坑完整拆一遍。这个题目适合正在做综合能源系统方向课题的研究生、需要复现论文算法的工程师以及想用 MatlabYalmip 快速搭建 MILP 调度模型的同学参考。下面直接进入正题。1. 这题真正要复现的是什么先看懂系统的能量流与耦合关系1.1 综合能源系统的典型结构与设备集合综合能源系统Integrated Energy System, IES之所以叫“综合”核心在于它打通了电、热、气三种能源载体。典型的园区级系统包含这些设备热电联产机组CHP最常见的是燃气轮机或内燃机同时产电和产热是电-热耦合的核心设备燃气锅炉GB补足热负荷缺口模型比较简单输入天然气、输出热量风力发电和光伏WT/PV作为清洁能源优先消纳在调度模型中通常处理为“负负荷”或带上限的出力变量电储能ES电池储能负责电负荷的削峰填谷蓄热罐HS热储能给热负荷增加调节柔性上级电网和天然气网能源买入口购电有分时电价购气有阶梯气价系统的能量流动可以概括为上级电网和天然气网作为外部能源来源经过 CHP、GB 等转换设备配合储能和新能源共同满足电、热、气三类负荷需求。很多复现论文还会加电转气P2G设备但那是另一条技术路线会显著增加模型复杂性如果是第一次复现建议先不加。建模的第一步是把这个系统图画成一张“能流表”每一列是一个设备每一行是电/热/气三种能源交叉点标注该设备对这种能源的“产出 1”或“消耗 -1”。这张表在后面写约束的时候非常管用因为所有平衡约束本质上都是“某时段某能源的所有流入等于所有流出”。1.2 综合需求响应不只是“削峰填谷”传统电力需求响应基本只针对电负荷比如峰谷电价引导用户错峰用电。综合需求响应Integrated Demand Response, IDR的关键区别在于“综合”二字电、热、气三类负荷同时具备可调柔性而且不同能源之间的替代关系也被利用起来。举个最直观的例子用户用燃气壁挂炉取暖在电力现货价格极高的时段可以削减一部分电负荷、同时增加一部分气负荷来补足相同的用能需求反过来在电价低谷时段用电采暖替代燃气供热。这种跨能源载体的负荷转移是综合需求响应区别于普通需求响应的核心价值。在建模层面我通常把 IDR 拆成三类动作可削减负荷在某一时段直接少用一部分能源代价是按固定单价支付补偿可转移负荷把某一时段的用能需求整体平移到其他时段需要保证一天内转移总量守恒可调节负荷负荷本身在一个区间内连续可调比如电采暖功率允许在额定值的 40%~100% 之间波动这三类动作分别对应不同的优化变量和约束后文会具体展开。1.3 碳交易机制在这类调度模型中的角色碳交易机制从字面上看像个“环保加分项”但在数学建模里它直接进目标函数。这意味着它不只是让结果更好看而是通过改变不同出力的经济成本从机制上促进系统偏向低碳运行。一般学术论文里对碳交易机制的建模是这套逻辑系统根据运行情况计算实际碳排放总量这个排放来自两部分一是从上级电网购电对应的间接排放用电网平均排放因子折算二是消耗天然气对应的直接排放按天然气排放因子折算。系统持有一定数量的免费碳排放配额如果实际排放低于配额可以把富余配额拿到市场上出售获得收益如果实际排放超过配额超额部分需要购买配额而且很多论文会设置阶梯式碳价——超额越多单价越高形成累进式惩罚。把这个机制放入优化模型后系统在做调度决策时会自动权衡多买电可能便宜但电网电的排放因子高碳交易成本随之上升多用气可能贵一些但天然气的单位排放相对更可控碳成本低。于是最优解会在“购能成本”和“碳交易成本”之间寻找平衡点。这也是为什么很多论文里只要加入碳交易机制CHP 和燃气锅炉的出力占比就会明显变化。2. 数学建模从论文公式到可直接编码的约束体系2.1 目标函数的三块成本购能、运维、碳交易优化调度的目标函数几乎都是总运行成本最小化常见写法是min C_total C_grid C_gas C_om C_co2 C_drC_grid向上级电网购电的成本等于各时段购电量乘以分时电价再累加C_gas购气成本单位是元/kWh 热价或元/m³取决于天然气的热值折算方式C_om设备运行维护成本通常是设备出力的线性函数比如 CHP 每发一度电收 x 元维护费C_co2碳交易成本或收益如果是负值说明卖出配额获得收益C_dr需求响应的补偿成本削负荷要付补偿转移负荷可能有激励成本有一点必须提醒目标函数各项的单位要完全统一。常见错误是把功率单位 MW 直接和能量单位 MWh 相加或者购气成本用了元/m³而设备模型里的热值折算用的是 MJ导致最终数量级差出好几个量级。后面第 4 章会专门讲单位问题。2.2 设备出力模型与耦合关系设备建模是整个优化模型的地基。这里给出最常用的线性模型CHP 热电联产机组P_chp(t) eta_chp * V_gas_chp(t) * LHV_gas % 电出力 H_chp(t) V_gas_chp(t) * LHV_gas - P_chp(t) % 热出力由能量平衡推出很多论文为了简化直接使用热电比耦合约束H_chp(t) K_chp * P_chp(t)K_chp 是热电比如果原论文没有特别说明建议用可调区间来建模而不是固定值H_chp_min H_chp(t) - K_chp * P_chp(t) H_chp_max这样既保留了 CHP“以热定电”或“以电定热”的灵活性又不会因为热电比参数不准确而把可行域限制死。燃气锅炉H_gb(t) eta_gb * V_gas_gb(t) * LHV_gas燃气锅炉只产热没有热电耦合建模最简单。电储能soc(t1) soc(t) (P_ch(t) * eta_ch - P_dis(t) / eta_dis) * delta_t soc_min soc(t) soc_max P_ch(t) P_ch_max * u_ch(t) P_dis(t) P_dis_max * u_dis(t) u_ch(t) u_dis(t) 1这里需要特别注意soc 是能量状态而功率乘以时段长度 delta_t 才是能量变化量。如果 delta_t 是 1 小时数值上功率和能量相等如果后面改成 15 分钟调度delta_t0.25忘记乘系数会直接导致储能状态计算错误。另外二元变量 u_ch 和 u_dis 的互斥约束是防止储能同时充放电的必要条件虽然惩罚项也能缓解但最可靠的还是显式约束。蓄热罐建模方式与电储能完全同构只是没有充放电互斥的约束那么严格因为热系统惯性更大允许一定程度的同时蓄放。2.3 碳交易的配额分配与阶梯碳价碳交易建模分两步先算排放再算成本。实际碳排放总量的计算E_total sum( P_grid(t) * EF_elec * delta_t ) sum( V_gas(t) * LHV_gas * EF_gas )其中 EF_elec 是电网购电的排放因子EF_gas 是天然气燃料的排放因子。注意这里的单位换算非常讨厌EF_elec 通常是 kgCO2/kWh而 EF_gas 可能是 kgCO2/MJ 或者 kgCO2/m³需要先统一到同一个能量单位体系里。免费配额常见做法是给定一个基准排放量乘以分配比例或者按设备额定容量折算E_quota quota_ratio * E_base超额排放量E_excess E_total - E_quota当 E_excess 为正成本为正当 E_excess 为负即排放低于配额未使用的配额可以出售碳交易成本为负。不过很多论文为了简化只考虑超额购买的情况复现的时候要看清原模型设计。阶梯碳价的线性化是这一节的重点放到第 4 章专门讲这里先只说概念阶梯碳价就是超额排放量越大边际碳价越高形成分段的线性成本函数。在 MILP 模型里这个函数需要用分段线性化技术处理不能直接在 Yalmip 里写 if-else。2.4 综合需求响应的价格弹性建模综合需求响应有两种建模路线。路线一价格弹性矩阵法。定义电负荷转移率与电价变化量的关系delta_P_e(t) / P_e_base(t) sum_j epsilon(t,j) * delta_price(j) / price_base(j)epsilon(t,j) 是自弹性系数tj或交叉弹性系数t≠j。这种方法理论上很漂亮但实际复现时有个大麻烦弹性矩阵的参数很难从论文或实测数据里得到很多时候论文作者自己都不清楚这些数是怎么标定出来的。我只建议在复现那些明确给出弹性矩阵数据的论文时采用这种方法。路线二比例边界法推荐。直接把需求响应能力表达为基线负荷的百分比并给转移负荷加守恒约束P_load(t) P_fix(t) P_trans(t) - P_cut(t) sum(P_trans) 0 % 转移负荷一天内总量守恒 -alpha * P_fix(t) P_trans(t) alpha * P_fix(t) 0 P_cut(t) beta * P_fix(t)相比弹性矩阵法这种方法参数更少、含义更直观——alpha 表示可转移比例beta 表示可削减比例——而且线性约束好求解。热负荷和气负荷的建模类似只是转移方向主要反映用能替代比如某时段削减了电采暖负荷同时增加了燃气供热需求。需求响应补偿成本根据削减量和削减单价计算转移部分如果没有明确的补偿单价约定可以简化为不进入目标函数只作为约束存在。2.5 可再生能源与平衡约束风电和光伏在调度模型里通常处理为“优先消纳”的出力变量0 P_wt(t) P_wt_pred(t) 0 P_pv(t) P_pv_pred(t)如果原论文考虑“弃风弃光”还可以在目标函数里加一项弃电惩罚但大多数标准模型默认全额消纳直接按预测值作为上限即可。系统平衡约束是最后收口的三条等式电平衡: P_grid(t) P_chp(t) P_wt(t) P_pv(t) P_dis(t) P_load(t) P_ch(t) 热平衡: H_chp(t) H_gb(t) H_hs_dis(t) H_load(t) H_hs_ch(t) 气平衡: V_gas(t) V_gas_grid(t) V_load(t)气网这边的平衡需要小心有的模型把气负荷分成“燃气锅炉用气”“CHP 用气”“居民气负荷”三类购气总量是三者之和有的模型把 CHP 和 GB 的直接购气已经折算进设备模型。复现的时候一旦发现约束数量对不上首先检查气平衡是否重复计算了。3. MatlabYalmip 实现变量定义、求解器配置与代码骨架3.1 为什么选 Yalmip 而不是手写优化器这类模型本质是一个混合整数线性规划MILP决策变量里有储能充放电状态这样的 0-1 变量。如果你打算用 Matlab 自带的linprog或者fmincon硬解配置约束时非常痛苦而且求解效率受限于算法选择。我的经验是直接用 Yalmip 做建模层。它是一个免费的 Matlab 工具箱核心作用是把优化变量、约束、目标函数用 Matlab 语法写出来然后后端对接商用或开源求解器。好处很明显不需要手写约束矩阵逻辑大幅简化求解器可切换Gurobi 不行就换 Cplex或者切回 Matlab 自带的intlinprog约束的增删改非常灵活做情景对比时尤其方便安装好 Yalmip 之后确认求解器能被识别即可yalmiptest3.2 决策变量的组织方式我用 24 时段调度为例展示一个完整的变量定义骨架。核心原则是同一类变量定义成列向量而不是逐时段写标量变量这样约束和目标函数都可以向量化代码干净且求解效率高。T 24; % 调度时段数典型是 24 小时 % 设备出力变量 P_chp sdpvar(T, 1); % CHP 电出力 H_chp sdpvar(T, 1); % CHP 热出力 H_gb sdpvar(T, 1); % 燃气锅炉热出力 P_wt sdpvar(T, 1); % 风电出力 P_pv sdpvar(T, 1); % 光伏出力 P_grid sdpvar(T, 1); % 上级电网购电量 % 电储能 P_ch sdpvar(T, 1); % 储能充电功率 P_dis sdpvar(T, 1); % 储能放电功率 u_ch binvar(T, 1); % 充电状态 u_dis binvar(T, 1); % 放电状态 soc sdpvar(T1, 1); % SOCT1 是为了包含初始时刻和结束时刻 % 需求响应变量 P_trans sdpvar(T, 1); % 可转移电负荷正为转入负荷 P_cut sdpvar(T, 1); % 可削减电负荷 % 碳交易相关 E_total sdpvar(1, 1); % 总碳排放量 C_co2 sdpvar(1, 1); % 碳交易成本binvar对应 0-1 变量sdpvar对应连续变量。储能变量的 T1 维度处理是很多新手容易忽略的——SOC 的状态转移约束会涉及 t1 和 tT1 两个边界时刻。3.3 目标函数与约束的代码化目标函数建议按成本项分块计算最后求和。这样后续调试时可以单独查看每一项的值方便定位异常。% 分时电价、气价、运维系数、碳排放参数示例数据 price_e [0.35; 0.5; ...; 0.3]; % 24x1元/kWh price_g 0.32; % 元/kWh天然气热值折算后 delta_t 1; % 时段长度 1h % 目标函数分项 C_grid price_e * P_grid * delta_t; % 购电成本 C_gas price_g * V_gas_total * delta_t; % 购气成本 C_om k_chp * (P_chp H_chp) k_gb * H_gb; % 运维成本 C_dr dr_price * P_cut * delta_t; % 需求响应补偿 objective C_grid C_gas C_om C_co2 C_dr;约束部分我用一个循环来写储能的状态转移约束这是最直观的写法Constraints []; % 储能SOC递推约束 for t 1:T Constraints [Constraints, soc(t1) soc(t) ... (P_ch(t) * eta_ch - P_dis(t) / eta_dis) * delta_t]; end % 初始和最终SOC约束 Constraints [Constraints, soc(1) soc_0]; Constraints [Constraints, soc(T1) soc_0]; % 周期性运行要求 % 储能充放电互斥 Constraints [Constraints, P_ch P_ch_max .* u_ch]; Constraints [Constraints, P_dis P_dis_max .* u_dis]; Constraints [Constraints, u_ch u_dis 1];注意等式约束用的是Yalmip 对浮点误差很敏感。如果你的参数导致矩阵病态经常会出现“不可行”或者求解结果明显违背物理规律优先检查是不是系数差了好几个数量级。其他约束比如功率平衡、设备上下限、需求响应变量边界写法上完全一致。把所有的Constraints累积起来最后调用求解器opts sdpsettings(solver, gurobi, verbose, 2); ops optimize(Constraints, objective, opts);求解完成之后取值用value()P_chp_opt value(P_chp); P_grid_opt value(P_grid); soc_opt value(soc);3.4 求解器选型与参数设置求解器选择直接影响能不能解出来、解多快。我按优先级排序求解器类型说明Gurobi商用MILP求解器有学术license首选求解快且数值稳定Cplex商用MILP求解器有学术license与Gurobi同级别项目要求时使用intlinprogMatlab自带MILP求解器免安装小规模模型够用规模大时偏慢SCIP开源MILP求解器实在没有商用license时的后备方案对于 24 时段的园区级 IES 调度模型决策变量数通常不到 500 个含 0-1 变量这属于小规模 MILPintlinprog也能在几秒到几十秒内求解。但如果你想做 96 时段15 分钟或者含多场景随机优化的扩展建议还是装 Gurobi。sdpsettings里比较实用的几个参数也提一下opts sdpsettings(solver, gurobi, ... verbose, 1, ... % 打印迭代日志 gurobi.MIPGap, 0.0001, ... % MIP相对间隙越小越精确 gurobi.TimeLimit, 300); % 300秒限制防止卡死4. 复现路上最容易翻车的五个细节4.1 单位与数量级错位这是我在帮别人排查代码时遇到最多的问题。论文里经常混用这几类单位功率kW/MW、能量kWh/MWh、热值MJ/m³、kWh/m³、排放因子kgCO2/kWh、kgCO2/MJ。一个典型事故现场购气价格是 2.8 元/m³而设备出力模型用的是 kW 和 kWh天然气的低位热值按 9.97 kWh/m³ 算。如果你直接把购气量乘气价但设备那边把 m³ 当成了 kWh就会把燃气能量算大或算小。调度结果会出现“CHP 疯狂发电、电储能完全不充电”这种明显不合理的现象。我的建议是在写代码之前把所有参数统一折算到一个单位体系里。我习惯统一使用 MW 和 MWh天然气能量统一用 MWh碳排放统一用 tCO2吨。在代码里专门建一个参数文件把原始参数和统一后的参数对照写清楚LHV_gas 9.97; % kWh/m3 LHV_gas_MWh LHV_gas / 1000; % MWh/m3 gas_price_per_kwh 2.8 / LHV_gas; % 元/kWh 热价 gas_price_per_mwh gas_price_per_kwh * 1000; % 元/MWh 热价最后检查目标函数的量纲所有成本项加起来应该单位一致都是元这样目标函数才有意义。4.2 储能 SOC 的初值/终值约束如果没有给 SOC 加初始值和终值约束会出现一个很搞笑的优化结果储能系统在最后一个时段把电全部放光SOC 直接降到下限因为“偷掉”这个能量可以减少购能成本。这在单日调度里是可行解但不是实际可运行的策略。解决办法是在模型里要求系统按日循环运行Constraints [Constraints, soc(1) 0.2]; % 初始SOC Constraints [Constraints, soc(T1) soc(1)]; % 周期回归如果论文中没有明确给出初始 SOC我一般取 SOC 上下限的中值然后强制终值等于初值。这样储能的作用是削峰填谷而不是额外提供“免费能量”。4.3 阶梯碳价的分段线性化很多论文里的碳交易机制是阶梯累进碳价超额排放量在第一个区间内的碳价是 c1超过一定量后进入第二个区间碳价跳到 c2以此类推。这种分段函数在 Yalmip 里不能用 if-else需要显式引入 0-1 变量来做区间指示。假设超额排放量 x 有三个区间[0,a]、[a,b]、[b,inf)对应累进碳价 c1c2c3可以这样线性化a 100; b 200; % 区间边界 c1 60; c2 90; c3 120; % 单位碳价元/吨 M 1e4; % 大M常数 u1 binvar(1,1); u2 binvar(1,1); u3 binvar(1,1); x1 sdpvar(1,1); x2 sdpvar(1,1); x3 sdpvar(1,1); Constraints [Constraints, x x1 x2 x3]; Constraints [Constraints, 0 x1 a * u1]; Constraints [Constraints, 0 x2 (b - a) * u2]; Constraints [Constraints, 0 x3 M * u3]; Constraints [Constraints, u1 u2 u3 1]; C_co2 c1 * x1 c2 * x2 c3 * x3;这样碳交易成本就变成了一个线性表达式可以安全放进 MILP。大M的值不要取得过大取略大于该变量物理上限的数量级即可太大会导致数值求解困难。4.4 需求响应变量的可行域综合需求响应建模时转移负荷最容易出问题。忘记写sum(P_trans) 0这个守恒约束优化器会直接在低电价时段疯狂“转入”负荷、在高电价时段“转出”负荷人为拼出无穷的能量来源最终得到一个失真且毫无意义的最优解。另外可削减负荷的上下限也要和基线负荷挂钩。如果没有显式约束优化器会把削减量开到无穷大来降低成本。正确写法是Constraints [Constraints, sum(P_trans) 0]; Constraints [Constraints, - alpha .* P_fix P_trans alpha .* P_fix]; Constraints [Constraints, 0 P_cut beta .* P_fix];这里的 alpha 和 beta 是论文给的负荷可调比例。如果论文没给常见范围是 0.1 到 0.2我一般先按 0.1 跑通模型再扫参看结果趋势。4.5 情景对比实验的设计复现论文还有一个隐藏任务复现论文里的实验结果图尤其是那几个对比情景的图表。如果只看最终最优解不做情景对比很难说明需求响应和碳交易机制各自的作用。标准的消融实验设计是四个情景情景是否考虑综合需求响应是否考虑碳交易机制情景1否否情景2是否情景3否是情景4是是每个情景除了这两个开关外其他所有参数设备参数、负荷数据、分时电价、气价、风电光伏预测必须完全一致。我在代码里会用一个标志位控制模型生成逻辑flag_idr 1; % 是否考虑综合需求响应 flag_co2 1; % 是否考虑碳交易机制 if ~flag_idr % 不添加负荷转移和削减的可调变量直接把负荷当成固定值 Constraints [Constraints, P_load P_fix]; else Constraints [Constraints, P_load P_fix P_trans - P_cut]; % ... 其他需求响应约束 end if flag_co2 % 目标函数里加入 C_co2 objective objective C_co2; end这样跑四次收集四组结果后面画对比图就非常顺手。5. 结果可视化与灵敏度分析该怎么写5.1 设备出力与功率平衡图复现类工作最直观的成果就是调度计划曲线很多人一到画图环节就把解锁代码写成一堆散点看起来乱。我通常用stairs画出阶梯状曲线因为 24 时段离散调度的本质就是阶梯式决策stairs比plot更真实。核心设备出力画在同一张图上figure; stairs(1:T, P_grid_opt, LineWidth, 1.5); hold on; stairs(1:T, P_chp_opt, LineWidth, 1.5); stairs(1:T, P_wt_opt P_pv_opt, LineWidth, 1.5); stairs(1:T, P_dis_opt - P_ch_opt, LineWidth, 1.5); legend(购电,CHP,新能源,储能净放电,Location,best); xlabel(时段/h); ylabel(功率/MW); grid on;再画一张功率平衡验证图把电负荷曲线和所有电源出力叠加确认每个时段加总都能对上。这一步虽然是验证步骤但很多人忽略导致论文里出现“平衡约束没满足”的低级错误。5.2 需求响应前后的负荷曲线对比需求响应的效果一般用原始负荷曲线和优化后的实际负荷曲线对比来呈现可以直观看到削峰填谷的效果figure; stairs(1:T, P_fix, --, LineWidth, 1.5); hold on; stairs(1:T, value(P_load), LineWidth, 1.5); legend(原始负荷,需求响应后负荷,Location,best); xlabel(时段/h); ylabel(负荷/MW); grid on;同时可以把四个情景的总成本、碳排放量、新能源消纳量整理成一个结果汇总表这是论文结论部分的标准呈现方式。用表格输出时四舍五入的精度要注意成本列保留两位小数碳排放列保留三位比较合适。5.3 碳交易价格与配额等参数敏感性分析灵敏度分析是这个题目里最容易出彩的部分。很多论文在算例分析部分都有这么一张图横轴是碳交易价格纵轴是系统总成本和碳排放量观察两条曲线如何相交或变化。我自己复现时会把核心参数封装成函数然后写一个循环扫描碳价carbon_prices 20:10:200; total_cost_list zeros(size(carbon_prices)); co2_list zeros(size(carbon_prices)); for i 1:length(carbon_prices) carbon_price carbon_prices(i); [total_cost_list(i), co2_list(i)] run_ies_model(carbon_price); end figure; yyaxis left; plot(carbon_prices, total_cost_list, o-, LineWidth, 1.5); ylabel(系统总成本/元); yyaxis right; plot(carbon_prices, co2_list, s-, LineWidth, 1.5); ylabel(碳排放量/吨); xlabel(碳价/(元/吨)); grid on;这种参数扫描做起来很快而且能讲清楚一个关键结论随着碳价升高系统的碳排放量下降但总成本不一定单调上升。因为需求响应和碳交易相互配合时系统会通过调整负荷曲线和能源结构把高碳出力的时段避开部分成本被负荷转移带来的购能成本下降抵消。如果在 5.3 之后要自然收尾我一般会加一小段个人操作体会。比如我在跑完碳价敏感性之后发现如果只开需求响应不开碳交易系统的总成本下降有限但负荷曲线明显变平如果只开碳交易不开需求响应碳排放降得很明显但成本上升更快。两个机制同时开启时才存在一个成本可以接受同时排放也有下降的甜点区间。这个观察不一定所有论文都明确写出来但对理解模型机理非常关键。在参数扫描时多花十几分钟把这类交叉现象找出来往往比复制论文的图更有价值。最后再分享一个小技巧这类模型调试时如果遇到“infeasible problem”不要急着去翻约束逐条看先用 Yalmip 的check(Constraints)找出哪些约束的值是 NaN 或 Inf然后专门盯着那几条约束的变量量纲查。绝大多数不可行问题根源都是单位不统一导致某条平衡约束的系数差了 10 的倍数。把这一关过了剩下的求解和出图通常就顺畅很多。