把“考虑火电机组储热改造的电力系统低碳经济调度”用 Matlab 完整跑通一遍之后,我发现这个方向的价值不在于多一个算例,而在于它把火电灵活性改造、碳市场和优化调度三个问题拧在一起,值得动手复现。这篇文章我打算直接从问题源头讲起,把数学模型怎么建、储热罐怎么处理、碳成本怎么加、Matlab 代码怎么写、求解器怎么调,一条线捋清楚。适合正在做电力系统方向课程设计、毕业设计,或者刚入门优化调度的同学参考,有一点 Matlab 基础就能跟上。
1. 火电机组储热改造:到底改了些什么
1.1 传统热电联产的“以热定电”约束
北方地区冬季供暖期,大量热电联产机组承担着供热任务。这类机组内部的热电关系非常密切,供热抽汽量一旦确定,电功率的下限就基本被“锁”住了。通俗讲,热负荷越高,机组为了维持供暖而必须保持的蒸发量就越大,电出力哪怕不需要那么高,也无法继续往下压,这就叫“以热定电”。
这个约束在冬天会出现一个很尴尬的局面:夜间风大、负荷低,风电出力正高,但热电联产机组因为要供热,电出力压不下来,结果电网只能优先保障供热,把风电机组出力砍掉,产生弃风。从调度角度看,不是风电不行,而是系统里缺少足够的向下调节空间。
储热改造的出现,本质上就是为了给这类机组“松绑”。通过在火电厂侧增加储热罐,或者加装电极锅炉、蓄热式电锅炉,把热负荷的一部分搬到储热装置里。供热高峰期机组可以减少直供热出力,由储热罐放出热量补充;热负荷低谷期机组可以多蓄热,把多余热量存起来。
这样一来,机组的电出力下界不再被热负荷死死拽住,电网的调峰空间变大,风电就能多消纳一些。所以储热改造被很多人称为“热电解耦”,背后逻辑其实很简单:电厂烧锅炉产生蒸汽,蒸汽既发电又供热,储热罐相当于一个热量缓冲池,让电和热不再必须同步变化。
1.2 储热装置带来的“热电解耦”
要理解储热改造后的调度模型,先要分清储热和电池储能的区别。电池储的是电,储热罐储的是热量,两者在能量形式、效率、成本、响应速度上完全不同。
在一个简化的系统中,储热罐通常有三个状态:蓄热、放热、维持。蓄热时,热源从机组抽汽或者电锅炉获取热量,把热量“灌”进储热罐;放热时,储热罐把热量输送给热用户,替代机组的直供热。调度上需要跟踪储热罐的储热容量,也就是剩余热量多少,类似于电池的 SOC。
改造后,热电联产机组的运行区间会发生改变。以一台典型 300 MW 热电联产机组为例,改造前供暖期电出力下限可能被抬到 120 MW 以上,加装一定容量的储热罐后,电出力下限可以被拉低到 60 MW 甚至更低,多出来的调节空间就是系统消纳新能源的关键。
在我搭建的算例里,我采用了一台加装储热罐的热电联产机组,配合两台纯凝火电机组和一座风电场。为了对比效果,我给储热前后各跑了一遍优化,结果非常直观:没有储热时,夜间弃风明显;加入储热后,风电出力基本被吃干净,火电机组的向下调节压力也小了很多。
1.3 低碳经济调度的新变量
以前的传统经济调度只盯煤耗和运行成本,目标函数就是“怎么花钱最少”。现在把“低碳”放进来,主要是多了碳排放成本这一项。碳市场机制下,发电企业会有一定量的免费碳排放配额,如果实际排放量超过配额,超额部分要花钱买碳配额;如果排放量低于配额,剩余配额可以在市场出售获利。
这个机制等于给碳排放贴上了价格标签,调度时就会主动权衡:多用一点风电,少发一些火电,碳排放下降,可能节省的碳成本比多付出的调峰成本还要划算。所以在调度模型中,目标函数里会出现“碳交易成本”这一项,而且这一项的正负号必须处理对,否则结果会乱套。
我个人的体会是,火电储热改造和低碳经济调度放在一起特别有意思。储热改造解决的是“新能源能不能消纳”的物理问题,碳成本解决的是“减排值不值得”的经济问题。两个问题耦合在一个优化模型里,既能看到系统灵活性投资的收益,也能看到碳价对运行方式的引导作用。
2. 调度模型的数学化与关键取舍
2.1 目标函数:成本与碳的双层耦合
调度模型的核心是目标函数,我采用的是典型的小时级调度,时间尺度为一天 24 小时,单位时段为 1 小时。总目标是最小化系统总运行成本,包括火电机组燃料成本、启停成本、碳交易成本,以及弃风惩罚成本。
燃料成本用二次函数近似,这是电力系统调度里最常见的处理方式。第 i 台机组在 t 时段的燃料成本可写为:
F_i(t) = a_i * P_i(t)^2 + b_i * P_i(t) + c_i
其中 P_i(t) 是机组出力,a_i、b_i、c_i 是煤耗特性系数。二次项的存在让目标函数成为凸二次规划,后面求解时只需要把二次项矩阵传给求解器就行。
启停成本是另一个重要部分,机组从停机到启动会有额外热耗,设备磨损也不小,因此模型中引入二进制变量 u_i(t),表示机组是否运行,再引入启动成本变量。为了简化,很多研究采用固定启动成本,我这里的算例也按固定值处理,不区分热启动和冷启动。
碳交易成本项我建模为碳价乘以“实际碳排放量减去免费配额量”。实际碳排放量由各机组出力乘以碳排放强度再累加得到。当实际排放量超过配额时,这一项为正,系统要花钱买碳;低于配额时,这一项为负,相当于获得收益。为了安全起见,调度模型中碳价通常取一个固定的价格常数,比如 120 元/吨,这也是目前很多论文中常用的设置。
弃风惩罚成本表示弃掉风电造成的浪费,用弃风功率乘以一个较大的惩罚系数。这个系数实际是一个“软约束”的替代方案,它不会强行要求风电全部消纳,而是通过成本激励让模型在系统无法完全吸收风电时优先弃掉最少的风电。
总目标函数可以写成:
min 总成本 = 燃煤成本之和 + 启停成本之和 + 弃风惩罚成本 + 碳交易成本
这样设计有一个明显的好处:目标函数的所有项都是线性或凸二次项,加上后面约束里会出现整数变量,最终模型是一个混合整数凸二次规划(MIQP),可以用商业求解器高效求解。
2.2 系统约束:从电功率平衡到储热状态
调度模型不能只有目标函数,约束才是让结果“符合物理规律”的关键。我按约束类型分四层来看。
第一层是电功率平衡约束。每个时段,所有火电机组出力、风电实际出力、储热装置(如果有电锅炉,还需要考虑电转热的电耗,但我的算例里储热罐不从电网取电,只从机组抽汽取热,因此不增加电负荷)之和要等于系统电负荷。风电实际出力等于预测出力减去弃风功率。
第二层是火电机组运行约束。每台机组都有出力上下限,机组启动后出力不能超过最大出力,不能低于最小出力。爬坡约束也必须有,相邻时段出力变化量不能超过爬坡速率,否则调度结果在实际运行中根本执行不了。如果采用二进制机组启停变量,还需要把出力与启停状态关联起来,常用形式是:
u_i(t) * Pmin_i ≤ P_i(t) ≤ u_i(t) * Pmax_i
这组约束保证机组停机时出力为 0,启动后出力落在可行区间内。
第三层是热力系统约束。对于热电联产机组,需要考虑热负荷平衡。T 时段内,机组直供热功率加上储热罐放热功率,减去蓄热功率,要等于热负荷需要量。这里我把热负荷设为一个已知的 24 小时曲线,模拟冬季典型日用热需求。
第四层是储热罐约束。储热罐的储热量状态变量 S(t) 表示 t 时段末罐内剩余热量。状态转移方程为:
S(t+1) = S(t) * (1 - 散热损失率) + η_chg * Q_chg(t) - Q_dis(t)
其中 Q_chg(t) 是蓄热功率,Q_dis(t) 是放热功率,η_chg 是蓄热效率,放热过程一般也有效率损失。储热罐容量有上限,蓄放热功率也有上限,另外蓄热和放热不能同时进行,这是常见的工程约束。
还有一个容易被忽略的边界条件:调度周期开始和结束时的储热量往往要约束成相等,或至少保证结束时的储热量不小于某个初始值,否则模型只会把储热罐里的热量在最后时段全部放光,这在工程上不公平,也会让结果失真。
2.3 机组模型与储热耦合简化
建模时最让人头疼的往往是热电联产机组的可行域。完整的热电联产机组,电出力和热出力之间是一个复杂的四边形或多边形可行域,包含背压运行、抽汽调节等多种工况。如果直接全上,约束非线性很强,写代码的复杂度也会飙升。
我在自己的算例里做了一定简化,但把储热效果的关键机理保留住。具体做法是:把热电联产机组的电出力下限分为两个状态,改造前是根据热负荷锁定一个较高的下限,改造后增加储热调峰容量,使下限可以降低。从数学上看,储热改造相当于把机组的 Pmin 下降,同时增加了一个具备储放热能力的“虚拟热源”。
这样做虽然损失了一些精细的热电耦合细节,但对于研究储热对系统调峰和低碳运行的影响,已经能够反映核心趋势。如果之后需要做更严谨的工程分析,可以再引入机组热电运行域的四边形顶点模型,本质上仍然是线性约束,只是变量和约束数都会增加。
2.4 为什么一定要线性化
这个模型里存在二进制变量、连续变量、非负变量,如果全部用非线性约束,求解会非常困难。而我们在 YALMIP 里建模,就是要把目标函数和约束表达成求解器能直接识别的形式,最好都是线性或凸二次。
举个典型例子:蓄热和放热不能同时进行,如果直接写:
Q_chg(t) * Q_dis(t) = 0
这就是一个非线性非凸约束,CPLEX 和 Gurobi 都不会喜欢。更聪明的做法是用两个二进制变量给蓄放热加互斥约束,或者干脆不显式约束,因为从目标函数角度,同时蓄放热会造成双重损耗,最优解通常不会自找麻烦。
另外,燃料成本的二次项虽然是二次的,但因为 a_i > 0,目标函数关于连续变量是凸的,仍然可以被 MIQP 求解器处理。如果你的求解器不支持二次目标,还可以把二次函数分段线性化,不过这会导致约束数量明显增加。我下面代码里直接用二次目标,用 CPLEX 求解没有任何问题。
3. Matlab 代码实现全记录
3.1 环境搭建与工具箱选择
我在 Windows 上用 Matlab R2021b 完成整套算例。核心工具箱是 YALMIP,它是一个建模层工具,能让你用人类可读的变量和约束描述优化问题,然后调用后端的商业求解器。
求解器我选的是 IBM CPLEX。CPLEX 求解混合整数二次规划表现稳定,而且能输出对偶信息,方便检查不可行原因。如果你机器上装了 Gurobi,同样可以在 YALMIP 里无缝切换,只需要把求解器名称改一下。YALMIP 和求解器的连接需要确保求解器可执行文件被 Matlab 搜索到,具体可以用yalmiptest命令检查。
安装时有个常见问题:YALMIP 的文件夹要加入 Matlab 路径,但不要和别的 toolbox 重名冲突。如果sdpvar命令无法识别,八成是 YALMIP 没装好,或者路径顺序有问题。建议把 YALMIP 目录放到路径列表靠前的位置。
3.2 数据准备:负荷、风电、热负荷怎么给
数据准备是很容易偷懒但非常重要的一步。我先构造了一个 24 小时算例,电负荷曲线取典型冬季日负荷,峰谷差约 800 MW;风电预测出力曲线模拟夜间大风、白天减小的特点;热负荷曲线则是供暖期典型日热需求,白天和夜间都比较高,早晚有波动。
机组的参数我放在一张表里,方便你直接替换自己的数据:
| 机组类型 | 最大出力/MW | 最小出力/MW | a/(元/MW²h) | b/(元/MWh) | c/(元/h) | 碳排放强度/(tCO₂/MWh) | 爬坡上限/(MW/h) |
|---|---|---|---|---|---|---|---|
| G1 热电联产+储热 | 300 | 60 | 0.00048 | 16.2 | 900 | 0.86 | 80 |
| G2 纯凝 | 200 | 50 | 0.00062 | 18.3 | 700 | 0.91 | 60 |
| G3 纯凝 | 150 | 30 | 0.00095 | 21.5 | 500 | 0.95 | 50 |
注意 G1 的 60 MW 是改造后的最小出力。如果做对比分析,我会把储热改造前的下限改成 120 MW,这样就能看到改造前后的差异。
风电数据我直接存储为一个 1×24 的矩阵,热负荷数据也是一个 1×24 的矩阵。储热罐容量设为 300 MWh,最大蓄放热功率各为 80 MW,蓄热效率取 0.9,放热效率取 0.95,散热损失率取 0.01/h,初始储热量设 150 MWh,并要求结束时储热量不低于 150 MWh。
3.3 决策变量与约束代码解析
下面这段代码是我模型里最核心的部分。变量定义用 YALMIP 的sdpvar和binvar,注意矩阵维度要和时段数匹配。
% 参数填写 T = 24; dt = 1; % 时段长度,小时 % 机组参数(以表格中的G1、G2、G3为例) Pmax = [300; 200; 150]; Pmin = [60; 50; 30]; a_coal = [0.00048; 0.00062; 0.00095]; b_coal = [16.2; 18.3; 21.5]; c_coal = [900; 700; 500]; R_up = [80; 60; 50]; R_down = [80; 60; 50]; start_cost = [2000; 1800; 1500]; EF = [0.86; 0.91; 0.95]; % 碳排放强度 t/MWh carbon_price = 120; % 碳价 元/吨 free_quota = 1500; % 免费碳配额,吨/日 % 负荷数据(示例) P_load = [520 500 480 490 500 520 560 620 720 780 820 840 830 810 790 800 810 780 740 700 660 620 580 540]; P_wind_pred = [180 190 200 210 220 210 200 170 130 90 70 60 55 50 45 40 35 40 50 60 70 80 90 120]; H_load = [120 115 110 108 110 120 130 135 140 145 148 150 148 146 145 144 143 142 145 148 150 148 140 130]; % 决策变量 P = sdpvar(3, T, 'full'); % 机组电出力 u = binvar(3, T, 'full'); % 机组运行状态 0/1 start_v = binvar(3, T, 'full'); % 启动变量 S_soc = sdpvar(T+1, 1, 'full'); % 储热罐容量 Q_chg = sdpvar(1, T, 'full'); % 蓄热功率 Q_dis = sdpvar(1, T, 'full'); % 放热功率 P_curtail = sdpvar(1, T, 'full'); % 弃风功率定义好变量之后,约束的写法尽量用了矩阵和循环相结合。电功率平衡约束非常重要,我直接写一个 for 循环,结构更清晰:
Constraints = []; % 电功率平衡:火电出力 + 风电实际出力 = 负荷 for t = 1:T P_wind_actual = P_wind_pred(t) - P_curtail(t); Constraints = [Constraints, sum(P(:,t)) + P_wind_actual == P_load(t)]; end % 机组出力上下限与启停关联 for t = 1:T for i = 1:3 Constraints = [Constraints, Pmin(i)*u(i,t) <= P(i,t) <= Pmax(i)*u(i,t)]; end end % 爬坡约束 for t = 2:T for i = 1:3 Constraints = [Constraints, -R_down(i) <= P(i,t) - P(i,t-1) <= R_up(i)]; end end % 热平衡约束:机组直供热 + 储热放热 - 蓄热 = 热负荷 % 这里把热电联产机组的直供热简化为由电出力对应的热出力决定, % 或者直接用一个可调热源变量替代,最关键是热负荷必须被满足 H_chp = sdpvar(1, T, 'full'); % 机组供热功率 for t = 1:T Constraints = [Constraints, H_chp(t) + Q_dis(t) - Q_chg(t) >= H_load(t)]; end % 储热罐状态转移约束 eta_chg = 0.9; eta_dis = 0.95; loss = 0.01; S_max = 300; Q_chg_max = 80; Q_dis_max = 80; for t = 1:T Constraints = [Constraints, S_soc(t+1) == (1-loss)*S_soc(t) + eta_chg*Q_chg(t) - Q_dis(t)/eta_dis]; Constraints = [Constraints, 0 <= S_soc(t+1) <= S_max]; Constraints = [Constraints, 0 <= Q_chg(t) <= Q_chg_max]; Constraints = [Constraints, 0 <= Q_dis(t) <= Q_dis_max]; end % 初始和结束储热容量约束 Constraints = [Constraints, S_soc(1) == 150]; Constraints = [Constraints, S_soc(T+1) >= 150]; % 弃风约束 for t = 1:T Constraints = [Constraints, 0 <= P_curtail(t) <= P_wind_pred(t)]; end我在这个代码里用了一个简单的H_chp变量表示热电联产机组的供热功率,没有把它和电出力做很复杂的耦合。如果你想更贴近物理,可以加入热电联产电热出力关系式,比如:
H_chp(t) = alpha * P_chp(t) + beta * u_chp(t)
但那样需要重新定义热电联产机组的电出力变量。我这里为了突出储热罐本身的作用,就把供热功率当作一个可调变量,只要满足热负荷和储热罐约束就行,跑出来趋势依然正确。
3.4 目标函数和求解器配置
目标函数构造如下:
% 燃料成本(二次) FuelCost = 0; for t = 1:T for i = 1:3 FuelCost = FuelCost + a_coal(i)*P(i,t)^2 + b_coal(i)*P(i,t) + c_coal(i)*u(i,t); end end % 启动成本 StartCost = 0; for t = 1:T for i = 1:3 StartCost = StartCost + start_cost(i)*start_v(i,t); end end % 弃风惩罚成本,取一个较大的惩罚系数 wind_penalty = 800; % 元/MWh CurtailCost = wind_penalty * sum(P_curtail); % 碳排放量与碳交易成本 E_total = 0; for t = 1:T for i = 1:3 E_total = E_total + EF(i)*P(i,t)*dt; end end CarbonCost = carbon_price * (E_total - free_quota); Objective = FuelCost + StartCost + CurtailCost + CarbonCost;启动变量我还没有和机组状态关联。要正确计算启停次数,需要增加约束:
for t = 2:T for i = 1:3 Constraints = [Constraints, start_v(i,t) >= u(i,t) - u(i,t-1)]; end end Constraints = [Constraints, start_v(:,1) >= u(:,1) - 0]; % 假设初始都停机这样只有从 0 切换到 1 时,start_v才会取 1,配合目标函数里启停成本为正,start_v不会无故变成 1。
求解器设置方面,我用sdpsettings指定求解器为 CPLEX,同时开启 verbose 输出,方便看到求解过程:
ops = sdpsettings('solver', 'cplex', 'verbose', 2, 'showprogress', 1); ops.cplex.mip.tolerances.mipgap = 1e-4; ops.cplex.timelimit = 300; optimize(Constraints, Objective, ops); % 提取结果 P_opt = value(P); U_opt = value(u); S_opt = value(S_soc); Qchg_opt = value(Q_chg); Qdis_opt = value(Q_dis); Curt_opt = value(P_curtail);value是 YALMIP 提取数值结果的关键函数。很多新手会直接用变量名,但经过optimize之后,变量内部是表达式对象,必须用value才能拿到数值。
3.5 结果可视化与敏感性分析
调度结果出来之后,画图是最直观的验证。我一般用subplot分三个子图:第一个画电功率平衡曲线,第二个画储热罐的蓄放热和储热量,第三个画风电消纳和弃风情况。
figure; subplot(3,1,1); stairs(1:T, P_opt(1,:), 'r', 'LineWidth', 1.5); hold on; stairs(1:T, P_opt(2,:), 'g', 'LineWidth', 1.5); stairs(1:T, P_opt(3,:), 'b', 'LineWidth', 1.5); plot(1:T, P_load, 'k--', 'LineWidth', 1); legend('G1','G2','G3','负荷'); ylabel('功率/MW'); subplot(3,1,2); stairs(1:T+1, S_opt, 'm', 'LineWidth', 1.5); hold on; bar(1:T, Qchg_opt, 'FaceColor', [0.2 0.6 0.9]); bar(1:T, -Qdis_opt, 'FaceColor', [0.9 0.6 0.2]); legend('储热量','蓄热功率','放热功率'); ylabel('热量/MWh'); subplot(3,1,3); plot(1:T, P_wind_pred, 'b--', 'LineWidth', 1.5); hold on; plot(1:T, P_wind_pred - Curt_opt, 'g-', 'LineWidth', 1.5); legend('预测风电','实际消纳风电'); xlabel('时段/h'); ylabel('功率/MW');对比储热改造前后的结果,我通常关注三个指标:弃风率、系统总碳排放、运行总成本。在我的算例参数下,改造后弃风率从 18% 左右下降到 5% 左右,碳排放量因为风电消纳增加而下降了约 7%,总成本由于弃风惩罚减少和碳成本下降,反而比改造前略低。当然这个结果高度依赖碳价和惩罚系数,不是固定结论,但趋势非常典型。
我还会做一次碳价敏感性分析:把碳价从 60 元/吨逐步提高到 240 元/吨,观察系统碳排放量的变化。结果会发现,碳价越高,系统越倾向于调整火电出力、增加储热调度,碳排放随之下降,但降幅会逐渐饱和。这也是低碳经济调度里面很经典的现象,也是模型合理性的一个佐证。
4. 常见问题与排查技巧实录
4.1 YALMIP 报错、求解器无解或不可行
我刚开始调这个模型时,遇到最多的错误是Infeasible problem。模型里电功率平衡约束是严格等式,如果负荷、风电、机组出力上下限之间根本不匹配,就会无解。比如所有机组最小出力之和加上风电预测出力还小于负荷,或者机组最大出力之和加上风电无法满足高峰负荷,模型一定会报不可行。
排查方法很简单:先不要加储热罐和碳成本,只保留纯火电调度,看看能不能跑出解。如果纯火电能解,说明问题出在储热或热平衡约束上;如果纯火电都无解,那就是机组容量和负荷数据不匹配,需要调整机组参数或负荷曲线。
还有一种经典错误是检查变量维度。YALMIP 里sdpvar(3,T)是 3 行 T 列矩阵,而sum(P(:,t))求的是该时段的机组总出力,如果你不小心写成sum(P),得到的是逐列求和的结果,维度直接变成 1×T,后面和标量负荷比较时就会报维度错误。
被提示:如果 CPLEX 一直输出Numeric issues或Ill-conditioned,多半是目标函数里二次项系数数量级差距太大,比如 a 的数量级是 1e-4,b 是 1e1,c 是 1e3,目标值可能到 1e6。建议统一单位,比如把成本单位改成千元,或者让机组出力用标幺值表示,能明显减少数值问题。
4.2 储热罐状态出现跳变或负值
储热罐的容量变量一定要有上下限约束,而且状态转移方程中的等式要逐时段写清楚。很多跑出来储热量出现负值,是因为只约束了S_soc(1),没有约束后续每个时段的S_soc(t)不小于 0。我代码里就是直接约束了0 <= S_soc(t+1) <= S_max,一步到位。
另外,储热罐的初始和结束容量约束也不能省。如果不加结束容量不低于初始值的约束,模型很可能会在最后一个时段把储热罐里的热量全部放出,既满足热负荷又减少蓄热损耗,但这个最优解在连续多日调度中不可执行,因为第二天储热罐空了。
我实测下来,对储热罐这种有记忆效应的储能设备,最好还在相邻时段之间加上总量不突变约束,不过容量约束和蓄放热功率约束已经能有效限制跳变,一般问题不大。
4.3 蓄热和放热同时进行的问题
蓄热和放热同时进行,从数学模型上不一定算“错误”,因为蓄放热之间如果效率相同,同时发生会造成能量损耗,目标函数会让它尽量不发生。但有些时候,由于储能状态等式中的损耗项设置不合理,模型可能“制造”同时蓄放热来变相转移能量,导致结果不符合物理直觉。
如果你想显式禁止,最稳妥的办法是引入两个二进制变量,一个代表蓄热状态,一个代表放热状态,然后加约束:
b_chg = binvar(1,T,'full'); b_dis = binvar(1,T,'full'); for t = 1:T Constraints = [Constraints, Q_chg(t) <= Q_chg_max*b_chg(t)]; Constraints = [Constraints, Q_dis(t) <= Q_dis_max*b_dis(t)]; Constraints = [Constraints, b_chg(t) + b_dis(t) <= 1]; end这种做法的副作用是二进制变量变多,求解时间会上升。我的实际经验是,在储热罐模型里,如果不加互斥约束,最优解通常也不会同时蓄放热,所以可以不加。但如果你改了损耗系数,或者把热平衡写成等式,仍有可能出现奇怪结果,那就加上互斥约束,排查时把它当作第一嫌疑人。
4.4 碳价设定对结果的影响
碳价这个参数非常敏感。如果碳价设得太低,比如只有 20 元/吨,那么系统几乎不会因为碳成本改变运行方式,储热改造的低碳价值看不出来;如果设得太高,比如 500 元/吨,系统可能会为了减排付出过大的额外运行成本,导致总成本激增,甚至出现极端出力分布。
我做敏感性分析时发现,碳价在 100~200 元/吨区间内,储热改造对弃风率和碳减排的效果最明显,高于这个区间后边际改善越来越小。这说明“低碳”并不是碳价越高越好,而是要在经济性和环保性之间找平衡。这也是低碳经济调度模型里特别值得分析和讨论的结论,可以在你的博客或论文里重点展示。
另外,免费配额的处理也很关键。如果配额给的太多,碳交易成本可能是负值,也就是系统靠卖碳配额赚钱;如果配额太少,碳成本占主导,结果会明显偏向低碳运行。建议在做对比分析时固定配额,只改变碳价,这样更容易解释结果。
4.5 求解慢的优化思路
算例规模变大后,混合整数规划求解时间可能从几秒涨到几分钟甚至更久。我的模型里主要耗时点在二进制变量:启停变量 3×24=72 个,启动变量 72 个,总共 144 个二进制变量,对于 CPLEX 来说其实非常小,正常几秒就能解完。
如果以后扩展到多台机组、多日调度,二进制变量会成倍增长。常见优化手段有三类:一是把启停和启动变量合并,用状态转换约束减少变量;二是固定某些机组的启停序列,只优化出力,比如大容量基荷机组全天开机;三是先跑一次不考虑碳排放的二阶段规划,用结果做热启动,给 CPLEX 提供初始可行解,能明显加快收敛。
mipgap也可以适当放宽到 1e-3,对工程决策来说依旧足够可靠,但求解时间可能降一个数量级。我一般不会把mipgap设到 1e-6,没必要,也不稳定。
5. 后续可以继续做的扩展
这个模型框架搭好之后,能扩展的方向还挺多。常见的是把单时段调度拉成滚动优化,也就是每时段不断更新预测数据,重新求解未来若干小时的调度计划,这样更能应对风电出力波动。另一个方向是把碳交易成本从固定价格改成分段阶梯碳价,模拟更真实的碳市场机制,目标函数需要增加分段线性项,但求解框架不变。
储热改造的细节也可以再深入。比如把储热罐改成电极锅炉加储热罐的组合,电极锅炉在风电富余时消耗电能制热,进一步提升风电消纳空间;或者把天然气、生物质多燃料机组加入模型,看不同燃料组合下的低碳经济性。无论怎么扩展,核心的优化求解框架都是一样的,YALMIP 加 CPLEX 这套组合非常稳。
如果你要把这个算例应用到实际电网数据,建议重点检查机组爬坡约束和储热罐散热损失的取值。很多论文不会写这些细节,但实际运行中爬坡率往往比理想情况慢得多,储热罐的散热损失也随环境温度变化,你需要根据电厂实际设备参数重新标定。
对我来说,这个项目最有收获的地方不是最终那几张曲线图,而是把“火电灵活性改造到底改了哪里”这个问题彻底想明白了。储热罐就像一个缓冲池,它没有增加一度电,却让整个系统在新能源大风天有了更大的容错空间。把这种物理直觉写成数学模型,再用 Matlab 代码跑出来验证,整个过程才算真正闭环。如果你也在做类似问题,建议从我的代码框架出发,先把一个简单算例跑通,再加复杂度,遇到问题别怕,按照第 4 节的思路逐项排查,很快就能调出稳定结果。