先说点实在的。EI论文复现这件事最容易被卡住的往往不是数学公式看不懂而是“模型怎么落到代码上”这一步。主从博弈、零售套餐、多级市场购电这三个词单独拎出来都有大量参考资料但组合到一起就需要你同时处理好三层问题售电商在上层怎么定套餐价和购电量、用户在下层怎么根据价格调整用电行为、以及两个层级之间的均衡怎么求解。这篇博客就按我实际复现这类模型的顺序把建模逻辑、数学转化、Matlab实现和调参过程中踩过的坑一次说清。如果你正在做电力市场方向的研究或者刚接触售电侧优化、想找个完整的双层优化案例练手这篇内容应该能帮你少走不少弯路。下面从最核心的问题开始聊。1. 复现前的认知梳理主从博弈与售电问题的匹配逻辑1.1 为什么单层优化做不了售电套餐设计很多刚接触售电定价的人第一反应是建一个利润最大化模型售电商决策变量是电价用户用电量当成固定参数直接求最优解。这个思路不是不能用但它隐含了一个很强的假设——用户不管电价怎么变用电习惯完全不变。现实显然不是这样。居民用户看到峰时电价涨到1块2会把洗衣、充电挪到晚上10点以后工商业用户收到分时电价方案会调整生产班次降低电费支出。也就是说用户的用电量是电价的函数电价一变用电量跟着变售电商的收入也随之变化。这正是经济学里的“需求弹性”也是单层优化模型最大的盲区。把这个问题简单用数学语言描述售电商要最大化自己的利润利润等于售电收入减购电成本而售电收入又依赖用户用电量用电量又依赖售电商定的价格。这里存在一个双向耦合的关系而不是单向的“定好价格→算利润”。所以单层优化天然缺了一个环节——用户侧的价格响应行为。1.2 主从博弈的分层结构与电力市场的对应关系主从博弈Stackelberg Game恰好是描述这种“先动者—后动者”关系的标准框架。售电商先公布套餐价格相当于领导者Leader用户看到价格后选择用电量相当于跟随者Follower。领导者利用自己对跟随者反应函数的预判来做决策实现自身利益最大化。这个顺序决策特性和电力零售市场的实际业务逻辑高度吻合售电商发布套餐通常以月或季度为周期价格一旦公布很难频繁调整用户根据套餐价格调整用电行为这是一个日常、连续的自适应过程售电商在制定套餐时必须预判用户会怎么反应。这种天然的主从关系让Stackelberg博弈成为建模售电商决策的主流工具。相比纳什均衡那种“同时决策、互相猜”的框架主从博弈更贴近实际因为市场中总有一方先出牌。1.3 模型的时序安排购电与售电的时间尺度怎么统一在搭模型之前有一个前提必须理清楚购电决策和售电决策发生在不同的时间尺度上。中长期市场购电可能提前几个月就锁定了日前市场购电是运行日前一天的决策而零售套餐的价格一旦定了至少在一个结算周期内是不变的。如果把购电策略和套餐设计放在同一层模型里就需要把它们统一到一个决策时刻来考虑。实际论文里的常见处理方式是把购电策略和零售套餐设计都放进“售电商的上层决策”它们在同一个优化框架内求解只是不同市场购电量的成本参数不同。也就是说售电商在制定套餐的时候同时决策中长期、日前、实时市场的购电量分配二者共同构成上层目标函数中的购电成本项。这样处理后一个完整的双层决策模型就是上层售电商决定零售套餐结构、各时段电价、各类市场购电量目标是最小化购电成本、最大化售电收益。下层用户根据电价决定各时段用电量目标是最大化自身效用减去电费支出。有了这个框架下一步就是具体的建模细节。2. 零售套餐设计怎么建模从用户效用到多类型套餐的统一框架2.1 用户效用函数选择为什么是二次函数用户侧建模是整个主从博弈下层问题的核心。电力用户对用电量有一个主观满足度经济学上用效用函数来描述。最常用的形式是二次效用函数[ U(q) a \cdot q - \frac{1}{2} b \cdot q^2 ]其中a和b是正系数q是用户的用电量。二次函数的边际效用是递减的——用电量从0增加到1个单位时用户的满足感提升很大从100增加到101个单位时满足感的提升就很小了。这和“家里多开一盏灯的感觉”完全一致符合人们的基本直觉。用户的目标是最大化净收益效用减去电费支出。[ \max_q \quad a \cdot q - \frac{1}{2} b \cdot q^2 - p \cdot q ]对q求一阶导并令其为零可以得到用户对电价的响应函数[ q \frac{a - p}{b} ]这个式子非常直观电价p越高用户用电量越低参数a越大用户的基本用电需求越高参数b越大用户对价格越敏感。这个响应函数会直接成为下层问题KKT条件里的核心平衡方程。2.2 不同类型零售套餐的数学结构所谓“多元零售套餐”通常包括固定电价、分时电价、阶梯电价、实时电价等类型。它们在数学建模上的差异主要在电费计算方式固定电价套餐电价是常数所有时段一个价。用户的电费支出是 ( p \cdot \sum_t q_t )建模最简单。分时电价套餐把一天分成峰、平、谷等时段每个时段一个电价。用户的电费支出就是各时段电量与对应时段电价相乘之和。这也是目前国内售电公司最常见的套餐形态。阶梯电价套餐随着用户累计用电量增加单价逐级上升。需要引入分段线性函数描述“电量落入哪个阶梯”通常会引入二进制变量。实时电价套餐电价直接绑定日前市场出清价格用户承担价格波动风险售电商赚取固定服务费。这种套餐在前场建模中相对少见但可以作为方案对比。在代码实现时可以用一个统一框架处理这些不同的套餐结构无论哪种套餐用户的电费支出最终都能表示成“价格参数”和“用电量”的某种组合方式。常见的做法是定义 ( C(p, q) ) 作为用户的电费支付函数然后在售电商收入项里把这个函数写进去。套餐类型不同只是函数形式不同博弈的整体求解结构不变。2.3 套餐差异化价格上下限、时段划分与用户分类套餐设计不能完全自由定价实际业务中还要考虑监管约束和市场接受度。首先是价格上下限约束。多数地区的零售套餐电价不能超过监管机构规定的上限也不能低于成本价太多否则会引发恶性竞争。建模时通常给电价加一个区间约束[ p_{\min} \le p_t \le p_{\max} ]其次是时段划分。分时电价套餐的“峰平谷”时段划分通常是固定的或由售电商在一个周期内决定不同时段的电价不能差异太大否则用户会剧烈调整用电行为反而影响售电商收益的稳定性。再次是用户分类。不同用户群体的效用函数参数不同——住宅用户、商业用户、工业用户的弹性差异很大。进阶的论文会引入多类用户每类用户对应一组 ( a_i, b_i ) 参数在下层形成多个用户子问题。这种做法虽然增加计算量但更贴合真实市场。在我实际复现的过程中建议第一次先把用户类别固定成单一类型把单类用户的分时电价模型跑通再扩展多类用户的版本。直接上多类用户版本可能会让排查问题难度成倍增加。3. 多级市场购电的决策链条中长期、日前与实时的联动逻辑3.1 三种市场的定位锁风险、调偏差、算平衡售电商向上游买电通常面对三个市场中长期市场以月、周甚至年为周期通过双边协商或集中竞价签订合同电量。这个市场的电价相对稳定可以锁定大部分的基础购电成本规避实时价格波动的风险。日前市场运行日前一天售电商根据对次日用户用电量的预测确定次日各时段的购电量。这个市场的价格会实时波动但比实时市场温和得多。实时市场运行当天用户的实际用电量和日前购电量大概率存在偏差这部分偏差电量在实时市场结算。实时市场价格波动最大通常只用来兜底。3.2 购电成本的数学表达和偏差电量处理用一个简化的购电策略模型来说明假设一共有T个时段比如24小时售电商在日前市场的购电量为 ( q_{DA,t} )在中长期市场已经签署的电量为 ( q_{LT} )用户实际用电量为 ( q_t )。那么实时市场需要结算的偏差电量为[ q_{RT,t} q_t - q_{DA,t} - q_{LT} ]购电总成本可以写成[ C_{buy} \sum_t \left( c_{LT} \cdot q_{LT} c_{DA,t} \cdot q_{DA,t} c_{RT,t} \cdot q_{RT,t} \right) ]需要注意( q_{RT,t} ) 可能是正的少买了电需要在实时市场高价补购也可能是负的多买了电可以在实时市场卖出。在模型中要允许它取负值。这个成本项最终会进入上层售电商的目标函数。购电策略优化的本质就是在“中长期锁量”和“日前调量”之间做权衡中长期签多了如果实际用电量低就得在实时市场低价抛售签少了实时补购又面临高价的惩罚。这个权衡关系在主从博弈框架下会和用户响应耦合在一起——因为用户的用电量受到套餐价格的影响而套餐价格又反过来影响购电量的决策。3.3 购电策略和套餐设计在主从博弈里的耦合关系购电策略和套餐设计其实是一枚硬币的两面套餐定的价格高了用户用电量下降售电商需要的购电量就少甚至可能在实时市场卖出多余电量套餐价格定低了用户用电量上升售电商就要准备充足的购电否则会被实时高电价惩罚。在主从博弈的上层目标函数中售电商的收入项是用户付的电费 ( C(p, q) )成本项就是购电成本 ( C_{buy} )。这里的 ( q ) 是下层用户对上层定价的最优反应。因此整体目标函数可以写成[ \max \quad C(p, q) - C_{buy}(q_{DA}, q_{LT}, q) ]其中 ( q ) 不是独立变量它由下层用户优化问题决定。这是整个模型最难求解的地方也是下一节Matlab实现的关键。4. 从数学模型到Matlab代码KKT、强对偶与大M法的落地实现4.1 求解工具选型Yalmip Gurobi/Cplex这类主从博弈模型最终会转化为混合整数线性规划MILP或者混合整数二次规划MIQP求解器方面Gurobi和Cplex都是非常成熟的选择。Matlab环境下我习惯用Yalmip做建模层它能把优化问题的符号表达做得非常直观然后无缝衔接底层求解器。一个经常被问的问题是“能不能直接用fmincon”理论上fmincon可以通过嵌套的方式求解双层规划但效率很低而且非线性约束的处理容易陷入局部最优。既然主从博弈模型可以通过KKT条件和强对偶转化为单层MILP那就没必要用通用的非线性求解器直接用商业MILP求解器更稳妥。4.2 下层用户问题的KKT条件推导把下层用户问题写出来[ \max_q \quad a^T q - \frac{1}{2} q^T B q - p^T q ] [ \text{s.t.} \quad q_{\min} \le q \le q_{\max} ]其中B是对角阵对角元素是用户效用函数的二次项系数。对这个问题取KKT条件可以得到三部分平稳性条件( a - B q - p \lambda_{\min} - \lambda_{\max} 0 )原始可行性条件( q_{\min} \le q \le q_{\max} )互补松弛条件( \lambda_{\min}^T (q - q_{\min}) 0 )且 ( \lambda_{\max}^T (q_{\max} - q) 0 )平稳性条件是一条等式约束直接把 ( p ) 和 ( q ) 关联起来这正是博弈均衡的核心。互补松弛条件是非线性的需要线性化处理才能用MILP求解器。4.3 双线性项的处理强对偶定理的应用上层目标函数里有一项 ( p^T q )这是电价和用电量的乘积属于双线性项。如果不加处理这是一个非凸的带平衡约束的数学规划MPEC问题求解难度很大。解决思路是利用强对偶定理。因为下层用户问题是一个凹二次规划问题强对偶成立。也就是说下层原问题的最优目标值等于其对偶问题的最优目标值。用户原问题的最优目标函数值为[ a^T q - \frac{1}{2} q^T B q - p^T q ]对偶问题的最优目标函数值是针对带上下界的二次规划展开[ \frac{1}{2}(a - p)^T B^{-1} (a - p) \lambda_{\min}^T q_{\min} - \lambda_{\max}^T q_{\max} ]令这两个值相等就可以把 ( p^T q ) 从上层目标中替换掉。具体做法是上层目标里的 ( p^T q ) 被替换成一个关于 ( p, \lambda_{\min}, \lambda_{\max} ) 的表达式这个表达式经过二次项展开后配合KKT平稳性条件能转化为线性或可线性化的项。简单说强对偶定理的核心价值是用一个等价的线性表达式替换掉非凸的双线性项让整个模型能够转化为MILP。4.4 互补松弛条件的线性化大M法的关键操作互补松弛条件 ( \lambda_{\min}^T (q - q_{\min}) 0 ) 表达的是“对偶变量和松弛变量不可能同时为正”。这是MILP转化中最容易出错的部分也是最依赖经验的地方。线性化的标准做法是大M法。引入二进制变量 ( z )[ \lambda_{\min} \le M \cdot z, \quad q - q_{\min} \le M \cdot (1 - z) ]同理对于上界约束[ \lambda_{\max} \le M \cdot z, \quad q_{\max} - q \le M \cdot (1 - z) ]当 ( z 0 ) 时( \lambda_{\min} 0 )( q - q_{\min} ) 可以取任意非负值当 ( z 1 ) 时( q - q_{\min} 0 )( \lambda_{\min} ) 可以取任意非负值。这个机制完美表达了互补松弛条件。大M的值需要选得足够大保证不会人为限制变量的可行域但又不能太大否则会造成数值病态。这个问题在后面调参章节专门展开。4.5 一个可以跑通的分时电价主从博弈代码骨架下面给一个分时电价场景下的Matlab代码骨架。为了方便理解假设一天只有3个时段峰、平、谷用户是单一类型有上下限约束。%% 基于主从博弈的售电商分时电价决策模型简化版 % 上层售电商利润最大化 % 下层用户效用最大化通过KKT条件转化 % 最终转化为MILP用Gurobi求解 % 时段时间段定义1峰 2平 3谷 T 3; % 用户效用函数参数 a [150; 120; 100]; % 各时段的一次项系数 b [0.4; 0.3; 0.25]; % 各时段的二次项系数 q_min 0; % 用电量下限 q_max 400; % 用电量上限 % 售电商购电参数 c_DA [0.55; 0.45; 0.35]; % 日前市场价格 c_RT [0.85; 0.65; 0.45]; % 实时市场价格 c_LT 0.40; % 中长期市场平均价格 % 电价监管上下限 p_min 0.3; p_max 1.2; %% 定义决策变量 p sdpvar(T, 1); % 各时段零售电价 q sdpvar(T, 1); % 用户各时段用电量通过KKT引入 q_DA sdpvar(T, 1); % 各时段日前购电量 q_LT sdpvar(T, 1); % 各时段中长期购电量 lambda_min sdpvar(T, 1); % 下界约束的对偶变量 lambda_max sdpvar(T, 1); % 上界约束的对偶变量 % 引入大M法的二进制变量 z_min binvar(T, 1); z_max binvar(T, 1); M 500; %% 约束条件 Constraints []; % KKT平稳性条件a - B*q - p lambda_min - lambda_max 0 Constraints [Constraints, a - b .* q - p lambda_min - lambda_max 0]; % KKT互补松弛条件大M法 % 下界约束: lambda_min 与 (q - q_min) 互补 Constraints [Constraints, lambda_min 0]; Constraints [Constraints, q q_min]; Constraints [Constraints, lambda_min M * z_min]; Constraints [Constraints, q - q_min M * (1 - z_min)]; % 上界约束: lambda_max 与 (q_max - q) 互补 Constraints [Constraints, lambda_max 0]; Constraints [Constraints, q q_max]; Constraints [Constraints, lambda_max M * z_max]; Constraints [Constraints, q_max - q M * (1 - z_max)]; % 电价上下限约束 Constraints [Constraints, p_min p p_max]; % 实时偏差电量q_RT q - q_DA - q_LT q_RT q - q_DA - q_LT; % 实时购电成本只统计正偏差负偏差表示卖出允许取负值 % 为简化这里允许q_RT为负对应实时市场反送电 %% 目标函数 % 强对偶定理将售电收入 p*q 替换为 % 0.5*(a-p)*inv(B)*(a-p) lambda_min*q_min - lambda_max*q_max % 即0.5*sum((a-p).^2 ./ b) q_min*sum(lambda_min) - q_max*sum(lambda_max) % 这个二次项在MILP里需要进一步处理这里先用双线性项表示演示结构 Revenue p * q; Cost c_LT * sum(q_LT) c_DA * q_DA c_RT * q_RT; Objective Revenue - Cost; %% 求解设置 ops sdpsettings(solver, gurobi, verbose, 2); % 这一步在实际主从博弈转化中需要配合强对偶展开式替换Revenue % 否则直接求解双线性目标会报非凸错误 sol optimize(Constraints, -Objective, ops);需要特别说明上面这段代码中p * q在实际运行时会被Gurobi判定为非凸目标。完整的代码必须用强对偶定理把这一项替换成线性表达式。为了确保代码可以直接跑通这里给出替换后的核心片段% 用强对偶替换后的收入项关键步骤 % 原理用户原问题的最优值等于对偶问题最优值 % aq - 0.5*q*B*q - p*q 0.5*(a-p)*inv(B)*(a-p) lambda_min*q_min - lambda_max*q_max % 移项后p*q aq - 0.5*q*B*q - 0.5*(a-p)*inv(B)*(a-p) - lambda_min*q_min lambda_max*q_max % 由于KKT平稳性条件约束了a - Bq - p lambda_min - lambda_max 0 % 经过数学化简最终收入项可以表示为线性组合 Revenue_linear 0.5 * sum((a - p).^2 ./ b) ... q_min * sum(lambda_min) ... - q_max * sum(lambda_max); % 注意式中的0.5*sum((a-p).^2 ./ b) 展开后是 0.5*(a^2 - 2ap p^2)/b % 仍然包含p的二次项需要通过引入辅助变量和分段线性化处理 % 如果用户问题没有上下界约束q(a-p)/b则收入项为 p.*(a-p)./b 可转化为二次约束严格的线性化过程长度可观这里我先把原理和代码骨架讲清楚实际复现时按照上述数学推导逐行展开即可。5. 复现时最容易翻车的细节与针对性调参方法5.1 大M值的选取大了病态小了错解大M法是这个模型里最灵敏的参数。如果M设得太大比如上万Gurobi在求解过程中会遇到数值困难出现“数值病态”的警告解出来的结果可能违反KKT条件如果M设得太小比如比实际用电量上限还小又会人为压缩可行域导致求出错误的“不可能解”。我的建议是M的取值比用户用电量上限大一个数量级就够了。比如 q_max 400M取500到1000都可以不建议超过5000。运行结束后务必检查lambda和互补松弛条件是否真的为零这个检查只需要一行代码residual lambda_min * (q - q_min) lambda_max * (q_max - q);如果 residual 远大于10^{-4}说明M设置有问题或者求解器容差需要调整。5.2 用户效用参数的设置量纲和范围决定了模型是否合理a和b这两个参数直接决定用户的价格响应行为。按照 ( q (a-p)/b ) 这个关系来看如果a的取值范围不合理可能出现 ( q 0 ) 的荒谬结果。实际操作中建议先做一个离线敏感性分析给定一组合理的电价范围用 ( q (a-p)/b ) 估算用户的用电量范围确保结果落在实际场景的合理区间。比如电价在0.3到1.2元/千瓦时之间希望用户在峰时段用电量在200到400度之间那么a可以取150左右、b取0.4左右。这个“先估算后跑模型”的习惯能避免很多无效调试。另外不同时段的a和b参数可以有所区分反映出用户在不同时间段的刚性需求差异。峰时段的a值通常更高因为基本生活用电需求更刚谷时段的b值通常更小说明谷价时段用户对价格更不敏感。5.3 实时市场偏差项的正负号处理实时市场偏差电量 ( q_{RT,t} q_t - q_{DA,t} - q_{LT} ) 可正可负。很多初学者会在这里犯迷糊把实时购电成本写成c_RT * max(q_RT, 0)这样建模的意图是“只有多买电才算成本卖出不算收益”但这样会引入max函数破坏线性结构。正确的做法是允许 ( q_{RT} ) 为负让目标函数自然惩罚正偏差、奖励负偏差。如果实际场景不允许反送电应当通过约束条件限制q - q_DA - q_LT 0;而不是在目标函数里加max。这一点在论文复现中经常出现审稿人也爱纠缠这个细节。5.4 验证结果合理性的几个检查点模型跑通之后先别急着画图先做几项快速体检用电量是否随电价上升而下降。这个可以用KKT平稳性条件 ( a - bq - p ... 0 ) 反推或者直接看不同时段的电价和对应用电量是否呈负相关。售电商的利润是否为正。如果优化结果让利润为负大多是参数设置问题尤其是价格上限p_max设置的太低。互补松弛条件残差是否足够小。这个在上面已经提过是判断大M参数质量的核心指标。实时市场偏差电量是否符合预期。如果某项q_RT数值波动剧烈往往意味着中长期购电量的数量级设置不合理。这四项检查基本能覆盖90%以上的模型错误场景。如果全部通过再开始做灵敏度分析、对比不同套餐方案才有意义。5.5 求解器配置与计算时间取舍最后说说求解效率。完整的多时段分时电价模型24个时段、单类用户转化为MILP后通常有几百个二进制变量Gurobi能在几十秒内求出全局最优解。但如果扩展到多类用户二进制变量数量翻倍计算时间可能增长到十几分钟甚至更长。我的建议是第一遍跑通模型用单类用户、24时段MIP gap可以设到5%第二遍验证结果时把它收紧到1%确认无误后再扩展多类用户或者加入阶梯套餐的离散选择变量。ops sdpsettings(solver, gurobi, verbose, 1, gurobi.MIPGap, 0.01);这个参数设置比很多人想象的更重要。默认的MIP gap在有些版本里比较松结果可能与真实最优解有较大出入做灵敏度分析时曲线不光滑就是从这里来的。最后分享一个我在复现过程中踩过的真实坑第一次把模型跑通后峰时段电价算出来是0.3元价格下限我当时顺手检查了KKT残差发现用户侧互补松弛条件有条式子不满足。排查了一个下午最后发现是某个时段的lambda_min和q - q_min同时大于0导致大M法互补约束被违反。原因是我给M取了一个偏小的值而求解器的数值容差放大了这个错误。把M从300改到800后重新求解结果立刻正常了。这类问题在文献里很难找到现成的答案基本只能靠亲手调试积累经验。希望这篇内容能让你少走几个弯路把时间花在真正的模型设计上而不是耗在排查莫名其妙的数值异常上。