Matlab实现电热综合能源系统日前经济调度模型

做电热综合能源系统日前经济调度研究时，我一直在琢磨一件事：怎么把可再生能源消纳的压力体现在优化模型里？这套用Matlab代码实现的电热综合能源系统日前经济调度模型，就是在这件事上折腾了几个月的结果。它不是那种只能跑通文档的玩具demo，而是能从负荷数据、装机参数一路算到24小时出力计划、弃风弃光指标的完整闭环。对于正在做综合能源系统优化、电力系统经济调度的研究生或者工程师，这套思路和代码可以直接拿来改、拿来用。

我先把话挑明：这个模型不是骡子式的“电/热两套系统各算各的”，而是把电锅炉、储热罐、热电联产机组放在同一个优化框架里，以24小时为周期做经济调度。目标是花钱最少，同时让风电、光伏尽量多上网。下面我把从模型设计到Matlab实现、再到调参踩坑的过程完整拆开讲。

1. 模型思路：电热耦合不是简单加个约束

1.1 为什么必须把“电”和“热”放在同一个优化框架里

如果你单独做电力调度，你会忽视供热机组的热出力约束；单做热力调度，你又看不到电网的实时电价信号。现实中的园区能源系统，最常见的场景是这样：冬季采暖季，热电联产机组（CHP）为了满足热负荷，锅炉必须维持较高出力，发电功率被“绑着”往上走。而这时候如果恰好是夜间风电大发，电网那侧消化不了这么多电，风机只能弃掉。这就是典型的“以热定电”导致弃风。

解决办法不是关掉热负荷，而是引入电锅炉、储热罐这些“解耦”设备。电锅炉让电力需求在供热侧多一个出口，储热罐则把多余热量存起来，等电价高的时候再放出来。这样一来，CHP机组的热出力不必紧跟热负荷，电出力也就有了下调空间。整个系统的运行灵活性上来了，可再生能源消纳的通道自然就打开了。所以，电热综合调度不是拍脑袋“为了综合而综合”，是电网灵活性不足倒逼出来的需求。

我这个模型把电、热两侧的平衡方程放在同一个线性规划里，调度对象包括常规发电机组、CHP机组、风电、光伏、电锅炉和储热装置。目标是在满足电热负荷的前提下，最小化总运行成本，同时用弃风弃光惩罚项引导模型优先消纳新能源。这种处理方式最直接，也最容易在论文里讲清楚。

1.2 日前经济调度到底在算什么事

日前调度，说白了就是提前一天决定每台设备在未来24小时每个时段的出力计划。我遇到不少初学者分不清“经济调度”和“机组组合”：经济调度假设机组启停状态已经给定，只优化各台设备的出力大小；而机组组合还要额外决定哪些机组在哪些时段开机、哪些停机，中间会引入0/1整数变量。我选择先做经济调度，是因为它能把问题本质讲清楚——所有约束都是线性的，目标函数也是线性的，调参、调试、可视化都非常顺手。等你把经济调度跑通，再往上加0/1变量做机组组合、加不确定场景做随机优化，就是在同一个骨架上扩展。

模型的时间分辨率取1小时，全天24个点。这个粒度对研究“消纳”问题足够了：风电低谷波峰、电价低谷时段都能捕捉到。如果你要做分钟级调度，模型规模会变大，但约束本质完全一样。我没在模型里考虑电网络潮流，默认整个园区是一个节点，所有设备挂在同一母线上。这样做的好处是省去节点电压、线路潮流这些非线性项，把注意力集中在“电-热耦合”和“可再生能源消纳”这两个核心问题上。等需要精细化网架分析，再叠加潮流约束也不迟。

2. 目标函数与约束条件建模：每个公式都对应一句人话

2.1 目标函数怎么写才合理

我的目标函数分四块：外购电成本、机组燃料成本、弃风弃光惩罚成本、设备运维成本。外购电是从上级电网买电的费用，用分时电价乘以购电功率逐时求和。机组燃料成本包括燃气轮机和CHP消耗天然气的费用，可以近似表示成出力的线性函数。经济调度研究里常用二次曲线，但实际参数标定很麻烦，线性近似在方案对比阶段完全够用。设备运维成本很简单，每台设备单位出力的维护单价乘上出力即可，不是影响结果的核心项。

最关键的是弃风弃光惩罚成本。风电光伏的预测出力是当天模型能消纳的上限，实际出力是决策变量，预测值和实际值之差就是弃电量。惩罚项的形式是：

C_pen * (P_wind_pred(t) - P_wind(t) + P_pv_pred(t) - P_pv(t))

这个惩罚系数要设得够大，比如500元/兆瓦时，这样模型“舍不得”弃风。如果系数设成10元/兆瓦时，而电价低谷时购电成本才0.3元/度，模型宁可让风机少发然后去电网买电，导致弃风率虚高，这显然不符合“消纳优先”的预期。

完整目标函数写成文字公式就是：

min Σt [ c_el(t)*P_buy(t) + c_gas*P_gt(t) + c_chp_el*P_chp(t) + c_chp_h*H_chp(t) + c_pen*(P_wind_pred(t)-P_wind(t)+P_pv_pred(t)-P_pv(t)) + c_om_eb*P_eb(t) ]

这里P_gt是燃气轮机电出力，P_chp、H_chp分别是CHP电、热出力，P_buy是购电功率，P_eb是电锅炉耗电功率。把这个目标写进Yalmip时，注意向量化和系数匹配：如果这些决策变量都是1×24的sdpvar向量，那么目标函数用sum(系数 .* 变量)这种形式，一行代码搞定，不需要写for循环。

2.2 电功率平衡和热功率平衡：能量守恒是硬约束

任何调度模型，最基本的约束就是功率平衡。电功率平衡如下：

P_gt(t) + P_chp(t) + P_wind(t) + P_pv(t) + P_buy(t) - P_eb(t) = P_load_e(t)

其中P_load_e是电负荷。这里P_buy是非负变量，假设系统只购电、不向电网倒送电。有些园区允许余电上网，那就额外加一个P_sell变量，等式右侧再减去它。热功率平衡则是：

H_chp(t) + H_eb(t) + H_dis(t) - H_char(t) = H_load(t)

H_eb是电锅炉产热功率，H_dis、H_char分别是储热罐的放热和充热功率。如果系统没有储热罐，这两个变量恒为0。

注意，热功率平衡里的H_eb和电平衡里的P_eb通过电锅炉效率耦合在一起：H_eb = eta_eb * P_eb，效率通常取0.9~0.95。电锅炉运行时会有热损耗，这部分损耗没有出现在热平衡里，而是体现在效率系数里。模型里不需要额外再扣一次热损失，别重复计入。

除了平衡约束，我还会给常规机组加爬坡约束：

-ramp_gt <= P_gt(t+1) - P_gt(t) <= ramp_gt

如果不加爬坡，模型会给出跳跃性的出力曲线，实际机组根本跟不上。但爬坡约束也不能设得太紧，否则可能直接导致无解，这一点后面调试部分会再谈。

2.3 热电联产机组的可行域：多边形约束怎么处理

CHP机组不是任意给定电出力和热出力都能运行。它的运行区间是一个多边形，由厂家给出的电-热运行范围决定。常见描述包括电出力上下限、热出力上下限、以及“电出力随热出力变化的斜率约束”。比如某台背压式CHP，电出力和热出力近似线性关系，可行域可以写成：

P_chp_min ≤ P_chp ≤ P_chp_maxH_chp_min ≤ H_chp ≤ H_chp_maxP_chp + k1 * H_chp ≥ C1P_chp + k2 * H_chp ≤ C2

具体系数由机组特性曲线标定。手动输入这些不等式很容易方向写反，我的做法是把多边形顶点坐标直接读进来，然后用凸包工具生成约束矩阵A和右侧向量b。Matlab里可以用convhull函数把顶点排序，再逐条边计算法向量，生成不等式。这样即使换一台不同参数的机组，只需要替换顶点矩阵就行，代码不用动。这个技巧在要多机组对比时特别管用。

2.4 电锅炉、储热罐等灵活性资源怎么建模

储热罐是模型里最有意思的部分。它的状态量是蓄热量S(t)，递推关系为：

S(t+1) = S(t) + eta_c * H_char(t) - H_dis(t) / eta_d - epsilon * S(t)

eta_c和eta_d分别是充热、放热效率，epsilon是散热损失系数。量纲要统一：S用MWh，H_char、H_dis用MW，时间步长是1小时，所以功率乘以1小时就是热量变化量。如果时间步长不是1小时，记得乘上对应间隔。

运行约束还包括蓄热量上下限、充放热功率上限：

0 ≤ S(t) ≤ S_max0 ≤ H_char(t) ≤ H_char_max0 ≤ H_dis(t) ≤ H_dis_max

还有一个日循环约束：为了场景能够重复执行（比如明天继续用同样策略），通常要求一个调度周期结束时的蓄热量等于初始值：

S(1) = S(25)（或S(end) = S(1)）

这个约束在夜间谷时段强制储热罐至少恢复到初始状态，否则模型会把储热罐里最后一点热量在峰时段全部放光，第二天无热可用。我一开始没加这个约束，结果是储热罐SOC在第24小时直接掉到下限，实际运行根本无法复现。加上周期约束后结果就正常多了。

3. Matlab实现：选对工具，代码能省一半时间

3.1 建模工具选型：Yalmip加求解器是黄金组合

Matlab里做优化模型的路子不少：手写矩阵调linprog、用CVX、用Gurobi自带接口，或者用Yalmip。我个人强推Yalmip。这玩意儿相当于一个“翻译官”，让你用接近数学表达式的语法定义变量和约束，然后自动翻译成求解器需要的标准形式。配合Gurobi或CPLEX，求解速度快到飞起。如果你暂时没有商业求解器，Yalmip默认调用Matlab的linprog，小规模LP也能跑，就是慢一些。我用的是Gurobi学术版，处理带大量稀疏矩阵的LP一直很稳，而且它在Yalmip里的配置就一行：sdpsettings('solver','gurobi')。

别用CVX做这件事。CVX对约束的拼接和指数锥问题的支持很优雅，但它的LP建模方式偏“规范型”，修改起来没有Yalmip自由。Yalmip的约束对象可以像搭积木一样往一个Constraints数组里塞，调试时还能用showconstr查看某条约束，体验舒服太多。

3.2 变量定义、约束组装的高频技巧

在Yalmip中，连续变量用sdpvar定义。除了矩阵还可以定义向量。假设系统有24时段，那么每个设备都定义成1×24的行向量：

P_wind = sdpvar(1, 24); % 风电实际出力 P_pv = sdpvar(1, 24); % 光伏实际出力 P_chp = sdpvar(1, 24); % CHP电出力 H_chp = sdpvar(1, 24); % CHP热出力 P_gt = sdpvar(1, 24); % 燃气轮机电出力 P_buy = sdpvar(1, 24); % 购电功率 P_eb = sdpvar(1, 24); % 电锅炉耗电 H_eb = sdpvar(1, 24); % 电锅炉产热 S = sdpvar(1, 24); % 储热罐蓄热量 H_char = sdpvar(1, 24); % 储热罐充热功率 H_dis = sdpvar(1, 24); % 储热罐放热功率

约束尽量用向量形式写，不要在for循环里一条条加。Yalmip对向量约束支持得很好。比如风电出力上下限，直接写：

C = [C, 0 <= P_wind <= P_wind_pred, 0 <= P_pv <= P_pv_pred];

这比循环里24次C = [C, P_wind(t) <= P_wind_pred(t)]要快得多，而且不容易维度出错。但如果你需要逐时段赋不同的系数（比如分时电价对应的购电成本），可以把系数向量乘上去，仍然是一行。

3.3 核心代码框架示例

给一个浓缩版主干，流程是：读数据、定义变量、写目标、写约束、求解、取结果。

%% 数据准备 T = 24; P_load_e = data.P_load_e; % 电负荷 1x24 H_load = data.H_load; % 热负荷 1x24 P_wind_pred = data.P_wind_pred;% 风电预测 1x24 P_pv_pred = data.P_pv_pred; % 光伏预测 1x24 c_el = data.c_el; % 分时电价 1x24 eta_eb = 0.95; % 电锅炉效率 P_eb_max = 1; % 电锅炉最大功率 MW S_max = 4; H_char_max = 1; H_dis_max = 1; % 储热罐参数 %% 变量定义 clear P_gt P_chp H_chp P_wind P_pv P_buy P_eb H_eb S H_char H_dis P_gt = sdpvar(1,T); P_chp = sdpvar(1,T); H_chp = sdpvar(1,T); P_wind = sdpvar(1,T); P_pv = sdpvar(1,T); P_buy = sdpvar(1,T); P_eb = sdpvar(1,T); H_eb = sdpvar(1,T); S = sdpvar(1,T); H_char = sdpvar(1,T); H_dis = sdpvar(1,T); %% 目标函数 c_chp_el = 0.35; c_chp_h = 0.15; % CHP电/热单位成本 c_gas = 0.22; c_pen = 500; % 燃气轮机成本、弃能惩罚 obj = sum(c_el .* P_buy) + c_gas*sum(P_gt) + c_chp_el*sum(P_chp) + c_chp_h*sum(H_chp) ... + c_pen*(sum(P_wind_pred - P_wind) + sum(P_pv_pred - P_pv)) ... + 0.01*sum(P_eb); %% 约束 C = []; % 电/热平衡 C = [C, P_gt + P_chp + P_wind + P_pv + P_buy - P_eb == P_load_e]; C = [C, H_chp + H_eb + H_dis - H_char == H_load]; % 新能源消纳 C = [C, 0 <= P_wind <= P_wind_pred, 0 <= P_pv <= P_pv_pred]; % 电锅炉 C = [C, H_eb == eta_eb * P_eb, 0 <= P_eb <= P_eb_max]; % CHP可行域 % 这里假设CHP电出力0.5~3MW，热出力0.2~2.5MW，并用两个耦合约束模拟多边形 C = [C, 0.5 <= P_chp <= 3, 0.2 <= H_chp <= 2.5]; C = [C, P_chp + 0.6*H_chp >= 0.8, P_chp + 0.8*H_chp <= 4.5]; % 燃气轮机约束 C = [C, 0.3 <= P_gt <= 2]; % 储热罐 C = [C, S(2:T) == S(1:T-1) + eta_c*H_char(1:T-1) - H_dis(1:T-1)/eta_d - epsilon*S(1:T-1)]; C = [C, S(end) == S(1)]; C = [C, 0 <= S <= S_max, 0 <= H_char <= H_char_max, 0 <= H_dis <= H_dis_max]; % 爬坡约束 C = [C, -0.5 <= diff(P_gt) <= 0.5]; %% 求解 ops = sdpsettings('solver','gurobi','verbose',1); optimize(C, obj, ops); %% 输出结果 P_chp_opt = value(P_chp);

注意代码里的S(2:T)和S(1:T-1)这种向量写法，Yalmip会自动展开成T-1条等式，非常简洁。diff(P_gt)是相邻时段的差值，也支持直接写进约束。如果报错，多半是向量方向问题，把列向量转成行向量就好。

3.4 结果可视化：把调度曲线画成图

优化完成后，用value()取出各变量数值，就可以画图了。我习惯画四张图：电功率平衡堆叠面积图、热功率平衡堆叠面积图、弃风弃光柱状图、储热罐SOC曲线。画图核心思路是用area函数把不同设备出力叠起来，和负荷曲线对比。

figure; area(1:T, [value(P_wind); value(P_pv); value(P_chp); value(P_gt); value(P_buy)]') hold on plot(1:T, P_load_e, 'k-', 'LineWidth', 2); legend('风电','光伏','CHP','燃气轮机','购电','电负荷');

注意面积图的输入矩阵是“时段×设备”，所以要对设备出力向量做转置。area会按列堆叠，顺序要和legend对应。储热罐SOC画折线就行。导出图片时，我一般用set(gcf,'Color','w')避免灰色底，再exportgraphics(gcf,'result.png','Resolution',300)，直接满足论文插图要求。

4. 实操过程与核心环节实现：从数据准备到结果复现

4.1 基础输入数据的准备

没有靠谱的数据，模型就是个空壳。我这儿给一套典型的参数，适合自己跑着玩：假设一个园区系统，电负荷峰值8MW，热负荷峰值5MW，风力装机4MW，光伏装机2MW，CHP装机电功率3MW/热功率2.5MW，燃气轮机装机2MW，电锅炉1MW，储热罐容量4MWh，最大充放热功率1MW。分时电价采用峰谷两段式：峰时段（10:00-15:00、18:00-21:00）0.8元/kWh，谷时段（23:00-7:00）0.3元/kWh，其他时段0.5元/kWh。风电预测曲线做成“夜间大、白天小”的形状，光伏则是“中午大、两端小”。这些数据我放在Excel里，Matlab用readtable读出来，再转成行向量。

设备参数可以整理成一张表：

设备	参数	数值
CHP机组	电出力下限/上限	0.5/3 MW
CHP机组	热出力下限/上限	0.2/2.5 MW
CHP机组	电成本系数	0.35 元/kWh
CHP机组	热成本系数	0.15 元/kWh
燃气轮机	出力下限/上限	0.3/2 MW
燃气轮机	单位燃料成本	0.22 元/kWh
电锅炉	电热转换效率	0.95
电锅炉	最大电功率	1 MW
储热罐	容量	4 MWh
储热罐	最大充/放热功率	1 MW
储热罐	充/放热效率	0.9
储热罐	散热损耗	0.05/h
弃风/弃光惩罚	单价	500 元/MWh

这些参数不用太较真，重点是量纲统一：功率用MW，热量也用MW（表示热功率），蓄热量用MWh，时间步长1小时。如果热负荷的单位是GJ/h，记得先除以3600换算成MW，否则等式平衡就完全对不上。单位出错是新手最容易掉进去的坑。

4.2 模型搭建中的几个关键调试节点

代码写完后，我先做“空跑测试”：把风电、光伏预测值全设为0，电锅炉和储热罐也都不能动，让模型退化成纯购电加燃气轮机加CHP的传统系统。这一步验证的是基础平衡约束和机组约束是否写得对。如果空跑都无解，那绝对不是新能源约束的问题，而是平衡等式、爬坡或可行域哪里“打架”了。第二步，把风电预测按比例放大，从0.5倍逐步调到2倍，观察弃风量变化。弃风量应当随预测增加而增加，并且增幅合理。如果预测功率明明很高，但弃风量为负，那一定是你把P_wind_pred - P_wind写反了。

第三步，临时把储热罐的周期约束去掉，看看SOC终值会掉到哪里。如果终端SOC和初值差距大，说明周期约束确实在起作用。这三个测试都通过，模型基本稳了，剩下的就是调惩罚系数和爬坡参数。

4.3 典型调度结果的解读

以冬季典型日为例，调度结果大概是这样的：深夜电价低谷、风电预测很高，此时CHP机组尽量压低电出力，但要满足热负荷，所以热出力维持下限；不足的热量由电锅炉在谷电时段补充，同时储热罐充电，把热量存起来。白天进入峰价时段后，CHP电出力回升，储热罐放热，替代部分电锅炉出力，降低购电成本。到晚上用电高峰，购电成本高，系统会启动燃气轮机发电，配合储热罐放热，实现“削峰填谷”。

从指标上看，如果没有电锅炉和储热罐，传统的“以热定电”模式下夜间弃风率可能到10%以上；加入电锅炉和储热罐后，弃风率能压到2%以内。你还可以特意跑一个“无电锅炉、无储热罐”的对照组，和完整模型对比，两张图放在一起，电-热灵活资源的价值一眼就能看出来。这种对比特别适合写论文的算例部分。

5. 常见问题与排查技巧实录

5.1 求解器返回无解或状态码inf

最常见的原因是约束集合不可行——某时段电热平衡根本无法同时满足。排查方法是把目标函数暂时屏蔽，只求解可行性问题：optimize(C, [], ops)，然后看check(C)的结果，哪条约束残差大于1e-4，问题就在哪。我遇到过几次无解，都是储热罐周期约束和SOC下限矛盾：比如初始SOC已经低到1MWh，而容量下限设置成2MWh，自然无解。把初值与容量下限统一，或者在周期约束中允许一定容差，问题就消失了。

5.2 弃风量异常大，先查惩罚系数

如果你把弃风惩罚系数设为10元/MWh，而峰时电价0.8元/kWh（折合800元/MWh），模型在“弃风少发电”和“购电满足负荷”之间做选择时，很可能选择弃风，因为弃风惩罚太便宜了。建议惩罚系数至少设为购电价上限的5到10倍，我默认设500元/MWh。另外，风电出力的上限向量必须是预测值，写成P_wind <= P_wind_pred，可别把预测值直接当成固定值赋给变量——那就不是在“消纳决策”了，而是在跑固定出力场景。

5.3 Yalmip报错“VARIABLE DOESN'T EXIST”或维度不匹配

这类问题十有八九是变量定义和作用域搞混了。确保所有sdpvar在同一个workspace里；如果写了函数封装，变量必须在函数内部重新定义。另一个高频坑是行列方向不一致：Yalmip对1×24和24×1非常敏感，从Excel读进来的列向量忘了转置，约束维度对不上，优化直接报错。我的习惯是读入数据后统一加一句data = data(:)';强制变行向量，后面所有变量定义都用行向量，基本不会再出维度问题。

5.4 求解时间太长，如何优化性能

24小时单节点模型，正常求解时间在毫秒到秒级。如果特别慢，先在sdpsettings里把'verbose',0关掉，减少输出。其次，检查是不是在约束里用了大量循环，尽量用向量化写法。Yalmip约束对象的数量本身会影响求解器前处理速度，把同类型约束一次性拼接，不要拆成几十条独立约束。最后，如果模型规模特别大，比如几百个节点、几千条约束，可以把求解器内部参数调一下，比如Gurobi的LPWarmStart开成2，或者把'ShowProgress'设为0。我做过对比：同样规模模型，从循环逐条加约束改成向量化，求解时间可以从200ms降到80ms，效果明显。

5.5 储热罐SOC曲线疯狂振荡

如果SOC曲线要么冲到上限要么贴下限，大概率是充放热功率约束没写对，或者效率和热损耗系数设得太极端。调试时先把散热损耗设为0，看SOC曲线是否平滑；再把充放热效率调成1，排除效率把状态量“吞掉”的可能。另外，如果储热罐同时充放正负交替，大概率是你没有禁止同步充放。LP模型里目标函数是成本最小化，通常不会出现既充又放的浪费行为，但如果惩罚项导致某种利润模式出现同步充放，就需要加0/1变量做互斥约束。我的模型由于目标函数只惩罚弃能、不奖励套利，没出现这个问题，如果你要做电价套利，就要小心这个坑。

5.6 电力平衡中有缺口，但检查所有约束都是对的

这时候去看看电锅炉。电锅炉耗电功率进了电平衡，产热功率进了热平衡，两者通过效率等式连接，这一步最容易漏。我见过好多次：电平衡右侧忘了加-P_eb，热平衡右侧又加了H_eb，结果电负荷缺口恰好等于电锅炉耗电，模型还“自我感觉良好”。把平衡等式拉出来，逐项核对变量正负号，尤其关注那些出现在两个平衡方程里的设备。

结尾

做完这套电热综合能源系统日前经济调度模型，我最大的体会是：建模过程最花时间的不是数学公式，而是怎么把实际物理设备的运行区间“翻译”成约束矩阵。CHP可行域、储热罐周期约束、弃风惩罚这些细节，一个地方没处理好，结果就会疯。解决的办法就是多做几个算例，多画几张图，对着结果看设备有没有“违规出工”。最后再分享一个小技巧：在Matlab里写大矩阵约束时，多用repmat和kron构造系数矩阵，再配合Yalmip的向量化约束，整个模型的调试体验会舒服很多。这个模型现在只是个确定性调度版本，如果后续要做随机优化或者分布鲁棒，把风电预测误差场景接入，本质上就是在这个骨架上加一堆场景复制，思路是相通的。