尧图网页设计 
简介面向能源工程研究者与智能电网开发者的技术资料围绕电动汽车与温控负载HVAC的灵活性刻画给出基于虚拟电池模型的日前优化调度Python实现涵盖参数设置、单台设备到集群的VB模型推导以及借助pulp库构建的线性规划提前计划策略。全包仅1个docx文件约24KB以论文复现笔记形式串联数学模型、代码实现与逐段解释配合96时段、15分钟步长的算例把EV充放电边界与HVAC舒适度约束转写为可求解的优化形式。已有72人学习适合具备中级以上编程能力、希望掌握需求响应资源建模与成本优化求解思路的读者通过图表化的逐时段结果呈现可直观理解负荷曲线稳定与备用服务保障的验证过程并为自行扩展规模、打磨模型细节提供可复用的脚本骨架与排错参照。1. 需求侧资源灵活性刻画日前调度里为什么先要把它压成一块虚拟电池日前优化调度里最先卡住的往往不是机组组合而是需求侧那一堆空调、热水器、储能、可平移负荷。它们各自有温度死区、启停次数、用户舒适度要求直接建模会变成上万变量加大量 0-1 变量求解器要么跑不动要么跑出来的解没人敢执行。虚拟电池模型Virtual Battery Model的价值就在这里把一批异质柔性负荷的可行域在功率-能量平面上外近似成一块等效电池——有充放电功率上下界、能量上下界、自放电系数和充放电效率。聚合商向上报的不再是一条曲线而是几个参数加一条能量轨迹调度侧用线性约束就能接纳。对做电力系统优化的 Python 工程师来说这套东西的收益很具体建模成本从 MILP 降到 LP求解时间从分钟级降到秒级而且参数还能直接对接日前市场的报价逻辑。适合两类人一类是要复现聚合灵活性论文里的模型一类是要在真实项目里把需求侧资源塞进日前出清。下面按模型推导、参数辨识、Python 落地、接入调度、结果校验的顺序走一遍。2. 虚拟电池模型的功率-能量边界与状态方程推导2.1 从温控负荷的一阶等效热参数模型起步空调、热泵这类温控负荷最常用的降阶模型是一阶等效热参数模型ETP离散形式写成T_in(k1) T_in(k) Δt/(R·C)·(T_out(k) − T_in(k)) − (Δt·η/C)·u(k)·P_rated其中 u(k) 取 0 或 1表示定频压缩机停或开。这套模型的参数在工程上是可以标定的常见做法是用历史温度曲线做最小二乘拟合而不是直接查铭牌。符号含义单位典型取值区间C等效热容房间家具J/K3e6 ~ 1e7R等效热阻的倒数关系此处为热阻K/W3e-3 ~ 1e-2η制冷能效比 COP—2.5 ~ 3.5P_rated额定电功率W800 ~ 2500Δt控制步长s90015 minδ温度死区宽度K0.5 ~ 2.0注意Δt 直接决定后面自放电系数 a 的大小。Δt 从 15 分钟改成 1 小时a 会从 0.97 掉到 0.88 附近能量约束的时段耦合强度差一个量级同一组虚拟电池参数不能跨步长复用。2.2 把温度状态平移成虚拟能量定义虚拟能量 E(k) C·(T_max − T_in(k)) / 3600单位 Wh其中 T_max T_set δ/2。这个量在物理上就是「当前室温距离允许温度上界还差多少热容量的能量」。把它代进 ETP 方程室温项会被消掉E(k1) E(k) − C/3600·[Δt/(R·C)·(T_out − T_in(k)) − Δt·η·u(k)·P_rated/C]再利用 T_out − T_in(k) (T_out − T_max) 3600·E(k)/C整理后得到E(k1) a·E(k) − gain b·u(k)其中三个系数分别定义为a 1 − Δt/(R·C)保持系数衡量一个步长内热量自然流失的比例b Δt·η·P_rated / 3600满功率运行一个步长注入的等效能量Whgain Δt·(T_out − T_max) / (3600·R)环境热增益折算成能量Wh对照通用虚拟电池形式 E(k1) a·E(k) η_ch·p_ch(k)·Δt − p_dis(k)·Δt/η_dis可以看出单台制冷设备的虚拟电池是「只能充电」的p_dis 恒为 0充电效率就是 COP。单机约束只有两组功率约束0 ≤ p(k) ≤ P_rated能量约束0 ≤ E(k) ≤ C·δ/3600这两个约束看起来简单但聚合之后可行域的形状完全不是矩形。2.3 聚合为什么不能只把上下界相加N 台设备的聚合可行域是各单机可行域的闵可夫斯基和。只把 P_rated 逐个相加得到的是保守内近似一定可行但会低估灵活性总量反过来把能量上下界逐个相加得到的区间里会混进物理不可达的状态——典型情形是每台设备都取能量上界但那个状态对应的功率方向不可能同时满足。所以标准做法分两步先把闵可夫斯基和的凸包边界算出来再用一个虚拟电池多面体把它包住。外近似的目标函数可以写成最小化多面体体积min (P_max − P_min) λ·(E_max − E_min)约束是让虚拟电池的可行域包含所有采样到的 (p, E) 可行点。这一步是纯 LP用 scipy 或 PuLP 都能解。λ 用来调节功率方向和能量方向的权重λ 取大意味着更看重压缩能量窗口通常在舒适度敏感的场景里这么设。3. Python 复现从单机参数到聚合虚拟电池参数3.1 环境准备与依赖清单Python 环境用 venv 就够不需要 conda。用 vscode 配置 python 开发环境时把解释器指到 .venv/bin/pythonWindows 是 .venv\Scripts\python.exe四件套依赖装齐python -m venv .venv source .venv/bin/activate # Windows: .venv\Scripts\activate pip install numpy pandas scipy matplotlib pulpnumpy 负责状态枚举和向量化pandas 管 24 时段的时间序列scipy 做拟合和插值pulp 写调度模型并调用 CBC 求解器matplotlib 只是最后画对比图用。3.2 单机虚拟电池边界的可达集枚举单台定频设备只有开停两档直接把可达能量格点做动态规划枚举比写解析边界更不容易出错import numpy as np def single_vb_params(C, R, eta, P_rated, T_out, T_set, delta, dt): 把 ETP 参数折算成单机虚拟电池的三个系数和容量上限。 T_max T_set delta / 2.0 T_min T_set - delta / 2.0 a 1.0 - dt / (R * C) b dt * eta * P_rated / 3600.0 # Wh满功率一步注入 gain dt * (T_out - T_max) / (3600.0 * R) # Wh环境热增益 E_cap C * (T_max - T_min) / 3600.0 # Wh能量窗口宽度 return dict(aa, bb, gaingain, E_capE_cap, P_ratedP_rated) def reachable_energy(params, horizon, levels80): 给定初始能量枚举每个时段可达的能量格点集合。 a, b, gain, E_cap params[a], params[b], params[gain], params[E_cap] grid np.linspace(0.0, E_cap, levels 1) tol E_cap / levels * 0.5 cur {0.0} # 初始能量由初始室温决定 history [sorted(cur)] for _ in range(horizon): nxt set() for e in cur: for u in (0.0, 1.0): # 定频机只有开 / 停 e_next a * e - gain b * u if -tol e_next E_cap tol: e_next min(max(e_next, 0.0), E_cap) nxt.add(round(e_next, 6)) cur nxt history.append(sorted(cur)) return grid, history逻辑上说明三点。a 越大初始能量对后续时段的影响衰减越慢24 步之后残余 a^24a0.97 时还有约 0.48所以日循环约束不是摆设。b 和 E_cap 的比值决定了单机的最小启停占空比示例参数下 b≈1.05 kWh、E_cap≈1.67 kWh占空比下限约 0.45低于这个值室温必然越死区。gain 是唯一和室外温度预测挂钩的项日前调度里它是不确定性的主要入口建议在第 5 章的校验里单独扫一遍。3.3 聚合采样与外近似拟合聚合层面用蒙特卡洛采样近似闵可夫斯基和再对每个时段取分位数作为边界。这样写代码短而且能直接处理异质设备def sample_aggregate(params, n_units, horizon, n_samples3000, seed0): 抽取聚合 (p, E) 可行点。p 单位 kWE 单位 kWh。 rng np.random.default_rng(seed) P np.zeros((n_samples, horizon)) E np.zeros((n_samples, horizon)) for s in range(n_samples): e_tot rng.uniform(0.0, params[E_cap] * n_units) for k in range(horizon): p_k 0.0 for _ in range(n_units): p_k rng.random() * params[P_rated] / 1000.0 # 连续松弛 e_tot params[a] * e_tot - params[gain] * n_units params[b] * p_k * 1000.0 e_tot np.clip(e_tot, 0.0, params[E_cap] * n_units) P[s, k], E[s, k] p_k, e_tot return P, E P, E sample_aggregate(single_vb_params(6e6, 5e-3, 2.8, 1500.0, 35.0, 25.0, 1.0, 900.0), n_units200, horizon24) P_min, P_max P.min(axis0), P.max(axis0) E_min, E_max E.min(axis0), E.max(axis0)外近似必须用真正的 min/max 而不是分位数否则边界会切掉真实可行点回代时必然越界。想收紧边界的话改用分位数之后再跑一次可行性回归把采样点逐个代入边界统计被切掉的比例控制在 1% 以内再接受。输出量含义单位示例量级200 台P_min(t)聚合最小充电功率kW0P_max(t)聚合最大充电功率kW300E_min(t)聚合能量下界kWh40 ~ 60E_max(t)聚合能量上界kWh300 ~ 3334. 虚拟电池接入日前优化调度目标函数与约束的线性化落地4.1 目标函数怎么设日前调度的目标是最小化购电成本加虚拟电池的等效退化惩罚min Σ_t [ c_DA(t)·P_grid(t)·Δt c_deg·(P_ch(t) P_dis(t))·Δt ]c_DA(t) 是 24 点分时电价c_deg 是舒适度或设备磨损的折算系数。c_deg 取 0 会让模型把充放电动作推满实际执行时室温抖得厉害一般取 0.01 到 0.05 元/kWh 量级就能压住频繁动作。这个写法是全线性的不需要任何二进制变量。4.2 约束集合与物理含义功率平衡P_grid(t) P_dis(t) L_base(t) P_ch(t)基线负荷由电网和放电共同供给能量动态E(t1) a·E(t) η_ch·P_ch(t)·Δt − P_dis(t)·Δt/η_dis能量边界E_min(t) ≤ E(t) ≤ E_max(t)功率边界0 ≤ P_ch(t) ≤ P_max(t)0 ≤ P_dis(t) ≤ P_dis_max(t)日循环E(0) E(T) E_init避免模型提前把能量抽干并网点容量0 ≤ P_grid(t) ≤ P_trafo温控负荷聚合的虚拟电池放电能力通常按可削减量单独算制冷设备本身不放能放电项来自负荷削减在功率平衡里的体现。4.3 PuLP 实现与求解import pulp import numpy as np T, dt 24, 1.0 price np.array([0.35, 0.32, 0.30] [0.55] * 5 [0.85] * 4 [0.70] * 4 [0.90] * 4 [0.60] * 3 [0.40]) # 示例 24 点电价 base np.full(T, 1800.0) # 基线负荷 kW? 见说明 base base / 1000.0 a, eta_ch, eta_dis 0.97, 0.93, 0.93 P_ch_max np.full(T, 300.0) P_dis_max np.full(T, 120.0) E_min np.full(T, 50.0) E_max_t np.linspace(300.0, 333.0, T) E_init, P_trafo, c_deg 180.0, 500.0, 0.02 m pulp.LpProblem(day_ahead_vb, pulp.LpMinimize) P_grid pulp.LpVariable.dicts(P_grid, range(T), lowBound0, upBoundP_trafo) P_ch pulp.LpVariable.dicts(P_ch, range(T), lowBound0) P_dis pulp.LpVariable.dicts(P_dis, range(T), lowBound0) E pulp.LpVariable.dicts(E, range(T 1), lowBound0) m pulp.lpSum(price[t] * P_grid[t] * dt c_deg * (P_ch[t] P_dis[t]) * dt for t in range(T)) for t in range(T): m P_grid[t] P_dis[t] base[t] P_ch[t] # 功率平衡 m P_ch[t] P_ch_max[t] m P_dis[t] P_dis_max[t] m E[t] E_min[t] m E[t] E_max_t[t] m E[t 1] a * E[t] eta_ch * P_ch[t] * dt - P_dis[t] * dt / eta_dis m E[T] E_min[T] m E[T] E_max_t[T] m E[0] E_init m E[T] E_init # 日循环 m.solve(pulp.PULP_CBC_CMD(msgFalse)) grid np.array([P_grid[t].value() for t in range(T)]) soc np.array([E[t].value() for t in range(T 1)])代码里有几处容易写错。第一E 的索引长度是 T1边界约束要覆盖 E[T]否则最后一个时段会无界。第二base 单位要和 P_ch_max 一致示例里统一用 kW。第三eta_ch 和 eta_dis 不能都取 1取 1 会让 LP 认为能量无损回代时 SOC 曲线和实际差一大截。第四日循环约束缺了的话模型会在最后几个时段把虚拟电池抽干来套利前 20 个时段的解都是假的。4.4 参数怎么调a 越接近 1跨时段耦合越强削峰效果越明显但对相邻时段功率变化率更敏感P_ch_max 一般取聚合额定功率的 100%P_dis_max 取 30% 到 40%因为可削减负荷远小于可增加负荷E_min 不能设 0留 15% 到 20% 的底能量否则室温贴着死区上界跑回代容易越界c_deg 从 0.01 起调观察到 P_ch 曲线出现高频抖动就加大5. 复现校验三个把虚拟电池调稳的参数与两个常见误用5.1 用回代仿真验证 LP 解是否真的可行LP 里满足能量边界不代表聚合后每台设备都不越死区因为外近似放大了可行域。最便宜也最有效的验收手段是把最优功率轨迹回代到单机模型def replay_single_machine(vb_traj, params, e_init): 把聚合功率按占比摊到单机检查能量是否越界。 e, violations, worst e_init, 0, 0.0 for p in vb_traj: e params[a] * e - params[gain] params[b] * p low, high 0.0, params[E_cap] if e low - 1e-6 or e high 1e-6: violations 1 worst max(worst, max(low - e, e - high)) e min(max(e, low), high) return violations, worstviolations 超过总步数的 10%就说明外近似太松需要把 E_max 按采样点收紧或者调大 c_deg 抑制激进动作。worst 反映越界幅度的量级控制在小数点后两位以内一般可以接受。5.2 三个最敏感的参数参数敏感方向调参建议a保持系数24 步残余 a^24直接决定耦合强度Δt 定了就由 R·C 标定别硬调η_ch / η_dis效率取 1 会系统性低估能量损耗制冷聚合 η_ch 取 0.85 ~ 0.95E_min / E_max窗口越宽灵活性越大回代越容易越界留 15% 底能量上界按采样 99 分位a 的敏感度可以在定步长前先扫一遍把 R·C 从 1e4 到 1e5 拉一条曲线看 a 从 0.88 到 0.99 变化时同样的削峰目标下 P_dis_max 需求差多少。差得大说明当前步长选得不合适应该把 Δt 改小而不是硬凑参数。5.3 两个常见误用第一个是把功率边界和能量边界独立聚合。正确做法是先算聚合可行域的凸包再用虚拟电池包住它中间那步闵可夫斯基和的采样不能省。第二个是把自放电项当 1 处理。储能的日自放电常被忽略温控负荷不能这么干——温差驱动的热增益每步都在扣能量忽略它会让 LP 高估可用灵活性实盘执行时削峰量直接缩水三成以上。5.4 用峰谷差和购电成本做验收指标最后拿两个数字收口一个对比基准场景虚拟电池功率全置 0一个对比接入场景cost float(np.sum(price * grid)) peak, valley grid.max(), grid.min() print(f日购电成本 {cost:.1f} 元峰谷差 {peak - valley:.1f} kW)我一般会把这个脚本挂在参数扫描外面a、c_deg、E_min 三个变量各取三档做个 27 组的小网格只看峰谷差降幅和 violations 两个指标。降幅不到 5% 的组合直接丢掉剩下的再按成本排序。这套流程跑下来虚拟电池参数就从「论文里抄的」变成「自己数据标定的」接进日前出清才站得住。本文还有配套的精品资源点击获取 需求侧柔性负荷虚拟电池建模与日前调度Python优化实践
简介面向能源工程研究者与智能电网开发者的技术资料围绕电动汽车与温控负载HVAC的灵活性刻画给出基于虚拟电池模型的日前优化调度Python实现涵盖参数设置、单台设备到集群的VB模型推导以及借助pulp库构建的线性规划提前计划策略。全包仅1个docx文件约24KB以论文复现笔记形式串联数学模型、代码实现与逐段解释配合96时段、15分钟步长的算例把EV充放电边界与HVAC舒适度约束转写为可求解的优化形式。已有72人学习适合具备中级以上编程能力、希望掌握需求响应资源建模与成本优化求解思路的读者通过图表化的逐时段结果呈现可直观理解负荷曲线稳定与备用服务保障的验证过程并为自行扩展规模、打磨模型细节提供可复用的脚本骨架与排错参照。1. 需求侧资源灵活性刻画日前调度里为什么先要把它压成一块虚拟电池日前优化调度里最先卡住的往往不是机组组合而是需求侧那一堆空调、热水器、储能、可平移负荷。它们各自有温度死区、启停次数、用户舒适度要求直接建模会变成上万变量加大量 0-1 变量求解器要么跑不动要么跑出来的解没人敢执行。虚拟电池模型Virtual Battery Model的价值就在这里把一批异质柔性负荷的可行域在功率-能量平面上外近似成一块等效电池——有充放电功率上下界、能量上下界、自放电系数和充放电效率。聚合商向上报的不再是一条曲线而是几个参数加一条能量轨迹调度侧用线性约束就能接纳。对做电力系统优化的 Python 工程师来说这套东西的收益很具体建模成本从 MILP 降到 LP求解时间从分钟级降到秒级而且参数还能直接对接日前市场的报价逻辑。适合两类人一类是要复现聚合灵活性论文里的模型一类是要在真实项目里把需求侧资源塞进日前出清。下面按模型推导、参数辨识、Python 落地、接入调度、结果校验的顺序走一遍。2. 虚拟电池模型的功率-能量边界与状态方程推导2.1 从温控负荷的一阶等效热参数模型起步空调、热泵这类温控负荷最常用的降阶模型是一阶等效热参数模型ETP离散形式写成T_in(k1) T_in(k) Δt/(R·C)·(T_out(k) − T_in(k)) − (Δt·η/C)·u(k)·P_rated其中 u(k) 取 0 或 1表示定频压缩机停或开。这套模型的参数在工程上是可以标定的常见做法是用历史温度曲线做最小二乘拟合而不是直接查铭牌。符号含义单位典型取值区间C等效热容房间家具J/K3e6 ~ 1e7R等效热阻的倒数关系此处为热阻K/W3e-3 ~ 1e-2η制冷能效比 COP—2.5 ~ 3.5P_rated额定电功率W800 ~ 2500Δt控制步长s90015 minδ温度死区宽度K0.5 ~ 2.0注意Δt 直接决定后面自放电系数 a 的大小。Δt 从 15 分钟改成 1 小时a 会从 0.97 掉到 0.88 附近能量约束的时段耦合强度差一个量级同一组虚拟电池参数不能跨步长复用。2.2 把温度状态平移成虚拟能量定义虚拟能量 E(k) C·(T_max − T_in(k)) / 3600单位 Wh其中 T_max T_set δ/2。这个量在物理上就是「当前室温距离允许温度上界还差多少热容量的能量」。把它代进 ETP 方程室温项会被消掉E(k1) E(k) − C/3600·[Δt/(R·C)·(T_out − T_in(k)) − Δt·η·u(k)·P_rated/C]再利用 T_out − T_in(k) (T_out − T_max) 3600·E(k)/C整理后得到E(k1) a·E(k) − gain b·u(k)其中三个系数分别定义为a 1 − Δt/(R·C)保持系数衡量一个步长内热量自然流失的比例b Δt·η·P_rated / 3600满功率运行一个步长注入的等效能量Whgain Δt·(T_out − T_max) / (3600·R)环境热增益折算成能量Wh对照通用虚拟电池形式 E(k1) a·E(k) η_ch·p_ch(k)·Δt − p_dis(k)·Δt/η_dis可以看出单台制冷设备的虚拟电池是「只能充电」的p_dis 恒为 0充电效率就是 COP。单机约束只有两组功率约束0 ≤ p(k) ≤ P_rated能量约束0 ≤ E(k) ≤ C·δ/3600这两个约束看起来简单但聚合之后可行域的形状完全不是矩形。2.3 聚合为什么不能只把上下界相加N 台设备的聚合可行域是各单机可行域的闵可夫斯基和。只把 P_rated 逐个相加得到的是保守内近似一定可行但会低估灵活性总量反过来把能量上下界逐个相加得到的区间里会混进物理不可达的状态——典型情形是每台设备都取能量上界但那个状态对应的功率方向不可能同时满足。所以标准做法分两步先把闵可夫斯基和的凸包边界算出来再用一个虚拟电池多面体把它包住。外近似的目标函数可以写成最小化多面体体积min (P_max − P_min) λ·(E_max − E_min)约束是让虚拟电池的可行域包含所有采样到的 (p, E) 可行点。这一步是纯 LP用 scipy 或 PuLP 都能解。λ 用来调节功率方向和能量方向的权重λ 取大意味着更看重压缩能量窗口通常在舒适度敏感的场景里这么设。3. Python 复现从单机参数到聚合虚拟电池参数3.1 环境准备与依赖清单Python 环境用 venv 就够不需要 conda。用 vscode 配置 python 开发环境时把解释器指到 .venv/bin/pythonWindows 是 .venv\Scripts\python.exe四件套依赖装齐python -m venv .venv source .venv/bin/activate # Windows: .venv\Scripts\activate pip install numpy pandas scipy matplotlib pulpnumpy 负责状态枚举和向量化pandas 管 24 时段的时间序列scipy 做拟合和插值pulp 写调度模型并调用 CBC 求解器matplotlib 只是最后画对比图用。3.2 单机虚拟电池边界的可达集枚举单台定频设备只有开停两档直接把可达能量格点做动态规划枚举比写解析边界更不容易出错import numpy as np def single_vb_params(C, R, eta, P_rated, T_out, T_set, delta, dt): 把 ETP 参数折算成单机虚拟电池的三个系数和容量上限。 T_max T_set delta / 2.0 T_min T_set - delta / 2.0 a 1.0 - dt / (R * C) b dt * eta * P_rated / 3600.0 # Wh满功率一步注入 gain dt * (T_out - T_max) / (3600.0 * R) # Wh环境热增益 E_cap C * (T_max - T_min) / 3600.0 # Wh能量窗口宽度 return dict(aa, bb, gaingain, E_capE_cap, P_ratedP_rated) def reachable_energy(params, horizon, levels80): 给定初始能量枚举每个时段可达的能量格点集合。 a, b, gain, E_cap params[a], params[b], params[gain], params[E_cap] grid np.linspace(0.0, E_cap, levels 1) tol E_cap / levels * 0.5 cur {0.0} # 初始能量由初始室温决定 history [sorted(cur)] for _ in range(horizon): nxt set() for e in cur: for u in (0.0, 1.0): # 定频机只有开 / 停 e_next a * e - gain b * u if -tol e_next E_cap tol: e_next min(max(e_next, 0.0), E_cap) nxt.add(round(e_next, 6)) cur nxt history.append(sorted(cur)) return grid, history逻辑上说明三点。a 越大初始能量对后续时段的影响衰减越慢24 步之后残余 a^24a0.97 时还有约 0.48所以日循环约束不是摆设。b 和 E_cap 的比值决定了单机的最小启停占空比示例参数下 b≈1.05 kWh、E_cap≈1.67 kWh占空比下限约 0.45低于这个值室温必然越死区。gain 是唯一和室外温度预测挂钩的项日前调度里它是不确定性的主要入口建议在第 5 章的校验里单独扫一遍。3.3 聚合采样与外近似拟合聚合层面用蒙特卡洛采样近似闵可夫斯基和再对每个时段取分位数作为边界。这样写代码短而且能直接处理异质设备def sample_aggregate(params, n_units, horizon, n_samples3000, seed0): 抽取聚合 (p, E) 可行点。p 单位 kWE 单位 kWh。 rng np.random.default_rng(seed) P np.zeros((n_samples, horizon)) E np.zeros((n_samples, horizon)) for s in range(n_samples): e_tot rng.uniform(0.0, params[E_cap] * n_units) for k in range(horizon): p_k 0.0 for _ in range(n_units): p_k rng.random() * params[P_rated] / 1000.0 # 连续松弛 e_tot params[a] * e_tot - params[gain] * n_units params[b] * p_k * 1000.0 e_tot np.clip(e_tot, 0.0, params[E_cap] * n_units) P[s, k], E[s, k] p_k, e_tot return P, E P, E sample_aggregate(single_vb_params(6e6, 5e-3, 2.8, 1500.0, 35.0, 25.0, 1.0, 900.0), n_units200, horizon24) P_min, P_max P.min(axis0), P.max(axis0) E_min, E_max E.min(axis0), E.max(axis0)外近似必须用真正的 min/max 而不是分位数否则边界会切掉真实可行点回代时必然越界。想收紧边界的话改用分位数之后再跑一次可行性回归把采样点逐个代入边界统计被切掉的比例控制在 1% 以内再接受。输出量含义单位示例量级200 台P_min(t)聚合最小充电功率kW0P_max(t)聚合最大充电功率kW300E_min(t)聚合能量下界kWh40 ~ 60E_max(t)聚合能量上界kWh300 ~ 3334. 虚拟电池接入日前优化调度目标函数与约束的线性化落地4.1 目标函数怎么设日前调度的目标是最小化购电成本加虚拟电池的等效退化惩罚min Σ_t [ c_DA(t)·P_grid(t)·Δt c_deg·(P_ch(t) P_dis(t))·Δt ]c_DA(t) 是 24 点分时电价c_deg 是舒适度或设备磨损的折算系数。c_deg 取 0 会让模型把充放电动作推满实际执行时室温抖得厉害一般取 0.01 到 0.05 元/kWh 量级就能压住频繁动作。这个写法是全线性的不需要任何二进制变量。4.2 约束集合与物理含义功率平衡P_grid(t) P_dis(t) L_base(t) P_ch(t)基线负荷由电网和放电共同供给能量动态E(t1) a·E(t) η_ch·P_ch(t)·Δt − P_dis(t)·Δt/η_dis能量边界E_min(t) ≤ E(t) ≤ E_max(t)功率边界0 ≤ P_ch(t) ≤ P_max(t)0 ≤ P_dis(t) ≤ P_dis_max(t)日循环E(0) E(T) E_init避免模型提前把能量抽干并网点容量0 ≤ P_grid(t) ≤ P_trafo温控负荷聚合的虚拟电池放电能力通常按可削减量单独算制冷设备本身不放能放电项来自负荷削减在功率平衡里的体现。4.3 PuLP 实现与求解import pulp import numpy as np T, dt 24, 1.0 price np.array([0.35, 0.32, 0.30] [0.55] * 5 [0.85] * 4 [0.70] * 4 [0.90] * 4 [0.60] * 3 [0.40]) # 示例 24 点电价 base np.full(T, 1800.0) # 基线负荷 kW? 见说明 base base / 1000.0 a, eta_ch, eta_dis 0.97, 0.93, 0.93 P_ch_max np.full(T, 300.0) P_dis_max np.full(T, 120.0) E_min np.full(T, 50.0) E_max_t np.linspace(300.0, 333.0, T) E_init, P_trafo, c_deg 180.0, 500.0, 0.02 m pulp.LpProblem(day_ahead_vb, pulp.LpMinimize) P_grid pulp.LpVariable.dicts(P_grid, range(T), lowBound0, upBoundP_trafo) P_ch pulp.LpVariable.dicts(P_ch, range(T), lowBound0) P_dis pulp.LpVariable.dicts(P_dis, range(T), lowBound0) E pulp.LpVariable.dicts(E, range(T 1), lowBound0) m pulp.lpSum(price[t] * P_grid[t] * dt c_deg * (P_ch[t] P_dis[t]) * dt for t in range(T)) for t in range(T): m P_grid[t] P_dis[t] base[t] P_ch[t] # 功率平衡 m P_ch[t] P_ch_max[t] m P_dis[t] P_dis_max[t] m E[t] E_min[t] m E[t] E_max_t[t] m E[t 1] a * E[t] eta_ch * P_ch[t] * dt - P_dis[t] * dt / eta_dis m E[T] E_min[T] m E[T] E_max_t[T] m E[0] E_init m E[T] E_init # 日循环 m.solve(pulp.PULP_CBC_CMD(msgFalse)) grid np.array([P_grid[t].value() for t in range(T)]) soc np.array([E[t].value() for t in range(T 1)])代码里有几处容易写错。第一E 的索引长度是 T1边界约束要覆盖 E[T]否则最后一个时段会无界。第二base 单位要和 P_ch_max 一致示例里统一用 kW。第三eta_ch 和 eta_dis 不能都取 1取 1 会让 LP 认为能量无损回代时 SOC 曲线和实际差一大截。第四日循环约束缺了的话模型会在最后几个时段把虚拟电池抽干来套利前 20 个时段的解都是假的。4.4 参数怎么调a 越接近 1跨时段耦合越强削峰效果越明显但对相邻时段功率变化率更敏感P_ch_max 一般取聚合额定功率的 100%P_dis_max 取 30% 到 40%因为可削减负荷远小于可增加负荷E_min 不能设 0留 15% 到 20% 的底能量否则室温贴着死区上界跑回代容易越界c_deg 从 0.01 起调观察到 P_ch 曲线出现高频抖动就加大5. 复现校验三个把虚拟电池调稳的参数与两个常见误用5.1 用回代仿真验证 LP 解是否真的可行LP 里满足能量边界不代表聚合后每台设备都不越死区因为外近似放大了可行域。最便宜也最有效的验收手段是把最优功率轨迹回代到单机模型def replay_single_machine(vb_traj, params, e_init): 把聚合功率按占比摊到单机检查能量是否越界。 e, violations, worst e_init, 0, 0.0 for p in vb_traj: e params[a] * e - params[gain] params[b] * p low, high 0.0, params[E_cap] if e low - 1e-6 or e high 1e-6: violations 1 worst max(worst, max(low - e, e - high)) e min(max(e, low), high) return violations, worstviolations 超过总步数的 10%就说明外近似太松需要把 E_max 按采样点收紧或者调大 c_deg 抑制激进动作。worst 反映越界幅度的量级控制在小数点后两位以内一般可以接受。5.2 三个最敏感的参数参数敏感方向调参建议a保持系数24 步残余 a^24直接决定耦合强度Δt 定了就由 R·C 标定别硬调η_ch / η_dis效率取 1 会系统性低估能量损耗制冷聚合 η_ch 取 0.85 ~ 0.95E_min / E_max窗口越宽灵活性越大回代越容易越界留 15% 底能量上界按采样 99 分位a 的敏感度可以在定步长前先扫一遍把 R·C 从 1e4 到 1e5 拉一条曲线看 a 从 0.88 到 0.99 变化时同样的削峰目标下 P_dis_max 需求差多少。差得大说明当前步长选得不合适应该把 Δt 改小而不是硬凑参数。5.3 两个常见误用第一个是把功率边界和能量边界独立聚合。正确做法是先算聚合可行域的凸包再用虚拟电池包住它中间那步闵可夫斯基和的采样不能省。第二个是把自放电项当 1 处理。储能的日自放电常被忽略温控负荷不能这么干——温差驱动的热增益每步都在扣能量忽略它会让 LP 高估可用灵活性实盘执行时削峰量直接缩水三成以上。5.4 用峰谷差和购电成本做验收指标最后拿两个数字收口一个对比基准场景虚拟电池功率全置 0一个对比接入场景cost float(np.sum(price * grid)) peak, valley grid.max(), grid.min() print(f日购电成本 {cost:.1f} 元峰谷差 {peak - valley:.1f} kW)我一般会把这个脚本挂在参数扫描外面a、c_deg、E_min 三个变量各取三档做个 27 组的小网格只看峰谷差降幅和 violations 两个指标。降幅不到 5% 的组合直接丢掉剩下的再按成本排序。这套流程跑下来虚拟电池参数就从「论文里抄的」变成「自己数据标定的」接进日前出清才站得住。本文还有配套的精品资源点击获取