Matlab+Yalmip实现综合能源系统低碳运行优化调度实战解析

“双碳”这两个字喊了也有几年了搞能源优化的人应该都深有体会。以前做调度盯着电费、气费把运行成本压下来就算交差现在不行了碳排放成了硬约束哪怕成本高一点点也要先把碳排降下来。很多做综合能源系统Integrated Energy SystemIES研究的朋友或者刚入行的工程师拿到一个包含光伏、风电、热电联产CHP、燃气锅炉、电储能这类设备的系统第一反应往往是不知道怎么把“低碳”和“经济性”同时塞进一个优化模型里。这篇博文就以我最近在Matlab里搭的一套低碳运行优化调度程序为例子把建模思路、目标函数设计、约束处理、求解实现和踩坑经历完整拆开讲一遍。内容偏实战适合正在写论文、做毕设、或者刚接触IES优化调度的同行参考直接照着改数据就能用。1. 程序整体设计与建模思路1.1 综合能源系统里都有哪些“家伙”拿到一个综合能源系统第一步不是急着写代码而是先把系统里的设备、能量流、信息流理清楚。我这套程序里默认的设备配置很典型光伏PV、风电WT、热电联产机组CHP、燃气锅炉GB、电储能ESS外加从上级电网购电、从天然气网购气这两个外部能源入口。这里要说清楚一个关键点CHP和燃气锅炉是“电-热”耦合的核心。CHP一边发电一边产热发电和产热之间存在强耦合关系燃气锅炉则只产热作为热力侧的补充。光伏和风电属于不可控可再生能源出力曲线基本靠预测数据给定调度模型里只能决定“用不用、用多少、弃不弃”。电储能是灵活性资源可以在电价低谷充电、高峰放电平抑可再生能源的波动。热负荷侧通常还包括一个热储能罐但为了控制程序复杂度我初期先用“电储能 直接热平衡”的方式跑通后期再叠加热储能。如果你手上的系统带了热储能建模思路和电储能几乎一样只是把能量状态换成热量状态而已。这套程序搭好之后设备增减就是改模块的事。1.2 调度模型的核心能量平衡与设备出力分配很多初学者最容易犯的错是拿起代码就写目标函数结果约束条件一加就冲突半天调不出可行解。我的习惯是先画能量平衡图再列数学表达式最后才写Matlab代码。系统的核心平衡关系分两个层面电力平衡光伏出力 风电出力 CHP发电 电储能放电 上级购电 电负荷 电储能充电热力平衡CHP产热 燃气锅炉产热 热负荷不要小看这两条平衡式整个优化模型都是围着它们转的。设备多了以后每条平衡式里还要加入损耗项、启停状态变量但那是后话。初版模型只要把这两个平衡等式写对后面的求解过程就顺了一半。这套系统的调度时间尺度我采用的是典型日调度时间分辨率为1小时总共24个调度时段。你也可以改成15分钟一个时段但那样决策变量数量会翻四倍求解时间相应拉长建议先1小时跑通再细化时间粒度。1.3 为什么选Matlab Yalmip这套组合工具选型这块我多说两句。目前IES优化调度领域的程序实现主流方案就那么几种Matlab亲生的优化工具箱、Yalmip Cplex/Gurobi、Python的Pyomo Gurobi还有人用GAMS。我最终选择Matlab纯粹是图它生态成熟、上手门槛低。Yalmip是一个建模语言工具包它最大的价值在于把“数学模型”和“求解器”解耦。你只需要用Yalmip的语法把变量、目标、约束描述出来底层调用Cplex还是Gurobi换个参数就切换了。这个特性在写论文时特别方便——审稿人质疑求解器你花两分钟就能换一个重新算。具体到版本我用的是Matlab R2022b Yalmip R20230630 Cplex 12.10。如果你手头没有Cplex用Gurobi也可以求解MILP混合整数线性规划都没问题。实在没有商业求解器先用Matlab内置的linprog跑一个去掉整数变量的LP版本也能验证模型逻辑。2. 设备建模与约束处理2.1 光伏与风电预测出力与弃用决策光伏和风电的建模是整套程序里最“简单但也最讲究”的部分。说简单是因为在调度层面我们基本不做物理建模直接采用预测出力曲线作为输入说讲究是因为弃光弃风决策处理不好模型要么过于乐观、要么约束过紧直接无解。光伏出力的简化模型是[ P_{PV}(t) P_{PV_rated} \times \frac{G(t)}{G_{STC}} \times [1 k \times (T(t) - 25)] ]其中(G(t))是实际光照强度(G_{STC}1000W/m^2)是标准测试条件光照(T(t))是光伏板温度(k)是温度系数通常取-0.004/°C左右。实际调度程序里更常见的做法是直接读入一个归一化的预测出力曲线再乘以装机容量。因为优化调度关心的是“在某时段能发多少电”光伏内部的物理细节不是关注重点。风电也类似读入预测出力曲线。但风电有个特殊性实际可用出力和调度允许出力之间存在一个“弃风”决策变量。这个变量的设置非常关键它在约束里长这样[ 0 \le P_{WT_dispatch}(t) \le P_{WT_forecast}(t) ]也就是调度后的风电出力可以在0到预测值之间任意取实际的弃风量就是预测值减去调度值。这种做法引入了“削峰填谷”的灵活性也是低碳调度的核心操作——当系统电负荷很低、但光伏风电大发时与其被强制降负荷不如主动弃掉一部分可再生能源换取整体运行成本不飙升。注意弃风弃光量一定要加进目标函数做惩罚惩罚系数不能太小否则模型会为了省成本疯狂弃风弃光但也不能太大否则等于硬性规定“必须全额消纳”模型灵活性没了。我一般取度电成本的1.5~2倍作为弃风弃光惩罚系数。2.2 热电联产机组的“电随热走”困境CHP是整个系统建模里最需要花心思的部件因为它同时横跨电力系统和热力系统而且不同技术路线的CHP运行约束完全不一样。我程序里默认采用的是抽凝式CHP它的可行运行区间可以描述为[ P_{CHP}(t) \in [P_{CHP_min}, P_{CHP_max}] ] [ H_{CHP}(t) \eta_{chp_h} \times F_{CHP}(t) ] [ P_{CHP}(t) \eta_{chp_e} \times F_{CHP}(t) ]这里的核心是热电比即产热和发电的比值。在简化模型里我们直接定义[ H_{CHP}(t) \alpha_{chp} \times P_{CHP}(t) ]这个(\alpha_{chp})就是热电比。对于背压式CHP热电比基本固定对于抽凝式CHP热电比在一定范围内可调。如果你的研究不需要做得太细固定热电比是最省事的做法如果要做深度分析可以用可行域多边形去约束P和H的耦合关系但程序复杂度会上一个台阶。CHP最麻烦的是爬坡约束和启停约束。爬坡约束是[ -P_{CHP_ramp_down} \le P_{CHP}(t) - P_{CHP}(t-1) \le P_{CHP_ramp_up} ]这个约束在Matlab里实现时需要把24小时的所有时段放在一起写用循环生成约束条件再拼接成矩阵也是最容易出错的地方。我的建议是先把一个机组的爬坡约束写对、测通再扩展到多机组。2.3 燃气锅炉与电储能灵活性担当燃气锅炉模型比CHP简单得多本质上就是一个“天然气→热量”的转换器[ H_{GB}(t) \eta_{gb} \times F_{GB}(t) ]其中(F_{GB}(t))是燃气锅炉消耗的天然气功率(\eta_{gb})是锅炉效率一般在0.85~0.95之间。燃气锅炉的唯一运行约束就是出力上下限[ 0 \le H_{GB}(t) \le H_{GB_max} ]它的作用很明确当CHP产热不足或电负荷要求CHP降出力时燃气锅炉顶上热力侧的缺口。电储能ESS是灵活性提升的关键设备建模用到经典的状态转移方程[ SOC(t1) SOC(t) \frac{P_{ESS_ch}(t) \times \eta_{ch} \times \Delta t}{E_{ESS_cap}} - \frac{P_{ESS_dis}(t) \times \Delta t}{E_{ESS_cap} \times \eta_{dis}} ]注意充放电是互斥的所以需要引入0-1整数变量[ P_{ESS_ch}(t) \le M \times z_{ch}(t) ] [ P_{ESS_dis}(t) \le M \times z_{dis}(t) ] [ z_{ch}(t) z_{dis}(t) \le 1 ]这里的(M)是一个足够大的常数在Yalmip里用约束生成的写法这三条不等式要循环写上24个小时。同时SOC本身有上下限约束通常工作在10%~90%区间不要让电池完全放空或充满这对实际寿命有利也对模型求解有利。2.4 碳排放约束双碳目标如何“落地”既然标题是“双碳目标下”的低碳调度碳排放的刻画自然不能含糊。我这里采用的是碳交易机制下的碳排放成本建模这也是目前IES低碳调度文献里最主流的写法。先算系统的实际碳排放量[ E_{CO2} E_{grid_buy} \times EF_{grid} (F_{CHP} F_{GB}) \times EF_{gas} ]其中(EF_{grid})是电网购电的碳排放因子单位是kgCO2/kWh(EF_{gas})是天然气的碳排放因子单位是kgCO2/kWh。然后引入免费碳排放配额[ E_{quota} (P_{grid_buy} P_{CHP}) \times EF_{quota} ]实际碳排放量与配额的差值就是需要去碳市场购买的碳配额量[ C_{carbon} \lambda_{CO2} \times (E_{CO2} - E_{quota}) ]当实际排放大于配额时这是一个正的碳交易成本当实际排放小于配额时这就是卖碳的收益。这个机制非常有意思它让系统有了“为减碳而减碳”的经济动机也让优化结果更有现实意义。这个模型在Matlab里实现时只要把碳交易价格(\lambda_{CO2})设成一个标量参数比如0.25元/kg碳成本就会自动进入目标函数。3. 优化模型与求解实现3.1 目标函数经济性与低碳性的博弈目标函数是整个程序灵魂所在。我最终采用的综合成本最小化目标函数包含四块购电购气成本、设备运行维护成本、碳交易成本、弃风弃光惩罚成本。[ \min \quad \sum_{t1}^{24} [C_{grid}(t) C_{gas}(t) C_{om}(t) C_{carbon}(t) C_{curtail}(t)] ]展开来看购电成本购电功率乘以分时电价电价曲线直接读入购气成本CHP和燃气锅炉的总耗气量乘以天然气价格这里要注意单位统一天然气的热值按kWh计算运行维护成本各设备出力乘以对应的单位运维成本系数光伏、风电的运维成本很低但也要算碳交易成本如2.4节所述这是一个线性表达式弃风弃光惩罚惩罚系数乘以弃风弃光电量在Yalmip里目标函数写起来非常直观。先用yalmip定义24维的优化变量然后按上面的表达式逐项叠加最后用optimize调用求解器。底层求解器会自动把问题转换成标准形式。我特别想强调的是低碳模型并不等于把所有成本都让位给碳成本。很多初学者把碳价格调得特别高导致模型为了减碳不惜一切代价最后结果一点都不经济。合理的做法是把碳价格设为实际碳市场的交易价格区间然后做敏感性分析——观察碳价从低到高变化时系统碳排放量和总成本的变化趋势这才是“双碳目标下低碳运行”研究的正确打开方式。3.2 Yalmip建模实操从变量定义到约束拼接我先贴一段精简的核心建模代码展示如何用Yalmip搭建这套优化模型。这里只展示框架完整代码建议自己动手补全补的过程中对约束逻辑的理解会深很多。%% 加载负荷与可再生能源预测数据 load(IES_data.mat); % 包含 P_load, H_load, P_pv, P_wt, 电价, 气价等 %% 定义决策变量 P_Grid sdpvar(1, 24, full); % 购电功率 F_CHP sdpvar(1, 24, full); % CHP购气功率 F_GB sdpvar(1, 24, full); % 燃气锅炉购气功率 P_CHP sdpvar(1, 24, full); % CHP发电功率 H_GB sdpvar(1, 24, full); % 燃气锅炉产热功率 P_ESS_ch sdpvar(1, 24, full); % 电储能充电功率 P_ESS_dis sdpvar(1, 24, full); % 电储能放电功率 SOC sdpvar(1, 25, full); % 电储能荷电状态25个点包含初末 z_ch binvar(1, 24, full); % 充电状态 z_dis binvar(1, 24, full); % 放电状态 % 弃风弃光变量也可以定义出来这里省略定义好变量后约束条件用循环或向量化方式拼接。下面这段是电储能SOC递推约束的写法Constraints []; for t 1:24 % SOC 状态递推 Constraints [Constraints, ... SOC(t1) SOC(t) (P_ESS_ch(t) * eta_ch - P_ESS_dis(t) / eta_dis) / ESS_cap]; % 充放电互斥 Constraints [Constraints, ... P_ESS_ch(t) 0, P_ESS_ch(t) ESS_cap * 0.5 * z_ch(t)]; Constraints [Constraints, ... P_ESS_dis(t) 0, P_ESS_dis(t) ESS_cap * 0.5 * z_dis(t)]; Constraints [Constraints, ... z_ch(t) z_dis(t) 1]; end这里把充放电功率上限设为了ESS容量的0.5倍对应2小时充满/放空倍率。实际参数要根据具体电池特性去改别照抄。SOC的初值设定为0.2末值要求回到0.2这是为了保证调度策略的可持续性不会出现一天结束电池电量被完全耗尽的情况。3.3 求解器配置与求解结果验证求解器配置这块也有不少坑。Yalmip调用Cplex的语法是ops sdpsettings(solver, cplex, verbose, 1, showprogress, 1); optimize(Constraints, Objective, ops);这里verbose设为1可以看到求解过程日志方便判断模型是否遇到数值问题。如果出现Infeasible problem大概率不是求解器的问题而是约束冲突需要回去检查约束条件。我的建议是刚开始调试时先用一个非常宽松的约束集合去跑比如去掉爬坡约束、把储能容量设得很大确认模型能出结果然后逐步收紧每加一组约束就重新求解一次这样能快速定位是哪个约束把模型“卡死”了。这个调试流程是我实测效率最高的方式比自己盯着代码猜要快得多。求解完成后一定要检查结果是否满足所有约束尤其是电功率平衡和热功率平衡这两条等式约束。我曾经因为负荷数据单位不一致一个用MW、一个用kW导致平衡等式有偏差结果曲线看起来合理实际电量差了好几兆瓦时。这种错误最坑因为不细看根本发现不了。建议在求解后增加一段校验代码residual_e P_PV P_WT P_CHP P_ESS_dis P_Grid - P_Load - P_ESS_ch; residual_h H_CHP H_GB - H_Load; fprintf(最大电功率不平衡量: %.4f kW\n, max(abs(residual_e))); fprintf(最大热功率不平衡量: %.4f kW\n, max(abs(residual_h)));这两行校验代码效率极高能帮你快速发现建模时手滑写错符号、漏加变量之类的问题。4. 结果分析与项目扩展方向4.1 典型日调度结果解读怎么看出模型的合理性求解完成后画图是必须的。我的惯例是画三张图第一张是电功率平衡图堆叠展示光伏、风电、CHP、储能、购电分别承担了多少负荷第二张是热功率平衡图看CHP和燃气锅炉的产热分配第三张是SOC曲线和碳成本曲线看储能充放电策略和碳排成本随时间的变化。从结果里重点观察三点分时电价引导的储能充放电策略是否合理、CHP在电负荷高峰和热负荷高峰时的出力是否合理、弃风弃光时段是否符合预期。举个例子如果电价低谷时段出现在凌晨而储能SOC在凌晨不涨反降那大概率是SOC初值或约束写错了。关于低碳效果评估我习惯引入两个指标碳减排率对比不引入碳交易机制时的碳排放量和综合成本变化率对比纯经济调度时的运行成本。这两个指标放在论文的结果分析里也好看评审一看就知道你的模型价值在哪里。我实测的一组典型数据里引入碳交易机制后碳排放量降低约12%而总成本上升约4%——这个权衡关系是符合预期的。4.2 常见问题与排查技巧实录前面说了很多经验这里集中整理一份常见问题速查表都是我在实际调试这套程序时真实遇到过的问题和对应的排查方法问题现象可能原因排查方法求解器报Infeasible约束过紧或设备参数冲突先用大范围约束测试最小模型再逐组加约束定位结果里SOC剧烈振荡SOC递推方向写反或eta取错检查递推等式的符号验证初末SOC约束是否生效弃风弃光量异常大惩罚系数太小将惩罚系数维持为度电成本的2倍重跑碳成本占目标函数比例过高碳价设置不合理做碳价敏感性分析画出碳价-碳排放曲线平衡校验最大残差1kW单位不统一或漏加设备项打印每个时间段的残差回溯到对应数据表检查求解时间过长整数变量太多或时间粒度太细先用1小时粒度跑通再考虑细化同时设定求解时间上限结果对初值过于敏感SOC初值设定不合理让SOC初末值相等或把初值放在可行域中间还有一个非常实用的技巧就是给每个设备单独跑一次“单机出力测试”。比如把所有其他设备关掉只保留电网购电和电负荷看模型能不能正确输出购电曲线再单独测试CHP单独供电供热时价格变化如何影响出力比例。这些测试能帮你快速把问题定位到单个模块而不是整个系统。4.3 从固定场景到多维扩展程序还能怎么改这套程序我用下来最大的感受是它扩展性很强。如果你不想止步于一个“跑通”的固定算例下面几个方向的改动难度都不高但效果会明显提升多场景调度把典型日换成夏冬两季、极端天气场景对比不同场景下的机组出力和碳排表现论文价值更高多能源耦合叠加加入氢能、地热、生物质等新设备模型框架不用大改只需新增对应的能量平衡方程和设备约束不确定性优化光伏风电预测本身就带误差把确定性优化改成两阶段鲁棒优化或随机优化程序整体结构不变但目标函数和约束要分层处理碳捕集装置给CHP加装碳捕集设备碳排放量变成可控变量碳成本模型会变得更丰富这些扩展方向在结构上都和现有程序兼容改动量取决于你研究的侧重点。我的建议是先把当前这套“预测-优化-校验-绘图”的闭环跑顺再决定往哪个方向深入。不要一上来就追热点加一堆新机制模型过于复杂后调试成本会指数级上升。最后分享一个我在实际运行中养成的习惯每次跑完优化把关键结果存成结构化数据文件包括目标函数各项成本明细、各设备逐时出力、SOC曲线、碳排放数据。积累几次不同场景的实验数据后横向对比分析就变得特别容易。这也为我后来做敏感性分析、写论文、给导师汇报省下了大把时间。