咱先拿简化过的「江南XX互通站」2021上半年到2023下半年的半年度数据凑凑活，真实数据不敢碰红线哈，但算法是一样的

MATLAB程序高速公路预测道路货运量和客运量预测货运量万吨52.1、61.3、70.8、68.5、76.2、82.7客运量万人次120.3、145.7、138.2、172.1、169.4、198.6朋友一开始找我的时候说他们领导突然要半年后的规划依据要“明年上半年可能会压垮哪段车道”这种粗略但有用的方向自己只会用Excel拖个散点图加趋势线偏差大到连他自己都不信。趋势线那种纯拟合线根本抓不住小样本时间序列可能带的“非平稳但有增长惯性”的味儿——毕竟咱高速不是纯自然科学实验有政策、有假期调休、有附近工厂搬进来搬走的小波动但核心趋势不会跑太远对吧这时候MATLAB里的灰色预测GM(1,1)简直是救星数据量少只要4个以上就行、不用正态分布假设、不用解释一堆复杂变量对这种“摸不准宏观完全变量但内部有小规律”的小场景太适配了。先上一段货客一起跑的MATLAB简化代码注释里我尽量用大白话别整那些“累减生成矩阵求逆”绕人的公式开头% 江南XX互通站凑活版预测GM(1,1) clear;clc;close all; % 先把后台擦干净强迫症友好 % 第一步放数据记住GM(1,1)喜欢列向量不管单变量多变量都竖着放 % 第一列是货运量第二列是客运量 raw_data [ 52.1, 120.3; 61.3, 145.7; 70.8, 138.2; 68.5, 172.1; 76.2, 169.4; 82.7, 198.6 ]; n size(raw_data, 1); % n是历史数据的时间点数量这里是6个半年度 m size(raw_data, 2); % m是预测的变量数2个货、客 future_num 2; % 想预测多久咱给朋友凑2个半年度的刚好明年全年的规划片段 % 第二步GM(1,1)核心预处理累加生成AGO % 为啥要累加因为把原始带波动的“锯齿线”拉成“平滑上升/下降线”灰色模型吃这个 ago_data cumsum(raw_data); % 插一句哈去年MATLAB好像更新了一个灰色预测工具箱函数gm11但我还是喜欢手写半核心不然连参数怎么调或者有没有明显调歪都不知道为啥说手写半核心踏实后面讲参数a发展系数和b灰作用量的时候你就懂了——如果发展系数a的绝对值超过0.5那GM(1,1)的长期预测就完全不准最多看个1-2个时间点刚好符合朋友的需求领导要的就是明年上半年的应急参考嘛。MATLAB程序高速公路预测道路货运量和客运量预测继续写参数计算和预测的部分% 第三步构造灰色模型的矩阵B和向量Y % 这个部分别纠结B矩阵为啥长这样记住这是邓聚龙老先生当年摸出来的“最优拟合平滑线的结构”就行 for i 1:m % 对货、客分别算单变量分开处理更稳 % 构造第i个变量的累减序列的中间值也叫紧邻均值生成 z 0.5*(ago_data(2:end,i) ago_data(1:end-1,i)); % 构造对应第i个变量的B矩阵和Y向量 B [-z, ones(n-1,1)]; Y raw_data(2:end,i); % 矩阵求逆找a和b公式其实就是最小二乘法拟合平滑AGO线的斜率截距变形 ab inv(B*B)*B*Y; a_i ab(1); b_i ab(2); % 第四步先预测累加序列AGO的未来值再累减还原成原始值别直接预测原始值会崩 x_ago_pred zeros(nfuture_num,1); x_ago_pred(1) raw_data(1,i); % 第一个累加值肯定等于原始第一个值 for k 1:nfuture_num-1 x_ago_pred(k1) (raw_data(1,i)-b_i/a_i)*exp(-a_i*k) b_i/a_i; end % 累减还原AGO的第k1个减第k个就是预测的原始第k1个值 x_pred(:,i) [raw_data(1,i); diff(x_ago_pred)]; % 第五步算个平均相对误差MAPE看看模型靠不靠谱别瞎给朋友 mape_i mean(abs((raw_data(2:end,i)-x_pred(2:n,i))./raw_data(2:end,i)))*100; fprintf(%s的MAPE是%.2f%%\n,{货运量,客运量}{i},mape_i); end好现在我们把这段代码粘到MATLAB里跑结果大概是这样货运量的MAPE是3.72%客运量的MAPE是4.15%这两个误差都在5%以内对于这种凑活的简化数据来说已经很够意思了领导绝对不会骂你瞎编。再上最后一段画图的代码可视化才是给小白哦不给领导看的神器% 第六步画图时间轴用半年度标清晰明了 time [2021.5, 2022, 2022.5, 2023, 2023.5, 2024, 2024.5]; % 前6个是历史后2个是预测 figure(Color,w,Position,[100,100,1000,400]); % 设个白色背景截图给领导好看 subplot(1,2,1); plot(time(1:n), raw_data(:,1), bo-, LineWidth,1.5, MarkerSize,8); hold on; plot(time(n:end), x_pred(n:end,1), r*--, LineWidth,1.5, MarkerSize,10); grid on; xlabel(年份半年度); ylabel(货运量万吨); title(江南XX互通站半年度货运量预测); legend(历史数据,预测数据,Location,northwest); subplot(1,2,2); plot(time(1:n), raw_data(:,2), go-, LineWidth,1.5, MarkerSize,8); hold on; plot(time(n:end), x_pred(n:end,2), m*--, LineWidth,1.5, MarkerSize,10); grid on; xlabel(年份半年度); ylabel(客运量万人次); title(江南XX互通站半年度客运量预测); legend(历史数据,预测数据,Location,northwest);跑出来的图应该是这样的货运量明年上半年大概能到89-90万吨下半年破95客运量明年上半年大概215-220万下半年235-240万。朋友拿到图直接给领导截了领导当天就安排了货车专用临时道的前期勘察完美收尾。哦对了再补个小提醒如果你的数据有明显的季节性波动比如春运五一国庆客运爆涨GM(1,1)可能不够用这时候可以加个季节性指数修正或者用MATLAB里的SARIMA工具箱函数但那又是另一个故事了。