企业官网建设流程全解析

做车辆稳定性分析的朋友，应该都绕不开“质心侧偏角-横摆角速度相平面”这张图。我记得第一次用 MATLAB 画出带鞍点和临界轨迹的二自由度车辆模型相平面时，才真正理解了什么叫“临界失稳”。这个项目不是什么论文级高深算法，而是一套非常经典的稳定性分析流程：把二自由度车辆模型写成状态方程，在相平面上找平衡点，识别鞍点，再追踪鞍点稳定流形画出临界轨迹，从而判断质心侧偏角和横摆角速度在什么组合下车辆会失控。它解决的问题很直接——车辆稳定边界在哪里，控制策略应该把状态约束在哪个区域内。适合正在做 ESP、ESC、四轮转向或底盘控制算法验证的工程师，也适合刚入门车辆动力学建模的研究生照着复现。

1. 项目背景与核心概念解读

1.1 二自由度车辆模型到底在描述什么

二自由度车辆模型，业内也常叫“自行车模型”，不是说车上只坐两个人，而是把车辆简化成前后两个轮胎侧向力作用点，车身只有一个横向平移自由度和一个横摆转动自由度。它忽略了悬架运动、侧倾、俯仰、纵向车速变化这些复杂因素，默认纵向车速恒定，只关心车辆在转向输入下的侧向运动和绕垂直轴的转动。

为什么这么极端简化还能用？因为车辆横向稳定性分析的核心矛盾，其实是轮胎侧向力与车身惯性力的平衡。只要把前后轴侧偏刚度、轴距、质心位置、质量、转动惯量这些关键参数保留住了，质心侧偏角β和横摆角速度ω的动态特性就能被刻画得足够准确。ESP之类的稳定性控制器做状态反馈、边界判断，用的也就是这两个状态量。

这个模型虽然简单，但它在做“相平面分析”时有个天然优势：系统只有两个状态，相平面是二维的，所有轨迹都能画在平面图上，平衡点、鞍点、极限环这些概念都可以直观地摆出来。换成三自由度甚至更高维度的模型，相平面就只能投影了，判断稳定边界就没这么直观。

1.2 相平面、鞍点、临界轨迹在说哪三件事

先说相平面。以质心侧偏角β为横轴、横摆角速度ω为纵轴，每一组状态(β, ω)都对应平面上的一个点。给定一个初始状态，车辆的运动方程会驱动这个点在平面上走出一条曲线，这就是相轨迹。整张相平面图就是无数条相轨迹组成的“流场”，能一眼看出不同初始状态下系统是收敛、发散还是绕圈。

再说鞍点。鞍点是系统的平衡点，也就是状态变化率为零的位置。车辆动力学模型通常不是线性的，尤其在轮胎力进入饱和区后，会出现多个平衡点。在这些平衡点里，有一种特殊的鞍点：在该点附近，相轨迹一边被吸引、一边被排斥，像马鞍一样。对车辆来说，鞍点往往对应着临界稳定的失控状态，两侧的相轨迹朝着截然不同的方向分流，一边回归稳定，一边滑向发散。

最后说临界轨迹。通俗点讲，临界轨迹就是鞍点的稳定流形，它是一条（或一组）特殊的相轨迹，把相平面分成两个区域：一侧是状态能回归稳定的“安全区”，另一侧是最终发散的“危险区”。很多控制策略的核心，就是保证系统状态不越过这条临界轨迹。所以把三者画在一张图上，等于同时回答了“哪里稳定、哪里临界、哪里必死”三个问题。

2. 模型搭建与仿真参数准备

2.1 运动微分方程推导与符号约定

二自由度车辆模型的微分方程，国内教材和 ISO 标准里的符号略有差异，但物理本质一样。我习惯用下面这套记号：

• m：整车质量，单位 kg
• Iz：绕质心竖直轴的转动惯量，单位 kg·m²
• u：纵向车速，单位 m/s
• δ：前轮转角，单位 rad
• a：质心到前轴距离，单位 m
• b：质心到后轴距离，单位 m
• Cf、Cr：前后轮侧偏刚度，单位 N/rad
• β：质心侧偏角，单位 rad
• ω：横摆角速度，单位 rad/s
• Fyf、Fyr：前后轴侧向力，单位 N

基于“纵向车速u恒定”的假设，侧向运动方程和横摆力矩方程可以写成：

m·u·(dβ/dt + ω) = Fyf·cosδ + Fyr

Iz·(dω/dt) = a·Fyf·cosδ - b·Fyr

这里的前后轮胎侧偏角为：

αf = δ - (β + a·ω/u)

αr = - (β - b·ω/u)

注意，这套式子默认侧偏角符号遵循SAE标准：侧偏角为正时，侧向力为负。如果你用ISO坐标系，符号方向要统一检查，否则画出来的鞍点位置会左右翻转。

如果进一步假设前后轮轮胎力处于线性区，即 Fyf = -Cf·αf，Fyr = -Cr·αr，那就得到线性二自由度模型。但线性模型在稳定性分析里有个致命局限：它没法反映轮胎饱和，系统最多只有一个平衡点，也不容易出现真正的鞍点结构。因此我在项目里用的是非线性轮胎模型，这样相平面才“有看头”，临界轨迹也才会出现。

2.2 参数选取与无量纲化

示例参数我直接用一台中级轿车的典型值，方便你复现时对照：

轮胎峰值侧向力需要结合垂直载荷和路面附着系数估算。前轴静态载荷约为 Fz_f = m·g·b/(a+b)，后轴约为 Fz_r = m·g·a/(a+b)。如果附着系数 μ=0.8，那么前轴峰值侧向力大约是 Fyf_peak = μ·Fz_f，后轴类似。

轮胎非线性部分我建议用双曲正切饱和函数近似，形式为：

Fy = -Fy_peak · tanh(C·α / Fy_peak)

这个形式的优点有两个。第一，小侧偏角时 tanh(x)≈x，退化成线性关系，和线性轮胎模型兼容；第二，大侧偏角时 Fy 趋近于 ±Fy_peak，自然限制在附着极限内，不会像线性模型那样无限增长。比起魔术公式，它没有那么多拟合系数，做平衡点求解和特征值分析时对初值不敏感，适合快速验证。

2.3 把状态方程写进 MATLAB

把上述公式转成 MATLAB 函数，我建议用结构体 p 存放车辆参数，这样后续做参数扫描时不用反复改函数签名。代码块示例：

function dz = vehicleDynamics(t, z, p) beta = z(1); r = z(2); % 默认前轮转角为零，便于分析无输入时的稳定性 delta = p.delta; % 前后轮侧偏角 alpha_f = delta - (beta + p.a * r / p.u); alpha_r = - (beta - p.b * r / p.u); % 非线性侧向力 Fyf = -p.Fyf_peak * tanh(p.Cf * alpha_f / p.Fyf_peak); Fyr = -p.Fyr_peak * tanh(p.Cr * alpha_r / p.Fyr_peak); % 状态方程 dz = zeros(2, 1); dz(1) = -r + (Fyf * cos(delta) + Fyr) / (p.m * p.u); dz(2) = (p.a * Fyf * cos(delta) - p.b * Fyr) / p.Iz; end

主脚本里给 p 赋初值：

p.m = 1500; p.Iz = 2500; p.a = 1.2; p.b = 1.6; p.Cf = 80000; p.Cr = 100000; p.u = 20; p.delta = 0; g = 9.81; Fz_f = p.m * g * p.b / (p.a + p.b); Fz_r = p.m * g * p.a / (p.a + p.b); p.mu = 0.8; p.Fyf_peak = p.mu * Fz_f; p.Fyr_peak = p.mu * Fz_r;

这里我故意把 p.delta 设成 0，因为项目标题讨论的是“临界状态下的稳定性”，也就是没有主动转向干预时车辆在什么状态下会失稳。如果你要看稳态圆周工况，把 delta 设成常数，平衡点会平移，鞍点位置也会变，但分析方法完全一致。

3. 相平面绘制与鞍点求解实操

3.1 相轨迹的数值积分方法

绘制相轨迹最直接的方法是在 β-ω 平面上铺一层初始状态网格，然后用数值积分工具求解每个初始点到未来一段时间的轨迹。MATLAB 里首选 ode45，因为这种车辆动力学方程一般不算刚性问题，四阶五阶变步长求解器足够用了。

网格密度要控制好。太稀看不到流场结构，太密图会变成一团黑线。我的做法是先用大步长铺出全局趋势，再在鞍点附近加密局部轨迹，这样既能看到完整相图，又不会让关键区域的细节被淹没。

一个常见的误区是：把所有相轨迹画出后，直接看“是否收敛到原点”来判断稳定域。这虽然直观，但对鞍点附近的轨迹很敏感，数值误差稍微大一点，临界状态可能被判成稳定或发散。所以相轨迹只能作为辅助观察，严格判定还得靠平衡点和流形分析。

绘制方向场并叠加部分相轨迹的代码如下：

beta_range = -0.25:0.01:0.25; omega_range = -0.7:0.02:0.7; [BETA, OMEGA] = meshgrid(beta_range, omega_range); U = zeros(size(BETA)); V = zeros(size(BETA)); for i = 1:numel(BETA) dz = vehicleDynamics(0, [BETA(i); OMEGA(i)], p); U(i) = dz(1); V(i) = dz(2); end % 归一化方向场，避免箭头长短不一影响视觉效果 quiver(BETA, OMEGA, ... U ./ sqrt(U.^2 + V.^2), ... V ./ sqrt(U.^2 + V.^2), 0.5, 'Color', [0.7 0.7 0.7]); hold on; % 绘制典型相轨迹 for beta0 = -0.2:0.02:0.2 for omega0 = -0.6:0.1:0.6 [~, X] = ode45(@(t,z) vehicleDynamics(t,z,p), [0 8], [beta0; omega0]); plot(X(:,1), X(:,2), 'b-', 'LineWidth', 0.4); end end

注意，正方向积分时间 tspan 不要给太长，否则大量轨迹发散到图外，满屏都是飞出去的线，反而看不出局部结构。一般给 5~10 秒就够判断趋势了。

3.2 鞍点位置的两种求解路径

找鞍点本质上是找非线性方程组的零点，再判断零点处的雅可比矩阵特征值。路径一：用 Symbolic Math Toolbox 做符号推导；路径二：直接用 fsolve 做数值搜索。我实际项目里用 fsolve 更多，因为车辆模型符号展开后很冗长，数值法更快。

方程组的平衡点条件是 dz = zeros(2,1)。fsolve 需要给初值，不同初值会收敛到不同平衡点，所以要配合多初值扫描：

fun = @(z) vehicleDynamics(0, z, p); eqs = zeros(0, 2); for beta0 = -0.3:0.05:0.3 for omega0 = -0.8:0.1:0.8 opt = optimoptions('fsolve', 'Display', 'off', ... 'FunctionTolerance', 1e-10, 'OptimalityTolerance', 1e-10); [eq_tmp, fval] = fsolve(fun, [beta0; omega0], opt); if norm(fval) < 1e-8 eqs = [eqs; eq_tmp.']; end end end % 合并重复平衡点 eqs = uniquetol(eqs, 1e-6, 'ByRows', true);

平衡点求出来后，对每个点计算雅可比矩阵。手推导太麻烦，我习惯用中心差分数值差分：

function J = numericalJacobian(fun, x) n = length(x); h = 1e-6; J = zeros(n, n); for i = 1:n xp = x; xm = x; xp(i) = xp(i) + h; xm(i) = xm(i) - h; J(:, i) = (fun(xp) - fun(xm)) / (2*h); end end

然后分别求特征值，若实部符号相反，该平衡点就是鞍点：

J = numericalJacobian(fun, eq_tmp); eigVals = eig(J); if prod(real(eigVals)) < 0 disp('该平衡点为鞍点'); end

3.3 完整绘制代码与效果解读

我把方向场、相轨迹、平衡点和鞍点画到一起，并额外把鞍点用红色星号标出来。下面这段代码可以直接跑：

% 1. 方向场与相轨迹 figure; hold on; beta_range = -0.25:0.01:0.25; omega_range = -0.7:0.02:0.7; [BETA, OMEGA] = meshgrid(beta_range, omega_range); U = zeros(size(BETA)); V = zeros(size(BETA)); for i = 1:numel(BETA) dz = vehicleDynamics(0, [BETA(i); OMEGA(i)], p); U(i) = dz(1); V(i) = dz(2); end quiver(BETA, OMEGA, U./sqrt(U.^2+V.^2), V./sqrt(V.^2+V.^2), 0.5, 'Color', [0.75 0.75 0.75]); % 2. 搜索平衡点并画符号 opts = optimoptions('fsolve', 'Display', 'off', 'FunctionTolerance', 1e-12); for beta0 = -0.25:0.05:0.25 for omega0 = -0.7:0.1:0.7 try eq_tmp = fsolve(fun, [beta0; omega0], opts); J_tmp = numericalJacobian(fun, eq_tmp); if norm(fun(eq_tmp)) < 1e-8 J_tmp = numericalJacobian(fun, eq_tmp); ev = eig(J_tmp); if all(real(ev) < 0) plot(eq_tmp(1), eq_tmp(2), 'ko', 'MarkerFaceColor', 'g', 'DisplayName', '稳定平衡点'); elseif prod(real(ev)) < 0 plot(eq_tmp(1), eq_tmp(2), 'kp', 'MarkerFaceColor', 'r', 'MarkerSize', 12, 'DisplayName', '鞍点'); else plot(eq_tmp(1), eq_tmp(2), 'k^', 'MarkerFaceColor', 'y', 'DisplayName', '不稳定平衡点'); end end catch continue; end end end

跑出来的图里，你会看到原点附近通常是一个稳定平衡点，代表车辆本身具备在无输入下恢复稳定的能力。远离原点的地方会出现一对鞍点，左右各一个，对称分布。鞍点附近的方向场非常有意思：箭头朝里的是稳定方向，朝外的是不稳定方向。从稳定方向延伸出去的轨迹，就是后面要画的临界轨迹。

4. 临界轨迹绘制与稳定域分析

4.1 稳定流形与临界轨迹的关系

临界轨迹本质上是鞍点的稳定流形，在二维非线性系统里，它由从鞍点出发的两条“臂”组成。稳定流形上的点随着时间推进会被鞍点“吸向”鞍点，但一旦偏离那么一点点，系统就会沿着不稳定方向飞出去。因此在相平面上，稳定流形就像一个分水岭，把收敛到稳定平衡点的区域和发散到无穷远的区域隔开。

要绘制稳定流形，最标准的数值方法是：首先求鞍点处的雅可比矩阵，找到在鞍点附近把轨迹“指向鞍点”的特征方向（实部为负的特征值对应的特征向量）。然后在该方向上给一个小扰动作为初始点，对原系统做负方向时间积分。因为负时间下，“指向鞍点”的稳定方向会变成从鞍点向外延伸的方向，所以你就能沿着整条稳定流形把它“拓”出来。

我写了个便于理解的小示例，假设已经找到鞍点坐标 eq_saddle 和稳定特征向量 v_s：

eps = 1e-5; Tmax = 20; for direction = [1, -1] x0 = eq_saddle + direction * eps * v_s; [~, X] = ode45(@(t,z) vehicleDynamics(t,z,p), [0 -Tmax], x0); plot(X(:,1), X(:,2), 'r-', 'LineWidth', 2, 'DisplayName', '临界轨迹'); end

画出来后会看到两条红色轨迹，从鞍点附近向两侧延伸，最终把整个相平面分割开。注意这里初始扰动 eps 不能太大，否则初始点已经偏离流形；也不能太小，否则会被数值精度吃掉。实际调试时通常取 1e-5 到 1e-3 之间。

4.2 稳定域的判定与可视化

临界轨迹画出来后，稳定域判定就清楚了：与鞍点稳定流形围成的包含稳定平衡点的区域，就是系统的稳定域。落在该区域内的初始状态，最终会稳定收敛到原点附近的稳定平衡点；落在区域外的状态，大概率会发散或跑到另一个稳定平衡点。

判断逻辑可以在仿真里用代码实现：给一批随机初始状态，分别积分 5 秒，看最终状态是否收敛到某个稳定平衡点邻域内。把所有能收敛的初始点画成散点，叠加在相平面图上，就能和临界轨迹互相验证。

N = 3000; beta_test = -0.25 + 0.5 * rand(N,1); omega_test = -0.7 + 1.4 * rand(N,1); T = 5; stable_flag = false(N,1); for i = 1:N [~, X] = ode45(@(t,z) vehicleDynamics(t,z,p), [0 T], [beta_test(i); omega_test(i)]); if abs(X(end,1)) < 0.02 && abs(X(end,2)) < 0.05 stable_flag(i) = true; end end scatter(beta_test(stable_flag), omega_test(stable_flag), 8, 'g', 'filled', 'Alpha', 0.4);

这种方法比单看相轨迹可靠，因为它是直接判终值，不会因为某条轨迹贴边而误判。注意阈值要结合控制精度需求设定，如果你做的是 ESP 标定，稳定域的保守性比精确性更重要，可以把阈值调得更严。

4.3 车速与附着系数的影响规律

做完基本流程后，强烈建议跑一遍参数扫描，因为临界轨迹不是固定不变的。车速 u 提高时，鞍点的位置会向内收缩，稳定域明显变小。原因很直观：高速下轮胎力提供的侧向加速度余量相对更小，状态稍微偏离，就可能突破附着极限，所以不安全区域扩大。

附着系数 μ 的影响更明显。μ 从 0.8 降到 0.4，相当于路面从干燥沥青变成湿滑路面，峰值侧向力减半，稳定域会大幅压缩，鞍点也会向原点靠拢。很多车辆稳定性控制器的逻辑，本质就是根据估算的 μ 来动态调整稳定边界，再决定是否介入制动或转向。

我在做批量仿真时，会把临界轨迹提取成一条边界曲线，再拟合成 μ 和 u 的函数，存成查找表给控制器用。这个方法虽然传统，但很稳，而且方便做硬件在环测试。

5. 高频问题与避坑心得

5.1 相轨迹乱飞：数值刚性与步长控制

刚开始画相轨迹时，最容易出现的问题是：某些初始点经过几秒后数值爆炸，轨迹一下子冲到图外。这不一定代表车辆真的失控，很多时候是数值积分步长太大导致的。普通乘用车动力学方程，尤其在轮胎饱和区，瞬时变化率可能相差几个数量级，ode45 会自动变步长，但偶尔还是需要调节精度参数。

我的经验是：把 odeset 里的 RelTol 设为 1e-6，AbsTol 设为 1e-7，尤其是两个状态量量纲不同（β 是弧度，ω 是 rad/s，数值范围差好几倍），AbsTol 必须设成向量：

opts = odeset('RelTol', 1e-6, 'AbsTol', [1e-7 1e-6]); [~, X] = ode45(@(t,z) vehicleDynamics(t,z,p), [0 T], z0, opts);

如果仍然发散，先检查状态方程里有没有符号错误，比如前轮转角 δ 的符号、侧偏角的符号，这些最容易阴人。把某个初始状态的导数打印出来，用手推公式核对几组值，比一味调积分器参数有效得多。

5.2 鞍点搜不到：初值、容差与算法选择

fsolve 本质是局部算法，初值离鞍点太远就会漏掉。我吃过不少亏，想象中鞍点应该在 β = ±0.15 附近，结果初值给到 0.2，fsolve 直接跑到奇异点。解决办法是“网格扫描 + 去重 + 类型判断”三步走，不要指望一次 fsolve 全找到。

另外，别忘了检查平衡点残差，有些点看似收敛但 fval 不为零，可能是数值奇异点而不是真正平衡点。我在代码里用norm(fun(eq_tmp)) < 1e-8做过筛，可以有效排除假平衡点。如果你用符号工具箱，vpasolve 能一次性把多个解列出来，但速度慢得多，模型复杂时不建议。

5.3 出图不直观：线宽、箭头与背景

相平面图信息密度大，最忌讳什么元素都堆在一起。方向场箭头本来就很多，如果相轨迹又用粗实线、红色，就容易和临界轨迹混淆。我的配色习惯是：

坐标轴范围也别贪大，围绕鞍点附近 0.5 rad/s 就够看清了。想让读者立刻抓住重点，还可以在图上用text标出“稳定域”“危险域”字样，比图例描述更直观。

5.4 一点个人实操建议

最后说点我自己的体会。做这种稳定性分析，别光盯着代码跑通，一定要先学会“读图”。拿到一张相平面图，先找原点附近的稳定平衡点，再看鞍点在哪，然后看临界轨迹穿过哪些区域。如果临界轨迹围成的稳定域形状和预期不符，先回头检查轮胎模型和参数，不要急着改控制器算法。

另一个建议是，把整个分析封装成一个脚本函数，输入参数是 u、mu、delta，输出是相平面图、鞍点坐标、临界轨迹数据。这样你后面换参数、出报告、做批处理都方便，也方便和其他同事分享。很多看似复杂的稳定性控制问题，到最后都会变成一张临界轨迹图的“平移”和“缩放”，你得把这张图的生成工具打磨利索。