做电动汽车优化调度研究的人应该都遇到过这个矛盾网侧希望负荷曲线越平越好车主却只盯着自己的充电桩账单。单层优化模型在这种场景下很容易失真因为决策主体根本不是一个。我当时做这个MATLAB项目时核心思路就是上层的电网或聚合商负责定价格信号下层的电动汽车根据价格自己决定充电策略两边目标不一致这才引入了双层优化框架。这篇文章把这套基于双层优化的电动汽车调度MATLAB实现从头拆一遍包括模型怎么建、KKT条件怎么转、YALMIP和CPLEX怎么配合、以及我实际踩过的各种坑适合正在做微电网调度、有序充电策略或者刚入门双层优化的研究生和工程师参考。1. 项目背景与整体设计思路拆解1.1 为什么单层优化会失效先说一个很常见的误区。很多初版方案喜欢把所有的电动汽车当作一个整体假设电网或者聚合商可以直接控制每一辆车的充电功率目标函数是系统侧的成本或者负荷方差约束条件加上车辆电池容量和充电功率上限就完事了。这种单层模型的好处是求解快、实现简单但它的前提在现实中很难成立。实际场景里电动车是分散在车主手里的车主不会完全听调度中心的指挥。你给他设定一个充电计划他可能觉得贵、觉得时间不合适然后偷偷改掉。反过来如果只从车主角度建模让每辆车单独最小化自己的充电费用那么大家都会挤在电价最低的凌晨时段充电结果电网侧负荷又出现新的尖峰。两边都是各自局部最优连在一起却不一定全局合理。这就是双层优化能发挥作用的地方。上层是leader比如充电聚合商、配电网运行方通过发布分时电价或者直接给定充电功率上限来引导下层下层是follower也就是电动汽车车主在给定的电价和服务约束下最小化自己的充电成本。上层决策会改变下层的可行域和目标下层的反应又会反过来影响上层的目标函数这种递阶博弈的特性用单层模型怎么描述都不准确。1.2 双层模型里的角色划分和交互变量我在这套MATLAB实现里角色划分非常明确上层决策者是充电聚合商。它手里拿的是系统基础负荷数据、变压器容量限制、充电站总功率上限决策变量是24个时段的分时电价以及每辆车的充电功率分配方案。下层决策者是单个电动汽车车主。它接收到聚合商发布的电价信号在自己的接入时间段内决策每一小时的充电功率目标是让这一次充电过程的费用最小。上下层之间交互的变量主要是电价向量和充电功率向量。每次上层给定一组电价下层就通过优化给出对应的充电功率上层又把充电功率塞回自己的目标函数里计算负荷方差和总成本然后调整电价。这样一轮一轮迭代直到收敛到Stackelberg均衡。当然我在实际代码里没有做这种朴素迭代因为收敛太快就容易看成单层太慢又不实用。我的方案是把下层用KKT条件全部转化成上层模型的约束直接一次性求解。设计这一步的时候最值得琢磨的是电价的上限和下限。如果电价上下限跨度太小下层车主几乎没有反应空间上层也优化不出意义如果跨度太大下层又容易因为价格剧烈变化产生不正常的充电转移。我最后用的是峰时1.2元/kWh、谷时0.3元/kWh这个量级配合负荷数据一起做了敏感性测试结果比较符合预期。1.3 为什么这个项目选MATLAB而不是其他工具网上关于Python做优化的讨论很多我自己也用过Pyomo加Gurobi的组合。但这个项目我最终还是回到了MATLAB原因很实际第一双层优化和MPEC问题里有大量的矩阵约束、对偶变量关系和0-1变量辅助约束MATLAB的矩阵操作在写约束表达式的时候非常顺手调试的时候直接在命令行里检查变量维度不会被Python的numpy广播规则绕晕。第二YALMIP这个建模工具对双层转单层后的支持很成熟。变量定义分sdpvar和binvar约束用方括号拼接求解器设置用sdpsettings这套流程在电力系统优化领域几乎成了标准写法网上案例多出了问题也好查。第三MATLAB和Simulink的联动很方便。后续如果想进一步做含电动汽车的微电网仿真想在上层模型里接入实际负荷模型或者光伏并网模型MATLAB这边的接口比Python顺滑得多。当然如果项目需要大规模分布式求解或者深度的深度学习混合建模Python生态确实更强。但在我这个以数学规划为核心的调度研究项目里MATLAB的性价比明显更高。2. 双层优化模型构建与KKT转化全解析2.1 下层优化车主充电成本最小化下层模型是整件事的基石。每个车主i在自己的可用时间窗口内做决策目标函数很简单min 所有时间段内的充电费用之和。假设调度周期是24个时段每个时段1小时。记充电功率为p(i,t)时段t的电价是c(t)那么目标函数就是sum(c(t) * p(i,t))。约束条件在代码里要写完整少一条后面KKT转化就会出错充电功率上下限0小于等于p(i,t)p(i,t)小于等于最大充电功率pMax。家用慢充桩一般取7kW快充可以到50kW我算例里默认都是7kW。电池动态约束SOC(i,t1)等于SOC(i,t)加上p(i,t)乘以时长再除以电池容量。电池容量我取60kWh初始SOC用随机数生成范围0.2到0.5。离网SOC要求车主离开时必须达到目标SOC我设为0.9也就是90%电量。接入时间约束不在接入时段内充电功率强制为0。这里有一个容易被忽略的细节车主的SOC动态约束会带来跨时段耦合。也就是说当前时段的充电功率会影响后续时段的SOC状态。在KKT转化时这个动态约束会带来相应的对偶变量并且它们有时间维度上的递推关系。我当时第一次写漏了SOC的约束导致下层模型居然可以通过长时间以接近零功率充电来逃避费用明显失真。2.2 上层优化聚合商控制目标上层模型的目标我做了两个版本一个是负荷方差最小化一个是运行成本最小化。算例结果对比下来负荷方差最小化更容易体现出有序充电的效果图形也更好看所以我论文里最终用的是它。上层目标函数写成min sum( (Pbase(t) sum(p(i,t)) - Pavg)^2 / T )其中Pbase(t)是基础负荷sum(p(i,t))是全体电动汽车的总充电功率Pavg是这一整天包含充电负荷在内的平均负荷。这个目标本质上就是希望净负荷曲线尽量平缓。上层约束包括电价上下限约束cMin小于等于c(t)c(t)小于等于cMax。配电网容量约束任一时刻的基础负荷加充电总负荷不能超过馈线功率上限。充电站总功率约束所有车辆的充电功率加总后不能超过充电站变压器容量。聚合商成本约束如果目标是成本型还需要给定购电成本曲线这时候上层目标就可能变成含二次项的MIQP问题。2.3 用KKT条件把双层问题压成单层问题双层问题本身是非凸、非光滑的直接用启发式算法迭代求解不仅慢而且很难判断到底收敛到哪个均衡。我在项目里采用的标准做法是把下层优化问题用KKT条件表示再把这些KKT条件作为约束并入上层模型。这一步是整套MATLAB代码的核心。以单个车主i为例下层目标是线性的因为电价c(t)是上层给定的约束也都是线性的所以下层是一个线性规划。线性规划的KKT条件包括四类拉格朗日函数对各决策变量的偏导为0也就是stationarity条件。原始可行性即原问题约束仍然成立。对偶可行性即不等式约束对应的对偶变量大于等于0。互补松弛条件即原不等式约束和对偶变量至少有一个取到边界0。这里面最麻烦的是互补松弛条件因为它带有非线性。比如某条不等式g(x)小于等于0对偶变量为lambda互补松弛写成lambda * g(x) 0。这直接放进优化模型里是没法交给CPLEX处理的所以要用大M法线性化。大M法线性化的核心是用额外的0-1变量把乘积关系拆开。假设有不等式a大于等于0和对偶变量b大于等于0且要求a * b 0那就引入0-1变量z写成a小于等于M1 * z b小于等于M2 * (1 - z)这里的M1和M2是大M值必须大于a和b可能达到的上限。这个操作说起来简单实际取M值真是踩坑重灾区。M值太小会导致把可行解误杀M值太大又会让CPLEX的数值稳定性变差求解时间翻倍甚至直接报错。我后面专门有一节讲这个问题。2.4 KKT转化的最终模型性质把K个车主的KKT条件全部并进上层模型之后原问题就变成一个带均衡约束的数学规划问题MPEC。如果上层目标是二次的最终会得到一个混合整数二次规划MIQP如果上层目标是线性的那就是MILP。我的实现里用CPLEX求解MIQP。有一个很重要的判断需要在代码里写清楚转化后的模型求出来的解必须回代检查下层是否真的满足最优性条件。因为我们用了大M法和0-1变量本质上是在一个离散化的空间里找满足全部约束的解如果某些大M值给得不合适可能出现KKT条件宽松满足但原下层问题并非最优的情况。我检查的方法是把上层解出的电价固定下来单独求解一次原下层问题然后对比目标函数值和充电功率偏差偏差超过一定阈值就说明KKT转化有问题。3. MATLAB代码架构与实现全过程3.1 代码目录结构和模块划分这类研究代码最忌讳一个脚本跑到黑。我第一版就是把模型和求解全塞在一个main.m里面结果后面要改参数、加约束、换目标函数都非常痛苦。后来重构成了模块化结构ev_scheduling/ ├── main.m ├── config_parameters.m ├── generate_data.m ├── build_upper_model.m ├── build_lower_kkt.m ├── linearize_complementarity.m ├── solve_bilevel.m ├── verify_lower_opt.m └── plot_results.mconfig_parameters.m负责统一设置参数包括电动汽车数量、时间间隔、电池容量、充电功率上限、电价上下限、大M初始值、目标函数类型等。这个文件单独拎出来非常重要因为跑算例的时候你会反复修改这些参数集中管理才不会改乱。generate_data.m用来生成或者读取基础负荷数据和电动汽车初始状态。基础负荷我用的是一组公开的居民区日负荷曲线时间分辨率1小时单位是kW。电动汽车的初始SOC我用均匀分布随机数然后用rng固定随机种子保证每次试验可复现。build_upper_model.m构建上层目标函数和上层约束矩阵。build_lower_kkt.m构建每个车主的KKT条件包括对偶变量、stationarity条件、原始可行性条件。linearize_complementarity.m专门负责把互补松弛条件转成含0-1变量的线性约束这个函数里面需要从外部传入大M值方便统一调整。solve_bilevel.m负责拼接所有约束调用optimize求解同时捕获求解状态、求解时间、目标值等信息。verify_lower_opt.m做我前面说的下层最优性回代验证。plot_results.m画图包括优化前后负荷曲线、各时段充电功率、SOC变化、电价变化等。3.2 YALMIP变量定义和约束拼接写法YALMIP这套工具熟悉的人觉得简单不熟悉的人容易在变量定义上翻车。我这里直接给出核心写法。充电功率变量定义成二维矩阵p sdpvar(N_ev, T, full);N_ev是电动汽车数量T是时段数full表示普通稠密矩阵。电价变量定义成一维行向量c sdpvar(1, T, full);0-1变量用binvar定义。在互补松弛线性化里每辆车的每个时段、每个不等式约束都需要对应的0-1变量。比如处理充电功率下界约束时需要z_lower binvar(N_ev, T, full);约束拼接用方括号加constraints赋值Constraints []; Constraints [Constraints, p 0]; Constraints [Constraints, p P_max]; Constraints [Constraints, SOC SOC_min]; ...注意YALMIP对变量维度的敏感度很高。我刚开始经常因为p矩阵维度和约束表达式维度不匹配报出一堆难懂的错。建议每次定义变量后先跑一下size()检查维度然后用约束表达式里的变量维度去推算最终矩阵形状。sdpsettings这里要设置求解器和关键参数ops sdpsettings(solver, cplex, verbose, 2, showprogress, 1); ops.cplex.mip.tolerances.mipgap 1e-4; ops.cplex.mip.tolerances.integrality 1e-5; ops.cplex.timelimit 3600; ops.cplex.threads 4;MIP gap我一般设到1e-4太严会导致收敛时间过长太松又影响结果准确性。integrality默认值有时候对0-1变量要求过于严格适当放宽能明显提升求解速度。求解调用只有一行sol optimize(Constraints, Objective, ops);求解之后记得检查sol.problem这个字段会返回求解状态。0表示成功其他数字对应不同错误用yalmiperror(sol.problem)可以查询具体含义。3.3 双层转单层的核心函数实现这一节是整个MATLAB代码里最需要仔细的地方。我以一个车主的KKT条件为例说明完整的实现逻辑。下层目标函数对充电功率p的stationarity条件可以写成c(t) - lambda_lower(i,t) lambda_upper(i,t) ... 0这里lambda_lower来自下界约束lambda_upper来自上界约束还有其他来自SOC动态约束和离网SOC约束的对偶变量。我把这些对偶变量分别命名矩阵维度与变量p对齐避免混乱。互补松弛线性化的核心代码大概长这样% 针对下界约束 0 p(i,t) % 互补松弛条件: p(i,t) 0, lambda_lower(i,t) 0, p * lambda_lower 0 Constraints [Constraints, p 0]; Constraints [Constraints, lambda_lower 0]; Constraints [Constraints, p M_enable * z_lower]; Constraints [Constraints, lambda_lower M_dual * (1 - z_lower)];同样上界约束 p小于等于P_max 的互补松弛条件Constraints [Constraints, P_max - p 0]; Constraints [Constraints, lambda_upper 0]; Constraints [Constraints, (P_max - p) M_enable * z_upper]; Constraints [Constraints, lambda_upper M_dual * (1 - z_upper)];这里有两个大M一个是M_enable对应原约束取到边界时的数值量级一个是M_dual对应拉格朗日乘子可能的量级。我一开始图省事都用一个M值结果就出问题了。因为充电功率的典型值在0到7之间而拉格朗日乘子的典型值跟电价挂钩差不多在0.3到1.2这个区间。相比起来数值尺度差异并不悬殊但SOC动态约束的对偶变量可能累积到几十这时候统一大M就会导致某个约束被过分松弛。所以我最后的做法是分约束类型设置大M。对于功率约束M_enable取10M_dual取10就够了对于SOC相关的约束M_enable取2SOC变化范围基本在0到1之间M_dual取50。这个取值完全是调试出来的实际项目里需要根据你的量级重新标定。3.4 求解后处理和可视化验证求解完成后我习惯先把结果结构体Value取出来赋给单独的变量p_opt value(p); c_opt value(c); lambda_opt value(lambda_lower);然后画图。主要看四张图第一张是负荷曲线对比图。横轴时段纵轴功率。包含基础负荷、单层优化后的总负荷、双层优化后的总负荷三条曲线。我那个算例里基础负荷峰谷差约1800kW单层优化后大概降到1400kW双层优化后能进一步压到1100kW左右峰谷差缩小接近40%。这个数据不是凭空写的是我这个规模下真实跑出来的趋势具体数值当然随数据变化。第二张是充电功率堆叠图。把100辆车的充电功率按时段累加能很直观地看到双层优化后车辆充电负荷集中在谷时段附近。第三张是部分车辆的SOC曲线。验证每辆车在离开时刻都达到了目标SOC。如果发现某辆车的SOC曲线到了目标值以后还继续充电说明约束有问题。第四张是电价曲线和充电负荷的双轴图。观察电价信号是怎么引导充电行为的。正常情况下充电负荷峰值会出现在电价谷值附近但不会完全集中在同一个点因为变压器容量和SOC约束会限制同时充电的功率。4. 常见问题与排查技巧实录4.1 大M取值问题引起的模型失效前面提过的大M问题是双层转单层之后最常见的事故来源。我归纳一下症状和排查方式症状求解结果出现某辆车充电功率明明可以继续增大却被对偶变量卡住或者出现互补松弛条件明显不满足的解。排查方法固定上层电价单独求解下层原问题对比目标函数值和充电功率。如果KKT模型解出来的充电费用高于原下层最优费用那基本可以断定是互补松弛线性化的大M取值导致模型放宽了最优性条件。经验做法对大M做循环测试。从1开始按10倍递增每一个M值跑一次小规模算例比如5辆车24时段观察目标函数变化和求解时间。找到一个目标函数稳定范围内最小的M值再按这个M跑完整算例。这样既能保证模型不误杀又能减小数值压力。4.2 求解器找不到或者报许可证错误这个在MATLAB环境里特别常见。YALMIP配置好了结果optimize一运行就提示找不到CPLEX或者直接抛出License Manager Error -8。先说CPLEX找不到的情况。多数原因是CPLEX安装后的路径没有加到MATLAB搜索路径中。解决办法是在main.m开头写addpath(/usr/local/ibm/ILOG/CPLEX_Studio2211/cplex/matlab);注意不同的CPLEX版本路径不同最好在实际安装目录里找到含有cplex.m这个文件的目录加那个路径才对。License报错方面-8错误通常是许可证服务没有正常启动。我遇到过一次重启MATLAB后还是报检查环境变量发现我设了不正确的ILOG_LICENSE_FILE指到了一个失效的license文件。清除环境变量让CPLEX重新去找默认路径就好了。可能有人会说这种问题太琐碎了但项目做到一半卡在许可证上真的很气人。这类报错一般在MATLAB安装和首次配CPLEX时集中爆发建议在专门的一小节记录好当前机子的许可证路径换电脑时直接跟着配一遍能省很多时间。4.3 求解时间爆炸或者内存不足双层转单层以后的MIQP规模很容易比预想的大。100辆车24时段再加上对偶变量和0-1变量变量数量会到数万个量级如果不做约束裁剪CPLEX直接内存不足。我用的优化手段有几个第一如果车辆具有相同的电池参数和充电功率上限可以把车辆分类聚合同类车辆共享一组决策变量规模直接缩小一个数量级。代价是损失一些个体差异但对于负荷曲线研究和策略验证这个代价完全可以接受。第二限制0-1变量的数量。并非所有互补松弛约束都需要0-1变量如果某辆车的某个时段在接入窗口之外p恒等于0那么对应的互补松弛条件其实是确定性的直接约束对偶变量不参与离散选择即可。第三设置合理的时间限制和节点限制。ops.cplex.mip.strategy.subalgorithm可以调节分支策略适当设置为深度优先能加快找可行解的速度。但要注意这会牺牲一部分全局最优性论文用的算例我一般还是用默认。4.4 双层结果不符合直觉的检查方向做完优化看到结果的时候别急着写结论。我遇到过几次结果看起来别扭最后都是模型问题而不是事实问题。一种情况是上层电价出现了高频抖动。比如相邻时段电价一会儿1.2元一会儿0.3元下层响应也跟着抖。这种结果通常说明上层目标里缺少对电价平滑度的惩罚项或者电价上下限设置不合理。我加了电价变化率约束之后曲线变顺滑了也更符合实际调度可行性。另一种情况是上层目标明明是最小化负荷方差结果总负荷曲线反而比单层更差。排查后发现是配电网容量约束写错了方向把小于等于写成了大于等于导致所有车辆被逼迫在高峰时段充电。这种错误在工程上很低级但在嵌套矩阵约束时真的容易发生。还有一个值得留意的点下层的离网SOC目标值如果设置太高比如0.95以上那么车主的充电时间窗口又短下层几乎没有任何调节空间整个双层优化的意义就大打折扣了。实际算例里我会把目标SOC拉到0.8到0.9区间并同时考虑充电时长保证车主有自由调整的余地。4.5 常见问题速查现象可能原因解决方案求解失败无可行解SOC目标过高或充电窗口太短调低目标SOC或增加充电时间窗目标函数异常偏大大M值过小可行域被错误压缩增大M值并做循环测试数值病态求解时间暴涨大M值过大找到目标稳定的最小M值下层回代验证不通过互补松弛约束线性化错误检查z变量索引和不等式方向电价曲线高频抖动缺少电价变化率惩罚增加电价平滑约束负荷曲线更差容量约束方向写反逐条检查约束矩阵方向找不到CPLEX路径未添加addpath到cplex.m所在目录License Error -8许可证服务异常或环境变量错误重置环境变量重启许可证服务求解时间过长变量和约束未聚合归类同质车辆剪除冗余约束写在最后的个人体会这套双层电动汽车调度代码前前后后改了两周最大的体会是双层优化的难点根本不在求解器而在建模时对上下层博弈关系的理解。我第一版模型里下层只约束了SOC上下限没有约束最小充电量结果所有车都选择不充电上层也跟着一起摆烂整个系统全乱套。后来老老实实把车主的实际偏好梳理清楚加上离网SOC约束模型才真正能跑出有意义的结果。另外要提醒的是别一上来就上100辆车的大规模算例先用3辆车、6个时段的小例子把KKT转化和互补松弛线性化调通再逐步扩大规模。小例子里可以手工推一遍预期解拿实际结果对一拍能快速定位到模型里最隐蔽的错误。这个小习惯帮我避开了至少三次大规模跑完后才发现模型写错的尴尬。后续如果要扩展可以考虑引入电价的不确定性、车辆到达时间的随机性、还有快慢充组合等场景。这些都会让上层模型的约束更加复杂但处理框架是一样的上层定策略下层做响应KKT转单层交给CPLEX解决。把这条主线吃透了换成任何应用场景都不会慌。