
简介面向能源工程领域研究者与工程师提供了一份论文复现资源聚焦需求侧资源电动汽车EV与温控负载HVAC的灵活性刻画及日前优化调度问题适合智能电网、灵活资源配置方向的学者和中级以上编程能力的技术开发者。资源以docx文档形式呈现共1个文件约24KB包含参数设置、单辆EV的虚拟电池模型、集群EV聚合、基于pulp的日前优化调度建模以及结果可视化等完整Python代码并附有详细注释与公式说明参数设置涵盖一天96个时段、15分钟步长以及日前电价、备用容量价格、EV与HVAC数量、电池容量、充放电功率限制、可调度时段和舒适温度范围等。已有72人浏览学习。通过具体算例可直观看到各资源在不同时段的功率与能量边界理解虚拟电池方法在降低系统运行成本、稳定负荷曲线及提供辅助服务方面的作用代码模块化组织便于在此基础上扩展更大规模场景或迁移到真实数据测试。对于希望掌握数学建模与Python求解能耗优化问题的读者尤为实用。1. 虚拟电池模型需求侧灵活性不是玄学是一组能进优化器的边界约束做电力系统日前优化调度的人最头疼的不是常规机组的出力分配而是需求侧资源怎么进模型。某园区上百台空调、几十台热泵和一排充电桩逐台建模会让混合整数规划规模爆炸直接忽略又会让削峰时段的响应量严重失真。虚拟电池模型Virtual Battery Model解决的就是这个问题把海量需求侧资源聚合成一台虚拟电池用一组线性不等式刻画它的功率上限、能量边界和自放电率嵌入日前优化调度只需几十行约束却保留了需求侧资源灵活性的主要特征。这篇博文我会从数学模型讲到 Python 实现覆盖参数辨识、调度嵌入、结果校验和典型踩坑适合正在复现相关论文或准备把需求响应策略落到优化调度代码里的开发者和研究者。2. 聚合成一台电池虚拟电池模型的数学刻画与参数体系2.1 功率边界、能量边界和自放电率一台虚拟电池只有五六个参数虚拟电池模型的思路一句话可以概括把需求侧资源的灵活性等效成一个储能系统。它有一个状态变量——等效能量 (E(k))有一组调节输入——功率 (P(k))。在日前调度的时间尺度上调度员不需要知道哪台空调在第几分钟转了只需要知道这个聚合体在每一小时能释放或吸收多少功率、能维持多久、恢复到什么状态。这正好是储能约束能表达的东西。离散时间下的核心演化方程写成[ E(k1) (1-\alpha) E(k) \eta_c P_c(k) \Delta t - \frac{P_d(k)}{\eta_d} \Delta t ]其中 (P_c(k) \ge 0) 是充电功率对应负荷侧吸收功率(P_d(k) \ge 0) 是放电功率对应负荷削减、向电网释放灵活性。(\alpha) 是自放电率温控负荷里它对应房间热量渗漏导致的自然回复(\eta_c)、(\eta_d) 是聚合充放电效率对空调集群来说这两个值通常接近 0.95 到 1.0。实际嵌入优化器时我更习惯把充放电合并成净功率 (P(k))约定为正是放电负荷削减随后模型退化成[ E(k1) (1-\alpha) E(k) - P(k) \Delta t ]配合的约束只有三组。功率边界 (P_{min} \le P(k) \le P_{max})能量边界 (E_{min} \le E(k) \le E_{max})以及末态约束 (E(K) \ge E_{target})。下面这张表把虚拟电池的全部参数列全参数符号典型取值范围物理含义功率下限(P_{min})负值MW聚合体允许吸收的最大功率功率上限(P_{max})正值MW聚合体允许释放的最大功率能量下限(E_{min})0 附近等效能量最低值能量上限(E_{max})数十 MWh等效能量最高值初始能量(E_0)(E_{min})~(E_{max})调度开始时等效 SOC自放电率(\alpha)0.01~0.05 / h热量损失和自然回复速率为什么说「盒式约束」够用因为温控负荷的开关周期在分钟级而日前调度步长是 1 小时在一个小时内聚合体的功率调节能力已经被热惯量平滑功率和能量之间的耦合可以用线性演化近似。这个近似不是严谨的充放电可行域但足够给出调度层面的安全边界。2.2 从单个温控负荷的 ETP 模型到聚合体虚拟电池参数从哪里来虚拟电池参数不是拍脑袋定的它来自对设备个体的物理建模。以空调和热泵这类温控负荷为例个体常用一阶等效热参数模型ETP描述[ C \frac{dT}{dt} \frac{T_{out} - T}{R} - s \cdot P_{nom} ](R) 是等效热阻单位 K/W(C) 是热容单位 J/K(T_{out}) 是室外温度(P_{nom}) 是设备额定电功率(s \in {0,1}) 是开关状态制冷模式下运行时取 1。温控逻辑是温度进入设定值附近的死区 ([T_{set}-\delta/2, T_{set}\delta/2]) 内就维持当前开关状态越过边界就翻转。离散化之后每个设备的状态更新写成[ T(k1) T(k) \frac{\Delta t}{R C}\left(T_{out}(k) - T(k)\right) - \frac{\Delta t}{C} s(k) P_{nom} ]这套模型的复杂之处在于 (s(k)) 是整数变量直接把几千台设备的 (s) 一起丢进日前调度是混合整数规划求解器分分钟翻车。常见做法是先用蒙特卡洛模拟一批设备参数跑出大量开关轨迹再做聚合。聚合的产物就是虚拟电池的功率边界和能量边界。模拟时要注意单位统一热阻热容用 SI 单位功率用 kW时间用小时或秒最后折算成 kWh。2.3 采样包络与凸包松弛为什么这套近似在工程上能用论文里标准的做法不是对每台设备建模而是对聚合可行域做凸包松弛。原论文的思路是给定一批随机生成的设备参数每个设备的状态被松弛成连续变量 ([0,1])再用线性约束逼近死区控制逻辑最后求解一系列线性规划把聚合功率投影到每个时段上得到功率包络。这个包络就是虚拟电池参数的来源。我在实际复现时发现完全照搬 LP 端点法对初学者不友好工程上可以先做采样包络蒙特卡洛生成几百组设备参数每组模拟一整天的开关轨迹同一时段聚合功率的 5% 分位数和 95% 分位数作为 (P_{min}(k)) 和 (P_{max}(k))累计能量的包络作为 (E_{min}) 和 (E_{max})。这个方法会略微低估可行域但胜在直观、代码量小且调度结果更安全。想要更精确再退回到论文的顶点枚举 LP我在下一章把两条路线都给出代码。3. Python 实现需求侧资源灵活性聚合数据生成与虚拟电池参数辨识3.1 生成需求侧资源簇设备参数按典型分布随机采样复现的第一步是先造数据。假设我们的对象是某园区内的空调集群500 台设备夏季制冷场景。设备参数分布参考常见文献取值等效热阻 (R) 在 0.5 到 1.5 K/W 之间均匀分布等效热容 (C) 在 2000 到 6000 J/K 之间额定功率 1.0 到 3.5 kW温度设定值围绕 24°C 正态分布死区带宽 0.5 到 1.5°C。代码如下import numpy as np import pandas as pd rng np.random.default_rng(42) n_dev 500 devices pd.DataFrame({ R: rng.uniform(0.5, 1.5, n_dev), # 热阻 K/W C: rng.uniform(2000, 6000, n_dev), # 热容 J/K P_nom: rng.uniform(1.0, 3.5, n_dev), # 额定功率 kW T_set: rng.normal(24.0, 1.0, n_dev), # 设定温度 degC deadband: rng.uniform(0.5, 1.5, n_dev), # 死区带宽 degC T0: rng.normal(24.0, 2.0, n_dev), # 初始温度 }) # 初始开关状态温度高于上限则开机制冷 devices[s] (devices[T0] devices[T_set]).astype(int) print(devices.head())这段代码里R和C决定设备热惯性P_nom决定单体功率上限T_set和deadband决定温控动作点。T0是初始温度如果初始温度高于设定值设备默认开机。要注意的是C单位是 J/K后面模拟时如果时间步长用秒、功率用 kW能量单位需要折算避免出现量级错误。seed42固定随机种子保证复现结果一致。3.2 蒙特卡洛模拟开关轨迹从个体启停到聚合功率曲线有了设备参数接下来做一整天的占空比模拟。步长取 15 分钟模拟精度偏低但对日前调度够用想要更精确可以缩到 1 分钟再重采样。模拟逻辑就是用上一章的 ETP 离散式逐设备推进温度判断是否越界翻转状态同时统计每个时刻的总功率K 24 * 4 # 15分钟步长一天96个点 dt_min 15.0 # 步长分钟 dt_sec dt_min * 60 # 步长秒 T_out np.full(K, 32.0) # 简化场景全天室外温度恒定32度 P_agg np.zeros(K) # 聚合功率 s_mat np.zeros((n_dev, K)) # 开关状态矩阵 for k in range(K): T devices[T0].values T_set devices[T_set].values db devices[deadband].values / 2.0 P_nom devices[P_nom].values R devices[R].values C devices[C].values s devices[s].values # ETP 离散更新 T_new T dt_sec / (R * C) * (T_out[k] - T) \ - (dt_sec / C) * s * P_nom * 1000.0 # W - kW devices[T0] T_new # 死区控制器翻转 upper T_set db lower T_set - db s_new s.copy() s_new[T_new upper] 1.0 s_new[T_new lower] 0.0 devices[s] s_new P_agg[k] np.sum(s * P_nom) # 开机设备的功率求和 s_mat[:, k] s这段模拟里温度更新公式中P_nom * 1000.0是把 kW 转成 W保证与热容 C 的单位一致。死区控制器是纯逻辑判断温度越过上限必须开机低于下限必须关机在区间内维持上一个状态。这就是占空比控制的本质。P_agg[k]是当前时刻所有制冷设备的聚合功率。把这个数组画出来能看到白天负荷高、设备同时开机率高的典型特征。真实场景里T_out应该用当地夏季典型日曲线替换恒定 32°C 只会得到偏保守的边界。3.3 采样包络法辨识虚拟电池参数功率边界与能量边界的提取一条P_agg轨迹只是可行域里的一个样本。要得到虚拟电池参数需要重复上面的模拟生成多组设备参数下的多条功率轨迹然后按分位数取包络。每组参数对应一个蒙特卡洛样本相当于在设备参数空间里采点n_scenario 200 P_samples np.zeros((n_scenario, K)) for i in range(n_scenario): # 重新生成一批设备参数并模拟P_agg存入P_samples[i] # 代码结构与上节相同这里省略重复的模拟逻辑 P_samples[i] simulate_aggregate_power( n_dev500, T_out_profileT_out, rng_seed1000 i ) # 功率边界取每个时段聚合功率的5%和95%分位数 P_min_k np.percentile(P_samples, 5, axis0) P_max_k np.percentile(P_samples, 95, axis0) # 能量边界对净功率做累计再取包络 E_samples np.cumsum(P_samples * dt_min / 60.0, axis1) # 累计 kWh E_min float(np.min(E_samples)) E_max float(np.max(E_samples)) # 虚拟电池参数汇总 vb_params { P_min: float(np.min(P_min_k)), P_max: float(np.max(P_max_k)), E_min: E_min, E_max: E_max, alpha: 0.02, # 待标定见避坑章节 E0: float(np.mean(E_samples[:, 0])), } print(vb_params)如果只看 5% 和 95% 分位数得到的是统计意义上的包络而非最坏情况边界。保守工程做法是把 5% 换成 0% 或 1%代价是可行域变小、调度成本上升。我一般会同时算两组参数一组用于日前优化另一组用于实时校验。E_samples的计算把功率对时间做累计得到等效能量轨迹它的全局最小值和最大值就是 (E_{min}) 和 (E_{max})。注意这个包络法忽略自放电率所以要先假定一个初值后续再用真实设备模拟去标定。3.4 用线性规划提取矩形可行域向论文方法再靠近一步采样包络法的缺点是不够严谨。论文里更规范的做法是把每个设备的开关状态松弛到 ([0,1])把死区控制近似成盒式温度约束然后解一个线性规划求每个时段的聚合功率极值。下面给一个简化版实现目标是最小化/最大化聚合功率约束是松弛后的设备动态from scipy.optimize import linprog K 24 n n_dev * (K 1) # 变量每台设备每个时段的状态 s_i(k) # 但全规模变量太大这里用单设备松弛近似先聚合再尺度化 # 简化策略把设备簇的功率占比视为常数只求总功率包络 c np.zeros(K) # 目标最大化总放电功率 c[:] -1.0 # linprog 默认最小化取负号求最大 A_ub [] # 松弛温度约束的线性表达 b_ub [] bounds [(0, 1)] * K # 开关状态松弛到[0,1] res_min linprog(np.ones(K), A_ubA_ub, b_ubb_ub, boundsbounds, methodhighs) res_max linprog(-np.ones(K), A_ubA_ub, b_ubb_ub, boundsbounds, methodhighs) print(最小聚合功率系数:, res_min.fun) print(最大聚合功率系数:, -res_max.fun)这段代码把完整 LP 简化成了只对总功率作尺度化估计真正的论文级实现需要把每台设备的 (T(k)) 和 (s(k)) 全部纳入变量并构建温度上下界的线性约束。工程上不必死磕这一步采样包络法足够支撑日前调度的决策LP 端点法更适合做离线校验和边界分析。如果你追求复现精度可以从scipy.optimize换到gurobipy或pulp把每台设备的状态矩阵完整建出来。4. 把虚拟电池嵌入日前优化调度目标函数、约束与求解实现4.1 日前调度建模常规机组与虚拟电池的联合优化调度模型的决策变量分为两块常规机组的出力 (P_g(k)) 和虚拟电池的净功率 (P(k))其中 (P(k)) 为正代表负荷削减、释放灵活性。目标函数是全天购电成本加调节成本[ \min \sum_{k0}^{K-1} \left( a_g P_g(k)^2 b_g P_g(k) c_{flex} \max(P(k), 0) \right) \Delta t ]实际代码里二次成本线性化成分段线性或者直接用linprog只保留一次项。约束包括功率平衡 (P_g(k) P(k) D(k))机组上下限、爬坡率以及虚拟电池自身的能量演化[ E(k1) (1-\alpha) E(k) - P(k) \Delta t ][ E_{min} \le E(k) \le E_{max}, \quad P_{min} \le P(k) \le P_{max} ]最后加终端约束 (E(K) \ge E_0)防止优化器把虚拟电池能量在最后时段榨干。这个终端约束是工程里的关键少了它调度结果会表现为最后几个小时灵活性全部耗尽真实设备根本没能力执行。4.2 用 Python 求解 24 小时调度一份可以直接跑的最小实现下面这份代码用scipy.optimize.linprog求解变量顺序是机组出力24 个 虚拟电池功率24 个 虚拟电池能量25 个。求解器选 HiGHS对中等规模 LP 足够快也避免依赖商业求解器import numpy as np from scipy.optimize import linprog K 24 dt 1.0 # 小时 # 虚拟电池参数来自第3章辨识结果 P_min, P_max -50.0, 40.0 # kW负为吸收功率正为释放功率 E_min, E_max 0.0, 120.0 # kWh alpha 0.02 E0 50.0 # 负荷曲线某园区夏季典型日虚构数据 D 300 80 * np.sin(np.linspace(0, 2 * np.pi, K)) 40 * np.cos(np.linspace(0, 4 * np.pi, K)) # 机组参数 Pg_min, Pg_max 100.0, 400.0 # kW b_g, c_flex 0.15, 0.08 # 成本系数 元/kWh # 变量: [Pg(0..23), P(0..23), E(0..24)] n_vars K K (K 1) # 目标机组购电成本 柔性调节成本P0 表示削减负荷需要补偿 c np.zeros(n_vars) c[0:K] b_g # 机组出力成本 c[K:2*K] c_flex * (P_max 0) # 简化放电全部给补偿更精细可用分段 c[2*K:] 0.0 # E 不计成本 # 功率平衡约束Pg(k) P(k) D(k) A_eq [] b_eq [] for k in range(K): row np.zeros(n_vars) row[k] 1.0 # Pg(k) row[K k] 1.0 # P(k) A_eq.append(row) b_eq.append(D[k]) # 虚拟电池能量演化E(k1) (1-alpha) E(k) - P(k) * dt for k in range(K): row np.zeros(n_vars) row[2 * K k 1] 1.0 row[2 * K k] -(1 - alpha) row[K k] dt A_eq.append(row) b_eq.append(0.0) # 不等式约束 A_ub [] b_ub [] # 机组上下限 for k in range(K): row np.zeros(n_vars); row[k] 1.0 A_ub.append(row); b_ub.append(Pg_max) row np.zeros(n_vars); row[k] -1.0 A_ub.append(row); b_ub.append(-Pg_min) # 虚拟电池功率上下限 for k in range(K): row np.zeros(n_vars); row[K k] 1.0 A_ub.append(row); b_ub.append(P_max) row np.zeros(n_vars); row[K k] -1.0 A_ub.append(row); b_ub.append(-P_min) # 能量上下限 for k in range(K 1): row np.zeros(n_vars); row[2 * K k] 1.0 A_ub.append(row); b_ub.append(E_max) row np.zeros(n_vars); row[2 * K k] -1.0 A_ub.append(row); b_ub.append(-E_min) # 终端能量约束E(K) E0 row np.zeros(n_vars); row[2 * K K] -1.0 A_ub.append(row); b_ub.append(-E0) res linprog(c, A_ubA_ub, b_ubb_ub, A_eqA_eq, b_eqb_eq, bounds[(None, None)] * n_vars, methodhighs) Pg res.x[0:K] Pvb res.x[K:2*K] E res.x[2*K:] print(求解成功:, res.success, 总成本:, res.fun) print(虚拟电池首时段调度功率(kW):, Pvb[0])这份代码的核心约束就三个功率平衡等式、能量演化等式、上下限不等式。注意目标函数里P0的补偿成本用了简化写法实际工程中应该引入两个非负变量拆分充放电否则linprog会对负的 (P) 也计成本导致优化结果偏向多吸收功率。拆分方法是P Pd - Pc其中 (Pd, Pc \ge 0)目标里只对 (Pd) 计成本约束里把P替换成Pd - Pc。4.3 从虚拟电池调度值到真实设备指令功率分配器与执行校验调度层给出的是聚合功率 (P(k))这个值最终要变成每台空调的开关指令。分配策略我常用两种按剩余调节容量比例分配或者按优先级顺序分配。比例分配实现最简单def dispatch_to_devices(P_target, s_current, devices): 把虚拟电池的调度功率目标分配给真实设备。 思路按每台设备可调节容量加权分配。 P_nom devices[P_nom].values s s_current.values # 可减载容量正在运行的设备可以关 can_reduce s * P_nom # 可增载容量未运行的设备可以开 can_increase (1 - s) * P_nom if P_target 0: # 需要减载/放电 weights can_reduce / np.sum(can_reduce 1e-6) delta np.minimum(weights * P_target, can_reduce) s_new s - (delta 0).astype(int) else: # 需要增载/充电 weights can_increase / np.sum(can_increase 1e-6) delta np.minimum(weights * (-P_target), can_increase) s_new s (delta 0).astype(int) return s_new分配到设备之后必须做执行校验把s_new重新输入 ETP 模拟器跑 15 分钟或 1 小时看真实聚合功率与调度目标的偏差。如果偏差超过 5%说明虚拟电池参数偏乐观需要回到第 3 章重新辨识。这一步是论文复现和工程落地的分水岭——只跑调度不解设备永远发现不了聚合松弛带来的误差。5. 论文复现避坑指南从模型退化为黑匣子的 5 个典型问题5.1 自放电率被忽略能量轨迹全天漂移现象虚拟电池的 SOC 曲线随时间单调下降明明调度设置了终端约束真实设备模拟出来却提前两个小时「没电了」。原因我把 (\alpha) 设成 0 跑通了第一版调度模型里能量只减不增。但真实温控负荷本身有热量渗漏房间在无人控制时温度会自然飘向室外温度这部分等效能量衰减就是自放电。忽略它模型高估了可用灵活性。解决用真实设备模拟器做开环实验——不给任何调度指令只统计聚合功率自然变化折算成等效能量衰减率标定 (\alpha)。对制冷场景(\alpha) 一般在 0.01 到 0.05 之间夏季室外温度越高数值越大。每次换场景都要重新标定不能一套参数用到底。5.2 充放电效率不对称导致虚拟电池凭空造出能量现象调度结果中虚拟电池的放出能量大于充入能量系统总成本为负优化器在「套利」储能效率而不是真正调度资源。原因效率参数放错位置。把演化方程写成 (E(k1)E(k)-\eta P(k)\Delta t) 时如果 (\eta1)放电过程反而增加等效能量模型出现能量永动机。这是复现虚拟电池最容易踩的坑尤其是参考不同论文时符号习惯不一致。解决严格使用充放电分离式充电效率 (\eta_c 1) 表示能量损失放电效率 (\eta_d 1) 同样表示损失。统一写成[ E(k1) (1-\alpha) E(k) \eta_c P_c(k)\Delta t - \frac{P_d(k)\Delta t}{\eta_d} ]并在求解前手动验算给定相同充放电量两轮循环之后 (E) 必须减少否则参数就是反的。5.3 矩形包络无条件外推真实设备根本執行不了现象调度优化给出连续 4 小时放电 35 kW设备级校验时前 2 小时能跟上第 3 小时开始聚合功率明显低于目标偏差超过 20%。原因采样包络法只在生成场景的室外温度、设备参数分布范围内有效。真实运行中如果设备处于不同热状态聚合可行域会收窄矩形包络过于乐观。尤其在傍晚空调自然需求上升可削减容量本身就减少。解决把虚拟电池参数按运行场景分段标定比如早晚峰各一套参数或者在校验环节引入「滚动预检」执行前 1 小时用当前设备状态重新做一次短周期模拟确认目标功率可达再下发指令。倾向保守时把功率边界的百分位数从 95% 下调到 90%损失一点经济性换来执行安全。5.4 调度步长和控制步长不一致SOC 在小时内失真现象日前调度步长 1 小时虚拟电池能量演化线性但真实设备开关周期 5 到 15 分钟聚合功率在一个小时内大幅波动。调度认为能量平滑变化实际却剧烈震荡终端约束形同虚设。原因模型的时间粒度与实际控制粒度不匹配。虚拟电池设计初衷是站在 1 小时粒度上的近似但设备层面的占空比波动无法被线性模型表达。解决调度模型保持 1 小时不变在功率分配器环节切换到 15 分钟步长把 (P(k)) 按小时目标分解成 4 个 15 分钟目标并在虚拟电池约束里把能量边界放宽 5% 作为波动裕度。这个方法保留了日前调度的计算规模又把执行偏差控制在可接受范围。5.5 初始 SOC 标定不对首个调度时段直接失衡现象调度结果给的第一时段灵活性很大但设备集群没有对应能力开局就偏差 30%。原因(E_0) 设置不反映真实设备状态。初始设备开关状态和温度分布决定了等效初始能量。用固定的 (E_00.5(E_{min}E_{max})) 是懒人做法真实场景里早晨的空调集群大部分处于停机状态等效能量偏上限。解决调度前做 30 分钟预模拟用当前设备温度和开关状态跑一次开环模拟统计等效能量作为 (E_0)。调度模型再把 (E_0) 所在区间收窄比如要求初始能量落在 ([E_0 - 5, E_0 5]) kWh既给了优化器自由度又锁住了真实性。6. 虚拟电池模型的验证方法与进阶方向验证一套虚拟电池参数是否合格最直接的办法是闭环回灌把日前调度的 (P^*(k)) 作为目标输入真实设备逐台模拟器对比每个时段的实际聚合功率与调度目标的偏差。偏差绝对值超过 5% 的时段占比超过 10%参数就要重标。更严格的做法是同时记录模拟中的等效能量轨迹检查它是否落在 ([E_{min}, E_{max}]) 区间内连续越界两台时就说明边界辨识失败。进阶方向上值得做的是三件事。第一把效率从常数改成线性函数对应不同功率水平下的效率差异第二加入 SOC 恢复区间约束做到「调度结束后能量回到指定带」避免资源被过度消耗第三多类负荷并联时把空调、热水器和储能各自聚合成一个虚拟电池再在优化器里并联而不是一锅烩成单台。三种方向都会增加约束数量但求解规模仍在可控范围。我自己的习惯是最后做一张图把虚拟电池调度曲线和逐设备模拟曲线叠在一起偏差大的时段高亮标红然后回去查那个时段是功率边界问题还是能量边界问题。这个方法帮我解决了至少三个看似是参数问题、实则是约束遗漏的bug。虚拟电池模型的优雅之处在于它把需求侧资源变成了一组可以放进任何优化器的线性边界但边界背后的物理意义永远要靠真实设备模拟来背书。希望这篇笔记对你复现和落地这个方向有帮助。本文还有配套的精品资源点击获取