我先说下这个项目的来龙去脉。这两年做综合能源系统优化的人越来越多尤其“电热综合能源系统”这个方向几乎是每个能源课题组都绕不开的题目。但很多新手拿到题后第一件事就是到处找代码、找模型找回来一堆文件却不知道怎么改、怎么跑、怎么看结果。这篇博文我就拿一个实际做过的主从博弈电热综合能源系统动态定价与能量管理项目来讲把MATLAB里从建模、求解到结果分析的完整链路拆开揉碎重点说清楚每个环节为什么这么做、代码哪里容易卡壳、调参有哪些门道。不管你是刚接触博弈论做优化的研究生还是已经跑过一些集中式调度想往分布式决策方向转的工程师这篇内容都值得花十分钟看完里面有大量实操细节是常规代码注释里根本不会写的。1. 项目概述与总体思路1.1 电热综合能源系统的核心场景所谓电热综合能源系统简单说就是把电力系统和热力系统放在一个框架里统一分析和优化。实际系统里通常包括上级电网、配电网、CHP热电联产机组、电锅炉、燃气锅炉、储能设备、热网管道以及楼宇、居民区、工业厂房等多种类型的用户。电和热在这里不是孤立存在的CHP一边发电一边产热电锅炉把多余的电转化为热热网再把热源送到用户端这就形成了强耦合关系。这种耦合带来的最大问题是你没法单独给电定价、单独给热定价因为用户用一度电可能导致热网负荷变化热网负荷变化又反过来影响电力系统的运行成本。传统的固定电价、固定热价模式在这种场景下效率很低峰时大家一起抢电谷时没人用热设备利用率上不去系统整体经济性也差。动态定价要解决的就是这个问题——让价格跟着供需关系走引导用户主动调整用能行为。1.2 为什么是主从博弈而不是集中优化很多第一次接触这个方向的人会问既然要优化直接用集中式调度不是更简单吗把所有设备、所有用户的数据收集上来建立一个统一的目标函数用MATLAB的linprog或者yalmip直接求解得到全局最优方案为什么非要引入主从博弈答案在于系统结构。集中式优化假设是一个“全能调度中心”能拿到所有信息并强制下发指令但在真实的电热综合能源系统中能源服务商和用户是不同利益主体。服务商想多赚利润用户想少花用能成本两边目标天然冲突。你直接替用户做决策用户不一定执行即使用户执行了也会因为缺乏价格激励而产生抵触。主从博弈Stackelberg game恰好是描述这种“一主多从”分层决策关系的标准框架。类比一下就很好理解主从博弈就像超市定价和消费者购物之间的关系。超市先定一个价格消费者根据价格决定买多少消费者买得多了超市发现定价太低下一轮涨价消费者又减少购买。来回博弈几次最终稳定在一个双方都能接受的价格和购买量上这就是纳什均衡。能源服务商就是超市用户就是消费者只不过这里交易的商品是电和热而且还要叠加复杂的设备运行约束。这类模型在MATLAB里不像单纯线性规划那样一个命令就能解完需要显式地构建上层价格决策变量、下层用户优化模型再设计迭代求解流程。核心工作量不在“解一个优化问题”而在“把两层优化问题正确耦合起来”。2. 主从博弈模型的数学构建2.1 上层能源服务商的动态定价模型上层的角色是能源服务商决策变量是各时段的购售电价和购售热价。目标函数为自身收益最大化主要包括售电收入、售热收入、向上级电网购电成本、燃气成本、设备运维成本这几块。用数学语言说假设调度周期为24小时时段集合记为T。上层目标可以写为% 上层目标函数示意非完整代码 % 收入电价*电负荷 热价*热负荷 % 成本购电成本 燃料成本 运维成本 Revenue sum( pi_e(t) .* P_load_total(t) pi_h(t) .* H_load_total(t) ); Cost sum( c_buy(t) .* P_grid(t) c_gas .* F_gas(t) c_om .* (P_chp(t) H_chp(t)) ); Profit Revenue - Cost;其中pi_e和pi_h是服务商制定的动态价格P_load_total和H_load_total是所有用户的总电负荷和总热负荷这两个量不是人为设定的而是下层用户优化后的结果。这就是主从博弈的关键上层只能通过调整价格间接影响负荷不能直接指定负荷。上层还需要满足一些运行约束比如CHP的发电功率与产热功率的可行域、电锅炉的爬坡约束、储能设备的充放电约束等。这些约束本质上刻画了服务商的“供给边界”价格不能定得脱离实际供给能力。2.2 下层用户侧能量管理模型下层是用户更准确的说法是“价格接受者”。每个用户在收到服务商发布的电价和热价后会根据自己的用电用热需求、设备特性、舒适度偏好以用能成本最小化为目标进行优化决定各个时段的购电量和购热量以及自家储能设备的充放电策略如果用户侧有光伏、储热罐等也可以一并纳入。一个典型的下层模型包含以下要素可转移电负荷比如电动汽车充电、洗衣机运行可以在一个时间窗内灵活安排但总用电量不变可削减热负荷供暖温度在一定范围内可调允许热功率在一定比例内波动用户舒适度惩罚项约束用户对温度偏离的容忍程度储能设备包括蓄电池和蓄热罐充放能约束、容量约束、爬坡约束。下层用户的优化问题是一个带约束的二次规划或线性规划目标函数一般为购能成本最小加上不适惩罚最小。这个子问题在MATLAB中可以用quadprog求解也可以用yalmip调用cplex或gurobi求解。2.3 博弈均衡与求解思路把上下层放在一起看整个问题是一个典型的双层优化问题。上层命令价格下层根据价格优化用能上层再根据用户的用能行为调整价格如此往复。理论上这个问题可以通过拉格朗日变换、KKT条件转成单层MPEC或MILP来求解但对新手来说这个方法推导复杂、变量爆炸而且求解器配置不当很容易出现数值问题。实际工程项目中更常用的做法是迭代式求解。流程如下设置初始电价和热价把当前价格传递给下层求解各用户的能量管理模型得到用电用热计划把所有用户的负荷汇总代入上层模型计算服务商收益并更新价格判断相邻两轮价格的差值是否小于收敛阈值如果是则停止否则回到第2步。这种做法的优势在于每一层都是独立的、易求解的优化问题MATLAB代码结构清晰调试方便。缺点是迭代次数可能较多、需要调收敛参数而且如果模型是非凸的无法保证收敛到全局最优。但对于动态定价研究来说迭代解出的均衡解已经足够说明问题。3. MATLAB工程实现从环境到代码3.1 环境准备MATLAB版本、YALMIP和求解器这个项目虽然用MATLAB但纯用内置的linprog、quadprog也能做只是代码写起来麻烦模型约束一多就非常容易出错。我强烈建议装上YALMIP工具箱再配一个商业求解器比如CPLEX或者Gurobi。YALMIP的作用是把优化模型用接近数学表达的方式写出来求解器负责真正计算。两者配合起来开发效率能提升一个量级。这里有个容易踩的坑不同MATLAB版本对YALMIP和求解器的兼容性差异很大。比如R2023b之前用旧版YALMIP配Gurobi经常出现“Unable to get licence”或者“Solver not found”的报错。如果你还在用旧版MATLAB不要追求最新版的YALMIP固件包去GitHub上下载与求解器版本匹配的历史release版本更稳。装好之后在MATLAB里运行yalmiptest测试一下能显示Optimizer成功才是真正配好了。另外提醒一句这个项目不需要深度学习工具箱热词里有人搜什么“bilstm代码matlab soc”那是另外的算法方向。这里用到的核心工具箱是Optimization Toolbox再加上YALMIP和外部求解器。3.2 代码整体架构我把整个项目的MATLAB代码拆成五个模块每个模块一个脚本或者函数文件文件/函数作用data_parameter.m输入基础数据24小时电负荷基线、热负荷基线、设备参数、电网分时电价、天然气价格price_init.m初始化电价和热价设置动态定价上下限user_optimize.m下层模型输入价格输出用户电热负荷曲线、储能策略、用户成本operator_optimize.m上层模型输入用户总负荷输出服务商收益、最优设备出力main_stackelberg.m主程序控制迭代流程判断收敛并输出结果这种模块化设计的好处是想改用户数量时只改user_optimize.m想加入新的储能设备时只改对应的约束块不会牵连其他部分。新手最常见的错误就是把所有代码堆在一个脚本里结果改一个参数要滚动半天才能找到受影响的地方。3.3 核心模块实现细节这里给出下层用户优化函数的关键代码骨架方便大家直接参考改造function [P_user, H_user, cost_user] user_optimize(pi_e, pi_h, params) % 输入电价序列 pi_e(24x1)热价序列 pi_h(24x1)结构体 params % 输出用户购电功率 P_user(24x1)购热功率 H_user(24x1)用户总成本 cost_user T 24; % 决策变量购电功率、购热功率、蓄电充电/放电功率、蓄热充/放热功率 P_buy sdpvar(T, 1); H_buy sdpvar(T, 1); P_ch sdpvar(T, 1); P_dis sdpvar(T, 1); H_ch sdpvar(T, 1); H_dis sdpvar(T, 1); % 储能荷电状态 SOC sdpvar(T1, 1); TEC sdpvar(T1, 1); % 储热罐等效温度类比SOC % 目标购能成本 舒适度惩罚 objective sum(pi_e .* P_buy pi_h .* H_buy) params.lambda * sum((H_buy - params.H_base).^2); constraints []; % 电功率平衡 constraints [constraints, P_buy P_dis - P_ch params.P_base]; % 热功率平衡允许一定的供热波动增加舒适度惩罚项 constraints [constraints, H_buy H_dis - H_ch params.H_set]; % 储能约束 for t 1:T constraints [constraints, SOC(t1) SOC(t) params.eta_ch*P_ch(t)/params.cap_elec - P_dis(t)/params.eta_dis/params.cap_elec]; constraints [constraints, 0.1 SOC(t1) 0.9]; constraints [constraints, 0 P_ch(t) params.max_ch]; constraints [constraints, 0 P_dis(t) params.max_dis]; constraints [constraints, P_ch(t)*P_dis(t) 0]; % 充放电互斥可用二进制变量替代 end % 调用求解器 options sdpsettings(solver, cplex, verbose, 0); optimize(constraints, objective, options); % 提取结果 P_user value(P_buy); H_user value(H_buy); cost_user value(objective); end注意这个代码里有几个我故意保留的关键点。第一电池充放电互斥约束用P_ch(t)*P_dis(t) 0这种写法在数值上非常不友好实际项目里应该引入两个二进制变量b_ch和b_dis约束P_ch(t) b_ch*max_ch、P_dis(t) b_dis*max_dis、b_ch b_dis 1。所以需要把模型声明为integer变量并用支持混合整数规划的求解器。第二热负荷目标函数里我加了二次惩罚项这是为了模拟用户对温度舒适度的容忍程度惩罚系数lambda越大热负荷波动越小用户越不愿参与需求响应。上层服务商模型同理约束CHP运行范围时可以使用一个简化的多边形可行域。很多论文里直接把CHP电热关系写成线性函数但在实际系统中这种线性化会遗漏很多运行边界至少要用三点式线性化把可行域近似成多边形。3.4 迭代收敛判定与参数配置迭代求解主从博弈时收敛判据是必须精心设计的。如果只看价格差小于阈值可能出现价格已经稳定但负荷还没稳定的情况如果只看负荷差又可能陷入价格微小振荡。我在实际项目中用的是双判据price_diff max(abs(pi_e_new - pi_e_old)) max(abs(pi_h_new - pi_h_old)); load_diff max(abs(P_load_total_new - P_load_total_old)) / max(P_load_total_old); if price_diff 1e-4 load_diff 1e-3 converged true; end价格阈值取1e-4负荷阈值取1e-3这个组合在大多数场景下都能在10到30次迭代内收敛。如果初始价格设置得太离谱比如电价初值远高于合理区间迭代次数会显著增加还可能出现震荡。建议用上级电网购电价加上一个微小的上浮百分比作为电价初值用天然气折算供热成本作为热价初值这样能让迭代起点落在合理区域内。另外要提一个很多论文里含糊带过的点价格更新时不要直接让新价格等于上层优化结果这样容易振荡。加入一个惯性系数alpha会更稳pi_e_update alpha * pi_e_new (1 - alpha) * pi_e_old; pi_h_update alpha * pi_h_new (1 - alpha) * pi_h_old;alpha取0.6到0.8之间比较合适。alpha太小收敛慢alpha太大容易发散。这是我跑了十几个算例之后得到的经验值。4. 动态定价与能量管理结果分析4.1 仿真场景与基础参数设置测试场景我采用的是典型北方园区一个CHP机组、一台电锅炉、一个蓄电池、一个蓄热罐下游接三个不同类型的用户——居民用户热负荷需求大、电负荷波动大、商业用户电负荷白天高、热负荷平稳、工业用户电热负荷相对恒定但总量大。调度周期24小时时间间隔1小时。上级电网的分时购电价设置峰时10:00-15:0018:00-21:001.2元/kWh平时7:00-10:0015:00-18:0021:00-23:000.8元/kWh谷时23:00-次日7:000.4元/kWh。天然气价格折算后为0.25元/kWh热值。CHP最大电出力500kW热电比固定为1.2电锅炉额定功率200kW效率0.95。4.2 24小时动态定价曲线解读跑完迭代后的动态电价和热价有一个非常明显的规律电价在上级电网购电价的基础上整体上浮但上浮幅度不是固定的而是跟随用户负荷弹性变化。用户弹性大的时段服务商不敢定太高的价怕把用户“吓跑”导致负荷骤降用户弹性小的时段服务商有更大的定价空间。热价的动态特性则不一样。由于CHP“以热定电”或者“以电定热”的运行特性热价曲线往往和电价曲线有一个错峰关系。比如夜间谷时电价低电锅炉更愿意在这种时段产热热价相应下调鼓励用户把热能需求转移到夜间白天峰时电价高CHP发电收益高热价会适度上升。这种联动机制是电热耦合最直观的体现也是动态定价研究的核心价值所在。4.3 能量管理效果对比我把动态定价模式下的结果和固定价格模式做个对比能非常清楚地看出能量管理的效果指标固定价格模式主从博弈动态定价模式变化幅度用户总购能成本元/日2186019640-10.1%服务商总收益元/日4820542612.6%系统峰谷差kW645498-22.8%电负荷峰值kW955892-6.6%用户舒适度惩罚元/日14.218.731.7%从这个表能看出来动态定价实现了服务商和用户的双赢服务商利润提高了约12%用户总成本反而下降了约10%。原因在于动态价格引导用户把一部分电负荷转移到谷时谷时电价低用户省了钱同时系统负荷峰谷差缩小服务商减少了高价购电和机组频繁调节的成本。代价是用户的舒适度损失略微增加体现为惩罚项从14.2涨到18.7但相对于省下的两千多元成本这点不舒适完全可以接受。储能设备在能量管理里也起到了关键作用。蓄电池的策略从“固定时段充放”变成了“追踪电价差”电价低时充电价高时放在价格差超过一定阈值时才动作避免了频繁充放造成的效率损失。蓄热罐则主要配合电锅炉在谷电时段储热在峰电时段放热实现了电热两个系统之间的时间转移。4.4 博弈均衡点的意义主从博弈迭代收敛后得到的解就是斯塔克伯格均衡。在这个均衡点上服务商没有单方面改变价格的动机因为任何价格调整都会导致用户负荷变化最终使得自身收益下降用户也没有单方面改变用能策略的动机因为在当前价格下任何调整都会增加成本。这是动态定价策略具有实操价值的基础。从MATLAB计算的角度判断是否真正达到均衡可以在迭代结束后额外做一次“偏离检测”把最终价格固定住重新求解一次用户模型对比用户成本是否比迭代过程中最后一轮的更小。如果用户成本没有因为模型变化而改变说明此时用户确实已经对价格做了最优响应。这一步虽然简单但很多代码实现的论文里根本没有做结果求解器明明收敛到一个次优解却当成均衡输出。5. 常见问题与排查技巧实录5.1 求解环境问题YALMIP配CPLEX/Gurobi的时候报Solver not found是最常见的问题。排查顺序如下先检查MATLAB路径里是否有YALMIP文件夹和求解器文件夹然后在命令行输入which cplex看能否定位到求解器路径如果定位不到把求解器的bin目录手动添加到MATLAB路径中最后运行yalmiptest检查求解器状态。还有不少人遇到sdpvar在旧版YALMIP下无法定义高维矩阵的问题这通常是因为YALMIP版本太老。建议直接从GitHub下载最新release版本不要用网上流传的几年前的安装包。另外运行前在代码开头加一句yalmip(clear)这行命令会清空之前遗留的YALMIP变量定义避免脚本重复运行buildup问题。5.2 模型求解中的数值问题我在调试时遇到最多的一个报错是NAN in solver input。这几乎都是因为参数里有Inf或者0/0造成的。比如设备效率参数如果从Excel读入时某些格是NaN传到约束里就炸了。建议在加载数据后统一加一步检查if any(isnan(params.cap_elec)) || any(isinf(params.max_ch)) error(输入数据包含无效数值请检查参数表); end还有一个非常隐蔽的问题储能SOC约束写成SOC(t1) SOC(t) ...但如果步长不一致比如从1小时改成15分钟所有率参数都要同步改。我见过一个项目把这个忘了结果储能容量越界形成的“最优解”实际上是不可行的。5.3 迭代收敛与性能调优迭代不收敛是最让人头疼的问题。我遇到过的场景主要有三种第一种价格上下限设置太宽。动态定价通常要给定一个合理的价格走廊比如电价的上下限设在0.3倍和2倍购电价之间。太宽会让迭代在几个极端点之间跳。解决办法是逐步收紧价格走廊比如先跑粗结果再在均衡点附近缩小范围重新迭代。第二种两方模型彼此“过度反应”。这种情况在价格更新系数alpha过大时容易出现上层一次涨价太多下层大幅削减负荷上层又大幅降价形成振荡。通过引入03节说的惯性系数可以稳定下来。第三种下层模型陷入了局部最优。如果用户储能模型包含二进制变量用某些求解器求解MILP时可能会因为初始解不好而停留在次优解。解决方案是给CPLEX/Gurobi设置更长的MIP gap容忍度或者适当放宽整数变量的数量比如把蓄热罐的连续模型当成线性模型处理先用松弛结果判断可行域然后再恢复整数约束求解。性能方面如果用户数量很大几十个用户每轮迭代都要跑几十个子问题时建议把user_optimize函数向量化用MATLAB的parfor并行求解各用户模型整个运行时间能减少一半以上。我在20个用户场景下实测串行跑25轮迭代约需8分钟并行跑只需要3分钟。5.4 代码调试的独家技巧最后分享几个我实际操作中总结的技巧第一单步调试上层和下层模型时不要一次把整个双层模型一起跑。初学者建议先把下层模型单独拎出来固定一组价格验证用户优化结果是否符合直觉比如谷时负荷增加峰时负荷减少。下层验证通过后再跑上层模型验证定价规则是否符合预期。分层验证的速度远远高于直接调试完整的博弈迭代。第二跟踪中间变量。在主循环里把每一轮的电价、热价、用户总成本、服务商总收益打印出来做成一个表格。我通常把迭代过程保存成iter_info.mat每次跑完看这个文件里各变量的轨迹能非常直观地判断收敛情况以及哪一轮出了问题。比如如果第8轮开始价格呈现周期性振荡那一定不是参数设置问题而是模型约束写错了。第三合理使用tic和toc记录每轮耗时。如果某一轮耗时比其他轮长很多说明求解器在该轮花费了大量时间处理整数变量大概率是因为该轮价格和约束边界产生了退化情况此时需要检查该时段的储能约束是不是被激活了。6. 后续扩展与个人体会这个项目的代码框架其实可以往好几个方向扩展。最简单的是把单用户换成多类型用户聚合把24小时扩展到288个时段做多时间尺度分析再复杂一点可以在下层加入光伏和电动汽车研究高比例可再生能源接入对价格策略的影响还可以把单阶段的静态主从博弈改成多阶段的动态博弈考虑设备运行状态对价格策略的约束。就我个人在这个项目里的体会主从博弈模型在MATLAB里的落地真正的难点从来不是数学推导而是工程细节。怎么把约束写得让求解器高效求解怎么设计迭代收敛判据让它既快又稳怎么从一堆输出曲线里提炼出有价值的结论这些才是决定一个代码能不能算出漂亮结果的胜负手。代码本身反而是最简单的一层有了清晰的模块结构任何一步出错都能快速定位。希望这篇内容能帮你少走一些弯路尤其在你第一次跑通上层和下层耦合迭代、看到价格曲线和负荷曲线逐渐收敛到均衡的那一刻会很有成就感。