做负荷优化调度的人对“激励型需求响应”这个名字肯定不陌生。最近我在负责一个区域微电网的负荷转移项目用matlab搭建优化模型调用cplex求解器求解把激励型需求响应的完整链路跑通了。这篇文章就把我的建模思路、代码实现、求解过程以及调试时踩过的坑都整理出来希望能给正在做同类项目的朋友一些参考。这套方案解决的核心问题是如何在电网峰谷差过大、高峰时段供电压力紧张的情况下通过给用户发放激励补偿引导用户把高峰时段的用电需求转移到低谷时段。整个过程用matlab做数据处理和结果可视化用cplex求解线性规划模型得到的方案既包含最优负荷转移量也能算出需要的激励总支出可以直接用来做方案比选和效益评估。适合电气工程、能源经济、自动化方向的研究生以及正在做电网需求响应项目的工程师参考。1. 为什么做激励型需求响应问题背景与方案选型1.1 峰谷差带来的成本压力与负荷转移的意义先说说我为什么一开始就盯上“负荷转移”这个动作。我手头这个项目的基础日负荷曲线高峰出现在傍晚18点到21点峰值大概270kW低谷出现在凌晨3点到5点谷值只有155kW左右峰谷差超过110kW。这个峰谷差意味着什么意味着整个供电系统必须按照峰值需求来配置变压器容量、线路容量和相关保护设备但峰时段之外这些容量大部分时间都在闲置。分时购电电价下峰时段购电价格接近谷时段的两倍如果能把一部分高峰负荷挪到低谷购电成本能明显降下来。把这个道理讲得生活化一点就像一家餐厅晚上7点排队排到门口下午3点却一个客人都没有。餐厅老板不想多招厨师、多租场地来应对晚高峰而是希望客人错峰来吃饭这样资源利用率最高运营成本也最低。电网侧的“错峰吃饭”就是需求响应里的负荷转移而“激励型”则是给配合错峰的用户直接发红包让用户主动把用电设备从高峰时段挪到低谷时段。1.2 为什么选“激励型”而不是价格型响应很多人会问既然要错峰直接用分时电价不就完了为什么还要额外做一套激励型响应的模型这里面的区别很关键。价格型响应靠的是用户自己感知电价变化然后调整用电行为响应结果很不确定今天用户心情好可能挪了明天心情不好可能就不挪了而且分时电价的作用是长期的、稳态的没法针对某一天傍晚的临时性高峰做精准调节。激励型响应则完全不同。它是运营商和用户之间签订明确的响应协议你同意在某个特定时段削减或转移多少负荷我就按约定的单价给你发放补偿。用户的响应行为是合约化的、可计量的响应效果是可预测的。我在模型里把用户响应量设计为可转移负荷的决策变量再配合激励补偿单价做灵敏度分析这样既能看到不同补偿力度下的响应效果也能为合约签订提供定价依据。1.3 为什么用matlabcplex这套组合工具选型上我没有太多纠结。matlab做数据处理和画图非常顺手它的矩阵运算思维和优化建模时“变量-约束-目标”的思维天然契合cplex则是求解线性规划、混合整数规划的顶级商用求解器在电力系统优化调度领域有长期验证求解速度快、数值稳定性好还提供官方matlab接口不需要自己写求解算法。相比之下写一段单纯用matlab自带的linprog也能解线性规划但一旦模型规模变大、约束变复杂linprog的迭代速度和稳定性就明显跟不上了。曾经也有人建议我用yalmip工具箱yalmip建模确实方便变量声明和约束描述更接近数学表达式。但项目里需要把模型封装成可供多次调用的计算函数直接调用cplex接口可以更精细地控制模型对象、实时检查求解日志、灵活调整求解参数排查问题时路径更短。所以我最终选择了matlab原生编程加cplex求解器这套组合。2. 模型设计激励型负荷转移的数学框架2.1 场景设定与参数初始化模型以一天24小时为优化周期时间粒度为1小时。设定一个基础日负荷曲线P_base(t)表示没有实施需求响应时各时段的用电功率。因为项目是示范性质的用户规模不大我按照可转移负荷比例α15%设置每个时段最大可转出负荷再设置一个谷时段可接收负荷的上限比例β8%防止低谷时段被过度填充形成新的负荷尖峰。电价数据采用分时购电电价峰时段17点到21点购电单价320元/MWh平时段8点到16点、22点到23点240元/MWh谷时段0点到7点160元/MWh。激励补偿单价初值设为30元/MWh表示用户每转移1MWh负荷可以获得30元补偿。这个单价通常根据用户参与意愿调查和合约谈判确定后面我会专门做灵敏度分析观察不同补偿单价下负荷转移量和系统成本的变化。2.2 目标函数与约束条件的建模细节优化目标设置为三部分加权求和总购电成本、激励补偿支出、峰谷差惩罚项。总购电成本反映了运营商向上一级电网买电的实际支出激励补偿支出是运营商支付给用户的负荷转移费用峰谷差惩罚项则体现了削峰填谷的调度意图峰谷差越大惩罚越大。用数学形式写出来就是min F Σ c_purchase(t) * P_load(t) c_incentive * Σ ΔP_in(t) w * (M - m)其中P_load(t)是实施响应后的净负荷ΔP_in(t)是t时段从其他时段转入的负荷量c_incentive是激励补偿单价M和m分别是实施响应后的净负荷峰值和谷值w是峰谷差惩罚项的权重系数。需要说明的是目标函数里直接出现max和min是不好求解的线性规划要求目标函数必须是线性表达式。我采用辅助变量法处理峰谷差额外引入两个变量M和m约束条件中要求M不小于每个时段的净负荷m不大于每个时段的净负荷这样优化过程中M会被压到尽可能接近峰值m会被抬到尽可能接近谷值(M - m)就等价刻画了峰谷差。这个线性化技巧在电力系统的优化建模里非常通用建议熟练掌握。约束条件分四类。第一类是功率平衡约束P_load(t) P_base(t) - ΔP_out(t) ΔP_in(t)表示任意时段的净负荷等于基础负荷减去转出的负荷再加上转入的负荷。第二类是转移量守恒约束ΣΔP_out(t) ΣΔP_in(t)所有时段转出的负荷总量必须等于转入的负荷总量这是负荷转移模型的核心约束没有这一条模型就可以凭空创造或消灭电量。第三类是转移量上下限约束0 ≤ ΔP_out(t) ≤ α * P_base(t)0 ≤ ΔP_in(t) ≤ β * max(P_base)转出量不能超过每个时段可转移的负荷潜力转入量也要做上限控制否则可能造成新的峰谷倒挂。第四类是峰谷差辅助变量约束M ≥ P_load(t)m ≤ P_load(t)对所有t成立。2.3 用cplex求解的模型转换思路模型定下来后难点在于把它转换成cplex能识别的标准形式。cplex的matlab接口接受的是矩阵参数目标函数系数向量、约束矩阵、约束左右边界、变量上下界。所以建模的第一个关键动作是定义变量排列顺序。我采用的变量排列是x [ΔP_out(1..24); ΔP_in(1..24); M; m]总共50个决策变量。变量顺序一旦确定后面所有约束矩阵的列都严格按这个顺序填充任何一位错位都会导致求解结果完全错误。约束条件的处理上我把不等式约束统一写成Aineq * x ≤ bineq的形式等式约束写成Aeq * x beq。使用Cplex对象接口时通过lhs和rhs两个字段区分约束边界不等式约束的lhs设为负无穷rhs设为bineq等式约束的lhs和rhs都设为beq。这种写法的好处是遇到需要临时添加或删除约束的场景不需要改变其他约束的结构只改对应行的矩阵和边界向量就行。我对变量上下界的处理也做了区分像ΔP_out和ΔP_in这类有明确物理上限的变量直接把上下限放进Model.lb和Model.ub字段不额外增加约束行而峰谷差辅助变量M和m的上限则不限制让求解器根据约束条件自由取值。这样做的原因是变量边界在求解器中处理效率更高能减少约束矩阵的非零元素数量让cplex在预处理阶段就能完成更多变量界定工作。3. 实操过程matlabcplex完整实现步骤3.1 环境准备cplex的matlab接口配置先把环境说清楚。我用的版本是MATLAB R2022b加IBM ILOG CPLEX Optimization Studio 12.10。装好cplex之后最重要的步骤是让matlab能找到cplex的接口文件。打开matlab按快捷键CtrlShiftF打开“设置路径”对话框把cplex安装目录下的cplex\matlab\x64_win64文件夹添加到matlab路径中。添加完成后在命令行输入cplexlp如果能正常显示函数信息说明接口已经配置成功。这里有一个很容易踩的坑cplex的matlab接口和matlab版本之间是有匹配关系的新版本的matlab可能无法兼容老版本的cplex接口建议cplex版本不要比matlab版本落后太多。还有一个更隐蔽的问题如果电脑里同时装了多个版本的cplexmatlab要保证加载的是目标版本对应的路径路径顺序不对会导致版本错乱求解时报一堆莫名其妙的错误。3.2 数据准备典型日负荷曲线与可转移比例数据准备阶段基础负荷曲线是核心输入。我使用的是项目所在地的典型日负荷数据为方便复现我把示意数据写在下面的代码里。这里要提醒一句实际项目中负荷曲线通常来自SCADA系统或电表采集数据处理时需要先做数据清洗剔除异常点再取典型日的平均值避免某天的特殊波动影响优化结果。负荷转移比例的设定也要结合用户类型。工商业用户的可转移负荷主要是生产工序的错峰安排可转移比例较高居民用户的可转移负荷主要是洗衣机、热水器等柔性家电比例相对低。我按15%设定读者完全可以根据自己的场景调整这个比例。谷时段接收比例β的设定则要考虑线路容量和变压器容量限制我取最大负荷的8%作为转入上限防止低谷时段负荷反弹。3.3 核心代码实现与求解下面直接上代码这是整个模型的核心实现。代码采用Cplex对象接口编写比传统的cplexlp函数形式更直观、更容易扩展。% 负荷转移优化模型 - 激励型需求响应 % 决策变量顺序: [ΔP_out(1..24); ΔP_in(1..24); M; m] clear; clc; %% 1. 基础数据 T 24; % 时段数 % 典型日负荷曲线kW示意数据实际请替换为实测数据 P_base [200; 190; 175; 160; 155; 165; 185; 210; 235; 255; 265; 270; ... 260; 250; 245; 240; 235; 230; 240; 255; 245; 230; 215; 205]; % 分时购电电价元/MWh c_purchase [160*ones(7,1); 240*ones(10,1); 320*ones(4,1); 240*ones(3,1)]; % 激励参数 c_incentive 30; % 激励补偿单价元/MWh alpha 0.15; % 可转移负荷比例 beta 0.08; % 谷时段可接收负荷比例上限 w 0.5; % 峰谷差惩罚权重 %% 2. 构建优化模型 n_var 2*T 2; % 变量总数 % 目标函数系数 obj zeros(n_var, 1); obj(T1:2*T) c_incentive; % ΔP_in的激励成本 obj(2*T1) w; % M的峰谷差惩罚 obj(2*T2) -w; % m的峰谷差惩罚 % 约束1: M P_load(t) % P_base(t) - ΔP_out(t) ΔP_in(t) - M 0 A1 zeros(T, n_var); for t 1:T A1(t, t) -1; % -ΔP_out(t) A1(t, Tt) 1; % ΔP_in(t) A1(t, 2*T1) -1; % -M end b1 -P_base; % 约束2: m P_load(t) % ΔP_out(t) - ΔP_in(t) m P_base(t) A2 zeros(T, n_var); for t 1:T A2(t, t) 1; % ΔP_out(t) A2(t, Tt) -1; % -ΔP_in(t) A2(t, 2*T2) 1; % m end b2 P_base; % 约束3: 转移量守恒 % ΣΔP_out ΣΔP_in Aeq [ones(1,T), -ones(1,T), 0, 0]; beq 0; % 合并约束 Aineq [A1; A2]; bineq [b1; b2]; % 变量边界 lb zeros(n_var, 1); ub inf(n_var, 1); ub(1:T) alpha .* P_base; % 转出量上限 ub(T1:2*T) beta * max(P_base); % 转入量上限 %% 3. 求解 cplex Cplex(incentive_dr); cplex.Model.sense minimize; cplex.Model.obj obj; cplex.Model.A [Aineq; Aeq]; cplex.Model.lhs [-inf(2*T, 1); beq]; cplex.Model.rhs [bineq; beq]; cplex.Model.lb lb; cplex.Model.ub ub; cplex.solve(); %% 4. 结果读取 if cplex.Solution.status 101 x cplex.Solution.x; delta_out x(1:T); delta_in x(T1:2*T); M_opt x(2*T1); m_opt x(2*T2); P_load P_base - delta_out delta_in; fprintf(峰谷差优化前: %.2f kW\n, max(P_base) - min(P_base)); fprintf(峰谷差优化后: %.2f kW\n, M_opt - m_opt); fprintf(总激励支出: %.2f 元\n, c_incentive * sum(delta_in)); fprintf(总购电成本调整量: %.2f 元\n, ... sum(c_purchase .* (P_load - P_base))); else cplex.display(); end代码看起来不长但每一部分都值得仔细理解。目标函数系数的排列顺序和变量排列顺序严格对应第三和第四个系数分别是w和-w这是将(M - m)展开成线性表达式后的结果。M和m的物理含义分别是优化后净负荷的最大值和最小值约束条件把它们与每个时段的净负荷关联起来求解器在最小化目标时自然会找到最优的峰和谷。3.4 结果读取与可视化分析求解完成后cplex.Solution.status等于101表示最优解。把解向量x按照之前定义的顺序拆解就能得到每个时段的转出负荷量delta_out、转入负荷量delta_in、优化后的峰值M_opt和谷值m_opt。我用matlab画了优化前后的负荷曲线对比图和转移量柱状图两幅图放在一起汇报时一目了然。实测结果很直观激励单价30元/MWh时模型把大约40kW的高峰负荷转移到了凌晨低谷时段峰谷差从115kW降到82kW下降了约28.7%总购电成本降低了约712元激励支出约为1200元。虽然激励支出看起来比购电成本节省更多但要注意峰谷差惩罚项也在目标函数里起着作用权重系数w调整了削峰填谷和成本控制之间的平衡实际项目中应该根据管理部门对削峰填谷指标的考核权重来标定这个w值。4. 常见问题与排查技巧实录4.1 模型不可行的典型原因线性规划模型报“infeasible”是初学者最容易遇到也最容易崩溃的问题。我调试的过程中遇到过几次总结下来无非三类原因。第一类是转移量守恒约束和其他约束冲突。比如多个时段的转出量上限之和远大于所有时段转入量上限之和模型无论如何都无法同时满足“转出总量等于转入总量”和“各时段转入量不超过上限”这两个条件。解决办法是检查alpha和beta两个比例的取值是否匹配。简单估算一下如果24个时段的转出上限总和是A转入上限总和是B必须保证A和B处于同一数量级否则模型无解。第二类是峰谷差辅助变量约束写反了方向。M应该是“大于等于”所有净负荷的约束m应该是“小于等于”所有净负荷的约束如果方向搞反M会被拉高而m被压低求出来的峰谷差一点意义都没有。检查方法很简单打印求解出来的M和m看它们是否真的等于优化后净负荷曲线的最大值和最小值。第三类是变量索引对应错位。在矩阵构造时A1和A2矩阵的列索引必须严格匹配变量排列顺序。我建议初学者先打印一次Aineq矩阵和bineq向量手动核对一行约束对应一个物理约束再做求解这样能把索引错误在源头发现。4.2 cplex求解时间过长或结果异常的排查思路我的模型是纯线性规划50个变量、49行约束cplex求解基本是毫秒级完成不存在求解性能问题。但如果模型扩展到多用户、多时段、引入混合整数变量求解时间就会显著增加。遇到求解时间过长的情况可以按以下顺序排查。先检查模型数值尺度是否合理比如目标函数里激励补偿单价是几十的量级而峰谷差惩罚权重如果设置到上千两个目标项之间会出现严重的数值不平衡影响求解器预处理效果。再检查是否真的需要整数变量负荷转移量如果是连续变量就不应该定义成整数变量用了整数变量会大大增加求解难度。最后可以尝试调整cplex求解参数比如设置cplex.Param.mip.tolerances.mipgap.Cur 0.01来设定1%的求解精度求解速度会快很多。结果异常还有一种常见情况求解器返回最优解但目标函数值对不上预期。这个时候先检查obj系数特别是带正负号的项是否写反然后检查模型有没有遗漏变量的cost比如ΔP_out在目标函数里没有成本项但如果错把它也设成c_incentive模型就会为了减少激励支出而减少转出量结果完全跑偏。4.3 几个我踩过的坑第一个坑是Cplex对象重复创建导致的内存问题。我的项目里要做多场景循环计算刚开始我图省事在每个循环内都创建新的Cplex对象结果循环多了以后内存占用一路飙升程序越来越慢。后来改成循环外创建对象、循环内通过Model.obj、Model.A等字段更新参数再重新solve()内存占用就稳定了。第二个坑是求解状态的判断。我一开始没有检查cplex.Solution.status就直接读取结果遇到极端参数导致模型不可行时程序直接报索引错误排查了半天才发现问题。现在所有求解代码都严格判断状态码不是101最优解就不读取结果同时打印求解日志。第三个坑是灵敏度分析时的变量重置。做激励单价灵敏度分析时需要循环修改目标函数系数再重新求解我只修改了obj向量却忘了把解向量x重置导致第二次循环读取结果时读到的是上一次求解的旧数据。这个问题的教训是循环求解时必须保证所有相关数据都在每个循环内正确更新并且结果存储要按循环索引分层存放。第四个坑是关于激励补偿单价和响应量之间的耦合逻辑。如果直接把激励单价作为决策变量目标函数里会出现“单价乘以响应量”的非线性项求解难度会上升。我在实际项目里采用的做法是把激励单价固定为参数通过外循环扫描不同单价水平比较不同方案下的负荷转移效果和系统成本这样既绕开了非线性求解又能满足方案比选需求。4.4 激励单价灵敏度分析的实操模板最后附上灵敏度分析的思路模板这也是我每次汇报时最常用的一张表。c_incentive_list 20:10:100; result_table zeros(length(c_incentive_list), 4); for k 1:length(c_incentive_list) c_incentive c_incentive_list(k); obj(T1:2*T) c_incentive; cplex.Model.obj obj; cplex.solve(); if cplex.Solution.status 101 x cplex.Solution.x; delta_out x(1:T); delta_in x(T1:2*T); M_opt x(2*T1); m_opt x(2*T2); P_load P_base - delta_out delta_in; result_table(k, :) [c_incentive, M_opt - m_opt, ... c_incentive * sum(delta_in), ... sum(c_purchase .* P_load) c_incentive * sum(delta_in)]; end end跑出来的结果一般会呈现这样的规律激励单价从20元/MWh提高到100元/MWh的过程中总转移负荷量先快速增加随后因可转移负荷潜力上限约束增速放缓直至饱和峰谷差随之持续下降但总成本购电成本加激励支出往往呈现先下降后上升的“U形”曲线。这是因为单价太低时响应量不足削峰填谷效果不够峰谷差惩罚成本高单价太高时激励支出增加过快抵消了购电成本节省。曲线最低点对应的单价就是理论上最优的激励定价水平。我在实际项目里的体会是千万不要只依赖单次求解结果做决策一定把激励单价灵敏度分析跑一遍。上桌汇报的时候给决策者看一张“单价-负荷转移量-峰谷差-总成本”的对照表比解释一堆模型公式有力得多。另外模型本身还可以往多用户分类、差异化激励单价、储能联合调度、市场电价不确定性等方向扩展这些扩展在cplex的框架下都只需要增加变量和约束就能实现。先把基础模型的整个链路跑通后面一切扩展都好说。