做车辆稳定性仿真的人迟早会碰上相平面这三个字。我第一次用 Matlab 画二自由度车辆的质心侧偏角—横摆角速度相平面是在研究 ESP 触发策略的时候。相平面能把系统的全局运动摊开在一张图上鞍点和临界轨迹直接告诉你车辆什么时候会突然从“能绕回来”变成“完全失控”。这套方法放在真实开发中非常实用尤其适合搞底盘控制、车辆动力学方向的研究生和工程师参考。这篇就把整套思路拆开讲清楚从二自由度模型建立、相平面绘制、鞍点定位到临界轨迹求解全部是能在 Matlab 里直接跑通的仿真流程。1. 相平面分析解决什么问题稳定边界的可视化1.1 为什么二自由度模型够用车辆工程里最常用的动力学模型是多体模型但做稳定性分析时二自由度模型反而是首选。它把整车抽象成只有质心侧偏角β和横摆角速度r两个自由度的系统忽略垂向载荷转移、悬架变形和空气动力学影响只保留侧向动力学中最核心的运动关系。这个简化听起来粗糙实际上却抓住了失稳的关键。车辆出现侧滑或者甩尾时最先变化的往往是质心侧偏角和横摆角速度这两个量直接反映车辆绕垂直轴转动和横向滑动的基本状态。用它们做坐标轴画相平面能把复杂整车状态压缩到二维平面上观察直观程度远超单纯看一条时域曲线。二自由度模型还有一个好处计算量小方便快速迭代。在做控制算法初步验证时不必每次都用 Carsim 或 Simulink 搭建十几自由度的整车模型直接用二自由度模型在 Matlab 里做上万次仿真也只要几秒。这也是为什么相平面方法被大量用于底盘控制策略开发的前期阶段包括稳定边界计算、控制阈值设定、路面附着估计等。1.2 相平面上的三类关键对象相平面不是随便画几条曲线就完了核心要看三个要素平衡点、鞍点、临界轨迹。平衡点是系统状态方程的零解即令质心侧偏角微分为零、横摆角速度微分为零时得到的点。线性二自由度模型通常只有一个稳定平衡点代表车辆最终进入稳恒转向状态。鞍点则是不稳定平衡点它像一座山峰和山口之间的连接处一部分轨迹被吸引过来另一部分轨迹又会被推开。在车辆稳定性问题里鞍点的存在意味着存在一个“临界状态”一旦车辆状态越过这个边界就会不可逆转地走向失稳。临界轨迹就是从鞍点附近延伸出的特殊轨线它把相平面划分为稳定域和不稳定域。只要车辆初始状态落在稳定域内即使有扰动也会收敛回稳定平衡点一旦落在稳定域外即使扰动消失车辆也会越走越偏直到完全失稳。三者结合就构成了一张车辆稳定性地图。这也是整篇文章要复现的核心结果。2. 模型搭建与状态方程推导先把数学摆清楚2.1 车辆参数和符号约定画相平面之前必须先把二自由度车辆模型的微分方程建立起来。这是后面的所有仿真基础参数定义错一个符号相平面就会完全变样。我常采用以下符号符号含义单位m整车质量kgVx纵向速度恒定m/slf质心到前轴距离mlr质心到后轴距离mCf前轮侧偏刚度取正值N/radCr后轮侧偏刚度取正值N/radIz横摆转动惯量kg·m²δ前轮转角radβ质心侧偏角radr横摆角速度rad/s侧偏刚度的符号坑过一次很多人。很多教科书把侧向力写成 Fy-Cα这时的 C 是正数也有教材直接写成 FyCα此时 C 为负数。统一采用正值的 Cf、Cr后续力表达式中要显式写出负号。下面的推导全部基于 Fy-Cα 的约定。2.2 完整状态方程推导根据牛顿第二定律和力矩平衡纵向速度 Vx 恒定情况下二自由度模型的动力学方程为m · Vx · (dβ/dt r) Fyf FyrIz · dr/dt lf · Fyf - lr · Fyr其中前后轴侧向力分别为Fyf -Cf · (β lf · r / Vx - δ)Fyr -Cr · (β - lr · r / Vx)把这两个式子代入整理成状态空间形式dβ/dt [ -Cf·(β lf·r/Vx - δ) - Cr·(β - lr·r/Vx) ] / (m·Vx) - rdr/dt [ -lf·Cf·(β lf·r/Vx - δ) lr·Cr·(β - lr·r/Vx) ] / Iz这一步用的是线性轮胎模型适合小侧偏角工况。实际车辆在极限工况下轮胎会饱和所以工程上画相平面时经常把线性侧偏刚度替换成双线性或魔术公式轮胎模型。替换的思路是先根据当前轮胎法向载荷和侧偏角计算饱和后的侧向力再代入同一套状态方程。下面仿真中如果希望结果更贴近实际可以把 Cf、Cr 替换为 Fy_f(αf) 和 Fy_r(αr)其中 αfβlf·r/Vx-δαrβ-lr·r/Vx。2.3 为什么用恒定纵向速度的假设一定要注意这个二自由度模型假设 Vx 恒定。实际车辆制动或加速时纵向速度变化会引入额外的状态变量使得相平面变成三维甚至更高维。但大多数稳定性分析中重点关心的是稳态转向或紧急工况下车辆是否失稳可以近似认为在很短时间内 Vx 变化不大所以保持恒定 Vx 是可接受的工程近似。如果做更精确的联合仿真可以把这套相平面分析输出给 Carsim 或者 dSPACE 模型让纵向速度动态变化再对比相平面结果。这样既保持了二自由度模型的快速计算能力又具备了高保真模型的准确度是工程上常见的两段式做法。3. Matlab 相平面绘制全流程从微分方程到矢量场3.1 建立状态方程函数先在 Matlab 里写一个状态方程函数输入为 t、状态量 x[β; r]、以及车辆参数输出为 [dβ; dr]。这里给出一段可直接用的参考代码function dx vehicle_states(t, x, param) beta x(1); r x(2); Vx param.Vx; delta param.delta; % 前、后轴侧偏角 alpha_f beta param.lf * r / Vx - delta; alpha_r beta - param.lr * r / Vx; % 线性侧向力模型高附着近似 Fyf -param.Cf * alpha_f; Fyr -param.Cr * alpha_r; dx zeros(2,1); dx(1) (Fyf Fyr) / (param.m * Vx) - r; dx(2) (param.lf * Fyf - param.lr * Fyr) / param.Iz; end如果使用非线性轮胎可以在vehicle_states中把 Fyf、Fyr 改成从查表函数中获得后面鞍点搜索和轨迹绘制不受影响。3.2 画相速度矢量场相平面里的每一个点 (β, r)代入状态方程后得到的速度向量 (dβ, dr)就是该点上的“风”。用 quiver 函数可以把整个速度场画出来这是相平面图最基础的骨架。beta_range -0.3:0.02:0.3; r_range -0.8:0.05:0.8; [Beta, R] meshgrid(beta_range, r_range); dBeta zeros(size(Beta)); dR zeros(size(R)); for i 1:numel(Beta) x [Beta(i); R(i)]; dx vehicle_states(0, x, param); dBeta(i) dx(1); dR(i) dx(2); end figure; quiver(Beta, R, dBeta, dR, b); xlabel(质心侧偏角 \beta (rad)); ylabel(横摆角速度 r (rad/s));网格步长要选合适。步长太大速度场看起来断断续续步长太小计算量成倍增加图片反而密集得看不清。一般 β 方向取 0.02 radr 方向取 0.05 rad/s已经能得到不错的视觉效果。如果关注某个局部区域再在该区间加密网格。3.3 从初值积分绘制相轨线矢量场只能给出每一个点的瞬时运动方向要看清具体演化路径还得从不同初值出发做轨迹积分。用 ode45 在给定初值上做正向积分即可t_span 0:0.01:5; x0_list [-0.1, -0.4; -0.1, 0.4; 0.05, 0.6; -0.2, 0.1]; hold on; for i 1:size(x0_list,1) x0 x0_list(i,:); [~, x_traj] ode45((t,x) vehicle_states(t,x,param), t_span, x0); plot(x_traj(:,1), x_traj(:,2), k-, LineWidth, 1.2); end这里有一个很容易踩的坑初始点如果正好落在临界轨迹上积分时间太长时轨迹会趋于鞍点而不是发散看起来像停在原地如果初始点离临界轨迹稍远轨迹才会明显朝稳定点或边界外方向走。所以画图时多选几组初值每组之间留一点间距才能真正看到稳定域和不稳定域的划分效果。4. 鞍点定位与临界轨迹绘制核心环节4.1 平衡点搜索与鞍点判别平衡点是状态方程中 dBeta0 且 dR0 的位置。线性模型可以通过解方程组直接求出但非线性模型只能数值搜索。Matlab 中可以用 fsolve 来寻找平衡点给定多个初值以避免漏解fun (x) vehicle_states(0, x, param); x_eq_candidates [-0.1 -0.1; 0.1 0.1; 0.2 0.5; -0.2 -0.5]; for i 1:size(x_eq_candidates,1) [x_eq, fval, exitflag] fsolve(fun, x_eq_candidates(i,:), optimoptions(fsolve,Display,off)); if exitflag 0 disp([平衡点: beta, num2str(x_eq(1)), , r, num2str(x_eq(2))]); end end找到候选平衡点后需要判断它是稳定点还是鞍点。方法是计算该点处系统的雅可比矩阵然后看特征值。若特征值实部均为负则是稳定平衡点若一个特征值实部为正、一个为负则是鞍点。雅可比矩阵可以用符号微分推导也可以用数值差分逼近。工程上我常用jacobianest这类数值工具箱或者直接在fsolve得到平衡点后用差分法求A zeros(2,2); delta_perturb 1e-6; for j 1:2 xp x_eq; xp(j) xp(j) delta_perturb; xm x_eq; xm(j) xm(j) - delta_perturb; A(:,j) (vehicle_states(0, xp, param) - vehicle_states(0, xm, param)) / (2*delta_perturb); end eig(A)注意差分步长别取得太大否则雅可比矩阵精度受损鞍点误判率会明显上升。1e-6 左右通常足够。4.2 沿特征向量反推临界轨迹临界轨迹其实就是鞍点的稳定流形和不稳定流形。对二维系统来说鞍点附近存在两条特征方向其中一条对应正实部特征值另一条对应负实部特征值。从鞍点出发沿着正特征方向向前积分得到不稳定流形沿着负特征方向反向积分得到稳定流形。把这两条流形画出来就是临界轨迹在相平面上的完整形态。实际操作如下先求出鞍点处的特征向量矩阵 V取对应于正实部特征值的特征向量 v_pos以及对应于负实部特征值的特征向量 v_neg。然后在鞍点位置基础上沿特征向量正负方向偏移一个很小的量比如 0.001 倍的特征向量长度作为积分初值。% 假设 x_saddle 是鞍点V(:,2) 对应正实部特征向量 disturb 0.001; t_plus 0:0.01:8; t_minus 0:-0.01:-8; x0_pos x_saddle disturb * V(:,2); traj_pos ode45_traj(x0_pos, t_plus); x0_neg x_saddle - disturb * V(:,2); traj_neg ode45_traj(x0_neg, t_plus);为什么偏移量要取这么小因为鞍点本身是不稳定平衡点初值如果偏离鞍点太远积分出来的轨迹就不是流形本身而是流形附近的普通轨迹。取一个很小扰动让初始点足够靠近特征向量方向这样积分出来的曲线会无限逼近临界轨迹。反向积分时要特别注意时间方向取负值。稳定流形上的点在正向时间下会收敛到鞍点但如果直接正向积分数值上几乎看不见这段轨迹因为误差会把轨线推开。反向积分反而能让轨迹从鞍点附近“回溯”出来得到完整的稳定流形。4.3 稳定边界与车辆稳定性区域解读当临界轨迹画出来后相平面上的稳定区域边界就清楚了。一般情况下稳定区域是围绕稳定平衡点的封闭或半封闭区域临界轨迹构成了区域的边界。车辆状态在稳定区域内受到扰动后会自动回到稳定点状态一旦跨出临界轨迹即使转向角不再变化也会进入失稳区。这个边界对应到实车上就是 ESP 或其他稳定性控制系统希望车辆不要越过的状态阈值。有了相平面图控制策略可以提前设定一条安全距离比如当质心侧偏角和横摆角速度的组合接近临界轨迹时系统提前限制发动机扭矩输出或施加单轮制动。还有一个重要现象是前轮转角 δ 增大时鞍点会逐渐靠近稳定平衡点稳定区域缩小车辆更容易失稳。通过绘制不同 δ 下的相平面图可以直观看到转向角对稳定裕度的影响。这也是相平面法在底盘开发中经常被用来做参数影响分析的原因。5. 仿真实战紧急变道工况下的稳定性分析5.1 工况参数设置这里给出一个真实的仿真实例。假设一辆乘用车参数如下参数取值单位m1500kgIz2500kg·m²lf1.2mlr1.5mCf60000N/radCr80000N/radVx20m/sδ0.1rad计算线性模型的平衡点时令 dBeta0、dR0有dβ/dt [ -60000·(β 1.2·r/20 - 0.1) - 80000·(β - 1.5·r/20) ] / (1500·20) - r 0dr/dt [ -1.2·60000·(β 1.2·r/20 - 0.1) 1.5·80000·(β - 1.5·r/20) ] / 2500 0解这个二元方程得到的平衡点大致在 β≈0.071 radr≈0.20 rad/s 附近。这个点是一个稳定焦点说明在这样的小转角和小侧偏条件下车辆最终能进入稳定的圆周运动状态。同时求解可以发现该模型中只存在一个平衡点没有鞍点。这并不意味着相平面没用了而是线性模型太简单只有在非线性轮胎模型下才会出现远处的鞍点。实际工程画相平面时几乎都使用非线性轮胎模型。5.2 非线性轮胎模型下的相平面结果把前后轮侧向力换成魔术公式的简化形式即饱和型 Fy -D·sin(C·atan(B·alpha))调节参数让侧偏刚度初始斜率与线性模型一致。再同样用 fsolve 搜索平衡点会发现除了原来的稳定焦点外左右两侧还各出现一个鞍点。这样画出来的相平面稳定区域不再是全平面而是被两条临界轨迹围成的一个有限区域。车辆如果初始状态处于区域内最终收敛到中心稳定焦点如果初始状态在区域外轨迹会向边界发散直观表现为质心侧偏角不断增大或横摆角速度剧烈震荡对应实车就是不可控的甩尾或掉头。我把不同初值积分结果画在同一张图上能够清楚看到三条特征曲线稳定点为吸引中心鞍点为分界点临界轨迹像两道屏障把图面分成几个不同的“流域”。这个画面就是车辆稳定性分析的“地图”。5.3 转向输入和车速对鞍点位置的影响用这个非线性模型分别设置 δ0.05、0.1、0.15 rad观察鞍点位置。δ 增大会使稳定焦点向更大侧偏角方向移动同时两侧鞍点向中心靠近。当 δ0.15 rad 时稳定焦点已经非常接近右侧鞍点稳定区域被压缩得很小这意味着车辆在较大转向输入下几乎没有安全裕度稍微一点扰动就会越过临界轨迹进入失稳区。纵向车速 Vx 的影响也很大同样 δ0.1 radVx 从 20 m/s 提升到 40 m/s鞍点与稳定点的距离变化明显稳定域快速收缩。这解释了为什么高速紧急变道比低速更容易失控——不是驾驶员反应慢了而是车辆稳定边界本身就在缩小。这种参数化分析在实际开发中有直接价值。比如标定 ESP 阈值时需要知道不同车速、不同路面附着下车辆还能承受多大的侧偏角组合。在 Matlab 里循环计算一批参数组合的相平面图生成一张“稳定裕度地图”比做几百次场地实验高效得多。6. 常见问题与踩坑记录相平面仿真的避坑指南6.1 轨迹发散与 ode45 求解器设置很多人在积分相轨线时发现轨迹一冲就飞出图边界数值直接变得极大。这是正常现象吗不完全是。如果系统真的不稳定轨迹发散是对的但有时是因为 ode45 默认误差容限下轨迹在接近鞍点时产生了不真实的跳跃或者因为状态方程在极端 β 值下存在数值溢出。我的建议是给状态方程加上角度保护。比如将输出限制在 β∈[-π/2,π/2] 范围内或者使用odeset设置更严格的相对误差options odeset(RelTol, 1e-8, AbsTol, [1e-8 1e-8]); [t, x_traj] ode45((t,x) vehicle_states(t,x,param), t_span, x0, options);RelTol 从默认的 1e-3 改成 1e-8 后靠近鞍点的轨迹精度会明显提升临界轨迹不会再动不动闪断。另外不要把仿真时间设得太长。在稳定域内轨迹会长期绕着稳定点转动画出来的曲线会堆积成一团看不清本来面目。先画 3~5 秒再根据收敛速度调整时间范围比一开始就跑到 20 秒靠谱得多。6.2 鞍点搜索失败的常见原因fsolve搜不到鞍点通常有三种原因第一搜索初值离真实鞍点太远。线性模型下鞍点可能根本不存在但非线性模型下鞍点往往出现在较大侧偏角位置比如 β±0.2 rad 甚至更远。初值范围给得太窄fsolve 容易停在局部平坦区。第二状态方程里用了分段线性或查表侧向力函数不平滑fsolve 默认算法容易失败。这种情况下先改用lsqnonlin做最小二乘搜索或者用MultiStart配合多个随机初值成功率会高很多。第三参数单位混了。比如 Cf 用 N/rad但 β 和 α 用的是 deg换算错误会导致搜索失败。统一用 rad检查所有角度单位。鞍点判定时还有一个细节用数值差分计算雅可比矩阵时扰动步长太大会导致特征值实部符号判断错误。可以先用线性模型解析解验证一下整个流程在线性模型里只存在稳定点没有鞍点如果程序报出鞍点说明数值方法有误换成非线性模型再出现鞍点才是正常现象。6.3 相平面图像的后处理与展示技巧相平面图画出来容易画得清晰有价值却需要经验。第一矢量场和轨迹不要抢视觉。矢量场建议用浅色、细箭头轨迹用深色实线平衡点用大圆点标出鞍点用方形或星形标出图例明确。第二临界轨迹一定要单独强调。推荐用红色粗线把两条流形画出来并在旁边标注“stable boundary”或“临界轨迹”。读者一眼就能看出稳定区域的范围。第三坐标范围要合理。β 范围通常取 -0.3~0.3 radr 范围取 -0.8~0.8 rad/s。如果只关心某一工况可以缩放局部区域但别把图缩得太窄否则看不到全局流形走向。第四把多个工况的相平面做成子图横向对比效果远好于把所有曲线叠加在一张图上。比如同一行放 δ0.05、0.1、0.15 三张图列放 Vx20、40 m/s 两组形成网格状对比图。这个布局在论文和汇报中特别能说明问题。最后别忘了把车辆参数和工况条件标注在图上或标题中。相平面结果对参数极其敏感不写清楚参数过两周自己再看都会怀疑图是怎么画出来的。我在实际项目中还有一个习惯先画线性模型的相平面做一次全流程验证再切换到非线性轮胎模型观察平衡点数量如何变化、稳定域如何收缩。这样既能发现代码里的低级错误也能对非线性带来的影响有更直观的感受。相平面这东西第一次看会觉得很玄一旦把鞍点、临界轨迹和稳定域三者对应到实车失控场景里你会发现它其实是一个特别接地气的分析工具。后续做 ESP 阈值标定、底盘失稳预警算法时这套图可以直接当输入条件用。