这两年做园区级综合能源项目最常被问到的就是风电光伏上得多了以后供暖期怎么把弃风弃光压下去我这套“考虑可再生能源消纳的电热综合能源系统日前经济调度模型Matlab代码实现”就是干这个的。它把电力调度和热力调度放到一个优化框架里提前24小时安排火电、热电联产、储热罐和风电光伏的出力计划。搞电力系统优化、综合能源规划或者写毕业论文的同学可以从这套模型里找到能直接改、直接跑、直接出曲线的思路。这类问题表面上是数学优化实际上核心是“电”和“热”两个系统之间的耦合怎么建模、怎么求解。很多初学者一开始只把热负荷当成一个简单的等式约束结果模型跑出来要么不收敛要么弃风量算出来是负的要么机组出力曲线抖得没法解释。这篇文章我把项目的建模思路、Matlab实现步骤、求解器选择和排查经验全部整理出来希望你能少走一点弯路。1. 这个模型到底在解决什么问题1.1 传统热电联产“以热定电”是如何挤占风电空间的先说一个最常见的场景。北方冬季供暖期间大量热负荷需要由热电联产机组CHP来满足。CHP机组有一个“以热定电”的特性只要它对外供热多少电出力就不是完全自由的存在一个由热电比决定的最小电出力。供热需求越大最小电出力就越高。这时候风电就尴尬了。晚上风电大发热负荷又高CHP机组为了供暖必须保持很高的最小电出力电网消纳不了那么多风电只能调度风电出力下降也就是弃风。这就像一个大客车占了整个车位后面的小车风电光伏再想停也停不进去。解决思路也不复杂给热力系统加一些灵活性资源比如储热罐、电锅炉、热泵把“以热定电”的紧耦合关系解开来。热负荷高峰时可以让CHP多发电、多供热多余热量存进储热罐风电大发时CHP可以少发电、少供热热负荷由储热罐放热补上。这样一来CHP的电出力下限不再被热负荷死死绑住风电的消纳空间自然就出来了。1.2 日前经济调度模型在调度时间尺度和目标上的定位“日前”指的是提前24小时通常以1小时为分辨率对未来一天各机组、各储能设备的出力计划做优化安排。它是整个调度体系里的“总指挥”日内滚动调度和实时调度会在它的基础上做修正。经济调度这个“经济”字眼决定了目标函数是成本类指标最常见的是运行成本最小化包括火电机组煤耗成本、CHP机组燃料成本、机组启停成本、向外部电网购电成本以及弃风弃光惩罚成本。有些人还会加入碳排放成本、设备运维成本这些都可以根据项目需要扩展。模型求解之后输出的不是一张简单的机组出力表而是一组完整的调度计划火电每个小时发多少、CHP每小时发多少电和多少热、储热罐每小时充放多少热、风电光伏实际消纳多少、弃了多少最终形成一系列负荷平衡曲线和成本数据。1.3 谁需要跑这套模型能拿来做什么我在实际接触中使用这类模型的人群大概分三类。第一类是电力系统或热能工程的研究生课题方向是综合能源系统优化需要用模型跑出案例数据来支持论文。第二类是综合能源服务公司的工程师做园区或区域级的源网荷储规划时需要用调度模型测算经济性和消纳提升效果。第三类是电网或热力公司的规划人员想评估新建储热罐、电锅炉到底能带来多大的灵活性收益。不管哪类读者跑通模型本身不是最终目的理解模型背后的物理逻辑才是关键。代码永远只是工具你只有知道每个约束在表达什么物理过程遇到问题才知道从哪里下手排查。2. 模型怎么建核心数学结构拆开讲2.1 电源侧与热源侧的统一建模思路建模第一步是明确系统里有哪些设备。最基础的配置包括常规火电机组、热电联产机组抽凝式或背压式、风电场、光伏电站、储热罐以及外部电网联络线。有条件的还可以加电锅炉、热泵、电储能。常规火电机组的出力范围就是最大最小技术出力加上爬坡约束P_f_min ≤ P_f(t) ≤ P_f_max -R_down ≤ P_f(t) - P_f(t-1) ≤ R_up风电和光伏则有一个预测可用出力值。实际调度中不一定全部消纳所以引入弃风变量0 ≤ P_w_curt(t) ≤ P_w_forecast(t) P_w_actual(t) P_w_forecast(t) - P_w_curt(t)CHP机组是模型的关键也是新手最容易建模错误的地方。抽凝式CHP机组的电出力和热出力不是独立变量它们落在热力可行域内。工程上常用多边形近似来描述这个可行域每个顶点对应一组“最小/最大电出力-供热量”的组合。采用多边形近似后约束变成一组线性不等式也就是用矩阵的形式表达可行域边界。2.2 热力系统的关键约束热平衡、储热罐和热网动态热力侧的第一个约束是热功率平衡。在简化的集总模型中可以写成Q_chp(t) Q_boiler(t) Q_hs_dis(t) - Q_hs_char(t) Q_load(t)其中 Q_chp 是CHP供热量Q_boiler 是电锅炉或燃气锅炉供热量Q_hs_char / Q_hs_dis 是储热罐充放热功率Q_load 是热负荷。储热罐模型与电储能类似包含储热容量上、下限、充放热功率限制以及为了保证调度周期可持续性的“始末状态一致”约束SOC(t1) SOC(t) (Q_hs_char(t)·η_c - Q_hs_dis(t)/η_d)·Δt SOC_min ≤ SOC(t) ≤ SOC_max SOC(0) SOC(T)实际项目里如果热网距离较长还需要考虑热网传输延迟和热惯性。热水的流动需要时间热负荷也不是即时反映到热源端这时可以引入时间延迟和一阶惯性环节。但这一部分会显著增加模型复杂度最开始的版本不建议直接上先把静态热平衡跑通再逐步加动态特性。2.3 弃风弃光惩罚项的建模技巧目标函数中弃风弃光惩罚成本的设计直接决定模型行为。惩罚项写作C_curt λ_w · Σ P_w_curt(t) λ_pv · Σ P_pv_curt(t)λ_w 和 λ_pv 是弃风弃光的单位惩罚系数。这个系数不能随便拍脑袋。设太小模型会觉得弃风也无所谓宁可让火电少调也不肯多消纳设太大模型会为了减少弃风而让火电和储能频繁调节甚至出现为了凌晨消纳一点风电白天多烧很多煤的反直觉结果。合理的做法是让惩罚系数略高于单位煤耗成本的某个比例再通过敏感性分析确定一个“既不弃风、又不导致过度调节”的临界值。这个我后面在排查章节会详细说。2.4 目标函数各项权重如何取舍目标函数如果只写一个煤耗成本模型可能选择弃风而少调火电因为弃风惩罚没写进去。如果只写弃风惩罚模型又会忽视火电成本。所以需要合理组合min F Σ C_fuel_f(t) Σ C_fuel_chp(t) Σ C_start(t) Σ C_buy(t) Σ λ · P_curt(t) Σ C_hs_op(t)一般情况下启停成本放在目标函数里可以避免机组频繁启停。但如果你刚开始调试建议先去掉0-1启停变量只做等值连续优化跑通后再加入启停。这样排查起来会简单很多。3. Matlab代码实现从手写矩阵到YALMIP建模3.1 环境准备工具箱和求解器选型我常用的配置是 Matlab R2021a 及以上版本配合 YALMIP 工具箱做建模然后用 Gurobi 或者 Cplex 求解。YALMIP 的好处是把优化问题写成接近数学表达的形式不用自己拼装大规模矩阵也不容易错位。如果你没有 Gurobi 或 Cplex可以用 Matlab 自带的linprog先跑一个小算例但要注意含机组启停之后问题变成混合整数线性规划MILP必须用intlinprog来解。intlinprog对中等规模问题还能应付但是一旦机组数量和时段数量上去求解速度会明显变慢。安装 YALMIP 很简单下载后把文件夹加入 Matlab 路径即可。Gurobi 需要有许可证一般学生和学术用户可以申请免费授权安装后记得在 YALMIP 中用sdpsettings(solver,gurobi)指定求解器。3.2 参数结构体管理数据让你的代码不失控模型数据量很大不建议用一堆散落的变量名。我在项目里习惯把数据分成三个结构体load_data电负荷、热负荷、风电预测出力、光伏预测出力曲线长度为 T。unit_data每个机组的技术参数包括出力上限、下限、爬坡率、煤耗系数。system_data储热罐参数、弃风惩罚系数、目标函数的权重系数。举个例子设定时间分辨率为1小时调度周期 T 24系统包含2台火电机组、2台CHP、1个风电场、1个储热罐T 24; n_f 2; % 火电机组数 n_chp 2; % CHP机组数 load_data.P_e [120 115 110 108 105 100 98 95 100 110 125 140 145 150 148 145 140 135 130 125 118 115 112 110]; % 电负荷,MW load_data.Q_h [80 78 75 73 70 68 66 70 75 85 90 95 96 94 92 90 88 85 83 82 80 79 78 76]; % 热负荷,MWth load_data.P_w_forecast [50 55 60 65 60 50 40 35 30 25 20 15 12 15 18 20 25 30 35 40 45 50 55 60]; % 风电预测,MW这样当你写约束时循环里面访问load_data.P_e(t)就非常清晰后面如果要改场景只需要替换对应的曲线值不需要改动约束代码。3.3 用 YALMIP 定义决策变量和约束决策变量的定义方式如下yalmip(clear) P_f sdpvar(n_f, T, full); % 火电出力 P_chp sdpvar(n_chp, T, full); % CHP电出力 Q_chp sdpvar(n_chp, T, full); % CHP热出力 Q_hs_char sdpvar(1, T, full); % 储热罐充热 Q_hs_dis sdpvar(1, T, full); % 储热罐放热 P_w_curt sdpvar(1, T, full); % 弃风 P_w_actual sdpvar(1, T, full);% 风电实际出力然后是等式和不等式的约束集合。我习惯把所有约束放在一个对象Constraints中不断追加Constraints []; % 火电出力上下限和爬坡约束 for t 1:T for i 1:n_f Constraints [Constraints, unit_data.P_f_min(i) P_f(i,t) unit_data.P_f_max(i)]; if t 1 Constraints [Constraints, P_f(i,t) - P_f(i,t-1) unit_data.ramp_up(i)]; Constraints [Constraints, P_f(i,t-1) - P_f(i,t) unit_data.ramp_down(i)]; end end end电力平衡约束是每一时刻所有电源出力之和等于电负荷注意风电出力是实际消纳值for t 1:T Constraints [Constraints, sum(P_f(:,t)) sum(P_chp(:,t)) P_w_actual(t) load_data.P_e(t)]; Constraints [Constraints, P_w_actual(t) load_data.P_w_forecast(t) - P_w_curt(t)]; Constraints [Constraints, P_w_curt(t) 0]; end热力侧约束对应如下% 储热罐SOC SOC sdpvar(1, T, full); eta_c 0.95; eta_d 0.95; for t 1:T Constraints [Constraints, sum(Q_chp(:,t)) Q_hs_dis(t) - Q_hs_char(t) load_data.Q_h(t)]; Constraints [Constraints, Q_hs_char(t) 0, Q_hs_dis(t) 0]; Constraints [Constraints, Q_hs_char(t) system_data.hs_char_max]; Constraints [Constraints, Q_hs_dis(t) system_data.hs_dis_max]; if t 1 Constraints [Constraints, SOC(t) system_data.hs_soc0 ... (eta_c*Q_hs_char(t) - Q_hs_dis(t)/eta_d)]; else Constraints [Constraints, SOC(t) SOC(t-1) ... (eta_c*Q_hs_char(t) - Q_hs_dis(t)/eta_d)]; end Constraints [Constraints, system_data.hs_soc_min SOC(t) system_data.hs_soc_max]; end % 周期始末一致 Constraints [Constraints, SOC(T) system_data.hs_soc0];3.4 目标函数、求解和结果可视化目标函数可以写成一个求和式Objective 0; for t 1:T for i 1:n_f % 火电煤耗成本使用二次函数简化 Objective Objective unit_data.alpha_f(i) * P_f(i,t)^2 unit_data.beta_f(i) * P_f(i,t); end % CHP煤耗成本 for i 1:n_chp Objective Objective unit_data.alpha_chp(i) * P_chp(i,t)^2 ... unit_data.gamma_chp(i) * Q_chp(i,t)^2 ... unit_data.delta_chp(i) * P_chp(i,t) * Q_chp(i,t); end % 弃风惩罚 Objective Objective system_data.lambda_w * P_w_curt(t); end注意含有机组出力平方项后问题变成二次规划QP如果再加0-1变量就是MIQP。为了提高求解稳定性很多项目会把二次煤耗成本分段线性化或者在YALMIP中直接使用Gurobi的MIQP求解能力。如果二次项导致求解很慢可以先改成线性成本系数C_fuel a_i · P_i b_i这也是工程上常用的近似手段。求解部分只有一行ops sdpsettings(solver, gurobi, verbose, 2); sol optimize(Constraints, Objective, ops);求解完成后检查sol.problem是否为0然后提取结果绘图。比如电力平衡图P_w_actual_value value(P_w_actual); P_f_value value(P_f); P_chp_value value(P_chp); figure; bar((1:T), [P_f_value, P_chp_value, P_w_actual_value], stacked); hold on; plot(load_data.P_e, k-, LineWidth, 2); xlabel(时段/h); ylabel(功率/MW); legend(火电出力,CHP电出力,风电实际出力,电负荷);类似地可以画出热平衡图、储热罐SOC曲线、弃风量曲线。曲线图在论文里非常有用而且能帮你直观地检查结果有没有问题。4. 跑通模型后我踩过的坑与排查实录4.1 求解器报“Infeasible”的常见原因这是最让人头疼的报错。模型变量多、约束多一不小心就无解。我遇到过几次经验是不要盯着全部约束看而是分层排查。第一步把热力侧约束全部注释掉只留电力平衡和火电约束检查是否有解。如果有解说明问题出在热力侧。第二步把储热罐SOC约束注释掉只留热平衡约束看是否有解。如果有了说明是SOC容量参数与系统调度需求不匹配。第三步检查储热罐始末SOC一致约束很多无解都是因为初始SOC设置太高或者储热容量太小导致无法做到24小时循环。还有一种隐蔽情况热负荷数据和电负荷数据本身不匹配。比如热负荷全天恒定而风电出力凌晨很大此时CHP无法把所有热负荷转移给储热罐因为储热罐容量不够。这个其实是物理不可行不是模型问题要调整储热罐容量参数。4.2 弃风量出现负值理论上弃风变量加了大于等于0的约束就不会为负。但如果你在写等式时漏了P_w_actual 0或者没有把P_w_actual定义为非负变量那么优化器为了凑平衡可能会让风电实际出力大于预测值从而产生负弃风这在物理上毫无意义。排查方法很简单看结果里P_w_actual是否超出预测出力如果超出检查变量定义P_w_actual sdpvar(1, T, full); % 必须追加 Constraints [Constraints, P_w_actual 0];另一种情况是弃风量绝对值很小但为负可能是数值精度问题。可以设置求解器的精度容差或在目标函数中对弃风量加一个极小的二次罚项帮助收敛。4.3 含0-1变量算得太慢怎么办加入机组启停变量后模型会变成大规模MILP。24个时段、10台机组对Gurobi来说不算难但如果约束写得不好求解时间可能从几秒变成几分钟。一个常见的加速手段是减少冗余约束。不要给所有时段、所有机组都写同样长度的爬坡约束有些机组在特定时段根本不会爬坡但约束越多分支定界的过程越慢。用binvar定义启停变量后可以加一些启动和停机动作的约束u_f binvar(n_f, T, full); % 机组最小启停时间约束 for i 1:n_f for t 2:T Constraints [Constraints, u_f(i,t) - u_f(i,t-1) u_f(i, max(1,t-MinUp(i)));]; end end如果项目只关心经济调度而不关心机组组合可以直接固定所有机组开机状态先不做0-1变量。在有些场景里机组组合已经由上一个模型确定当前模型只做负荷分配那用连续变量就够了求解速度会快非常多。4.4 弃风惩罚系数怎么调才合理罚系数调参是我觉得最考验经验的部分。我一般采取两步走。第一步把惩罚系数设为0跑一次模型得到基础的弃风量。第二步逐步增加惩罚系数观察弃风量变化和总成本变化。你会看到这样一条曲线惩罚系数从0增大时弃风量快速下降总成本小幅上升但继续增大后弃风量下降缓慢总成本却加速上升说明模型为了消纳最后那一点风电付出了过高代价。在实际项目中我通常取“弃风量-成本曲线”的拐点位置作为惩罚系数。如果系统里还有储热罐还要看储热罐SOC是否经常触及上限。如果SOC频繁顶到上限说明储能容量约束在限制消纳这时候调大惩罚系数没有意义不如增加储热容量或提高充放热功率限制。4.5 可视化结果里的“锯齿形”出力怎么解释有时机组出力曲线会每小时大幅波动看起来非常不自然。这通常有两个原因。一是数据曲线本身不平滑比如电负荷预测有突变火电机组为了跟随负荷会频繁调整。二是目标函数里没有加调节成本或惩罚优化器会在不同时段间“零成本”地来回调整从而形成锯齿。处理办法是给出力变化量加一个小的惩罚项或者在目标函数中加入出力波动项的二次惩罚。比如火电出力相邻时段差值平方乘以一个小系数可以显著平滑曲线。但注意系数不能太大否则机组响应能力会被过度约束经济性变差。5. 从这套模型出发还能往哪里扩展这套日前经济调度模型是一个非常好的基础框架我在实际项目中基于它做过几个方向的扩展。最典型的是把日前模型和日内滚动模型衔接日前先安排机组组合和储热计划日内每15分钟滚动修正可以应对风电预测误差带来的不确定性。另一个方向是引入碳排放约束和碳交易价格。在目标函数中加入碳配额成本和碳交易收益后模型会自动倾向于让CHP机组多供热、少发电再配合电锅炉消纳新能源碳成本降低的效果非常明显。这个方向对政策研究和综合能源项目投资分析都很有价值。如果系统规模更大还可以把网络约束考虑进来区分多个电力节点和多个热力节点。那时模型会变成多区域耦合优化变量规模大幅增加求解效率就成了新的挑战。建议先把本文这套单节点的模型完全跑透再去加网络潮流否则很容易被复杂矩阵结构劝退。6. 总结一下我这套代码的调试心得模型本身并不难难的是让模型结果符合物理直觉。我在实际调试中最大的体会是永远不要相信第一次的求解结果一定要把约束拆开来验证。比如先单独看电力平衡是否成立再看热平衡是否成立最后看储热罐SOC是否超限。这些数据用value()提取出来后直接画图对比一秒钟就能看出问题。还有一个小技巧在跑完整模型之前先跑一个2时段或4时段的简化版本。如果2时段的模型都调不通直接跑24时段只会更难排查。小规模算例里你能手工算出来每个变量大概应该在什么范围然后和求解结果互相印证比一直盯着报错信息有效率得多。最后再分享一个关于代码组织的经验。模型跑通之后一定要把参数和结果通过结构体保存下来比如result.P_f value(P_f)方便后续处理。我见过不少同学模型跑通但没保存原始变量后面想改个图或者重新统计某个指标只能重新跑一遍优化非常浪费时间。把这些细节做好你手里的这套调度模型就能真正变成可以反复使用的分析工具。