Matlab联合优化调度：破解热电联产机组风电消纳难题

1. 项目背景与核心问题拆解

1.1 热电联产机组为什么会“堵住”风电消纳

做风电、火电联合调度的朋友应该都遇到过这个场景：冬季夜间，风电场满发，但全网用不完这些电，于是调度只能下令弃风。明明是可再生能源，白白扔掉，看着就心疼。问题出在哪？很大程度出在那批“以热定电”的热电联产机组（CHP）身上。

北方城市供热季是靠热电厂带基础热负荷的。供热管网需要稳定的热源，而CHP机组的运行模式是“以热定电”——只要供热量确定了，发电出力就被机组的电热特性曲线绑死了。冬季夜间恰好是供热的峰值时段，CHP机组必须高负荷运转，发出来的电却正好撞上电网负荷低谷和风电大发期。电多、用不掉，风电就得让路。这个矛盾我前几年在北方某省份的实际调度数据里看得非常直观：三月供热季结束前后，弃风率能差出五六个百分点。

传统思路是上电锅炉、蓄热罐这类设备来解耦热电耦合，但这涉及一台台机组怎么搭配运行的问题。单纯加设备，不解决机组之间的协同调度，投资打了水漂，风电照样弃。这篇项目做的就是一个“联合优化控制”框架：把多台CHP机组、电加热装置、电能负荷、热负荷和风电出力放在同一个优化模型里，用Matlab求解出未来若干个时段内每台机组的出力计划，目标就是在保证供热和供电的前提下，最大程度把风电消纳掉。

1.2 这个优化控制的本质是在做什么

剥开“联合优化控制”这层外壳，本质就是一个带约束的多时段经济调度问题。我更愿意把它理解为“给一整片电热联合系统做一个未来24小时的动作编排计划”：

• 风电出力多出来的时段，怎么压CHP电出力、怎么用储热来顶上供热缺口；
• 风电出力小的时段，怎么让CHP多发电、给储热罐蓄热，给后续时段留出调节空间；
• 所有动作都要满足电网功率平衡、热网平衡、机组爬坡速率、出力上下限等硬约束。

项目实施路径分成四步：先梳理CHP机组的电热运行可行域，再搭目标函数和约束的数学表达式，然后通过Matlab调用规划求解器完成优化计算，最后做多场景仿真验证模型有效性。下面我把每一个环节的实操细节都展开说清楚，包括我踩过的坑。

2. 数学模型搭建方法

2.1 目标函数怎么选才符合工程直觉

目标函数是优化模型的“指挥棒”，定得对不对直接决定结果是否符合调度预期。这个项目里有一个容易被忽略的细节：最大化风电消纳和最小化系统运行成本，表面上方向一致，实际建模时处理方式不同。

“最大化风电消纳”最直接的写法是让目标函数里包含一个风电出力最大化项，比如最大化每个时段风电场并网出力之和。但在实际工程调度里，单纯追求风电消纳会导致求解器利用一切手段压CHP出力，哪怕某些CHP机组的运行效率很高、成本很有优势，也被一刀切地压到最低，这不符合经济性要求。

我见过不少初学的朋友在这里翻车：目标函数只写风电消纳最大化，求解出来的方案是“不管成本多高都要弃煤保风”。等对上实际考核指标（比如煤耗、供热保障率）时方案根本没法落地。

这个项目采用的替代做法是把目标函数改为最小化系统运行总成本，同时在成本的各项组成中引入“弃风惩罚成本”：

• 弃风成本系数设置得比常规机组发电成本高一个量级，这样求解器在权衡“压CHP出力节省的煤耗成本”和“弃风惩罚成本”时，会主动选择优先消纳风电；
• 但当风电已经无法全额消纳时，系统才会选择弃风，避免无限制牺牲经济性。

这种间接建模方式叫“惩罚系数法”，是工程里最实用、最不容易出问题的目标表达方式。具体写作：

弃风惩罚成本 = K_penalty * sum( Wind_available(t) - Wind_used(t) )

其中 K_penalty 取值要根据系统具体情况调整，通常设定在500~1000元/MWh量级，远高于火电边际煤耗成本。热负荷侧如果也追求经济性，目标函数里再加一个恒温控制惩罚项，防止为了省电费导致供热温度波动过大。

2.2 弃风弃在哪里？风电并网约束的表达

风电消纳问题在模型里的核心约束是风电场并网功率约束：

0 <= P_wind_used(t) <= P_wind_available(t) P_wind_used(t) <= P_load(t) + P_electric_boiler(t) - P_chp_sum(t)

第二条不等式实际上是电网功率平衡的变体，它表达了一个直观逻辑：风电能用的最大量，取决于“用电负荷＋电锅炉额外吃掉的电－CHP必须发的电”。风电用不完，就是因为右边的可用空间被CHP的“以热定电”出力给挤占了。

电锅炉在这里扮演的角色很关键。它本质上是把风电转化成热量供给供热管网，相当于给风电“找个出口”。在约束表达里，电锅炉出力算作电负荷的一部分，同时计入热功率平衡方程：

• 电功率平衡：风电出力 + CHP电出力 = 电负荷 + 电锅炉耗电
• 热功率平衡：CHP热出力 + 电锅炉产热 + 储热罐放热 = 热负荷 + 储热罐蓄热

这个“电转热”（power-to-heat，P2H）的解耦思路，是北方供热地区提高风电消纳能力最成熟的手段之一。

2.3 CHP机组电热可行域的建模细节

CHP机组不是随便一个“电出力和热出力”组合都能运行。它有一个电热运行可行域，主要由以下几段构成：

• 背压工况线：热电比固定，电出力随热出力线性变化，通常是可行域的上边界；
• 纯凝工况线：热出力为零时的最小和最大电出力区间;
• 抽汽调节区间：电出力可以在一定范围内独立于热出力调节。

把这四条边界围出来的多边形区域用一组线性不等式约束表示，就是CHP机组的可行域。项目里我用的是最常见的四边形近似法：

P_chp_min <= P_chp(t) <= P_chp_max 0 <= H_chp(t) <= H_chp_max P_chp(t) - c1 * H_chp(t) >= P_min_pure_condensing P_chp(t) - c2 * H_chp(t) <= P_max_back_pressure

其中 c1、c2 是机组的电热特性系数，由厂家热力试验数据拟合得到。这个四边形近似在工程精度足够，但要注意：如果实际机组是抽凝式且运行区间包含多个工作模式，四边形会高估调节能力，需要做更精细的可行域聚类。这一点我后面在常见问题部分还会再提。

3. Matlab代码实现要点

3.1 工具箱选型：我为什么选YALMIP+Gurobi

Matlab平台做这类优化调度，可选方案有好几套：手写线性规划直接调linprog、用fmincon做非线性规划、用YALMIP建模后调Gurobi/CPLEX。这个项目由于约束多、变量维度大（每台机组、每个时段都有变量），linprog手动拼矩阵容易拼错且难维护，fmincon处理大规模线性约束又不够稳，最终我选的是YALMIP+Gurobi组合。

原因是这个模型本质上是混合整数线性规划（MILP）——机组启停状态的0-1变量和功率连续变量混在一起，Gurobi在求解速度上有碾压级优势。举个例子，24时段、3台CHP机组、1台电锅炉、1个储热罐的模型，变量数大约800个、约束约1500条，Gurobi默认参数下一般十几秒就能收敛到0.1%的MIP gap。用linprog手搓基本要在建模阶段耗一半时间，fmincon则可能因为非凸的可行域陷入局部最优。

当然，没有Gurobi license的朋友可以改用免费的CBC求解器，YALMIP底层可以直接调用。速度和稳定性略差，但中小规模场景完全够用。我建议初学先用CBC跑通流程，再换Gurobi提升性能。

3.2 数据结构的组织：从杂乱数据到规范矩阵

Matlab代码实现这个模型，最大难点不是求解公式，而是数据的组织。我强烈建议一开始就把所有输入数据封装成结构体（struct），不要散落着写在脚本里。项目里我的数据组织方式如下：

% 输入数据结构 sys.H = 24; % 调度时段数 sys.nt = 3; % CHP机组数 % 风电预测数据 wind.P_avail = [34 32 30 ...]; % 各时段可用风电出力，MW % 负荷数据 load.P_demand = [120 115 ...]; % 各时段电负荷，MW load.H_demand = [210 205 ...]; % 各时段热负荷，MW % CHP机组参数（3台机组分别定义） chp(1).P_max = 120; chp(1).P_min = 40; chp(1).H_max = 120; chp(1).c1 = 0.25; chp(1).c2 = 0.75; chp(1).ramp_rate = 20; % MW/h，爬坡速率 % 电锅炉参数 eb.P_max = 60; % 最大耗电功率，MW eb.eta = 0.98; % 电转热效率 % 储热罐参数 ts.V_max = 300; % 最大储热量，MWh ts.charge_rate_max = 30; % 最大充/放热功率，MW

把这些数据集中定义之后，模型构建阶段反复引用同一个结构体字段，不容易发生“变量名打错但语法合法”这类隐错。我见过太多同行把数据散落在一堆变量里，调试的时候根本理不清哪些属于同一台机组。

3.3 核心代码框架逐段拆解

模型构建的核心代码分五大块：变量定义、目标函数、等式约束、不等式约束、求解与结果回读。

第一步：定义决策变量

% 定义一个27x24的优化问题（3CHP+电锅炉+储热+风电+功率平衡相关） P_chp = sdpvar(3, 24); % CHP电出力 H_chp = sdpvar(3, 24); % CHP热出力 P_wind = sdpvar(1, 24); % 风电并网出力 P_eb = sdpvar(1, 24); % 电锅炉耗电功率 E_ts = sdpvar(1, 25); % 储热罐储热量（25个节点，覆盖首尾） H_ts_charge = sdpvar(1, 24);% 储热罐蓄热功率 H_ts_release = sdpvar(1,24);% 储热罐放热功率 U_chp = binvar(3, 24); % 机组启停状态

特别注意这里储热罐的储热量 E_ts 定义成了25个点而不是24个，这是为了表达储能状态在相邻时段之间的递推关系：E_ts(t+1) = E_ts(t) + 充热 - 放热。很多新手直接定义24个点，写递推约束时发现 t=24 时刻的状态没法闭环，后期还得回头改维度，不如一开始就预留一个节点。

第二步：目标函数

% 成本项 objective = 0; for t = 1:24 % 煤耗成本：线性近似，也可以取二次函数的线性化分段 objective = objective + sum(coal_cost_coef .* P_chp(:,t)); % 启停成本 objective = objective + sum(start_cost .* max(0, U_chp(:,t)-U_chp(:,t-1))); % 弃风惩罚成本（核心） objective = objective + K_penalty * (wind.P_avail(t) - P_wind(t)); end

第三步：约束条件组装

constraints = []; % 电功率平衡 for t = 1:24 constraints = [constraints, sum(P_chp(:,t)) + P_wind(t) == ... load.P_demand(t) + P_eb(t)]; end % 热功率平衡 for t = 1:24 constraints = [constraints, sum(H_chp(:,t)) + eb.eta * P_eb(t) ... + H_ts_release(t) - H_ts_charge(t) == load.H_demand(t)]; end % CHP可行域约束 for k = 1:3 for t = 1:24 constraints = [constraints, chp(k).P_min * U_chp(k,t) <= P_chp(k,t) <= chp(k).P_max * U_chp(k,t)]; constraints = [constraints, 0 <= H_chp(k,t) <= chp(k).H_max * U_chp(k,t)]; constraints = [constraints, P_chp(k,t) - chp(k).c1 * H_chp(k,t) >= chp(k).P_min]; constraints = [constraints, P_chp(k,t) - chp(k).c2 * H_chp(k,t) <= chp(k).P_max]; end end

第四步：求解与回读

ops = sdpsettings('solver', 'gurobi', 'verbose', 2, 'mipgap', 0.001); optimize(constraints, objective, ops); % 回读结果 P_chp_opt = value(P_chp); H_chp_opt = value(H_chp); P_wind_opt = value(P_wind); P_eb_opt = value(P_eb); E_ts_opt = value(E_ts);

3.4 爬坡约束和储热罐约束：两类容易漏掉的物理约束

除了上述基础约束，还必须给机组加爬坡约束。CHP机组有功出力不能瞬间跳变，通常用：

for k = 1:3 for t = 2:24 constraints = [constraints, ... -chp(k).ramp_rate <= P_chp(k,t) - P_chp(k,t-1) <= chp(k).ramp_rate]; end end

这里我吃过亏：不加爬坡约束的模型优化结果很“漂亮”，风电消纳率高得惊人，但拿到现场一执行，机组响应不了那么快，计划全部作废。加上爬坡约束后消纳率会下降两三个百分点，这才是真实可信的结果。

储热罐方面有两组约束：一是储热量上下限和初始/末尾能量约束（防止调度周期结束把热全部放空），二是充放热功率互补约束（同一时段不能同时充和放）。

% 储热罐能量递推 for t = 1:24 constraints = [constraints, E_ts(t+1) == E_ts(t) + ... (H_ts_charge(t) - H_ts_release(t))]; end % 储热罐容量约束 constraints = [constraints, E_ts(1) == ts.E_initial]; constraints = [constraints, E_ts(25) == ts.E_final]; constraints = [constraints, 0 <= E_ts(2:25) <= ts.V_max]; % 充放热不能同时进行（用二进制变量互斥） C_ts = binvar(1, 24); constraints = [constraints, H_ts_charge(t) <= ts.charge_rate_max * C_ts(t)]; constraints = [constraints, H_ts_release(t) <= ts.charge_rate_max * (1-C_ts(t))];

末能量约束这一条特别重要。如果把 E_ts(25) 放开不管，求解器会把所有热都留在罐子里（等于白给的成本），或者把储热量压到零（相当于调度吃光了所有储能余量），两种都不是运营商想要的。通常做法是让始末储热量相等，保证日循环可持续性。

4. 案例仿真与结果分析

4.1 典型场景设置

我用一个基础算例来验证模型。场景设定如下：

• 3台CHP机组，容量分别为120MW、100MW、80MW；
• 电锅炉1台，最大耗电60MW，效率0.98；
• 储热罐有效容量300MWh，最大充放热功率30MW，初始储热量100MWh；
• 调度周期24小时；
• 风电预测数据选用北方某风电场冬季典型日，夜间23:00-07:00风电大发，白天出力回落；
• 热负荷全天高位运行，尤其夜间更大；电负荷白天高、夜间低。

这个场景基本复现了“风电大发时段恰好供热需求高、电负荷低”的典型困境。

4.2 没有联合优化时发生了什么

为了说明效果，先跑一版约束“CHP按以热定电方式运行、无电锅炉无储热罐”的基准场景。结果一目了然：夜间23:00至次日07:00，CHP为了带供热，总电出力最低仍有160MW左右，而电负荷只有120MW，硬生生占掉了风电的上网空间。风电可用功率在夜间峰值达到100MW，实际并网只有30MW左右，弃风率高达60%以上。这数据看着触目惊心，但在北方供热期是真实存在的普遍现象。

4.3 引入联合优化后的效果

加入电锅炉和储热罐、并运行联合优化模型后，结果发生了明显变化：

• 风电消纳比例从约52%提升到了约90%以上；
• 弃风集中时段明显缩短，从原来夜间连续7小时弃风缩减到凌晨两三个小时的轻度弃风；
• CHP总电出力曲线不再是完全跟着热负荷走，而是在风电大发的时段主动压低出力，把供热任务部分转移给电锅炉和储热罐放热；
• 系统总煤耗成本比基准场景有所上升（因为电锅炉耗电增加），但加上弃风惩罚成本后的“综合成本”明显下降。

我额外做了一个 CHP+储热罐（无电锅炉）的对照算例，发现储热罐单独使用只能把风电消纳率提高到70%左右。原因在于储热罐只能平移热负荷的时序，不能真正扩大用电空间；电锅炉则直接把风电转化为热负荷，相当于“新增了用电需求”，两者配合才能最大化消纳效果。这个结论和业内多篇文献的结论一致，也验证了模型的合理性。

提示：实际工程中不要指望消纳率能到100%。受机组最小技术出力、网络约束等因素影响，轻微弃风往往是经济上最优的。模型的价值是把弃风控制到合理区间，而不是数字上做到最好看。

4.4 不同惩罚系数下的敏感性分析

消纳效果和成本之间有一个平衡点，而这个平衡点主要由弃风惩罚系数 K_penalty 控制。我跑了一组敏感性分析：

K_penalty (元/MWh)	风电消纳率	系统总成本（万元/日）	机组启停次数
200	72%	156	6
500	86%	162	8
800	92%	164	12
1200	94%	166	18

注意看，惩罚系数从800提高到1200，消纳率只提高2个百分点，成本却多出2万元，机组启停次数也大幅增加——频繁启停对设备寿命非常不利。也就是说，盲目追求“零弃风”并不理性。实际工程里我一般把这个系数定在600~900元/MWh之间，让模型在“尽量消纳”和“避免过度调节”之间取一个平衡。

5. 常见问题与调试心得

5.1 求解器报“数值鲁棒性”问题怎么办

YALMIP/Gurobi求解时经常遇到提示：“Numerical problems / Ill-conditioned data”。这通常不是模型写错了，而是量纲差异太大。比如储热罐蓄热量单位是MWh（数值在几百），而弃风成本单位是元（数值在几万甚至几十万），目标函数里同时出现不同量级的数值，数值求解器容易出状况。

我的处理方式是统一量纲：所有功率量用MW，能量量用MWh，费用项全部除以1e3变成“千元”；或者干脆把目标函数整体缩小1000倍，保持求解数值在1e-3到1e5区间内。这类细节在大规模模型中影响更大，数值扰动会导致收敛速度骤降甚至报错。

还有个细节，储热罐的储热量维度（MWh）和充放热功率维度（MW）在递推约束里相乘要匹配：E_ts(t+1) = E_ts(t) + [充放热功率 × 1小时]，这里每个时段长度是1小时，所以功率数值可以直接累加。如果改成15分钟一个时段，递推式里必须乘以0.25，否则能量不守恒的bug非常隐蔽，查半天查不出来。

5.2 风电消纳率算出来明显偏高？检查约束有没有漏

我发现一个很常见的排查路径：模型刚建完跑出来风电几乎全额消纳，高到不真实。这时候优先排查三类约束：

• 是否遗漏了CHP机组出力下限和可行域下边界约束？如果 P_chp 的边界写成了 0 <= P_chp <= P_max，那么求解器可以把 CHP 出力压到零来腾空间给风电，这在物理上完全不合理；
• 是否遗漏了爬坡约束？夜间风电大发时，求解器可能让CHP从晚高峰的高出力直接跳到夜间的低出力，一个时段内跨了太大幅度；
• 是否给电锅炉设置了过大的容量？如果电锅炉容量超过热负荷需求，模型会无脑把电锅炉拉到满发来消耗风电，热平衡被“强行”满足，这也不是真实可行方案。

这三条我全踩过。排查手段也不复杂：把优化得到的各变量曲线打印出来，逐一核对每个时段功率平衡是否严格满足、设备出力是否落在可行域内，半小时就能定位问题。

5.3 整数变量导致的求解时间剧增

模型里有机组启停状态和储热罐互斥状态两类二进制变量，如果再加电锅炉的开关状态，总整数变量会达到“机组数×时段数×2~3”个。以10台机组、24时段算，就有至少480个整数变量，直接跑默认参数可能要好几分钟。

优化手段有三招：

• 第一，把非必要整数变量替换成连续变量+边界约束。比如电锅炉如果运行模式是“深度调峰期间才开”，可以预判时段，提前固定开关状态，减少整数变量；
• 第二，给Gurobi设置合理的MIP gap，比如 0.001 或 0.005。对调度计划而言，1%以内的次优解完全可接受，没必要追求数学上的绝对最优；
• 第三，用“滚动优化”替代全局24小时优化，例如每4小时滚动一次，每次只优化未来4小时并冻结后20小时的预测信息，模型规模大幅下降，求解速度从分钟级降到秒级。

5.4 CHP可行域多边形过简化的后果

最后专门提一个模型精度问题。前面用四边形近似CHP可行域，其实隐性地“缩减”了机组实际调节空间。在实际项目中，我遇到过抽汽式CHP机组的实际可行域接近五边形甚至六边形（考虑低压缸最小冷却流量约束），四边形近似可能导致：

• 计算结果偏向保守，系统调节能力被低估，风电消纳率低于实际可达到水平；
• 极端情况下优化的运行点在四边形内但不在真实可行域内，下发给现场后机组无法跟踪执行。

解决方法是把机组厂家的热力特性试验数据找出来，用凸包法（convhull）拟合真实可行域的多边形顶点，然后转化成多面体约束：

% 用凸包顶点构建可行域约束（YALMIP支持polygon约束） Poly = [P1 H1; P2 H2; P3 H3; ...]; % 从试验数据得到的顶点 constraints = [constraints, ismember([P_chp(k,t); H_chp(k,t)], Poly)];

如果没有试验数据，至少也要把四边形的关键顶点坐标标定准确，不要随便拍脑袋定系数。C1、C2这类系数偏差10%，最终调度计划就能差出一个档次。

6. 一些值得坚持的工程习惯

项目做到这个程度，模板代码其实已经完成核心验证。如果后续要往真实现场走，我建议再补齐三个动作：

• 把优化结果做24小时功率曲线可视化分析，每一张图都要能解释“为什么在这个时刻CHP压低出力、储热罐放热”；
• 做一个“风电预测误差”影响分析，在风电预测值±20%扰动下重复求解，观察消纳率的波动范围；
• 预留与SCADA/AGC系统的数据接口，把优化结果输出成标准格式的调度指令文件。

我个人在实际操作中有一个习惯，值得分享：模型跑通之后第一件事不是调参数刷高消纳率，而是先把“基准场景”的每一条约束都对一遍物理意义，再放开对比。Matlab里用YALMIP建模的最大便利是约束可读性高，逐条打印出来就能当报告用。这套代码我后来又扩充过多区域版本、加了日前实时两阶段调度逻辑，核心框架基本没动过。前期把模型和数据结构写规范，后面改起来会顺手很多。