做车辆稳定性仿真这些年相平面是绕不开的一个工具。尤其是当你手里没有CarSim、没有dSPACE那套昂贵设备只想用Matlab把二自由度车辆模型跑透的时候质心侧偏角-横摆角速度相平面就是最直观、最便宜的一扇窗。最近我在整理一套完整的仿真流程从模型搭建到鞍点识别再到临界轨迹绘制踩了不少坑也总结出一些可以复现的经验正好分享出来。这篇文章记录的是我用Matlab实现“二自由度车辆相平面分析 鞍点识别 临界轨迹分界线绘制”的完整过程。手把手讲清楚模型怎么建、鞍点怎么找、临界轨迹怎么画以及一路上踩过的坑。适合正在做车辆稳定性控制、底盘集成、ESP算法验证的朋友也适合刚接触相平面分析、想快速上手复现的初学者。核心关键词就是这六个字二自由度、相平面、鞍点、临界轨迹、Myelination——啊不好意思说顺嘴了是“临界轨迹、Matlab仿真”。如果你已经跑过二自由度模型的时序响应但总觉得“看曲线看不出名堂”这篇文章就是给你准备的。1. 项目要做的事为什么相平面能替你说清楚“稳定不稳定”1.1 相平面到底画的是什么相平面本质上是一个状态空间投影。对于二自由度车辆模型我们关心的两个状态就是质心侧偏角β单位rad和横摆角速度γ单位rad/s。以β为横轴、γ为纵轴把不同初始状态β₀, γ₀出发的系统运动轨迹全部画在一张图上就得到了相平面图。为什么看相平面比看时序曲线更爽因为时序响应回答的是“某个初始状态下系统怎么走”而相平面回答的是“所有初始状态下系统分别往哪里去”。车辆稳定性不是单一工况问题它依赖初始状态、车速、路面附着、前轮转角等因素。相平面把这种“状态依赖”直接可视化收敛到原点的区域是安全区远离原点的发散区就是失控区二者之间的边界——临界轨迹就是稳定域的几何边界。我之前做过一个ESP阈值标定项目单纯看横摆角速度时序判断车辆是否失稳很容易误报。后来把质心侧偏角-横摆角速度相平面画出来配合鞍点位置去定阈值误报率明显下降。因为相平面能提前告诉你从某个状态出发1秒后车辆是救得回来还是救不回来。1.2 二自由度模型够用且不复杂的首选二自由度车辆模型又叫“自行车模型”它把左右两侧车轮合并为一个等效车轮只保留两个自由度侧向运动和横摆运动。纵向速度假设恒定不考虑侧倾、俯仰、空气动力学。听起来简化得很粗暴但大量研究表明在侧向加速度不超过0.4g的普通行驶工况下二自由度模型的精度足够用于稳定性分析和控制设计。真正的车辆稳定性研究核心矛盾是轮胎侧偏力的非线性。二自由度模型加上非线性轮胎模型Fiala或魔术公式就能抓住车辆失稳的本质轮胎先在侧偏角较小时近似线性、在大侧偏角或低附着路面上饱和导致横摆力矩无法继续抵抗扰动。这个非线性行为在相平面上会呈现出鞍点、稳定域边界等丰富的结构。1.3 鞍点和临界轨迹在工程上意味着什么鞍点是相平面上最特殊的平衡点。它的雅可比矩阵特征值一正一负也就是说在这个点的邻域内一个方向的状态被吸引向它另一个方向的状态却背离它。工程上鞍点往往被解释为“临界失稳点”——状态一旦朝不稳定方向越过鞍点车辆就会进入不可控区域。临界轨迹则是从鞍点出发、将相平面划分为稳定域与不稳定域的分界线。它可以被直观地理解为“最后一道防线”初始状态落在分界线以内车辆状态会收敛到稳定平衡点落在分界线以外状态迅速发散、侧偏角急剧增大此时驾驶员很难通过转向操作救回车辆。这就是ESP、ESC等稳定性控制系统需要实时估计稳定域的理论基础。2. 模型与方程从自行车模型到相平面状态方程2.1 自由度与基本假设二自由度模型的设定如下车辆纵向速度Vx恒定不考虑纵向动力学。两个自由度质心侧偏角β或侧向速度Vy、横摆角速度γ。前轮转角δf直接作为系统输入不考虑转向系统惯量和延迟。左右车轮合并为等效前轮和等效后轮位于车辆中轴线上。忽略侧倾、俯仰、轮胎回正力矩、纵向力耦合。这些假设决定了模型的使用边界只适合分析前轮转角激励下的横向动力学和横摆响应不适合大纵向加减速、严重侧倾等极限工况。但在稳定性和相平面分析这个范畴内它是最经典、最高效的起点。2.2 侧向动力学方程的推导拿一个质量为m、横摆转动惯量为Iz的车辆模型前轴到质心距离Lf后轴到质心距离Lr。前轮侧偏力Fyf后轮侧偏力Fyr。侧向力平衡方程m·Vx·(β̇ γ) Fyf Fyr横摆力矩平衡方程Iz·γ̇ Lf·Fyf - Lr·Fyr第一个方程左边为什么是m·Vx·(β̇γ)因为质心侧偏角β表示速度方向相对车身纵轴的夹角当车辆既有横向速度变化β̇又有横摆运动γ时质心绝对加速度在车身横向的分量近似为Vx·(β̇γ)。这是车辆动力学中很基础的几何关系推导过程可以这样想车身以角速度γ旋转速度矢量方向也随之旋转旋转带来的横向加速度贡献恰好是Vx·γ再加上β自身的增长率贡献Vx·β̇。前后轮侧偏角的几何关系也很重要αf β (Lf·γ)/Vx - δf αr β - (Lr·γ)/Vx这个式子成立的物理背景是前轮中心的速度方向相对于车身纵轴有β Lf·γ/Vx的偏角车身旋转带来的附加侧向速度再减去转向角δf就得到轮胎实际运动方向与轮胎对称面之间的夹角——侧偏角。2.3 轮胎侧偏特性线性与非线性Fiala线性轮胎模型把轮胎侧偏力简化成Fy -C·α相当于一根线性弹簧。它能解释很多基础现象但有一个致命缺陷侧偏角增大后轮胎力并不会无限增大而是趋于饱和。饱和之后二自由度线性模型就完全失效了——这也是为什么很多初学者用线性模型画相平面发现永远只有原点一个平衡点看不到鞍点。我推荐用Fiala轮胎模型它参数少、实现简单但抓住了轮胎饱和这个关键非线性特性当|α| ≤ α_sl时 Fy -C·tan(α) (C²/(3·μ·Fz))·tan(α)·|tan(α)| - (C³/(27·(μ·Fz)²))·tan³(α)当|α| α_sl时 Fy -μ·Fz·sign(α)其中峰值侧偏角α_sl atan(3·μ·Fz/C)。别看这一堆多项式吓人它的物理意义很清晰第一项是线性侧偏力主导项第二、第三项是胎面局部滑移后对线性关系的修正最终让轮胎力在α接近α_sl时达到峰值μ·Fz之后保持饱和。与魔术公式相比Fiala模型只需要侧偏刚度C、垂直载荷Fz和附着系数μ三个参数在相平面分析这种定性定量结合的场合非常合适。2.4 参数取值一组常用的仿真标定下面这组参数取自某典型中型轿车的等效单轨模型也是我实际调试中用得比较顺的一组参数含义数值m整车质量1270 kgIz横摆转动惯量1536 kg·m²Lf质心到前轴距离1.016 mLr质心到后轴距离1.460 mCf前轮等效侧偏刚度44000 N/radCr后轮等效侧偏刚度66000 N/radμ路面附着系数0.85Vx纵向车速25 m/s90 km/hδf前轮转角0可调注意Cf小于Cr也就是前轴侧偏刚度小于后轴车辆呈现适当的不足转向特性。这个设定符合常见民用车的操纵倾向。如果反过来设成过度转向线性临界速度会很低相平面形态会完全不同这一点后面在实操部分还会展开。前后轴垂直载荷按静载荷分配计算 Fzf m·g·Lr/(LfLr)Fzr m·g·Lf/(LfLr)带入数值大致是Fzf ≈ 7343 NFzr ≈ 5107 N。这两个值直接影响轮胎饱和点非常重要。3. 鞍点和临界轨迹的原理与数值求法3.1 平衡点求解不止一个解相平面分析的第一步是找到系统所有平衡点也就是满足状态方程右侧同时为零的点Fyf Fyr m·Vx·γ Lf·Fyf - Lr·Fyr 0对于线性轮胎模型这个方程组只有一个解直线行驶工况下就是原点。但引入非线性轮胎模型之后状况完全不同轮胎力在饱和区呈现出“非单调”但饱和的特性大侧偏角组合下也会出现满足力矩平衡的状态点。于是系统可能有多个平衡点——除了原点附近那个稳定平衡点还存在一到两个鞍点。具体有几个平衡点、它们在哪里取决于Vx、δf、μ和刚度参数。我的经验是高速20m/s以上、低附着μ0.5或大转角工况下鞍点出现的概率显著增加。所以如果你跑低速工况发现只有一个平衡点别急着怀疑代码先调高速试试。数值上我不会只依赖fsolve从一个初始猜测出发因为多解方程非常容易漏根。更稳妥的做法是在β-γ平面上做一次粗网格扫描在每个网格点计算状态方程的值找出符号发生变化的区域作为潜在根区域然后在这些区域附近用fsolve精化。这样可以最大概率拿到所有平衡点。3.2 平衡点分类稳定焦点、鞍点与不稳定结点找到平衡点之后要在该点做线性化得到2×2雅可比矩阵J。然后看J的特征值λ1、λ2两个特征值实部都为负 → 稳定焦点或稳定结点状态收敛两个特征值实部都为正 → 不稳定结点或焦点状态发散一个实部为正、一个实部为负 → 鞍点特征值为共轭复数 → 焦点实部符号决定稳定与否。在车辆相平面分析中最常见的组合是一个稳定焦点通常就在原点或稳态转向点附近加一个或两个鞍点。稳定焦点对应车辆线性稳定工作点鞍点对应失稳边界上的临界状态。判断特征值只是第一步真正有用的是鞍点处对应的特征向量。特征向量给出了局部稳定流形和不稳定流形的方向稳定特征向量方向指向上游状态被吸引到鞍点不稳定特征向量方向指向下游状态被鞍点推向远方。这两个方向就是我们绘制临界轨迹的出发参考。3.3 临界轨迹的数值积分思路临界轨迹数学上叫鞍点的稳定流形和不稳定流形的组合。在工程相平面图上我们通常重点绘制不稳定流形因为它是稳定域边界的主要组成部分。数值求法的基本原理是从鞍点出发沿着不稳定特征向量的方向偏移一个很小的距离ε然后把这个偏移点作为初始点用ode45正向积分。积分得到的轨迹就会沿着不稳定流形向后延伸最终成为划分稳定域与失稳域的分界线。注意一个关键细节不能直接从鞍点本身出发积分。因为鞍点是一个平衡点从它本身出发的数值积分会一直停留在那里根本不会沿着流形移动。偏移距离ε的选择是个经验活太小积分器会在鞍点附近徘徊很久才“找到”流形方向太大初始点已经偏离真实流形轨迹会在后续积分中明显失真。我的经验是取状态幅值的0.1%~0.5%对于β量级0.1~0.3 rad的系统ε1e-3左右比较合适。同样稳定流形的绘制需要对系统时间反向积分也就是tspan设为[0 -5]。这里有个数值按摩反向积分时原来的稳定特征值会变成不稳定数值误差会被指数放大非常容易发散。所以绘制稳定流形时通常需要更严的误差容差或者分段推进。如果只为了工程上界定稳定域可以只绘制不稳定流形分支效果也够用。3.4 稳定域边界工程上的“可救区”鞍点的不稳定流形在相平面上画出来后会形成一个很直观的几何结构两条临界轨迹从鞍点出发向两侧延伸把它们所在区域的流动“分隔”开来。落在边界内侧的初始状态最终收敛到稳定平衡点落在边界外侧的初始状态β或γ持续增大车辆进入不可控状态。这个边界在工程上的意义就是“可救区”。ESP系统如果真的把车辆状态点跟踪到相平面上一旦发现状态点逼近临界轨迹就提前施加横摆力矩控制比单纯依靠阈值判断侧偏角大小要提前得多、可靠得多。这也是我为什么坚持做鞍点临界轨迹分析——它能直接把一个“是否可救”的定性结论画成一张几何边界图。4. Matlab仿真全流程实现4.1 总体流程与代码架构整个仿真分五步定义参数和状态方程 → 网格化初始状态并积分得到相平面流线 → 搜索平衡点 → 识别鞍点 → 从鞍点出发绘制临界轨迹。整体代码组织成一个主脚本加三个局部函数。Matlab R2016b以后支持在脚本末尾直接定义局部函数所以可以放在一个.m文件里直接跑。完整脚本可以分成几段我这里逐步讲解最后你可以拼起来运行。先说代码架构clearvars; clc; close all; % 车辆参数结构体 p.m 1270; p.Iz 1536; p.Lf 1.016; p.Lr 1.460; p.Cf 44000; p.Cr 66000; p.Vx 25; p.mu 0.85; p.delta_f 0; p.g 9.81; p.Fzf p.m * p.g * p.Lr / (p.Lf p.Lr); p.Fzr p.m * p.g * p.Lf / (p.Lf p.Lr);参数结构体p是所有函数共享的这样后续改参数非常方便不用到处改。4.2 第一步定义状态方程和轮胎力函数状态方程函数接收时间t虽然方程不显含时间但ode45要求第一个参数是t、状态向量x和参数结构体p。x的约定是x(1)βx(2)γ。function dx vehicle_state(t, x, p) beta x(1); gamma x(2); alpha_f beta p.Lf * gamma / p.Vx - p.delta_f; alpha_r beta - p.Lr * gamma / p.Vx; Fyf tire_force(alpha_f, p.Cf, p.mu, p.Fzf); Fyr tire_force(alpha_r, p.Cr, p.mu, p.Fzr); dx zeros(2, 1); dx(1) (Fyf Fyr) / (p.m * p.Vx) - gamma; dx(2) (p.Lf * Fyf - p.Lr * Fyr) / p.Iz; end轮胎力函数实现Fiala模型function Fy tire_force(alpha, C, mu, Fz) alpha_sl atan(3 * mu * Fz / C); ta tan(alpha); if abs(alpha) alpha_sl Fy -C * ta (C^2 / (3 * mu * Fz)) * ta * abs(ta) ... - (C^3 / (27 * (mu * Fz)^2)) * ta^3; else Fy -mu * Fz * sign(alpha); end end注意Matlab的所有三角函数都默认使用弧度别在别处换算成角度再喂进来后面绘图标签也要写清楚单位。4.3 第二步绘制基础相平面在β∈[-0.25, 0.3]、γ∈[-1.0, 1.0]的范围内取16×16个网格点作为初始状态从每个点出发用ode45积分4秒把轨迹画出来。beta_span [-0.25, 0.3]; gamma_span [-1.0, 1.0]; n_beta 16; n_gamma 16; betas linspace(beta_span(1), beta_span(2), n_beta); gammas linspace(gamma_span(1), gamma_span(2), n_gamma); T_end 4; options odeset(RelTol, 1e-6, AbsTol, 1e-8); figure(Color, w); hold on; grid on; box on; for i 1:n_beta for j 1:n_gamma x0 [betas(i); gammas(j)]; [~, x] ode45((t, x) vehicle_state(t, x, p), [0 T_end], x0, options); plot(x(:, 1), x(:, 2), LineWidth, 0.5, Color, [0.7, 0.7, 0.9]); end end xlabel(质心侧偏角 \beta [rad], FontSize, 12); ylabel(横摆角速度 \gamma [rad/s], FontSize, 12); title(sprintf(二自由度车辆相平面 Vx%.0f m/s, \\delta_f%.2f rad, p.Vx, p.delta_f));网格密度我建议16×16到20×20之间。太密了轨迹线重叠严重整个图变成一团浆糊太疏了看不清流场的整体走向。轨迹颜色统一用浅蓝细线平衡点和临界轨迹用深色粗线覆盖视觉层次最好。每次网格仿真要执行256次ode45积分我实测在普通笔记本上大概需要几十秒到一两分钟建议加个进度提示比如每隔若干格disp一下免得干等。4.4 第三步搜索平衡点并识别鞍点平衡点求解我提供两套思路。第一套是快速脚本方案直接用fsolve搭配多组初始猜测。第二套是稳稳妥妥的网格扫描方案先粗定位再精化。下面给出第一套实现init_guesses [0, 0; 0.1, 0.3; -0.1, -0.3; 0.2, 0.6; -0.2, -0.6; ... 0.25, 0.9; -0.25, -0.9; 0.05, 0.15; -0.05, -0.15]; opts optimoptions(fsolve, Display, off, ... FunctionTolerance, 1e-12, ... OptimalityTolerance, 1e-12); equilibria []; for k 1:size(init_guesses, 1) x0 init_guesses(k, :); [xeq, ~, exitflag] fsolve((x) vehicle_state(0, x, p), x0, opts); if exitflag 0 xeq xeq; if isempty(equilibria) equilibria xeq; elseif min(vecnorm(equilibria - xeq, 2, 2)) 1e-4 equilibria [equilibria; xeq]; end end end得到平衡点集合后对每个平衡点计算数值雅可比矩阵并求特征值判断鞍点saddle_point []; for k 1:size(equilibria, 1) xeq equilibria(k, :); J numerical_jacobian((x) vehicle_state(0, x, p), xeq); lambda eig(J); fprintf(平衡点(%.4f, %.4f) 特征值: %.3f, %.3f\n, ... xeq(1), xeq(2), real(lambda(1)), real(lambda(2))); if real(lambda(1)) * real(lambda(2)) 0 saddle_point xeq; end end数值雅可比用中心差分实现比起直接调用符号工具箱速度快且没有解析求导的烦恼function J numerical_jacobian(fun, x0) h 1e-6; n length(x0); J zeros(n, n); for i 1:n xp x0; xp(i) xp(i) h; xm x0; xm(i) xm(i) - h; J(:, i) (fun(xp) - fun(xm)) / (2 * h); end end中心差分步长h选1e-6比较合适。太大截断误差明显太小浮点舍入误差主导。1e-6在大多数车辆动力学方程尺度下都很安全。如果你用网格扫描方式粗定位大致是这样把β-γ平面分成50×50网格计算每个网格点的状态变化率dx找dx(1)0与dx(2)0的等值线交点然后用这些交点作fsolve的初始猜测。这个方法能抓到的平衡点数量通常比固定猜测多尤其在鞍点位置偏离常规认知时。4.5 第四步从鞍点出发绘制临界轨迹识别出鞍点后计算该点的雅可比矩阵找到正实部特征值对应的特征向量作为不稳定流形方向。从鞍点偏移±ε·v后正向积分得到两条临界轨迹分支。如果是直线行驶对称工况你会发现两条分支几乎对称如果加了前轮转角两条分支就是非对称的很有意思。if ~isempty(saddle_point) J_s numerical_jacobian((x) vehicle_state(0, x, p), saddle_point); [V, D] eig(J_s); lambda diag(D); [~, idx_plus] max(real(lambda)); v_unstable V(:, idx_plus); v_unstable v_unstable / norm(v_unstable); eps_offset 1e-3; t_span_fwd [0 5]; options_strict odeset(RelTol, 1e-8, AbsTol, 1e-10); x_start1 saddle_point eps_offset * v_unstable; x_start2 saddle_point - eps_offset * v_unstable; [~, x_branch1] ode45((t, x) vehicle_state(t, x, p), t_span_fwd, x_start1, options_strict); [~, x_branch2] ode45((t, x) vehicle_state(t, x, p), t_span_fwd, x_start2, options_strict); plot(x_branch1(:, 1), x_branch1(:, 2), r-, LineWidth, 2.2); plot(x_branch2(:, 1), x_branch2(:, 2), r-, LineWidth, 2.2); plot(saddle_point(1), saddle_point(2), rp, MarkerSize, 15, ... MarkerFaceColor, r, MarkerEdgeColor, k); end我特别强调一下积分容差。临界轨迹的绘制对数值误差极敏感尤其是从鞍点出发这种极度不稳定的流形。RelTol我直接用1e-8而不是默认的1e-3误差容差设到1e-10。别嫌慢这条轨迹的稳定性直接决定了你画出来的边界是光滑曲线还是锯齿线。如果你还想绘制稳定流形即从鞍点反向积分得到的另一部分边界方法类似但tspan用[0 -5]初始偏移沿稳定特征向量方向。前面说过反向积分数值上更容易发散建议每隔一步就检查一下轨迹是否已经冲到离谱范围必要时用分段积分把上一段的终点作为下一段的起点摸着石头过河。4.6 第五步结果可视化与解读最终出图时我习惯把相平面流线、平衡点、鞍点、临界轨迹四类信息叠在一张图上。为了好读我还会加图例标注清楚哪些是流线、哪些是平衡点、哪些是临界轨迹。如果同一张图里有多组速度或转角工况用subplot排开逐张对比。一个等效可视化的技巧除了画轨迹线还可以对网格上的每个初始状态点做“最终收敛分类”——看它落在哪个平衡点的吸引域内然后用散点颜色填充稳定域和不稳定域。这种“相平面热图”比线条更直观尤其做汇报展示的时候特别好用。实现起来也不复杂就是在原有的16×16网格循环里判断最后时刻状态离哪个平衡点最近再用scatter打点。5. 实操中的常见问题与避坑经验5.1 相平面杂乱难辨问题出在网格和线宽很多新手第一次画出相平面惊呼“怎么全是一团线”。这通常不是代码错而是网格太密、线宽太宽、颜色太深。16×16网格配0.5磅线宽属于比较能看的组合。如果你想更清爽可以把网格降到12×12或者把每条轨迹的透明度调低Color第四个分量。还有一个我自己常用的技巧不画完整4秒轨迹只画后2秒的轨迹片段。因为相平面图里初始段的状态快速移动线条密集看了容易头晕而末段的轨迹慢慢接近吸引子或发散方向信息量反而最集中。5.2 找不到鞍点线性轮胎模型的局限如果你用线性轮胎模型Fy-C·α那整套分析会非常“干净”——干净到没有鞍点可画。线性系统只可能有一个平衡点稳定就是不稳定的收敛到它不稳定就是发散而没有任何边界。所以相平面里出现不了鞍点是因为你还没把非线性轮胎模型引进来。切换到Fiala模型后鞍点会从“无到有”。但注意不是所有参数下都有鞍点。低速、高附着、小转角时整个相平面像一个巨大的漩涡所有轨迹都向原点收敛鞍点消失。你要做的就是提高车速、降低附着系数或者加大前轮转角把失稳边界逼出来。5.3 临界轨迹“跑偏”积分设置与初始扰动的平衡临界轨迹跑偏我前前后后遇到过三种原因。第一种是偏移量ε太大初始点没有落在真实不稳定流形上。结果画出来的“临界轨迹”是条随机的发散曲线完全不能当边界用。这时把ε从1e-3降到1e-4重新跑。第二种是积分容差太松。默认RelTol是1e-3对普通轨迹没问题对鞍点流形就差太远了。我建议至少设到1e-6画分界线时直接上1e-8。第三种是积分时间过长。临界轨迹从鞍点出发后快速延伸如果T_end太大轨迹早就跑到β超过0.5甚至1.0的远场去了图上看起来就是一条冲出边界的线。把T_end限制在3~5秒看它把稳定域围起来就停。5.4 参数敏感换一组车辆参数后的注意事项车辆参数变化会显著影响相平面形态我吃过亏所以特别提醒三点。第一前后轴刚度比例变化会直接改变鞍点位置和临界轨迹的倾斜方向。Cf和Cr互换一下临界轨迹可能从“左上右下”变成“右上左下”。第二速度Vx的影响最大。Vx从15m/s变到30m/s稳定域可能从几乎覆盖整个图面缩小到只有原点周围一小圈。扫速度时我建议每次增加5m/s连续画3~5张图你就能很清楚地看到稳定域收缩的过程。第三μ附着系数控制着轮胎饱和点。μ越低饱和点越早出现鞍点越靠近原点。雨雪路面μ0.3左右的稳定域比干燥路面缩小一半以上这是ESP算法最关注的一个场景。5.5 平衡点搜索的漏解问题fsolve多初始猜测方案虽然简单但确实容易漏根。如果你发现相平面上某个区域的轨迹明显“绕着某个点转”但平衡点列表里没有那个点那大概率就是漏解了。我的改进方法前面提过的网格扫描法先用50×50网格计算dx(1)、dx(2)的符号场找它们分别从正变负的等值线交点把这些交点作为fsolve的初始猜测。这样基本能保证把所有平衡点一网打尽。条件允许的话还可以用pde或者continuation方法比如对Vx做参数扫描从已知平衡点逐步追踪来跟踪鞍点随车速的变化轨迹但那是另一个层面的工作一般项目用不着。还有个小细节如果系统有多个很接近的平衡点去重阈值不能设得太大。我前面代码里用1e-4够用。如果设成1e-2可能把物理上不同的平衡点错误合并导致信息丢失。一点个人经验整套流程跑通之后我最深的体会是相平面看起来只是一堆轨迹线但配合鞍点和临界轨迹它其实把车辆的稳定裕度直接画在了你面前。每次调ESP控制阈值、评估某个底盘参数改动带来的稳定性变化我都会先拉一组相平面看看特定工况下稳定域的形态心里先有个底再去做控制参数整定。最后分享一个我一直在用的小技巧不要只画单一速度下的相平面用一个subplot把15、20、25、30 m/s四个速度的相平面排在一起对比。你很快就会发现临界轨迹的“收拢”过程特别直观——稳定域从开阔地变成一条窄走廊车辆的安全余量到底有多大一目了然。这个方法我推荐给每一个做车辆稳定性的朋友比纯看时序响应靠谱多了。 二自由度车辆相平面分析：Matlab鞍点识别与临界轨迹绘制
做车辆稳定性仿真这些年相平面是绕不开的一个工具。尤其是当你手里没有CarSim、没有dSPACE那套昂贵设备只想用Matlab把二自由度车辆模型跑透的时候质心侧偏角-横摆角速度相平面就是最直观、最便宜的一扇窗。最近我在整理一套完整的仿真流程从模型搭建到鞍点识别再到临界轨迹绘制踩了不少坑也总结出一些可以复现的经验正好分享出来。这篇文章记录的是我用Matlab实现“二自由度车辆相平面分析 鞍点识别 临界轨迹分界线绘制”的完整过程。手把手讲清楚模型怎么建、鞍点怎么找、临界轨迹怎么画以及一路上踩过的坑。适合正在做车辆稳定性控制、底盘集成、ESP算法验证的朋友也适合刚接触相平面分析、想快速上手复现的初学者。核心关键词就是这六个字二自由度、相平面、鞍点、临界轨迹、Myelination——啊不好意思说顺嘴了是“临界轨迹、Matlab仿真”。如果你已经跑过二自由度模型的时序响应但总觉得“看曲线看不出名堂”这篇文章就是给你准备的。1. 项目要做的事为什么相平面能替你说清楚“稳定不稳定”1.1 相平面到底画的是什么相平面本质上是一个状态空间投影。对于二自由度车辆模型我们关心的两个状态就是质心侧偏角β单位rad和横摆角速度γ单位rad/s。以β为横轴、γ为纵轴把不同初始状态β₀, γ₀出发的系统运动轨迹全部画在一张图上就得到了相平面图。为什么看相平面比看时序曲线更爽因为时序响应回答的是“某个初始状态下系统怎么走”而相平面回答的是“所有初始状态下系统分别往哪里去”。车辆稳定性不是单一工况问题它依赖初始状态、车速、路面附着、前轮转角等因素。相平面把这种“状态依赖”直接可视化收敛到原点的区域是安全区远离原点的发散区就是失控区二者之间的边界——临界轨迹就是稳定域的几何边界。我之前做过一个ESP阈值标定项目单纯看横摆角速度时序判断车辆是否失稳很容易误报。后来把质心侧偏角-横摆角速度相平面画出来配合鞍点位置去定阈值误报率明显下降。因为相平面能提前告诉你从某个状态出发1秒后车辆是救得回来还是救不回来。1.2 二自由度模型够用且不复杂的首选二自由度车辆模型又叫“自行车模型”它把左右两侧车轮合并为一个等效车轮只保留两个自由度侧向运动和横摆运动。纵向速度假设恒定不考虑侧倾、俯仰、空气动力学。听起来简化得很粗暴但大量研究表明在侧向加速度不超过0.4g的普通行驶工况下二自由度模型的精度足够用于稳定性分析和控制设计。真正的车辆稳定性研究核心矛盾是轮胎侧偏力的非线性。二自由度模型加上非线性轮胎模型Fiala或魔术公式就能抓住车辆失稳的本质轮胎先在侧偏角较小时近似线性、在大侧偏角或低附着路面上饱和导致横摆力矩无法继续抵抗扰动。这个非线性行为在相平面上会呈现出鞍点、稳定域边界等丰富的结构。1.3 鞍点和临界轨迹在工程上意味着什么鞍点是相平面上最特殊的平衡点。它的雅可比矩阵特征值一正一负也就是说在这个点的邻域内一个方向的状态被吸引向它另一个方向的状态却背离它。工程上鞍点往往被解释为“临界失稳点”——状态一旦朝不稳定方向越过鞍点车辆就会进入不可控区域。临界轨迹则是从鞍点出发、将相平面划分为稳定域与不稳定域的分界线。它可以被直观地理解为“最后一道防线”初始状态落在分界线以内车辆状态会收敛到稳定平衡点落在分界线以外状态迅速发散、侧偏角急剧增大此时驾驶员很难通过转向操作救回车辆。这就是ESP、ESC等稳定性控制系统需要实时估计稳定域的理论基础。2. 模型与方程从自行车模型到相平面状态方程2.1 自由度与基本假设二自由度模型的设定如下车辆纵向速度Vx恒定不考虑纵向动力学。两个自由度质心侧偏角β或侧向速度Vy、横摆角速度γ。前轮转角δf直接作为系统输入不考虑转向系统惯量和延迟。左右车轮合并为等效前轮和等效后轮位于车辆中轴线上。忽略侧倾、俯仰、轮胎回正力矩、纵向力耦合。这些假设决定了模型的使用边界只适合分析前轮转角激励下的横向动力学和横摆响应不适合大纵向加减速、严重侧倾等极限工况。但在稳定性和相平面分析这个范畴内它是最经典、最高效的起点。2.2 侧向动力学方程的推导拿一个质量为m、横摆转动惯量为Iz的车辆模型前轴到质心距离Lf后轴到质心距离Lr。前轮侧偏力Fyf后轮侧偏力Fyr。侧向力平衡方程m·Vx·(β̇ γ) Fyf Fyr横摆力矩平衡方程Iz·γ̇ Lf·Fyf - Lr·Fyr第一个方程左边为什么是m·Vx·(β̇γ)因为质心侧偏角β表示速度方向相对车身纵轴的夹角当车辆既有横向速度变化β̇又有横摆运动γ时质心绝对加速度在车身横向的分量近似为Vx·(β̇γ)。这是车辆动力学中很基础的几何关系推导过程可以这样想车身以角速度γ旋转速度矢量方向也随之旋转旋转带来的横向加速度贡献恰好是Vx·γ再加上β自身的增长率贡献Vx·β̇。前后轮侧偏角的几何关系也很重要αf β (Lf·γ)/Vx - δf αr β - (Lr·γ)/Vx这个式子成立的物理背景是前轮中心的速度方向相对于车身纵轴有β Lf·γ/Vx的偏角车身旋转带来的附加侧向速度再减去转向角δf就得到轮胎实际运动方向与轮胎对称面之间的夹角——侧偏角。2.3 轮胎侧偏特性线性与非线性Fiala线性轮胎模型把轮胎侧偏力简化成Fy -C·α相当于一根线性弹簧。它能解释很多基础现象但有一个致命缺陷侧偏角增大后轮胎力并不会无限增大而是趋于饱和。饱和之后二自由度线性模型就完全失效了——这也是为什么很多初学者用线性模型画相平面发现永远只有原点一个平衡点看不到鞍点。我推荐用Fiala轮胎模型它参数少、实现简单但抓住了轮胎饱和这个关键非线性特性当|α| ≤ α_sl时 Fy -C·tan(α) (C²/(3·μ·Fz))·tan(α)·|tan(α)| - (C³/(27·(μ·Fz)²))·tan³(α)当|α| α_sl时 Fy -μ·Fz·sign(α)其中峰值侧偏角α_sl atan(3·μ·Fz/C)。别看这一堆多项式吓人它的物理意义很清晰第一项是线性侧偏力主导项第二、第三项是胎面局部滑移后对线性关系的修正最终让轮胎力在α接近α_sl时达到峰值μ·Fz之后保持饱和。与魔术公式相比Fiala模型只需要侧偏刚度C、垂直载荷Fz和附着系数μ三个参数在相平面分析这种定性定量结合的场合非常合适。2.4 参数取值一组常用的仿真标定下面这组参数取自某典型中型轿车的等效单轨模型也是我实际调试中用得比较顺的一组参数含义数值m整车质量1270 kgIz横摆转动惯量1536 kg·m²Lf质心到前轴距离1.016 mLr质心到后轴距离1.460 mCf前轮等效侧偏刚度44000 N/radCr后轮等效侧偏刚度66000 N/radμ路面附着系数0.85Vx纵向车速25 m/s90 km/hδf前轮转角0可调注意Cf小于Cr也就是前轴侧偏刚度小于后轴车辆呈现适当的不足转向特性。这个设定符合常见民用车的操纵倾向。如果反过来设成过度转向线性临界速度会很低相平面形态会完全不同这一点后面在实操部分还会展开。前后轴垂直载荷按静载荷分配计算 Fzf m·g·Lr/(LfLr)Fzr m·g·Lf/(LfLr)带入数值大致是Fzf ≈ 7343 NFzr ≈ 5107 N。这两个值直接影响轮胎饱和点非常重要。3. 鞍点和临界轨迹的原理与数值求法3.1 平衡点求解不止一个解相平面分析的第一步是找到系统所有平衡点也就是满足状态方程右侧同时为零的点Fyf Fyr m·Vx·γ Lf·Fyf - Lr·Fyr 0对于线性轮胎模型这个方程组只有一个解直线行驶工况下就是原点。但引入非线性轮胎模型之后状况完全不同轮胎力在饱和区呈现出“非单调”但饱和的特性大侧偏角组合下也会出现满足力矩平衡的状态点。于是系统可能有多个平衡点——除了原点附近那个稳定平衡点还存在一到两个鞍点。具体有几个平衡点、它们在哪里取决于Vx、δf、μ和刚度参数。我的经验是高速20m/s以上、低附着μ0.5或大转角工况下鞍点出现的概率显著增加。所以如果你跑低速工况发现只有一个平衡点别急着怀疑代码先调高速试试。数值上我不会只依赖fsolve从一个初始猜测出发因为多解方程非常容易漏根。更稳妥的做法是在β-γ平面上做一次粗网格扫描在每个网格点计算状态方程的值找出符号发生变化的区域作为潜在根区域然后在这些区域附近用fsolve精化。这样可以最大概率拿到所有平衡点。3.2 平衡点分类稳定焦点、鞍点与不稳定结点找到平衡点之后要在该点做线性化得到2×2雅可比矩阵J。然后看J的特征值λ1、λ2两个特征值实部都为负 → 稳定焦点或稳定结点状态收敛两个特征值实部都为正 → 不稳定结点或焦点状态发散一个实部为正、一个实部为负 → 鞍点特征值为共轭复数 → 焦点实部符号决定稳定与否。在车辆相平面分析中最常见的组合是一个稳定焦点通常就在原点或稳态转向点附近加一个或两个鞍点。稳定焦点对应车辆线性稳定工作点鞍点对应失稳边界上的临界状态。判断特征值只是第一步真正有用的是鞍点处对应的特征向量。特征向量给出了局部稳定流形和不稳定流形的方向稳定特征向量方向指向上游状态被吸引到鞍点不稳定特征向量方向指向下游状态被鞍点推向远方。这两个方向就是我们绘制临界轨迹的出发参考。3.3 临界轨迹的数值积分思路临界轨迹数学上叫鞍点的稳定流形和不稳定流形的组合。在工程相平面图上我们通常重点绘制不稳定流形因为它是稳定域边界的主要组成部分。数值求法的基本原理是从鞍点出发沿着不稳定特征向量的方向偏移一个很小的距离ε然后把这个偏移点作为初始点用ode45正向积分。积分得到的轨迹就会沿着不稳定流形向后延伸最终成为划分稳定域与失稳域的分界线。注意一个关键细节不能直接从鞍点本身出发积分。因为鞍点是一个平衡点从它本身出发的数值积分会一直停留在那里根本不会沿着流形移动。偏移距离ε的选择是个经验活太小积分器会在鞍点附近徘徊很久才“找到”流形方向太大初始点已经偏离真实流形轨迹会在后续积分中明显失真。我的经验是取状态幅值的0.1%~0.5%对于β量级0.1~0.3 rad的系统ε1e-3左右比较合适。同样稳定流形的绘制需要对系统时间反向积分也就是tspan设为[0 -5]。这里有个数值按摩反向积分时原来的稳定特征值会变成不稳定数值误差会被指数放大非常容易发散。所以绘制稳定流形时通常需要更严的误差容差或者分段推进。如果只为了工程上界定稳定域可以只绘制不稳定流形分支效果也够用。3.4 稳定域边界工程上的“可救区”鞍点的不稳定流形在相平面上画出来后会形成一个很直观的几何结构两条临界轨迹从鞍点出发向两侧延伸把它们所在区域的流动“分隔”开来。落在边界内侧的初始状态最终收敛到稳定平衡点落在边界外侧的初始状态β或γ持续增大车辆进入不可控状态。这个边界在工程上的意义就是“可救区”。ESP系统如果真的把车辆状态点跟踪到相平面上一旦发现状态点逼近临界轨迹就提前施加横摆力矩控制比单纯依靠阈值判断侧偏角大小要提前得多、可靠得多。这也是我为什么坚持做鞍点临界轨迹分析——它能直接把一个“是否可救”的定性结论画成一张几何边界图。4. Matlab仿真全流程实现4.1 总体流程与代码架构整个仿真分五步定义参数和状态方程 → 网格化初始状态并积分得到相平面流线 → 搜索平衡点 → 识别鞍点 → 从鞍点出发绘制临界轨迹。整体代码组织成一个主脚本加三个局部函数。Matlab R2016b以后支持在脚本末尾直接定义局部函数所以可以放在一个.m文件里直接跑。完整脚本可以分成几段我这里逐步讲解最后你可以拼起来运行。先说代码架构clearvars; clc; close all; % 车辆参数结构体 p.m 1270; p.Iz 1536; p.Lf 1.016; p.Lr 1.460; p.Cf 44000; p.Cr 66000; p.Vx 25; p.mu 0.85; p.delta_f 0; p.g 9.81; p.Fzf p.m * p.g * p.Lr / (p.Lf p.Lr); p.Fzr p.m * p.g * p.Lf / (p.Lf p.Lr);参数结构体p是所有函数共享的这样后续改参数非常方便不用到处改。4.2 第一步定义状态方程和轮胎力函数状态方程函数接收时间t虽然方程不显含时间但ode45要求第一个参数是t、状态向量x和参数结构体p。x的约定是x(1)βx(2)γ。function dx vehicle_state(t, x, p) beta x(1); gamma x(2); alpha_f beta p.Lf * gamma / p.Vx - p.delta_f; alpha_r beta - p.Lr * gamma / p.Vx; Fyf tire_force(alpha_f, p.Cf, p.mu, p.Fzf); Fyr tire_force(alpha_r, p.Cr, p.mu, p.Fzr); dx zeros(2, 1); dx(1) (Fyf Fyr) / (p.m * p.Vx) - gamma; dx(2) (p.Lf * Fyf - p.Lr * Fyr) / p.Iz; end轮胎力函数实现Fiala模型function Fy tire_force(alpha, C, mu, Fz) alpha_sl atan(3 * mu * Fz / C); ta tan(alpha); if abs(alpha) alpha_sl Fy -C * ta (C^2 / (3 * mu * Fz)) * ta * abs(ta) ... - (C^3 / (27 * (mu * Fz)^2)) * ta^3; else Fy -mu * Fz * sign(alpha); end end注意Matlab的所有三角函数都默认使用弧度别在别处换算成角度再喂进来后面绘图标签也要写清楚单位。4.3 第二步绘制基础相平面在β∈[-0.25, 0.3]、γ∈[-1.0, 1.0]的范围内取16×16个网格点作为初始状态从每个点出发用ode45积分4秒把轨迹画出来。beta_span [-0.25, 0.3]; gamma_span [-1.0, 1.0]; n_beta 16; n_gamma 16; betas linspace(beta_span(1), beta_span(2), n_beta); gammas linspace(gamma_span(1), gamma_span(2), n_gamma); T_end 4; options odeset(RelTol, 1e-6, AbsTol, 1e-8); figure(Color, w); hold on; grid on; box on; for i 1:n_beta for j 1:n_gamma x0 [betas(i); gammas(j)]; [~, x] ode45((t, x) vehicle_state(t, x, p), [0 T_end], x0, options); plot(x(:, 1), x(:, 2), LineWidth, 0.5, Color, [0.7, 0.7, 0.9]); end end xlabel(质心侧偏角 \beta [rad], FontSize, 12); ylabel(横摆角速度 \gamma [rad/s], FontSize, 12); title(sprintf(二自由度车辆相平面 Vx%.0f m/s, \\delta_f%.2f rad, p.Vx, p.delta_f));网格密度我建议16×16到20×20之间。太密了轨迹线重叠严重整个图变成一团浆糊太疏了看不清流场的整体走向。轨迹颜色统一用浅蓝细线平衡点和临界轨迹用深色粗线覆盖视觉层次最好。每次网格仿真要执行256次ode45积分我实测在普通笔记本上大概需要几十秒到一两分钟建议加个进度提示比如每隔若干格disp一下免得干等。4.4 第三步搜索平衡点并识别鞍点平衡点求解我提供两套思路。第一套是快速脚本方案直接用fsolve搭配多组初始猜测。第二套是稳稳妥妥的网格扫描方案先粗定位再精化。下面给出第一套实现init_guesses [0, 0; 0.1, 0.3; -0.1, -0.3; 0.2, 0.6; -0.2, -0.6; ... 0.25, 0.9; -0.25, -0.9; 0.05, 0.15; -0.05, -0.15]; opts optimoptions(fsolve, Display, off, ... FunctionTolerance, 1e-12, ... OptimalityTolerance, 1e-12); equilibria []; for k 1:size(init_guesses, 1) x0 init_guesses(k, :); [xeq, ~, exitflag] fsolve((x) vehicle_state(0, x, p), x0, opts); if exitflag 0 xeq xeq; if isempty(equilibria) equilibria xeq; elseif min(vecnorm(equilibria - xeq, 2, 2)) 1e-4 equilibria [equilibria; xeq]; end end end得到平衡点集合后对每个平衡点计算数值雅可比矩阵并求特征值判断鞍点saddle_point []; for k 1:size(equilibria, 1) xeq equilibria(k, :); J numerical_jacobian((x) vehicle_state(0, x, p), xeq); lambda eig(J); fprintf(平衡点(%.4f, %.4f) 特征值: %.3f, %.3f\n, ... xeq(1), xeq(2), real(lambda(1)), real(lambda(2))); if real(lambda(1)) * real(lambda(2)) 0 saddle_point xeq; end end数值雅可比用中心差分实现比起直接调用符号工具箱速度快且没有解析求导的烦恼function J numerical_jacobian(fun, x0) h 1e-6; n length(x0); J zeros(n, n); for i 1:n xp x0; xp(i) xp(i) h; xm x0; xm(i) xm(i) - h; J(:, i) (fun(xp) - fun(xm)) / (2 * h); end end中心差分步长h选1e-6比较合适。太大截断误差明显太小浮点舍入误差主导。1e-6在大多数车辆动力学方程尺度下都很安全。如果你用网格扫描方式粗定位大致是这样把β-γ平面分成50×50网格计算每个网格点的状态变化率dx找dx(1)0与dx(2)0的等值线交点然后用这些交点作fsolve的初始猜测。这个方法能抓到的平衡点数量通常比固定猜测多尤其在鞍点位置偏离常规认知时。4.5 第四步从鞍点出发绘制临界轨迹识别出鞍点后计算该点的雅可比矩阵找到正实部特征值对应的特征向量作为不稳定流形方向。从鞍点偏移±ε·v后正向积分得到两条临界轨迹分支。如果是直线行驶对称工况你会发现两条分支几乎对称如果加了前轮转角两条分支就是非对称的很有意思。if ~isempty(saddle_point) J_s numerical_jacobian((x) vehicle_state(0, x, p), saddle_point); [V, D] eig(J_s); lambda diag(D); [~, idx_plus] max(real(lambda)); v_unstable V(:, idx_plus); v_unstable v_unstable / norm(v_unstable); eps_offset 1e-3; t_span_fwd [0 5]; options_strict odeset(RelTol, 1e-8, AbsTol, 1e-10); x_start1 saddle_point eps_offset * v_unstable; x_start2 saddle_point - eps_offset * v_unstable; [~, x_branch1] ode45((t, x) vehicle_state(t, x, p), t_span_fwd, x_start1, options_strict); [~, x_branch2] ode45((t, x) vehicle_state(t, x, p), t_span_fwd, x_start2, options_strict); plot(x_branch1(:, 1), x_branch1(:, 2), r-, LineWidth, 2.2); plot(x_branch2(:, 1), x_branch2(:, 2), r-, LineWidth, 2.2); plot(saddle_point(1), saddle_point(2), rp, MarkerSize, 15, ... MarkerFaceColor, r, MarkerEdgeColor, k); end我特别强调一下积分容差。临界轨迹的绘制对数值误差极敏感尤其是从鞍点出发这种极度不稳定的流形。RelTol我直接用1e-8而不是默认的1e-3误差容差设到1e-10。别嫌慢这条轨迹的稳定性直接决定了你画出来的边界是光滑曲线还是锯齿线。如果你还想绘制稳定流形即从鞍点反向积分得到的另一部分边界方法类似但tspan用[0 -5]初始偏移沿稳定特征向量方向。前面说过反向积分数值上更容易发散建议每隔一步就检查一下轨迹是否已经冲到离谱范围必要时用分段积分把上一段的终点作为下一段的起点摸着石头过河。4.6 第五步结果可视化与解读最终出图时我习惯把相平面流线、平衡点、鞍点、临界轨迹四类信息叠在一张图上。为了好读我还会加图例标注清楚哪些是流线、哪些是平衡点、哪些是临界轨迹。如果同一张图里有多组速度或转角工况用subplot排开逐张对比。一个等效可视化的技巧除了画轨迹线还可以对网格上的每个初始状态点做“最终收敛分类”——看它落在哪个平衡点的吸引域内然后用散点颜色填充稳定域和不稳定域。这种“相平面热图”比线条更直观尤其做汇报展示的时候特别好用。实现起来也不复杂就是在原有的16×16网格循环里判断最后时刻状态离哪个平衡点最近再用scatter打点。5. 实操中的常见问题与避坑经验5.1 相平面杂乱难辨问题出在网格和线宽很多新手第一次画出相平面惊呼“怎么全是一团线”。这通常不是代码错而是网格太密、线宽太宽、颜色太深。16×16网格配0.5磅线宽属于比较能看的组合。如果你想更清爽可以把网格降到12×12或者把每条轨迹的透明度调低Color第四个分量。还有一个我自己常用的技巧不画完整4秒轨迹只画后2秒的轨迹片段。因为相平面图里初始段的状态快速移动线条密集看了容易头晕而末段的轨迹慢慢接近吸引子或发散方向信息量反而最集中。5.2 找不到鞍点线性轮胎模型的局限如果你用线性轮胎模型Fy-C·α那整套分析会非常“干净”——干净到没有鞍点可画。线性系统只可能有一个平衡点稳定就是不稳定的收敛到它不稳定就是发散而没有任何边界。所以相平面里出现不了鞍点是因为你还没把非线性轮胎模型引进来。切换到Fiala模型后鞍点会从“无到有”。但注意不是所有参数下都有鞍点。低速、高附着、小转角时整个相平面像一个巨大的漩涡所有轨迹都向原点收敛鞍点消失。你要做的就是提高车速、降低附着系数或者加大前轮转角把失稳边界逼出来。5.3 临界轨迹“跑偏”积分设置与初始扰动的平衡临界轨迹跑偏我前前后后遇到过三种原因。第一种是偏移量ε太大初始点没有落在真实不稳定流形上。结果画出来的“临界轨迹”是条随机的发散曲线完全不能当边界用。这时把ε从1e-3降到1e-4重新跑。第二种是积分容差太松。默认RelTol是1e-3对普通轨迹没问题对鞍点流形就差太远了。我建议至少设到1e-6画分界线时直接上1e-8。第三种是积分时间过长。临界轨迹从鞍点出发后快速延伸如果T_end太大轨迹早就跑到β超过0.5甚至1.0的远场去了图上看起来就是一条冲出边界的线。把T_end限制在3~5秒看它把稳定域围起来就停。5.4 参数敏感换一组车辆参数后的注意事项车辆参数变化会显著影响相平面形态我吃过亏所以特别提醒三点。第一前后轴刚度比例变化会直接改变鞍点位置和临界轨迹的倾斜方向。Cf和Cr互换一下临界轨迹可能从“左上右下”变成“右上左下”。第二速度Vx的影响最大。Vx从15m/s变到30m/s稳定域可能从几乎覆盖整个图面缩小到只有原点周围一小圈。扫速度时我建议每次增加5m/s连续画3~5张图你就能很清楚地看到稳定域收缩的过程。第三μ附着系数控制着轮胎饱和点。μ越低饱和点越早出现鞍点越靠近原点。雨雪路面μ0.3左右的稳定域比干燥路面缩小一半以上这是ESP算法最关注的一个场景。5.5 平衡点搜索的漏解问题fsolve多初始猜测方案虽然简单但确实容易漏根。如果你发现相平面上某个区域的轨迹明显“绕着某个点转”但平衡点列表里没有那个点那大概率就是漏解了。我的改进方法前面提过的网格扫描法先用50×50网格计算dx(1)、dx(2)的符号场找它们分别从正变负的等值线交点把这些交点作为fsolve的初始猜测。这样基本能保证把所有平衡点一网打尽。条件允许的话还可以用pde或者continuation方法比如对Vx做参数扫描从已知平衡点逐步追踪来跟踪鞍点随车速的变化轨迹但那是另一个层面的工作一般项目用不着。还有个小细节如果系统有多个很接近的平衡点去重阈值不能设得太大。我前面代码里用1e-4够用。如果设成1e-2可能把物理上不同的平衡点错误合并导致信息丢失。一点个人经验整套流程跑通之后我最深的体会是相平面看起来只是一堆轨迹线但配合鞍点和临界轨迹它其实把车辆的稳定裕度直接画在了你面前。每次调ESP控制阈值、评估某个底盘参数改动带来的稳定性变化我都会先拉一组相平面看看特定工况下稳定域的形态心里先有个底再去做控制参数整定。最后分享一个我一直在用的小技巧不要只画单一速度下的相平面用一个subplot把15、20、25、30 m/s四个速度的相平面排在一起对比。你很快就会发现临界轨迹的“收拢”过程特别直观——稳定域从开阔地变成一条窄走廊车辆的安全余量到底有多大一目了然。这个方法我推荐给每一个做车辆稳定性的朋友比纯看时序响应靠谱多了。