简介这套MATLAB程序包与《交通流建模》一书配套面向交通工程、智能交通领域的初学者和研究者旨在借助运行代码消除交通流理论与仿真实践之间的鸿沟适合课程设计、毕业设计阶段动手验证三参数关系与不同模型特点。压缩包内含16个m文件整体仅12KB属于轻量级代码集按功能划分出基本图绘制、宏观LWR与微观跟驰仿真、元胞自动机规则、参数初始化、坐标后处理及扰动加载等模块可独立运行或配合教材逐段调试。目前已有282人在线学习浏览。通过实际操作学习者既能掌握流量-密度-速度三参数关系的可视化方法也能对比宏观模型、元胞自动机模型和微观跟驰模型的适用场景同时培养数据预处理、仿真结果分析与信号配时等控制策略验证能力可作为后续研究或课程设计的起点。1. 交通流建模 MATLAB 程序包宏观守恒方程与微观跟驰模型的同台复现赶过晚高峰的人几乎都见过导航地图上的深红色路段明明没有事故车流却走走停停最后又莫名恢复。这类现象正是交通流建模要解释的核心。Traffic Flow Modelling 是交通工程里把宏观连续模型、微观跟驰模型讲得比较系统的经典书这套 MATLAB 程序就是它的配套代码。压缩包内有 17 个 .m 文件ESM 是场景入口CASEPARS、MODELPARS、NUMPARS、PLOTPARS 四组参数文件分别描述路网、模型、数值格式和绘图配置TrafficFlowSimulation.m 是主循环TimeStepMacro.m 和 TimeStepMicro.m 分别负责宏观与微观推进。也就是说拿到代码后不用从零搭仿真器只要读参数、理顺主循环就能复现 LWR 密度波也能看到微观跟驰模型里的走走停停。适合正在啃交通流理论的学生、要做仿真的智能交通从业者以及习惯用 MATLAB 复现论文的工程师。2. 从 CASEPARS 到 TrafficFlowSimulation参数文件与仿真主循环的装配方式2.1 四组参数文件各自管什么这套程序的第一道门槛不是公式而是参数文件。ESM 没有 .m 后缀第一次读代码时容易被忽略。我的习惯是先看 ESM再看四组 PARS 文件最后打开 TrafficFlowSimulation.m。参数文件管什么典型字段举例CASEPARS.m场景与边界条件roadLength, numLanes, initialDensity, boundaryType, bottleneckPositionMODELPARS.m车辆行为模型modelType, freeSpeed, rhoJam, s0, tau, maxAccelNUMPARS.m离散化与数值格式numCells, dx, dt, simTime, sampleEveryPLOTPARS.m输出与可视化plotMode, colorRange, movingRefVCASEPARS 管“在哪里仿”MODELPARS 管“车怎么跑”NUMPARS 管“算得稳不稳”PLOTPARS 管“看得出看不出”。把场景和模型参数分开后同一个路网可以直接切换宏观、微观模型不需要改道路长度和边界条件。NUMPARS.dt 是宏观程序里最敏感的参数如果freeSpeed * dt / dx超过 1守恒格式就会失稳微观程序则不受这个约束主要看每辆车每个时间步的位移是否超过前车间距。2.2 ESM 如何拼装参数并触发初始化ESM 的等效脚本片段如下。它没有像函数一样被call而是按脚本方式把四组参数放进当前工作区再调用初始化函数。% ESM 场景装载脚本伪代码按该文件实际工作方式整理 CASEPARS; % 载入场景参数 MODELPARS; % 载入模型参数 NUMPARS; % 载入数值参数 PLOTPARS; % 载入绘图参数 state Initialization(CASEPARS, MODELPARS, NUMPARS); [state, logs] TrafficFlowSimulation(state, CASEPARS, MODELPARS, NUMPARS, PLOTPARS);这段代码展示的是一条完整数据流四组参数进入工作区Initialization 根据宏观还是微观创建初始状态TrafficFlowSimulation 返回更新后的 state 和日志 logs。实际函数签名可能和我写的不完全一致但数据走向是这条线。四个 PARS 文件共用一套结构体命名省去了传递大量单独变量的麻烦代价是读代码时要随时注意字段来源比如CASE.numVehicle和NUM.numCells都可能存在前者用于微观初始化后者用于宏观网格。Initialization 里最常做的事情是构造一个非均匀初始密度。下面这段等价逻辑可以用来制造初始排队段% Initialization.m 中可能出现的密度初始化 x0 linspace(0, CASE.roadLength, NUM.numCells); rho0 CASE.initialDensity * ones(NUM.numCells, 1); % 在道路 30%~50% 区间把密度翻倍模拟初始拥堵 idx round(NUM.numCells * 0.3) : round(NUM.numCells * 0.5); rho0(idx) rho0(idx) * 2; state.x x0; state.rho rho0; state.v MODEL.freeSpeed * ones(NUM.numCells, 1);逻辑说明把道路离散成NUM.numCells个格子再在中间段设置一个高密度块。这样宏观程序一开始就会向右释放消散波、向左释放冲击波不需要额外做扰动。微观初始化更简单一般是等间距放置车辆再给同一个初始速度。参数说明initialDensity的单位要和roadLength、numCells对齐如果写成辆/km 而道路长度用米算出来的车辆数会差三个量级这是最常见的初始化错误。2.3 TrafficFlowSimulation.m 的主循环骨架宏观和微观的推进函数差别很大但主循环骨架是统一的初始化、步进、按采样周期记录日志。% TrafficFlowSimulation.m 主循环骨架 for it 1 : NUM.simTime / NUM.dt if strcmp(MODEL.type, macro) state TimeStepMacro(state, CASEPARS, MODELPARS, NUMPARS); else state TimeStepMicro(state, CASEPARS, MODELPARS, NUMPARS); end if mod(it, PLOTPARS.sampleEvery) 0 logs.t(end 1) it * NUM.dt; logs.rho(end 1, :) state.rho; logs.v(end 1, :) state.v; end end逻辑说明TimeStepMacro和TimeStepMicro都接收 state 和全部参数返回更新后的 state。PLOTPARS.sampleEvery决定日志密度设为 1 表示每个时间步都记录适合短时仿真对微观程序来说由于时间步长通常远小于宏观建议设 5 或 10否则内存里会堆满大矩阵。这里有一个常见边界问题NUM.simTime / NUM.dt不一定是整数直接用for it 1 : N会把尾步丢掉。我一般会改成while t NUM.simTime每步更新t it * NUM.dt保证最后一段时刻也被记录。3. FundamentalDiagram.m 与 TimeStepMacro.mLWR 守恒方程的数值落地3.1 三参数关系与 FundamentalDiagram.m交通流三参数关系是判断程序是否跑对的第一道标准。宏观模型里最基本的满足关系是q ρ * v也就是流量等于密度乘以速度。常见线性速度-密度关系为v vf * (1 - ρ / ρjam)代入后得到抛物线的流量-密度关系。实际高速路基本图更像三角形因为自由流段和拥堵段的斜率不同但教科书配套程序往往先用线性关系讲原理。% FundamentalDiagram.m 中 Q-K 关系绘制核心部分 rho linspace(0, 200, 200); % 密度轴辆/km vf 80; % 自由流速度km/h rhoJam 200; % 阻塞密度辆/km v vf * (1 - rho / rhoJam); % 线性速度-密度关系 q rho .* v; % 流量 密度 * 速度 plot(rho, q, LineWidth, 1.5); xlabel(密度 K (辆/km)); ylabel(流量 Q (辆/h)); grid on;逻辑说明当密度为零时流量为零当密度等于阻塞密度时速度为零、流量也为零最大流量出现在rho rhoJam / 2处对应临界密度。程序里可以把这个理论曲线保存下来再把仿真输出的密度-流量散点叠上去基本图上出现明显偏离就说明初始场、边界或数值格式有问题这是最直接的自检方式。三参数关系表达式关键含义速度-密度 V-Kv v(ρ)密度升高速度下降流量-密度 Q-Kq ρ v(ρ)临界密度处流量最大流量-速度 Q-Vq ρ(v) v高速低流、低速高流两个分支3.2 TimeStepMacro.m 的迎风差分实现LWR 守恒方程可以写作∂ρ/∂t ∂q/∂x 0。离散时不能直接用中心差分因为交通波有明确方向信息从上游传向下游所以要用迎风格式取上风侧流量。% TimeStepMacro.m 中的一阶迎风更新片段 function state TimeStepMacro(state, CASE, MODEL, NUM) rho state.rho; v Speed(rho, MODEL); % 平衡速度 q rho .* v; % 当前流量场 % 迎风取上游格点流量作为当前格点边界流量 q_flux [q(1); q(1:end-1)]; dqdx (q - q_flux) / NUM.dx; rho_new rho - NUM.dt * dqdx; % 物理约束密度不能越界 rho_new max(0, min(CASE.rhoJam, rho_new)); state.rho rho_new; state.v Speed(rho_new, MODEL); end逻辑说明q_flux把每个格点的上游邻居当作入流流量这是典型的一阶迎风。rho_new rho - dt * dqdx就是守恒方程显式时间推进。密度截断保证数值解始终落在物理区间内。这个格式能复现激波但会有数值耗散激波前沿会被抹平一点要减小耗散可以把q_flux换成 Godunov 精确黎曼解或者加 Minmod 斜率限制本质都是更准确地计算网格边界处的流量。参数说明稳定性要求满足 CFL 条件也就是max(v) * dt / dx 1。这里max(v)直接取MODEL.freeSpeed即可不需要在运行时算最大值。如果看到密度出现负值或激波前方振荡先看NUM.dt是否偏大再看边界条件是否写对。闭环道路必须做周期边界也就是最后一个格点的上游邻居要回到最后一个或第一个格点不能简单补零开环道路则要给定入口流量否则车辆在边界处会凭空消失。3.3 Speed.m 与 Position2Spacing.m 的角色差别Speed.m在宏观和微观两个分支里的含义不同。宏观分支输入密度输出平衡速度微观分支输入车头间距输出期望速度。Position2Spacing.m则是把位置向量转成车头间距数组为微观模型准备数据。下面这种写法把两类模型统一到一个接口里% Speed.m 的等价实现 function v Speed(rhoOrSpacing, MODEL) if strcmp(MODEL.type, macro) v MODEL.freeSpeed * max(0, 1 - rhoOrSpacing / MODEL.rhoJam); else v MODEL.freeSpeed * (1 - exp(-rhoOrSpacing / MODEL.s0)); end end逻辑说明宏观分支里max(0, ...)保证密度超过阻塞密度后速度也不会变负微观分支用指数松弛表示间距越大越接近自由流速度。两个分支都依赖MODELPARS里的参数所以同一个Speed.m可以被宏观和微观两个推进函数调用。Position2Spacing.m则不同如果输入是宏观格点位置输出平均间距如果输入是微观车辆位置输出相邻车辆的车头间距。宏观里平均间距约等于1/ρ微观里需要减掉车长才是有效间距这个差别很容易在跨模型对比时踩坑。4. TimeStepMicro 与 DoDisturbance微观跟驰模型如何复现走走停停4.1 微观状态量更新从间距到速度再到位置微观模型的状态量很简单每辆车有位置x、速度v车辆之间通过车头间距互相影响。TimeStepMicro.m每次推进分两步先算期望速度或加速度再更新速度和位置。下面是一段等效的松弛跟驰模型实现。% TimeStepMicro.m 的车辆循环等效代码 function state TimeStepMicro(state, CASE, MODEL, NUM) x state.x; % 单车位置向量 v state.v; % 单车速度向量 nv length(x); % 车辆数 spacing Position2Spacing(x, nv); % 车头间距 v_new zeros(nv, 1); for i 1 : nv s spacing(i); vDes MODEL.freeSpeed * (1 - exp(-s / MODEL.s0)); % 期望速度 a (vDes - v(i)) / MODEL.tau; % 松弛加速度 v_new(i) max(0, v(i) a * NUM.dt); % 速度更新 end state.v v_new; state.x x v_new * NUM.dt; % 位置更新 state.spacing Position2Spacing(state.x, nv); end逻辑说明这里用的是松弛跟驰模型核心是期望速度公式vDes vf * (1 - exp(-s / s0))。当间距很小时期望速度接近零车辆减速当间距足够大时期望速度接近自由流速度。MODEL.tau是松弛时间决定驾驶员把速度调整到期望速度的快慢。注意位置更新用的是更新后的v_new相当于一种半隐式处理比用旧速度更不易出现车辆重叠。参数说明s0通常取 5~20 mtau取 1~3 s。tau太小时车辆起步过猛tau太大时车流重新加速很慢正好对应走停状态里“消散慢”的过程。微观程序的NUM.dt通常取 0.01~0.1 s比宏观小一个量级所以总步数会多很多。4.2 DoDisturbance.m在车队中部注入一次速度扰动微观仿真的难点不是让车流跑起来而是制造出真实的扰动。DoDisturbance.m在某个时间窗口内压低某辆车速度相当于人为制造一次刹车灯波。% DoDisturbance.m 的等效实现 function state DoDisturbance(state, it, CASE, MODEL, NUM) t it * NUM.dt; if t CASE.disturbTimeStart t CASE.disturbTimeEnd idx round(length(state.x) * CASE.disturbPositionRatio); % 把第 idx 辆车限制到低速模拟前车短暂减速 state.v(idx) min(state.v(idx), CASE.disturbSpeed); end end在主循环里调用时要和TimeStepMicro配合if strcmp(MODEL.type, micro) state DoDisturbance(state, it, CASE, MODEL, NUM); state TimeStepMicro(state, CASE, MODEL, NUM); end逻辑说明扰动位置disturbPositionRatio通常取 0.5~0.7不要放在头车。放在头车只会让单车减速放在车队中部会让后车因为间距变小而逐辆减速同时前车逐渐拉开形成局部密度升高并向上游传播。若扰动窗口持续 3~5 s模拟的是刹车灯波若持续到仿真结束模拟的是持续瓶颈。注意DoDisturbance之后要立刻执行一次TimeStepMicro不要连续调用两次扰动否则同一时刻速度会被压两次。4.3 固定坐标与移动坐标两种看拥堵演化的方式微观输出是每个时刻的所有车辆位置直接看轨迹图容易眼花。这套程序提供了两个观测视角固定坐标和移动坐标。PlotFixedCoordinates.m用于固定路面位置观察密度或速度随时间变化PlotMovingCoordinates.m把坐标系以参考速度v_ref平移方便跟踪密度波传播。% PreprocessPlotMovingCoordinates.m 的核心变换 function xm PreprocessPlotMovingCoordinates(x, t, v_ref) xm x - v_ref * t; % 伽利略变换到移动坐标系 end然后可以用散点图绘制移动坐标下的轨迹% PlotMovingCoordinates.m 绘制移动坐标下的轨迹点 scatter(xm, t, 6, v, filled); % 颜色表示速度 xlabel(移动坐标 x - v_{ref} t (m)); ylabel(时间 t (s)); colorbar;逻辑说明v_ref选在自由流速度附近时自由行驶的车辆轨迹接近竖直走停状态下的减速波会变成明显斜线因为波速通常低于自由流速度。如果所有轨迹都向同一个方向倾斜说明v_ref选得不合适可以取扰动窗口内波速的平均值重新变换。固定坐标绘图则更简单先把车辆位置插值到规则网格上再用imagesc画时空图PreprocessPlotFixedCoordinates.m做的主要就是这种插值网格间距太大会抹掉短时拥堵太小又会放大插值噪声。5. ESM 场景切换与守恒检验把程序改造成参数实验平台5.1 用 ESM 做批量场景扫描ESM 最大的价值是场景可复用。把 ESM 放在一个批量脚本里循环执行就能对反应时间、自由流速度、扰动窗口等参数做扫描。注意 MATLAB 的脚本循环里变量会残留跑下一组场景前要重新加载四组 PARS。% runExperiments.m 简单批量扫描 reactTimes [0.5, 1.0, 2.0]; for i 1 : length(reactTimes) ESM; % 载入默认场景 MODELPARS.tau reactTimes(i); % 覆盖反应时间 state Initialization(CASEPARS, MODELPARS, NUMPARS); [state, logs] TrafficFlowSimulation(state, CASEPARS, MODELPARS, NUMPARS, PLOTPARS); results(i).tau reactTimes(i); results(i).meanFlow mean(logs.q(end - 200 : end)); end这里ESM每次都会重新执行变量名相同所以会覆盖上一次的MODELPARS这是脚本式参数装配的副作用。如果要跑更严格的实验建议把 ESM 改成函数[CASE, MODEL, NUM, PLOT] SetupScenario(scenarioID)返回参数结构体而不是往工作区塞变量。5.2 用守恒律验证数值解闭环道路下车辆数应该守恒。宏观程序里车辆总数等于密度对空间的积分也就是sum(rho) * dx。每跑完一步都检查这个值如果漂移超过千分之一说明边界或差分格式有问题。% 闭环边界下车辆数守恒检查 totalVehicle sum(state.rho) * NUM.dx; if abs(totalVehicle - totalVehicle0) totalVehicle0 * 1e-3 warning(车辆数不守恒: %.4f - %.4f, totalVehicle0, totalVehicle); end微观模型没有这个守恒量问题因为车辆是离散对象只要不超越前车总数就不变但微观程序要检查最小间距是否小于车长否则说明发生了碰撞。5.3 用“扰动窗口提前 20 秒”验证冲击波速度这套程序里最有意思的验证方式是只改CASE.disturbTimeStart其他参数不动。把扰动窗口提前 20 s再分别用固定坐标和移动坐标画两张时空图对比两次仿真的密度波前沿位置。你会看到拥堵排队头的反向传播速度基本不变这个速度在 LWR 模型里对应特征速度在排队论里对应冲击波速度。若用移动坐标系把v_ref设成这个反向传播速度的绝对值密度波前沿会变成一条接近竖直的线说明坐标变换和波速测量都对上了。这一步跑通之后这套交通流 MATLAB 程序就不再是“别人的代码”而是你可以自由改参数、做实验、出图写结论的仿真平台。本文还有配套的精品资源点击获取 交通流建模MATLAB解析：LWR宏观守恒与微观跟驰模型
简介这套MATLAB程序包与《交通流建模》一书配套面向交通工程、智能交通领域的初学者和研究者旨在借助运行代码消除交通流理论与仿真实践之间的鸿沟适合课程设计、毕业设计阶段动手验证三参数关系与不同模型特点。压缩包内含16个m文件整体仅12KB属于轻量级代码集按功能划分出基本图绘制、宏观LWR与微观跟驰仿真、元胞自动机规则、参数初始化、坐标后处理及扰动加载等模块可独立运行或配合教材逐段调试。目前已有282人在线学习浏览。通过实际操作学习者既能掌握流量-密度-速度三参数关系的可视化方法也能对比宏观模型、元胞自动机模型和微观跟驰模型的适用场景同时培养数据预处理、仿真结果分析与信号配时等控制策略验证能力可作为后续研究或课程设计的起点。1. 交通流建模 MATLAB 程序包宏观守恒方程与微观跟驰模型的同台复现赶过晚高峰的人几乎都见过导航地图上的深红色路段明明没有事故车流却走走停停最后又莫名恢复。这类现象正是交通流建模要解释的核心。Traffic Flow Modelling 是交通工程里把宏观连续模型、微观跟驰模型讲得比较系统的经典书这套 MATLAB 程序就是它的配套代码。压缩包内有 17 个 .m 文件ESM 是场景入口CASEPARS、MODELPARS、NUMPARS、PLOTPARS 四组参数文件分别描述路网、模型、数值格式和绘图配置TrafficFlowSimulation.m 是主循环TimeStepMacro.m 和 TimeStepMicro.m 分别负责宏观与微观推进。也就是说拿到代码后不用从零搭仿真器只要读参数、理顺主循环就能复现 LWR 密度波也能看到微观跟驰模型里的走走停停。适合正在啃交通流理论的学生、要做仿真的智能交通从业者以及习惯用 MATLAB 复现论文的工程师。2. 从 CASEPARS 到 TrafficFlowSimulation参数文件与仿真主循环的装配方式2.1 四组参数文件各自管什么这套程序的第一道门槛不是公式而是参数文件。ESM 没有 .m 后缀第一次读代码时容易被忽略。我的习惯是先看 ESM再看四组 PARS 文件最后打开 TrafficFlowSimulation.m。参数文件管什么典型字段举例CASEPARS.m场景与边界条件roadLength, numLanes, initialDensity, boundaryType, bottleneckPositionMODELPARS.m车辆行为模型modelType, freeSpeed, rhoJam, s0, tau, maxAccelNUMPARS.m离散化与数值格式numCells, dx, dt, simTime, sampleEveryPLOTPARS.m输出与可视化plotMode, colorRange, movingRefVCASEPARS 管“在哪里仿”MODELPARS 管“车怎么跑”NUMPARS 管“算得稳不稳”PLOTPARS 管“看得出看不出”。把场景和模型参数分开后同一个路网可以直接切换宏观、微观模型不需要改道路长度和边界条件。NUMPARS.dt 是宏观程序里最敏感的参数如果freeSpeed * dt / dx超过 1守恒格式就会失稳微观程序则不受这个约束主要看每辆车每个时间步的位移是否超过前车间距。2.2 ESM 如何拼装参数并触发初始化ESM 的等效脚本片段如下。它没有像函数一样被call而是按脚本方式把四组参数放进当前工作区再调用初始化函数。% ESM 场景装载脚本伪代码按该文件实际工作方式整理 CASEPARS; % 载入场景参数 MODELPARS; % 载入模型参数 NUMPARS; % 载入数值参数 PLOTPARS; % 载入绘图参数 state Initialization(CASEPARS, MODELPARS, NUMPARS); [state, logs] TrafficFlowSimulation(state, CASEPARS, MODELPARS, NUMPARS, PLOTPARS);这段代码展示的是一条完整数据流四组参数进入工作区Initialization 根据宏观还是微观创建初始状态TrafficFlowSimulation 返回更新后的 state 和日志 logs。实际函数签名可能和我写的不完全一致但数据走向是这条线。四个 PARS 文件共用一套结构体命名省去了传递大量单独变量的麻烦代价是读代码时要随时注意字段来源比如CASE.numVehicle和NUM.numCells都可能存在前者用于微观初始化后者用于宏观网格。Initialization 里最常做的事情是构造一个非均匀初始密度。下面这段等价逻辑可以用来制造初始排队段% Initialization.m 中可能出现的密度初始化 x0 linspace(0, CASE.roadLength, NUM.numCells); rho0 CASE.initialDensity * ones(NUM.numCells, 1); % 在道路 30%~50% 区间把密度翻倍模拟初始拥堵 idx round(NUM.numCells * 0.3) : round(NUM.numCells * 0.5); rho0(idx) rho0(idx) * 2; state.x x0; state.rho rho0; state.v MODEL.freeSpeed * ones(NUM.numCells, 1);逻辑说明把道路离散成NUM.numCells个格子再在中间段设置一个高密度块。这样宏观程序一开始就会向右释放消散波、向左释放冲击波不需要额外做扰动。微观初始化更简单一般是等间距放置车辆再给同一个初始速度。参数说明initialDensity的单位要和roadLength、numCells对齐如果写成辆/km 而道路长度用米算出来的车辆数会差三个量级这是最常见的初始化错误。2.3 TrafficFlowSimulation.m 的主循环骨架宏观和微观的推进函数差别很大但主循环骨架是统一的初始化、步进、按采样周期记录日志。% TrafficFlowSimulation.m 主循环骨架 for it 1 : NUM.simTime / NUM.dt if strcmp(MODEL.type, macro) state TimeStepMacro(state, CASEPARS, MODELPARS, NUMPARS); else state TimeStepMicro(state, CASEPARS, MODELPARS, NUMPARS); end if mod(it, PLOTPARS.sampleEvery) 0 logs.t(end 1) it * NUM.dt; logs.rho(end 1, :) state.rho; logs.v(end 1, :) state.v; end end逻辑说明TimeStepMacro和TimeStepMicro都接收 state 和全部参数返回更新后的 state。PLOTPARS.sampleEvery决定日志密度设为 1 表示每个时间步都记录适合短时仿真对微观程序来说由于时间步长通常远小于宏观建议设 5 或 10否则内存里会堆满大矩阵。这里有一个常见边界问题NUM.simTime / NUM.dt不一定是整数直接用for it 1 : N会把尾步丢掉。我一般会改成while t NUM.simTime每步更新t it * NUM.dt保证最后一段时刻也被记录。3. FundamentalDiagram.m 与 TimeStepMacro.mLWR 守恒方程的数值落地3.1 三参数关系与 FundamentalDiagram.m交通流三参数关系是判断程序是否跑对的第一道标准。宏观模型里最基本的满足关系是q ρ * v也就是流量等于密度乘以速度。常见线性速度-密度关系为v vf * (1 - ρ / ρjam)代入后得到抛物线的流量-密度关系。实际高速路基本图更像三角形因为自由流段和拥堵段的斜率不同但教科书配套程序往往先用线性关系讲原理。% FundamentalDiagram.m 中 Q-K 关系绘制核心部分 rho linspace(0, 200, 200); % 密度轴辆/km vf 80; % 自由流速度km/h rhoJam 200; % 阻塞密度辆/km v vf * (1 - rho / rhoJam); % 线性速度-密度关系 q rho .* v; % 流量 密度 * 速度 plot(rho, q, LineWidth, 1.5); xlabel(密度 K (辆/km)); ylabel(流量 Q (辆/h)); grid on;逻辑说明当密度为零时流量为零当密度等于阻塞密度时速度为零、流量也为零最大流量出现在rho rhoJam / 2处对应临界密度。程序里可以把这个理论曲线保存下来再把仿真输出的密度-流量散点叠上去基本图上出现明显偏离就说明初始场、边界或数值格式有问题这是最直接的自检方式。三参数关系表达式关键含义速度-密度 V-Kv v(ρ)密度升高速度下降流量-密度 Q-Kq ρ v(ρ)临界密度处流量最大流量-速度 Q-Vq ρ(v) v高速低流、低速高流两个分支3.2 TimeStepMacro.m 的迎风差分实现LWR 守恒方程可以写作∂ρ/∂t ∂q/∂x 0。离散时不能直接用中心差分因为交通波有明确方向信息从上游传向下游所以要用迎风格式取上风侧流量。% TimeStepMacro.m 中的一阶迎风更新片段 function state TimeStepMacro(state, CASE, MODEL, NUM) rho state.rho; v Speed(rho, MODEL); % 平衡速度 q rho .* v; % 当前流量场 % 迎风取上游格点流量作为当前格点边界流量 q_flux [q(1); q(1:end-1)]; dqdx (q - q_flux) / NUM.dx; rho_new rho - NUM.dt * dqdx; % 物理约束密度不能越界 rho_new max(0, min(CASE.rhoJam, rho_new)); state.rho rho_new; state.v Speed(rho_new, MODEL); end逻辑说明q_flux把每个格点的上游邻居当作入流流量这是典型的一阶迎风。rho_new rho - dt * dqdx就是守恒方程显式时间推进。密度截断保证数值解始终落在物理区间内。这个格式能复现激波但会有数值耗散激波前沿会被抹平一点要减小耗散可以把q_flux换成 Godunov 精确黎曼解或者加 Minmod 斜率限制本质都是更准确地计算网格边界处的流量。参数说明稳定性要求满足 CFL 条件也就是max(v) * dt / dx 1。这里max(v)直接取MODEL.freeSpeed即可不需要在运行时算最大值。如果看到密度出现负值或激波前方振荡先看NUM.dt是否偏大再看边界条件是否写对。闭环道路必须做周期边界也就是最后一个格点的上游邻居要回到最后一个或第一个格点不能简单补零开环道路则要给定入口流量否则车辆在边界处会凭空消失。3.3 Speed.m 与 Position2Spacing.m 的角色差别Speed.m在宏观和微观两个分支里的含义不同。宏观分支输入密度输出平衡速度微观分支输入车头间距输出期望速度。Position2Spacing.m则是把位置向量转成车头间距数组为微观模型准备数据。下面这种写法把两类模型统一到一个接口里% Speed.m 的等价实现 function v Speed(rhoOrSpacing, MODEL) if strcmp(MODEL.type, macro) v MODEL.freeSpeed * max(0, 1 - rhoOrSpacing / MODEL.rhoJam); else v MODEL.freeSpeed * (1 - exp(-rhoOrSpacing / MODEL.s0)); end end逻辑说明宏观分支里max(0, ...)保证密度超过阻塞密度后速度也不会变负微观分支用指数松弛表示间距越大越接近自由流速度。两个分支都依赖MODELPARS里的参数所以同一个Speed.m可以被宏观和微观两个推进函数调用。Position2Spacing.m则不同如果输入是宏观格点位置输出平均间距如果输入是微观车辆位置输出相邻车辆的车头间距。宏观里平均间距约等于1/ρ微观里需要减掉车长才是有效间距这个差别很容易在跨模型对比时踩坑。4. TimeStepMicro 与 DoDisturbance微观跟驰模型如何复现走走停停4.1 微观状态量更新从间距到速度再到位置微观模型的状态量很简单每辆车有位置x、速度v车辆之间通过车头间距互相影响。TimeStepMicro.m每次推进分两步先算期望速度或加速度再更新速度和位置。下面是一段等效的松弛跟驰模型实现。% TimeStepMicro.m 的车辆循环等效代码 function state TimeStepMicro(state, CASE, MODEL, NUM) x state.x; % 单车位置向量 v state.v; % 单车速度向量 nv length(x); % 车辆数 spacing Position2Spacing(x, nv); % 车头间距 v_new zeros(nv, 1); for i 1 : nv s spacing(i); vDes MODEL.freeSpeed * (1 - exp(-s / MODEL.s0)); % 期望速度 a (vDes - v(i)) / MODEL.tau; % 松弛加速度 v_new(i) max(0, v(i) a * NUM.dt); % 速度更新 end state.v v_new; state.x x v_new * NUM.dt; % 位置更新 state.spacing Position2Spacing(state.x, nv); end逻辑说明这里用的是松弛跟驰模型核心是期望速度公式vDes vf * (1 - exp(-s / s0))。当间距很小时期望速度接近零车辆减速当间距足够大时期望速度接近自由流速度。MODEL.tau是松弛时间决定驾驶员把速度调整到期望速度的快慢。注意位置更新用的是更新后的v_new相当于一种半隐式处理比用旧速度更不易出现车辆重叠。参数说明s0通常取 5~20 mtau取 1~3 s。tau太小时车辆起步过猛tau太大时车流重新加速很慢正好对应走停状态里“消散慢”的过程。微观程序的NUM.dt通常取 0.01~0.1 s比宏观小一个量级所以总步数会多很多。4.2 DoDisturbance.m在车队中部注入一次速度扰动微观仿真的难点不是让车流跑起来而是制造出真实的扰动。DoDisturbance.m在某个时间窗口内压低某辆车速度相当于人为制造一次刹车灯波。% DoDisturbance.m 的等效实现 function state DoDisturbance(state, it, CASE, MODEL, NUM) t it * NUM.dt; if t CASE.disturbTimeStart t CASE.disturbTimeEnd idx round(length(state.x) * CASE.disturbPositionRatio); % 把第 idx 辆车限制到低速模拟前车短暂减速 state.v(idx) min(state.v(idx), CASE.disturbSpeed); end end在主循环里调用时要和TimeStepMicro配合if strcmp(MODEL.type, micro) state DoDisturbance(state, it, CASE, MODEL, NUM); state TimeStepMicro(state, CASE, MODEL, NUM); end逻辑说明扰动位置disturbPositionRatio通常取 0.5~0.7不要放在头车。放在头车只会让单车减速放在车队中部会让后车因为间距变小而逐辆减速同时前车逐渐拉开形成局部密度升高并向上游传播。若扰动窗口持续 3~5 s模拟的是刹车灯波若持续到仿真结束模拟的是持续瓶颈。注意DoDisturbance之后要立刻执行一次TimeStepMicro不要连续调用两次扰动否则同一时刻速度会被压两次。4.3 固定坐标与移动坐标两种看拥堵演化的方式微观输出是每个时刻的所有车辆位置直接看轨迹图容易眼花。这套程序提供了两个观测视角固定坐标和移动坐标。PlotFixedCoordinates.m用于固定路面位置观察密度或速度随时间变化PlotMovingCoordinates.m把坐标系以参考速度v_ref平移方便跟踪密度波传播。% PreprocessPlotMovingCoordinates.m 的核心变换 function xm PreprocessPlotMovingCoordinates(x, t, v_ref) xm x - v_ref * t; % 伽利略变换到移动坐标系 end然后可以用散点图绘制移动坐标下的轨迹% PlotMovingCoordinates.m 绘制移动坐标下的轨迹点 scatter(xm, t, 6, v, filled); % 颜色表示速度 xlabel(移动坐标 x - v_{ref} t (m)); ylabel(时间 t (s)); colorbar;逻辑说明v_ref选在自由流速度附近时自由行驶的车辆轨迹接近竖直走停状态下的减速波会变成明显斜线因为波速通常低于自由流速度。如果所有轨迹都向同一个方向倾斜说明v_ref选得不合适可以取扰动窗口内波速的平均值重新变换。固定坐标绘图则更简单先把车辆位置插值到规则网格上再用imagesc画时空图PreprocessPlotFixedCoordinates.m做的主要就是这种插值网格间距太大会抹掉短时拥堵太小又会放大插值噪声。5. ESM 场景切换与守恒检验把程序改造成参数实验平台5.1 用 ESM 做批量场景扫描ESM 最大的价值是场景可复用。把 ESM 放在一个批量脚本里循环执行就能对反应时间、自由流速度、扰动窗口等参数做扫描。注意 MATLAB 的脚本循环里变量会残留跑下一组场景前要重新加载四组 PARS。% runExperiments.m 简单批量扫描 reactTimes [0.5, 1.0, 2.0]; for i 1 : length(reactTimes) ESM; % 载入默认场景 MODELPARS.tau reactTimes(i); % 覆盖反应时间 state Initialization(CASEPARS, MODELPARS, NUMPARS); [state, logs] TrafficFlowSimulation(state, CASEPARS, MODELPARS, NUMPARS, PLOTPARS); results(i).tau reactTimes(i); results(i).meanFlow mean(logs.q(end - 200 : end)); end这里ESM每次都会重新执行变量名相同所以会覆盖上一次的MODELPARS这是脚本式参数装配的副作用。如果要跑更严格的实验建议把 ESM 改成函数[CASE, MODEL, NUM, PLOT] SetupScenario(scenarioID)返回参数结构体而不是往工作区塞变量。5.2 用守恒律验证数值解闭环道路下车辆数应该守恒。宏观程序里车辆总数等于密度对空间的积分也就是sum(rho) * dx。每跑完一步都检查这个值如果漂移超过千分之一说明边界或差分格式有问题。% 闭环边界下车辆数守恒检查 totalVehicle sum(state.rho) * NUM.dx; if abs(totalVehicle - totalVehicle0) totalVehicle0 * 1e-3 warning(车辆数不守恒: %.4f - %.4f, totalVehicle0, totalVehicle); end微观模型没有这个守恒量问题因为车辆是离散对象只要不超越前车总数就不变但微观程序要检查最小间距是否小于车长否则说明发生了碰撞。5.3 用“扰动窗口提前 20 秒”验证冲击波速度这套程序里最有意思的验证方式是只改CASE.disturbTimeStart其他参数不动。把扰动窗口提前 20 s再分别用固定坐标和移动坐标画两张时空图对比两次仿真的密度波前沿位置。你会看到拥堵排队头的反向传播速度基本不变这个速度在 LWR 模型里对应特征速度在排队论里对应冲击波速度。若用移动坐标系把v_ref设成这个反向传播速度的绝对值密度波前沿会变成一条接近竖直的线说明坐标变换和波速测量都对上了。这一步跑通之后这套交通流 MATLAB 程序就不再是“别人的代码”而是你可以自由改参数、做实验、出图写结论的仿真平台。本文还有配套的精品资源点击获取