简介这份压缩包围绕考虑用户舒适度的冷热电多能互补综合能源系统优化调度问题提供了一套基于Matlab 2014/2019a的完整仿真代码与配套资料面向电气工程、能源动力与控制方向的本科生、硕士生及科研人员适合用于课程设计、论文复现或作为综合能源系统优化调度研究方向的基础参考。压缩包共24个文件总大小约4.95MB核心为4个.m脚本覆盖经济性优化与排放优化两种目标的主程序并配有12张运行结果图、3个程序思路说明文档、2个Excel数据文件以及2篇关于冷热电气多能互补微能源网鲁棒优化调度的PDF论文和1个代码说明Word文档便于结合图、表和文献快速复核运行细节。资源目前已有265人学习下载。借助这份资料不仅可以系统学习用户舒适度约束下冷热电联供系统的多能互补建模方法、优化目标与运行策略还可以直接修改运行、观察不同目标函数对调度结果的影响为深入掌握智能优化算法在电力系统中的应用打下基础。1. 考虑用户舒适度的冷热电多能互补系统优化调度到底在优化什么冷热电多能互补综合能源系统的优化调度表面上是把一个园区或建筑群的冷、热、电三种能源需求分摊给光伏、储能、燃气轮机、吸收式制冷机、电制冷机等设备去满足让总运行成本最低。但真正让调度模型从“能跑”变成“能落地”的往往不是设备效率曲线拟合得有多准而是用户舒适度约束怎么进模型。一个只盯着电费单的调度策略很可能在夏季午后把室内温度压到 22℃在冬季清晨把热水温度顶到 60℃账面上省了钱用户侧投诉却把收益全部吃掉。这个问题在工程上被拆成两层第一层是能量平衡即每一时刻冷、热、电的供需必须守恒第二层是舒适度边界即室内温度、相对湿度、热水温度等用户体验指标必须落在可接受区间内。前者是等式约束后者是不等式约束两者叠加后调度问题的可行域被大幅收窄求解难度也随之上升。本文按“建模→目标→求解→MATLAB 实现→验证”这条线把一套可复现的优化调度方案讲清楚适合正在做微网优化调度、综合能源系统日前调度的工程师也适合拿 MATLAB 做毕设或竞赛比如电工杯微电网调度类题目的学生参考。2. 冷热电多能互补系统的设备模型与能量平衡约束2.1 系统架构与能量流先画图再写方程典型的冷热电多能互补系统包含三条能量母线电力母线、热力母线、冷力母线。电源侧由光伏、市电、燃气轮机或燃气锅炉组成热力侧由燃气轮机余热回收、燃气锅炉、储热罐组成冷力侧由吸收式制冷机、电制冷机、蓄冷罐组成。用户侧则同时存在电负荷、热负荷和冷负荷。调度的本质是在每一个调度时段通常取 1 小时内决定各设备的出力值使得三条母线上的供需平衡。在 MATLAB 中建模的第一步不是写目标函数而是先把设备出力定义为决策变量。常见做法是构造一个二维决策变量矩阵行索引对应设备列索引对应时段。例如对一个含光伏、燃气轮机、电制冷机、吸收式制冷机、燃气锅炉、储热罐、蓄冷罐的系统决策变量可以组织为% 设备出力矩阵维度7个设备 x 24个时段 nDevice 7; nPeriod 24; x sdpvar(nDevice, nPeriod, full); % 行1: 光伏出力; 行2: 燃气轮机发电; 行3: 电制冷机耗电; % 行4: 吸收式制冷机制冷量; 行5: 燃气锅炉产热; 行6: 储热罐充放热(正放负充); % 行7: 蓄冷罐充放冷(正放负充)这里使用 YALMIP 的sdpvar声明优化变量full表示变量矩阵无特殊结构。之所以把充放能设备放在同一行变量里用正负号区分是为了减少变量数量、加快求解速度。如果拆成两个变量虽然物理含义更直白但会引入额外的二进制变量来防止同时充放模型规模会明显膨胀。2.2 能量平衡约束等式约束怎么写才不漏项电平衡约束要求光伏出力 燃气轮机发电 市电购电 蓄冷罐放冷等效耗电如果有 储热罐相关耗电 电负荷 电制冷机耗电 其他电耗。写成 MATLAB 约束如下% 电平衡光伏 燃气轮机 购电 电负荷 电制冷耗电 Constraints [Constraints, ... x(1, :) x(2, :) x_powerGrid P_e_load x(3, :)]; % 热平衡燃气锅炉产热 燃气轮机余热 储热罐放热 热负荷 Constraints [Constraints, ... x(5, :) eta_hr * x(2, :) x(6, :) P_h_load]; % 冷平衡吸收式制冷 电制冷 蓄冷罐放冷 冷负荷 Constraints [Constraints, ... x(4, :) COP_ec * x(3, :) x(7, :) P_c_load];燃气轮机的余热回收量通常与发电量成正比用eta_hr表示余热回收系数电制冷机的制冷量等于耗电量乘以能效比COP_ec。注意热平衡中x(6, :)的正负号——正值代表放热负值代表蓄热。这种写法下储热罐自身损耗可以放在容量约束里考虑简化模型但保留主要物理特征。2.3 设备出力上下限与爬坡约束每台设备的出力不能超过额定容量这组约束直接写成变量边界即可。燃气轮机和燃气锅炉还有爬坡约束即相邻两个时段的出力变化不能过于剧烈否则会加速设备磨损甚至触发保护停机。% 设备出力上下限 Constraints [Constraints, ... 0 x(1, :) P_pv_max ... % 光伏0到额定容量 0 x(2, :) P_gt_max ... % 燃气轮机发电上限 0 x(3, :) P_ec_max ... % 电制冷机耗电上限 0 x(4, :) P_ac_max ... % 吸收式制冷机冷量上限 0 x(5, :) P_gb_max]; % 燃气锅炉产热上限 % 燃气轮机爬坡约束相邻时段出力变化不超过额定容量的30% ramp_rate 0.3 * P_gt_max; Constraints [Constraints, ... -ramp_rate diff(x(2, :)) ramp_rate];diff(x(2, :))计算出 24 个时段中相邻两个时段的出力差是一个长度为 23 的向量。爬坡约束对燃气轮机这类响应速度较慢的设备是必须的但对光伏和电制冷机一般可以省略因为前者出力由光照决定、不可控后者响应速度快、爬坡限制远小于其物理能力。3. 用户舒适度约束的建模方法PMV 指标与温度区间3.1 为什么舒适度必须进约束而不是进目标有些文献把舒适度作为目标函数的一项和运行成本加权求和。这种做法在数学上更灵活却有一个工程上的尴尬成本和舒适度的量纲不同权重系数很难标定。权重给大了调度结果会为了 0.1℃ 的舒适度提升多花几百元电费权给小了舒适度约束形同虚设。我一般把舒适度写成硬约束——设定一个可接受的舒适度区间限制条件内让成本自由优化。这样做还有一个额外的好处模型仍然是线性约束不引入非线性项可以直接用线性规划或混合整数线性规划求解器处理求解速度和收敛性都有保障。3.2 PMV 指标的线性化近似最常用的舒适度指标是 PMVPredicted Mean Vote预测平均投票数取值在 -3 到 3 之间0 代表中性最舒适正值为偏热负值为偏冷。PMV 的原始计算式高度非线性包含人体代谢率、服装热阻、空气流速等参数直接放进优化模型会变成非线性约束。工程上的做法是在给定室内风速和服装热阻的条件下把 PMV 近似为室内温度和相对湿度的线性函数% PMV线性近似T_in为室内温度, RH为相对湿度 % 系数a1, a2, c由实验数据回归得到典型值约 a10.21, a20.03, c-7.5 PMV a1 * T_in a2 * RH c; % 舒适度约束PMV在[-0.5, 0.5]区间内 Constraints [Constraints, ... -0.5 PMV 0.5];这里室内温度T_in不是决策变量而是由建筑热动态模型决定的中间变量。最常见的一阶热动态模型为% 一阶热动态室内温度变化 供热/供冷功率 - 围护结构散热 % C_build为建筑热容, T_out为室外温度, R_wall为围护结构热阻 T_in(:, t1) T_in(:, t) (Q_hvac(:, t) - (T_in(:, t) - T_out(:, t)) / R_wall) / C_build;Q_hvac 是暖通设备这里即冷热联供设备提供给室内的热/冷量。注意这个方程在 MATLAB 中要用逐时段递推的方式展开成约束而不是写成循环赋值——优化问题中变量之间的时序关系要靠约束表达不能靠循环计算。3.3 舒适度约束与设备出力的耦合舒适度约束引入后设备出力不再是独立变量。电制冷机的耗电量、吸收式制冷机的制冷量会直接影响室内温度的变化轨迹进而影响 PMV 是否越界。这意味着电平衡、热平衡、冷平衡三条母线方程和舒适度约束是强耦合的。求解器需要在满足所有耦合约束的前提下找到成本最低的调度方案。这一耦合关系是问题难度的主要来源之一。如果系统不含蓄冷蓄热设备那么每个时段的冷热出力只影响该时段的温度24 个时段可以解耦成 24 个小问题。一旦加入储热罐和蓄冷罐温度就成了跨时段的“记忆量”今天下午蓄的冷可能在晚间释放整个时间维度必须作为一个整体求解模型规模变为原来的约 24 倍这也是为什么这类问题通常需要商业求解器而不是手写梯度下降法。4. 目标函数与求解策略最低成本还是最低碳排放4.1 运行成本目标函数的完整构成优化调度的目标函数一般取系统日运行成本最小包括燃料成本、购电成本、设备启停成本和运维成本四部分。燃料成本针对燃气轮机和燃气锅炉购电成本针对市电交互。若按分时电价购电成本表达式如下% 目标函数总运行成本最小 % 购电成本分时电价 x 购电量 cost_power sum(price_buy .* x_powerGrid); % 燃料成本燃气轮机 燃气锅炉单位热值价格 x 消耗功率 fuel_gt sum(gas_price * x(2, :) / eta_gt); % eta_gt为发电效率 fuel_gb sum(gas_price * x(5, :) / eta_gb); % eta_gb为锅炉效率 % 运维成本按出力比例估算系数一般为0.01-0.05元/kWh cost_om sum(om_gt * x(2, :)) sum(om_ec * abs(x(3, :))) sum(om_gb * x(5, :)); % 总目标 Objective cost_power fuel_gt fuel_gb cost_om;分时电价向量price_buy长度为 24可直接从当地电价政策表写入。燃气轮机的燃料成本按发电量除以发电效率折算成天然气消耗量再乘以气价。运维成本通常远小于燃料成本但不可省略因为它影响了设备间的出力分配——两台设备效率接近时运维成本低的优先出力。4.2 求解器选择CPLEX、Gurobi 还是默认求解器目标函数是线性的约束也是线性的这是一个标准的线性规划问题。但如果引入设备启停状态即某台设备最小运行时间的约束就需要引入二进制变量问题退化为混合整数线性规划。求解器支持问题类型许可方式MATLAB 接入方式CPLEXLP / MILP / QP商业授权yalmip(solver,cplex)GurobiLP / MILP / QP商业授权yalmip(solver,gurobi)GLPKLP / MILP开源yalmip(solver,glpk)MATLAB linprogLP随 MATLAB 提供直接用 linprogintlinprogMILP随 MATLAB 提供直接用 intlinprog对于 24 时段、7 个连续变量加若干二进制变量的模型linprog或intlinprog完全够用求解时间通常在秒级。只有把规模扩大到数百个节点、分钟级调度粒度时才需要 CPLEX 或 Gurobi。YALMIP 的优势是建模层与求解层解耦换求解器只需改一行配置。4.3 用 YALMIP 组合完整模型的骨架把前面的模型片段拼起来一个完整的日前优化调度骨架如下% 初始化YALMIP yalmip(clear); % 决策变量7个设备x24时段 购电量 x sdpvar(7, 24, full); x_powerGrid sdpvar(1, 24, full); % 约束集合 Constraints []; % 能量平衡约束代码见2.2节 % 设备出力约束代码见2.3节 % 舒适度约束代码见3.2节 % 储热/蓄冷罐容量约束 Constraints [Constraints, ... SOC_tank_min cumsum(x(6, :)) / cap_tank SOC_tank_max ... SOC_cold_min cumsum(x(7, :)) / cap_cold SOC_cold_max]; % 目标函数代码见4.1节 Objective cost_power fuel_gt fuel_gb cost_om; % 求解 ops sdpsettings(solver, linprog, verbose, 0); result optimize(Constraints, Objective, ops); % 结果检查 if result.problem 0 disp(调度方案求解成功); x_opt value(x); else disp([求解失败: result.info]); end储热罐容量约束用cumsum累加各时段的充放热量表示罐内剩余能量始终处在上下限之间。这段代码直接拷贝到 MATLAB 中配上参数定义即可运行。注意value(x)必须放在optimize成功后调用否则拿到的是初值。5. 冷热电联供微网调度在 MATLAB 中的实现细节5.1 数据准备负荷曲线与气象数据怎么组织调度模型的输入包括三类数据用户电、热、冷负荷曲线室外温度曲线以及光伏出力预测曲线。这三类数据在 MATLAB 中统一组织为 24 维行向量便于直接参与矩阵运算。% 夏季典型日负荷曲线 (kW) P_e_load [230 215 205 200 210 240 280 320 345 360 350 330 ... 320 335 340 350 365 360 340 310 290 270 250 235]; P_h_load [80 75 70 65 60 55 50 45 40 35 30 35 ... 40 45 50 55 45 40 35 30 25 30 45 60]; % 热负荷系数乘以面积 P_c_load [50 45 40 35 30 60 95 130 155 170 165 155 ... 150 155 160 165 175 180 170 140 110 85 70 60]; % 冷负荷 % 室外温度 (℃) T_out [26 25 24 23 23 24 27 29 31 33 34 35 ... 36 36 35 35 34 33 32 30 29 28 27 26]; % 光伏出力预测 (kW)归一化后乘以装机容量 pv_norm [0 0 0 0 0.02 0.08 0.15 0.25 0.38 0.50 0.62 0.72 ... 0.78 0.75 0.65 0.52 0.38 0.22 0.10 0.02 0 0 0 0]; P_pv 300 * pv_norm; % 光伏装机300kW冷负荷的峰值出现在 16:00 到 18:00这是因为建筑围护结构的蓄热效应使得冷负荷峰值滞后于室外温度峰值约 2 到 3 小时。热负荷则在早晚出现双峰对应人员上下班后的热水需求。5.2 参数表与变量命名规范模型参数集中管理在脚本开头避免在约束代码中直接写魔法数。设备参数建议按以下方式组织%% 设备参数 P_gt_max 150; % 燃气轮机额定发电容量 (kW) eta_gt 0.35; % 燃气轮机发电效率 eta_hr 0.45; % 余热回收效率 P_gb_max 200; % 燃气锅炉额定产热容量 (kW) eta_gb 0.85; % 燃气锅炉效率 P_ec_max 120; % 电制冷机额定耗电功率 (kW) COP_ec 3.5; % 电制冷机能效比 P_ac_max 180; % 吸收式制冷机额定冷量 (kW) COP_ac 1.2; % 吸收式制冷机能效比热驱动 cap_tank 300; % 储热罐容量 (kWh) cap_cold 250; % 蓄冷罐容量 (kWh) %% 舒适度参数 C_build 800; % 建筑等效热容 (kWh/℃) R_wall 0.08; % 围护结构等效热阻 (℃/kW) a1 0.21; % PMV线性化温度系数 a2 0.03; % PMV线性化湿度系数 c -7.5; % PMV线性化常数项 %% 价格参数 price_buy [0.48 0.48 0.48 0.48 0.48 0.58 0.68 0.78 0.88 0.88 0.88 0.88 ... 0.88 0.88 0.88 0.78 0.68 0.78 0.88 0.88 0.88 0.68 0.58 0.48]; gas_price 2.8; % 天然气价格 (元/m³)按热值折算后为元/kWh om_gt 0.02; % 燃气轮机运维成本 (元/kWh) om_gb 0.015; % 燃气锅炉运维成本 (元/kWh) om_ec 0.01; % 电制冷机运维成本 (元/kWh)参数的选取直接影响调度结果。燃气轮机效率取 0.35 时发电成本约为 0.8 元/kWh高于低谷电价但低于高峰电价因此燃气轮机通常在电价高峰时段满发低谷时段停运。电制冷机 COP 取 3.5 时制冷成本约为 0.14 元/kWh远低于吸收式制冷机由天然气驱动成本约 0.33 元/kWh所以只要电价不是极端高系统倾向于优先用电制冷机。5.3 运行结果的可视化与数据导出求解完成后把结果画成调度堆叠图是最直观的验证方式。堆叠图能明显看出各设备在不同时段的出力分配是否合理以及舒适度约束是否真的被激活。% 调度结果绘图 figure; subplot(3,1,1); area([x_opt(1,:) x_opt(2,:) x_powerGrid], LineWidth, 1.2); legend(光伏, 燃气轮机, 市电购电, Location, northwest); title(电力平衡调度结果); ylabel(功率 (kW)); subplot(3,1,2); area([x_opt(5,:) eta_hr*x_opt(2,:) x_opt(6,:)], LineWidth, 1.2); legend(燃气锅炉, 余热回收, 储热罐, Location, northwest); title(热力平衡调度结果); ylabel(功率 (kW)); subplot(3,1,3); area([x_opt(4,:) COP_ec*x_opt(3,:) x_opt(7,:)], LineWidth, 1.2); legend(吸收式制冷, 电制冷, 蓄冷罐, Location, northwest); title(冷力平衡调度结果); ylabel(功率 (kW));一个值得注意的现象是如果储热罐的容量约束SOC_tank_min设置得过紧例如要求任何时刻罐内剩余能量不低于 50%燃气轮机在夜间低谷时段也不得不出力为储热罐充热导致成本升高。这时候需要检查储热罐约束是否恰当而不是怀疑求解器出了问题。5.4 求解失败时的排查顺序optimize返回的result.problem值不为 0 时按以下顺序排查先移除舒适度约束看问题是否变可行——如果移除后仍然无解说明能量平衡或设备容量约束有矛盾例如冷负荷峰值大于所有制冷设备总容量如果移除后可行再逐个加回舒适度约束找出是哪条约束导致无解。另一个高频错误是变量维度不匹配sdpvar(7, 24)与长度为 24 的行向量相乘时如果其中一个用了列向量MATLAB 会报维度错误或直接得到错误结果。建议在optimize之前先检查所有Constraints中涉及的变量规模。6. 冷热惯性利用用舒适度区间换调度弹性的一个验证技巧讨论冷热电多能互补系统的舒适度约束最容易被忽视的一点是舒适度区间本身是一种储能。室内温度允许在 24℃ 到 27℃ 之间波动时就等于给了系统一个容量为C_build * 3的热储能只是这个“储能”的功率上限受空调/暖通设备出力限制容量上限受温度区间限制。实际运行中利用建筑热惯性把部分冷负荷从电价高峰时段平移到低谷时段是降低运行成本最有效的手段之一。验证这个结论的方法是做一组对比实验把 PMV 约束从[-0.5, 0.5]放宽到[-1, 1]分别运行优化模型对比两次调度的总成本和设备出力曲线。用 MATLAB 脚本可以自动完成这组对比% 对比不同舒适度区间下的运行成本 pmv_ranges {[-0.5 0.5], [-1 1]}; total_cost zeros(1, length(pmv_ranges)); for i 1:length(pmv_ranges) % 重新构建约束PMV上下限取自pmv_ranges Constraints build_constraints(pmv_ranges{i}); Objective build_objective(); ops sdpsettings(solver, linprog, verbose, 0); result optimize(Constraints, Objective, ops); total_cost(i) value(Objective); % 计算舒适度约束的拉格朗日乘子 lambda dual(Constraints); end % 输出对比结果 fprintf(PMV区间[-0.5,0.5]: 总成本 %.2f 元\n, total_cost(1)); fprintf(PMV区间[-1,1]: 总成本 %.2f 元\n, total_cost(2));dual(Constraints)返回各约束的对偶变量也就是约束的影子价格。影子价格最大的那条约束就是限制成本下降的瓶颈——在夏季场景中通常是 PMV 上限约束在冬季则变成 PMV 下限约束。这个信息直接告诉你如果要进一步降本应该放宽哪个约束或者加装什么设备而不是盲目增大所有设备的容量。另外这套模型稍作修改就可以从日前调度扩展为日内滚动调度把 24 时段改为 96 时段15 分钟粒度每 15 分钟用最新负荷预测数据重新求解一次只执行第一个时段的调度指令。这就是模型预测控制在微网优化调度中的典型应用方式。改动量不大但要注意把舒适度约束改成基于当前实测室内温度的递推形式否则预测误差会随着滚动次数累积导致温度频繁越界。最后提一个从工程现场反馈回来的细节舒适度模型的参数——C_build建筑等效热容和R_wall围护结构热阻——不能从设计图纸直接抄。这两个参数需要根据历史运行数据做参数辨识常见做法是用最小二乘法拟合实测的室内温度变化曲线。参数偏小调度结果会过于乐观频繁启停设备参数偏大温度变化看起来迟钝真实运行时却经常越限。把辨识后的参数代入模型才是一次真正能交给现场运行人员使用的优化调度方案。本文还有配套的精品资源点击获取