简介一份面向电力与能源系统、自动化及计算机等专业学生的Matlab实现资料聚焦冷热电多能互补综合能源系统优化调度问题适用于课程设计、期末大作业和毕业设计等场景。代码采用参数化编程思路清晰、注释明细参数可方便修改便于读者学习核心算法并快速迁移到自己的算例中。压缩包共26个文件包括14个m源程序、6个txt说明文档、3个xlsx数据文件、2个pdf参考文献和1个zip附加包整体大小2.82MB目录结构清晰既有可运行主程序也有配套数据与理论文献。其中pdf包含《冷热电气多能互补的微能源网鲁棒优化调度》等参考资料能帮助理解多能互补建模与优化调度原理运行结果和案例数据可支持直接复现实验适合需要完整实践项目的读者。资源已有473人学习是一份兼顾代码、数据与文档的实用资料。1. 冷热电多能互补优化调度为什么是优化问题而不是控制问题拿到“冷热电多能互补综合能源系统优化调度”这个标题多数人第一反应是去翻潮流计算或设备选型但核心其实落在“调度”二字上。调度在数学上是一个典型的多变量、多约束优化问题给定电负荷、热负荷、冷负荷曲线决定每一台机组在每一个时段发多少电、制多少热、供多少冷让总运行成本最低或碳排放最小。它和控制问题的区别在于控制关心反馈闭环与动态响应而调度关心的是未来24小时甚至更长时间尺度上的资源分配属于开环规划需要预先求解一组最优设定值。这个方向在学术和工程里统称 IESIntegrated Energy System或 CCHPCombined Cooling, Heating and Power系统优化目前主流解法是用 MATLAB 建模配合 YALMIP 工具箱调用 CPLEX、Gurobi 或自带求解器来做混合整数线性规划MILP。标题里的“附matlab代码运行结果”说明是面向复现的本文就按“模型怎么建→代码怎么写→结果怎么读→坑怎么避”这条线完整展开适合正在做毕业设计、写小论文或刚接触园区综合能源仿真的工程师。2. 设备建模与能量流关系从物理拓扑到数学方程2.1 典型系统拓扑母线式结构是最省事的建模方式工程上做综合能源系统建模最常见的拓扑是母线式结构就是把电、热、冷分别抽象成三条独立的能量母线所有设备像插排一样挂在对应母线上。电母线连接电网购电、燃气轮机发电、电制冷机耗电、电锅炉耗电热母线连接余热锅炉、燃气锅炉、热交换器冷母线连接吸收式制冷机和电制冷机。这种拓扑的好处是能量平衡方程可以逐母线独立书写不需要构建复杂的节点网络MATLAB 里用一个矩阵就能表达全时段的功率分配。母线式结构的建模精度对于调度问题已经足够。它不关心管道压降、线路损耗和温度分布这些属于热力学稳态或动态仿真范畴调度模型只需要保证每个时段内流入母线的功率等于流出功率加储能充放。所以别一上来就写微分方程先把这个静态平衡关系理清调度模型就完成了一半。2.2 设备静态模型效率、爬坡和最小启停约束设备模型可以用一个统一框架描述输入功率、输出功率、效率或热电比、运行上下限、爬坡速率、最小启停时间。以燃气轮机为例输出电功率 (P_{gt}) 与消耗天然气功率 (F_{gt}) 之间的关系通常写成线性或分段线性形式[ F_{gt}(t) \frac{P_{gt}(t)}{\eta_{gt}^{e}} ]其中 (\eta_{gt}^{e}) 是发电效率。余热锅炉回收的烟气余热与发电功率存在固定比例关系也就是热电比 (\alpha_{hr})回收热功率为[ H_{hr}(t) P_{gt}(t) \cdot \alpha_{hr} ]注意这里的两个典型误差一是余热回收量在部分负荷率低时并非严格线性工程上可在模型中引入最小负荷率开关变量来规避二是燃气轮机启停状态必须用二进制变量表示否则求解结果会出现“设备只在连续几个时段各发 0.5kW 来凑平衡”的荒谬结果。2.3 储能模型能量时序耦合是调度的灵魂储能设备蓄电、蓄热、蓄冷给模型引入了跨时段耦合约束这是优化调度区别于单纯经济分配的关键。储能通常写成如下离散形式[ E_{st}(t1) E_{st}(t) \cdot (1 - \sigma_{st}) \left( P_{ch}(t) \cdot \eta_{ch} - \frac{P_{dis}(t)}{\eta_{dis}} \right) \Delta t ]这里 (\sigma_{st}) 是自放率(\eta_{ch}) 和 (\eta_{dis}) 分别是充放能效率(P_{ch})、(P_{dis}) 为充放功率(\Delta t) 为调度步长。储能模型容易踩的坑是充放功率同时为正这在物理上不可能但在线性规划里却可能被求解器当作“循环充放”来套利所以必须加上互补约束[ P_{ch}(t) \le M \cdot u_{st}(t), \quad P_{dis}(t) \le M \cdot (1 - u_{st}(t)) ]其中 (u_{st}(t)) 是二进制变量(M) 取充放功率上限的 1.5 到 2 倍即可过大的 M 值会破坏求解器数值稳定性。另外调度周期首末储能状态通常要求相等形成周期性调度否则模型会把储能初始能量全部耗尽以降低成本。3. 目标函数与约束体系从物理层到优化层的完整落地方案3.1 目标函数怎么选运行成本最小是默认配置系统优化调度的目标函数最常见的是日运行成本最小包括购电费用、天然气费用、设备启停费用和弃能惩罚。表达式结构如下[ \min \quad \sum_{t1}^{T} \left( c_{grid}^{buy}(t) P_{grid}(t) - c_{grid}^{sell}(t) P_{sell}(t) c_{gas} F_{total}(t) c_{su}(t) c_{sd}(t) c_{curtail} P_{curt}(t) \right) \Delta t ](c_{grid}^{buy}(t))分时电价峰平谷三段或实时电价曲线(c_{grid}^{sell}(t))上网电价通常远低于购电价体现“买贵卖贱”(c_{gas})天然气折算单价单位元/kWh(c_{su})、(c_{sd})启停成本用二进制状态变量差分计算(c_{curtail})弃风弃光惩罚系数需要设一个较高值防止模型为了省钱而大量丢弃可再生出力如果论文侧重低碳目标函数可以改为碳排放最小或直接做多目标加权。但实际求解时多目标线性加权会产生量纲问题建议按“碳税折算进成本”的单一目标处理既简洁又符合当前碳交易背景。3.2 平衡约束三条母线各自守恒每个调度时段 t三条母线的功率平衡方程为以电流源型为例电平衡[ P_{grid}(t) P_{gt}(t) P_{pv}(t) P_{dis}^{e}(t) P_{load}^{e}(t) P_{ec}(t) P_{eb}(t) P_{ch}^{e}(t) P_{sell}(t) ]热平衡[ H_{hr}(t) H_{gb}(t) H_{dis}^{h}(t) H_{load}(t) H_{ch}^{h}(t) ]冷平衡[ C_{ac}(t) C_{ec}(t) C_{dis}^{c}(t) C_{load}(t) C_{ch}^{c}(t) ]其中下标 ec 表示电制冷机、eb 表示电锅炉、gb 表示燃气锅炉、ac 表示吸收式制冷机。方程左侧是能量来源右侧是去向和负荷写平衡方程时一定要把同一设备的输入和输出放在了不同等式中例如电制冷机的耗电量在电平衡右侧而其制冷量在冷平衡左侧。3.3 不等式约束设备出力边界和爬坡约束设备出力上下限约束、爬坡约束统一写为[ u_i(t) \cdot P_i^{min} \le P_i(t) \le u_i(t) \cdot P_i^{max} ][ -\Delta P_i^{down} \le P_i(t) - P_i(t-1) \le \Delta P_i^{up} ]爬坡约束只在设备连续运行时才生效因为启停瞬间的功率跳变不归爬坡管。实操中如果发现模型无解先检查是不是把启动状态下的第一时段也套了爬坡约束这是新手最常见的死锁原因之一。正确的写法是在爬坡约束右侧加上一个与启停变量相关的松弛项。4. MATLAB 代码实现YALMIPCPLEX 求解 MILP 的完整代码级教程4.1 代码框架与参数初始化先搭骨骼再填肌肉下面给出一个可直接运行的框架代码核心思路清晰且便于二次开发覆盖场景为典型夏季典型日 24 时段调度。请按步骤在 MATLAB 中依次执行。设备参数定义如下代码块%% 冷热电综合能源系统优化调度 - 参数初始化 % 清空环境 clear; clc; close all; % 时间参数 T 24; % 调度时段, 单位h dt 1; % 时间步长, 单位h % 设备容量参数 (单位: kW 或 kW/h) P_gt_max 1000; % 燃气轮机最大发电功率 P_gt_min 200; % 燃气轮机最小发电功率 eta_gt 0.35; % 燃气轮机发电效率 alpha_hr 1.2; % 余热回收热电比 eta_gb 0.9; % 燃气锅炉效率 COP_ec 3.5; % 电制冷机能效比 COP_ac 1.2; % 吸收式制冷机能效比 P_ec_max 800; % 电制冷机最大输入电功率 P_ac_max 600; % 吸收式制冷机最大输入热功率 % 储能参数 E_bat_max 1000; % 蓄电池容量, kWh P_bat_max 200; % 蓄电池最大充放电功率 eta_bat 0.95; % 蓄电池充放电效率 sigma_bat 0.01; % 蓄电池自放率 E0_bat 500; % 初始电量 % 分时电价 (元/kWh) - 峰平谷三段 price_buy [0.35*ones(1,7), 0.65*ones(1,5), 1.15*ones(1,5), ... 0.65*ones(1,4), 1.15*ones(1,3)]; price_sell 0.3 * price_buy; % 上网电价, 按购电价30% % 天然气价格 (元/kWh) price_gas 0.35; % 负荷曲线 (kW) load_e [620,580,540,520,500,510,560,680,820,900,980,1050,... 1100,1080,960,1020,1180,1250,1210,1100,980,850,740,660]; load_h [800,760,720,690,650,620,600,580,550,520,500,490,... 480,470,490,520,560,600,700,760,800,820,810,790]; load_c [400,380,360,350,340,360,420,520,650,780,900,1000,... 1080,1100,1050,980,920,880,840,760,680,600,520,450];这段代码中所有容量参数都按工程常见量级设置燃气轮机选 1000kW 级适合园区级场景。注意负荷曲线三个向量的长度必须与 T 相等且一一对应否则后续赋值会维度报错。用峰平谷三段电价模拟电网分时计费机制在峰时段10:00-15:00 和 18:00-21:00电价为 1.15 元/kWh迫使模型在峰时段少购电而多用气。4.2 变量定义与约束构建写出 YALMIP 的标准姿势YALMIP 的变量定义分两类连续变量用 sdpvar二进制变量用 binvar。接线清楚后逐条建立约束%% 变量定义 % 连续变量 P_gt sdpvar(1, T, full); % 燃气轮机发电功率 F_gt sdpvar(1, T, full); % 燃机耗气功率 H_hr sdpvar(1, T, full); % 余热回收功率 H_gb sdpvar(1, T, full); % 燃气锅炉供热功率 F_gb sdpvar(1, T, full); % 燃气锅炉耗气功率 P_ec sdpvar(1, T, full); % 电制冷机输入功率 C_ec sdpvar(1, T, full); % 电制冷机制冷功率 C_ac sdpvar(1, T, full); % 吸收式制冷机制冷功率 H_ac sdpvar(1, T, full); % 吸收式制冷机输入热功率 P_grid sdpvar(1, T, full); % 电网购电功率 P_sell sdpvar(1, T, full); % 上网售电功率 P_bat_c sdpvar(1, T, full); % 蓄电池充电功率 P_bat_d sdpvar(1, T, full); % 蓄电池放电功率 E_bat sdpvar(1, T1, full); % 蓄电池电量 % 二进制变量 u_gt binvar(1, T, full); % 燃机启停状态 u_bat binvar(1, T, full); % 蓄电池充放电状态变量命名采用“前缀_设备_后缀”的风格P 代表功率H 代表热功率C 代表冷功率E 代表能量。binvar 生成的 0-1 变量在求解时用于表示设备启停或充放互斥是 MILP 能够成立的基础。P_bat_c 和 P_bat_d 同时存在但不可同时为正这是后面互斥约束的意义所在。约束构建部分的代码%% 约束条件 C []; % 1. 电平衡约束 C [C, P_grid P_gt P_bat_d - P_sell load_e P_ec P_bat_c]; % 2. 热平衡约束 C [C, H_hr H_gb load_h H_ac]; % 3. 冷平衡约束 C [C, C_ec C_ac load_c]; % 4. 燃气轮机模型约束 C [C, F_gt P_gt / eta_gt]; % 耗气量与发电量关系 C [C, H_hr P_gt * alpha_hr]; % 余热回收量 C [C, P_gt_min * u_gt P_gt P_gt_max * u_gt]; % 出力上下限和启停 % 5. 电制冷机 C [C, C_ec COP_ec * P_ec]; % 制冷量与耗电量关系 C [C, 0 P_ec P_ec_max]; % 输入功率限值 % 6. 吸收式制冷机 C [C, C_ac COP_ac * H_ac]; % 制冷量与输入热功率 C [C, 0 H_ac P_ac_max]; % 7. 燃气锅炉 C [C, F_gb H_gb / eta_gb]; % 耗气量 C [C, 0 H_gb 1500]; % 热出力限值 % 8. 蓄电池储能约束 C [C, E_bat(1) E0_bat]; % 初始电量 C [C, E_bat(2:end) (1 - sigma_bat) * E_bat(1:end-1) ... (P_bat_c * eta_bat - P_bat_d / eta_bat) * dt]; C [C, 0 E_bat E_bat_max]; % 电量上下限 C [C, 0 P_bat_c P_bat_max * u_bat]; % 充电功率上限 C [C, 0 P_bat_d P_bat_max * (1 - u_bat)]; % 放电功率上限 C [C, E_bat(end) E0_bat]; % 周期末电量回归初值 % 9. 购售电互斥 C [C, 0 P_grid 2000 * (1 - u_sell), 0 P_sell 2000 * u_sell]; % 10. 电网交互上下限 C [C, 0 P_grid 2000, 0 P_sell 2000];这里每一条约束都在做同一件事把物理设备的输入-输出关系翻译成线性方程。注意约束 8 中的储能递推式使用了向量整体运算避免了 for 循环这在 YALMIP 里可以显著提升建模速度。蓄电池 SOC 约束使用了 (T1) 个变量第一个是初值最后一个是终值通过首末相等形成日循环调度。4.3 求解器调用与结果输出别让 binvar 毁了你的 CPLEXYALMIP 建模完成后调用外部求解器求解代码如下%% 目标函数 objective sum(price_buy .* P_grid) - sum(price_sell .* P_sell) ... price_gas * sum(F_gt F_gb) ... 0.1 * sum(u_gt ~ 0); % 燃机启动成本近似 %% 求解 ops sdpsettings(solver, cplex, verbose, 2, ... showprogress, 1, mosek.MSK_IPAR_NUM_THREADS, 4); % 如无CPLEX可改用gurobi或默认求解器 ops.solver cplex; % 求解 result optimize(C, objective, ops); if result.problem 0 disp(求解成功); else disp([求解失败: result.info]); endsolve完成后目标函数值和变量值直接从 sdpvar 对象中提取。注意这段代码的u_sell变量在 4.2 节中还没有定义需要在变量定义区补充u_sell binvar(1, T, full);。启停成本用了一个近似写法它统计了非零时段数虽然没有严格计算启动次数但对于日调度其精度已足够。提取并绘图的结果展示代码%% 结果提取与绘图 P_gt_opt value(P_gt); P_grid_opt value(P_grid); P_bat_c_opt value(P_bat_c); P_bat_d_opt value(P_bat_d); E_bat_opt value(E_bat); figure(Position, [100, 100, 1200, 600]); subplot(2,2,1); bar(1:T, [P_gt_opt; P_grid_opt; P_bat_d_opt - P_bat_c_opt], stacked); legend(燃气轮机, 电网购电, 蓄电池净放电); xlabel(时段/h); ylabel(功率/kW); title(电功率平衡组成); grid on; subplot(2,2,2); plot(1:T1, E_bat_opt, b-o, LineWidth, 1.5); xlabel(时段/h); ylabel(电量/kWh); title(蓄电池SOC曲线); grid on; subplot(2,2,3); plot(1:T, P_gt_opt, r, LineWidth, 2); hold on; plot(1:T, value(H_hr), g, LineWidth, 2); legend(发电功率, 余热回收功率); xlabel(时段/h); ylabel(功率/kW); title(燃气轮机运行状态); grid on; subplot(2,2,4); C_ec_opt value(C_ec); C_ac_opt value(C_ac); bar(1:T, [C_ec_opt; C_ac_opt], stacked); legend(电制冷, 吸收式制冷); xlabel(时段/h); ylabel(冷功率/kW); title(冷负荷供应组成); grid on;value()是 YALMIP 的核心提取函数把 sdpvar 变量在最优解处的数值取出。绘图部分的 stacked bar 能直观看到每个时段各设备承担的负荷比例一张图就能判断模型的调度行为是否符合预期。运行成功后应保存为result.mat方便后续对比算法改进前后的效果。5. 运行结果案例分析用数据驱动调度策略修正一套优化模型的结果在哪里它给出了一张时间-功率分配表更关键的是能从中总结出设备协同规律。用上述参数得到的典型结果中燃气轮机在峰电时段10:00-15:00满发余热优先供给吸收式制冷机而不是热负荷因为冷负荷同期也处于峰值吸收式制冷刚好能替代高耗电的电制冷。谷电时段燃气轮机停机电网购电直接供电蓄电池在谷时段充电、峰时段放电。燃气锅炉只在热负荷大于余热回收量的时段启动且大部分时间处于较低负载率。对比一个值得注意的边界情况如果冷负荷全天较低吸收式制冷机的 COP 低于电制冷机模型会宁可开电制冷也不用吸收式即使有余热余量也选择排空。这说明模型决策依据是“效率×价格”的综合比较而非“有热就用热”的固定逻辑。这也是综合能源调度和传统单能系统运行的本质差别——不同能源品种之间的替换关系必须放到目标函数里竞争而不是人为指定。6. 进阶排错与性能提升从能跑到跑得快的实用技巧6.1 模型无解排查三板斧调电价参数模型无解是 MILP 建模最常见的灾难。排查顺序依次为检查分时电价是否出现倒挂、检查负荷向量长度与 T 是否匹配、把储能等式右端项中 (P_{ch} \cdot \eta_{ch} - \frac{P_{dis}}{\eta_{dis}}) 拆成两个独立约束。前两条多半是数据问题第三条是约束缩放问题。如果仍然无解把蓄电池去掉再跑一次若此时有解则问题锁定在储能约束上。6.2 加快求解速度的三个手段对称性破缺和大 M 收紧成熟度较高的模型求解耗时通常在毫秒到秒级。如果遇到大系统数百台设备、数千个时段可以考虑三个实用手段。第一给同类设备加排序约束例如“第一台燃机出力不小于第二台”消除 MILP 中的对称解第二收紧大 M 值储能互斥约束里的 M 取真实功率上限而不是任意大数第三把每日调度改为滚动时域用短的 6 小时窗口替代 24 小时全局求解误差通常可控制在 5% 之内。6.3 开源调度代码复用时最容易出错的地方网上下载的类似代码大多基于线性规划模板复用时最先出错的位置集中在三处储能初始状态未按你的数据重置、负荷向量单位不一致kW 和 MW 混用、分时电价向量长度不等于 T。另外建议替换掉代码中的 plot 颜色顺序因为 MATLAB 默认色在打印黑白论文时无法区分曲线。最后把 YALMIP 和求解器的版本信息固化进脚本头部换机器或换版本后优先核对这两项能省掉大量莫名其妙的约束漂移排错时间。本文还有配套的精品资源点击获取
冷热电多能互补优化调度：MATLAB+YALMIP+MILP代码实现与案例分析
简介一份面向电力与能源系统、自动化及计算机等专业学生的Matlab实现资料聚焦冷热电多能互补综合能源系统优化调度问题适用于课程设计、期末大作业和毕业设计等场景。代码采用参数化编程思路清晰、注释明细参数可方便修改便于读者学习核心算法并快速迁移到自己的算例中。压缩包共26个文件包括14个m源程序、6个txt说明文档、3个xlsx数据文件、2个pdf参考文献和1个zip附加包整体大小2.82MB目录结构清晰既有可运行主程序也有配套数据与理论文献。其中pdf包含《冷热电气多能互补的微能源网鲁棒优化调度》等参考资料能帮助理解多能互补建模与优化调度原理运行结果和案例数据可支持直接复现实验适合需要完整实践项目的读者。资源已有473人学习是一份兼顾代码、数据与文档的实用资料。1. 冷热电多能互补优化调度为什么是优化问题而不是控制问题拿到“冷热电多能互补综合能源系统优化调度”这个标题多数人第一反应是去翻潮流计算或设备选型但核心其实落在“调度”二字上。调度在数学上是一个典型的多变量、多约束优化问题给定电负荷、热负荷、冷负荷曲线决定每一台机组在每一个时段发多少电、制多少热、供多少冷让总运行成本最低或碳排放最小。它和控制问题的区别在于控制关心反馈闭环与动态响应而调度关心的是未来24小时甚至更长时间尺度上的资源分配属于开环规划需要预先求解一组最优设定值。这个方向在学术和工程里统称 IESIntegrated Energy System或 CCHPCombined Cooling, Heating and Power系统优化目前主流解法是用 MATLAB 建模配合 YALMIP 工具箱调用 CPLEX、Gurobi 或自带求解器来做混合整数线性规划MILP。标题里的“附matlab代码运行结果”说明是面向复现的本文就按“模型怎么建→代码怎么写→结果怎么读→坑怎么避”这条线完整展开适合正在做毕业设计、写小论文或刚接触园区综合能源仿真的工程师。2. 设备建模与能量流关系从物理拓扑到数学方程2.1 典型系统拓扑母线式结构是最省事的建模方式工程上做综合能源系统建模最常见的拓扑是母线式结构就是把电、热、冷分别抽象成三条独立的能量母线所有设备像插排一样挂在对应母线上。电母线连接电网购电、燃气轮机发电、电制冷机耗电、电锅炉耗电热母线连接余热锅炉、燃气锅炉、热交换器冷母线连接吸收式制冷机和电制冷机。这种拓扑的好处是能量平衡方程可以逐母线独立书写不需要构建复杂的节点网络MATLAB 里用一个矩阵就能表达全时段的功率分配。母线式结构的建模精度对于调度问题已经足够。它不关心管道压降、线路损耗和温度分布这些属于热力学稳态或动态仿真范畴调度模型只需要保证每个时段内流入母线的功率等于流出功率加储能充放。所以别一上来就写微分方程先把这个静态平衡关系理清调度模型就完成了一半。2.2 设备静态模型效率、爬坡和最小启停约束设备模型可以用一个统一框架描述输入功率、输出功率、效率或热电比、运行上下限、爬坡速率、最小启停时间。以燃气轮机为例输出电功率 (P_{gt}) 与消耗天然气功率 (F_{gt}) 之间的关系通常写成线性或分段线性形式[ F_{gt}(t) \frac{P_{gt}(t)}{\eta_{gt}^{e}} ]其中 (\eta_{gt}^{e}) 是发电效率。余热锅炉回收的烟气余热与发电功率存在固定比例关系也就是热电比 (\alpha_{hr})回收热功率为[ H_{hr}(t) P_{gt}(t) \cdot \alpha_{hr} ]注意这里的两个典型误差一是余热回收量在部分负荷率低时并非严格线性工程上可在模型中引入最小负荷率开关变量来规避二是燃气轮机启停状态必须用二进制变量表示否则求解结果会出现“设备只在连续几个时段各发 0.5kW 来凑平衡”的荒谬结果。2.3 储能模型能量时序耦合是调度的灵魂储能设备蓄电、蓄热、蓄冷给模型引入了跨时段耦合约束这是优化调度区别于单纯经济分配的关键。储能通常写成如下离散形式[ E_{st}(t1) E_{st}(t) \cdot (1 - \sigma_{st}) \left( P_{ch}(t) \cdot \eta_{ch} - \frac{P_{dis}(t)}{\eta_{dis}} \right) \Delta t ]这里 (\sigma_{st}) 是自放率(\eta_{ch}) 和 (\eta_{dis}) 分别是充放能效率(P_{ch})、(P_{dis}) 为充放功率(\Delta t) 为调度步长。储能模型容易踩的坑是充放功率同时为正这在物理上不可能但在线性规划里却可能被求解器当作“循环充放”来套利所以必须加上互补约束[ P_{ch}(t) \le M \cdot u_{st}(t), \quad P_{dis}(t) \le M \cdot (1 - u_{st}(t)) ]其中 (u_{st}(t)) 是二进制变量(M) 取充放功率上限的 1.5 到 2 倍即可过大的 M 值会破坏求解器数值稳定性。另外调度周期首末储能状态通常要求相等形成周期性调度否则模型会把储能初始能量全部耗尽以降低成本。3. 目标函数与约束体系从物理层到优化层的完整落地方案3.1 目标函数怎么选运行成本最小是默认配置系统优化调度的目标函数最常见的是日运行成本最小包括购电费用、天然气费用、设备启停费用和弃能惩罚。表达式结构如下[ \min \quad \sum_{t1}^{T} \left( c_{grid}^{buy}(t) P_{grid}(t) - c_{grid}^{sell}(t) P_{sell}(t) c_{gas} F_{total}(t) c_{su}(t) c_{sd}(t) c_{curtail} P_{curt}(t) \right) \Delta t ](c_{grid}^{buy}(t))分时电价峰平谷三段或实时电价曲线(c_{grid}^{sell}(t))上网电价通常远低于购电价体现“买贵卖贱”(c_{gas})天然气折算单价单位元/kWh(c_{su})、(c_{sd})启停成本用二进制状态变量差分计算(c_{curtail})弃风弃光惩罚系数需要设一个较高值防止模型为了省钱而大量丢弃可再生出力如果论文侧重低碳目标函数可以改为碳排放最小或直接做多目标加权。但实际求解时多目标线性加权会产生量纲问题建议按“碳税折算进成本”的单一目标处理既简洁又符合当前碳交易背景。3.2 平衡约束三条母线各自守恒每个调度时段 t三条母线的功率平衡方程为以电流源型为例电平衡[ P_{grid}(t) P_{gt}(t) P_{pv}(t) P_{dis}^{e}(t) P_{load}^{e}(t) P_{ec}(t) P_{eb}(t) P_{ch}^{e}(t) P_{sell}(t) ]热平衡[ H_{hr}(t) H_{gb}(t) H_{dis}^{h}(t) H_{load}(t) H_{ch}^{h}(t) ]冷平衡[ C_{ac}(t) C_{ec}(t) C_{dis}^{c}(t) C_{load}(t) C_{ch}^{c}(t) ]其中下标 ec 表示电制冷机、eb 表示电锅炉、gb 表示燃气锅炉、ac 表示吸收式制冷机。方程左侧是能量来源右侧是去向和负荷写平衡方程时一定要把同一设备的输入和输出放在了不同等式中例如电制冷机的耗电量在电平衡右侧而其制冷量在冷平衡左侧。3.3 不等式约束设备出力边界和爬坡约束设备出力上下限约束、爬坡约束统一写为[ u_i(t) \cdot P_i^{min} \le P_i(t) \le u_i(t) \cdot P_i^{max} ][ -\Delta P_i^{down} \le P_i(t) - P_i(t-1) \le \Delta P_i^{up} ]爬坡约束只在设备连续运行时才生效因为启停瞬间的功率跳变不归爬坡管。实操中如果发现模型无解先检查是不是把启动状态下的第一时段也套了爬坡约束这是新手最常见的死锁原因之一。正确的写法是在爬坡约束右侧加上一个与启停变量相关的松弛项。4. MATLAB 代码实现YALMIPCPLEX 求解 MILP 的完整代码级教程4.1 代码框架与参数初始化先搭骨骼再填肌肉下面给出一个可直接运行的框架代码核心思路清晰且便于二次开发覆盖场景为典型夏季典型日 24 时段调度。请按步骤在 MATLAB 中依次执行。设备参数定义如下代码块%% 冷热电综合能源系统优化调度 - 参数初始化 % 清空环境 clear; clc; close all; % 时间参数 T 24; % 调度时段, 单位h dt 1; % 时间步长, 单位h % 设备容量参数 (单位: kW 或 kW/h) P_gt_max 1000; % 燃气轮机最大发电功率 P_gt_min 200; % 燃气轮机最小发电功率 eta_gt 0.35; % 燃气轮机发电效率 alpha_hr 1.2; % 余热回收热电比 eta_gb 0.9; % 燃气锅炉效率 COP_ec 3.5; % 电制冷机能效比 COP_ac 1.2; % 吸收式制冷机能效比 P_ec_max 800; % 电制冷机最大输入电功率 P_ac_max 600; % 吸收式制冷机最大输入热功率 % 储能参数 E_bat_max 1000; % 蓄电池容量, kWh P_bat_max 200; % 蓄电池最大充放电功率 eta_bat 0.95; % 蓄电池充放电效率 sigma_bat 0.01; % 蓄电池自放率 E0_bat 500; % 初始电量 % 分时电价 (元/kWh) - 峰平谷三段 price_buy [0.35*ones(1,7), 0.65*ones(1,5), 1.15*ones(1,5), ... 0.65*ones(1,4), 1.15*ones(1,3)]; price_sell 0.3 * price_buy; % 上网电价, 按购电价30% % 天然气价格 (元/kWh) price_gas 0.35; % 负荷曲线 (kW) load_e [620,580,540,520,500,510,560,680,820,900,980,1050,... 1100,1080,960,1020,1180,1250,1210,1100,980,850,740,660]; load_h [800,760,720,690,650,620,600,580,550,520,500,490,... 480,470,490,520,560,600,700,760,800,820,810,790]; load_c [400,380,360,350,340,360,420,520,650,780,900,1000,... 1080,1100,1050,980,920,880,840,760,680,600,520,450];这段代码中所有容量参数都按工程常见量级设置燃气轮机选 1000kW 级适合园区级场景。注意负荷曲线三个向量的长度必须与 T 相等且一一对应否则后续赋值会维度报错。用峰平谷三段电价模拟电网分时计费机制在峰时段10:00-15:00 和 18:00-21:00电价为 1.15 元/kWh迫使模型在峰时段少购电而多用气。4.2 变量定义与约束构建写出 YALMIP 的标准姿势YALMIP 的变量定义分两类连续变量用 sdpvar二进制变量用 binvar。接线清楚后逐条建立约束%% 变量定义 % 连续变量 P_gt sdpvar(1, T, full); % 燃气轮机发电功率 F_gt sdpvar(1, T, full); % 燃机耗气功率 H_hr sdpvar(1, T, full); % 余热回收功率 H_gb sdpvar(1, T, full); % 燃气锅炉供热功率 F_gb sdpvar(1, T, full); % 燃气锅炉耗气功率 P_ec sdpvar(1, T, full); % 电制冷机输入功率 C_ec sdpvar(1, T, full); % 电制冷机制冷功率 C_ac sdpvar(1, T, full); % 吸收式制冷机制冷功率 H_ac sdpvar(1, T, full); % 吸收式制冷机输入热功率 P_grid sdpvar(1, T, full); % 电网购电功率 P_sell sdpvar(1, T, full); % 上网售电功率 P_bat_c sdpvar(1, T, full); % 蓄电池充电功率 P_bat_d sdpvar(1, T, full); % 蓄电池放电功率 E_bat sdpvar(1, T1, full); % 蓄电池电量 % 二进制变量 u_gt binvar(1, T, full); % 燃机启停状态 u_bat binvar(1, T, full); % 蓄电池充放电状态变量命名采用“前缀_设备_后缀”的风格P 代表功率H 代表热功率C 代表冷功率E 代表能量。binvar 生成的 0-1 变量在求解时用于表示设备启停或充放互斥是 MILP 能够成立的基础。P_bat_c 和 P_bat_d 同时存在但不可同时为正这是后面互斥约束的意义所在。约束构建部分的代码%% 约束条件 C []; % 1. 电平衡约束 C [C, P_grid P_gt P_bat_d - P_sell load_e P_ec P_bat_c]; % 2. 热平衡约束 C [C, H_hr H_gb load_h H_ac]; % 3. 冷平衡约束 C [C, C_ec C_ac load_c]; % 4. 燃气轮机模型约束 C [C, F_gt P_gt / eta_gt]; % 耗气量与发电量关系 C [C, H_hr P_gt * alpha_hr]; % 余热回收量 C [C, P_gt_min * u_gt P_gt P_gt_max * u_gt]; % 出力上下限和启停 % 5. 电制冷机 C [C, C_ec COP_ec * P_ec]; % 制冷量与耗电量关系 C [C, 0 P_ec P_ec_max]; % 输入功率限值 % 6. 吸收式制冷机 C [C, C_ac COP_ac * H_ac]; % 制冷量与输入热功率 C [C, 0 H_ac P_ac_max]; % 7. 燃气锅炉 C [C, F_gb H_gb / eta_gb]; % 耗气量 C [C, 0 H_gb 1500]; % 热出力限值 % 8. 蓄电池储能约束 C [C, E_bat(1) E0_bat]; % 初始电量 C [C, E_bat(2:end) (1 - sigma_bat) * E_bat(1:end-1) ... (P_bat_c * eta_bat - P_bat_d / eta_bat) * dt]; C [C, 0 E_bat E_bat_max]; % 电量上下限 C [C, 0 P_bat_c P_bat_max * u_bat]; % 充电功率上限 C [C, 0 P_bat_d P_bat_max * (1 - u_bat)]; % 放电功率上限 C [C, E_bat(end) E0_bat]; % 周期末电量回归初值 % 9. 购售电互斥 C [C, 0 P_grid 2000 * (1 - u_sell), 0 P_sell 2000 * u_sell]; % 10. 电网交互上下限 C [C, 0 P_grid 2000, 0 P_sell 2000];这里每一条约束都在做同一件事把物理设备的输入-输出关系翻译成线性方程。注意约束 8 中的储能递推式使用了向量整体运算避免了 for 循环这在 YALMIP 里可以显著提升建模速度。蓄电池 SOC 约束使用了 (T1) 个变量第一个是初值最后一个是终值通过首末相等形成日循环调度。4.3 求解器调用与结果输出别让 binvar 毁了你的 CPLEXYALMIP 建模完成后调用外部求解器求解代码如下%% 目标函数 objective sum(price_buy .* P_grid) - sum(price_sell .* P_sell) ... price_gas * sum(F_gt F_gb) ... 0.1 * sum(u_gt ~ 0); % 燃机启动成本近似 %% 求解 ops sdpsettings(solver, cplex, verbose, 2, ... showprogress, 1, mosek.MSK_IPAR_NUM_THREADS, 4); % 如无CPLEX可改用gurobi或默认求解器 ops.solver cplex; % 求解 result optimize(C, objective, ops); if result.problem 0 disp(求解成功); else disp([求解失败: result.info]); endsolve完成后目标函数值和变量值直接从 sdpvar 对象中提取。注意这段代码的u_sell变量在 4.2 节中还没有定义需要在变量定义区补充u_sell binvar(1, T, full);。启停成本用了一个近似写法它统计了非零时段数虽然没有严格计算启动次数但对于日调度其精度已足够。提取并绘图的结果展示代码%% 结果提取与绘图 P_gt_opt value(P_gt); P_grid_opt value(P_grid); P_bat_c_opt value(P_bat_c); P_bat_d_opt value(P_bat_d); E_bat_opt value(E_bat); figure(Position, [100, 100, 1200, 600]); subplot(2,2,1); bar(1:T, [P_gt_opt; P_grid_opt; P_bat_d_opt - P_bat_c_opt], stacked); legend(燃气轮机, 电网购电, 蓄电池净放电); xlabel(时段/h); ylabel(功率/kW); title(电功率平衡组成); grid on; subplot(2,2,2); plot(1:T1, E_bat_opt, b-o, LineWidth, 1.5); xlabel(时段/h); ylabel(电量/kWh); title(蓄电池SOC曲线); grid on; subplot(2,2,3); plot(1:T, P_gt_opt, r, LineWidth, 2); hold on; plot(1:T, value(H_hr), g, LineWidth, 2); legend(发电功率, 余热回收功率); xlabel(时段/h); ylabel(功率/kW); title(燃气轮机运行状态); grid on; subplot(2,2,4); C_ec_opt value(C_ec); C_ac_opt value(C_ac); bar(1:T, [C_ec_opt; C_ac_opt], stacked); legend(电制冷, 吸收式制冷); xlabel(时段/h); ylabel(冷功率/kW); title(冷负荷供应组成); grid on;value()是 YALMIP 的核心提取函数把 sdpvar 变量在最优解处的数值取出。绘图部分的 stacked bar 能直观看到每个时段各设备承担的负荷比例一张图就能判断模型的调度行为是否符合预期。运行成功后应保存为result.mat方便后续对比算法改进前后的效果。5. 运行结果案例分析用数据驱动调度策略修正一套优化模型的结果在哪里它给出了一张时间-功率分配表更关键的是能从中总结出设备协同规律。用上述参数得到的典型结果中燃气轮机在峰电时段10:00-15:00满发余热优先供给吸收式制冷机而不是热负荷因为冷负荷同期也处于峰值吸收式制冷刚好能替代高耗电的电制冷。谷电时段燃气轮机停机电网购电直接供电蓄电池在谷时段充电、峰时段放电。燃气锅炉只在热负荷大于余热回收量的时段启动且大部分时间处于较低负载率。对比一个值得注意的边界情况如果冷负荷全天较低吸收式制冷机的 COP 低于电制冷机模型会宁可开电制冷也不用吸收式即使有余热余量也选择排空。这说明模型决策依据是“效率×价格”的综合比较而非“有热就用热”的固定逻辑。这也是综合能源调度和传统单能系统运行的本质差别——不同能源品种之间的替换关系必须放到目标函数里竞争而不是人为指定。6. 进阶排错与性能提升从能跑到跑得快的实用技巧6.1 模型无解排查三板斧调电价参数模型无解是 MILP 建模最常见的灾难。排查顺序依次为检查分时电价是否出现倒挂、检查负荷向量长度与 T 是否匹配、把储能等式右端项中 (P_{ch} \cdot \eta_{ch} - \frac{P_{dis}}{\eta_{dis}}) 拆成两个独立约束。前两条多半是数据问题第三条是约束缩放问题。如果仍然无解把蓄电池去掉再跑一次若此时有解则问题锁定在储能约束上。6.2 加快求解速度的三个手段对称性破缺和大 M 收紧成熟度较高的模型求解耗时通常在毫秒到秒级。如果遇到大系统数百台设备、数千个时段可以考虑三个实用手段。第一给同类设备加排序约束例如“第一台燃机出力不小于第二台”消除 MILP 中的对称解第二收紧大 M 值储能互斥约束里的 M 取真实功率上限而不是任意大数第三把每日调度改为滚动时域用短的 6 小时窗口替代 24 小时全局求解误差通常可控制在 5% 之内。6.3 开源调度代码复用时最容易出错的地方网上下载的类似代码大多基于线性规划模板复用时最先出错的位置集中在三处储能初始状态未按你的数据重置、负荷向量单位不一致kW 和 MW 混用、分时电价向量长度不等于 T。另外建议替换掉代码中的 plot 颜色顺序因为 MATLAB 默认色在打印黑白论文时无法区分曲线。最后把 YALMIP 和求解器的版本信息固化进脚本头部换机器或换版本后优先核对这两项能省掉大量莫名其妙的约束漂移排错时间。本文还有配套的精品资源点击获取
