企业官网建设流程全解析

我最早被热电联产机组和风电消纳绑在一起“折磨”，是在一次供热期的调度方案复盘上。当时风电出力并不低，系统却必须压掉一部分风电场出力，原因不是线路堵，而是热电厂为了保证居民供暖，机组电出力降不下来。这个现象在业内人士口中叫“以热定电”，它直接导致冬季夜间风电越高、弃风越猛。后来我做了大半年联合优化控制方向的工作，把热电联产机组、蓄热罐、电锅炉放进同一个模型里协调，才真正把弃风率压下去。这篇博文就把整套思路和Matlab代码实现完整拆开讲，适合正在做电力系统优化调度、热电联产灵活性改造、新能源消纳方向课题的研究生和工程师参考。

1. 为什么供热季的风电会被“挤”出去——热电联产的刚性约束

1.1 供热期“以热定电”是如何压缩风电空间的

热电联产机组和普通纯凝火电机组最大的区别在于，它发电的同时还要对外供热。根据机组形式不同，热和电的耦合方式有两种：

• 背压式机组：排汽全部用于供热，电出力与热出力完全刚性绑定，热发多少就产生多少电，基本没有独立调节能力。
• 抽汽式机组：从汽轮机中间级抽一部分蒸汽去供热，剩余蒸汽继续发电。这类机组可以在一定范围内调整抽汽量，因此电出力有一个围绕热负荷变化的可行域，但整体上仍然是“热需求越大，最小电出力越高”。

在供热期，系统调度员面对的第一约束就是热负荷必须满足。只要热负荷降不下来，所有带供热的机组就有一个无法突破的最低电出力。风电消纳空间可以简单写成：

风电消纳空间 = 系统总负荷 - 常规机组最小技术出力 - 热电联产机组最小电出力 - 联络线净受入功率

这个公式看起来普通，却是理解一切问题的基础。冬季夜间是典型低谷负荷时段，同时热负荷达到全天最高，此时热电联产机组因为供热被“压”在很高的电出力水平上，留给风电的调节空间就很小。风电场出力越高，超出空间的那部分就只能被弃掉。所以会出现一个反直觉的现象：风越大，弃风越严重。

1.2 联合优化控制到底“联合”了什么

单台热电联产机组的“以热定电”约束是物理规律，没法直接消除。但把多台机组和多种热源放到一起看，情况就不一样了。所谓联合优化控制，是把以下对象统一纳入一个优化模型，在满足热负荷和电负荷的前提下，找到成本最低、弃风最小的运行方式：

• 多台热电联产机组之间的电热出力分配
• 蓄热罐的充放热计划
• 电锅炉等电转热设备的启停与功率
• 常规纯凝机组的出力
• 风电场出力与弃风量
• 联络线交换功率（如果涉及区域级调度）

本质上是把“热”从“电”的刚性束缚中部分解放出来。比如白天电负荷高、热负荷低时，让热电联产机组多发电并多抽汽供热，把热量储存在蓄热罐里；夜间风电出力高、电负荷低时，用蓄热罐放热替代一部分机组供热，让机组电出力降下去，把电网空间让给风电。再加上电锅炉直接把风电转为热能，整个系统的风电接纳能力会有明显改善。

1.3 灵活性手段的建模价值

做优化控制，不能只讲“改造”，要讲“可控性”。蓄热罐、电锅炉这类灵活性资源，价值全在于它们给模型增加了可调的决策变量和储能动态约束：

• 蓄热罐相当于给热负荷增加了一个“时间平移缓冲区”，约束表现为相邻时段的储热量递推关系。
• 电锅炉相当于给电负荷增加了一个“可调节耗电设备”，约束表现为功率上下限与爬坡限制。
• 电转热过程相当于把风电消纳空间转换为热能存储，提高了系统对风电预测误差的包容度。

有了这些变量，联合优化控制就从“调一台机组”升级成了“调一个系统”。我不会在模型里幻想彻底消除以热定电，而是通过协调手段把刚性约束的影响最小化。这也是下文建模部分所有数学表达的核心出发点。

2. 联合优化控制模型的搭建：目标函数与约束怎么落笔

2.1 目标函数怎么写：经济成本与弃风惩罚

优化控制首先要明确“优化什么”。对这个题目，最自然的目标是最小化系统运行总成本，同时在目标函数里加入弃风惩罚项。为什么不能直接把目标设成“最大化风电消纳量”？因为电负荷、热负荷、机组爬坡、热网时间延迟都有限制，强行消纳可能让某些机组运行在极端低负荷区间，导致煤耗和损耗上升，反而得不偿失。

我使用的目标函数形式为：

min sum( C_fuel * F(P_chp, H_chp) ) + sum( C_start * start_flag ) + sum( C_eb * P_eb ) + sum( lambda_wind * P_wind_curtail )

其中：

• F(P_chp, H_chp)是热电联产机组的燃料消耗函数，本文用线性或分段线性近似。
• C_fuel是燃料价格，单位为元/吨标准煤。
• C_start是机组启动成本，单位元/次。
• C_eb是电锅炉用能成本系数。
• lambda_wind是弃风惩罚系数，单位元/MWh。
• P_wind_curtail是弃风功率，也就是预测可发风电减去实际并网风电的部分。

弃风惩罚系数必须设置得比常规发电成本高，否则优化结果会出现“经济上合理、物理上浪费”的弃风。我习惯取煤耗成本的2到3倍，具体要看项目对弃风率的考核权重。

时间尺度上，通常采用日前计划或者日内滚动优化，时间间隔取1小时，优化窗口为24小时。更短的时间尺度（比如15分钟）可以在滚动预测环节使用，但建模框架完全相同，只需把热网延迟处理得更细。

2.2 约束条件逐条拆解

约束是模型能不能跑出可信结果的关键。我把约束分为六类，每一类都对应实际物理系统中的一个刚性条件：

电功率平衡

sum(P_chp) + sum(P_thermal) + P_wind + sum(P_eb) = P_load + P_loss

其中P_load是系统电负荷，P_loss是网损，小规模算例可以忽略，工程上通常按负荷的一定百分比计入。

热功率平衡

sum(H_chp) + H_tank_dis + H_eb = H_load + H_tank_chg

蓄热罐的充放热在这个式子里体现为热负荷的时间平移。热量可以从机组侧“存”到罐里，也可以从罐里“放”到热网。

热电联产机组可行域约束

抽汽式机组的热电关系是一个凸包络区域，常用一组线性不等式描述：

P_chp(i,t) + alpha1 * H_chp(i,t) >= P_min_ref(i) P_chp(i,t) + alpha2 * H_chp(i,t) <= P_max_ref(i) P_chp(i,t) >= 0 0 <= H_chp(i,t) <= H_max(i)

这里的alpha1和alpha2是根据机组热特性曲线拟合得到的系数，代表背压工况和纯凝工况之间的转换斜率。工程上如果直接用一个线性关系刻画，误差会很大，建议至少用两个斜率分段逼近。

机组爬坡约束

-delta_down * dt <= P_chp(i,t) - P_chp(i,t-1) <= delta_up * dt

这个约束防止模型给出物理上无法实现的“瞬跳”。我见过不少初学模型因为漏了爬坡约束，结果最优解的出力曲线在相邻时段直接跳了几十兆瓦，调度员根本执行不了。

蓄热罐动态约束

E_tank(t+1) = E_tank(t) + eta_chg * H_tank_chg(t) - H_tank_dis(t) / eta_dis 0 <= E_tank(t) <= E_tank_max

蓄热罐不是一个简单的“可以随便用”的电源，它有容量上下限和充放热效率。eta_chg和eta_dis通常取0.9到0.98之间；数值上看着小，但长时间尺度下累加效应非常明显。

弃风与风电出力约束

P_wind(t) = P_wind_forecast(t) - P_wind_curtail(t) 0 <= P_wind(t) <= P_wind_forecast(t) P_wind_curtail(t) >= 0

弃风功率在这里是一个松弛变量，它的存在保证模型在所有风电预测偏差下都可解。最终目标函数里的惩罚项会驱动求解器把这个变量压到尽可能小。

2.3 决策变量与参数表

为了方便Matlab建模代码与后文对应，把模型中的符号统一如下。

参数部分，每个机组的P_min、P_max、H_max、爬坡率、煤耗斜率都需要从机组出厂热力试验报告或历史运行数据中拟合出来。如果用公开测算数据，一定要注明数据来源或至少说清楚是算例假设，不能直接当真实工程数据用。

3. Matlab求解实现：从数据准备到YALMIP建模与求解

3.1 建模工具和求解器选择

Matlab下做优化建模，我优先推荐YALMIP加Gurobi的组合。YALMIP最大的优势是把变量定义、约束叠加、目标函数表达三件事句法化，写起来和论文公式几乎一一对应，排错快。求解器方面，Gurobi或CPLEX对混合整数线性规划的处理能力很强，但如果你的环境没有商业求解器许可，intlinprog同样能完成这个任务。

我在实际项目中两者的用法差异不大：先在YALMIP里用sdpvar定义变量，用optimize()调用求解器，唯一区别是optimize()里的求解器参数。如果模型规模不大，比如机组数量不超过5台、时间窗口不超过48小时，intlinprog完全够用；规模再大，才需要考虑Gurobi。

3.2 代码结构怎么搭

代码不追求花哨，重点是让模型逻辑清晰。我的工程目录一般这样组织：

case_data.m # 所有算例参数：负荷、热负荷、风电、机组参数 build_model.m # 定义变量、目标函数、约束 solve_case.m # 主脚本，调用build_model并求解 plot_result.m # 绘制出力曲线、热平衡曲线、弃风率

case_data.m返回一个结构体data，字段名和模型符号保持对应，这样从参数到代码基本不需要查表翻译。

3.3 核心建模代码示例

下面给出一个精简但可运行的核心片段，对应2.2节的约束框架。这个片段省略了部分边界细节，但结构完整，可以直接扩展。

%% solve_case.m 主脚本片段 data = case_data(); I = data.I; % 热电联产机组台数 T = data.T; % 时间窗口长度 dt = data.dt; % 时间间隔，小时 %% 定义决策变量 P_chp = sdpvar(I, T, 'full'); H_chp = sdpvar(I, T, 'full'); P_th = sdpvar(1, T, 'full'); P_wind = sdpvar(1, T, 'full'); P_curtail = sdpvar(1, T, 'full'); H_tank_chg = sdpvar(1, T, 'full'); H_tank_dis = sdpvar(1, T, 'full'); E_tank = sdpvar(1, T+1, 'full'); P_eb = sdpvar(1, T, 'full'); u_chp = binvar(I, T); % 机组启停状态，0-1变量 %% 目标函数 obj = 0; for t = 1:T for i = 1:I % 煤耗成本，线性化近似 obj = obj + data.fuel_cost(i) * (data.a(i) * P_chp(i,t) + data.b(i) * H_chp(i,t)); obj = obj + data.start_cost(i) * max(0, u_chp(i,t) - ... (t == 1 ? 0 : u_chp(i,t-1))); end obj = obj + data.c_eb * P_eb(t); obj = obj + data.lambda_wind * P_curtail(t); end %% 约束 Constraints = []; % 电功率平衡 for t = 1:T Constraints = [Constraints, sum(P_chp(:,t)) + P_th(t) + P_wind(t) + P_eb(t) == data.P_load(t)]; end % 热功率平衡 for t = 1:T Constraints = [Constraints, sum(H_chp(:,t)) + H_tank_dis(t) + data.eta_eb * P_eb(t) == data.H_load(t) + H_tank_chg(t)]; end % 热电联产可行域：以两段线性包络为例 for i = 1:I for t = 1:T Constraints = [Constraints, P_chp(i,t) + data.alpha1(i) * H_chp(i,t) >= data.P_min(i) * u_chp(i,t)]; Constraints = [Constraints, P_chp(i,t) + data.alpha2(i) * H_chp(i,t) <= data.P_max(i) * u_chp(i,t)]; Constraints = [Constraints, 0 <= H_chp(i,t) <= data.H_max(i) * u_chp(i,t)]; Constraints = [Constraints, P_chp(i,t) >= 0]; end end % 爬坡约束 for i = 1:I for t = 2:T Constraints = [Constraints, P_chp(i,t) - P_chp(i,t-1) <= data.delta_up(i) * dt]; Constraints = [Constraints, P_chp(i,t-1) - P_chp(i,t) <= data.delta_down(i) * dt]; end end % 蓄热罐动态 Constraints = [Constraints, E_tank(1) == data.E_tank_init]; for t = 1:T Constraints = [Constraints, E_tank(t+1) == E_tank(t) + data.eta_chg * H_tank_chg(t) - H_tank_dis(t) / data.eta_dis]; Constraints = [Constraints, 0 <= E_tank(t) <= data.E_tank_max]; Constraints = [Constraints, 0 <= H_tank_chg(t) <= data.H_tank_max]; Constraints = [Constraints, 0 <= H_tank_dis(t) <= data.H_tank_max]; end % 风电约束 for t = 1:T Constraints = [Constraints, P_wind(t) + P_curtail(t) == data.P_wind_forecast(t)]; Constraints = [Constraints, 0 <= P_wind(t) <= data.P_wind_forecast(t)]; Constraints = [Constraints, P_curtail(t) >= 0]; end %% 求解 options = sdpsettings('solver', 'gurobi', 'verbose', 1, 'dualize', 0); sol = optimize(Constraints, obj, options); if sol.problem == 0 fprintf('求解成功，目标值为 %.2f\n', value(obj)); else disp('求解失败，请检查约束和求解器日志'); end

这段代码需要配合case_data.m才能跑通，但核心信息都在这了。需要注意的一点是，我在约束里用了u_chp这个0-1变量，这意味着我们求解的是混合整数线性规划。对于仅做研究的场景，也可以去掉启停变量，把所有机组的u_chp固定为1，模型退化为线性规划，求解速度会快很多。

3.4 从求解结果到调度指令

求解完成后，把value()的内容取出来，写入结构体，方便后续绘图和分析：

result.P_chp = value(P_chp); result.H_chp = value(H_chp); result.P_wind = value(P_wind); result.P_curtail = value(P_curtail); result.E_tank = value(E_tank); result.H_tank_chg = value(H_tank_chg); result.H_tank_dis = value(H_tank_dis); result.P_eb = value(P_eb);

我习惯把弃风率定义为：

弃风率 = sum(P_curtail) / sum(P_wind_forecast) * 100%

然后在主脚本末尾直接打印。如果优化结果里弃风率仍然很高，先检查是不是约束给得太紧，比如热平衡里的H_load是否被硬性卡死，或者蓄热罐容量是不是设成了0。

4. 一个典型供热日场景的仿真结果

4.1 算例设置

为了验证模型效果，我设置了一个典型的供热日场景。包含3台热电联产机组、1台纯凝机组、1个风电场、1台电锅炉和1个蓄热罐。具体参数如下。

热负荷曲线按“早晚高、午间低”设置，电负荷曲线按“午间高、夜间低”设置，风电预测曲线按“夜间高、白天低”设置。这个组合对应了北方供热期最典型的“风热矛盾”场景。

基准方案采用传统的“以热定电”调度：热电联产机组按满足热负荷的最低电出力运行，蓄热罐和电锅炉不参与优化。对比方案则使用本文的联合优化控制模型，所有设备都在优化窗口内自由协调。

4.2 结果对比

联合优化方案以煤耗增加约90吨标准煤为代价，把弃风率从18.6%压到了5.2%，多消纳了151MWh风电。从成本角度看，弃风惩罚项在大幅下降，虽然煤耗成本上升，但总目标值仍然更优。

这个结果符合预期：风电多发的代价不是零，而是以一部分煤耗成本去换风电消纳空间。如果目标函数里不设置弃风惩罚，优化器很可能为了省煤而不去消纳那部分低价值风电。所以我在实际项目里非常看重惩罚系数的敏感性分析。

4.3 出力曲线的变化逻辑

观察最优调度曲线，最明显的特征是夜间时段热电联产机组电出力被明显压低。这主要来自三个机制共同作用：

• 蓄热罐在午后和傍晚蓄热，夜间放热，替代了部分机组供热。
• 电锅炉在夜间风电峰值时段启动，把富余风电直接转化为热能。
• 热电联产机组在白天电负荷较高时多发电、多蓄热，晚上则降低电出力。

如果没有蓄热罐和电锅炉，仅仅依靠机组间的分配优化，弃风率大概只能从18.6%降到12%左右。也就是说，灵活性资源的引入是联合优化控制里最关键的变量，机组间的“小优化”和跨时段的“大优化”效果差距明显。

我在结果分析里习惯额外输出一个“机组最小出力包络线”图，把所有CHP机组在每个时段允许的最小总电出力画出来。这条包络线在联合优化方案里会比基准方案低一大截，风电空间的变化一眼就能看出来。绘图用常规的plot和area就能完成，不需要额外工具箱。

5. 工程化落地中的几个坑和经验补充

5.1 热负荷数据不能直接用“给定曲线”

很多初版模型把H_load当成常数序列直接输入，这是最大的坑之一。实际热网有热惯性，热媒从热源到用户侧有数小时延迟，供热管道本身也有储热能力。用瞬时热负荷建模，会让优化结果高估蓄热罐的调节能力，甚至出现“热源出力与用户需求在时间上错位”的假最优。

我在模型里补充的处理方法是加入简化热网延迟模型：把热负荷按时间序列平移若干时段，或者用一个一阶惯性环节表示热网储能。如果项目阶段还处理不了热网模型，至少要在结果分析中给热平衡加一个允许的误差带，不能把热平衡约束卡得太死。

5.2 热电联产可行域线性化必须保证凸性

机组热力特性曲线通常是非线性的，实际可行域要取凸包络。如果直接用离散点连线，很容易得到凹多边形，求解器在凹可行域上做线性规划，最优解可能落在物理不可行的点上。这个问题在调试时特别隐蔽，因为目标函数值差距不大，但出力点实际无法实现。

我的建议是：拿到每个工况点的实际电出力、热出力和煤耗后，先画出热电关系散点图，然后人工判断边界形状。如果要做多段线性化，每段的斜率必须使整个包络保持凸性。YALMIP对线性约束组成的集合不强制凸性检查，所以这个把关必须自己做。

5.3 风电预测误差怎么兜底

日前优化用的是预测曲线，日内实际风速往往偏差不小。如果模型只有一个确定性的风电序列，实际执行时要么弃风超标，要么功率不足导致电平衡被破坏。兜底手段有三层：

• 在电平衡约束中加入旋转备用约束，比如sum(P_max) - sum(P_chp) >= reserve_up(t)。
• 把弃风惩罚系数调到一个合理水平，让模型在没有完全把握时不至于冒险抢占过多风电空间。
• 日内滚动更新：每1小时或每4小时用最新风电预测重算一次优化，只执行第一步的动作。

第三层是工程上最有效的手段。滚动优化对计算时间的要求更高，但Matlab环境下配合Gurobi，T=24的模型单次求解通常在几秒内完成，完全可以支撑小时级滚动。

5.4 Matlab求解性能优化的小经验

模型规模变大后，求解时间会明显上升。我调优的顺序是先看整数变量数量。很多时候机组启停变量并不是必需的，尤其是调度周期内所有CHP机组都处于连续运行状态时，直接把u_chp固定为1，问题从混合整数线性规划退化为纯线性规划，求解时间可能缩短一个数量级。

如果启停变量确实需要保留，可以给求解器设置一个合理的MIP gap，例如2%。Gurobi的MIPgap参数设为0.02，通常能在牺牲极少精度的情况下大幅缩短求解时间。还有一个容易被忽略的点：sdpsettings里把solver指定为具体求解器，比让YALMIP自动选择更稳定，否则遇到复杂的约束组合时，YALMIP可能选到一个处理不了的问题类型。

5.5 调试时一定要打印的几个诊断量

最后分享一个我自己总结的调试检查清单。每次优化跑完，除了看弃风率，我固定检查以下五个量：

• 热平衡松弛量：每个时段的sum(H_chp) + H_tank_dis + eta_eb*P_eb - H_load - H_tank_chg，应该严格为0。
• 蓄热罐SOC终值：如果结束时刻罐内还有大量热量，且第二天初始值固定为低位，模型可能在第一天“故意”多蓄热来规避惩罚，需要加入末状态约束。
• 每台CHP的出力落点：看是否在可行域内部，而不是极限边界。
• 弃风惩罚项是否在目标函数中起主导作用：如果惩罚项占了总目标50%以上，说明弃风惩罚系数可能过高。
• 机组爬坡约束是否被激活：如果某些时段反复被激活，说明热负荷或风电变化太陡，模型在硬顶着限制走。

这套检查清单帮我发现过好几次模型“看起来收敛、实际矛盾”的情况。最经典的一次是蓄热罐SOC约束写反了方向，导致夜间放热逻辑变成夜间蓄热，优化器把热电联产机组电出力抬得更高，弃风率不降反升。不打印SOC逐时曲线根本发现不了。

这套联合优化控制模型我已经在两个不同规模的项目里复用过，大框架基本不变，变的主要是热网延迟的精细程度和约束的松紧尺度。对你手头的具体问题，建议先把基准算例跑通，确认热平衡、电平衡、蓄热罐动态这三个核心约束没有矛盾，再逐步加复杂度。等到弃风率数字开始随参数变化而规律性变化时，模型基本就靠谱了。