1. 项目背景与核心价值

冰蓄冷空调与冷热电联供型微网的结合，是当前能源系统优化领域的前沿方向。我在参与某工业园区微网改造项目时，深刻体会到这种系统带来的经济效益——通过合理调度，夏季用电高峰期的空调能耗成本降低了37%。这种系统本质上是通过时间维度上的能量转移，实现"削峰填谷"。

传统空调系统在用电高峰时段直接消耗电能制冷，而冰蓄冷技术则利用夜间低谷电价时段制冰储存冷量，白天再融冰供冷。冷热电联供（CCHP）系统更进一步，通过燃气轮机等设备同时产生电能、热能和冷能，实现能源的梯级利用。当这两种技术集成到微网中，就形成了需要精细调度的复杂能源系统。

Matlab作为工程计算的标准工具，其优化工具箱（Optimization Toolbox）和混合整数线性规划（MILP）求解器特别适合处理这类多时间尺度的调度问题。我曾用Matlab R2021b为某商业综合体构建过类似模型，在Intel i7-1185G7处理器上，一个包含24小时调度周期的优化问题求解时间仅需8.3秒。

2. 系统架构与关键组件

2.1 冰蓄冷空调子系统

冰蓄冷系统的核心是双工况制冷机组和蓄冰槽。在制冰模式下，制冷机在-5~-3℃的蒸发温度下运行，COP（性能系数）约为2.8-3.2；而在常规制冷模式下，蒸发温度升至2~4℃，COP可提升至4.5-5.3。实际项目中，我们通常选用乙二醇溶液作为载冷剂，其质量浓度控制在25%-30%之间，既能保证低温流动性，又避免过度黏滞。

蓄冰槽的数学模型需要考虑显热和潜热两部分。以某项目使用的 encapsulated ice storage为例，单个冰球（直径77mm）的潜热储存量约为175kJ，显热容量为1.8kJ/℃。Matlab建模时需要建立分段函数：

function Q = ice_storage(T, m_ice) if T >= 0 % 仅显热 Q = m_ice * 2.09 * (T - 0); % 水比热容2.09 kJ/(kg·K) else % 显热+潜热 Q = m_ice * (2.09*(0-T) + 334); end end

2.2 冷热电联供子系统

燃气内燃机是CCHP系统的核心，其热电关系可通过以下经验公式描述：

P_elec = η_elec * Q_fuel P_heat = η_heat * Q_fuel * (1 - η_elec)

其中η_elec通常为35%-45%，η_heat可达40%-50%。在Matlab中，我们用线性不等式约束表示这种耦合关系。某项目使用的颜巴赫JMS 612 GS-N.L机组，在75%负荷时实测电效率41.2%，热回收效率43.8%。

吸收式制冷机的性能用COP_abs表示，通常为0.7-1.2，其冷量输出与热输入的关系为：

Q_cool = COP_abs * Q_heat_in;

3. 多时间尺度优化框架

3.1 时间尺度划分

我们采用三级时间尺度架构：

• 长期调度层（24小时）：以1小时为间隔，决定机组启停和蓄冰策略
• 短期调度层（1小时）：以15分钟为间隔，调整设备出力
• 实时控制层（15分钟）：以1分钟为间隔，处理功率波动

在Matlab中，这通过分层优化实现。主问题（长期调度）使用intlinprog求解启停计划，子问题（短期调度）用fmincon进行连续变量优化。

3.2 目标函数构建

目标函数包含三项：

f = w1*Cost + w2*Emission + w3*Comfort;

其中：

• 经济成本Cost包括购电费用、燃气费用和维护成本
• 排放Emission折算CO2当量（电网电：0.85kg/kWh，燃气：0.19kg/kWh）
• 舒适度Comfort用PMV指标偏差的平方和表示

权重系数需根据项目需求调整。某办公楼项目采用w1=0.6, w2=0.3, w3=0.1，通过AHP（层次分析法）确定。

4. Matlab实现关键代码解析

4.1 混合整数规划建模

使用optimproblem建立问题框架：

prob = optimproblem('ObjectiveSense','minimize'); % 定义二进制变量表示设备状态 u_CHP = optimvar('u_CHP',24,'Type','integer','LowerBound',0,'UpperBound',1); % 添加约束 prob.Constraints.powerBalance = sum(P_gen) == sum(P_load) + sum(P_charge);

4.2 分段线性化处理

对于非线性设备特性，采用分段线性逼近。以燃气轮机为例：

% 划分5个运行区间 breakpoints = [0 0.3 0.5 0.8 1.0]*P_max; slopes = [0.28 0.32 0.35 0.38]; % 各段效率斜率 intercepts = [0 0.015 0.022 0.03]; % 截距

4.3 并行计算加速

利用parfor加速多场景计算：

parfor i = 1:num_scenarios [xopt(i), fval(i)] = solve(prob,'Options',options); end

5. 实际应用中的经验技巧

5.1 初始解生成策略

好的初始解能显著缩短求解时间。我们开发了基于规则的启发式方法：

1. 优先使用谷电时段（23:00-7:00）制冰
2. CCHP机组在电价高峰时段（10:00-12:00, 18:00-21:00）至少运行50%负荷
3. 吸收式制冷机在环境温度>28℃时优先启用

在Matlab中通过x0参数传递初始解：

options = optimoptions('intlinprog','InitialPoint',x0);

5.2 模型简化技巧

通过灵敏度分析确定可简化的部分。某项目中，我们发现蓄电池的充放电效率对总成本影响<0.5%，遂将其简化为固定效率（90%）模型，使求解时间从324秒降至187秒。

5.3 求解器参数调优

关键参数设置建议：

options = optimoptions('intlinprog',... 'RelativeGapTolerance',0.01,... % 允许1%的gap 'MaxTime',300,... % 限制求解时间 'Heuristics','advanced',... % 使用高级启发式 'CutGeneration','advanced'); % 生成更多切割平面

6. 典型问题排查指南

6.1 不可行解分析

当模型返回infeasible时，按以下步骤排查：

1. 检查功率平衡约束是否过紧（可暂时放宽5%测试）
2. 验证设备爬坡速率约束是否合理（特别是燃机启动阶段）
3. 使用irreducible inconsistent subsystem (IIS)分析：

[~,infeasibleConstraints] = iis(prob);

6.2 数值不稳定处理

遇到"Numerical instability"警告时：

1. 对变量进行归一化处理（如将功率单位从W改为MW）
2. 调整约束容差：

options.ConstraintTolerance = 1e-6;

3. 避免不同数量级系数混合（如将0.0001改为1e-4）

7. 性能优化实测数据

在某医院项目中（建筑面积8.2万㎡），不同调度策略对比：

求解性能：

• 模型规模：1,248变量（其中96个整数变量），892约束
• 求解时间：平均47秒（ThinkPad P15v, i7-11800H）
• 最优间隙：0.83%