主从博弈在电动汽车有序充电中的Matlab实现与调试指南

想象一下这个场景小区地下车库里停着三四十辆电动车晚上六点下班回来插上枪早上七点拔枪走人。表面上看每辆车只是“充个电”但站在小区能源代理商的角度这几十辆车的充电负荷叠加起来完全能把配变容量顶到红线。更麻烦的是所有车主都习惯在晚上七八点插上枪就开始充正好撞在电网晚高峰上代理商的购电成本直接被拉满。怎么定价、怎么引导这些车错峰充电还能保证自己赚钱、用户也愿意配合这就是一个非常典型的双层决策问题——也就是主从博弈。这个课题我前后折腾了小半年从纯理论建模到Matlab跑通踩了不少坑也积累了一套可以复用的实现框架。如果你正在做智能小区、需求响应、电动汽车有序充电方向或者准备拿这个方向做毕业设计这篇文章应该能帮你少走一大半弯路。我不打算讲太多云里雾里的理论推导重点放在模型怎么建、代码怎么写、迭代为什么不收敛、怎么排查这类实操问题上争取你看完就能自己搭出一版能跑的仿真。1. 这个课题到底在做什么先说清楚整个问题的主角是谁。这里其实有四方参与者电网、代理商、电动汽车用户还有小区里的常规负荷空调、照明、冰箱这些。电网只管卖电给代理商执行的是分时电价或者实时电价。代理商夹在中间从电网批量买电再以零售价卖给小区用户赚的是价差。电动汽车用户呢响应代理商的零售电价决定自己什么时候充、充多少功率目标是充电成本尽量低、用车需求尽量满足。常规负荷就当作固定基线不参与博弈。这三方的利益是冲突的。用户希望电价越低越好最好谷电价时段全部覆盖自己充电窗口代理商希望零售电价定得高一点来扩大利润但价格定太高用户会减少充电量或者干脆不在你这里充边际收入反而下降。这种“我决定一个策略你再根据我的策略做最优响应而我预先就知道你会这样响应”的先后决策关系就是主从博弈也叫Stackelberg博弈。用术语说代理商是领导者先公布各时段充电价格用户是跟随者观察到价格后优化自己的充电计划。代理商在做定价决策时必须把用户会怎么反应纳入考量这就要求把下层的用户充电优化问题嵌入上层的定价问题里构成双层结构。这个结构用Matlab实现主要分成三块建模、求解、仿真分析。建模就是把上下两层目标函数和约束写成数学表达式求解是设计算法把双层问题算出来仿真分析是生成用户数据、跑结果、画曲线、验证均衡。后面几章我分别展开。2. 数学建模上下两层各自治什么2.1 上层代理商定价模型代理商在上层的决策变量是未来24小时或者更短的调度周期内的零售电价。假设时间轴按1小时划分成24个时段用 (t1,2,\dots,24) 表示电价记作 (p_t)单位是元/kWh。代理商的目标是最大化自己的净利润也就是售电收入减去购电成本。这里可以加一个运营成本项但一般情况下运营成本相对固定可以省略。表达式是[ \max_{p_t} \sum_{t1}^{24} p_t \cdot L_t(p) - \sum_{t1}^{24} c_t \cdot L_t(p) ]其中 (L_t(p)) 是所有电动汽车在第 (t) 时段的总充电功率它是电价向量 (p) 的函数因为用户会根据电价调整充电计划。(c_t) 是代理商从电网买电的成本电价由分时购电曲线给定。这个目标函数看起来很简洁但它不是普通优化问题因为 (L_t(p)) 无法写成 (p) 的显式解析表达式。它是下层用户优化问题的解。所以上层问题没法直接用梯度下降或者二次规划求解器硬解必须把下层问题一起考虑进来。上层还需要加约束。最常见的约束包括零售电价不能超过某个上限比如阶梯电价的上限不能低于购电价太多否则亏本有时还会加一条总利润不能为负的要求防止出现“赔本赚吆喝”的局部解。这些约束在实现时都写成线性不等式即可。2.2 下层电动汽车用户充电决策下层有多个独立用户每个用户 (i) 有自己的充电需求。假设用户 (i) 的车在某个时段接入电网插枪某个时段离开拔枪接入时间 (a_i) 到离开时间 (d_i) 就是它的可充电窗口。用户 (i) 的决策变量是在每个接入时段 (t \in [a_i, d_i]) 的充电功率 (P_{i,t})单位kW。它的目标函数要考虑两个方面一是充电费用尽量少二是离开时电池电量尽量达到预期目标。这里如果用严格的“满意度”建模会很复杂大多数论文里更常用的做法是在目标函数里加上一个惩罚项用来惩罚实际离开电量与期望电量的偏差。目标函数可以写成[ \min_{P_{i,t}} \sum_{ta_i}^{d_i} p_t \cdot P_{i,t} \Delta t \lambda_i \left( SOC_i^{dep} - SOC_i^{target} \right)^2 ]其中 (\Delta t) 就是时间粒度这里是1小时(SOC_i^{dep}) 是离开时的SOC(SOC_i^{target}) 是用户期望的SOC(\lambda_i) 是用户对电量完成度的重视程度不同用户可以设置不一样体现“有些人怕亏电有些人无所谓”的差异化。约束条件有以下几类充电功率上下限(0 \le P_{i,t} \le P_i^{max})(P_i^{max}) 一般是慢充桩的额定功率7kW是常见值。电池SOC递推方程(SOC_{t1} SOC_t \eta \cdot P_{i,t} \Delta t / B_i)其中 (\eta) 是充电效率(B_i) 是电池容量。SOC不能越限(SOC_{min} \le SOC_t \le SOC_{max})一般SOC范围是0.1到0.9。离开时SOC要满足期望(SOC_i^{dep} \ge SOC_i^{target})。注意这个约束可以替代目标函数里的惩罚项但两者同时用效果更好因为硬约束容易导致“为了达标不惜一切代价”的充电行为在特殊时段可能推高负荷。仔细看下层问题目标函数包括电价乘以功率线性项和SOC偏差平方二次项加上线性约束这是一个标准的二次规划问题用Matlab自带的quadprog就能直接求解。这一点是整个实现里最让人省心的地方因为二次规划求解非常成熟几乎不会遇到算法不收敛的问题。2.3 谁是领导者、谁是跟随者博弈均衡与存在性主从博弈的均衡解就是一组电价 (p^) 和一组用户的充电计划 (P^)使得给定 (p^)每个用户的 (P_i^) 都是各自下层问题的最优解给定所有用户都会按 (P_i^) 响应(p^) 是上层代理商的最优定价。换句话说上层不能通过单方面改变电价变得更有利下层任何一个用户也不能在电价不变的情况下单方面改变充电计划获得更多收益。这正是Stackelberg均衡。理论上只要上层价格可行集是紧的、下层目标函数是凸的Stackelberg均衡的存在性是有保障的。因为下层QP是严格凸的最优响应函数对电价是连续的上层连续函数在紧集上有最优解。这类保证对工程仿真来说已经足够实际代码里更多要考虑的还是迭代算法的收敛和实现稳定性。3. 求解算法设计从纸面到可跑代码这一章是项目里最核心的部分也是很多人从论文走向Matlab时最容易卡住的地方。把双层问题写成硬啃的数学表达式和变成能跑出结果的代码中间隔着一层非常厚的方法论。我介绍三种常见求解思路说说它们各自的适用场景和坑。3.1 方案一KKT条件转单层MPEC这是论文里最“正统”的做法。既然下层是凸QP那可以直接写出它的KKT条件把下层问题变成了上层的约束于是双层问题转化为一个带互补约束的单层数学规划问题也就是MPEC。再通过Fortuny-Amat变换大M法把互补松弛条件转成混合整数线性约束整个问题就变成了一个MILP/MIQP可以交给cplex、gurobi这些商用求解器。这个方案的优点是理论上可以得到精确解不需要调迭代步长也不存在不收敛的问题。缺点是建模工作量非常大。你要手推每个用户QP的拉格朗日函数把KKT条件里的所有互补约束都展开处理大量的辅助变量和线性化。一个用户就有二三十个约束10个用户就是几百条约束写错一个符号排查起来非常要命。而且商用求解器不是免费的很多学生手里只有Matlab自带工具箱这个问题就很现实。我的经验是如果你只是做小规模验证、时间充足可以尝试这条路如果你希望快速复现一个完整案例、还要改参数做对比实验不建议一上来就啃KKT。3.2 方案二双层迭代法最实用的工程方案实际工程项目里绝大多数人会选择迭代法。思路非常直接代理商先给定一组初始电价用户根据这个电价各自求解自己的充电QP得到充电计划代理商统计所有用户的充电负荷计算自己的利润代理商根据利润信号调整电价比如梯度上升或价格搜索重复2到4步直到电价和负荷不再变化。这个方案的好处是每个子问题都很简单。下层是若干个小QP用quadprog秒解上层本质上是一个无梯度或者梯度数值化的优化问题用最直接的搜索方法也能跑通。核心代码量可以压缩到两百行以内。缺点是顶层迭代的收敛性全看细节处理。最常见的坑是价格振荡电价一会涨一会跌系统在那来回震荡就是不收敛。解决办法我在第五章细讲。我在实践里推荐的迭代框架是内层用matlab的quadprog解用户充电计划外层用带自适应步长的梯度上升或者简单的二分搜索更新电价。这种结构稳定、好调试、也方便扩展比如把15分钟粒度改成5分钟粒度。3.3 几种方案的对比与工程取舍为了让你快速选型我把三种方案整理成一张表方案精度建模复杂度计算速度商用求解器依赖推荐场景KKT转MPEC MILP精确高中需要论文理论证明、小规模算例双层迭代 QP近似低快不需要工程仿真、快速复现、做对比实验元启发式 QP近似中慢不需要上层函数难处理、需要全局搜索我在这个项目里跑通的是第二种方案。说实话对于智能小区代理商定价这类问题上层变量是24维电价向量下层是若干个QP迭代法的精度已经足够支撑分析结论了。没有必要为了“精确”去背负一堆MPEC的建模包袱。4. Matlab实现的关键细节Matlab这套代码的完整逻辑我会按“场景生成 → 下层QP求解 → 上层价格更新 → 结果输出”四步来说。这也是照着写就能跑通的顺序。4.1 场景生成50辆EV的一条龙参数设计做博弈仿真不能只用三五辆车太少体现不出负荷聚合效应太多代码又调得慢。我建议用50辆EV跑一遍全24小时场景速度比较合适。每辆EV的参数包括接入时间、离开时间、初始SOC、期望SOC、电池容量、最大充电功率。这些参数不能拍脑袋定要有一定的统计依据。我的做法是用正态分布随机生成并用rng固定随机种子保证每次实验可复现。rng(2024); N 50; % 电动车数量 H 24; % 调度时段数 dt 1; % 时间粒度1小时 EV.arrive round(normrnd(18, 2, N, 1)); % 接入时间均值18点标准差2小时 EV.leave round(normrnd(7, 1, N, 1)); % 离开时间均值7点 EV.leave(EV.leave EV.arrive) 8; % 防止时间倒挂 EV.soc_init min(max(normrnd(0.3, 0.1, N, 1), 0.1), 0.6); % 初始电量 EV.soc_target min(max(normrnd(0.9, 0.05, N, 1), 0.7), 1); % 期望电量 EV.capacity normrnd(60, 5, N, 1); % 电池容量kWh EV.pmax 7 * ones(N, 1); % 最大充电功率kW EV.eta 0.9; % 充电效率这里有几个细节要提醒。第一生成时间后一定要检查逻辑错误比如离开时间早于接入时间这种粗心错误会让后面所有索引都错乱。第二normrnd生成的SOC可能超出合理范围要记得clip到物理可行区间。第三电池容量建议在55到65之间浮动所有车都用60反而显得不真实。4.2 下层QP建模用quadprog还是YALMIP下层用户充电优化是个二次规划Matlab下面最好用的是quadprog。因为下层每个用户规模不大不需要装YALMIP直接用quadprog就够了。YALMIP在求解MILP、需要调用gurobi那些场景才有优势这里用反而是过度设计。quadprog的标准形式是[ \min_x \frac{1}{2} x^T H x f^T x ]对用户 (i) 来说决策变量 (x) 是它可充电时段内的充电功率向量。假设它有 (m_i) 个可充电时段则 (x) 是 (m_i \times 1) 的向量。目标函数需要处理两部分充电费用这是线性项SOC偏差惩罚这是二次项。把SOC偏差用线性变换代入功率变量可以整理出H和f。如果你不想手推矩阵也可以用YALMIP的符号建模牺牲一点效率换建模清晰度。手写quadprog矩阵有个地方特别容易翻车SOC递推方程里的系数矩阵。因为SOC是分段累加充电功率的效果约束矩阵是一个下三角形式但你要是把时段索引搞错约束就变成了“看到未来”。这里给一个YALMIP写法的参考虽然实际运行稍慢一点但它不容易错用来对照验证是很方便的x sdpvar(m, 1); soc soc_init; for t 1:m soc soc eta * x(t) * dt / capacity; constraints [constraints, soc soc_max, soc soc_min]; end constraints [constraints, soc soc_target, 0 x pmax]; objective price_win * x * dt lambda * (soc_target - soc)^2; optimize(constraints, objective);你只需要知道价格向量price_win怎么映射到用户的充电窗口上。很多新手卡在这一步原因是没搞清楚“用户的可充电时段”和“全局24时段”之间的索引关系。我建议单独写一个函数把全局时段索引转换成语义清晰的充电窗口索引宁肯多花十分钟写好这个映射关系也不要在一个几百行的主脚本里硬找索引。4.3 上层电价更新的迭代框架上层电价更新是整段代码最需要耐心调试的部分。我最终采用的是一个带约束的梯度上升法。代理商利润对电价的梯度理论上是负荷加上价格对需求的影响项但这个影响项没法解析得到所以实际代码里我用数值差分来估计梯度。也就是把 (p_t) 微调一下看利润变化量再决定往哪个方向调。% 初始化价格 price ones(H, 1) * 0.8; % 初始统一价 price_lb 0.3 * ones(H, 1); % 价格下限 price_ub 1.5 * ones(H, 1); % 价格上限 step 0.05; tol 1e-4; for k 1:200 % step1: 已知price求解所有用户充电计划 load_total zeros(H, 1); for i 1:N x_opt solve_user_qp(price, EV(i)); load_total load_total map_to_global(x_opt, EV(i)); end % step2: 计算当前利润 profit(k) sum((price - buy_price) .* load_total) * dt; % step3: 数值梯度 grad zeros(H, 1); for t 1:H price_eps price; price_eps(t) price_eps(t) eps_p; % 重新求解所有用户的QP计算新价格下的利润 load_eps zeros(H, 1); for i 1:N x_eps solve_user_qp(price_eps, EV(i)); load_eps load_eps map_to_global(x_eps, EV(i)); end profit_eps sum((price_eps - buy_price) .* load_eps) * dt; grad(t) (profit_eps - profit(k)) / eps_p; end % step4: 梯度上升 投影到可行区间 price_new min(max(price step * grad, price_lb), price_ub); % step5: 收敛判断 if norm(price_new - price) / norm(price) tol price price_new; break; end price price_new; end这段代码逻辑上是通的但效率不高因为每更新一个 (p_t) 都要重新求解所有用户的QP。我在实际实现里还会加一个批量并行层把不同时刻的数值梯度计算用parfor替代for速度能有明显提升。另外电价更新步长的选择很考验手感我一般从0.05开始如果发现振荡就不断减半。4.4 结果展示与画图仿真跑完重点输出三类图第一是电价收敛过程。横轴是迭代次数纵轴是各时段电价或某几个典型时段的电价值这个图能直观看到算法是否收敛、有没有振荡。第二是最优电价和充电负荷的时间曲线。把代理商最终给出的24小时电价曲线、电动汽车总充电负荷、常规基础负荷画在一张图里对比峰谷时段能否有效错开。第三是SOC箱线图。把50辆车离开时的SOC分布画成箱线图检查是否都达到期望值、有没有极端低的离群车辆。这类图在实验报告和论文里非常出效果也能验证模型约束是否被正确执行。画图代码建议用subplot排布三张图一页保存成PDF矢量图后续写论文直接用。5. 常见问题与调试记录5.1 迭代振荡问题这是最常见的故障表现是价格在收敛区间来回跳动利润曲线像锯齿一样。出现这种问题十有八九是上层步长太大或者电价上下限范围太宽没有物理约束。我的处理思路是分级排查。首先把步长减半看是否改善如果减半后振荡减轻说明是步长问题就逐步调试出合适的步长。其次是做价格平滑每次迭代不是直接采用新价格而是新旧价格加权平均price_new alpha * price_candidate (1 - alpha) * price_old;这里的alpha取0.3到0.5之间能让迭代过程稳定很多。付出代价是收敛速度下降但工程上完全可接受。还有一种振荡来自下层用户数量带来的非光滑响应。想象一下电价只上涨了一点点但用户对这个价格的反应是“我干脆从晚上8点改到凌晨2点充”这是一个跳变的反应导致负荷突变、利润突变、梯度突变。这种非光滑行为不是步长能解决的得靠加平滑项或者使用用户需求曲线的分段线性近似来缓解。5.2 quadprog报错和约束矩阵排查quadprog报错百分之九十是约束矩阵维度不对。我在调试时习惯把每个矩阵的size单独打印出来检查disp([Aeq size: , num2str(size(Aeq)), x_size: , num2str(size(x))]);还有一个容易被忽略的点quadprog默认要求H是半正定。如果你的QP目标函数形式给得不恰当H可能出现负特征值导致求解器直接拒绝求解。如果出现这个报错优先检查SOC惩罚项的二次项系数符号这一项很容易因为手写公式时符号错误变成负的。5.3 计算慢和瓶颈定位50辆车、24个时段、200次迭代如果用数值梯度法内层循环次数非常惊人外层200次迭代每次梯度更新要H次重新求解所有用户QP就是 (200 \times 24 \times 50 240000) 次QP求解。这个规模在普通笔记本上可能要跑很久。实际优化思路有几个下层QP不是每次都需要完整初始化所有变量。复用一个计算好的初始点可以省掉一部分求解时间。引入并行计算。Matlab的parfor在用户循环和梯度循环都可以并行前提是每个求解函数里尽可能减少共享变量。减少外层迭代的上限。实际收敛一般二十次左右就差不多了不用硬跑200次。最省事的是把数值梯度改成解析近似或者用无梯度方法。我的经验是用一次迭代采集整个利润曲线的走势来判断各时段电价该升还是该降虽然粗糙但收敛非常快。5.4 结果验证的几个土办法算完博弈均衡怎么确认结果真的可信我一般用三个土办法验证。第一看下层KKT残差。既然下层是凸QP直接调用quadprog得到的解就应该满足最优性。用Matlab的输出信息exitflag检查一下当exitflag为1时说明求解正常。第二上层用多起点验证。跑三次以上每次给不同的初始电价看最终收敛的均衡价格是否一致。如果三次都一样说明至少是局部最优点如果差别很大说明上层问题有多个局部最优需要进一步用随机搜索找全局最优。第三做参数敏感性分析。比如把用户的 (\lambda) 值整体调大也就是大家更在乎电量充足度理论上最优电价应该上涨因为用户对电价的敏感度下降了如果仿真结果方向相反大概率是模型哪里写错了。这几个土办法不需要额外工具纯靠Matlab现有反馈就能实现强烈建议每个项目都过一遍。6. 个人经验与教训这个项目做下来我对主从博弈在能源管理里的定位有了一些更新的认识。很多人一听到博弈就觉得很玄其实从工程实现角度讲主从博弈更像是一个“嵌套决策”的分层求解框架。代理商先走一步用户随后最优响应两边通过价格和负荷互相耦合。想清楚这一点Matlab代码的架构就自然地出来了外层价格循环、内层用户QP、中间用迭代串起来。如果让我重头再做一次我会在最初就给下层用户分类而不是所有50辆车一把抓。实际场景里用户充电习惯可以归类为“下班即充型”、“夜间错峰型”、“紧急补电型”等几类同类内部参数差异不大聚合后再博弈计算量可以降一个数量级而且结论反而更容易解释。这也是后续可以做扩展的方向比如把不同的用户群视作下层多个子博弈者代理商对他们制定差异化套餐电价。另外还有一个很实用的小技巧在迭代过程中把每一轮的电价向量、负荷向量、利润值全部存下来画成动态图或者做回放。这比只看最终结果要直观得多特别是在给别人讲方案的时候非常加分。你只要在循环体里加三行赋值语句最后用animatedline或者记录数组画出来就行。最后说一句实在话主从博弈模型能不能产出有价值的结论关键不在博弈求解器跑得有多精而在上层的定价目标函数和下层的用户效用函数建得是否符合实际。函数建得贴近真实行为哪怕迭代法近似求解也有工程指导意义函数建得天花乱坠但脱离实际再精确的均衡解也只是数字游戏。这一点无论是做仿真还是将来做真实项目的落地都值得时刻记住。