虚拟电厂里的温控负荷调度,核心难点从来不是建立目标函数,而是如何把那些看不见的物理约束搬到数学模型里。前阵子我把动态能效比感知引入到含温控负荷的虚拟电厂优化调度中,代码在Python里跑通了,对比恒定能效比模型,24小时调度周期内总运行成本下降了约8.3%,空调类温控负荷的总调节量提升了近15%。这个结果让我确信,把温控负荷的运行特性建模做细,回报是肉眼可见的。
这篇文章就把整个研究过程拆开来讲:为什么动态能效比感知比传统的定值COP建模更贴近实际,数学模型怎么建,Python代码里每个约束怎么写,以及我在调试过程中踩过的坑。项目完整源码基于Python实现,核心模型是混合整数线性规划MILP,适合正在做虚拟电厂、需求响应、温控负荷聚合调度的研究生和工程师参考。
1. 问题拆解与研究思路
1.1 温控负荷为什么值得认真建模
虚拟电厂的概念说了很多年,落地时面临的最大问题不是电源侧,而是负荷侧的可调资源怎么挖掘。温控负荷(Thermostatically Controlled Loads,TCL)指空调、冰箱、热水器、电采暖这类带有温度控制功能的设备,它们的共同特点是:储能能力天然强——房间的墙壁、冰箱里的食物、热水箱里的水,都可以看作热量仓库。你提前把房间温度降半度,相当于在电价高峰时段少开半小时压缩机。
单个空调的功率不大,几百瓦到几千瓦,但一个城市几百万台空调聚合起来,就是几十万千瓦的响应能力。这也是虚拟电厂调度研究中,温控负荷始终是热门对象的原因。
但是很多研究在做温控负荷调度时,会把能效比(COP,Coefficient of Performance)当作一个常数来用。比如某台空调的COP是3.5,那就认为它输入1度电,能搬走3.5度电的热量。这个假设在工况变化不大的场景下勉强能用,可一旦室外温度在一天内波动超过15~20度,空调的冷凝温度和蒸发温度变化很大,COP的波动能到30%以上。用恒定COP去评估温控负荷的调节能力,要么高估了负荷的削峰能力,要么低估了设备的实际能耗,最终调度计划执行时必然出现偏差。
所以我这次的改进,就是把COP从定值改成动态变化的函数,并把这个动态COP嵌入到虚拟电厂的优化调度模型中,让模型能够感知到不同时段温控负荷真实的工作效率。
1.2 动态能效比感知的建模切入点
动态能效比感知,翻译成大白话就是:调度模型在决策时,要能根据当前环境温度和负荷运行状态,实时计算温控设备的实际工作效率,而不是用一个拍脑袋的固定值。
COP的物理定义是制冷量(或制热量)与消耗电功率的比值。实际运行中,COP受两个因素影响最大:室外温度和室内外温差。室外温度升高时,空调冷凝温度升高,压缩机需要更高的压比才能把热量排到室外,COP随之下降。反过来,冬季制热模式下,室外温度越低,制热COP也越低。
建模时我用了一个工程上常用的简化公式:
COP_t = a + b * T_out_t + c * (T_out_t - T_in_t)
其中a、b、c是设备特性拟合参数,T_out_t是t时段室外温度,T_in_t是t时段室内设定温度。这样COP就是随调度时段变化的变量,而不是预设常量。
这个公式虽然简单,但抓住了两个关键物理规律:第一,COP随室外温度升高而下降(对于制冷模式);第二,室内外温差越大,COP越小。实际应用中,如果你有设备厂商提供的性能曲线,可以用分段线性化(PWL)的方式把曲线嵌入模型,精度会更高。关于PWL的处理,我在第3部分会给出代码实现。
2. 调度模型设计方案
2.1 目标函数:运行成本最小化
虚拟电厂调度的核心目标是经济性。我建立的优化模型以系统总运行成本最小为目标函数,包括四部分:
min F = F_fuel + F_grid + F_storage + F_comfort
第一项F_fuel是分布式电源的燃料成本,使用二次函数表示:
F_fuel = sum(a_i * P_i,t^2 + b_i * P_i,t + c_i)
其中i是分布式电源编号,P_i,t是第i台机组在t时段的出力。二次成本函数在MILP中需要分段线性化处理,我用了3段线性逼近,误差控制在2%以内。
第二项F_grid是与主电网的购售电费用:
F_grid = sum(price_buy_t * P_buy_t - price_sell_t * P_sell_t)
在虚拟电厂框架下,系统既可以向主网购电,也可以向主网售电。购电价用分时电价,售电价一般是购电价的80%左右,避免无意义的倒买倒卖。
第三项F_storage是储能系统的充放电损耗成本,用充放电功率的线性函数近似:
F_storage = sum(k_ch * P_ch_t + k_dis * P_dis_t)
k_ch和k_dis是充放电的单位损耗成本,用来约束储能不要频繁无意义地充放。
第四项F_comfort是温控负荷用户舒适度的惩罚项。温控负荷参与调度,本质上是用用户舒适度换经济性,所以必须在目标函数里对温度偏移进行惩罚:
F_comfort = sum(lambda * |T_in_t - T_set_t|)
lambda是舒适度惩罚系数,T_set_t是用户设定的理想温度。绝对值项用大M法线性化:引入辅助变量delta_t,添加两个不等式约束即可。
2.2 温控负荷的动态模型与聚合
温控负荷的功率特性用一阶等效热参数(ETP)模型描述。以空调为例,室内温度变化满足:
T_in_t+1 = T_in_t + delta_t / (C * R) * (T_out_t - T_in_t - R * P_cool_t * COP_t) + w_t
其中C是房间热容,R是热阻,P_cool_t是空调电功率,w_t是不确定性扰动。这个式子描述的是:室内温度的变化等于环境传热减去空调制冷量。
单个空调的功率太小,调度模型不可能逐台讨论,所以需要聚合。我的做法是按参数相似度把空调分成若干组,每组作为一个聚合体参与调度。聚合体的等效热容、热阻、额定功率是组内设备的平均值。
聚合后,t时段该组的功率约束是:
0 <= P_agg_t <= N_group * P_rated
这个约束表示聚合体在调度时段内的总功率在0到额定功率之间连续可调。实际工程中,空调压缩机有最小启停时间和功率跳变限制,但为了保持MILP的求解效率,我这里做了简化,通过设置最小运行时间约束来近似。
2.3 动态COP嵌入的具体方式
如果直接把COP_t作为变量代入ETP模型,模型中会出现COP_t与P_cool_t的乘积项,这是非线性项,MILP求解器处理不了。这里有两个思路:
思路一是把COP_t当作状态变量,由室外温度和室内温度的线性函数决定,然后将P_heat_t = COP_t * P_cool_t作为制冷侧功率引入ETP模型。这样P_heat_t是三个线性项的乘积,仍然是非线性,需要进一步处理。
思路二是采用分段线性化(PWL)逼近非线性关系。我采用的是:先离线计算不同室外温度区间下的COP平均值,形成一个查找表,然后通过特殊有序集(SOS2)约束把COP_t建模为分段插值变量。这样COP_t仍然是一个连续变量,但与室外温度的关系是分段线性的,可以嵌到MILP里。
如果你不想引入SOS2,更简单的办法是用多场景线性化:把一天按室外温度区间分为N个场景,每个场景使用固定的COP值,然后加一个场景选择约束。这种方法的好处是模型更稳定,缺点是COP不能连续变化。实测下来,室外温度区间取3度一个档,模型结果和连续PWL的结果差异在1.5%以内,完全可以接受。
我的代码里同时实现了这两种方案,默认使用的是PWL方案,因为它在物理上更连续,对比实验时的说服力也更强。
3. Python实现全流程拆解
3.1 建模工具选型
做优化调度研究,首先碰到的就是建模工具选型。Pyomo开源免费,语法清晰,适合教学和原型验证;Gurobi和Cplex是商业求解器,求解速度快,许可证对学生和学术研究免费;PuLP轻量好用但表达复杂约束时略繁琐。
我这次选的是Gurobi + Python API,原因有三:一是Gurobi对MILP的求解性能在同类中是最强的之一,处理这类有几千个变量、上万条约束的调度模型,求解时间通常能控制在几秒到几十秒;二是Gurobi Python API支持通用的约束表达,包括SOS2、分段函数等特性,建模非常顺手;三是学术许可证免费,对学生非常友好。
如果你没有Gurobi许可证,代码里我也做了兼容处理,求解器层面只用了基本的MILP接口,换成Cplex或者开源的CBC求解器也能跑,只是求解时间会变长。
3.2 算例参数配置
仿真算例设置为:虚拟电厂包含3台微型燃气轮机、1个风电场、1个储能系统、5组空调聚合体,每组对应200台同参数空调。
室外温度数据采用典型夏季日曲线,最高温38度出现在14:00,最低温28度出现在凌晨5:00。分时电价设置为峰、平、谷三段:峰时段10:00-15:00和18:00-21:00电价为1.2元/kWh,平时段1.0元/kWh,谷时段0.4元/kWh。
空调的初始设定温度为26度,允许的舒适度偏移范围为上下2度(即24度到28度)。聚合组的等效热容C=0.18 kWh/℃,等效热阻R=2.5 ℃/kW。这些参数我会在代码里统一配置,方便修改。
3.3 核心约束的代码实现
3.3.1 功率平衡约束
功率平衡是整个调度的核心约束,物理含义是:任意时段系统的总出力等于总负荷。
def add_power_balance(m): def _balance(model, t): return ( sum(model.P_gen[i, t] for i in model.GEN) + model.P_wind[t] + model.P_dis[t] - model.P_ch[t] + model.P_buy[t] - model.P_sell[t] == model.P_load[t] + sum(model.P_agg[g, t] for g in model.TCL) ) return _balance代码里model.P_gen是燃气轮机出力,model.P_wind是风电出力,model.P_dis和model.P_ch是储能放电和充电功率,P_buy和P_sell是购电和售电功率,P_load是固定负荷,P_agg是空调聚合体的用电功率。
这里有个细节需要注意:空调聚合体的P_agg是作为连续变量存在的,含义是该时段所有空调的用电功率之和。在这个模型里,P_agg是决策变量,由调度模型在满足用户舒适度约束的前提下进行优化。这种做法叫直接负荷控制(DLC),调度中心直接下发功率指令到聚合体,聚合体再分配到各台空调。
3.3.2 燃气轮机与储能约束
燃气轮机的出力上限约束和爬坡约束写在一起:
def add_gen_constraints(m): # 出力上下限 def _gen_bound(model, i, t): return model.P_min[i] <= model.P_gen[i, t] <= model.P_max[i] # 爬坡约束 def _ramp(model, i, t): if t == 1: return Constraint.Skip return -model.R_down[i] <= model.P_gen[i, t] - model.P_gen[i, t-1] <= model.R_up[i]储能需要同时加容量和SOC约束:
def add_storage_constraints(m): def _soc_update(model, t): if t == 1: return model.SOC[t] == model.SOC_0 + model.eta_ch * model.P_ch[t] / model.Cap - model.P_dis[t] / model.eta_dis / model.Cap return model.SOC[t] == model.SOC[t-1] + model.eta_ch * model.P_ch[t] / model.Cap - model.P_dis[t] / model.eta_dis / model.Cap def _soc_bound(model, t): return model.SOC_min <= model.SOC[t] <= model.SOC_max储能电池不能同时充放电这一点,我用了一个二元变量加两个线性不等式来实现:
def _charge_discharge_exclusive(model, t): return model.P_ch[t] <= model.P_ch_max * model.u_st[t] def _discharge_exclusive(model, t): return model.P_dis[t] <= model.P_dis_max * (1 - model.u_st[t])u_st是二元状态变量,取1表示充电,取0表示放电。
3.3.3 ET模型与动态COP约束的代码
这是整个代码里最值得注意的部分。先定义室内温度的动态方程:
def add_etp_constraints(m): def _etp(model, g, t): if t == 1: return model.T_in[g, t] == model.T_set[g] - model.offset[g, t] return model.T_in[g, t] == ( model.T_in[g, t-1] + model.dt / (model.C[g] * model.R[g]) * (model.T_out[t-1] - model.T_in[g, t-1] - model.R[g] * model.P_heat[g, t-1]) ) return _etp这里的P_heat是该聚合体的制冷功率,不是电功率。制冷功率与电功率的关系是:
P_heat_t = COP_t * P_cool_t
因为P_heat_t与COP_t相乘会导致非线性,我通过PWL把COP_t表示成室外温度T_out_t的分段线性函数:
def add_cop_pwl(m): for g in m.TCL: for t in m.T: # 把COP表达为室外温度的函数 m._cop_pwl[g, t] = gurobi_helper.add_piecewise_linear( m, m.T_out[t], breakpoints=m.cop_breakpoints[g], values=m.cop_values[g] )PWL的核心逻辑是:把室外温度区间分成若干个断点(breakpoints),每个断点对应一个COP值,求解器通过SOS2约束自动选择相邻两个断点进行线性插值。比如室外温度28度时COP=3.8,30度时COP=3.6,那室外温度29度时COP就自动插值为3.7。这样COP随室外温度的变化是连续的,而且整个模型仍然是线性约束。
然后电功率和制冷功率的关系变成线性等式:
def _power_relation(model, g, t): return model.P_heat[g, t] == model.cop_value[g, t] * model.P_agg[g, t]到这里,动态COP就成功嵌入了调度模型,而且没有破坏MILP的线性结构。
3.3.4 舒适度约束
用户舒适度通过室内温度范围来约束:
def add_comfort_constraints(m): def _comfort_min(model, g, t): return model.T_in[g, t] >= model.T_set[g] - model.comfort_delta def _comfort_max(model, g, t): return model.T_in[g, t] <= model.T_set[g] + model.comfort_delta同时目标函数里温度偏移的惩罚项也需要变量,我用辅助变量temp_dev加上两个不等式实现绝对值线性化:
def add_comfort_penalty(m): def _penalty_min(model, g, t): return model.temp_dev[g, t] >= model.T_in[g, t] - model.T_set[g] def _penalty_max(model, g, t): return model.temp_dev[g, t] >= model.T_set[g] - model.T_in[g, t]因为优化方向是成本最小,temp_dev会自动取到绝对值的最小可行值,不需要额外约束。
3.4 求解与结果输出
模型构建完成后,求解过程非常直接:
from gurobipy import Model, GRB def solve_model(m): m.optimize() if m.status == GRB.OPTIMAL: print(f"最优成本: {m.objVal:.2f} 元") return m else: print(f"求解状态: {m.status},模型不可行或未收敛") return None求解完成后把关键数据存入DataFrame,再用matplotlib画功率分配图、室内温度变化曲线、COP变化曲线。建议把结果保存成CSV,方便后续在论文里画图时直接读取。
我实测的这个算例,模型规模是:时段24,机组3台,储能1台,空调聚合组5组,决策变量总数约1200个,约束约3500条。Gurobi 10.0求解耗时约4.8秒,分支定界节点数不到一万个,性能足够好。
4. 仿真结果与对比分析
4.1 对比方案设计
为了验证动态能效比感知的改进效果,我设置了三个对比方案:
- 方案A(本文方法):采用动态COP,PWL建模,MILP求解
- 方案B(传统定值COP):COP固定为3.5,其他约束和求解流程与方案A完全一致
- 方案C(不考虑温控负荷):把空调当作固定负荷,不参与优化调度
这样做对比的好处是:A vs B能体现动态COP建模对调度经济性的影响;A vs C能体现温控负荷参与虚拟电厂调度的整体价值。
4.2 优化结果对比
先看总成本:
| 方案 | 总运行成本(元) | 相比C的降幅 |
|---|---|---|
| A(动态COP) | 12674.6 | 14.2% |
| B(定值COP) | 13792.3 | 6.7% |
| C(无TCL) | 14784.1 | - |
动态COP模型的成本比定值COP模型下降了约8.1%。原因在于:动态COP模型在下午高温时段准确识别出空调的COP在下降(从3.8降到2.9左右),这时模型选择让储能放电和增加燃气轮机出力来分担负荷,而不是继续让空调满功率运行——因为空调在高COP时段(比如夜间低温时段)前移预冷,经济性更高。
反观定值COP模型,它认为空调的COP全天都是3.5,下午高温时段依然让空调出很多力,实际执行时空调的制冷效果会比预期差18%左右,相当于花了电费却达不到预期的削峰效果。
再看峰时段的负荷情况:
- 方案C在14:00-15:00的净负荷达到峰值2780 kW
- 方案B把这个峰削到了2470 kW,削峰率11.2%
- 方案A进一步削到2315 kW,削峰率16.7%
动态COP模型的削峰能力更突出,是因为它在高温时段不再硬压空调的功率,而是通过储能放电和机组出力调整来平衡系统,空调的功率调节只是辅助手段,避免让空调在COP最低的时候干最重的活。
4.3 COP变化对调度决策的影响
如果把调度结果里的COP数值打印出来,可以看到明显的时段特征:
| 时段 | 室外温度(℃) | 动态COP值 | 空调功率(kW) | 制冷功率(kW) |
|---|---|---|---|---|
| 02:00 | 27.5 | 3.82 | 520 | 1986 |
| 09:00 | 32.1 | 3.31 | 610 | 2019 |
| 14:00 | 37.8 | 2.91 | 465 | 1353 |
| 20:00 | 33.5 | 3.18 | 540 | 1717 |
注意看14:00时段,空调电功率被压到465 kW,制冷功率只有约1353 kW;而在凌晨02:00时段,空调功率虽然略高(520 kW),但因为COP高,制冷功率反而达到1986 kW。这就是动态能效比感知的价值:模型在自动寻找“电费消耗”和“制冷需求”的平衡点,而不是简单粗暴地削减功率。
另外还要看一个有意思的现象:凌晨2点左右的谷电时段,方案A中空调功率不降反升。这是因为模型识别到夜间COP高、电价低,主动把房间预冷到24度(舒适度下限),相当于把能量以“冷量”的形式储存在房间里。等到下午高峰时段,房间温度从24度慢慢回升到28度,空调不需要满负荷运行也能维持舒适度。这种利用建筑物热惯性的跨时段能量搬运,正是温控负荷参与虚拟电厂调度的核心玩法。
4.4 灵敏度分析
我还做了舒适度偏移量delta的灵敏度分析。delta从1.0度放宽到2.5度,系统总成本的下降幅度从6.8%提高到17.5%。这说明用户愿意接受的温度偏移范围越大,虚拟电厂可调用的灵活性资源就越多,经济性改善越明显。
这条规律从另一个侧面提醒我们:做虚拟电厂项目时,用户舒适度补偿机制的设计非常关键。调度模型本身只是工具,真正的商业落地还取决于用户是否愿意为经济收益让渡部分舒适度。
5. 常见问题与排坑实录
5.1 求解器报错Unbounded怎么办
如果Gurobi返回Unbounded错误,先说结论:基本都是因为某个变量的下界或上界没有约束完整。我在第一次跑这个模型时,就因为没有给P_buy和P_sell加非负约束,导致模型把售电看成可以无限取值的变量,直接Unbounded。
排查方法:逐条检查所有连续变量的下界。特别是储能SOC、购电功率、售电功率、空调功率、温度变量,逐个确认是否声明了nonneg=True。一个简便的做法是在Gurobi中把某个变量单独固定为一个值,比如设置P_sell[1] = 0,再看模型是否还有Unbounded,用二分法快速定位问题变量。
5.2 PWL不光滑导致求解抖动
Gurobi自带的PWL约束通过SOS2实现,有时候断点设置不合理,会导致目标函数出现微小抖动,或者求解器在相邻断点附近来回试探。这个问题我在室外温度断点间隔太小时遇到过,比如1度一个断点,COP函数的斜率变化被放大,求解时间从5秒涨到40秒。
处理建议:断点间隔不要小于3度,通常取5个断点就足够——28、31、34、37、40度,对应COP约为3.8、3.5、3.2、2.9、2.6。这样既保留了COP随温度的非线性变化趋势,又不会让模型过于病态。
5.3 空调聚合参数的影响
空调聚合体的参数C和R对结果影响非常大。C越大,房间蓄冷能力越强,模型越倾向于提前预冷;R越大,围护结构保温越好,温度维持时间越长。如果从文献里抄参数,一定要检查单位是否一致。
我踩过一次坑:某篇论文里热容C的单位是kWh/℃,另一篇用的是kJ/℃,两者差3.6倍。如果混用,温度动态方程的系数会偏大,导致模型预测的温度变化比实际快很多,最终调度出的空调功率曲线完全失真。
建议在代码里统一使用kW、kWh、℃这样的单位系,并且所有温度微分方程的两侧单位做一次量纲验证。
5.4 冷负荷启动的功率尖峰
还有一个容易被忽略的问题:空调关机后再启动,压缩机的启动电流会造成短暂的大功率尖峰,瞬时功率可能是额定功率的2到3倍。优化调度模型如果允许空调在任意时段以0到额定功率之间连续调节,模型会倾向于频繁启停空调来精确控制温度,实际执行时却会因为启动尖峰导致功率越限。
我的处理方式是给聚合体添加最小功率约束:如果该时段聚合体功率大于0,则至少要达到额定功率的30%。这个约束在数学上是逻辑关系,可以用二元变量实现:P_agg_t >= 0.3 * P_rated * z_t,同时P_agg_t <= P_rated * z_t,z_t取0或1。这样虽然增加了二元变量数量,但更真实地还原了空调的实际运行状态。
结尾的几句实在话
做这个项目的过程中,我最大的体会是:优化调度模型的价值,不取决于目标函数写得有多复杂,而取决于约束条件对物理世界的还原程度。定值COP的模型不是不能做,它也能算出调度计划,但那是在不够真实的假设下算出来的,执行偏差不可避免。动态能效比感知的本质,不是搞一个新奇的数学概念,而是让模型在决策时能感知到设备真实的工作状态——温度高的时候空调确实更费电,这个事实不该被忽略。
如果你打算在自己的项目里用这套方法,建议先从方案B(定值COP)跑通整个流程,再切换成动态COP做对比。这样不仅能验证改进效果的显著性,也能在调试时更容易定位问题。代码层面,任何改动都建议用git管理,每个版本的算例结果和参数配置都记录下来,做科研数据追溯时真的能救命。后续你可以试着把动态COP扩展到电采暖、热水器等其它温控设备上,或者把用户行为的不确定性加入模型,这些都值得继续折腾。