我先交代一下背景。这个项目是我在实际课题里被问到最多的一类问题——多主体综合能源系统、需求响应、电能交互、主从博弈四个关键词堆在一起看着像四座山但真正落地成Matlab代码时核心就一句话谁先出招谁后应对怎么收敛到均衡点。这篇博文我尽量不把通用建模教科书复述一遍而是把整套思路、建模细节、Matlab实现的关键节点以及我实际跑代码时踩过的坑都摊开来讲清楚。1. 项目概述与整体设计思路1.1 为什么选择主从博弈框架先说说这类问题最本质的背景。传统电力系统调度默认假设所有参与者都是听从调度指令的单一整体用集中式优化就能拿到全局最优解。但综合能源系统发展到多主体阶段后情况就不一样了园区里有综合能源运营商楼下有商业用户、工业用户、居民用户不同主体都有自己的利益诉求。运营商想通过制定能源价格让自己的收益最大化用户呢则会根据运营商给出的价格调整自己的用电、用热行为——这就是需求响应。两者之间存在明显的先后顺序和利益冲突。运营商先定价格用户后调整用能计划这就天然构成一个Stackelberg博弈也就是主从博弈。主从博弈的好处是逻辑清晰领导者运营商拥有先动优势但同时必须考虑跟随者用户的理性响应曲线找到能让自己收益最大化的定价方案。相比集中式优化主从博弈更贴近实际的多主体交易场景也更符合当前电力市场改革后售电公司终端用户的竞争格局。从科研角度看这类模型的好处也显而易见可以写出双层优化形式上层是运营商的收益最大化下层是用户的用能成本最小化整个问题有明确的KKT条件可以去处理数学上非常成熟审稿人也买账。1.2 需求响应与电能交互在博弈模型中的定位这个题目里有两个交互值得单独强调。第一个是需求响应它连接的是运营商的价格策略和用户的用能行为第二个是电能交互它连接的是多个用能主体之间的功率流动。很多人做多主体优化时把需求响应看成负荷约束就完事了但真正的博弈模型里需求响应是下层用户问题的一部分。运营商定价之后用户的负荷曲线、储能充放电策略、甚至可转移负荷的启停时间都会跟着变变完之后又反过来影响运营商的收益和系统供需平衡。所以需求和响应不是一条静态曲线而是整个博弈的动态反馈环节。电能交互则是多主体场景特有的耦合。多个用户/微电网/能源站之间可以通过公共联络线进行电能买卖。交互功率的大小和方向受制于线路容量和运行约束而且交互电价同样可以由运营商制定或通过市场机制形成。加入电能交互的收益在于系统可以从各自为政变成互联互济高负荷时段能够互相支援整体运行成本自然下降。再往下拆这两块内容就会直接反映在数学模型的变量和约束中具体怎么建我在下一节详细展开。2. 核心数学模型拆解2.1 上层领导者综合能源运营商的调度模型我先给出运营商的典型设备构成。大多数学术文章和实际工程里的IES运营商主要设备包括CHP热电联产机组、燃气锅炉、电锅炉、储能电池、蓄热罐还可能配一部分分布式光伏和风机。运营商的目标函数一般是自己的净收益最大化常见形式是max F_leader Σ(售电收入 售热收入 售气收入 电能交互收入) - Σ(购电成本 购气成本 运行维护成本 储能折旧成本)其中售电/售热/售气收入是用户侧交易得到的电能交互收入是与外部配电网交易得到的。这个目标的本质是卖出去的钱减去买进来的钱和运行损耗的成本非常贴合售电公司的实际运营逻辑。约束条件要从多个维度搭能量平衡约束任意时段电力均满足CHP发电 PV输出 储能放电 外购电力 交互功率 用户电负荷 电锅炉用电 储能充电。热力平衡同理CHP余热 燃气锅炉 蓄热罐放热 热负荷 蓄热罐蓄热。设备运行约束CHP存在电出力上下限、爬坡速率约束同时热电比耦合燃气锅炉有热出力上下限储能电池有SOC递推方程和充放功率限制而且要避免同时充放电。交互功率约束与配电网的购售电功率有上限且为了安全考虑通常不允许在同一时段既买入又卖出。这些约束用MatlabYALMIP表达并不复杂最需要注意的是储能SOC的递推约束要写对时序以及CHP的热电耦合约束如果引入二进制变量问题规模会翻倍要有心理准备。2.2 跟随者用户侧需求响应模型下层用户的目标函数是自身用能成本最小约束条件则体现需求响应能力。在实际代码中常见的用户侧模型分三类我在实现时建议分开写清楚后面Debug会方便很多。第一类价格型需求响应Price-based DR这一类基于需求价格弹性矩阵。核心思路是用户在电价升高时降低用电需求电价降低时增加用电需求。写成约束时负荷变化量与电价变化量之间的关系为ΔL_i E_ii * (Δp_i / p_i0) * L_i0 Σ(E_ij * (Δp_j / p_j0) * L_i0)其中E_ii是自弹性系数E_ij是交叉弹性系数。弹性矩阵一般由统计数据或实测数据得到代码里可以直接写成常数矩阵然后作为参数传入。需要注意价格型需求响应一般只改变用户的用电量曲线形状不保证总用电量不变所以如果做电费中性假设无差异曲线分析需要额外加一个总用电量守恒约束。这一点很多初学的同行容易漏。第二类激励型需求响应Incentive-based DR这类需求响应通常指可中断负荷IL和直接负荷控制DLC。运营商与用户签订协议在系统高峰时段可以切除一部分负荷同时要支付相应的补偿费用。建模时引入二进制变量表示负荷是否被中断并加中断次数限制。最终在下层目标函数中用户会权衡削减负荷拿到的补偿与舒适度损失这个平衡其实靠价格和补偿系数调节。第三类可转移负荷Shiftable Load典型场景是电动汽车充电、洗衣机和消毒柜这类可延时负荷。约束条件为负荷必须在一个允许的时间窗口内完成并且有连续运行要求。建模时对每个时间窗口设置开始时间变量用二进制变量表示t时段是否启动。这类约束的线性化形式非常成熟代码写起来就是用big-M方法。我个人习惯把三类需求响应拆成三个函数文件——calc_pdr.m、calc_idr.m、calc_slr.m输入是电价序列和原始负荷参数输出是响应后的负荷序列。这样方便单独调试也可以灵活决定算例里到底激活哪几类需求响应。2.3 电能交互与耦合约束电能交互是多主体问题最核心的耦合部分也是最容易在编程时出bug的地方。先明确主体的边界。假设系统中有N个用户主体和一个运营商主体。用户之间通过公共母线建立联系运营商作为中间调控方组织交互功率。电能交互的建模通常用以下几条约束交互功率限值约束任意两个主体i和j之间的交互功率P_ex_ij(t)不能超过联络线允许容量P_ex_max。交互平衡约束总交互功率的代数和为零即ΣP_ex_ij(t) 0。这是保证多主体之间电能守恒的关键。交互功率方向约束若引入二进制变量区分买入/卖出状态就需要加入互斥约束避免同一时段既买入又卖出这既不合常理也容易导致求解无解。交互价格约束如果交互电价由运营商制定价格要在购电成本和售电价格之间以保证运营商有利可图且用户不会亏本买卖。通常设为购电成本的1.1~1.3倍。从编程角度说交互变量是N×T维矩阵维度不大但会带来强烈的耦合性。如果问题是集中式求解这些约束直接堆到总约束里就行如果走分布式求解则要把交互功率看成边界变量通常会在目标函数里加一个二次惩罚项来加速收敛。2.4 主从博弈求解的两条技术路线这道题看着复杂实际上业界主流的求解路线分两种我把各自的特点讲透。路线一KKT条件强对偶定理MPEC方法思路很简单下层问题是用户的最小化问题如果它是线性规划LP或凸二次规划QP就可以写出它的KKT条件然后用KKT条件替换下层优化问题将双层问题转化为单层MPEC带均衡约束的数学规划。但此时目标函数中出现电价×负荷这类双线性项需要借助强对偶定理或者利用KKT互补条件做线性化最终转化为混合整数线性规划MILP调用Cplex/Gurobi求解。这条路线胜在求解效率高、最优性有保证学术文章里用得最多。缺点是推导过程繁琐涉及大量补充松弛变量手推很容易出错。我的经验是用Matlab的符号工具箱做推导辅助或者先用小规模样例验证KKT条件正确性再扩展到完整算例。路线二分布式迭代方法启发式嵌套优化这条路线是直接把博弈写成迭代形式外层用启发式算法粒子群、遗传算法、差分进化搜索上层变量价格、交互电价内层用Cplex/Yalmip求解下层用户的优化问题然后把用户的响应结果反馈到上层计算上层目标函数。迭代收敛到的平衡点就是Stackelberg均衡。这条路线好写好理解但收敛性会让人头疼尤其是粒子群参数没调好时价格迭代曲线经常震荡。后面第4章我会详细讲我怎么处理这个问题的。总体建议如果对求解精度有硬要求或者被审稿人质问是否全局最优尽量走KKT路线如果只是为了快速搭Demo、验证策略有效性分布式迭代够用。3. Matlab仿真实现全流程3.1 代码总体框架一个标准的主从博弈综合能源系统Matlab工程我推荐按以下目录组织IES_Stackelberg/ ├── main.m ├── data/ │ ├── load_data.m % 负荷与可再生能源曲线 │ ├── price_data.m % 分时电价、气价参数 │ └── system_params.m % 设备参数 ├── models/ │ ├── upper_leader.m % 上层运营商优化模型 │ ├── lower_user.m % 下层用户需求响应模型 │ └── interaction.m % 电能交互与耦合约束 ├── solvers/ │ ├── solve_kkt.m % KKT单层化求解 │ └── solve_iter.m % 分布式迭代求解 └── utils/ ├── plot_results.m % 绘图 └── check_converge.m% 收敛性检测main.m是整个程序的主入口统一负责数据加载、模型组装、求解和结果输出。把数据和模型分离是很多Matlab新手容易忽略的一点实际项目一跑起来你会发现改参数、换数据是常态不分离的后果就是改一个参数要翻半天代码。3.2 参数定义与数据初始化我在代码里直接用结构体管理参数这比散落的变量清晰得多。以机组和储能为例%% 系统参数设置 devices struct(); devices.CHP.e_min 100; % CHP电出力下限kW devices.CHP.e_max 800; % CHP电出力上限kW devices.CHP.hpr 1.3; % 热电比实际项目中可取0.9~1.5 devices.CHP.eta_e 0.35; % 发电效率 devices.CHP.eta_h 0.45; % 供热效率 devices.CHP.ramp 200; % 爬坡速率kW/h devices.GB.h_min 0; devices.GB.h_max 500; % 燃气锅炉热出力上限kW devices.GB.eta 0.9; % 热效率 devices.ES.cap 1000; % 储能电池容量kWh devices.ES.p_max 200; % 最大充放电功率kW devices.ES.eta_ch 0.95; % 充电效率 devices.ES.eta_dis 0.95; % 放电效率 devices.ES.soc_min 0.1; % SOC下限 devices.ES.soc_max 0.9; % SOC上限 devices.ES.soc_0 0.2; % 初始SOC时间尺度方面我用的标准调度周期是24小时步长为1小时这是这个领域最常见的配置。如果你做的是日内滚动优化也可以改成96个时段每15分钟一个点这时候变量维度会直接乘以4求解时间呈指数上涨要有心理准备。3.3 需求响应模型的编码实例以价格型需求响应为例这里我直接展示核心代码编写思路。弹性矩阵可以直接从外部数据文件读入function L_new price_dr(L0, p, p_ref, E) % 输入 % L0 - 基准负荷1×T向量 % p - 实际电价1×T向量 % p_ref- 参考电价基准情景电价1×T向量 % E - 需求价格弹性矩阵T×T % 输出 % L_new- 需求响应后的负荷1×T向量 T length(L0); delta_p_ratio (p - p_ref) ./ p_ref; % 电价变化率向量 L_new L0; for i 1:T flexible 0; for j 1:T flexible flexible E(i,j) * delta_p_ratio(j) * L0(i); end L_new(i) L0(i) * (1 flexible); end % 可选总用电量守恒约束修正保证用户总用电量不变 % L_new L_new * sum(L0) / sum(L_new); end这里有个很容易踩的细节弹性矩阵E的主对角线元素自弹性通常为负值交叉弹性非对角线元素通常为正值这代表自身电价上升削减自身负荷其他时段电价上升时转移负荷到本时段。如果你写反了符号那响应出来的负荷曲线会完全反直觉甚至会放大峰谷差。可中断负荷激励型DR的代码我用YALMIP的binvar来定义中断状态变量T 24; I_il binvar(1, T); % 1表示负荷被中断 P_il_max 50; % 每个时段可中断负荷上限(kW) N_il_max 4; % 最多中断的时段数 Constraints [Constraints, 0 P_il I_il * P_il_max]; Constraints [Constraints, sum(I_il) N_il_max];注意你有没有想过为什么要限制N_il_max如果运营商给的中断补偿足够高用户恨不得每个时段都被拉闸启动次数限制本质上是模拟用户舒适度和设备寿命的现实约束不做限制的话模型过于理想化。3.4 博弈迭代主循环实现分布式迭代是理解整套博弈计算流程的最好入口代码逻辑也很直白我贴一个框架% 初始化 iter 0; lambda 0.5 * (c_buy s_sell); % 初始电价取购售电价的中间值 alpha 0.3; % 价格更新步长需要调参 tol 1e-4; max_iter 30; err inf; while err tol iter max_iter iter iter 1; % 第一步已知运营商电价求解下层用户问题 [L_user, cost_user] lower_user_solve(lambda, price_data, device_data); % 第二步将用户响应结果传回上层求解运营商优化问题 [P_operator, revenue_operator] upper_leader_solve(L_user, price_data, device_data); % 第三步更新价格这里用最经典的迭代公式 lambda_new lambda alpha * (L_user - L_target); % 第四步检测收敛性 err max(abs(lambda_new - lambda)); lambda lambda_new; % 记录每一轮的结果方便画收敛曲线 history.lambda(iter, :) lambda; history.revenue(iter) revenue_operator; history.cost(iter) cost_user; end这里最核心也最玄学的部分就是第三步的价格更新。我一开始用固定步长时价格迭代曲线像心电图一样疯狂震荡后来换了阻尼步长效果立竿见影。所谓阻尼步长就是随着迭代次数增加逐渐缩小更新幅度相当于每次价格调整得越来越谨慎alpha alpha0 / sqrt(iter); % alpha0为初始步长再进阶一点的做法是用平均值法每次取前几轮的平均值作为新价格或者借鉴次梯度类算法的平滑技巧。具体怎么选往往要看你的电价目标曲线是什么形状这个没有通用的万能答案建议多看几个典型场景的收敛曲线再定。3.5 结果绘图与分析Matlab出图是整套代码最后也是最能直观体现成果的环节。我一般画这么几张图第一张是系统电功率平衡图。用堆叠面积图展示CHP、光伏、储能放电、外购电、交互功率随时间的变化直观检验任何时候功率都是平衡的这个图在写论文、做汇报时是必备的。第二张是价格与负荷的响应曲线对比图。把无需求响应时的负荷曲线和有需求响应时的负荷曲线画在一起叠加分时电价作为右侧坐标轴一眼就能看出需求响应是否削峰填谷、用户是否做到合理避峰。第三张是储能SOC和充放电功率图。储能调度策略是否合理看这个图最直观SOC曲线应该是平滑过渡的不能有突兀跳变。如果SOC从0.1瞬间跳到0.9说明程序里的充放电功率约束或者时序关系写错了。第四张是博弈收敛曲线。把每一轮迭代的运营商收益或电价画出来呈现单调收敛趋势这是审稿人和答辩老师最喜欢看到的一张图也是验证模型求解可靠性的直接证据。4. 频率踩坑记录与排查技巧4.1 迭代不收敛价格震荡的真相分布式迭代最常见的翻车现场就是价格震荡不收敛。我踩过最深的一次坑是在做多主体电能交互时直接把交互功率增量乘以固定步长叠加到电价上结果50次迭代下来电价不是收敛到均衡而是呈现周期为2的振荡——一会儿极度抬高、一会儿极度压低完全没法用。后来我梳理清楚原因了固定步长太大时价格更新会越过均衡点步长太小时又收敛得太慢。而且多主体交互让响应函数带有一定的惯性价格更新必须考虑前一时刻的交互功率方向。我的解决办法分两步。第一步改步长策略采用衰减步长方案让前期快速逼近、后期精细收敛。第二步加低通滤波即每次更新的价格取当前迭代计算值和上一轮价格的加权平均权重可以取0.3/0.7效果立竿见影。lambda_new 0.3 * lambda_calc 0.7 * lambda_prev;有人管这叫价格平滑本质上就是抑制高频振荡思路跟控制理论里的阻尼器差不多。4.2 求解器报无可行解的排查清单无论你用的是YalmipCplex还是Gurobi遇到infeasible problem都是家常便饭。我总结了一套排查顺序能解决90%的问题第一优先级查等式约束。能量守恒约束其实最容易被写错比如储能充放电同时发生的隐性bug会导致等式约束永远无法满足因为同一个时刻既要充电又要放电超出物理可能。这时候给储能加一个充放互斥约束logical constraint或二进制变量就能解决。第二优先级查SOC初始条件。储能初始SOC设为0.2但终态SOC要求等于0.2时如果系统净电量不足或者充放电效率损失被忽略很容易导致无解。建议先去掉SOC终值约束跑一遍如果解出来了说明是能量不平衡问题再调设备参数或者加惩罚项。第三优先级查交互功率约束。多主体同时买卖会导致功率平衡约束被破坏。一个快速定位方法是把交互约束全部注释掉如果问题可解就逐渐加上限制二分定位到出问题的约束。4.3 运行时间爆炸的处理思路综合能源系统的维度一旦上去Cplex求解时间很容易从秒变成小时。我实测过一个24节点、24时段的算例直接求解MILP用了近40分钟这在参数调试阶段根本无法接受。优化手段从易到难排列如果上层变量是多维连续变量先用线性规划松弛版本测速大量使用big-M时M值不宜取得过大过大会让松弛的可行域过大求解器分支定界效率暴跌工程上取设备最大功率的2倍左右足够如果可以接受精度损失将原问题分解为多个子问题依次求解用热启动参数传入初始解Cplex能大幅减少分支切割的工作量。4.4 结果的合理性自查代码能跑出结果不等于结果是正确的。我每次跑完程序至少会做三件事第一检查设备出力是否在上下限边界内重点看CHP和GB有没有超过容量输出。如果某设备长期运行在最大值而系统仍然满足不了负荷那说明设备容量配置不合理或者需求响应参数设置太保守。第二检查储能SOC曲线是否符合物理规律。SOC变化趋势应该与充放电功率曲线一致充多放少则SOC上升反之下降。如果SOC曲线毛刺很多大概率是时序索引写错了比如用t和t-1时出了边界。第三核对总账。把运营商的总收入、总成本、净收益分别打印出来检查利润是否合理。如果运营商亏本运行大概率是交互电价或售能价格设置不合理如果用户成本比基准情景还高那需求响应就失去意义了。5. 实操经验与个人心得这套主从博弈优化调度模型我从写第一个简陋单用户版本到现在加入多主体交互前后折腾了大半年。最大的体会是数学模型写出来再漂亮落到Matlab代码上都是从一个个细节堆起来的。弹性矩阵的正负号、SOC递推的索引、价格迭代的步长任何一个地方出错结果看起来都会跟合理差那么一点但又很难一眼看出错在哪。这里分享一个我后来固定下来的调试习惯先跑一个最小的2时段、2主体诱导弹例把每一个约束的残差打印出来手算一遍验证结果完全对得上再扩展到完整24时段算例。虽然前期多花一小时但后面Debug省出来的时间可能是一天。再有一个建议是主从博弈模型的验证不能只看最终目标函数值一定要画迭代过程的收敛曲线。收敛曲线既证明了算法有效性也方便观察是否存在局部振荡。那些看着数字合理但收敛曲线乱跳的结果在答辩和审稿阶段很容易被打回来。如果后续你想在这个方向上继续扩展可以考虑把单层博弈拓展为主从-对等混合博弈即运营商与用户之间是主从关系但多个用户之间是对等关系需要同时求解Stackelberg均衡和Nash均衡这也是现在比较热门的交叉方向。模型复杂度会提升一个量级求解方法也需要调整但底层这套Matlab框架完全可以复用变量定义和约束组织方式不需要推倒重来。 多主体综合能源系统主从博弈优化调度Matlab实现
我先交代一下背景。这个项目是我在实际课题里被问到最多的一类问题——多主体综合能源系统、需求响应、电能交互、主从博弈四个关键词堆在一起看着像四座山但真正落地成Matlab代码时核心就一句话谁先出招谁后应对怎么收敛到均衡点。这篇博文我尽量不把通用建模教科书复述一遍而是把整套思路、建模细节、Matlab实现的关键节点以及我实际跑代码时踩过的坑都摊开来讲清楚。1. 项目概述与整体设计思路1.1 为什么选择主从博弈框架先说说这类问题最本质的背景。传统电力系统调度默认假设所有参与者都是听从调度指令的单一整体用集中式优化就能拿到全局最优解。但综合能源系统发展到多主体阶段后情况就不一样了园区里有综合能源运营商楼下有商业用户、工业用户、居民用户不同主体都有自己的利益诉求。运营商想通过制定能源价格让自己的收益最大化用户呢则会根据运营商给出的价格调整自己的用电、用热行为——这就是需求响应。两者之间存在明显的先后顺序和利益冲突。运营商先定价格用户后调整用能计划这就天然构成一个Stackelberg博弈也就是主从博弈。主从博弈的好处是逻辑清晰领导者运营商拥有先动优势但同时必须考虑跟随者用户的理性响应曲线找到能让自己收益最大化的定价方案。相比集中式优化主从博弈更贴近实际的多主体交易场景也更符合当前电力市场改革后售电公司终端用户的竞争格局。从科研角度看这类模型的好处也显而易见可以写出双层优化形式上层是运营商的收益最大化下层是用户的用能成本最小化整个问题有明确的KKT条件可以去处理数学上非常成熟审稿人也买账。1.2 需求响应与电能交互在博弈模型中的定位这个题目里有两个交互值得单独强调。第一个是需求响应它连接的是运营商的价格策略和用户的用能行为第二个是电能交互它连接的是多个用能主体之间的功率流动。很多人做多主体优化时把需求响应看成负荷约束就完事了但真正的博弈模型里需求响应是下层用户问题的一部分。运营商定价之后用户的负荷曲线、储能充放电策略、甚至可转移负荷的启停时间都会跟着变变完之后又反过来影响运营商的收益和系统供需平衡。所以需求和响应不是一条静态曲线而是整个博弈的动态反馈环节。电能交互则是多主体场景特有的耦合。多个用户/微电网/能源站之间可以通过公共联络线进行电能买卖。交互功率的大小和方向受制于线路容量和运行约束而且交互电价同样可以由运营商制定或通过市场机制形成。加入电能交互的收益在于系统可以从各自为政变成互联互济高负荷时段能够互相支援整体运行成本自然下降。再往下拆这两块内容就会直接反映在数学模型的变量和约束中具体怎么建我在下一节详细展开。2. 核心数学模型拆解2.1 上层领导者综合能源运营商的调度模型我先给出运营商的典型设备构成。大多数学术文章和实际工程里的IES运营商主要设备包括CHP热电联产机组、燃气锅炉、电锅炉、储能电池、蓄热罐还可能配一部分分布式光伏和风机。运营商的目标函数一般是自己的净收益最大化常见形式是max F_leader Σ(售电收入 售热收入 售气收入 电能交互收入) - Σ(购电成本 购气成本 运行维护成本 储能折旧成本)其中售电/售热/售气收入是用户侧交易得到的电能交互收入是与外部配电网交易得到的。这个目标的本质是卖出去的钱减去买进来的钱和运行损耗的成本非常贴合售电公司的实际运营逻辑。约束条件要从多个维度搭能量平衡约束任意时段电力均满足CHP发电 PV输出 储能放电 外购电力 交互功率 用户电负荷 电锅炉用电 储能充电。热力平衡同理CHP余热 燃气锅炉 蓄热罐放热 热负荷 蓄热罐蓄热。设备运行约束CHP存在电出力上下限、爬坡速率约束同时热电比耦合燃气锅炉有热出力上下限储能电池有SOC递推方程和充放功率限制而且要避免同时充放电。交互功率约束与配电网的购售电功率有上限且为了安全考虑通常不允许在同一时段既买入又卖出。这些约束用MatlabYALMIP表达并不复杂最需要注意的是储能SOC的递推约束要写对时序以及CHP的热电耦合约束如果引入二进制变量问题规模会翻倍要有心理准备。2.2 跟随者用户侧需求响应模型下层用户的目标函数是自身用能成本最小约束条件则体现需求响应能力。在实际代码中常见的用户侧模型分三类我在实现时建议分开写清楚后面Debug会方便很多。第一类价格型需求响应Price-based DR这一类基于需求价格弹性矩阵。核心思路是用户在电价升高时降低用电需求电价降低时增加用电需求。写成约束时负荷变化量与电价变化量之间的关系为ΔL_i E_ii * (Δp_i / p_i0) * L_i0 Σ(E_ij * (Δp_j / p_j0) * L_i0)其中E_ii是自弹性系数E_ij是交叉弹性系数。弹性矩阵一般由统计数据或实测数据得到代码里可以直接写成常数矩阵然后作为参数传入。需要注意价格型需求响应一般只改变用户的用电量曲线形状不保证总用电量不变所以如果做电费中性假设无差异曲线分析需要额外加一个总用电量守恒约束。这一点很多初学的同行容易漏。第二类激励型需求响应Incentive-based DR这类需求响应通常指可中断负荷IL和直接负荷控制DLC。运营商与用户签订协议在系统高峰时段可以切除一部分负荷同时要支付相应的补偿费用。建模时引入二进制变量表示负荷是否被中断并加中断次数限制。最终在下层目标函数中用户会权衡削减负荷拿到的补偿与舒适度损失这个平衡其实靠价格和补偿系数调节。第三类可转移负荷Shiftable Load典型场景是电动汽车充电、洗衣机和消毒柜这类可延时负荷。约束条件为负荷必须在一个允许的时间窗口内完成并且有连续运行要求。建模时对每个时间窗口设置开始时间变量用二进制变量表示t时段是否启动。这类约束的线性化形式非常成熟代码写起来就是用big-M方法。我个人习惯把三类需求响应拆成三个函数文件——calc_pdr.m、calc_idr.m、calc_slr.m输入是电价序列和原始负荷参数输出是响应后的负荷序列。这样方便单独调试也可以灵活决定算例里到底激活哪几类需求响应。2.3 电能交互与耦合约束电能交互是多主体问题最核心的耦合部分也是最容易在编程时出bug的地方。先明确主体的边界。假设系统中有N个用户主体和一个运营商主体。用户之间通过公共母线建立联系运营商作为中间调控方组织交互功率。电能交互的建模通常用以下几条约束交互功率限值约束任意两个主体i和j之间的交互功率P_ex_ij(t)不能超过联络线允许容量P_ex_max。交互平衡约束总交互功率的代数和为零即ΣP_ex_ij(t) 0。这是保证多主体之间电能守恒的关键。交互功率方向约束若引入二进制变量区分买入/卖出状态就需要加入互斥约束避免同一时段既买入又卖出这既不合常理也容易导致求解无解。交互价格约束如果交互电价由运营商制定价格要在购电成本和售电价格之间以保证运营商有利可图且用户不会亏本买卖。通常设为购电成本的1.1~1.3倍。从编程角度说交互变量是N×T维矩阵维度不大但会带来强烈的耦合性。如果问题是集中式求解这些约束直接堆到总约束里就行如果走分布式求解则要把交互功率看成边界变量通常会在目标函数里加一个二次惩罚项来加速收敛。2.4 主从博弈求解的两条技术路线这道题看着复杂实际上业界主流的求解路线分两种我把各自的特点讲透。路线一KKT条件强对偶定理MPEC方法思路很简单下层问题是用户的最小化问题如果它是线性规划LP或凸二次规划QP就可以写出它的KKT条件然后用KKT条件替换下层优化问题将双层问题转化为单层MPEC带均衡约束的数学规划。但此时目标函数中出现电价×负荷这类双线性项需要借助强对偶定理或者利用KKT互补条件做线性化最终转化为混合整数线性规划MILP调用Cplex/Gurobi求解。这条路线胜在求解效率高、最优性有保证学术文章里用得最多。缺点是推导过程繁琐涉及大量补充松弛变量手推很容易出错。我的经验是用Matlab的符号工具箱做推导辅助或者先用小规模样例验证KKT条件正确性再扩展到完整算例。路线二分布式迭代方法启发式嵌套优化这条路线是直接把博弈写成迭代形式外层用启发式算法粒子群、遗传算法、差分进化搜索上层变量价格、交互电价内层用Cplex/Yalmip求解下层用户的优化问题然后把用户的响应结果反馈到上层计算上层目标函数。迭代收敛到的平衡点就是Stackelberg均衡。这条路线好写好理解但收敛性会让人头疼尤其是粒子群参数没调好时价格迭代曲线经常震荡。后面第4章我会详细讲我怎么处理这个问题的。总体建议如果对求解精度有硬要求或者被审稿人质问是否全局最优尽量走KKT路线如果只是为了快速搭Demo、验证策略有效性分布式迭代够用。3. Matlab仿真实现全流程3.1 代码总体框架一个标准的主从博弈综合能源系统Matlab工程我推荐按以下目录组织IES_Stackelberg/ ├── main.m ├── data/ │ ├── load_data.m % 负荷与可再生能源曲线 │ ├── price_data.m % 分时电价、气价参数 │ └── system_params.m % 设备参数 ├── models/ │ ├── upper_leader.m % 上层运营商优化模型 │ ├── lower_user.m % 下层用户需求响应模型 │ └── interaction.m % 电能交互与耦合约束 ├── solvers/ │ ├── solve_kkt.m % KKT单层化求解 │ └── solve_iter.m % 分布式迭代求解 └── utils/ ├── plot_results.m % 绘图 └── check_converge.m% 收敛性检测main.m是整个程序的主入口统一负责数据加载、模型组装、求解和结果输出。把数据和模型分离是很多Matlab新手容易忽略的一点实际项目一跑起来你会发现改参数、换数据是常态不分离的后果就是改一个参数要翻半天代码。3.2 参数定义与数据初始化我在代码里直接用结构体管理参数这比散落的变量清晰得多。以机组和储能为例%% 系统参数设置 devices struct(); devices.CHP.e_min 100; % CHP电出力下限kW devices.CHP.e_max 800; % CHP电出力上限kW devices.CHP.hpr 1.3; % 热电比实际项目中可取0.9~1.5 devices.CHP.eta_e 0.35; % 发电效率 devices.CHP.eta_h 0.45; % 供热效率 devices.CHP.ramp 200; % 爬坡速率kW/h devices.GB.h_min 0; devices.GB.h_max 500; % 燃气锅炉热出力上限kW devices.GB.eta 0.9; % 热效率 devices.ES.cap 1000; % 储能电池容量kWh devices.ES.p_max 200; % 最大充放电功率kW devices.ES.eta_ch 0.95; % 充电效率 devices.ES.eta_dis 0.95; % 放电效率 devices.ES.soc_min 0.1; % SOC下限 devices.ES.soc_max 0.9; % SOC上限 devices.ES.soc_0 0.2; % 初始SOC时间尺度方面我用的标准调度周期是24小时步长为1小时这是这个领域最常见的配置。如果你做的是日内滚动优化也可以改成96个时段每15分钟一个点这时候变量维度会直接乘以4求解时间呈指数上涨要有心理准备。3.3 需求响应模型的编码实例以价格型需求响应为例这里我直接展示核心代码编写思路。弹性矩阵可以直接从外部数据文件读入function L_new price_dr(L0, p, p_ref, E) % 输入 % L0 - 基准负荷1×T向量 % p - 实际电价1×T向量 % p_ref- 参考电价基准情景电价1×T向量 % E - 需求价格弹性矩阵T×T % 输出 % L_new- 需求响应后的负荷1×T向量 T length(L0); delta_p_ratio (p - p_ref) ./ p_ref; % 电价变化率向量 L_new L0; for i 1:T flexible 0; for j 1:T flexible flexible E(i,j) * delta_p_ratio(j) * L0(i); end L_new(i) L0(i) * (1 flexible); end % 可选总用电量守恒约束修正保证用户总用电量不变 % L_new L_new * sum(L0) / sum(L_new); end这里有个很容易踩的细节弹性矩阵E的主对角线元素自弹性通常为负值交叉弹性非对角线元素通常为正值这代表自身电价上升削减自身负荷其他时段电价上升时转移负荷到本时段。如果你写反了符号那响应出来的负荷曲线会完全反直觉甚至会放大峰谷差。可中断负荷激励型DR的代码我用YALMIP的binvar来定义中断状态变量T 24; I_il binvar(1, T); % 1表示负荷被中断 P_il_max 50; % 每个时段可中断负荷上限(kW) N_il_max 4; % 最多中断的时段数 Constraints [Constraints, 0 P_il I_il * P_il_max]; Constraints [Constraints, sum(I_il) N_il_max];注意你有没有想过为什么要限制N_il_max如果运营商给的中断补偿足够高用户恨不得每个时段都被拉闸启动次数限制本质上是模拟用户舒适度和设备寿命的现实约束不做限制的话模型过于理想化。3.4 博弈迭代主循环实现分布式迭代是理解整套博弈计算流程的最好入口代码逻辑也很直白我贴一个框架% 初始化 iter 0; lambda 0.5 * (c_buy s_sell); % 初始电价取购售电价的中间值 alpha 0.3; % 价格更新步长需要调参 tol 1e-4; max_iter 30; err inf; while err tol iter max_iter iter iter 1; % 第一步已知运营商电价求解下层用户问题 [L_user, cost_user] lower_user_solve(lambda, price_data, device_data); % 第二步将用户响应结果传回上层求解运营商优化问题 [P_operator, revenue_operator] upper_leader_solve(L_user, price_data, device_data); % 第三步更新价格这里用最经典的迭代公式 lambda_new lambda alpha * (L_user - L_target); % 第四步检测收敛性 err max(abs(lambda_new - lambda)); lambda lambda_new; % 记录每一轮的结果方便画收敛曲线 history.lambda(iter, :) lambda; history.revenue(iter) revenue_operator; history.cost(iter) cost_user; end这里最核心也最玄学的部分就是第三步的价格更新。我一开始用固定步长时价格迭代曲线像心电图一样疯狂震荡后来换了阻尼步长效果立竿见影。所谓阻尼步长就是随着迭代次数增加逐渐缩小更新幅度相当于每次价格调整得越来越谨慎alpha alpha0 / sqrt(iter); % alpha0为初始步长再进阶一点的做法是用平均值法每次取前几轮的平均值作为新价格或者借鉴次梯度类算法的平滑技巧。具体怎么选往往要看你的电价目标曲线是什么形状这个没有通用的万能答案建议多看几个典型场景的收敛曲线再定。3.5 结果绘图与分析Matlab出图是整套代码最后也是最能直观体现成果的环节。我一般画这么几张图第一张是系统电功率平衡图。用堆叠面积图展示CHP、光伏、储能放电、外购电、交互功率随时间的变化直观检验任何时候功率都是平衡的这个图在写论文、做汇报时是必备的。第二张是价格与负荷的响应曲线对比图。把无需求响应时的负荷曲线和有需求响应时的负荷曲线画在一起叠加分时电价作为右侧坐标轴一眼就能看出需求响应是否削峰填谷、用户是否做到合理避峰。第三张是储能SOC和充放电功率图。储能调度策略是否合理看这个图最直观SOC曲线应该是平滑过渡的不能有突兀跳变。如果SOC从0.1瞬间跳到0.9说明程序里的充放电功率约束或者时序关系写错了。第四张是博弈收敛曲线。把每一轮迭代的运营商收益或电价画出来呈现单调收敛趋势这是审稿人和答辩老师最喜欢看到的一张图也是验证模型求解可靠性的直接证据。4. 频率踩坑记录与排查技巧4.1 迭代不收敛价格震荡的真相分布式迭代最常见的翻车现场就是价格震荡不收敛。我踩过最深的一次坑是在做多主体电能交互时直接把交互功率增量乘以固定步长叠加到电价上结果50次迭代下来电价不是收敛到均衡而是呈现周期为2的振荡——一会儿极度抬高、一会儿极度压低完全没法用。后来我梳理清楚原因了固定步长太大时价格更新会越过均衡点步长太小时又收敛得太慢。而且多主体交互让响应函数带有一定的惯性价格更新必须考虑前一时刻的交互功率方向。我的解决办法分两步。第一步改步长策略采用衰减步长方案让前期快速逼近、后期精细收敛。第二步加低通滤波即每次更新的价格取当前迭代计算值和上一轮价格的加权平均权重可以取0.3/0.7效果立竿见影。lambda_new 0.3 * lambda_calc 0.7 * lambda_prev;有人管这叫价格平滑本质上就是抑制高频振荡思路跟控制理论里的阻尼器差不多。4.2 求解器报无可行解的排查清单无论你用的是YalmipCplex还是Gurobi遇到infeasible problem都是家常便饭。我总结了一套排查顺序能解决90%的问题第一优先级查等式约束。能量守恒约束其实最容易被写错比如储能充放电同时发生的隐性bug会导致等式约束永远无法满足因为同一个时刻既要充电又要放电超出物理可能。这时候给储能加一个充放互斥约束logical constraint或二进制变量就能解决。第二优先级查SOC初始条件。储能初始SOC设为0.2但终态SOC要求等于0.2时如果系统净电量不足或者充放电效率损失被忽略很容易导致无解。建议先去掉SOC终值约束跑一遍如果解出来了说明是能量不平衡问题再调设备参数或者加惩罚项。第三优先级查交互功率约束。多主体同时买卖会导致功率平衡约束被破坏。一个快速定位方法是把交互约束全部注释掉如果问题可解就逐渐加上限制二分定位到出问题的约束。4.3 运行时间爆炸的处理思路综合能源系统的维度一旦上去Cplex求解时间很容易从秒变成小时。我实测过一个24节点、24时段的算例直接求解MILP用了近40分钟这在参数调试阶段根本无法接受。优化手段从易到难排列如果上层变量是多维连续变量先用线性规划松弛版本测速大量使用big-M时M值不宜取得过大过大会让松弛的可行域过大求解器分支定界效率暴跌工程上取设备最大功率的2倍左右足够如果可以接受精度损失将原问题分解为多个子问题依次求解用热启动参数传入初始解Cplex能大幅减少分支切割的工作量。4.4 结果的合理性自查代码能跑出结果不等于结果是正确的。我每次跑完程序至少会做三件事第一检查设备出力是否在上下限边界内重点看CHP和GB有没有超过容量输出。如果某设备长期运行在最大值而系统仍然满足不了负荷那说明设备容量配置不合理或者需求响应参数设置太保守。第二检查储能SOC曲线是否符合物理规律。SOC变化趋势应该与充放电功率曲线一致充多放少则SOC上升反之下降。如果SOC曲线毛刺很多大概率是时序索引写错了比如用t和t-1时出了边界。第三核对总账。把运营商的总收入、总成本、净收益分别打印出来检查利润是否合理。如果运营商亏本运行大概率是交互电价或售能价格设置不合理如果用户成本比基准情景还高那需求响应就失去意义了。5. 实操经验与个人心得这套主从博弈优化调度模型我从写第一个简陋单用户版本到现在加入多主体交互前后折腾了大半年。最大的体会是数学模型写出来再漂亮落到Matlab代码上都是从一个个细节堆起来的。弹性矩阵的正负号、SOC递推的索引、价格迭代的步长任何一个地方出错结果看起来都会跟合理差那么一点但又很难一眼看出错在哪。这里分享一个我后来固定下来的调试习惯先跑一个最小的2时段、2主体诱导弹例把每一个约束的残差打印出来手算一遍验证结果完全对得上再扩展到完整24时段算例。虽然前期多花一小时但后面Debug省出来的时间可能是一天。再有一个建议是主从博弈模型的验证不能只看最终目标函数值一定要画迭代过程的收敛曲线。收敛曲线既证明了算法有效性也方便观察是否存在局部振荡。那些看着数字合理但收敛曲线乱跳的结果在答辩和审稿阶段很容易被打回来。如果后续你想在这个方向上继续扩展可以考虑把单层博弈拓展为主从-对等混合博弈即运营商与用户之间是主从关系但多个用户之间是对等关系需要同时求解Stackelberg均衡和Nash均衡这也是现在比较热门的交叉方向。模型复杂度会提升一个量级求解方法也需要调整但底层这套Matlab框架完全可以复用变量定义和约束组织方式不需要推倒重来。