蒙特卡洛模拟在电动汽车充电负荷计算中的应用与MATLAB实现

做配电网规划和充电设施布局的同行应该都有同感电动汽车充电负荷最大的难点不是算不清楚而是它本身就带着很强的随机性。你问一辆车几点开始充电不同的人答案完全不一样你问它出发前电池还剩多少电不同车型、不同驾驶习惯差异更大。传统的确定性负荷计算在这种场景下很容易失真而蒙特卡洛模拟恰好就是处理这类随机性问题的通用工具。这篇文章把我用MATLAB做大规模电动汽车充电负荷计算的全过程拆开讲一遍从概率建模、参数设置到代码实现和结果验证最后把实际跑仿真时踩过的坑也列出来适合正在做充电负荷预测、配电网承载力分析、充电桩规划或者写相关论文的朋友参考。1. 先搞清楚要算什么充电负荷问题的数学本质1.1 目标函数与建模思路计算大规模电动汽车充电负荷本质上要回答一个问题一个区域内成千上万辆电动汽车在一天24小时里每个时刻总共会从电网取多少电。写成数学形式并不复杂。把一天按15分钟一个点切成96个时段定义第 i 辆车在第 t 个时段的充电功率为 p_i(t)那么整个区域的总充电负荷就是P(t) Σ p_i(t)i 1, 2, ..., Nt 1, 2, ..., 96这里的难点不在求和而在 p_i(t) 本身是随机变量。一辆车什么时候开始充电、充多久、功率多大取决于起始充电时间、起始SOC荷电状态、充电功率、电池容量、充电效率等一系列因素。这些因素不是固定的数值而是服从某种概率分布的随机变量。所以充电负荷计算本质上是一个“大量随机单体行为叠加”的问题适合用蒙特卡洛模拟来解决先对每辆车的各项参数进行随机抽样生成一辆“虚拟电动车”的完整充电行为得到一个单体的负荷序列然后把成千上万个单体负荷按时间对齐叠加最终得到总负荷曲线。1.2 为什么偏偏用蒙特卡洛模拟有人可能会问能不能用解析法直接算期望曲线理论上当然可以但实践上非常麻烦。每个随机变量都有自己的概率密度函数要算总负荷的期望需要把所有随机变量的分布做卷积。起始充电时间如果服从正态分布、日行驶里程如果服从对数正态分布它们组合之后的总负荷分布往往没有简洁的解析表达式。变量一多解析法基本就走不通了。蒙特卡洛的思路完全不同。它不追求推导出完美的数学公式而是靠“大量抽样 统计平均”逼近真实分布。你只需要做两件事一是知道每个随机变量满足什么分布二是会从这些分布里抽样。剩下的事情交给大数定律抽样的样本量足够大时样本均值会收敛到数学期望样本分布会逼近真实分布。这种处理方式特别适合现在做充电负荷研究。原因有三点模型扩展方便。今天只考虑私家车明天要加上公交车、出租车、网约车只需要增加对应的分布和抽样规则不用重推公式。对“边界条件”友好。比如电池SOC不能超过100%、低于某个阈值就不充电这些条件在解析法里处理起来很别扭但放在蒙特卡洛的循环里就是几个if判断。可解释性强。直接能看到每辆车的充电行为输出结果也是直观的负荷曲线评审和写报告都容易讲清楚。蒙特卡洛的代价也很明显计算量大而且结果是统计估计而非精确值。误差大致和 1/√N 成正比N是仿真车辆数。这意味着想要误差减半样本量需要变成原来的4倍。实际工程里一般取几千到几万辆车就能得到比较平滑的曲线完全在可接受范围内。2. 单体充电行为建模每个随机变量都要有依据2.1 起始充电时间决定峰值位置的第一个分布起始充电时间对负荷曲线形态的影响最大它直接决定了充电高峰出现在一天中的什么位置。大家最常用的是正态分布模型。这是因为实际出行统计数据比如美国的NHTS全国家庭出行调查显示私家车下班回家后接入充电桩的时间集中在傍晚到夜间。一个常见参数是均值19点标准差2.4小时也就是t_start ~ N(19, 2.4²)但是直接用一个标准的正态分布会有一个实际问题抽样出来的时间可能出现在白天比如上午10点或者下午2点。对于“私家车晚上回家充电”这个场景白天充电显然不合理。所以必须对分布做截断处理只在17点到24点、以及次日0点到8点这个区间内抽样。这里的8点界限也不是随手写的它对应的是“大部分车主在早上出门前拔枪”的时间窗口。如果你觉得某个地区早上6点以后几乎没人充电把上限改成6点也完全可以。这恰恰体现了这类仿真研究的特点分布形式可以有通用模板但具体参数必须结合实际场景调整。在MATLAB里做截断抽样最简单的办法是用while循环配合判断条件t_start normrnd(19, 2.4); while t_start 24 || t_start 0 || (t_start 8 t_start 17) t_start normrnd(19, 2.4); end有人觉得用while循环一直抽效率太低其实不用担心。正态分布在均值19附近概率密度最高落在17到24区间的概率超过80%落到非法区间的概率不高重采样几次就出来了。真正要注意的是别把截断条件写反否则循环可能卡死。2.2 日行驶里程与起始SOC一对绑定的随机变量车主每天行驶多少公里直接决定了电池在接入充电桩时还剩多少电。大量统计研究倾向于用对数正态分布来描述日行驶里程概率密度函数是f(x) 1 / (x·σ·√(2π)) · exp( -(ln x - μ)² / (2σ²) )经验参数一般取 μ 3.2σ 0.88此时日均行驶里程的中位数大约在 e^3.2 ≈ 24.5 km均值会高一些整个分布向右拖尾少数车主每天跑上百公里。这个拖尾特性很重要因为正是这部分“高里程用户”对充电负荷的贡献最大。有了里程起始SOC就可以估算SOC_init max(0.1, 1 - mileage / R)R是车辆满电续航里程比如300km。这个式子的逻辑是假设车主每天出发时都是满电一天消耗的电量等于行驶里程除以续航里程那么到家时剩余电量就是1减去消耗比例。加上 max(0.1, ...) 是为了防止极端情况下算出负的SOC对电池保护下限做个兜底。当然真实情况比这个线性关系复杂比如有的车主习惯电量低于30%才充有的车子还受温度影响。但对于区域级负荷计算这个线性近似的精度已经足够。很多公开发表的论文也是这么处理的关键是你要在论文里把这个假设写清楚。2.3 充电功率、电池容量和充电效率充电功率取决于充电桩类型和车载充电机能力常见分两类交流慢充功率3.5kW或7kW主要出现在居民小区、单位停车场充电时间长。直流快充功率30kW到120kW不等主要出现在公共充电站充电时间短。实际计算时要按场景配置比例。如果研究某居民区夜间充电负荷全部用7kW慢充就合理如果研究城区的公共充电站就要把快充桩功率设成30kW、60kW甚至120kW并按占比抽样。电池容量同样不是固定值。不同车型电池容量可能从40kWh到100kWh不等可以按一个离散分布来抽样也可以简化为取典型值60kWh。取典型值的好处是计算快、结果好解释代价是忽略了车型差异。如果论文需要更高的精细度建议按比例混合几种典型车型这样能看出不同车型数量比例对负荷的影响。充电效率在计算中容易被忽略但它直接影响电网侧取电量。电池侧充进去1kWh电网侧实际要供大约1.1kWh因为充电过程有AC/DC变换损耗和电池内阻损耗。工程上取 η 0.9 比较常见也就是说电网取电量 电池充电所需电量 / 0.9。2.4 一个典型车主的充电时序推导为了把上面的参数串起来我举个具体的例子。假设一辆车电池容量 E_bat 60kWh续航里程 R 300km慢充功率 P 7kW充电效率 η 0.9。某次抽样的日行驶里程是 m 40km那么起始SOC大约是SOC_init 1 - 40/300 ≈ 0.867这里就出现了一个边界情况如果目标SOC上限设为0.9那么这辆车需要的充电量是(0.9 - 0.867) × 60 2kWh折算到电网侧是2/0.9 ≈ 2.22kWh。按7kW功率计算充电时间只有0.32小时大约19分钟。这个例子说明并非所有车接入充电桩都会长时间充电。日行驶里程比较短的车主可能插上枪很快就充满就拔了。真正对负荷曲线产生持续贡献的是那些日行驶里程在几十公里以上、回家时SOC比较低的车。再看一个更典型的情况某车主日行驶里程60kmSOC_init 1 - 60/300 0.8仍低于0.9所以需要充电。所需充电量为0.2 × 60 12kWh电网侧取电12/0.9 ≈ 13.33kWh充电时间约为13.33/7 ≈ 1.9小时也就是大约2小时。这个时长就非常接近实际生活中“晚上下班充两三个小时”的体验了。映射到15分钟一个时段大约占8个时段。3. MATLAB代码从框架到落地能直接跑的版本3.1 程序主框架与核心参数初始化我先把一个可以直接运行的版本放出来后面逐行解释。这个版本针对的是“私家车夜间慢充”场景车辆数设1万辆一天分成96个时段每个时段15分钟。为了让论文结果可复现程序开头设置了随机数种子。% 蒙特卡洛电动汽车充电负荷计算 - 私家车夜间慢充场景 clear; clc; rng(2024); % 参数配置 N 10000; % 仿真车辆数 T 96; % 一天分成96个15分钟时段 dt 15 / 60; % 每个时段长度小时 P_charge 7; % 交流慢充功率kW E_bat 60; % 电池容量kWh SOC_end 0.9; % 目标充电SOC上限 eta 0.9; % 充电效率 R 300; % 满电续航里程km % 预分配负荷矩阵 Load zeros(N, T); % 蒙特卡洛主循环 for i 1:N % 1. 抽样日行驶里程服从对数正态分布 mileage lognrnd(3.2, 0.88); mileage max(mileage, 0.1); % 防止出现极端小值 % 2. 估算起始SOC SOC_init max(0.1, 1 - mileage / R); % 3. 判断当天是否需要充电 if SOC_init SOC_end continue; end % 4. 抽样起始充电时间限制在17-24点和0-8点 t_start normrnd(19, 2.4); while t_start 24 || t_start 0 || (t_start 8 t_start 17) t_start normrnd(19, 2.4); end % 5. 计算充电时长并折算成时段数 t_charge_h (SOC_end - SOC_init) * E_bat / (P_charge * eta); n_charge ceil(t_charge_h / dt); % 6. 将充电事件映射到96点负荷序列 idx max(1, ceil(t_start / dt)); for k 1:n_charge Load(i, idx) P_charge; idx idx 1; if idx T idx 1; % 跨日循环 end end end % 汇总统计 TotalLoad sum(Load, 1); % 各时段总充电负荷单位kW hours (0.5:T) * dt; % 每个时段对应的中心时刻 plot(hours, TotalLoad, LineWidth, 1.5); xlabel(时刻 (h)); ylabel(总充电负荷 (kW)); grid on;这个程序跑完后会得到一条24小时的充电负荷曲线纵轴是总负荷横轴是时刻。3.2 单体抽样与负荷序列生成逐行解释下面把主循环里的关键步骤拎出来说清楚。第一步抽样日行驶里程。lognrnd(3.2, 0.88) 是MATLAB内置的对数正态分布抽样函数。抽出来的数据分布特征非常符合实际大多数车在10到40km之间少数车超过100km。加一行 mileage max(mileage, 0.1) 是为了防止抽样到极其接近0的数值那会让SOC_init约等于1计算上没有问题但是偶尔会出现0后面算SOC时不好看。第二步估算SOC_init。直接用线性公式同时用max(0.1, ...)做下限保护。这一步顺便给出了“是否需要充电”的判断依据。如果起始SOC已经达到0.9就认为当天不需要充电continue跳过。这个逻辑看似简单但对最终结果影响很大。如果漏掉这个判断、让所有车都充电负荷峰值会虚高一截论文数据就不真实了。第三步抽样起始充电时间。while循环的作用前面说过是为了把抽样结果限制在“傍晚到深夜、凌晨到清晨”这个合理窗口内。这里有一点值得注意t_start 24的情况ceil(24/dt) 96正好映射到当天最后一个时段。而 t_start 0.1 时idx ceil(0.1/0.25) 1对应0点到0点15分。整个时间映射是连续的。第四步计算充电时长。这里用的是t_charge_h (SOC_end - SOC_init) * E_bat / (P_charge * eta)分子是“从起始SOC充到目标SOC需要充入电池的能量”单位kWh分母是“电网侧实际消耗功率”因为充电效率小于1所以电网侧功率要按P_charge × eta来算。再用ceil向上取整得到需要占用的时段数。注意我用了ceil而不是round或floor因为充电时间只要超过一个时段的起始点就该把这个时段算进去否则会低估负荷。第五步把充电事件写入负荷矩阵。这一步是程序里最容易出错的地方。核心逻辑是从起始时段idx开始连续给n_charge个时段赋充电功率。如果时间跨过24点就把索引绕回到1模拟“次日凌晨继续充电”的行为。有人问为什么要用Load(i, idx)这样先存单车负荷再求和而不直接累加到TotalLoad上主要好处是方便后续做统计比如每辆车的平均充电时长、总充电电量写论文时这些都是很有价值的数据。3.3 从单群到多群私家车、公交、出租混合建模上面的程序只考虑了私家车夜间慢充现实中一个区域往往同时存在私家车、公交车、出租车网约车等多种车型。它们的充电行为差异很大私家车慢充为主接入时间集中在傍晚和夜间功率7kW左右。公交车白天运行间隙快充晚上收班后在场站集中充电功率60kW到120kW。出租车行驶里程大、充电频次高白天随机快充为主功率30kW到60kW。混合建模的思路很简单按车辆类型占比拆分总车辆数分别跑一次蒙特卡洛最后把各类型负荷曲线叠加。% 车辆类型配置示例 car_ratio [0.7, 0.15, 0.15]; % 私家车、公交车、出租车占比 P_level [7, 90, 60]; % 各类车典型充电功率 % 对每一类车辆分别调用子函数生成对应负荷曲线具体实现时建议把“生成单车充电事件”封装成一个子函数输入是车辆类型输出是一个长度为T的负荷序列。这样主程序就只负责按比例循环调用子函数和累加结果结构清晰后面对某个类型的参数做敏感性分析时也方便。3.4 加速技巧向量化、parfor和随机数流管理当车辆数从1万增加到10万甚至更多时逐车循环会变得很慢。我实际测试过96点时段的程序1万辆车在普通PC上大概需要几秒到几十秒但10万辆就可能需要好几分钟。如果每个时段粒度升级到1440个点按分钟运行时间会进一步增长。有几个实用加速手段。第一尽量用向量化。如果想避免内层k循环可以先记录每辆车的起始时段和持续时段数然后生成一个“充电事件表”最后一次性累加到负荷矩阵上。这个方法代码稍微复杂但速度提升明显。第二用parfor并行。把for i 1:N改成parfor i 1:N前提是需要Parallel Computing Toolbox。但这里有个隐藏坑parfor的每个worker如果共享同一个全局随机数流不同worker可能生成极其相似的结果导致仿真可信度下降。最简单的处理方式是提前为每个worker设置不同的随机种子% 在parfor循环内部开头加一行 s RandStream(mlfg6331_64, Seed, 2024 10 * i); RandStream.setGlobalStream(s);这样每个循环独立种子的做法虽然牺牲了一点效率但保证了并行结果可复现。第三如果只是做参数敏感性分析不需要每次都用新的随机样本可以先一次性生成所有随机数矩阵再分块喂给主循环也能省去大量重复初始化开销。4. 仿真结果长什么样典型曲线与场景对比4.1 无序充电的典型日负荷曲线按照上面程序跑出来1万辆纯私家车慢充场景的结果大致是这样的凌晨0点到6点负荷逐渐下降到接近零早晨7点前后几乎没有充电行为白天有小幅波动到了傍晚18点开始快速拉升19点到20点左右出现峰值然后缓慢回落凌晨1点后逐步归零。峰值出现在19到20点这非常符合直觉。因为大多数车主下班到家后接入充电桩7kW的慢充功率会持续好几个小时但大家接入时间高度集中在18到21点之间于是形成了明显的“晚高峰叠加”。这个结果对配电网来说不是什么好消息现在的电网本来就有一个晚间生活用电高峰大约在19点到22点电动汽车无序充电的负荷高峰正好和它重叠。如果电动汽车渗透率提升配电网的晚间峰值压力会非常明显。这也是为什么很多研究都盯着“有序充电”和“错峰充电”策略的原因。4.2 不同渗透率下峰值负荷的变化规律把上面的程序稍微改一下加上区域车辆总数和渗透率两个变量就能做场景对比。我列一个示意结果具体数值取决于参数设定但规律是一致的。电动汽车渗透率EV数量充电负荷峰值峰值出现时刻相比基荷的增幅0%08000 kW20:00-10%30008800 kW20:0010%30%900011350 kW19:3042%50%1500014200 kW19:3078%最值得注意的不是数值本身而是增幅并不是线性的。原因在于相当一部分车当天根本不需要充电或者充电时长很短它们对峰值负荷的贡献远小于“车辆数 × 充电功率”的上限估计。这也提醒我们做配电网规划时不能用“所有EV同时以最大功率充电”这种最保守假设否则会造成巨大的设备浪费但也不能因此就对影响过于乐观因为随着渗透率升高晚高峰负荷的重叠效应会越来越明显。4.3 结果到底准不准三种验证思路蒙特卡洛仿真的结果不能只看曲线“顺眼”还要做验证。说三种我实际用过的思路。第一和实测数据对比。如果你手头有某个充电站的真实负荷记录可以把同场景参数代入仿真然后对比仿真负荷曲线和实测曲线的形状与峰值。两者峰值通常不会完全相等因为很多车主行为参数都是估计的但趋势和峰值时段应当高度一致。第二和解析解对比。在只保留一个随机变量、其他参数全部固定的简化条件下总负荷的期望可以用卷积直接算出来。把蒙特卡洛结果和解析结果放在一起对比如果偏差在几个百分点以内至少说明主循环逻辑没有写错。第三重复实验看置信区间。保持参数不变换不同的随机种子跑20次统计每个时段负荷的均值和标准差。如果峰值负荷的变异系数标准差/均值小于5%说明样本量足够结果稳定如果变异系数很大就需要增大N。5. 跑仿真时踩过的坑常见问题与排查实录5.1 现象、原因与解决办法速查表这个表里整理的问题基本都是我在实际运行这类仿真时真遇到过的包括自己写错的和帮别人改代码时发现的问题。现象可能原因解决办法充电峰值出现在白天中午起始充电时间分布没有截断抽样到了白天用while循环限制在17-24点与0-8点负荷曲线毛刺多不光滑车辆数太小或日行驶里程的方差过大产生离群值增大N到1万以上对里程设置合理上限峰值负荷高得离谱没有判断SOC_init SOC_end所有车都充电在抽样后加continue条件每次运行结果都不一样没有设置随机数种子程序开头加 rng(2024)凌晨点出现尖峰负荷跨日循环的索引处理有误idx绕回1时多写了时段检查idx T时的处理和for循环计数出现NaN或Inflognrnd参数不对或里程出现0对里程做max(mileage, 0.1)下限保护parfor结果和for结果差异很大不同worker的随机数流未同步每个迭代内部单独设置RandStream种子程序跑得很慢内层循环过多矩阵未预分配预分配Load矩阵记录事件表后一次性累加这里再展开说两个最容易“中招”的细节。第一个是continue的位置。很多人会把“是否需要充电”的判断放在抽样时间之后结果没问题但是白白浪费了抽样时间还有人把判断写在SOC_init计算结果之前直接用预设值做判断等于所有车都充电或者都不充电曲线完全失真。正确做法是先算SOC_init再判断再抽起始时间。第二个是时段索引的边界。例如t_start 24ceil(24/0.25) 96对应23:45到24:00这个时段这没问题。但如果是t_start 23.9ceil(23.9/0.25) ceil(95.6) 96也一样落在最后一个时段。如果这里改用round就可能把23.9映射到第95时段等于把充电开始时间提前了15分钟虽然单辆车影响很小但一万辆车的开始时间整体偏移峰值位置就会失真。5.2 提升可复现性的几个小习惯做仿真研究最怕的是自己过两个礼拜回来看代码已经搞不清楚当时用的什么参数。我的习惯是把所有参数打包成一个结构体para每次运行前自动保存一份带时间戳的结果文件。para.N N; para.T T; para.dt dt; para.P_charge P_charge; para.E_bat E_bat; para.SOC_end SOC_end; para.eta eta; para.R R; % ... 其他参数 save([sim_, datestr(now, yyyymmdd_HHMMSS), .mat], TotalLoad, para);同理跑参数敏感性分析时每次只改动一个参数把其他参数固定结果对比才有意义。改参数的时候不要直接改代码里的数字而是把参数放进for循环里批量跑比如分别计算渗透率10%、30%、50%的场景最后画在一张图上对比。还有一个实用小技巧如果终端输出需要进度条就在主循环里加waitbar(i/N, Simulating...)。虽然会拖慢一点运行速度但当你跑10万辆车的时候能知道大概还要等多久心里踏实很多。跑完仿真后建议把总负荷的均值、峰值、峰值时刻、总充电电量都算出来做成一张结果汇总表。这些数字就是论文结论或者项目报告的核心论据不要等到最后再回去翻Load矩阵那样非常容易出错。最后再分享一点个人体会。用蒙特卡洛做充电负荷计算我走过一段弯路一开始沉迷于调“更好看”的分布函数和更复杂的模型结果发现对最终结论影响最大的反而是参数取值的合理性。你把起始充电时间的均值从19点改成18点晚高峰峰值就可能提前一小时出现你把充电效率从0.9改成0.85总充电电量就会增加大约6%。这些参数背后都应该有实际调研数据支撑而不是拍脑袋定的。所以做这类研究我建议先把精力花在收集本地充电行为数据上哪怕只是从几篇公开论文或行业报告里整理出参数范围也比直接抄一套“万能参数”要靠谱得多。代码只是把模型跑起来的工具真正的功夫全在建模时的选择和取舍上。