最近在做一个含冰蓄冷装置的冷电联供型微网经济优化运行项目要用MATLAB把优化调度模型完整跑通。最开始我以为是常规的冷热电联供模型改改参数就行真正动手才发现冰蓄冷装置的引入让整个优化问题的复杂度抬高了一大截——不仅要同时平衡电负荷和冷负荷两条能量流还得处理蓄能设备带来的时序耦合约束。本文就把我用MATLAB从建模到求解的完整过程梳理一遍包括核心代码思路、典型日结果数据和调试中踩过的坑给正在做微网经济调度或者综合能源系统优化的朋友做个参考。1. 冷电联供微网为什么要配冰蓄冷装置1.1 系统结构冷、电、蓄三条能量流怎么耦合我接的这个项目场景是一个园区级微网主要负荷是冷负荷和电负荷。供能侧的核心是一台燃气轮机发出来的电优先满足自身负荷不够的部分从电网购入。燃气轮机的余热被余热锅炉回收驱动溴化锂吸收式制冷机供冷同时还有一台电制冷机作为冷负荷的补充和调节手段。关键就在于系统里加了一台冰蓄冷装置——夜间电价低谷时段让电制冷机多出力把冷量以冰的形式存起来白天电价高峰时再融冰放冷。这个系统的能量流比单纯的电热联供要复杂。因为冷和电两条能量流在电制冷机上耦合电制冷机既可以用电直接制冷供冷负荷也可以把制冷量灌进蓄冰槽。这样一来优化调度的问题就不只是设备开多少出力而是要同步决定燃气轮机发多少电、电网买多少电、电制冷机产生的冷量分多少给负荷、分多少给冰槽、冰槽什么时候蓄、什么时候放。每一步决策都会影响下一步的可行域这也是这类模型最需要花心思的地方。1.2 冰蓄冷的移峰填谷经济账冰蓄冷的核心是利用相变潜热储能1kg水变成冰要释放约334kJ的热量这个能量密度比水的显热蓄冷大得多所以在同样的蓄冷量需求下冰蓄冷槽的体积可以做得很小。工程上选它而不是大水罐核心就是看中这一点。经济账怎么算假设当地分时电价是峰时0.95元/kWh、谷时0.35元/kWh电制冷机的综合COP按3.5算。夜间花1kWh谷电制冰理论上能产出约3.5kWh的冷量即使考虑制冰工况下COP会降到3.0左右、蓄冰与融冰环节的综合效率按0.75算最后也能从冰槽取出约2.2kWh的冷量成本只有0.35元。这2.2kWh的冷量如果白天直接用高峰电价电制冷来补需要耗电约2.2/3.50.63kWh电费约0.6元。一来一回每kWh谷电用于制冰能省下约0.25元。一天蓄个几百kWh日运行成本降一两百块钱是现实的。当然实际优化中的蓄冰策略不是简单按峰谷电价差套利就完事。燃气轮机发电的同时会产生余热余热驱动的吸收式制冷边际成本很低天然气价格、购电价格、机组效率、冷负荷曲线之间是耦合的。所以最后到底什么时候蓄冰、蓄多少冰得靠优化模型去算而不是拍脑袋定。这也是本文要讲的核心。1.3 这类课题的典型应用场景带冰蓄冷的冷电联供微网最典型的应用场景就是白天冷负荷大、峰谷电价差距明显的园区——大型商业综合体、数据中心、食品冷链加工厂都属于这一类。这些地方空调或工艺制冷负荷占大头而且负荷曲线和电价曲线高度重叠白天越热越要制冷偏偏又赶上电价高峰。冰蓄冷正好能把制冷需求和电力消耗在时间上解耦。对做课题或者前期方案论证的人来说经济优化运行的意义在于设备都选型好了、参数都知道了怎么安排24小时内的出力计划让运行成本最低。这一步算得准后面无论是写可研报告还是做能量管理系统EMS的调度策略都有直接价值。2. 经济优化运行的数学模型怎么搭2.1 目标函数运行成本到底包含哪几项优化目标很直接一个调度周期通常取24小时内的总运行成本最小。我把成本拆成四项燃气轮机的燃料成本由发电功率和效率反推燃料消耗量再乘天然气价格。如果用线性化效率就是c_gas × (P_GT/η_GT) × Δt。电网购电成本按分时电价c_buy(t) × P_buy(t) × Δt计算。如果允许余电上网再减去售电收益。设备运行维护成本燃气轮机、电制冷机、蓄冰槽都按出力或启停次数计提系数一般是几厘到几分钱每kWh。碳排放成本可选如果课题要求考虑碳排放可以在目标函数里加一个CO2排放惩罚项。这部分对结果有一定影响但本文没把它放进主模型免得把问题绕大。需要特别提醒的是成本函数里所有功率×时间才是能量千万别忘了乘调度步长Δt。很多人第一次跑模型结果莫名其妙查到最后就是单位问题。2.2 约束条件的核心功率平衡、蓄能连续性、设备界限模型的约束我大致分成三类。第一类是瞬时平衡。每个时刻都要满足电功率平衡P_GT(t) P_buy(t) P_E_load(t) P_EC(t)其中P_EC是电制冷机耗电。冷功率平衡Q_AC(t) Q_direct(t) Q_dis(t) Q_C_load(t)其中Q_AC是吸收式制冷量Q_direct是电制冷机直接供冷量Q_dis是融冰放冷量。第二类是这类系统最关键的部分——蓄能连续性。蓄冰槽的状态方程是S(t1) S(t) η_ch × Q_ch(t) - Q_dis(t) / η_dis其中Q_ch是蓄冰功率进入冰槽的冷量η_ch是蓄冷效率η_dis是融冰效率。这个方程必须逐时段写进约束让蓄冰量S始终在0到槽容量之间同时给定初始值和末值约束。第三类是设备出力界限和逻辑约束。燃气轮机有出力上下限、最小运行状态电制冷机有制冷量上限蓄放冷功率有上限电制冷机的总制冷量Q_ET Q_direct Q_ch必须同时满足直接供冷和蓄冰两部分的需要。2.3 为什么这是一个MILP问题如果你把上面所有约束列出来会发现这个模型其实是线性的——功率、冷量、蓄冰量都是连续变量成本函数也是线性的。但问题出在设备启停逻辑上燃气轮机要么开要么停这个要么需要0-1变量蓄冰和融冰不能同时进行也需要0-1变量来做互斥约束。有0-1变量之后问题就从线性规划LP变成了混合整数线性规划MILP。很多初学者上来就用fmincon这类非线性求解器把0-1变量当成连续变量去优化出来的结果要么是0.37这种半开机状态要么就是违反物理逻辑的方案。我建议直接把模型规整成MILP交给成熟的商业求解器去解。24时段规模的MILP对Gurobi/CPLEX来说基本是秒解完全没必要在算法层面自己造轮子。3. MATLAB建模与求解的完整流程3.1 工具箱选型与环境配置MATLAB上做MILP求解我常用的组合是Yalmip Gurobi。Yalmip是一个建模层让约束和目标的写法接近数学模型本身代码可读性好很多Gurobi是底层求解器处理MILP的效率要比MATLAB自带的intlinprog高不少。如果你没有Gurobi的许可证用intlinprog也能跑只是规模大或者二进制变量多的时候求解时间会明显拉长。我实测下来24时段的冷电联供优化模型连续变量几十个、二进制变量二十几个intlinprog大概几秒到几十秒能出结果同样的模型让Gurobi去算基本是秒开。如果只是验证算法linprog/intlinprog够用如果要做多场景扫描或全年8760小时的规划建议直接上Gurobi或Cplex。另外新版本MATLAB对Yalmip的兼容性没问题但装工具箱的时候注意管理员权限和路径设置别装完发现Yalmip找不到solver。3.2 核心代码框架从参数到求解逐步拆解下面给一个可以直接改着用的核心框架。为了演示我把天然气价格按热值折算成了元/kWh热值写入参数电价和负荷数据留成占位符实际使用时用实测数据填充。首先定义基础参数T 24; dt 1; % 调度周期与步长(小时) P_GT_max 500; P_GT_min 100; % 燃气轮机出力界限 kW eta_GT 0.30; eta_rec 0.35; % 发电效率、余热回收效率 COP_AC 1.2; COP_EC 3.5; % 制冷机性能系数 Q_ICS_max 1000; Q_ch_max 250; % 蓄冰槽容量kWh、蓄冷功率上限kW Q_dis_max 250; eta_ch 0.85; % 放冷功率上限、蓄冷效率 eta_dis 0.90; c_gas 0.27; % 融冰效率、燃料成本(元/kWh热值) % c_buy、P_E_load、Q_C_load 根据实测数据预先赋值均为1xT数组然后定义变量和约束。蓄冰槽状态用1×(T1)的变量方便写递推关系同时用两个0-1变量分别控制机组启停和蓄/放冰互斥P_GT sdpvar(1, T); P_EC sdpvar(1, T); Q_ET sdpvar(1, T); Q_direct sdpvar(1, T); Q_ch sdpvar(1, T); Q_dis sdpvar(1, T); Q_AC sdpvar(1, T); P_buy sdpvar(1, T); S_ice sdpvar(1, T1); u_GT binvar(1, T); u_ch binvar(1, T); Constraints []; % 蓄冰槽状态转移 for t 1:T Constraints [Constraints, ... S_ice(t1) S_ice(t) eta_ch*Q_ch(t) - Q_dis(t)/eta_dis]; Constraints [Constraints, 0 S_ice(t) Q_ICS_max]; Constraints [Constraints, Q_ch(t) u_ch(t)*Q_ch_max]; Constraints [Constraints, Q_dis(t) (1-u_ch(t))*Q_dis_max]; end % 循环调度蓄冰槽当天回到零点也可改成末值初值 Constraints [Constraints, S_ice(1) 0, S_ice(T1) S_ice(1)]; % 电、冷平衡 Constraints [Constraints, P_GT P_buy P_E_load P_EC]; Constraints [Constraints, Q_AC Q_direct Q_dis Q_C_load]; Constraints [Constraints, Q_ET Q_direct Q_ch]; % 电制冷机与吸收式制冷 Constraints [Constraints, Q_ET COP_EC * P_EC, 0 Q_ET 400]; Constraints [Constraints, 0 Q_AC eta_rec*(1-eta_GT)*(P_GT/eta_GT)*COP_AC]; % 燃气轮机与电网购电上限 Constraints [Constraints, P_GT_min*u_GT P_GT P_GT_max*u_GT]; Constraints [Constraints, 0 P_buy 800];最后是目标函数和求解调用。注意所有成本项都要乘dtObjective sum(c_gas*(P_GT/eta_GT)*dt ... % 燃料成本 c_buy.*P_buy*dt ... % 购电成本 0.02*P_GT*dt ... % 燃气轮机运维 0.01*Q_ET*dt ... % 电制冷机运维 0.005*(Q_chQ_dis)*dt); % 蓄冰槽运维 options sdpsettings(solver,gurobi,verbose,0); sol optimize(Constraints, Objective, options);这段代码只是按示例参数写的实际项目里要根据设备手册把效率曲线、启停成本、最小运行时间补进去。另外如果燃气轮机余热还要分摊一部分去供热吸收式制冷的可用热量上限需要再乘以一个供热分配系数大家按自己的系统结构调整即可。3.3 结果后处理与可视化求解完之后用value()把变量取出来。先别急着画图第一步是检查约束是否闭合把每一时刻的电平衡、冷平衡、蓄冰状态递推都算一遍看误差在不在10的负6次方以内。这一步能挡住绝大多数看起来收敛但结果很蠢的情况。画图方面我通常画三张第一张是电功率平衡堆叠图把燃气轮机、购电、电制冷耗电、电负荷几根线叠在一起第二张是冷功率平衡图吸收式制冷、电制冷直接供冷、融冰放冷三条线叠一块第三张是蓄冰槽的容量变化曲线最直观地看出夜间蓄冰、白天放冷的调度效果。用MATLAB的area函数加legend就够了导出图片用exportgraphics新版MATLAB的矢量图导出质量会好很多。4. 有蓄冰和没蓄冰运行策略差在哪4.1 对比方案设置模型跑通之后我做了三组对比A方案不加蓄冰槽电制冷机直接供冷冷负荷全部靠吸收式制冷加电制冷满足。B方案蓄冰槽接入系统但采用固定策略——夜间23点到次日7点强制蓄冰、白天10点到18点强制放冷。C方案蓄冰槽接入系统蓄/放冷时机完全交给优化器决定就是上面那个MILP模型的最优解。三组用同一套设备参数和负荷数据只改变蓄冰是否参与、蓄冰策略是否优化。这样对比最有说服力先看蓄冰本身有没有价值再看加蓄冰但不优化策略和加蓄冰且优化策略之间差多少。4.2 典型日仿真结果解读我用一组某园区典型日的数据负荷数值做过脱敏处理跑了优化结果大致如下方案购电成本(元)燃气成本(元)运维成本(元)总成本(元)高峰购电峰值(kW)A 无蓄冰178022601604200680B 固定策略蓄冰145023802104040560C 优化蓄冰132023401803840470这是示例性结果不同电价条件和负荷形态下数值会有变化但趋势是一致的。先说结论从总成本看C方案比A方案大约低了8.6%效果很可观B方案虽然也有改善但明显不如C——这说明蓄冰装置不是装了就完事运行策略不当照样浪费设备潜力。从购电峰值看C方案的高峰购电功率比A方案低了200多kW对园区配电容量来说意义很大。再看机组出力曲线。C方案里燃气轮机在夜间低谷时段满负荷或接近满负荷运行多发的电除了满足夜间电负荷主要就是富余给电制冷机制冰白天高峰时段反而是吸收式制冷加融冰放冷在扛冷负荷电制冷机基本不在高电价时段开。这个把燃气轮机发电和电制冷造冰都挪到夜间、白天靠蓄冷放冷的策略就是冰蓄冷在冷电联供系统里的核心价值。4.3 蓄冰槽运行特性解读与末状态约束蓄冰槽的容量曲线有很强的时段特征夜间从零点开始库存爬升早晨7点左右达到峰值接近900kWh白天10点以后开始下降到晚上20点左右基本清空。这也符合直觉——优化器就是在夜间用低价电制冰、白天用冰代替高价电制冷之间做选择只不过它同时还考虑了燃气轮机的余热出力和电网购电上限所以蓄冰速率不是恒定值。提示蓄冰槽的末状态约束非常重要。如果只做单日优化而不加约束优化器会在第24小时把冰全部放光因为最后时刻剩余的冰没有任何收益。这样算出来的是单日内的最优却不是可持续的调度方案。我自己的做法是把末状态设成等于初始状态S(T1) S(1)这样每一天的蓄冰量是循环衔接的如果有明确的第二天初始库存数据也可以用固定值代替。5. 实际调试中踩过的坑5.1 求解返回无解或非数值结果的排查链路第一次跑MILP大概率会遇到infeasible或者Yalmip报NaN的报错。我总结了一个排查顺序基本能解决90%的问题第一步去掉所有0-1变量把模型降成纯LP先跑。如果LP都无解说明不是整数变量的问题而是约束相互矛盾。第二步逐条检查输入数据向量电价、负荷曲线里有没有NaN、Inf、负数特别是手填的Excel数据经常带空单元格。第三步检查蓄冰槽的递推约束——这是最容易出错的地方。S(1)0和S(T1)S(1)同时写的时候如果蓄冰/融冰能力小于负荷需求导致一天内无法清空约束冲突就会出现。解决办法是先把末状态约束放宽成S(T1)S(1)看有没有解。第四步用value()把每个约束的残差打印出来定位是哪一行约束不满足。这套流程走下来大多数无解问题半小时内就能定位。最怕的是不做排查直接调求解器参数调半天还是无解。5.2 非线性约束与物理可行性问题我之前有一版模型图省事写了一个Q_ch.*Q_dis0的约束表示蓄冰和融冰不能同时进行结果Yalmip直接报非凸约束Gurobi拒绝求解。后来改成用0-1变量u_ch来互斥问题就消失了。这种同一时刻动作方向互斥在储能建模里非常常见建议一律用二进制变量线性化不要写乘积约束。另外还要注意效率参数的方向。蓄冰槽的状态方程里蓄冷要乘效率、放冷要除以效率——如果两个方向都乘或者都除会导致系统越蓄越多的假象求解器会利用这个bug无限套利算出来的成本低到离谱。遇到成本结果低得不合理时优先查这类能量不守恒的约束。5.3 从24小时到全年的扩展思路做完单日优化很多人自然会想扩展到全年8760小时。这时候直接用MILP求解一整年变量规模会爆炸Gurobi也未必扛得住。我的做法是两步走第一步用K-means把全年负荷和电价曲线聚成3到5个典型日场景对典型日分别做调度优化得到的结果乘以该场景出现天数估算全年运行成本第二步如果要求更高再做滚动调度——比如每24小时窗口滚动一次只用未来48小时的数据既保证实时性又降低求解规模。这种从离线全局优化到在线滚动调度的思路恰恰是蓄能设备发挥作用的关键。蓄冰槽的核心价值就在于跨时段转移能量所以调度算法必须能处理昨天蓄的冰今天还能用这类窗口边缘问题滚动时要把当前蓄冰量作为初值传给下一轮。6. 给想复现这个课题的同学几点实在建议6.1 新手入门的搭建顺序不要一上来就写完整MILP。建议分四步走先跑固定电价加无蓄冰的纯LP再加分时电价再加蓄冰槽暂不加0-1互斥最后补启停逻辑和互斥变量。每一步都能验证结果合理性出问题也容易定位。我见过不少同学一步到位写完整模型结果无解了完全不知道从哪查起。6.2 容易出错的两个细节第一个是效率参数和单位统一。我见过太多人把kW和kWh混着用把效率乘反了结果优化器在凭空生能量。动手前先列一个单位表所有功率用kW、能量用kWh、电价用元/kWh、时间步长用小时能避免90%的低级错误。第二个是结果要人工检查。优化器给出的解未必物理上好看拿到结果后手动模拟一遍蓄冰量变化看看是否符合直觉有没有出现同时蓄冰又融冰空转启停这类现象。如果有说明约束漏了。6.3 这个模型还能怎么扩展这个方向再往后做可以考虑引入光伏或者风电的随机出力、需求侧响应、电池储能与蓄冰槽的配合优化。每一步加上去模型复杂度都会涨一截但底层框架还是本文这套功率平衡加蓄能连续性加设备界限加MILP求解。先把基础版本跑通再谈扩展。我在实际项目中最大的体会就是建模思路越贴近物理过程调试时越省心优化器给出的结果也越可信。 
长河君 · 2026/10/3 18:19:45 · 阅读 1.2k 冷电联供微网含冰蓄冷装置的经济优化运行：MATLAB建模与求解实践
最近在做一个含冰蓄冷装置的冷电联供型微网经济优化运行项目要用MATLAB把优化调度模型完整跑通。最开始我以为是常规的冷热电联供模型改改参数就行真正动手才发现冰蓄冷装置的引入让整个优化问题的复杂度抬高了一大截——不仅要同时平衡电负荷和冷负荷两条能量流还得处理蓄能设备带来的时序耦合约束。本文就把我用MATLAB从建模到求解的完整过程梳理一遍包括核心代码思路、典型日结果数据和调试中踩过的坑给正在做微网经济调度或者综合能源系统优化的朋友做个参考。1. 冷电联供微网为什么要配冰蓄冷装置1.1 系统结构冷、电、蓄三条能量流怎么耦合我接的这个项目场景是一个园区级微网主要负荷是冷负荷和电负荷。供能侧的核心是一台燃气轮机发出来的电优先满足自身负荷不够的部分从电网购入。燃气轮机的余热被余热锅炉回收驱动溴化锂吸收式制冷机供冷同时还有一台电制冷机作为冷负荷的补充和调节手段。关键就在于系统里加了一台冰蓄冷装置——夜间电价低谷时段让电制冷机多出力把冷量以冰的形式存起来白天电价高峰时再融冰放冷。这个系统的能量流比单纯的电热联供要复杂。因为冷和电两条能量流在电制冷机上耦合电制冷机既可以用电直接制冷供冷负荷也可以把制冷量灌进蓄冰槽。这样一来优化调度的问题就不只是设备开多少出力而是要同步决定燃气轮机发多少电、电网买多少电、电制冷机产生的冷量分多少给负荷、分多少给冰槽、冰槽什么时候蓄、什么时候放。每一步决策都会影响下一步的可行域这也是这类模型最需要花心思的地方。1.2 冰蓄冷的移峰填谷经济账冰蓄冷的核心是利用相变潜热储能1kg水变成冰要释放约334kJ的热量这个能量密度比水的显热蓄冷大得多所以在同样的蓄冷量需求下冰蓄冷槽的体积可以做得很小。工程上选它而不是大水罐核心就是看中这一点。经济账怎么算假设当地分时电价是峰时0.95元/kWh、谷时0.35元/kWh电制冷机的综合COP按3.5算。夜间花1kWh谷电制冰理论上能产出约3.5kWh的冷量即使考虑制冰工况下COP会降到3.0左右、蓄冰与融冰环节的综合效率按0.75算最后也能从冰槽取出约2.2kWh的冷量成本只有0.35元。这2.2kWh的冷量如果白天直接用高峰电价电制冷来补需要耗电约2.2/3.50.63kWh电费约0.6元。一来一回每kWh谷电用于制冰能省下约0.25元。一天蓄个几百kWh日运行成本降一两百块钱是现实的。当然实际优化中的蓄冰策略不是简单按峰谷电价差套利就完事。燃气轮机发电的同时会产生余热余热驱动的吸收式制冷边际成本很低天然气价格、购电价格、机组效率、冷负荷曲线之间是耦合的。所以最后到底什么时候蓄冰、蓄多少冰得靠优化模型去算而不是拍脑袋定。这也是本文要讲的核心。1.3 这类课题的典型应用场景带冰蓄冷的冷电联供微网最典型的应用场景就是白天冷负荷大、峰谷电价差距明显的园区——大型商业综合体、数据中心、食品冷链加工厂都属于这一类。这些地方空调或工艺制冷负荷占大头而且负荷曲线和电价曲线高度重叠白天越热越要制冷偏偏又赶上电价高峰。冰蓄冷正好能把制冷需求和电力消耗在时间上解耦。对做课题或者前期方案论证的人来说经济优化运行的意义在于设备都选型好了、参数都知道了怎么安排24小时内的出力计划让运行成本最低。这一步算得准后面无论是写可研报告还是做能量管理系统EMS的调度策略都有直接价值。2. 经济优化运行的数学模型怎么搭2.1 目标函数运行成本到底包含哪几项优化目标很直接一个调度周期通常取24小时内的总运行成本最小。我把成本拆成四项燃气轮机的燃料成本由发电功率和效率反推燃料消耗量再乘天然气价格。如果用线性化效率就是c_gas × (P_GT/η_GT) × Δt。电网购电成本按分时电价c_buy(t) × P_buy(t) × Δt计算。如果允许余电上网再减去售电收益。设备运行维护成本燃气轮机、电制冷机、蓄冰槽都按出力或启停次数计提系数一般是几厘到几分钱每kWh。碳排放成本可选如果课题要求考虑碳排放可以在目标函数里加一个CO2排放惩罚项。这部分对结果有一定影响但本文没把它放进主模型免得把问题绕大。需要特别提醒的是成本函数里所有功率×时间才是能量千万别忘了乘调度步长Δt。很多人第一次跑模型结果莫名其妙查到最后就是单位问题。2.2 约束条件的核心功率平衡、蓄能连续性、设备界限模型的约束我大致分成三类。第一类是瞬时平衡。每个时刻都要满足电功率平衡P_GT(t) P_buy(t) P_E_load(t) P_EC(t)其中P_EC是电制冷机耗电。冷功率平衡Q_AC(t) Q_direct(t) Q_dis(t) Q_C_load(t)其中Q_AC是吸收式制冷量Q_direct是电制冷机直接供冷量Q_dis是融冰放冷量。第二类是这类系统最关键的部分——蓄能连续性。蓄冰槽的状态方程是S(t1) S(t) η_ch × Q_ch(t) - Q_dis(t) / η_dis其中Q_ch是蓄冰功率进入冰槽的冷量η_ch是蓄冷效率η_dis是融冰效率。这个方程必须逐时段写进约束让蓄冰量S始终在0到槽容量之间同时给定初始值和末值约束。第三类是设备出力界限和逻辑约束。燃气轮机有出力上下限、最小运行状态电制冷机有制冷量上限蓄放冷功率有上限电制冷机的总制冷量Q_ET Q_direct Q_ch必须同时满足直接供冷和蓄冰两部分的需要。2.3 为什么这是一个MILP问题如果你把上面所有约束列出来会发现这个模型其实是线性的——功率、冷量、蓄冰量都是连续变量成本函数也是线性的。但问题出在设备启停逻辑上燃气轮机要么开要么停这个要么需要0-1变量蓄冰和融冰不能同时进行也需要0-1变量来做互斥约束。有0-1变量之后问题就从线性规划LP变成了混合整数线性规划MILP。很多初学者上来就用fmincon这类非线性求解器把0-1变量当成连续变量去优化出来的结果要么是0.37这种半开机状态要么就是违反物理逻辑的方案。我建议直接把模型规整成MILP交给成熟的商业求解器去解。24时段规模的MILP对Gurobi/CPLEX来说基本是秒解完全没必要在算法层面自己造轮子。3. MATLAB建模与求解的完整流程3.1 工具箱选型与环境配置MATLAB上做MILP求解我常用的组合是Yalmip Gurobi。Yalmip是一个建模层让约束和目标的写法接近数学模型本身代码可读性好很多Gurobi是底层求解器处理MILP的效率要比MATLAB自带的intlinprog高不少。如果你没有Gurobi的许可证用intlinprog也能跑只是规模大或者二进制变量多的时候求解时间会明显拉长。我实测下来24时段的冷电联供优化模型连续变量几十个、二进制变量二十几个intlinprog大概几秒到几十秒能出结果同样的模型让Gurobi去算基本是秒开。如果只是验证算法linprog/intlinprog够用如果要做多场景扫描或全年8760小时的规划建议直接上Gurobi或Cplex。另外新版本MATLAB对Yalmip的兼容性没问题但装工具箱的时候注意管理员权限和路径设置别装完发现Yalmip找不到solver。3.2 核心代码框架从参数到求解逐步拆解下面给一个可以直接改着用的核心框架。为了演示我把天然气价格按热值折算成了元/kWh热值写入参数电价和负荷数据留成占位符实际使用时用实测数据填充。首先定义基础参数T 24; dt 1; % 调度周期与步长(小时) P_GT_max 500; P_GT_min 100; % 燃气轮机出力界限 kW eta_GT 0.30; eta_rec 0.35; % 发电效率、余热回收效率 COP_AC 1.2; COP_EC 3.5; % 制冷机性能系数 Q_ICS_max 1000; Q_ch_max 250; % 蓄冰槽容量kWh、蓄冷功率上限kW Q_dis_max 250; eta_ch 0.85; % 放冷功率上限、蓄冷效率 eta_dis 0.90; c_gas 0.27; % 融冰效率、燃料成本(元/kWh热值) % c_buy、P_E_load、Q_C_load 根据实测数据预先赋值均为1xT数组然后定义变量和约束。蓄冰槽状态用1×(T1)的变量方便写递推关系同时用两个0-1变量分别控制机组启停和蓄/放冰互斥P_GT sdpvar(1, T); P_EC sdpvar(1, T); Q_ET sdpvar(1, T); Q_direct sdpvar(1, T); Q_ch sdpvar(1, T); Q_dis sdpvar(1, T); Q_AC sdpvar(1, T); P_buy sdpvar(1, T); S_ice sdpvar(1, T1); u_GT binvar(1, T); u_ch binvar(1, T); Constraints []; % 蓄冰槽状态转移 for t 1:T Constraints [Constraints, ... S_ice(t1) S_ice(t) eta_ch*Q_ch(t) - Q_dis(t)/eta_dis]; Constraints [Constraints, 0 S_ice(t) Q_ICS_max]; Constraints [Constraints, Q_ch(t) u_ch(t)*Q_ch_max]; Constraints [Constraints, Q_dis(t) (1-u_ch(t))*Q_dis_max]; end % 循环调度蓄冰槽当天回到零点也可改成末值初值 Constraints [Constraints, S_ice(1) 0, S_ice(T1) S_ice(1)]; % 电、冷平衡 Constraints [Constraints, P_GT P_buy P_E_load P_EC]; Constraints [Constraints, Q_AC Q_direct Q_dis Q_C_load]; Constraints [Constraints, Q_ET Q_direct Q_ch]; % 电制冷机与吸收式制冷 Constraints [Constraints, Q_ET COP_EC * P_EC, 0 Q_ET 400]; Constraints [Constraints, 0 Q_AC eta_rec*(1-eta_GT)*(P_GT/eta_GT)*COP_AC]; % 燃气轮机与电网购电上限 Constraints [Constraints, P_GT_min*u_GT P_GT P_GT_max*u_GT]; Constraints [Constraints, 0 P_buy 800];最后是目标函数和求解调用。注意所有成本项都要乘dtObjective sum(c_gas*(P_GT/eta_GT)*dt ... % 燃料成本 c_buy.*P_buy*dt ... % 购电成本 0.02*P_GT*dt ... % 燃气轮机运维 0.01*Q_ET*dt ... % 电制冷机运维 0.005*(Q_chQ_dis)*dt); % 蓄冰槽运维 options sdpsettings(solver,gurobi,verbose,0); sol optimize(Constraints, Objective, options);这段代码只是按示例参数写的实际项目里要根据设备手册把效率曲线、启停成本、最小运行时间补进去。另外如果燃气轮机余热还要分摊一部分去供热吸收式制冷的可用热量上限需要再乘以一个供热分配系数大家按自己的系统结构调整即可。3.3 结果后处理与可视化求解完之后用value()把变量取出来。先别急着画图第一步是检查约束是否闭合把每一时刻的电平衡、冷平衡、蓄冰状态递推都算一遍看误差在不在10的负6次方以内。这一步能挡住绝大多数看起来收敛但结果很蠢的情况。画图方面我通常画三张第一张是电功率平衡堆叠图把燃气轮机、购电、电制冷耗电、电负荷几根线叠在一起第二张是冷功率平衡图吸收式制冷、电制冷直接供冷、融冰放冷三条线叠一块第三张是蓄冰槽的容量变化曲线最直观地看出夜间蓄冰、白天放冷的调度效果。用MATLAB的area函数加legend就够了导出图片用exportgraphics新版MATLAB的矢量图导出质量会好很多。4. 有蓄冰和没蓄冰运行策略差在哪4.1 对比方案设置模型跑通之后我做了三组对比A方案不加蓄冰槽电制冷机直接供冷冷负荷全部靠吸收式制冷加电制冷满足。B方案蓄冰槽接入系统但采用固定策略——夜间23点到次日7点强制蓄冰、白天10点到18点强制放冷。C方案蓄冰槽接入系统蓄/放冷时机完全交给优化器决定就是上面那个MILP模型的最优解。三组用同一套设备参数和负荷数据只改变蓄冰是否参与、蓄冰策略是否优化。这样对比最有说服力先看蓄冰本身有没有价值再看加蓄冰但不优化策略和加蓄冰且优化策略之间差多少。4.2 典型日仿真结果解读我用一组某园区典型日的数据负荷数值做过脱敏处理跑了优化结果大致如下方案购电成本(元)燃气成本(元)运维成本(元)总成本(元)高峰购电峰值(kW)A 无蓄冰178022601604200680B 固定策略蓄冰145023802104040560C 优化蓄冰132023401803840470这是示例性结果不同电价条件和负荷形态下数值会有变化但趋势是一致的。先说结论从总成本看C方案比A方案大约低了8.6%效果很可观B方案虽然也有改善但明显不如C——这说明蓄冰装置不是装了就完事运行策略不当照样浪费设备潜力。从购电峰值看C方案的高峰购电功率比A方案低了200多kW对园区配电容量来说意义很大。再看机组出力曲线。C方案里燃气轮机在夜间低谷时段满负荷或接近满负荷运行多发的电除了满足夜间电负荷主要就是富余给电制冷机制冰白天高峰时段反而是吸收式制冷加融冰放冷在扛冷负荷电制冷机基本不在高电价时段开。这个把燃气轮机发电和电制冷造冰都挪到夜间、白天靠蓄冷放冷的策略就是冰蓄冷在冷电联供系统里的核心价值。4.3 蓄冰槽运行特性解读与末状态约束蓄冰槽的容量曲线有很强的时段特征夜间从零点开始库存爬升早晨7点左右达到峰值接近900kWh白天10点以后开始下降到晚上20点左右基本清空。这也符合直觉——优化器就是在夜间用低价电制冰、白天用冰代替高价电制冷之间做选择只不过它同时还考虑了燃气轮机的余热出力和电网购电上限所以蓄冰速率不是恒定值。提示蓄冰槽的末状态约束非常重要。如果只做单日优化而不加约束优化器会在第24小时把冰全部放光因为最后时刻剩余的冰没有任何收益。这样算出来的是单日内的最优却不是可持续的调度方案。我自己的做法是把末状态设成等于初始状态S(T1) S(1)这样每一天的蓄冰量是循环衔接的如果有明确的第二天初始库存数据也可以用固定值代替。5. 实际调试中踩过的坑5.1 求解返回无解或非数值结果的排查链路第一次跑MILP大概率会遇到infeasible或者Yalmip报NaN的报错。我总结了一个排查顺序基本能解决90%的问题第一步去掉所有0-1变量把模型降成纯LP先跑。如果LP都无解说明不是整数变量的问题而是约束相互矛盾。第二步逐条检查输入数据向量电价、负荷曲线里有没有NaN、Inf、负数特别是手填的Excel数据经常带空单元格。第三步检查蓄冰槽的递推约束——这是最容易出错的地方。S(1)0和S(T1)S(1)同时写的时候如果蓄冰/融冰能力小于负荷需求导致一天内无法清空约束冲突就会出现。解决办法是先把末状态约束放宽成S(T1)S(1)看有没有解。第四步用value()把每个约束的残差打印出来定位是哪一行约束不满足。这套流程走下来大多数无解问题半小时内就能定位。最怕的是不做排查直接调求解器参数调半天还是无解。5.2 非线性约束与物理可行性问题我之前有一版模型图省事写了一个Q_ch.*Q_dis0的约束表示蓄冰和融冰不能同时进行结果Yalmip直接报非凸约束Gurobi拒绝求解。后来改成用0-1变量u_ch来互斥问题就消失了。这种同一时刻动作方向互斥在储能建模里非常常见建议一律用二进制变量线性化不要写乘积约束。另外还要注意效率参数的方向。蓄冰槽的状态方程里蓄冷要乘效率、放冷要除以效率——如果两个方向都乘或者都除会导致系统越蓄越多的假象求解器会利用这个bug无限套利算出来的成本低到离谱。遇到成本结果低得不合理时优先查这类能量不守恒的约束。5.3 从24小时到全年的扩展思路做完单日优化很多人自然会想扩展到全年8760小时。这时候直接用MILP求解一整年变量规模会爆炸Gurobi也未必扛得住。我的做法是两步走第一步用K-means把全年负荷和电价曲线聚成3到5个典型日场景对典型日分别做调度优化得到的结果乘以该场景出现天数估算全年运行成本第二步如果要求更高再做滚动调度——比如每24小时窗口滚动一次只用未来48小时的数据既保证实时性又降低求解规模。这种从离线全局优化到在线滚动调度的思路恰恰是蓄能设备发挥作用的关键。蓄冰槽的核心价值就在于跨时段转移能量所以调度算法必须能处理昨天蓄的冰今天还能用这类窗口边缘问题滚动时要把当前蓄冰量作为初值传给下一轮。6. 给想复现这个课题的同学几点实在建议6.1 新手入门的搭建顺序不要一上来就写完整MILP。建议分四步走先跑固定电价加无蓄冰的纯LP再加分时电价再加蓄冰槽暂不加0-1互斥最后补启停逻辑和互斥变量。每一步都能验证结果合理性出问题也容易定位。我见过不少同学一步到位写完整模型结果无解了完全不知道从哪查起。6.2 容易出错的两个细节第一个是效率参数和单位统一。我见过太多人把kW和kWh混着用把效率乘反了结果优化器在凭空生能量。动手前先列一个单位表所有功率用kW、能量用kWh、电价用元/kWh、时间步长用小时能避免90%的低级错误。第二个是结果要人工检查。优化器给出的解未必物理上好看拿到结果后手动模拟一遍蓄冰量变化看看是否符合直觉有没有出现同时蓄冰又融冰空转启停这类现象。如果有说明约束漏了。6.3 这个模型还能怎么扩展这个方向再往后做可以考虑引入光伏或者风电的随机出力、需求侧响应、电池储能与蓄冰槽的配合优化。每一步加上去模型复杂度都会涨一截但底层框架还是本文这套功率平衡加蓄能连续性加设备界限加MILP求解。先把基础版本跑通再谈扩展。我在实际项目中最大的体会就是建模思路越贴近物理过程调试时越省心优化器给出的结果也越可信。