企业官网建设流程全解析

最近在做一个含冰蓄冷装置的冷电联供型微网经济优化运行项目，要用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.5=0.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(t+1) = 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×(T+1)的变量，方便写递推关系；同时用两个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, T+1); u_GT = binvar(1, T); u_ch = binvar(1, T); Constraints = []; % 蓄冰槽状态转移 for t = 1:T Constraints = [Constraints, ... S_ice(t+1) == 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(T+1) == 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];

最后是目标函数和求解调用。注意所有成本项都要乘dt：

Objective = 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_ch+Q_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 典型日仿真结果解读

我用一组某园区典型日的数据（负荷数值做过脱敏处理）跑了优化，结果大致如下：

这是示例性结果，不同电价条件和负荷形态下数值会有变化，但趋势是一致的。先说结论：从总成本看，C方案比A方案大约低了8.6%，效果很可观；B方案虽然也有改善，但明显不如C——这说明蓄冰装置不是装了就完事，运行策略不当照样浪费设备潜力。从购电峰值看，C方案的高峰购电功率比A方案低了200多kW，对园区配电容量来说意义很大。

再看机组出力曲线。C方案里燃气轮机在夜间低谷时段满负荷或接近满负荷运行，多发的电除了满足夜间电负荷，主要就是富余给电制冷机制冰；白天高峰时段反而是吸收式制冷加融冰放冷在扛冷负荷，电制冷机基本不在高电价时段开。这个"把燃气轮机发电和电制冷造冰都挪到夜间、白天靠蓄冷放冷"的策略，就是冰蓄冷在冷电联供系统里的核心价值。

4.3 蓄冰槽运行特性解读与末状态约束

蓄冰槽的容量曲线有很强的时段特征：夜间从零点开始库存爬升，早晨7点左右达到峰值接近900kWh；白天10点以后开始下降，到晚上20点左右基本清空。这也符合直觉——优化器就是在"夜间用低价电制冰、白天用冰代替高价电制冷"之间做选择，只不过它同时还考虑了燃气轮机的余热出力和电网购电上限，所以蓄冰速率不是恒定值。

提示：蓄冰槽的末状态约束非常重要。如果只做单日优化而不加约束，优化器会在第24小时把冰全部放光，因为"最后时刻剩余的冰"没有任何收益。这样算出来的是单日内的最优，却不是可持续的调度方案。我自己的做法是把末状态设成等于初始状态（S(T+1) = S(1)），这样每一天的蓄冰量是循环衔接的；如果有明确的"第二天初始库存"数据，也可以用固定值代替。

5. 实际调试中踩过的坑

5.1 求解返回"无解"或"非数值结果"的排查链路

第一次跑MILP，大概率会遇到infeasible或者Yalmip报NaN的报错。我总结了一个排查顺序，基本能解决90%的问题：

第一步，去掉所有0-1变量，把模型降成纯LP先跑。如果LP都无解，说明不是整数变量的问题，而是约束相互矛盾。

第二步，逐条检查输入数据向量：电价、负荷曲线里有没有NaN、Inf、负数，特别是手填的Excel数据经常带空单元格。

第三步，检查蓄冰槽的递推约束——这是最容易出错的地方。S(1)==0和S(T+1)==S(1)同时写的时候，如果蓄冰/融冰能力小于负荷需求导致一天内无法清空，约束冲突就会出现。解决办法是先把末状态约束放宽成S(T+1)>=S(1)，看有没有解。

第四步，用value()把每个约束的残差打印出来，定位是哪一行约束不满足。

这套流程走下来，大多数"无解"问题半小时内就能定位。最怕的是不做排查直接调求解器参数，调半天还是无解。

5.2 非线性约束与物理可行性问题

我之前有一版模型图省事，写了一个Q_ch.*Q_dis==0的约束表示"蓄冰和融冰不能同时进行"，结果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求解"。先把基础版本跑通，再谈扩展。我在实际项目中最大的体会就是：建模思路越贴近物理过程，调试时越省心，优化器给出的结果也越可信。