冰蓄冷空调多时间尺度优化调度：MATLAB三层滚动控制实现

简介面向电力系统优化调度研究人员这份含冰蓄冷空调的多时间尺度优化调度Matlab源码与配套数据兼顾日前调度与日内电热冷协调调度两大场景可直接用于冷热负荷参与电网调峰的建模与仿真实验。压缩包共15个文件以Matlab脚本为主配以Excel数据表、MAT数据文件以及运行说明文档整体约127KB目录按日前调度、日内电热冷调度、图形表格输出和场景削减等功能模块清晰划分。代码包含主调度程序、图形与表格生成程序、场景削减程序及配套数据并针对日前与日内子模块的运行顺序、Matlab重启要求等常见问题给出提示可帮助学习者快速复现结果并进行二次开发。目前已有106人学习下载适合具备一定Matlab编程基础、希望研究冰蓄冷空调多时间尺度优化调度方法的读者。1. 冰蓄冷空调的多时间尺度优化调度到底在调什么电力现货市场把峰谷价差拉到3倍以上的城市冰蓄冷空调不再只是“夜间制冰、白天放冷”的设备而是建筑冷源侧最灵活的电负荷。但真正跑过调度模型的人都清楚电价套利只是结果不是目标。蓄冰槽容量有限、主机不能同时制冰和供冷、冷负荷预测在24小时内的误差可能超过20%这些约束叠加在一起单靠一个日前线性规划根本扛不住。多时间尺度优化调度要做的事就是把“明天怎么排”拆成“未来4小时怎么调”和“接下来几分钟怎么修正”让每一层只承担自己能承受的预测误差。这篇东西从物理模型、MATLAB源码结构、三层滚动调度实现一直讲到验证方法。适合暖通、源网荷储、建筑能源优化方向的工程师也适合想用MATLAB把MPC落到实际数据的同学。2. 含冰蓄冷空调的物理模型与优化目标先建对变量再谈调度2.1 制冷机组与蓄冰槽的状态变量怎么选冰蓄冷系统建模时最常见的错误是变量定义得太粗。比如只把“蓄冰量”当成一个百分数或者把主机耗电功率当成控制量这样的模型在优化后很难换算成设备指令。我一般分四类变量连续状态量蓄冰槽剩余冷量 $S(k)$单位 kWh有时也用 kJ但必须全系统统一连续控制量主机供冷功率 $Q_{cool}(k)$、主机产冰功率 $Q_{ice}(k)$、融冰供冷功率 $Q_{melt}(k)$0-1 模式变量$m(k)$1 表示主机处于制冷工况0 表示制冰工况时变输入冷负荷预测 $Q_{load}(k)$、分段电价 $price(k)$、室外温度 $T_{out}(k)$。把主机供冷功率和产冰功率分开列是因为双工况主机的 COP 差异很大。有些源码只用一个总功率和一个总制冷量再用比例系数分配最后做模式互斥时非常麻烦。蓄冰槽的离散状态递推方程一般写成$$ S(k1) S(k) \left( \eta_{ice} Q_{ice}(k) - \frac{Q_{melt}(k)}{\eta_{melt}} - \theta S(k) \right) \Delta t $$其中 $\eta_{ice}$ 是制冰效率$\eta_{melt}$ 是融冰效率$\theta$ 是蓄冰槽自损失系数。对于 8 小时蓄冰、白天放冰的系统$\theta$ 取 0.02~0.05/小时如果保温做得好直接取 0 在 24 小时时域内误差不大。但规划时域超过 24 小时时自损失项不能省否则会低估蓄冰槽的补充需求。2.1.1 双工况主机的效率模型制冷机组在制冷和制冰两种工况下 COP 差别很大不能用一个常数。常见做法是给一个分段线性或二次修正% 双工况主机效率模型Q 是负载率对应的实际制冷量 COP_cool (Q) 4.5 - 0.8 * (Q / Q_max) - 0.3 * (Q / Q_max).^2; % 制冷工况 COP_ice (Q) 3.6 - 0.9 * (Q / Q_max) - 0.2 * (Q / Q_max).^2; % 制冰工况 P_e Q_cool / COP_cool(Q_cool) Q_ice / COP_ice(Q_ice);逻辑说明这里把耗电量拆成制冷耗电和制冰耗电两项。实际运行中两项二选一但优化模型里如果直接相加求解器可能让Q_cool和Q_ice同时为正得到一个物理上不存在的“假运行点”所以后面必须加互斥约束。参数说明Q_max是主机额定制冷量单位 kWCOP 表达式里的系数按典型离心机取负载率越低 COP 越低。若源码里只写了固定 COP调度结果会偏乐观白天电价高峰时段主机出力可能被高估真正的瓶颈到执行阶段才会暴露。2.2 目标函数电费最小只是第一层最朴素的优化目标是全周期电费最低$$ \min \sum_{k1}^{N} price(k) \cdot P_e(k) \Delta t $$但直接把这个目标丢给求解器结果往往会出现“每小时都在剧烈启停”的抖动因为电价的每个拐点都会让主机开关一次。电费最小只是骨架工程上要在目标里加三个修正项主机启停惩罚相邻时刻制冷量变化量 $\lambda |Q_{cool}(k1)-Q_{cool}(k)|$弃冷惩罚如果冷负荷预测偏高但蓄冰不足需要切负荷惩罚系数要远高于最高电价峰段主机直供惩罚很多地区的需求响应要求高峰时段主机直供受限可以把峰时段Q_cool的费率乘 1.5~2.0。目标函数写成$$ \min \sum_{k1}^{N} \left( price(k) P_e(k) \lambda |\Delta Q_{cool}(k)| \mu \cdot e_{shed}(k) \right) \Delta t $$其中 $e_{shed}(k)$ 是切冷量辅助变量满足 $e_{shed}(k) \ge Q_{load}(k) - Q_{cool}(k) - Q_{melt}(k)$ 且 $e_{shed}(k) \ge 0$。MATLAB YALMIP 里这样写% 目标函数和惩罚项定义N 为当前优化时域 e_shed sdpvar(N, 1); % 切冷量非负 constraints [constraints, e_shed 0]; constraints [constraints, e_shed Q_load - Q_cool - Q_melt]; objective sum(price .* P_e) * dt ... lambda * sum(abs(diff(Q_cool))) ... mu * sum(e_shed) * dt;逻辑说明e_shed是松弛变量它吸收冷负荷供需不平衡惩罚系数mu要设得足够高否则求解器会在电费便宜时偷偷少供冷。我一般设mu 最高峰时电价 × 2这样切负荷只会在极端情况下出现。参数说明lambda的单位是元/kWdiff(Q_cool)得到的是功率变化量所以lambda的大小要结合dt看。若dt1小时lambda取 0.2~0.5 元/kW 可以平滑波动太小没约束太大会让主机完全不出力蓄冰槽夜间蓄不满白天全靠融冰最后切负荷惩罚反而更高。2.3 多时间尺度的核心理由预测误差分期消纳如果你只有一组精确的 24 小时负荷预测用单层优化就够。现实里冷负荷预测在 24 小时内的误差可达 10%30%尤其夏季午后雷暴、大型会议临时加场这类突变一次性压到日前计划里要么蓄冰不足要么过度蓄冰。多时间尺度调度把决策分成三层每一层面对不同的预测误差水平控制层级预测误差控制周期主要应对手段日前计划10%30%1h确定总蓄冰量和 SOC 参考曲线日内滚动3%8%15min滚动重算未来 4 小时实时反馈1%3%30s60s小幅修正主机功率消除脉动偏差用一句话概括日前负责定“总预算”日内负责把预算重新分配实时负责保证执行不偏轨。所以多时间尺度优化不只是一个线性规划问题而是一个滚动优化 反馈修正问题。MATLAB 里做这个框架的优势是矩阵运算快、YALMIP 建模方便、仿真循环写起来直观源码加数据很快能跑出对比曲线。3. 用MATLAB搭建多时间尺度优化调度的最小可运行框架3.1 数据文件长什么样CSV里的单位决定你后面要不要乘1000一套含冰蓄冷空调的调度数据文件里时间序列是绝对主角。我一般把一天的数据按小时存成 CSV表头如下字段单位说明timestamp-日期-小时load_forecastkW冷负荷预测值load_actualkW实际冷负荷做回测用price元/kWh峰谷分时电价temp℃干球温度用于修正 COP注意很多源码里电价单位是“分/kWh”负荷单位是 MW。如果不统一目标函数数值会差 1000 倍结果就是约束全部失效。加载后第一步做单位校验data readtable(schedule_input.csv, PreserveVariableNames, true); % 统一单位到 kW 和 元/kWh if max(data.load_forecast) 100 data.load_forecast data.load_forecast * 1000; end if max(data.price) 10 data.price data.price / 100; end逻辑说明读入后先看数据范围用启发式规则统一单位能避免“为什么优化结果全是 0”这类问题。PreserveVariableNames在较新的 MATLAB 版本里保留表头原名老版本用VariableNamingRule,preserve否则下划线会被改写成驼峰。参数说明如果数据文件是 15 分钟间隔max(data.load_forecast)的阈值判断要改prices 单位是元/kWh 时通常在 0.1~1.5 之间大于 10 大概率是“分”单位。3.2 用YALMIP把优化问题写成分层MPCYALMIP 是 MATLAB 下的建模工具箱不是求解器。先装好 YALMIP再配一个求解器纯线性问题用 MATLAB 自带的linprog有 0-1 变量时用intlinprog或 Gurobi、CPLEX。基础模型包含三类约束冷负荷平衡、设备出力上下限、蓄冰槽状态递推。% 最小多时间尺度调度 - 单窗口 N 24; % 规划时域小时 dt 1; % 采样间隔小时 Q_cool sdpvar(N, 1); % 主机供冷功率 Q_ice sdpvar(N, 1); % 主机产冰功率 Q_melt sdpvar(N, 1); % 融冰供冷功率 S sdpvar(N1,1);% 蓄冰量状态 constraints [S(1) S0, S 0, S S_max]; for k 1:N % 冷负荷平衡直供 融冰 负荷 constraints [constraints, Q_cool(k) Q_melt(k) load_forecast(k)]; % 出力上下限 constraints [constraints, 0 Q_cool(k) Q_ch_max, 0 Q_ice(k) Q_ch_max]; % 状态递推 constraints [constraints, S(k1) S(k) (eta_ice*Q_ice(k) - Q_melt(k)/eta_melt) * dt]; end逻辑说明constraints先收集最后交给optimize。这里刻意没写模式互斥因为概念验证阶段可以先把主机供冷量和产冰量分开后面再加 0-1 变量。状态变量S的长度是 N1因为末尾时刻的状态也要有边界约束否则蓄冰槽可能在终端被放空或大量剩冷。参数说明S_max是蓄冰槽总蓄冷容量常见取法按第二天峰值冷负荷的 40%60% 配置太小没法转移太多负荷太大夜间制冰时间不够而且蓄冰槽占地面积和成本都会上升。3.3 调用求解器从linprog到整数变量的切换前面模型如果目标函数是线性的且没有引入binvar可以交给linprog求解% 用 YALMIP 求解线性目标 objective sum(price .* P_e) * dt; options sdpsettings(solver, linprog, verbose, 0); result optimize(constraints, objective, options); % 结果回读 Q_cool_opt value(Q_cool); S_opt value(S);逻辑说明optimize的第三个参数是目标函数第二个参数是约束集。求解完成后用value()把 sdpvar 转成数值。调试时把verbose设为 1可以看到求解器收敛信息调参阶段可以设 0避免控制台刷屏。参数说明如果模型里用到了binvar求解器要换成intlinprog或 Gurobi。MATLAB 自带的intlinprog对小规模 24 小时模型没问题但日内滚动层要跑 96 次变量多了以后建议用 Gurobi 或 CPLEX速度差距能达到 10 倍以上。求解器选择参考下表模型复杂度变量类型推荐求解器适用场景纯线性连续变量linprog可行性验证、大时域线性 二次代价连续变量quadprog带软约束惩罚项线性 / 混合整数有 0-1 变量intlinprog、Gurobi、CPLEX模式切换、启停变量非线性 COP 耦合非线性约束fmincon、IPOPT效率曲线非凸时工程上我更推荐把 COP 曲线分段线性化或凸二次化保证每个时域的求解稳定可复现。直接用非线性 COP 加了 fmincon结果可能因初始点不同而漂移源码给别人时很容易被质疑。4. 日前-日内-实时三层调度代码结构、滚动窗口与关键参数4.1 日前计划层24小时小时级计划日前层要做的不是精确执行而是规划边界。它用凌晨的电价预测和冷负荷预测为全天每个小时算出蓄冰槽 SOC 参考曲线、主机直供计划、每小时融冰量上限。判断一个日前计划好不好不是看总电费最低而是看它给日内层留了多少调整空间。如果 SOC 曲线在白天每一段都顶到上限日内一遇到负荷偏高就只能切负荷没有任何余量。常见做法是先解 24 小时开环问题% 日前层求解24小时计划 [Q_cool_day, Q_ice_day, Q_melt_day, S_ref] solve_schedule(N_horizon24); % 保存到结构体供日内调用 plan.S_ref S_ref; plan.Q_melt_ref Q_melt_day;逻辑说明solve_schedule内部就是第 3 章那个模型但约束里要额外加“末期蓄冰量”约束比如S(N1) S0否则求解器会把蓄冰槽在最后一个小时全部放空导致第二天早上没有冷量可用。实际运行中第二天凌晨还要制冰所以末期蓄冰量应该设定为次日零点前的目标值。参数说明N_horizon取 24 是标准做法。如果求解时间吃紧可以缩到 16 小时但覆盖不到次日早高峰容易失衡。日前计划频率低一天只需算一次对求解速度不敏感。4.2 日内滚动层15分钟窗口重算未来4小时日内层是三层里最值得调试的一层。每 15 分钟用最新的短期负荷预测重算未来 4 小时然后只执行第一个控制步。MATLAB 里用循环实现每个时域问题的初值必须跟随当前实际状态。% 日内滚动优化主循环 H 16; % 4小时 / 15分钟 16个控制点 mpc_step 0.25; % 控制周期小时 for t 1:Horizon_points % 取最新的预测切片 load_pred load_forecast_short(t:tH-1); price_slice price(t:tH-1); % 求解当前窗口 [u, S_mpc] solve_mpc_window(load_pred, price_slice, S0, S_actual(t-1)); % 只执行第一个动作 u_apply(t, :) u(:,1); % 用实际负荷更新系统状态 S_actual(t) system_step(S_actual(t-1), u_apply(t,:), load_actual(t)); end逻辑说明solve_mpc_window内部与第 3 章相同但约束里要加“终点牵引项”把窗口末尾状态S_mpc(H1)尽量拉回日前参考S_ref(tH)。不加的话每个窗口都会为了省电费把蓄冰量在窗口末尾压到最低长时间滚动下来蓄冰槽放得过早。S0必须从仿真状态读不能用上一窗口的预测值否则滚动闭环会漂移Horizon_points是全天控制点数15 分钟步长就是 96预测窗口 4 小时覆盖下午电价高峰和负荷拐点计算量适中。如果只做 2 小时日内层会变得近视分不清晚高峰是短时脉冲还是持续负荷。4.3 实时反馈层把蓄冰槽SOC拉回计划日内层以 15 分钟为周期重算但两次计算之间如果有冷负荷脉冲主机得不到命令。实时层就是在这 15 分钟之间用简单规则或快速 MPC 修正主机功率把 SOC 偏差往回调。MATLAB 源码里这一层通常写成回调函数或者 Simulink 里的 MATLAB Function 块。实时修正的典型逻辑是比例控制% 实时反馈比例修正SOC偏差 SOC_error (S_actual - S_ref_current) / S_max; dP Kp * Q_max * SOC_error; % Kp一般取0.1~0.3 u_corrected u_plan dP;逻辑说明SOC 偏差为正说明蓄冰过多实时层可以让主机少产冷、多融冰偏差为负说明冰不够则限制融冰速率并提高主机功率。Kp 过大会造成主机功率抖动过小则 15 分钟内的偏差积累无法消除。离线仿真时可以在每个控制周期内用简单状态方程推演这一小段的内部状态。参数说明SOC_error是无量纲量0~1Q_max是主机额定制冷量所以dP也有功率单位可以直接叠加到主机指令上。Kp 取 0.2 左右比较稳如果系统模型里主机有爬坡限制实时修正量还要再过一层爬坡约束。5. MATLAB源码里最常被改坏的三个模块约束矩阵、SOC更新、模式切换5.1 约束矩阵的行列对应关系不用 YALMIP、直接拼linprog矩阵时最容易犯的错误是行列不对齐。变量按[Q_cool; Q_ice; Q_melt; S(2:end)]排布时Aeq的列数必须等于变量总个数行数等于等式约束数。很多源码里S(k1) - S(k) - eta_ice*Q_ice Q_melt/eta_melt的符号或者索引写错白天蓄冰槽越界优化器还在输出“求解成功”。调试步骤我一般这样走% 校验约束矩阵尺寸 assert(size(Aeq, 2) length(x0), 约束列数不等于决策变量个数); % 用已知可行解检查等式约束残差 resid Aeq * x_init - beq; max(abs(resid)) 1e-6逻辑说明用一个已知可行解去测约束矩阵能在调优化器之前把拼矩阵的错误找出来。比如全部变量为 0同时把初值S(1)设为一个正数看看beq是否对得上。参数说明还有一种隐蔽错误是循环内用Aeq(total,:) ...拼接但total递增写错了位置导致某行被覆盖。检查nnz(Aeq)是否等于预期非零元数量可以快速暴露这种问题。5.2 蓄冰量SOC的离散递推不要写错符号SOC 更新公式是最常被改坏的地方。有人从文献抄代码时会误把制冰加项写成减项而且因为夜间制冰时负荷低、目标函数引导蓄冰即使符号反了优化器也可能靠融冰负数来凑约束结果看起来“结果合理”实际物理上不可能。正确写法是% 正确的 SOC 递推单位换算成 kWh S(k1) S(k) (eta_ice * Q_ice(k) - Q_melt(k) / eta_melt) * dt; % 常见错误1符号反了 % S(k1) S(k) (- eta_ice * Q_ice(k) Q_melt(k) / eta_melt) * dt; % 常见错误2漏了dtQ 是功率S 是能量 % S(k1) S(k) eta_ice * Q_ice(k) - Q_melt(k) / eta_melt;逻辑说明如果dt是 1 小时漏掉dt看不出来控制周期改成 15 分钟后能量误差立刻放大 4 倍。所以dt必须单独设成变量所有模型常量都乘它避免在约束矩阵里硬编码1。参数说明还要确认Q_ice和Q_melt使用同一个功率单位。有些源码用冷吨 RT1 RT 约等于 3.517 kW不换算的话 SOC 会偏离。电机驱动冷机时制冰工况的耗电量要按制冰工况 COP 折算不能沿用制冷工况 COP。5.3 制冷/制冰模式切换不用整数变量会得到“假冰”双工况主机在同一时刻要么制冷要么制冰。如果只写Q_cool Q_ice Q_max求解器很可能让两者各出力 50%得到一个实际不存在的多点运行状态。要避免这个得引入 0-1 模式变量% 模式互斥m1制冷m0制冰 m binvar(N, 1); M Q_ch_max; % 大M数取额定容量稍大一点 constraints [constraints, Q_cool M .* m]; constraints [constraints, Q_ice M .* (1 - m)]; constraints [constraints, Q_cool 0, Q_ice 0];逻辑说明当m1时Q_ice被强制为 0当m0时Q_cool被强制为 0。M要够大但不能过大过大会降低混合整数规划求解速度也不利于数值稳定性。工程上M取1.1 * max(Q_cool_max, Q_ice_max)即可。参数说明加了 0-1 变量后求解器要从linprog切换成intlinprog。如果手拼矩阵要把二进制变量展平进x并在intcon参数里指定索引。只做模式互斥还不够日内滚动层可能在相邻两个窗口让模式来回切。常见的做法是在目标函数里加切换惩罚% 切换次数惩罚sw |diff(m)| sw binvar(N-1, 1); constraints [constraints, sw diff(m), sw -diff(m)]; objective objective lambda_switch * sum(sw);逻辑说明diff(m)在模式变化时为 ±1sw肯定取 1模式不变时为 0。lambda_switch设得足够大求解器会自动减少无意义启停。代价是目标函数里绝对值项被线性化问题变成混合整数线性规划求解时间增加。参数说明lambda_switch的量纲与目标函数一致代表一次启停折算成多少电费成本。我一般先取当天最高电价的一半再根据启停次数曲线微调如果白天模式切换仍然频繁就再放大到一倍。6. 验证调度效果与仿真边界在把源码交给别人之前先做这四件事多时间尺度调度源码跑通只是起点我一般验收时做四个检查能挡住绝大多数“看起来合理结果”的脏数据。6.1 开环回代与闭环滚动对拍用同一组历史数据把日前计划直接一次性求解得到的 SOC 轨迹与日内滚动仿真得到的实际 SOC 轨迹画在一张图上。如果两者差值超过 20%说明日前计划的边界太紧日内层无法跟随。这个检查能直接暴露是否需要加长日内预测时域plot(t, S_day, LineWidth, 1.5); hold on; plot(t, S_actual, --, LineWidth, 1.5); legend(日前计划, 日内实际); ylabel(蓄冰量/kWh);6.2 预测误差灵敏度分析把负荷预测误差从 0% 调整到 ±30%观察总电费波动。如果 10% 预测误差就导致电费上升超过 5%说明调度方案过度依赖预测精度需要增加蓄冰槽余量或加强实时反馈。这个测试在 MATLAB 里就是循环改load_forecast后重复跑调度。6.3 检查约束违反和求解状态每次求解后记录sol.info.problem或exitflag用表格汇总。如果连续多次求解失败优先放宽 SOC 边界或提高惩罚系数。MATLAB 中可以用sol optimize(constraints, objective, options); if sol.info.problem ~ 0 fprintf(求解失败 code%d\n, sol.info.problem); end6.4 用能量守恒指标验证全系统最后检查全天的能量守恒total_chilled sum(Q_cool * dt) sum(Q_melt * dt) S0 - S_end; total_load sum(Q_load * dt); if abs(total_chilled - total_load) 1e-3 warning(能量不平衡: %.4f kWh, total_chilled - total_load); end这四个检查做完数据里的坑基本都能翻出来。剩下的问题通常就出现在预测数据本身和约束边界取值上。本文还有配套的精品资源点击获取