基于主从博弈的售电商零售套餐与多级市场购电策略Matlab复现

电力市场方向的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^contractQ_t^day是日前市场购电量价格c_t^dayQ_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小时调度周期为例时段数T24。用户数量初始设为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再配合随机种子的固定结果差异很小速度却快了很多。如果你也要跑大规模参数建议先从S50开始调通再逐步增加。另外中长期合约价格我设为一段固定均值只在部分时段允许签约模拟实际合约约束。合约电量上限按用户总负荷的某个比例设置比如最高能覆盖总电量的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 调试小技巧与复现建议调试这类模型时最忌讳一上来就调全规模参数。我的调试顺序是先设N2、T4、S1把模型跑通。T4相当于只模拟凌晨、上午、下午、晚上四个时段N2代表两类用户S1就是确定性场景跳过所有随机因素。这样模型体量小任何错误都能从约束或目标里快速定位。跑通确定性案例后再逐步增加T、N和S。另外一个非常实用的技巧单独验证用户下层问题是否正确。先固定一组套餐价格单独求解一个用户的优化问题看看最优电量是否等于解析解q(a-p)/b。如果解析解和数值解对不上说明用户问题有约束写错或者双线性项加错地方。这个检查比直接跑整个博弈靠谱得多因为下层问题错漏往往会传导到最终的MPEC里让人查无可查。调试时可以把所有约束目标写成文本导出用稀疏矩阵可视化查看非零元分布能快速找到异常链接。我见过很多次所谓“模型不可行”其实就是某个时段电量平衡约束里漏乘了一个用户数量缩放系数导致供求永远对不齐。最后再分享一个实用建议这类复现项目不要执着于“参数和论文一模一样”。顶刊论文里的数值结果往往经过了大量标定你复现时只需要抓到核心模型和趋势。先把机理层跑通再根据你自己的数据调出合理曲线这样写论文或做汇报时更有说服力。我接下来准备把多类用户对套餐的自选择行为用更严格的离散选择模型表达再对比一下和当前连续近似的差距这个方向有进展的话会继续更新。