搞电动汽车充电调度的同学一定遇到过这种尴尬：运营商想靠高电价多赚钱，用户看到电价高了就磨蹭到半夜再充，两边对着干，最后要么充电桩利用率惨淡，要么负荷尖峰把小区变压器压得喘不过气。我之前用MATLAB做了一套基于主从博弈的智能小区充电管理代码，核心思路就是让代理商当"领导者"先出动态电价，电动汽车用户当"跟随者"再决定怎么充，把这种对抗关系用Stackelberg博弈变成可求解的数学模型。这个项目从模型推导、Yalmip建模到CPLEX求解全部落地，代码量不大，但里面的推导细节和调参经验挺值得展开说说，就想写出来分享给同样做电力系统仿真的朋友。
内容方向偏智能电网、电动汽车充电管理方向的仿真科研，也适合正在学双层优化和博弈论应用的人参考。如果你手头要复现类似的定价策略，或者刚接触Yalmip和CPLEX，这篇应该能帮你少走不少弯路。
1. 这个项目在解决什么问题
1.1 小区充电博弈的痛点
先说你实际会遇到的场景。一个普通居民小区，晚上六七点开始，住户陆续回家把车插上充电桩。如果不加任何管理，所有人都在这个时间段抢功率，小区变压器容量很快就到顶了，甚至可能出现跳闸风险。而充电代理商（可以理解成运营充电桩的物业公司或运营商）不希望变压器扩容，那样成本太高，更希望用户错峰充电，把负荷平摊到半夜。
问题在于用户没有动力配合。用户在意的是充电费用和第二天早上电量够不够，哪个时段便宜就哪个时段充，并不关心变压器是否过载。代理商如果只用一个固定电价，那所有用户都会涌向同一个低价时段，形成新的负荷尖峰；如果电价定太高，用户又可能直接不在小区充了，代理商赚不到钱。
这本质上是一个利益冲突问题。代理商有价格制定权，用户有充电量选择权，双方各自优化自己的目标，互相影响。常规的单层优化没法处理这种互动关系，因为它默认所有主体听统一调度。现实中用户不会完全听调度，你让他十点再充，他偏要七点就插上。所以需要一个能刻画"你定规则我响应"的建模框架，主从博弈就是这个场景下的首选模型。
1.2 为什么是主从博弈
主从博弈也叫Stackelberg博弈，核心特征是一方先行动，另一方观察到行动后再响应。先行动的一方叫领导者，后响应的一方叫跟随者。在这个项目里，代理商是领导者，它先公布未来24小时每个时段的充电电价；电动汽车用户是跟随者，看到电价后决定自己每个时段充多少功率。
容易混淆的是，这不是一个动态视频博弈，而是一个"同时决策但带有先后层级"的静态优化问题。代理商在制定电价时就要预判用户的反应，它心里有数：电价提高一毛钱，夜间时段的充电需求会下降多少；电价压低五分钱，高峰时段的负荷能转移多少。借助下层用户的优化模型，可以把这种"预判"数学化为下层问题的最优解函数，然后代回上层目标函数。这样一来，代理商的定价决策就不再是拍脑袋，而是建立在对用户理性响应的准确建模之上。
为什么不用传统的双层规划？其实主从博弈在数学形式上就是双层优化，只是博弈论给了它更清晰的经济学解释：上层目标最大化代理商利润，下层目标最小化用户充电费用。两者通过价格和电量耦合。用双层优化的角度理解更容易上手，用博弈论的角度理解更容易设计均衡分析，二者不冲突。
1.3 项目边界与适用场景
这套代码针对的是"单代理商、多用户、24小时静态定价"的典型场景。代理商管理一个小区内的若干充电桩，所有用户都有夜间充电需求，代理商统一向电网购电再转售给用户。时间粒度设为1小时，视作24个决策时段。用户车辆参数（电池容量、初始SOC、到达离开时间）已知且各异，但代理商不做实时滚动优化，而是提前一天制定全天电价。
它适合用来做策略验证和学术研究，比如分析不同用户规模下均衡价格的变化、对比博弈定价与分时电价的削峰填谷效果、观察价格上限对利润和社会福利的影响。如果要做在线实时控制，还需要把模型改成模型预测控制的滚动框架，交易间隔从1小时压缩到15分钟甚至5分钟，这是后续扩展的方向。
2. 主从博弈模型设计与数学推导
2.1 系统架构与前提假设
建模之前必须先定清楚系统里有哪些主体、哪些变量、哪些参数。这个项目里有两个核心主体：充电代理商和电动汽车用户。代理商从电网以购电价买电，同时向用户收取充电服务费，最终用户支付的总价格就是代理商制定的动态电价，购电价则参考当地分时电价结构，在代码里作为已知参数输入。
时间轴设置为1天24个时段，每个时段长度为1小时。这个粒度对小区慢充场景完全够用，因为慢充本身持续几个小时，1小时的决策精度已经能反映错峰转移的效果。如果后续想做得更细，改成96个15分钟时段，模型结构基本不用动，只是变量数量翻四倍。
基本的简化和假设有四条。第一，忽略充电桩的功率损耗和电池老化成本，充电效率统一设为1。第二，每辆车的充电功率连续可调，范围从0到额定功率，模拟慢充桩的调功特性。第三，用户必须满足最低充电需求，也就是离开时SOC要达到自己设定的目标值，这是硬约束。第四，传输线路损耗和电压约束暂时不考虑，只把变压器容量作为总负荷上限。
这些假设看似简化得很厉害，但对策略层研究是合理的。充电管理研究的核心是电价引导下的负荷转移规律，不是配电网络的潮流细节。如果谁一上来就把三相不平衡、谐波、线损全塞进去，模型复杂到根本没法分析博弈关系，得不偿失。
2.2 上层代理商动态定价模型
上层模型里，代理商需要决策的就是24个时段的充电价格，记为向量，每个元素是对应时段的单位电价，单位是元/kWh。代理商的目标是最大化全天利润，利润等于售电收入减去购电成本。
售电收入是所有用户在各个时段充电量的加权和，权重就是对应时段价格。购电成本是总充电负荷乘以对应时段的购电价。这里的总负荷不是自由变量，而是下层用户优化后的结果，由所有用户最优充电策略叠加得到。也就是说，代理商的利润函数里，隐含了"价格影响用户行为"这一层关系。
约束条件方面，代理商不能随意定价。一方面要设置价格上下限，避免极端高价损害用户利益，这在现实中也对应政府对充电服务费的监管；另一方面可以加一个利润率约束，比如要求平均售价不低于购电成本的某个比例，确保代理商可持续运营。实际操作中我倾向于只设价格上下限，其余约束不加，因为约束越多越容易出现下层可行性被破坏的问题。
这里还要解释一个关键点：为什么上层目标不是简单地在低谷时段降价、高峰时段涨价？因为用户对高价的响应不是无限度的，如果某个时段价格高到用户宁可少充电也不充，代理商的总收入反而下降。所以定价问题本质上是寻找一个"价格-需求"平衡点，这正是建模的价值所在。
2.3 下层用户充电决策模型
下层模型中，每个用户是一个独立的小型优化问题。用户的决策变量是自己在24个时段内的充电功率曲线。用户的目标很朴素：在自己接入充电桩的时间窗口内，最小化总充电费用。
每个用户的约束条件包括：充电功率不能超过充电桩额定功率；不在接入时间窗口内的充电功率强制为0；离开时必须达到目标SOC对应的充电量需求。用数学语言描述，就是把充电功率限制在0和额定功率之间，并要求累计充电量不小于目标充电量。
需要注意一个看似矛盾的点：目标函数是最小化费用，而目标充电量约束是"不小于"某个值，所以用户理论上可能超额充电，但超额充电只会增加费用，不会带来任何收益，因此最优解一定会精确停留在目标充电量上，不会多充。这样"不小于"约束在实际解中等式自动成立。
每个用户的到达时间、离开时间、电池容量、初始SOC、目标SOC、最大充电功率都是独立设定的。比如有的车电池容量40kWh，初始SOC 30%，目标SOC 85%，那么最少需要补充22kWh电量；另一辆车可能只需要10kWh，它的充电灵活性就大得多，可以在更便宜的时段轻松完成充电。这种用户之间的异质性是后面分析价格弹性差异的基础。
2.4 双层问题如何转化为单层可解模型
这个双层问题不能直接丢给求解器，因为上下层各自有独立的目标函数和决策变量，嵌套在同一个问题里。最经典的求解思路是把下层优化问题用KKT条件替换，变成一个数学上等价的单层约束系统。
下层用户优化问题是一个线性规划，目标函数和约束都是线性的，因此满足强对偶条件，KKT条件是原问题最优解的充分必要条件。KKT条件包含四类内容：拉格朗日函数对决策变量的导数为零（即平稳性条件）；原始可行性条件；对偶变量非负条件；互补松弛条件。
平稳性条件把上层电价变量和下层拉格朗日乘子关联起来，这是上下层耦合的核心。互补松弛条件处理的是"如果某个不等式约束不起作用，则对应乘子必须为0"这类逻辑关系。由于互补松弛条件是非线性的（乘子乘以约束边界等于0），需要引入二进制辅助变量并用大M法线性化。
最终得到的单层问题是一个混合整数线性规划（MILP），包含连续变量（价格、充电功率、拉格朗日乘子）和二进制变量（互补松弛条件指示）。求解这个MILP就能得到Stackelberg均衡解。代价是模型规模变大，但对几十辆车的规模，CPLEX完全能胜任。
3. MATLAB代码实现与调试细节
3.1 环境准备与工具箱选择
MATLAB代码的核心依赖是Yalmip优化建模工具箱和商用求解器。Yalmip的作用是把优化模型用接近数学公式的方式写出来，避免手写大规模的稀疏矩阵，大幅降低建模出错的概率。求解器方面，推荐CPLEX或Gurobi，两者都能高效求解混合整数线性规划。
我用的是R2020a版本的MATLAB配CPLEX 12.10。Yalmip安装很简单，把解压后的文件夹加入MATLAB路径即可。启动时可以用yalmiptest命令检查求解器是否被正确识别。常见的坑是装了Yalmip但MATLAB找不到CPLEX接口，大概率是因为求解器安装时间在Yalmip之后，重新启动MATLAB或者重跑savepath就能解决。
如果没有商用求解器，也可以先用MATLAB自带的intlinprog凑合跑小规模算例。Yalmip对intlinprog是支持的，但求解速度比CPLEX慢很多，而且数值稳定性要差一些，几十辆车的算例可能会跑得很痛苦。做严肃的科研复现，还是建议直接上CPLEX或Gurobi。
3.2 代码结构与核心变量定义
这个项目我把代码拆成几个文件：主脚本负责整体流程，参数文件集中定义所有车辆和电网数据，模型构建函数负责在Yalmip中创建优化模型，绘图函数负责结果可视化。这样拆的好处是改参数时不用翻主脚本，跑敏感性分析也方便。
核心变量的定义部分，价格用p = sdpvar(1, T)表示，其中T是时段数24。每辆车的充电功率因为接入时段不同，建议用元胞数组：P{i} = sdpvar(1, T)，这样每辆车都有完整的24维变量向量。虽然各车有大量未接入时段变量最终被约束为0，但用元胞数组写代码最直观，后续优化性能时再改成只对有效时段建变量。
拉格朗日乘子也对应定义：下层约束中充电功率下限对应的乘子、充电功率上限对应的乘子、目标充电量约束对应的乘子，分别用sdpvar创建。特别注意这些乘子必须加非负约束，否则KKT条件就错了。
二进制变量用于互补松弛条件的线性化，数量大约等于车辆数乘以时段数再乘以2，再加上车辆数。举例：50辆车、24个时段，光功率上界和下界的互补松弛就要引入2400个二进制变量，这还能扛；如果加到200辆车，MILP规模会显著上升，需要走聚合或分解的路线。
3.3 关键代码段的实现思路
先展示KKT平稳性条件的核心代码。对每个时段、每辆车，价格乘以时间段长度与充电功率、对应乘子之间的关系要满足平稳性条件。用Yalmip写的话就是：
% 平稳性条件：c - lam_low + lam_up - mu = 0 % c = p(t) * dt for i = 1:N for t = 1:T Model = [Model, p(t)*dt - lam_low{i}(t) + lam_up{i}(t) - mu(i) == 0]; end end这段代码的量级很小，但它是上下层耦合的关键。lam_low和lam_up分别对应充电功率下界和上界约束的乘子，mu对应目标电量约束的乘子，三者都和该时段的充电价格直接挂钩。
再来看互补松弛条件怎么用大M法线性化。以充电功率上界约束为例，原始条件是：如果某个时段的充电功率没有达到额定功率上限，则上界对应的乘子必须为0。用二进制变量和两个不等式来模拟这个逻辑：
% 互补松弛：lam_up * (P - Pmax) = 0 % 用 z_up 变量指示是否达到上限 M_power = 100; % 功率类大M for i = 1:N for t = 1:T Model = [Model, lam_up{i}(t) >= 0]; Model = [Model, lam_up{i}(t) <= M_power * z_up{i}(t)]; Model = [Model, P{i}(t) - Pmax(i) <= M_power * (1 - z_up{i}(t))]; Model = [Model, P{i}(t) - Pmax(i) >= 0]; % 这里其实是 P 的上界，作为原始约束一部分 end end这里有个容易写错的地方：P{i}(t) - Pmax(i) >= 0其实是上界约束本身，它不能只说一半。在单层化之后，所有原始约束都必须保留。所以写代码时要把原始下界约束、原始上界约束、原始电量约束都加进去，再加上KKT条件。少任何一条原始约束，解都可能偏离真实可行域。
目标函数是整个模型最需要花心思的地方。如果直接把售电收入写成价格乘充电功率，目标函数里会出现双线性项，问题变成非凸的，CPLEX没法保证全局最优。这里用强对偶定理做替换：因为下层是线性规划，原问题最优值等于对偶问题最优值，于是可以把下层目标函数中的双线性项替换成对偶变量表达式。具体代码如下：
% 目标函数：上层利润最大化 % 收入项用强对偶替换后的形式表达 revenue_expr = 0; for i = 1:N revenue_expr = revenue_expr + Ereq(i) * mu(i) - Pmax(i) * sum(lam_up{i}); end cost_expr = sum(c_buy .* sum_power) * dt; Objective = revenue_expr - cost_expr;这里的sum_power是所有车辆各时段功率之和，c_buy是分时购电价。通过这个替换，目标函数完全线性化，模型变成标准MILP，求解器可以稳定求出全局最优解。很多初学双层优化的人卡在这一步，模型建好却跑不动，就是因为没有意识到双线性项需要先处理。
3.4 参数初始化怎么设最合理
参数初始化直接决定结果合不合理，也决定模型是否有可行解。我常用的参数设置如下：车辆数量默认50辆，电池容量在40到60kWh之间均匀分布，最大充电功率统一取7kW，对应单相慢充桩。每辆车的初始SOC取20%到50%均匀分布，目标SOC统一设90%。到达时间集中在18:00到21:00随机生成，离开时间集中在7:00到9:00随机生成。
购电价参考典型分时电价结构：白天高峰时段（比如10:00-15:00和18:00-21:00）电价较高，夜间低谷时段（比如23:00-次日6:00）电价较低。代理商售电价格的上下限我通常设为购电价的0.9倍到1.3倍，这样既能保证代理商有利润空间，又不至于让用户费用高得离谱。
有一个特别需要注意的参数：用户电量需求不能设得太极端。比如初始SOC只有10%却要求离开时达到95%，这辆车需要充50多度电，在7kW额定功率下至少要充8个小时，如果它的接入时间窗口只有6个小时，模型直接无解。所以初始化时最好写一段校验逻辑，自动检查每辆车的最大可充电量是否满足目标需求，不满足就调整目标SOC，而不是让求解器返回infeasible。
4. 仿真结果分析与对比
4.1 均衡结果的判读
模型求解完成后，第一个要看的是价格曲线和总负荷曲线的对应关系。正常的结果应该是：夜间低谷时段电价相对较低，对应充电负荷的集中区域；高峰时段电价较高，对应充电负荷的谷底。代理商通过这种价差引导用户把充电需求从晚间高峰转移到凌晨低谷。
再看代理商的利润构成。收入端里，目标电量对应的乘子乘以需求的累加项是大头，而功率上界对应的乘子通常只在充电达到上限的时段才起作用。如果某个时段的乘子处处为0，说明该时段的功率约束没有收紧，充电功率没有触顶，这个时段还有继续提高价格的余地。
还要检查用户的充电策略是否合理。理论上每个用户都会在自己接入窗口内选择电价最低的时段充满需求电量。如果某个用户最终计费时段分布完全不符合这个规律，大概率是KKT条件或者互补松弛条件写错了。我调试时经常把某辆车的功率曲线单独画出来，和电价曲线叠在同一张图上，一眼就能看出用户是否理性选择低价时段充电。
4.2 用户规模变化时的表现
模拟用户规模从30辆逐步增加到120辆，观察均衡价格和负荷曲线的变化，会发现很有意思的规律。用户数量少时，单辆车的充电需求对总负荷贡献大，代理商定价空间相对受限，不敢定太高，因为少一个用户就会显著影响收入。用户规模增大后，个体用户的议价能力变弱，代理商可以小幅抬高整体价格水平，利润增长曲线比车辆数量增长得更快。
另一个明显变化是负荷集中度。用户数量多时，每个用户的充电窗口叠加在一起，几乎覆盖全天，代理商可以通过价格信号把负荷平滑地分散到多个时段，变压器峰值负载率反而可能下降。但如果用户到达时间高度集中，比如全部挤在19:00到20:00之间接入，即使数量再多，可转移的时段也有限，削峰效果会打折扣。
这个结果对实际运营有直接意义：小区充电运营商的定价策略不能只看自己这一侧的成本，还要看用户的到达规律和充电需求刚性。用户越集中、需求越刚性，价格的调节空间越小；用户越分散、需求越弹性，价格工具的威力越大。
4.3 与固定电价和分时电价的对比
我做了三组对比实验：固定电价全天不变、常规分时电价、本文的主从博弈动态定价。固定电价下，用户全都在自己接入后的最低价时段充电，但因为全天价格一样，大家没有错峰动机，结果就是晚间接入高峰时段的负荷冲得很高，变压器很危险。
分时电价好一些，夜间低谷电价能引导部分错峰，但问题是分时电价是静态的，它不会考虑当前小区实际负荷水平。某天如果大量车辆都在低谷时段同时充电，甚至低谷时段本身也可能出现过载。分时电价只管"时段的价差"，不管"小区变压器当前的剩余容量"。
主从博弈动态定价的优势在于价格直接响应所有用户的聚合需求。代理商在知道所有车辆参数和充电需求的情况下，提前算出一个价格曲线，这个曲线天然考虑了变压器容量约束和用户响应行为。三组实验对比下来，博弈定价方案的负荷峰谷差最小，同时代理商利润也没有比理想的分时电价场景低太多。也就是说，用户视角看，充电费用增幅有限；运营商视角看，削峰效果好、运营压力小；电网视角看，负荷曲线更友好。三方利益达到一个相对平衡的状态。
5. 踩过的坑与排查经验
5.1 不可行解与求解器报错
跑MILP最怕的就是求解器直接返回"Problem is infeasible"（不可行）。排查不可行问题，第一步不是改约束，而是把模型拆开：先固定上层价格，单独求解下层用户模型，看下层模型本身是否有可行解。如果单车跑就有问题，那就是参数初始化出错，比如前面说的充电需求超过可充电时长；如果单车没问题，那就是上下层耦合时加了过多的价格约束，把下层可行域挤没了。
另一个常见原因是互补松弛条件里的大M取值太小，导致二进制变量的逻辑约束把原本可行的方案剪掉了。遇到这种问题可以先把互补松弛条件整体注释掉，只保留原始约束和KKT平稳性条件，再求解一次。如果能得到合理的连续松弛解，说明瓶颈就在互补松弛条件的建模细节上。
5.2 大M参数的选择技巧
大M法的核心是选一个足够大的常数，让二进制变量取0时对应的不等式自然满足，取1时强制对应等式成立。M太小会误伤可行域，M太大则会导致数值病态，求解器在判定整数变量和连续变量关系时可能出现精度问题。
功率相关约束的大M不需要很大，因为充电功率本身有上界，取一个比最大充电功率大10倍左右的值就行。我在50辆车、7kW额定功率的场景下用100完全够用。电量相关约束的大M与电池容量、目标SOC相关，我习惯取1000，因为充电电量也就是几十千瓦时量级，1000给了足够松弛度，没有出现过数值问题。
一个实用经验是：M能用精确边界算出来，就不要偷懒取固定大数。比如功率上界互补松弛条件里的M可以写成Pmax(i)本身，因为它是该变量理论上限。精确边界能提升求解速度，还能减少数值警告。
5.3 目标函数里的双线性项怎么消
这点前面提过，实际操作时是让很多初学者崩溃的地方。初始建模时很自然地写出"价格乘以功率"作为售电收入，丢给Yalmip后发现报错说目标函数非凸，或者CPLEX直接警告无法保证最优解。这个双线性项必须处理。
处理方式是强对偶替换，前提是下层优化是线性规划。因为线性规划强对偶成立，原问题目标最优值等于对偶问题目标最优值，于是把下层目标函数中用连续变量和价格构成的乘积项替换成对偶变量表达式。替换之后整个上层目标函数变成线性函数，模型从非凸问题变成MILP。
这里要特别提醒：如果下层模型里添加了二次惩罚项，比如为了模拟用户舒适度而加入充电功率平方项，下层就不是线性规划而是二次规划，强对偶替换依然可以部分使用，但对偶问题形式会更复杂。为了保证代码简单和求解稳定，我的项目把下层目标保持为纯线性。
5.4 大规模算例提速方案
当车辆数量超过100辆时，MILP的求解时间可能从几十秒飙升到几分钟甚至更久。提速方案有很多，我最常用的是用户聚合。先把用户按电池容量、初始SOC、到达离开时间做聚类，每一类用户用一个聚合充电变量表示，最后再把单辆车的功率约束改写成聚合功率约束，相当于在数学上把多个相同特征的用户合并成一个"代表用户"。
聚合后模型变量数量大幅下降，求解速度可以提升一个数量级，但代价是用户异质性刻画变差。一个折中方案是保持每个用户的原始模型，但只把每辆车接入窗口外的变量不创建，而不是先创建再约束为0。这样相当于把很多一直为0的变量从模型里删掉，规模也能显著减小。
另外可以试试给求解器设置合理的运行参数。比如CPLEX的mipgap从默认值放松到1%，求解时间经常能缩短一半以上，而结果差异在可接受范围内。科研场景下如果只需观察趋势而非追求极致精度，这个方法性价比很高。
最后再分享一个我调试时常用的小技巧：先把模型里所有整数变量清空，把它当成线性规划求解一次，把松弛解打印出来。这个连续松弛解往往能暴露模型里约束方向写反、变量符号定义错误这类低级问题，而且排查速度比直接看复杂MILP结果快得多。等你把连续松弛调顺了，再恢复整数变量求精确解，思路会清晰很多。这套项目如果后续想继续扩展，可以考虑加入V2G双向充电、屋顶光伏和储能系统，让代理商在定价时同时调度放电功率和储能充放策略，博弈框架不变，只是下层用户的决策空间更丰富。先把这个基础模型吃透，后面的扩展都是水到渠成的事。
