前阵子我一直在啃一篇关于售电商套餐设计与购电策略的EI论文课题全称是“基于主从博弈的售电商多元零售套餐设计与多级市场购电策略”附带的Matlab代码我前后调了一周多才完全跑通。一开始我以为这只是一个普通的双层优化问题把KKT条件往下一推丢给求解器就完事。真正动手之后才发现模型框架看着清爽实操层面的坑却比公式里的符号还多。这篇内容适合正在复现电力市场类论文的研究生、想理解主从博弈怎么落地的工程师也适合被“EI复现”这几个字劝退的人。我会把模型的建模思路、Matlab实现骨架、参数设计、常见报错和排查经验都摊开讲尽量让你少走几个月的弯路。1. 项目到底在解决什么问题1.1 售电商的两头受气与主从博弈的登场售电商这个角色很有意思它夹在批发市场和零售用户之间两头都得伺候好。购电侧它要在中长期合约、日前市场、实时平衡市场之间分配电量每个市场的价格波动规律不一样买多了或买少了都会产生偏差成本售电侧它要给用户设计零售套餐定高了没人签定低了自己要亏而且用户会根据电价调整用电行为。这两个决策是互相嵌套的你定的套餐价格会影响用户用电量用户用电量反过来又决定你要在市场上买多少电。如果用单层优化建模很容易忽略用户的理性响应最终给出的套餐价格往往是“自嗨”。主从博弈天然契合这个场景售电商先报套餐用户再决定用电量售电商在做决策时已经把用户的反应考虑进去了。你可以把它理解成房东定价、租客决定租几间房的博弈房东不会傻到定一个把租客全部吓跑的价格。1.2 EI论文复现到底在复现什么很多论文把模型写得云里雾里公式推导一段带过参数表给得含含糊糊。所谓“复现”不只是把代码跑出几张图而是要还原完整的建模思路上层是售电商的套餐定价和购电决策下层是用户的效用最大化负荷响应然后用KKT条件把双层问题压成单层最后交给求解器算。这套东西复现下来你会对“主从博弈”有完全不一样的理解。很多人在论文里看到Stackelberg博弈觉得很高深其实落到代码上就是几组约束和变量的组合。这里有一个很现实的问题需要说清楚文献里往往只给核心模型不会给完整参数也不会有调试过程。比如用户效用函数的系数、偏差惩罚成本、价格上限这些关键参数很多论文要么放在附录小字里要么压根不提复现时只能根据算例结果反推或者自己设定这也是为什么同一篇论文不同人复现出的图形会差很多。2. 主从博弈模型的核心框架拆解2.1 上层问题套餐定价与购电优化售电商的决策大致可以拆成两块一块是怎么给用户设计零售套餐另一块是怎么在不同的批发市场里买电。在大多数EI论文里这两个决策是放在同一个优化问题里联立求解的。上层决策变量包括零售电价向量、中长期购电量、日前市场各时段购电量等。目标函数很直观就是最大化售电利润。设用户在某时段的负荷为d_t零售电价为λ_t则售电收入为∑λ_t·d_t购电成本需要考虑中长期合约、日前市场购电和实时平衡电量三个部分。实时平衡电量通常用偏差惩罚来处理即实际负荷与已购电量之差乘以惩罚系数。约束条件里比较关键的是零售电价必须落在用户可接受范围内中长期购电量有上限日前市场购电量不能超过最大申报容量。另外还有一个隐含约束容易被忽略——用户参与约束也就是用户在参与套餐后获得的效用至少要高于他不参与任何套餐、直接用默认电价时的效用否则用户不会签这个套餐。这个约束在数学上表现为一个不等式但它是保证套餐能被市场接受的前提。这里有个常见的误区有人会把套餐设计和购电策略拆成两个独立问题分开求解先定电价再根据电价结果去买电。这种做法在严格意义上是次优的因为购电策略会影响售电成本而售电成本又会影响你能够承受的套餐价格。只有放在同一个模型里联立优化才能体现出“主从博弈”的全局最优性。2.2 下层问题用户效用最大化的需求响应用户作为跟随者面对售电商给定的零售电价会通过调整用电量来最大化自己的净效用。这里的净效用等于用电带来的满意程度减去电费支出。常见的用户效用函数形式为 U(d_t) a_t·d_t - 0.5·b_t·d_t²其中a_t和b_t是与用户偏好相关的系数。用户的决策问题可以写成 max ∑[U(d_t) - λ_t·d_t]对d_t求导并令其等于零就能得到用户的最优用电量响应函数 d_t (a_t - λ_t) / b_t这个响应函数是整个主从博弈建模的关键。它说明用户的负荷并不是固定不变的而是电价的线性函数。a_t可以理解为用户在该时段的用电意愿上限b_t反映用户对价格的敏感程度。b_t越大用户对电价越敏感微小的价格调整就能引起明显的负荷变化。在复杂一点的模型里用户不只有一个而是分成若干类别每一类用户有不同的效用函数参数。这种情况下下层问题就要对每一类用户分别写最优性条件然后全部带入上层问题。还有些模型会考虑负荷的时移特性比如用户可以把部分可转移负荷从高峰时段挪到低谷时段这就需要引入跨时段的耦合约束下层问题的KKT条件会变得更复杂。2.3 从Stackelberg问题到可计算的MILP/MIQP主从博弈的直接求解并不容易因为上下层决策是嵌套的。但如果我们把下层的优化问题替换成它的KKT最优性条件就可以把双层问题转化成单层问题。这个思路在电力市场文献里非常常见也是复现代码的核心。具体来说下层用户的优化问题相对简单它的最优性条件就是前面推导出的电价-负荷响应关系。把这个关系代入上层目标函数之后问题就变成了含有双线性项的单层优化问题。双线性项主要来自目标函数里的λ_t·d_t。好在d_t是λ_t的线性函数代进去之后λ_t·d_t会变成关于λ_t的二次函数因此目标函数仍然是一个二次函数。如果问题里没有整数变量可以用二次规划求解如果涉及套餐类型选择的0-1变量或者阶梯电价的整数分段就变成混合整数二次规划。在实际复现中很多论文用的就是这类MIQP模型求解器直接选Gurobi或者Cplex小算例几秒钟就能出结果。还有一类更复杂的处理方式是用KKT条件互补松弛约束把用户的带约束优化问题完全转写再用大M法把互补松弛条件线性化。这种做法的好处是用户可以带约束条件比如负荷调整上下限模型表达能力更强。代价是引入了额外的二元变量求解规模会变大。我复现的这篇论文采用的是第二种方案因为它的下层用户模型里带有负荷转移约束不能只靠一阶条件简单代换。3. Matlab环境配置与代码骨架3.1 环境准备版本、求解器与工具箱复现这类电力市场优化模型Matlab环境配置其实比想象中重要。我建议直接用R2022b或者更新一点的稳定版本不一定非要追最新版。我自己用R2023a跑Yalmip加Gurobi的体验最稳新出的版本有时候第三方工具箱还没来得及适配反而容易出兼容性问题。安装的时候要注意Optimization Toolbox必装这是Matlab自带优化功能的基础。如果模型是纯MILP/MIQP用自带intlinprog也能跑但性能上比商业求解器差不少。我自己实测下来同样一个算例intlinprog跑两分钟才收敛换Gurobi十几秒就出来了。学术用途的Gurobi可以申请免费licenseCplex也有类似的学术授权配置起来都不算麻烦。如果你只用Yalmip它自带一个编译器检测功能安装完Gurobi后记得在Matlab里运行一下yalmiptest确认求解器能被正确识别。另外并行计算工具箱建议顺手装上后面做多场景对比实验时能省不少时间。3.2 模块化代码结构与核心函数实现整个复现代码我建议按功能拆成几个模块不要在脚本里堆两三百行。我自己的习惯是这样划分数据参数模块、用户响应模块、博弈模型求解模块、购电策略计算模块、结果输出与可视化模块。用户响应模块是最简单的因为下层问题已经可以解析表达。如果你用KKT重构方式需要写一个函数把用户的效用系数、电价向量转成最优负荷向量。这里有一个细节值得注意如果用户有负荷调整上下限约束响应函数就不是简单的线性表达式而是一个带饱和区的分段函数写代码的时候要用min和max把上下界钳住。核心的求解模块用Yalmip建模。Yalmip的语法对这类双层转单层的问题非常友好可以直接用sdpvar定义决策变量用optimize一键求解。我在代码里会把上层决策变量如各时段零售电价、各市场购电量和下层用户的响应变量一起定义出来再把KKT条件转成约束加进去最后设置求解器选项。3.3 求解主循环与均衡结果导出主从博弈模型转成单层之后就不需要迭代求解了直接一次性调用求解器就能得到Stackelberg均衡解。但为了验证算法的正确性通常还会写一个后验模块把求出来的套餐价格重新代回用户的优化问题确认用户的最优响应和你模型里用的负荷一致。我在代码里会输出以下几个关键量最优套餐价格曲线、各类用户的负荷响应曲线、各市场购电量、售电商总利润、用户总福利。这些量在论文图标里基本都是标准输出提前准备好了后面画图也顺手。求解完成后用value()取出变量值再写个结构体保存结果方便后面多个场景之间做对比。% Yalmip求解主从博弈模型核心框架示意 P sdpvar(T,1); % 零售电价上层决策变量 Qday sdpvar(T,1); % 日前市场购电量 Qf sdpvar(1,1); % 中长期合约购电量 D sdpvar(T,1); % 用户负荷经KKT条件与P关联 Constraints [...]; % 含KKT条件、价格上下限、购电容量约束等 Objective ...; % 售电商利润最大化 ops sdpsettings(solver,gurobi,verbose,2); optimize(Constraints, -Objective, ops); P_opt value(P); D_opt value(D); Qday_opt value(Qday);这里我习惯把目标函数取负号因为Yalmip默认是求最小值写成最小时直接optimize(Constraints, Objective)也可以但工程上我更习惯把最大化问题统一换算成最小化避免符号混淆。4. 参数设置、场景设计与你应该跑哪些对比实验4.1 用户、套餐与市场价格参数怎么定参数是复现的重头戏也是论文里最容易被“略过”的部分。我复现时采用了一组典型算例参数24小时一个优化周期用户分三类居民用户、商业用户、工业用户。每类用户的效用函数系数不同工业用户价格敏感度高居民用户价格敏感度低。这个设定比较贴近实际而且三类用户的负荷曲线叠加起来能明显看出峰谷差。套餐设计上我用的是多套零售套餐并存的方式包括固定电价套餐、分时电价套餐和阶梯电价套餐。固定电价适合价格敏感度低、怕波动的用户分时电价鼓励用户在低谷多用电阶梯电价则对高用电量用户起到约束作用。每类用户根据自身效用最大化原则选择套餐这又是一个嵌套在里面的二元选择问题。批发市场侧我设定了一条长期合约价格曲线一条日前市场价格曲线和一条实时平衡价格曲线。日前市场价格在峰时段拉高谷时段压低实时平衡价格在偏差出现时产生惩罚成本。这里的参数直接影响最终结果建议参考相关期刊论文的算例参数去取值别自己随意拍脑袋。比如偏差惩罚系数如果设得太高模型会把所有电量都尽量在中长期和日前市场买齐实时平衡几乎不用结果反而没有参考价值惩罚太低模型又会过度依赖实时市场风险成本被低估。我最后用的是惩罚系数约为日前均价的1.5倍这个数值在不少文献里都能找到依据。4.2 对照实验设计三个必跑的场景判断你的代码和模型有没有实现到位最直接的办法是跑对照组。我建议至少跑三个场景。第一个场景是基准场景售电商只提供单一固定电价套餐购电侧只走日前市场。这个场景等价于传统售电模式指标是后面对比的基线。第二个场景是套餐优化场景售电商提供多元零售套餐但购电策略仍然只走日前市场。这个场景用来单独衡量套餐设计对利润和用户福利的影响。第三个场景是完整场景多元零售套餐加多级市场购电策略全部打开。这个场景对应论文的完整模型也是最终要呈现的结果。把三个场景的利润、负荷曲线、峰谷差、用户福利放在一张表里你就能清楚看到每一层优化分别贡献了多少收益。我在自己的复现结果里场景三相比场景一利润提升了大概18%其中套餐优化的贡献占了大头购电策略优化提供了进一步的改善这个量级也和文献报道一致。如果你的结果提升幅度过大或者过小就要回头检查参数和约束条件是否合理。4.3 如何判断你的结果是对的这里说一个最容易被忽视的问题模型跑通了结果也出来了但你怎么确定这个结果是“对”的光看目标函数值没有意义需要做几个验证。首先是均衡验证。把求出来的套餐价格固定住单独求解用户的最优负荷响应看和模型里用户负荷是否一致。如果两个结果不一致说明KKT条件转写或互补松弛处理有问题。其次是单边偏离检验。在最优套餐价格基础上给某个时段的价格加一个小扰动重新求解用户响应和售电商利润。如果扰动后的利润比原结果低说明原结果是局部最优的候选解如果反而更高说明模型或求解器出了问题。这个检验简单又有效我每次跑完新算例都会做一遍。第三是补松弛校验。对于用到互补松弛条件的场景检查一下乘积项是否在容差范围内趋近于零。如果残差很大通常是大M参数取得不合适需要调大或者改换其他的线性化方案。5. 复现过程中的坑与排查手册5.1 模型层的坑大M取值、双线性项与不收敛先说大M取值。互补松弛条件线性化时的大M参数非常敏感M取得太大数值计算会出现病态求解器收敛慢甚至给出错误的最优解M取得太小又可能把可行域截掉导致结果偏离真实最优解。我自己的经验是按照问题物理边界来估算M值比如电价上限乘以负荷上限再乘一个裕度系数算出来多少就填多少不要偷懒直接填一个很大的数。双线性项处理是另一个重灾区。如果模型里含有连续变量乘积比如购电量和实时价格相乘求解器会直接报非凸或者不收敛。处理办法通常有三种一是利用下层响应函数消元把双线性项变成单变量二次项二是用McCormick包络做松弛但会有松弛误差三是引入辅助变量配合大M法做精确线性化。具体用哪种取决于你的模型结构。优先尝试第一种因为它最干净。还有一类不收敛问题来自目标函数数级差异过大。比如售电收入是百万级别而惩罚成本在千级别求解器在数值上会忽略小量级项。出现这种情况我一般会对目标函数的各项做归一化处理或者给各项加上合理的权重系数。别小看这个操作很多时候模型在理论上没问题跑起来结果离谱就是因为数值尺度不一致。5.2 软件层的坑license、版本与Yalmip配置Matlab软件层面我遇到的坑也不少。最典型的是license问题。这类优化模型经常要在实验室服务器上跑不少人用远程桌面连服务器时发现Matlab打不开或者启动时报mathworks licensing error 9。这个错误通常是license文件和当前机器的hostid绑定不一致导致的。解决办法是检查当前机器的MAC地址是否和license文件里记录的一致如果不一致需要重新激活或者在Matlab启动脚本里显式指定环境变量MLM_LICENSE_FILE指向正确的license文件。关于版本密钥我的建议是不要用来路不明的所谓密钥轻则激活失败重则被官方拉黑。学校有校园授权就用校园版没有就申请官方试用版完全够用。另外重装Matlab时如果遇到“删除不干净”的问题记得把环境变量和用户目录下的MathWorks残留配置一起清理掉否则新的license激活会被旧配置干扰。Yalmip版本和求解器版本不匹配也很常见。旧版Yalmip可能不认识新版本Gurobi的接口导致yalmiptest通不过或者求解时直接报错。遇到这类问题升级Yalmip到最新版通常能解决或者去Gurobi官网下载对应Matlab接口文件手动配置。我还在代码里遇到过Nonconvex quadratic的警告这一般是模型里出现了未线性化的双线性项需要回到5.1的三种处理办法里排查。5.3 数据与可视化数组操作、循环画图与图件导出后处理阶段的坑虽然没有那么致命但真的很浪费时间。第一是数据导入导出。如果你的电价数据存在Excel里读进来之后时间列常常是datetime类型直接拼接字符串做横轴标签会报错建议先datestr()转格式或者用string()转换再用datetime统一管理。第二是数组索引的坑。处理24小时数据时经常要提取特定时段比如峰时段8到11点用s(:, 8:11)这种列取法很方便但要注意行和列的顺序我经常因为索引方向搞反导致画出来的曲线完全对不上。建议在处理之前先用size()确认维度或者用reshape把所有数据统一成列向量能省不少事。第三是画图循环。多场景对比图比如基准场景和优化场景的负荷曲线画在同一张图里要用循环统一设置颜色和线型别手写三遍plot再手敲三个legend。legend可以用cell数组动态生成避免每改一个场景都要改图例。我自己习惯用legend({基线负荷,优化后负荷}, Location,best)这种写法配合set(gca,FontName,Times New Roman,FontSize,11)统一字体导出图片时用exportgraphics(gcf,xxx.png,Resolution,300)导师要的清晰度基本都能满足。下面是我遇到的高频问题速查表按严重程度排了一下现象可能原因解决建议求解器报不能处理二次约束模型含有非凸双线性项优先消元其次考虑分段线性化目标值偏离常识很远目标函数各项数量级差异过大对目标各项做量纲归一化用户负荷响应与KKT代换结果不一致用户约束漏写或大M法参数不当检查下层约束调整M值图例位置重叠或文字过小直接无脑用默认设置用legend指定位置统一FontSizeMatlab远程桌面打不开/license报错license与hostid不匹配检查MAC地址重设环境变量6. 后续还能往哪些方向扩展6.1 从单售电商到多售电商竞争我复现的模型是单一售电商作为领导者、多个用户作为跟随者的结构。现实中一个区域往往有多家售电商在竞争用户可以选择签约其中任意一家。这种情况下的均衡就是多个领导者之间的博弈属于均衡约束均衡问题EPEC求解难度比单层Stackelberg高出很多。如果你想在这个方向做工作一个常见的做法是把多售电商博弈处理成迭代过程每一轮固定其他售电商的策略求解单个售电商的主从博弈然后循环更新直到收敛。但这种迭代方式不保证收敛到唯一均衡对初值敏感需要配合小步长更新或者松弛技巧。复现完基础模型之后再往这个方向扩展你会对博弈模型的适用范围有更清楚的认识。6.2 不确定性、动态与数据驱动方向目前这个模型假设市场价格和用户参数是确定性已知的。实际上日前市场价格、实时平衡价格都有很强的不确定性用户负荷也存在随机波动。把不确定性纳入模型可以考虑鲁棒优化、分布鲁棒优化或者随机规划。这些扩展会让模型从MIQP变成更复杂的结构求解难度成倍增加。如果你对数据驱动感兴趣可以把用户历史负荷数据拿来做聚类用聚类结果直接标定不同类型用户的效用函数参数而不是像我前面那样手工设定a_t和b_t。Matlab里聚类工具箱可以直接上手把负荷曲线聚类之后每一类的响应系数可以通回归估计出来。另外如果要对动态定价过程做仿真拿离散时间状态方程来描述用户负荷变化Matlab里ode45或者ss这类工具也都能直接配合优化模型使用。6.3 我复现完之后最想说的一句话整个项目做下来我最深的感受是复现EI论文不是对着公式敲代码而是在还原作者每一步建模思考。很多关键的约束条件、参数取值在论文里可能只是半行符号代码实现时却决定了整个模型能不能算出符合直觉的结果。你如果能坚持把每个约束为什么存在、每个参数为什么取这个量级都搞清楚这篇论文就算“吃透”了。另外一个小建议代码注释要写详细特别是每个约束对应的论文公式编号。我的习惯是在每条Yalmip约束后面加一行注释标明它来自论文的第几个公式或者哪一段描述。等过两个月再回头调试或者改参数的时候你会感谢当时的自己。这次的分享就到这里希望对正在做或者准备做售电商博弈模型复现的同学有点参考价值。