非均衡动态交通分配模型MATLAB实现：从原理到代码拆解

简介非均衡动态交通分配模型的MATLAB实现是一套面向交通规划、智能交通及大学相关专业学生的完整代码资源。它通过基础模型I、II及其IS变体结合动态网络加载脚本揭示了非均衡条件下交通流时空调配与路径选择的分析方法为实现拥堵预测与设施优化提供了有效工具。压缩包共13个文件其中8个.m脚本实现模型核心算法与主程序5个.mat文件提供SiouxFalls等路网拓扑、OD需求及路径流量数据包体仅1.19MB便于精准获取与本地运行。代码采用参数化编程注释清晰支持MATLAB 2014/2019a/2024a等版本用户可直接运行附赠案例数据也可按需调整流量、时间、容量等参数以适配不同模拟场景。对于计算机、电子信息工程等专业的大学生这份资源可支撑课程设计、期末大作业或毕业设计对于相关领域研究者与工程师又提供了可复现的实验基础与对比参考。目前已有60人学习该资源具备实践借鉴价值。1. 为什么动态交通分配模型要区分“均衡”与“非均衡”读完城市交通分配相关文献你会发现大多论文都在解 Wardrop 均衡用户都选最短路径最终谁也不能通过对换路径降低自己的出行成本。可落到工程上均衡模型的问题很明显——求解需要的路径枚举和固定点迭代在中等规模路网上就非常慢更别说是动态场景时段切得越细变量规模涨得越快。非均衡模型不追求严格均衡用增量加载、连续平均法这类启发式逼近算得快、结构清晰还能在滚动预测、动态信号优化这类对实时性敏感的场合里直接用。这套 MATLAB 实现把非均衡动态交通分配DTD的骨架拆得很清楚Base_Model_I 与 Base_Model_II 两套基准模型配合 SiouxFalls 标准路网和 OD 数据跑起来就能看到流量在时轴上逐步传播的过程。如果你是做课程设计、毕设刚开始接触 DTA或者工程师想快速验证路网改造方案这份代码能帮你省掉两周搭框架的时间。接下来我会按模型原理、基准脚本逻辑、数据文件和实际改路网四个层面把它拆透。2. 非均衡 DTA 的模型拆解从 UE/SUE 说起2.1 均衡模型的理论底座与不适用场景经典静态交通分配里用户平衡User Equilibrium, UE要求同一 OD 对间所有被使用路径的阻抗相等且等于最短阻抗。系统最优SO则要求总阻抗最小。UE 有严格数学定义对应变分不等式的解Frank-Wolfe 算法是标准解法。可一旦放到动态环境Dynamic Traffic Assignment, DTA路径行程时间变成时变函数需要把时间离散成区间UE 的求解就变成了复杂的高度非线性问题收敛性很难保证。这里说“非均衡”不是指模型完全忽略用户选择而是放弃严格的均衡条件。工程上常见的做法有两种一是把 OD 需求按比例分到几条预生成的路径上不做迭代收敛二是用迭代近似比如连续平均法MSA每步用旧的路径费用估算新路径用递减步长合并流量。项目里的 Base_Model_I 和 Base_Model_II 基本可以对应这两条思路。2.2 非均衡路径选择规则增量加载与迭代式加载增量加载Incremental Assignment的做法是把 OD 矩阵分成若干份每次加载一小部分根据当前路网流量计算行程时间再更新路径费用再把下一份流出需求加进路网。所以早加载的流量会影响后来者的选路决策这比一次在全空网上加载要接近现实。另一种是迭代比例分配类似 MSA先按自由流时间算最短路径并全量加载然后按当前流量重新算最短路径再算新的路径流量用 1/k 这样的步长去和上一次结果做平均。k 是迭代次数这样能在固定迭代步数内稳定输出。项目中的 Base_Model_I_IS.m 和 Base_Model_II_IS.m 里的 “IS” 我常用作 “Incremental Strategy” 或 “Iterative Strategy” 理解前两个脚本对应基础分配逻辑后两个脚本在同一个框架上加了策略控制。看代码时你会发现两类脚本共用同一份 OD_info.mat 和 Network_planning_parameters.mat说明数据层与算法层是解耦的。2.3 MATLAB 里如何表示路网和 OD 矩阵打开 SiouxFalls6180_pp.mat里面通常不是 MATLAB 的 graph 对象而是交通规划软件里常见的节点-路段结构。我一般先加载用 whos 看一下变量名和维度load(SiouxFalls6180_pp.mat); whos % 查看所有变量名、尺寸、类型假设其中有net和od两个结构体原数据文件名带 6180说明有 6180 个 OD 对节点数和路段数视网络版本而定典型的是% 路网结构体一般包含 % net.node: N x 2节点编号和坐标 % net.link: M x 4[起点编号 终点编号 自由流时间 容量] % OD 矩阵full matrix 或 sparse matrixOD(i,j) 表示 i 到 j 的出行需求 net struct(... node, [1 1 0; 2 2 1; 3 3 0], ... % 示例实际以文件为准 link, [1 2 3.5 100; 2 3 4.0 120], ... od, [0 10; 5 0]);真实数据里Network_planning_parameters.mat约束了模型的时间片长度、仿真周期数、拥挤系数等OD_info.mat里存 OD 矩阵和需求时变曲线。加载后务必检查isfield再进入主循环否则动态加载时读不到字段会直接报Index exceeds matrix dimensions。3. Base_Model_I 与 Base_Model_II两个基准模型的实现逻辑3.1 Base_Model_I固定路径集合下的流量加载Base_Model_I.m 负责做静态或准动态的非均衡分配先把 OD 矩阵按比例切块例如把总量分为 5 份20%、20%、20%、20%、20%每一轮选择当前行程时间下的最短路径把这一份需求放到路径上然后更新该路径经过路段的流量从而增加通行时间。这个逻辑在代码里的判断条件就是是否按轮次循环% Base_Model_I.m 核心循环伪结构示例实际以源码为准 numSlice size(OD_matrix, 1); % OD矩阵维度 nParts 5; % 分割份数参数化后可改 linkFlow zeros(nLinks, 1); % 路段流量 for k 1:nParts partOD OD_matrix / nParts; % 每个切片的需求量 [pathSet, travelTime] shortpath(net, partOD); % 当前费用下最短路 linkFlow updateLinkFlow(linkFlow, pathSet, partOD); net.link(:,3) BPR(net.link(:,3), linkFlow); % 更新阻抗 end这里的shortpath是典型的 Dijkstra 或 yen 算法项目里可能内嵌在某个函数里。前半句说明每份切片都按当前路网状态选路所以前一份流量造成拥挤后后续流量会被引导到替代路径这个过程模拟了交通流的动态转移。参数说明nParts分割份数越大结果越接近均衡但计算量线性增加net.link(:,3)是自由流时间BPR函数用 BPR 曲线t t0 * (1 alpha * (q/c)^beta)计算拥挤后的阻抗默认 alpha0.15beta4。实际项目里这两参数也放在 Network_planning_parameters.mat 中。3.2 逐步迭代的非均衡分配Base_Model_II 的迭代框架Base_Model_II.m 走的是迭代平均路线。它与 I 型不同不是在一次分配里切分 OD 量而是做多轮加载每轮结束后用 MSA 步长平均上一轮的路径流量与当前轮的结果。这样即使路径费用变化输出流量也不会震荡太剧烈。% Base_Model_II.m 迭代框架示例结构 maxIter 30; % 最大迭代次数 rho 0.15; alpha 4; % BPR 参数 X_old zeros(nLinks, 1); for iter 1:maxIter [pathNum, pathTime] shortestPath4All(net, OD_matrix); X_new assignODToLinks(pathNum, OD_matrix); % 全量加载 step 1 / (iter 1); % MSA 步长 X_msa (1 - step) * X_old step * X_new; % 更新路段时间 net.link(:, 3) freeTime .* (1 rho * (X_msa ./ capacity).^alpha); % 收敛判断前后两轮流量相对误差小于阈值则停止 if norm(X_msa - X_old, inf) / norm(X_old, inf) 1e-4 break; end X_old X_msa; end每次迭代的全量加载算是“盲目加载”但由于步长递减最终输出在理论上可以靠拢到一个稳定状态。注意1/(iter1)是 MSA 的关键如果步长不衰减后期流量会一直在真值附近抖动。项目里 Base_Model_II_IS.m 增加了对迭代次数的控制逻辑比如进行 15 次迭代后强制输出中间结果方便观察收敛过程。这节再给出一个实际运行中提取收敛指标的片段% 记录每轮最大流量差绘制收敛曲线 convergence(iter) norm(X_msa - X_old, inf); if iter 1 semilogy(1:iter, convergence); xlabel(迭代轮次); ylabel(最大流量差对数坐标); grid on; pause(0.02); end这个图能直观判断 MSA 是否在收敛。非均衡模型的收敛阈值不像严格 UE 那样有理论保证看曲线进入平台期就可以停止。3.3 参数化编程哪些参数可以直接改这套代码的优点是参数集中。下面把最常改的参数和取值范围列出来方便你调整场景。参数名称对应文件/位置说明推荐值/调整方向时间切片数Network_planning_parameters.mat动态模型中一个仿真周期划分成多少段6-24片段越多越贴近连续流但内存涨快仿真总时长Network_planning_parameters.mat模拟多少分钟的路网运行状态60-120 分钟视 OD 时变曲线覆盖范围而定OD 分割份数Base_Model_I.m 内变量增量加载的轮次4-8大于 10 效果接近均衡但耗时上升最大迭代数Base_Model_II.m 内变量 maxIterMSA 循环上限20-50看网络规模BPR alphaNetwork_planning_parameters.mat拥堵程度对时间的影响系数0.15 是标准值城市密集区可取 0.2BPR betaNetwork_planning_parameters.mat拥堵函数的非线性指数2-5beta 越大越敏感收敛阈值Base_Model_II.m前后两轮流量差的最大允许值1e-4 到 1e-3课程设计用 1e-3 可接受我在实际跑的时候会把 OD 分割份数和最大迭代次数拆出来做两层循环先粗调后细调。.mat里的参数建议用脚本统一写入不要手动改二进制文件比如param.timeSlice 12; param.simHorizon 60; param.bprAlpha 0.15; param.bprBeta 4; save(Network_planning_parameters.mat, param);这样就避免了在 MATLAB 2014 和 2024 之间切换时旧版保存的 struct 字段顺序影响加载。4. SiouxFalls 6180_pp 数据集与动态网络加载4.1 数据文件清单与网络拓扑资源包里几个.mat文件对应不同的数据层。SiouxFalls6180_pp.mat 与 SiouxFalls6180_pp_68.mat 是从标准 Sioux Falls 网络裁出来的两个变体前者含 24 个节点、76 条单向路段不同版本略有差异后者标注了 68 个路径或 OD 相关的信息。我的判断是_68这个版本给的是预设的 68 条候选路径集用来做固定路径分配的避免每次重新生成路径。文件数据内容用途OD_info.matOD 矩阵、各 OD 的时变需求模型输入Network_planning_parameters.mat仿真时间、切片数、BPR 参数模型参数配置Path_flow_data.mat路径与路段关联矩阵、路径流量结果校验SiouxFalls6180_pp.mat路网拓扑、路段属性、节点坐标基础路网SiouxFalls6180_pp_68.mat带 68 条候选路径的路网路径预生成场景跑通主程序前先确认这几个文件都在同一个目录并检查which路径fnames {OD_info.mat,Path_flow_data.mat,Network_planning_parameters.mat, ... SiouxFalls6180_pp.mat,SiouxFalls6180_pp_68.mat}; for i 1:length(fnames) fprintf(%s: %s\n, fnames{i}, which(fnames{i})); end如果输出显示空白说明文件不在搜索路径里要在脚本开头用addpath加上相对路径或绝对路径。4.2 DYNAMIC_NETWORK_LOADING.m 的工作机制DYNAMIC_NETWORK_LOADING.m 是动态加载的核心函数。它把 OD 需求按时间片拆分随着时间推进流量进入路网经过一个个路段在路段出口形成流量分布。常见做法是用了一个流出率函数比如根据 Smith 模型或 CTM元胞传输模型的简化版本把路段流量按时间推进function [linkState, linkTravel] DYNAMIC_NETWORK_LOADING(OD_series, net, param) % linkState: 每个时间段每个路段的累积流量 % linkTravel: 每个路段在每个时间片的实际行驶时间 nLinks size(net.link, 1); nTime param.timeSlice; linkState zeros(nLinks, nTime); linkTravel zeros(nLinks, nTime); for t 1:nTime inflow OD_series(:, :, t); % 该时段进入路网的出行需求 linkState(:, t) propagateTraffic(linkState(:, t-1), inflow, net); % 计算该时段车流经历的实际时间包含排队延迟与拥挤影响 linkTravel(:, t) net.link(:, 3) .* (1 param.bprAlpha * ... (linkState(:, t)./net.link(:, 4)).^param.bprBeta); end end这段是骨架代码与项目中将net.link(:,4)作为容量字段的假设一致。动态加载的意义在于同一条路段在不同时间段进入的车流离开时经历的时间是不同的。如果只关心最终总流量直接用静态加载即可但你想知道“早高峰出门的车多快能到收费站”这类时变问题就必须依赖 linkTravel 这个输出。4.3 从 MATLAB 控制台到结果可视化运行非均衡模型通常从主脚本入口启动。比如你的工作目录里有Base_Model_I.m直接在命令窗口执行cd(你的解压目录/DTD-master); % 假设 DTD-master 是主文件夹 run(Base_Model_I.m);程序运行期间如果用到了dispstat.m那是用来输出动态进度信息的小工具它能在命令行同一行刷新百分比进度条不影响计算逻辑。如果你在无桌面环境运行脚本会跳过显示部分只保存结果。跑完后把动态加载的流量拿出来画时空图可以直观看到流量是否随时间向网络深处推进% 假设 linkState 是 DYNAMIC_NETWORK_LOADING 输出的 nLinks x nTime 矩阵 [t_axis, link_index] meshgrid(1:size(linkState,2), 1:size(linkState,1)); figure; surf(t_axis, link_index, linkState); xlabel(时间片); ylabel(路段编号); zlabel(流量); grid on; title(动态交通加载流量时空分布);图上如果某几个路段永远是同一时刻出现尖峰说明 OD 加载粒度太粗或时间切片划分和网络长度不匹配。适当增加param.timeSlice或缩小每个时间片的 OD 量能平滑曲线。5. 把模型改到自己路网上数据准备、标定与常见坑5.1 自定义路网的数据格式转换要换到自己的路网最关键的一步是把 CAD 或 GIS 里的道路网络转成 MATLAB 结构体。我一般先转成两个 CSVnode_file.csv节点号、x、y和link_file.csv起点、终点、自由流时间、容量。然后在 MATLAB 里批量导入。% 从 CSV 构建 net 结构体字段与原始数据保持一致 nodes readtable(node_file.csv); links readtable(link_file.csv); net.node table2array(nodes); net.link [links.begin_node, links.end_node, ... links.free_flow_time, links.capacity]; % 检查路段方向如果 link_file 里同时有双向路段确认 a-b 和 b-a 都存在 if any(links.begin_node links.end_node) warning(存在反向路段缺失请检查 link_file 是否包含双向边); end转换完成后用 0 流量状态跑一次 Base_Model_I看路网里是否每条边都有自由流时间。常见错误是把free_flow_time当成 0导致模型里所有最短路都选一跳完成提前结束时序模拟。% 检查是否有零通行时间的路段 bad find(net.link(:,3) 0); if ~isempty(bad) fprintf(警告%d 条路段自由流时间 0\n, numel(bad)); net.link(bad, 3) 1; % 保底值后面按实际再标定 end5.2 非均衡分配结果的一致性验证模型能跑通不等于结果的数值正确。我最常用的是“传播守恒检查”从发出端累计算每个时间切片进入路网的总车数减去到达端累计离开的车数差值应该等于停留在路网内部的车数且不小于 0。% 守恒检查所有时间片总流入 总流出 路网内滞留 totalIn sum(sum(inflow_series)); % inflow_series 为 nOdPairs x nTime 的进入量 totalOut sum(sum(outflow_series)); % outflow_series 为 nLinks x nTime 的出口流量 staying totalIn - totalOut; if staying 0 || staying 0.05 * totalIn error(流量不守恒检查 DYNAMIC_NETWORK_LOADING 的流出率是否与流入率匹配); else fprintf(守恒验证通过路网内滞留车辆 %g\n, staying); end另外拿 Base_Model_I 和 Base_Model_II 的表层结果对比也有价值。因为两种算法思想不同最终路段流量会有偏差但偏差超过 30% 时说明 OD 分割份数太小或迭代步长太大应该调整参数方向。5.3 关于 .mat 文件版本、路径依赖与避开乱码MATLAB 2014 保存的 mat 文件是 v7 格式2024a 仍然兼容但反过来 2024a 高版本保存的 v7.3 格式文件在旧版 MATLAB 2014 上可能读不出来。如果打算在两个版本间交叉验证保存数据时明确指定版本。save(SiouxFalls6180_pp.mat, net, od, -v7); % 低版本最高兼容另外项目解压路径里别带中文和空格。MATLAB 的 m 文件编码默认跟随系统 locale如果文件里含中文注释在 Linux 或西文系统上打开可能乱码导致脚本无法执行。稳妥做法是打开脚本后用“另存为”编码选 UTF-8再在运行前用edit Base_Model_I.m检查第一行注释是否正常。最后一个技巧用dispstat.m输出棋盘格日志时它会把光标固定在某一列刷新遇到 MATLAB Desktop 卡顿是正常的不用管它。如果嫌刷屏太快把它换成fprintf每五行打一次进度即可if mod(iter, 5) 0 fprintf(Iteration %03d, 最大流量误差 %g\n, iter, convergence(iter)); endmkdocs 类的日志格式对对比迭代很有用直接复制到命令行窗口就能一眼定位到哪轮迭代开始收敛比盯进度条高效得多。本文还有配套的精品资源点击获取