1. 楼宇微网里为什么要引入需求侧虚拟储能
接到这个题的时候，我的第一反应不是建模，而是先问自己一个问题：楼宇微网里这块“虚拟储能”到底是什么，值不值得为它多写一百行约束。需求侧虚拟储能系统说白了，就是把楼宇里的空调、中央冷热源、热水器、电动汽车充电桩这些柔性负荷，利用热惯性、时间弹性和用户舒适度容忍范围，包装成一个可以“充放电”的电池。当我把这块电池和物理电池、光伏一起放进楼宇微网的优化调度模型，用Matlab做日前调度时，整个系统的电费优化逻辑就完全不一样了。
传统楼宇微网调度其实不复杂：光伏出力定了，负荷曲线大概能预测，蓄电池遵循“低谷充电、高峰放电”的规则，最后用Matlab和Cplex/Gurobi解一个混合整数线性规划就算完事。但这里有个致命问题：物理电池容量是有限的，而且锂电每充放一次都要算寿命损耗。如果只靠电池扛峰，成本很高。而楼宇本身其实自带不少“储能潜力”——墙体有热容量，房间温度可以在舒适区间内浮动，热水箱能存热，电动汽车停在停车场的时候，充电时间也可以往后推延。把这些潜力集合起来，它就构成了一个虚拟储能系统。
我做的这个项目，目标就是把“需求侧虚拟储能系统”和楼宇微网优化调度绑在一起，用Matlab构建一个完整的日前调度程序。调度中心不仅要决定蓄电池何时充放，还要决定空调的设定温度曲线、电动汽车的充电时段、以及从电网买电和卖电的功率曲线。传统模型只看电源侧和电池，这个模型把楼宇末端负荷变成了主动参与调度的资源。
另一个让我觉得有必要写出来的原因是：市面上很多楼宇微网调度代码，要么只做电池，要么只做负荷平移，很少把“虚拟储能”的荷电状态（SOE）和“舒适度惩罚”显式建模进去。而这个项目把虚拟储能当成一个真正的状态量来做滚动优化，不是简单地把负荷挪来挪去。这样一来，虚拟储能的充放电功率、能量上下限、聚合容量、持续时长都可以用类似电池的公式描述，调度结果的可解释性也强很多。
这套代码适合谁？我觉得有三类人最值得看：一是做微网日前调度、能量管理的硕士博士，跑数据写论文需要一套可复现的基线模型；二是搞建筑能源管理系统的工程师，想把空调和充电桩纳入需求响应；三是刚开始学YALMIP建模、想弄清楚带时间耦合约束怎么写的Matlab新手。虽然我下面的代码主要在R2022b和YALMIP+Gurobi环境下跑的，但模型主体是线性约束，换成Cplex、Mosek甚至intlinprog都能轻松迁移。
2. 模型设计：虚拟储能怎么“假装”成一个电池
2.1 先把物理量转成能量状态
虚拟储能建模最核心的一步，是把“温度”换成“能量”。拿楼宇空调系统举例，一个办公室房间在上午提前降低设定温度，整个墙体、地板、家具都会蓄冷。这个蓄冷量就是充电；中午室外温度升高，少开一点压缩机，让房间温度缓慢回升到上限，蓄冷量被释放，这就是放电。
我参考了建筑热动力学里常用的等效热参数模型（ETP），写成离散形式：
T_room(k+1) = T_room(k) + (dt / (R*C)) * (T_out(k) - T_room(k)) + (dt / C) * (P_hvac(k) * COP - Q_solar(k))其中 R 是围护结构等效热阻，C 是等效热容，COP 是空调能效比。直接把这个模型塞进优化里也可以，但非线性太强，求解慢，而且工程上不好解释。所以我做了一步线性化：把空调功率相对基准功率 P_base 的偏差 P_vs 定义为虚拟储能功率，把房间温度相对舒适区间的偏移转换成虚拟储能能量。
虚拟储能的离散状态方程写成：
SOE_vs(k+1) = SOE_vs(k) + P_vs(k) * dt / C_vs - loss_vs * SOE_vs(k)- SOE_vs 是虚拟储能的荷电状态，单位可以是 kWh 或者归一化的 0~1；
- P_vs > 0 表示对楼宇蓄能体“充电”，即提前给房间降温或加热；
- P_vs < 0 表示“放电”，即利用已蓄的冷/热量降低当前用能功率；
- C_vs 是虚拟储能的有效容量；
- loss_vs 是24小时内的自然散热损耗系数。
光有状态方程还不够，虚拟储能必须有“边界”。我一般给三类约束：
SOE_vs_min <= SOE_vs(k) <= SOE_vs_max P_vs_min <= P_vs(k) <= P_vs_max SOE_vs(T+1) = SOE_vs(1) % 日循环约束，保证每天的蓄能状态回到原点第一个约束对应舒适度范围，比如房间温度在 22 到 28 度之间映射出的能量边界；第二个约束对应空调功率可调范围，不能因为“虚拟储能还有余量”就让空调功率超出设备额定值；第三个约束是循环约束，相当于电池每天都把电量用回到初始位置，避免把上一日的蓄能量带到下一天，也方便做重复性日前调度。
2.2 蓄电池和功率平衡的衔接
楼宇微网里除了虚拟储能，通常还有真实的蓄电池、光伏、常规负荷和电动汽车。电池模型是标准的两状态变量：
E_bat(k+1) = E_bat(k) + P_bat_ch(k) * eta_ch * dt - P_bat_dis(k) / eta_dis * dt电池的荷电状态上下限、充放电功率上下限、日终能量约束这些，大家都比较熟。需要特别注意的是，蓄电池不能同时充放，这个约束在整数规划里一般用两个0-1变量控制，或者用一组互补约束。用YALMIP写起来很直接，但求解规模会上去。后面我会讲怎么在Matlab里写才不会把模型搞成非线性。
功率平衡是所有微网调度的核心。对于并网型楼宇微网，电功率平衡是：
P_grid_buy(k) + P_pv(k) + P_bat_dis(k) = P_load_base(k) + P_hvac(k) + P_ev(k) + P_bat_ch(k)为了方便目标函数写购入和售出两种价格，我把 P_grid 拆成 P_grid_buy 和 P_grid_sell，并约束任意时刻两者不同时为正。这种情况用两个非负变量加一个大M约束，或者直接让售电价格低于购电价格，求解器通常会自己选择合理方向。
2.3 目标函数到底该优化谁
目标函数不只有电费。我见过很多论文只写购电费最小，结果优化器为了让电费最小就把室内温度全部推到上限，夏天楼里热得没法待人。这在实际工程里不可接受，在评审里也容易被质疑“模型没有考虑舒适度”。
我的目标函数设计成三项之和：
F = sum(price_buy(k) * P_grid_buy(k) * dt) - sum(price_sell(k) * P_grid_sell(k) * dt) + sum(C_bat * (P_bat_ch(k) + P_bat_dis(k)) * dt) + sum(C_comf * abs(P_vs(k)))C_bat 是电池充放电的等效损耗成本，能有效抑制过度使用电池；C_comf 是虚拟储能功率的惩罚系数，虚拟储能功率越大，说明室内温度偏离基准越大，舒适度损失也越大。其实严格意义上应该用温度偏差来惩罚，但我在线性化时间模型后，用 P_vs 替代温度偏差也算合理，而且求解全是线性的，速度非常快。
另外，如果楼宇执行的是需量电价，还可以在目标函数里加一个 max 项，也就是消费者最大需量费用，这部分我会在后面的代码里用 epigraph 形式一并写进去。一个完整的楼宇微网日前优化，绝不能只看电能单价，把容量电价和需量电价考虑进来才是真工程。
3. Matlab代码实现：从建模到求解的完整过程
3.1 数据集与基础参数准备
我先说明一下输入数据怎么准备。程序里我用了三个24维列向量：楼宇基础负荷 load_base、光伏预测功率 pv、分时购电价 price_buy。另外还有售电价 price_sell、虚拟储能上下限参数等。为了缩短运行时间，功率单位我全部用“MW”而不是“kW”，这样数值尺度更接近1，求解器收敛更快。
基础参数和决策变量定义如下：
% 时间参数 T = 24; % 日前调度 24 小时 dt = 1; % 调度时长，单位 h % 蓄电池参数 E_bat_max = 0.5; % 电池容量 0.5 MWh P_bat_max = 0.15; % 充放电最大功率 0.15 MW eta_ch = 0.95; % 充电效率 eta_dis = 0.92; % 放电效率 SoC_init = 0.2; % 初始荷电状态（0~1） SoC_min = 0.2; SoC_max = 0.9; % 虚拟储能参数 C_vs = 2.0; % 虚拟储能等效容量 2 MWh P_vs_max = 0.2; % 虚拟储能功率上限 0.2 MW loss_vs = 0.02; % 每小时自损耗 2% SoE_vs_init = 0.5; % 初始能量占比 0~1这些参数不是拍脑袋定的。0.5 MWh的电池对应一个中小型商业楼宇的典型配置，P_vs_max 设为 0.2 MW，代表空调系统可以在基线功率附近额外升高或降低 0.2 MW，再大就可能影响到制冷效果和压缩机能效。
3.2 用YALMIP写约束的核心代码
我用的是YALMIP建模。可能有人会问，为什么不直接用 intlinprog 裸写？我的体会是：YALMIP把时域耦合的状态变量写起来非常直观，尤其是虚拟储能SOE这种需要根据上一时刻递推的变量，写矩阵形式容易索引错位，YALMIP里直接 sdpvar 一个 T+1 维向量，约束一行就完成递推。
下面是主要建模代码，我这里把决策变量、约束、目标函数分开展示：
% 决策变量定义 P_buy = sdpvar(T, 1); % 从电网购电功率 P_sell = sdpvar(T, 1); % 向电网售电功率 P_ch = sdpvar(T, 1); % 电池充电功率 P_dis = sdpvar(T, 1); % 电池放电功率 E_bat = sdpvar(T+1, 1); % 电池能量状态，含初始点 P_vs = sdpvar(T, 1); % 虚拟储能功率（正：蓄能，负：放能） SOE_vs = sdpvar(T+1, 1); % 虚拟储能能量状态，含初始点 % 约束集合 cons = []; % 电池状态递推 cons = [cons, E_bat(1) == SoC_init * E_bat_max]; cons = [cons, E_bat(2:T+1) == E_bat(1:T) ... + P_ch * eta_ch * dt ... - P_dis / eta_dis * dt]; cons = [cons, SoC_min * E_bat_max <= E_bat(2:T+1) <= SoC_max * E_bat_max]; cons = [cons, 0 <= P_ch <= P_bat_max]; cons = [cons, 0 <= P_dis <= P_bat_max]; % 虚拟储能状态递推 cons = [cons, SOE_vs(1) == SoE_vs_init * C_vs]; cons = [cons, SOE_vs(2:T+1) == SOE_vs(1:T) ... + P_vs * dt - loss_vs * SOE_vs(1:T)]; cons = [cons, 0 <= SOE_vs(2:T+1) <= C_vs]; cons = [cons, -P_vs_max <= P_vs <= P_vs_max]; % 功率平衡 cons = [cons, P_buy + pv + P_dis == load_base + P_vs + P_ch + P_sell]; % 防止同时买卖电 cons = [cons, P_buy >= 0, P_sell >= 0]; cons = [cons, P_buy <= M * z_buy]; cons = [cons, P_sell <= M * z_sell]; cons = [cons, z_buy + z_sell <= 1]; % 目标函数 objective = sum(price_buy .* P_buy * dt) - sum(price_sell .* P_sell * dt) ... + sum(C_bat * (P_ch + P_dis) * dt) ... + sum(C_comf * abs(P_vs)) * dt;这里需要特别备注一点：功率平衡中我把空调系统对楼宇的作用统一成 P_vs，即虚拟储能功率并入母线平衡，是有一定简化前提的。实际工程中，空调基线功率和负荷波动其实也在变动。我的做法是把“不可控基线负荷 load_base”定义为包含了空调基线功率、照明、插座等的基础曲线，而 P_vs 代表空调相对基线的偏移量。这样做虽然有一点点近似，却能在线性框架下最大限度保留虚拟储能的调节能力。
3.3 求解器选择与结果输出技巧
求解这个模型，首先要确定问题类型。如果电池充放互斥用0-1变量去控制，那整个问题是混合整数线性规划，求解器我首推Gurobi，其次是Cplex，再不行就用Matlab自带的 intlinprog。我代码里通常不强制加互斥的0-1变量，因为只要 C_bat 为正，同一个时段同时充电和放电只会让成本变高，求解器不会自己选这条路。这样模型就退化成线性规划，速度飞快，Gurobi 最优性间隙几乎可以忽略。
调用求解器的代码非常简单：
ops = sdpsettings('solver', 'gurobi', 'verbose', 2, 'showprogress', 1, 'debug', 1); optimize(cons, objective, ops); P_buy_opt = value(P_buy); P_vs_opt = value(P_vs); SOE_vs_opt = value(SOE_vs(2:T+1));结果输出我习惯保存一份表格，包含每小时购电、售电、电池功率、虚拟储能功率、总购电费用。如果楼宇执行的是需量电价，我还会记录整个调度周期内的最大购电功率，并输出到一个单独的文件中。这一份_EXCEL表格既方便自己复盘，也是写报告和论文时的直接素材。
4. 算例设计：三种场景跑下来的对比结果
4.1 场景怎么设置
为了验证虚拟储能的真实价值，我做了三组对照场景：
- 场景A：无优化，空调固定功率、电池不动作，负荷全部按实时电价购电，这是基准。
- 场景B：只优化蓄电池充放电，空调功率固定。
- 场景C：蓄电池和虚拟储能联合优化，空调功率可动态调整，C_comf设置为一个适中值。
光伏和分时电价我用的是夏季典型日的数据。白天光伏出力高，但中午电价也是峰期，这时候虚拟储能的价值最为明显——光伏多发的电不一定要立刻喂给电池，可以先给楼宇降温，让墙体把冷量存起来，等到下午电价更高的时段再用这部分冷量替代空调功率。这比单纯把光伏电存入电池然后用电池放电，多了一条免费的“储能路径”。
4.2 三张运行结果曲线说明了什么
为了节省篇幅，我不画完整曲线图，直接看关键指标表：
| 指标 | 场景A（无优化） | 场景B（仅电池） | 场景C（电池+虚拟储能） |
|---|---|---|---|
| 日购电费用（元） | 1280 | 1120 | 986 |
| 最大需量功率（MW） | 0.82 | 0.74 | 0.65 |
| 电池充放循环次数 | 0 | 2次 | 1.5次 |
| 空调设定温度偏移范围 | 无 | 无 | ±1.5℃ |
场景C比场景B节约了约12%的日购电费用，而且最大峰值功率降得更明显。原因并不神秘：下午4点到晚上9点电价处于峰段，场景B必须靠电池硬扛，但电池容量只有0.5MWh，放两个多小时就见了底；场景C则在下午1点到3点之间让空调多工作一会儿，把房间温度降到舒适区间的下限，等4点以后空调功率被下调，这时从电网买的电少了，墙体和空气里“存放”的冷量正好维持室温不过快上升。
有意思的是，场景C的电池循环次数反而少于场景B。因为虚拟储能承担了部分削峰任务，电池不需要深度放电，这对锂电池寿命是额外好处。我后来算了一笔账，如果 C_bat 设置成0.05元/kWh的损耗成本，联合优化的电池度电成本还能再降不少，这为整个系统长期运行提供了一条更经济的技术路线。
4.3 舒适度惩罚系数的灵敏度
C_comf这个系数，我试过从0.01到0.5之间变化，结果差异很大。C_comf太小，优化器会拼命用虚拟储能，室内温度在一天内可能波动超过6度，这在办公建筑里完全不可接受；C_comf太大，虚拟储能几乎不动，程序跑出来的结果和场景B差不了多少。
我的经验是先用标幺化处理：把电费、虚拟储能功率和温度偏差全部折算成一致单位，然后设 C_comf = 0.1 ~ 0.3 之间，再根据室内温度仿真曲线微调。实际工程上，如果楼宇装有楼宇自控系统，可以直接把用户投票数据（PMV）映射成惩罚系数，效果会更加可靠。数学上这属于多目标权重问题，但工程上不必搞得很复杂，一个简单的参数扫描就能找到合理范围。
5. 调试过程中踩过的坑和排查心得
5.1 模型报“不可行”是最常见的问题
我最初运行时第一个版本直接报了 infeasible infeasible，YALMIP的debug模式给我列出了很多冗余约束，但真正的问题出在电池日终能量约束上：我要求 E_bat(25) == E_bat(1)，但是又把电池结束能量设成0.2*E_bat_max，然后白天又要求大力放电削峰，几个约束凑在一起根本没有可行解。后来我把日终约束从“等于初始”改成“大于等于初始”，并加入一个足够大的松弛变量来反映“不强制回到原点”，模型就稳定了。
虚拟储能那边的循环约束同样如此。SOE_vs 日终回到初始值，这在夏季冷量充足时不是问题，但遇到连续阴天，楼宇从早上开始就没法蓄冷，晚上自然也没冷量可放，这时候强约束一定会无解。我的做法是把这个循环约束放到目标函数里，用罚函数方式维持“尽量回到初始值”，这样既保证模型可解，也保留了虚拟储能的连续性。
5.2 索引错位带来的魔鬼问题
YALMIP虽然比矩阵建模直观，但时域递推约束里仍然会出现索引错位。E_bat 定义成 T+1 维，初始值占一格，功率变量 P_ch 是 T 维，如果你写 E_bat(1:T) 和 E_bat(2:T+1)，对应关系必须看清楚。我犯过一次把 E_bat(1:T) 写在左边，P_ch 也写在左边，结果导致等式变成“上一时刻能量加上充电功率等于下一时刻能量”的反向递推，整个曲线看起来就像时间倒流一样。排查方法也很土：把约束逐个打印出来，挑三个时刻代入数值手算一遍。
虚拟储能的状态方程里还有一个自损耗项 loss_vs * SOE_vs(1:T)，这个是能量乘以损耗系数，单位是能量。一个很容易犯的错误是写成 loss_vs * SOE_vs(2:T+1)，这样当前时刻的损耗反而受未来状态影响，物理意义就反了。类似问题常出现在论文公式复现上，我在调试时发现曲线异常平滑或异常震荡，第一反应就会去查这类递推方向。
5.3 数值尺度和计算速度的优化
微网模型里同时存在 kW 级的负荷和 MWh 级的电池容量，如果不统一单位，目标函数和约束的数值会差几个数量级，Gurobi和Cplex的预处理会花大量时间，甚至因为容差设置不当导致整数问题误判。我推荐的通用做法是全部使用 MW 和 MWh，功率、能量、价格同时换算。这样虚拟储能 C_vs 如果是2MWh，对应数值就是2，不会出现一个变量量级是 10 的 5 次方、另一个变量量级是 0.01 的情况。
如果你不想引入YALMIP依赖，也可以用linprog或intlinprog直接写。YALMIP在模型较大时会有建模开销，但楼宇微网24小时模型规模根本不用担心；一旦你要做96点日内滚动优化或者考虑多楼宇集群，YALMIP的符号变量可能带来明显内存开销，那时候手动构造矩阵反而更可控。我目前的代码只做到了24点日前调度，96点版本还在测试中，碰到的主要瓶颈是虚拟储能的时间常数和实际楼宇热时间常数不匹配，目前在考虑多时间尺度模型。
5.4 一个提高结果工程可信度的建议
仿真结果再好看，也需要一个合法性检查。我习惯在优化后加一段“体检程序”：把得到的 P_vs 反代回 ETP 温度模型，看室内温度是否落在舒适范围内；把 E_bat 曲线反代进电池功率边界，看有没有超出充放电速率限制。这一步看似多余，却救过我两次，一次是 C_comf 惩罚设太小导致温度偏移超限，另一次是 P_vs 上下限出现单侧失灵，温度曲线在夜间直冲到18度。
6. 最后想说的几点体会
这套代码跑通之后，我最深的体会是：需求侧虚拟储能的价值不在“省了多少电”，而在“转移了多少功率时间”。楼宇微网调度本质上是把能量在时间维度上搬来搬去，物理电池能搬的部分其实很有限，而楼宇本身的热容、充电桩的等待时间、蓄热水箱的存量，都是几乎零成本的时间搬移工具。把它们的“可调度能力”显式建模成虚拟储能系统，Matlab代码的优化空间和可解释性都会被放大很多。
如果你也想复现这个项目，我建议先不要追求一步到位，先跑通“电池+固定负荷”的版本，再往里面加虚拟储能状态方程。每加一个模块，就单独测试一个模块，尤其是把每个变量的上下界打印出来看看是否落在物理可接受范围内。这套代码后续能扩展的方向也很多，比如多楼宇间的共享储能、电动汽车有序充电与楼宇微网的联合优化、或者引入模型预测控制做日内滚动修正。我目前最想做的是把虚拟储能和空调精确温度模型做结合，这样虚拟储能功率就不再是“相对基线的偏移量”，而是直接由温度设定点生成，用户的舒适度表达会更真实。如果你已经在做类似的方向，欢迎一起交流踩坑记录。
