做车辆稳定性分析的朋友应该都绕不开“质心侧偏角-横摆角速度相平面”这张图。我记得第一次用 MATLAB 画出带鞍点和临界轨迹的二自由度车辆模型相平面时才真正理解了什么叫“临界失稳”。这个项目不是什么论文级高深算法而是一套非常经典的稳定性分析流程把二自由度车辆模型写成状态方程在相平面上找平衡点识别鞍点再追踪鞍点稳定流形画出临界轨迹从而判断质心侧偏角和横摆角速度在什么组合下车辆会失控。它解决的问题很直接——车辆稳定边界在哪里控制策略应该把状态约束在哪个区域内。适合正在做 ESP、ESC、四轮转向或底盘控制算法验证的工程师也适合刚入门车辆动力学建模的研究生照着复现。项目背景与核心概念解读1.1 二自由度车辆模型到底在描述什么二自由度车辆模型业内也常叫“自行车模型”不是说车上只坐两个人而是把车辆简化成前后两个轮胎侧向力作用点车身只有一个横向平移自由度和一个横摆转动自由度。它忽略了悬架运动、侧倾、俯仰、纵向车速变化这些复杂因素默认纵向车速恒定只关心车辆在转向输入下的侧向运动和绕垂直轴的转动。为什么这么极端简化还能用因为车辆横向稳定性分析的核心矛盾其实是轮胎侧向力与车身惯性力的平衡。只要把前后轴侧偏刚度、轴距、质心位置、质量、转动惯量这些关键参数保留住了质心侧偏角β和横摆角速度ω的动态特性就能被刻画得足够准确。ESP之类的稳定性控制器做状态反馈、边界判断用的也就是这两个状态量。这个模型虽然简单但它在做“相平面分析”时有个天然优势系统只有两个状态相平面是二维的所有轨迹都能画在平面图上平衡点、鞍点、极限环这些概念都可以直观地摆出来。换成三自由度甚至更高维度的模型相平面就只能投影了判断稳定边界就没这么直观。1.2 相平面、鞍点、临界轨迹在说哪三件事先说相平面。以质心侧偏角β为横轴、横摆角速度ω为纵轴每一组状态(β, ω)都对应平面上的一个点。给定一个初始状态车辆的运动方程会驱动这个点在平面上走出一条曲线这就是相轨迹。整张相平面图就是无数条相轨迹组成的“流场”能一眼看出不同初始状态下系统是收敛、发散还是绕圈。再说鞍点。鞍点是系统的平衡点也就是状态变化率为零的位置。车辆动力学模型通常不是线性的尤其在轮胎力进入饱和区后会出现多个平衡点。在这些平衡点里有一种特殊的鞍点在该点附近相轨迹一边被吸引、一边被排斥像马鞍一样。对车辆来说鞍点往往对应着临界稳定的失控状态两侧的相轨迹朝着截然不同的方向分流一边回归稳定一边滑向发散。最后说临界轨迹。通俗点讲临界轨迹就是鞍点的稳定流形它是一条或一组特殊的相轨迹把相平面分成两个区域一侧是状态能回归稳定的“安全区”另一侧是最终发散的“危险区”。很多控制策略的核心就是保证系统状态不越过这条临界轨迹。所以把三者画在一张图上等于同时回答了“哪里稳定、哪里临界、哪里必死”三个问题。模型搭建与仿真参数准备2.1 运动微分方程推导与符号约定二自由度车辆模型的微分方程国内教材和 ISO 标准里的符号略有差异但物理本质一样。我习惯用下面这套记号m整车质量单位 kgIz绕质心竖直轴的转动惯量单位 kg·m²u纵向车速单位 m/sδ前轮转角单位 rada质心到前轴距离单位 mb质心到后轴距离单位 mCf、Cr前后轮侧偏刚度单位 N/radβ质心侧偏角单位 radω横摆角速度单位 rad/sFyf、Fyr前后轴侧向力单位 N基于“纵向车速u恒定”的假设侧向运动方程和横摆力矩方程可以写成m·u·(dβ/dt ω) Fyf·cosδ FyrIz·(dω/dt) a·Fyf·cosδ - b·Fyr这里的前后轮胎侧偏角为αf δ - (β a·ω/u)αr - (β - b·ω/u)注意这套式子默认侧偏角符号遵循SAE标准侧偏角为正时侧向力为负。如果你用ISO坐标系符号方向要统一检查否则画出来的鞍点位置会左右翻转。如果进一步假设前后轮轮胎力处于线性区即 Fyf -Cf·αfFyr -Cr·αr那就得到线性二自由度模型。但线性模型在稳定性分析里有个致命局限它没法反映轮胎饱和系统最多只有一个平衡点也不容易出现真正的鞍点结构。因此我在项目里用的是非线性轮胎模型这样相平面才“有看头”临界轨迹也才会出现。2.2 参数选取与无量纲化示例参数我直接用一台中级轿车的典型值方便你复现时对照参数数值单位m1500kgIz2500kg·m²a1.2mb1.6mCf80000N/radCr100000N/radu20m/s轮胎峰值侧向力需要结合垂直载荷和路面附着系数估算。前轴静态载荷约为 Fz_f m·g·b/(ab)后轴约为 Fz_r m·g·a/(ab)。如果附着系数 μ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.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(该平衡点为鞍点); end3.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.^2V.^2), V./sqrt(V.^2V.^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.1 稳定流形与临界轨迹的关系临界轨迹本质上是鞍点的稳定流形在二维非线性系统里它由从鞍点出发的两条“臂”组成。稳定流形上的点随着时间推进会被鞍点“吸向”鞍点但一旦偏离那么一点点系统就会沿着不稳定方向飞出去。因此在相平面上稳定流形就像一个分水岭把收敛到稳定平衡点的区域和发散到无穷远的区域隔开。要绘制稳定流形最标准的数值方法是首先求鞍点处的雅可比矩阵找到在鞍点附近把轨迹“指向鞍点”的特征方向实部为负的特征值对应的特征向量。然后在该方向上给一个小扰动作为初始点对原系统做负方向时间积分。因为负时间下“指向鞍点”的稳定方向会变成从鞍点向外延伸的方向所以你就能沿着整条稳定流形把它“拓”出来。我写了个便于理解的小示例假设已经找到鞍点坐标 eq_saddle 和稳定特征向量 v_seps 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.1 相轨迹乱飞数值刚性与步长控制刚开始画相轨迹时最容易出现的问题是某些初始点经过几秒后数值爆炸轨迹一下子冲到图外。这不一定代表车辆真的失控很多时候是数值积分步长太大导致的。普通乘用车动力学方程尤其在轮胎饱和区瞬时变化率可能相差几个数量级ode45 会自动变步长但偶尔还是需要调节精度参数。我的经验是把 odeset 里的 RelTol 设为 1e-6AbsTol 设为 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.2fsolve 直接跑到奇异点。解决办法是“网格扫描 去重 类型判断”三步走不要指望一次 fsolve 全找到。另外别忘了检查平衡点残差有些点看似收敛但 fval 不为零可能是数值奇异点而不是真正平衡点。我在代码里用norm(fun(eq_tmp)) 1e-8做过筛可以有效排除假平衡点。如果你用符号工具箱vpasolve 能一次性把多个解列出来但速度慢得多模型复杂时不建议。5.3 出图不直观线宽、箭头与背景相平面图信息密度大最忌讳什么元素都堆在一起。方向场箭头本来就很多如果相轨迹又用粗实线、红色就容易和临界轨迹混淆。我的配色习惯是元素颜色线型方向场浅灰细箭头普通相轨迹蓝色0.4pt 细线临界轨迹红色2pt 粗线稳定平衡点绿色实心圆点鞍点红色实心五角星坐标轴范围也别贪大围绕鞍点附近 0.5 rad/s 就够看清了。想让读者立刻抓住重点还可以在图上用text标出“稳定域”“危险域”字样比图例描述更直观。5.4 一点个人实操建议最后说点我自己的体会。做这种稳定性分析别光盯着代码跑通一定要先学会“读图”。拿到一张相平面图先找原点附近的稳定平衡点再看鞍点在哪然后看临界轨迹穿过哪些区域。如果临界轨迹围成的稳定域形状和预期不符先回头检查轮胎模型和参数不要急着改控制器算法。另一个建议是把整个分析封装成一个脚本函数输入参数是 u、mu、delta输出是相平面图、鞍点坐标、临界轨迹数据。这样你后面换参数、出报告、做批处理都方便也方便和其他同事分享。很多看似复杂的稳定性控制问题到最后都会变成一张临界轨迹图的“平移”和“缩放”你得把这张图的生成工具打磨利索。