1. 项目概述与核心需求解析做电热综合能源系统调度研究的朋友应该都有过这种体会单独做电力系统调度吧风电一多弃风率就压不下来单独做热力系统调度吧热负荷的强约束又让机组运行方式僵得要命。真正把电和热放在一个模型里做日前调度能直观看到可再生能源消纳空间被“掰开揉碎”地优化出来那个成就感确实不一样。这篇分享就是围绕“考虑可再生能源消纳的电热综合能源系统日前经济调度模型”展开的核心工具是Matlab重点是建模思路和代码实现的完整落地。这个课题解决的问题一句话概括就是在一个同时包含供电、供热、储热和可再生能源接入的区域综合能源系统里如何提前24小时日前制定各机组的出力计划让总运行成本最低同时最大程度消纳风电/光伏。它适合三类人看一是正在做电热综合能源系统方向毕业设计或论文研究的学生二是刚接触YALMIP/CPLEX调度建模、想做骨干模型打底的工程师三是想把手里的纯热力调度模型升级为电热联合优化模型的研究人员。从技术链条来说这个模型的核心动作可以拆成几个环节构建典型电热系统拓扑确定目标函数为运行成本煤耗 购电 弃风惩罚等确定约束条件为电功率平衡、热功率平衡、机组出力上下限、爬坡约束、热电联产机组运行域约束、储能/蓄热罐状态约束、电网交互功率约束然后对风电出力不确定性做场景处理或旋转备用处理最后调用求解器得到各时段出力计划和总成本。标题里“可再生能源消纳”这个表述其实对应了模型里的两类硬性指标一是弃风弃光电量占比二是为消纳可再生能源而配置的调度手段如电锅炉、蓄热罐、热网延时发挥的实际作用。这篇博客会从“为什么要把电和热放在一起做日前调度”讲起一直到把Matlab代码的关键片段、参数配置、求解器和常见坑全部展开你可以直接照着搭出一个能跑的模型。1.1 为什么电和热必须联立建模很多刚入门的朋友会问电和热有必要放在一个模型里吗分开调度不是更简单我这里直接给结论分不开。原因在于热电联产机组CHP是天然的电热耦合设备它的电出力和热出力被“以热定电”的运行方式绑定在一起。冬天热负荷上来的时候CHP被迫高负荷运行电出力跟着爬高这时候风电稍微一发力系统就消纳不了只能弃风。如果不把热网、蓄热罐、电锅炉这些“热侧资源”纳入调度模型你就只能干瞪眼看着弃风没有任何调节手段。反过来如果把热网管道的蓄热特性、蓄热罐的充放热能力、电锅炉的灵活调节能力都放进一个优化模型CHP就可以在热负荷高峰期少发电、多抽汽供热同时用电锅炉补足电负荷缺口把弃风降下来。这就是电热联立建模的价值所在——用热侧的灵活性给电侧腾出消纳空间。从数学上看电热联合调度是一个混合整数线性规划MILP问题因为机组启停状态是0/1变量而CHP运行域约束是线性不等式组。MILP模型在Matlab里最通用的解法就是YALMIPRmpex/CPLEX/Gurobi这也是这个项目选用Matlab实现的根本原因——Matlab生态里YALMIP建模非常成熟矩阵化约束写法能把模型代码压缩到几百行而且调试方便可视化也顺手。1.2 模型的整体架构和调度流程这个模型的调度流程可以拆成“数据输入 → 场景生成 → 模型构建 → 求解 → 结果分析”五个阶段。数据输入包括负荷曲线电/热、风电出力预测曲线、机组参数容量、煤耗系数、爬坡速率、CHP运行域顶点、热网参数管道延时、蓄热罐容量、电网购电电价与联络线功率上限。场景生成阶段如果按确定性模型做就用预测曲线直接算如果你想把风电不确定性也纳入考量可以用蒙特卡洛抽样生成若干风电出力场景模型变成随机规划。这里我建议第一次搭模型时先用确定性版本跑通后续再加场景否则调试难度会陡增。模型构建阶段就是把目标函数和约束写成YALMIP格式这里有几个关键设计决策会在后文详细展开CHP运行域的凸包建模方式、热网延时的时间平移近似、蓄热罐的离散时间状态方程。求解阶段推荐CPLEX或GurobiMILP求解速度在毫秒到秒级完全满足日前调度的计算需求。最后是结果分析。需要输出的结果包括各机组各时段电出力/热出力曲线、弃风电量及弃风率、系统总运行成本、蓄热罐充放热状态、联络线交换功率、CHP运行点轨迹。这些结果既能用于论文图表展示也能验证“热侧灵活性释放了多少风电消纳空间”这一核心结论。2. 电热综合能源系统建模思路拆解2.1 系统拓扑哪里来电、哪里来热、哪里储能我搭的典型电热综合能源系统包含以下元件常规火电机组纯电出力热电联产机组可同时出电和热这里用抽汽式CHP风电场出力不确定、优先消纳电锅炉用电转热作为弃风消纳手段蓄热罐热的时移容器以及外部电网通过联络线进行功率交换。热负荷由CHP和电锅炉联合供给电负荷由火电、CHP、风电和购电联合满足。这里需要注意一个建模决策热网管道本身的动态特性水流延时、管道蓄热是否要建模。如果按完整动态模型做需要用偏微分方程描述热网温度场复杂度很高如果按“稳态蓄热系数修正”的简化方式可以在调度模型中用一组线性约束近似管道热惯性。做研究的话推荐先用简化模型把调度问题本身做透再逐步增加热网动态细节。系统内的能量流关系一句话可以描述清楚电能从“火电CHP风电电网”流入母线一部分直接供电负荷一部分被电锅炉变成热能送入热网热能从“CHP电锅炉”流入热网经过蓄热罐的充放调节后满足热负荷。这个拓扑对应到Matlab代码里就是各时段的功率平衡等式约束。2.2 目标函数经济性与消纳指标如何权衡日前经济调度的目标函数是最小化系统总运行成本通常包含四项常规机组煤耗成本通常用二次函数MILP中可线性化、CHP机组煤耗成本与电出力和热出力都相关、向外部电网购电成本、弃风惩罚成本。在需要考虑启停的场景中还需要加入机组启停成本项。写目标函数之前有一个必须先想清楚的问题弃风惩罚成本设多少合适如果设得太小模型宁可弃风也不调整热侧如果设得太大模型会不计代价地让电锅炉满发或让蓄热罐疯狂充放热可能在物理上不合理。参考常见文献做法弃风惩罚系数一般取煤耗成本的数倍比如500-1000元/MWh同时给弃风量设置一个上限比例约束比如弃风率不超过5%防止优化结果过于极端。另一个需要权衡的点是经济性和消纳目标在目标函数里的关系。有些模型会把“弃风率最小”作为第二个优化目标做多目标优化但在日前经济调度这个场景下更常见的做法是把弃风惩罚并入运行成本让模型自行权衡“用热侧灵活性消纳风电节省的燃料成本”与“调节热侧带来的额外运维/电耗成本”。这种方式实现简单而且结果就是单一的最优成本方案不用做帕累托前沿分析更适合工程落地。2.3 约束体系功率平衡、热力平衡和CHP运行域约束体系是建模的重头戏我把它分为四层来拆解。第一层是电功率平衡约束。每个时段内火电出力加CHP电出力加风电消纳量加购电量必须等于电负荷加电锅炉耗电量。这里“风电消纳量”不是系统预测出力而是可变量预测值与消纳量的差值就是弃风量。这样做的好处是模型内部自然计算出弃风指标不用事后处理。第二层是热功率平衡约束。CHP热出力加电锅炉热出力加蓄热罐放热功率减去充热功率等于热负荷。注意蓄热罐的充放热可以不同时进行需要引入充放热状态二值变量来避免同时充放导致的无效循环。第三层是CHP运行域约束。这是整个模型最核心、也是出错率最高的部分。抽汽式CHP机的运行域是一个凸多边形顶点集合可以用线性不等式组描述电出力和热出力必须落在该多边形内。在YALMIP里可以很方便地把“点是否在凸包内”这个几何约束写成线性不等式组但前提是你要正确列出顶点坐标。第四层是爬坡约束和状态约束。火电和CHP都存在单位时间出力变化上限爬坡率蓄热罐有容量上下限和充放热功率上限联络线有传输功率上限这些都要一一起约束。还需要注意机组出力上下限约束不能仅写容量区间——对CHP来说电出力上限和热出力上限是耦合的必须用运行域多边形来约束而不是分别写两个独立区间。2.4 热网延时与蓄热罐的建模细节热网输配过程中存在水流延时热量从热源输出到用户端需要一定时间。在日前调度模型中最常见的处理方式是近似延时模型把热网分成若干管道段每段的传输延时是固定的那么时段的调度计划里热源在当前时段输出的热量会在tL时段才到达负荷侧。这种“纯延时”近似在工程上可接受且能在MILP框架内用简单的时段平移实现。蓄热罐的建模相对直接但有几个细节我会特别提醒。蓄热罐的状态方程是SOC(t1) SOC(t) η_ch * P_ch(t) - P_dis(t) / η_dis其中SOC是蓄热量η是充放热效率。由于充放热过程存在效率损失模型天然会避免同时充放如果加了充放状态二值变量还可以严格杜绝这种无效操作。蓄热罐的初始SOC一般设为满罐的一半避免优化结果依赖初始条件。关于热网延时的参数取值我做过的算例里区域供热网典型延时在30-60分钟量级对应到1小时分辨率调度模型里就是L1的时段平移。如果调度时间分辨率设为15分钟延时参数就变成L2到L4。这里提醒一点延时参数不是拍脑袋定的应当根据供热管道长度、管径和水流速度估算或者用热网动态仿真软件标定。3. Matlab代码实现的关键环节3.1 数据准备典型算例的机组参数和负荷曲线先说算例数据怎么准备因为这部分决定了模型能不能跑出合理结果。我的建议是先用一个“能出图、能解释”的小型算例验证模型逻辑再逐步扩大到更复杂的系统。举一个典型配置系统包含1台纯凝火电机组、1台抽汽式CHP机组、1个风电场装机容量较大比如300MW、1台电锅炉、1个蓄热罐、1条与外部电网的联络线。时间范围为24小时分辨率为1小时。纯凝火电机组参数最小出力100MW最大出力300MW煤耗系数a0.0002元/MW²h、b0.15元/MWh、c0.01元/h这里的系数只做示例实际应查对应机组的热耗率曲线爬坡率为30MW/h。CHP机组参数电出力范围100-200MW热出力范围0-150MW运行域用顶点表示比如(-20,0)、(200,0)、(200,150)、(20,150)坐标含义分别是电出力和热出力其中有一个坐标含负电出力原因在于“以热定电”时电出力下方区域不是矩形而是斜线边界负数顶点是几何抽汽区边界的自然延伸代码中必须处理为允许的数学区域。电锅炉参数额定功率50MW热电转换效率0.95运行时最小出力为0即可完全关停。蓄热罐参数容量200MWh最大充放热功率50MW充热效率0.9放热效率0.9初始SOC为100MWh。负荷曲线和风电预测曲线我习惯用一段正弦为基础叠加随机波动。负荷曲线基础值取典型冬季热负荷高的形态电负荷40-60之间的场景不好做建议电负荷峰值取300MW左右热负荷峰值取120MW左右。风电预测曲线取夜间大、白天小的反调峰形态这样的算例能直观展示弃风和消纳的矛盾。3.2 求解器选型与YALMIP环境配置Matlab环境下做MILP求解首先需要安装YALMIP工具包然后准备一个支持MILP的商业或开源求解器。我实测过的组合里最推荐 YALMIP Gurobi其次 YALMIP CPLEX两者在中小规模MILP上都很快。如果没有商业求解器授权可以用 YALMIP SCIP 或 YALMIP GLPK但求解速度会明显变慢大规模算例不建议碰。安装YALMIP很简单去GitHub下载源码把根目录添加到Matlab路径即可。Gurobi则需要申请学术授权或商业授权安装完成后通过gurobi_setup命令配置。然后测试求解器是否可用yalmiptest这个命令会列出可用的求解器。如果输出里出现了Fmincon仅限非线性、Gurobi、CPLEX等条目说明环境配置成功。建模过程中需要反复提到“线性化”。纯凝火电机组的煤耗成本如果是二次函数需要先用分段线性化处理才能进MILP模型。YALMIP里可以用binvar和辅助变量实现或者直接用implies描述分段区间。实操上经常偷懒用一次函数近似或者把二次函数分成4段线性近似误差已经非常小。3.3 核心代码目标函数与约束的实现要点这里给出模型的核心代码骨架展示YALMIP的矩阵化写法。完整可跑的项目代码篇幅长但关键模式可以在这里完整呈现。假设变量定义如下% 时间范围 T 24; % 决策变量 Pg sdpvar(1, T); % 纯凝火电出力 Pchp sdpvar(1, T); % CHP电出力 Hchp sdpvar(1, T); % CHP热出力 Pw sdpvar(1, T); % 风电消纳量 Peb sdpvar(1, T); % 电锅炉耗电功率 Qeb sdpvar(1, T); % 电锅炉产热功率 Pbuy sdpvar(1, T); % 联络线购电功率 SOC sdpvar(1, T1); % 蓄热罐蓄热量 Hch sdpvar(1, T); % 蓄热罐充热功率 Hdis sdpvar(1, T); % 蓄热罐放热功率 u binvar(1, T); % 蓄热罐充放状态变量1为充放然后写目标函数% 成本系数 c_coal 0.15; % 煤耗成本系数元/MWh c_penalty 800; % 弃风惩罚成本元/MWh c_buy [30 30 30 30 30 35 40 45 50 55 55 50 45 45 50 55 55 50 45 40 35 30 25 25]; % 分时电价 % 目标函数 Objective sum(c_coal * (Pg Pchp)) sum(c_buy .* Pbuy) c_penalty * sum(Pw_max - Pw);其中Pw_max是风电预测最大可用出力。这里特别注意弃风量写成Pw_max - Pw而Pw是模型决策变量这样目标函数才能引导模型尽量多消纳风电。约束部分功率平衡和热平衡Constraints []; % 电功率平衡 Constraints [Constraints, Pg Pchp Pw Pbuy Pd Peb]; % 热功率平衡 Constraints [Constraints, Hchp Qeb Hdis Hd Hch];CHP运行域约束用凸包顶点构造% CHP运行域顶点电出力, 热出力 % 这里采用典型的抽汽式CHP凸多边形 V [100 0; 200 0; 200 150; 20 150]; A_chp []; b_chp []; % 利用YALMIP的con2vert或手动构造半空间表达式 % 推荐直接写成多项式不等式 % 具体凸包形式需根据V计算这里给出直线段约束的核心 % 下边界斜线H -1.5 * (P - 100) 0 % 更稳妥做法是用顶点坐标求半空间表达 Constraints [Constraints, Pchp 100, Pchp 200, Hchp 0, Hchp 150]; Constraints [Constraints, Hchp 1.5*(200 - Pchp) 0]; % 根据运行域几何确定这里提醒一个最常见的坑CHP运行域的顶点坐标怎么来的它不应当拍脑袋。正确的做法是从CHP机组的实际技术参数出发——纯凝工况最大电出力对应一个顶点最大抽汽工况对应边界最小凝汽工况对应边界。文献里常用的数值可以直接引用但自己的算例里如果换了机组顶点必须同步换。爬坡约束和蓄热罐% 爬坡约束 Ramp_up 30; Ramp_down 30; for t 2:T Constraints [Constraints, Pg(t) - Pg(t-1) Ramp_up]; Constraints [Constraints, Pg(t-1) - Pg(t) Ramp_down]; Constraints [Constraints, Pchp(t) - Pchp(t-1) 20]; Constraints [Constraints, Pchp(t-1) - Pchp(t) 20]; end % 蓄热罐状态方程 eta_ch 0.9; eta_dis 0.9; for t 1:T Constraints [Constraints, SOC(t1) SOC(t) eta_ch * Hch(t) - Hdis(t)/eta_dis]; Constraints [Constraints, SOC(t) 200, SOC(t) 20]; Constraints [Constraints, Hch(t) 50*u(t), Hdis(t) 50*(1-u(t))]; end蓄热罐特别要注意初始SOC设置SOC(1) 100。如果初始SOC不设模型会利用特殊的初始状态获利比如一开始蓄热罐满罐直接放热满足负荷严重失真。求解部分ops sdpsettings(solver, gurobi, verbose, 0); result optimize(Constraints, Objective, ops); Pw_opt value(Pw); Pg_opt value(Pg); Pchp_opt value(Pchp); Hchp_opt value(Hchp); SOC_opt value(SOC);一个非常关键的细节运行域的下边界推导。如果CHP运行域顶点是(100,0)、(200,0)、(200,150)、(20,150)下边界其实就是从(100,0)到(200,0)的水平线但抽汽式CHP的下边界通常是从(100,0)到(20,150)的斜线表示热出力越大、最小电出力越小。在代码中必须严格按几何求半空间表达式否则会出现“模型优化出CHP在巨大热出力下极小电出力”的物理上不可行结果——虽然数学上不违和但完全不符合CHP特性。3.4 场景扩充把确定性模型改成考虑风电不确定性的随机模型基础模型跑通之后如果想把它扩充为考虑风电不确定性的版本可以在预测功率曲线上叠加场景集。做法是生成N个风电出力场景每个场景对应一组功率平衡约束目标函数改为期望运行成本这就是经典的两阶段随机规划。具体实现Nscen 20; % 生成N个风电场景基于预测值加扰动的正态分布并截断到0到装机容量之间 Pw_scen zeros(Nscen, T); for s 1:Nscen for t 1:T Pw_scen(s,t) min(max(Pw_max(t) randn*10, 0), 300); end end随机模型里决策变量分为两类第一阶段变量机组出力计划、蓄热罐计划不依赖场景第二阶段变量实际风电消纳量、调整后的功率平衡依赖场景。由于YALMIP天然支持多下标变量代码改动并不大。但要注意如果保留机组启停变量随机MILP的求解规模会随场景数线性增长20个场景一般Gurobi还能轻松处理100个场景就开始吃力。3.5 结果可视化关键图表的Matlab绘制方法调度结果不能只看一堆数字必须画图才能看出规律。我按论文常用的输出图梳理一遍。第一张图电功率平衡堆叠图。横轴是时间纵轴是功率用area函数把火电、CHP、风电、购电分别填色堆叠总高度要等于电负荷加电锅炉耗电。这张图能直接看出风电的消纳情况以及CHP在热负荷约束下的“压电”效应。第二张图热功率平衡图。画出CHP热出力、电锅炉热出力、蓄热罐放热/充热、热负荷曲线重点看蓄热罐在什么时段充热、什么时段放热——一般会在晚上风电大时充热白天热负荷高峰期放热。第三张图弃风量柱状图。把每个时段的弃风量画成柱状配合风电预测曲线、实际消纳曲线一起看能明显看出弃风时段集中在夜间以及加了蓄热罐/电锅炉之后弃风削减的幅度。第四张图蓄热罐SOC曲线。横轴时间、纵轴蓄热量看SOC是否在上下限内平滑变化可以直观判断蓄热罐约束是否被合理激活。绘图代码示例figure; bar(1:T, Pw_max - Pw_opt, r); hold on; plot(1:T, Pw_opt, b-o, LineWidth, 1.5); plot(1:T, Pw_max, k--, LineWidth, 1.2); legend(弃风功率, 实际消纳风电, 风电预测); xlabel(时段/h); ylabel(功率/MW);4. 常见问题与排查技巧实录4.1 求解器报错或返回不可行解怎么办做MILP调度模型遇到最频繁的问题就是Infeasible problem。我的排查顺序是这样的先检查等式约束特别是电功率平衡和热功率平衡——很多不可行都出现在“某些时段所有可调资源加起来不够满足负荷”这时模型确实无解。为了快速定位不可行原因有一个非常实用的技巧把目标函数暂时改成常数0然后逐个约束取消测试或者用YALMIP的check命令检查每条约束的残差check(Constraints)输出里会标明哪条约束残差最大、哪个变量越界。这个方法比我反复注释代码快得多。另一个常见原因是蓄热罐的SOC范围设置有误。比如容量200MWh你却把SOC约束写在0 SOC 200而初始SOC正好100但充放热功率上限500.9的效率一折算可能出现“时间顺序上无论如何也跳不过某约束”的隐性矛盾。这种检查脚本很难发现需要实物推演一遍。4.2 结果明显不合理的四类典型情况如果调试后模型有解但结果看着不对那大概率是以下四类问题。第一弃风量等于0。这不一定算结果“错”但如果你设计的算例里风电高发时段和热负荷高峰重叠弃风却完全为零就要怀疑是不是弃风惩罚系数设太大导致模型不敢弃风逼迫机组运行在极端工况。检查一下CHP运行点是否跑到了边界外、爬坡是否勉强卡上限如果是把弃风惩罚系数适度调低。第二蓄热罐状态曲线出现“锯齿”。充热、放热每时段都切换说明目标函数里缺少对充放热切换的经济惩罚或逻辑限制模型在用蓄热罐做“套利”。这时候加一条充放热状态切换成本或者把SOC的平滑性纳入目标函数锯齿就能消除。第三CHP运行点跑出了运行域。这基本可以断定为运行域约束写错了。最常见的手误是把CHP上下限分开写成Pchp 100, Pchp 200, Hchp 0, Hchp 150但这根本不能刻画“热出力越大最小电出力越小”的耦合关系。必须用多边形半空间约束完整描述。第四结果里电锅炉整日满发。这说明热负荷主要由电锅炉满足CHP热出力被压低。这种情况虽然也能消纳风电但不经济。需要检查电锅炉的运行成本是否过低、CHP煤耗系数是否设得过高。调参本质上权衡热侧灵活性与经济性也是调度的核心。4.3 求解速度慢的处理办法中小规模算例比如24时段、5台机组、20个场景Gurobi一般一秒内出结果。如果求解时间超过10秒先检查模型里是不是有大量不必要的二进制变量。一个常见坑是给每个时段的蓄热罐充放热状态都定义了0/1变量其实24小时每个时段都充放切换本身就不合理完全可以引入“最小充放持续时间”约束减少模型的组合爆炸。或者如果只关心SOC的最终变化可以把充放热状态变量做成连续变量加互补约束虽然理论上仍是非凸但实际求解中收敛性尚可。另外YALMIP会把约束逐条展开如果你把约束写在for循环里变量维度很大时YALMIP的解析开销会明显增加。建议对循环约束用矩阵形式一次生成比如爬坡约束用索引矩阵乘法能快出好几倍。4.4 参数灵敏度分析的实用思路做研究的时候单跑一种算例说明不了问题通常需要分析关键参数变化对结果的影响。我碰到的项目里经常需要做三类灵敏度分析。一是弃风惩罚系数的灵敏度。从200到2000元/MWh分档扫描画出弃风率变化曲线能看到明显的阈值效应。二是蓄热罐容量的灵敏度。从0到500MWh分档扫描看蓄热罐容量每增加多少MWh弃风率下降多少这个结果特别适合写进论文的“消纳能力提升分析”章节。三是电锅炉功率上限的灵敏度。考察电锅炉从20MW到80MW变化时系统总成本和弃风量的Pareto变化趋势。实现上就是“参数循环 重新求解 记录关键指标”核心代码不超过30行但能产出的图表远远多于一票到底的单一算例。4.5 一套准确的数据校验方法最后分享一个我自己踩过很多坑后总结出的校验流程。拿到任意一组调度结果先按公式手动检查几个关键时段挑热负荷最大的那个时段手动计算CHP热出力加上电锅炉产热、蓄热罐放热是否严格等于热负荷挑风电最大的那个时段手动核对电功率平衡里风电消纳量是否被其他电源挤占再检查蓄热罐SOC的连续性——用上一时段的SOC代入状态方程能否算到当前时段的SOC值。这套校验不用花多少时间但能避免大量的“结果看起来合理但实际上是错的”尴尬情况。尤其是投论文或者做项目汇报前务必做一遍这个检查。我个人的实际感受是电热综合能源系统的日前经济调度模型难的不是某一个约束的数学表达而是把物理逻辑完整地翻译成约束集合并保证翻译过程不丢失关键耦合信息。顺着“为什么电热要联合” → “模型怎么搭” → “代码怎么落” → “问题怎么查”这条路走下来你收获的不只是一个能跑的模型还有一套能迁移到后续研究方向比如加入碳捕集、需求响应、配电网潮流的方法论。