做车辆稳定性控制的人，手机相册里大概率存过这样一张图：横轴是质心侧偏角β，纵轴是横摆角速度r，整个平面上被几条曲线硬生生切开，一侧代表能稳住，一侧代表彻底失控。这张图就是二自由度车辆模型下的β-r相平面图。很多课程和论文都会提到它，但真到自己用MATLAB复现时，才发现从模型搭建、鞍点求解到临界轨迹绘制，处处是细节。这篇文章把我从零开始做这套仿真的完整思路和代码写出来，包括我踩过的坑和验证技巧，给正在做车辆稳定性研究、底盘控制算法验证以及相关课程设计的朋友参考。
1. 为什么用相平面分析车辆稳定性
1.1 二自由度模型：抓住主矛盾，丢掉细枝末节
先问一个问题：为什么几乎所有车辆稳定性教材都从二自由度模型开始？因为车辆横向动力学在低频段的主矛盾，其实就是两个自由度：侧向平移和横摆转动。其它自由度，比如侧倾、俯仰、垂向跳动，对车身稳定性的第一层影响没那么直接，在做控制策略验证时完全可以在后续再补。二自由度模型把车辆简化为一辆刚性自行车，前轮负责输入转角δ（也就是驾驶员或上层控制器给的转向角），后轮随动，车身通过质心侧偏角β和横摆角速度r两个状态量描述。这样描述的不是车辆在二维路面上的位置，而是车辆当前“姿态”的变化率关系。
模型方程本身不复杂。对纵向速度恒定的情况，侧向力与横摆力矩平衡方程可以写成：
m * V * (dβ/dt + r) = Fyf + Fyr Iz * dr/dt = a * Fyf - b * Fyr其中m是整车质量，Iz是绕垂向轴的转动惯量，a是质心到前轴距离，b是质心到后轴距离，Fyf和Fyr是前后轴等效侧偏力。前后轮侧偏角与状态量的几何关系为：
αf = β + a*r/V - δ αr = β - b*r/V注意δ是前轮转角，V是纵向车速。当轮胎工作在线性区时，侧偏力Fy = -C * α，C就是侧偏刚度。把线性关系代入，就得到经典的二自由度线性状态方程。线性模型的优势是有解析解，可以做传递函数、特征值分析，很多控制器的初始参数都从它推导。
1.2 β-r相平面的坐标系选择逻辑
为什么偏偏选β和r作为相平面坐标，而不是侧向速度vy和横摆角速度r？主要原因有两条。第一，质心侧偏角直接决定了车身实际运动方向与车头方向的夹角，驾驶员对β的感知非常敏锐，在后轮侧滑失稳时，β往往最先快速增大。第二，轮胎侧偏角αf、αr都用β和r的线性组合表达，画在β-r平面里，侧偏角的等值线是直线，物理含义一目了然。如果换成vy-r坐标系，同一张图上的稳定边界会扭曲，不够直观。
相平面的横纵轴量程一般怎么定？常规沥青路面、车速20m/s左右时，β变化范围通常在±0.2~±0.3 rad左右，r大约在±0.8~1.0 rad/s。超过这个范围基本已经进入严重侧滑状态，继续画下去意义不大。画图之前先把量程确定好，后面找平衡点、追踪临界轨迹都会少很多麻烦。
2. 鞍点与临界轨迹：稳定边界的本质
2.1 平衡点的“好”与“坏”：稳定焦点 vs 鞍点
相平面图上有一种特殊位置，叫作平衡点。所谓平衡点，就是状态导数都为零的位置，即dβ/dt=0、dr/dt=0同时满足。直观理解，就是车辆在固定方向盘转角下，质心侧偏角和横摆角速度都不再变化，进入一种稳态转弯状态。低速小转角工况下，这样的平衡点通常只有一个，对应稳定的稳态转向；但在高速大转角或低附着路面，由于轮胎力饱和，侧偏力不再是α的线性函数，平衡点可能不止一个。平衡点的“好坏”，需要看它附近的流线走向：如果把系统状态微微推离平衡点，扰动后还能自动回到这里，就是稳定平衡点或稳定焦点；如果微微推离后状态反而越走越远，就是不稳定的。
怎么精确区分？在平衡点处对系统做泰勒展开，保留一阶项，得到2×2的雅可比矩阵J。算J的特征值，就有了判据。
| 特征值分布 | 平衡点类型 | 相平面表现 |
|---|---|---|
| 两个实部都为负 | 稳定节点/稳定焦点 | 状态向平衡点收敛 |
| 一个正、一个负 | 鞍点 | 状态沿特定方向被放大 |
| 至少一个正实部 | 不稳定节点/不稳定焦点 | 状态远离平衡点 |
鞍点在车辆动力学里不是抽象的数学名词。它对应的状态往往就是车辆进入后轴侧滑的临界状态：沿一个方向，车辆还能勉强保持姿态；沿另一个方向，车辆会迅速甩尾。物理上可以想象成一个山坳口——山坳口本身是山谷里的一个点，但两侧都是山坡，站上去稍微偏一点，就会滑向不同的山沟。鞍点的β、r组合值越高，说明当前稳定裕度越差。
2.2 临界轨迹就是鞍点区域延伸出的分界线
有了鞍点，临界轨迹就有着落了。临界轨迹本质上是鞍点的不变流形，也就是从鞍点附近出发、在非线性微分方程驱动下形成的特殊轨道。把系统在鞍点处线性化，能得到两个特征向量方向：一个对应正特征值，代表状态最容易增大的方向（不稳定方向）；另一个对应负特征值，代表状态容易被拉回的方向（稳定方向）。从鞍点出发，沿不稳定方向取一个很小的偏移量作为初始点，再对完整的非线性车辆模型做数值积分，这条轨迹就会沿着鞍点的“脊线”延伸出去，把相平面分隔成不同区域。
这条过鞍点的轨迹，就是通常说的临界轨迹，也叫相平面分界线。分界线一侧的状态最终会收敛到稳定平衡点，另一侧则会发散或者滑向另一个吸引域。对车辆稳定性分析而言，临界轨迹围成的区域可以看作车辆的安全稳定域：状态点在这个区域内，车辆有自行恢复的趋势；一旦越界，只靠轮胎本身的非线性特性很难自救，必须依赖ESP之类的主动控制介入。
这里有个细节：临界轨迹不是随便从鞍点拉两条直线就算完事。鞍点附近的局部线性化只决定了初始方向，真正的分界线要在完整非线性模型下积分延伸，因为远离鞍点后，轮胎力饱和特性会显著改变轨道的曲率。这也是为什么代码实现里，必须把非线性轮胎模型写进被积函数。
2.3 从鞍点到“鞍点规划”的控制视角
近两年讨论比较多的一个词叫鞍点规划。听起来玄乎，放到车辆稳定性语境下并不复杂：既然鞍点是不稳定边界上的关键位置，那么在做轨迹规划或控制目标设计时，就要主动让参考状态远离鞍点，而不是让控制器把状态往鞍点附近推。比如路径跟踪算法规划出期望β和r轨迹，如果这条轨迹要穿过鞍点附近区域，那么系统受扰动后正特征值方向上的误差会被指数放大，控制增益再高也很难压住。合理做法是在规划阶段引入稳定域约束，基于离线相平面图把鞍点邻域标记为禁区，相当于给控制器增加了一条安全边界。这个思想对于漂移控制、极限工况避障、冰雪路面驾驶辅助都有应用价值。
3. MATLAB仿真完整实现流程
3.1 车辆参数与轮胎模型选择
这里以一辆典型乘用车为例，所有参数只是量级参考，真要用到实车上，需要通过整车试验或轮胎辨识获得。
| 参数 | 符号 | 取值 | 单位 |
|---|---|---|---|
| 整车质量 | m | 1500 | kg |
| 横摆转动惯量 | Iz | 2500 | kg·m² |
| 质心到前轴距离 | a | 1.2 | m |
| 质心到后轴距离 | b | 1.4 | m |
| 前轮侧偏刚度 | Cf | 80000 | N/rad |
| 后轮侧偏刚度 | Cr | 120000 | N/rad |
| 路面附着系数 | μ | 0.8 | - |
| 纵向车速 | V | 20 | m/s |
| 前轮转角 | δ | 0.1 | rad |
前后轴垂直载荷按照静载分配计算：
Fzf = m * g * b / (a + b) Fzr = m * g * a / (a + b)轮胎力模型是这一步最重要的选择。如果直接用最简单的线性模型Fy = -C * α，那么给定一个δ，系统只有唯一平衡点，根本不会出现鞍点，也就没有临界轨迹可画。要出现多平衡点和鞍点，轮胎模型必须带有饱和甚至回落特性。我常用的Fiala简化模型是一个典型的工程选择，它在侧偏角较小时表现与线性轮胎一致，在侧偏角增大后侧偏力先饱和再回落，和真实轮胎在大侧偏角下的趋势接近。简化模型的目标不是百分之百复现真实轮胎曲线，而是保证相平面分析能展现出结构性的边界特征。
3.2 从状态方程到平衡点扫描
把前面推导写成MATLAB函数，最核心的就是状态方程和轮胎力两部分：
function dq = vehicle_dynamics(q, delta, p) beta = q(1); r = q(2); V = p.V; % 前后轮侧偏角 alpha_f = beta + p.a * r / V - delta; alpha_r = beta - p.b * r / V; % Fiala 简化轮胎模型 Fyf = tire_fiala(alpha_f, p.Cf, p.mu, p.Fzf); Fyr = tire_fiala(alpha_r, p.Cr, p.mu, p.Fzr); % 二自由度车身方程 dbeta = (Fyf + Fyr) / (p.m * V) - r; dr = (p.a * Fyf - p.b * Fyr) / p.Iz; dq = [dbeta; dr]; end function Fy = tire_fiala(alpha, C, mu, Fz) alpha_s = atan(3 * mu * Fz / C); if abs(alpha) < alpha_s Fy = -C * alpha ... + (C^3 / (3 * mu * Fz)) * abs(alpha) * alpha ... - (C^5 / (27 * (mu * Fz)^2)) * alpha^3; else Fy = -mu * Fz * sign(alpha); end end有了状态方程，下一步找平衡点。用fsolve对多个网格初值做扫描是一种直观做法：
beta_scan = linspace(-0.3, 0.3, 21); r_scan = linspace(-0.7, 0.7, 21); eq_points = []; for bi = beta_scan for ri = r_scan [q, fval, exitflag] = fsolve(@(q) vehicle_dynamics(q, delta, p), [bi, ri], opts); if exitflag > 0 && norm(fval) < 1e-6 eq_points = [eq_points; q(:)']; end end end % 按行去重，容差1e-4 eq_points = uniquetol(eq_points, 1e-4, 'ByRows', true);扫描初值的网格不能太稀。以β步长0.05、r步长0.1为例，一般来说够用；但碰到平衡点恰好落在网格之间时，容易漏掉。更稳的做法是结合等高线图：先用contour画出dβ/dt=0和dr/dt=0两条零线，从零线的交点数预计平衡点个数，再针对每个交点附近取初值做fsolve。我经常两种方法结合，先用contourc提取零线上的点云，再用点云作为初值。
3.3 鞍点识别与临界轨迹数值追踪
得到所有平衡点后，逐个计算雅可比矩阵。手推雅可比解析式容易出错，代码里用中心差分最省事：
function J = numerical_jacobian(q, delta, p) h = 1e-6; f0 = vehicle_dynamics(q, delta, p); J = zeros(2, 2); for i = 1:2 qp = q; qm = q; qp(i) = qp(i) + h; qm(i) = qm(i) - h; J(:, i) = (vehicle_dynamics(qp, delta, p) - vehicle_dynamics(qm, delta, p)) / (2 * h); end end筛选鞍点的逻辑就是找特征值一正一负的平衡点：
for i = 1:size(eq_points, 1) q = eq_points(i, :)'; J = numerical_jacobian(q, delta, p); lambda = real(eig(J)); if prod(lambda) < 0 saddle_point = q; break; end end筛选后要可视化确认：在相平面图上把每个平衡点点出来，稳定平衡点用圆圈，鞍点用三角。如果分类结果和矢量场走势对不上，先怀疑差分步长h是否合适，再检查轮胎力公式的符号。对于轮胎力这种带饱和的分段函数，h太小，函数值变化会被舍入噪声淹没，我常用1e-6到1e-5之间的值。
画分界线时，先求鞍点处雅可比矩阵的特征向量，取对应正特征值的特征向量u，归一化后得到局部不稳定方向。然后设一个很小偏移量eps0，分别以鞍点加减eps0*u为初值，调用ODE45对非线性状态方程积分：
[V, D] = eig(J); lambda = diag(D); iu = find(real(lambda) > 0); u = real(V(:, iu(1))); u = u / norm(u); eps0 = 1e-3; tspan = [0 20]; options = odeset('RelTol', 1e-6, 'AbsTol', 1e-8, 'Events', @(t, q) event_out_of_range(t, q)); [t1, y1] = ode45(@(t, q) vehicle_dynamics(q, delta, p), tspan, saddle_point + eps0 * u, options); [t2, y2] = ode45(@(t, q) vehicle_dynamics(q, delta, p), tspan, saddle_point - eps0 * u, options);积分初值的偏移量eps0很有讲究。太小了，ODE45在鞍点附近会因为数值误差而偏离正确流形；太大了，初值点已经落到分界线以外，积分出来的轨迹就是普通轨道，不是临界轨迹。经验上，对于β、r这种量级分别在0.1和0.5左右的系统，eps0取1e-3到1e-2比较合适。如果发现轨迹怪怪地绕圈，先减小eps0试试。
3.4 把相平面画得清楚可读
矢量场绘制我推荐用quiver加自归一化。网格上每个点的方向导数除以自身幅值，让所有箭头等长，图面干净，能清晰显示流线走向；但这样会丢失速度大小信息，所以可以再叠加一个颜色背景表示状态变化快慢。整理后的绘图代码大致是：
beta_grid = linspace(-0.3, 0.3, 30); r_grid = linspace(-0.8, 0.8, 30); [B, R] = meshgrid(beta_grid, r_grid); DB = zeros(size(B)); DR = zeros(size(B)); for i = 1:numel(B) dq = vehicle_dynamics([B(i); R(i)], delta, p); DB(i) = dq(1); DR(i) = dq(2); end % 归一化箭头长度 magnitude = sqrt(DB.^2 + DR.^2); DBn = DB ./ magnitude; DRn = DR ./ magnitude; figure; hold on; quiver(B, R, DBn, DRn, 0.5, 'k'); % 零线 contour(B, R, DB, [0 0], 'b', 'LineWidth', 1.5); contour(B, R, DR, [0 0], 'g', 'LineWidth', 1.5); % 平衡点与临界轨迹 plot(y1(:,1), y1(:,2), 'r-', 'LineWidth', 2); plot(y2(:,1), y2(:,2), 'r-', 'LineWidth', 2); plot(eq_stable(1), eq_stable(2), 'go', 'MarkerSize', 10); plot(saddle_point(1), saddle_point(2), 'r^', 'MarkerSize', 10); axis([-0.3 0.3 -0.8 0.8]); xlabel('质心侧偏角 \beta [rad]'); ylabel('横摆角速度 r [rad/s]');若嫌quiver参数太多，可以用MATLAB自带streamslice直接画流线，视觉效果更接近论文插图。不过streamslice对网格较敏感，建议先用粗网格试跑，确认矢量场方向没问题再加密。
临界轨迹画完之后一定要做验证，否则分界线位置对不对根本没人知道。验证方法很简单：在分界线两侧各取几个初始点，分别做时域仿真；如果一侧全部收敛到稳定平衡点，另一侧全部发散或进入另一个吸引域，说明临界轨迹基本可信。如果发现穿帮，优先检查轮胎模型参数还是积分方向选错了。
4. 相平面绘制中的常见问题与实用技巧
4.1 线性模型画不出鞍点，这是最大的坑
先说最大的坑。很多初学者上来直接套线性模型，结果无论怎么改车速和转角，相平面图永远只有一个稳定平衡点，找半天也看不到鞍点。原因很简单：线性轮胎模型下状态方程是线性的，平衡点唯一，不存在特征值一正一负的可能。鞍点本质上是非线性系统在多个平衡点时才会出现的边界点，所以必须用带饱和特性的非线性轮胎模型。如果只是想快速出图，把Fiala模型塞进去，立刻就能看到平衡点分裂、鞍点出现。
这里还要提醒一句：有些文章里的线性模型稳定分析用的是“横摆角速度增益”或者特征值实部，那是线性稳定性分析，和相平面上的鞍点不是一回事。线性稳定性分析能判断单一平衡点是否稳定，但画不出稳定域边界。
4.2 特征向量方向、初始偏移量与ODE积分稳定性
第二个常见问题是ODE45在临界轨迹附近疯狂报错，或者算出来的轨迹乱跳。原因往往是鞍点初始偏移方向拿错了，或者eps0取得过大。另外，用ODE45和ODE15s对同样初值在长时间积分后可能会有明显差异，尤其在鞍点附近这种数值敏感区域。画分界线时建议用带严格误差控制的ODE45，并把RelTol设成1e-6，AbsTol设成1e-8。不要为了图快省略误差控制参数。
还有一个问题是fsolve返回的平衡点不唯一时去重难。可以用uniquetol设置1e-5的容差，按行去重。要注意，不同初值可能收敛到同一点，也可能直接收敛到命令行窗口报错。需要在循环里判断exitflag是否大于0，并且检查残差norm(fval)是否足够小，否则很容易把不在零线上的点当作平衡点。
4.3 参数敏感性：车速、转角和附着系数如何改变稳定域
相平面稳定域随参数变化很敏感。我跑过的几个典型工况对比：
- 车速从20m/s提高到30m/s，稳定域明显向原点收缩，鞍点位置离原点更近，临界轨迹也更加弯曲。
- 前轮转角从0.05rad增加到0.15rad，稳定平衡点整体向β负方向移动，鞍点几乎贴近稳定平衡点，两者距离越小，说明可用稳定域越窄。
- 附着系数从0.8降到0.3，由于轮胎饱和点提前，分界线内缩幅度非常大，哪怕在中等车速下也可能出现大范围失稳。
这些趋势和实车感觉完全对得上，用来横向验证模型是否合理很有用。比如参数调完发现高速小转角下稳定域比低速大转角还大，那多半是轮胎力公式或载荷分配出问题了。
下面整理一个快速排查表，遇到问题可以对着查：
| 问题现象 | 可能原因 | 解决办法 |
|---|---|---|
| 相图上找不到鞍点 | 用了线性轮胎模型 | 换成饱和型非线性轮胎模型 |
| 平衡点扫描漏点 | fsolve初值太少或分布不均 | 先用零线交点判断个数，再扫多个初值 |
| 临界轨迹乱跳 | eps0太大、方向不对 | 减小eps0，重新检查特征向量方向 |
| 鞍点位置和流线不符 | 差分步长h不合适 | 调h到1e-6~1e-5之间重算雅可比 |
| 分界线两侧初始点时域仿真不收敛 | 模型参数或积分方向不对 | 验证初始点选取，核对状态方程符号 |
5. 从相平面到稳定性控制的实战思路
5.1 用分界线做稳定域评估
相平面图的最大价值就是能给出一个直觉的稳定域边界。工程上常把临界轨迹围成的区域面积作为稳定性指标：面积越大，车辆在当前工况下的稳定裕度越大。即使不做形式化的稳定性证明，离线跑几百组工况、画出不同δ和V下的稳定域，也能直接指导底盘调校，比如在某个车速区间ESP介入阈值应该放宽还是收紧。更进一步，可以把临界轨迹拟合成一个包络函数，在线判断当前β、r是否在边界内，作为触发稳定控制的开关信号。
5.2 控制器在线使用鞍点信息的思路
很多稳定性控制器的核心逻辑其实可以统一成一句话：保证状态点不越过临界轨迹，并在越过时施加干预让其回退。比如基于前馈加反馈的主动前轮转向、后轮转矩分配，控制目标可以是让状态点落到稳定平衡点附近，同时远离鞍点。由于临界轨迹是离线算好的，在线运行时只需要计算当前点到分界线的带符号距离，计算成本极低，完全可以在实时控制器里跑。这也是我喜欢先用相平面做预分析的原因：它把复杂的稳定性问题压缩成一张二维地图，方便控制逻辑设计和标定。
5.3 后续扩展方向
这套MATLAB仿真还可以继续扩展。轮胎模型升级为魔术公式后，临界轨迹会更精确，但参数标定工作量会上一个台阶；纵向速度V不再恒定，改成与纵向动力学耦合的联合仿真，就能研究制动或驱动过程中的稳定性切换。更进一步，把状态方程封装成S-Function，接入实时硬件环境做HIL测试，或者与车辆动力学软件联合验证控制效果，都是很成熟的工程路径。不过无论外围怎么扩展，相平面分析和鞍点识别的核心思路不会变。
我个人在实际操作中的体会是：画出一张漂亮的相平面图不难，难的是确保这张图经得起验证。每一次参数调整后，我都会在分界线两侧成对撒点做时域仿真，看轨迹是否真的按边界分开；只要有一次穿界，我就会重新检查轮胎力公式和数值积分设置，而不是直接改线去凑预期效果。最后再分享一个小技巧：在临界轨迹上多取几个点，把该点的状态导数方向画成短箭头，你会看到所有箭头都和轨迹相切，这个几何性质可以作为快速检查画得对不对的经验方法。
