简介面向电气工程与能源领域毕业设计与课题研究的电-气-热综合能源系统耦合调度/优化调度源码包针对新能源大规模接入导致主网下网功率波动加剧的问题提供计及新能源出力不确定性的协同优化建模与求解实现。程序采用动态场景法刻画出力不确定目标为最小化系统运行成本与主网下网功率波动并考虑温控负荷调节能力、配电网交流潮流及天然气网潮流约束经分段线性化与二阶锥松弛转化为混合整数二阶锥规划附夏季、冬季算例仿真验证气网惯性平抑波动和温控负荷降低成本的综合效果。压缩包共25个文件以20个Matlab源程序为主体包含模型与数据说明文本、可参考论文、Word解析及数据表格整体5.12MB结构清晰便于对照运行与二次开发。已有170人学习下载适合需要快速掌握综合能源耦合调度模型、开展仿真验证的本科生、研究生或相关从业者。1. 电-气-热综合能源系统耦合调度、优化调度源程序先复现模型再谈跑通拿到“电-气-热综合能源系统耦合调度、优化调度”这类源程序我的第一反应不是直接打开主函数按 F5而是先把配套论文里的系统拓扑图找出来。原因很简单这类源程序的价值不在“能跑出一个成本数字”而在数字背后电、气、热三个网络是怎么被耦合到一起的。它属于综合能源系统优化调度里最典型的模板——用 CHP、P2G、燃气锅炉、储热罐把三种能源载体连起来做日前 24 小时或 96 时段的协同经济调度附带论文和 WORD 解析的作用就是把模型讲清楚让代码不再是黑匣子。适合正在写这个方向论文的研究生也适合想把手头程序与论文图表逐点核对、做场景扩展的从业者。前提是你愿意把源程序当成模型来改而不是当成一个只出结果的工具。2. 耦合调度模型先立住能量枢纽映射、三类网络约束与目标函数怎么选2.1 能量枢纽建模CHP、P2G、燃气锅炉、储热罐在程序里的接口关系翻开源程序时第一眼看到的不是物理网络而是一组接口变量。无论程序写成什么样设备层的耦合关系都是固定的CHP 消耗气功率输入同时输出电和热P2G 消耗电输出气燃气锅炉消耗气输出热储热罐负责热功率的时序平移。把这四个设备的接口关系列出来程序骨架就清楚了设备输入变量输出变量效率/容量接口CHPF_chp气功率P_chp电、Q_chp热eta_ge、eta_ghP2GP_p2g电G_p2g气eta_p2g燃气锅炉F_gb气Q_gb热eta_gb储热罐H_sto充放热E_sto能量状态sto_cap、max功率我一般会先把论文设备参数表转成一个结构体统一管理而不是让参数散落在代码各个角落% 把论文设备参数表转成结构体统一管理避免散落各处的魔法数字 para.chp_eta_ge 0.35; % CHP 电效率来自论文设备表 para.chp_eta_gh 0.45; % CHP 热效率 para.p2g_eta 0.60; % P2G 电转气效率 para.gb_eta 0.85; % 燃气锅炉热效率 para.sto_cap 60; % 储热罐容量MWh para.sto_max_ch 10; % 最大充放热功率MW这里的效率耦合本质上是在做“单位统一后的功率折算”不是能量守恒的直接形式。尤其要注意气侧程序里变量 F_chp、F_gb、G_p2g 大多统一用“MW 热值”而不是 m³/h否则后面购气成本会差一个热值系数。固定效率是多数论文程序的默认做法虽然实际机组变工况效率会有变化但先把固定效率版本跑通再替换成线性化效率曲线是更稳妥的路径。2.2 电网-气网-热网三类约束变量边界与量纲处理是共同难点耦合调度程序里最难的不是设备而是三类网络约束各自有一套变量和单位逻辑。电网约束普遍用节点功率平衡加直流潮流近似变量是母线相角和线路潮流但不少论文简化成“从上级电网购电”的单节点模型程序里只有一个 P_buy 变量。如果代码里出现了相角变量和线路潮流矩阵说明它真的在算网架这时要注意线路容量约束的编号是否与论文拓扑图一致。气网约束的核心是节点流量平衡加 Weymouth 方程表达的是管道流量与两端气压平方差的关系。这是整个程序里最容易让求解器翻车的地方因为平方项直接放进 MIP 是非线性的后面第 4 章会专门讲处理办法。热网稍微复杂一点水力方程、温度降落方程、节点温度混合都能写但多数论文源程序只保留热功率平衡和储热动态把供热系统当“热母管”处理。所以先看清楚你手上这段代码做了哪个粒度的建模再决定要不要补热网细节。三类网络的单位不同是调试时最容易踩的坑建议所有加减运算之前先折算到公共单位载体常用原始单位折算到公共单位电kW / MW1 MW 1000 kW天然气m³/h、万 m³/h按热值约 35.7 MJ/Nm³ 折算为 MW热GJ/h、Gcal/h1 GJ/h ≈ 0.2778 MW2.3 目标函数与日前调度框架成本型还是碳型决定程序改动量目标函数决定这个程序做出来是拿来算什么。最常见的配套论文是“运行成本最小”公式大概是购电成本加购气成本加设备运维成本也有一批论文写的是“碳排放最小”或“弃风弃光最小”。程序上的差异就是目标函数第一项和第二项的区别写起来像这样% 成本型与碳型目标函数骨架单位统一到 MW 和 元/MWh 后直接相加 base_cost sum(1000 * c_ele .* P_buy) sum(1000 * c_gas .* G_buy); carbon_cost carbon_price * (sum(P_buy) * ele_co2 sum(F_chp F_gb) * gas_co2); obj base_cost carbon_cost; % 碳价权重从论文参数表抄需要注意 c_ele 和 c_gas 的单位不同前者按元/kWh 给后者可能按元/m³ 给。如果论文里的购气价写的是 2.5 元/m³要先用热值折算成元/MWh再进目标函数。改动量最大的不是目标函数本身而是为了配合“碳型目标”额外引入的碳排放系数和设备启停变量。先确定你手上程序是哪类目标再动代码否则你会在调试时发现结果总是跟论文对不上。3. 把论文算例变成可运行程序参数表落地、最小调度代码与求解器调用3.1 先按这个顺序读论文填参数拓扑图、参数表、负荷曲线拿到源程序和配套论文我建议按三步走顺序不要反。先看系统拓扑图确定电网几个节点、气网几个节点、热网是不是成网这决定了程序里变量矩阵的大小。再看设备参数表把 CHP 容量、效率、P2G 容量、储热罐容量抄进 para 结构体。最后看日负荷曲线注意横坐标到底是 24 个小时还是 96 个 15 分钟断面很多程序跑出来曲线对不上问题不在算法在这里。论文里负荷数据常常只有图没有表格这时需要从图上手工数字化。我的习惯是把电、热、气三条曲线分别存成 CSV数据文件内容我习惯的命名load_ele.csv24 或 96 点电负荷单位 MWload_heat.csv热负荷曲线单位 MWload_gas.csv气负荷曲线折算成 MW 热值如果论文实在没给热负荷曲线常见处理是按论文给定的热电比从电负荷缩放出来。不要自己凭空造一条平滑曲线后面你做图表复现时负荷形状对不上会非常难排查。3.2 最小可跑的日前调度代码以 CHP P2G 储热为例跑通一个能量枢纽这份代码没有建电网潮流、气网管道和热网水力只保留能量枢纽级的功率平衡是这个方向源程序的最小能跑骨架。绝大多数带网络约束的论文程序本质上就是在骨架上加变量、加约束行% iehs_dispatch_demo.m % 电-气-热综合能源系统日前调度最小示例MATLAB YALMIP Gurobi/CPLEX clear; clc; yalmip(clear); %% 1 参数区负荷与分时价格 T 24; Pd [80 75 70 68 72 80 92 105 120 130 132 128 ... 125 127 130 132 128 115 105 95 88 82 78 75]; % 电负荷 MW Qd [45 44 43 41 42 44 46 49 53 56 55 53 ... 52 52 53 55 52 48 45 43 42 41 40 39]; % 热负荷 MW Gd [60 58 56 52 50 48 55 67 80 90 94 90 ... 87 84 85 86 83 78 70 62 58 55 53 52]; % 气负荷 MW热值 c_ele 550 * ones(1, T); % 基础购电价元/MWh c_ele(1:6) 350; c_ele(12:18) 850; % 谷段 0-6点峰段 11-18点 c_gas 320 * ones(1, T); % 购气价元/MWh热值 %% 2 设备效率与容量参数 eta_ge 0.35; eta_gh 0.45; % CHP 电/热效率 eta_p2g 0.60; % P2G 效率 eta_gb 0.85; % 燃气锅炉效率 P_chp_max 120; P_chp_min 20; % CHP 电出力上下限 P_p2g_max 50; % P2G 最大耗电功率 Q_gb_max 80; % 燃气锅炉最大热出力 E_sto_max 60; E_sto_0 30; % 储热罐容量与初始能量 h_sto_max 10; % 最大充放热功率 %% 3 决策变量 F_chp sdpvar(1, T); % CHP 输入气功率 P_chp sdpvar(1, T); % CHP 电出力 Q_chp sdpvar(1, T); % CHP 热出力 P_p2g sdpvar(1, T); % P2G 耗电 G_p2g sdpvar(1, T); % P2G 产气 F_gb sdpvar(1, T); % 燃气锅炉耗气 Q_gb sdpvar(1, T); % 燃气锅炉产热 P_buy sdpvar(1, T); % 向上级电网购电 G_buy sdpvar(1, T); % 从气网购气 E_sto sdpvar(1, T1); % 储热罐能量状态 H_sto sdpvar(1, T); % 放热为正充热为负 %% 4 约束 C []; C [C, P_chp eta_ge * F_chp, Q_chp eta_gh * F_chp]; C [C, P_chp_min P_chp P_chp_max, 0 Q_chp 120]; C [C, G_p2g eta_p2g * P_p2g, 0 P_p2g P_p2g_max]; C [C, Q_gb eta_gb * F_gb, 0 Q_gb Q_gb_max]; C [C, E_sto(2:T1) E_sto(1:T) - H_sto]; % 储热动态 C [C, 0 E_sto E_sto_max]; C [C, E_sto(1) E_sto_0, E_sto(T1) E_sto_0]; % 日循环约束 C [C, -h_sto_max H_sto h_sto_max]; C [C, Pd P_p2g P_buy P_chp]; % 电功率平衡 C [C, Qd Q_chp Q_gb H_sto]; % 热功率平衡 C [C, Gd G_buy G_p2g - F_chp - F_gb]; % 气功率平衡 C [C, 0 P_buy, 0 G_buy]; %% 5 目标购电成本 购气成本储热动作尽量平缓 obj sum(1000 * c_ele .* P_buy) sum(1000 * c_gas .* G_buy) ... 1e-3 * sum(abs(H_sto)); %% 6 求解与结果 ops sdpsettings(solver, gurobi, verbose, 1); sol optimize(C, obj, ops); if sol.problem ~ 0 warning(求解异常%s, sol.info); end P_buy_v value(P_buy); P_chp_v value(P_chp); Q_chp_v value(Q_chp); P_p2g_v value(P_p2g); fprintf(总运行成本: %.2f\n, value(obj));这段代码里最关键的是三条平衡方程。电平衡里 P2G 耗电被当成额外的“电负荷”由购电和 CHP 供电共同承担热平衡里 H_sto 放热为正相当于热负荷的直接供给方气平衡写成“购气 P2G 产气 气负荷 机组耗气”符号正好和电、热平衡相反容易写反要特别盯住。目标函数最后加的1e-3 * sum(abs(H_sto))是一个很小的正则项目的是让储热罐不要频繁充放数值别超过运行成本的千分之一否则结果会偏离最小成本。这套骨架要扩展成论文里的完整模型动三个地方就行。第一把 T 从 24 改成 96负荷曲线换成 15 分钟采样储热罐的充放功率上限对应的能量步长要乘以 0.25。第二加上网络约束时在每个节点写各自的功率平衡而不是一个母线平衡。第三加二进制变量做机组启停或 P2G 分段运行。扩展顺序按“先设备层再网络层最后整型变量”来每加一层就用一次这个骨架验证可解性。3.3 求解器调用与结果导出YALMIP 参数确认与检查求解状态YALMIP 写模型只是前半段求解器配置不对会浪费大量时间。代码里的sol.problem字段是第一个要检查的0 表示正常1 表示不可行2 表示无界。不可行十有八九是约束冲突或单位不一致无界则常见于目标函数符号写反或变量缺下界。排查时用check(C)看每条约束的残差比人肉读代码快得多。导出结果我用writematrix把电、热、气三个网络的调度结果和一维变量一起落盘方便后面画图和论文图表对照result [(1:T), Pd, Qd, Gd, P_buy_v, P_chp_v, Q_chp_v, P_p2g_v]; writematrix(result, result_dispatch.csv);如果本机装的是 CPLEX 或 MOSEK把sdpsettings(solver, gurobi)里的求解器名换掉即可纯 LP 模型也可以换成linprog。YALMIP 的好处是求解器接口统一模型不用重写。4. 线性化与收敛参数让气网热网约束可解、可调、不翻车4.1 Weymouth 方程怎么进求解器分段线性化才是论文源程序的主流气网管道流量和气压的关系通常用 Weymouth 方程表达管道流量正比于两端气压平方差的开方符号由压差方向决定。这个方程直接写进优化模型是非线性的Gurobi 不会接。常见做法是换成“压差平方”为自变量做分段线性化把平方关系分成若干段用二进制变量选择当前落在哪一段% Weymouth 增量分段线性化片段以压差平方 D2 为自变量 K 8; % 分段数 d2_bound linspace(0, D2max, K1); % 压差平方分段点 seg d2_bound(2:K1) - d2_bound(1:K); % 各段区间宽度 F_end Kij * sqrt(d2_bound); % 各分段端点流量 slope (F_end(2:K1) - F_end(1:K)) ./ seg; % 各段斜率 xk sdpvar(1, K); % 各段内偏移量 zk binvar(1, K); % 该段是否被激活 D2 sdpvar(1, 1); % 当前支路压差平方 C [C, D2 d2_bound(1) sum(xk)]; % 压差平方由各段偏移叠加 C [C, sum(zk) 1]; % 只能落在其中一段 C [C, 0 xk seg .* zk]; % 偏移量被所在段夹住 F_ij_lin F_end(1) sum(slope .* xk); % 线性化后的管道流量 C [C, F_ij F_ij_lin]; % 用 F_ij_lin 替换原 F_ij这是一个典型的增量线性化写法比简单的“折线连接”更可靠因为它用二进制变量保证了解点一定落在同一条分段内不会出现跨段外推。代价是每段都会引入一个 0-1 变量支路数量多时整数变量数量会明显增加。分段数 K 选 4 时模型跑得快但压差误差可能到 10%选 10 精度好但求解时间翻倍。我一般先 8 段跑通确认结果合理再决定要不要加段数。如果发现所有节点压力都顶在边界值上多半是 Kij 的单位和流量单位对不上要先做支路流量与压差平方的散点拟合确认系数在同一数量级再调整段数。4.2 热网温度怎么处理外层迭代回水温度内层求解调度热网的水力方程和温度降落方程如果全部写进 MIP模型规模会爆炸。论文源程序里常见做法是准稳态顺序迭代调度模型求解时只保留热功率平衡回水温度在外面用一个循环迭代更新供水温度固定回水温度影响热负荷对应的质量流量。代码骨架如下% 热网回水温度与调度模型解耦的迭代框架 T_ret 40 * ones(1, T); % 回水温度初值℃ for iter 1:30 sol optimize(C, obj, ops); % 内层求解调度模型 Q_heat value(Q_chp) value(Q_gb) value(H_sto); m_dot Q_heat ./ (Cp * (T_supply - T_ret)); % 质量流量 T_ret_cal T_supply - 1.03 * Q_heat ./ (Cp * max(m_dot, 1e-3)); % 回水温度 T_ret_new 0.7 * T_ret 0.3 * T_ret_cal; % 阻尼更新防震荡 if max(abs(T_ret_new - T_ret)) 0.05 T_ret T_ret_new; break; end T_ret T_ret_new; end这里max(m_dot, 1e-3)是防止供热功率接近 0 的时段出现除零。系数 1.03 表示管网散热造成的等效热损失程序里经常给一个略大于 1 的常数。阻尼系数 0.7/0.3 的搭配是我试过比较稳的配置太激进容易在 40 和 70 度之间来回跳太保守要迭代 30 次以上。如果外层迭代一直发散先怀疑初值把 T_ret 初值设成论文给的回水温度附近而不是从 0 度开始猜。4.3 求解器收敛参数从 MIPGap 到 TimeLimit 的常规设置综合能源调度模型一旦带上网络约束和二进制变量最容易出问题的不是可解性而是整数爆炸。求解器参数设置不是越高精度越好而是要匹配调试阶段的目标。首次跑通时我通常这样设置ops sdpsettings(solver, gurobi); ops sdpsettings(ops, gurobi.MIPGap, 1e-4, ... % 最优间隙 gurobi.TimeLimit, 600, ...% 单次求解时间上限 gurobi.MIPFocus, 1, ... % 优先找可行解 verbose, 1);MIPGap 设成 1e-4 对论文算例足够96 时段加网络约束时先放宽到 5e-3拿到可行解再收紧。TimeLimit 设 600 秒是为了防止模型卡死在一棵树上。MIPFocus 设为 1 代表优先找可行解适合首次跑通如果已经确定模型没问题、想证明最优性再改成 2。调参时看求解器日志里的 gap 变化如果 gap 从 100% 一路往下掉说明模型没问题只是慢如果 gap 长时间纹丝不动多半是线性化约束写错或 Big-M 取值过大。Big-M 常见的坑是为了“绝对没问题”取 10000结果数值病态各支路流量全往边界上跑。合理做法是把 Big-M 压到支路流量上限的 1.1 倍。5. 避坑电-气-热源程序调试中的五个常见现场5.1 P2G 在谷电时段疯狂耗电目标函数里漏了“气平衡”符号现象结果里 P2G 满发但对应产气量没有体现在气负荷那边或者气负荷没变却多出一大笔购气。原因气平衡方程里 G_p2g 符号写反了或者 P2G 效率被填成大于 1谷电价格低时耗电变成“无本套利”。解决先打印约束残差用 YALMIP 的check(C)看每一行约束的残差再重点查 P2G 的效率和变量单位。给购气变量加上0 G_buy 购气上限避免出现负购气这种假解。5.2 气网节点压力全是边界值Weymouth 松弛参数给得太宽现象所有节点压力顶在上限或下限管道流量却和论文结果相差很大。原因压差平方项被过度松弛可行域被放大压力约束形同虚设。解决缩小压差平方的上限 D2max检查 Kij 单位是否与流量单位一致。用 4.1 的分段线性化之前先画一画支路流量与压差平方的散点图确认 Kij 在同一数量级。模型对压力不敏感说明你松弛的不是精度而是物理约束本身。5.3 热网回水温度迭代发散不要动求解器先改初值现象外层迭代的 T_ret 在 40 和 70 度之间反复横跳内层求解一直正常但没有稳定解。原因回水温度更新补偿系数太激进或者初值远离论文设定的运行工况供热功率接近 0 的时段出现了除零。解决初值改成论文给的回水温度上下 2 度以内调试期把阻尼系数改成对半开比如0.5 * T_ret 0.5 * T_ret_cal再逐步加大比例质量流量加max(m_dot, 1e-3)下限。热网这个部分是最需要耐心调的因为它不是求解器问题是初值问题。5.4 Gurobi 报“二次等式”或“整数变量”错误先查约束里混入了平方项现象求解器报出 Quadratic equality 或 Q 1 之类的错误模型直接不求解。原因约束里把气压平方直接写成等式或者两个二进制变量相乘被写进了约束。解决用 YALMIP 的export(C, obj, ops)把模型导出成文件打开看报错的具体行。平方等式改成 4.1 的分段线性化或二阶锥松弛两个二进制变量相乘用辅助变量替换。这是这个方向最典型的血泪经验多数报错不是算法问题而是模型表达不合规。5.5 同一段程序换个电脑跑不通环境版本差异比模型更常见现象自己在熟悉的机器上跑得好好的拷到另一台机器报 YALMIP undefined function、Gurobi license 错误或内存不足。原因YALMIP 版本和 Gurobi API 版本不匹配MATLAB 路径里没有包含 YALMIP 目录求解器路径写死了。解决在程序开头加环境准备代码addpath(genpath(yalmip_path))用yalmiptest检查当前可用的求解器再跑主程序。把环境检查写进脚本头部而不是指望每台电脑都配好这是省时间最划算的一步。6. 从跑通到改进用论文图表复现和场景扩展给程序上保险6.1 用论文的日调度曲线复现画图对照比看成本数更快程序跑通后第一件事不是看总成本数字而是画堆叠图对照论文里的日调度曲线。把电出力、热出力、P2G 耗电和购电画在同一张图里形状相似但幅值差一点多半是单位折算问题形状完全不一样先看三条平衡方程的符号再看负荷数据有没有对齐。绘图代码不需要复杂figure; plot(1:T, P_buy_v, b, 1:T, P_chp_v, r, 1:T, P_p2g_v, m, LineWidth, 1.5); legend(购电,CHP电出力,P2G耗电); xlabel(时段/h); ylabel(功率/MW);6.2 三类改参数的验证极端负荷、检修停运、碳价灵敏度验证程序能不能用于写论文我一般会做三组改动。第一把电负荷整体乘 1.3看购电和 P2G 是否按预期响应这一步能揭示平衡方程是否写反。第二把 CHP 的可用状态置 0模拟检修停运看系统是否自动转向燃气锅炉和储热放热。第三把碳价从 0 扫到 100画出碳排放量和总成本的灵敏度曲线这是评审最喜欢问的部分。这三组改动都不需要动模型结构只改参数就能看出程序的行为是否符合物理直觉。6.3 一个多年调试形成的习惯我的习惯是拿到程序先花半小时把一个脚本写好把所有参数、论文公式出处、单位换算过程以注释形式写进代码头部然后才跑第一遍。这个脚本不产生任何结果但能保证你在第 20 次改模型时还能知道自己当初为什么把某个效率填成 0.45。综合能源调度的源程序调试绝大多数时间浪费在“记不清参数从哪来”上而不是算法本身。这个习惯帮我躲过了很多返工希望帮到你。本文还有配套的精品资源点击获取
