做电力系统调度的朋友这两年应该没少被“火电灵活性改造”和“双碳目标”这两个词轮番轰炸。火电机组储热改造通俗点说就是给传统火电机组加装一个巨大的“热水瓶”——在新能源大发、电力过剩的时候把多余的电能转换成热能存起来等到晚高峰电力紧张的时候再把储存的热能释放出来发电或供热。这样一来火电机组就不再是那个“要么满发、要么停机”的笨重大块头而是能灵活调节、配合新能源消纳的调节性电源。我最近正好在Matlab里把这套“储热改造”的模型完整跑通结合低碳经济调度目标做了优化今天把整个思路、模型搭建过程、以及我用Matlab实现时的经验和踩过的坑一次性整理出来分享给正在做相关课题的研究生或者刚入行的工程师朋友。这篇文章不仅适合做纯理论研究的读者也适合那些想把模型落地成可运行代码的实操派。看完之后你至少能收获一套明确的目标函数和约束条件的数学表达、一份能直接参考的Matlab代码框架以及几个在结果调试阶段非常容易卡壳的问题排查思路。1. 为什么储热改造能同时兼顾“低碳”和“经济”先说一个很多人容易混淆的概念。提到火电灵活性改造很多人第一反应是“降低最小出力”也就是让机组在低负荷下也能稳定运行。但单纯压负荷是有物理极限的锅炉的稳燃条件、脱硝系统的烟温窗口都限制了机组不能无限往下压。储热改造的核心思路则完全不同它不是在“发电”这个环节上硬扛而是把能量在“电”和“热”两个形态之间做一个缓冲和转移。1.1 火电厂的“热水瓶”是如何参与调度的具体的工程实现上主流的路线有两种一种是抽汽蓄热在机组运行期间从汽轮机中间级抽出一部分蒸汽加热储热介质导热油、熔盐或者高压热水把热能储存在罐体里另一种是电极锅炉直接用多余的电能加热水或熔盐。前者的能量转化效率高抽汽本身是热能的直接转移后者在“深度调峰”场景下更常用——因为电极锅炉完全可以把机组在低谷时段变成一个大功率用电设备硬生生在厂内多造出一块负荷。两种方式殊途同归都会在系统层面留下一个关键特征火电机组的功率可以被“分时转移”。低谷时段多发电存起来或者用电制热存起来高峰时段释放。对于调度模型而言我们关心的不是罐体内部怎么对流换热而是它在一天24小时里每个时段充了多少热、放了多少热、罐内剩余多少热。这就是抽象的“储热罐模型”。1.2 低碳经济调度到底在优化什么传统的经济调度目标是让系统总发电成本最小本质上是燃料成本的最小化。有了碳交易机制之后成本函数里多了一项碳配额成本——机组排放多少二氧化碳就得购买相应的碳排放配额。这样原来“烧煤便宜”的算计就变了如果一个机组碳排放强度高它的综合运行成本会显著上升低碳机组和可再生能源在调度排序中的优先级就会自然提前。所以这里的“低碳”和“经济”不是两个割裂的目标而是通过碳价这个信号把碳排放量“货币化”之后统一纳入经济调度这一个目标函数里。这就是低碳经济调度的核心思想以含碳成本的总运行成本为目标以系统安全稳定运行为约束求解各机组的出力计划。对于储热改造机组我们还要额外考虑储热罐的充放热计划让它在碳约束下发挥最大的调峰和减排价值。2. 数学模型的完整搭建目标函数与三层约束体系这一部分是整篇博文的灵魂。模型搭得清不清晰直接决定后面Matlab代码是“一次跑通”还是“改了三天还在报错”。我在建模时把它拆成三块目标函数的构建、常规电力系统约束、储热改造特有的运行约束。2.1 目标函数燃料成本、碳成本、弃风惩罚的加权博弈低碳经济调度的目标函数可以写成min F Σ(F_fuel F_carbon) F_wind_penalty逐项拆开来说。第一项是燃料成本也就是煤耗成本。工程上通常用二次函数近似描述机组出力P与煤耗量之间的关系F_fuel a × P² b × P c其中a、b、c是煤耗特性系数。这个二次项在后面的Matlab求解时会被线性化处理一会展开讲。第二项是碳成本。引入碳交易机制后每台机组根据实际排放量M和免费配额E_free的差值产生购买或出售配额的成本。排放量M和出力之间存在线性关系M δ × P这里的δ就是机组的碳排放强度系数t/MWh。碳成本的计算式是F_carbon λ × (M - E_free)如果实际排放低于配额括号内是负值意味可以卖配额赚钱这就是低碳机组的经济激励。第三项是弃风惩罚。我习惯在目标函数里加一个单位弃风惩罚成本如果不加模型可能出现“为了省钱而大量弃风”的极端结果——这在碳价较低、煤价也比较低的时候经常发生。加了惩罚之后模型会在“减少燃煤”和“多消纳新能源”之间自动找平衡。2.2 功率平衡与机组运行约束最容易被忽略的细节这一类约束属于电力系统调度模型的标配但我还是要强调几个容易出错的点。首先是功率平衡约束。电力系统在任何时刻都必须保持发电和负荷的实时平衡即所有机组出力之和等于负荷需求。这个约束我在建模初期犯过一个低级错误——忘了把储热系统的耗电功率加进负荷侧。加了电极锅炉之后厂内负荷增加功率平衡的等号右边不是简单的系统负荷而是“系统负荷储热系统用电功率-储热放热供出的电功率”如果储热系统是热电联产形式则还需要耦合热负荷。其次是机组出力上下限约束。对于改造后的机组出力的上限不变但下限值得仔细讨论。纯粹通过降低最小出力改造的机组下限可以从50%额定降到30%甚至20%但如果是靠储热实现调峰的机组下限其实并没有降低多少——它是在低负荷时段把多出来的电转成热蓄起来等效上相当于“用了更少的电、但还是烧着同样多的煤”。这一点在模型里要小心我见过有人直接在代码里把储热机组的最小出力改成了20%结果目标函数算出一个非常“乐观”却不现实的结果。还有爬坡约束。机组在相邻时段之间的出力变化不能超过爬坡速率限制。储热系统的充放热切换会导致机组电出力变化率变大如果爬坡约束设置不当求解器会很难收敛。2.3 储热罐模型连续变量 状态耦合约束储热罐是本模型的特殊之处也是Matlab代码中最容易出bug的地方。它的状态可以用四个变量来描述SOC(t)t时段末罐内剩余热量MWh对应储能系统的荷电状态H_in(t)t时段的充热功率MWH_out(t)t时段的放热功率MWS(t)0-1变量表示储热系统整体处于充热还是放热状态用于防止同时充放核心的状态转移方程是SOC(t) SOC(t-1) (H_in × η_in - H_out / η_out) × Δt其中η_in是充热效率、η_out是放热效率。注意一个关键细节放热侧的效率放在分母上。也就是说如果想要在放热侧输出1MWh的热量储热罐内部其实要释放出1/η_out的热量这是工程上的定义差异直接写代码时容易搞反。同时还要约束充放热功率上限和罐容上限0 ≤ H_in(t) ≤ H_in_max × S(t) 0 ≤ H_out(t) ≤ H_out_max × (1 - S(t)) SOC_min ≤ SOC(t) ≤ SOC_maxS(t)是0-1整数变量它的作用是禁止储热罐在同一个时段既充又放——虽然现实中这种操作没有物理意义但如果不加这个变量约束求解器就可能在数学上做出“既充又放”的荒谬结果来凑最优值。2.4 碳排放约束配额与阶梯碳价的建模碳交易的建模方式直接影响目标函数的线性化难度。最基础的做法是固定碳价碳成本是线性的但实际碳市场中配额价格往往是阶梯式的——超过一定排放量之后超出部分价格更高即阶梯碳价。在Matlab实现里我用了一个分段线性化的方法来处理这个阶梯函数代价是增加了几组辅助变量和额外的约束但换来了更贴合实际的价格信号。阶梯碳价函数可以表达为基于同一个Class同一类型or同一默契设定下它是分段函数排放量M在配额区间内价格为λ1超过配额后价格为λ2λ2 λ1。我把它写成等价的两段线性函数通过一组0-1变量约束每一段的选择这在Matlabyalmip下的实现比较直观我会在代码部分展示。3. 求解流程与混合整数线性规划(MILP)问题的转化的完整思路模型搭建完成之后不要急着把所有约束typeset进Matlab就开始求解。我需要先问自己一个问题这个模型是线性规划(LP)、非线性规划(NLP)还是混合整数线性规划(MILP)选错了求解器后面所有工作都是白费。3.1 非线性项的处理方法与理论依据在我搭建的目标函数中燃料成本a×P²是一个典型的二次项属于非线性规划。但是如前面所说实际求解时我将其分段线性化处理——把机组出力的可行区间分成若干段每一段用一个线性函数近似引入对应的0-1变量。这样一来目标函数和约束全部变成线性形式可以用成熟的混合整数线性规划求解器高效求解。为什么不直接保留二次项用非线性求解器呢因为电力系统调度模型动辄几十台机组乘上24个时段非线性求解器的求解时间可能呈指数增长收敛性和全局最优性都没有保证。MILP的求解器经过几十年的发展分支定界算法已经非常成熟在几秒到几分钟内就可以获得全局最优解。实践下来分段线性化造成的精度损失非常小完全在工程接受范围内。3.2 为什么阶梯碳价让问题变成了MILP阶梯碳价中引入了0-1变量选择“处于哪个碳价段”这就是混合整数变量的来源。储能系统的充电/放电状态S(t)是另一个0-1变量。这些整数变量的引入把原本的线性规划问题升维成了MILP问题。MILP的求解时间随整数变量数量上升而增长所以在编写Matlab代码时一个非常关键的优化习惯是尽量减少不必要的整数变量。能通过连续变量边界约束解决的就不要引入0-1变量。3.3 求解器的选型与Matlab环境配置建议常见的Matlab求解器搭配有两种。一种是Matlab自带的linprog和intlinprog前者解决线性规划后者解决混合整数线性规划另一种是Matlabyalmip工具箱再外接cplex或者gurobi。我个人强烈建议使用yalmip作为建模语言哪怕你最终用intlinprog求解yalmip的符号建模方式也能大幅减少调试时间——它的约束表达几乎和数学公式一一对应代码可读性比纯粹用矩阵拼装高出两个量级。求解器方面cplex和gurobi的求解性能在MILP方向上比intlinprog好很多。对于小规模24时段、几台机组的问题intlinprog勉强够用但如果要做24小时×数十台机组或者进行多场景分析还是建议用cplex或者gurobi。4. Matlab代码实现的架构设计与关键模块解读终于到了代码部分。我按照“数据准备–模型搭建–求解–结果可视化”四个步骤来组织代码下面的内容都是基于这个框架展开的。4.1 数据准备环节的细节参数与数据结构为了方便修改和测试不同场景我用一个脚本文件load_data.m来统一管理所有基础数据。这里面包含系统负荷曲线24个时段的有功负荷需求MW风电预测出力曲线24个时段的可用风电出力MW火电机组参数额定容量、技术最小出力、煤耗系数、爬坡速率、碳排放强度储热罐参数罐容上限、充热效率、放热效率、最大充/放热功率、初始SOC碳交易参数碳价、免费配额比例或者绝对配额量数据准备阶段最容易被忽略的是各参数的物理单位一致性。我就吃过亏负荷和机组出力单位用MW但储热罐容量用了MWh在状态转移方程里计算SOC时需要乘以Δt单位小时才能把功率换算成能量。有时候一不留神单位混了模型结果看起来正常但数值解读完全错误。4.2 基于yalmip的约束构建核心代码下面是模型搭建的核心框架我用yalmip语法直接声明变量并添加约束代码逻辑跟2.3小节的数学公式一一对应。% 时间段定义 T 24; dt 1; % 每个时段时长小时 % 变量声明 P sdpvar(1, T); % 火电机组出力 P_wind sdpvar(1, T); % 实际风电消纳功率 P_wind_curtail sdpvar(1, T); % 弃风功率 H_in sdpvar(1, T); % 储热充热功率 H_out sdpvar(1, T); % 储热放热功率 SOC sdpvar(1, T); % 储热罐状态 S_state binvar(1, T); % 充/放状态 % 功率平衡约束注意储热系统耗电体现在负荷侧 Constraints []; for t 1:T Constraints [Constraints, P(t) P_wind(t) Load(t) H_in(t)]; Constraints [Constraints, P_wind(t) P_wind_curtail(t) P_wind_forecast(t)]; end % 机组出力上下限约束 Constraints [Constraints, P_min P P_max]; % 储热罐状态转移约束 SOC(1) SOC_initial (H_in(1) * eta_in - H_out(1) / eta_out) * dt; for t 2:T Constraints [Constraints, ... SOC(t) SOC(t-1) (H_in(t) * eta_in - H_out(t) / eta_out) * dt]; end % 储热罐容量约束 Constraints [Constraints, SOC_min SOC SOC_max]; % 充放热功率约束含互斥状态 Constraints [Constraints, 0 H_in H_in_max * S_state]; Constraints [Constraints, 0 H_out H_out_max * (1 - S_state)]; % 目标函数燃料成本 碳成本 弃风惩罚 Objective sum(fuel_coeff_a * P.^2 fuel_coeff_b * P fuel_coeff_c) ... carbon_price * sum(carbon_intensity * P - carbon_quota) ... wind_penalty * sum(P_wind_curtail);第一版跑通之后我需要额外说明一点关于燃料成本的二次项处理。在上面的Objective中直接写了P.²这个在yalmip中是二次规划(QP)的问题描述方式在某些求解器下可以被接受。但为了让模型完全线性化、能够用纯MILP求解器处理我实际代码里用的是对P分段线性化的做法——把P的可行区间分成K段用sdpvar的二进制辅助变量组实现每段线性近似。这个改动虽然增加了代码量但求解稳定性明显更好。4.3 结果可视化与关键图表输出求解完成后我用Matlab自带的plot命令绘制几类关键图表机组出力曲线与负荷曲线的对比直观看到储热低谷时段“增负荷”的效果储热罐SOC随时间变化曲线检查能量是否守恒、是否超出罐容上限弃风率柱状图对比改造前后的消纳效果成本构成堆叠图燃料/碳/弃风惩罚展示低碳改造带来的成本结构变化可视化不只是为了写论文插图更是排查模型错误的重要工具。我记得第一次跑出来的SOC曲线竟然在凌晨2点和凌晨6点之间出现了突然掉到负数的情况排查了半天才意识到是初始SOC设置和约束顺序导致的状态转移方程逻辑错误。类似这种问题只盯着一堆数字很难感觉到异常画成图一眼就看出来了。5. 实际算例测试从仿真结果看储热改造的价值光有模型和代码还不够我用一个简化的算例来验证模型的有效性。算例规模不大但是足够说明储热改造在低碳经济调度中的核心价值。5.1 算例场景与参数设定我设定了一个含一台600MW火电机组、200MW风电场的简化系统。机组额定出力600MW技术最小出力不含储热时为300MW煤耗系数取600MW机组的典型值碳排放强度0.8t/MWh碳价取50元/吨。储热系统容量设定为300MWh最大充热功率100MW充热效率0.9放热效率0.85。负荷曲线采用典型冬季日负荷曲线白天有两个高峰午高峰和晚高峰夜间出现低谷约350MW。风电出力设定为夜间大发的典型反调峰特性——夜间风电出力高、负荷低这恰恰是弃风最严重的时段。5.2 改造前后的对比结果我把“无储热改造”传统机组最小出力约束为300MW和“有储热改造”两种场景分别求解得到了以下几组关键对比数据。对比项无储热改造有储热改造机组低谷最小出力300MW320MW实际电出力夜间弃风电量68MWh12MWh日碳排放量3648t3396t系统总运行成本含碳61.2万元58.7万元第一眼看去改造后机组在低谷时段的“电出力”反而从300MW升到了320MW这不科学其实这正是储热改造的作用机理——夜间机组维持320MW的电出力比原来还多了20MW但同时开启了80MW的充热功率等效来看机组从电网侧“净吸收”的电力只有240MW相当于腾出了60MW的新能源消纳空间。这就是“增加电出力、同时增加储热负荷”实现深度调峰的等效原理。5.3 结果背后的机理分析在我设定的算例里储热改造带来的收益主要来自三个部分一是弃风电量从68MWh降到12MWh多消纳的风电替代了部分燃煤发电二是机组在夜间处于更低等效负荷煤耗率下降三是因为碳排放总量下降需要购买的碳配额减少碳成本降低。有意思的是如果我把碳价从50元/吨提高到120元/吨改造的经济优势会更明显总成本的差异从2.5万元扩大到4万元以上。这也侧面印证了一个结论碳价越高储热改造的投资回报越有吸引力。反过来如果碳价过低单纯靠储热改造省下来的碳成本可能无法覆盖改造和运维费用这也是很多电厂在实际投资决策中观望的原因。6. 在这套模型实战中踩过的坑与解决经验最后分享几个我在调试这套Matlab代码时真实遇到的问题都是外面论文和教程里不太会写的细节但实操中极其折磨人。6.1 储热罐SOC状态约束的初始化陷阱第一个坑是SOC的初始化。solve之前的sdpvar变量没有初始值求解器通常会从全零开始迭代。但储能系统的初始SOC在现实中往往不是0通常设定为50%左右这个初值对前几个时段的充放热决策影响巨大。我在第一版代码里直接把SOC(1) SOC_initial作为一个硬约束但忽略了SOC_initial和SOC(24)之间的一致性关系——调度周期末的SOC应该尽量回到初始值附近否则等于在“挪用”下一周期的储能容量。后来我加了一项终端SOC约束或者目标函数里加终端SOC偏差惩罚项结果才合理。6.2 yalmip对多个约束添加方式导致的性能差异yalmip支持在循环里逐条添加约束也支持一次传入约束数组。很多人习惯在for循环里写成for t 1:T Constraints [Constraints, P(t) P_wind(t) Load(t) H_in(t)]; end这种做法在T24时问题不大但如果扩展到8760个小时的全年调度甚至只是T9615分钟一个点循环内反复拼接数组会让模型构建时间暴增。更好的写法是用矩阵表达式一次性构建P P_wind Load H_in; % 向量化约束对所有t同时成立这不仅仅是代码风格问题。yalmip内部对向量约束的处理比标量循环高效得多我之前一个24时段的算例改成向量化后模型构建时间从17秒降到了0.8秒。6.3 求解器选择不当导致的“伪最优解”使用intlinprog求解时默认的整数可行性容差IntegerTolerance和最优性容差RelativeGapTolerance如果保持默认值有时会返回一个事实上次优的解尤其在储能SOC这种连续性很强的变量上表现为SOC曲线出现小幅抖振。我的做法是显式设定求解器的容差参数options optimoptions(intlinprog, RelativeGapTolerance, 0.001, ... IntegerTolerance, 1e-5);如果用yalmip可以在solve命令中的options结构里传入cplex或gurobi的mipgap参数。这个细节能让结果在论文复现和对比分析时更可信。6.4 防止同时充放热约束推导的细节最后再说说防止同时充放热的0-1约束。我在2.3小节给出的做法是0 H_in H_in_max * S_state; 0 H_out H_out_max * (1 - S_state);第一版我在运行时发现求解器偶尔会给出某个时段H_in和H_out同时大于0但数值很小的结果因为如果数值足够小对目标函数的扰动和SOC状态的影响都微乎其微求解器在容差范围内就“无所谓”了。解决方法是在目标函数里加一个非常小的惩罚项例如0.001 × sum(H_in H_out)或者直接给H_in设置一个下限阈值当S_state1时H_in必须大于某个最小充热功率避免求解器用接近于零的充热流量来凑约束。这种问题在实际工程中不会出现但在数学模型里如果不处理就会成为论文评审时被质疑的逻辑漏洞。回头来看火电机组储热改造的低碳经济调度模型在数学上并不算一个新提出的复杂模型它的价值更多在于“储热系统与传统机组运行约束的耦合建模”以及“如何在碳交易机制下让改造过的机组在调度中真正发挥价值”。Matlab和yalmip的组合让模型的搭建和迭代非常顺手——我甚至觉得用yalmip描述储能约束时的直观程度比手推矩阵形式要高出一个维度调试成本大幅降低这也是我为什么愿意把这一套方案完整分享出来的主要原因。如果你的课题也需要在Matlab里做类似的电力系统优化可以试着先从我的这个框架出发替换成你自己的机组参数和负荷曲线应该能很快跑出第一版结果。 火电机组储热改造的低碳经济调度模型与Matlab实现
做电力系统调度的朋友这两年应该没少被“火电灵活性改造”和“双碳目标”这两个词轮番轰炸。火电机组储热改造通俗点说就是给传统火电机组加装一个巨大的“热水瓶”——在新能源大发、电力过剩的时候把多余的电能转换成热能存起来等到晚高峰电力紧张的时候再把储存的热能释放出来发电或供热。这样一来火电机组就不再是那个“要么满发、要么停机”的笨重大块头而是能灵活调节、配合新能源消纳的调节性电源。我最近正好在Matlab里把这套“储热改造”的模型完整跑通结合低碳经济调度目标做了优化今天把整个思路、模型搭建过程、以及我用Matlab实现时的经验和踩过的坑一次性整理出来分享给正在做相关课题的研究生或者刚入行的工程师朋友。这篇文章不仅适合做纯理论研究的读者也适合那些想把模型落地成可运行代码的实操派。看完之后你至少能收获一套明确的目标函数和约束条件的数学表达、一份能直接参考的Matlab代码框架以及几个在结果调试阶段非常容易卡壳的问题排查思路。1. 为什么储热改造能同时兼顾“低碳”和“经济”先说一个很多人容易混淆的概念。提到火电灵活性改造很多人第一反应是“降低最小出力”也就是让机组在低负荷下也能稳定运行。但单纯压负荷是有物理极限的锅炉的稳燃条件、脱硝系统的烟温窗口都限制了机组不能无限往下压。储热改造的核心思路则完全不同它不是在“发电”这个环节上硬扛而是把能量在“电”和“热”两个形态之间做一个缓冲和转移。1.1 火电厂的“热水瓶”是如何参与调度的具体的工程实现上主流的路线有两种一种是抽汽蓄热在机组运行期间从汽轮机中间级抽出一部分蒸汽加热储热介质导热油、熔盐或者高压热水把热能储存在罐体里另一种是电极锅炉直接用多余的电能加热水或熔盐。前者的能量转化效率高抽汽本身是热能的直接转移后者在“深度调峰”场景下更常用——因为电极锅炉完全可以把机组在低谷时段变成一个大功率用电设备硬生生在厂内多造出一块负荷。两种方式殊途同归都会在系统层面留下一个关键特征火电机组的功率可以被“分时转移”。低谷时段多发电存起来或者用电制热存起来高峰时段释放。对于调度模型而言我们关心的不是罐体内部怎么对流换热而是它在一天24小时里每个时段充了多少热、放了多少热、罐内剩余多少热。这就是抽象的“储热罐模型”。1.2 低碳经济调度到底在优化什么传统的经济调度目标是让系统总发电成本最小本质上是燃料成本的最小化。有了碳交易机制之后成本函数里多了一项碳配额成本——机组排放多少二氧化碳就得购买相应的碳排放配额。这样原来“烧煤便宜”的算计就变了如果一个机组碳排放强度高它的综合运行成本会显著上升低碳机组和可再生能源在调度排序中的优先级就会自然提前。所以这里的“低碳”和“经济”不是两个割裂的目标而是通过碳价这个信号把碳排放量“货币化”之后统一纳入经济调度这一个目标函数里。这就是低碳经济调度的核心思想以含碳成本的总运行成本为目标以系统安全稳定运行为约束求解各机组的出力计划。对于储热改造机组我们还要额外考虑储热罐的充放热计划让它在碳约束下发挥最大的调峰和减排价值。2. 数学模型的完整搭建目标函数与三层约束体系这一部分是整篇博文的灵魂。模型搭得清不清晰直接决定后面Matlab代码是“一次跑通”还是“改了三天还在报错”。我在建模时把它拆成三块目标函数的构建、常规电力系统约束、储热改造特有的运行约束。2.1 目标函数燃料成本、碳成本、弃风惩罚的加权博弈低碳经济调度的目标函数可以写成min F Σ(F_fuel F_carbon) F_wind_penalty逐项拆开来说。第一项是燃料成本也就是煤耗成本。工程上通常用二次函数近似描述机组出力P与煤耗量之间的关系F_fuel a × P² b × P c其中a、b、c是煤耗特性系数。这个二次项在后面的Matlab求解时会被线性化处理一会展开讲。第二项是碳成本。引入碳交易机制后每台机组根据实际排放量M和免费配额E_free的差值产生购买或出售配额的成本。排放量M和出力之间存在线性关系M δ × P这里的δ就是机组的碳排放强度系数t/MWh。碳成本的计算式是F_carbon λ × (M - E_free)如果实际排放低于配额括号内是负值意味可以卖配额赚钱这就是低碳机组的经济激励。第三项是弃风惩罚。我习惯在目标函数里加一个单位弃风惩罚成本如果不加模型可能出现“为了省钱而大量弃风”的极端结果——这在碳价较低、煤价也比较低的时候经常发生。加了惩罚之后模型会在“减少燃煤”和“多消纳新能源”之间自动找平衡。2.2 功率平衡与机组运行约束最容易被忽略的细节这一类约束属于电力系统调度模型的标配但我还是要强调几个容易出错的点。首先是功率平衡约束。电力系统在任何时刻都必须保持发电和负荷的实时平衡即所有机组出力之和等于负荷需求。这个约束我在建模初期犯过一个低级错误——忘了把储热系统的耗电功率加进负荷侧。加了电极锅炉之后厂内负荷增加功率平衡的等号右边不是简单的系统负荷而是“系统负荷储热系统用电功率-储热放热供出的电功率”如果储热系统是热电联产形式则还需要耦合热负荷。其次是机组出力上下限约束。对于改造后的机组出力的上限不变但下限值得仔细讨论。纯粹通过降低最小出力改造的机组下限可以从50%额定降到30%甚至20%但如果是靠储热实现调峰的机组下限其实并没有降低多少——它是在低负荷时段把多出来的电转成热蓄起来等效上相当于“用了更少的电、但还是烧着同样多的煤”。这一点在模型里要小心我见过有人直接在代码里把储热机组的最小出力改成了20%结果目标函数算出一个非常“乐观”却不现实的结果。还有爬坡约束。机组在相邻时段之间的出力变化不能超过爬坡速率限制。储热系统的充放热切换会导致机组电出力变化率变大如果爬坡约束设置不当求解器会很难收敛。2.3 储热罐模型连续变量 状态耦合约束储热罐是本模型的特殊之处也是Matlab代码中最容易出bug的地方。它的状态可以用四个变量来描述SOC(t)t时段末罐内剩余热量MWh对应储能系统的荷电状态H_in(t)t时段的充热功率MWH_out(t)t时段的放热功率MWS(t)0-1变量表示储热系统整体处于充热还是放热状态用于防止同时充放核心的状态转移方程是SOC(t) SOC(t-1) (H_in × η_in - H_out / η_out) × Δt其中η_in是充热效率、η_out是放热效率。注意一个关键细节放热侧的效率放在分母上。也就是说如果想要在放热侧输出1MWh的热量储热罐内部其实要释放出1/η_out的热量这是工程上的定义差异直接写代码时容易搞反。同时还要约束充放热功率上限和罐容上限0 ≤ H_in(t) ≤ H_in_max × S(t) 0 ≤ H_out(t) ≤ H_out_max × (1 - S(t)) SOC_min ≤ SOC(t) ≤ SOC_maxS(t)是0-1整数变量它的作用是禁止储热罐在同一个时段既充又放——虽然现实中这种操作没有物理意义但如果不加这个变量约束求解器就可能在数学上做出“既充又放”的荒谬结果来凑最优值。2.4 碳排放约束配额与阶梯碳价的建模碳交易的建模方式直接影响目标函数的线性化难度。最基础的做法是固定碳价碳成本是线性的但实际碳市场中配额价格往往是阶梯式的——超过一定排放量之后超出部分价格更高即阶梯碳价。在Matlab实现里我用了一个分段线性化的方法来处理这个阶梯函数代价是增加了几组辅助变量和额外的约束但换来了更贴合实际的价格信号。阶梯碳价函数可以表达为基于同一个Class同一类型or同一默契设定下它是分段函数排放量M在配额区间内价格为λ1超过配额后价格为λ2λ2 λ1。我把它写成等价的两段线性函数通过一组0-1变量约束每一段的选择这在Matlabyalmip下的实现比较直观我会在代码部分展示。3. 求解流程与混合整数线性规划(MILP)问题的转化的完整思路模型搭建完成之后不要急着把所有约束typeset进Matlab就开始求解。我需要先问自己一个问题这个模型是线性规划(LP)、非线性规划(NLP)还是混合整数线性规划(MILP)选错了求解器后面所有工作都是白费。3.1 非线性项的处理方法与理论依据在我搭建的目标函数中燃料成本a×P²是一个典型的二次项属于非线性规划。但是如前面所说实际求解时我将其分段线性化处理——把机组出力的可行区间分成若干段每一段用一个线性函数近似引入对应的0-1变量。这样一来目标函数和约束全部变成线性形式可以用成熟的混合整数线性规划求解器高效求解。为什么不直接保留二次项用非线性求解器呢因为电力系统调度模型动辄几十台机组乘上24个时段非线性求解器的求解时间可能呈指数增长收敛性和全局最优性都没有保证。MILP的求解器经过几十年的发展分支定界算法已经非常成熟在几秒到几分钟内就可以获得全局最优解。实践下来分段线性化造成的精度损失非常小完全在工程接受范围内。3.2 为什么阶梯碳价让问题变成了MILP阶梯碳价中引入了0-1变量选择“处于哪个碳价段”这就是混合整数变量的来源。储能系统的充电/放电状态S(t)是另一个0-1变量。这些整数变量的引入把原本的线性规划问题升维成了MILP问题。MILP的求解时间随整数变量数量上升而增长所以在编写Matlab代码时一个非常关键的优化习惯是尽量减少不必要的整数变量。能通过连续变量边界约束解决的就不要引入0-1变量。3.3 求解器的选型与Matlab环境配置建议常见的Matlab求解器搭配有两种。一种是Matlab自带的linprog和intlinprog前者解决线性规划后者解决混合整数线性规划另一种是Matlabyalmip工具箱再外接cplex或者gurobi。我个人强烈建议使用yalmip作为建模语言哪怕你最终用intlinprog求解yalmip的符号建模方式也能大幅减少调试时间——它的约束表达几乎和数学公式一一对应代码可读性比纯粹用矩阵拼装高出两个量级。求解器方面cplex和gurobi的求解性能在MILP方向上比intlinprog好很多。对于小规模24时段、几台机组的问题intlinprog勉强够用但如果要做24小时×数十台机组或者进行多场景分析还是建议用cplex或者gurobi。4. Matlab代码实现的架构设计与关键模块解读终于到了代码部分。我按照“数据准备–模型搭建–求解–结果可视化”四个步骤来组织代码下面的内容都是基于这个框架展开的。4.1 数据准备环节的细节参数与数据结构为了方便修改和测试不同场景我用一个脚本文件load_data.m来统一管理所有基础数据。这里面包含系统负荷曲线24个时段的有功负荷需求MW风电预测出力曲线24个时段的可用风电出力MW火电机组参数额定容量、技术最小出力、煤耗系数、爬坡速率、碳排放强度储热罐参数罐容上限、充热效率、放热效率、最大充/放热功率、初始SOC碳交易参数碳价、免费配额比例或者绝对配额量数据准备阶段最容易被忽略的是各参数的物理单位一致性。我就吃过亏负荷和机组出力单位用MW但储热罐容量用了MWh在状态转移方程里计算SOC时需要乘以Δt单位小时才能把功率换算成能量。有时候一不留神单位混了模型结果看起来正常但数值解读完全错误。4.2 基于yalmip的约束构建核心代码下面是模型搭建的核心框架我用yalmip语法直接声明变量并添加约束代码逻辑跟2.3小节的数学公式一一对应。% 时间段定义 T 24; dt 1; % 每个时段时长小时 % 变量声明 P sdpvar(1, T); % 火电机组出力 P_wind sdpvar(1, T); % 实际风电消纳功率 P_wind_curtail sdpvar(1, T); % 弃风功率 H_in sdpvar(1, T); % 储热充热功率 H_out sdpvar(1, T); % 储热放热功率 SOC sdpvar(1, T); % 储热罐状态 S_state binvar(1, T); % 充/放状态 % 功率平衡约束注意储热系统耗电体现在负荷侧 Constraints []; for t 1:T Constraints [Constraints, P(t) P_wind(t) Load(t) H_in(t)]; Constraints [Constraints, P_wind(t) P_wind_curtail(t) P_wind_forecast(t)]; end % 机组出力上下限约束 Constraints [Constraints, P_min P P_max]; % 储热罐状态转移约束 SOC(1) SOC_initial (H_in(1) * eta_in - H_out(1) / eta_out) * dt; for t 2:T Constraints [Constraints, ... SOC(t) SOC(t-1) (H_in(t) * eta_in - H_out(t) / eta_out) * dt]; end % 储热罐容量约束 Constraints [Constraints, SOC_min SOC SOC_max]; % 充放热功率约束含互斥状态 Constraints [Constraints, 0 H_in H_in_max * S_state]; Constraints [Constraints, 0 H_out H_out_max * (1 - S_state)]; % 目标函数燃料成本 碳成本 弃风惩罚 Objective sum(fuel_coeff_a * P.^2 fuel_coeff_b * P fuel_coeff_c) ... carbon_price * sum(carbon_intensity * P - carbon_quota) ... wind_penalty * sum(P_wind_curtail);第一版跑通之后我需要额外说明一点关于燃料成本的二次项处理。在上面的Objective中直接写了P.²这个在yalmip中是二次规划(QP)的问题描述方式在某些求解器下可以被接受。但为了让模型完全线性化、能够用纯MILP求解器处理我实际代码里用的是对P分段线性化的做法——把P的可行区间分成K段用sdpvar的二进制辅助变量组实现每段线性近似。这个改动虽然增加了代码量但求解稳定性明显更好。4.3 结果可视化与关键图表输出求解完成后我用Matlab自带的plot命令绘制几类关键图表机组出力曲线与负荷曲线的对比直观看到储热低谷时段“增负荷”的效果储热罐SOC随时间变化曲线检查能量是否守恒、是否超出罐容上限弃风率柱状图对比改造前后的消纳效果成本构成堆叠图燃料/碳/弃风惩罚展示低碳改造带来的成本结构变化可视化不只是为了写论文插图更是排查模型错误的重要工具。我记得第一次跑出来的SOC曲线竟然在凌晨2点和凌晨6点之间出现了突然掉到负数的情况排查了半天才意识到是初始SOC设置和约束顺序导致的状态转移方程逻辑错误。类似这种问题只盯着一堆数字很难感觉到异常画成图一眼就看出来了。5. 实际算例测试从仿真结果看储热改造的价值光有模型和代码还不够我用一个简化的算例来验证模型的有效性。算例规模不大但是足够说明储热改造在低碳经济调度中的核心价值。5.1 算例场景与参数设定我设定了一个含一台600MW火电机组、200MW风电场的简化系统。机组额定出力600MW技术最小出力不含储热时为300MW煤耗系数取600MW机组的典型值碳排放强度0.8t/MWh碳价取50元/吨。储热系统容量设定为300MWh最大充热功率100MW充热效率0.9放热效率0.85。负荷曲线采用典型冬季日负荷曲线白天有两个高峰午高峰和晚高峰夜间出现低谷约350MW。风电出力设定为夜间大发的典型反调峰特性——夜间风电出力高、负荷低这恰恰是弃风最严重的时段。5.2 改造前后的对比结果我把“无储热改造”传统机组最小出力约束为300MW和“有储热改造”两种场景分别求解得到了以下几组关键对比数据。对比项无储热改造有储热改造机组低谷最小出力300MW320MW实际电出力夜间弃风电量68MWh12MWh日碳排放量3648t3396t系统总运行成本含碳61.2万元58.7万元第一眼看去改造后机组在低谷时段的“电出力”反而从300MW升到了320MW这不科学其实这正是储热改造的作用机理——夜间机组维持320MW的电出力比原来还多了20MW但同时开启了80MW的充热功率等效来看机组从电网侧“净吸收”的电力只有240MW相当于腾出了60MW的新能源消纳空间。这就是“增加电出力、同时增加储热负荷”实现深度调峰的等效原理。5.3 结果背后的机理分析在我设定的算例里储热改造带来的收益主要来自三个部分一是弃风电量从68MWh降到12MWh多消纳的风电替代了部分燃煤发电二是机组在夜间处于更低等效负荷煤耗率下降三是因为碳排放总量下降需要购买的碳配额减少碳成本降低。有意思的是如果我把碳价从50元/吨提高到120元/吨改造的经济优势会更明显总成本的差异从2.5万元扩大到4万元以上。这也侧面印证了一个结论碳价越高储热改造的投资回报越有吸引力。反过来如果碳价过低单纯靠储热改造省下来的碳成本可能无法覆盖改造和运维费用这也是很多电厂在实际投资决策中观望的原因。6. 在这套模型实战中踩过的坑与解决经验最后分享几个我在调试这套Matlab代码时真实遇到的问题都是外面论文和教程里不太会写的细节但实操中极其折磨人。6.1 储热罐SOC状态约束的初始化陷阱第一个坑是SOC的初始化。solve之前的sdpvar变量没有初始值求解器通常会从全零开始迭代。但储能系统的初始SOC在现实中往往不是0通常设定为50%左右这个初值对前几个时段的充放热决策影响巨大。我在第一版代码里直接把SOC(1) SOC_initial作为一个硬约束但忽略了SOC_initial和SOC(24)之间的一致性关系——调度周期末的SOC应该尽量回到初始值附近否则等于在“挪用”下一周期的储能容量。后来我加了一项终端SOC约束或者目标函数里加终端SOC偏差惩罚项结果才合理。6.2 yalmip对多个约束添加方式导致的性能差异yalmip支持在循环里逐条添加约束也支持一次传入约束数组。很多人习惯在for循环里写成for t 1:T Constraints [Constraints, P(t) P_wind(t) Load(t) H_in(t)]; end这种做法在T24时问题不大但如果扩展到8760个小时的全年调度甚至只是T9615分钟一个点循环内反复拼接数组会让模型构建时间暴增。更好的写法是用矩阵表达式一次性构建P P_wind Load H_in; % 向量化约束对所有t同时成立这不仅仅是代码风格问题。yalmip内部对向量约束的处理比标量循环高效得多我之前一个24时段的算例改成向量化后模型构建时间从17秒降到了0.8秒。6.3 求解器选择不当导致的“伪最优解”使用intlinprog求解时默认的整数可行性容差IntegerTolerance和最优性容差RelativeGapTolerance如果保持默认值有时会返回一个事实上次优的解尤其在储能SOC这种连续性很强的变量上表现为SOC曲线出现小幅抖振。我的做法是显式设定求解器的容差参数options optimoptions(intlinprog, RelativeGapTolerance, 0.001, ... IntegerTolerance, 1e-5);如果用yalmip可以在solve命令中的options结构里传入cplex或gurobi的mipgap参数。这个细节能让结果在论文复现和对比分析时更可信。6.4 防止同时充放热约束推导的细节最后再说说防止同时充放热的0-1约束。我在2.3小节给出的做法是0 H_in H_in_max * S_state; 0 H_out H_out_max * (1 - S_state);第一版我在运行时发现求解器偶尔会给出某个时段H_in和H_out同时大于0但数值很小的结果因为如果数值足够小对目标函数的扰动和SOC状态的影响都微乎其微求解器在容差范围内就“无所谓”了。解决方法是在目标函数里加一个非常小的惩罚项例如0.001 × sum(H_in H_out)或者直接给H_in设置一个下限阈值当S_state1时H_in必须大于某个最小充热功率避免求解器用接近于零的充热流量来凑约束。这种问题在实际工程中不会出现但在数学模型里如果不处理就会成为论文评审时被质疑的逻辑漏洞。回头来看火电机组储热改造的低碳经济调度模型在数学上并不算一个新提出的复杂模型它的价值更多在于“储热系统与传统机组运行约束的耦合建模”以及“如何在碳交易机制下让改造过的机组在调度中真正发挥价值”。Matlab和yalmip的组合让模型的搭建和迭代非常顺手——我甚至觉得用yalmip描述储能约束时的直观程度比手推矩阵形式要高出一个维度调试成本大幅降低这也是我为什么愿意把这一套方案完整分享出来的主要原因。如果你的课题也需要在Matlab里做类似的电力系统优化可以试着先从我的这个框架出发替换成你自己的机组参数和负荷曲线应该能很快跑出第一版结果。