简介一套面向智能电网电动汽车充电管理场景的MATLAB实现方案基于YALMIPCPLEX构建主从博弈模型适用于电力系统优化、充电负荷调控方向的科研与工程人员。方案模拟智能小区代理商作为领导者制定电价电动汽车用户作为跟随者调整充电行为通过Stackelberg均衡求解实现供需平衡与成本优化。资源以zip压缩包提供大小约406KB压缩包内当前文件数量显示为0具体文件类型未标注。已有4410人学习下载在同类博弈优化资源中具有一定参考热度。内容主要涉及主从博弈的数学模型定义、YALMIP优化问题构建、CPLEX求解器调用以及不同场景下的数据输入与结果分析脚本便于读者理解定价策略对用户充电决策的影响机制。借助该程序可进一步掌握智能电网环境下电动汽车有序充电的建模与求解思路为扩展多主体协同优化研究提供基础。1. 主从博弈电动汽车管理先让上层亮价再让下层回话一个 30 个车位的公共充电站白天电价便宜时桩前冷清晚上下班后所有车挤在同一时段充电。把电价改成峰谷分时后负荷峰削下去一点用户成本也降了但站里利润反而掉了。原因是“管理策略”和“用户响应”是两件事管理方先出牌用户看到电价后再调整充电计划用户的调整反过来又改变管理方的收益。这正是主从博弈Stackelberg game要解决的问题。这篇博文讲的是在 MATLAB 中用 YALMIP 建立主从博弈模型再用 CPLEX 求解转化后的单层 MILP完成电动汽车充电定价与功率分配。适合做充电聚合、有序充电、V2G 仿真的工程师以及想把双层优化从论文落到代码的研究生。2. 先写上下层模型为什么单层优化在这里会失效2.1 单层优化为什么预判不了用户管理方如果只以“全站负荷最平”为目标直接给出每辆车的充电曲线这是一种集中式单层调度前提是用户完全服从指令。实际场景中用户按电价和便利性做自己的最优决策集中调度给出的曲线往往执行不下去。反过来如果完全让用户自由选时段所有人会往同一个低电价时段挤价格信号立刻被“挤没”。主从博弈把决策拆成两层上层是充电站或聚合商先公开分时电价下层是 EV 用户在收到的电价下做成本最小的充电决策。上层做决策时已经知道下层会如何响应所谓 Stackelberg 均衡就是“我先出牌并预判你如何应对”的均衡。2.2 上层电价决策和下层充电响应的数学表达上层问题写为max Σ (p_t - c_t) · x_t(p) s.t. p_min ≤ p_t ≤ p_max x_t(p) D_t ≤ C_max其中 c_t 是管理方的购电成本D_t 是背景负荷C_max 是变压器容量。关键在 x_t(p)它不是自由的调度变量而是下层用户问题的最优解。下层问题对给定电价 p 写为min Σ p_t · x_t s.t. 0 ≤ x_t ≤ P_max Σ x_t E这里故意把下层写成线性规划。一方面是 LP 的 KKT 条件好写强对偶严格成立另一方面是转成 MILP 后 CPLEX 可以直接处理。如果换成带电池损耗的二次目标转化思路完全相同但目标里会多出凹二次项对求解器的凸性判断更敏感新手容易被这类细节卡住。参数含义典型值T调度时段数24P_max单桩最大充电功率7 kWE单车日充电需求40 kWhp_min / p_max电价下/上限0.2 / 1.5 元/kWhD_t背景负荷实测或模拟数据C_max接入变压器容量50 kW2.3 双层问题直接交给 CPLEX 会卡在哪CPLEX 的强项是 LP、QP 和 MILP但双层优化即使上下层都是 LP整体可行域通常非凸。把上层目标与下层最优性条件放在一起时会出现 p_t · x_t 这种两个决策变量相乘的项构成非凸二次项。CPLEX 默认不接受非凸二次目标直接求解会报错或给出不可信的结果。标准做法分两步第一步对下层写 KKT 条件把“x 是 p 的函数”改写成一组约束第二步用强对偶把上层目标中的双线性项 px 替换成对偶变量的线性表达式得到一个单层 MILP。在论文里这叫单层化在实际工程里就是“让模型能被 MILP 求解器接受”的常规动作。建议建模时先自查目标表达式T 24; p sdpvar(T, 1); x sdpvar(T, 1); obj_naive (p - 0.3) * x; % 直接写上层目标 if degree(obj_naive) 1 disp(linear objective); else disp(nonlinear objective, need reformulation); enddegree返回表达式最高次数。上面的 obj_naive 是 p 与 x 的乘积次数为 2说明上层目标中存在双线性项不能直接交给 CPLEX。这一步检查在 YALMIP 里很便宜却能提前暴露“模型能不能被求解器接受”的根本问题。提示下层模型必须是凸的KKT 才是充要条件。本文的 LP 下层满足要求如果你把下层改成非凸整数规划就不能直接套下面的 KKT 单层化。3. 用KKT和强对偶把双层拍平YALMIP建模的核心步骤3.1 对下层写KKT四组条件缺一不可针对 2.2 节的下层问题写出拉格朗日函数L Σ [p_t·x_t λ_ub_t·(x_t - P_max) - λ_lb_t·x_t] ν·(Σx_t - E)其中 λ_ub_t ≥ 0 对应 x_t ≤ P_maxλ_lb_t ≥ 0 对应 -x_t ≤ 0ν 自由对应等式约束 Σx_t E。KKT 条件分四组平稳性∂L/∂x_t p_t λ_ub_t - λ_lb_t ν 0原始可行性0 ≤ x_t ≤ P_maxΣx_t E对偶可行性λ_ub_t ≥ 0λ_lb_t ≥ 0互补松弛λ_ub_t·(P_max - x_t) 0λ_lb_t·x_t 0平稳性、原始可行性和对偶可行性都是线性约束只有互补松弛是乘积为 0 的非线性条件。YALMIP 提供了kkt函数可以自动生成下层 KKT但互补条件仍需手工线性化自己手写反而更容易控制乘子符号。这里推荐手写因为符号方向错了强对偶表达式会差一个负号最终解出的“均衡”完全错误。3.2 互补松弛的大M线性化二进制变量到底怎么加互补条件 λ·(约束剩余量) 0 等价于“要么乘子为 0要么约束取等号”。对每个时段引入两个二进制变量用大 M 约束实现二选一M 500; b_ub binvar(T,1); % 控制 x P_max 是否激活 b_lb binvar(T,1); % 控制 x 0 是否激活 comp_cons [Pmax - x M * b_ub, lambda_ub M * (1 - b_ub), ... x M * b_lb, lambda_lb M * (1 - b_lb)];逐条看当 b_ub 0 时Pmax - x ≤ 0配合原始可行性 x ≤ Pmax推出 x Pmax此时 λ_ub 可以大于 0。当 b_ub 1 时λ_ub ≤ 0配合对偶可行性 λ_ub ≥ 0推出 λ_ub 0此时 x 可以小于 Pmax。b_lb 同理负责处理 λ_lb_t·x_t 0。大 M 的取值直接影响 CPLEX 分支效率。M 太小会把最优解排除在可行域外M 太大则让 LP 松弛过松分支定界要探索大量节点。对本文参数500 是一个安全的起点如果对偶变量可能到几千就相应放到 5000但不要无脑设成 1e6。3.3 用强对偶消掉上层目标里的px上层目标里 Σ(p_t - c_t)·x_t 含 p_t·x_t。消除它的办法是利用下层 LP 的强对偶。对 2.2 节的下层问题对偶问题的最优目标等于原问题最优目标px -P_max·Σλ_ub_t - E·ν注意 ν 前面的符号。这里的 ν 对应拉格朗日里的 ν·(Σx_t - E)不是 ν·(E - Σx_t)。如果把等式约束写成 E - Σx_t 0强对偶表达式就变成 px -P_max·Σλ_ub_t E·ν符号完全不同。这是手写 KKT 时最容易出错的点。替换后的上层目标只包含线性项revenue -Pmax * sum(lambda_ub) - E * nu; obj -(revenue - c * x);YALMIP 默认做最小化所以对上层最大化目标取负号。c * x 是购电成本revenue 是用户充电费用两者之差是聚合商利润。经过替换目标里不再有 p 和 x 的乘积模型从非凸双线性问题变成 MILP。3.4 合并约束集CPLEX能解的单层MILP把上层约束、KKT 约束和互补线性化约束合并成一个约束集cons [p pmin, p pmax, ... x D Cmax, ... % 上层容量约束 x 0, x Pmax, sum(x) E, ... % 原始可行性 p lambda_ub - lambda_lb nu 0, ... % 平稳性 lambda_ub 0, lambda_lb 0, ... % 对偶可行性 comp_cons]; % 线性化互补模型变量包括 p、x、λ_ub、λ_lb、ν以及二进制 b_ub、b_lb。所有约束都是线性等式或不等式CPLEX 用分支切割法求解。求解结果同时给出上层电价和下层充电功率两侧都不需要迭代求解一次求解就是 Stackelberg 均衡。注意 x D Cmax 是逐时段约束x 和 D 都必须是同维度向量。如果 D 在 MATLAB 里是行向量代码会触发隐式扩展把约束变成矩阵约束YALMIP 不会立刻报错但模型规模会悄悄变大求解结果也不对。这种隐蔽的维度 bug 在双层模型里比求解器报错更难排查。4. 可跑的MATLABYALMIPCPLEX代码和关键参数4.1 先确认CPLEX被YALMIP正常调用模型写得再对求解器没接上也是白搭。首次运行前先做一次空模型测试ops sdpsettings(solver, cplex, verbose, 1); sol optimize([], [], ops); if sol.problem 0 disp(CPLEX ready); else disp(sol.info); endsol.problem 0 表示求解成功。如果在optimize阶段提示没有找到求解器先检查yalmiptest输出或执行which cplex确认 MATLAB 能访问 CPLEX 的可执行文件。YALMIP 对 CPLEX 的调用是通过 MATLAB 接口完成的常见问题不是 CPLEX 本身没装好而是 MATLAB 路径里同时存在多个求解器YALMIP 默认选了别的求解器。4.2 单桩24时段完整代码下面是一段可以直接复制的完整代码。背景负荷用正弦函数模拟只是为了演示实际使用应替换成实测数据。%% 参数设置 T 24; Pmax 7; % 单桩最大功率kW E 40; % 单日需求kWh pmin 0.2; pmax 1.5; % 电价上下限 Cmax 50; % 变压器容量kW D 20 10 * sin((0:T-1) / T * 2 * pi); D D(:); % 强制转成列向量避免维度隐患 c 0.3 0.1 * rand(T, 1); % 购电成本元/kWh %% YALMIP变量 p sdpvar(T, 1); % 电价上层决策 x sdpvar(T, 1); % 充电功率下层决策 lambda_ub sdpvar(T, 1); % 对偶乘子x Pmax lambda_lb sdpvar(T, 1); % 对偶乘子x 0 nu sdpvar(1, 1); % 对偶乘子sum(x) E b_ub binvar(T, 1); % 互补松弛开关 b_lb binvar(T, 1); M 500; % 大M值 %% 约束 cons [p pmin, p pmax, ... x D Cmax, ... x 0, x Pmax, sum(x) E, ... p lambda_ub - lambda_lb nu 0, ... lambda_ub 0, lambda_lb 0, ... Pmax - x M * b_ub, lambda_ub M * (1 - b_ub), ... x M * b_lb, lambda_lb M * (1 - b_lb)]; %% 上层目标强对偶替换 px revenue -Pmax * sum(lambda_ub) - E * nu; obj -(revenue - c * x); %% 求解 ops sdpsettings(solver, cplex, verbose, 1); ops.cplex.mip.tolerances.mipgap 1e-4; ops.cplex.mip.limits.timelimit 60; sol optimize(cons, obj, ops); %% 结果 if sol.problem 0 figure; subplot(2,1,1); stairs(value(p)); title(分时电价); subplot(2,1,2); bar(value(x)); hold on; plot(D, r); legend(EV充电功率, 背景负荷); else disp(sol.info); end代码里每个变量都对应 3.2 和 3.3 节中的数学符号。revenue是用强对偶算出的用户充电费c * x是购电成本obj是聚合商利润的负值。求解后value(p)是均衡电价value(x)是用户充电功率value(lambda_ub)和value(lambda_lb)是对偶乘子可用于后续校验。注意D D(:)这行不能省。MATLAB 2016b 之后支持隐式扩展行向量 D 和列向量 x 相加会生成矩阵YALMIP 会把矩阵约束当成逐元素约束模型维数膨胀结果无法解释。4.3 三个直接影响解质量的CPLEX参数CPLEX 参数很多对这类单层化 MILP 模型影响最大的三个是 mipgap、时间限制和线程数。参数设置路径建议值作用MIP 相对间隙mip.tolerances.mipgap1e-4控制最优性误差太小会显著增加分支节点求解时间上限mip.limits.timelimit60避免多车场景陷入长时间搜索并行线程数threads4加速求解共享计算环境设置过高反而变慢代码中的设置方式ops sdpsettings(solver, cplex, verbose, 1); ops.cplex.mip.tolerances.mipgap 1e-4; ops.cplex.mip.limits.timelimit 60; ops.cplex.threads 4;mipgap 是相对误差不是绝对误差。如果目标函数本身数值很大1e-4 可能过于严格建议先观察一个简单场景的求解时间再决定是否放宽到 5e-4。timelimit到点后 CPLEX 会返回当前最好整数解但 sol.problem 仍然可能是 0需要额外检查sol.solvertime和sol.info确认解是收敛还是超时截断。4.4 从单桩到N辆EV的矩阵化改造实际项目里车不只一辆把 x、lambda_ub、lambda_lb 从向量改成矩阵即可。N 辆车每辆车有自己的总需求 E_n 和最大功率 P_max,n时段索引仍为 T。N 5; E_i [40; 50; 30; 45; 60]; % 每辆车需求 Pmax_i [7; 7; 22; 7; 11]; % 每辆车最大功率 x sdpvar(N, T, full); lambda_ub sdpvar(N, T, full); lambda_lb sdpvar(N, T, full); b_ub binvar(N, T, full); b_lb binvar(N, T, full); nu sdpvar(N, 1); % 每辆车一个等式乘子约束按维度对齐。聚合功率约束写成sum(x,1) D Cmax平稳性约束写成cons [repmat(p, N, 1) lambda_ub - lambda_lb repmat(nu, 1, T) 0, ... sum(x, 2) E_i];目标函数把每辆车的 revenue 相加购电成本按时段聚合revenue sum(-Pmax_i .* sum(lambda_ub, 2) - E_i .* nu); obj -(revenue - c * sum(x, 1));这里sum(x,1)是 1×T 的聚合功率向量转置后与 c 做内积得到总购电成本。多车模型的变量数量和二进制数量线性增长CPLEX 的求解时间通常不会是线性增长的mipgap 和 timelimit 的意义在这里才真正体现。5. 均衡质量校验互补残差和固定电价基准5.1 互补残差确认KKT没有被大M带偏求解完成后第一件事不是画图而是检查互补条件是否真的满足。KKT 互补条件是 λ_ub_t·(Pmax - x_t) 0 和 λ_lb_t·x_t 0。用 CPLEX 求出的数值解互补残差一般应该在 1e-4 量级。res_ub value(lambda_ub) .* (Pmax - value(x)); res_lb value(lambda_lb) .* value(x); max_res max([res_ub(:); res_lb(:)]); if max_res 1e-3 warning(互补残差过大检查M取值或CPLEX容差); end残差过大的原因一般有两个一是大 M 取值偏小把真正的最优互补关系排除在可行域外二是 CPLEX 的整数可行性容差太松二进制变量没有严格取到 0 或 1。处理方法分别是增大 M 和把mip.tolerances.integrality从默认值调到 1e-6 或更小。M 不能无限大否则 LP 松弛质量变差分支定界会变得很慢。5.2 固定电价基准区分“博弈解”和“拍脑袋解”一个很有效的对照实验把电价固定成某个常数解下层 LP得到用户响应和聚合商利润。这个结果可以当作“没有博弈”的基准。p_fixed 0.6 * ones(T, 1); x_lp sdpvar(T, 1); cons_lp [x_lp 0, x_lp Pmax, sum(x_lp) E]; optimize(cons_lp, p_fixed * x_lp, sdpsettings(solver, cplex, verbose, 0)); profit_fixed (p_fixed - c) * value(x_lp);把 profit_fixed 与主从博弈的聚合商利润做对比。如果主从解利润低于固定电价基准说明 KKT 或强对偶表达式符号写反了或者互补条件线性化写错了。这个基准不需要额外调参几行代码就能完成是排查“模型写着写着就错了”的最快路径。5.3 更贴近实际的三个改动方向第一把下层 LP 换成含电池损耗的 QP 下层目标为 Σ(p_t·x_t α·x_t²)。强对偶会引入凹二次项CPLEX 需要把optimalitytarget设置为 3 才能处理求解速度会下降。第二加入每辆车的到达和离开时段窗口把窗口外的 x_{n,t} 固定为 0甚至用二进制变量表示“是否在充电”代价是模型从 MILP 变成 MIQP。第三把单车需求 E 改成概率分布用多个场景做蒙特卡洛仿真观察电价策略在不同需求样本下的鲁棒性。差分进化、粒子群这类元启发式算法可以处理双层模型但在小规模确定性模型上尽快换成 CPLEX 求解 MILP 才是更可控的做法。本文还有配套的精品资源点击获取
