简介:围绕电力系统日前优化调度中的需求侧资源灵活性刻画,资源以论文复现形式提供基于虚拟电池(VB)模型的完整Python实现,面向能源工程研究者、智能电网从业者及具备中级编程能力的技术人员。压缩包共1个docx文档,大小24KB,文档按参数设置、EV单辆与集群VB建模、HVAC建模、pulp线性规划调度求解等模块组织,并附代码逐段解释与结果图形。已有72人学习下载。读者可掌握从单个设备到集群的虚拟电池模型转换方法,理解日前电价、备用容量价格及EV/HVAC运行约束如何在Python中表达,并复现日前计划策略,观察需求响应资源参与后负荷曲线与运行成本的变化。文中数据为理论模拟样本,适合作为方法验证、教学演示及后续研究的基础;在更大规模或真实环境部署前,仍需进一步打磨模型细节。
1. 虚拟电池模型实操:从公式到可跑的调度代码
做电力系统优化的人大多遇到过这类尴尬:论文里把电动汽车和空调的灵活性写得头头是道,公式推了一整页,但到了自己动手建模时,连“这个变量到底是边界还是状态”都要纠结半天。这篇论文复现笔记直接聚焦虚拟电池(Virtual Battery, VB)模型在日前调度中的完整落地过程,从单辆EV的功率轨迹推导到集群聚合、再到基于pulp的线性规划求解,全部用Python实现并给出可运行代码。适合已经有电力系统优化基础、想快速把需求侧资源灵活性模型落到工程项目里的工程师,也适合正在复现论文、需要参照代码做二次开发的研究者。文中会明确标出哪些参数是示例值、哪些公式做了简化处理,避免你把教学代码当成生产级解决方案。
2. 单辆EV的VB模型推导:为什么电能边界是时域耦合的
2.1 三个时间参数决定充放电能力
EV的灵活性本质上由并网时段、电池容量和充放电功率共同决定。single_ev_vb_model函数接受T_i_o(并网时段)和T_i_d(离网时段),这两个参数直接划定了EV参与调度的可行域。代码中每个时段计算电能上限时用了min((E_i_To + P_i_max * (t - T_i_o) * time_step), E_i_max),这里E_i_To是并网初始能量,P_i_max * (t - T_i_o) * time_step表示从并网到现在以最大功率充电累积的能量增量,整体取min是因为物理上电池容量存在硬上限。
很多初学者会把E_i_max当成一个常数约束直接写进优化模型,但在VB框架里它需要通过逐时段的min操作转化为动态边界,两种做法在数学上等价,但后者可以直接嵌入调度模型的约束矩阵,省掉额外的状态变量。
注意P_i_min是负值,代表放电功率上限。E_i_t_EV_min[t] = max((E_i_ex + P_i_min * (T_i_d - t) * time_step), E_i_min)这行的物理含义更关键:为了保证离网时能量能达到E_i_ex,从当前时段到离网时段这段时间内,最保守的放电策略是让能量始终不低于E_i_ex减去剩余时间内的最大放电量。这种“倒推”约束是VB模型区别于普通储能模型的核心。
2.2 功率限制的隐式耦合:能不能充,取决于下一时段
函数里有一段不太显眼的逻辑:P_i_t_EV_max[t] = min((E_i_t_EV_max[t + 1] - E_i_t_EV_min[t]) / time_step, P_i_max)。这行代码做的是把电能上下界转换为功率上下界,但转换方式有讲究。如果直接从E_i_max - E_i_min除以时间步长,得到一个宽泛的功率范围,但VB模型要求逐时段考虑:本时段的最大功率充入后,下一时段的电能上界不能越限,所以用E_i_t_EV_max[t + 1] - E_i_t_EV_min[t]作为动态约束。
这也解释了为什么论文里强调“功率限制不是常数”。EV的充放电功率上限虽然由硬件决定,但在调度模型中实际可用功率还受到当前能量状态下可充可放空间的限制。代码中P_i_t_EV_min[t] = max((E_i_t_EV_min[t + 1] - E_i_t_EV_max[t]) / time_step, P_i_min)同理,只是换成了下界推导。如果跳过这个时域耦合直接建模,优化的结果很可能出现“功率设备约束满足、但电池能量连续性破坏”的违规解。
2.3 函数边界条件与索引陷阱
代码里for t in range(T)循环内对t < T_i_d - 1的检查并不影响函数内部的电能计算,但功率限制计算依赖t+1时段的电能边界,如果t = T - 1时会索引越界。因此在实现时,我给函数加了前置检查:assert T_i_d <= T,并且最后一段功率不做计算直接置零,避免数组访问异常。这种边界条件在论文里通常一两句话带过,但落地时最容易翻车。
提示:输入
T_i_o、T_i_d一般给整数时段编号(0到95),不要传时间戳或datetime对象,函数内部用了大量整数索引运算。
3. HVAC的VB建模:热力学参数的辨识与简化
3.1 一阶等效热参数模型:a和b从哪里来
single_hvac_vb_model函数里用了a和b两个系数,对应论文中的alpha = exp(-time_step / (R * C))和beta = R * (1 - alpha),其中R是等效热阻,C是等效热容。如果直接按热力学公式推导,HVAC的室内温度演化满足theta_t = alpha * theta_{t-1} + (1 - alpha) * (theta_ambient - R * P_cooling),当制冷功率为0时室温趋近环境温度,制冷功率开启时室温下降。
代码里没有显式计算R和C,而是直接把a、b作为输入参数传入。theta_t = theta_t_prev + a * (theta_0 - theta_t_prev) - b * P_HVAC_max * (P_t_HVACVB[t - 1] > 0)这行其实是一个简化:b * P_HVAC_max合并了热阻、COP和功率的乘积。实际项目中我会先用最小二乘法辨识a和b:给一台空调给定开启/关闭序列,记录室内温度变化曲线,然后用np.linalg.lstsq拟合出最优a和b。
3.2 温度死区与功率状态的强耦合
HVAC建模的难点在于温度死区(epsilon)的存在使得功率不再是连续决策变量。if theta_t < theta_c - epsilon / 2: P_t_HVACVB[t] = P_HVAC_max表示温度低于舒适下限时强制满功率制冷——注意这里逻辑是反的,制冷时功率为正表示消耗电功率,温度应该下降,如果温度已经低于死区下限,说明需要减少制冷量,但代码里写成了满功率,这是一个典型的伪代码矛盾点。
论文原文的处理方式是通过优化求解器在死区内自由选择功率值,即不是在每个时段硬性判定开或关,而是把温度约束转化为不等式的约束条件。我在复现时采用了另一种方式:用连续变量表示功率,在死区内允许0到P_max之间的任意值,温度演化按theta_t = alpha * theta_{t-1} + (1 - alpha) * (theta_ambient - R * P_cooling[t])线性计算,这样就不需要if判断,整体可以嵌入线性规划框架。
这种连续松弛的代价是损失了设备开关状态的物理真实性,但换来的是pulp可以直接求解。对于集群级调度来说,单台设备的开关状态精度无关紧要,集群层面的聚合功率才是决策变量关注的对象。
3.3 电能边界与基准功率的选取
代码里P_t_base = (theta_0 - theta_c) / (hvac_cop * 1)计算的是维持设定温度所需的基准功率,E_t_HVACVB[t] = (1 - a) * E_t_HVACVB[t - 1] + (P_t_HVACVB[t] - P_t_base) * time_step是在迭代计算相对于基准功率的能量偏移。E_t_HVACVB不是电池的绝对能量,而是一个虚拟的能量偏移量,可正可负,代表热储能相对于设定点的盈余或亏欠。
这里的hvac_cop * 1中的热阻R被硬编码为1,实际项目中需要根据房间面积、墙体材料等参数估算,建议用R = 1 / (UA)计算,其中UA是建筑的围护结构传热系数。如果R给得太小,P_t_base会偏大,虚拟电池的容量边界计算失真。
4. 集群VB模型与聚合边界:从单设备到虚拟电厂的跨越
4.1 集群EV的累加逻辑与并网离网修正项
cluster_ev_vb_model函数做的事非常简单:把每辆EV的功率上下界逐时段累加,得到集群的可调功率范围。但E_t_EVVB_exchange这个变量的计算容易被忽略——E_t_EVVB_exchange[t] += E_i_t_EV_max[t + 1] * (W_i_t_plus_1 ** 2 - W_i_t ** 2)。
这行代码的实际意义是:EV在并网瞬间(W从0变为1)和离网瞬间(W从1变为0),其电池能量相对于VB模型的等效能量发生跳变。物理上,EV插枪瞬间电池以当前SOC接入系统,这对集群的能量连续性是一种“注入”或“抽离”。数学上加W_i_t_plus_1 ** 2 - W_i_t ** 2是为了处理边界时刻的符号:接入时刻产生正跳变,离开时刻产生负跳变,平方操作让两个时刻都有定义。
项目中我通常把E_t_EVVB_exchange作为调度模型中的外部能量注入项写进约束:E_t_VB[t] = E_t_VB[t-1] + (P_t_VB[t] - P_t_base[t]) * time_step + E_t_EVVB_exchange[t],这样能量连续性约束才能自洽。
4.2 聚合边界的松弛问题
把200辆EV、200台HVAC直接线性累加得到的集群边界,在数学上是可行的,但会忽略一个关键事实:单设备的功率限制在聚合层面不是同步发生的。比如200辆EV的充电功率上限都是6.5kW,但同一时段内所有车都同时满充的可能性很低,因为每辆车的SOC、并网时段、离网时段各不相同。
TB模型论文里通常用闵可夫斯基求和来收紧聚合边界,代码仓库里没有实现这一步,只用简单的逐点求和,这会导致调度结果偏保守——系统在备用容量申报时用了过低的灵活性上限。如果要追求更经济的调度结果,建议在集群模型中加入一个容量系数,按同时率经验值(EV取0.7到0.9,HVAC取0.8到1.0)对峰值功率做折减。
4.3 集群HVAC的优化边界计算:简化与实用的折中
cluster_hvac_vb_model里用scipy.optimize.minimize的伪代码只给出了框架,实际跑起来会有问题。目标函数只考虑了温度偏差平方和,没有加入功率变化惩罚,优化结果容易出现相邻时段功率剧烈波动。而且power_constraint返回的是布尔值而不是可微的约束向量,如果直接用于序列最小二乘规划求解器,可能无法满足约束的雅可比矩阵要求。
我在复现时改用以下方案:把温度演化方程写成线性等式约束嵌入pulp模型,决策变量是每台HVAC每时段的功率值,目标函数是温度偏差的绝对值之和加上功率变化惩罚。这个模型规模是:num_hvac * T个变量,num_hvac * T个温度约束,对pulp来说完全可解。初始可行解可以全部置为P_t_base,让求解器自行优化。
4.4 为什么基荷功率不能直接相加
如果每台HVAC的基准功率不同(因为热阻R和设定温度不同),集群基准功率应该是逐时段求和。但代码示例里用了[20] * num_hvac作为环境温度、[22] * num_hvac作为设定温度,如果所有HVAC参数完全一致,逐时段基准功率就是常数,调度模型会得到对称解。真实场景中环境温度随时间变化,不同的朝向、楼层导致热负荷曲线不同,基准功率会随时间波动,VB模型的能量边界也会跟着偏移。聚合时建议先算每台设备的基准功率曲线再求和,不要先平均参数再算基准功率,因为热力学模型对参数均值的非线性使得两者结果不一致。
5. 日前优化调度模型:pulp的建模细节与约束处理
5.1 决策变量范围与VB边界的映射
LpVariable.dicts("P_t_EVVB", range(T), lowBound=P_t_EVVB_min, upBound=P_t_EVVB_max)这段代码直观地把VB模型的功率上下界映射为变量的上下界。但需要注意,lowBound和upBound传入的是numpy数组,pulp能正常处理按数组下标取值的边界吗?实际验证结果是:pulp对数组型的边界支持不友好,运行时会报类型错误。
正确的做法是逐个时段创建变量:LpVariable(f"P_t_EVVB_{t}", lowBound=P_t_EVVB_min[t], upBound=P_t_EVVB_max[t])。这个细节处理虽然啰嗦,但能避免pulp内部类型混用导致的不稳定问题。变量数量越大,用字典推导式构造的效率优势越明显。
5.2 旋转备用容量约束:备用申报值如何进入目标函数
目标函数里R_t[t] * reserve_price[t]表示备用容量的收益,但代码中的约束prob += R_t[t] <= 200只是硬编码的常数上限,没有体现“备用可用性”的物理约束。现实中申报的备用容量必须能在系统需要时被调用,也就是说EV和HVAC在当前调度点附近应有足够的可调空间。
在复现时我给约束模型增加了如下关系:P_t_EVVB[t] + R_t_EV[t] <= P_t_EVVB_max[t]和P_t_EVVB[t] - R_t_EV[t] >= P_t_EVVB_min[t],其中R_t_EV[t]是EV提供的备用容量。这样备用申报不会超出VB模型的功率可行域,避免了系统在调用备用时越限的隐性风险。R_t在代码里是单独变量,建议拆分为R_t_EV和R_t_HVAC两部分,各自受对应设备集群的边界约束,再加总得到总备用。
5.3 功率平衡约束中的符号约定
目标函数(P_t_EVVB[t] + P_t_HVACVB[t] + P_load[t] - P_pv[t]) * day_ahead_price[t]把EV和HVAC的功率视为购电成本项,P_load是固定负荷,P_pv是光伏出力。这个模型里EV和HVAC都是以消耗功率为正方向,如果EV处于放电状态(功率为负),实际上是在向电网卖电,目标函数会自动把它算作负成本,即收益,这在数学表达上是对的。
但约束P_t_EVVB[t] + P_t_HVACVB[t] + P_load[t] - P_pv[t] <= P_limit要求净负荷不超过配变容量限值。这里隐含的假设是光伏出力和EV放电可以抵消负荷,代码里没有加入弃光约束。如果光伏出力大于负荷加充电功率,这个不等式天然满足,但实际工程中需要允许一定比例的弃光率,在模型中应加P_pv_curtail[t] >= 0和一个弃光惩罚项。
5.4 代码版本兼容性检查
pulp的API在2.7版本之后把LpProblem的名称参数改为关键字形式,prob = LpProblem("Day_Ahead_Optimization", LpMaximize)在旧版本能用,但新版建议改成prob = LpProblem(name="Day_Ahead_Optimization", sense=LpMaximize)。另外pulp库默认使用CBC求解器,线性规划模型规模不大时性能没问题,但整数变量和二次目标会让求解变慢,建议纯线性目标搭配连续变量。
6. 验证与踩坑记录:算例跑通后的五个高频问题
6.1 画图前先检查能量守恒
可视化结果之前,先验证VB模型内部的能量连续性。计算E_t_VB[t]的数值应该等于上一时段能量加上充放电量再减去基准功率的累积量,误差超过1e-6说明积分步长或索引计算有误。常见错误是time_step没有转换成小时导致能量单位不匹配——充电功率6.5kW乘以15分钟应得1.625kWh,如果直接用6.5 * 15会得到97.5kWh的错误结果。建议所有功率乘以时间都通过time_step换算,并在函数入口打印能量总量比对一次。
6.2 pulpy的LpAffineExpression不支持numpy数组广播
构建目标函数时如果用objective += np.sum(vars * prices)会直接抛异常,因为pulp的表达式对象不支持与numpy数组的逐元素乘法广播。正确写法是循环累加:objective += pulp.lpSum([vars[t] * prices[t] for t in range(T)])。遇到大模型必须用lpSum构造迭代器,直接sum()也可用但效率略低。
6.3 伪代码的功率约束会返回布尔值导致求解失败
power_constraint函数返回True或False,在pulp里这种约束会报类型错误。正确做法是把这个约束转化为pulp.LpAffineExpression不等式:prob += (P[t * num_hvac + z] <= P_HVAC_max_list[z]),逐台设备、逐时段地添加。这个循环可能产生数千条约束,但CBC求解器处理线性约束的效率完全可以应对。
6.4 不要直接使用论文的固定温度死区值做人舒适约束
hvac_temperature_deadband = [-0.5, 0.5]如果直接作为约束边界写入模型,温度会在死区边缘反复振荡,导致功率变量也剧烈跳变。更符合实际的做法是把死区宽度作为软约束,目标函数中加入温度越限惩罚项。设定温度theta_c作为决策变量或调整为恒温器设定值,允许优化模型在舒适范围内调整,这样调度结果更平滑且费用更低。
6.5 集群聚合前各设备的时标必须对齐
多辆EV的并网时段可能不同,集群模型按固定时段索引累加时,如果某辆EV的T_i_d小于当前时段t,该设备的功率边界对集群无贡献,但函数里没有自动跳过。建议在single_ev_vb_model内部循环前判断t < T_i_o or t >= T_i_d时所有输出置零,避免在聚合时把已离网设备的边界错误累加进来。没有这个判断,集群总功率上限在凌晨时可能被离网车辆抬高几十千瓦。
从那以后,我每次跑调度优化前都会强制走一遍能量守恒检查和设备在线状态检查,这两个小习惯能拦住80%以上的建模错误。希望这份笔记能让你在复现VB模型时少踩几个坑。
本文还有配套的精品资源,点击获取