我最早接触含冰蓄冷的冷电联供微网是在一个制冷季的能源审计项目里。白天办公楼冷负荷爆表正好撞上电价高峰段夜里空调没需求电价反而跌到谷段的三分之一。当时脑子里蹦出来的第一个念头就是要是能在夜里把电变成冰存起来白天再用这些冰扛冷负荷电费不就省下来了等真正把它写进微网经济优化模型里才发现冰蓄冷远不止“搬移冷量”这么简单——它会牵动燃气轮机出力、电网购电、电制冷机投停、蓄冰槽调度的一连串连锁反应。这篇文章就是把我在MATLAB里完整跑通“含冰蓄冷装置的冷电联供型微网经济优化运行”的过程、模型思路和踩过的坑整理出来适合正在做综合能源、微电网优化方向的研究生以及做园区能源规划的工程师参考。这类课题的核心其实可以拆成三个问题系统里有哪些设备、冷和电怎么流动经济优化的目标函数和约束条件怎么用数学语言表达以及如何用MATLAB把模型建起来、求解出来、把结果讲清楚。下面我就按这个逻辑一步步展开最后把代码实现和典型日算例结果一起贴出来保证你照着做能跑出图。1. 先把这个系统的能量流彻底捋清楚1.1 微网里都有谁谁在供电谁在供冷含冰蓄冷的冷电联供微网典型的设备清单大概是这么一套燃气轮机或者燃气内燃机、余热制冷机组通常是溴化锂吸收式、电压缩制冷机、冰蓄冷装置再加上电负荷和冷负荷这两个“必须满足的口子”。有些配置里还会加燃气锅炉作为备用热源但这个课题如果只关心冷电联供燃气锅炉可以先不放进模型。燃气轮机烧天然气发电电优先供给微网内的电负荷发完电剩下的烟气余热不能浪费送进溴化锂机组去制冷这部分冷量就是“联供”的核心产出。电制冷机的角色比较灵活白天冷负荷高的时候它可以直接制冷补缺口夜里电价便宜的时候它又可以把冷量以冰的形式存在蓄冰槽里。而冰蓄冷装置本身不产冷它只是把“夜间便宜的电”转化成“白天能用的冷”相当于在时间轴上做能量的搬运。把这套系统画成能量流来看燃气轮机的出口有两条支路一条是电母线一条是烟气余热母线电母线上挂着电负荷、电制冷机以及一个可选的“向电网购电/售电”接口冷母线上挂着溴化锂机组、电制冷机、融冰出口和冷负荷。优化调度做的事情就是每个调度时段决定燃气轮机发多少电、电网买多少电、电制冷机产多少冷且多少直接供冷多少去蓄冰、蓄冰槽融多少冰出来最终让电和冷两条母线都刚好平衡。这里我建议初学者先把“冷母线”和“电母线”分开画两个口子的平衡约束分开写。很多第一次做这个课题的人容易把电制冷机的耗电功率和它的产冷量混在一起一写约束就乱套。记住一条电平衡只算功率冷平衡只算冷量电制冷机就是它们之间的“转换器”用COP把两者勾起来。1.2 冰蓄冷为什么能省钱它的“天时”在哪冰蓄冷省钱的根本逻辑不是设备效率而是时间套利。理解这件事特别重要因为它决定了你后面积分时电价、冷负荷曲线时该怎么思考。举个生活的例子假设你所在的地区白天电价1.2元/度夜间谷段电价0.4元/度白天开空调制冷需要的电费是晚上的三倍。冰蓄冷相当于你在晚上用电把水冻成冰白天冰块融化的冷量来吹空调等于把白天的用电量“挪”到了晚上交费。表面上总耗电量没怎么变但电费结构变了这就是它经济性的第一来源。在含燃气轮机的冷电联供微网里经济性会比“单纯峰谷套利”复杂一层。因为燃气轮机烧燃气发电也要算燃料成本它余热制冷还能顶掉一部分电制冷机的耗电。所以系统最优调度不是简单地“谷段蓄冰、峰段融冰”而是要全局打量某个时段到底是买电网的电更划算还是自己烧气发电更划算发出来的电是给电负荷用还是给电制冷机用冷负荷缺口是用溴化锂余热供、还是电制冷直供、还是把蓄冰槽里的冰放出来。这些都是优化模型需要解决的问题。还有一点容易被忽略并不是所有冷负荷都适合用冰蓄冷来扛。冰蓄冷装置有一个蓄冰容量上限蓄满了就是蓄满了放空了就是放空了。而且融冰出力的速率也有限制不能指望一个小时内把整夜的冰全部融完。所以它更像一个“蓄水池”有容积约束、有进出水速率约束调度模型必须把它当储能设备来建模而不是一个可以任意使用的“免费冷源”。2. 经济优化模型究竟在优化什么2.1 目标函数一天下来究竟花多少钱经济优化运行的“经济”二字必须落到一个可量化的函数上。最常见的做法是取一个典型日把全天分成T个调度时段通常取24个时段每小时一个点目标函数是这24个小时里微网的总运行成本最小化。总运行成本一般包括三块从电网购电的费用、燃气轮机消耗天然气的燃料费、以及各设备的运维费用。写出来大概是这样的形式min ∑(购电电价 × 购电功率 气价 × 燃气耗量 运维系数 × 出力)其中购电电价是分段函数峰、平、谷各时段电价不一样这在MATLAB里可以用一个长度24的向量直接表示每一个元素就是对应时段的电价。燃气耗量跟燃气轮机出力之间的关系可以用一个线性的热耗率表达式来近似燃气耗量 燃气轮机出力 / 发电效率再换算成燃料费。运维费用通常按设备的运行出力乘以一个很小的单位运维成本系数来估算比如电制冷机每产出1kW冷量收0.01元、燃气轮机每发1kW电收0.02元这部分费用虽然数额不大但加入目标函数之后能让优化结果更贴近实际避免出现设备出力频繁跳变的毛刺。有些模型还会考虑微网向电网售电的收入也就是如果燃气轮机发多了电、用不完可以卖给电网赚钱。这种情况下目标函数里就要减掉卖电收入但大多数课题为了聚焦冰蓄冷的调度逻辑会先假设微网不从电网卖电或者只允许从电网买电不允许反向送电。建议你入门阶段先做成“只买不卖”把主框架跑通了再放开。2.2 约束条件电力平衡、冷量平衡和蓄冰槽的“脾气”目标函数定了之后约束条件就是模型的骨架。缺了约束的优化就是一个无约束最大化利润问题很荒谬。对这套系统来说至少有以下几类约束必须写清楚。第一是电功率平衡约束。每个时段燃气轮机出力加上从电网购电的功率必须等于电负荷加上所有用电设备的耗电功率。用电设备包括电制冷机注意它既可能直接制冷也可能蓄冰两部分耗电都要计入以及其他辅助设备。这条约束是等式约束物理含义是“微网内部电能不能凭空产生或消失”。第二是冷功率平衡约束。每个时段溴化锂机组产冷量、电制冷机直供冷量、蓄冰槽融冰供冷量三者之和必须等于冷负荷。同样也是等式约束。这里溴化锂机组的产冷量跟燃气轮机的余热息息相关一般用热电比或者余热回收效率把它和燃气轮机出力耦合成一个线性关系。第三是蓄冰槽的动态约束。这是整个模型中最有“储能味”的一条。蓄冰槽的蓄冷状态可以用SOCState of Charge荷冷状态来表示它像一个水桶本时段结束时桶里的水量等于上一时段剩下的水量加上本时段蓄进去的冷量蓄冰减去本时段放出来的冷量融冰。写成离散形式就是SOC(t1) SOC(t) (蓄冰功率 - 融冰功率) × 时段时长 / 蓄冰槽容量SOC不能小于0也不能大于1这对应蓄冰槽不能“透支”也不能“溢出”。另外蓄冰槽的蓄冰速率和融冰速率都各自有上限取决于制冷机和融冰换热器的实际能力。对于完整的日调度通常还会要求SOC在一天开始和一天结束时保持一致代表这个调度策略可以循环运行。第四是设备出力上下限约束和爬坡约束。燃气轮机有最小稳定出力和最大额定出力电制冷机的制冷能力也有上下限。爬坡约束描述的是设备从一个时段的出力跳到下一个时段时变化速度不能太快比如燃气轮机每15分钟最多只能爬升百分之多少这对燃机实际运行很重要。第五是互斥约束。同一时刻蓄冰槽在蓄冰和在融冰是两种矛盾状态理论上不应该同时发生。这需要引入二进制变量比如用0-1变量y表示蓄冰工况y1时只能蓄冰、不能融冰y0时只能融冰、不能蓄冰。加了二进制变量整个规划问题就变成了混合整数线性规划MILP求解规模会变大但这也是这类课题的标准做法。2.3 关键参数的计算与单位陷阱参数处理是实操中第一个“翻车高发区”。我见过太多同学把冷负荷单位搞混结果优化出来的结果全是“天文数字”。这里必须提一下几个最常见的单位坑。冷量单位有kW和RT冷吨两种二者之间的换算是1RT约等于3.517kW。很多文献里冷负荷曲线用的是RT电负荷用的是kW如果直接放进同一个模型里算冷平衡约束就废了。我习惯的做法是在数据读入MATLAB之后统一把冷负荷单位转换成kW再参与计算所有设备参数也都转成kW制最后出图时再根据需要标回RT。蓄冰槽的容量常见有两种标法一种是用蓄冷量表示比如“××RTh”或“××kWh”另一种是用蓄冰量表示比如“××吨冰”。如果是后者换算更麻烦一点因为冰的融解热是334kJ/kg1吨冰全部融化大约可以吸收93kWh左右的热量。实际建模时我建议直接以蓄冰槽的冷量容量为基准定义SOC这样最简单。COP的取值也要分场景。同样是那台电制冷机直接制冷的COP和制冰工况的COP通常不一样制冰工况的冷凝温度更低、效率大约会打个八折。模型里最好拆成两个参数COP_direct和COP_ice蓄冰功率除以COP_ice才是它的耗电功率。很多简化模型干脆只用一个COP误差其实不小尤其当蓄冰量占比较大时算出来的夜间耗电量会偏低省电费效果会被高估。还有一个容易忽略的是分时电价边界时段的处理。比如峰段从8:00开始那么7:59和8:00这两个时段如果分别落在谷段和峰段电价会出现一个突跳。优化模型会强烈倾向于把大量蓄冰安排在谷段末尾这是合理的但如果蓄冰速率上限设置得过大模型可能安排极端出力出现实际设备无法响应的调度指令。所以蓄冰速率上限要按设备的真实能力来约束不能拍脑袋给一个很大的数。3. MATLAB求解路线怎么选3.1 方案一YALMIP Cplex 精确求解做混合整数线性规划我在MATLAB里最常用的组合是YALMIP加外部求解器。YALMIP是一个建模工具箱它最大的价值是让你把优化问题写成接近数学公式的代码而不用关心底层求解器的具体调用方式。Cplex或者Gurobi扮演的是求解引擎的角色负责真正把MILP解出来。两者配合起来写代码的效率远高于手写linprog矩阵。为什么推荐这条路因为这个课题的模型里天然有二进制变量用YALMIP可以非常自然地写出“蓄冰和融冰互斥”这类逻辑约束。如果不用YALMIP你自己用ma tlab的intlinprog来建模就必须把所有的约束条件手动整理成矩阵形式Ax≤b而且要把二进制变量和连续变量拼成一个长向量变量一多下标映射能把人绕晕。YALMIP的典型写法是T 24; x sdpvar(T, 1); % 连续变量 y binvar(T, 1); % 二进制变量 Constraints []; for t 1:T Constraints [Constraints, x(t) 0, x(t) 100]; end Cost sum(x); optimize(Constraints, Cost, sdpsettings(solver, cplex));跑完用value()函数取出各变量的数值然后直接plot非常顺手。需要提醒的是Cplex和Gurobi都是商业求解器但如果是在校学生可以通过学校许可证或者求解器官方申请学术版授权。实在装不上也可以用YALMIP内置支持的免费求解器如CBC、GLPK。免费的求解器速度会差一些但这个课题的规模不大一般也能在可接受的时间内跑完。3.2 方案二启发式算法兜底有些朋友的环境里装不了外部求解器或者导师坚持要求用智能优化算法这时会考虑粒子群、遗传算法这类启发式方法。我不拦着但必须先说清楚它的代价。启发式算法最麻烦的是处理约束。MILP里那些等式约束电平衡、冷平衡和不等式约束在粒子群算法里不能像YALMIP那样直接写进去。常见的做法是把约束转成罚函数加到目标函数里或者设计满足约束的编码方式。罚函数法说白了就是“如果违反约束就在目标函数里加一个很大的惩罚值”但惩罚系数怎么取是个玄学取小了结果违反约束取大了惩罚项直接盖过真实成本优化过程变成在瞎找。我见过不少代码罚函数系数拍脑袋设了个10000结果目标函数全是罚出来的数完全失去经济意义。另一个问题是重复性差。启发式算法有随机性每次跑的最终结果都不一样除非设定随机数种子。写论文做对比实验时“优化结果不稳定”这个缺点非常恼人。所以我的建议是如果条件允许优先用精确求解器跑MILP启发式算法可以作为对比方法在论文里验证智能算法的有效性但不建议一开始就陷进去。3.3 核心代码结构与关键片段下面这段是我整理的一套可运行的YALMIP关键代码骨架覆盖了变量定义、约束书写和求解调用。参数部分我统一写成一个结构体方便后续修改。%% 数据准备假设每个变量都是长度为T的列向量 T 24; dt 1; % 调度时段数和步长 P_load [...]; % 24h电负荷kW Q_load [...]; % 24h冷负荷kW注意单位统一 c_buy [谷段电价; 峰段电价; ...]; % 24h购电价 % 设备参数 cap_ice 500; % 蓄冰槽容量kWh rate_charge 150; % 蓄冰速率上限kW rate_melt 200; % 融冰速率上限kW COP_ec 3.5; % 电制冷机直供COP COP_ice 2.8; % 电制冷机制冰工况COP COP_br 1.2; % 溴化锂余热制冷系数 eta_gt 0.35; % 燃气轮机发电效率 gas_price 2.6; % 天然气价元/m3 P_gt_max 500; P_gt_min 50; % 燃气轮机出力上下限 Q_ec_max 600; % 电制冷机最大产冷量 %% 变量定义 P_gt sdpvar(T,1); % 燃气轮机出力 P_buy sdpvar(T,1); % 购电功率 Q_ec_direct sdpvar(T,1); % 电制冷机直供冷量 Q_ec_charge sdpvar(T,1); % 电制冷机用于蓄冰的冷量 Q_melt sdpvar(T,1); % 融冰供冷量 soc sdpvar(T1,1); % 蓄冰槽SOC y binvar(T,1); % 蓄冰/融冰互斥标志 %% 约束 Constraints []; for t 1:T % 电平衡 Constraints [Constraints, P_gt(t) P_buy(t) P_load(t) ... Q_ec_direct(t)/COP_ec Q_ec_charge(t)/COP_ice]; % 冷平衡 Constraints [Constraints, COP_br*P_gt(t) Q_ec_direct(t) Q_melt(t) Q_load(t)]; % 设备上下限 Constraints [Constraints, P_gt_min P_gt(t) P_gt_max]; Constraints [Constraints, 0 P_buy(t) 1000]; Constraints [Constraints, 0 Q_ec_direct(t) Q_ec_max]; Constraints [Constraints, 0 Q_ec_charge(t) rate_charge]; Constraints [Constraints, 0 Q_melt(t) rate_melt]; % 互斥约束 Constraints [Constraints, Q_ec_charge(t) rate_charge*y(t)]; Constraints [Constraints, Q_melt(t) rate_melt*(1-y(t))]; % 蓄冰槽动态 Constraints [Constraints, soc(t1) soc(t) ... (Q_ec_charge(t) - Q_melt(t))*dt/cap_ice]; Constraints [Constraints, 0 soc(t1) 1]; end % 循环周期约束 Constraints [Constraints, soc(1) soc(T1)]; %% 目标函数与求解 Cost sum(c_buy.*P_buy*dt ... gas_price*P_gt*dt/eta_gt/10 ... % 天然气耗量换算注意气热值单位 0.02*P_gt 0.01*(Q_ec_directQ_ec_charge)/COP_ec*dt); ops sdpsettings(solver,cplex,verbose,1); result optimize(Constraints, Cost, ops); if result.problem 0 P_gt_opt value(P_gt); P_buy_opt value(P_buy); soc_opt value(soc); else disp(求解失败检查约束和求解器设置); end这段代码我尽量保持了精简但核心要素都齐全。运行起来之后你可以把P_gt_opt、soc_opt这些结果画成柱状图或阶梯图观察各设备的出力时序。如果遇到具体报错后面第5节我会分析常见问题。4. 典型日算例结果怎么看4.1 场景与参数设定为了让你对结果有个直觉我说一个自己调过的典型算例。假设一个商业园区冷负荷高峰出现在下午14:00到17:00峰值大约是800kW电负荷跟冷负荷走势相近但耦合不强峰值约600kW。分时电价按常见的三段式谷段0:00-8:00电价为0.4元/kWh平段8:00-11:00和17:00-21:00为0.8元/kWh峰段11:00-17:00和21:00-24:00为1.2元/kWh。燃气轮机最大出力500kW最小稳定出力50kW发电效率35%溴化锂余热制冷系数取1.2电制冷机最大产冷量600kWCOP取3.5制冰工况COP取2.8蓄冰槽容量500kWh蓄冰速率上限150kW融冰速率上限200kW。天然气价按2.6元/m³再按天然气热值折算单位发电成本。这个场景我跑出来的典型优化结果有几个鲜明特征。第一蓄冰槽的SOC曲线呈现出非常典型的“夜间蓄、白天放”形态从凌晨1点左右开始SOC一路爬升到早上8点前达到接近1然后从上午10点左右开始持续下降到傍晚冷负荷高峰期基本放空。第二燃气轮机并非全天满发而是在电价较高的峰段加大出力在谷段则压低出力甚至接近最小出力因为谷段买电网电更划算。第三电制冷机在谷段大量运行一部分冷量直供夜间的基数冷负荷大部分冷量用于制冰蓄存到了白天峰段电制冷机的直供冷量会被压缩冷负荷主要靠融冰和余热制冷顶上。4.2 优化结果怎么解读对比数据才有说服力光看一个优化结果没啥说服力做研究必须要有对比。最直接也最常用的对比方案是设一个“无冰蓄冷”的对照组把冰蓄冷相关变量全部置0蓄冰槽容量设为0其他参数保持不变重新求解一个不含冰蓄冷的冷电联供微网优化模型。两个模型跑出来的总运行成本一对比冰蓄冷的经济价值就出来了。在我这个算例里含冰蓄冷方案的全天总运行成本大约是4800元无冰蓄冷方案大约是5500元左右节约比例在10%到15%。这个数字会让你直观感受到“移峰填谷”的能量有多大。你还可以继续做敏感性分析拉大峰谷电价差看节约比例怎么变或者把蓄冰槽容量从300改到800看经济性提升是否边际递减。这类分析放在论文里就是非常扎实的章节。结果图的呈现上我建议至少画三张图。第一张是蓄冰槽SOC曲线一眼看出蓄放规律第二张是冷负荷供需堆叠图把溴化锂供冷、电制冷直供、融冰供冷三个部分堆叠起来跟冷负荷曲线对整齐谁在什么时段出力一目了然第三张是电功率平衡图把燃气轮机出力、电网购电、电负荷和电制冷耗电画在一起。这几张图画完之后整个系统的调度策略就非常立体了。5. 实操中常见的问题与排查技巧5.1 求解器与MATLAB代码层面的报错我在调试这套模型时遇到过不少报错挑几个高频的写在这里供你排查。“No suitable solver found”是YALMIP最常见的报错。它表示求解器没有正确配置到YALMIP的路径里。解决方法是检查Cplex或Gurobi的MATLAB接口是否安装然后用yalmiptest命令测试一下哪些求解器可用。对于MILP问题YALMIP需要有支持整数规划的求解器如果你只有linprog那确实解不了带binvar的模型。一个应急办法是把binvar去掉假设蓄冰和融冰可以同时发生模型退化成LP用自带linprog也能跑但结果在物理上不严格。“Index exceeds array bounds”这种报错大多是变量维度不匹配。YALMIP里如果soc定义成(T1,1)而循环里误写成soc(t)就很容易越界。另外约束中两个变量相加时要注意它们都是列向量还是行向量。MATLAB对1×24和24×1的隐式扩展有时能算但结果会莫名其妙。我的习惯是所有时间序列变量统一定义为T×1列向量绝不混用。“Infeasible problem”提示无可行解是最让人头疼的。我第一次跑就遇上了后来发现是SOC的周期约束和初始SOC设置矛盾了我既设了soc(1)0.5又要求soc(T1)soc(1)这本来没问题问题是我在凌晨安排了很大的蓄冰量导致soc(1)0.5的约束限制了初值偏高无法在24小时内完成一次完整的蓄放循环。解决方法是先跑一次不带SOC周期约束的模型观察末尾SOC自然落点再来决定周期约束怎么写或者干脆让soc(1)soc(T1)0表示一天从空槽开始也相对合理。5.2 模型与工程层面的坑单位不一致是最隐蔽的坑。冷负荷曲线如果是从RT标定的而蓄冰槽容量是用kWh标定的你必须在数据准备阶段就统一成kW和kWh。我建议在代码开头加一个单位转换参数比如RT_TO_KW 3.517所有冷量数据乘完再进模型同时用注释写明转换关系免得过几天自己都忘了。蓄冰槽的自损耗也不要完全忽略。实际蓄冰槽一天下来会有5%到10%的冷量损失如果你完全不建模优化结果会过于乐观。最简单的方法是在SOC动态约束里加一个衰减系数比如SOC(t1) 0.98×SOC(t) 入流 - 出流虽然只是粗略近似但结果的可信度提升不少。电制冷机同时承担直供和制冰时COP分开取值的细节前面已经强调过。这里再补充一个实操技巧如果你发现优化结果中电制冷机蓄冰速率频繁触及上限说明蓄冰槽容量相对冷负荷峰值偏大经济上出现了“为了蓄冰而蓄冰”的边际效益递减。此时可以做一个蓄冰槽容量敏感性扫描找出边际收益骤降的拐点这个拐点就是合理的容量配置。最后提醒一点关于爬坡约束。如果燃气轮机的爬坡约束设置得太紧求解时间会明显变长因为二进制变量和连续变量之间的耦合更紧了。我一般先不加入爬坡约束快速验证整个模型框架无误后再把它加回去做完整仿真。这样能大幅缩短前期调试时间。我个人做完这个课题后最大的体会是冰蓄冷装置的引入真正改变了微网调度的“相位”——把电力需求从白天挪到黑夜把冷量供给从即时生产变成了“先储后放”。MATLAB在这个课题里最大的价值是借助YALMIP把模型和算法解耦让研究者的精力集中在物理建模和结果分析上而不是陷入求解器的矩阵地狱。最后再分享一个小技巧拿到任何冷电联供微网优化模型先跑一个不蓄冰的纯电制冷场景把电平衡、冷平衡和负荷曲线彻底验证对了再一步步加入余热制冷、冰蓄冷、互斥约束。按这个顺序递进你会发现自己调试模型的效率至少翻一倍。