最近我在做电动汽车有序充放电的仿真研究被问得最多的一个问题就是用MATLAB的Yalmip建模、CPLEX求解时最小化总负荷峰谷差这个目标到底该怎么落地表面上看就是一句 min(max - min)但真正把模型写出来才发现到处都是细节——SOC的单位换算、充放电效率怎么进递推公式、max和min怎么转成CPLEX能处理的约束、以及最坑的充了马上放伪循环。这篇博文就把我完整跑过的一条路线摊开讲从目标函数数学建模到Yalmip代码实现到CPLEX求解调参再到结果校验和常见坑的排查。适合刚接触有序充电仿真的同学也适合那些已经能跑简单线性规划、想往混合整数规划再走一步的人。1. 写目标函数前先把最小化峰谷差翻译成数学语言1.1 峰谷差为什么值得作为优化目标峰谷差指的是一个调度周期内总负荷的最大值减去最小值。这个指标在电网侧和用户侧都很敏感对电网来说峰谷差越大备用容量和调峰压力越大对商业楼宇或小区配电房来说变压器容量和基本电费通常按最大需量计费峰值压不下去一年到头都要为那几分钟的尖峰买单。电动汽车充电负荷天然具备时间灵活性所以标题里用最小化总负荷峰谷差作为目标是特别贴近工程语义的选型。不过这个目标函数有两个隐藏脾气值得先说清楚。第一它只关心最大点和最小点不关心中间曲线怎么走。两段曲线可能拥有完全相同的峰谷差但一段平滑、一段剧烈波动目标函数的值却一样。这意味着最优解往往不唯一求解器会随机给你一个可行方向中间那段可能并不好看。第二max 和 min 在数学上是不可微算子。如果直接把max(L)-min(L)丢给求解器Yalmip 虽然能做一定程度的自动改写但为了稳健和可控业内标准做法是引入辅助变量把它转成一组线性约束。这个方法不需要任何非线性求解器CPLEX 能直接啃下来。1.2 用辅助变量把 max-min 目标线性化假设我们总共有 T 个时段每时段总负荷记为 L(t)它等于基础负荷加上所有车的净充电功率。我们要最小化的是目标 max_{t} L(t) - min_{t} L(t)这个表达式没法直接进线性规划。于是引入两个标量变量 L_max 和 L_min分别表示总负荷的上包络和下包络L_max L(t), 对每个 t L_min L(t), 对每个 t因为目标是最小化L_max - L_min而约束逼着 L_max 至少大于所有 L(t)、L_min 至多小于所有 L(t)所以在最优解处L_max 会自动压到最大值附近L_min 会自动抬到最小值附近。这一步是整套建模里最关键的转换我见过不少新手在这里用minimize(max(L)-min(L))直接跑然后被求解器一堆警告绕晕最后绕回来改用辅助变量世界立刻清净了。1.3 总负荷是基础负荷和充电负荷叠加出来的别只看车端曲线还有一点特别容易忽略总负荷不是电动汽车充电曲线本身而是基础负荷曲线逐时段叠加净充电功率。基础负荷就是小区、园区在没有任何电动车时原本的用电曲线高峰通常在白天和傍晚深夜低谷。电动车充电负荷叠加进去之后可能出现原本深夜低谷被填平、但傍晚高峰更高的后果。所以建模时目标函数里的 L(t) 必须写成L(t) P_base(t) sum_i P_ch(i,t) - sum_i P_dis(i,t)其中 P_ch 是充电功率P_dis 是放电功率二者单位都是 kW。把这个叠加关系作为辅助变量计算式放进约束或后处理里优化器才会真正考虑基础负荷的时序影响。2. 约束条件搭建一辆电动汽车的物理边界2.1 功率上下限与充放电互斥binvar 的用法每一辆车的充电功率和放电功率都不是任意值它们受制于车载充电机、电池管理系统和充电桩的最大功率。另外同一辆车在同一时段不能同时充电和放电——这在物理上不成立也会毁电池。所以每辆车每个时段需要引入一个二元变量 y(i,t)y(i,t) 1 表示该时段充电 y(i,t) 0 表示该时段放电或空闲然后约束写成0 P_ch(i,t) P_max(i) * y(i,t) 0 P_dis(i,t) P_max(i) * (1 - y(i,t))这两条约束同时限制了功率上限和充放互斥。当 y1 时P_dis 被压到 0当 y0 时P_ch 被压到 0。实测下来这种写法比逐个 if 分支判断要干净得多求解器也喜欢这种结构。2.2 SOC 递推公式效率、单位、时长三者缺一不可电池 SOC荷电状态是描述车辆能量状态的唯一变量。递推公式写错一个系数结果就会整体跑偏。最常见的形式是SOC(i,t1) SOC(i,t) ( eta_ch * P_ch(i,t) - P_dis(i,t) / eta_dis ) * dt / E(i)这里每个量各司其职dt 是一个时段的时长单位小时。比如 15 分钟一个断面dt 0.251 小时一个断面dt 1。功率单位是 kW乘完 dt 得到的是 kWh。eta_ch 是充电效率通常在 0.9~0.95。含义是从电网侧取 1 kW 电真正进电池的只有 eta_ch 份。eta_dis 是放电效率同理。电池侧放出功率经过逆变器损耗后实际送到电网侧的是 P_dis所以电池内部消耗是 P_dis / eta_dis。E(i) 是电池容量单位 kWh。SOC 是无量纲比例所以最终用能量增量除以电池容量。我建议在模型里先写一个注释块把单位标清楚。等到调结果的时候你会感谢当年的自己。2.3 接入时段边界非接入时刻必须钳零很多车辆不是全天在线。每辆车有自己的到达时段 arrival 和离开时段 departure。在尚未接入、或已经离开的时段充放电功率必须严格为 0否则优化器会把充电安排在你根本没插枪的时间——我管这叫幽灵充电。对应约束很简单P_ch(i, 1:arrival(i)-1) 0 P_ch(i, departure(i):T) 0 P_dis 同理同时SOC 在非接入时段保持不变。只要功率为 0递推公式自然会维持 SOC 恒定不需要额外约束。离开时段的目标则必须写在 departure 这个截点而不是随便某个中间时刻。2.4 系统层约束变压器容量和充电桩数量如果仿真的是一个小区或园区通常还要加系统级约束。最常用的是每时段总负荷不超过变压器或进线容量P_base(t) sum_i P_ch(i,t) - sum_i P_dis(i,t) S_max这个约束其实和峰谷差目标高度相关它限制的是绝对峰值能保证求解结果不过载。另一种是充电桩数量约束。假设站里只有 K 个桩那同一时刻正在充电的车数不能超过 Ksum_i y(i,t) K注意这里 y 是充电状态指示变量如果车辆插着但不充电会有偏差不过大多数场景直接用正在充电作为占用桩的近似已经够用。3. Yalmip代码实战从车辆参数表到可运行脚本3.1 先造一份合理的测试数据建模之前先把测试场景数据准备好。我常用的是 96 个时段每个时段 15 分钟模拟一整天车辆数量从 20 到 100 台不等。每台车需要这些字段电池容量 E、最大充电/放电功率 P_max、到达时段 arrival、离开时段 departure、初始 SOC、目标 SOC。典型值如下表参数典型值说明T96一天 24 小时15 分钟一个断面N100参与优化的车辆数E50~75 kWh不同车型差异明显P_max7 kW常规交流慢充上限eta_ch / eta_dis0.92 / 0.92充放电双向效率SOC_arr0.2~0.5到达时剩余电量随机分布SOC_target0.9离开时用户期望电量arrival早高峰后至傍晚模拟通勤规律departurearrival 后 8~24 时段夜间长时间停留用随机数生成时务必让departure - arrival留出足够的充电窗口否则可能出现无解后面第 5 章会专门说。3.2 变量声明与约束装配核心代码拆开讲Yalmip 里建模我习惯分四步走。第一步声明变量第二步把约束垒成一个 cell 或直接拼接第三步写目标第四步调用 optimize。%% 变量声明 P_ch sdpvar(N, T, full); % 充电功率单位 kW P_dis sdpvar(N, T, full); % 放电功率单位 kW SOC sdpvar(N, T1, full); % SOC第 T1 列对应时段末 y binvar(N, T, full); % 充放状态1 充0 放 L_max sdpvar(1, 1); % 辅助变量总负荷上包络 L_min sdpvar(1, 1); % 辅助变量总负荷下包络第二步是约束装配。我习惯用cons []然后不断追加虽然循环装配在规模大时稍慢但可读性最好cons []; %% 初始 SOC 和功率上下限 cons [cons, SOC(:, 1) SOC_arr(:)]; cons [cons, 0 P_ch P_max .* y]; cons [cons, 0 P_dis P_max .* (1 - y)]; %% SOC 递推 for i 1:N for t 1:T cons [cons, SOC(i, t1) SOC(i, t) ... (eta_ch * P_ch(i, t) - P_dis(i, t) / eta_dis) * dt / E(i)]; end end %% 离开时目标 SOC for i 1:N cons [cons, SOC(i, departure(i)) SOC_target(i)]; end %% 非接入时段功率钳零 for i 1:N if arrival(i) 1 cons [cons, P_ch(i, 1:arrival(i)-1) 0, ... P_dis(i, 1:arrival(i)-1) 0]; end if departure(i) T cons [cons, P_ch(i, departure(i):T) 0, ... P_dis(i, departure(i):T) 0]; end end第三步把总负荷算出来然后接辅助变量L_total P_base sum(P_ch, 1) - sum(P_dis, 1); for t 1:T cons [cons, L_max L_total(t)]; cons [cons, L_min L_total(t)]; end %% SOC 上下限保护 cons [cons, 0.1 SOC 0.95];3.3 求解调用与结果提取求解核心代码其实就一行options sdpsettings(solver, cplex, verbose, 1, ... cplex.mipgap, 0.01); sol optimize(cons, L_max - L_min, options);sdpsettings里的verbose控制求解器输出等级调试时用 2正式跑用 1 或 0。cplex.mipgap是 CPLEX 的混合整数规划相对间隙阈值我后面会展开讲。求解完成后所有变量用value()提取P_ch_val value(P_ch); P_dis_val value(P_dis); SOC_val value(SOC); L_val value(L_total); fval value(L_max - L_min);出图时把基础负荷、无序充电总负荷、有序充电总负荷三条曲线叠一起一眼就能看出削峰填谷效果。无序充电曲线一般直接用车到即充、充满停的逻辑仿真得到。4. CPLEX求解规模、参数和结果校验4.1 这个MILP到底有多大先心里有个谱100 辆车、96 个时段二元变量数量是 N×T 9600 个连续变量包括 P_ch、P_dis、SOC大约 3 万个约束数量在 4 万左右。这个规模对 CPLEX 来说属于中等偏小通常几秒到几分钟能出结果不需要特别焦虑。真正影响求解时间的是二元变量和连续变量的耦合方式。上面这种功率乘以二元变量的约束是标准的大 M 结构CPLEX 处理得很成熟。如果所有车的参数都一样求解器内部还能做对称性检测速度更快。4.2 求解器参数配置我的实测推荐值我一般不会用默认参数直接跑会按实际需求调这四类mipgap默认值通常很严比如 1e-4。对工程问题来说峰谷差差 1% 根本无感把 mipgap 放到 0.01 或 0.005求解时间经常能降一个数量级。timelimit给求解器一个时间上限比如 600 秒避免陷入马拉松式搜索。nodeslimit限制分支节点数配合 mipgap 用防止内存爆炸。integrality tolerance一般保持默认除非碰到半整数解导致的异常现象。对应的sdpsettings写法是options sdpsettings(solver, cplex, ... cplex.mipgap, 0.01, ... cplex.timelimit, 600, ... cplex.nodeslimit, 200000, ... verbose, 1);另外提醒一句Yalmip 的模型本身是求解器无关的。如果没有 CPLEX 授权换成开源的 CBC、SCIP只要把 solver 字段改掉其余代码基本不动。区别只在求解速度和大规模问题的稳定性上。4.3 求解完成后必须做的三类校验求解器说最优解找到了不代表结果一定物理合理。我每跑完一个 case 都会做三层校验缺一不可。第一层数值校验。手动按 SOC 递推公式重新算一遍和 value(SOC) 逐时段比较误差应该到 1e-6 级别。如果对不上多半是效率公式写错或索引错位。第二层边界校验。检查每辆车每个时段的功率是否落在上下限内、是否有同一车同一时段 P_ch 和 P_dis 同时非零SOC 是否始终在上下限之内。通常把违反项单独打印出来定位是哪台车、哪个时段。第三层结果合理性。把有序充电的峰谷差和无序充电的峰谷差做个对比看改善比例。如果改善幅度异常大或者异常小回头查基础负荷和车辆参数是否合理。5. 我在这个项目里踩过的五个坑5.1 忽略充放电效率导致电池无限循环最早期版本我的 SOC 递推里没有效率系数默认 eta1。结果求解器给出的最优解里大量车辆在关联时段反复充、放、充、放——从电网侧看能量没有任何损耗峰谷差被这些循环电流压得很低但电池在现实里早就报废了。加上 eta_ch 和 eta_dis 之后循环不再划算因为每次往返都有能量损失求解器自然会减少无意义的充放切换。不过如果目标压得太狠解里仍可能出现个别车辆的多次切换建议在统计结果时把每辆车的充放切换次数算出来超过阈值就要注意了。5.2 随机测试数据不能随便生成可达性校验用随机函数生成到达、离开时段和目标 SOC 时很容易造出物理不可达的数据。比如一辆车到达时 SOC 只有 0.15目标是 0.9但接入窗口只有 6 个小时充电功率又只有 3.5 kW怎么算都充不满。我在代码里加了一个预处理函数生成每辆车之后先做可达性判断needed_energy (SOC_target(i) - SOC_arr(i)) * E(i); available_energy_max P_max(i) * eta_ch * (departure(i) - arrival(i)) * dt; if needed_energy available_energy_max % 丢弃或调整目标 end不满足就调低目标 SOC或者直接把这辆车从算例里拿掉。这个检查只要几行代码却能省掉后面排查无解原因的大量时间。5.3 削峰不填谷的退化解还有一次结果让我盯着曲线看了半天峰谷差确实被压到很小但代价是傍晚高峰削平了深夜低谷却没怎么填——因为目标函数只要求最大值减最小值最小不区分是削峰还是填谷做到的。两辆车、三个候选方案可能峰谷差完全一样但工程上我们通常更希望同时削峰和填谷。我的处理办法是在目标里加一个很小的二次项或斜率惩罚项。但要注意如果直接加sum((L_total - mean(L_total)).^2)模型就从 MILP 变成 MIQPCPLEX 也能解但速度会慢。更稳妥的做法是保留界面L_max - L_min lambda * sum(abs(L_total(t1) - L_total(t)))绝对值在 Yalmip 里可以用 asb 函数或增加辅助变量转成线性约束。比如引入正负偏差变量用一对约束表示相邻时段负荷差的绝对上限。这样目标仍是线性的求解器非常开心。5.4 时段粒度从24切到96求解时间非线性暴涨从 24 时段改成 96 时段变量数是 4 倍但我实测求解时间涨了不止 10 倍。原因很简单整数变量变多后分支定界树的搜索空间是超线性增长的。如果你是新手或者只是验证算法我建议分三步走先用 24 时段、20 辆车跑通逻辑确认目标函数和约束没问题然后扩到 96 时段、50 辆车观察求解时间最后再上 100 辆车。如果时间依然不可接受优先放宽 mipgap或者考虑滚动时域优化——把一天的调度拆成几个重叠的短窗口每个窗口只优化 8 到 12 个小时。滚动时域虽然牺牲了全局最优但工程上往往足够好。5.5 结果曲线中间抖动怎么压平目标只盯最大最小中间段经常会出现类似锯齿的上下抖动。比如负荷曲线在 800 kW 到 850 kW 之间反复横跳而峰谷差早就锁死在目标值附近了。这种曲线拿去汇报非常难看。我的经验是把最大允许的单位时段功率变化率设成一个约束。比如限制相邻时段总负荷变化不超过 30 kW-30 L_total(t1) - L_total(t) 30如果约束过紧可能导致无解就把这个斜率限制放宽或者在目标里加一个非常小的惩罚系数。这种处理让曲线平滑很多而且对峰谷差的影响通常很小。要提醒的是斜率约束本质上限制了电动汽车功率调节的爬坡率这在现实中确实存在只是值需要结合具体充电桩特性来设置。整个项目走下来我个人最大的体会是Yalmip 和 CPLEX 的语法都不难真正的门槛在建模之前的翻译环节——把工程语义换成数学表达再把数学表达换成一组对求解器友好的线性约束。遇到无解或怪异解时不要急着怀疑求解器先把 SOC 递推公式的单位、效率、索引对齐一遍大多数问题都出在这些不起眼的地方。