去年年底我接了一个评估任务两家充电桩运营商要在同一个区域电网里共同参与电力市场甲方想知道各家电量该怎么报、报完之后会不会出现恶性竞争。老实说纯做单主体的优化调度MATLAB里几分钟就能出结果但一旦引入“互相影响、互相牵制”问题就从优化变成了博弈。为了把这件事验证清楚我在MATLAB里搭了一套两阶段充电桩市场投标测试环境核心是把电动汽车调度潜力揉进投标约束再用纳什均衡作为整个策略链的收敛目标。这篇文章就把从建模、代码架构到调试细节的全部过程写出来适合正在复现同类课题的研究生以及想评估多充电站联合入市方案的工程师参考。1. 调度潜力的量化口径从单辆车到充电站聚合的可调边界1.1 单台EV的可调度域到底长什么样投标量不是拍脑袋定的它受物理边界约束。一台电动汽车插到充电桩上能参与调度的能力由四件事决定电池容量上限、充放电功率上限、接入时段长度、用户离开时要求的最低电量。我习惯把调度潜力想象成一个带时间锁的仓库电池容量是仓库空间充电功率是进货口速度放电功率是出货口速度用户离网最低SOC是“无论如何都要留下的保底库存”。白天车停在那儿不急着走的那几个小时就是仓库能被市场利用的时间窗。在测试环境里单个EV在每个调度时段的可行充放电功率带可以这样估计% get_flex_powerband: 计算单台EV的可行充放电功率带 [kW] % flex_up 0 表示可充电功率上限flex_low 0 表示可放电功率下限 function [flex_low, flex_up] get_flex_powerband(ev, t_horizon) nT length(t_horizon); flex_low zeros(nT, 1); flex_up zeros(nT, 1); for t 1:nT % 当前电量 soc_now ev.soc(t); % 可充电上限受额定快充功率与“别把电池充爆”共同限制 cap_up (ev.e_batt_cap - soc_now) / ev.eta_c / ev.dt; flex_up(t) min(ev.pc_max, cap_up); % 可放电下限受额定放电功率与“离网保底电量”限制 cap_lo (soc_now - ev.e_dep_min) * ev.eta_d / ev.dt; flex_low(t) -min(ev.pd_max, cap_lo); end end这里的e_dep_min是用户设定的离网最低电量比如60 kWh的电池、要求离开时不低于80%那就是48 kWh。eta_c和eta_d分别是充放电效率通常取0.9到0.95别小看这几个百分点的效率它直接影响投标量的经济边界。需要注意这里为了测试环境简洁我用功率带近似可行域。真实工程里还要加温度修正、SOC区间非线性折减、电池健康状态衰减等细节但做博文这套测试已经够了。1.2 从单桩到集群聚合不是简单的加减法单台车算完之后要让整个充电站参与市场还必须把几十上百台EV的功率带做时间对齐后的加总。为什么说“不是简单的加减法”因为每台EV的接入时间、离开时间、初始SOC都不一样。直接用全站额定功率之和当可调度容量结果是严重虚高优化器会把调度潜力用在一个根本不存在的边界上。正确的聚合方式是在每个调度时段只把那些“此刻已在网内”的EV的上下限分别累加for k 1:length(station.ev_list) [lo, up] get_flex_powerband(station.ev_list(k), t_horizon); station.flex_low station.flex_low lo; station.flex_up station.flex_up up; end这个逐时段累加逻辑很简单但很容易踩坑。我得提醒一句聚合后的功率带只是“各行段独立可行”的集合如果把多个时段的电量加起来作为跨时段能量约束还要额外校验能量连续约束。比如某台EV在17:00到18:00可放电50 kWh但它在18:00就要走这时候该时段的放电能力带就不能平移到后面时段。测试环境的简化处理是只允许调度器在EV在线时段内使用其功率带跨时段能量转移严格遵循电池状态递推方程。聚合出的flex_low和flex_up就是两阶段投标模型里最核心的“调度潜力”输入第一阶段用它约束投标量的上限第二阶段用它决定实际功率分配。1.3 调用潜力的代价电池退化成本与用户体验惩罚免费调用电动汽车电池会让模型做出极其激进、实际不可行的策略。用户把车停在充电站不是为了让运营商拿电池去市场套利。所以收益函数里必须有两笔账第一笔是电池退化成本。充放电循环次数和放电深度相关我采用线性近似每释放1 kWh造成约0.05元到0.15元的循环寿命损耗具体系数按磷酸铁锂或三元锂电池经验取值。这个系数在测试环境里做成参数方便灵敏度分析。第二笔是用户体验惩罚。如果为了市场收益车离开时总是贴着最低SOC走用户下回一上车发现续航焦虑下次就不会再来这家站。我在模型里对低于目标SOC的部分施加线性惩罚系数实际离网电量与用户目标电量差得越多惩罚越大。这个惩罚本质上是在收益函数里加一个带符号的二次项保证了内部优化问题仍是凸的后面求解会顺很多。把这三节汇总调度潜力的量化从“单台车物理边界”到“全站可调度功率带”再到“调用成本”形成了一套完整的输入数据。这组数据直接喂给投标模型下一步进入博弈框架。2. 两阶段投标策略的博弈框架参与人、策略空间与纳什均衡落脚点2.1 两阶段的市场时序日前投标与实时平衡测试环境模拟的是最常见的现货市场时序。第一阶段是日前市场每个运营商在T1日到来之前把24小时或96个时段的分时购电计划报给市场。这个计划一旦中标就锁定了该运营商在日前市场的购电量和购电成本。第二阶段是实时市场真实运行当天电动汽车实际接入情况与预估值有偏差运营商需要在实时市场买入或卖出偏差电量。实时价格与日前价格不同而且会受全体市场参与者总偏差的影响这就形成了博弈耦合的关键点。在这个框架下一个“两阶段投标策略”指的是第一阶段做出的投标量要能够最大化期望收益同时充分保留第二阶段用调度潜力去纠偏的空间。换句话说第一阶段的决策不是单纯的“买多少电”而是“买多少电才能让自己在第二阶段有足够的调整空间且不被罚”。2.2 参与人、策略空间与收益函数博弈参与人是同一个配电网区域内的I个充电桩运营商程序里记为nSta。每个运营商的策略是分时投标电量向量[ q_i (q_{i,1}, q_{i,2}, \dots, q_{i,T})^\top ]约束条件是0到聚合调度潜力上限之间[ 0 \le q_{i,t} \le \overline{P}_{i,t} \cdot \Delta t ]收益函数我写成了四段式[ U_i(q_i, q_{-i}) \text{零售收入} - \text{日前购电成本} - \text{实时平衡成本} - \text{调度代价} ]零售收入对外售电单价×实际售电量这是相对固定的部分。日前购电成本日前出清价×投标量。实时平衡成本最关键如果实际用电量超过日前投标量就要以实时价格补购差额反之如果买多了就要以实时价格卖出剩余。这导致每个运营商的实际收益不仅取决于自己报了多少还取决于大家都在实时市场买还是卖。为了在测试环境里保持凸性我把实时价格近似成基准价格加上一个与全体偏差量成正比的线性项[ \lambda_{rt}(\omega) \lambda_{rt}^{base}(\omega) \gamma \cdot \sum_{j1}^{I} (D_j - q_j) ]这里的(D_j)是运营商j在实时阶段的实际用电量(\gamma)是价格弹性系数。这个近似的物理解释是大家都报少了、实时都在补电实时价格会被买压抬高大家都报多了、实时都在卖电实时价格会被卖压压低。正是这个耦合项让不同运营商的投标策略真正变成了一个博弈问题。2.3 纳什均衡作为“无人愿意主动变卦”的判据有了参与人、策略空间和收益函数接下来要回答“最后会稳定在什么地方”。在非合作博弈里最自然的解概念就是纳什均衡。纳什均衡定义为一组策略(q^(q_1^, q_2^, \dots, q_I^))使得任何单一运营商i在其他人策略不变时都无法通过单方面改变自己的投标量来获得更高收益。换句话说均衡状态下每个人手里的策略都是对其他所有人的最优回应谁先偷偷改谁就吃亏。测试环境里我用这个判据还有一个工程原因充电桩运营商之间不存在强制契约谁都可以随时改报价策略。比起“大家一起联合最优”这种需要信任和分钱机制的方案纳什均衡是更现实的预测结果。联合最优解在收益分配谈不拢时根本站不住脚——这个区别后面实验部分还会专门验证。3. MATLAB测试环境架构按“对象-模型-求解-场景”切分的工程实现3.1 为什么用面向对象而不是一摞脚本MATLAB里做优化模型的人常犯一个毛病所有状态都塞进全局变量脚本从上到下跑一遍。这个习惯在一次性分析里没问题但想系统性地测试不同场景、不同运营商数量、不同价格曲线就会变成噩梦。我这套测试环境按“对象-模型-求解-场景”四个层次切分核心是五个类EVFlex描述单台车的动态参数ChargingStation描述充电站的聚合调度潜力MarketModel描述日前与实时市场价格规则NashEqSolver负责外层均衡迭代CaseBuilder负责根据场景配置生成对象。这样每个部件都能单独改、单独测试。3.2 核心类的设计与关键属性拿ChargingStation举例它对外只暴露两个核心数据聚合后的flex_low和flex_up外加一个站级配置batt_degrade_cost。这样的设计让上层博弈求解器不需要关心站内到底停了几台特斯拉还是几台比亚迪它只需要知道“这个站能提供的可调度功率带是多少”。classdef ChargingStation handle properties id % 运营商编号 ev_list % EVFlex对象数组 flex_low % 聚合后可放电功率下限 [kW] flex_up % 聚合后可充电功率上限 [kW] batt_degrade_cost % 每kWh调用对应的退化成本 retail_price % 面向用户的售电单价 end methods function obj aggregate_flex(obj, t_horizon) nT length(t_horizon); obj.flex_low zeros(nT, 1); obj.flex_up zeros(nT, 1); for k 1:length(obj.ev_list) [lo, up] get_flex_powerband(obj.ev_list(k), t_horizon); obj.flex_low obj.flex_low lo; obj.flex_up obj.flex_up up; end end end endMarketModel类负责生成价格曲线和实时清算规则。它的输入是基准日前的分时出清价lambda_da、各时段的实时基准价lambda_rt_base、价格弹性系数gamma输出是给定全体策略后每个人收益所需的成本和偏差。这一类我建议做得越薄越好尽量只做计算不做决策方便后面替换成更精细的市场清算逻辑。3.3 求解器选型内部优化与博弈外层的分工测试环境里有两层求解任务选型逻辑完全不同。外层是均衡迭代我采用“对角化最佳响应”算法理由后面专门讲。内层是单个运营商在固定对手策略下的最优投标问题我一开始用YALMIP建模、调MATLAB自带的quadprog求解因为收益函数取负号后是凸二次规划。后来为了减少第三方依赖直接手写二次规划的quadprog接口调用效果一样部署也省心。为什么不直接用KKT条件把所有运营商的最优性条件联立成混合互补问题MCP一次性求解理论上可以实际上调试痛苦。多运营商场景下KKT系统规模随时段数和运营商数暴涨而且对初始点极度敏感一旦商业求解器没装或者许可证有问题整个测试环境就瘫痪了。最佳响应迭代则不一样每一轮只求解I个中等规模的二次规划即使中途某一轮发散了也能从日志里定位是哪个运营商的子问题出了问题。3.4 代码目录与运行流程目录结构如下目录/文件作用config/场景参数、价格曲线、电池参数等配置文件model/EVFlex、ChargingStation、MarketModel等类定义solver/NashEqSolver、最佳响应求解器、Nash gap计算case/基准场景构建脚本CaseBuilderplot/收敛曲线、均衡结果、最佳响应曲线绘图脚本run_case.m主入口加载场景→求解均衡→输出结果random_restart_test.m多初始点验证均衡稳定性的辅助脚本运行流程上我坚持一个原则run_case.m只做三件事——调用CaseBuilder生成对象、调用NashEqSolver求均衡、调用plot脚本画图。任何想加的新实验场景都写成独立的配置文件而不是在主脚本里堆条件分支。这样跑“运营商数量从2变到6”这类灵敏度测试时只需要改一个参数再重新运行同一个入口。4. 核心求解算法最佳响应迭代、阻尼更新与均衡间隙判定4.1 单运营商最佳响应问题怎么求解在对手策略固定的情况下运营商i需要求解[ q_i^* \arg\max_{q_i} U_i(q_i, q_{-i}) ]我把它整理成标准二次规划形式。目标函数取负号变成凸最小化约束条件都是线性的投标量介于0和调度潜力之间。这里有三个必须处理的细节第一实时平衡成本里包含(q_i)和(q_{-i})的乘积项由于对手策略在外层是常数乘开后对(q_i)来说仍是二次项没问题。第二调度潜力约束在市场化前已经由aggregate_flex算好直接作为ub传给quadprog。第三用户SOC惩罚项会引入关于(q_i)的二次项我把它合并进同一个Hessian矩阵。4.2 外层迭代对角化最佳响应外层算法采用最经典的高斯-赛德尔式对角化最佳响应每一轮里依次更新每个运营商的策略更新第i个运营商时其他运营商一律用本轮最新值。全部更新完后再用阻尼因子做一次松弛。q repmat(q0, nSta, 1); % nSta行每行是一个运营商的T时投标量 alpha 0.4; % 阻尼系数 for iter 1:maxIter q_old q; for i 1:nSta q_other q; q_other(i,:) []; q(i,:) solve_best_response(i, q_other, params, q(i,:)); end % 阻尼更新防止振荡 q alpha * q (1 - alpha) * q_old; % 计算纳什间隙 gap compute_nash_gap(q, params); if gap tol break; end end为什么q_other(i,:) []之后还能正确索引我这份代码里用的是“移除第i行”的直接写法逻辑清楚。如果运营商数量少这样做没问题运营商数量多了建议改用掩码索引避免反复复制矩阵带来的开销。solve_best_response内部就是调quadprog把目标函数、线性约束、上下界组装好function q_i solve_best_response(i, q_other, params, q_init) H build_hessian(i, params); % 凸二次项 f build_gradient(i, q_other, params); % 包含对手策略的线性项 Aineq []; bineq []; lb zeros(params.T, 1); ub params.flex_up(i, :); options optimoptions(quadprog, Display, off, ... OptimalityTolerance, 1e-8); q_i quadprog(H, f, Aineq, bineq, [], [], lb, ub, q_init, options); end这里ub用的是该运营商第一阶段可申报的最大电量由聚合调度潜力决定。如果希望更贴近市场规则还可以把ub设成变压器容量限制和调度潜力的较小值测试环境里我留了这样的参数接口。4.3 收敛判定用Nash gap而不是看迭代差只看策略变化量收敛很容易被假象骗了——策略变化很小不代表每个运营商都已经在自己最优点上。我采用标准的均衡间隙Nash gap判据计算所有运营商“如果单独换成最佳响应能额外多赚多少收益”的总和。[ \text{gap} \frac{1}{I} \sum_{i1}^{I} \left[ U_i(BR_i(q_{-i}), q_{-i}) - U_i(q_i, q_{-i}) \right] ]当gap小于阈值测试环境里取1e-4的相对基准时判定收敛。这个判据的意义在于它直接度量“偏离均衡的诱惑”有多大。gap收敛到0说明没有任何一个运行商有动机单方面改变策略。对应代码也很短function gap compute_nash_gap(q, params) nSta size(q, 1); total_gain 0; for i 1:nSta q_other q; q_other(i,:) []; br_i solve_best_response(i, q_other, params, q(i,:)); total_gain total_gain calc_utility(br_i, q_other, params) ... - calc_utility(q(i,:), q_other, params); end gap total_gain / nSta; end4.4 内层不收敛或非凸时的备选方案如果测试中发现迭代振荡、gap下降不到阈值我会依次检查三件事阻尼系数、初始点、凸性前提。阻尼系数的经验值是0.3到0.5。太大收敛快但容易振荡太小则要跑很多轮。初始点不要全取同一个值否则各运营商策略完全对称迭代很容易停在鞍点附近我习惯从q0 ub * rand(nSta, T)开始。如果模型里加了太复杂的非线性市场清算规则导致内层问题非凸quadprog会直接失败。备选方案有两种一是把实时价格弹性项线性化在均衡点附近做泰勒展开二是改用fmincon配sqp算法求解内层问题外层照常用最佳响应只是每一轮慢一些。测试环境的第一版我坚持保持凸性只有跑通基础流程后才尝试扩展。5. 测试场景设计均衡存在性验证、灵敏度分析与扩展性压力测试5.1 基准场景与关键参数测试环境里的基准场景设置如下参数取值运营商数量2每站EV数量20EV电池容量60 kWh单桩最大充电功率60 kW单桩最大放电功率30 kW用户离网最低SOC0.8调度时段长度15 min日前价格谷0.30元/kWh峰1.20元/kWh实时价格弹性系数gamma0.02元/kWh²电池退化成本0.08元/kWh两组运营商的差异主要在外接EV的到达时间分布A站接入时间集中在9:00-17:00B站集中在17:00-21:00。这样设计是为了让两个运营商的调度潜力曲线形成“错峰互补”博弈会更有意思。5.2 三种模式的结果对比在这个场景下我对比了三种模式第一种是“单边最优”每个运营商完全忽略对手策略只按自己的预测投标。结果是A站拿到约3.2万元日收益B站只有1.8万元看起来A站赚大了但把A站自己的策略代入B站视角检验B站只要把投标量提高15%就能额外多赚约4200元。这说明单边最优根本不是稳定点。第二种是“纳什均衡”两个运营商都按最佳响应迭代收敛后的策略投标。收益分别约为2.71万元和2.66万元。差距很小任何一方单方面改变策略最多只能增加不到300元不足以抵消市场风险和违约风险。这正是测试环境想得到的结论。第三种是“联合最优”假设两个运营商合并成一家所有EV统一调度总收益最高约5.55万元比纳什均衡高出约1800元。如果把收益均分每方能拿2.775万元确实比纳什均衡高。问题在于没有外部约束每一方都有动机虚报成本、在分钱时多拿最后联合方案往往谈崩。这个对比说明纳什均衡不是效率最高的解但它是不需要信任契约时最现实的稳定解。5.3 均衡点的稳定性多初始点与最佳响应曲线验证为了验证求到的均衡点不是“碰巧撞上的”我做了一个最朴素也最有效的检验随机初始化20组策略分别跑完整的最佳响应迭代观察收敛点是否一致。结果表明20组初始点全部收敛到同一个策略附近投标量的最大差异小于0.3%。这给测试环境一个很有力的信心在当前凸性假设下纳什均衡解是数值稳定的。对于两个运营商的情况我还画了最佳响应曲线固定玩家A的策略为横轴玩家B的最佳响应落在纵轴上反过来再画另一条曲线。两条曲线的交点就是纳什均衡。这个可视化对甲方特别有说服力——比直接给一行“gap1e-5”直观得多。5.4 扩展性压力测试测试环境最终要验证的另一个问题是算法规模能不能撑住真实场景。我做了运营商数量2/4/6、EV总规模从40到240的递进测试。运营商数EV总数单轮迭代耗时均衡所需轮数总耗时2400.2 s183.6 s4800.9 s2522 s62403.1 s33102 s单轮迭代耗时基本随运营商数和EV数线性增长轮数增长更温和。6个运营商、240辆车的情况下总耗时约100秒对测试环境完全可接受。如果继续升级到几十个运营商我会把内层quadprog换成并行parfor并且把外层改为异步更新的分布式最佳响应思路不变。6. 调参经验与代码调试中的关键坑位6.1 迭代振荡先查阻尼、再查初始点、最后查数据量纲我第一次跑出来的gap曲线是一个“锯齿形”——下降几步又弹回去完全不收敛。排查后发现不是算法逻辑问题而是实时价格弹性系数gamma设得太大导致收益函数对q_i的二次项曲率过大最佳响应更新幅度远超可行域。解决方式是组合拳把gamma从0.02改成0.005作为基础场景阻尼因子从0.5降到0.35初始点改成各运营商差异化策略。改完之后gap曲线变成单调下降大约15轮内达到阈值。数据量纲的问题也在这里暴露出来收益是元量级SOC是无量纲小数混合建模时如果忘记换算成同一个基准Hessian矩阵里会出现10^4和10^0数量级并存的病态矩阵。我的处理方式是全流程以kWh和元/kWh为基准单位最后展示结果时再折成一天总收益。6.2 双线性项的凸化处理是建模阶段的胜负手两阶段问题里最麻烦的是第一阶段投标量与第二阶段实际电量偏差共同决定实时成本这会产生双线性交叉项(q_i \cdot (D_i - q_i))。虽然(D_i)是一段表达式而不是常数但在测试环境里我做了关键简化将第二阶段内部调度问题的解看作投标量的线性函数。这个简化基于一个合理假设内部调度问题本身是线性规划在最优基点处对(q_i)的灵敏度是固定的。于是实时成本近似为关于(q_i)的二次函数整个内层问题保持凸性。这个技巧在学术文献里叫“线性化伴随约束”在工程博客里说人话就是如果你的调度计划是跟着投标量线性调整的那模型就还能用二次规划一次求出来。6.3 验证均衡时的“多初始点”习惯必须养成只跑一次收敛到gap很小还不够。我在测试环境里写了一个辅助脚本random_restart_test.m专门做多初始点检验。每次运行都记录收敛点、gap、迭代轮数最后对比不同初始点下的策略差异和收益差异。这个习惯帮我揪出过一次问题第三版模型加了用户SOC惩罚项之后初始点取全零和取随机数收敛点差了大约5%。原因是SOC惩罚项让目标函数在可行域边界出现了很小的非凸区域影响了部分子问题的求解。后来我把SOC惩罚项的二次系数调小并给每台EV的e_dep_min加了±2%的随机扰动才重新恢复全局收敛。6.4 关于“能跑通”和“可信”之间的差距测试环境跑通均衡只是第一步可信度要靠后处理补。我每次出结果都会额外做一个单方扰动检验把某个运营商的投标量在上上下下各偏移1%到10%重新计算它自己的收益画出一条“抛物线”确认最优值确实落在均衡点处。这个检验单独写进报告甲方看完基本不会质疑算法有没有算对。另外还有一个容易被忽略的点均衡结果对不同场景参数的敏感性。不能只在基准价格曲线下求一次均衡就下结论。我后来把电价曲线替换成光伏大发日、阴雨天两种场景gamma和退化成本也做了±30%的扫描发现均衡投标量的变化完全在预期范围内这才敢把结论往外放。最后分享一个我做这个项目时最深的体会两阶段博弈模型的工程难点不在于纳什均衡这个概念本身而在于把物理约束、市场规则、经济成本揉进同一个可求解的框架。MATLAB的好处是面向对象和自带优化工具箱能快速搭出完整测试链路坑都在细节里功率带聚合要对齐时段双线性项要尽早凸化数值量纲要统一均衡验证要多初始点。把这套测试环境跑通之后后续换市场规则、换调度潜力模型、换求解器都只是替换对应模块的事。如果读者想把模型再往前推一步我建议可以往两个方向扩展一是把实时价格的不确定性建模成场景集合并加CVaR风险约束二是把单轮迭代改成异步分布式求解去适应更大规模的运营商集群。