电力市场方向的Matlab复现坑有多深，估计踩过的人都懂。我最近把一篇基于主从博弈的售电商零售套餐设计与多级市场购电策略的EI论文复现完了，从建模到代码调通用了大概两个星期。这篇文章把里面的模型设计、数学转化和Matlab实现细节一次性说清楚，希望能帮你省掉大把试错时间。
先说这个项目是干嘛的。售电公司（也叫售电商）从批发市场买电，再打包成各种零售套餐卖给终端用户。批发侧有中长期市场、日前市场、实时市场，价格波动大；零售侧又要面对不同类型用户，有人对电价敏感，有人根本不在乎。怎么在套餐里设计价格信号，让用户主动削峰填谷，同时自己又在多级市场里把购电成本压到最低，这就是主从博弈要解决的问题。适合看这篇的，基本是电力市场、能源经济方向的研究生和刚转做售电业务的工程师，需要对优化理论和Matlab有基础，但不用太深。
1. 项目整体设计与思路拆解
1.1 售电商在电力市场中的角色与痛点
售电商本质上是做“批发-零售”价差生意的中间商。传统模式下，售电商只做一件事：按固定价格从上级购电，再按另一个固定价格卖给用户，赚一个固定价差。这套玩法在管制时代行得通，因为批发侧成本相对稳定，用户侧没有选择权。市场放开之后，问题就来了：批发侧价格会随供需实时波动，尤其是日前和实时市场价格可能一天之内差出好几倍，而零售侧还是固定价的话，价差风险就全部压在售电商身上。
再往深一层看，用户用电习惯天然是刚性的，但并非完全不可调。如果零售侧永远是一个固定价格，用户没有任何动力改变用电行为，高峰时段该用还是用，售电商被迫在高价时段去多买电，成本只能自己扛。如果能设计一套零售套餐，让用户在高峰时段主动压负荷、低谷时段多用电，售电商的购电曲线就能被削峰填谷，整体采购均价也能降下来。这才是“零售套餐设计”的真正意义。
这个项目里要解决的痛点可以拆成两个：一是怎么设计套餐，让用户愿意跟着价格信号走；二是怎么在不一样的市场里购电，使得最终的期望收益最大、风险可控。两个问题不是孤立的，套餐设计决定了售电商与用户之间的账单结算，购电决策决定了成本端，二者在主从博弈框架里天然耦合。
1.2 为什么选主从博弈而不是集中优化
很多人拿到这个题目会先想到：把所有用户和售电商放一起做集中优化，直接最大化全网或全公司的总收益不就行了？我一开始也觉得集中优化简单，但复现之后发现这和市场机制根本对不上。
集中优化隐含一个假设：用户是“可控单元”，让你怎么用就怎么用。可现实里的用户有自己的效用函数，面对电价会自主决定用多少电。你强制压低某时段用电，用户会不干，因为影响了他自身的舒适度和生产需求。想要准确刻画这种“你定规则、我选行为”的交互关系，Stackelberg主从博弈是更合适的人选。
在主从博弈框架下，售电商是领导者（leader），先发布零售套餐价格；用户是跟随者（follower），看到价格后根据自身效用最大化决定用电量。售电商在预期到用户会这样响应的情况下，再反过来优化套餐和购电策略，最终达到一个均衡点：在这个点上，售电商单方面改套餐不会增加利润，用户单方面改用电量也不会增加效用。这个均衡就是论文里说的“均衡解”。
这里有个关键点不能忽略：主从博弈得到的是纳什均衡的一个特例，而不是全局最优。售电商作为领导者之所以能占便宜，是因为他先动，可以把用户的响应行为当作内生变量放进去优化。而集中优化把用户当牵线木偶，虽然在数学上更容易解，但结果在真实市场里不可执行。这也是评审和读者最看重“建模是否符合市场物理结构”的原因。
1.3 多级市场购电策略的设计逻辑
“多级市场”指的不是一个平台，而是不同时间尺度上连续衔接的电力交易市场。大致可以分成三层：中长期合约市场、日前市场和实时市场。
中长期合约市场很早就要签订，价格相对稳定，但电量不一定能完全匹配实际负荷，因为几个月前没法精确预测用户到底用多少电。日前市场在运行日前一天开市，可以根据更准确的负荷预测去买缺口电量，价格比中长期灵活。实时市场离实际运行最近，用来处理最后的电量偏差，但价格波动最剧烈，可能高得离谱。
购电策略设计就是在三个市场上做电量分配。如果负荷预测准，就可以多在日前买、少在实时买；如果预测不稳，就需要靠中长期合约兜底，或者承受一定的实时偏差风险。再进一步，还可以在优化目标里加入风险度量，比如CVaR（条件风险价值），防止出现极端价格场景下收益大幅恶化的情况。
把这个多级购电决策和零售套餐设计放在同一个博弈模型里，是这篇论文比较“顶”的地方。套餐价格会影响用户用电量，用电量又决定需要在各个市场上买多少电，而各市场的价格又直接影响售电商收益，最终反过来影响套餐定价。三件事必须联动优化，缺一块都不完整。
2. 核心模块解析与实操要点
2.1 零售套餐的数学建模：分时、阶梯、固定型怎么选
实际项目中，零售套餐的类型有很多种，但论文里常用的大概有三类：固定价格套餐、分时电价套餐和阶梯电量套餐。
固定价格套餐最简单，所有时段、所有电量都执行同一个单价。数学上就是一个常数p_fixed，用户不管什么时候用电，结算单价不变。这种套餐对用户友好，但售电商没有调节负荷的工具，通常只在用户基数大、负荷曲线平稳时才划算。
分时电价套餐把一天分成峰、平、谷时段，分别给不同的电价。设时段集合为T，每个时段t的执行电价为p_t，用户在时段t的用电量为q_t，那么用户的电费就是Σ_t p_t q_t。这种套餐能让人直观感知价格信号，实践中最常用。规划设计核心是怎么划分时段、怎么设定峰谷价差。价差太小，用户没反应；价差太大，用户可能反弹，甚至导致负荷大幅转移，反过来影响售电商的购电成本。
阶梯电量套餐是按总用电量分档计价，比如月用电量小于X度时按低价，超过X度后按高价。数学上通常用分段线性函数表示，会引入整数变量或者特殊顺序集合（SOS2）来处理分档，实现复杂度比分时高一些。在博弈模型里，阶梯套餐会让用户问题变成一个带离散选择的优化问题，对求解器要求更高。
我复现时选用的是“分时套餐为主、固定套餐为辅”的组合。理由有两点：第一，分时电价本身是连续变量，能让KKT转化后的模型保持结构；第二，论文里强调“多元零售套餐”，意味着用户可以在多种套餐里选择，细分用户群体。我把用户按价格弹性分成几类，每一类只能选择对应的一种套餐，避免用户同时拥有多套餐的选择权而引入复杂的离散博弈。
2.2 用户响应模型与需求弹性
用户响应是整个博弈模型的“灵魂”。如果用户对价格没反应，那售电商怎么优化套餐都没用。刻画用户响应的常用方法是给出一个效用函数。假设用户在时段t用电q_t，获得效用U_t(q_t)，同时支付电费p_t q_t，用户的目标就是最大化净效用：
max Σ_t [U_t(q_t) - p_t q_t]
为了保证用户问题的凸性、确保KKT条件成立，效用函数一般取凹函数。我用的最常见形式是二次效用函数：
U_t(q_t) = a_t q_t - (b_t / 2) q_t^2
其中a_t和b_t是用户偏好参数，a_t可以理解为基础意愿支付价格，b_t是效用衰减系数。对q_t求导并令为零，能得到线性需求函数：
q_t = (a_t - p_t) / b_t
这就是经典的价格-电量线性响应关系。1/b_t就是价格弹性系数，b_t越小，用户对价格越敏感，电价一波动用电量就大幅变化；b_t越大，用户越“不差钱”，怎么涨价都不太改习惯。实际设定参数时，我会参考论文里的典型数值范围，再结合负荷基线数据微调。
需要特别提醒的是，用户并不只有电量这一个决策维度。真实场景中用户还可能有储能、电动汽车、可转移负荷等资源，这些都会在后续优化模型中增加约束。为了控制复现复杂度，第一版我只保留弹性负荷约束：每个时段用电量有上下限，一天总用电量保持在一个允许范围内。这样既能体现用户响应，又不会让下层问题变成难解的混合整数问题。
2.3 多级市场购电模型与风险度量
购电策略的建模核心是“总量平衡”。一天24小时，每个时段售电给用户的电量必须等于从各个市场买到的电量之和。我按论文的设定，把购电分为三部分：中长期合约电量、日前市场购电量和实时市场平衡电量。设时段t的总售电负荷为Q_t（等于所有用户该时段用电量之和），则：
Q_t = Q_t^contract + Q_t^day + Q_t^real
其中Q_t^contract是中长期合同电量，执行合同电价c_t^contract；Q_t^day是日前市场购电量，价格c_t^day；Q_t^real是实时市场购电量，价格c_t^real。这里有一个隐含假设：实时电量差可能为正也可能为负，如果实际负荷低于购电量，售电商可以在实时市场卖出多余电量，这时候Q_t^real为负，收益为正。
难点在于日前和实时价格是不确定的。为了刻画不确定性，我用场景法：根据历史价格数据，通过蒙特卡洛抽样生成S个价格场景，每个场景有出现概率ρ_s。那么售电商的期望购电成本就是所有场景下购电成本的平均值。不过期望值只能描述平均水平，碰上极端高价场景可能亏到怀疑人生，所以论文里通常都会加一个风险度量。
最常用的是CVaR。说人话就是：在所有可能出现的情景里，把最差的那5%或10%情景挑出来，计算这些最差情景的平均收益或平均损失，把它作为一个风险指标纳入优化目标。数学上可以用一个辅助变量加上一组约束线性化表示。加入CVaR后，购电策略会变得更加保守：风险厌恶系数越高，售电商越愿意多签中长期合约、多买日前电量，减少对实时高价市场的依赖。
我在代码里用了一个通用写法：收益的期望值减去β × CVaR，β就是风险厌恶系数。β等于0时是风险中性，只追求平均收益；β越大，越偏向低风险高确定性的购电组合。
2.4 主从博弈的KKT条件与强对偶重构
这是整个复现里最烧脑的部分，也是论文“数学含金量”所在。主从博弈原本是一个双层优化问题：上层售电商决策套餐和购电量，下层用户决策用电量。双层问题直接扔给求解器是解不了的，必须先把下层问题转化为单层约束，才能变成一个单层数学规划。
转化核心是“用用户的KKT条件替换用户优化问题”。因为用户问题是一个凸二次规划，满足正则条件，所以它的最优解一定满足KKT条件。我们把这些KKT条件（包括站点条件、对偶可行性、互补松弛等）全部写成约束，附加到上层问题中，原来的双层问题就变成了一个带互补约束的单层优化问题（也就是MPEC）。
KKT条件里最麻烦的是互补松弛条件，形式是“两个非负变量相乘等于0”。这种条件用传统线性规划解不了，需要用大M法做线性化。具体做法是：引入布尔变量，当其中一个乘子大于0时强制另一个为0。大M法里M的选取特别关键，这个我后面在常见问题里详细说。
另一个转化路线是“强对偶定理”。用户问题是凸的，满足强对偶，所以上层问题中可以把下层目标函数替换成它的对偶目标函数，再结合用户的原始约束和对偶约束，同样能得到单层问题。两条路线殊途同归，实际代码中我建议先用KKT条件，因为对YALMIP建模的人来说更直观，但要注意最后形成的MPEC问题中会产生双线性项，需要用M法处理。
如果你不想一上来就啃MPEC，也可以用“差分进化+内层优化”的迭代方法：外层由售电商更新价格，内层求解用户问题，然后根据用户的反应再修正价格。迭代法的优点是实现简单，缺点是不保证收敛到真正的Stackelberg均衡，遇到非凸情况容易振荡。论文复现如果想严谨，还是优先走单层MPEC路线。
3. 实操过程与核心环节实现
3.1 代码整体框架与文件组织
Matlab代码结构我建议按模块拆分，别把所有东西塞在一个脚本里。我复现时的工程目录大概长这样：
project/ ├── main.m % 主程序入口 ├── config/ │ └── parameters.m % 所有参数集中设置 ├── data/ │ ├── load_profile.xlsx % 用户负荷基线 │ └── price_history.xlsx % 历史市场价格 ├── scenarios/ │ └── generate_scenarios.m % 生成日前/实时价格场景 ├── model/ │ ├── build_user_model.m % 构建用户下层优化问题 │ ├── build_retailer_model.m % 构建售电商上层优化问题 │ └── kkt_linearization.m % KKT条件与大M线性化辅助函数 ├── solver/ │ └── solve_mpec.m % 调用商业求解器解MPEC └── results/ ├── plot_results.m % 绘图脚本 └── save_outputs.m % 保存结果main.m 的逻辑很简单：加载参数，生成场景，构建模型，求解，画图。实际开发时我会先把parameters.m里的所有参数确认好，再跑模型。因为参数一旦中途改动，后面所有结果都乱套。
3.2 关键参数与场景生成
参数设置直接决定结果的合理性。我以24小时调度周期为例，时段数T=24。用户数量初始设为10个，分成3类：高弹性、中弹性和低弹性。每个用户的负荷基线从某个典型负荷曲线缩放而来，避免所有用户曲线都一样，否则博弈结果会退化。
价格场景生成我建议不要用太复杂的随机过程。第一版用历史电价的均值向量和协方差矩阵，直接调用Matlab的mvnrnd生成正态分布场景。代码如下：
mu_day = mean(day_price_history, 1); % 24×1 sigma_day = cov(day_price_history); % 24×24 scenarios_day = mvnrnd(mu_day, sigma_day, S);要注意的是，正态分布会生成负价格，虽然现货市场偶尔也有负电价，但大部分论文场景里不希望出现太多负值。我会在采样后做截断处理，把小于0的数值统一设为0.01。更精细的做法是用Copula或历史场景重采样，但那会增加复现工作量，性价比不高。
场景数量S我最初设了500，结果MPEC模型太大，求解慢得让人抓狂。后来降到200，再配合随机种子的固定，结果差异很小，速度却快了很多。如果你也要跑大规模参数，建议先从S=50开始调通，再逐步增加。
另外，中长期合约价格我设为一段固定均值，只在部分时段允许签约，模拟实际合约约束。合约电量上限按用户总负荷的某个比例设置，比如最高能覆盖总电量的60%，剩下的必须走日前和实时市场。
3.3 主从博弈求解的具体流程
如果走MPEC路线，整体流程可以概括为：
- 初始化所有参数，生成场景。
- 用YALMIP定义上层决策变量：
- 套餐价格向量
p_t（可能有多类套餐）； - 各市场购电量
Q_t^contract,Q_t^day_s,Q_t^real_s； - 用户用电量
q_{i,t}； - 辅助风险变量
VaR,z_s。
- 套餐价格向量
- 写出售电商目标函数：收益减去购电成本减去风险惩罚。
- 加入下层用户的KKT条件作为约束，包括：
- 用户一阶最优条件；
- 对偶变量非负约束；
- 互补松弛条件（用大M法线性化）。
- 调用求解器求解单层MPEC问题。
- 检查均衡条件，输出结果。
这里给出一个核心的YALMIP建模片段（伪代码风格）：
% 变量定义 p = sdpvar(T, 1); % 套餐价格 Qc = sdpvar(T, 1); % 中长期购电 Qd = sdpvar(T, S); % 日前购电，场景相关 Qr = sdpvar(T, S); % 实时购电 q = sdpvar(N, T); % 每个用户的用电量 lambda = sdpvar(N, T); % 用户问题的一阶乘子 mu_lb = sdpvar(N, T); % 下界乘子 mu_ub = sdpvar(N, T); % 上界乘子 % 目标函数 retail_revenue = sum(p .* sum(q, 1)); purchase_cost = sum(Qc .* c_contract, 1) + mean(sum(Qd .* c_day + Qr .* c_real, 2)); risk_cost = beta * cvar_term; objective = - (retail_revenue - purchase_cost - risk_cost); % 用户KKT一阶条件 constraints = []; for i = 1:N for t = 1:T constraints = [constraints, a(i,t) - b(i,t) * q(i,t) - p(t) + lambda(i,t) - mu_lb(i,t) + mu_ub(i,t) == 0]; end % 互补松弛条件用大M法线性化，这里省略具体代码 end注意上面目标函数里我用了mean(sum(...))，这是把所有场景的购电成本做了平均，对应期望值。风险项我没有展开，因为CVaR线性化需要额外变量和约束，网上有很多现成模板，直接套用即可。
求解器配置是另一个重点。YALMIP只是建模语言，真正解MPEC还需要底层的商业求解器。我用的商业求解器支持二次约束和整数变量。如果求解器不支持非凸二次项，可能需要把双线性项通过M法线性化。模型规模比较大时，求解时间会很长，建议在求解器参数里开启并发计算，并设置最大求解时间。
3.4 结果输出与敏感性分析
结果输出至少要包含四张图：零售套餐价格曲线、用户总负荷响应前后对比、三个市场购电量分配、迭代收敛曲线（如果用迭代法）。我在实际跑出来之后通常会先看套餐价格曲线，判断峰谷价差是否在合理区间。如果峰时段价格是谷时段的2到3倍，说明套餐设计已经能有效引导用户避峰；如果价差太小，用户响应不明显；价差太大，用户可能把负荷全部转到低谷，造成新的“谷峰倒置”，反而增加日前市场购电压力。
购电分配图上能看到一个明显规律：风险厌恶系数β较小时，日前市场购电占比高，实时市场也有一定比例，因为这样可以抓住低价机会；β变大后，中长期合约占比稳定升高，实时市场几乎被“压”到只能平衡偏差，体现出风险规避特征。
敏感性分析我做了三个维度：用户弹性系数、价格场景数量和风险厌恶系数。用户弹性系数整体放大时，同一个套餐价格下用户响应幅度更大，售电商的峰谷价差会适当缩小，因为不用订太高的价差用户就动了。价格场景数量增加到一定值后，购电策略基本稳定，说明场景数量不是越多越好。风险厌恶系数的影响我刚才说了，就是让购电组合向确定性市场倾斜。
我的一个体会是：论文里的结果图看着很理想，但自己跑下来的结果往往没有论文里那么平滑。主要是因为参数没有做精细的标定。你可以通过微调用户效用参数a_t和b_t，让响应后的负荷曲线不至于太夸张，再去对齐论文中的趋势图。
4. 常见问题与排查技巧实录
4.1 YALMIP变量定义与求解器选择容易踩的坑
第一个常见的坑是：把普通数值矩阵和sdpvar变量运算混在一起，经常报“Conversion to double from sdpvar not possible”。这是因为YALMIP的重载不支持某个操作，或者变量维度不匹配。我的经验是：只要涉及优化变量，所有矩阵运算都放到YALMIP环境中，不要试图先计算再转换。
第二个坑是求解器路径没有配置好。YALMIP本身不是一个求解器，它只是把你的模型翻译成求解器能吃的标准形式。如果你在命令行直接输入optimize(constraints, objective)而没有任何商业求解器，YALMIP会调用内置的linprog或fmincon，遇到混合整数和二次约束就崩了。解决方法是提前在Matlab里设置求解器路径，并在optimize函数中明确指定用哪个求解器。如果没有可用的商业求解器，至少要用sdpt3这类免费求解器先做可行性测试。
第三个坑和模型规模有关：YALMIP内部的变量展开会导致求解器内存爆掉。尤其是场景数量S几百、用户数量N几十的时候，变量数很容易过万。建议给每个变量加注释，构建约束时用稀疏方式而不是循环拼接。循环里每写一次constraints = [constraints, ...]都会增加内部复杂度，对性能影响很大。最好先构建多个约束块，最后一次性合并。
4.2 迭代不收敛与振荡问题
如果你用的是迭代法而不是单层MPEC，最常遇到的是价格在两次迭代间来回跳，就是“振荡”。最典型的现象是：峰时段价格从0.8跳到1.2，下一次又返回0.8，反复几次不收敛。
振荡的根源在于下层用户问题可能存在多解，或者上层目标对价格过于敏感。一个很有效的修复方法是“加步长阻尼”：每次迭代更新价格时不要直接取最新解，而是：
p_new = p_old + α (p_candidate - p_old)
其中α取0.2到0.5之间。比如设α=0.3，价格变化就平滑得多。刚开始我为了图快，α直接设1，结果振荡到怀疑人生。加了阻尼之后，大概20轮左右能稳定下来。
另一个手段是给用户问题加正则项，比如在用户目标里加一个很小的二次惩罚ε * ||q||^2，把原问题变成强凸问题。ε不能太大，否则会改变原博弈的均衡解；一般取10^-6量级就够了。这个做法在数学上是加了一个小扰动，但会显著改善数值稳定性。
如果用了单层MPEC还是出现“解不稳定”，那很可能大M法中的M取值不合适，导致互补约束没有被正确满足。这时可以检查Kkt互补松弛残差，如果没有趋近于0，就调整M的值，或者把互补条件改写成带小阈值的约束。
4.3 数值问题：大M法取值与尺度缩放
大M法在MPEC转化里几乎是必须的，但M的取值是一个技术活。M太小，比如设为10，某些乘子可能超过M，原本应为0的乘积被错误地允许为正，导致满足约束的解根本不是真正的互补解。M太大，比如设成10^9，商业求解器内部数值稳定性变差，容易出现“数值病态”，报错不可行或求解缓慢。
我的经验是，M和变量数量级相关。电价一般在几十到几百元/MWh，电量在几十到几百MWh，乘子大小跟目标函数梯度有关，通常也在10^2到10^4量级。因此M先在1000到10000之间试，然后检查互补松弛条件的实际残差，再逐步调整。不要一开始就放一个超大的M。
尺度缩放是另一个容易被忽略的点。售电商一天的收入可能是几万块，而某个乘子可能是0.0001，两个数量级差太大，求解器误差容限就失效了。我建议先把所有变量做归一化，比如把电量价格统一除以基准值，让目标函数与约束的量级都在1附近，求解成功后反归一化再展示结果。这一步能减少很多无解的困扰。
4.4 调试小技巧与复现建议
调试这类模型时，最忌讳一上来就调全规模参数。我的调试顺序是：先设N=2、T=4、S=1，把模型跑通。T=4相当于只模拟凌晨、上午、下午、晚上四个时段，N=2代表两类用户，S=1就是确定性场景，跳过所有随机因素。这样模型体量小，任何错误都能从约束或目标里快速定位。跑通确定性案例后再逐步增加T、N和S。
另外一个非常实用的技巧：单独验证用户下层问题是否正确。先固定一组套餐价格，单独求解一个用户的优化问题，看看最优电量是否等于解析解q=(a-p)/b。如果解析解和数值解对不上，说明用户问题有约束写错或者双线性项加错地方。这个检查比直接跑整个博弈靠谱得多，因为下层问题错漏往往会传导到最终的MPEC里，让人查无可查。
调试时可以把所有约束目标写成文本导出，用稀疏矩阵可视化查看非零元分布，能快速找到异常链接。我见过很多次所谓“模型不可行”，其实就是某个时段电量平衡约束里漏乘了一个用户数量缩放系数，导致供求永远对不齐。
最后再分享一个实用建议：这类复现项目不要执着于“参数和论文一模一样”。顶刊论文里的数值结果往往经过了大量标定，你复现时只需要抓到核心模型和趋势。先把机理层跑通，再根据你自己的数据调出合理曲线，这样写论文或做汇报时更有说服力。我接下来准备把多类用户对套餐的自选择行为用更严格的离散选择模型表达，再对比一下和当前连续近似的差距，这个方向有进展的话会继续更新。
