魔术公式轮胎模型MATLAB实现与参数辨识指南

做整车动力学仿真的朋友，迟早会和魔术公式轮胎模型打交道。它是一套用三角函数表达轮胎纵向力、侧向力与回正力矩的半经验模型，配合MATLAB代码实现，是底盘、悬架和稳定性控制仿真里最常见的轮胎表达方式。这篇文章不打算从教科书推导开始，也不直接丢一个完整工程包出来，而是想把我自己从“套公式”到“能落地”这个过程中觉得最关键的东西拆开讲：公式结构、参数含义、代码写法，以及用真实数据做参数辨识时那些文档里不会写的坑。

如果你只是需要一条“看着像轮胎”的曲线，那网上现成的代码已经很多；但如果你想用魔术公式做参数研究、把轮胎模型接进整车模型，或者想把试验台架数据拟合成可用参数，这篇会比较对路。我会尽量用直白的话，配合可以直接跑的MATLAB示例，把每一步为什么这么做交代清楚。

1. 轮胎模型在整车仿真里的位置：为什么非要魔术公式

1.1 整车仿真的瓶颈常常是轮胎

做过整车操纵稳定性仿真的人可能都有这种感受：悬架模型、转向模型、车身模型都有比较成熟的建模套路，但一到轮胎，问题就来了。轮胎不是简单弹簧，也不是线性阻尼，它在大侧偏角、大滑移率下表现出强烈的非线性，而且不同载荷、不同路面下特性完全不一样。可以说，整车里最“不线性”的部件就是轮胎。

底盘控制类项目里，比如ABS、TCS、ESC的算法验证，最终作用对象就是轮胎力。轮胎模型给得不准，控制器算法做得再精细，仿真结论也没有说服力。所以在项目前期，我一般会花不少时间在轮胎模型上，先把轮胎这条曲线“立住”，再谈后面的整车验证。这个投入很值得，因为轮胎模型一旦偏差大，后面所有东西都会被带偏。

1.2 纯理论模型vs半经验模型：工程上的取舍

轮胎模型大体可以分两类：一类是纯物理模型，从橡胶材料、胎体结构、接地印迹应力分布去推导；另一类是半经验模型，用特定数学函数去逼近试验数据。纯物理模型精度高、外推能力强，但参数多，很多参数在工程阶段根本拿不到，而且求解成本高，不适合做实时的整车控制仿真。

魔术公式属于半经验模型，它不追求从物理机理解释每一个力是怎么产生的，而是用一组带明确形状控制的三角函数，把轮胎的稳态力特性拟合出来。优点很明显：参数数量少、计算速度快、拟合精度足够、数值上也稳定。所以它在工程仿真里特别流行。代价就是你必须有可靠的试验数据去做参数辨识，否则模型只是“看起来像”。

2. 魔术公式的数学骨架：B、C、D、E到底在控制什么

2.1 一条公式看懂魔术公式

魔术公式的基础形式并不复杂，核心就是下面这个函数：

y = D sin(C arctan(Bx - E(Bx - arctan(Bx)))) + S_v

其中，x = X + S_h。X 是自变量，对纵向力来说通常是滑移率，对侧向力来说通常是侧偏角。D、C、B、E分别负责控制曲线的峰值、形状、刚度和曲率；S_h 和 S_v 是水平偏移和垂直偏移，用来处理试验数据里曲线不过原点的情况。

这个公式初看有点绕，但拆开就很简单：先让自变量 x 乘上刚度因子 B，然后用 arctan 把大输入压缩到有限区间，再用 C 控制形状，最后用 D 控制峰值，E 调节过渡段的弯曲程度。整条曲线天然被限幅在 D 附近，不会出现力无限增长的情况，这符合轮胎力最终饱和的物理特性。

2.2 四个系数的物理映射

我把这几个系数的物理含义整理成了一张表，方便对照理解：

注意，这个公式里的角度自变量必须用弧度，或者把 B 做对应缩放。我在早期犯过把角度直接按度数代进去的错，结果 B 值怎么拟合都不对。后来统一成弧度之后，参数辨识稳定多了。这是很小的细节，但影响很大。

2.3 载荷对参数的影响：先别急着写死

魔术公式的另一个关键点是：参数并不是常数，而是垂直载荷 F_z 的函数。轮胎受到的垂直载荷越大，侧向力峰值 D 通常也越大，刚度因子 B 也会跟着变化。如果只测了一个载荷下的数据，那拟合出来的参数只能用于那个载荷，仿真时一旦载荷变化，轮胎力就会明显失真。

所以工程上常见的做法是做多组不同载荷下的试验，分别拟合 B、C、D、E、S_h、S_v，再把这些参数做成随载荷变化的插值表。如果数据量不够，也可以对 D 做线性化，让它随 F_z 近似线性增长，其余参数只保留若干离散点。实现上没有统一标准，但原则很明确：不要把一个载荷下拟合出来的参数当成全局模型用。

3. 动手写MATLAB之前，先把这三件事想清楚

3.1 数据从哪来：试验、仿真还是文献

写代码前，我建议先确定参数来源。最常见的是轮胎转鼓试验台或整车道路试验的数据；没有试验条件时，也可以用成熟的车辆动力学软件先导出一组参考曲线，再用魔术公式去拟合。这个思路在项目预研阶段特别常用，虽然数据精度不如真实试验，但可以先把模型框架跑通。

如果是研究性质，没有试验数据，那也可以用文献里常见轮胎参数先占位。我不反对直接用经典参数，但心里要有数：那是某个特定尺寸轮胎的参数，不一定代表你关心的轮胎。在用MATLAB跑通代码之后，第一件事就是检查模型在常用滑移率范围内的稳态力是否符合物理直觉，比如侧向力应该随着侧偏角增大先线性上升后饱和，纵向力应该随滑移率出现峰值后回落，这些都是基本判断。

3.2 单位统一，尤其是角度

单位问题看起来基础，但很多人栽跟头。我见过有人在同一个脚本里，前一个模块用弧度，后一个模块用角度，拟合出来的 B 值互相矛盾的案例。魔术公式里，滑移率是无量纲的，一般 0 表示纯滚动，1 表示完全抱死；而侧偏角的单位要统一为弧度。

如果你从试验数据里拿到的侧偏角是度数，那导入MATLAB时先除以180再乘以pi，或者用 deg2rad 函数转一下。拟合完成之后，验证阶段也要保持同样的单位习惯。否则同样一组参数，画出来的曲线形态会差很多，而且你会误以为模型不准。

3.3 确定输入输出方向：纯工况还是联合工况

魔术公式的经典形式是针对纯纵向滑移和纯侧偏工况的。所谓纯工况，就是要么只制动/驱动，要么只转向，不把两者混在一起。实际车辆行驶时往往是纵向力和侧向力同时存在，比如边制动边转向，这时候需要联合工况模型。

联合工况的建模逻辑通常是先分别算纯工况下的纵向力、侧向力和回正力矩，再用相似理论或权重组合进行耦合。标题里的“魔术公式轮胎模型研究”如果定位在基础研究，我建议先把纯工况做扎实，再接联合修正。不要一上来就想一个函数搞定所有工况，那样参数辨识会和真实数据打架，很难收敛。

4. MATLAB实现：从公式到可复用的函数

4.1 函数接口设计

写MATLAB代码时，我习惯把核心公式封装成一个函数，输入是自变量、垂直载荷和参数结构体，输出是力。这样后面做参数拟合、批量扫描、整车集成都方便。如果直接把公式散写在脚本里，参数一改就得改多处，非常容易出错。

接口上，我推荐用结构体传参数，而不是六个零散输入。结构体字段名直接叫 B、C、D、E、Sh、Sv，可读性好，也方便和拟合结果对接。垂直载荷单独传参，是为了后续扩展成载荷插值表。

4.2 一个能跑的纵向力示例

下面这个代码是纵向力纯工况的最小可用版本，核心就是魔术公式标准形式。为了简化，这个示例先假设参数在某个固定载荷下已经拟合好：

function Fx = mf_fx(kappa, Fz, p) % mf_fx 魔术公式纵向力计算（纯滑移工况） % 输入: % kappa - 滑移率，无量纲，0表示自由滚动 % Fz - 垂直载荷，单位N % p - 参数结构体，包含 B C D E Sh Sv % 输出: % Fx - 纵向力，单位N x = kappa + p.Sh; Bx = p.B * x; Fx = p.D * sin(p.C * atan(Bx - p.E * (Bx - atan(Bx)))) + p.Sv; end

这个函数不长，但已经把模型核心表达清楚了。调用方式也很简单：

p.B = 10; p.C = 1.65; p.D = 3200; p.E = 0.4; p.Sh = 0; p.Sv = 0; kappa = linspace(-1, 1, 200); Fz = 4000; Fx = mf_fx(kappa, Fz, p); plot(kappa, Fx); xlabel('滑移率'); ylabel('纵向力/N'); grid on;

跑完之后你会看到一条从负到正、先增长后饱和的曲线。如果你手里的参数不是这组值，曲线形状可能不同，但整体趋势应该是符合轮胎特性的。

4.3 扩展成侧向力与载荷插值

侧向力的实现逻辑完全一样，只是把自变量从滑移率换成侧偏角 alpha，参数换成侧向力参数。要注意的是，侧偏角很小的时候，侧向力几乎线性增长，这个线性段的斜率就是“侧偏刚度”，很多人喜欢直接用线性模型近似。魔术公式的优势在于，它能表达出大侧偏角下的饱和特性，所以在极限工况仿真里不能丢掉。

载荷插值也很容易在MATLAB里实现。你可以在一个数据结构里保存多个载荷点下的拟合参数，然后用 interp1 对参数做线性插值：

function [Fy] = mf_fy_load(alpha, Fz, paramTable) % paramTable.FzList: 载荷点数组 % paramTable.P: 结构体数组，每个元素对应一个载荷点的参数 idx = find(paramTable.FzList <= Fz, 1, 'last'); if isempty(idx), idx = 1; end if idx >= length(paramTable.FzList), idx = length(paramTable.FzList) - 1; end Fz0 = paramTable.FzList(idx); Fz1 = paramTable.FzList(idx + 1); w = (Fz - Fz0) / (Fz1 - Fz0); p0 = paramTable.P(idx); p1 = paramTable.P(idx + 1); p.B = p0.B + w * (p1.B - p0.B); p.C = p0.C + w * (p1.C - p0.C); p.D = p0.D + w * (p1.D - p0.D); p.E = p0.E + w * (p1.E - p0.E); p.Sh = p0.Sh + w * (p1.Sh - p0.Sh); p.Sv = p0.Sv + w * (p1.Sv - p0.Sv); Fy = mf_fy(alpha, Fz, p); end

这个函数只是线性插值的示例，真正的工程代码还要加边界检查、载荷超出范围时的外推策略。但思路就是：载荷变了，参数跟着变，而不是换一套完全不同的公式。

5. 参数辨识：让仿真曲线贴合真实轮胎数据

5.1 最小二乘拟合的目标函数

拿到真实试验数据之后，我们需要让模型输出尽量贴近测量值。这在MATLAB里通常用 lsqcurvefit 来完成。目标函数不需要自己手动写残差平方和，只要把模型函数传进去，优化算法会自动调整参数，让模型曲线和试验数据差最小。

纵向力的拟合代码大概长这样：

% data_x: 试验滑移率 % data_Fx: 试验纵向力 modelfun = @(p, x) mf_fx_fit(p, x, Fz_fixed); p0 = [10, 1.65, 3000, 0.3, 0, 0]; lb = [0, 0.5, 0, -1, -0.1, -100]; ub = [50, 3, 10000, 2, 0.1, 100]; options = optimoptions('lsqcurvefit', 'Display', 'iter', 'MaxFunctionEvaluations', 5000); p_fit = lsqcurvefit(modelfun, p0, data_x, data_Fx, lb, ub, options);

这里的 mf_fx_fit 是固定垂直载荷下的核心公式，输入是一维参数向量 p，内部把它映射到 B、C、D、E、Sh、Sv。不要把结构体直接丢给 lsqcurvefit，优化器只认数值向量。

5.2 初值、边界和归一化的经验

参数辨识最容易翻车的地方就是初值和边界。单独用 lsqcurvefit 并不能保证找到全局最优，初值给得离谱，结果可能就掉进局部极小值。

我的习惯是先肉眼估计几个关键量：D 的初值看试验曲线峰值，直接取峰值力附近；B 的初值看原点附近斜率除以 D，再适当缩小；C 的初值侧向力取1.3左右，纵向力取1.6左右；E 的初值取0.5以内。边界设定上，C 一般限制在0.5到2.5之间，E 限制在-1到2之间，D 不能为负。这些边界不需要太紧，但能避免优化器跑到明显不合理区域。

另外，如果试验数据和模型在数值量级上差距太大，比如力是几千牛，滑移率是0到0.3，可以考虑先做归一化。把力除以峰值，把自变量做缩放，这样优化过程更稳定。我之前拟合高速轮胎数据时，直接拿原始数值跑，优化器总是说迭代不收敛；后来把数据归一化到0到1附近，再拟合，很快就收敛了。拟合完再反归一化回物理量，效果完全一致。

5.3 拟合结果怎么评价

拟合完成别只看 R² 或者均方根误差一个指标。我一般会画出模型曲线和试验数据点叠在一起的图，并重点看三个区域：原点附近、峰值附近、大滑移率回落段。原点附近决定刚度能不能对上，峰值附近决定最大附着能力的表达，回落段决定制动工况下的曲线形状。

如果原点附近拟合不好，往往说明 B 或者 S_h 不对；如果峰值附近差，多半是 D 和 E 的耦合关系没调好。还要注意有没有过拟合单条曲线。比如你只拟合了一个载荷下的小数据集，可能会得到一条完美穿过所有点的曲线，但换个载荷就完全不能用。所以做参数辨识时，一定要预留一部分工况数据不参与拟合，用独立数据做验证。

6. MATLAB实现与参数辨识中容易踩的坑

6.1 稳态模型解决不了瞬态问题

魔术公式本质上是稳态轮胎模型，它描述的是轮胎在固定滑移率、固定侧偏角下的力。但真实车辆转向、制动时，轮胎力是有建立过程的，也就是所谓的松弛长度效应。如果你把魔术公式直接接入整车模型，并且模型步长很小，你会看到轮胎力跳变瞬间跟上输入变化，这其实不符合物理。

解决思路有两个：一种是在模型前端加一阶惯性环节，对侧偏角或滑移率做滤波，让轮胎状态逐步变化；另一种是直接使用带瞬态特性的轮胎模型。后面这种复杂度会高不少。对大多数底盘控制仿真的前期研究来说，用一阶惯性处理已经够了。这也是魔术公式研究中容易被忽略的一个点，因为纯公式看起来没有时间项。

6.2 参数相关性与多组工况的矛盾

参数辨识时，B、E、C 之间其实存在很强的相关性。很多时候你会发现，B 调大一点、E 跟着调大一点，得到的曲线差别不大。这不是优化器的问题，而是模型本身的参数冗余。出现这种情况时，我一般会固定其中一个参数，只拟合其余几个，或者采用分步拟合：先用低频段粗略确定 B，再固定 B 去拟合 E 和 C。

多组工况拟合时还有个麻烦：不同载荷下拟合出的参数可能并不平滑。这可能是因为轮胎本身在不同载荷下特性有突变，也可能是试验数据里有噪声。处理办法是先看原始试验数据，确认是否每条曲线的趋势都一致；然后把参数随载荷的变化做成平滑插值，而不是直接用分散点。不要强行让所有载荷下的参数落在一条光滑曲线上，那不一定符合真实轮胎。

6.3 我的参数管理习惯

最后分享一点个人经验。我在做这类模型研究时，会把参数管理的核心放在“可复现”上。每个载荷点的拟合参数、拟合用的数据文件、以及拟合的初值边界设置，都会存成独立的MATLAB脚本或者结构体文件，而不是只保留最终结果。因为过了两周你再回来看代码，如果只有一条最终曲线，根本不知道它对应什么工况、什么轮胎、什么试验条件。

另外，所有拟合参数我都会写上一个简短的来源备注。这个备注不是给机器看的，而是给未来的自己看的。比如“来自某预研项目的高速工况数据”或者“来自某试验台架的干路面数据”这种简短说明，能省去很多重新摸索的时间。魔术公式本身不复杂，复杂的是它周围那一堆数据和参数的真实性管理，而这恰恰决定了你的仿真结果到底有没有参考价值。