做科研的人拿到一个EI复现代码最尴尬的是啥不是代码跑不通而是代码能跑通但完全看不懂它在干什么。尤其是“主从博弈”这种名字一听就头大的优化模型如果只是把代码当作黑盒跑一遍出几张图那这论文基本算白复现了。这篇文章我就把“基于主从博弈的售电商多元零售套餐设计与多级市场购电策略”这个项目彻底拆开从模型结构到每一段代码背后的数学逻辑再到调参踩坑全部过一遍。先给还没接触过这个方向的朋友说清楚这个项目解决的是售电公司在电力市场环境下的两个核心决策问题——怎么给终端用户设计零售套餐电价、电量、时段差异以及怎么在上游批发市场中长期、日前、实时合理分配购电比例让公司收益最大、风险可控。这两个问题不是独立的用户会根据你定的套餐调整用电你的套餐又反过来影响你在市场上的购电成本和偏差考核费用所以要用主从博弈把两边同时建模求解。你如果正在做电力市场、售电公司决策、需求响应、Stackelberg博弈相关的课题或者需要一份能直接跑通出图的Matlab毕业设计代码那这篇文章正好合适。1. 主从博弈的整体建模思路谁先动谁后应1.1 售电商和用户其实是一个“先出牌”和“后接招”的关系主从博弈的学术名字叫Stackelberg博弈它的核心在于“领导者先行动跟随者根据领导者的行动作出最优反应领导者再预测跟随者的反应来优化自己的策略”。放到这个项目里售电商是领导者用户是跟随者。先捋一下决策顺序第一步售电商制定零售套餐套餐里包含固定电价、分时电价、阶梯折扣等条目。第二步用户看到套餐之后根据自己的用电习惯和用电舒适度偏好调整各时段用电量目标是让自己效用最大化也就是“用得舒服”和“花得少”之间的平衡。第三步售电商知道用户会对套餐作出这种反应于是回头优化套餐参数使得扣除购电成本之后自己的净收益最大。看到没有这就是一个典型的“你先定规则、我按规则行动、你再根据我的行动改规则”的循环。这个逻辑用数学写出来就是一个双层优化问题Bilevel Optimization Problem上层是售电商的利润最大化下层是用户的效用最大化。代码里如果只用单层优化那等于把用户的反应忽略了算出来的“最优套餐”在真实场景里根本站不住脚。1.2 双层模型具体在优化什么目标函数和关键变量这个项目的数学模型我直接把它拆开讲比看论文里的公式舒服得多。上层问题售电商视角售电商的收益 卖电收入 - 批发市场购电成本 - 偏差考核费用 - 套餐运营成本其中卖电收入就是所有用户在各时段按照零售套餐支付的电费总和购电成本就是在中长期市场、日前市场、实时市场买入电量的加权成本。偏差考核费用更关键因为用户实际用电和你当初在市场上申报的购电量不可能一模一样偏差大了电网会收惩罚性费用所以购电策略必须考虑用户响应的不确定性。下层问题用户视角用户的目标函数是效用最大化效用 用电带来的满意感 - 电费支出。满意感通常用二次效用函数表示因为边际效用递减——电用得越多每多一度电带来的幸福感提升越小。用户通过调整各时段用电量来最大化效用而用电量就是这个下层问题的决策变量。上下层之间的耦合变量就是售电商定的零售电价和用户反馈的用电量。电价一变用户用电量就变用电量一变售电商收入就变购电偏差也就变了。这就是博弈的“耦合”所在。模型必须同时求解上下层而不是分开算。1.3 为什么必须用博弈而不是单层优化很多人第一次做这个题会问我直接让售电商的利润函数里面带一个用户用电量的表达式然后一起优化不行吗理论上行但实践中有两个问题。第一个问题是用户用电量并不仅仅是电价的简单函数它由用户自己的优化问题决定。用户是“理性主体”他会自己衡量用电效用和电费支出你直接把他的用电行为硬编码成电价的一次函数模型就失真了。主从博弈的优势在于用下层优化问题来内生刻画用户行为不需要人工假设所谓的“需求价格弹性曲线”。第二个问题是如果售电商直接决定用户用电量那就意味着售电商拥有对用户用电的完全控制权这在真实市场里不成立。用户有选择权有不配合的自由。主从博弈尊重跟随者的决策空间让双方都在自己的约束下做到最优这才是市场机制的本质。所以这个项目如果用单一优化模型去复现审稿人一眼就能挑出毛病需求响应被你建模成了指令控制不是市场激励。这就是为什么顶刊论文普遍采用主从博弈框架的原因它把“激励机制设计”这个核心问题嵌进了数学模型里。2. 多元零售套餐设计不是拍脑袋定价格是求解出来的2.1 套餐的四个核心维度电价水平、峰谷划分、时段比例、优惠策略项目标题里强调“多元零售套餐”这四个字不是营销话术。在实际建模里套餐的可变参数至少包含以下四个维度。电价水平是基准典型做法是阶梯电价用得越多单价越高或者反过来用得越多折扣越大取决于售电商想让用户多用还是少用。峰谷划分是时间维度把一天切成峰段、平段、谷段分别定价引导用户把用电从高峰挪到低谷。时段比例是各时段的持续小时数比如峰段从8点到12点还是从8点到11点会直接影响用户的转移空间。优惠策略包括固定折扣率、用电量达到某个门槛后打折、或者套餐绑定某种负荷特性的定向优惠。这些参数在代码里通常被定义为一个结构体数组或者一个矩阵每一行对应一种套餐方案每一列对应一个参数维度。初始化的时候参数取值范围的合理性直接决定求解器能不能搜到可行解。我见过不少复现代码把电价上限设成实际市场价的几十倍求解器一跑全是Inf就是因为没做参数合理性检查。2.2 用户的效用函数怎么定二次函数里的系数有什么讲究用户侧建模是主从博弈里最容易出问题的地方。最常见的写法是用户效用 a * 用电量 - b * 用电量^2 - 电费支出其中a和b是效用偏好系数a越大代表用户越喜欢用电比如夏天家里空调多b越大代表边际效用递减越快用够了就真的够了。这两个系数不是随便填的它们决定了用户对价格变化的敏感程度。代码里如果只是设个a1、b0.1然后祈祷模型收敛那你大概率会得到一组荒诞结果——用户用电量在电价暴涨时几乎没变化因为a太大b太小用户“离不开电”。正确做法是根据实际负荷数据反推。比如某用户平均用电量是10度平均电价是0.5元那么用户的边际效用应该在电价附近波动也就是说2a用电量大概是0.5这个量级。用这个关系反推a和b的数量级自然就合理了。2.3 套餐参数的求解细节KKT条件与强对偶转换下层用户问题如果是个凸优化问题二次效用函数对应的目标函数通常是凹函数最大化凹函数等价于最小化凸函数所以是凸的那就可以用KKT条件把下层问题转成上层问题的约束条件从而把双层问题变成单层问题。具体来说对用户的目标函数求导令其等于零再补上用户用电量的上下限约束和对应的互补松弛条件原来的下层优化问题就被替换成了一组方程和不等式。这样上层售电商的问题就变成了一个带非线性约束的单层优化问题求解器可以直接处理。但这里有个坑下层问题的约束如果不是线性或者不是凸的KKT条件就不是充要条件转化之后可能丢失最优解。用户用电量如果是连续变量、目标函数是二次可导凹函数那没问题。但如果你给用户加了0-1变量比如某些电器开/关下层问题变成混合整数规划KKT转换就会失效。这个项目里如果用户只调整连续用电量那KKT是安全的如果论文里用到了离散设备那就得考虑用强对偶条件或者直接枚举用户策略空间了。看到代码里出现KKT相关约束时务必确认下层问题的凸性假设成立否则整个模型的理论基础就是错的。3. 多级市场购电策略中长期、日前、实时怎么分配3.1 三级市场的角色分工为什么不能只在一个市场买电售电商购电不是一次性把事情办完的。电力市场的层级结构决定了购电必须分阶段进行而且每个市场的功能完全不同。中长期市场是“保基本盘”的地方提前几个月甚至一年签订电量合同价格锁定但灵活性差。日前市场是“细化调整”的地方提前一天根据负荷预测确定次日各时段购电量价格比中长期贵一点但更贴近真实需求。实时市场是“查漏补缺”的地方实际运行前几小时甚至实时买卖价格波动大、风险高只在偏差出现时用。一个合理的购电策略是中长期锁定基荷的60%到70%日前市场买掉剩余需求的20%到30%实时市场只处理5%以内的偏差。这个比例不是固定的需要根据市场风险偏好和预测精度动态优化。代码里如果这三个市场的购电量都是决策变量那么约束条件要特别留意——实时市场的购电量通常可以买也可以卖即允许正负但中长期和日前市场一般只允许正向电量。3.2 购电决策变量怎么定义两个阶段之间的传递逻辑购电策略的核心决策变量是各市场各时段的购电量。结构上通常这样定义中长期购电矩阵维度是[时段数量]因为中长期合约基本都是分时段电量曲线日前购电矩阵同样是[时段数量]实时购电矩阵维度略有不同因为实时市场可以15分钟甚至5分钟为一个出清周期。代码实现里最容易被忽略的是“电量传递关系”用户最终用电量 中长期购电量 日前购电量 实时购电量这个等式不是自动成立的需要作为约束条件显式写进模型。很多复现代码跑出来结果很奇怪购电总量和售电总量根本对不上一看就是少了这个平衡约束。还有个细节是偏差考核。实际运行中售电商申报的总购电量和用户实际用电量之间会有偏差偏差超过一定阈值就要交考核费用。这个费用在目标函数里通常表示为偏差电量乘以惩罚系数的二次函数或线性函数。惩罚系数设多大很影响结果——设小了售电商会懒得做精细化预测全指望实时市场兜底设大了售电商会过度保守中长期买得太多失去市场灵活性。一般做法是让惩罚成本略高于日前市场最高电价逼着模型平衡预测和偏差。3.3 不确定性怎么处理场景法还是鲁棒优化用户用电行为有随机性市场电价也有随机性购电策略如果只基于一个确定的预测值那跟赌徒没什么区别。代码里常见两种处理方式。场景法最直观生成一批典型场景比如用蒙特卡洛从历史数据里抽样得到100组负荷和电价场景模型的目标函数改成“所有场景下的期望收益最大化”。好处是模型直观、容易实现坏处是计算量大而且场景生成的质量直接决定结果质量。如果只生成10个场景求解结果会很“碎”生成100个以上求解时间又可能受不了。鲁棒优化则是考虑“最坏情况下的最优策略”也就是不管市场怎么波动保证在最差情形下收益不低于某个底线。这种模型求解效率高但结果偏保守可能在普通场景下损失不少利润。这个项目如果想做“高级一点”的复现可以在原来确定性模型基础上加上场景生成模块。先跑确定性版本验证代码正确性再上场景版本做结果对比这样论文里能多出一张“确定性vs随机性购电策略对比图”工作量不大但效果很明显。4. Matlab代码实现架构从数学公式到可运行代码的路线图4.1 为什么用MatlabYALMIP这个组合在电力市场里的地位电力系统优化领域Matlab配YALMIP再加一个商业求解器CPLEX或Gurobi几乎就是事实标准。YALMIP的作用是把数学建模语言翻译成求解器能懂的格式你不用手写单纯形法或者内点法只需要把变量、目标函数、约束条件按照语法写清楚就行。为什么不用Python不是不能用Python的Pyomo也能做类似工作但对于很多IEEE/EI论文里的复现代码Matlab版本更常见因为老一辈研究者用Matlab用了几十年代码生态已经沉淀下来了。你做EI论文复现直接用原论文的Matlab环境是最省事的。YALMIP的Symbolic建模方式非常接近论文里的数学公式写法调试起来很直观这是其他工具很难替代的。4.2 代码结构怎么组织六个模块的职责划分一份能看的复现代码绝对不是把所有代码堆在一个脚本里。这个项目建议按下述六个模块拆分。主程序模块负责定义问题规模用户数量、时段数量、市场阶段数量、调用各子模块、汇总结果并出图。参数初始化模块集中管理所有可调参数——用户效用系数、市场电价曲线、惩罚系数、套餐参数上下限等方便基准测试时统一调参。零售套餐建模模块负责生成套餐候选方案或者套餐参数的决策变量结构。用户响应求解模块用于求解下层用户最优用电问题通常是调用YALMIP直接求解。购电优化模块承载上层主问题包括购电平衡约束、偏差考核约束和目标函数。结果处理与可视化模块把优化出来的套餐价格、用户用电量、购电分配比例、各主体收益画成图。每个模块建议写成函数不要写成脚本。函数的好处是接口清晰、可复用你可以只改参数传入值就能跑不同场景而不需要每次复制粘贴大段代码。我见过一些复现代码通篇是脚本数据存在一堆变量里换个用户数量就得改好几处这种代码看着跑得通实际上扩展性为零。4.3 KKT转化在代码里具体长什么样一个最小示例这个项目最关键、也最容易被写错的部分就是把下层用户问题用KKT条件转成上层约束。这里我给出一个极简框架展示代码长什么样。假设只有一个用户t个时段变量是各时段用电量P_load(t)零售电价是price(t)。用户问题为% 用户侧目标最大化效用即最小化负效用 User_obj sum(0.5 * a_ut * P_load.^2 - b_ut * P_load) sum(price .* P_load); % 上面的形式是“效用 - 电费”取负所以在求解器里是“最小化”表达式写成KKT条件后在主问题里的样子是% KKT 条件替代用户问题 % stationarity: 一阶导为零 Stationarity_con [a_ut .* P_load - b_ut price - lambda_lo lambda_up 0]; % 其中 lambda_lo 是下限约束的对偶变量lambda_up 是上限约束的对偶变量 % complementarity: 松弛互补条件 Comp_lo [lambda_lo 0, P_load P_min, lambda_lo .* (P_load - P_min) 0]; Comp_up [lambda_up 0, P_load P_max, lambda_up .* (P_max - P_load) 0];注意这里的互补松弛条件涉及“变量×变量等于零”这种非线性等式YALMIP默认是不喜欢这种约束的因为非线性等式容易导致求解困难。实际处理时常用大M法把互补条件线性化% 引入0-1辅助变量 z_lo 和 z_up把互补条件拆成线性不等式 M 1e4; % 足够大的常数 Linearized_comp_lo [lambda_lo M * z_lo, P_load - P_min M * (1 - z_lo)]; Linearized_comp_up [lambda_up M * z_up, P_max - P_load M * (1 - z_up)];这就把原问题变成了一个混合整数线性/二次规划CPLEX这类求解器处理起来就稳定多了。4.4 如果没有CPLEX/Gurobi许可证怎么办替代方案实测对比商业求解器许可证确实是个拦路虎。Matlab自带的是linprog和intlinprog但处理混合整数非线性问题MINLP时表现不太行。我的实测经验是如果问题规模不大用户数不超过5个、时段24个用Matlab自带的ga遗传算法配合内层YALMIP求解也能跑出结果但稳定性和速度都不如CPLEX。另一个可行路径是SCIP它对外提供免费的学术许可证YALMIP可以直接对接SCIP。实测下来对于这个项目的模型规模SCIP的求解时间比CPLEX慢大约3到5倍但至少能稳定跑出结果。做毕设或者课程设计的话SCIP完全够用。如果论文要求严格的求解效率对比实验再去解决CPLEX/Gurobi的许可证问题。5. 实操过程与核心环节实现从空文件夹到出结果全流程5.1 环境安装与验证别在第一步就卡住这步是劝退最多人的环节没有之一。建议按下面这个顺序一次性搞定。先安装Matlab然后确认版本至少要R2020a以上因为老版本对YALMIP某些语法支持有问题。然后下载YALMIP把整个文件夹放到Matlab的toolbox目录下在Matlab里执行addpath(genpath(你的YALMIP路径)); savepath;验证安装的关键命令是yalmiptest这个命令会跑一系列测试用例如果绝大部分显示成功那YALMIP就OK了。别只看“No problems found”这种提示重点看“sdpvar”“milp”这些关键测试项有没有通过。然后安装求解器。如果临时没有CPLEX先用Matlab自带的linprog凑合跑线性版本但遇到整数变量时大概率会报错不用慌那是因为求解器不支持整数规划不是你的代码错了。环境配置完后强烈建议先跑一个小的测试算例。比如设置24时段、3个用户、无中长期合约只测日前和实时购电。这个规模跑出来非常快你可以确认“模型能不能解”而不是“模型解得好不好”。等基础算例能跑通再逐步加复杂约束和更多用户。5.2 参数初始化的经验数值第一版就能收敛的配置这里直接给一组我反复验证过的初始参数照着设基本不会翻车。用户效用系数方面a_ut取值在[0.8, 1.2]区间b_ut取值在[0.01, 0.05]区间这样用户对电价的敏感度比较接近真实水平。零售电价初始值方面谷段电价在市场价的0.8倍左右平段约1.0倍峰段约1.3倍。购电比例初始值方面中长期电量占比60%日前30%实时10%偏差惩罚系数设为日前市场电价的1.5倍。时段划分方面峰段8小时9点到17点平段8小时谷段8小时。这里最关键的原则是所有初始值必须让模型中所有约束在可行域内部而不是刚好在边界上。比如P_load的初始值如果取的是上下限平均值那第一轮迭代每个时刻都是松弛的求解器更容易找到方向。如果初始值给定在边界上很多求解器会因为活动约束太多而夭折。5.3 运行结果怎么看四个关键指标是否合理出图之后别急着高兴先用四个指标验证结果合理性。用户在峰段的用电量是否低于平段如果峰段电价明显高于平段而用户峰段用电量还更高基本可以确定下层模型配置有问题。售电商总购电量和总用电量的平衡误差是否趋近于零如果购电量比售电量多出一大截说明平衡约束写错了。中长期购电量在总购电量中的占比是否处于50%到70%区间如果中长期占比超过90%说明现货价格波动惩罚设得太高售电商过分保守了。售电商最终利润是否为正且量级与售电收入匹配如果利润为负但收入很大说明购电成本或者惩罚项的计算有符号错误。这四个指标是代码是否复现正确的“体检报告”比盯着迭代次数收敛曲线管用多了。结果合理之后再去调整套餐参数作敏感性分析这时候才进入正常的“科研环节”。6. 常见问题与排查技巧实录六个高频大坑6.1 求解器报“Infeasible problem”或“No solution found”这是复现这个项目遇到最多的问题。排查思路按下面顺序来。先查约束是否自相矛盾。这个项目最常见的矛盾是P_load的上下限跨度过小同时购电平衡约束严格要求和售电量相等结果任何一组解都无法同时满足。松弛方法是把平衡约束的等式改成等式加容忍度比如左侧允许在[-1e-4, 1e-4]之间波动。再查有没有忘记定义某个变量。YALMIP的报错有时候不会精确指出未定义变量只会笼统说Infeasible。检查代码里每个出现在约束中的变量是否都在前面用sdpvar定义了。还要查大M常数是否过大。我习惯从1e3开始如果不可行就逐级检查1e4、1e5。M值太大可能导致数值不稳定太小则可能导致KKT条件不成立。6.2 YALMIP报“Cannot use the solver”或者“No suitable solver”遇到这种报错先检查是否装了求解器以及YALMIP是否识别到了求解器路径。在Matlab里执行yalmiptest看CPLEX或者Gurobi是否出现在可用列表里。如果没出现通常是solver的bin文件夹没有添加到Matlab路径。另外模型本身如果因为非线性等式比如互补松弛条件没做线性化被YALMIP判定成MINLP问题而你没有安装MINLP求解器比如Bonmin也会报这个错。处理方法就是把所有互补条件都用大M法转成线性约束不要让YALMIP识别成MINLP。6.3 求解时间爆炸从10分钟到跑了2小时还没出结果这个问题在用户数超过10个以后非常常见。两层优化如果用了复杂的KKT条件非线性项会显著拖慢求解速度。我的建议是分级降规模。先把24时段聚合成6个特征时段、算试算再把每个售电商对应的用户群进行聚合比如按负荷特性分成3类每类代表一群用户然后固定部分决策变量比如先用启发式算法固定购电分配比例只优化零售套餐确认没问题之后再两阶段联合优化。还有一招是冷启动先用上次保存的最优解作为初始点传给求解器YALMIP支持assign和guess等初始解赋值函数能大幅缩短求解时间。6.4 结果图“曲线满天飞”但看不出规律画图前先做统计复现代码跑完直接plot一堆曲线经常是一团乱麻。我习惯的做法是先对结果做统计分析再画图展示关键规律。比如把所有场景下的购电比例均值算出来做一个饼图把各时段套餐电价和用户用电量的变化曲线画在同一张图里双纵轴展示把不同套餐方案下用户总用电量和售电商利润的散点图做出来加趋势线。这样出来的图不仅自己看得懂放进论文里也立得住。6.5 用户用电量出现负值或者超出物理极限约束根本没起作用这几乎可以肯定是约束定义范围有误。检查P_load的上下限是否在初始化时写错了顺序比如把最大值写成了最小值。另外检查KKT转化时是否漏写了上界或下界对应的互补约束。最简单的方法是先去掉KKT转化直接用fmincon解用户问题做对照看看用电量是否落在合理区间内。如果fmincon的结果正常而YALMIP的主问题结果异常那就是互补条件的线性化出了问题。6.6 代码换了一台电脑就跑不出同样结果数值稳定性问题不同电脑的浮点运算结果会有细微差异但对于优化问题来说细微差异可能导致求解路径完全不同。这就是为什么保存随机种子、保存求解器选项、固定初值如此重要。代码开头建议加上rng(2024); % 固定随机种子并且把求解器的容差参数在代码中显式设置好比如ops sdpsettings(solver, cplex, verbose, 2, cplex.mip.tolerances.mipgap, 0.001);这样换机器跑只要求解器版本一致结果基本可以复现。7. 这组代码还能往哪些方向扩展这个项目做完之后如果精力允许有几个方向可以进一步扩展对发论文或者做毕设加分会很明显。第一个方向是考虑储能设备。售电商本身如果配置了储能那就可以在电价低谷充电、高峰放电套利的同时还能平衡偏差多一级决策会让模型更好看代码里只需要新增储能容量约束和充放电状态变量。第二个方向是多售电商竞争。目前只有一个售电商作领导者如果扩容成两个售电商同时定价用户可以在两者之间选择博弈就从单领导者变成多领导者——这已经是博弈论前沿课题了做出来比较容易出成果。第三个方向是把需求响应机制做得更落地比如加入可中断负荷、电动汽车充电调度、温控负荷聚合让用户侧的响应行为更有物理解释而不是纯靠效用函数驱动。我个人在实际复现这个项目时最深的体会是主从博弈模型最难的地方从来不是数学推导而是让上下层之间的变量传递不出错、让KKT条件和求解器稳定配合。YALMIP给了我们一个非常方便的声明式建模环境但查错的时候还是得回到物理直觉——电价高了用户就该少用电购电比例就该追求降低偏差成本如果结果违反了这种直觉那模型多半有bug而不是现实世界出了问题。最后再分享一个小技巧当你自己完全看不懂某个约束为什么会带来这个结果时把这个约束单独拿出来测试一下找一个最简场景单独求解通常立刻就能定位问题所在。别在大模型里反复试那是性价比最低的调试方式。