简介:面向智能电网中的电动汽车充电管理问题,这份Matlab程序(基于Yalmip和Cplex求解器实现)复现了《基于主从博弈的智能小区代理商定价策略及电动汽车充电管理》一文的核心模型,适合电力系统、优化调度方向的研究生及工程师借鉴学习。该模型将代理商与电动汽车车主各自追求利益最大化的行为构造为主从博弈,并利用Karush-Kuhn-Tucker最优性条件和线性规划对偶定理,将双层博弈转化为混合整数线性规划问题求解,最终获得全局最优的定价策略,对研究电价响应下的电动汽车需求侧管理具有直接参考价值。压缩包为zip格式,大小5.16MB,程序主体为Matlab源码;已有929人学习。读者可从中梳理主从博弈建模到混合整数线性规划求解的完整实现路径,掌握代理商店内定价与电动汽车充电调度的建模思路,并可通过修改参数分析不同电价场景下的负荷变化。
1. 一个 zip 里的 Stackelberg 博弈:EV 充电管理为什么先定价再响应
晚上 7 点,你运营的充电站门口停了 10 台等着充电的车。拍一个固定电价最简单,但价格定高了,理性车主只充“刚需”那几度电,站里的营收反而上不去;定低了,你连购电成本都收不回来。EV 充电管理里的 Stackelberg Game(斯坦科尔伯格博弈)就是处理这种“谁先出牌”的问题:充电站作为领导者先发布电价,车主作为跟随者根据电价决定充多少。标题里这个 zip 包做的正是这样一件事——把上层定价和下层充电量响应建成一个双层优化模型,再用 KKT 条件压成单层问题交给求解器解出均衡。它解决的是充电价格怎么定、车主响应怎么建模、均衡收益怎么算这一整条链路。适合做 V2G 调度、充电站运营、需求响应和电力市场方向的人直接拿来当起点,也适合想搞懂双层规划怎么落地的算法工程师。
2. 先把博弈模型立住:充电站领导者、EV 跟随者与双层数学结构
2.1 领导者和跟随者分别优化什么:决策变量、目标函数与约束
Stackelberg 博弈和普通 Nash 博弈最大的区别在“先后”。充电站先公布电价 p,每台 EV 看到 p 之后再去优化自己的充电量 x_i。EV 的决策是在充电站的价格之后做出的,所以充电站必须“预知”车主会怎么响应,这就形成了上层套下层的结构。
上层是充电站(领导者),决策变量就是充电电价 p。它的收益是卖电收入减去购电成本:
收益 = (p - c0) * Σ x_i其中 c0 是充电站从电网买电的度电成本,x_i 是第 i 台车的充电量。注意这里有个隐含前提:充电站不能直接命令每台车“你必须充多少”,它只能通过价格去引导。如果它能直接调度每台车,那就是集中式优化,不是博弈了——这一点到第 6 章验证时会再对比。
下层是每台 EV(跟随者),决策变量是充电量 x_i。车主有自己的效用函数,学术界最常用的是二次型:
效用_i = a_i * x_i - 0.5 * b_i * x_i^2 - p * x_ia_i 可以理解为车主的“支付意愿”,代表他这趟充电有多迫切;b_i 是效用饱和斜率,充电量越接近满电,额外再充一度电带来的边际满足感越低。车主要最大化自己的效用,约束只有两个:
0 <= x_i <= Pmax_iPmax_i 是这台车当前允许的最大充电功率,由电池 SOC、充电桩功率、车载充电机共同决定。
把上下层放在一起看,就得到标准的双层规划:上层先定 p,下层根据 p 求出最优 x_i(p),上层再拿这个响应去评估自己的收益。正因为下层问题有解析结构,这个博弈才能往下推进。
2.2 把双层问题压成单层:KKT 条件与大 M 互补线性化
双层问题不能直接扔给常规求解器,因为变量 p 同时出现在上下层。一类代码包的做法是“枚举价格 + 反复解下层”,只适合小算例;工程上更通用的是把下层问题用 KKT 条件替换,压成一个单层问题。
下层是凸优化(二次目标加线性约束),所以 KKT 条件是充要条件。对第 i 台车的下层问题写 Lagrange 函数,得到最优性条件:
a_i - b_i * x_i - p + λ_i - μ_i = 0其中 λ_i 对应 x_i >= 0 的对偶变量,μ_i 对应 x_i <= Pmax_i 的对偶变量。除了这条平稳条件,还要带上互补松弛:
λ_i >= 0, λ_i * x_i = 0 μ_i >= 0, μ_i * (Pmax_i - x_i) = 0互补条件是非线性的,不能直接进 MIP。常见做法是引入二元变量 u_i、v_i 和大 M 常数做线性化,这一点是这一整类代码包的核心,也是后面最容易翻车的地方。另一种特例:当下层只有单个变量且没有耦合约束时,反应函数可以直接解析写出来:
x_i(p) = clamp((a_i - p) / b_i, 0, Pmax_i)但这个解析式只在“每台车之间互不影响”时才成立。一旦有配网容量、变压器上限这类耦合约束,就必须老老实实走 KKT 路线。
2.3 最小算例的参数表:价格、效用斜率、功率上限与量纲
写代码之前,先定一套能跑通的最小参数。下面这张表是后续 MATLAB 算例的基础,所有量纲统一成 kW、元/kWh:
| 参数 | 符号 | 取值 | 说明 |
|---|---|---|---|
| EV 数量 | n | 5 | 最小算例只放 5 台车 |
| 支付意愿 | a_i | 12 ± 2 | 车主愿意承受的价位的上限基准 |
| 效用斜率 | b_i | 0.4 ± 0.1 | 决定响应曲线的陡峭程度 |
| 最大充电功率 | Pmax_i | 10 ± 3 kW | 每台车的充电上限 |
| 购电成本 | c0 | 5 元/kWh | 充电站从电网侧买电的成本 |
| 价格下限 | pmin | 0.5 元/kWh | 防止负电价 |
| 价格上限 | pmax | 15 元/kWh | 防无界的关键参数 |
这里最容易埋雷的是量纲。很多论文里 p 是 0.8 元/kWh,但 a 写的是 12,两者一除反应函数直接出现负值。要么全用元/kWh,要么全用元/MWh,绝不能混用。第 5 章我会把这条单独拎出来讲。
3. 用 MATLAB + YALMIP 跑通最小算例:KKT 转换与 MIQP 代码全解
3.1 从 zip 包到工作区:目录结构、数据文件和求解器准备
这类代码 zip 解压后通常是一组.m文件加一两个数据或说明文档,常见结构是main.m(主算例)、solve_.m(求解封装)、data_case.m(参数定义)。先把所有文件放到同一个工作目录,然后确认三件事:
第一,安装 YALMIP,并把 Gurobi 或 Cplex 配到 MATLAB 路径里。本算例的求解器选择是 Gurobi,因为目标函数是双线性的,后面要开 NonConvex 参数。第二,运行yalmiptest确认 YALMIP 能正常调用求解器。第三,注意解压路径里不要带中文和空格,MATLAB 对路径里的中文支持时好时坏,这是复现阶段最低级也最常见的翻车点。
如果只装了免费的求解器,比如linprog那种,MIQP 是跑不了的。可以暂时把目标函数简化为线性扫描,或者装 Gurobi 学术版。别在环境上浪费时间,这个模型的价值在博弈结构本身,不在求解器选型。
3.2 完整可跑代码:下层 KKT 一次性写进约束
下面这段代码完整实现了 2.2 节的 KKT 转换。核心思想是:把下层每台车的 KKT 条件作为约束,和上层的价格变量 p、目标函数放同一个模型里,交给 Gurobi 做 MIQP。
%% 参数定义 clear; clc; n = 5; % EV 数量 rng(1); % 固定随机种子,方便复现 a = 12 + 2 * rand(1, n); % 车主的支付意愿 b = 0.4 + 0.1 * rand(1, n); % 效用曲线斜率 Pmax = 10 + 3 * rand(1, n); % 每台车最大充电功率(kW) c0 = 5; % 充电站购电成本(元/kWh) pmin = 0.5; pmax = 15; % 价格上界,防止无界 %% 决策变量 p = sdpvar(1, 1); % 上层决策:充电电价 x = sdpvar(1, n); % 下层决策:每台车充电量 lambda = sdpvar(1, n); % 对偶变量,对应 x >= 0 mu = sdpvar(1, n); % 对偶变量,对应 x <= Pmax u = binvar(1, n); % u(i)=1 表示 x(i) 不在下边界 v = binvar(1, n); % v(i)=1 表示 x(i) 不在上边界 M = 30; % 大 M 常数,取值原则见 5.4 %% 约束:下层 KKT 条件 + 上层价格边界 Cons = []; for i = 1:n Cons = [Cons, x(i) >= 0]; Cons = [Cons, x(i) <= Pmax(i)]; % KKT 平稳条件 Cons = [Cons, a(i) - b(i)*x(i) - p + lambda(i) - mu(i) == 0]; % 互补条件线性化:x=0 或 x=Pmax 两个边界用二元变量标记 Cons = [Cons, x(i) <= Pmax(i)*u(i)]; % u=0 时 x=0 Cons = [Cons, x(i) >= Pmax(i) - Pmax(i)*v(i)]; % v=0 时 x=Pmax Cons = [Cons, lambda(i) >= 0, lambda(i) <= M*u(i)]; Cons = [Cons, mu(i) >= 0, mu(i) <= M*v(i)]; Cons = [Cons, u(i) + v(i) >= 1]; % 至少一个边界未激活 end Cons = [Cons, p >= pmin, p <= pmax]; %% 目标函数:上层收益最大化 Obj = -(p - c0) * sum(x); % YALMIP 默认最小化,加负号 %% 求解:双线性目标需要打开 NonConvex ops = sdpsettings('solver', 'gurobi', 'gurobi.NonConvex', 2, 'verbose', 1); optimize(Cons, Obj, ops); %% 输出 p_star = value(p); x_star = value(x); fprintf('均衡价格: %.2f 元/kWh\n', p_star); fprintf('总充电量: %.2f kW\n', sum(x_star)); fprintf('充电站收益: %.2f 元\n', (p_star - c0) * sum(x_star));逻辑说明:第 18 行的平稳条件把上层价格 p 和下层充电量 x 耦合在一起,等价于“每台车在自己的约束内对 p 做出了最优响应”。第 20 到 23 行的二元变量组合是互补条件的线性化替代,u 和 v 分别标记 x 是否贴在两个边界上。正常情况下三选一:内点时 u=1、v=1;x=0 时 u=0、v=1;x=Pmax 时 u=1、v=0。
参数说明:M 取 30 不是随便拍的。a 的上限在 14 左右,pmax 是 15,价格和支付意愿的最大量级刚好 30 出头。大 M 太大会破坏数值精度,太小又会把可行域切掉,这个平衡在第 5.4 节具体展开。目标函数里(p - c0) * sum(x)是双线性项,Gurobi 默认不处理非凸二次目标,必须显式打开NonConvex=2,否则求解器直接报错。
3.3 结果怎么读:价格响应、收益组成与补写曲线
跑完先看 p_star 落在什么位置。一个合理的均衡价格应该高于购电成本 c0,同时低于车主的支付意愿上限。如果 p_star 贴着 pmax 走,说明上层吃定了车主的响应弹性,价格机制没有起到“引导”作用,这往往是 a_i 设得普遍偏高;如果 p_star 贴近 pmin,说明车主对价格极其敏感,稍微涨价大家就不充了。
接着画一张充电量分布图,能直观看到哪台车贴了边界、哪台车在内点:
figure; bar(x_star, 0.5); hold on; yline(p_star, '--r', '均衡价格'); xlabel('EV 编号'); ylabel('充电量 (kW)'); legend('充电量', '均衡价格', 'Location', 'best');一块块柱子整齐贴着 Pmax 说明用电意愿很强烈;如果柱子和 Pmax 之间留出明显空隙,说明价格已经压到了部分车主的边际效用以下,博弈正在起作用。收益组成也要拆开看:总收益 = 电量收益 - 购电成本,把这两项分别打印出来,能判断是量的贡献大还是价的贡献大,这对第 4 章的敏感性分析很有用。
4. 三个必调参数与敏感性分析:价格上限、车辆规模、配网容量
4.1 价格上限 pmax:领导者收益的拐点在哪里
很多初学者以为 pmax 只是用来避坑的“安全阀”,实际它是一个重要的政策参数。一些地区的充电服务费有上限,pmax 就是那个上限值;上限越低,充电站的定价空间越窄,收益曲线形态完全不同。
把 3.2 的求解过程封装成函数solve_stackelberg(n, pmax, Cap),返回均衡价格、总充电量和收益,然后做一次扫描:
pmax_list = 8:0.5:20; rev_list = zeros(size(pmax_list)); for k = 1:length(pmax_list) [~, ~, rev] = solve_stackelberg(5, pmax_list(k), inf); rev_list(k) = rev; end plot(pmax_list, rev_list, '-o'); xlabel('价格上限 pmax (元/kWh)'); ylabel('充电站收益 (元)');观察两点:收益峰值出现在哪个 pmax;峰值之后收益是走平还是回落。走平说明上限已经不再约束均衡价格,回落则说明过高的价格自由度反而让上层“贪心”地把价格抬过头,吓退了一批车主,总电量损失超过了单价提升的收益。这就是 Stackelberg 博弈里“权力不总是越大越好”的直观体现。
实际代码包里这类循环通常会配一个for加save,把每次运行的均衡结果存成.mat,方便后续画曲线上限对比图。注意每次求解都要重新调用 Gurobi,5 台车小算例很快,但如果后面规模放大到 200 台,这个循环会非常慢,最好改成并行parfor。
4.2 EV 规模从 5 辆扩到 200 辆:求解时长与均衡形态
把 n 从 5 拉到 200,变量数和二元变量数同时增长。MIQP 的复杂度主要来自二元变量,这台模型里每台车两个二元变量,200 台就是 400 个 0-1 变量,配合双线性目标,求解时长会明显上升但通常还能接受。
n_list = [5, 10, 20, 50, 100, 200]; time_list = zeros(size(n_list)); for k = 1:length(n_list) t0 = tic; [p_star, x_star, ~] = solve_stackelberg(n_list(k), 15, inf); time_list(k) = toc(t0); end跑之前先想清楚一个问题:当 n 变大,a_i 和 b_i 还是按12 + 2*rand生成吗?如果随机分布不变,规模上来以后“总充电量”一定增长,但均衡价格可能变化不大。更有信息量的做法是固定 a、b 的分布,只增加车辆数量,看单位车辆的平均收益是否被摊薄——这是充电站扩建时最关心的指标。
如果跑 200 台时求解器明显变慢,先把verbose关掉,再确认是不是NonConvex=2导致的 MIP 节点数爆炸。有时候把二元变量去掉,改成第 6 章的枚举法,反而更快。200 台车的场景不是非要 MIQP 不可。
4.3 配网容量约束收紧:削峰后的收益损失怎么评估
真实充电站不可能无限供电,变压器容量、线路载流量都会限制总充电功率。设配电网给的容量上限为 Cap,给上层加一条物理约束:
sum(x) <= Cap注意这条约束加在上层还是下层,结果完全不同。它本质是配电网对充电站的物理限制,不是车主自己能感知的效用约束,所以常见做法是加到上层;如果加到下层,意味着车主之间要抢容量,下层就必须引入耦合约束,KKT 条件要从头再写一遍,很多代码包刻意回避了这一层复杂性。
加约束后重新求解:Cap 足够大时结果和原来一致;Cap 开始卡住总充电量后,均衡价格会怎么动?这是这个模型最有意思的现象。当容量受限,充电站没有动力再用低价吸引更多电量,反而可能抬高价格,优先服务支付意愿最高的那批车主。收益曲线会在某个 Cap 出现拐点,这个拐点就是在“增容成本”和“收益损失”之间做投资决策的参考点。
需要警惕的是,如果 Cap 过紧,模型可能显示不可行。这时要检查 x_i 的下界——如果你额外加了强制充电量,而下层的 KKT 条件不允许亏本充电,不可行是正常的,说明这个场景下价格机制本身就撑不起需求。
5. 复现 zip 代码的 5 条高频踩坑记录与排查清单
5.1 现象:求解器报错 Nonconvex,双线性目标引发的第一次翻车
现象是 Gurobi 启动后直接弹错,提示二次目标非凸,甚至 YALMIP 连求解器都选不出来。原因就是上层目标(p - c0) * sum(x)里 p 和 x 都是变量,问题本质是非凸 MIQP,Gurobi 默认按凸问题检查。解决方法是显式打开:
sdpsettings('solver', 'gurobi', 'gurobi.NonConvex', 2)如果换成 Cplex 或 Mosek,注意各自的开关写法不同,Mosek 对非凸二次目标支持很差,建议直接用 BMIBNB 或 BONMIN 这一类全局求解器兜底。我自己遇到最多的情况是,明明代码包里提供了这个设置,但用户单独复制了某一段求解代码而漏掉了设置行,所以排查第一步永远是检查sdpsettings是否完整。
5.2 现象:模型显示无界,价格上限丢失的典型没收尾
现象是求解器返回Infeasible or unbounded,或者 p_star 打印出来是个离谱的数值。九成原因是 pmax 没设或设得太大,上层可以无限抬高价格;另一半原因是目标函数符号写反了。YALMIP 的optimize默认最小化,代码里必须对收益取负号。我见过一个案例,代码里写的是Obj = (p - c0) * sum(x),结果求解器把价格压到 pmin、充电量拉到最大,收益变成了负数,人还一脸懵。
排查步骤固定:先检查约束列表里有没有p >= pmin, p <= pmax;再检查目标函数负号;最后把 Obj 打印出来对着量级看一眼。这三步能滤掉 80% 的无界问题。
5.3 现象:数值和论文对不上,量纲没统一结果差一个量级
你复现了一段代码,算出的 p_star 是 9.2,论文里写 0.92,十倍关系。这不是代码写错了,是量纲不一致。论文里常用元/MWh,而代码里随手写的是元/kWh;购电成本 5 元和 0.5 元,差一个数量级,反应函数计算结果当然对不上。
解决只有一条路:全模型统一量纲。我习惯先在参数表里强制加一列“单位”,每个参数后面标注,再写一行断言检查量纲一致性。不要相信记忆,一个月后重跑这个 zip,量纲就是第一个要查的地方。这类代码包最常见的数据文件里,a_i 和 p 的量纲经常混着来,尤其在把多个论文的代码拼在一起时。
5.4 现象:大 M 一开大就瞎:互补条件线性化的玄学区间
大 M 看起来是个简单参数,实际是 KKT 线性化最容易踩坑的地方。M 设成 1e6,求解结果经常出现“x=0 但对偶变量 λ 也等于 0”这种违背互补条件的伪解;M 设成 10,又可能把真正可行的边界解切掉。
原因是大 M 同时进入平稳约束和互补约束,它必须大于所有可能出现的最优对偶变量,但又不能大到破坏求解器数值精度。经验做法是估算对偶变量的上界:对这台模型,λ 的量级不会超过 max(a_i) + pmax,μ 的量级同理,所以 M 取 30 左右足够,不建议超过 100。调试时可以把 M 从大到小扫一遍,看均衡价格和充电量是否稳定在一个区间。如果 M 的影响显著,说明模型本身有数值病态,先修量纲,不要和大 M 较劲。
5.5 现象:全部充电量为 0,支付意愿参数 a 和价格区间配合失误
跑出来的 x_star 全 0,充电站收益为 0 或为负。背后的数学原因是:当电价 p 高于车主的支付意愿 a_i 时,反应函数(a_i - p)/b_i为负,截断到 0,所有车都不充了。也就是说,pmax 设得比所有 a_i 都高,而上层只要稍微把价格抬起来一点,就会把所有车主赶走,最后均衡只能落在“没人充电”这个糟糕的点上。
解决方式不是去约束上层的行为,而是调整参数分布区间。让a_i不低于pmax的 60% 左右,保证在 pmax 附近仍有部分车主愿意充电;同时把 c0 保持在 pmin 之上,否则上层没有经济动机去压低价格。排查时打印一行min(a), max(a), pmin, pmax,四个数放在一起,几乎立刻就能判断参数区间是否自洽。
6. 收尾验证技巧:用枚举和集中式最优给 Stackelberg 解做体检
6.1 枚举价格网格,校验 KKT 解的金标准
KKT 路线本身也可能出问题,尤其是在大 M、二元变量、非凸求解器这“三件套”组合下。最可靠的验证方法,是利用下层没有耦合约束时反应函数可以解析求解这一点,直接枚举价格:
p_grid = linspace(pmin, pmax, 2000); rev_grid = zeros(size(p_grid)); for k = 1:length(p_grid) pk = p_grid(k); xk = min(max((a - pk) ./ b, 0), Pmax); rev_grid(k) = (pk - c0) * sum(xk); end [rev_enum, idx] = max(rev_grid); p_enum = p_grid(idx);把枚举得到的 p_enum 和 rev_enum 与 3.2 节 KKT 求解结果对比。偏差在 1% 以内基本可以放心;偏差大就要回头查大 M 取值和互补条件。这个方法还有一个额外的好处:枚举法对求解器版本零依赖,任何环境都能跑,是复现 zip 代码时最扎实的“后悔药”。
6.2 和集中式最优比效率损失,才知道你建的是不是博弈
最后做一个集中式对照:假设充电站直接把每台车的充电量当决策变量,追求全系统收益最大化,不通过价格中介:
xc = sdpvar(1, n); Obj_c = -sum(a .* xc - 0.5 * b .* xc.^2 - c0 * xc); optimize([xc >= 0, xc <= Pmax], Obj_c, sdpsettings('solver', 'gurobi'));集中式收益一定不差于 Stackelberg 均衡收益,两者的差距就是价格机制带来的效率损失。差距小,说明充电站的定价引导很高效;差距大,说明车主效用曲线的峰值区域和购电成本之间错位明显,可能需要考虑分时段定价而不是单一电价。我自己做这类项目时的习惯是先跑枚举、再跑集中式对照,最后才信 KKT 路线的结果——三个数能对上,才算真的把博弈模型复现通了。希望这套思路对你手里的这个 zip 也有帮助。
本文还有配套的精品资源,点击获取