最近我把二自由度车辆稳定性分析又重新做了一遍重点放在质心侧偏角β和横摆角速度r组成的相平面上并把鞍点、临界轨迹的绘制流程完整跑通了。这个工作很多人会用“相平面法”四个字一笔带过但真正动手时模型怎么简化、轮胎力怎么取非线性、平衡点怎么搜、临界轨迹从哪个特征方向出发每一步都会影响结果。这次我用MATLAB把整条链路串起来整理成一篇偏向实操的记录适合正在做车辆稳定性控制、ESP策略预研、或者课程里刚好卡在相平面仿真的朋友直接参考。下面所有代码我都按“能复制到本地跑出图”的标准写参数也给了具体值但更希望大家看明白每一步背后的原理。先给结论只要前、后轴轮胎力带有饱和甚至侧偏力下降段β-r相平面里常常同时存在一个稳定平衡点和一个鞍点。鞍点的稳定流形就是我们要找的临界轨迹它像一条分水岭把“还能回稳”和“即将失控”两种状态分开。这句话说起来简单真正在MATLAB里画对需要把状态方程、非线性轮胎模型、平衡点搜索、特征向量积分这几个环节都处理好。1. 为什么是 β-r 相平面稳定性边界问题的切入点1.1 两个状态量怎么来的常用的二自由度车辆模型也叫自行车模型只保留车辆侧向运动和横摆运动两个自由度。整车被简化成前后两个车轮侧向力集中到前轴和后轴上。车辆纵向速度V设为恒定状态量取质心侧偏角β和横摆角速度r。这两个状态量下的微分方程可以写成m·V·(dβ/dt r) Fy_f Fy_r Iz·dr/dt a·Fy_f - b·Fy_r其中m是整车质量Iz绕铅垂轴的转动惯量a、b分别是质心到前轴、后轴的距离Fy_f和Fy_r是前、后轴侧向力。第一个式子本质上就是“侧向力改变速度方向”所以dβ/dt的表达式里会带出一个横摆角速度r第二个式子则是横摆力矩平衡决定了车辆如何旋转。之所以用β和r作为相平面状态量是因为横摆角速度r可以直接通过陀螺仪测到质心侧偏角β则直观反映车身的“漂移程度”这两个量在ESP、车辆稳定性控制里都极端重要。实际仿真时先把车辆纵向速度固定下来。比如取V20m/s就相当于假设车辆在某个中高速工况下做准稳态横摆运动短时间内车速变化不大。这个假设对稳定性分析够用但要注意每条相轨迹、每张相平面图只对应一个固定车速。如果你想研究不同车速下的稳定性边界不是在同一张图里把线乱画而是每个车速单独算一张图再叠加对比。1.2 相平面中的鞍点与临界轨迹线性轮胎模型下给定一个前轮转角δ系统只有一个平衡点相轨迹要么全部收敛到该点要么全部发散不存在“部分区域稳定、部分区域失控”的复杂结构。实际问题不是这样。大侧偏角时轮胎力饱和甚至下降会让微分方程出现多个平衡点通常是一个稳定焦点或稳定结点加上一个鞍点。鞍点的名字很形象就像山谷垭口一样。它在相平面中同时有一个“吸引方向”和一个“排斥方向”。在鞍点附近沿某一个特征方向轨迹会向鞍点靠拢沿另一个特征方向轨迹会远离鞍点。真正决定稳定域边界的是鞍点的稳定流形也就是沿着稳定特征方向能够最终收敛到鞍点的所有状态点的集合。这条线两侧的轨迹行为完全不同一侧会被拉回稳定平衡点另一侧会发散。这条稳定流形就是标题里说的“临界轨迹”。我常用一个生活化类比稳定平衡点是山谷底部鞍点是山脊上的最低点临界轨迹就是山脊本身。水往两边流一滴水如果刚好落在山脊线上会顺着山脊去鞍点稍微偏一点点就会滑向其中一侧。汽车稳定性的相平面边界本质上就是这个山脊线的数学投影。2. 非线性轮胎与二自由度模型建模2.1 二自由度车辆模型和滑移角推导先约定正方向。根据右手坐标系车辆向左转向为正前轮转角δ取正值时车辆左转。前轮滑移角α_f等于车轮实际指向和车轮速度方向之间的差后轮滑移角α_r同样定义。推导后得到两个关键表达式α_f δ - β - a·r / V α_r -β b·r / V前轮因为位置在质心前方质心侧偏角β和横摆角速度r都会让车轮速度方向变化所以α_f里同时出现a·r/V这一项后轮在质心后方横摆角速度项前面的符号相反。这里所有角度单位都是弧度千万别用角度制直接代入否则力算出来全是错的。整车参数我按常见中型轿车标定取值如下m1573kgIz2873kg·m²a1.10mb1.60mCf80000N/radCr90000N/radμ0.85V20m/sδ4°。前、后轴垂向载荷按静态分配计算Fzfm·g·b/(ab)Fzrm·g·a/(ab)这里g取9.81m/s²。算出来的载荷用于轮胎饱和模型。不要把前轴侧偏刚度和后轴侧偏刚度搞混这两个车轮的载荷不同、侧偏特性也不同一旦写反平衡点位置和稳定域形状会明显变差。2.2 Fiala 饱和轮胎模型的选择相平面里要出现鞍点轮胎力必须非线性。线性轮胎模型只有一个平衡点画出来的结果就是所有向量都指向同一个方向看不出临界轨迹。我这次用简化的Fiala模型也叫抛物线侧向力模型。它的特点是小侧偏角时侧向力和α近似线性侧偏角超过某个峰后侧向力不再增大甚至回落并最终饱和到滑动摩擦力。正是侧向力“先升后降”的特性给系统引入了鞍点。简化Fiala模型的具体表达式在MATLAB里这样写function Fy fiala_approx(alpha, C, mu, Fz) alpha_sl 3*mu*Fz / C; % 近似饱和滑移角单位rad if abs(alpha) alpha_sl Fy -C*alpha C^2/(3*mu*Fz) * abs(alpha).*alpha ... - C^3/(27*(mu*Fz)^2) * alpha.^3; else Fy -mu*Fz * sign(alpha); end end这个式子里的C就是侧偏刚度mu是路面附着系数Fz是轴载荷。当alpha非常小的时候后面高次项几乎可以忽略Fy约等于-C·alpha和线性模型一致当alpha增大到alpha_sl时表达式正好给出Fy-mu·Fz也就是达到附着极限再往后就保持饱和。选择这个模型而不使用复杂的Magic Formula是因为它参数少、曲线光滑、便于求雅可比矩阵而且足够表现出相平面分析需要的非线性特征。需要说明的是真实轮胎在小侧偏角下侧向力通常也有松弛效应和残余侧偏力但二自由度相平面分析更关注定常侧向力特性所以这里的静态Fiala映射已经够用。如果你后面要接整车动力学再替换成魔术公式也不迟。2.3 MATLAB 函数封装把状态方程封装成一个独立的函数文件后面计算平衡点、画流场、积分相轨迹都调用它。输入x[β; r]输出xdot[dβ/dt; dr/dt]function xdot bicycle_beta_r(~, x, p) beta x(1); r x(2); alpha_f p.delta_rad - beta - p.a * r / p.Vx; alpha_r -beta p.b * r / p.Vx; Fyf fiala_approx(alpha_f, p.Cf, p.mu, p.Fzf); Fyr fiala_approx(alpha_r, p.Cr, p.mu, p.Fzr); xdot zeros(2,1); xdot(1) (Fyf Fyr) / (p.m * p.Vx) - r; xdot(2) (p.a * Fyf - p.b * Fyr) / p.Iz; end这里第一个参数写成~是为了直接兼容ode45的标准接口因为我们的微分方程本身不显式依赖时间t。参数p用结构体传递方便在main脚本里一次改车重、轴距、车速、附着系数不需要重复修改函数签名。这样封装之后除了画相平面还可以直接对某一组初始状态做时间历程仿真。3. 相平面、鞍点与临界轨迹的MATLAB绘制3.1 搜索平衡点网格预扫描 fsolve 精化平衡点的条件很简单dβ/dt0且dr/dt0。直接对整个状态空间用fsolve乱找绝对不行因为非线性方程组有多个解fsolve从一个初值出发只能收敛到一个零点而且很容易跑到数值溢出。我建议先做一次粗网格扫描把那些“函数值足够小”的点挑出来当作fsolve初值再精化求解。具体做法是在β和r的合理范围内打网格每个点调用一次bicycle_beta_r计算状态导数的模长。只要模长小于某个阈值就记下这个网格点作为“种子点”。然后对每个种子点做fsolve收敛后再做去重防止多个初值收敛到同一个平衡点上。阈值不能设得太严比如0.8~1.0左右比较可靠因为网格点一般不会正好落在平衡点附近。核心代码如下beta_vec linspace(-0.6, 0.2, 60); r_vec linspace(-0.4, 1.0, 60); seeds []; for b beta_vec for rr r_vec x [b, rr]; f bicycle_beta_r(0, x, p); if norm(f) 0.8 seeds(end1,:) x; %#okAGROW end end end all_eqs zeros(0,2); for i 1:size(seeds,1) x0 seeds(i,:); [xeq, fval, exitflag] fsolve((x) bicycle_beta_r(0,x,p), x0, ... optimoptions(fsolve,Display,off)); if exitflag 0 norm(fval) 1e-8 if isempty(all_eqs) || min(sum((all_eqs - xeq).^2, 2)) 1e-4 all_eqs(end1,:) xeq; %#okAGROW end end end实际跑下来大部分工况都会筛出两个平衡点。一个靠近原点附近另一个在质心侧偏角比较大的地方。前者通常实部为负数是稳定平衡点后者实部一正一负就是鞍点。如果某组参数下只找到一个平衡点多半是路面附着系数很高、转角很小非线性效果还没出来这时候需要适当减小μ或增大δ再看。3.2 鞍点的雅可比判断与特征方向找到平衡点后不能手工判断它是不是鞍点要算雅可比矩阵。雅可比矩阵在平衡点处的特征值符号决定了局部稳定性。对于二维系统两个特征值实部都为负是稳定点一正一负是鞍点两个实部都为正则是不稳定点。如果特征值是复数实部为负就是稳定焦点实部为正则是不稳定焦点。相平面分析里我们最关心的是鞍点因为它决定了稳定域边界。数值雅可比用中心差分计算简单可靠J zeros(2,2); h 1e-6; for i 1:2 xp xeq; xp(i) xp(i) h; fp bicycle_beta_r(0, xp, p); xm xeq; xm(i) xm(i) - h; fm bicycle_beta_r(0, xm, p); J(:,i) (fp - fm) / (2*h); end [V, D] eig(J); ev diag(D);如果isreal(ev)为真且ev(1)*ev(2) 0就说明这是鞍点。这时雅可比矩阵的特征向量列V里对应负特征值的那一列就是稳定特征方向。在画临界轨迹之前先把归一化后的特征向量取出来。特征向量的方向不需要过度纠结因为后面的积分会用正负两个方向同时追踪。3.3 绘制流场与多条相轨迹流场用quiver画。为了保证图上箭头长度不会互相遮挡最好把向量模做归一化。另外要注意单位问题状态方程里β是弧度但相平面图的横坐标习惯标成度看起来直观。把β转成deg时导数的x分量也要乘以180/pi否则箭头方向会在图上发生扭曲。建议的做法是在网格上用弧度计算导数随后把β和dβ/dt同时转换角度再用quiver绘制。beta_grid linspace(-0.6, 0.2, 25); r_grid linspace(-0.4, 1.0, 25); [Beta, R] meshgrid(beta_grid, r_grid); dBeta zeros(size(Beta)); dR zeros(size(R)); for i 1:numel(Beta) xdot bicycle_beta_r(0, [Beta(i), R(i)], p); dBeta(i) xdot(1); dR(i) xdot(2); end dBeta_deg dBeta * 180 / pi; speed sqrt(dBeta_deg.^2 dR.^2); quiver(Beta*180/pi, R, dBeta_deg./speed, dR./speed, 0.45, b);除了流场还要叠加若干条真实的相轨迹。可以从不同初值出发用ode45积分一段时间。积分时间终止条件要限制在状态范围内否则发散轨迹会一路跑出去图的范围控制不住。最简单的方法是指定一个积分终止事件或者直接取0到3秒先看趋势。对稳定区域内的初值轨迹会绕圈并最终收紧到稳定平衡点对不稳定区域内的初值r和β会一起增长轨迹一路冲向图边界。opts odeset(RelTol, 1e-6, AbsTol, 1e-8); init_list [-0.05 0.2; -0.15 0.35; -0.25 0.5; -0.05 0.7; -0.35 0.15]; hold on for i 1:size(init_list,1) [t, y] ode45((t,x) bicycle_beta_r(t,x,p), [0 3], init_list(i,:), opts); plot(y(:,1)*180/pi, y(:,2), Color, [0 0.4 0.8]); end初始状态最好不要直接从平衡点开始。从平衡点附近出发轨迹几乎不动看不出全局行为。多取几个初值点覆盖不同象限才能把流场的整体趋势展现出来。3.4 从鞍点出发绘制临界轨迹关键步骤临界轨迹不是随便取一个鞍点附近的点积分出来的而是需要沿着鞍点的稳定特征方向用极小的偏移量作为初值再对时间反向积分。这个步骤很多人容易做错。如果直接在鞍点处积分微分方程右端为零轨迹根本不动如果偏移量太大初值已经进入了非线性区积出来的轨迹会偏移真正的不变流形看起来歪歪扭扭。稳定特征向量对应负实部特征值沿这个方向的轨迹从远处向鞍点收敛。要画出整条稳定流形需要把初值放在“离鞍点超近、又恰好踩在稳定特征方向上”的位置然后让时间倒退从局部上等价于沿稳定流形远离鞍点。两个方向的正负偏移刚好覆盖稳定流形的两支。function [V_s, V_u] saddle_directions(xs, p) % 计算平衡点处雅可比、特征向量假设xs是鞍点 h 1e-6; J zeros(2,2); for i 1:2 xp xs; xp(i) xp(i) h; fp bicycle_beta_r(0, xp, p); xm xs; xm(i) xm(i) - h; fm bicycle_beta_r(0, xm, p); J(:,i) (fp - fm) / (2*h); end [V, D] eig(J); ev diag(D); [~, idx] sort(real(ev)); V_s V(:, idx(1)); % 对应最小实部通常是稳定方向 if length(idx) 2 V_u V(:, idx(2)); else V_u V(:, idx(1)); end V_s V_s / norm(V_s); V_u V_u / norm(V_u); end绘制稳定流形时从鞍点xs出发取x0 xs s * eps0 * V_s其中s取1和-1eps0取1e-5左右。积分区间用[0 -10]表示时间倒退。为了不让轨迹跑出画图范围可以给ode45加一个事件函数一旦β或r越界就终止积分。eps0 1e-5; opts odeset(RelTol, 1e-6, AbsTol, 1e-8, Events, event_out_of_range); for s [1 -1] x0 xs s * eps0 * V_s; [t, y] ode45((t,x) bicycle_beta_r(t,x,p), [0 -10], x0, opts); plot(y(:,1)*180/pi, y(:,2), r-, LineWidth, 2); end事件函数里判断β是否小于某个下限或大于上限以及r是否超出范围。代码可以写成function [value,isterminal,direction] event_out_of_range(~,x) value [x(1)0.8; 0.5-x(1); x(2)1.0; 1.8-x(2)]; isterminal ones(4,1); direction zeros(4,1); end这组事件值的范围按图窗范围写就行原则是只要轨迹跑出观察窗口就停不让ode45在白区域外空转。积分时间的长度决定了临界轨迹延伸多远。太短流形画不完整太长某些远离鞍点的轨迹数值误差会累积导致线的末段抖动。我的经验是先用[0 -10]试跑看曲线是否已经到达图边界再决定要不要缩短或延长。除了稳定流形还可以顺手把鞍点的不稳定流形画出来作为辅助。不稳定流形从鞍点沿V_u方向正向积分通常会向稳定平衡点方向卷绕。它不构成稳定域的边界但能帮助理解相平面结构。画的时候用虚线区分for s [1 -1] x0 xs s * eps0 * V_u; [t, y] ode45((t,x) bicycle_beta_r(t,x,p), [0 5], x0, opts); plot(y(:,1)*180/pi, y(:,2), m--, LineWidth, 1.5); end最终图里红色实线是白色稳定流形紫色虚线是不稳定流形蓝线是一般相轨迹。稳定流形两侧的蓝线走向截然不同很容易看出临界轨迹的分离作用。3.5 仿真结果怎么看在我给的参数工况下大致会出现两个平衡点。稳定平衡点落在β约等于-4.1°、r约等于0.31rad/s附近鞍点落在β约等于-13.5°到-14°、r约等于0.22到0.25rad/s附近。这个位置不是固定的当你把δ从4°改为6°或者把μ从0.85降到0.4稳定平衡点和鞍点都会移动临界轨迹的形状也会明显变化。看图时先找鞍点然后看两条红色稳定流形。如果观察窗口覆盖了它们可以看到稳定流形从鞍点两侧延伸出去大致像一张轻轻张开的嘴把稳定平衡点包围在中间。稳定平均点附近的蓝线会旋绕收紧说明车辆最终进入稳态圆周行驶而红色稳定流形外侧的蓝线要么β越来越负要么r持续变大说明车辆姿态正在失控。这组参数下稳态横摆角速度大约在0.3rad/s左右质心侧偏角只有几度对应一个正常的稳态转向工况。鞍点对应的β值更大说明此时车身侧偏已经很严重接近后轴饱和边界。临界轨迹恰恰落在两者之间。工程上可以把这个临界轨迹作为稳定包络的某个截面而不是简单设定“β小于某个固定角度就安全”。4. 实操经验、常见问题与工程扩展4.1 高频问题排查做仿真时最头疼的不是公式推导而是结果对不上。下面这张表是我调试过程中遇到的典型问题基本覆盖了大多数情况现象常见原因处理办法fsolve找不到鞍点初值离解太远或网格范围太小先粗网格扫描生成种子点后再fsolve最后去重临界轨迹画出来歪歪扭扭鞍点附近初值偏移量太大eps0取1e-4到1e-6之间方向严格沿稳定特征向量轨迹积分到一半NaN状态量越界轮胎模型进入饱和后数值突变加事件终止函数限定β和r范围降低RelTolquiver箭头方向看起来明显不对β轴用角度显示但向量的x分量没有换算将dβ/dt乘以180/pi后再归一化显示改变δ或V后平衡点找不到不同工况下解的位置差很多每次单独重新扫描不要沿用上一次的搜索区间结果只有一个平衡点轮胎非线性没被激发β和r范围太小增大δ到4°以上或降低μ到0.4~0.6再试ode45误差导致临界轨迹末端抖动默认容差对长时间积分不够设置RelTol1e-6AbsTol1e-8必要时使用事件终止第4条经常被忽略。很多教程直接用弧度画β出图后又把刻度标签改成度结果相轨迹箭头方向完全不对还以为是模型错了。我个人的习惯是所有计算都保留弧度画图前统一转成度并把导数中的β分量同步转换这样字体标注和数值都正常。还有一个容易被坑的地方Fiala模型中的alpha_sl不是固定的。它和轴载荷、侧偏刚度、附着系数都相关。换一组车型参数后若前、后轴载荷分配不同alpha_sl也会不同不能拿前一组的饱和角度套到后一组上。我每次改参数都会先打印一下Fy曲线和不同alpha下的力值确认没有出现数据溢出。4.2 工程中基于临界轨迹的稳定域应用纯画图当然不够更实际的价值在于给稳定性控制算法提供边界。工程上常常会把质心侧偏角-横摆角速度相平面中的临界轨迹拟合成简化表达式比如用双曲线、椭圆弧或多段直线包络作为车辆稳定性的安全包络线。控制器一旦检测到当前状态点越过临界轨迹就主动调节驱动力矩或施加修正横摆力矩。这种做法的优点是直观、计算量小适合实时估计。缺点是临界轨迹随车速V、前轮转角δ、路面附着系数μ都在移动如果直接用一张图的固定边界很容易误判。比较好的做法是把不同工况离线做好参数表把临界轨迹的关键点位置存为二维插值表运行时根据当前车速和转角实时选取最近的一组临界状态。从相平面的角度β和r的轨迹呈现明显方向转变的地方往往是车辆后轴侧向力到达饱和的时刻。后轴先饱和车辆就会趋于甩尾前轴先饱和车辆则趋于推头。二自由度模型虽然不能反映复杂悬架和载荷转移但用来判断车辆何时进入不稳定区依然是低成本、高效率的方法。我之前在真车上验证过β-r相平面边界和实车侧滑前兆趋势上吻合得很好。最后说一个我自己的习惯画临界轨迹时我几乎不会只跑一条流形曲线。我会在鞍点附近取多个微小偏移方向围绕稳定特征向量正负两侧加上几个小角度偏转然后把多条积分结果叠在一起看包络。原因是相平面的临界轨迹对数值误差很敏感单条曲线可能只是视觉上贴边多条线相互印证后才敢确认当前参数下的临界边界可信。这个习惯帮我排掉过很多看起来漂亮、实际错误的“假分界线”。如果你后续想把这套流程扩展得更深可以往两个方向走一是把二自由度模型升级为带垂向载荷转移的四轮模型在Fiala力里加上轴荷变化二是引入驾驶员闭环模型把前轮转角δ从常数改成有驾驶员反馈的时间函数。这样相平面分析就从“开环车辆固有不稳定性判断”走向“人-车-路闭环稳定性判断”思路一下子就能打开。