1. 为什么做这个模型问题拆解与研究动机先说个实际场景。过去电网调度只看电热力调度只看热两个系统各管各的。但随着风电、光伏装机量越来越大尤其是北方地区冬季供暖期电网低谷时段的风电出力恰恰是最猛的时候而热电机组因为要顶着供热任务没法压低电出力结果就是“弃风”——风机眼睁睁看着风在吹却不能发电上网。这是典型的“以热定电”约束造成的系统刚性。电热综合能源系统的核心思路就是把电网和热网放在一个框架里统一调度。通过热电联产机组CHP、电锅炉、储热罐这些耦合元件让电和热可以互相转化、互相补偿。夜间风电多的时候让电锅炉把多余的电变成热存起来或者让储热罐吸收热电机组多产的热量这样热电机组就能放心压低电出力给风电让出空间。白天负荷高的时候再把存储的热能释放出来。这一套逻辑就是整个模型的灵魂所在。所谓“日前经济调度”指的是在当天运行前一天基于日前预测的风电、光伏出力曲线、电负荷曲线、热负荷曲线制定未来24小时或者96个时间段每个时段15分钟内各设备的出力计划。它属于“前瞻性调度”和实时调度、日内滚动调度是不同时间尺度的配合关系。日前调度做得好不好直接关系到第二天能不能既安全又经济地把系统跑下来。这个课题适合谁来研究一是电气工程、热能工程、控制科学方向的研究生尤其是做综合能源系统、微电网、智能电网调度的课题组二是电网调度中心、综合能源服务公司的工程师需要做一些仿真验证和方案比选三是对优化建模、Matlab编程感兴趣的科研入门者。这篇文章会从问题的数学建模开始逐步讲到Matlab代码实现最后附上我踩过的坑和调试经验尽量让你能照着落地。2. 模型整体设计目标函数与约束条件体系2.1 设备模型系统里都有哪些元件在写目标函数之前必须先搞清楚系统里有哪些设备、每个设备的运行边界是什么。我设计的这个典型电热综合能源系统包括风电、光伏可再生能源出力作为系统的输入具有一定的预测误差但在日前调度中通常按预测值处理只考虑出力上限约束。热电联产机组CHP既能发电又能发热是电热耦合的关键环节。它的运行区间不是矩形而是一个多边形或四边形可行域电出力和热出力之间是强耦合的。纯凝汽式火电机组只发电不发热作为系统的常规电源。电锅炉用电制热是消纳过剩可再生能源的“海绵”。储热罐热储能设备可以存热、放热起到时间平移作用。电负荷和热负荷系统的刚性需求必须实时满足。这些元件里面CHP和储热罐是调度的核心变量枢纽。理解了这两个设备的特性模型就搞懂了一半。2.2 目标函数经济性从哪些成本来衡量日前经济调度的目标是在满足所有运行约束的前提下让系统一天的总运行成本最小。这个成本通常包括以下几项购电成本从上级电网购电的费用可以用分时电价表示峰谷价差越大对调度策略的影响越明显。购气成本燃气轮机或燃气锅炉的燃料成本由天然气消耗量和气价决定。CHP机组的燃料成本这是最复杂的一项。CHP的燃料成本同时取决于电出力和热出力不是简单的关于单一变量的函数。弃风弃光惩罚成本为了体现消纳可再生能源的优先级可以在目标函数里加入一个惩罚项当系统被迫削减风电、光伏出力时产生惩罚成本。这个惩罚系数要设置得足够大但也不能太大否则会造成数值病态以保证系统优先调度可再生能源。运维成本各设备的启停和运行维护成本通常是出力的线性函数。目标函数写成数学形式类似min F 购电成本 燃料成本 弃风弃光惩罚 运维成本这里有一个关键设计决策要不要把碳排放成本考虑进去我建议基础版本先不加因为碳约束会增加变量和求解复杂度而且很多入门文献也不含这部分。如果你想做扩展可以在约束里加一个总的碳排放上限作为可选的“碳约束模式”后面我会在扩展方向里详细谈。2.3 约束条件系统的“规矩”清单约束条件是调度模型的骨架我把它们分成几个子模块1电功率平衡约束每个时刻所有电源出力之和加上从上级电网购电必须等于电负荷加上电锅炉的用电量。这是最基本的平衡方程。2热功率平衡约束CHP热出力、电锅炉热出力、储热罐放热之和必须等于热负荷加上储热罐充热。注意充放热的符号定义别把方向搞反了。3CHP机组可行域约束这是整个模型里最核心、也最容易写错的地方。CHP的发电功率和供热功率不是独立调节的而是落在图中的一个多边形区域内。典型可以用四个或五个线性不等式描述例如电出力上下限随热出力变化热出力上下限最大电出力随热出力增加而下降的斜线最小电出力随热出力变化的约束这个区域如果描述不准要么调度结果不切实际比如发电功率和供热功率组合落在了矩形“外挂”的角落要么根本无解。后面我会给出一个标准的约束不等式组。4电锅炉模型电锅炉的输入是电力输出是热量转换效率一般在0.95以上。它的约束包括用电功率的上限以及产热量与用电量之间的线性关系乘以效率系数。5储热罐模型储热罐的状态转移方程储热量(t1) 储热量(t) 充热功率 × 效率 − 放热功率 ÷ 放热效率这是一个跨时段耦合约束也是让模型“记住”历史状态的关键。另外还要避免同时充放热现实中很少有设备一边充一边放满发状态下这样做也是浪费。通常引入一个二值变量来锁定充热/放热状态但这会引入整数变量增加求解复杂度。6机组爬坡约束CHP和纯凝机组的电出力变化速率有上限相邻两个时段的出力差不能太大。这在日前调度里尤其重要如果忽略生成的调度计划可能在分钟内根本无法执行。7可再生能源出力约束风电和光伏出力必须在0和预测值之间。目标函数里的弃风惩罚项会让优化器“尽量”让它们跑满但受系统消纳能力和安全约束限制必要时仍然允许削减。我把主要的约束整理成一张表方便对照检查约束类型数学描述简写说明电平衡ΣP_gen P_buy P_load P_eb每个时段都成立热平衡H_chp H_eb H_dis H_load H_chg注意充放热方向CHP可行域A·[P_chp; H_chp] ≤ b线性不等式组核心电锅炉H_eb η_eb × P_eb0 ≤ P_eb ≤ P_max效率约0.95储热罐S(t1) S(t) Q_chg·η_c − Q_dis/η_d跨时段耦合储能容量S_min ≤ S(t) ≤ S_max上下限约束同时充放热u_chg u_dis ≤ 1二值变量约束爬坡−Ramp ≤ P(t1) − P(t) ≤ Ramp相邻时段限制可再生上限0 ≤ P_wind ≤ P_wind_forecast预测值上限每一条约束都对应到实际系统的物理限制写代码的时候务必一条一条地核对公式和变量的方向我做过的项目里约有一半的“无解故障”最终都是方向写反了或符号写错了。3. 模型求解思路与工具链选型3.1 从非线性问题到可求解问题线性化是怎么玩的电热综合能源系统的调度问题本来带有非线性因素最典型的是CHP燃料成本曲线是一个凸二次函数。直接求解混合整数二次规划MIQP在系统规模小的时候没问题但一旦机组数量增加、时段粒度变细求解时间会非常难看。在实际工程里大家通用的做法是分段线性化。把二次燃料成本曲线划分为若干段每段用线性函数近似然后再引入二进制变量表示当前激活哪个段。这样做之后整个模型就从MIQP转化成了混合整数线性规划MILP。MILP的求解技术已经非常成熟Cplex和Gurobi这类商业求解器可以在毫秒到秒级的时间内搞定大部分实际问题。这里我给你的建议是基础版直接用分段线性化不要贪心去保留非线性精度。原因有二一是调度问题的精度本来就受限于预测误差没必要在成本曲线上追求千分之一的精度二是MILP的鲁棒性和求解速度对于跑参数扫掠、多场景对比这些操作体验完全不是一个级别。3.2 为什么用Matlab Yalmip Cplex这套组合工具链的选择直接决定了你是在做“原型验证”还是在做“工程落地”。我的习惯组合是Matlab Yalmip 求解器Cplex或Gurobi。这套组合的好处在于Matlab矩阵运算和数据处理方便适合快速搭建模型、做可视化分析。Yalmip一个建模工具箱让你用接近数学公式的写法来描述优化问题不需要手动把约束展开成矩阵形式调试效率提升非常大尤其适合模型还处于频繁调整阶段的场景。Cplex / Gurobi商业求解器性能强悍接口支持好。学术用户通常可以申请免费授权Gurobi的学生版和试用版也很好获取。很多人一上来就用Cplex的纯API或者MATLAB的intlinprog写大模型结果光是矩阵装配就要写几百行不说稍作修改就得重写一大片代码。Yalmip把这些事情的繁琐程度降低了一个数量级。一个最小化的Yalmip模型框架长这样% 定义决策变量 P_chp sdpvar(24, 1); % CHP电出力 H_chp sdpvar(24, 1); % CHP热出力 P_eb sdpvar(24, 1); % 电锅炉用电 H_storage sdpvar(24, 1); % 储热罐充放热正为放热负为充热 S sdpvar(25, 1); % 储热罐状态量 u_chg binvar(24, 1); % 充热状态二值变量 u_dis binvar(24, 1); % 放热状态二值变量 % 约束条件以电平衡为例 constraints []; constraints [constraints, P_chp P_wind P_buy P_load P_eb]; % 目标函数 objective sum(fuel_cost(P_chp, H_chp)) sum(price_buy .* P_buy) ... penalty * sum(P_wind_forecast - P_wind); % 求解 optimize(constraints, objective, sdpsettings(solver, cplex, verbose, 2));从上面这段就能看出Yalmip的好处是代码几乎和数学公式一一对应读代码的人不需要“翻译”就能理解模型的含义。后面我会把这个框架扩展到完整版。3.3 求解参数设置别让求解器默认配置拖后腿很多新手忽略求解参数设置直接调用optimize然后吐槽“怎么跑了几分钟还没结果”。我常用的调整项有MIP gap相对间隙Cplex默认是0强求最优是好事但代价可能是额外的求解时间。实践中我通常把它设置为0.011%也就是允许次优解偏差1%求解时间能减少一个数量级。调度问题对成本的微小偏差并不敏感但求解时长影响的是你迭代方案的效率。时间上限设置一个合理的时间上限比如60秒跑不完就返回当前最优可行解。这样在批量跑场景时不会因为某一个场景卡住而全军覆没。数值容差如果模型scale差异很大比如储热罐存储量上千GJ而电出力只有几十MW可以开启求解器的scaling选项或者手动对变量做归一化减少数值问题导致的不可行判定。Yalmip里设置方法options sdpsettings(solver, cplex, ... cplex.mip.tolerances.mipgap, 0.01, ... cplex.mip.limits.timelimit, 60, ... verbose, 2);4. 核心代码实现从数学模型到Matlab落地4.1 参数准备数据从哪里来、怎么生成做日前调度需要的数据包括24小时电负荷预测曲线典型日为双峰曲线早峰和晚峰。24小时热负荷预测曲线供暖期特点是夜间高、白天相对低基本跟气温变化反着来。风电、光伏预测出力曲线风电夜间强光伏白天有出力。分时电价峰谷平三个时段价格差异明显。设备参数CHP爬坡速率、储热罐容量、电锅炉功率上限等。如果你手头没有实测数据先自己生成“人工典型日”数据跑通代码这是最务实的做法。实测数据可以后续替换。这里我强调一个技巧把参数写成结构体或单独建一个case_data.m文件。我之前见过不少人把数据全写在主脚本里后面参数一多脚本变得非常难维护。建议的代码结构是main.m——主调度脚本case_data.m——所有系统参数和预测数据的定义build_model.m——构建目标函数和约束返回约束集合和目标函数solve_model.m——调用求解器返回结果结构体plot_results.m——结果可视化这么做的好处是后期做参数敏感性分析、多场景对比时只需要修改case_data.m或循环调用前三个函数就行。4.2 核心建模代码详解我直接贴一段经过简化的核心建模代码然后逐段解释关键点%% 1. 加载数据和参数 run(case_data.m); %% 2. 定义变量以24时段为例 T 24; P_chp sdpvar(T, 1); % CHP电出力 H_chp sdpvar(T, 1); % CHP热出力 P_wind sdpvar(T, 1); % 风电并网功率 P_pv sdpvar(T, 1); % 光伏并网功率 P_buy sdpvar(T, 1); % 购电功率 P_eb sdpvar(T, 1); % 电锅炉电功率 H_eb sdpvar(T, 1); % 电锅炉热功率 H_sto_chg sdpvar(T, 1); % 储热罐充热功率 H_sto_dis sdpvar(T, 1); % 储热罐放热功率 S_sto sdpvar(T1, 1); % 储热罐储热量含初始时刻 % 二值变量充放热状态互斥 u_chg binvar(T, 1); u_dis binvar(T, 1); %% 3. 约束条件 C []; % 3.1 电功率平衡 C [C, P_chp P_wind P_pv P_buy P_load P_eb]; % 3.2 热功率平衡 C [C, H_chp H_eb H_sto_dis H_load H_sto_chg]; % 3.3 CHP可行域以典型四边形可行域为例 % P_chp_min k_pm * H_chp P_chp ... % P_chp_max - k_pp * H_chp P_chp ... C [C, P_chp P_chp_min km_chp .* H_chp]; % 最小电出力约束 C [C, P_chp P_chp_max - kp_chp .* H_chp]; % 最大电出力约束 C [C, H_chp H_chp_min]; % 最小热出力 C [C, H_chp H_chp_max]; % 最大热出力 % 3.4 电锅炉模型 C [C, H_eb eta_eb .* P_eb]; C [C, 0 P_eb P_eb_max]; % 3.5 储热罐模型 C [C, S_sto(2:end) S_sto(1:end-1) eta_chg .* H_sto_chg - H_sto_dis ./ eta_dis]; C [C, S_min S_sto S_max]; C [C, 0 H_sto_chg H_sto_max .* u_chg]; C [C, 0 H_sto_dis H_sto_max .* u_dis]; C [C, u_chg u_dis 1]; % 3.6 可再生能源约束 C [C, 0 P_wind P_wind_forecast]; C [C, 0 P_pv P_pv_forecast]; %% 4. 目标函数 % CHP燃料成本分段线性化假设可分为3段 % 这里简化写法实际需要根据分段点构造成本系数 C_chp_fuel chp_fuel_cost(P_chp, H_chp); C_buy sum(price_buy .* P_buy); C_penalty_wind wind_penalty * sum(P_wind_forecast - P_wind); C_penalty_pv pv_penalty * sum(P_pv_forecast - P_pv); C_om sum(om_chp .* P_chp) sum(om_eb .* P_eb); objective C_chp_fuel C_buy C_penalty_wind C_penalty_pv C_om; %% 5. 求解 options sdpsettings(solver, cplex, verbose, 2, ... cplex.mip.tolerances.mipgap, 0.01, ... cplex.mip.limits.timelimit, 60); sol optimize(C, objective, options); if sol.problem 0 % 提取结果 P_chp_opt value(P_chp); H_chp_opt value(H_chp); % ... 其余变量类似 else disp(求解失败请检查约束和参数); disp(sol.info); end4.3 几个关键实现细节1CHP可行域的坑上面代码里CHP可行域简化成了四边形。但实际中CHP可行域可能是五边形甚至更复杂包含背压工况、抽汽工况等模式。简单起见可以先从四边形做起跑通了再增加复杂度。重要的是理解CHP的电出力上下限都和热出力相关而不是固定的常数。新手最容易犯的错误就是把CHP写成“电出力在一定范围内、热出力在一定范围内”的矩形约束这样得到的结果在物理上不可行。2储热罐状态变量的维度我在代码里把S_sto定义成了T1维从第1个到第25个时刻其中第1个是调度起始时刻的初始储热量第25个是调度结束时刻的终态储热量。如果不加循环约束优化器可能会把储热罐吸干在最后一个时段把热量全放掉这一点在现实中是不允许的或者至少不鼓励。加一个周期边界约束S_sto(1) S_sto(end)会好很多表示一天调度结束时储热量恢复到初始水平。3二值变量与“大M”技巧同时充放热约束用了两个二值变量u_chg和u_dis。如果储能设备比较多二值变量的数量会线性增长求解时间随之上升。另一种方案是让充放热功率共享一个变量正值为放热负值为充热这样就不需要二值变量了。代价是这种模型没法区分充热效率和放热效率如果两者不同。工程上如果你不需要区分充放效率尽量用带符号的单一变量方案求解速度更快。5. 常见问题与排查技巧实录这一节是我最想写给实操者的部分。理论模型看起来清清爽爽但只要一跑到代码里各种妖魔鬼怪就会冒出来。5.1 求解器返回“不可行”infeasible这是最常见的问题。模型从数学上看天衣无缝求解器却直接说“无解”。我自己的排查套路按顺序检查第一步检查变量范围和数据类型。有没有变量被定义成了整数变量但实际上是连续的有没有变量漏掉了取值范围最典型的情况是忘记给P_buy加非负约束Yalmip里如果默认变量是自由变量优化器可能通过负购电即售电来强行满足平衡。第二步逐步注释约束定位问题。把约束群逐一注释观察哪一组约束去掉后问题变可解。这是一个二分查找式的调试流程我做过最难的案例逐条排查了2个多小时最后发现是CHP可行域里两个不等式符号反了。第三步检查数据和量纲。电负荷单位是kW热负荷单位是GJ满满的都是量纲不一致的坑。储热罐容量和时段放热功率之间的单位必须统一换算。第四步看松弛解relaxation。先用solvesdp把整数变量全都松弛成连续变量看看松弛问题是否可解。如果松弛问题都不可行那一定是约束本身冲突如果松弛可解但整数问题不可行那就是整数约束太紧的问题。5.2 求解很慢或者卡住常见原因一MIP问题规模太大。比如把96个时段、多个机组、多个储能设备全都考虑进去二值变量可能有几百个。对策是先跑24时段版本验证逻辑再扩展96时段MIP gap从0放宽到0.01或0.05减少不必要的二值变量。常见原因二数值病态。设备参数之间数量级差距太大比如热负荷以GW计而风电只有几MW矩阵条件数很差。对策是对变量做归一化或者统一量纲到同一数量级比如全部用MW和MWh。常见原因三求解器参数没配置。这是最可惜的浪费。同一模型Cplex默认参数下可能要5分钟设置好MIP gap和时间限制后可能30秒就出结果。我这里做一个速查表方便大家对照排查症状可疑原因排查方法解决建议直接报不可行约束符号方向反了逐条注释约束定位核对数学推导直接报不可行量纲不统一检查单位和系数数量级统一单位归一化求解时间过长二值变量太多统计整数变量数量改用单变量储能方案放宽MIP gap解不合理储热罐空放缺少周期边界约束检查终态储热量加S(1)S(end)约束CHP出力超出可行域可行域不等式写错画出可行域验证检查各顶点坐标是否满足约束风电明明充裕却被弃惩罚系数太小查看弃风惩罚值占比调大惩罚系数求解器不收敛数值病态查看KKT残差和scaling信息归一化数据和变量5.3 结果违反直觉时的“拍脑袋检查法”有一次我的模型跑出来凌晨3点让电锅炉满功率运行但同时CHP还在高电出力状态。乍一看觉得不合理但仔细一分析凌晨3点风大、电价低让电锅炉产热存起来是正确的CHP高电出力是因为热负荷在那段时间确实高而且爬坡约束限制它不能快速降下来。这个结果其实是系统的优化选择不是模型出错了。碰到反直觉的结果先不要急着改模型而是做三件事画出各个设备出力的时序图看相互之间是否满足平衡关系手工算几个关键时段的功率平衡验证结果的自洽性对比基准工况比如没有储热罐、没有电锅炉的调度结果看储能设备加入后是否改善了经济性。5.4 一个我复现过的真实调试案例去年我帮一位师弟调试他的模型卡了整整两周现象是加了储热罐之后系统总成本反而上升了。这个结果理论上不可能因为储热罐给了调度更多的自由度成本只可能不增或降低。最终定位到问题出在储热罐的初始储热量设置上。他给初始储热量设了一个很大的值导致优化器发现“已经有很多热量可以用了为什么还要烧锅炉”于是热出力显著下降但为了让储热罐在调度周期结束时回到初始值最后一个时段需要大量充热。充热来自电锅炉电锅炉的用电又推高了购电成本一来一回反而让总成本变高了。虽然模型本身没问题但这提醒我们初始储热量的选取本质上相当于把前一天到当天交接的“历史状态”作为输入。如果你关注的是储热罐对调度策略的影响初始储热量应该设置成让它有调节空间的中间值而不是随便填一个数。这个坑我敢说90%的新手都会踩。6. 仿真结果分析与扩展方向6.1 结果可视化应该看哪些图、关注什么调度模型跑完后结果的呈现方式其实很重要。我自己常用的绘图包括第一张图电功率平衡堆叠图。横轴是24小时纵轴是功率风电、光伏、CHP、购电、电锅炉用电分别用不同颜色堆叠电负荷曲线作为总需求叠加展示。这张图能最直观地看出风电在夜间是否被消纳。第二张图热功率平衡堆叠图。类似地展示CHP热出力、电锅炉热出力、储热罐充放热的时序变化。第三张图储热罐储热量曲线。看储热罐一天内的充放节奏正常情况应该是夜间充热、白天放热的“削峰填谷”规律。我的习惯是写一个plot_results.m函数把所有变量结果传进去自动生成这几张图。调试的时候每跑一次模型只要扫一眼这三张图模型的合理性就心里有数了。6.2 模型扩展方向这个框架能怎么“长”基础版模型跑通之后你可以根据自己的研究需要往以下方向扩展方向一考虑不确定性。风电、光伏的预测误差是实际调度无法回避的问题。可以在日前调度的基础上叠加一个日内滚动修正环节形成“日前计划日内调整”的两阶段架构。更学术的做法是改成随机规划或鲁棒优化模型用场景法把多个可能的预测情景塞进来。方向二多区域互联。如果多个园区或多个微网之间可以交换电能和热能模型就要加入联络线功率约束、跨区域交易成本等。这个方向适合做分布式协同调度的课题。方向三加入碳排放约束。在目标函数中增加碳排放成本或总碳排上限会让调度决策在“经济最优”和“低碳最优”之间做一个权衡。这个方向很适合当前双碳背景下的研究选题。方向四用户侧需求响应。把一部分负荷变成可削减或可转移的弹性负荷比如电动汽车充电、蓄热电采暖、可中断工业负荷这样系统调度又多了可控资源。方向五设备模型精细化。电锅炉的启停次数、寿命损耗储热罐的散热损失CHP的变工况效率曲线等都可以从恒定效率改成随运行状态变化的非线性模型不过这会明显增加求解难度。6.3 给新手的两个实用建议建议一从抄作业开始但要“抄懂”。刚入门的时候找一篇公开代码的论文把代码下载下来一行一行地运行、理解、修改。不要急着从零写因为调度模型的初值、量纲、约束细节太容易出错了有参考代码能节省大量“调bug”的时间。但抄完之后一定要做一次“无参考复现”——关掉原代码自己从头写一遍这个过程才能检验你是不是真懂了。建议二用基线结果做对照。在你调模型的过程中时刻保留一个最简单的“对照组”比如只有CHP和电锅炉没有储热罐。每一次加约束、加设备都和对照组的结果对比看看成本怎么变、风电消纳怎么变。这不仅是验证模型合理性的手段也是论文里写“敏感性分析”的好素材。我个人对这个方向最大的体会是电热综合能源系统的调度模型本质上是一个“多能源品种协同的问题”。你写的每一行约束背后都是一笔物理账优化器给出的每一个决策背后都是一笔经济账。把这个账算清楚模型就离实际可用不远了。另外补充一个自己在实际项目中常用的调试技巧跑场景对比的时候用save把每次的结果文件保存下来命名带上日期和场景标识比如result_20250110_wind_high.mat。这样过几天你回头看还能记起来当时跑的是哪组配置。别问我为什么强调这一点——我已经数不清有多少次因为结果文件被覆盖被迫重跑一整天场景的日子了。
