做电热综合能源系统日前经济调度研究时我一直在琢磨一件事怎么把可再生能源消纳的压力体现在优化模型里这套用Matlab代码实现的电热综合能源系统日前经济调度模型就是在这件事上折腾了几个月的结果。它不是那种只能跑通文档的玩具demo而是能从负荷数据、装机参数一路算到24小时出力计划、弃风弃光指标的完整闭环。对于正在做综合能源系统优化、电力系统经济调度的研究生或者工程师这套思路和代码可以直接拿来改、拿来用。我先把话挑明这个模型不是骡子式的“电/热两套系统各算各的”而是把电锅炉、储热罐、热电联产机组放在同一个优化框架里以24小时为周期做经济调度。目标是花钱最少同时让风电、光伏尽量多上网。下面我把从模型设计到Matlab实现、再到调参踩坑的过程完整拆开讲。1. 模型思路电热耦合不是简单加个约束1.1 为什么必须把“电”和“热”放在同一个优化框架里如果你单独做电力调度你会忽视供热机组的热出力约束单做热力调度你又看不到电网的实时电价信号。现实中的园区能源系统最常见的场景是这样冬季采暖季热电联产机组CHP为了满足热负荷锅炉必须维持较高出力发电功率被“绑着”往上走。而这时候如果恰好是夜间风电大发电网那侧消化不了这么多电风机只能弃掉。这就是典型的“以热定电”导致弃风。解决办法不是关掉热负荷而是引入电锅炉、储热罐这些“解耦”设备。电锅炉让电力需求在供热侧多一个出口储热罐则把多余热量存起来等电价高的时候再放出来。这样一来CHP机组的热出力不必紧跟热负荷电出力也就有了下调空间。整个系统的运行灵活性上来了可再生能源消纳的通道自然就打开了。所以电热综合调度不是拍脑袋“为了综合而综合”是电网灵活性不足倒逼出来的需求。我这个模型把电、热两侧的平衡方程放在同一个线性规划里调度对象包括常规发电机组、CHP机组、风电、光伏、电锅炉和储热装置。目标是在满足电热负荷的前提下最小化总运行成本同时用弃风弃光惩罚项引导模型优先消纳新能源。这种处理方式最直接也最容易在论文里讲清楚。1.2 日前经济调度到底在算什么事日前调度说白了就是提前一天决定每台设备在未来24小时每个时段的出力计划。我遇到不少初学者分不清“经济调度”和“机组组合”经济调度假设机组启停状态已经给定只优化各台设备的出力大小而机组组合还要额外决定哪些机组在哪些时段开机、哪些停机中间会引入0/1整数变量。我选择先做经济调度是因为它能把问题本质讲清楚——所有约束都是线性的目标函数也是线性的调参、调试、可视化都非常顺手。等你把经济调度跑通再往上加0/1变量做机组组合、加不确定场景做随机优化就是在同一个骨架上扩展。模型的时间分辨率取1小时全天24个点。这个粒度对研究“消纳”问题足够了风电低谷波峰、电价低谷时段都能捕捉到。如果你要做分钟级调度模型规模会变大但约束本质完全一样。我没在模型里考虑电网络潮流默认整个园区是一个节点所有设备挂在同一母线上。这样做的好处是省去节点电压、线路潮流这些非线性项把注意力集中在“电-热耦合”和“可再生能源消纳”这两个核心问题上。等需要精细化网架分析再叠加潮流约束也不迟。2. 目标函数与约束条件建模每个公式都对应一句人话2.1 目标函数怎么写才合理我的目标函数分四块外购电成本、机组燃料成本、弃风弃光惩罚成本、设备运维成本。外购电是从上级电网买电的费用用分时电价乘以购电功率逐时求和。机组燃料成本包括燃气轮机和CHP消耗天然气的费用可以近似表示成出力的线性函数。经济调度研究里常用二次曲线但实际参数标定很麻烦线性近似在方案对比阶段完全够用。设备运维成本很简单每台设备单位出力的维护单价乘上出力即可不是影响结果的核心项。最关键的是弃风弃光惩罚成本。风电光伏的预测出力是当天模型能消纳的上限实际出力是决策变量预测值和实际值之差就是弃电量。惩罚项的形式是C_pen * (P_wind_pred(t) - P_wind(t) P_pv_pred(t) - P_pv(t))这个惩罚系数要设得够大比如500元/兆瓦时这样模型“舍不得”弃风。如果系数设成10元/兆瓦时而电价低谷时购电成本才0.3元/度模型宁可让风机少发然后去电网买电导致弃风率虚高这显然不符合“消纳优先”的预期。完整目标函数写成文字公式就是min Σt [ c_el(t)*P_buy(t) c_gas*P_gt(t) c_chp_el*P_chp(t) c_chp_h*H_chp(t) c_pen*(P_wind_pred(t)-P_wind(t)P_pv_pred(t)-P_pv(t)) c_om_eb*P_eb(t) ]这里P_gt是燃气轮机电出力P_chp、H_chp分别是CHP电、热出力P_buy是购电功率P_eb是电锅炉耗电功率。把这个目标写进Yalmip时注意向量化和系数匹配如果这些决策变量都是1×24的sdpvar向量那么目标函数用sum(系数 .* 变量)这种形式一行代码搞定不需要写for循环。2.2 电功率平衡和热功率平衡能量守恒是硬约束任何调度模型最基本的约束就是功率平衡。电功率平衡如下P_gt(t) P_chp(t) P_wind(t) P_pv(t) P_buy(t) - P_eb(t) P_load_e(t)其中P_load_e是电负荷。这里P_buy是非负变量假设系统只购电、不向电网倒送电。有些园区允许余电上网那就额外加一个P_sell变量等式右侧再减去它。热功率平衡则是H_chp(t) H_eb(t) H_dis(t) - H_char(t) H_load(t)H_eb是电锅炉产热功率H_dis、H_char分别是储热罐的放热和充热功率。如果系统没有储热罐这两个变量恒为0。注意热功率平衡里的H_eb和电平衡里的P_eb通过电锅炉效率耦合在一起H_eb eta_eb * P_eb效率通常取0.9~0.95。电锅炉运行时会有热损耗这部分损耗没有出现在热平衡里而是体现在效率系数里。模型里不需要额外再扣一次热损失别重复计入。除了平衡约束我还会给常规机组加爬坡约束-ramp_gt P_gt(t1) - P_gt(t) ramp_gt如果不加爬坡模型会给出跳跃性的出力曲线实际机组根本跟不上。但爬坡约束也不能设得太紧否则可能直接导致无解这一点后面调试部分会再谈。2.3 热电联产机组的可行域多边形约束怎么处理CHP机组不是任意给定电出力和热出力都能运行。它的运行区间是一个多边形由厂家给出的电-热运行范围决定。常见描述包括电出力上下限、热出力上下限、以及“电出力随热出力变化的斜率约束”。比如某台背压式CHP电出力和热出力近似线性关系可行域可以写成P_chp_min ≤ P_chp ≤ P_chp_maxH_chp_min ≤ H_chp ≤ H_chp_maxP_chp k1 * H_chp ≥ C1P_chp k2 * H_chp ≤ C2具体系数由机组特性曲线标定。手动输入这些不等式很容易方向写反我的做法是把多边形顶点坐标直接读进来然后用凸包工具生成约束矩阵A和右侧向量b。Matlab里可以用convhull函数把顶点排序再逐条边计算法向量生成不等式。这样即使换一台不同参数的机组只需要替换顶点矩阵就行代码不用动。这个技巧在要多机组对比时特别管用。2.4 电锅炉、储热罐等灵活性资源怎么建模储热罐是模型里最有意思的部分。它的状态量是蓄热量S(t)递推关系为S(t1) S(t) eta_c * H_char(t) - H_dis(t) / eta_d - epsilon * S(t)eta_c和eta_d分别是充热、放热效率epsilon是散热损失系数。量纲要统一S用MWhH_char、H_dis用MW时间步长是1小时所以功率乘以1小时就是热量变化量。如果时间步长不是1小时记得乘上对应间隔。运行约束还包括蓄热量上下限、充放热功率上限0 ≤ S(t) ≤ S_max0 ≤ H_char(t) ≤ H_char_max0 ≤ H_dis(t) ≤ H_dis_max还有一个日循环约束为了场景能够重复执行比如明天继续用同样策略通常要求一个调度周期结束时的蓄热量等于初始值S(1) S(25)或S(end) S(1)这个约束在夜间谷时段强制储热罐至少恢复到初始状态否则模型会把储热罐里最后一点热量在峰时段全部放光第二天无热可用。我一开始没加这个约束结果是储热罐SOC在第24小时直接掉到下限实际运行根本无法复现。加上周期约束后结果就正常多了。3. Matlab实现选对工具代码能省一半时间3.1 建模工具选型Yalmip加求解器是黄金组合Matlab里做优化模型的路子不少手写矩阵调linprog、用CVX、用Gurobi自带接口或者用Yalmip。我个人强推Yalmip。这玩意儿相当于一个“翻译官”让你用接近数学表达式的语法定义变量和约束然后自动翻译成求解器需要的标准形式。配合Gurobi或CPLEX求解速度快到飞起。如果你暂时没有商业求解器Yalmip默认调用Matlab的linprog小规模LP也能跑就是慢一些。我用的是Gurobi学术版处理带大量稀疏矩阵的LP一直很稳而且它在Yalmip里的配置就一行sdpsettings(solver,gurobi)。别用CVX做这件事。CVX对约束的拼接和指数锥问题的支持很优雅但它的LP建模方式偏“规范型”修改起来没有Yalmip自由。Yalmip的约束对象可以像搭积木一样往一个Constraints数组里塞调试时还能用showconstr查看某条约束体验舒服太多。3.2 变量定义、约束组装的高频技巧在Yalmip中连续变量用sdpvar定义。除了矩阵还可以定义向量。假设系统有24时段那么每个设备都定义成1×24的行向量P_wind sdpvar(1, 24); % 风电实际出力 P_pv sdpvar(1, 24); % 光伏实际出力 P_chp sdpvar(1, 24); % CHP电出力 H_chp sdpvar(1, 24); % CHP热出力 P_gt sdpvar(1, 24); % 燃气轮机电出力 P_buy sdpvar(1, 24); % 购电功率 P_eb sdpvar(1, 24); % 电锅炉耗电 H_eb sdpvar(1, 24); % 电锅炉产热 S sdpvar(1, 24); % 储热罐蓄热量 H_char sdpvar(1, 24); % 储热罐充热功率 H_dis sdpvar(1, 24); % 储热罐放热功率约束尽量用向量形式写不要在for循环里一条条加。Yalmip对向量约束支持得很好。比如风电出力上下限直接写C [C, 0 P_wind P_wind_pred, 0 P_pv P_pv_pred];这比循环里24次C [C, P_wind(t) P_wind_pred(t)]要快得多而且不容易维度出错。但如果你需要逐时段赋不同的系数比如分时电价对应的购电成本可以把系数向量乘上去仍然是一行。3.3 核心代码框架示例给一个浓缩版主干流程是读数据、定义变量、写目标、写约束、求解、取结果。%% 数据准备 T 24; P_load_e data.P_load_e; % 电负荷 1x24 H_load data.H_load; % 热负荷 1x24 P_wind_pred data.P_wind_pred;% 风电预测 1x24 P_pv_pred data.P_pv_pred; % 光伏预测 1x24 c_el data.c_el; % 分时电价 1x24 eta_eb 0.95; % 电锅炉效率 P_eb_max 1; % 电锅炉最大功率 MW S_max 4; H_char_max 1; H_dis_max 1; % 储热罐参数 %% 变量定义 clear P_gt P_chp H_chp P_wind P_pv P_buy P_eb H_eb S H_char H_dis P_gt sdpvar(1,T); P_chp sdpvar(1,T); H_chp sdpvar(1,T); P_wind sdpvar(1,T); P_pv sdpvar(1,T); P_buy sdpvar(1,T); P_eb sdpvar(1,T); H_eb sdpvar(1,T); S sdpvar(1,T); H_char sdpvar(1,T); H_dis sdpvar(1,T); %% 目标函数 c_chp_el 0.35; c_chp_h 0.15; % CHP电/热单位成本 c_gas 0.22; c_pen 500; % 燃气轮机成本、弃能惩罚 obj sum(c_el .* P_buy) c_gas*sum(P_gt) c_chp_el*sum(P_chp) c_chp_h*sum(H_chp) ... c_pen*(sum(P_wind_pred - P_wind) sum(P_pv_pred - P_pv)) ... 0.01*sum(P_eb); %% 约束 C []; % 电/热平衡 C [C, P_gt P_chp P_wind P_pv P_buy - P_eb P_load_e]; C [C, H_chp H_eb H_dis - H_char H_load]; % 新能源消纳 C [C, 0 P_wind P_wind_pred, 0 P_pv P_pv_pred]; % 电锅炉 C [C, H_eb eta_eb * P_eb, 0 P_eb P_eb_max]; % CHP可行域 % 这里假设CHP电出力0.5~3MW热出力0.2~2.5MW并用两个耦合约束模拟多边形 C [C, 0.5 P_chp 3, 0.2 H_chp 2.5]; C [C, P_chp 0.6*H_chp 0.8, P_chp 0.8*H_chp 4.5]; % 燃气轮机约束 C [C, 0.3 P_gt 2]; % 储热罐 C [C, S(2:T) S(1:T-1) eta_c*H_char(1:T-1) - H_dis(1:T-1)/eta_d - epsilon*S(1:T-1)]; C [C, S(end) S(1)]; C [C, 0 S S_max, 0 H_char H_char_max, 0 H_dis H_dis_max]; % 爬坡约束 C [C, -0.5 diff(P_gt) 0.5]; %% 求解 ops sdpsettings(solver,gurobi,verbose,1); optimize(C, obj, ops); %% 输出结果 P_chp_opt value(P_chp);注意代码里的S(2:T)和S(1:T-1)这种向量写法Yalmip会自动展开成T-1条等式非常简洁。diff(P_gt)是相邻时段的差值也支持直接写进约束。如果报错多半是向量方向问题把列向量转成行向量就好。3.4 结果可视化把调度曲线画成图优化完成后用value()取出各变量数值就可以画图了。我习惯画四张图电功率平衡堆叠面积图、热功率平衡堆叠面积图、弃风弃光柱状图、储热罐SOC曲线。画图核心思路是用area函数把不同设备出力叠起来和负荷曲线对比。figure; area(1:T, [value(P_wind); value(P_pv); value(P_chp); value(P_gt); value(P_buy)]) hold on plot(1:T, P_load_e, k-, LineWidth, 2); legend(风电,光伏,CHP,燃气轮机,购电,电负荷);注意面积图的输入矩阵是“时段×设备”所以要对设备出力向量做转置。area会按列堆叠顺序要和legend对应。储热罐SOC画折线就行。导出图片时我一般用set(gcf,Color,w)避免灰色底再exportgraphics(gcf,result.png,Resolution,300)直接满足论文插图要求。4. 实操过程与核心环节实现从数据准备到结果复现4.1 基础输入数据的准备没有靠谱的数据模型就是个空壳。我这儿给一套典型的参数适合自己跑着玩假设一个园区系统电负荷峰值8MW热负荷峰值5MW风力装机4MW光伏装机2MWCHP装机电功率3MW/热功率2.5MW燃气轮机装机2MW电锅炉1MW储热罐容量4MWh最大充放热功率1MW。分时电价采用峰谷两段式峰时段10:00-15:00、18:00-21:000.8元/kWh谷时段23:00-7:000.3元/kWh其他时段0.5元/kWh。风电预测曲线做成“夜间大、白天小”的形状光伏则是“中午大、两端小”。这些数据我放在Excel里Matlab用readtable读出来再转成行向量。设备参数可以整理成一张表设备参数数值CHP机组电出力下限/上限0.5/3 MWCHP机组热出力下限/上限0.2/2.5 MWCHP机组电成本系数0.35 元/kWhCHP机组热成本系数0.15 元/kWh燃气轮机出力下限/上限0.3/2 MW燃气轮机单位燃料成本0.22 元/kWh电锅炉电热转换效率0.95电锅炉最大电功率1 MW储热罐容量4 MWh储热罐最大充/放热功率1 MW储热罐充/放热效率0.9储热罐散热损耗0.05/h弃风/弃光惩罚单价500 元/MWh这些参数不用太较真重点是量纲统一功率用MW热量也用MW表示热功率蓄热量用MWh时间步长1小时。如果热负荷的单位是GJ/h记得先除以3600换算成MW否则等式平衡就完全对不上。单位出错是新手最容易掉进去的坑。4.2 模型搭建中的几个关键调试节点代码写完后我先做“空跑测试”把风电、光伏预测值全设为0电锅炉和储热罐也都不能动让模型退化成纯购电加燃气轮机加CHP的传统系统。这一步验证的是基础平衡约束和机组约束是否写得对。如果空跑都无解那绝对不是新能源约束的问题而是平衡等式、爬坡或可行域哪里“打架”了。第二步把风电预测按比例放大从0.5倍逐步调到2倍观察弃风量变化。弃风量应当随预测增加而增加并且增幅合理。如果预测功率明明很高但弃风量为负那一定是你把P_wind_pred - P_wind写反了。第三步临时把储热罐的周期约束去掉看看SOC终值会掉到哪里。如果终端SOC和初值差距大说明周期约束确实在起作用。这三个测试都通过模型基本稳了剩下的就是调惩罚系数和爬坡参数。4.3 典型调度结果的解读以冬季典型日为例调度结果大概是这样的深夜电价低谷、风电预测很高此时CHP机组尽量压低电出力但要满足热负荷所以热出力维持下限不足的热量由电锅炉在谷电时段补充同时储热罐充电把热量存起来。白天进入峰价时段后CHP电出力回升储热罐放热替代部分电锅炉出力降低购电成本。到晚上用电高峰购电成本高系统会启动燃气轮机发电配合储热罐放热实现“削峰填谷”。从指标上看如果没有电锅炉和储热罐传统的“以热定电”模式下夜间弃风率可能到10%以上加入电锅炉和储热罐后弃风率能压到2%以内。你还可以特意跑一个“无电锅炉、无储热罐”的对照组和完整模型对比两张图放在一起电-热灵活资源的价值一眼就能看出来。这种对比特别适合写论文的算例部分。5. 常见问题与排查技巧实录5.1 求解器返回无解或状态码inf最常见的原因是约束集合不可行——某时段电热平衡根本无法同时满足。排查方法是把目标函数暂时屏蔽只求解可行性问题optimize(C, [], ops)然后看check(C)的结果哪条约束残差大于1e-4问题就在哪。我遇到过几次无解都是储热罐周期约束和SOC下限矛盾比如初始SOC已经低到1MWh而容量下限设置成2MWh自然无解。把初值与容量下限统一或者在周期约束中允许一定容差问题就消失了。5.2 弃风量异常大先查惩罚系数如果你把弃风惩罚系数设为10元/MWh而峰时电价0.8元/kWh折合800元/MWh模型在“弃风少发电”和“购电满足负荷”之间做选择时很可能选择弃风因为弃风惩罚太便宜了。建议惩罚系数至少设为购电价上限的5到10倍我默认设500元/MWh。另外风电出力的上限向量必须是预测值写成P_wind P_wind_pred可别把预测值直接当成固定值赋给变量——那就不是在“消纳决策”了而是在跑固定出力场景。5.3 Yalmip报错“VARIABLE DOESNT EXIST”或维度不匹配这类问题十有八九是变量定义和作用域搞混了。确保所有sdpvar在同一个workspace里如果写了函数封装变量必须在函数内部重新定义。另一个高频坑是行列方向不一致Yalmip对1×24和24×1非常敏感从Excel读进来的列向量忘了转置约束维度对不上优化直接报错。我的习惯是读入数据后统一加一句data data(:);强制变行向量后面所有变量定义都用行向量基本不会再出维度问题。5.4 求解时间太长如何优化性能24小时单节点模型正常求解时间在毫秒到秒级。如果特别慢先在sdpsettings里把verbose,0关掉减少输出。其次检查是不是在约束里用了大量循环尽量用向量化写法。Yalmip约束对象的数量本身会影响求解器前处理速度把同类型约束一次性拼接不要拆成几十条独立约束。最后如果模型规模特别大比如几百个节点、几千条约束可以把求解器内部参数调一下比如Gurobi的LPWarmStart开成2或者把ShowProgress设为0。我做过对比同样规模模型从循环逐条加约束改成向量化求解时间可以从200ms降到80ms效果明显。5.5 储热罐SOC曲线疯狂振荡如果SOC曲线要么冲到上限要么贴下限大概率是充放热功率约束没写对或者效率和热损耗系数设得太极端。调试时先把散热损耗设为0看SOC曲线是否平滑再把充放热效率调成1排除效率把状态量“吞掉”的可能。另外如果储热罐同时充放正负交替大概率是你没有禁止同步充放。LP模型里目标函数是成本最小化通常不会出现既充又放的浪费行为但如果惩罚项导致某种利润模式出现同步充放就需要加0/1变量做互斥约束。我的模型由于目标函数只惩罚弃能、不奖励套利没出现这个问题如果你要做电价套利就要小心这个坑。5.6 电力平衡中有缺口但检查所有约束都是对的这时候去看看电锅炉。电锅炉耗电功率进了电平衡产热功率进了热平衡两者通过效率等式连接这一步最容易漏。我见过好多次电平衡右侧忘了加-P_eb热平衡右侧又加了H_eb结果电负荷缺口恰好等于电锅炉耗电模型还“自我感觉良好”。把平衡等式拉出来逐项核对变量正负号尤其关注那些出现在两个平衡方程里的设备。结尾做完这套电热综合能源系统日前经济调度模型我最大的体会是建模过程最花时间的不是数学公式而是怎么把实际物理设备的运行区间“翻译”成约束矩阵。CHP可行域、储热罐周期约束、弃风惩罚这些细节一个地方没处理好结果就会疯。解决的办法就是多做几个算例多画几张图对着结果看设备有没有“违规出工”。最后再分享一个小技巧在Matlab里写大矩阵约束时多用repmat和kron构造系数矩阵再配合Yalmip的向量化约束整个模型的调试体验会舒服很多。这个模型现在只是个确定性调度版本如果后续要做随机优化或者分布鲁棒把风电预测误差场景接入本质上就是在这个骨架上加一堆场景复制思路是相通的。