在底盘稳定性分析这个领域里“相平面”这三个字出现的频率特别高尤其是做ESP标定、四轮转向控制、或者是写车辆动力学毕业论文的时候。我最早接触二自由度车辆模型配相平面是因为当时要评估一台车在低附着路面上的失稳边界。项目需求很简单给定车速和前轮转角把车辆在质心侧偏角-横摆角速度坐标系下的全部运动趋势画出来再把鞍点和临界轨迹标出来。做完之后你会发现车辆是否稳态、失控边界在哪、控制器阈值该定在什么范围全都清清楚楚写在这张图上。这篇文章就把这套“二自由度车辆相平面鞍点临界轨迹”的MATLAB仿真方法完整拆开讲一遍从模型推导到代码实现再到调试坑点尽量让读者能照着复现。1. 项目到底在解决什么问题为什么要画β-r相平面1.1 二自由度模型的“二”在哪里整车运动如果仔细展开有纵向、侧向、垂向、侧倾、俯仰、横摆等多个自由度轮子还有旋转自由度分析起来非常复杂。但如果我们只关心一件事——车辆在高速行驶时会不会甩尾失控那么真正起决定性作用的是两个状态量质心侧偏角β和横摆角速度r。这就是二自由度模型的核心思想把车辆想象成一辆“自行车”前后轴分别等效成一个轮子纵向车速恒定忽略侧倾和垂向运动只保留侧向平移和横摆转动这两个自由度。这样的简化不是偷懒而是工程上的刻意取舍。稳定控制关注的是车辆在极限工况下的横向动态纵向速度变化往往比横向失稳慢得多可以把它当作常数处理。侧倾运动会对轮胎载荷有影响但在一阶分析阶段先把主要矛盾β和r的动态耦合看清楚再逐步细化也不迟。MATLAB仿真中用二自由度模型跑相平面本质上就是用最少的参数描述车辆最核心的失稳机理。1.2 相平面上的鞍点和临界轨迹对应了什么物理现象相平面简单说就是把状态空间铺开成一个平面横轴是质心侧偏角β纵轴是横摆角速度r。车辆当前处于什么状态就是平面上的一个点状态怎么变化就是平面上该点处的“箭头”。在这个平面上箭头指向会形成几种典型结构。箭头都指向某个点说明这里是稳定的平衡点车辆状态最终会收敛到这里箭头都从某个点向外发散说明这是不稳定平衡点还有一种特殊的点一部分方向流入、另一部分方向流出这就是鞍点。鞍点在车辆动力学里非常重要因为鞍点的稳定流形——也就是“临界轨迹”——恰好把相平面分成了两个区域临界轨迹内侧的状态最终会收敛到稳态外侧则会发散失稳。拿一个生活化的比喻想象一个山谷盆地的分水岭。盆地里的小雨滴最终都会汇入最低处稳定平衡点但分水岭上的水滴如果稍微偏向一侧就会流向完全不同的方向。临界轨迹就是这条“分水岭”鞍点就是分水岭上的一个鞍部两边是山坡两端是山脊。理解了这层物理含义再看相平面图就不会觉得它只是一堆箭头和曲线了。1.3 这套仿真能力的工程用途我在实际项目里用这套仿真解决过三类问题。第一类是失稳边界预测在某一车速和转角下车辆能承受的最大质心侧偏角到底是多少、横摆角速度超过多少就救不回来直接看临界轨迹的位置就能估算。第二类是控制器触发阈值设计ESP类控制器的介入时机本质上就是判定“当前状态是否接近临界轨迹”把相平面上的稳定裕度换算成触发阈值比拍脑袋定经验值可靠得多。第三类是开环特性评估在开发四轮转向或者LQR转向控制器之前先画出原车的相平面看看系统的固有阻尼、稳定域大小、鞍点位置控制器的设计目标就非常明确了。这套工具看起来只是“画图”但它是整条底盘稳定性控制链路上的第一步也是沟通车辆动力学理论到工程标定之间的一座重要桥梁。2. 模型建立与关键数学细节2.1 状态方程推导与符号约定二自由度模型的推导并不复杂关键在于符号约定要统一。整车质量m绕质心的横摆惯量Iz质心到前轴距离a质心到后轴距离b前轮侧偏刚度Cf后轮侧偏刚度Cr纵向车速Vx恒定前轮转角δ作为输入。基于牛顿第二定律侧向力和横摆力矩平衡可以得到侧向力方程mVx(β̇ r) Fyf·cosδ Fyr横摆力矩方程Iz·ṙ a·Fyf·cosδ - b·Fyr其中Fyf是前轴等效侧偏力Fyr是后轴等效侧偏力。轮胎侧偏角在小角度假设下可以近似为前轮侧偏角αf δ - β - a·r/Vx后轮侧偏角αr -β b·r/Vx把侧偏角代到轮胎模型里就能得到关于β和r的两个一阶微分方程正好组成状态方程。这里需要注意Vx不能取得太小否则侧偏角公式里的ar/Vx会非常大模型本身失效。一般Vx在10m/s以上使用这个简化模型比较合理。2.2 为什么必须上非线性轮胎模型很多入门资料用的是线性轮胎模型也就是把轮胎侧偏力直接写成Fy -C·α。这在线性区没问题做稳态响应或者模态分析足够用。但只要转向角稍大或者路面附着系数偏低轮胎进入非线性饱和区后线性模型给出的相平面就是一个单一稳定焦点根本不会有鞍点也就画不出临界轨迹。要观察失稳边界就必须让轮胎力随着侧偏角增大而饱和。我用的是Fiala轮胎模型它用三段式曲线描述侧偏力当|α| α_sl时Fy -C·α C²/(3μFz)·|α|·α - C³/(27(μFz)²)·α³当|α| ≥ α_sl时Fy -μFz·sign(α)这里的α_sl 3μFz/C是轮胎线性区结束的临界侧偏角。这个模型很实用参数只有侧偏刚度C、附着系数μ、垂直载荷Fz三个却能把轮胎从线性到饱和再到完全滑移的过渡描述出来足够用于相平面定性分析。2.3 平衡点、鞍点和稳定流形的数学定义相平面上任意一个“箭头归零”的点就是平衡点数学上满足β̇0和ṙ0。由于轮胎模型是非线性的通常有多个平衡点。为了判断每个平衡点的性质需要对状态方程在平衡点处求雅可比矩阵再看特征值。雅可比矩阵的特征值如果实部全为负平衡点是稳定的车辆状态会趋向它如果实部一正一负就是鞍点——这正是临界轨迹的“源头”如果实部全为正则是完全不稳定的平衡点。鞍点的稳定流形就是沿着对应负实部特征值特征向量方向的一簇轨迹这簇轨迹在相平面上绘制出来就是临界轨迹即稳定域边界。在实际计算中雅可比矩阵可以用数值差分法近似得到不需要推出解析表达式这对非线性轮胎模型来说非常方便。3. MATLAB实现从参数到相平面图3.1 参数初始化与状态方程函数我常用的一组示例参数如下表对应一台中型乘用车在干沥青路面上的工况参数符号取值单位整车质量m1500kg横摆惯量Iz2500kg·m²质心到前轴a1.20m质心到后轴b1.40m前轮侧偏刚度Cf80000N/rad后轮侧偏刚度Cr100000N/rad附着系数μ0.85-纵向车速Vx20m/s前轮转角δ0.02rad对应的MATLAB参数结构体可以这样写p.m 1500; p.Iz 2500; p.a 1.20; p.b 1.40; p.Cf 80000; p.Cr 100000; p.mu 0.85; Vx 20; % 纵向车速 m/s delta 0.02; % 前轮转角 rad状态方程函数是整套代码的核心建议单独写成函数文件便于相平面扫描和ODE积分复用function [dbeta, dr] vehicle_2dof(beta, r, p, Vx, delta) % 二自由度车辆模型用于相平面分析 % beta : 质心侧偏角 [rad] % r : 横摆角速度 [rad/s] % p : 车辆参数结构体 % Vx : 纵向车速 [m/s] % delta: 前轮转角 [rad] % 静态垂直载荷分配 Fzf p.m * 9.81 * p.b / (p.a p.b); Fzr p.m * 9.81 * p.a / (p.a p.b); % 轮胎侧偏角 alpha_f delta - beta - p.a * r / Vx; alpha_r -beta p.b * r / Vx; % Fiala 轮胎模型 Fyf tire_fiala(alpha_f, p.Cf, p.mu, Fzf); Fyr tire_fiala(alpha_r, p.Cr, p.mu, Fzr); % 整车动力学方程 dbeta (Fyf * cos(delta) Fyr) / (p.m * Vx) - r; dr (p.a * Fyf * cos(delta) - p.b * Fyr) / p.Iz; end轮胎模型函数单独写function Fy tire_fiala(alpha, C, mu, Fz) % Fiala 轮胎模型简化形式 alpha_sl 3 * mu * Fz / C; % 线性区临界侧偏角 alpha_abs abs(alpha); if alpha_abs alpha_sl Fy -mu * Fz * sign(alpha); % 完全滑移区 else Fy -C * alpha ... C^2 / (3 * mu * Fz) * alpha_abs .* alpha ... - C^3 / (27 * (mu * Fz)^2) * alpha.^3; % 过渡区 end end这里建议所有角度都用弧度制避免后面度/弧度混用造成的各种诡异结果。我最早被坑过一次把δ按角度传进去平衡点位置完全对不上。3.2 相平面网格与向量场绘制画相平面第一步是生成网格然后在每个网格点调用状态方程函数计算方向向量用quiver画出箭头。beta_g linspace(-0.25, 0.25, 17); r_g linspace(-0.6, 0.6, 17); [BGR, RGR] meshgrid(beta_g, r_g); dB zeros(size(BGR)); dR zeros(size(BGR)); for i 1:numel(BGR) [dB(i), dR(i)] vehicle_2dof(BGR(i), RGR(i), p, Vx, delta); end figure(Color, w, Position, [120 120 720 600]); hold on; grid on; box on; quiver(BGR, RGR, dB, dR, 0.8, Color, [0.85 0.85 0.85], LineWidth, 0.6);网格密度取17×17比较合适。太密箭头挤在一起看不清太疏又看不出流场走向。quiver的缩放因子0.8是我试过比较舒服的显示效果箭头长度太长会盖住轨迹太短则流场走向不明显。不过只靠箭头还是不够直观最好再叠加几条从不同初始状态出发的实际积分轨迹。用ode45从典型初始点出发让系统自己“走”几步轨迹线会清晰展示收敛和发散的行为tspan [0 5]; opt odeset(RelTol, 1e-6, AbsTol, 1e-8, ... Events, (t, x) event_out_of_domain(t, x)); for b0 [-0.2 -0.1 0 0.1 0.2] for r0 [-0.35 -0.15 0 0.15 0.35] x0 [b0; r0]; [~, xout] ode45((t, x) state_ode(t, x, p, Vx, delta), tspan, x0, opt); plot(xout(:,1), xout(:,2), Color, [0 0.45 0.75], LineWidth, 1.0); end end这里的state_ode是一个适配ode45的包装函数function dx state_ode(~, x, p, Vx, delta) dx zeros(2,1); [dx(1), dx(2)] vehicle_2dof(x(1), x(2), p, Vx, delta); end事件函数的作用是限制轨迹不会飞到根本不可能出现的远处同时也能避免长时间积分浪费时间function [value, isterminal, direction] event_out_of_domain(~, x) value(1) 0.5 - abs(x(1)); % |beta| 0.5 rad value(2) 1.2 - abs(x(2)); % |r| 1.2 rad/s isterminal [1; 1]; direction 0; end3.3 搜索平衡点并判定鞍点平衡点是相平面上f(x)0的点。由于非线性方程组没有解析解我用两层策略先在粗网格上找“场强”局部极小点作为初值再用fsolve精求解。场强的定义是每个网格点处状态变化量的模长即sqrt(dβ² dr²)。平衡点附近场强接近零所以找局部极小值就能得到很好的初值候选beta_vec linspace(-0.30, 0.30, 61); r_vec linspace(-0.8, 0.8, 61); [BG, RG] meshgrid(beta_vec, r_vec); fmag zeros(size(BG)); for i 1:numel(BG) [db, dr] vehicle_2dof(BG(i), RG(i), p, Vx, delta); fmag(i) sqrt(db^2 dr^2); end % 找局部极小值作为 fsolve 初值 mask islocalmin(fmag, MinProminence, 0.02, FlatSelection, center); [row, col] find(mask); options optimoptions(fsolve, Display, off, ... TolFun, 1e-12, TolX, 1e-12); eq_points []; for k 1:length(row) x0 [BG(row(k), col(k)); RG(row(k), col(k))]; xeq fsolve((x) state_ode(0, x, p, Vx, delta), x0, options); % 过滤未收敛的结果 if max(abs(xeq - x0)) 0.3 % 去重如果与已有平衡点太近则跳过 if isempty(eq_points) || min(sqrt(sum((eq_points - xeq.).^2, 2))) 1e-3 eq_points [eq_points; xeq.]; end end end fprintf(找到 %d 个平衡点\n, size(eq_points, 1)); disp(eq_points);接下来对每个平衡点计算数值雅可比矩阵和特征值判断类型function J num_jacobian(f, x) h 1e-6; J zeros(2, 2); for k 1:2 xp x; xm x; xp(k) xp(k) h; xm(k) xm(k) - h; J(:, k) (f(xp) - f(xm)) / (2*h); end end function type classify_equilibrium(xeq, p, Vx, delta) f (x) state_ode(0, x, p, Vx, delta); J num_jacobian(f, xeq(:)); ev eig(J); if all(real(ev) 0) type stable; elseif all(real(ev) 0) type unstable; elseif any(real(ev) 0) any(real(ev) 0) type saddle; else type center; end end在这组参数下我跑出来的典型结果是一个稳定焦点位于中央偏左下方两个鞍点分布在它的两侧。稳定焦点对应车辆的正常稳态转向鞍点则对应临界失稳状态。如果你跑出来的平衡点只有稳定焦点而没有鞍点大概率是前轮转角太小或者路面附着太高轮胎基本处于线性区此时可以适当增大δ或者降低μ再看。3.4 临界轨迹绘制鞍点稳定流形的数值求法有了鞍点接下来画临界轨迹。关键点在于鞍点处的雅可比矩阵有一正一负两个实特征值负实部特征值对应的特征向量方向就是稳定流形的切方向。在鞍点附近沿着这个方向取一个极小的偏移点然后让系统“倒着走”——也就是逆时间积分——轨迹就会沿着稳定流形向外延伸画出完整的分界线。% 在相平面上叠加临界轨迹 h_crit []; for k 1:size(eq_points, 1) type classify_equilibrium(eq_points(k,:), p, Vx, delta); if ~strcmp(type, saddle) continue; end xeq eq_points(k,:); f (x) state_ode(0, x, p, Vx, delta); J num_jacobian(f, xeq(:)); [V, D] eig(J); [~, idx] min(real(diag(D))); % 取负实部最小的特征值 v_stable V(:, idx); v_stable v_stable / norm(v_stable); for sgn [-1 1] x0 xeq(:) sgn * 1e-3 * v_stable; tback [0 -8]; opt_back odeset(RelTol, 1e-8, AbsTol, 1e-10, ... Events, (t, x) event_out_of_domain(t, x)); [~, xcrit] ode45((t, x) state_ode(t, x, p, Vx, delta), tback, x0, opt_back); h_crit(end1) plot(xcrit(:,1), xcrit(:,2), -, ... Color, [0.85 0.33 0.1], LineWidth, 2.2); end end这里有个细节要重点提醒逆时间积分时积分器容差要设得比普通轨迹严格很多。临界轨迹本身是一组“中性稳定”的曲线如果RelTol太宽松数值误差会让轨迹从流形上“滑下去”导致画出来的分界线向内或向外偏移看起来歪歪扭扭的。我一般设RelTol1e-8AbsTol1e-10。完整的绘图脚本最后再加上坐标轴标签和图例xlabel(质心侧偏角 \beta [rad]); ylabel(横摆角速度 r [rad/s]); title(sprintf(二自由度车辆相平面特性 V_x%.0f m/s \\delta%.3f rad, Vx, delta)); legend(h_crit(1), {临界轨迹鞍点稳定流形}, Location, northeastoutside); axis([-0.3 0.3 -0.7 0.7]);跑出来的图上你应该能看到中央蓝色轨迹线快速收敛到稳定焦点周围箭头全都指向那个焦点红色临界轨迹从两个鞍点向两侧延伸像两条边界线一样把相平面分成了内外两个区域。落在红色曲线内部的状态点无论初始怎么偏离最终都会被拉回稳态一旦跨过红色曲线状态量就会迅速发散β急剧增大r也无法收敛——对应到实际车辆上就是甩尾、侧滑、彻底失去控制。4. 实际运行中的常见问题与调试心得4.1 鞍点“消失”了参数变化对稳定域的影响很多读者照着代码跑了一遍发现鞍点找不着图上一片箭头都指向中央完全没有分界线。这种情况绝大多数是因为参数取在了“安全区”里。前轮转角很小、车速适中、附着系数很高时轮胎几乎都工作在线性区非线性饱和效应不明显系统在相平面上只有一个全局稳定焦点自然没有鞍点。想让鞍点出现有两个常用办法。一是增大前轮转角比如从0.02rad逐渐加到0.05rad以上二是降低附着系数比如把μ从0.85降到0.4或0.3模拟雨天或冰雪路面。你会发现随着参数变化两个鞍点会从远处逐渐向中央移动临界轨迹围出的稳定域越来越小。等到参数超过某一临界值稳定焦点和鞍点相遇合并平衡点数量骤减系统彻底失去稳定平衡点。这个过程在动力学里叫鞍结分岔在实车上对应的就是“到了某个车速和转角组合下车辆无论如何都救不回来”的现象。4.2 单位不一致与收敛性坑这听起来像基础问题但我在帮别人看代码时遇到最多。侧偏角公式里δ、β、r/Vx的单位必须一致。MATLAB三角函数默认输入是弧度所以δ、β、α都用弧度横摆角速度用rad/s。如果习惯用角度制要么所有输入统一转成弧度要么在轮胎模型里单独转换切忌混用。还有fsolve初值问题。直接给一个随机初值往往收敛不到想要的平衡点甚至跑偏到远方的错误解。网格扫描局部极小值粗定位fsolve精求解这套流程虽然多写几行代码但实际效率反而更高。记得加收敛判断和去重逻辑否则同一平衡点会被重复计算很多次。4.3 临界轨迹向外漂移的问题画临界轨迹最常见的困扰是红线上半段规则、下半段却突然朝某个方向飘出去。这通常是因为逆时间积分的种子点偏离鞍点太远或者积分容差太松。种子点离鞍点不能太远一般取1e-3量级的偏移太大会脱离稳定流形附近的线性区域轨迹一开始就歪了。如果已经将容差收到1e-8还是漂移可以考虑把逆时间积分拆成多段每段从上一段终点附近重新校正到流形方向但这一步通常用不上保持密切观察即可。另外事件函数防止轨迹飞出绘图范围是必要的推荐保留。没有事件终止时逆时间积分可能会跑出β±5rad这种毫无物理意义的状态白白浪费时间还可能触发数值异常。4.4 ODE积分器选择的经验相平面系统绝大多数情况下是光滑的用ode45就行。但在某些参数组合下比如μ非常低、前轮转角很大的时候状态方程会变“硬”ode45可能需要极小的步长速度很慢。这时候换成ode15s试试往往能明显提速。判断“硬不硬”的经验是看同一条轨迹ode45和ode15s的时间消耗差异超过两倍就往ode15s切换。瞬时积分速度本身并不是问题但大范围扫参数时累计差距非常可观。还有ODE选项里的MaxStep也建议设一下比如0.2秒避免在某些区域步长自动变得太大导致轨迹细节丢失。个人实操体会与扩展建议把这套相平面工具在手里跑顺之后我最直观的感触是它把车辆稳定性从“玄学般的驾驶感受”变成了一张可以度量的图。临界轨迹围出的区域大小直接就是当前工况下的稳定裕度鞍点离当前状态点的距离可以作为ESP介入的提前量指标。我后来在项目中做过的最有效扩展是批量扫不同的Vx和δ组合把每种工况下的鞍点位置和稳定域面积记录成MAP表然后再叠加车辆状态点实时判断是否需要控制介入。这套离线相平面在线查表的思路比单纯用固定阈值做稳定性判断要可靠得多。如果后续想深入还可以在这个基础上做三件事一是把Fiala模型换成魔术公式对比轮胎模型差异对临界轨迹的影响二是把司机转向输入从恒定δ改成随时间变化观察瞬态工况下相平面结构的迁移三是把相平面工具封装成函数批量生成不同参数下的对比图做出类似“鞍点规划”的决策辅助图。只要基础设施到位每一条扩展路径都是水到渠成的事。