从零实现魔术公式轮胎模型：Matlab拟合与车辆动力学仿真应用

做车辆动力学仿真和底盘控制开发的朋友，对"魔术公式轮胎模型"这个词一定不陌生。这套由荷兰学者Pacejka提出的半经验轮胎模型，核心是一组三角函数，用几个参数就把轮胎的纵向力、侧向力和回正力矩表达成滑移率、侧偏角和垂向载荷的函数。在Matlab里把模型和代码实现打通，是车辆工程专业学生、底盘工程师和ADAS算法人员绕不开的一步。我最早接触它是在做EPS力矩标定项目的时候，手头只有一套台架试验数据，手册上的参数又不匹配，硬是花了两个多星期才把曲线拟合到工程可用的程度。这篇文章就把我从零实现魔术公式的经验完整梳理一遍，从公式原理、Matlab代码结构、参数拟合方法到常见坑，一次讲清楚。

这篇文章适合谁看？如果你正在做底盘控制、ADAS横向/纵向控制、整车操稳仿真，或者只是在准备课程大作业和数据拟合项目，下面的内容都能直接落地。我会把参数含义、函数实现、拟合流程和避坑经验全部拆开讲，保证你看完能动手写代码，也能理解自己在写什么。

1. 魔术公式到底在表达什么

1.1 一条轮胎曲线背后的数学逻辑

魔术公式轮胎模型的基本形式是一句话能写完的：

Y = D * sin(C * arctan(B * x - E * (B * x - arctan(B * x)))) + Sv

其中 x = X + Sh。X 是输入量（滑移率或者侧偏角），Y 是输出量（纵向力、侧向力或者回正力矩）。B、C、D、E、Sh、Sv 这六个字母就是"魔术"的来源，也是做参数拟合时要去标定的未知数。

先别急着跳公式，我用最直白的方式解释一下每个参数干了什么。D 是峰值因子，直接决定了曲线最高点的高度，你把D改了，整个输出的峰值就跟着变，它是最容易从数据里读出来的参数。C 是形状因子，控制曲线整体形态偏向正弦还是偏向一条带饱和的S形，范围通常在 1 到 2 之间。B 是刚度因子，它和 C、D 三者相乘得到的 BCD 正好等于曲线在原点附近的斜率，也就是说它控制的是小输入下力的增长快慢，这一条对整车模型的影响非常直接：转向刚开始转动时，车辆横摆响应快不快，基本就是它说了算。E 是曲率因子，它单独控制曲线峰值附近弯曲下降的幅度，让曲线从"圆滑的拱顶"变成"带尖角的峰值"。

Sh 和 Sv 则是一对偏移量，分别沿输入轴和输出轴平移整条曲线。为什么需要平移？因为真实轮胎在零侧偏角时，由于胎体不对称、残余内应力、滚动阻力和测量系统误差，输出力不一定是零，这时候靠 Sh 和 Sv 就能把曲线校准回试验数据。

从数学形状上看，这个公式可以理解成在 arctan 外层再套一层 arctan 再乘上 sin。arctan 本身是一条先近似线性、后逐渐饱和的曲线，两层组合之后，正好能表达轮胎在小滑移区域近似线性、在大滑移区域出现附着极限并缓慢回落的形态。如果用多项式拟合这种数据，往往在数据范围边缘出现不可控外推，而魔术公式由于 arctan 天然饱和，外推表现稳定很多，这也是它在工程中被广泛采用的重要原因。

1.2 三张曲线面孔：纵滑、侧偏与回正力矩

同一个魔术公式，在三种不同的受力场景下呈现三张面孔，只是参数不同，输入变量不同。

第一张是纵向力与滑移率的关系，也就是 Fx-kappa 曲线。滑移率是无量纲量，制动时取负值，驱动时取正值，典型仿真范围在 -0.3 到 0.3 之间。滑移率很小时，纵向力随滑移率近似线性增长；滑移率到 0.1 到 0.2 附近，轮胎进入附着饱和区，纵向力达到峰值；再继续增大滑移率，力反而会缓慢下降。这个形态在ABS开发中非常关键，因为控制器需要判断当前到底是处在线性区还是饱和区。

第二张是侧向力与侧偏角的关系，也就是 Fy-alpha 曲线。输入侧偏角一般用度数表示，但在计算函数里必须转成弧度。曲线在小角度（约 2 到 3 度以内）近似线性，之后逐渐饱和，到 10 度到 15 度左右达到极限。这条曲线决定了车辆稳态转向特性，也是ESP和主动转向控制的重要基础。

第三张是回正力矩与侧偏角的关系，Mz-alpha 曲线。回正力矩的形态比较特殊：随侧偏角先增大，然后很快到达峰值，再下降穿越零点变成负值，呈现"先扬后抑"的形态。别小看这条曲线，EPS系统的方向盘手感、自动回正控制、转向系统残余力矩补偿，全都依赖它。记得我做EPS标定时，最头疼的就是回正力矩模型参数不对，导致低速回正仿真结果总是差了半拍。

另外必须提一下载荷依赖问题。同一套轮胎参数在 2000N 和 8000N 的垂向载荷下是完全不同的曲线，峰值会随载荷整体上移。工程上的标准做法是取多个载荷工况分别标定一组参数，或者在模型里用参数化方程描述 D、B 等随 Fz 的变化，比如 D = a1 * Fz^2 + a2 * Fz。前者实现简单，在固定载荷仿真中完全够用；后者适合全工况仿真，也是商用车辆动力学软件（如CarSim、TruckSim）内部的做法。用Matlab做代码实现时，我最推荐先做多载荷点分别拟合，验证思路通了再考虑参数化，一上来就啃MF 5.2完整公式容易被一堆系数淹没。

2. 用Matlab实现一套可复用的计算函数

2.1 参数结构体的设计思路

写Matlab代码的第一步，是先把参数组织起来。我见过不少人把B、C、D、E这些参数散落在脚本里，画图时手动改数字，拟合时又复制一遍，改一处漏一处，最后对不上号。正确的做法是定义一个结构体，把同一套工况的所有参数集中放进去。

% 纵向力参数示例，对应某乘用车轮胎在Fz=4000N附近的标定结果 params_fx = struct( ... 'B', 11.3, ... 'C', 1.78, ... 'D', 5830, ... 'E', 0.98, ... 'Sh', 0.01, ... 'Sv', -35);

这样的一组参数意味着什么？在滑移率零点附近，初始斜率 K = BCD = 11.3 * 1.78 * 5830，大约 117000 N 每单位滑移率。也就是说，滑移率增加 0.01，纵向力约增加 1170N，这个量级对单条轮胎来说是合理的。D 是 5830N，说明这条轮胎在 4000N 垂向载荷下最多能提供约 5800N 的纵向力，附着系数已经接近 1.45，这是偏向高性能轮胎的数据。Sh 和 Sv 比较小，说明试验曲线基本过原点，装车后由于制造误差产生的偏移很小。

侧向力的参数自然是另一套，写法完全一样。这里要提醒一句：B 因子在使用侧偏角时，单位必须对应弧度制。如果一组侧向力参数的 B 是 6.9，那它的物理意义是每弧度对应 6.9 个单位的"等效曲率增长"，不是每度。很多初次接触的人在拟合时直接把角度传进去，结果得到的 B 比正常值大了近 57 倍，曲线形变到无法直视。

2.2 核心计算函数代码实现

有了参数结构体，计算函数可以写得很干净。我给出一套我常用的纯纵滑和纯侧偏计算函数，它们完全基于魔术公式原始形态，没有任何花哨扩展，适合理解、修改和后期扩展。

function Fx = magic_formula_fx(kappa, Fz, p) % kappa: 滑移率,无量纲,驱动为正,制动为负 % Fz: 垂向载荷,单位N % p: 参数结构体,包含 B C D E Sh Sv x = kappa + p.Sh; arg = p.B * x - p.E * (p.B * x - atan(p.B * x)); Fx = p.D * sin(p.C * atan(arg)) + p.Sv; end

function Fy = magic_formula_fy(alpha_deg, Fz, p) % alpha_deg: 侧偏角,单位deg % Fz: 垂向载荷,单位N % p: 参数结构体,包含 B C D E Sh Sv alpha = deg2rad(alpha_deg); % 统一转为弧度,与B因子的单位匹配 x = alpha + p.Sh; arg = p.B * x - p.E * (p.B * x - atan(p.B * x)); Fy = p.D * sin(p.C * atan(arg)) + p.Sv; end

你会发现两个函数的主体几乎一样，只是输入单位处理不同。这正是魔术公式的优点：同一套数学框架，换参数就是另一条曲线。函数里 Fz 参数当前只做占位，因为你已经针对固定载荷标定了 D、B 等参数；如果你后续要引入载荷参数化，只需在函数内部根据 Fz 去更新 p 里的 D 和 B，比如写一个 load_dependent_params(Fz) 函数返回更新后的参数结构体。

2.3 主脚本绘图与曲线合理性检查

函数写完之后，主脚本就是装配逻辑了。定义参数、生成输入范围、调用函数、画图，一气呵成。

% 画纵向力-滑移率曲线 kappa = linspace(-0.3, 0.3, 300); Fx = magic_formula_fx(kappa, 4000, params_fx); figure('Color', 'w'); plot(kappa, Fx, 'LineWidth', 2); xlabel('滑移率 \kappa'); ylabel('纵向力 Fx (N)'); grid on; title('纵向力与滑移率特性');

把这一段跑出来之后，先别急着继续写侧向力，花一分钟做三个主观检查：第一，曲线是否连续光滑，有没有突变或NaN；第二，曲线在零点附近是否单调上升，初始斜率是否和 BCD 的估算一致；第三，峰值出现在哪里，后面有没有回落趋势。这三个点的物理合理性，比任何代码注释都重要。如果曲线在滑移率 0.1 之前就一路冲到天上不回头，那大概率是参数单位写错了；如果曲线看起来像一条直线，根本不饱和，那可能是 D 设置得过大，把弧度制下的有效范围全压在了线性段。

3. 从试验数据到模型参数：拟合流程全记录

3.1 参数初值估算：先用眼睛读曲线

很多教程会直接摆一段 lsqnonlin 代码让你跑，但闭着眼睛调参数的结果，往往是残差下降了一点，输出曲线却长得和试验数据完全两码事。原因很简单：非线性最小二乘对初值非常敏感，魔术公式的参数里还有明显的相关性。所以我的习惯是先做一轮手动的初值估算，把每个参数的范围压到合理区间，再交给优化算法去精调。

初值估算最朴素的方法是"看图说话"。拿到试验数据后，按这几个步骤走：

1. 找峰值，D 的初值直接取曲线最大输出值，再减掉 Sv 的估计值。
2. 在零点附近画一条切线，读初始斜率 K。然后假定 C 在 1.3 到 1.8 之间取一个值，用 B = K / (C * D) 算出 B。
3. 看曲线尾部。如果峰值之后明显下探，E 取正值（通常在 0 到 1 之间）；如果曲线到达峰值后基本平着走，E 取接近 0 的值或小幅负值。
4. 看曲线的零点偏移。如果输入为 0 时输出不为 0，它的值就是 Sv 的初值；如果峰值在横轴上不在原点对称位置，说明需要 Sh 来平移。
5. 最后检查一下 B*x 的量级是否合理。如果 B 的初值算出来是 50 以上，别急着用，回去检查斜率 K 的量纲，很可能单位又混了。

这套"几何法"最妙的地方，是每一步都有明确的物理对应，不需要任何优化算法就能得到一个能跑通的初值。我通常用平滑后的数据（比如滑动平均）来做这一步，避免把试验噪声当成曲线特征。

3.2 lsqnonlin非线性最小二乘拟合实战

初值给定之后，就可以上 lsgnonlin 了。这是Matlab自带的非线性最小二乘求解器，专门干这种"给定模型和残差，找参数让残差平方和最小"的活。

% 试验数据（示意，请替换成您的台架数据） kappa_data = [-0.2 -0.1 -0.05 -0.02 0 0.02 0.05 0.1 0.2]; Fx_data = [-5200 -4800 -3900 -1800 0 1900 4100 4900 5100]; % 固定C,减少参数耦合 C_fixed = 1.8; % 待拟合参数: [B, D, E, Sh, Sv] p0 = [11, 5600, 0.9, 0, 0]; model = @(p, kappa) p(2) * sin(C_fixed * atan( ... p(1)*(kappa+p(4)) - p(3)*(p(1)*(kappa+p(4)) - atan(p(1)*(kappa+p(4)))))) + p(5); resid = @(p) model(p, kappa_data) - Fx_data; opts = optimoptions('lsqnonlin', 'Display', 'iter', ... 'MaxIterations', 500, ... 'FunctionTolerance', 1e-8, ... 'StepTolerance', 1e-8); lb = [0, 0, -5, -0.1, -500]; ub = [50, 10000, 5, 0.1, 500]; p_fit = lsqnonlin(resid, p0, lb, ub, opts);

为什么先固定 C？因为 B、C、D 三个参数在 BCD 这个乘积里高度耦合，如果三个一起自由变化，优化器会沿着"增加B减小C"这条等高线滑来滑去，收敛很慢甚至不收敛。先把 C 固定在一个合理值，让优化器集中火力去调 B 和 D，等曲线主体形态对了，再放开 C 做最终精修，会稳定得多。

每次拟合结束都要读一下优化器输出的残差范数，并画出拟合曲线和试验数据对比图。如果肉眼可见的偏差已经消失，再算一个R方，确认方差解释率在0.95以上。工程里过度追求小数点后四位没有意义，能抓住趋势就够用了。

3.3 多载荷点拟合的组织方式

实际项目里不可能只在单一载荷下做标定，整车在弯道制动时载荷转移非常普遍，所以至少要在 2000N、4000N、8000N 三个载荷点分别重复拟合流程。很多教程建议把所有载荷数据混在一起统一拟合，这个做法要慎重。纯魔术公式六参数里没有显式包含 Fz，把不同载荷的数据混合在一起，拟合器会在"高载荷高峰值"和"低载荷低峰值"之间取一个折中的 D，结果哪条曲线都拟合不好。

我更推荐先把每个载荷点看成独立的数据集，分别跑一遍上面说的"初值估算+lsqnonlin"流程，得到三组参数。然后再看 D、B 随 Fz 的变化趋势：如果 D 基本随 Fz 线性增长，就可以用二次多项式去拟合 D(Fz)、B(Fz)、E(Fz)，把结果写进参数化函数里，这样仿真时就能连续插值出任意载荷下的参数。这个分两个阶段的做法，比单一阶段"混沌拟合"好调试得多，也方便检查哪条曲线出了问题。

4. 常见问题与排查技巧实录

4.1 单位、符号与数值陷阱

算是最常见的三类坑，我按频率排序：

第一，侧偏角单位混用。函数里写的是角度转弧度，但拟合时如果直接拿角度值做输入，B 因子会被严重低估。检查方法很简单：看拟合出的 B，如果侧偏角工况下 B 的数量级在 0.1 到 1 之间，基本是弧度制；如果在几到十几之间，大概率是角度制混进去了。

第二，滑移率正负号约定不一致。AGV、车辆动力学仿真软件、试验台架对驱动/制动的符号约定可能有差异，纯代码层面没有对错，但一旦前后不一致，拟合出的 Sh 会凭空出现很大的偏移，掩盖真实物理偏移。建议在代码注释里明确写上"驱动为正、制动为负"或相反，保持全局统一。

第三，极限输入下的数值风险。虽然 arctan 天然饱和，但当你输入的 kappa 超出 ±1，Bx 可能达到上百，atan(Bx) 趋近 π/2，这没问题。问题是 if E 取得太极端（大于 10），括号里的 Bx - arctan(Bx) 会变成很大的数，再乘以 E，可能导致三重括号里的值异常波动，进而输出振荡。处理办法是限制 E 的范围为 [-5, 5]，并在画图时扫一下输入范围边界，确认无异常。

4.2 拟合不收敛与参数越界

拟合不收敛和越界的本质原因是一样的：参数空间太大，而目标函数表面太平坦。具体到魔术公式，B、C、D 之间强相关，会导致优化器在"陡峭山谷"里反复振荡，迟迟不收敛。

解决思路分三步。第一步，固定 C，先拟合 B、D，这招能把收敛难度降低一大半。第二步，把数据范围截取到足够覆盖饱和区，至少包含峰值前后各几个点，否则 E 和 Sh 根本被激发不出来，优化器自然瞎跑。第三步，给上下界约束，上面代码里的 lb、ub 不是摆设，它们把搜索空间限制在物理可解释的范围内，既防止越界，也加速收敛。

如果还是收敛失败，我建议用全局优化工具箱的 particleswarm 或 MultiStart 配合 lsqnonlin 做多起点搜索。先用粒子群大致扫一遍参数空间，再用扫到的点作为 lsqnonlin 初值精修，实测在难度较高的回正力矩拟合中能稳定收敛。

4.3 曲线形态异常检查清单

经验多了之后，几乎不用看数字，扫一眼曲线形状就能定位问题。下面这份清单是我总结的速查表，建议存一份：

这里特别说一句回正力矩曲线的拟合。它先升后降再穿越零点，这种非单调形态对初值要求极高，我的经验是先固定 C 和 D，手动把 E 从 0 慢慢增大，观察曲线尾部下探趋势是否和试验一致，然后才交给优化器做精细调整。这一步一旦偷懒，优化器很容易把曲线拟合出一条"一直在下降"的形态，残差还不小，让人一头雾水。

最后再分享一个非常实用的经验：写一个交互式可视化函数，把参数结构体和数据都传进去，快速画出拟合曲线，再用Matlab的实时脚本或者 App Designer 做个滑块控件，直接拖动 B、C、D、E 观察曲线变化。这样调参虽然粗糙，但比盲跑优化器高效得多，能在十秒内锁定一个接近最优的初值范围。我做车辆动力学预研项目时，基本都是靠这套"滑块预调+lsqnonlin精修"的组合拳来处理标定数据，效果稳定，调试效率也高。如果你也在跟魔术公式模型较劲，不妨先把这个工具搭起来，再回头碰数据拟合，会顺很多。