关键词：MATLAB、数据建模、曲线拟合、最小二乘、polyfit、fitlm、系统辨识、ARX
实验课上测回来一堆散点，怎么把它变成一条能预测的曲线？本文沿着"曲线拟合 → 多元回归 → 系统辨识"这条工程上最常用的主线，把 MATLAB 数据建模的核心工具链串一遍：从最小二乘的正规方程推导，到 polyfit/fitlm/arx 的完整实战代码，再到一份"踩坑清单"。
一、数据建模在工程中的位置
1.1 机理建模 vs 数据驱动建模
工程上建立数学模型大体有两条路：
- 机理建模（白箱）：从物理定律出发推导方程，比如牛顿冷却定律、电路的基尔霍夫方程。优点是外推可靠、可解释性强；缺点是复杂系统（化工反应、电池老化）机理往往写不全。
- 数据驱动建模（黑箱/灰箱）：不问机理，直接从输入输出数据中学习映射关系 y = f(x)。拟合、回归、系统辨识、机器学习都属于这一类。
| 维度 | 机理建模 | 数据驱动建模 |
|---|---|---|
| 知识来源 | 物理/化学定律 | 实验数据 |
| 可解释性 | 强，参数有物理意义 | 弱，参数多为数学符号 |
| 外推能力 | 较好 | 差，训练域外不可靠 |
| 建模成本 | 高（需领域专家） | 低（有数据就能做） |
| 典型工具 | Simulink 物理建模 | Curve Fitting / System Identification Toolbox |
实际工程中两者常常混合使用：机理定结构、数据定参数（所谓"灰箱建模"）。本文聚焦数据驱动这条线。
1.2 为什么是 MATLAB
在工程建模领域，MATLAB 的地位来自三点：矩阵即一等公民的语法让最小二乘这类线性代数操作一行搞定；Curve Fitting、Statistics and Machine Learning、System Identification 三个 Toolbox 覆盖了从静态拟合到动态辨识的完整链条；拟合结果可以直接进 Simulink 做仿真。对工科生而言，它是"从实验数据到控制器"最短的路径。
二、曲线拟合：最小二乘的数学与工具
2.1 最小二乘原理与正规方程
给定 n 个观测点，假设模型为线性参数形式（注意：对参数线性，对 x 可以非线性，多项式就满足）：
写成矩阵形式，其中设计矩阵 X 的第 i 行为
。最小二乘的目标是最小化残差平方和：
对求导并令其为零，得到正规方程（Normal Equation）：
这就是所有拟合工具的数学内核。需要提醒的是：实际计算中 MATLAB 并不会真的去求逆（inv(X'*X)数值上不稳定），而是用 QR 分解或 SVD 求解——这也是反斜杠运算符\（即beta = X\y）优于手写正规方程的原因。
2.2 polyfit / polyval：最轻量的组合
polyfit(x, y, n)做 n 阶多项式拟合，返回按降幂排列的系数向量；polyval(p, x)用系数求值。R2020+ 推荐带三输出的用法以获得中心化与缩放（mu），消除病态：
[p, S, mu] = polyfit(x, y, 3); % mu 含均值和标准差，内部做 (x-mu(1))/mu(2) yhat = polyval(p, x, S, mu); % 求值时必须把 mu 传回去2.3 fit 函数与 Curve Fitting Toolbox
fit是更通用的入口，支持上百种模型库（多项式、指数、傅里叶、高斯、样条、自定义方程），并直接返回置信区间：
[fo, gof] = fit(x, y, 'poly2'); % 二次多项式 ci = confint(fo, 0.95); % 95% 参数置信区间 gof.rsquare % 拟合优度 R²交互式工具cftool可以拖拽调模型、实时看残差，适合探索阶段。
polyfit与fit怎么选？日常经验是：快速画图、嵌进脚本做批量处理，用polyfit——它返回纯数值向量，不依赖额外 Toolbox（基础 MATLAB 自带）；需要置信区间、拟合优度报表、非多项式模型（指数、有理式、自定义方程），用fit——它返回cfit对象，plot、coeffvalues、predint（预测区间）一条龙。两者底层同为最小二乘，数值结果一致。
2.4 拟合优度：R²、RMSE 与置信区间
判断拟合好坏不能只看"线穿得漂不漂亮"，常用三个定量指标：
| 指标 | 公式 | 含义与经验判据 |
|---|---|---|
| RMSE | 与 y 同量纲，越小越好；应小于测量噪声水平 | |
| R² | 解释方差比例，越接近 1 越好；但随阶数单调不减，不能单独用来选阶 | |
| 参数置信区间 | 区间跨过 0 说明该参数不显著，模型可能过度复杂 |
2.5 过拟合：阶数不是越高越好
n 个点总能被一个 n-1 阶多项式精确穿过，但那只是记住了噪声。过拟合的典型征兆：训练 RMSE 极小、高阶项系数巨大且置信区间很宽、曲线在数据点之间剧烈振荡（Runge 现象）。正确做法是划分训练集/验证集，用验证集误差选模型复杂度——下面的实战会完整演示。
三、代码实战一：温度传感器标定
场景：一只 NTC 热敏电阻温度传感器，标定实验在 0~100 ℃ 内测得 21 组数据（标准温度计读数t_ref为横轴，传感器电压换算的读数t_sen为纵轴），需要建立标定曲线并评估不同阶数的多项式。
%% 温度传感器标定：多项式阶数选择实战 clear; clc; close all; rng(42); % 固定随机种子，结果可复现 %% 1. 构造带噪声的标定数据（真实物理近似为二次关系） t_ref = linspace(0, 100, 21)'; % 标准温度计：0~100℃，21 个点 t_true = 0.002*t_ref.^2 + 0.9*t_ref + 1; % 假设的真实标定曲线 t_sen = t_true + 0.8*randn(size(t_ref)); % 叠加 σ=0.8℃ 的测量噪声 %% 2. 划分训练集 / 验证集（隔一个取一个，保证覆盖全温区） idx_val = 2:2:numel(t_ref); % 偶数下标做验证 x_tr = t_ref; x_tr(idx_val) = []; y_tr = t_sen; y_tr(idx_val) = []; x_va = t_ref(idx_val); y_va = t_sen(idx_val); %% 3. 对比 1 / 2 / 3 / 9 阶多项式 orders = [1 2 3 9]; figure; tiledlayout(2, 2, 'Padding', 'compact'); results = zeros(numel(orders), 3); % 存 [阶数, 训练RMSE, 验证RMSE] for k = 1:numel(orders) n = orders(k); [p, ~, mu] = polyfit(x_tr, y_tr, n); % 带中心化的拟合，数值更稳 r_tr = y_tr - polyval(p, x_tr, [], mu); % 训练残差 r_va = y_va - polyval(p, x_va, [], mu); % 验证残差 rmse_tr = sqrt(mean(r_tr.^2)); rmse_va = sqrt(mean(r_va.^2)); results(k, :) = [n, rmse_tr, rmse_va]; nexttile; plot(x_tr, y_tr, 'bo', x_va, y_va, 'rs'); hold on; fplot(@(x) polyval(p, x, [], mu), [0 100], 'k-'); title(sprintf('%d 阶: RMSE_{tr}=%.2f, RMSE_{va}=%.2f', n, rmse_tr, rmse_va)); legend('训练集', '验证集', '拟合曲线', 'Location', 'northwest'); end %% 4. 残差分析：残差应像白噪声，无趋势、无周期 [p2, ~, mu2] = polyfit(x_tr, y_tr, 2); % 选定的 2 阶模型 resid = y_tr - polyval(p2, x_tr, [], mu2); figure; subplot(2,1,1); stem(x_tr, resid, 'filled'); yline(0); title('二阶模型残差'); ylabel('残差 / ℃'); subplot(2,1,2); qqplot(resid); % QQ 图检验残差正态性 title('残差 QQ 图'); %% 5. 输出对比表 fprintf('阶数 训练RMSE 验证RMSE\n'); fprintf('%3d %8.3f %8.3f\n', results.');运行结果与分析（数值因噪声实现略有浮动，规律稳定）：
| 阶数 | 训练 RMSE / ℃ | 验证 RMSE / ℃ | 现象 |
|---|---|---|---|
| 1（直线） | ≈ 3.6 | ≈ 3.5 | 欠拟合，残差呈明显抛物线趋势 |
| 2 | ≈ 0.7 | ≈ 0.8 | 与噪声水平（σ=0.8）相当，最佳 |
| 3 | ≈ 0.7 | ≈ 0.9 | 训练误差不再下降，验证误差开始回升 |
| 9 | ≈ 0.5 | ≫ 2.0 | 训练误差虚低，验证误差爆炸——典型过拟合 |
三条可复用的经验：
- 训练误差降到噪声水平就停——再低就是在拟合噪声；
- 残差图比 R² 更有信息量：残差里只要有趋势或周期结构，就说明模型缺项；
- 9 阶拟合 11 个训练点，设计矩阵 X^TX 条件数极大，必须做
polyfit的中心化（mu），否则系数完全不可信。
四、回归建模与系统辨识
单变量拟合只是起点。工程问题里输出往往受多个因素影响（回归），或者系统本身有惯性、记忆（动态系统辨识）。
4.1 多元线性回归：regress 与 fitlm
经典接口regress(y, X)返回系数、置信区间、残差与统计量，其中X需自行添加全 1 列作为截距。现代接口fitlm（Statistics and Machine Learning Toolbox）更推荐——直接返回模型对象，系数检验、ANOVA、诊断图一应俱全：
mdl = fitlm(X, y); % X 为 n×p 矩阵，无需加截距列 mdl.Coefficients % 系数表：估计值、标准误、t 统计量、p 值 plotDiagnostics(mdl, 'cookd'); % Cook 距离，找强影响点fitlm还支持 Wilkinson 公式写法：fitlm(tbl, 'y ~ x1 + x2 + x1:x2')中x1:x2表示交互项，x1*x2等价于x1 + x2 + x1:x2。
4.2 逐步回归：stepwiselm
特征很多但不知道哪些有用时，stepwiselm按 p 值自动进/出变量，是快速的特征筛选手段：
mdl = stepwiselm(X, y, 'interactions', ... % 从纯交互模型出发逐步筛选 'Criterion', 'bic', ... % 用 BIC 而非 p 值，更防过拟合 'Upper', 'quadratic'); % 最多考虑到二次项注意：逐步回归是贪婪搜索，不保证全局最优，且对共线性敏感，结果应结合领域知识解读。
4.3 动态系统辨识：ARX 模型
如果数据是时间序列且系统有动态特性（如电机转速对电压的响应、室温对加热功率的响应），静态回归失效——此刻的输出还依赖过去的输入输出。System Identification Toolbox 中的 ARX（AutoRegressive with eXogenous input）模型描述为：
即
核心工作流三步：
data = iddata(y, u, Ts); % 打包：输出、输入、采样周期 model = arx(data, [na nb nk]); % na/nb: 多项式阶次；nk: 纯延迟 compare(data_val, model); % 在验证数据上对比仿真输出配套命令：iddata打包数据，arx估计参数，present(model)查看参数与不确定度，resid(model, data)做残差白性检验，compare用验证数据算 FIT 百分比。更复杂的系统可换ssest（状态空间）、nlarx（非线性 ARX）。
4.4 时间序列预测的思路
纯时间序列（无外生输入）可视为 ARX 的特例——AR 模型，思路一致：定阶（AIC/BIC 或交叉验证）→ 估计 → 残差白性检验（残差还含相关性说明信息没榨干）→ 多步预测时滚动回代。MATLAB 中可用ar、或 Econometrics Toolbox 的arima做 ARIMA 类模型。
定阶是最容易被敷衍的一步。除了网格搜索比验证 FIT，aic(model)可以直接给出信息准则值；实践中常把 AIC、BIC、验证误差三个判据放在一起看——三者指向同一个阶次时才敢拍板。另一个细节：做时间序列划分时必须按时间先后切，前段训练、后段验证，绝不能随机抽样，否则未来信息泄漏进训练集，验证误差会虚假地好看。
五、代码实战二：加热炉温度的 ARX 辨识
场景：电加热炉，输入u为加热功率（0~100%），输出y为炉温（℃），采样周期 2 s。用前 70% 数据辨识 ARX 模型，后 30% 验证。
%% 加热炉 ARX 系统辨识实战 clear; clc; close all; rng(7); %% 1. 生成仿真数据：一阶惯性 + 纯延迟 + 噪声（模拟真实采集） N = 600; Ts = 2; % 600 个样本，采样 2 s u = 50 + 30*randn(N,1); % 功率在 50% 附近随机扰动（保证激励充分） u = max(0, min(100, u)); y = zeros(N,1); for t = 3:N % y(t) = 0.9*y(t-1) + 0.08*u(t-2) + 噪声 y(t) = 0.9*y(t-1) + 0.08*u(t-2) + 0.3*randn; end %% 2. 打包并划分估计段 / 验证段 data = iddata(y, u, Ts, 'InputName', '功率', 'OutputName', '炉温'); Ne = round(0.7*N); data_est = data(1:Ne); % 前 70% 用于估计参数 data_val = data(Ne+1:N); % 后 30% 用于验证 %% 3. 网格搜索定阶：na, nb ∈ {1,2,3}, nk ∈ {1,2} best.fit = -Inf; for na = 1:3, for nb = 1:3, for nk = 1:2 m = arx(data_est, [na nb nk]); [~, fit] = compare(data_val, m); % 验证段 FIT 百分比 if fit > best.fit best.fit = fit; best.orders = [na nb nk]; end end, end, end fprintf('最优阶次 [na nb nk] = [%d %d %d], 验证 FIT = %.1f%%\n', ... best.orders, best.fit); %% 4. 用最优阶次重新估计并诊断 model = arx(data_est, best.orders); present(model); % 参数值 ± 标准差 figure; compare(data_val, model); % 验证段实测 vs 仿真 figure; resid(model, data_val); % 残差白性：自相关系数应落在置信带内典型运行结果：最优阶次收敛到[1 1 2]（与生成数据的真实结构一致），验证段 FIT 约 90% 以上，残差自相关系数基本落在 99% 置信带内——说明动态信息已被模型充分提取。若残差检验不过，应回去加大阶次，而不是直接宣布建模完成。
六、常见坑与优化清单
| 症状 | 可能原因 | 对策 |
|---|---|---|
polyfit警告 "Polynomial is badly conditioned" | x 量纲大、阶数高，X^TX 条件数爆炸 | 用[p,S,mu]=polyfit(...)中心化；或先把 x 归一化 |
| 系数巨大、正负交替、置信区间极宽 | 过拟合或共线性 | 降阶；验证集选阶；fitlm查 VIF，剔除/合并相关特征 |
| 训练 R²≈0.99 但换批数据就崩 | 过拟合 / 数据划分泄漏 | 严格划分训练/验证（时间序列按时间切，不能乱序抽样） |
| 拟合曲线在数据范围外发散 | 多项式外推失效（天性如此） | 只在标定范围内使用；需要外推时改用机理模型或渐近形式合理的函数（指数、有理式） |
regress结果与fitlm差一列 | regress需手动加全 1 截距列 | X = [ones(n,1), X]；或直接用fitlm |
| 公式写法报错 | 交互项写成x1:x2之外的形式 | 记住'y ~ x1*x2'= 主效应 + 交互项；分类变量自动虚拟化 |
| ARX 辨识结果 FIT 极低 | 输入激励不充分（u 恒定不变） | 实验设计阶段就要让输入充分变化（PRBS 伪随机信号是标准做法） |
| 残差检验不通过 | 阶次偏低 / 存在未建模动态 | 升阶；考虑armax（加噪声模型）或ssest |
| 训练/验证误差都很高 | 模型结构根本不对（欠拟合） | 回到残差图看结构；别急着加数据 |
七、小结
- 数据建模三阶梯：静态单变量用
polyfit/fit，多因素用fitlm/stepwiselm，动态系统用iddata+arx系工具； - 最小二乘的内核是正规方程 \hat\beta=(X^TX)^{-1}X^Ty，但工程代码请交给
\、polyfit（带mu）这类数值稳定的实现； - 验证集误差是选模型复杂度的唯一可靠标准，训练误差和 R² 都会骗人；
- 残差分析贯穿始终：静态看趋势，动态看白性，残差里有结构 = 模型里有遗漏；
- 数据质量决定上限：标定点要覆盖使用范围，辨识实验的输入激励要充分——模型永远超不过数据。
参考资料
- MathWorks.Curve Fitting Toolbox User's Guide. R2024b 版官方文档. https://www.mathworks.com/help/curvefit/
- MathWorks.System Identification Toolbox User's Guide. https://www.mathworks.com/help/ident/
- MathWorks.fitlm / stepwiselm 函数参考（Statistics and Machine Learning Toolbox）. https://www.mathworks.com/help/stats/fitlm.html
- Ljung, L.System Identification: Theory for the User(2nd Edition). Prentice Hall, 1999.（系统辨识领域的经典教材）
- Montgomery, D. C., Peck, E. A., Vining, G. G.Introduction to Linear Regression Analysis(5th Edition). Wiley, 2012.
- 姜启源, 谢金星, 叶俊. 《数学模型（第五版）》. 高等教育出版社, 2018.（中文建模入门经典）
