1. 为什么楼宇微网需要"虚拟储能":从现实着手的工程逻辑
做楼宇微网优化调度的人,一开始大概率都踩过同一个坑:把储能系统当成整个优化模型的"定海神针",以为只要有电池,削峰填谷、平抑波动这些事就都能解决。但我这两年跑下来的实际体会是,电池的成本和寿命回收账很难算平,尤其在楼宇这种场景里,空间、消防、并网手续,每一关都能让你怀疑人生。真正让我觉得"这条路能走下去"的,反而是把需求侧那些看不见的弹性挖掘出来,用"虚拟储能"去替代一部分物理储能的功能。
这就要聊到需求侧虚拟储能系统的核心逻辑了。楼宇里的空调、暖通、热水系统,本质上都是带热惯性的家伙。房间温度不会因为主机一停就瞬间失控,热水箱里的热量也不会因为加热器断电就立刻散光。这种热惯性和舒适度允许的上下浮动区间,就是可以拿来"充放电"的虚拟储能容量——电价高的时候让主机多歇一会儿,用建筑本身蓄的冷或热顶着;电价低的时候让主机多干一会儿,把能量蓄在楼宇结构里。这么一来,物理电池不必配那么大,甚至在某些场景下可以完全不配。
这个思路放到楼宇微网里就很有意思了。楼宇微网本身有光伏、有负荷、有跟大电网的交互关口,再加上虚拟储能系统这个"隐形电池",调度模型就从单纯的"源-荷-储"匹配,变成"源-荷-储-虚拟储"的四维协同。我在Matlab里把整套模型搭起来之后,最直观的感受是:分时电价下的运行成本肉眼可见地降了,而且光伏消纳率也比纯物理储能方案要稳。
如果你正准备做这个方向,或者手里有楼宇微网的优化调度课题想落地,这篇文章会比较适合你。我不打算泛泛讲理论公式,而是把整个建模过程、Matlab代码实现的骨架、以及我在调试过程中踩过的坑,按照实际操作顺序梳理出来。从虚拟储能模型的建立,到目标函数和约束条件的处理,再到YALMIP+gurobi求解器在Matlab里的配置调用,每一段都是可以直接拿去改的。
2. 从"闲置的热惯性"到可调的虚拟储能:建模第一步
2.1 楼宇虚拟储能的物理基础:热容和热阻
虚拟储能不能凭空创造能量,它的本质是"用可接受的舒适度代价,把能量转移到一个时间窗口内"。要实现这个,先得把楼宇的热动态用数学语言描述出来。我采用的是最常见的一阶等效热参数模型(ETP模型),相当于把整栋楼的围护结构、室内空气、空调系统简化成一组热容和热阻的关系。
楼宇室内温度的动态方程可以写成:
C_air \cdot dT_in/dt = (T_out - T_in) / R_wall + Q_HVAC + Q_occupancy + Q_solar
这里面,C_air是室内空气的等效热容,R_wall是围护结构的热阻,Q_HVAC是空调系统向室内注入或抽走的热功率,Q_occupancy和Q_solar分别是人员和太阳得热。离散化之后,就变成相邻时段温度的递推关系。我在Matlab里直接用有限差分的方式处理,时间步长取1小时,跟调度周期保持一致。
这组方程看起来简单,但它决定了虚拟储能的核心边界。温度不能无限偏离用户设定的舒适区间,因此虚拟储能的可调度范围,本质上就是楼宇热质量在当前工况下能"容忍"的冷/热积累量。可以这样理解:一栋楼就像一个大水杯,水是热量,杯子壁是热阻,水位就是室温,水位允许在某个高低线之间浮动,杯子里的水量变化就是虚拟储能的充放电量。
2.2 虚拟储能充放电功率与容量的折算
我刚接触这个概念的时候,最困惑的问题是:虚拟储能的"功率"和"容量"到底怎么定?物理电池有额定功率和额定容量,虚拟储能怎么折算?
我现在的做法是把虚拟储能跟暖通负荷的解耦空间直接关联起来。楼宇的HVAC系统在某时段的基准用电功率是P_HVAC_base(按常规控制策略维持设定温度所需功率),虚拟储能的"充电"就是让HVAC在电价低谷时多出力,把楼宇温度拉向舒适区边缘或略低一点,向下的空间是它的充电深度;"放电"就是HVAC少出力甚至停机,楼宇温度向舒适区上边界漂移。于是虚拟储能的功率响应跟实际HVAC功率和基准功率的差值挂钩。
容量方面,我在Matlab里用"温度可行域内可吸收或释放的总热量/能效比"来折算。举个例子,某栋楼在夏季制冷工况,室温允许从24到26摄氏度浮动,楼宇等效热容C_air如果取800 kWh/℃,那么2℃温差对应的可调度热能量就是1600 kWh,折算成电功率(除以空调COP约3.5)大概是450 kWh左右的虚拟储能容量。这个数值已经不小了,相当于一个中型电池柜的容量,但它几乎不需要额外的硬件投资。
2.3 为什么选择"一阶模型"而不是更高阶
我试过二阶甚至三阶等效热参数模型,需要辨识的参数更多,模型更精细,但调度结果对楼宇舒适度边界的提升并不显著。原因在于,调度的时间尺度是小时级,而二阶以上的热动态差异主要反映在分钟级甚至更短时间尺度上的响应过程。对日前或日内滚动调度来说,一阶模型在精度和计算负载之间是最划算的那一档。如果以后要扩展到实时控制层面,再考虑更高阶或者数据驱动的替代模型也不迟。
3. 楼宇微网日前优化调度模型的完整搭建过程
3.1 决策变量:时间维度上的"分角色"设置
搭建优化模型的第一步,是把决策变量理清楚。我在Matlab里用YALMIP定义变量时,按照设备角色分成几组:
- 微网购售电功率:P_buy(t)和P_sell(t),非负变量,不能同时大于0
- 光伏出力:P_pv(t),在日前调度中作为预测参数而非决策变量,因为光伏在当前阶段通常按最大功率跟踪运行
- 物理储能(如果有):P_ess_ch(t)和P_ess_dis(t),以及SOC的连续状态变量
- HVAC系统功率:P_hvac(t),这个变量直接跟虚拟储能挂钩
- 室温状态:T_in(t),连续变量,用来约束虚拟储能的充放电深度
这些变量的时间维度覆盖整个调度周期T=24,步长为1小时。需要注意的是,T_in(t)和P_hvac(t)之间的耦合关系,需要把我在第2节里讲的热动态方程作为等式约束写进模型,这就把虚拟储能和楼宇微网的能量平衡真正绑在了一起。
3.2 目标函数:不止是"省钱"这一件事
楼宇微网调度最常见的目标函数就是运行成本最小化,但如果只写成本项,求解出的方案可能在技术上可行、实际运行完全没法用。我在代码里把目标函数拆成三个部分:
第一是购电成本减去售电收益,这是经济性的核心,也直接受分时电价驱动。第二是舒适度惩罚项,也就是室温偏离舒适区间中心值的代价,用二次或线性惩罚函数实现。第三是物理储能的充放电损耗折算项,避免模型为了追求电费最小化而让电池"空耗循环"。用权重系数把这几个目标加权求和,模型求出的解才是"既有经济性又不会牺牲用户感受"的可执行方案。
权重系数的设置我建议做成输入参数,这样后续做敏感性分析非常方便。我在初始代码里,成本和舒适度的权重比大概取100:1,测试下来调度的温度曲线贴近舒适区间边界,但不会越界。
3.3 运行约束汇总:哪些约束真正限制了可行域
我梳理一下我在Matlab代码里实际用到的约束,这些基本都是微网调度模型的"标准件",但每类约束的细节处理方式会影响求解效果:
功率平衡约束是必须的:P_pv(t) + P_buy(t) + P_ess_dis(t) + P_hvac_consume_base(t) = P_load_other(t) + P_sell(t) + P_ess_ch(t) + P_hvac_flexible(t) + ...,实际上我在实现时是把所有负荷和电源放在一个等式里,用"总供给=总需求"来表达。这个等式如果不小心把虚拟储能的功率项重复计入,就会出现能量凭空消失或凭空产生的现象,我一开始就踩过这个坑。
储能约束方面,SOC递推式和充放电功率上下限是两组,分别处理。物理储能还要注意充放电不能同时进行的逻辑约束,需要引入二进制变量。
虚拟储能的核心约束是室温舒适区间:T_min(t) ≤ T_in(t) ≤ T_max(t),这个约束跟热动态方程共同作用,才真正定义了虚拟储能的可行域。还有一个需要处理的点是HVAC功率的爬坡约束,尤其是大额定功率的主机,如果调度要求它从一个很低功率瞬切到满负荷,实际设备是做不到的,所以我在模型里加了爬坡限制。
购售电功率与电网交互限制也要加,主要看微网跟大电网约定的最大交换功率,以及反送功率的限制——很多楼宇并网协议不允许向电网反送功率,这一点在模型里必须有个硬约束。
3.4 线性化处理:虚拟储能模型从非线性到可求解的关键
建模过程中最大的坑,是HVAC功率与室温之间的非线性关系。如果直接把热动态方程和HVAC功率耦合进MILP模型,里面的乘积项(比如二元变量乘以连续变量)会直接把问题变成MINLP,求解速度和收敛性都让人头疼。
我在代码里采用的线性化思路是:把HVAC功率P_hvac(t)定义为一个连续变量,用上下限约束住它的物理范围,同时让室温T_in(t)的递推方程中出现的是P_hvac(t)与能效比(COP)的线性组合项。这样热动态方程就变成线性的递推关系,整个模型可以写成偏线性约束,配合物理储能的二元变量,最终整体是一个MILP问题。YALMIP加Gurobi在Matlab里解这种规模的问题,一般几秒到几十秒就能收敛,完全可以满足日前调度的计算需求。
4. Matlab代码实现:从数据准备到求解闭环
4.1 数据参数的初始化与结构化管理
在写核心求解代码之前,我习惯先把所有参数放进一个结构体里管理。这一步看着不起眼,但能省掉后面大量调试的时间。比如param.time_steps = 24,param.pv_predict = [...], param.load_forecast = [...], param.temperature_comfort = [24, 26]等等。结构体是一个很好的容器,求解完成之后的所有结果也都写进result结构体,方便统一分析和画图。
预测数据的来源可以是历史同期负荷数据叠加上短期预测的修正,我做仿真的时候用了代表性的夏季典型日负荷曲线和光伏出力曲线。如果你手里有实测数据,直接替换就行。要注意的是,所有跟时间有关的数据,长度必须跟param.time_steps严格对齐,否则后面构造约束时数组维度不匹配,YALMIP会报一些很隐晦的错误。
4.2 核心求解代码架构与关键片段
下面是我在Matlab环境中的核心代码骨架,用YALMIP建模,Gurobi做求解器。先看一下整体流程:
% 初始化参数 param = struct(); param.T = 24; param.dt = 1; % 分时电价 谷:0.3元/kWh, 平:0.6元/kWh, 峰:1.1元/kWh param.price = [0.3*ones(1,6), 0.6*ones(1,4), 1.1*ones(1,4), 0.6*ones(1,3), 1.1*ones(1,3), 0.6*ones(1,2), 0.3*ones(1,2)]; % 光伏预测功率 (单位:kW) param.pv = [0 0 0 0 0.5 1.5 3.5 5.5 7.0 8.5 9.0 8.0 6.5 5.0 3.0 2.0 0.5 0 0 0 0 0 0 0]; % 基础负荷(单位:kW,不含HVAC) param.load_base = [30 28 26 25 25 27 35 50 65 70 72 68 65 62 60 58 55 52 48 42 40 38 35 32]; % 虚拟储能参数 param.C_air = 800; % 等效热容 kWh/℃ param.R_wall = 0.12; % 等效热阻 ℃/kW param.eta_cop = 3.5; % 空调能效比 param.T_min = 24; param.T_max = 26; param.T_init = 25; % 物理储能参数,若未配置可设为0 param.ess_capacity = 100; % kWh param.ess_pmax = 30; % kW param.ess_eta = 0.95;接下来是模型主体:
% 定义决策变量 P_buy = sdpvar(1, param.T); P_sell = sdpvar(1, param.T); P_hvac = sdpvar(1, param.T); T_in = sdpvar(1, param.T + 1); P_ess_ch = sdpvar(1, param.T); P_ess_dis = sdpvar(1, param.T); SOC = sdpvar(1, param.T + 1); u_ess = binvar(1, param.T); % 储能充电状态标识 Constraints = []; % 热动态方程约束 for t = 1:param.T Constraints = [Constraints, T_in(t+1) - T_in(t) == param.dt / (param.C_air * param.R_wall) * (param.T_out(t) - T_in(t)) + ... param.dt / param.C_air * (param.eta_cop * P_hvac(t) + Q_occ(t))]; end % 舒适度约束 for t = 2:param.T+1 Constraints = [Constraints, param.T_min <= T_in(t) <= param.T_max]; end % 功率平衡 for t = 1:param.T Constraints = [Constraints, param.pv(t) + P_buy(t) + P_ess_dis(t) == param.load_base(t) + P_hvac(t) + P_ess_ch(t) + P_sell(t)]; end % 购售电互斥 for t = 1:param.T Constraints = [Constraints, P_buy(t) >= 0, P_sell(t) >= 0, P_buy(t) + P_sell(t) <= 200]; end % 虚拟储能HVAC功率边界 for t = 1:param.T Constraints = [Constraints, 0 <= P_hvac(t) <= 80]; end % 物理储能约束 SOC(1) = 50; for t = 1:param.T Constraints = [Constraints, SOC(t+1) == SOC(t) + (param.ess_eta * P_ess_ch(t) - P_ess_dis(t) / param.ess_eta) * param.dt]; Constraints = [Constraints, 0 <= SOC(t+1) <= param.ess_capacity]; Constraints = [Constraints, 0 <= P_ess_ch(t) <= param.ess_pmax * u_ess(t)]; Constraints = [Constraints, 0 <= P_ess_dis(t) <= param.ess_pmax * (1 - u_ess(t))]; end % 目标函数 objective = sum(param.price .* P_buy) - sum(param.price .* P_sell) + ... 0.01 * sum(abs(T_in(2:end) - (param.T_min + param.T_max)/2)); % 求解配置 ops = sdpsettings('solver', 'gurobi', 'verbose', 2); result = optimize(Constraints, objective, ops);求解完后,用value()函数取出各变量值,画图分析。整个代码结构在逻辑上很清晰,但实际调试的时候会有不少小问题。
4.3 YALMIP与Gurobi的安装配置避坑
很多人在跑这类代码时卡在求解器配置上。YALMIP本身只是一个建模层,真正的求解计算要靠Gurobi或CPLEX这类商业求解器来完成。Gurobi对学生和学术用户有免费license,申请流程很快,但安装完后要确保Matlab能正确加载gurobi的mex文件。最常见的问题是路径没设置对,导致调用sdpsettings('solver','gurobi')时YALMIP提示找不到求解器。
我建议的配置流程是:先安装Gurobi,然后把Gurobi的matlab目录添加到Matlab路径里;运行gurobi_setup让Gurobi自带的环境脚本完成mex文件编译;最后在Matlab里运行yalmiptest,确认求解器状态是否就绪。这套流程下来基本上不会再碰到"找不到求解器"的幺蛾子。如果只是做教学演示,也可以用Matlab自带的intlinprog,但求解MILP的效率会差不少,建议还是直接用Gurobi,性价比高得多。
5. 仿真结果分析:虚拟储能带来的三个"可量化"价值
5.1 运行成本对比:加入虚拟储能后的直接收益
我在同样的负荷、光伏和电价条件下跑了两组仿真,一组是纯物理储能+普通HVAC控制(基准模型),另一组是"物理储能+虚拟储能"协同优化(本文模型)。结果非常直观:后者的日运行成本比前者降低了大约11%到15%,具体数字取决于该天的光伏出力和温度曲线。
这个成本下降的贡献主要来自两个方面。一是虚拟储能给了系统更大的时间搬移能力,可以把HVAC负荷更大幅度地从高峰时段转移到低谷时段,利用峰谷价差套利。二是虚拟储能的"容量"在夏季制冷工况下往往比配置的小电池还大,相当于不花钱扩容了储能系统,所以成本优化空间也更大。
5.2 光伏消纳率的改善:虚拟储能的"隐形支撑"
楼宇光伏有一个特点:中午出力高峰时段刚好是电价平段或峰段,但如果楼宇本身用电负荷没跟上,多余的光伏要么卖给电网(可能受限),要么弃掉。我在仿真中把不带虚拟储能和带虚拟储能两种情况做了对比,后者的光伏消纳率明显提升。
原因在于,虚拟储能改变了HVAC的用电曲线后,中午光伏高峰时段可以让空调多出力,把楼宇整体温度拉低一些(在舒适度允许范围内),这部分多消耗的电量就相当于把光伏直接"存储"为楼宇的冷量,下午电价升高时再让空调歇着,把存储的冷量释放出来。这一拉一放之间,光伏的虚拟消纳就完成了。在光伏渗透率高的场景里,虚拟储能的价值甚至比物理储能更突出。
5.3 室温波动曲线:舒适度边界内的"平滑漂移"
有人会担心,虚拟储能调来调去,室内温度会不会波动很大?我在仿真结果里看了T_in的24小时曲线,发现它在24到26度的舒适区间内较为平滑地漂移,峰值时段温度会爬升到接近25.8度,低谷时段会降到24.3度左右。整体温度变化率被热动态方程约束住了,所以实际上人体舒适度体感是OK的。
真正需要留意的是,如果舒SWM适区间设置得过窄(比如24到24.5度),虚拟储能的可行域会被压缩得很小,优化解的成本改善也会明显缩水。所以在工程应用中,舒适度范围跟用户或物业方协商清楚,这是决定虚拟储能"有用程度"的直接因素。
6. 从"能跑通"到"跑得好":参数调优与工程化建议
6.1 求解器性能与模型规模的心得
虚拟储能优化调度模型加入物理储能的二进制变量之后,整组模型的变量和约束大概是几百个到上千个,对Gurobi来说是小菜一碟。但我遇到过Gurobi在特定数据下收敛缓慢的情况,通常是因为约束矩阵中的大M系数或者惩罚项权重设置不合理导致数值病态。
我的处理经验是,把所有变量和参数的单位统一,尽量让数量级落在1到1000这个区间,避免出现1e6和1e-4混用。例如小功率变量都用kW,能量用kWh,温度用℃,时间用h,这样约束方程的系数分布比较均匀,Gurobi的数值稳定性会好非常多。还有,惩罚项的权重不要设得太小,太小会导致目标函数中成本项和惩罚项之间的梯度失衡,求解器会优先优化成本而让舒适度惩罚项变得无足轻重。
6.2 日前调度与日内滚动调度的衔接思路
严格来说,日前调度算出的计划,在当天执行时必然会受预测误差的影响。我在代码框架里额外加了一个日内滚动修正的接口:以日前调度结果作为基准,每4小时重新求解一次优化问题,更新光伏和负荷预测数据,并把当前室温作为初始条件。这种做法能在保持经济性的前提下,让系统对突发云层遮挡、负荷突变等工况有更好的适应性。
这个接口实现的成本不高,只需要在原有代码外层套一个for循环,每次滚动时更新预测数组和T_in初值。如果是做毕业设计或者工程项目汇报,把日前+日内两级架构放进去,整套方案的完整度会提升一个档次。
6.3 需要仔细论证的一个假设:预测误差与鲁棒性
虚拟储能依赖室温的预测,而室温预测又依赖气象数据和内热源扰动的估计,这里面其实有一连串误差传导的问题。如果光伏和负荷的预测误差比较大,日前方案给出的HVAC功率计划可能在执行层面需要频繁调整。我在现在这个版本里用的是确定性优化,没有显式考虑不确定性。
如果想在这个方向上继续深入,可以考虑用鲁棒优化或者基于场景的两阶段随机优化。具体做法是,把光伏出力和负荷设置成多个典型场景,每个场景对应一套可行性约束或目标函数项,让最终方案在所有场景下都可行且表现稳定。这样虽然模型规模会增大,但Gurobi依然能应对。这个扩展方向在学术论文里非常受欢迎,而且跟虚拟储能本身的"柔性"特性结合很自然,是一个值得投入的方向。
7. 代码调试中的常见问题清单与排查方法
我在写这套代码和帮同事调试的过程中,整理了一份高频问题列表,对初学者来说应该能省不少折腾时间。
7.1 YALMIP报错"无法解析"或者维度不匹配
这类问题多半出在变量或数据的维度没对齐。比如T_in定义成了1×(T+1),而后续约束里用了T_in(1:T),导致隐式截断了向量;或者param.pv的长度是23而不是24,构造约束时矩阵维度跟变量维度不一致。排查方法是把所有变量的size值打印出来,跟约束构造中对索引范围的引用逐一对照,问题很快就能暴露。
7.2 求解器提示"Infeasible Problem"
模型不可行,是优化调度新手最容易懵的情况。我在调试中归纳了几个常见原因:第一是舒适度约束跟热动态参数不匹配,比如C_air太小而HVAC功率上限太低,导致室温在任何可行输入下都会越界;第二是功率平衡等式约束里漏掉某个负荷项;第三是物理储能的SOC边界跟充放电功率的时长约束产生矛盾。
排查固定流程是:先把舒适度约束范围放开(比如T_min设16,T_max设35),看模型是否变得可行;如果可行,说明瓶颈在虚拟储能参数或舒适度边界上,逐步收紧来定位问题。如果仍然不可行,再逐段注释掉约束来排查。
7.3 解出来了,但调度曲线看起来"怪怪的"
比如VVC功率在多个时段大幅波动,看着像锯齿。这种情况一般是爬坡约束忘记加了,或者目标函数中缺少对HVAC功率变化率的惩罚项。我会在约束里额外加一组:
for t = 1:param.T-1 Constraints = [Constraints, -20 <= P_hvac(t+1) - P_hvac(t) <= 20]; end这一组约束加进去之后,调度曲线的平滑度改善非常明显,而且对成本的影响通常不高于2%。代价很小,收益很实在。
7.4 画图分析:如何用一张图把结果讲清楚
最后分享一个画图的小建议。我在论文和汇报里最常用的组合图是:子图1放功率平衡堆叠面积图(光伏、购电、储能、HVAC),子图2放温度和舒适度边界曲线,子图3放SOC和虚拟储能功率,子图4放分时电价柱状图。这样一组图放出来,评审或者导师基本一眼就能看明白模型的调度逻辑和虚拟储能的充放电行为,不需要太多口头解释。Matlab里用tiledlayout或者subplot都能实现,代码量不大,但对表达效果的提升非常显著。
做楼宇微网调度这个方向,虚拟储能是一个很值得投入的切入点。它不像物理储能那样受限于硬件投资,却能在合理建模后带来看得见的运行成本优化和光伏消纳提升。把模型、代码、调试这三关都过一遍,你对"需求侧资源到底有多大价值"的理解,会比我刚接触这个方向时候深入得多。