做电力市场方向的复现，最花时间的往往不是看懂公式，而是把主从博弈这种双层决策结构改写成能求解的Matlab代码。这篇基于主从博弈的售电商多元零售套餐设计与多级市场购电策略，正好是这类模型里比较典型的一篇：上层售电商定套餐价，下层用户根据价格改变用电行为，中间还要跨中长期、日前、实时多个市场做购电分配。很多人一拿到这类EI论文就开始找代码，结果卡在模型转化上。这篇博客我不讲虚的，直接按我复现时的思路，把问题从业务逻辑拆到数学公式，再落到Yalmip代码，最后把常见报错和处理方法一并列出来。
适合正在做电力零售套餐、需求响应或者购电策略方向研究的人，也适合刚接触双层优化与MPEC/MILP复现的初学者。你不需要提前懂太多电力市场机制，但最好有最基础的优化理论和Matlab操作经验。读完这篇，你能搞清楚这类论文的核心建模主线，也能照着这个思路把一个双层博弈模型从公式推到可运行代码。
1. 复现前必须先想清楚：这个模型在解决售电公司的什么问题
复现这类论文，第一步不是打开代码，而是把业务逻辑讲成数学语言。业务逻辑不清，后面的公式和代码都是空中楼架。
1.1 一桩生意拆成两层决策
售电公司的生意其实像一家“买进卖出”的批发商：左手在批发市场采购电力，右手把电包装成不同套餐卖给终端用户。批发侧，购电不是一次性行为，而是跨多个市场的组合决策——中长期合约锁定基础电量，现货市场用于调整偏差，实时市场兜底。零售侧，用户不是被动接受价格，而是会根据套餐价格调节自己的用电行为。这两侧决策互相影响：套餐定得贵，用户少用电，购电量就得跟着调；购电成本高了，套餐价格又不得不涨。
这就是为什么论文要把问题写成双层结构。上层是售电公司的利润最大化问题，决策变量是套餐价格和各时段购电量；下层是所有用户各自独立地最大化自己的净效用，决策变量是选择哪个套餐、用多少电。上层做决策时必须预见下层会怎么反应，而不是拍脑袋定价。这种“我先走、你跟着最优反应”的博弈结构，正是主从博弈。
初看起来，双层规划只是“两个优化问题叠在一起”，实际建模时远比这麻烦。下层用户的优化目标在上层模型里不是约束条件，而是一个“必须是最优解”的隐式条件。你不能直接要求用户用电量等于某个值，只能通过KKT条件或者强对偶关系把它转化为等价约束。这个转化是整个复现工作的核心，也是大多数代码出错的根源。
1.2 主从博弈的“主”和“从”分别是谁
主从博弈中，先动的一方叫领导者，后动的一方叫跟随者。在这篇论文的语境里，领导者是售电公司，跟随者是各类用户，进一步还可以把多级市场的出清机制看成中间层的市场规则。
具体来说：售电公司先公布零售套餐的价格、时段划分和约束条件；用户看到套餐后，计算哪种套餐对自己最划算，决定签约和用电量；售电公司再把用户的总需求拿去市场上采购。因为售电公司了解用户的效用函数和选择逻辑，所以它可以提前“模拟”用户的下层问题，把用户的反应写进自己的优化模型。这正是主从博弈最吸引人的地方：它把“用户会怎样反应”这个看似不可控的因素，变成了一种可以计算的最优响应。
复现时有一个容易误解的点：用户是不是具有完全理性？论文通常假设用户是价格响应型，追求自身效用最大化，这个假设在小用户场景里成立，对大型工业用户未必完全成立。但作为复现，我们跟着论文的假设走就可以了。不要试图在复现阶段引入行为经济学模型，那会让求解复杂度直接失控。
1.3 多元套餐与多级市场购电怎么耦合
多元零售套餐不是简单的“峰谷电价”三个字。常见的套餐结构包括：固定电价套餐、分时电价套餐、阶梯折扣套餐、套餐电量上下限、违约金条款等。每种套餐对应一组价格参数和约束条件，用户在不同套餐之间做0-1选择，在同一套餐内部决定逐时段用电量。
购电侧的多级市场则对应多条成本曲线：中长期合约有约定电量和固定价格，日前现货有预测价格曲线，实时平衡市场有更高的惩罚性价格。售电公司在零售套餐定完、用户需求已知之后，才能将总负荷拆解分配到各个购电市场，目标是让总购电成本最小，同时抵御价格不确定性风险。
耦合点就在这里：用户的用电量取决于套餐价格，套餐价格又必须覆盖购电成本，而购电成本又取决于总用电量。这三者互相缠绕，形成一个均衡问题。如果不用主从博弈框架，你就只能固定其中一个变量，比如先算出一个假想的负荷曲线再去买电，那策略性能就差很多。这也是这类论文能够发到高水平期刊的根本原因：它的均衡概念强，变量联动充分，结果比直观方法更有说服力。
复现时建议先画一张“决策时间线”：售电公司先定什么、用户后定什么、市场再出清什么，然后在图上标注每个环节的决策变量和参数。这张图画清楚，公式建模就成功了一半。
2. 数学模型拆解：从目标函数到底层约束
下面以这类论文最常见的通用框架为例，具体符号和数值以你手上的目标论文为准。我尽量把每个式子背后的经济含义讲清楚，而不是只堆公式。
2.1 上层售电公司：利润最大化的完整表达式
上层目标函数看起来复杂，本质就是“卖电收入减买电成本，再减各项惩罚项”。以连续时段模型为例，设 t 表示时段，I 表示用户类型，J 表示套餐类型。售电公司的决策变量包括：套餐价格 p_{j,i,t}（或分时价格向量）、中长期购电量 x_f、现货购电量 x_s。
收入项是 Σ_{i,t} (用户i的套餐电价 × 用户i的实际用电量)。注意这里不能直接用用户的申报负荷，而要用下层求解出来的用电量 q_i,t。这就让收入项变成了“上层变量与下层变量耦合的乘积项”，也就是整个模型的非线性来源。
购电成本项则按市场分块。中长期购电是固定价 P_f，现货购电按日前价格 P_s，实时平衡价格带惩罚系数。很多论文还会额外加一个风险项，比如CVaR，用来衡量现货价格波动对利润的冲击。加风险项以后，问题就从确定性优化变成随机优化，场景数量直接决定求解规模。
写目标函数时，建议按“时段×用户×市场”三个维度逐项列出来，不要图省事写成一个总括公式。在Matlab里实现时，用矩阵运算把三个维度压在二维里，反而容易漏项。宁可先写成几个小sdpvar变量的加和，再合并成总的obj，这样后续检查每一项的经济含义都方便。
2.2 下层用户：效用最大化与套餐选择
下层用户的问题可以拆成两步：先挑套餐，再定用电量。很多初学者把这两步混成一个大优化问题，结果变量爆炸，求解极其困难。
常规处理方式是这样的：对每个套餐 j，用户在自己预算和套餐约束下选择一个最优用电量向量 q_j，得到一个最大净收益 U_j；然后用户在所有套餐里选 U_j 最大的那个，也就是引入0-1选择变量 y_ij 满足 Σ_j y_ij = 1。这样下层问题既包含连续变量 q_i,j,t，又包含整数变量 y_ij。
效用函数通常用二次型：U(q) = a q - 0.5 b q^2。为什么用二次型？因为它能让边际效用递减，曲线平滑，而且配上网电费以后求KKT条件很简单。a 可以理解为用户对电的基础需求系数，b 是需求弹性系数。参数设置时 a 要大于套餐边际电价，否则用户最优用电量会落到下界，模型退化。
这个双层-整数结构是复现时最大的拦路虎。处理思路一般有两条：一条是把整数选择变量留在上层，下层只保留连续用电决策；另一条是用大M法把下层整数问题的逻辑约束全部线性化，再取KKT。两条路都能走通，但第一条在代码实现上更直观，避免对混合整数问题直接求KKT带来的理论瑕疵。
2.3 双层问题单层化：KKT条件与互补松弛
把双层问题交给求解器前，必须先把下层优化问题转化为等价的单层约束。最常用的工具就是KKT条件。对连续的下层凸优化问题，KKT条件是必要且充分的。
设下层问题为：给定套餐价格 p，用户最大化 U(q) - p q，约束下限 q ≥ 0 和上限 q ≤ q_max。写拉格朗日函数，对 q 求导得到驻点条件。同时加上原问题可行、对偶变量非负、互补松弛条件，三者合起来就是下层的完整最优性条件。
KKT条件在Matlab里并不神奇，就是把几个等式和不等式写进约束列表。真正麻烦的是互补松弛条件：它形如“λ ≥ 0，g(x) ≥ 0，且 λ g(x) = 0”，是一个非光滑的非线性约束，不能直接丢给线性求解器。这时要用大M法把它离散化：引入0-1变量 z，把 λ g(x) = 0 拆成 λ ≤ M z 和 g(x) ≤ M(1-z) 两条线性不等式。M的取值要足够大，但不能大到让模型数值失真，这一点后面会专门讲。
2.4 非线性项线性化：大M法与辅助变量
除了互补松弛，模型里还有几处非线性需要处理。最典型的是上层收入项里的 p×q 乘积。套餐价格 p 是变量，用电量 q 是下层变量，两者相乘得到双线性项。如果 p 是连续变量，这类问题本质上还是非凸的，只能通过迭代或者特殊情况处理。
如果套餐价格被局限在几个离散档位上，问题就变简单了：引入0-1变量，把乘积项表示成“0-1变量×连续变量”，再用大M法拆成线性约束。很多EI论文为了可解性，都会把价格设计成离散候选集，然后整个问题退化为MILP。复现时最稳妥的做法，是沿用论文的离散化假设，不要擅自改成连续价格。
对用户选择变量 y 与用电量 q 的乘积，也可以用同样的辅助变量法线性化。一句话总结：把所有非凸双线性项，统统用0-1变量和上界常量去拆。模型看起来约束变多了，但求解器反而更愿意接。
2.5 参数标定：先跑小规模再上全规模
复现论文，最忌讳一上来就复现全部场景。我建议先构造一个最小可复现案例：1个售电公司、2类用户、3种套餐、24个时段、1个购电市场。用这个规模把模型跑通，验证目标函数值、套餐定价和各时段购电量的合理性，再把用户类型、套餐数量、场景数量逐步加上去。
参数标定有一些常见的经验值可以参考：居民用户的基础负荷系数 a 可以设在0.8~1.2元/kWh之间，b 设在0.01~0.05之间，这样24小时总用电量大致落在几十到几百千瓦时量级；现货价格曲线可以用典型日负荷曲线的反比关系生成，中长期价格略低于现货均价。注意这些参数要配合目标值量级做调整，而不是照抄。
很多复现失败其实是参数量级的问题：目标函数值动不动就上亿，约束的右端项只有几十，求解器的数值容差直接罢工。建议所有价格单位统一用“元/MWh”，电量用“MWh”，收益用“万元”，把量级压缩到10^3以内，求解稳定性会好很多。
3. Matlab代码实现：从零到出结果的几个关键环节
复现这类双层规划论文，代码结构比你想的重要。你可以直接用脚本堆，也可以按模块拆开，我倾向于模块化。
3.1 代码整体结构与文件划分
我一般按数据、模型、求解、结果四层来组织：
- data_xxx.m：生成所有输入数据，包括用户参数、套餐候选集、市场分段、预测价格曲线；
- build_single_level.m：把上层问题和下层KKT条件拼成单层模型；
- solve_main.m：调用Yalmip建模，设置求解器参数，求解并保存结果；
- post_process.m：把解还原成业务指标，比如各套餐签约用户数、总利润、购电结构占比；
- utils/：存放线性化辅助函数，把双线性项和互补约束打包，避免主文件太久太长。
这个分层的好处是调试范围可控。上层模型写错了，你只需要检查 build_single_level.m；数据错了，不用翻模型文件。很多初学者把所有代码挤在一个scripts里，300行以后连自己在哪一行都找不到，复现成本翻倍。
3.2 场景生成与数据预处理
多级市场购电策略的核心输入是各个市场的价格场景。常见做法是用历史现货价格聚类得到若干个典型场景，每个场景给一个概率。生成场景时要注意不同市场之间的相关性：日前价格高的时候，实时平衡价格往往也高，直接用独立随机数生成会让模型高价买入又高价卖出，出现看起来异常的最优解。
数据预处理里还有一个容易忽略的点：用户用电弹性不是跨类型统一的。居民用户在峰时段对电价敏感，工业用户对整体电价水平敏感。如果把所有用户都设成同一套弹性系数，模型结果会严重偏向某一类套餐，结论就失真了。复现时宁可先按论文的表把参数逐个录入，也不要偷懒用全统一参数。
我在处理时习惯把一切输入数据先画个图：负荷曲线、价格曲线、套餐候选集都在图上过一眼。数据进入模型之前先做可视化检查，比求解完发现结果荒谬再去回头查数据，效率高十倍。
3.3 Yalmip建模：目标函数和约束的写法
Yalmip是这类复现的首选建模工具。它把线性规划、整数规划、二次规划统一成同一套语法，改约束比直接写求解器API高效得多。核心代码骨架是这样的：
% 变量定义 p_pkg = sdpvar(nPackage, nTime, 'full'); % 套餐价格 x_f = sdpvar(nBlock, nTime, 'full'); % 中长期分块购电量 x_s = sdpvar(1, nTime); % 现货购电量 y_ij = binvar(nUserType, nPackage, 'full'); % 用户套餐选择 q_use = sdpvar(nUserType, nPackage, nTime, 'full'); % 各套餐下用电量 % 约束 C = []; % 每个用户只能选一个套餐 C = [C, sum(y_ij, 2) == ones(nUserType, 1)]; % 套餐电量上下限 C = [C, repmat(q_min, 1, nTime) <= q_use <= repmat(q_max, 1, nTime)]; % 市场购电平衡 C = [C, sum(sum(sum(q_use, 2), 3)) - sum(x_s, 2) == sum(x_f, 2)]; % 目标函数（最大化利润 → 最小化负利润） profit = sum(sum(sum(p_pkg .* q_use))) ... - sum(sum(x_f .* P_forward)) ... - sum(x_s .* P_spot); optimize(C, -profit, sdpsettings('solver', 'gurobi', 'verbose', 2));这段代码只是一个示意图架，真正跑之前还需要把KKT条件、互补松弛线性化、非线性项的辅助变量全部补充进来。但骨架能帮你确认一件事情：所有变量的维度是否对齐。维度对不对，是新手最常见的错误，Yalmip不会帮你自动广播维度，错了就是报错或NaN。
3.4 互补约束的代码写法
用大M法实现互补松弛，建议封装成一个函数，避免在主模型里反复写M和二进制变量。一个典型写法如下：
function C = complementarity(C, lambda, g, M) % lambda >= 0, g(x) >= 0, lambda * g = 0 % 引入0-1变量，把乘积约束拆成两条线性约束 z = binvar(size(g,1), size(g,2)); C = [C, g >= 0, lambda >= 0]; C = [C, g <= M * z]; C = [C, lambda <= M * (1 - z)]; end使用时注意 g 和 lambda 的维度要一致。M的取值一般取约束可能达到的上界，我的习惯是将M设成该约束实际最大值的1.5到2倍，而不是凭感觉给个10^6。M太大，求解器会出现数值问题；M太小，会错误地截断可行域。给M设值前，先跑一版没有互补约束的松弛解，看一下g和lambda的量级，再确定M。
3.5 求解器配置与结果校验
这类模型最终基本是MILP或MIQP，我默认用商用求解器。Yalmip配置Gurobi只有一行：sdpsettings('solver', 'gurobi')。如果机器上没有Gurobi，也可以用Cplex或者SCIP，但性能差距在整数变量多的时候非常明显。
求解器参数上，MIPGap默认是万分之一，复现初期建议放宽到1%甚至2%，先把结果跑出来。这里有个反直觉的点：对双层均衡模型，过于严格的MIPGap会让求解器在最优解附近反复转圈，迟迟不退出。放宽后快速拿到近似最优解，再用这个解作为热启动去求更严格解，反而更快。
结果校验是最容易被跳过的环节。拿到解以后，我最少做三件事：第一，把套餐价格和购电量代回下层用户模型，单独重新求解一遍，看用户的选择是否和上一级模型的 y_ij 一致；第二，检查KKT互补松弛约束的残差，看所有二阶乘项是否都接近0；第三，对比不同场景下的购电结构占总负荷的百分比，看中长期占比是否落在合理区间。这三道检查都过了，才敢说复现成功。
4. 复现过程中踩过的坑与排查方法
双层模型求解失败，绝大多数情况不是求解器能力不行，而是建模细节没到位。这里把我遇到的高频问题列出来，每条都附上排查思路。
4.1 问题一：模型报无解，不一定是模型错
求解器返回Infeasible的时候，第一反应别去改约束，先看是哪个约束导致不可行。Yalmip可以调用 diagnose(C) 或者通过求解器的 IIS 功能报告矛盾约束组。Gurobi支持找到最小不可行子集，Yalmip里把sdpsettings的iis开关打开就行。
我遇到最多的情形是购电平衡约束写反了方向。比如把“售电量 = 购电量”误写成“售电量 ≥ 购电量”，一开始不会报错，但加上现货平衡条件后就会出现无解。另一个高频原因是套餐选择约束和用户数量不匹配：某类用户在任何套餐下都不盈利，模型只能让它不选，但程序要求每个用户必须选一个套餐，于是不可行。这时候要么删掉这类用户，要么增加一个“不消费”选项。
检查不可行问题的顺序我建议是：先查约束的符号和方向，再查维度是否匹配，最后查M的上界是否太小。千万不要动辄把约束删掉来“试”，那样只会把模型越改越错。
4.2 问题二：互补松弛约束的大M怎么取
前面提过M的取值问题，这里展开说下。互补约束 λ g = 0 里的 g 可以是价格出清相关的等式残差，量级可能到10^3；而 λ 是对偶变量，量级可能小到10^-3。如果M统一取1e6，小量级的变量会被淹没，求解器认为 λ 和 g 的乘积“天然为0”，导致本该满足的均衡条件被跳过，结果看起来合理但其实是错的。
经验做法是分变量设置M。对每一组互补约束，估计自身g和λ的上界，然后取 mid = 1.5*max(ub_g, ub_lambda)。如果两个上界差了好几个数量级，建议把约束的单位换算到同一量级后再做。这个操作在复现前就要想好，因为M一多，查起来特别痛苦。
4.3 问题三：量纲不一致导致数值病态
量纲问题前面已经埋了伏笔。最常见的病态场景是：负荷以MWh计（几百到几万），电价以元/MWh计（几百），购电费就是10^6量级；而套餐价格如果写成元/kWh，和负荷一乘又变成10^7量级；再加上变量的上下界从小数到十万，数值范围横跨五个数量级，求解器的预求解很难处理。
我的习惯是统一用“元/MWh”和“MWh”，所有收益和成本以“万元”显示。这样目标函数的数值在10^3左右，约束右侧在10^2左右，求解器跑起来稳定很多。另外一个细节是不要用过大或过小的惩罚系数，比如10^7这种，会让线性规划条件数变大，预处理失效。
4.4 问题四：求解时间爆炸
双层模型变成MILP后，解空间会比普通MILP大得多，求解时间爆炸是常态。首先检查模型里是不是不小心产生了太多二进制变量。每个互补约束引入1个0-1变量，如果互补约束本身是向量，那0-1变量数就是维度乘积，很快突破几十万个。
应对方法有几条：场景削减，用聚类把50个价格场景压到5个典型场景；边界压缩，把套餐价格候选集从100档压到20档；热启动，先用小规模算出一个解，固定部分变量再求精细解。另外，MIPGap放宽到1%往往能省下一半时间，对复现研究来说精度完全够。
4.5 常见问题与解决方案速查表
| 现象 | 常见原因 | 解决思路 |
|---|---|---|
| 求解器报Infeasible | 购电平衡方向写反、套餐选择约束过死 | 用IIS定位矛盾约束组 |
| 目标值出现NaN | 变量维度错位、稀疏矩阵有NaN初始值 | 逐项打印收入的每个加和项 |
| 套餐价格被压到上界 | 购电成本项漏了或多算了一倍 | 检查购电成本中是否重复计入负荷 |
| 互补约束残差过大 | M值太大或太小 | 分约束设置M并校核量级 |
| 长时间不出最优解 | 整数变量过多、MIPGap过紧 | 削减场景、放宽MIPGap、热启动 |
| 结果购电结构异常 | 场景生成忽略了市场价格相关性 | 用聚类或copula生成场景 |
这张表是我在实际复现中反复用到的自查清单，每个问题都至少遇到过。建议把这几个问题挂在代码文件的注释里，跑挂了照着查一遍。
5. 从复现到自己的研究：几个值得扩展的方向
复现论文不是终点。把代码跑通之后，更值钱的事情是思考怎么改模型、加场景、换机制，把它变成自己的研究增量。
5.1 扩展方向一：引入更多层级的市场互动
原论文的多级市场购电如果是三层，你可以考虑引入日前与实时两市场的决策顺序。现实中，日前市场先出清，售电商根据日前结果调整实时补量。把这个时序放进去，模型就从静态Stackelberg变成带时序的MPEC，虽然求解更难，但结论更贴近实际。
5.2 扩展方向二：把单售电商扩展为多售电商竞争
原模型是单一售电商面对用户的Stackelberg博弈。如果你感兴趣的是零售侧竞争，可以把上层改成多个售电商同时决策，下层用户仍然做最优选择，这样模型从单领导者变成多领导者竞争，属于EPEC问题。EPEC的求解比MPEC更麻烦，经常要用对角化或启发式方法，但做成了就是很好的创新点。
5.3 扩展方向三：把确定性模型写成随机规划
原论文如果在价格和负荷上用的是确定性场景，那么很自然可以扩展成两阶段随机规划：第一阶段定套餐价格，第二阶段根据场景出清购电并结算。场景削减、机会约束、CVaR风险度量的加入都会显著提升论文的工作量和现实意义。
5.4 给复现新手的几条实操建议
第一条建议：把“跑通代码”和“复现论文”分开。跑通代码只是起点，能在你的机器上重复出论文的关键数值和趋势图，才算复现成功。第二条建议：养成改一个参数就记录一组结果的习惯。双层模型的参数敏感性很强，某个系数差10%，结果就可能从“全买中长期”变成“全买现货”，没有日志你根本无从排查。第三条建议：多和别人讨论KKT转化那一步。我发现大部分做这类问题的同学卡住，不是卡在代码，而是卡在“为什么能用互补松弛”这个概念上。把这个概念讲清楚一次，后面整个复现都会顺很多。
我个人的体会是，复现这类EI论文最大的收获不是代码本身，而是你终于理解了“博弈均衡”在真实交易里是怎么被算出来的。论文里的曲线和数值都只是结果，真正值钱的是那套把两层决策压缩成一组线性约束的思维方法。这套方法学会了，后面做电动汽车充电定价、微电网内部调度、需求响应激励机制，都会发现全是同一个套路。
