企业官网建设流程全解析

简介：面向能源工程研究者与智能电网开发者的技术资料，围绕电动汽车与温控负载（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(k+1) = T_in(k) + Δt/(R·C)·(T_out(k) − T_in(k)) − (Δt·η/C)·u(k)·P_rated

其中 u(k) 取 0 或 1，表示定频压缩机停或开。这套模型的参数在工程上是可以标定的，常见做法是用历史温度曲线做最小二乘拟合，而不是直接查铭牌。

注意：Δ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(k+1) = 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(k+1) = a·E(k) − gain + b·u(k)

其中三个系数分别定义为：

• a = 1 − Δt/(R·C)，保持系数，衡量一个步长内热量自然流失的比例
• b = Δt·η·P_rated / 3600，满功率运行一个步长注入的等效能量（Wh）
• gain = Δt·(T_out − T_max) / (3600·R)，环境热增益折算成能量（Wh）

对照通用虚拟电池形式 E(k+1) = 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/python（Windows 是 .venv\Scripts\python.exe），四件套依赖装齐：

python -m venv .venv source .venv/bin/activate # Windows: .venv\Scripts\activate pip install numpy pandas scipy matplotlib pulp

numpy 负责状态枚举和向量化，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(a=a, b=b, gain=gain, E_cap=E_cap, P_rated=P_rated) def reachable_energy(params, horizon, levels=80): """给定初始能量，枚举每个时段可达的能量格点集合。""" 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^24，a=0.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_samples=3000, seed=0): """抽取聚合 (p, E) 可行点。p 单位 kW，E 单位 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_units=200, horizon=24) P_min, P_max = P.min(axis=0), P.max(axis=0) E_min, E_max = E.min(axis=0), E.max(axis=0)

外近似必须用真正的 min/max 而不是分位数，否则边界会切掉真实可行点，回代时必然越界。想收紧边界的话，改用分位数之后再跑一次可行性回归：把采样点逐个代入边界，统计被切掉的比例，控制在 1% 以内再接受。

4. 虚拟电池接入日前优化调度：目标函数与约束的线性化落地

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(t+1) = 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), lowBound=0, upBound=P_trafo) P_ch = pulp.LpVariable.dicts("P_ch", range(T), lowBound=0) P_dis = pulp.LpVariable.dicts("P_dis", range(T), lowBound=0) E = pulp.LpVariable.dicts("E", range(T + 1), lowBound=0) 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(msg=False)) 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 的索引长度是 T+1，边界约束要覆盖 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, worst

violations 超过总步数的 10%，就说明外近似太松，需要把 E_max 按采样点收紧，或者调大 c_deg 抑制激进动作。worst 反映越界幅度的量级，控制在小数点后两位以内一般可以接受。

5.2 三个最敏感的参数

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% 的组合直接丢掉，剩下的再按成本排序。这套流程跑下来，虚拟电池参数就从「论文里抄的」变成「自己数据标定的」，接进日前出清才站得住。

本文还有配套的精品资源，点击获取