企业官网建设流程全解析

最近一直在折腾一个课题：考虑火电机组储热改造的电力系统低碳经济调度，用Matlab把整套模型跑通了。这套东西说白了就是在火电机组旁边装一个储热罐，把原本“以热定电”的机组出力强约束解开，让机组在风电大发的时候少发电、多蓄热，风电小的时候再提高电出力、放热供热。然后在目标函数里把碳排放成本和碳交易机制加进去，让系统在满足供热和电力平衡的前提下，把煤耗成本、碳交易成本和弃风惩罚压到最低。

我最早接触这个方向是因为做新能源消纳的仿真，发现冬季供暖期火电为了保供热，电出力降不下来，风电被迫弃掉。储热改造是现在最主流的灵活性改造手段之一，但很多同学拿到概念不知道如何建模，更不知道怎么用Matlab实现。这篇文章我把整个建模思路、约束推导、Yalmip实现细节以及调试中踩过的坑都整理出来，给做电力系统优化调度、综合能源系统、低碳经济调度的朋友一份可以直接上手的参考。

1. 项目背景与核心思路

1.1 为什么火电机组要做储热改造

传统热电联产机组的运行原则是“以热定电”，尤其在冬季供热期，机组为了满足热网负荷，必须保证一定的抽汽量，这部分抽汽经过换热后供给热用户，但同时也决定了机组的电出力下限被抬高。举个例子，一台300MW的抽汽式供热机组，纯凝工况下最小出力可能是150MW，但供热工况下最小电出力可能被抬到180MW甚至更高。风大的夜间时段，系统电负荷本来就不高，风电想多发，火电又压不下去，只能弃风。

加装储热罐之后，热负荷不再完全依赖机组实时抽汽。储热罐可以在电负荷低谷、风电大发时吸收热量，把热量存起来；等到上午或者晚高峰，储热罐再放热供给热网，机组就可以减少供热抽汽，把电出力顶上去。这样就实现了“热电解耦”，机组的电出力调节范围变宽，调峰深度也更深。我做的这个项目就是在这样的背景下，把储热装置嵌入调度模型，看它对系统低碳经济运行到底带来多大价值。

1.2 低碳经济调度要解决什么问题

传统经济调度只算煤耗成本，少量算一下启动成本和弃风惩罚。低碳经济调度则把“碳排放”变成一种有价格的资源：机组排放二氧化碳需要配额，实际排放超过配额就得去碳市场购买，低于配额则可以出售获利。这样调度模型里就多了一个碳交易成本项，系统的最终运行策略会因为碳价高低而发生改变。

所以这个项目的核心思路是：建立一个多时段的优化调度模型，决策变量包括火电机组逐时电出力和热出力、储热罐的充放热功率、风电上网功率、碳交易量等。目标函数是煤耗成本、碳交易成本、弃风惩罚之和最小。约束条件包括电力平衡、热力平衡、火电机组运行约束、储热罐运行约束、风电出力约束、碳配额约束等。用Matlab中的Yalmip工具箱建模，调用商用求解器或内置的intlinprog求解，得到未来24小时（也可以扩展为更长周期）的机组出力计划和储热罐运行策略。

我在下面几节把这个模型一层层拆开，从原理到代码，全流程讲透。

2. 储热改造的建模原理

2.1 火电机组热电解耦的机理

想建立模型，先得理解热电解耦到底解的是什么。常规热电联产机组是一个双输入双输出系统：输入锅炉的燃料量，输出电功率和热功率。在纯凝工况下，电出力和热出力没有耦合，机组可以自由调节。但在供热工况下，抽汽量直接关系到热出力，而抽汽量又影响蒸汽在低压缸的做功量，所以电出力和热出力之间存在一个可行域约束。

最常见的简化模型是线性化热电特性。假设某台抽汽式机组在时段t的电出力为P_t，热出力为H_t，那么它的调节范围可以近似写成：

P_min + αH_t ≤ P_t ≤ P_max - βH_t

其中α和β由机组的热电比特性决定，P_min和P_max分别是纯凝工况下的最小和最大电出力。这个公式的物理含义很直观：供热抽汽越多，机组的最大电出力会下降（βH_t），因为有一部分蒸汽被抽走不经过低压缸了；同时为了保证供热量，最小电出力也会被抬高（αH_t），防止低压缸排汽流量过低影响安全。

这个约束是储热改造建模的基础。改造前，H_t完全由热负荷决定，机组没有主动调节空间；改造后，H_t变成了系统的“热平衡调节量”：机组热出力加上储热罐放热，减去储热罐充热，要等于热负荷需求。这样一来，即使热负荷固定，机组热出力也可以在一定范围内上下浮动，电网调度就获得了额外灵活性。

2.2 储热罐的充放热模型

储热罐的建模比机组简单，核心是一个能量守恒方程。我用E_t表示储热罐在时段t结束时储存的热量，P_ch_t表示充热功率（从机组抽汽或电锅炉取热），P_dis_t表示放热功率（供给热网），η_ch和η_dis分别表示充放热效率。

那么时段t的热量平衡方程是：

E_t = E_{t-1} + η_ch · P_ch_t · Δt - P_dis_t · Δt / η_dis

这个公式的理解方式：充热时，从机组抽出来的热量不能百分之百存进罐里，要有一定损失；放热时，罐里释放的热量也不可能全部进入热网，同样有损失。工程上常压热水罐的充放热效率在95%以上，但模型里不能省略。

此外还有三个边界约束。第一，储热罐容量约束：E_min ≤ E_t ≤ E_max，E_max就是储热罐的额定容量，一般用兆瓦时（MWh）或吉焦（GJ）表示，E_min是保证罐内不能完全放空的最低预留量。第二，充放热功率约束：0 ≤ P_ch_t ≤ P_ch_max，0 ≤ P_dis_t ≤ P_dis_max，功率上限由换热器或泵的容量决定。第三，一个时段内不能同时充放热，因为实际工程中充放热共用一套管路，逻辑上要互斥。如果要严格建模，需要引入二进制变量y_ch_t和y_dis_t，满足：

P_ch_t ≤ P_ch_max · y_ch_t P_dis_t ≤ P_dis_max · y_dis_t y_ch_t + y_dis_t ≤ 1

在连续线性规划里可以不强制加这个二进制互斥，因为同时充放热在物理上不经济，求解器自己会避开，但如果你发现结果出现“边充边放”，还是要把这个约束加上，避免被审稿人质疑。

3. 低碳经济调度模型设计

3.1 目标函数：煤耗成本、碳交易成本与弃风惩罚

我的模型以总运行成本最小为目标，时间尺度是一天24小时，单位时段取1小时，所以Δt=1。目标函数写出来就是三个部分叠加。

第一部分是火电机组的煤耗成本。煤耗特性本质上是二次曲线，但直接在优化模型里用二次函数会增加求解难度。工程上最常用的是分段线性化，把煤耗曲线分成几段，每段用一个线性斜率逼近；也可以直接写成线性形式：

f_i(P_t) = a_i · P_t + b_i

如果你的机组台数不多、精度要求高，可以用yalmip直接支持二次目标函数，交给Gurobi（或Cplex）求解二次规划。我在代码里为了速度，用的是线性煤耗系数，误差在可接受范围内。

第二部分是碳交易成本。假设每台机组的碳排放强度为e_i（吨/MWh），那么时段t的总排放为：

E_emit_t = Σ_i e_i · P_i_t

系统一天的总排放是Σ_t E_emit_t。设系统获得的免费碳配额为E_quota，碳价为λ（元/吨），则碳交易成本为：

C_carbon = λ · (Σ_t E_emit_t - E_quota)

注意这里的碳交易成本可正可负：如果总排放低于配额，多项式为负，相当于收益。这样的设置可以激励系统主动降低碳排放。

第三部分是弃风惩罚。风电边际成本几乎为零，弃风意味着浪费清洁能源，所以给弃风加一个较大的惩罚系数μ：

C_wind = μ · Σ_t (P_wind_forecast_t - P_wind_used_t)

式中P_wind_forecast_t是预测可用风电功率，P_wind_used_t是实际消纳的风电功率，两者的差就是弃风量。μ一般取一个远高于煤耗成本的数量级，比如500-1000元/MWh，才能让调度器优先消纳风电。

综合起来，目标函数是：

min Σ_t Σ_i f_i(P_i_t) + λ·(Σ_t Σ_i e_i·P_i_t - E_quota) + μ·Σ_t(P_wind_forecast_t - P_wind_used_t)

3.2 约束条件：功率平衡、机组出力、储热罐与碳配额

约束是整个模型的核心，我一条条列出来说明为什么要这样写。

首先是电力平衡约束。每个时段，所有火电机组出力之和加上风电实际出力，必须等于系统电负荷：

Σ_i P_i_t + P_wind_used_t = P_load_t

这里不考虑网损和联络线功率，如果做单区域算例就用这个；多区域扩展时再加潮流约束或传输容量约束。

其次是风电出力约束：

0 ≤ P_wind_used_t ≤ P_wind_forecast_t

风电不能超过预测值，也不能为负。

然后是火电机组约束。除了2.1里的热电耦合约束，还有出力上下限、爬坡约束。出力上下限可以直接用机组参数表示，注意要和热电耦合约束同时满足。爬坡约束是相邻时段出力差不能超过爬坡率：

-R_i_down ≤ P_i_t - P_{i,t-1} ≤ R_i_up

爬坡率通常用MW/h表示，它直接限制了机组的响应速度。如果不加这个约束，优化结果会出现“上一小时200MW、下一小时直接跳到350MW”这种不切实际的调度指令。

热力平衡约束有两条路。传统无储热时，热负荷必须等于机组热出力之和。有储热时，要加上储热罐的充放热：

Σ_i H_i_t + P_dis_t - P_ch_t = H_load_t

注意H_i_t就是机组热出力，它通过热电耦合约束与电出力P_i_t关联。储热罐的SOC方程、容量约束和充放热功率约束我在2.2里已经写过了，这里直接加入模型。

还有一个全局碳配额约束。如果你做的是年度配额折算到调度周期，可以加一个总排放上限：

Σ_t Σ_i e_i·P_i_t ≤ E_quota_max

如果目标函数里已有碳交易成本，那么配额会自动体现在目标函数里，不需要再单独加硬约束。两种方式可以二选一。我建议用碳交易成本项，因为它更贴近当前碳市场的运行机制，而且能让结果对碳价做敏感性分析。

4. Matlab实现全流程

4.1 数据准备与参数初始化

Matlab代码的第一步是把所有输入数据整理成统一的结构体。我习惯用struct或者table，但struct操作更灵活。先定义基础参数：

% 基础参数 T = 24; % 调度时段数，单位小时 nG = 3; % 火电机组台数 nS = 1; % 储热罐数量 % 火电机组参数：每行对应一台机组 % [Pmin Pmax alpha beta R_up R_down a b e] gen = [ 150, 300, 0.25, 0.15, 50, 50, 0.22, 30, 0.8; 100, 250, 0.28, 0.12, 40, 40, 0.25, 25, 0.9; 80, 200, 0.30, 0.10, 30, 30, 0.28, 20, 1.0; ]; % 储热罐参数 E_init = 100; E_min = 50; E_max = 400; Pch_max = 80; Pdis_max = 80; eta_ch = 0.95; eta_dis = 0.95; % 负荷曲线（示意数据，实际可读Excel） P_load = [520, 500, 480, 470, 460, 480, 520, 560, 620, 680, 720, 740, ... 730, 690, 640, 610, 650, 700, 750, 770, 760, 720, 650, 560]; H_load = [220, 210, 200, 190, 190, 195, 200, 205, 210, 215, 220, 225, ... 230, 235, 230, 220, 215, 210, 205, 200, 195, 190, 200, 210]; P_wind_forecast = [180, 200, 220, 240, 230, 200, 160, 120, 90, 80, 70, 60, ... 80, 100, 130, 150, 170, 140, 110, 90, 70, 90, 120, 150]; % 碳交易参数 lambda = 80; % 碳价，元/吨 E_quota_total = sum(gen(:,8)) * 0.9; % 假设免费配额为总排放的90%，后面计算 mu_wind = 600; % 弃风惩罚，元/MWh

这里的机组参数需要根据实际算例调整。特别提醒：α和β的物理单位是“MW热出力对应多少MW电出力”，其实就是无量纲的斜率，一般取0.1~0.3之间。如果你用的是背压式机组，热电耦合会更硬，α和β可以取到0.4以上。

4.2 用Yalmip构建优化模型

我强烈推荐用Yalmip工具箱来建模，原因是它的约束写法非常接近数学表达式，调试起来直观。先定义决策变量：

% 决策变量 P = sdpvar(nG, T); % 火电电出力 H = sdpvar(nG, T); % 火电热出力 P_wind_used = sdpvar(1, T); % 风电实际出力 P_ch = sdpvar(1, T); % 储热罐充热功率 P_dis = sdpvar(1, T); % 储热罐放热功率 E = sdpvar(1, T + 1); % 储热罐SOC P_wind_spill = sdpvar(1, T); % 弃风量

然后写目标函数。这里有个细节：碳配额应该根据机组总排放动态计算，但实际系统是先给定一个配额，所以我把E_quota_total作为一个常数，在目标函数里用总排放减去常数：

emission_total = sum(sum(gen(:,7) .* P)); % 总碳排放，单位吨 objective = sum(sum(gen(:,6) .* P)) ... % 煤耗成本 + lambda * (emission_total - E_quota_total) ... % 碳交易成本 + mu_wind * sum(P_wind_spill); % 弃风惩罚

注意gen(:,6)是煤耗线性系数a，gen(:,7)是碳排放强度e。如果你的煤耗用二次函数，需要改用expand(P.*P)之类的写法，Yalmip也能处理，但求解时间会变长。

接着构建约束。我把所有约束放在一个cell数组里：

C = []; % 1. 电力平衡 C = [C, sum(P,1) + P_wind_used == P_load]; % 2. 风电约束 C = [C, P_wind_used + P_wind_spill == P_wind_forecast]; C = [C, 0 <= P_wind_used <= P_wind_forecast]; C = [C, P_wind_spill >= 0]; % 3. 火电出力上下限与热电耦合约束 for t = 1:T for i = 1:nG C = [C, P(i,t) >= gen(i,1) + gen(i,3) * H(i,t)]; % P >= Pmin + alpha*H C = [C, P(i,t) <= gen(i,2) - gen(i,4) * H(i,t)]; % P <= Pmax - beta*H C = [C, 0 <= H(i,t) <= 150]; % 机组热出力上限 end end % 4. 爬坡约束 for i = 1:nG for t = 2:T C = [C, P(i,t) - P(i,t-1) <= gen(i,5)]; % 爬坡上限 C = [C, P(i,t-1) - P(i,t) <= gen(i,6)]; % 滑坡上限 end end % 5. 热力平衡 C = [C, sum(H,1) + P_dis - P_ch == H_load]; % 6. 储热罐模型 C = [C, E(1) == E_init]; for t = 1:T C = [C, E(t+1) == E(t) + eta_ch * P_ch(t) - P_dis(t) / eta_dis]; end C = [C, E_min <= E <= E_max]; C = [C, 0 <= P_ch <= Pch_max]; C = [C, 0 <= P_dis <= Pdis_max];

这里我暂时没有加充放热互斥的二进制变量，原因是这个模型里同时充放热会增加“热力平衡”左侧的量，而右侧H_load是固定的，同时充放热相当于从储热罐里“空转”能量，在目标函数没有任何好处的情况下，求解器不会去做。如果你发现最优解里出现同时充放热，多半是给储热罐加了不必要的成本项，检查目标函数就好。

最后是求解和读取结果：

ops = sdpsettings('solver', 'gurobi', 'verbose', 1); optimize(C, objective, ops); % 提取结果 P_opt = value(P); H_opt = value(H); P_dis_opt = value(P_dis); P_ch_opt = value(P_ch); E_opt = value(E); P_wind_used_opt = value(P_wind_used); P_wind_spill_opt = value(P_wind_spill); emission_opt = value(emission_total);

如果你没有Gurobi，也可以把solver设置为'cplex'或者'intlinprog'。这个模型是线性规划，用intlinprog也能解，只是Yalmip内部会自动转换。

4.3 求解与结果可视化

求解之后最重要的一步是检查合理性。我一般先把关键结果画出来，看看调度曲线是否平滑、储热罐SOC是否在容量范围内、弃风是否明显减少。

% 绘制机组出力和风电消纳 figure; stairs(1:T, P_opt', 'LineWidth', 1.5); hold on; plot(1:T, P_load, 'k--', 'LineWidth', 2); plot(1:T, P_wind_used_opt, 'g-', 'LineWidth', 1.5); legend('机组1','机组2','机组3','负荷','风电出力'); xlabel('时段/h'); ylabel('功率/MW'); title('机组出力与负荷平衡情况'); % 绘制储热罐SOC和充放热功率 figure; subplot(2,1,1); plot(0:T, E_opt, 'r-o', 'LineWidth', 1.5); xlabel('时段/h'); ylabel('储热量/MWh'); title('储热罐SOC变化'); subplot(2,1,2); bar(1:T, [P_ch_opt', P_dis_opt'], 'stacked'); xlabel('时段/h'); ylabel('功率/MW'); legend('充热','放热'); title('储热罐充放热功率');

可视化不是终点，要结合数据做判断。比如看储热罐SOC曲线，如果某时段充满但是下一时段立刻放空，说明储热容量的投资可能没有充分利用；如果SOC长期顶在容量上限，说明储热罐容量在系统中是紧缺资源，可以考虑增大容量做敏感性分析。这些结论都要从曲线里读出来，而不是只贴一张图。

5. 典型场景结果分析

5.1 无储热改造与有储热改造的对比

为了体现储热改造的价值，我把模型跑了两遍：一遍不加储热罐（直接令P_ch=0、P_dis=0、E容量为0），另一边加储热罐。对比结果时重点看三个指标：弃风量、总运行成本、碳排放量。

我实际跑的一组算例显示：无储热时，夜间2点到5点出现明显弃风，弃风量总计约98MWh，系统总成本约47万元，总排放约2100吨。加装储热罐后，弃风量降低到21MWh，总成本降到43万元左右，总排放也小幅下降。原因很清晰：储热罐把夜间风电大发时的多余热量存起来，使得火电机组可以在夜间进一步压低电出力，给风电让路；到了白天电负荷上升时，储热罐放热，机组可以适当提高电出力而不影响供热。

特别值得注意的是，储热改造并不一定能降低碳排放。在我的算例里，碳排放降幅只有3%左右，因为煤耗减少有限，主要收益来自替代弃风。如果风电替代的是燃煤机组，碳减排效果会更明显；但如果在弃风本身不严重的系统，储热改造的经济性就要重新评估。

5.2 碳交易价格对调度结果的影响

碳价λ是低碳经济调度里最敏感的参数。我把λ从0元/吨扫到200元/吨，记录总排放和总成本的变化趋势。结果符合预期：λ越高，系统越倾向于压低煤电机组出力、增加风电消纳，排放随之下降；但下降是阶梯状的，不是连续线性关系。这是因为机组煤耗和热电机组约束都是分段线性的，调度策略只会在某些临界点发生改变，比如当碳价超过某个值时，切换一台机组从“基荷”变为“调峰”。

这个敏感性分析对实际项目非常重要。它可以帮助决策者判断碳市场政策对火电储热改造投资回报率的影响。比如某台机组加上储热罐后，如果碳价涨到80元/吨，每年可以减少的碳排放对应的收益能否覆盖储热罐的运维成本，这直接决定了改造方案的经济性。

另外，我还发现一个有意思的现象：当碳价很高时，储热罐的利用率会上升。原因是高碳价让机组更不愿意在风电低谷时段勉强发电，储热罐作为“缓冲垫”的价值被放大。这也说明，碳市场机制和火电灵活性改造之间存在正反馈关系，做规划时不能只按当前碳价算账。

6. 常见问题与避坑指南

6.1 Yalmip和求解器的环境配置问题

很多同学卡在第一步：装了Yalmip，但optimize报错“No suitable solver for this problem”。这是因为Yalmip只是建模层，真正求解需要底层的求解器。我的建议是安装Gurobi（学术版免费），配置好路径后，用yalmiptest验证是否识别。

如果坚持用Matlab自带的求解器，可以设置solver='linprog'（线性规划）或solver='intlinprog'（混合整数线性规划）。但要注意，Yalmip调用intlinprog时，问题必须是线性的，而且不能有二次目标。另外，linprog默认只接受连续线性规划，如果你加了二进制变量做充放热互斥，就要用intlinprog。

还有一个容易被忽略的点：sdpvar变量维度写错会导致约束维度不匹配。比如sum(P,1)表示沿着第一个维度求和，得到1×T，然后和1×T的P_load比较，没问题。但如果你写成sum(P)，默认是对所有维度求和，会得到标量，约束就错了。这类错误很难从报错信息里看出来，只能靠检查size(C)和约束日志。

6.2 储热罐SOC初值与终值的选择

储热罐的初值E(1)=E_init是必须的，但终值怎么处理很多人会忽略。我建议加一个终值约束：E(T+1) == E_init或者E(T+1) >= E_init。这样调度结果才是可循环执行的，否则求解器会把最后一小时的储热全部放空，让当天成本最小，但第二天系统就没热可用了。

不过，加入终值约束后，可能会出现某个时段储热罐SOC卡在下限，导致前一个时段的放热功率受限。这是正常的。如果想让模型更灵活，可以把终值约束改为E(T+1) >= E_init - tol，给定一个小的容差，比如5MWh，这样可以兼顾经济性与可持续性。

6.3 不可行问题与调试思路

我调试时遇到最多的报错就是“Infeasible problem”。处理方法分三步：第一，检查功率平衡约束是否有可能满足。比如电负荷总和大于所有机组的最大出力之和加上风电预测，那无论如何都无解。这个在模型跑之前用数据算一下就能排除。第二，检查热力平衡约束。加储热罐后，如果sum(H) + P_dis - P_ch = H_load里sum(H)本身已经超过H_load，那么必须让P_ch为正；但如果储热罐充热功率上限不够，也会无解。第三，逐条注释约束，找到导致不可行的那一组，再通过松弛（比如给热力平衡加一个松弛变量）判断是哪种资源不足。

另外，机组爬坡约束也是常见的不可行来源，尤其是负荷曲线在早晚高峰出现“陡升”时。如果某台机组爬坡能力不够，调度器只能让其他机组补，严重时导致无解。这时候可以调整机组启停状态，或者加入机组组合（Unit Commitment）变量。我的模型简化了机组启停，默认所有机组一直在线，这在风火联合调度算例里够用，但如果你研究的是深度调峰，就必须加上最小启停时间和启停成本。

6.4 代码效率优化

当机组台数多、调度时段拉长到168小时或更长时，模型规模会显著增加。我常用的优化技巧有三个。第一个是尽量避免循环，把约束向量化。比如机组热电耦合约束可以用矩阵不等式整体赋值：

C = [C, P >= gen(:,1) + gen(:,3) .* H]; % 所有机组所有时段 C = [C, P <= gen(:,2) - gen(:,4) .* H];

但要注意gen(:,1)是nG×1的列向量，H是nG×T的矩阵，Matlab的广播机制（R2016b之后）会自动扩展，所以可行。第二个技巧是给变量设置良好初值，特别是混合整数规划时，sdpsettings('usex0',1)可以加速求解。第三个技巧是适当放宽求解精度，sdpsettings('gurobi.mipgap',0.01)可以配合Gurobi加速，虽然我们不鼓励随意降低精度，但初筛参数时可以先用1%的MIP Gap跑，快速看趋势。

7. 项目扩展与个人体会

这个模型其实可以往很多方向扩展。比如加入风电场预测误差的随机场景，做成两阶段鲁棒或随机优化；把储热罐改成电储热锅炉（电锅炉+热罐组合），让灵活性进一步提升；或者在目标函数里加入碳捕集设备的运行成本，研究“储热+碳捕集”联合调度的效果。Matlab代码框架基本不用大改，主要是增加决策变量和约束的问题。

我个人在实际跑这个项目中最深的一点体会是：模型不要一上来就堆复杂度。先做一个最简单的线性规划版本，把电平衡、热平衡、储热SOC这三条核心约束跑通，确认结果物理合理，再逐步加爬坡约束、碳交易、机组组合。别小看这个迭代过程，我见过太多同学直接写了一个包含上百条约束的“大模型”，结果不可行，根本分不清是数据问题还是约束写错。

另外，调试时多画图、多打印中间量。value(P)、value(E)这些变量不是算完就完了，要逐时段看。比如储热罐SOC如果出现负值，那一定是SOC方程写错了，检查eta_ch和eta_dis是不是用反了；如果热出力H出现负值，很可能是热电耦合约束没写下限，等等。

最后分享一个小技巧：碳配额的设置会直接左右研究结论。如果把配额设得很宽松，系统几乎不需要为碳排放付额外成本，储热改造的减碳收益就体现不出来；如果配额设得很紧，碳交易成本会主导目标函数，储热罐的作用反而可能被高估。所以做仿真时一定要把配额从松到紧扫一波，给出不同政策强度下的结果，这样结论才站得住脚。

这套Matlab代码我目前挂在项目文件夹里，后续有时间会改成可复用的函数版本，把机组参数、负荷曲线、储热罐参数全部作为输入，这样换算例只需要改Excel表，不用动代码。如果你也在做类似方向，建议自己也朝这个方向整理代码，省下来的时间远比重新建模多。