简介面向冷热电多能互补综合能源系统优化调度场景这套 Matlab 源码及运行结果包提供了完整的仿真实现。代码支持 matlab2014/2019a/2021a 运行采用参数化编程思路清晰且注释详细方便用户按需修改冷热电负荷、设备容量等参数并内置案例数据可直接复现实验。压缩包共 26 个文件以 14 个 .m 脚本为核心辅以 6 个 .txt 说明文档、3 个 .xlsx 数据表格、2 篇对应 pdf 论文参考及 1 个嵌套压缩包整体大小 2.82MB结构紧凑便于课程设计、期末大作业或毕业设计快速上手。资源同时涵盖经济性调度与排放优化两种目标模式参考了舒适度与冷热电气多能互补的综合能源网络等文献可帮助研究者理解微能源网鲁棒优化调度的建模与求解思路。目前已有 473 人学习使用适合电气工程、能源系统及自动化相关专业学生作为算法仿真与论文复现的参考。1. 冷热电多能互补综合能源系统优化调度先把物理过程变成数学方程对以燃气轮机为核心的园区型综合能源系统电、热、冷三条母线通过余热回收、电制冷机和吸收式制冷机耦合在一起任何一个负荷波动都会沿这条能量链传导到其他母线上。冷热电多能互补综合能源系统优化调度要做的是在逐时负荷、分时电价和光伏预测出力都已知的前提下算出未来 24 小时每一台设备的最优出力计划让运行成本、碳排放等指标尽量小。这个问题的本质接近一个带大量等式不等式约束的混合整数优化问题变量是各设备出力目标函数是经济性约束是设备爬坡范围、储能状态与母线平衡。下面给出的 MATLAB 建模思路不依赖某个特定压缩包的数据结构你可以把它改写到自己手头的园区数据上几个关键点分别是能源枢纽建模、目标函数线性化、求解器选型、结果自校验。2. 冷热电多能互补系统建模把设备出力写进三条母线方程2.1 电、热、冷三条母线各只有一个平衡方程调度模型里最常见的抽象方式是“能量母线”每一种能量形式对应一条母线母线之间由设备连接。电母线上连着燃气轮机、光伏、电网购电、电制冷机、电储能热母线上连着燃气轮机的余热、燃气锅炉、蓄热槽同时向吸收式制冷机供热水冷母线上连着电制冷机、吸收式制冷机和蓄冷槽。这样做的好处是设备内部的热力学细节被压缩成一个输入输出的代数关系平衡方程的数量和能量形式一一对应。典型的三条母线平衡方程可以写成下面这组形式t 表示调度时段所有功率单位取 kW电母线 P_gt(t) P_pv(t) P_buy(t) P_edis(t) L_e(t) P_ec(t) P_ech(t) 热母线 Q_gt(t) Q_gb(t) Q_hs_dis(t) L_h(t) Q_ac_in(t) Q_hs_ch(t) 冷母线 COP_ec * P_ec(t) COP_ac * Q_ac_in(t) L_c(t) Q_cs_ch(t) - Q_cs_dis(t)其中 P_gt 是燃气轮机电出力Q_gt 是它回收的余热P_pv 是光伏出力P_buy 是向电网购电P_ec 是输入电制冷机的电功率Q_ac_in 是输入吸收式制冷机的热功率P_ech/P_edis 是电储能充放电功率Q_hs_ch/Q_hs_dis 是蓄热槽充放热Q_cs_ch/Q_cs_dis 是蓄冷槽充放冷。COP_ec 与 COP_ac 是两类制冷机的能效比。把冷负荷与热负荷分开写是因为两者温度品位不同在工程上一般也不混用而 AC 的存在让热转冷成为可能这是“冷热电互补”里最关键的一个耦合点。2.2 设备出力模型与常用参数表母线方程解决的是“能量从哪来、到哪去”设备模型解决的是“输入输出怎么换算”。调度用的设备模型不需要精确到燃烧室温度或压缩机压比线性或分段线性关系就足够。下面这些关系在 MATLAB 代码里可以直接作为等式或不等式约束出现。设备输入→输出关系出力范围典型值备注燃气轮机 CHPQ_gt (η_h/η_e) * P_gtP_gt ∈ [0, 60] kW热电比固定时是线性关系燃气锅炉 GBQ_gb η_gb * F_gb * LHVQ_gb ∈ [0, 100] kWF_gb 为燃气流量电制冷机 ECQ_c COP_ec * P_ecP_ec ∈ [0, 40] kWCOP_ec 取 3.5吸收式制冷机 ACQ_c COP_ac * Q_ac_inQ_ac_in ∈ [0, 80] kWCOP_ac 取 1.2电储能 ESSSOC 递推方程容量 50 kWh充放效率约 0.95蓄热槽 TSS同上容量 80 kWh散热损失率 2%蓄冷槽 CSS同上容量 100 kWh不适合频繁启停一个常见的坑是燃气轮机的热电比。取 η_e0.35、η_h0.45 时Q_gt (0.45/0.35) P_gt意味着发 10 kW 电的同时回收约 12.86 kW 热。很多入门代码把电和热当成两个独立变量随意约束结果求解器给出的出力组合在物理上根本不存在。处理办法很简单让 Q_gt 用 P_gt 的等式表达不要另设成独立变量或者采用线性化可行域来近似变热电比运行区。2.3 储能 SOC 递推方程与日循环约束电、热、冷三类储能虽然介质不同递推方程完全一样E(t1) (1 - σ) * E(t) η_ch * P_ch(t) - P_dis(t) / η_disσ 是自损率η_ch 是充电效率η_dis 是放电效率。对冷热电系统来说这个方程要写三份分别对应电储能 SOC、蓄热槽储热量和蓄冷槽储冷量。在 MATLAB 里通常写成向量化约束% 电储能SOC方程T为24E_e为SOC变量P_ech/P_edis为充放电功率 sigma_e 0.02; eta_ech 0.95; eta_edis 0.95; prob.Constraints.soc_e E_e(2:T) ... (1 - sigma_e) * E_e(1:T-1) eta_ech * P_ech(1:T-1) ... - P_edis(1:T-1) / eta_edis; % 一周循环初末SOC相等避免调度把储能“用完” prob.Constraints.init_end_e E_e(1) E_e(T);这里的写法用到了 MATLAB 的 problem-based 优化接口E_e、P_ech、P_edis 都是 T 维优化变量。E_e(2:T) 和 E_e(1:T-1) 构成逐时递推避免写 24 条 for 循环求解性能会好很多。init_end_e 这条约束很关键否则优化器会倾向于在最后一个时段把储能放空得到一个不可持续的计划。给储能变量设置上下界时E_e 的区间是 [0, 50]充放电各设 [0, 20]这些数字来自设备铭牌而不是随意取的罚函数参数。3. 目标函数与约束条件把调度问题做成标准优化形态3.1 经济目标函数的线性化写法运行成本一般由购电费、燃气费、设备维护费三部分组成。分时电价是逐时变化的购电费写成 price(t) * P_buy(t) 的累加燃气费由燃气轮机和燃气锅炉消耗的天然气共同决定按流量折成费用维护费通常按设备出力乘以一个很小的系数比如 0.02 元/kWh。于是目标函数写成min Σ [ C_buy(t) * P_buy(t) C_gas * ( P_gt(t)/η_e Q_gb(t)/η_gb ) / LHV C_m * ( P_gt(t) Q_gb(t) P_ec(t) ) ]C_buy 是分时电价C_gas 是天然气单价LHV 是天然气低位热值取 9.7 kWh/m³C_m 是维护费系数。所有项都是决策变量的线性组合这决定了整个问题可以用线性规划或混合整数线性规划求解。如果目标是碳排放最小只要把购电和燃气的碳排放系数加进去替换价格项即可问题形态不变。% 目标函数购电费 燃气费 维护费 C_gas 3.5; % 天然气单价元/m³ LHV 9.7; % 低位热值kWh/m³ C_m 0.02; % 单位运维成本元/kWh prob.Objective sum(price_e .* P_buy) ... C_gas / LHV * (sum(P_gt)/eta_e sum(Q_gb)/eta_gb) ... C_m * (sum(P_gt) sum(Q_gb) sum(P_ec));需要留意单位一致性P_gt、Q_gb 的单位是 kW时间粒度为 1 小时所以 sum(P_gt) 本身就表示全天发电量 kWh电价单位是元/kWh两者相乘直接得到元。若调度间隔改成 15 分钟sum 前要乘以 0.25。3.2 约束清单范围、爬坡、储能与耦合约束条件按作用对象可以分成五层。第一层是设备出力上下限也就是给每个变量设置 LowerBound 和 UpperBound。第二层是 CHP 热电耦合和制冷机 COP 关系用等式表达。第三层是爬坡约束燃气轮机每分钟只能增加或减少一定出力折算到小时级别就是 |P_gt(t) - P_gt(t-1)| ≤ R_gt。爬坡约束需要用两个不等式实现直接写成 P_gt(t) - P_gt(t-1) ≤ R_gt 和 P_gt(t-1) - P_gt(t) ≤ R_gt。第四层是储能 SOC 范围与充放电功率上限。第五层是母线平衡等式。约束项数学表达典型参数CHP 热电耦合Q_gt (η_h/η_e) P_gt0.45/0.35燃气轮机爬坡|P_gt(t)-P_gt(t-1)| ≤ 2020 kW/h购电上限P_buy ≤ 200并网容量 200 kW储能 SOC 上下限0.2C ≤ E ≤ 0.9C防止过充过放蓄冷槽充放Q_cs_ch, Q_cs_dis ≤ 30受限于制冷机余量储能 SOC 上下限一般不放成全 0 到全容量建议留 10% 到 20% 的安全裕量尤其是蓄冷槽。蓄冷槽在电费低的夜间蓄冷、白天融冰供冷如果允许 SOC 到 0第二天遇到极端热负荷时没有应对余地加入上下限之后求解结果会自然把 SOC 维持在合理区间。3.3 求解器选择linprog、intlinprog 还是 Yalmip调度模型到底用什么工具取决于变量类型。全部是连续变量用 MATLAB 自带的 linprog 就够出现 0-1 启停变量或储能同时充放互斥约束要用 intlinprog模型里存在非线性项比如变热电比或罚函数二次项才值得引入 Yalmip 配合商业求解器。在实际代码中我一般先用连续模型跑通确认平衡方程没有写错再加 0-1 变量处理启停和储能互斥。原因是连续模型的求解和调试都要快得多。对 24 时段 × 10 个设备的规模优化工具箱里的 linprog 和 intlinprog 都能在几秒内解完完全不需要一开始就上 Gurobi。只有在做 8760 小时全年仿真、单次求解时间超过分钟级时才考虑换求解器和紧缩变量数量。求解路径建模工作量适用规模依赖linprog/intlinprog中千级变量以内MATLAB 优化工具箱Yalmip Gurobi低大规模 MILPYalmip 工具箱 第三方求解器自编梯度/启发式高非线性强需要写大量回调MATLAB R2023b 及之后版本里优化工具箱可以直接在附加功能管理中安装装好后用 which intlinprog 验证路径。如果只做线性目标加线性约束problem-based 接口比 solver-based 接口更适合反复改约束因为它把变量名和方程写法都保留成了可读的表达式。4. MATLAB 代码落地从数据准备到结果曲线绘制4.1 第一段代码把 24 小时负荷、分时电价和光伏出力组织成列向量调度代码的第一步是数据区。电、热、冷负荷曲线、分时电价、光伏预测出力全部整理成 24×1 的列向量顺序和优化变量的 t 索引对齐。负荷单位统一用 kW电价用元/kWh。以下是一组典型园区数据组织方式T 24; % 电负荷工作日上午和傍晚两个峰 L_e [62 58 55 54 56 62 78 95 105 110 116 118 ... 112 108 104 116 118 124 120 108 92 80 70 65]; % 热负荷夜间采暖白天部分生活热水 L_h [80 84 86 82 78 72 60 45 38 32 28 26 ... 25 24 26 28 30 36 48 62 74 82 84 82]; % 冷负荷中午前后最高夜间接近零 L_c [0 0 0 0 0 8 30 62 85 96 102 92 ... 84 92 96 90 78 52 30 12 0 0 0 0]; % 光伏预测出力中午峰值 P_pv_fc [0 0 0 0 0 2 12 25 35 46 52 55 ... 54 48 38 26 14 5 0 0 0 0 0 0]; % 分时电价峰 1.2 / 平 0.8 / 谷 0.4 price_e 0.4 * ones(T,1); price_e(8:11) 1.2; price_e(18:21) 1.2; price_e(12:17) 0.8;这段数据只是示例实际使用时要替换成自己系统的预测值。注意热负荷和冷负荷可能存在交叠时段这在夏季尤其明显代码里 L_h 与 L_c 是独立预测的最后通过吸收式制冷机的热输入在平衡方程中产生耦合。4.2 第二段代码用 problem-based 接口构造并求解调度问题下面给出一个可直接运行的连续变量调度模型。把决策变量、约束、目标函数依次写入 prob然后调用 solve。代码如下prob optimproblem(ObjectiveSense, minimize); % 决策变量燃气轮机电出力、锅炉热出力、两类制冷机输入、购电 P_gt optimvar(P_gt, T, 1, LowerBound, 0, UpperBound, 60); Q_gb optimvar(Q_gb, T, 1, LowerBound, 0, UpperBound, 100); P_ec optimvar(P_ec, T, 1, LowerBound, 0, UpperBound, 40); Q_ac optimvar(Q_ac, T, 1, LowerBound, 0, UpperBound, 80); P_buy optimvar(P_buy, T, 1, LowerBound, 0, UpperBound, 200); % 电储能变量 E_e optimvar(E_e, T, 1, LowerBound, 0, UpperBound, 50); P_ech optimvar(P_ech, T, 1, LowerBound, 0, UpperBound, 20); P_edis optimvar(P_edis, T, 1, LowerBound, 0, UpperBound, 20); % 吸收式制冷机驱动的热来自CHP余热因此先用等式表达CHP热电耦合 Q_gt (eta_h / eta_e) * P_gt; % 三条母线平衡 prob.Constraints.ele P_gt P_pv_fc P_buy P_edis ... L_e P_ec P_ech; prob.Constraints.heat Q_gt Q_gb L_h Q_ac; prob.Constraints.cool COP_ec * P_ec COP_ac * Q_ac L_c; % 电储能SOC递推 prob.Constraints.soc E_e(2:T) ... (1 - sigma_e) * E_e(1:T-1) eta_ech * P_ech(1:T-1) ... - P_edis(1:T-1) / eta_edis; prob.Constraints.soc_init E_e(1) E_e(T); % 目标函数购电费燃气费维护费 prob.Objective sum(price_e .* P_buy) ... C_gas / LHV * (sum(P_gt)/eta_e sum(Q_gb)/eta_gb) ... C_m * (sum(P_gt) sum(Q_gb) sum(P_ec)); % 求解 [sol, fval] solve(prob);这段代码里最有必要说明的是 Q_ac 的语义它表示输入吸收式制冷机的热功率最终在冷母线方程里乘以 COP_ac 变成冷量。因为制热和制冷有可能同时发生热母线和冷母线方程必须分开写否则会出现一部分热被用来制冷、制冷量又超过冷负荷的冗余最优解。如果不希望出现这种情况可以再加一条不等式 Q_ac ≤ L_h Q_gt Q_gb - Q_hs_out但一般平衡等式已经足够。求解器默认会从变量类型推断使用 linprog。所有变量都是连续变量所以调用的是线性规划算法solve 返回的 sol 是一个结构体sol.P_gt、sol.P_buy 分别对应各设备的逐时出力。fval 是全天总运行成本。4.3 第三段代码提取结果、逐时堆叠图与数据表求解之后把优化变量转成普通数组再画图。常见做法是画三张图设备逐时出力堆叠图、各母线电量平衡图、储能 SOC 变化曲线。堆叠图最能看出“多能互补”的实际含义代码如下% 从求解结果中提取调度计划 P_gt_sol value(sol.P_gt); P_buy_sol value(sol.P_buy); P_ec_sol value(sol.P_ec); E_e_sol value(sol.E_e); % 电母线出力堆叠图 t (1:T); figure(Color,w); bar(t, [P_buy_sol, P_gt_sol, P_pv_fc], stacked); hold on; plot(t, L_e, k-o, LineWidth, 1.5); xlabel(时刻/h); ylabel(功率/kW); legend(购电,燃气轮机,光伏,电负荷,Location,best); grid on; % SOC曲线确认初末相等 figure(Color,w); plot(t, E_e_sol, r-o, LineWidth, 1.5); xlabel(时刻/h); ylabel(SOC/kWh); title(电储能SOC日循环曲线);bar 的 stacked 参数把购电、燃气轮机、光伏三段堆起来叠加出来的上边界应该和电负荷曲线在多数时段重合两条线之差就是储能充放造成的净功率。SOC 曲线首尾相等说明日循环约束起作用了如果不等回去检查 init_end 约束是否写成两点比较而不是完整时段比较。4.4 解不出来或结果明显异常时的排查顺序这一步是实际改动代码时最容易花时间的地方。排查顺序按三个方向走先看数据再看约束最后看求解器状态。数据方面L_e、L_h、L_c 不要有 NaN 或负值光伏出力不能超过设备上限电价向量长度必须等于 T。约束方面用 show(prob.Constraints.ele) 检查每个等式的变量顺序重点看有没有把 P_ech 和 P_edis 的符号写反。求解器状态方面exitflag 为 1 是收敛为 -2 是无可行解这时应放宽某一条灵敏约束比如购电上限从 200 提到 250确认是不是容量不够导致无解。通常压缩包里附带“运行结果”目录时里面会有对应的 .mat 文件保存变量快照方便你把自己的数据和示例结果做差分对比。5. 运行结果的验证技巧从曲线合法性到参数灵敏度5.1 平衡残差与 SOC 初末一致性校验调度结果拿到手不要先看目标函数值先做两层机械校验。第一层是母线平衡残差把求解出的出力代回三条平衡方程残差应小于 1e-6。MATLAB 里直接算res_ele P_gt_sol P_pv_fc value(sol.P_buy) value(sol.P_edis) ... - L_e - P_ec_sol - value(sol.P_ech); disp([电母线最大残差, num2str(max(abs(res_ele)))]); res_soc E_e_sol(2:end) - 0.98*E_e_sol(1:end-1) ... - 0.95*value(sol.P_ech(1:end-1)) ... value(sol.P_edis(1:end-1))/0.95; disp([SOC递推最大残差, num2str(max(abs(res_soc)))]);如果残差在 1e-3 量级通常是单位混了把 kW 与 kWh 混用如果残差在几百多半是等式约束没写进 prob检查一下变量名是否被覆盖。第二层校验是看储能 SOC 是否始终处于上下限之内、初末是否相等。一个常见错误是 E_e 初始值被固定成某个常数而不是变量这会让优化器牺牲第一个时段的可行性去迁就初始状态。正确做法是把 E_e(1) 也作为优化变量用 init_end 约束把首尾绑在一起。5.2 灵敏度实验分时电价与光伏出力对调度决策的影响调度代码跑通后最值得做的一次实验是电价灵敏度。把峰时电价从 1.2 提到 1.5观察两个现象燃气轮机在峰时段的出力是否上升电储能在低谷充电、峰时放电的功率是否变高。如果峰时电价提高了储能出力曲线却没有变化大概率是 SOC 上下限或充放功率上限设得过于宽松导致目标函数中维护费占主导这时把维护费系数调低一个数量级再看曲线。同样的方法也可以验证光伏出力把 P_pv_fc 整体乘以 0.5冷热电各自的购电量应同时发生变化尤其是中午时段。这种灵敏度分析在论文和工程报告中都是必要的。它证明了模型不是一条写死的规则表而是真正在响应价格信号。压缩包里的“运行结果”如果只有一张总成本对比图通常不够我把调度计划按峰谷平三个时段分解成表格分别列出购电、燃气轮机发电、电制冷机耗电三列这样才能看出互补行为发生的具体时段。5.3 把结果导出为规范化格式最后一步是把调度结果和控制指令对齐。用 writetable 把逐时调度计划导出成 CSV列名建议是 hour、P_gt、P_gb、P_ec、P_ac、P_buy、P_sell、ESS_SOC、TSS_SOC、CSS_SOC。导出前把所有数值保留两位小数避免给出虚假精度。对连续变量模型导出的功率值直接可用如果是带 0-1 启停变量的版本还要单独导出一列 on_off 状态供后续执行层读取。这样跑出来的结果既能在 MATLAB 里画图复现也能进报表系统作日调度复盘。本文还有配套的精品资源点击获取