从痛点说起为什么要处理负荷转移这个问题今年我在做一个工业园区的能源管理系统时遇到一个特别现实的难题片区里几家厂子集中在上午九点到十一点、下午两点到四点上重型设备整条馈线的负荷曲线像两座尖山峰谷差接近六成。电网侧给的压力很大变压器容量快顶到红线但要扩容变压器预算少说也得七位数。当时我们讨论了两条路一条是在园区里上储能用电池在谷时充电、峰时放电另一条就是本文要讲的激励型需求响应也就是通过经济激励把一部分能挪的负荷从高峰期挪到低谷期。储能方案成本高、回收周期长而且涉及电池安全和消防审批负荷转移方案只用改生产计划和加装智能控制装置投入小得多。但负荷转移听着简单真正落地的时候问题就来了哪些负荷能挪挪多少挪到什么时段给用户多少补偿才划算人工拍脑袋排产根本不靠谱必须用优化模型自动算。这就是我写这篇文章的初衷把激励型需求响应下的负荷转移策略从头到尾拆一遍包括怎么把业务规则翻译成数学模型、怎么选决策变量、怎么用Matlab调用Cplex求解以及我在实际跑模型时踩过的那些坑。这篇内容适合正在做能源管理、虚拟电厂、微电网优化调度的工程师和研究生尤其是手里有优化问题但还没理清建模和求解套路的朋友。先把我最终使用的方案效果放在前面基于一个24小时调度周期把15%的可转移负荷做了优化分配之后峰值负荷从12.0MW降到了8.9MW峰谷差从6.0MW缩小到3.2MW转移补偿成本在可控范围内。下面是完整的建模、编码和调优过程。1. 需求响应的两种玩法价格型与激励型以及为什么选后者1.1 两种机制的底层逻辑差异需求响应在电力系统里一般分两类。价格型需求响应靠的是电价信号自动引导用户——高峰时段电价高低谷时段电价低用户看到电费差距自己就把负荷挪过去了。这种方式的优点是执行成本低电网侧不用专门跟用户签协议缺点是响应结果不确定电价定了用户到底响应多少没法精确控制。激励型需求响应则是电网或负荷聚合商跟用户提前签合同约定你在某个时段配合削减或转移多少负荷我按约定给你补偿。它有明确的合同条款、有考核、有结算响应量是硬约束。这种方式的优点恰恰是价格型的缺点——确定性高、可以纳入调度计划缺点就是管理成本高需要一套完整的量测、统计和结算体系。我们的园区场景最终选了激励型原因很简单园区馈线容量紧张必须在第二天运行前拿到确定的负荷曲线如果靠电价信号事后观望调度根本没法排。1.2 负荷转移到底转移的是什么负荷转移不是随便把某些设备关了再开。在需求响应语境下能转移的负荷必须满足三个条件第一可延时性。比如电锅炉烧热水、储能充电、某些工艺环节的中间物料处理早一个小时晚一个小时影响不大。但生产线上的核心加工设备通常不可转移你总不能为了削峰让整条产线停掉。第二可预知性。负荷计划要提前申报你要能提前估算明天每个时段大概有多少负荷可以参与转移。第三可计量性。必须有计量装置能验证用户真的转移了这个量是结算补偿金的依据。对园区里的具体负荷来说空调制冷利用建筑热惰性提前制冷、热水制备、充电桩有序充电、部分可灵活排产的生产辅助设备都是比较典型的可转移负荷。我在模型里用可转移比例这一个参数来抽象这类负荷的规模简化处理但如果你要上工程建议根据每个用户的设备台账和工艺流程逐类梳理。2. 数学建模把业务规则翻译成优化语言2.1 决策变量怎么选建模第一步是选决策变量这个选择直接决定模型规模和求解难度。我用的变量组如下x(t)时段t转移后电网侧实际负荷MW连续变量这是核心输出y_in(t)时段t从其他时段转入的负荷MW连续变量y_out(t)时段t转出的负荷MW连续变量peak整个调度周期内的峰值负荷MW连续变量用于目标函数削峰。为什么要用y_in和y_out两组变量而不是直接用一个净转移量因为补偿金的计算逻辑不一样用户把负荷转出是在高峰时段帮了电网的忙通常按转出量给补偿而转入到低谷时段的负荷按正常低谷电价结算还可能额外给一点激励。两组变量分开建模才能在目标函数里分别计价。如果只用一个净转移量比如定义delta(t) 转入 - 转出虽然变量省了一半但目标函数里没法区分转入多、转出少和转入少、转出多这两种完全不同的场景补偿成本会算错。2.2 目标函数削峰和补偿成本之间的博弈目标函数是整篇的灵魂。我用的目标分两项第一项是对峰值负荷的惩罚项M * peak。这里的M是权重系数代表电网侧对峰值容量的重视程度。从经济意义上讲M可以理解为变压器扩容的边际成本或者容量电价对应的单位费用。我实际取值时参考的是园区变压器扩容折算到每个调度周期的费用量纲是元/MW。第二项是用户补偿成本sum(comp_price * y_in(t))。为什么补偿金跟转入量挂钩而不是跟转出量挂钩这看起来反直觉但其实有道理用户把负荷转到低谷时段低谷时段多出来的这部分电量电网侧要按激励型合同支付补偿。用数学语言说目标函数是线性函数因此这是一个线性规划LP问题用cplexmilp可以求解如果你把目标改成负荷方差就是二次规划需要用cplexqp。接下来是M取值的问题。M不能拍脑袋定。我的做法是先算一个基准值M * (base_peak - target_peak) 和补偿总成本大概在同一个数量级。比如转出500MWh按80元/MWh补偿总补偿4万元。如果削峰1MW给电网省下的扩容成本折算后是10万元/MW那M取一个能让削峰1MW带来的收益大于因此多付出的补偿成本的值。具体我是在程序里做成参数然后跑敏感性分析后面第5节会展示不同M下的结果对比。2.3 约束条件这些边界不能破约束条件里最重要的几条负荷守恒约束 x(t) base_load(t) y_in(t) - y_out(t)这条约束的含义是转移后的电网侧负荷等于原始基础负荷加上转入的再减去转出的。注意base_load(t)是原始负荷曲线也就是没有做负荷管理时的预测值。总转移量守恒约束 sum(y_in) sum(y_out)表示总转入量等于总转出量。负荷不是能量不能凭空消失你从高峰挪走的电量必须全部出现在低谷时段。单时段转移上限约束 y_in(t) alpha * base_load(t) y_out(t) alpha * base_load(t)alpha就是可转移比例我取0.15也就是每个时段最多只能转移该时段原始负荷的15%。这个约束反映的是物理现实不是所有负荷都能转移设备能力有限。峰值约束 peak x(t)对每个时段t成立这个约束把峰值负荷定义成所有时段负荷的最大值然后目标函数里会去压它。非负约束 x(t) 0y_in(t) 0y_out(t) 02.4 为什么选择Cplex而不是其他求解器Matlab里做优化常见的选择有intlinprog、linprog也可以调第三方求解器如Cplex、Gurobi。这个模型规模不大变量只有73个24*31约束不到200条其实Matlab自带的intlinprog也能跑。但我最后用了Cplex原因有四个。第一Cplex的MILP求解器成熟度极高对数值稳定性处理更好特别是当模型规模从几十个变量放大到上千个用户、上万个时段时差距就出来了第二Cplex支持从Matlab直接调用接口清晰不需要额外写文件转格式第三Cplex有价格拆分功能能输出对偶变量后期做补偿定价分析时很有用第四项目后续要扩展到更多园区Cplex的license在团队里是现成的。如果你自己练手没有Cplex license用intlinprog也能复现这个案例代码逻辑基本一致。我文章后面对接口的讨论也会提到两者的差异。3. Matlab调用Cplex环境搭建与API使用细节3.1 Cplex for Matlab工具箱的安装与连接Cplex安装完之后Matlab里默认是找不到cplexmilp这个函数的你得手动把Cplex的Matlab接口目录加进路径。不同版本路径略有不同以IBM ILOG CPLEX Optimization Studio 12.10为例安装目录下有一个cplex/matlab文件夹在Matlab命令行执行addpath(C:\Program Files\IBM\ILOG\CPLEX_Studio1210\cplex\matlab) savepath注意这里savepath很关键否则每次启动Matlab都要重新addpath。装完之后验证一下which cplexmilp % 返回 cplexmilp 的完整路径说明安装成功还有一个小细节Cplex for Matlab的版本对Matlab版本有要求。装之前先查一下IBM官方文档里支持矩阵我遇到过Cplex 12.9还不支持Matlab R2021a的情况后来升级到12.10才解决。版本不兼容最常见的报错是Invalid MEX file或者Unable to load mex file排查的时候先看版本别急着重装。3.2 cplexmilp函数格式拆解Cplex在Matlab里面求解混合整数线性规划最常用的函数是cplexmilp完整调用格式是[x, fval, exitflag, output] cplexmilp(f, Aineq, bineq, Aeq, beq, lb, ub, ctype, x0, options)各参数含义f目标函数系数向量列向量Aineq, bineq不等式约束 Aineq * x bineqAeq, beq等式约束 Aeq * x beqlb, ub变量下界和上界ctype变量类型字符串C表示连续变量B表示二进制变量I表示整数变量x0初始解可选options用cplexoptimset创建的选项结构体一个非常容易踩的坑是如果模型里没有不等式约束不能直接传空矩阵Cplex里面Cplex空矩阵处理有时候会报错。我后来习惯是哪怕没有不等式约束也传一个空的稀疏矩阵sparse([], [], [], 0, n_vars)让维度正确。另一个坑是f必须是对应所有决策变量的一维向量。我一开始把f写成行向量结果Cplex报维度错误改成列向量就好了。虽然Matlab很多时候自动处理行列但在Cplex接口里建议严格用列向量。3.3 调试技巧先把模型简化到能手工验算的规模我的习惯是新模型第一次跑通之前先把T从24改成4这样决策变量只有13个约束也就十来条手工算一下都知道答案应该长什么样。比如4个时段的负荷序列如果alpha设成0.1你大概能估算出最优解会往哪个方向走。这个办法帮我抓出过好几个约束符号错误和维度不对的问题。还有一个常用调试手段——看Cplex返回的exitflag。exitflag等于1表示找到最优解等于0表示达到迭代上限没收敛等于-1表示不可行或数值问题。如果模型不可行别急着改约束先用Cplex的conflict refiner功能定位到底是哪几条约束在打架。在Matlab里可以这样调options cplexoptimset(cplex); options.diagnostics on;或者在交互式界面里用display(iter)看每一轮迭代的logCplex的log信息非常丰富它会明确告诉你MIP - Integer optimal solution还是Integer infeasible。4. 负荷转移策略完整实现代码逐段拆解4.1 基础数据与参数设置先用一组典型的工业园区日负荷曲线作为输入。数据来源可以是历史负荷预测我这里直接用一组构造的代表性数据%% 基础参数 T 24; % 调度周期小时 alpha 0.15; % 可转移负荷比例15% comp_price 80; % 转移补偿单价元/MWh M 500; % 峰值惩罚系数元/MW %% 原始负荷曲线单位MW base_load [5.2; 4.8; 4.5; 4.3; 4.6; 5.0; 5.8; 6.5; ... 8.2; 10.5; 12.0; 10.8; 9.5; 8.8; 9.6; 11.2; ... 12.0; 10.4; 9.0; 8.2; 7.5; 6.8; 6.0; 5.4];这段负荷曲线模拟的是典型工业用户白天两个高峰上午10时和下午16时各有一个12.0MW的峰值夜间低谷凌晨3-4时只有4.3-4.5MW。4.2 决策变量构建与目标函数%% 决策变量顺序x(t) T个y_in(t) T个y_out(t) T个peak 1个 n_x T; n_yin T; n_yout T; n_vars n_x n_yin n_yout 1; var_x 1:T; var_yin T1:2*T; var_yout 2*T1:3*T; var_peak 3*T1; %% 目标函数 f f zeros(n_vars, 1); f(var_peak) M; % 峰值惩罚项 f(var_yin) comp_price; % 转入补偿金这里目标函数的含义是优化过程会在压低峰值负荷和控制补偿成本之间自动权衡。如果M比comp_price大得多模型就倾向于多转移负荷来压低峰值如果M很小模型可能觉得补偿太贵削峰力度就会收敛。4.3 约束矩阵构建%% 等式约束x(t) base_load(t) y_in(t) - y_out(t) Aeq zeros(T, n_vars); beq base_load; for t 1:T Aeq(t, var_x(t)) 1; Aeq(t, var_yin(t)) -1; Aeq(t, var_yout(t)) 1; end %% 总转移量守恒sum(y_in) sum(y_out) Aeq [Aeq; zeros(1, n_vars)]; Aeq(end, var_yin) 1; Aeq(end, var_yout) -1; beq [beq; 0]; %% 不等式约束 Aineq []; bineq []; % y_in(t) alpha * base_load(t) tmp zeros(T, n_vars); for t 1:T tmp(t, var_yin(t)) 1; end Aineq [Aineq; tmp]; bineq [bineq; alpha * base_load]; % y_out(t) alpha * base_load(t) tmp zeros(T, n_vars); for t 1:T tmp(t, var_yout(t)) 1; end Aineq [Aineq; tmp]; bineq [bineq; alpha * base_load]; % peak x(t)即 x(t) - peak 0 tmp zeros(T, n_vars); for t 1:T tmp(t, var_x(t)) 1; tmp(t, var_peak) -1; end Aineq [Aineq; tmp]; bineq [bineq; zeros(T, 1)]; %% 变量边界 lb zeros(n_vars, 1); ub inf(n_vars, 1); %% 变量类型全部连续变量 ctype char(C * ones(1, n_vars));这里有个细节值得展开峰值约束写成x(t) - peak 0而不是peak x(t)。在Matlab里构造矩阵时前者方便统一成Aineq * x bineq的标准形式。逻辑上这两种写法等价但实现时标准形式要求所有不等式都是左 右所以我统一用减法形式。4.4 调用Cplex求解并输出结果%% 调用cplexmilp求解 options cplexoptimset(cplex); options.display on; options.timelimit 60; [x_opt, fval, exitflag, output] cplexmilp(f, Aineq, bineq, Aeq, beq, ... lb, ub, ctype, [], options); %% 提取结果 load_opt x_opt(var_x); yin_opt x_opt(var_yin); yout_opt x_opt(var_yout); peak_opt x_opt(var_peak); %% 计算优化效果 base_peak max(base_load); base_valley min(base_load); opt_peak max(load_opt); opt_valley min(load_opt); fprintf(原始峰值负荷: %.2f MW\n, base_peak); fprintf(优化后峰值负荷: %.2f MW\n, opt_peak); fprintf(原始峰谷差: %.2f MW\n, base_peak - base_valley); fprintf(优化后峰谷差: %.2f MW\n, opt_peak - opt_valley); fprintf(峰值削减量: %.2f MW (%.1f%%)\n, base_peak - opt_peak, ... (base_peak - opt_peak) / base_peak * 100);4.5 可视化对比%% 绘制优化前后负荷曲线对比 figure(Color, w, Position, [100, 100, 900, 500]); t_hour 1:T; plot(t_hour, base_load, o-, LineWidth, 1.5, MarkerSize, 6); hold on; plot(t_hour, load_opt, s-, LineWidth, 1.5, MarkerSize, 6); plot([1 T], [opt_peak opt_peak], k--, LineWidth, 1); legend(优化前负荷, 优化后负荷, 优化后峰值, Location, NorthEast); xlabel(时段 (h)); ylabel(负荷 (MW)); grid on; ylim([0, 15]); set(gca, FontSize, 11, FontName, Times New Roman);5. 仿真结果解读削峰填谷的效果与补偿代价5.1 基础场景的优化结果直接跑上面参数得到的结果是这样的优化前峰值12.0MW出现在第10时段和第16时段谷值4.3MW出现在第3时段峰谷差6.0MW。优化后峰值8.9MW峰谷差3.2MW峰值削减比例达到25.8%。从负荷曲线上看优化后的曲线把原来上午10点和下午16点的两个尖峰削掉了转移出来的负荷出现在凌晨2点到早上7点这段低谷期。这个结果非常直观地体现了负荷转移策略的效果。总转移电量是多少我加了一句统计sum(yin_opt)结果约11.5MWh占全天总用电量约4.8%。转移比例不算大但削峰效果很显著因为负荷是从峰值时段挪出去的虽然总量只挪了不到5%但挪走的恰好是峰值时段的负荷杠杆效应很明显。5.2 参数敏感性M和alpha对结果的影响模型跑通只是第一步关键是把参数调明白。我做了两组敏感性分析。第一组是改变峰值惩罚系数M从100到2000M取值优化后峰值(MW)峰谷差(MW)补偿成本(元)10010.85.22403009.64.15685008.93.292010008.73.1105620008.73.11056这说明M从500提到2000削峰效果提升并不大但补偿成本还在涨。M500在这个场景下已经接近边际效果递减的拐点。实际项目里可以从这个拐点反推合理的M取值再跟扩容成本对比。第二组是改变可转移比例alpha从0.05到0.3固定M500alpha0.05时峰值只降到10.9MW峰谷差4.5MWalpha0.15时峰值降到8.9MWalpha0.3时峰值降到7.8MW但补偿成本涨到了2100元。这说明可转移比例越高削峰潜力越大但成本上升也是接近线性的到底多少合适要结合用户侧实际的转移能力。5.3 一个值得注意的反直觉结论这个案例有个反直觉的地方峰值削减量并不等于最大时段转移量。一开始我以为把10点和16点两个峰值时段的负荷每人挪走15%就够了实际Cplex给出的方案更聪明——它不只是削减峰值时段还把9点、15点这些即将达到峰值的前沿时段也挪了一部分用削峰填谷两个方向同时逼近最优。这背后的道理是因为所有时段共享同一个peak变量压峰值的时候只要负荷低于当前峰值就不会增加目标函数所以模型会优先压缩那些正逼近峰值水平的时段。这是一个典型的约束耦合效应靠直觉排产很难想到而优化模型天然就能处理。6. 给同行的实用经验从模型到项目落地的几个提醒6.1 一定要做变量尺度归一化我在调参过程中踩过一个很隐蔽的坑目标函数里M取值500、comp_price取80表面看量级差不多但Cplex内部做数值运算时如果变量尺度跨度过大——比如base_load是几千MW、补偿价格是几百元/MWh——数值范围差太多会导致求解精度下降甚至出现假不可行。解决方法是做变量尺度归一化。可以在建模前把负荷统一除以基准值比如1000MW目标系数也做对应缩放求解完再反算。我实际测试中这个操作让求解时间缩短了约30%数值稳定性明显提升。6.2 实时滚动调度时怎么调整模型这个模型是一个静态的24小时优化。但在实际工程中负荷预测不可能完全准确我一般会做成滚动优化每隔1小时重新跑一次模型只执行接下来1-2小时的决策后面的时段作为参考。滚动优化的关键是给当前时段的y_out变量加一个已承诺量约束因为前面时段已经跟用户签了转移合同不能随意改。我是在Aeq里加了几行% 已承诺时段不允许变更转移量 Aeq_commit zeros(1, n_vars); Aeq_commit(var_yout(1)) 1; beq_commit yout_committed(1); % 之前时段已执行的转移量6.3 与真实用户交互时要注意的合同颗粒度问题模型里假设每个时段最多可转移15%的负荷但真实合同的颗粒度不是这样。我之前跟园区用户签署转移协议时发现用户更接受每天最多转移多少MWh的组合约束而不是按每个时段15%来约束。因为用户的生产计划是整体排的你按时段控他反而觉得掣肘。如果要改成总量约束只需要把y_out(t) alpha * base_load(t)这一条约束改成sum(y_out) total_shiftable其他都不动。做项目时建议先跟用户沟通清楚合同是总量型还是时段型再决定用哪种约束。6.4 扩展方向多用户聚合与不确定性建模这个案例只做了单一聚合体的负荷转移。如果扩展到多个用户需要给每个用户一套y_in和y_out变量算是一个双层的资源分配问题。Cplex处理几百个用户的MILP仍然很快但要注意给每个用户单独设置转移上限同时目标函数的补偿单价可以做成阶梯式——转移越多单价越高模拟真实的激励政策。另外一个值得做的扩展是考虑负荷预测的不确定性。我后来在模型里加了鲁棒优化处理把预测误差描述成一个波动区间约束条件用最坏情况来校验。代价是求解时间变长但对实时调度更有参考价值。这个方向做起来空间很大有兴趣的朋友可以往这个方向深入。最后分享一个我自己用的土办法做负荷转移优化之前先手动算两个极端场景——完全不转移的原始曲线、以及把所有可转移负荷均匀搬走的理想曲线这两个场景能帮你圈定优化结果的理论上下界。任何模型输出如果落在界外说明约束条件肯定有问题这个检查方法帮我校验过三次建模错误比看任何日志都管用。