1. 冷热电联供微网为什么需要“冰蓄冷”这个角色先聊聊这个课题的来龙去脉。做综合能源微网优化调度的朋友应该都有体会冷热电联供CCHPCombined Cooling, Heating and Power系统最大的难点不是设备建模有多复杂而是三种能量形式的时序耦合让整个优化问题变得格外棘手。传统的电力调度里电负荷是即时的——发多少用多少发电和用电在同一时刻必须平衡。但到了冷热电联供这个场景事情就不一样了燃气轮机或内燃机发电时会产生大量余热这部分余热可以用来供热也可以通过吸收式制冷机转化为冷量。问题是电负荷高峰、热负荷高峰、冷负荷高峰往往不在同一个时刻出现。比如夏季傍晚电负荷因为空调还在高位但太阳下山后冷负荷已经开始回落余热驱动的吸收式制冷机出力就显得过剩到了正午冷负荷峰值压得人喘不过气电制冷机又得拼命耗电把微网的电平衡逼得很紧。冰蓄冷空调在这里的角色非常特殊——它本质上是一个冷量的储能装置。制冷主机在电价低谷或者余热充足的时候先制冰储存冷量等到冷负荷高峰再融冰释放把原本必须即时平衡的冷负荷变成了可以跨时段调度的灵活资源。这就像在电力系统里加了一台储能电池只不过存的是“冷”而不是“电”。加了冰蓄冷之后微网的调度自由度会明显上升燃气轮机的出力安排、电网购电策略、电制冷机和吸收式制冷机的配比都有了更大的调整空间。这个项目标题里的“多时间尺度优化调度”还有“Matlab代码实现”实际指向的是一套从日前计划到日内滚动修正再到实时调整的完整调度框架。我最初接触这个课题的时候以为所谓多时间尺度就是把三个优化问题分别写一遍然后串起来跑做完了才发现里面的坑远比我预想的多冰蓄冷的状态变量如何在时间尺度之间传递、冷负荷平衡方程里蓄冰放冷怎么折算、不同时间尺度下预测精度差异如何影响约束的松弛策略——这些问题随便拎一个出来都够折腾几天的。这篇文章就把我从模型搭建到代码实现再到排错调试的完整经验做一个整理对正在做CCHP微网调度、冰蓄冷空调、或者需要在Matlab里实现多时间尺度优化的读者应该会比较有参考价值。2. 系统视角下的谋划逻辑2.1 冰蓄冷不是单纯的储能设备而是微网的“冷量调节器”在建模之前需要先把冰蓄冷空调在微网中的角色想清楚。我在第一版模型里把它当成了普通的储能设备——有蓄冰功率、融冰功率、蓄冷量上下限看着跟电池储能模型的套路差不多。但真正跑起来之后发现这种简化忽略了冰蓄冷与冷负荷之间的特殊关系。普通电池储能是向电网发功率的独立电源而冰蓄冷是与电制冷机、吸收式制冷机共同服务冷负荷的协同参与方。冷负荷平衡方程里电制冷机出力、吸收式制冷机出力、融冰放冷量必须加在一起严格等于冷负荷需求制冷工况下。蓄冰工况下方程形式又不同——制冷主机一边供冷一边制冰冷量流向需要分清楚。如果只把冰蓄冷写成一个单输入单输出的储能模块冷量平衡方程会出现明显的物理错误蓄冰时冷量去了哪里融冰时冷量从哪来这些必须用同一套物理逻辑把状态跑通。这是我的第一条经验冰蓄冷模型的核心是明确的工况划分。蓄冰工况Charge、融冰工况Discharge、停机工况Idle之间的切换不能靠逻辑判断生硬实现而要借助状态变量和约束方程的数学关系来反映切换条件。这在下面的模型部分我会详细展开代码思路。2.2 微网能量流电、热、冷三者如何耦合冷热电联供微网的物理结构简要说来是这样一个流程燃气轮机或内燃机燃烧天然气发电输出电功率给微网母线发电产生的余热被余热回收装置捕获一部分直接供给热负荷另一部分进入吸收式制冷机产出冷量冷负荷还可以由电驱动制冷机常规压缩式制冷来满足电制冷机的耗电则进一步增加了电网侧的购电压力冰蓄冷空调作为可转移的冷量来源既可以在谷电时段蓄冰也可以在冷负荷高峰时段融冰供冷。电、热、冷三条能量流在这里形成了交叉耦合电网购电、燃气轮机发电、电制冷耗电、电负荷需求构成电力平衡余热回收、补燃锅炉如果有、热负荷需求构成热力平衡吸收式制冷、电制冷、冰蓄冷放冷、冷负荷需求构成冷量平衡。一个时段内燃气轮机出力调高意味着发电多了、余热也多了供热供冷的潜力同步上升但天然气的燃料成本也在上涨——每个决策变量的变动都会牵动全局。这就是为什么这类优化调度要建数学模型靠人工经验来权衡多能量流的配合关系在维度和规模上完全不现实。Matlab代码实现的时候第一步就是把这种物理耦合关系转化为约束方程组的数学表达之后的优化计算才有依据可言。3. 多时间尺度优化调度从日前计划到日内滚动的递进逻辑3.1 为什么单一时间尺度的调度策略行不通先讲一个反直觉的结论不少论文里把日前优化做得特别精细96个时段逐点优化、约束条件写到极致、目标函数安排得天衣无缝但到了实际运行你会发现这方案就像天气预报的“预测明日气温25摄氏度”——看着精确实际中午到了30度怎么办这里的根本矛盾在于预测信息的时间尺度与不确定性。日前调度依赖的是提前24小时的负荷与可再生能源预测这个时间跨度下预测误差不容小觑。离实际运行还有十几个小时冷负荷预报与实际值差个10%~15%冰蓄冷的蓄冰计划就可能完全偏离最优状态——要么冰蓄多了消化不掉要么蓄少了高峰期能放的全放光了也没顶住冷负荷。所以工业界和学术界不约而同地走向了分段决策的思路不同的决策时间尺度用不同精度的预测信息同时用短时间尺度的决策来修正长时间尺度决策的预见性不足整套调度通过“预测-决策-修正-再预测-再决策”的闭环推进。这就是多时间尺度调度的底层原理。3.2 日前调度层24小时宏观计划日前调度的决策粒度一般是1小时规划未来24小时内各设备的出力基准值。这一层的核心任务包括根据冷热电负荷预测曲线确定燃气轮机的启停计划和出力基线制定冰蓄冷空调的蓄冰/放冰策略的宏观时机比如夜间谷电时段完成主要蓄冰安排与上级电网的购电协议在分时电价机制下确定各时段购电量为后续日内调度提供一个可供追踪的参照轨迹。日前层的优化模型通常是一个混合整数线性规划MILP问题整型变量主要用于刻画设备启停状态和冰蓄冷工况切换。注意冰蓄冷在该层必须有足够时间尺度的“远视”——因为蓄冰是跨多小时的能量搬运只看局部时段必然做出短视决策。3.3 日内滚动调度层15分钟粒度跟踪修正到了日内层时间粒度紧缩到15分钟甚至5分钟。这一层的角色定位是基于最新更新的超短期预测数据对日前计划进行局部修正。日内层的核心逻辑和日前层不同它不是全盘推翻而是沿着日前轨迹做收敛式调整。也就是说日内优化会在目标函数中加入对日前计划偏离的惩罚项——维持主要设备的出力模式稳定小幅调整出力设定值来响应新的预测信息。冰蓄冷在日内层的灵活性会进一步得到体现。如果日内发现冷负荷比预期低不需要马上停蓄冰机组或调燃气轮机出力可以先降低融冰放冷速率把多余的冷量留在冰槽里供晚高峰继续使用这就是储能介质对不确定性天然具备的缓冲能力。3.4 实时调整层应对分钟级波动实时层的任务更纯粹一些以更短的控制周期跟踪日内计划处理实际运行中出现的分钟级偏差。在这一层并需要做大规模全局优化——时间根本不够。工程上常用的做法是控制周期内只解除和微网安全运行密切相关的约束比如母线电压越限、频率偏差做局部快速调整——调整电制冷出力、微调蓄冰或融冰速率、调节吸收式制冷机热水阀门开度让系统实际状态贴近日内计划值。三层之间的信息传递关系可以这样理解日前层给出基准轨迹日内层根据新预测修正轨迹实时层把实际运行拉回轨迹。三层配合的基础就是“预测信息精度随时间尺度缩短而逐步提高”这个基本事实。4. 核心数学模型目标函数与约束条件的搭建细节4.1 目标函数全周期运行成本最小化这个系统的经济性优化目标归结起来是让整个调度周期内微网的总运行成本最小。成本项一般包括这几个部分$$C_{total} \sum_{t} \left( C_{fuel}(t) C_{grid}(t) C_{om}(t) C_{switch}(t) \right)$$其中$C_{fuel}(t)$ 是燃气轮机消耗天然气的燃料成本$C_{grid}(t)$ 是向电网购电/售电费用分时电价下按时段电价结算$C_{om}(t)$ 是设备运行维护成本通常折算为与出力成比例的系数$C_{switch}(t)$ 是设备启停切换的惩罚项可选用于避免机组频繁启停。建模中容易被忽略的一个问题是冰蓄冷系统本身也存在运行维护成本且蓄冰工况和融冰工况下的损耗系数不同需要分开标定。我这里第一版目标函数漏了这个分项结果优化器给出的方案是频繁切换蓄冰/融冰工况来套利——从经济账看毫无意义因为忽略了工况切换带来的机械损耗。判断这类行为的本质还是模型经济项不全面导致的。4.2 能量平衡约束电、热、冷三条主线每条能量流的平衡方程都是硬约束方程形式如下电力平衡微网母线 $$P_{gt}(t) P_{grid}(t) P_{pv}(t) P_{load}(t) P_{ec}(t) P_{ice_charge}(t)$$热力平衡热网母线 $$Q_{recovery}(t) Q_{boiler}(t) Q_{load}(t) Q_{ac_abs}(t)$$冷量平衡冷网母线 $$Q_{ec}(t) Q_{ac_abs}(t) Q_{ice_discharge}(t) Q_{cool_load}(t)$$这三个方程是整个模型的骨架写代码时一定反复检查。我调试时发现的最多的错误就出在这三个方程上。4.3 冰蓄冷模型工况切换与状态转移冰蓄冷是本模型中最有特色也最容易出错的部分。首先蓄冰和融冰被建模为互斥事件$$0 \le P_{ice_charge}(t) \le u_{charge}(t) \cdot P_{charge}^{max}$$ $$0 \le Q_{ice_discharge}(t) \le u_{discharge}(t) \cdot Q_{discharge}^{max}$$ $$u_{charge}(t) u_{discharge}(t) \le 1$$其次冰槽的蓄冷量状态随工况转移 $$S_{ice}(t1) S_{ice}(t) \eta_{charge} \cdot P_{ice_charge}(t) - \frac{Q_{ice_discharge}(t)}{\eta_{discharge}}$$其中 $S_{ice}(t)$ 是 t 时段末的蓄冷量$\eta$ 是蓄/放冷效率。蓄冰和融冰效率互为倒数关系——这个细节也是调试重点我后面会专门细说。还有一个关键细节蓄冰工况下电制冷机的耗电 $P_{ice_charge}(t)$ 其实同时也在产出冷量制冷机一边制冰一边供冷所以蓄冰工况下电制冷机的冷量收益不仅要进冰槽还要同时满足一部分冷负荷。如果你的模型把 $P_{ec}$ 和 $P_{ice_charge}$ 当成两个完全独立的变量冷量平衡就会算错。4.4 设备出力约束与爬坡约束燃气轮机的运行区间一般设置在一个可运行出力上下限之间 $$u_{gt}(t) \cdot P_{gt}^{min} \le P_{gt}(t) \le u_{gt}(t) \cdot P_{gt}^{max}$$热回收装置出力与燃气轮机出力成比例约束因为余热的量取决于发电功率 $$Q_{recovery}(t) \alpha \cdot P_{gt}(t)$$吸收式制冷机的出力也有区间和速率限制电制冷机的爬坡能力相对较强但也不宜写成上下限完全不设防的形式否则容易在优化里出现剧烈的出力甩动。4.5 多时间尺度之间的衔接约束这部分我要特别强调因为很多实现失败都发生在各层模型衔接时状态边界条件设置混乱上。日前调度结束时刻的蓄冰量 $S_{ice}^{day-ahead}(T_{end})$ 是什么状态日内滚动启动时的初始条件是不是把它带进去了冰蓄冷的一个调度周期通常是跨天的白天融冰、夜间蓄冰回满如果各层之间只是简单把决策结果当作初始条件而忽略了冰槽状态本身的慢动态特性蓄冷需要数小时持续放冷才能完成那么日内优化在面对高峰期冷负荷时会错误地认为冰槽还有余量实际却早放空了。推荐的衔接方式是日内层以日前层给出的各时段终态蓄冷量为参考基准加入对该基准值偏离的软约束目标函数中加惩罚项而不是硬性固定相等。这样既保持了计划的连续性又给实时修正留了空间。5. Matlab代码实现从单层优化到多层嵌套的工程落地5.1 代码的模块化设计思路Matlab作为这一场景的开发环境代码的组织方式直接决定了后续调试和维护的效率。我采用的结构如下data/存放负荷曲线、电价曲线、设备参数的Excel或.mat数据文件models/设备数学模型的函数文件燃气轮机模型、电制冷机模型、冰蓄冷模型等optimizers/日前、日内、实时三层的调度求解函数utils/公共工具函数数据读取、结果格式化、可视化输出。建模时最推荐的表达是采用MatlabYALMIP工具箱。YALMIP的符号化建模方式能极大降低建模工作量——你只需要对照数学模型逐行写约束方程YALMIP会负责把优化模型转换成求解器可识别的标准形式。以日前调度层为例核心骨架大致如下% 日前调度层24小时1小时步长MILP求解 T 24; P_gt sdpvar(1, T, full); % 燃气轮机出力 u_gt binvar(1, T, full); % 燃气轮机启停状态 P_grid sdpvar(1, T, full); % 电网购电 P_ec sdpvar(1, T, full); % 电制冷机耗电 Q_ice sdpvar(1, T, full); % 冰蓄冷放冷量 S_ice sdpvar(1, T1, full); % 蓄冷量状态 u_charge binvar(1, T, full); % 蓄冰状态 u_discharge binvar(1, T, full);% 放冰状态 Constraints []; % 设备出力上下限约束、平衡约束、冰蓄冷互斥约束... Constraints [Constraints, ...]; % 目标函数总运行成本最小化 Objective sum(fuel_cost .* P_gt) sum(grid_price .* P_grid) ...; optimize(Constraints, Objective, sdpsettings(solver, gurobi, verbose, 0));5.2 求解器选择为什么我推荐Gurobi或CPLEXMatlab自带的linprog和intlinprog可以作为验证用但真正跑全套模型日用层含整型变量的问题规模达到数千个优化变量、数千条约束时intlinprog的效率与稳定性明显不足。个人实测下来Gurobi在本类问题的求解速度上大概是intlinprog的10~20倍而且对大规模MILP的求解稳定性好得多。求解器相关注意点YALMIP只是建模层真正算还是靠求解器如果没有正版许可可以先使用intlinprog验证正确性再切换到Gurobi/CPLEX做压力测试SDPSettings里建议关闭求解器输出verbose0否则终端会刷出一大片日志影响调试效率。5.3 参数配置与数据读取的预处理Matlab代码的可复现性在很大程度上取决于参数整理是否规范。建议把设备参数集中管理而不是散落在各个脚本里。比较有效的做法是写一个init_params.m返回结构化数组params struct(); params.gt.efficiency 0.35; % 燃气轮机发电效率 params.gt.pmin 100; % 最小出力 kW params.gt.pmax 1000; % 最大出力 kW params.ice.capacity 2000; % 冰槽总容量 kWh params.ice.charge_rate 300; % 最大蓄冰功率 kW params.ice.discharge_rate 500; % 最大放冷功率 kW params.ice.eta_charge 0.90; % 蓄冰效率 params.ice.eta_discharge 0.95; % 放冷效率 params.grid.pbuy [0.3*ones(1,8), ...]; % 分时购电价参数文件独立的直接好处是换数据或换设备型号时不需要改动模型和优化代码只改配置文件即可。这对毕业设计需要做多工况对比比如几种典型日负荷曲线下调度效果对比尤其重要。5.4 从单层到多层的代码改造策略不要一上来就写三层模型那样调试起来会很痛苦。我建议用三层递进的方式先跑通单时段单点优化固定某一时刻的负荷和电价验证设备模型与基本的能量平衡约束正确性。扩展到日前单层优化验证跨时段约束冰蓄冷状态转移、爬坡约束、日内负荷曲线趋势输出日前计划曲线。在日前层基础上加日内滚动把日前层结果作为日内层的基准轨迹输入同时把预测数据替换为日内短时预测系列。最后接实时调整层以日内计划为基准做局部修正输出设备实际运行曲线。每一步的输入数据、输出结果存档要设计好方便随时回溯验证。我自己的经验是在每个阶段的输出目录按结果标识_时间戳保存拥挤的数据空间也比找不到之前的运行结果要好得多。5.5 典型参数与求解规模的参考下表是我在几组合适的数据下记录的求解规模和耗时参考可以作为你验证模型的对照标准调度层时间粒度规划时长决策变量规模求解耗时参考日前层1 h24 h约500个变量含80个整型变量5~20 sGurobi日内层15 min4 h滚动约400个变量含80个整型变量3~10 sGurobi实时层5 min15 min约100个变量少或无常量1 slinprog即可如果你在跑自己数据时发现求解耗时异常优先检查是不是把连续型变量的尺度弄得太离谱比如功率单位是kW却把制冷量写成kW·h数值尺度差异大大会让求解器收敛变慢甚至数值困难。6. 踩坑实录与参数整定冰蓄冷模型中最容易出错的三类问题6.1 蓄冷量与放冷量、耗电量之间的单位折算混乱这个坑我在调试时踩了整整三天。冰蓄冷系统里涉及三个本质不同但数值上容易混淆的量电制冷机的耗电功率电单位kW制冷机制取的冷量冷单位kW或RT制冷吨冰槽储存的蓄冷量冷能单位kWh。这三者在能量转换关系上有明确的系数耗电乘以制冷能效比COP得到制冷量制冷量的一部分进入冰槽储存。如果代码里直接在平衡方程中把三者相加而未做折算物理意义就全乱了。最终的解决办法是建立一套统一的单位转换约定并写进参数文件里作为一个显式的转换参数params.conversion.cop_ec 3.5; % 电制冷机能效比 COP params.conversion.kwh_rt 3.517; % 1 USRT 3.517 kW然后在冷量平衡方程中严格区分电制冷机产冷量 电制冷机耗电量 × COP冰蓄冷的融冰放冷量直接用冷量功率单位kW参与平衡蓄冷量状态的单位固定为kWh。6.2 蓄冰/融冰互斥约束在多时间尺度下的隐式冲突按照日前层写互斥约束日内层也写实时层也写乍看很合理。但在某次调试中我发现日内层的解在某个时刻出现了“先蓄冰再融冰”这种毫无经济意义的现象——不是代码错误而是互斥约束写得太强了。原因是这样日内层滚动优化时时间窗从当前时刻往前推而日前层给出的冰槽蓄冷量基准值可能已经预告了未来的蓄/放冰计划。日内层的优化变量在窗口起点到未来某时刻可能既要满足蓄冰状态又要满足基准轨迹约束两个约束打架时求解器在边界处被迫做了状态反复跳变的妥协。解决思路是放宽日内层中冰蓄冷工况切换约束允许状态在相邻时段发生切换但加上切换次数上限约束total variation penalty或者切换成本让求解器在冲突时刻选择最经济的做法而不是给出无意义的来回切换。我在代码里增加了对一天内蓄/放工况切换总次数的限制参数实测效果明显改善。6.3 可复现性的“莫测”故障同一份代码换了机器结果不同还有一次比较困扰的经验是关于结果可复现性的同样的代码在A电脑跑出来的结果总是和B电脑有细微出入。排查了很久最终发现是YALMIP生成模型时内部随机数的种子问题求解器若使用启发式初始解需要随机种子。解决办法很直接在调用求解器之前显式设置随机数种子并在SDPSettings里指定求解算法的详细参数rng(20240501); % 固定随机种子 options sdpsettings(solver,gurobi,verbose,0, ... gurobi.MIPGap, 0.01, ... % 设置MIP间隙 gurobi.TimeLimit, 300, ... % 设置求解时间上限 gurobi.Seed, 42); % Gurobi自己的随机种子这种细节在发布代码供后续测试复现时非常重要。如果不固定种子别人下载你的代码跑出来的结果是无法验证的这也是学术代码分享时必须注意的问题。6.4 关于预测数据千万不要在模型中对预测误差掉以轻心最后一个值得讲的坑是关于预测数据的输入方式。很多参考代码直接用一个负荷曲线文件喂给日前层和日内层而且假定日内层的预测和实际完全相同——这在概念上就错了完全失去了多时间尺度的意义。多时间尺度调度的核心价值就在于短时间尺度能获取更准确的预测信息所以日内层的输入应该是带误差、需要更新的短时预测曲线而不是直接使用实际值。最直接的做法是用蒙特卡洛方法生成若干场景下的预测误差分布或者至少用“实测值叠加正态扰动”的方式模拟短期预测数据的偏差测试调度策略的鲁棒性。如果所有层都用同一份完美预测数据跑多时间尺度的优势根本显现不出来——从那种结果里你会得出一个错误结论多时间尺度模型和普通单层模型一样。这一点在论文对比实验设计上尤其要留意。7. 从复现到扩展三组实用改进方向7.1 增加储电/蓄热模块后的协调控制如果你的研究需要把冰蓄冷扩展到更完整的储充体系可以考虑加入电储能电池和蓄热罐。混合储能系统的高价值在于不同储能介质的物理特性差异非常大电池响应快、循环次数有限、适合短时功率调节蓄冰槽响应慢、容量大、适合长时间冷量转移蓄热罐介于两者之间。对各种储能介质在时间尺度、功率密度、容量密度上的差异建立差异化的调度模型比随便组装多个储能方程要有意义得多。7.2 用算法对比体现多时间尺度的优势实现多时间尺度调度后值得对比的是简化模型——例如同样一套参数跑单层日前优化即多时间尺度的三层框架压缩成一层加上实际误差代入结果与多时间尺度的全套模型作经济性对比这样能直观量化多时间尺度带来的收益比例。这类对比不仅说服力强而且能帮助你从底层理解这套框架的价值来源对论文润色有很大帮助。7.3 从MILP到鲁棒/随机优化的进阶方向考虑进一步深入的话可以将求解器能力和Matlab环境扩展到鲁棒优化或随机优化框架。冰蓄冷系统的核心不确定性主要来自冷负荷预测误差而在多时间尺度框架下日内层会不断获得新的预测信息来修正恰好可以用滚动时域控制理论和随机模型预测控制Stochastic MPC来做更严谨的表述。如果时间宽裕这是我推荐的深挖方向——它对工程直觉的培养和论文水平的提升都很有好处。整套项目从模型建立到代码实现到最后的调试排错回头再看最核心的收获其实不是“会调Gurobi”“会写YALMIP约束”而是理解了多时间尺度这个框架的物理意义它是在跟不确定性打交道。不同精度的预测信息有不同的时效窗口调度决策也要跟着这个时效窗口去分层而不是试图用一个完美模型把所有问题一次性解决。冰蓄冷空调在这个框架里起到了很妙的“缓冲器”作用——它让冷负荷不再是分秒必争的硬需求而变成了可以按经济性安排的灵活资源。如果你正在做类似课题建议先把冷负荷预测曲线摸透把冰蓄冷的物理特性吃透再动手写Matlab代码这样后面每步推起来会顺很多。 冰蓄冷与多时间尺度优化调度：CCHP微网Matlab实现
1. 冷热电联供微网为什么需要“冰蓄冷”这个角色先聊聊这个课题的来龙去脉。做综合能源微网优化调度的朋友应该都有体会冷热电联供CCHPCombined Cooling, Heating and Power系统最大的难点不是设备建模有多复杂而是三种能量形式的时序耦合让整个优化问题变得格外棘手。传统的电力调度里电负荷是即时的——发多少用多少发电和用电在同一时刻必须平衡。但到了冷热电联供这个场景事情就不一样了燃气轮机或内燃机发电时会产生大量余热这部分余热可以用来供热也可以通过吸收式制冷机转化为冷量。问题是电负荷高峰、热负荷高峰、冷负荷高峰往往不在同一个时刻出现。比如夏季傍晚电负荷因为空调还在高位但太阳下山后冷负荷已经开始回落余热驱动的吸收式制冷机出力就显得过剩到了正午冷负荷峰值压得人喘不过气电制冷机又得拼命耗电把微网的电平衡逼得很紧。冰蓄冷空调在这里的角色非常特殊——它本质上是一个冷量的储能装置。制冷主机在电价低谷或者余热充足的时候先制冰储存冷量等到冷负荷高峰再融冰释放把原本必须即时平衡的冷负荷变成了可以跨时段调度的灵活资源。这就像在电力系统里加了一台储能电池只不过存的是“冷”而不是“电”。加了冰蓄冷之后微网的调度自由度会明显上升燃气轮机的出力安排、电网购电策略、电制冷机和吸收式制冷机的配比都有了更大的调整空间。这个项目标题里的“多时间尺度优化调度”还有“Matlab代码实现”实际指向的是一套从日前计划到日内滚动修正再到实时调整的完整调度框架。我最初接触这个课题的时候以为所谓多时间尺度就是把三个优化问题分别写一遍然后串起来跑做完了才发现里面的坑远比我预想的多冰蓄冷的状态变量如何在时间尺度之间传递、冷负荷平衡方程里蓄冰放冷怎么折算、不同时间尺度下预测精度差异如何影响约束的松弛策略——这些问题随便拎一个出来都够折腾几天的。这篇文章就把我从模型搭建到代码实现再到排错调试的完整经验做一个整理对正在做CCHP微网调度、冰蓄冷空调、或者需要在Matlab里实现多时间尺度优化的读者应该会比较有参考价值。2. 系统视角下的谋划逻辑2.1 冰蓄冷不是单纯的储能设备而是微网的“冷量调节器”在建模之前需要先把冰蓄冷空调在微网中的角色想清楚。我在第一版模型里把它当成了普通的储能设备——有蓄冰功率、融冰功率、蓄冷量上下限看着跟电池储能模型的套路差不多。但真正跑起来之后发现这种简化忽略了冰蓄冷与冷负荷之间的特殊关系。普通电池储能是向电网发功率的独立电源而冰蓄冷是与电制冷机、吸收式制冷机共同服务冷负荷的协同参与方。冷负荷平衡方程里电制冷机出力、吸收式制冷机出力、融冰放冷量必须加在一起严格等于冷负荷需求制冷工况下。蓄冰工况下方程形式又不同——制冷主机一边供冷一边制冰冷量流向需要分清楚。如果只把冰蓄冷写成一个单输入单输出的储能模块冷量平衡方程会出现明显的物理错误蓄冰时冷量去了哪里融冰时冷量从哪来这些必须用同一套物理逻辑把状态跑通。这是我的第一条经验冰蓄冷模型的核心是明确的工况划分。蓄冰工况Charge、融冰工况Discharge、停机工况Idle之间的切换不能靠逻辑判断生硬实现而要借助状态变量和约束方程的数学关系来反映切换条件。这在下面的模型部分我会详细展开代码思路。2.2 微网能量流电、热、冷三者如何耦合冷热电联供微网的物理结构简要说来是这样一个流程燃气轮机或内燃机燃烧天然气发电输出电功率给微网母线发电产生的余热被余热回收装置捕获一部分直接供给热负荷另一部分进入吸收式制冷机产出冷量冷负荷还可以由电驱动制冷机常规压缩式制冷来满足电制冷机的耗电则进一步增加了电网侧的购电压力冰蓄冷空调作为可转移的冷量来源既可以在谷电时段蓄冰也可以在冷负荷高峰时段融冰供冷。电、热、冷三条能量流在这里形成了交叉耦合电网购电、燃气轮机发电、电制冷耗电、电负荷需求构成电力平衡余热回收、补燃锅炉如果有、热负荷需求构成热力平衡吸收式制冷、电制冷、冰蓄冷放冷、冷负荷需求构成冷量平衡。一个时段内燃气轮机出力调高意味着发电多了、余热也多了供热供冷的潜力同步上升但天然气的燃料成本也在上涨——每个决策变量的变动都会牵动全局。这就是为什么这类优化调度要建数学模型靠人工经验来权衡多能量流的配合关系在维度和规模上完全不现实。Matlab代码实现的时候第一步就是把这种物理耦合关系转化为约束方程组的数学表达之后的优化计算才有依据可言。3. 多时间尺度优化调度从日前计划到日内滚动的递进逻辑3.1 为什么单一时间尺度的调度策略行不通先讲一个反直觉的结论不少论文里把日前优化做得特别精细96个时段逐点优化、约束条件写到极致、目标函数安排得天衣无缝但到了实际运行你会发现这方案就像天气预报的“预测明日气温25摄氏度”——看着精确实际中午到了30度怎么办这里的根本矛盾在于预测信息的时间尺度与不确定性。日前调度依赖的是提前24小时的负荷与可再生能源预测这个时间跨度下预测误差不容小觑。离实际运行还有十几个小时冷负荷预报与实际值差个10%~15%冰蓄冷的蓄冰计划就可能完全偏离最优状态——要么冰蓄多了消化不掉要么蓄少了高峰期能放的全放光了也没顶住冷负荷。所以工业界和学术界不约而同地走向了分段决策的思路不同的决策时间尺度用不同精度的预测信息同时用短时间尺度的决策来修正长时间尺度决策的预见性不足整套调度通过“预测-决策-修正-再预测-再决策”的闭环推进。这就是多时间尺度调度的底层原理。3.2 日前调度层24小时宏观计划日前调度的决策粒度一般是1小时规划未来24小时内各设备的出力基准值。这一层的核心任务包括根据冷热电负荷预测曲线确定燃气轮机的启停计划和出力基线制定冰蓄冷空调的蓄冰/放冰策略的宏观时机比如夜间谷电时段完成主要蓄冰安排与上级电网的购电协议在分时电价机制下确定各时段购电量为后续日内调度提供一个可供追踪的参照轨迹。日前层的优化模型通常是一个混合整数线性规划MILP问题整型变量主要用于刻画设备启停状态和冰蓄冷工况切换。注意冰蓄冷在该层必须有足够时间尺度的“远视”——因为蓄冰是跨多小时的能量搬运只看局部时段必然做出短视决策。3.3 日内滚动调度层15分钟粒度跟踪修正到了日内层时间粒度紧缩到15分钟甚至5分钟。这一层的角色定位是基于最新更新的超短期预测数据对日前计划进行局部修正。日内层的核心逻辑和日前层不同它不是全盘推翻而是沿着日前轨迹做收敛式调整。也就是说日内优化会在目标函数中加入对日前计划偏离的惩罚项——维持主要设备的出力模式稳定小幅调整出力设定值来响应新的预测信息。冰蓄冷在日内层的灵活性会进一步得到体现。如果日内发现冷负荷比预期低不需要马上停蓄冰机组或调燃气轮机出力可以先降低融冰放冷速率把多余的冷量留在冰槽里供晚高峰继续使用这就是储能介质对不确定性天然具备的缓冲能力。3.4 实时调整层应对分钟级波动实时层的任务更纯粹一些以更短的控制周期跟踪日内计划处理实际运行中出现的分钟级偏差。在这一层并需要做大规模全局优化——时间根本不够。工程上常用的做法是控制周期内只解除和微网安全运行密切相关的约束比如母线电压越限、频率偏差做局部快速调整——调整电制冷出力、微调蓄冰或融冰速率、调节吸收式制冷机热水阀门开度让系统实际状态贴近日内计划值。三层之间的信息传递关系可以这样理解日前层给出基准轨迹日内层根据新预测修正轨迹实时层把实际运行拉回轨迹。三层配合的基础就是“预测信息精度随时间尺度缩短而逐步提高”这个基本事实。4. 核心数学模型目标函数与约束条件的搭建细节4.1 目标函数全周期运行成本最小化这个系统的经济性优化目标归结起来是让整个调度周期内微网的总运行成本最小。成本项一般包括这几个部分$$C_{total} \sum_{t} \left( C_{fuel}(t) C_{grid}(t) C_{om}(t) C_{switch}(t) \right)$$其中$C_{fuel}(t)$ 是燃气轮机消耗天然气的燃料成本$C_{grid}(t)$ 是向电网购电/售电费用分时电价下按时段电价结算$C_{om}(t)$ 是设备运行维护成本通常折算为与出力成比例的系数$C_{switch}(t)$ 是设备启停切换的惩罚项可选用于避免机组频繁启停。建模中容易被忽略的一个问题是冰蓄冷系统本身也存在运行维护成本且蓄冰工况和融冰工况下的损耗系数不同需要分开标定。我这里第一版目标函数漏了这个分项结果优化器给出的方案是频繁切换蓄冰/融冰工况来套利——从经济账看毫无意义因为忽略了工况切换带来的机械损耗。判断这类行为的本质还是模型经济项不全面导致的。4.2 能量平衡约束电、热、冷三条主线每条能量流的平衡方程都是硬约束方程形式如下电力平衡微网母线 $$P_{gt}(t) P_{grid}(t) P_{pv}(t) P_{load}(t) P_{ec}(t) P_{ice_charge}(t)$$热力平衡热网母线 $$Q_{recovery}(t) Q_{boiler}(t) Q_{load}(t) Q_{ac_abs}(t)$$冷量平衡冷网母线 $$Q_{ec}(t) Q_{ac_abs}(t) Q_{ice_discharge}(t) Q_{cool_load}(t)$$这三个方程是整个模型的骨架写代码时一定反复检查。我调试时发现的最多的错误就出在这三个方程上。4.3 冰蓄冷模型工况切换与状态转移冰蓄冷是本模型中最有特色也最容易出错的部分。首先蓄冰和融冰被建模为互斥事件$$0 \le P_{ice_charge}(t) \le u_{charge}(t) \cdot P_{charge}^{max}$$ $$0 \le Q_{ice_discharge}(t) \le u_{discharge}(t) \cdot Q_{discharge}^{max}$$ $$u_{charge}(t) u_{discharge}(t) \le 1$$其次冰槽的蓄冷量状态随工况转移 $$S_{ice}(t1) S_{ice}(t) \eta_{charge} \cdot P_{ice_charge}(t) - \frac{Q_{ice_discharge}(t)}{\eta_{discharge}}$$其中 $S_{ice}(t)$ 是 t 时段末的蓄冷量$\eta$ 是蓄/放冷效率。蓄冰和融冰效率互为倒数关系——这个细节也是调试重点我后面会专门细说。还有一个关键细节蓄冰工况下电制冷机的耗电 $P_{ice_charge}(t)$ 其实同时也在产出冷量制冷机一边制冰一边供冷所以蓄冰工况下电制冷机的冷量收益不仅要进冰槽还要同时满足一部分冷负荷。如果你的模型把 $P_{ec}$ 和 $P_{ice_charge}$ 当成两个完全独立的变量冷量平衡就会算错。4.4 设备出力约束与爬坡约束燃气轮机的运行区间一般设置在一个可运行出力上下限之间 $$u_{gt}(t) \cdot P_{gt}^{min} \le P_{gt}(t) \le u_{gt}(t) \cdot P_{gt}^{max}$$热回收装置出力与燃气轮机出力成比例约束因为余热的量取决于发电功率 $$Q_{recovery}(t) \alpha \cdot P_{gt}(t)$$吸收式制冷机的出力也有区间和速率限制电制冷机的爬坡能力相对较强但也不宜写成上下限完全不设防的形式否则容易在优化里出现剧烈的出力甩动。4.5 多时间尺度之间的衔接约束这部分我要特别强调因为很多实现失败都发生在各层模型衔接时状态边界条件设置混乱上。日前调度结束时刻的蓄冰量 $S_{ice}^{day-ahead}(T_{end})$ 是什么状态日内滚动启动时的初始条件是不是把它带进去了冰蓄冷的一个调度周期通常是跨天的白天融冰、夜间蓄冰回满如果各层之间只是简单把决策结果当作初始条件而忽略了冰槽状态本身的慢动态特性蓄冷需要数小时持续放冷才能完成那么日内优化在面对高峰期冷负荷时会错误地认为冰槽还有余量实际却早放空了。推荐的衔接方式是日内层以日前层给出的各时段终态蓄冷量为参考基准加入对该基准值偏离的软约束目标函数中加惩罚项而不是硬性固定相等。这样既保持了计划的连续性又给实时修正留了空间。5. Matlab代码实现从单层优化到多层嵌套的工程落地5.1 代码的模块化设计思路Matlab作为这一场景的开发环境代码的组织方式直接决定了后续调试和维护的效率。我采用的结构如下data/存放负荷曲线、电价曲线、设备参数的Excel或.mat数据文件models/设备数学模型的函数文件燃气轮机模型、电制冷机模型、冰蓄冷模型等optimizers/日前、日内、实时三层的调度求解函数utils/公共工具函数数据读取、结果格式化、可视化输出。建模时最推荐的表达是采用MatlabYALMIP工具箱。YALMIP的符号化建模方式能极大降低建模工作量——你只需要对照数学模型逐行写约束方程YALMIP会负责把优化模型转换成求解器可识别的标准形式。以日前调度层为例核心骨架大致如下% 日前调度层24小时1小时步长MILP求解 T 24; P_gt sdpvar(1, T, full); % 燃气轮机出力 u_gt binvar(1, T, full); % 燃气轮机启停状态 P_grid sdpvar(1, T, full); % 电网购电 P_ec sdpvar(1, T, full); % 电制冷机耗电 Q_ice sdpvar(1, T, full); % 冰蓄冷放冷量 S_ice sdpvar(1, T1, full); % 蓄冷量状态 u_charge binvar(1, T, full); % 蓄冰状态 u_discharge binvar(1, T, full);% 放冰状态 Constraints []; % 设备出力上下限约束、平衡约束、冰蓄冷互斥约束... Constraints [Constraints, ...]; % 目标函数总运行成本最小化 Objective sum(fuel_cost .* P_gt) sum(grid_price .* P_grid) ...; optimize(Constraints, Objective, sdpsettings(solver, gurobi, verbose, 0));5.2 求解器选择为什么我推荐Gurobi或CPLEXMatlab自带的linprog和intlinprog可以作为验证用但真正跑全套模型日用层含整型变量的问题规模达到数千个优化变量、数千条约束时intlinprog的效率与稳定性明显不足。个人实测下来Gurobi在本类问题的求解速度上大概是intlinprog的10~20倍而且对大规模MILP的求解稳定性好得多。求解器相关注意点YALMIP只是建模层真正算还是靠求解器如果没有正版许可可以先使用intlinprog验证正确性再切换到Gurobi/CPLEX做压力测试SDPSettings里建议关闭求解器输出verbose0否则终端会刷出一大片日志影响调试效率。5.3 参数配置与数据读取的预处理Matlab代码的可复现性在很大程度上取决于参数整理是否规范。建议把设备参数集中管理而不是散落在各个脚本里。比较有效的做法是写一个init_params.m返回结构化数组params struct(); params.gt.efficiency 0.35; % 燃气轮机发电效率 params.gt.pmin 100; % 最小出力 kW params.gt.pmax 1000; % 最大出力 kW params.ice.capacity 2000; % 冰槽总容量 kWh params.ice.charge_rate 300; % 最大蓄冰功率 kW params.ice.discharge_rate 500; % 最大放冷功率 kW params.ice.eta_charge 0.90; % 蓄冰效率 params.ice.eta_discharge 0.95; % 放冷效率 params.grid.pbuy [0.3*ones(1,8), ...]; % 分时购电价参数文件独立的直接好处是换数据或换设备型号时不需要改动模型和优化代码只改配置文件即可。这对毕业设计需要做多工况对比比如几种典型日负荷曲线下调度效果对比尤其重要。5.4 从单层到多层的代码改造策略不要一上来就写三层模型那样调试起来会很痛苦。我建议用三层递进的方式先跑通单时段单点优化固定某一时刻的负荷和电价验证设备模型与基本的能量平衡约束正确性。扩展到日前单层优化验证跨时段约束冰蓄冷状态转移、爬坡约束、日内负荷曲线趋势输出日前计划曲线。在日前层基础上加日内滚动把日前层结果作为日内层的基准轨迹输入同时把预测数据替换为日内短时预测系列。最后接实时调整层以日内计划为基准做局部修正输出设备实际运行曲线。每一步的输入数据、输出结果存档要设计好方便随时回溯验证。我自己的经验是在每个阶段的输出目录按结果标识_时间戳保存拥挤的数据空间也比找不到之前的运行结果要好得多。5.5 典型参数与求解规模的参考下表是我在几组合适的数据下记录的求解规模和耗时参考可以作为你验证模型的对照标准调度层时间粒度规划时长决策变量规模求解耗时参考日前层1 h24 h约500个变量含80个整型变量5~20 sGurobi日内层15 min4 h滚动约400个变量含80个整型变量3~10 sGurobi实时层5 min15 min约100个变量少或无常量1 slinprog即可如果你在跑自己数据时发现求解耗时异常优先检查是不是把连续型变量的尺度弄得太离谱比如功率单位是kW却把制冷量写成kW·h数值尺度差异大大会让求解器收敛变慢甚至数值困难。6. 踩坑实录与参数整定冰蓄冷模型中最容易出错的三类问题6.1 蓄冷量与放冷量、耗电量之间的单位折算混乱这个坑我在调试时踩了整整三天。冰蓄冷系统里涉及三个本质不同但数值上容易混淆的量电制冷机的耗电功率电单位kW制冷机制取的冷量冷单位kW或RT制冷吨冰槽储存的蓄冷量冷能单位kWh。这三者在能量转换关系上有明确的系数耗电乘以制冷能效比COP得到制冷量制冷量的一部分进入冰槽储存。如果代码里直接在平衡方程中把三者相加而未做折算物理意义就全乱了。最终的解决办法是建立一套统一的单位转换约定并写进参数文件里作为一个显式的转换参数params.conversion.cop_ec 3.5; % 电制冷机能效比 COP params.conversion.kwh_rt 3.517; % 1 USRT 3.517 kW然后在冷量平衡方程中严格区分电制冷机产冷量 电制冷机耗电量 × COP冰蓄冷的融冰放冷量直接用冷量功率单位kW参与平衡蓄冷量状态的单位固定为kWh。6.2 蓄冰/融冰互斥约束在多时间尺度下的隐式冲突按照日前层写互斥约束日内层也写实时层也写乍看很合理。但在某次调试中我发现日内层的解在某个时刻出现了“先蓄冰再融冰”这种毫无经济意义的现象——不是代码错误而是互斥约束写得太强了。原因是这样日内层滚动优化时时间窗从当前时刻往前推而日前层给出的冰槽蓄冷量基准值可能已经预告了未来的蓄/放冰计划。日内层的优化变量在窗口起点到未来某时刻可能既要满足蓄冰状态又要满足基准轨迹约束两个约束打架时求解器在边界处被迫做了状态反复跳变的妥协。解决思路是放宽日内层中冰蓄冷工况切换约束允许状态在相邻时段发生切换但加上切换次数上限约束total variation penalty或者切换成本让求解器在冲突时刻选择最经济的做法而不是给出无意义的来回切换。我在代码里增加了对一天内蓄/放工况切换总次数的限制参数实测效果明显改善。6.3 可复现性的“莫测”故障同一份代码换了机器结果不同还有一次比较困扰的经验是关于结果可复现性的同样的代码在A电脑跑出来的结果总是和B电脑有细微出入。排查了很久最终发现是YALMIP生成模型时内部随机数的种子问题求解器若使用启发式初始解需要随机种子。解决办法很直接在调用求解器之前显式设置随机数种子并在SDPSettings里指定求解算法的详细参数rng(20240501); % 固定随机种子 options sdpsettings(solver,gurobi,verbose,0, ... gurobi.MIPGap, 0.01, ... % 设置MIP间隙 gurobi.TimeLimit, 300, ... % 设置求解时间上限 gurobi.Seed, 42); % Gurobi自己的随机种子这种细节在发布代码供后续测试复现时非常重要。如果不固定种子别人下载你的代码跑出来的结果是无法验证的这也是学术代码分享时必须注意的问题。6.4 关于预测数据千万不要在模型中对预测误差掉以轻心最后一个值得讲的坑是关于预测数据的输入方式。很多参考代码直接用一个负荷曲线文件喂给日前层和日内层而且假定日内层的预测和实际完全相同——这在概念上就错了完全失去了多时间尺度的意义。多时间尺度调度的核心价值就在于短时间尺度能获取更准确的预测信息所以日内层的输入应该是带误差、需要更新的短时预测曲线而不是直接使用实际值。最直接的做法是用蒙特卡洛方法生成若干场景下的预测误差分布或者至少用“实测值叠加正态扰动”的方式模拟短期预测数据的偏差测试调度策略的鲁棒性。如果所有层都用同一份完美预测数据跑多时间尺度的优势根本显现不出来——从那种结果里你会得出一个错误结论多时间尺度模型和普通单层模型一样。这一点在论文对比实验设计上尤其要留意。7. 从复现到扩展三组实用改进方向7.1 增加储电/蓄热模块后的协调控制如果你的研究需要把冰蓄冷扩展到更完整的储充体系可以考虑加入电储能电池和蓄热罐。混合储能系统的高价值在于不同储能介质的物理特性差异非常大电池响应快、循环次数有限、适合短时功率调节蓄冰槽响应慢、容量大、适合长时间冷量转移蓄热罐介于两者之间。对各种储能介质在时间尺度、功率密度、容量密度上的差异建立差异化的调度模型比随便组装多个储能方程要有意义得多。7.2 用算法对比体现多时间尺度的优势实现多时间尺度调度后值得对比的是简化模型——例如同样一套参数跑单层日前优化即多时间尺度的三层框架压缩成一层加上实际误差代入结果与多时间尺度的全套模型作经济性对比这样能直观量化多时间尺度带来的收益比例。这类对比不仅说服力强而且能帮助你从底层理解这套框架的价值来源对论文润色有很大帮助。7.3 从MILP到鲁棒/随机优化的进阶方向考虑进一步深入的话可以将求解器能力和Matlab环境扩展到鲁棒优化或随机优化框架。冰蓄冷系统的核心不确定性主要来自冷负荷预测误差而在多时间尺度框架下日内层会不断获得新的预测信息来修正恰好可以用滚动时域控制理论和随机模型预测控制Stochastic MPC来做更严谨的表述。如果时间宽裕这是我推荐的深挖方向——它对工程直觉的培养和论文水平的提升都很有好处。整套项目从模型建立到代码实现到最后的调试排错回头再看最核心的收获其实不是“会调Gurobi”“会写YALMIP约束”而是理解了多时间尺度这个框架的物理意义它是在跟不确定性打交道。不同精度的预测信息有不同的时效窗口调度决策也要跟着这个时效窗口去分层而不是试图用一个完美模型把所有问题一次性解决。冰蓄冷空调在这个框架里起到了很妙的“缓冲器”作用——它让冷负荷不再是分秒必争的硬需求而变成了可以按经济性安排的灵活资源。如果你正在做类似课题建议先把冷负荷预测曲线摸透把冰蓄冷的物理特性吃透再动手写Matlab代码这样后面每步推起来会顺很多。