1. 先搞清楚:低碳调度为什么必须同时考虑这两个因素
1.1 从“成本最小”到“低碳+经济”的转变
搞综合能源系统调度的人都知道,传统经济调度模型的核心只有一件事:在满足负荷需求的前提下让发电成本最低。但当你开始做低碳调度课题,会发现事情没那么简单——目标里多了一个碳排放维度,运行方式里多了碳捕集电厂这类新元素,负荷侧也不再是写死的数据,而是具备主动调节能力的响应资源。
我最早接触这个模拟项目X的时候,需求非常简短:“用MATLAB建一个含碳捕集电厂和需求响应的低碳综合能源系统调度模型。”这话听起来不复杂,真正动手以后才发现,光是“碳捕集电厂到底怎么参与调度”这一个点,就能推翻好几版方案。
核心的认知转变在于:碳捕集电厂不是传统意义上的电源。它在发电的同时还要消耗一部分电能去捕集自己产生的CO2,所以对外提供的净出力不是锅炉燃烧出来的毛出力,而是毛出力减去捕集能耗。这个耦合关系如果做不透,后面所有约束都失真。而需求响应又是一个维度——以前调度把负荷曲线当刚性边界,现在负荷侧通过价格信号或者激励政策,可以按计划平移和削减,源荷之间的互动关系彻底变了。
1.2 整体框架:源-网-荷-储-碳五要素
我把这个模型的整体框架概括为“源-网-荷-储-碳”五要素。源侧有常规燃煤机组、带碳捕集改造的煤电机组、燃气热电联产机组、风电场和光伏;网侧在这个课题里简化为单母线汇集,不展开网络拓扑;荷侧除了刚性电负荷和热负荷之外,新增价格型和激励型两类需求响应资源;储能侧包含电储能电池和储热罐;碳侧则把机组碳排放、捕集量、捕集能耗和净出力耦合在一起。
在这个框架下,调度模型的优化任务是:在24个时段内,给出每台机组出力、碳捕集率、储能充放电功率、需求响应削减量或转移量等所有决策变量的最优值,使得总费用最低——总费用同时包含燃料成本、运行维护成本、启停成本、碳交易成本和需求响应补偿成本——同时还要满足电力平衡、热力平衡、机组爬坡、储能SOC、碳捕集耦合等约束。
这套东西说穿了不神秘:先建立各单元的数学模型,再统一到同一个优化框架里用求解器去算。但真的动手建的时候,细节非常多。这篇文章就把我的建模思路、MATLAB实现过程和踩过的坑完整写出来,给做同类方向的同行一个可以参考的底稿。
2. 碳捕集电厂建模:核心是“捕集能耗耦合”
2.1 捕集单位CO2要付出多少能量代价
碳捕集电厂的物理过程是:燃烧产生的烟气进入吸收塔,通过化学吸收剂(典型如醇胺溶液)吸收CO2,富液被送进再生塔用蒸汽加热解吸。整个过程要消耗大量热能,同时泵和压缩机也需要电能。在调度模型里,完全仿真这个热力过程是不现实的,工程上普遍采用一种简化但有效的近似——把单位捕集能耗折算成统一的电功率损耗。
具体到我建的这个模型,设机组i在时段t的毛出力为P(i,t),那么这个时段机组产生的CO2总量为:
E(i,t) = e_i × P(i,t)
其中e_i是机组单位发电量对应的CO2排放强度,单位是t/MWh。我的算例里燃煤机组取0.8左右,燃气机组取0.4左右,Hot这部分数值可以根据煤质和燃气热值标定。
碳捕集系统的捕集量记为:
C(i,t) = ε(i,t) × E(i,t)
ε(i,t)是捕集率,这是调度模型的决策变量。捕集率不能乱设,运行范围一般在0.6到0.95之间,太低没有减碳效果,太高再生能耗剧烈上升。捕集能耗按下式计算:
D(i,t) = λ_i × C(i,t)
λ_i是单位捕集能耗,典型值在0.2-0.4 MWh/tCO2之间。我的是取0.3,这个数值直接影响净出力计算。
于是碳捕集机组的净出力为:
P_net(i,t) = P(i,t) − D(i,t)
把上面的式子代进去就会发现一个重要事实:
P_net(i,t) = P(i,t) × (1 − λ_i × ε(i,t) × e_i)
也就是说,净出力等于毛出力乘上一个与捕集率相关的系数。捕集率越高,净出力损失越大。如果捕集率是决策变量,这个式子就引入了决策变量的乘积,直接使用会让模型出现非线性。
2.2 捕集率调节范围与工程约束
捕集率连续可调听起来很理想,实际工程中再生塔的负荷调整是缓慢的,相邻时段捕集率不能剧烈跳变。我在模型里加了这么一条爬坡约束:
|ε(i,t) − ε(i,t−1)| ≤ Δε_max
一般Δε_max取0.1左右。加这个约束在初学时容易被忽略,但很重要——否则求解器会让捕集率在相邻时段从0.9直接跳变到0.6,数值上可行,物理上根本做不到。
处理捕集率与毛出力乘积的非线性,我在项目X里实际用了两种方案。
第一种方案是把捕集率离散成几个档位,例如0.6、0.7、0.8、0.9四档,用一组二进制变量选择档位。这样P_net(i,t)和P(i,t)的关系在每个档位内是线性关系,整体变成混合整数线性约束,交给MILP求解器很友好。
第二种方案更简单,第一版实现时我用的是它——把捕集率固定为常数ε=0.85,只保留“这台机组是否开启碳捕集”这个0-1变量。先把整体流程跑通,确认调度框架没有问题,再回头把捕集率档位放进去。
另外还有一个容量问题需要注意:捕集系统本身有处理上限。即使机组满发产生大量CO2,捕集系统也未必能全部吸收。所以要加这个约束:
C(i,t) ≤ C_cap_i
即捕集量不能超过捕集装置的设计处理能力,在实际工程里这是很常见却被很多模型忽略的限制。
3. 需求响应建模:让负荷侧也具备调节能力
3.1 价格型需求响应:弹性矩阵怎么设置
需求响应分两类,一类是价格型,一类是激励型。价格型DR的核心是用户根据电价变化自主调整用电行为,峰时少用、谷时多用。建模最常用的工具是需求价格弹性矩阵。
自弹性描述本时段电价变化对本时段负荷的影响,交叉弹性描述本时段电价变化对其他时段负荷的影响。公式可以写成:
ΔP_load(t) = P_base(t) × Σ E(t,τ) × (ρ(τ) − ρ_base(τ)) / ρ_base(τ)
其中E(t,τ)是弹性矩阵元素。自弹性一般为负,取值在−0.2到−0.5之间;交叉弹性一般为正,取值在0.1到0.3之间。设置完弹性矩阵后,负荷形状会因为电价变化而变化,从而起到削峰填谷效果。
不过我这里要坦白说一句:价格型DR建模虽然被论文写得多,真正落地到MATLAB优化模型里,它和峰谷电价、弹性系数之间是很强的耦合关系,参数选不好容易失真。我在项目X里做第一版时,先实现了激励型DR,因为它的约束结构更简单,效果也更直观可控——直接告诉用户可以削减多少负荷,每削减一兆瓦给多少补偿,边界清清楚楚,也方便做线性化。
3.2 激励型需求响应:可削减与可转移负荷
激励型DR在我的模型里有两种形式:可削减负荷和可转移负荷。
可削减负荷,就是系统在高峰时段通过协议直接切除一部分非关键负荷,比如商业空调、工业辅助设备。建模时引入削减量δ_down(i,t),满足:
0 ≤ δ_down(t) ≤ δ_down_max(t)
其中δ_down_max(t)是时段t可削减负荷的上限,通常取基础负荷的5%-10%。每次削减,系统需要支付补偿费用,这要进目标函数,否则模型会为了省钱无节制地削减,显然不合理。
可转移负荷更复杂一点,它对应那些可以在时间上灵活安排的负荷,比如电动汽车充电、洗衣机、热水器等。这类负荷的特点是总用电量不变,只是在时间轴上移动。建模时引入“转出量”和“转入量”两个非负变量,它们满足:
转出总量 = 转入总量
这两个条件保证转移前后总用电量守恒。我最开始写的时候漏了总量守恒约束,结果优化出来的结果用电量凭空增加了很多,结果对不上。
不管哪类DR,都别忘了加上响应后的负荷不能为负这条硬约束。也要考虑用户侧的舒适度限制——削减量太大会影响正常使用,因此可以再设一个满意度约束,比如单时段最大削减比例不能超过某个阈值。我在算例里把可削减比例上限设在8%,意思是峰时最多削掉负荷的8%,这个数值通常来自历史负荷分析或者合同条款。
4. 调度优化模型的完整数学表达
4.1 目标函数:成本和碳排放怎么放进同一句话里
低碳调度的目标函数不能只看燃料成本,也不能只盯碳排放,现实中大多是通过碳交易机制把碳排放量折算成经济成本。这样做有个好处:模型始终是单目标优化,求解难度小,又能在公平的比较框架下评估碳排放的经济价值。
我的目标函数包含五部分,写成:
min F = C_fuel + C_om + C_startup + C_co2 + C_dr
C_fuel是燃料成本,燃煤和燃气分别乘上各自的耗量系数和单价。C_om是运行维护成本,通常用单位出力成本系数乘以出力来近似。C_startup是机组启停成本,每台机组每次启动要花一笔固定费用,这一项会让模型自动合理安排启停,避免频繁启停带来的经济损耗。C_co2是碳交易成本,计算公式是:
C_co2 = λ_price × (实际总碳排放 − 免费配额)
实际总碳排放包含所有火电机组的排放;免费配额是给定参数。如果排放超过配额就要买碳配额,反之可以卖,这就是“碳价信号”。当碳价提高时,模型会自动增加碳捕集系统的使用率,也会更多调用可再生机组。C_dr是需求响应补偿成本,包括可削减负荷补偿和可转移负荷补偿。补偿价格比正常电价略低,这样模型是在“不得不削减”时才削减,而不是把DR当免费手段滥用。
4.2 约束体系:从功率平衡到碳捕集耦合
约束体系是模型的地基,我把它们分成五组逐一来写。
第一组是系统功率平衡。电力平衡方程如下:
Σ P_net(i,t) + P_wind(t) + P_pv(t) + P_bat_discharge(t) + P_CHP_e(t) + δ_in(t) = P_load_base(t) + P_bat_charge(t) + P_electric_boiler(t) + δ_out(t) + δ_down(t) + P_pump_CCS(t)
这里要特别提醒:等式左边的输入要带净出力,同时负荷右侧要加上电锅炉耗电、储能充电功率、转出负荷和削减负荷。还有P_pump_CCS(t)——碳捕集系统的泵和压缩机额定耗电。因为综合能源系统里还有热负荷,需要补充热力平衡方程,燃气热电联产机组的热出力和电锅炉产热、储热罐放热共同供给热负荷,这个平衡约束和电力平衡并列。
第二组是常规机组约束。每台火电机组有功出力要落在上下限之间,相邻时段爬坡速率要受限制,启停时间和出力相关也可以加,但第一版我把启停逻辑简化了,只保留了基本的出力上下限和爬坡约束。
第三组是碳捕集机组的耦合约束,就是前面第二节里的那组公式。如果采用捕集率档位离散方案,还需要把档位选择变量和捕集率关联起来。
第四组是储能约束。电储能的SOC更新方程为:
SOC(t+1) = SOC(t) + η_ch × P_ch(t) − P_dis(t) / η_dis
SOC要限制在0.1-0.9之间,充放电功率分别有上限,而且同一时段不能同时充放电——这个用两个二进制变量和互斥约束来实现。储热罐的建模类似,只是热储能损耗略高一点。
第五组是需求响应约束,前面第三节已经写了,可削减量上限、可转移量上下限、转移总量守恒、响应后负荷非负。
最后还需要考虑风电和光伏出力上限。可再生机组在调度模型一般可以简化为可弃风可弃光——出力可以在预测值以内任意取值,这样模型在负荷低谷时可以选择弃掉一部分可再生能源,换取系统更好的经济性,这个和实际风电场的运行方式也对得上。
5. MATLAB代码实现:从建模到求解的详细流程
5.1 求解器选型与YALMIP配置
MATLAB建调度模型,最常用的方式不是自己写求解算法,而是用YALMIP加一个商业求解器。YALMIP是一个建模工具箱,它的价值在于把优化问题描述语言翻译成求解器能识别的标准输入。求解器我选的是Gurobi,因为MILP问题它的求解速度在同类里面非常靠前。
安装好工具箱之后,调用方式很简单,核心代码框架如下:
% 定义决策变量 P = sdpvar(n_gen, T, 'full'); % 各机组出力 eps_ccs = sdpvar(1, T, 'full'); % 碳捕集率 u = binvar(n_gen, T, 'full'); % 机组启停状态 delta_down = sdpvar(T, 1, 'full'); % 需求响应削减量 P_ch = sdpvar(1, T, 'full'); % 储能充电功率 P_dis = sdpvar(1, T, 'full'); % 储能放电功率 SOC = sdpvar(1, T+1, 'full'); % 储能荷电状态 % 目标函数 F = sum(sum(fuel_cost_coeff .* P)) + ... co2_price * sum(E_total) + ... dr_price * sum(delta_down); % 具体每一项要根据参数展开 % 求解器设置 ops = sdpsettings('solver', 'gurobi', 'verbose', 2, 'gurobi.TimeLimit', 300); sol = optimize(constraints, F, ops);这里有一个小细节:YALMIP的sdpvar和binvar定义变量时,第三单引号里面的'full'不是必须的,但显式写出来能让维度更清楚。我做项目的时候遇到过因为维度没对齐导致约束写错而检查半天的情况,后面统一用'full'就没再出问题。
5.2 碳捕集耦合约束代码实现
碳捕集机组的耦合约束,是整个模型里面最有技术含量的部分。下面是捕集率离散档位的实现方式。假设碳捕集机组索引是1,我们把捕集率分为三档,每档对应一个二进制变量z1、z2、z3。
% 碳捕集机组耦合约束 % 排放总量 = 排放强度 * 毛出力 E_ccs = e_coal * P(1, :); % 1 x T 行向量 % 捕集量 = 捕集率 * 排放总量 C_captured = sum(eps_levels .* z, 1); % 1 x T % 其中 z 是 3 x T 的二进制矩阵,eps_levels 是 3x1 向量 % 捕集能耗 D_ccs = lambda_ccs * C_captured; % 净出力 = 毛出力 - 捕集能耗 P_net_ccs = P(1, :) - D_ccs; % 净出力上下限约束 constraints = [constraints, P_net_ccs >= P_min_ccs * sum(z, 1)]; constraints = [constraints, P_net_ccs <= P_max_ccs * sum(z, 1)]; % 每个时段只能选一个捕集率档位 constraints = [constraints, sum(z, 1) == ones(1, T)];这里有个容易出错的点:eps_levels .* z是向量和矩阵逐列相乘,我一开始写成一行一行加,总是报维度错误。后来统一用矩阵乘法和sum(...,1)做,代码简洁也不容易错。
捕集率爬坡约束再加一条:
constraints = [constraints, abs(eps_ccs(:, 2:T) - eps_ccs(:, 1:T-1)) <= 0.1];注意如果用了档位离散化,eps_ccs实际上不是连续变量,而是由z二进制变量计算出来的,这句爬坡约束也要改成对sum(eps_levels .* z, 1)做差分。
5.3 需求响应约束与场景生成
需求响应的约束实现相对简单,但需要额外引入转移变量。我定义的转出量为P_out_t、转入量为P_in_t,它们都是T维非负变量。
% 可削减负荷上限 constraints = [constraints, 0 <= delta_down <= delta_down_max]; % 可转移负荷:转出量和转入量 constraints = [constraints, P_out_t >= 0, P_out_t <= P_trans_max]; constraints = [constraints, P_in_t >= 0, P_in_t <= P_trans_max]; % 转移总量守恒 constraints = [constraints, sum(P_out_t) == sum(P_in_t)]; % 延迟约束:某些时段不能转入 % 比如凌晨负荷低,不允许再转入增加负担 constraints = [constraints, P_in_t(1:6) == 0];最后一个P_in_t(1:6) == 0是我根据实际场景加的,凌晨时段本身负荷低,这时候不允许负荷再转入,否则谷电时段又被填出一个小峰,反而违背了削峰填谷的初衷。这个细节在实际项目中很有用。
初始场景数据要处理好,我写了一个数据加载脚本,把24小时的负荷数据、风电预测、光伏预测、电价、热负荷等参数汇总成一个结构体,然后统一传给优化模型。这样后面做敏感性分析、改参数,都不用动模型部分的代码,只改数据文件就行。
6. 仿真算例与效果对比
6.1 算例系统设计与参数设置
为了验证模型,我设计了一个小型的综合能源测试系统:3台火电机组,其中机组1是煤电并改造为碳捕集,机组2是普通煤电,机组3是燃气热电联产机组,同时加上1座风电场、1个电储能装置和1套电锅炉加储热罐。负荷曲线取自某个典型冬季日的公开测试数据,风电出力按“预测值可弃”的方式接入系统。
碳交易价格我设置了多个档位做敏感性分析,基准值取40元/tCO2。碳配额按历史排放的80%发放,如果模型总排放低于配额,富余配额可以按碳价出售冲抵成本。需求响应补偿价格定为每MWh 120元,比正常峰时电价低一点,这样模型不会滥用DR。这些参数我都是按行业平均水平取的,不同的价格组合会影响结果,但整体趋势一致。
6.2 结果对比与分析
我跑了三个场景做对比:场景A是无碳捕集且无需求响应,相当于传统经济调度,作为基准;场景B加入碳捕集电厂但无需求响应;场景C同时加入碳捕集电厂和需求响应。
在我的参数设置下,场景B相比场景A,系统总碳排放量下降约14%,同时运行总成本小幅上升约6%。这一部分成本增加主要来自碳捕集能耗和捕集装置的运行维护。这还是碳价40元/t时的情况,碳价翻倍后,碳捕集电厂被调用得更多,碳排还会进一步降低,但总成本上升也更明显。
场景C在场景B的基础上加入需求响应后,碳排放量又降低了约5%-7%,并且最关键的指标——系统最大净负荷——降低了约9%。峰谷差缩小后,机组2和机组3不再需要频繁启停,系统启停成本下降,可再生能源消纳率反而提升了。这里面的机制值得说一句:需求响应削减了高峰负荷,碳捕集电厂在高峰时段就不用被迫满发,净出力裕度更大,整个系统更有余量去接纳风电。
有意思的是我对比发现,需求响应和碳捕集并不是独立起作用的。需求响应削峰之后,碳捕集设备的捕集率可以维持在较高水平,因为不需要为了短时顶峰而压捕集率,这是两个技术叠加后的联动效应。模型中把这两者的耦合关系建进去,才能看到这种效果;如果只是在两个独立模型里分别算,这个结论是得不出来的。
7. 调试过程中踩过的坑
7.1 维度问题:最隐蔽也最耗时
YALMIP建模最大的坑就是维度不匹配。我第一版跑出来的功率平衡值对不上,排查了两天才发现是sdpvar(1, T, 'full')和sdpvar(T, 1, 'full')混用了,行向量和列向量在矩阵运算里出错,导致某一行约束被错误地赋值。
后来的经验是:所有变量全部统一成行向量或者全部统一成列向量,不做混用。我统一用行向量,sdpvar(1, T, 'full'),然后所有约束和矩阵运算都按行向量来写,再也没有出现过类似问题。还有一招是写完每个约束后,用size(constraints(end))快速检查约束矩阵尺寸是否符合预期。
7.2 求解无解时的排查路径
优化无解是另一个高频问题。跑sol = optimize(...)后返回sol.problem == 1,最开始我第一反应是约束写错,但其实问题往往出在数据上。我的排查顺序是这样的:
先检查所有平衡约束的等式是否真的平衡——把各组件的出力求和,和负荷对比。很多无解是因为负荷数据里包含了DR削减量,但我在平衡方程里又减了一次,重复计算。然后检查SOC更新和储热罐的动态方程,有没有出现初始SOC和终止SOC矛盾的情况。最后检查碳捕集耦合约束里有没有把P_net_ccs和某个常数上限比较,导致净出力上限小于下限的荒唐情况。
我的建议是,当模型无解时先把碳捕集率和DR变量固定为常数,比如固定ε=0.85、δ=0,看一个纯火电加储能的调度模型能不能解。如果能,说明问题在耦合部分;如果还不能,就是基础约束的问题。
7.3 模型合理性怎么验证
模型能求解不等于结果合理。我做合理性验证有几条经验:一是把碳价设为0、DR补偿设为极大值,结果应该退化为几乎不用DR的常规调度模型;二是把风电预测设为0,系统必须完全靠火电和储能满足负荷,看是否出现违反能量守恒的结果;三是调高碳价,模型应该逐渐增加碳捕集率,减少碳排放,如果趋势反了,多半是目标函数里符号写反了。
还有一个最简单有效的验证方法:把优化后的机组出力、负荷平衡、SOC变化曲线画出来,逐时段人工核对。比如负荷曲线在某个时段大幅下降,储能电池在前后时段是否出现了合理的充放电变化。眼睛对齐图形和数据往往能发现公式里看不出来的问题。
8. 一点真实的个人体会
这个项目做下来,我最大的体会是:低碳调度模型真正难的不是数学公式,而是把设备特性和工程约束用正确的方式翻译成优化语言。碳捕集电厂不是简单的“出力减去能耗”,需求响应不是简单的“峰值削减”,它们之间的互动有非常强的耦合效应。把模型里的这些关系处理好,得到的调度结果才能真的指导实际运行。
另外也想对刚开始做这类课题的同行说一句:先跑通一个简单版本,别一上来就追求面面俱到。我第一次做的时候想把所有细节都建模进去,光是碳捕集系统的储液罐就设计了六种运行模式,结果约束数量翻倍,求解器计算时间从几十秒涨到十几分钟,而且数值稳定性变得很差。后来我把模型砍到最精简的版本,把核心逻辑验证清楚,再一层一层往里面加功能,后面每一步都能控制住质量。
如果你准备在自己的论文或者项目里做类似的方向,建议从最基础的“火电+风电+储能”调度模型起步,跑通之后再逐步加入碳捕集、需求响应、多能互补这些模块。每加一个模块就做一次完整性检查,别一口气堆上去。这个顺序看起来慢,但踩坑最少,代码质量反而最高。最后把碳捕集率档位、DR补偿价格这些参数都做成可配置项,你会发现在后面做敏感性分析时,这样改造的收益远超预期。