虚拟电厂温控负荷调度：动态能效比感知的MILP建模与Python实现

虚拟电厂里的温控负荷调度核心难点从来不是建立目标函数而是如何把那些看不见的物理约束搬到数学模型里。前阵子我把动态能效比感知引入到含温控负荷的虚拟电厂优化调度中代码在Python里跑通了对比恒定能效比模型24小时调度周期内总运行成本下降了约8.3%空调类温控负荷的总调节量提升了近15%。这个结果让我确信把温控负荷的运行特性建模做细回报是肉眼可见的。这篇文章就把整个研究过程拆开来讲为什么动态能效比感知比传统的定值COP建模更贴近实际数学模型怎么建Python代码里每个约束怎么写以及我在调试过程中踩过的坑。项目完整源码基于Python实现核心模型是混合整数线性规划MILP适合正在做虚拟电厂、需求响应、温控负荷聚合调度的研究生和工程师参考。1. 问题拆解与研究思路1.1 温控负荷为什么值得认真建模虚拟电厂的概念说了很多年落地时面临的最大问题不是电源侧而是负荷侧的可调资源怎么挖掘。温控负荷Thermostatically Controlled LoadsTCL指空调、冰箱、热水器、电采暖这类带有温度控制功能的设备它们的共同特点是储能能力天然强——房间的墙壁、冰箱里的食物、热水箱里的水都可以看作热量仓库。你提前把房间温度降半度相当于在电价高峰时段少开半小时压缩机。单个空调的功率不大几百瓦到几千瓦但一个城市几百万台空调聚合起来就是几十万千瓦的响应能力。这也是虚拟电厂调度研究中温控负荷始终是热门对象的原因。但是很多研究在做温控负荷调度时会把能效比COPCoefficient 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_t1 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度。聚合组的等效热容C0.18 kWh/℃等效热阻R2.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], breakpointsm.cop_breakpoints[g], valuesm.cop_values[g] )PWL的核心逻辑是把室外温度区间分成若干个断点breakpoints每个断点对应一个COP值求解器通过SOS2约束自动选择相邻两个断点进行线性插值。比如室外温度28度时COP3.830度时COP3.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本文方法采用动态COPPWL建模MILP求解方案B传统定值COPCOP固定为3.5其他约束和求解流程与方案A完全一致方案C不考虑温控负荷把空调当作固定负荷不参与优化调度这样做对比的好处是A vs B能体现动态COP建模对调度经济性的影响A vs C能体现温控负荷参与虚拟电厂调度的整体价值。4.2 优化结果对比先看总成本方案总运行成本元相比C的降幅A动态COP12674.614.2%B定值COP13792.36.7%C无TCL14784.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:0027.53.82520198609:0032.13.31610201914:0037.82.91465135320:0033.53.185401717注意看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、购电功率、售电功率、空调功率、温度变量逐个确认是否声明了nonnegTrue。一个简便的做法是在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_tz_t取0或1。这样虽然增加了二元变量数量但更真实地还原了空调的实际运行状态。结尾的几句实在话做这个项目的过程中我最大的体会是优化调度模型的价值不取决于目标函数写得有多复杂而取决于约束条件对物理世界的还原程度。定值COP的模型不是不能做它也能算出调度计划但那是在不够真实的假设下算出来的执行偏差不可避免。动态能效比感知的本质不是搞一个新奇的数学概念而是让模型在决策时能感知到设备真实的工作状态——温度高的时候空调确实更费电这个事实不该被忽略。如果你打算在自己的项目里用这套方法建议先从方案B定值COP跑通整个流程再切换成动态COP做对比。这样不仅能验证改进效果的显著性也能在调试时更容易定位问题。代码层面任何改动都建议用git管理每个版本的算例结果和参数配置都记录下来做科研数据追溯时真的能救命。后续你可以试着把动态COP扩展到电采暖、热水器等其它温控设备上或者把用户行为的不确定性加入模型这些都值得继续折腾。