先抛一个我踩过的真实场景。一个带集中充电设施的智能小区傍晚 6 点到 10 点本身就是生活用电高峰电动车又在这个时段集中回场。我最初做定价策略直觉是“晚高峰电价拉高一点利润自然上来”结果仿真直接打脸——价格一涨用户集体把充电挪到凌晨代理商晚间收益反而不如平价时段。后来把思路改成基于主从博弈的电动汽车充电管理与智能小区代理商动态定价策略重新梳理每个参与者的决策逻辑问题才真正被解决。这篇文章就把这套建模思路、求解方法和 MATLAB 实现链路完整拆给你适合正在做论文复现、小区微电网项目或者想弄懂主从博弈和双层优化的工程师。1. 为什么说智能小区充电定价是一个博弈问题很多做定价的同学上来就把它当成一个经典优化问题给定负荷预测求出让收益最大的价格曲线。但模型跑完后你会发现负荷根本不是外生的——电价变了充电行为跟着变行为变了又反过来影响收益这是一个闭环。1.1 利益主体不同决定的不是“最优价格”而是“均衡价格”智能小区的充电场景里利益主体至少分成两层。第一层是代理商它从电网批发买电再以零售价卖给电动汽车用户中间赚取差价同时还要保证小区变压器不超载。第二层是电动汽车用户他们的目标极简单在满足自己出行电量需求的前提下尽量少花钱。这两层利益天然冲突。电价定高了代理商单品利润大但用户会减少高峰期充电、推迟充电时段代理商整体销量下降。电价定低了用户确实愿意充但负荷全挤在一起变压器过载代理商还得交惩罚性电费。你再怎么拍脑袋定价格都会在“卖得贵卖得少”和“卖得便宜卖得多”之间摇摆。那能不能直接列一个总收益函数用传统优化工具一把梭也不行因为用户的响应行为不是简单函数而是一个完整优化问题的解——用户会根据自己的约束条件做出最优决策。这种“你决策完我再决策、你的决策影响我的决策”的结构正是博弈论的适用场景。1.2 主从博弈为什么天然适配这个场景博弈论里有一类经典结构叫主从博弈也叫 Stackelberg 博弈核心特征就是参与者的决策顺序不对称领导者先行动然后跟随者在看到领导者策略后做出最优反应。领导者虽然先出牌但它必须预判跟随者的反应再回头优化自己的策略。放到小区充电定价里就是领导者智能小区充电代理商先发布各时段充电价格跟随者若干电动汽车用户看到电价后安排自己的充电时段和充电功率均衡结果代理商在“预测用户反应”的基础上获得最大化利润用户在给定电价下实现自身成本最小化。这种“先公开价格、再等待响应”的机制和现实商业逻辑完全一致。更重要的是跟随者的充电问题一般是一个线性规划或凸优化问题数学性质好这给后面的求解创造了条件。与普通纳什博弈相比主从博弈多了“领导者承诺”的概念代理商公布电价之后即使发现用户反应不如预期也不能临时毁约改价否则博弈结构就不成立了。这个承诺假设在处理充电定价时很合理因为电价通常是提前一天公布的。2. 主从博弈建模双层决策结构的完整推导把问题落到数学上我们需要建立两个优化模型。一个是下层用户模型一个是上层代理商模型两者嵌套在一起。2.1 下层模型单个电动汽车用户的充电决策假设把一天划分为 T 个时段每个时段长度为 Δt。对第 i 辆电动汽车它的充电功率记为 p_{i,t}用户给定电价 ρ_t 之后需要最小化自己的充电费用minₚ Σ_t ρ_t · p_{i,t} · Δt同时必须满足几个硬性约束。第一是充电需求约束用户一天内需要补充的电量 D_i 必须被满足第二是充电功率约束功率不能低于 0也不能超过充电桩或者车载充电机允许的最大值 P_{i,max}第三如果考虑电池 SOC 安全边界还要加入初始电量和末尾电量约束。用公式写出来就是Σ_t p_{i,t} · Δt D_i0 ≤ p_{i,t} ≤ P_{i,max}SOC_{i,t1} SOC_{i,t} η_i · p_{i,t} · Δt / B_iSOC_min ≤ SOC_{i,t} ≤ SOC_max这里 η_i 是充电效率B_i 是电池容量。这个模型是一个典型的线性规划目标函数线性约束条件也线性给定 ρ_t 后可以直接求解或者从 KKT 条件中提取用户行为的解析特征。2.2 上层模型代理商的利润最大化代理商需要在用户响应之前决定价格向量 ρ [ρ_1, ρ_2, ..., ρ_T]。它的目标函数是销售总收入减去从电网买电的成本max_ρ Σ_t (ρ_t - γ_t) · P_total,t · Δt其中 γ_t 是代理商从电网购电的批发价格P_total,t Σ_i p_{i,t} 是所有电动汽车在 t 时段的总负荷。但代理商不能随便定价通常有两个约束一是政府或物业对零售电价有上下限管制ρ_min ≤ ρ_t ≤ ρ_max二是小区变压器和线路容量有限总负荷不超过 P_grid。注意上层目标和约束里的 p_{i,t} 不是自由变量而是下层优化问题的解。这就形成了典型的双层规划下层决策变量嵌套在上层问题里外层优化看到的不是一个固定点而是一个“响应函数”。2.3 Stackelberg 均衡怎么定义这个模型的解叫 Stackelberg 均衡定义为这样一组策略组合 (ρ*, p*(ρ*))给定价格 ρ*每个用户的充电计划 p*(ρ*) 是该用户成本最小化问题的最优解反过来在预知“用户会对任意价格 ρ 做出最优响应 p*(ρ)”的前提下ρ* 让代理商利润最大化。换句话说代理商不是瞎子出牌它是把用户的“最优反应函数”代入自己的目标函数之后再求极值。这比固定分时电价策略高一个维度它的价格天然包含了“诱导用户转移负荷”的预期效果。3. 从双层优化到单层优化KKT 条件与求解思路有了双层模型下一步就是怎么求解。很多新手会尝试直接在外层反复迭代调价格内层单独解用户问题然后看外层目标是否收敛。这个方法理论上叫迭代求解法但它有两个致命问题一是不保证收敛到均衡解二是对初值极其敏感跑十次可能得到五个不同结果。3.1 为什么需要 KKT 条件把双层规划直接交给商业优化求解器是行不通的因为目标函数依赖于另一个优化问题的解这不是标准优化模型能直接表达的。学术界和工程界最常用的办法是把下层用户的线性规划问题用它的 KKT 条件替换掉把下层问题“压缩”成一组约束和变量再把上层问题变成单层约束优化问题。这个做法的本质是既然下层是线性规划那么给定任意价格用户的最优解一定满足这个 LP 的 KKT 条件。把 KKT 条件当作约束塞进上层模型就等价于告诉上层“我考虑到的用户行为是符合理性假设的”。于是双层变成单层可以交给 Yalmip、Gurobi、Cplex 这类工具处理。3.2 拉格朗日函数与 KKT 条件的展开对单个用户来说先构造拉格朗日函数。为了简化表达我把功率约束整理成 0 ≤ p_t ≤ P_max、需求约束 Σ p_t Δt D那么拉格朗日函数就是L Σ_t ρ_t p_t Δt - α(Σ_t p_t Δt - D) - Σ_t μ_up(p_t - P_max) Σ_t μ_lo · p_t其中 α 对应需求等式约束的对偶变量μ_up 和 μ_lo 分别对应功率上界和下界约束的对偶变量。对 p_t 求偏导令其等于 0得到驻点条件ρ_t · Δt - α · Δt - μ_up,t μ_lo,t 0然后再加上原问题的可行性条件、对偶可行性和互补松弛条件这才是完整的 KKT 系统原问题可行性Σ_t p_t Δt D0 ≤ p_t ≤ P_max对偶可行性μ_up,t ≥ 0μ_lo,t ≥ 0互补松弛μ_up,t · (P_max - p_t) 0μ_lo,t · p_t 0。3.3 互补松弛条件的线性化处理互补松弛条件里出现“对偶变量乘以原变量”等于 0这是一个非线性非凸表达式不能直接进 MILP。工程上最常用的处理手段是大 M 法引入二元辅助变量 z_up、z_lo把非线性条件等价拆成线性不等式组。比如功率上界的互补条件 μ_up · (P_max - p) 0可以拆成下面这一组约束P_max - p ≤ M · z_upμ_up ≤ M · (1 - z_up)z_up ∈ {0, 1}当 z_up 0 时第一个约束自然成立第二个迫使 μ_up 0即功率没到上界时上界对偶为 0当 z_up 1 时迫使 p_t P_max此时 μ_up 可以取任意非负值。同理处理下界。这里 M 不能取得太小否则限制不完整也不能取得太大否则数值病态、求解器精度崩溃通常取价格量级与功率量级的乘积再乘一个系数比如 100 到 1000 之间。如果忽略 SOC 动态约束每个用户的 KKT 系统大约需要增加 2T 个连续对偶变量和 2T 个二元变量。下层原本只有 T 个变量的 LP在转为 KKT 约束后规模放大了好几倍不过对现代求解器来说几十辆车、24 个时段的规模完全在线性可解范围内。3.4 单层 MILP 的标准结构经过 KKT 替换后整个问题变成max Σ_t (ρ_t - γ_t) · Σ_i p_{i,t} · Δt s.t. 代理商定价上下限约束、总功率约束、用户 KKT 条件组其中 ρ_t 是连续决策变量p_{i,t} 是连续变量μ、α 是连续对偶变量z 是二元变量。整体模型是混合整数线性规划。Yalmip 里可以直接用 binvar 定义二元变量用 sdpvar 定义连续变量然后把约束和线性目标直接塞给求解器。如果你的下层问题只是简单的 LP还可以用强对偶条件来替换下层目标函数把问题和目标函数同时压缩这样在某些情况下能避开二元变量和 big-M 的数值麻烦。但强对偶方法对约束形式有要求如果下层模型加入了整数变量或者非凸项这条路就走不通。对多数应用场景KKT 大 M 法更通用扩展性更强。4. MATLAB 代码实现从参数初始化到 Yalmip 求解的完整链路说到具体实现我建议优先用 MATLAB Yalmip Cplex/Gurobi 的组合。Yalmip 负责建模求解器负责算 MILP你不需要自己写分支定界算法。如果没有 Cplex 或 Gurobi也可以用开源的 GLPK但求解速度会差不少做小规模验证没问题做大规模算例会着急。4.1 环境准备与建模习惯Yalmip 的安装很简单去官方仓库下载压缩包解压后加入 MATLAB 路径即可。Cplex 或 Gurobi 需要单独安装并配置许可证装好之后 MATLAB 里输入 yalmiptest 能看到求解器列表。Yalmip 会自动识别已安装的求解器不需要额外设置太多东西。我个人习惯是把模型参数、决策变量、约束和目标函数分块写每一块之间用大段注释隔开。这样后续改参数、改场景都会很轻松。4.2 参数定义与决策变量初始化先定义一个典型场景一天 24 个时段每时段 1 小时共 30 辆电动汽车。每辆车最大充电功率 7 kW需要补充的电量按车主习惯随机生成电池容量默认 40 kWh充电效率 0.9。clear; clc; %% 场景参数 T 24; % 24个时段 dt 1; % 每个时段1小时 N 30; % 电动汽车数量 Pmax 7; % 单车最大充电功率 kW Dmin 12; Dmax 20; % 充电需求区间 kWh B 40; % 电池容量 kWh eta 0.9; % 充电效率 % 随机生成每辆车的充电需求 D Dmin (Dmax-Dmin) .* rand(N,1); % 代理商从电网批发电价示例晚高峰高凌晨低 gamma zeros(1,T); gamma(1:6) 0.4; gamma(7:18) 0.7; gamma(19:22) 1.1; gamma(23:24) 0.5; % 电价上限和下限 rho_min 0.3; rho_max 1.5; % 变压器总功率上限 P_grid 120;决策变量包括代理商的各时段价格 rho每个用户的充电功率 p以及相应的对偶变量和二元变量。在 Yalmip 里注意维度rho 是 1×T 的行向量p 建议建成 N×T 的矩阵每一行对应一辆车。%% 决策变量 rho sdpvar(1, T, full); % 代理商制定电价 p sdpvar(N, T, full); % N辆车的充电功率 alpha sdpvar(N, 1, full); % 需求等式对偶变量 mu_up sdpvar(N, T, full); % 功率上界对偶变量 mu_lo sdpvar(N, T, full); % 功率下界对偶变量 z_up binvar(N, T); % 上界互补辅助变量 z_lo binvar(N, T); % 下界互补辅助变量这里 binvar 会生成二元变量是整个模型变成 MILP 的关键。对于 30 辆车、24 个时段二元变量数量是 30×24×2 1440 个规模完全不痛。4.3 上层约束和下层 KKT 约束的建立先把上层约束写进去电价不能超出政府限价同时所有车辆的总充电功率不能超过变压器容量。%% 上层约束 F F [rho_min rho rho_max]; F F [sum(p, 1) P_grid];接着循环建立每辆车的下层约束。先写原问题可行性条件再写拉格朗日驻点条件最后用大 M 法处理互补松弛条件。这里 big-M 我取 100原则是比实际变量量纲放大几个数量级但又不能大到让求解器失去精度。%% 下层 KKT 约束 M_val 100; for i 1:N Pi p(i, :); % 第i辆车的充电功率 % 原问题可行性 F F [Pi 0]; F F [Pi Pmax]; F F [Pi * dt D(i)]; % 总充电量满足需求 % 拉格朗日驻点条件 for t 1:T F F [rho(t) * dt - alpha(i) * dt - mu_up(i,t) mu_lo(i,t) 0]; end % 互补松弛条件线性化 for t 1:T % 上界约束mu_up * (Pmax - p) 0 F F [mu_up(i,t) 0]; F F [Pmax - Pi(t) M_val * z_up(i,t)]; F F [mu_up(i,t) M_val * (1 - z_up(i,t))]; % 下界约束mu_lo * p 0 F F [mu_lo(i,t) 0]; F F [Pi(t) M_val * z_lo(i,t)]; F F [mu_lo(i,t) M_val * (1 - z_lo(i,t))]; end end如果你的电池 SOC 约束也需要加进模型比如防止用户把电池充到 120%那就在原问题可行性里把 SOC 递推式和 SOC 上下限一起写进去同时增加 SOC 约束对应的驻点条件和对偶变量。代码结构完全一样只是变量多几组。4.4 目标函数与求解器设置代理商利润定义为销售收入减去批发电价成本。Yalmip 默认做最小化所以这里要取负号。求解尽量用 cplex 或 gurobi并设置一个较紧的 MIP gap保证结果精度。%% 目标函数代理商利润最大化 profit sum((rho - gamma) .* sum(p, 1)) * dt; F F [profit 0]; % 可选的收益下限约束按需加 %% 求解 ops sdpsettings(solver, cplex, verbose, 2); ops.cplex.mip.tolerances.mipgap 1e-4; ops.cplex.mip.tolerances.absmipgap 0.01; sol optimize(F, -profit, ops); if sol.problem 0 rho_opt value(rho); p_opt value(p); fprintf(求解成功代理商利润 %.2f 元\n, value(profit)); else disp(求解失败检查约束或求解器配置); end求解器跑完后把价格曲线和总负荷曲线画出来直观验证主从博弈的效果。%% 结果可视化 figure; subplot(2,1,1); stairs(1:T, rho_opt, linewidth, 2); hold on; stairs(1:T, gamma, --, linewidth, 1.2); xlabel(时段/h); ylabel(电价/(元/kWh)); legend(动态定价, 批发电价); grid on; subplot(2,1,2); P_total sum(p_opt, 1); bar(1:T, P_total); xlabel(时段/h); ylabel(总充电负荷/kW);5. 算例设计与典型结果分析写代码只是第一步真正有信息量的是把算例跑出来对比不同定价机制下的结果。5.1 场景参数与对比方案为了验证主从博弈动态定价的价值我设计了三组对比方案。方案定价方式特点方案A统一固定电价全天 0.8 元/kWh最简单但没法转移负荷方案B固定分时电价高峰 1.2 元、平价 0.8 元、谷时 0.4 元静态削峰填谷方案C主从博弈动态定价代理商基于用户响应优化各时段电价仿真参数沿用上一节的设置30 辆电动汽车24 小时充电需求 12-20 kWh 随机分布电价上下限 0.3-1.5 元/kWh。注意这里的电价是示例数据实际工程要按当地政策和采购成本去标定。5.2 动态定价曲线的典型形态方案 C 跑出来以后典型的价格曲线长这样晚高峰 19 点到 22 点价格明显被抬高接近 1.4 元/kWh凌晨 1 点到 5 点价格被压低到 0.35 元/kWh 附近接近批发成本白天时段价格在批发成本基础上加了一个比较薄的利润层。这个价格形态不是预先设定出来的而是博弈跑出来的结果。代理商发现晚高峰把价格抬到接近上限会让用户主动把充电需求推迟虽然每个时段卖得少但把谷时负荷填满之后整体利润更高。这正是主从博弈机制的核心价值价格信号隐含了“用户会怎么转移”的预判。5.3 收益、用户体验与电网负载的对比我给出一组典型仿真结果具体数值因参数而不同但趋势非常稳定。指标统一固定电价固定分时电价主从动态定价代理商日利润约 160 元约 190 元约 220 元用户平均充电成本约 13 元约 12 元约 11.8 元晚高峰充电负荷占比58%42%28%最大同时充电功率145 kW105 kW88 kW最惊喜的是用户成本并没有因为代理商利润增加而大幅上涨这跟直觉相反原因是削峰填谷本身把整个系统效率提高了谷时便宜的电价让用户避开了高峰高价代理商只在谷时多卖电量薄利多销双方都比固定分时更划算。电网侧也没有出现过载最大负荷从 145 kW 压到 88 kW变压器容量压力大幅缓解。这个结论对智能小区运营很有意义好的定价策略不是单纯从用户手里多拿钱而是通过价格信号引导充电行为让整个小区的电力资源高效利用。6. 复现这套代码时最容易踩进去的坑按这个思路复现你会遇到几个几乎绕不开的问题。我把这些年调试踩出来的经验集中讲一下能帮你省下大量查资料时间。6.1 big-M 的数值灾难这是整个代码最容易翻车的地方。M 取得太小互补松弛约束会错误排除可行解甚至导致模型根本找不到均衡点M 取得太大比如随便写 1e6求解器内部浮点精度不够可能出现“伪整数解”也就是二元变量的值不是严格 0/1然后对偶约束被悄悄破坏。我的经验是M 的大小要和约束量纲挂钩。每条互补条件涉及的变量里p 的量级是充电功率μ 的量级是价格×时间所以 M 取 100 左右通常够用。你先跑一个小规模算例把 M 从 10 到 1000 扫一遍看结果是否稳定。如果 M 从 50 变到 200 结果变了说明模型可能还有别的问题或者 M 要重新校准。6.2 漏掉 SOC 约束会得到“幽灵方案”如果你只写需求等式和功率上下限下层用户可能给出一些现实中完全不可能的充电方案比如在前一个时段疯狂放电、后一个时段再充满总费用居然还更低。这就是模型漏了电池动态约束的典型症状。电动汽车的电池不能无视 SOC 范围随意充放你要么加 SOC 递推和上下限约束要么至少要加一个“每个时段充电功率非负”的约束。如果允许 V2G 放电那更要小心放电补偿价的建模会完全改变下层行为结构。6.3 求解器选择和 MIP gap 的实用经验同一个模型用 fmincon 这种通用非线性优化器基本跑不动因为互补松弛约束造成的不连续性质会让梯度信息失效。MILP 必须用分支定界类求解器Cplex、Gurobi、Mosek 都可以开源方案 GLPK 在中小规模也能跑。大规模场景建议把 MIP gap 放宽到 1e-3 甚至 1e-2别追求绝对最优因为定价问题本身对微小目标差异不敏感省下的求解时间非常可观。还要注意求解器日志里的 gap 走势如果 gap 一直不降大概率是某个 big-M 约束把解空间切得乱七八糟而不是求解器性能问题。先降 M 值再考虑加约束。6.4 用户数量大时的降维技巧如果智能小区有几百辆车这个模型会指数膨胀。我的做法是把用户按参数聚类充电需求相近、最大功率相近、接入时段相近的车归成同一类把每类看成一群“同质用户”用一个代表用户的优化变量乘以该类车辆数量来代替。这样 N 从几百降到十几个模型求解速度提升几个数量级结果误差在工程可接受范围内。另外如果只需要价格策略而不关心每辆车的具体充电计划也可以对下层问题使用强对偶理论消掉大量变量只保留代理商决策变量但建模难度会提高适合进阶折腾。最后再说一个扩展思路。当前模型里代理商只管卖电价格信号是唯一的控制手段。如果小区将来接入光伏储能代理商手里就多了储能充放电这个“调节旋钮”可以把主从博弈扩展成“价格 储能调度”的联合优化下层照样是用户充电决策上层加储能 SOC 状态变量整个框架不用推翻重来只需要在现有 MATLAB 代码上再加一组变量和约束。这也是我认为这套主从博弈模型最具迁移价值的地方——框架稳定扩展省事前期把边界条件和数值细节调稳后面换场景基本就是换参数的事。