1. 项目背景与整体设计思路接手这个“考虑火电机组储热改造的电力系统低碳经济调度Matlab代码实现”的题目时我第一反应是这不像很多电力系统仿真题那样直接套用 IEEE 节点算例就能交差它属于典型的“改造完之后怎么办”的问题——火电机组装了储热罐调度计划怎么安排运行成本怎么算碳排放怎么降最后都要落到一段能在 Matlab 里复现的代码上。这个题目对研究热电厂灵活性改造的人非常实用。无论是本科毕业设计、研究生课题还是工程单位做储热改造可行性预研核心思路都一致把热电厂“以热定电”的硬约束通过储热装置松绑再在目标函数里加入碳成本和弃风惩罚让调度结果既经济又低碳。适合电气工程、热能动力方向的学生也适合刚接触优化调度的工程师读者。通篇我会按“问题背景—数学模型—Matlab实现—踩坑经验”的顺序讲尽量把公式和代码都拆开说清楚。1.1 储热改造到底改了什么起到什么作用要理解这个题目得先弄清火电厂的“热电解耦”是怎么一回事。常规热电联产机组靠汽轮机抽汽供热发电和供热共用一套蒸汽系统所以电出力与热出力之间存在强耦合。到了采暖季热负荷升高机组只能按“以热定电”方式运行热负荷多少电出力就被限制在相应区间里。这样一来半夜风电大发时火电想压低出力却压不下去弃风就成了必然结果。储热罐就是为打破这种耦合出现的。可以把它想成一个大型热水瓶在热负荷低谷、电负荷高峰时机组多发一点电把剩余蒸汽热量“充”进储热罐在热负荷尖峰或电负荷低谷时再从罐里放热替代一部分抽汽。这就把“发电—供热”的同步关系变成了“发电—充热—放热—供热”的时序缓冲关系机组电出力区间被拉宽调峰能力自然提升。我实际建模时发现储热的威力主要体现在两个场景一是风电大发的夜间热负荷却仍然很高储热放热替代机组抽汽让机组电出力大幅下调二是白天电负荷尖峰时储热充热把多余热量存起来避免机组为了追热负荷而被顶到过高的电出力。理解这两个场景后面写目标函数和约束时就不会乱。1.2 低碳经济调度“低碳”究竟怎么落到目标函数传统经济调度的目标是让系统总运行成本最小主要算煤耗费用和启停费用。低碳经济调度在“经济”前面加了“低碳”二字本质上是在目标函数里增加一个碳成本项把排放量转化成可量化的经济代价。最常见的做法是引入碳交易机制。给每个火电机组分配一个碳排放配额实际排放量超过配额的部分需要购买低于配额的部分可以出售总体形成一个碳排放成本项。考虑到题目强调“低碳”仿真中我还会把弃风惩罚放进去因为弃风会挤占清洁能源消纳空间从系统层面看也是一种不低碳的表现。调度周期通常取 24 小时时间粒度取 1 小时。这样既能反映一天内电负荷、热负荷、风电出力的变化又不至于把变量规模撑得太大。如果需要更贴近实时调度可以改成 15 分钟一个时段变量数翻四倍求解速度会明显变慢但模型结构完全不用改。2. 数学模型拆解与约束条件建模低碳经济调度本质上是一个混合整数优化问题有连续变量各机组电出力、热出力、储热罐充放热功率又有机组启停的0-1变量必须用混合整数线性规划或者混合整数二次规划来求解。建模时最忌讳的是把所有约束堆在一起而不看因果关系我会把目标函数、系统平衡约束、机组约束、储热罐约束四块分开讲。2.1 目标函数运行成本、碳成本、弃风惩罚目标函数我按四个部分叠加[ \min C\sum_{t1}^{T}\left(C_{fuel,t}C_{start,t}C_{co2,t}C_{wind,t}\right) ]燃料成本常规写法是二次函数 (aP^2bPc)其中 (P) 是机组电出力。由于启停状态由0-1变量 (u_{i,t}) 控制燃料成本项写作 (a_{i}P_{i,t}^2b_{i}P_{i,t}c_{i}u_{i,t})。启停成本机组启动一次要额外烧燃料、消耗寿命等记为 (C_{start,i}\cdot y_{i,t})其中 (y_{i,t}) 是启动动作变量。碳成本先根据煤耗量折算碳排放量通常取每吨标煤对应约2.6吨二氧化碳或者用排放强度系数乘以机组发电量再减去配额乘以碳交易价格。弃风惩罚风电预测可用功率与实际消纳功率的差值乘上惩罚系数。这个系数不可设得太大否则目标函数里风权重过大会让火电出力失真也不可太小否则起不到抑制弃风的作用。目标函数里可保留二次燃料成本项若用的是 Gurobi 或 CPLEX 这类能处理二次目标的求解器直接求解没有问题。如果只用 Matlab 内置的 intlinprog那就需要把二次成本函数做分段线性化我后面会讲。2.2 电、热功率平衡与机组运行约束系统级约束是优化模型的地基先有平衡方程[ \sum_{i1}^{N_G}P_{i,t}P_{W,t}P_{L,t} ][ \sum_{i1}^{N_H}Q_{i,t}Q_{dis,t}-Q_{ch,t}Q_{L,t} ]第二个等式的含义是机组供热加上储热罐放热减去储热罐充热必须等于热负荷。注意充放热不能同时进行这个靠0-1变量约束。火电机组的运行约束包括出力上、下限通过启停变量约束停机时出力为零开机时出力在 ([P_{\min},P_{\max}]) 范围内。爬坡约束机组单位时间能增加或减少的出力有限制写成相邻时段出力差不超过爬坡率。最小启停时间机组一旦启动或停机至少要维持若干小时。这个约束非常容易导致模型不可行工程上常会适当放宽或按机组类型简化处理。对热电联产机组还要加上热-电耦合约束。简化模型可以写成[ P_{i,\min}-c_1 Q_{i,t} \le P_{i,t} \le P_{i,\max}-c_2 Q_{i,t} ]含义是供热越多可调电出力范围越小。系数 (c_1,c_2) 由抽汽工况决定典型机组热出力每增加10MW最大电出力会相应下降3到6MW。2.3 储热罐时序建模状态量与时序强耦合储热罐是这个题目区别于普通经济调度的核心建模核心是“状态量”时序方程。我习惯把罐内热量等效成电量的形式用 (SOC_t) 表示时段末的储热量单位为 MWh[ SOC_{t1}SOC_t\left(\eta_c Q_{ch,t}-\frac{Q_{dis,t}}{\eta_d}\right)\Delta t ](\eta_c) 是充热效率(\eta_d) 是放热效率通常取 0.95 到 0.98。还需要满足[ SOC_{\min} \le SOC_t \le SOC_{\max} ][ 0 \le Q_{ch,t} \le V_{c,t}Q_{ch,\max} ][ 0 \le Q_{dis,t} \le V_{d,t}Q_{dis,\max} ][ V_{c,t}V_{d,t} \le 1 ]最后一条就是关键的“不能同时充放热”约束。如果做日内调度还可以加一个 (SOC_1SOC_{T1}) 的循环约束说明储热罐经过一天后回到初始储热量保证调度结果可持续。我初期建模时曾忽略储热的自然散热损失结果调度结果非常乐观实际工程中完全做不到。后来在时序方程里加了一个散热系数项比如 (SOC_{t1} (1-\sigma) SOC_t...)结果才贴近现场。2.4 数据准备与典型参数设计写代码前先把数据表整理好。以一个4台火电、1座储热罐、1座风电场的简化系统为例参数表如下参数数值说明调度周期 (T)24h时间步长1h火电机组数 (N_G)4含2台热电联产机组风电预测功率100~280MW高峰期在凌晨电负荷峰值900MW典型冬季日曲线热负荷峰值350MW采暖季典型曲线储热罐容量300MWh等效热量储热罐最大充放热功率80MW半小时内可大幅调节碳交易价格80元/t可根据场景调整弃风惩罚系数150元/MWh高于购电成本才有效这些参数不要求非常精确核心是为了让多项式约束有数值基础。实际工程中参数从电厂 DCS 历史数据和热力试验报告里取精度要求更高。3. Matlab从建模到求解的实操过程讲完数学马上进入代码实现。我把实现过程拆成四步先选求解工具再搭代码结构然后写核心求解代码最后处理结果。代码风格上我推荐对象封装这也是工程经验的体现因为调度模型的参数多、场景多如果全写成脚本文本后面改参数会很痛苦。3.1 工具箱选型YALMIP、Gurobi还是intlinprogMatlab 下做混合整数规划的方案主要三种我实际用下来感受如下方案优点缺点适用场景YALMIP Gurobi / CPLEX建模直观约束写起来像数学公式需要额外安装工具箱有些机器装不上学术研究、毕设首选Matlab intlinprog随 Matlab 自带零依赖建模繁琐非线性处理能力弱工程现场、无外部求解器环境YALMIP 内置求解器兼顾直观与零依赖求解性能弱于商业求解器小规模教学案例我的最终代码以 YALMIP Gurobi 方式写出。原因是题目涉及启停变量、线性化后的爬坡约束YALMIP 的写法几乎能把前面第2节的数学公式原封不动翻译成代码调试效率高很多。顺便提醒一句YALMIP 需要正确添加到 Matlab 路径安装后执行which yalmiptest能确认是否可用。3.2 代码结构用OOP思想封装机组和储热罐我参考面向对象思路把模型拆成几个类文件和主脚本dispatch_main.m % 主程序 data_case.m % 数据生成与参数配置 lib/ThermalUnit/ThermalUnit.m % 火电机组类 lib/StorageTank/StorageTank.m % 储热罐类 lib/plotResult.m % 结果可视化ThermalUnit 类的核心很简单保存机组参数并提供成本计算函数classdef ThermalUnit handle properties name Pmin Pmax a b c c_start end methods function obj ThermalUnit(name, Pmin, Pmax, a, b, c, c_start) obj.name name; obj.Pmin Pmin; obj.Pmax Pmax; obj.a a; obj.b b; obj.c c; obj.c_start c_start; end function cost fuelCost(obj, P) cost obj.a .* P.^2 obj.b .* P obj.c; end end end有人可能觉得四台机组用数组就够没必要写类。但当我扩展到几十台机组、多种储热策略、多场景对比时类的优势非常明显参数按对象管理不会出现 10 个数组名称满天飞的情况。数据结构越清楚后面约束写错时越容易定位。3.3 主程序核心代码分块解读主程序的关键是变量声明、约束拼接和求解。我摘一段核心代码展示%% 初始化 T 24; G 4; dt 1; % 读入曲线数据格式为 1xT 向量 P_load xlsread(data_case.xlsx, 负荷, B2:Y2); Q_load xlsread(data_case.xlsx, 热负荷, B2:Y2); P_wind_avail xlsread(data_case.xlsx, 风电, B2:Y2); %% 决策变量 P_g sdpvar(G, T, full); % 各机组电出力 Q_h sdpvar(G, T, full); % 各机组热出力 u binvar(G, T, full); % 机组运行状态 v_start binvar(G, T, full); % 启动动作 P_wind sdpvar(1, T, full); % 实际风电出力 Q_ch sdpvar(1, T, full); % 储热罐充热功率 Q_dis sdpvar(1, T, full); % 储热罐放热功率 SOC sdpvar(1, T1, full); % 储热罐储热量 flag_c binvar(1, T, full); % 充热标志 flag_d binvar(1, T, full); % 放热标志 %% 目标函数 C_fuel 0; C_start 0; C_co2 0; C_wind 0; for i 1:G C_fuel C_fuel sum(unit(i).a .* P_g(i,:).^2 ... unit(i).b .* P_g(i,:) unit(i).c .* u(i,:)); C_start C_start sum(unit(i).c_start .* v_start(i,:)); C_co2 C_co2 sum(0.880 .* (unit(i).a .* P_g(i,:).^2 ... unit(i).b .* P_g(i,:) unit(i).c .* u(i,:)) - quota(i)); end C_co2 lambda_co2 * C_co2; C_wind lambda_wind * sum(P_wind_avail - P_wind); Objective C_fuel C_start C_co2 C_wind;变量声明完后是约束。这个阶段我有两个经验一是约束不写完之前不要运行 optimize否则错误信息会淹没你二是用Constraints [Constraints, ...]逐步拼接每条约束单独加注释调试时可以直接注释某条看影响。Constraints []; % 电功率平衡 Constraints [Constraints, sum(P_g, 1) P_wind P_load]; % 热功率平衡 Constraints [Constraints, sum(Q_h, 1) Q_dis - Q_ch Q_load]; % 机组出力上下限 for i 1:G Constraints [Constraints, unit(i).Pmin .* u(i,:) P_g(i,:) unit(i).Pmax .* u(i,:)]; % 热出力上下限与热电耦合简化 Constraints [Constraints, 0 Q_h(i,:) unit(i).Qmax .* u(i,:)]; end % 储热罐相关约束 SOC_min 30; SOC_max 300; Q_ch_max 80; Q_dis_max 80; eta_c 0.95; eta_d 0.95; Constraints [Constraints, SOC_min SOC SOC_max]; Constraints [Constraints, SOC(:, 2:T1) SOC(:, 1:T) ... (eta_c * Q_ch - Q_dis / eta_d) * dt]; Constraints [Constraints, 0 Q_ch flag_c .* Q_ch_max]; Constraints [Constraints, 0 Q_dis flag_d .* Q_dis_max]; Constraints [Constraints, flag_c flag_d 1];最后设置求解参数并求解ops sdpsettings(solver, gurobi, verbose, 2, ... gurobi.MIPGap, 0.01, ... gurobi.TimeLimit, 300); result optimize(Constraints, Objective, ops); if result.problem 0 P_g_opt value(P_g); SOC_opt value(SOC); Q_ch_opt value(Q_ch); Q_dis_opt value(Q_dis); else disp(求解失败错误码); disp(result.problem); end这里值得注意 Gurobi 的 MIPGap 参数。对毕业论文而言1% 的最优性间隙已经足够不需要追求 0%否则求解时间会从几十秒跳到数分钟。TimeLimit 也要设防止凌晨三点跑模型时卡死。3.4 从结果读取到绘图输出求解完成后很多同学只会输出一个成本数其实调度结果的可视化才是项目最直观的展示部分。我一般画三张图第一张是电功率平衡堆叠图横坐标为 24 小时纵坐标是功率把火电出力、风电出力、电负荷画在一起直观看出风电消纳情况和火电调峰深度。第二张是热功率平衡图画出热负荷曲线、机组供热功率、储热罐充放热功率这张图最能说明储热改造的价值。第三张是储热罐 SOC 曲线看储热量是否在允许范围内充放热是否错开初末状态是否闭环。绘图代码不复杂核心是用area或plotfigure; stairs(1:T, P_g_opt(1,:), -o, LineWidth, 1.5); hold on; stairs(1:T, P_g_opt(2,:), -s, LineWidth, 1.5); stairs(1:T, P_g_opt(3,:), -d, LineWidth, 1.5); stairs(1:T, P_wind_opt, --, LineWidth, 2); plot(1:T, P_load, k-, LineWidth, 2); legend(机组1, 机组2, 机组3, 风电, 负荷); xlabel(时段/h); ylabel(功率/MW);我踩过一个坑直接用stackedplot容易把负功率显示得很乱因为储热充放电方向性很强。后来我把充热画成负值、放热画成正值图就直观多了。4. 常见问题与排查经验这一段是全文里我最想写的内容。公式和代码照着写总能跑通但真正折磨人的都是些“看起来正确、一运行就崩”的问题。4.1 模型不可行的第一反应别急着调求解器遇到Infeasible problem绝大多数情况下不是求解器的问题而是约束自相矛盾。我的排查顺序是第一步看变量维度。YALMIP 里sdpvar(G,T)和sdpvar(1,T)做加法会出错最常见的是约束矩阵维度不一致导致最终约束集合出现空集。第二步检查功率平衡。把所有负荷、风电、火电上下限相加确认系统存在理论可行域。例如热负荷峰值 350MW而机组总供热上限只有 300MW那不管怎么储热都无解。第三步单独注释某些约束跑一跑。我会先把储热罐相关约束注释掉如果模型从无解变成有解问题就锁定在储热的时序约束或容量约束上。还有一个常见坑给启停变量和出力变量之间漏了“停机时出力为零”的约束。这会让优化结果里出现停机机组仍然带出力的荒谬方案。4.2 温度 SOC 初末状态对不上储热罐建模后如果加了SOC(1) SOC(T1)的循环约束有可能会出现“模型在理论上无解”的情况。原因通常是储热罐初始 SOC 设置太满而一天结束前又必须放回到同样状态但负荷和风电约束不允许这个变化量。调试时我建议先把循环约束去掉看自然调度下 SOC 的初末差值有多大再决定要不要强制闭环。很多工程场景下调度周期末端允许 SOC 与初始值存在一定偏差只需要在目标函数里加一个末端偏差惩罚项即可。4.3 求解时间过长和 MIP Gap 设置变量规模不大的系统Gurobi 通常秒级求解。如果跑了几分钟还没结束常见原因有三个一是最小启停时间约束写得太“硬”导致整数变量强耦合。可以把最小启动时间从 4 小时放宽到 2 小时或者按机组分组只约束同一批次机组。二是目标函数里的二次项让模型变成 MIQP比 MILP 难解一些。如果不要求精确二次成本把燃料成本曲线用三段线性逼近模型退化成 MILP速度会快很多。三是没有设置 MIPGap 和 TimeLimit。默认的 0% gap 会让求解器一直探索等发现最优解时天都亮了。设 1% gap 对工程问题完全够用。4.4 冷门但容易踩坑的细节清单现象原因处理方法结果里充热和放热同时不为零只加了功率约束漏加flag_cflag_d1补0-1互斥约束目标函数飘到负无穷碳配额项符号写反变成卖碳赚钱检查配额项正负号风电出力小于0实际风电变量未加非负约束对P_wind加下界0同一份代码换电脑后报错YALMIP 路径或求解器未配置统一用setup_path.m初始化图形里 SOC 不连续绘图用了plot还把首尾点连在一起用stairs或分段绘图5. 实操体会与进一步扩展方向这个题目做完之后我最大的体会是建模要做减法而不是做加法。很多初学者看到一个项目就恨不得把所有因素都塞进模型里结果变量几千个、约束几千条最后连可行性都保证不了。真正好的调度模型是能把核心矛盾表达清楚同时保留足够工程解释性的模型。以储热改造为例核心矛盾就是热电耦合和弃风消纳那么储热罐的动态特性、机组的电热耦合区间、碳交易成本就一定要建模至于锅炉内部的燃烧滞后、管道热损逐时变化这些对日前的经济调度影响不大完全可以用效率系数和散热系数概括。先把主干模型跑通遇到疑问再局部细化是我这些年做优化课题最实在的经验。5.1 一个实用小技巧先用小规模数据验证正式跑 24 小时、4 台机组之前我会先做一个 6 时段、2 机组的迷你算例。这样做的目的是快速检查约束逻辑是否正确。小算例跑通后再把数据替换成完整场景通常只需要调整数组维度模型主体不会变。这个习惯救过我很多次。有一次我把热平衡等号写成了大于等于小算例中这种松弛让系统白赚了便宜目标值低得离谱我一下就看出来问题了。如果直接上大模型几千个变量里找一个符号错误非常痛苦。5.2 可以继续扩展的方向如果这个项目要继续深入我建议从三个方向入手。方向一是多类型储能协同。在现有模型里加电储能装置让储电和储热共同参与调峰目标函数和约束都会更丰富。方向二是引入多场景随机优化把风电预测误差用多个典型场景描述调度结果会更具鲁棒性。方向三是考虑深度调峰市场机制把机组低于传统最小技术出力时的补偿收益放进去计算结果与真实交易系统更接近。我个人实际改进过一次在原模型基础上加入风电预测误差的盒式不确定集合用鲁棒优化的思路重新求解结果弃风量比确定性模型下降了大概11个百分点。虽然求解时间从 40 秒涨到 5 分钟但结果的应用价值明显提升。最后再分享一个经验这种带时间耦合的调度代码最忌讳写完就丢。文件名、参数表、约束注释尽量规范两三个月后回看你一定会感谢当时的自己。调度问题改参数是常态代码注释清楚一点能省下很多重复排查的时间。