做车辆稳定性分析的人早晚都要跟相平面打交道。我第一次用MATLAB画车辆相平面是在做操纵稳定性课题的时候——看横摆角速度的时间响应曲线看得头大转向输入一给曲线一颤一颤的根本说不清这车到底是“稳”还是“不稳”。后来改成在相平面上看轨迹状态到底会不会跑偏、能不能绕回平衡点一眼就清楚。这篇文章就聊聊怎么用MATLAB把车辆稳定性相平面画出来包括状态方程怎么搭、程序怎么组织、稳定域边界怎么叠加以及我在绘图调试中踩过的几个坑。适合正在做车辆动力学仿真、底盘控制策略验证、或者写毕业论文需要出相平面图的朋友照着抄就能跑出第一批结果。1. 为什么车辆稳定性分析总绕不开相平面1.1 只看时间响应曲线很难判断“全局稳不稳”很多刚接触车辆稳定性分析的人第一步习惯是给某个初始扰动然后看β(t)、r(t)曲线是否收敛。这个做法没错但有几个实际局限。第一时间响应只代表你选的那一个初始状态。车辆实际运行中路面激励、驾驶员误操作带来的初始状态千差万别单条曲线只能说明“这个初值下系统是收敛的”不能说明“另一个更大的侧偏角下还收不收敛”。第二β和r是相互耦合的看两条一维曲线很难直观感受耦合关系——β回稳了但r还在振荡这种现象在曲线上容易误判。第三车辆动力学模型本质上是非线性的非线性系统在相平面上的行为特征比如鞍点、分界线、稳定域边界在时间域曲线里几乎无法直接读出。相平面最大的价值是把“状态值”和“状态变化率”同时放在一张二维图上。纵轴是状态导数横轴是状态本身系统从某个状态出发下一时刻怎么走全由轨迹的切线方向决定。你看相平面图就像看一张地形图哪里有谷底稳定平衡点、哪里有山脊分界线、哪里会滑落一目了然。1.2 车辆稳定性相平面到底看哪两个状态车辆侧向动力学里常用的是二自由度模型状态变量通常取质心侧偏角β和横摆角速度r。但画相平面时有两种主流选择β-r平面横轴质心侧偏角纵轴横摆角速度适合分析整车横摆运动与侧滑运动的耦合关系。β-βdot平面横轴质心侧偏角纵轴质心侧偏角变化率这是论文里最常见的组合因为βdot直接反映侧向加速度趋势工程上可以由侧向加速度传感器推算试验标定方便。两种平面结论等价但β-βdot平面在控制层更好理解后面我会专门讲程序的切换方式。相平面图的核心信息可以拆成四个要素平衡点状态导数等于0的点系统可能长期停留的工况。稳定平衡点轨迹围绕它收敛一般是焦点或节点对应车辆稳定行驶状态。鞍点一个方向吸引、另一个方向排斥的平衡点是整个相平面里最“危险”的位置稳定域边界通常锚定在鞍点上。分界线把收敛区域和发散区域隔开的轨迹工程上叫稳定域边界。读数的时候记住一条经验轨迹如果从初始点一路走向稳定平衡点说明车辆在该初始状态下能恢复稳定如果轨迹越过鞍点附近的分界线跑向发散区域说明车辆已经进入失稳工况需要ESC等控制系统介入。相平面方法之所以在底盘稳定性控制里被广泛采用就是因为它能把这个“稳/不稳”的判定变成可计算的区域划分问题。2. 二自由度参考模型推导与MATLAB状态方程构建2.1 模型假设为什么二自由度就够用做相平面分析不需要上复杂的19自由度整车模型二自由度自行车模型是标准起点。假设如下车辆纵向速度u恒定不考虑纵向动力学左右轮合并前后轴各用一个等效轮胎表示侧偏力与侧偏角在工作点附近近似线性轮胎侧偏刚度取常数忽略空气阻力、侧倾、俯仰影响。这套假设在侧向加速度不超过0.4g左右的工况下精度足够。相平面分析关注的主要矛盾是“质心侧偏-横摆运动”的耦合关系二自由度模型正好抓住了这个核心。参数太多反而干扰判断——你先用简单模型把相平面的逻辑跑通再往模型里加非线性轮胎、侧倾自由度都来得及。2.2 从受力平衡到状态方程的完整推导设车辆质量m横摆转动惯量Iz质心到前轴距离a到后轴距离b轴距Lab前轮转角δ质心侧偏角β横摆角速度r。侧向力平衡方程m u (βdot r) Fyf Fyr横摆力矩平衡方程Iz rdot a Fyf - b Fyr这里Fyf和Fyr是前后轴等效侧偏力。前轮侧偏角近似为αf δ - (β a r / u)后轮侧偏角近似为αr -(β - b r / u)在线性轮胎假设下Fyf Cf αfFyr Cr αr。其中Cf、Cr取正数整套符号约定保证方向盘左转输入时侧向力向左物理自洽。代入整理后得到状态方程βdot [Cf (δ - β - a r/u) Cr (-β b r/u)] / (m u) - rrdot [a Cf (δ - β - a r/u) - b Cr (-β b r/u)] / Iz写成矩阵形式取δ0[βdot; rdot] A [β; r]A矩阵四个元素分别为A(1,1) -(Cf Cr)/(m u) A(1,2) -1 - (a Cf - b Cr)/(m u²) A(2,1) -(a Cf - b Cr)/Iz A(2,2) -(a² Cf b² Cr)/(Iz u)这个矩阵形式主要用来算平衡点和特征值判断局部稳定性后面会用到。2.3 车辆参数与状态方程MATLAB函数我在建模型阶段惯用参数表如下是一组典型轿车的参考值你换成自己项目里的参数即可参数符号数值单位整车质量m1470kg横摆转动惯量Iz2350kg·m²质心到前轴距离a1.05m质心到后轴距离b1.41m前轴等效侧偏刚度Cf80000N/rad后轴等效侧偏刚度Cr80000N/rad纵向速度u20m/s前轮转角δ0rad状态方程函数写成MATLAB文件方便ode45直接调用function dx vehicle_single_track(t, x, P) % x [beta; r] beta x(1); r x(2); delta P.delta; alpha_f delta - (beta P.a * r / P.u); alpha_r -(beta - P.b * r / P.u); Fyf P.Cf * alpha_f; Fyr P.Cr * alpha_r; beta_dot (Fyf Fyr) / (P.m * P.u) - r; r_dot (P.a * Fyf - P.b * Fyr) / P.Iz; dx [beta_dot; r_dot]; end参数用结构体P统一管理后面做变参数扫描比如改变车速u、路面附着μ非常方便不用改函数内部代码。注意侧偏刚度的符号定义不同教材有差异关键是保证自己整套公式自洽——你只要发现相轨迹的收敛方向反了把α表达式里的正负号整体反过来就行。3. MATLAB相平面绘制核心程序从单条轨迹到轨迹场3.1 单条相轨迹先跑通一条线万事开头难第一条轨迹跑通后面就是批处理。给定初始状态β0、r0调用ode45P.m 1470; P.Iz 2350; P.a 1.05; P.b 1.41; P.Cf 80000; P.Cr 80000; P.u 20; P.delta 0; x0 [8*pi/180; 0.3]; % 初始质心侧偏角 8°初始横摆角速度 0.3 rad/s options odeset(RelTol, 1e-6, AbsTol, 1e-8); [t, x] ode45((t, x) vehicle_single_track(t, x, P), [0 10], x0, options); figure; plot(x(:,1)*180/pi, x(:,2), b-, LineWidth, 1.5); xlabel(质心侧偏角 \beta (deg)); ylabel(横摆角速度 r (rad/s)); grid on;注意单位问题状态方程里β统一用弧度计算画图时乘以180/pi换成度。横摆角速度单位rad/s保留。这一步跑完你应该看到一条从初值点出发、螺旋收敛到原点的轨迹。3.2 网格初值批量生成相轨迹场单条轨迹只能说明一个初始状态要画出论文里那种漂亮的相轨迹场就得在β-r平面上铺一片初值网格对每个初值积分一条轨迹beta0_list -12*pi/180:1.5*pi/180:12*pi/180; r0_list -0.6:0.05:0.6; figure; hold on; for i 1:length(beta0_list) for j 1:length(r0_list) x0 [beta0_list(i); r0_list(j)]; [~, x] ode45((t, x) vehicle_single_track(t, x, P), [0 10], x0, options); plot(x(:,1)*180/pi, x(:,2), k-, LineWidth, 0.3); end end hold off;这段代码跑出来就是一张密密的轨迹场图。网格间距怎么选我实测的经验是β方向步长1°到2°r方向步长0.04到0.06 rad/s比较合适。太疏看不出分界线的走向太密整张图糊成一团黑色而且计算时间成倍增加。如果电脑性能一般先用粗网格验证逻辑再加密出图。3.3 平衡点标注与局部稳定性验证轨迹图上只有曲线还不够要把平衡点标出来才有参考价值。δ0时线性模型的平衡点可以通过fsolve求eq_fun (x) vehicle_single_track(0, x, P); x_eq fsolve(eq_fun, [0; 0]); plot(x_eq(1)*180/pi, x_eq(2), ro, MarkerSize, 8, MarkerFaceColor, r);对线性模型原点就是唯一平衡点。判断局部稳定性就求A矩阵特征值A [-(P.CfP.Cr)/(P.m*P.u), -1-(P.a*P.Cf-P.b*P.Cr)/(P.m*P.u^2); -(P.a*P.Cf-P.b*P.Cr)/P.Iz, -(P.a^2*P.CfP.b^2*P.Cr)/(P.Iz*P.u)]; eig(A)如果特征值实部都在左半平面原点局部稳定相轨迹会螺旋收进原点。这一步不要跳过先确认模型整车参数在这个车速下的稳定性趋势再去看轨迹场心里才有底。4. 稳定域边界叠加从线性模型到非线性轮胎模型4.1 为什么线性模型画不出分界线很多人在论文里看到的相平面图都有一个像“鱼嘴”一样的稳定域边界边界内轨迹收敛到平衡点边界外轨迹发散。但你可能已经发现上面用线性轮胎模型跑出来的轨迹场所有轨迹最后都收敛到原点根本没有什么“边界”。原因在于线性模型只有一个全局平衡点系统要么全局稳定要么全局发散不存在“区域性的稳定”。真实车辆之所以有大侧偏角下失稳的现象是因为轮胎侧偏力在大侧偏角时饱和了——侧偏角再大力也上不去了。这个非线性是产生鞍点和稳定域边界的物理根源。4.2 轮胎饱和模型与状态方程改造在不引入魔术公式之前可以先用一个“线性限幅”的简化饱和模型把现象复现出来。侧偏力线性增长但达到附着极限μFz之后保持恒定function Fy tire_limited(C, alpha, mu, Fz) % 线性限幅饱和轮胎力模型 F_lin C * alpha; F_max mu * Fz; Fy max(min(F_lin, F_max), -F_max); end前后轴垂直载荷按下式估算Fzf m g b / (a b) Fzr m g a / (a b)改造后的状态方程函数function dx vehicle_single_track_sat(t, x, P) beta x(1); r x(2); delta P.delta; alpha_f delta - (beta P.a * r / P.u); alpha_r -(beta - P.b * r / P.u); Fyf tire_limited(P.Cf, alpha_f, P.mu, P.Fzf); Fyr tire_limited(P.Cr, alpha_r, P.mu, P.Fzr); beta_dot (Fyf Fyr) / (P.m * P.u) - r; r_dot (P.a * Fyf - P.b * Fyr) / P.Iz; dx [beta_dot; r_dot]; endP里增加P.mu 0.85P.Fzf和P.Fzr按上面公式算好。用这套模型重新跑网格初值你会发现三个平衡点原点附近区域轨迹收敛两侧出现两个鞍点鞍点向外方向轨迹发散。鞍点连接起来的那条线就是稳定域分界线。4.3 多初值法绘制稳定域包络鞍点位置可以通过fsolve分别搜索左右两侧的解x_eq_left fsolve((x) vehicle_single_track_sat(0, x, P), [-0.25; 0]); x_eq_right fsolve((x) vehicle_single_track_sat(0, x, P), [0.25; 0]);精确绘制分界线需要求鞍点的稳定流形在程序上要把状态方程反向积分操作有一定门槛。工程上更常用的做法是“多初值法”把相平面区域划分成网格对每个网格点积分判断终点是否落在稳定平衡点附近然后用等值线画出发散/收敛的分界。beta0_grid linspace(-0.4, 0.4, 41); r0_grid linspace(-0.8, 0.8, 41); stable zeros(length(beta0_grid), length(r0_grid)); x_eq_stable [0; 0]; % 原点稳定平衡点 tol 0.005; % 判定收敛半径 for i 1:length(beta0_grid) for j 1:length(r0_grid) x0 [beta0_grid(i); r0_grid(j)]; [~, x] ode45((t, x) vehicle_single_track_sat(t, x, P), [0 10], x0, options); xf x(end, 1:2); stable(i, j) norm(xf - x_eq_stable) tol; end end figure; contourf(beta0_grid*180/pi, r0_grid, stable, [0.5 0.5], k-, LineWidth, 2);这段代码跑完后你能清晰地看到稳定域在β维度的边界大约在正负十几度附近超过这个范围状态就收不回来了。这个稳定域包络线就是ESC、ABS等底盘控制策略设计时最关心的“安全边界”。注意contourf的坐标方向stable矩阵第一维对应beta0_grid第二维对应r0_grid转置后传入才能保证轴方向正确。别问我怎么知道的我第一次画出来稳定域横竖颠倒调了半小时才发现是维度顺序问题。5. 相平面程序实战几个“看起来对但结果错”的坑5.1 仿真时间太短轨迹中途截断像发散用ode45跑相轨迹时最隐蔽的坑是仿真时间设太短。比如[t, x] ode45(..., [0 3], ...)3秒可能只够轨迹绕半圈还没收敛到平衡点。图上的轨迹末端显得很“开放”如果不注意会误判为发散。我当时调试一个高速工况明明特征值显示稳定相轨迹图却密密麻麻像个爆炸图。后来把一条单条轨迹的时间拉长到20秒才看清它其实在慢慢螺旋收敛只是低速衰减很慢而已。解决思路有两个。一是合理设置积分时长一般取10秒起步画完检查轨迹末端是否集中在平衡点附近二是用事件函数提前终止积分既省时间又不影响结果function [value, isterminal, direction] stopEvent(t, x, eq, tol) value norm(x - eq) - tol; isterminal 1; direction 0; end调用时把事件函数加进odesetoptions odeset(RelTol, 1e-6, AbsTol, 1e-8, ... Events, (t, x) stopEvent(t, x, [0;0], 0.001));这样轨迹一旦进入平衡点邻域就停止积分批量跑网格的效率能提升好几倍。5.2 网格密度与绘图次序的博弈批量画轨迹场时网格步长和线宽非常影响成图质量。β方向间隔太小、r方向间隔太小400个初值以下还行上千条轨迹挤在一起Matlab的绘图开销和内存占用会明显上升而且图面全是黑线绞成一团。我现在的经验是分两轮出图第一轮用粗网格β步长2°r步长0.06快速验证趋势第二轮用中等网格β步长1°r步长0.04出最终图。线宽用0.2到0.3颜色用黑色或深灰色网格太细密时可以用LineColorAlpha之类的透明度属性或者隔一个点连绘一条保证图面清爽。还有一个绘图次序的细节先画所有离散轨迹最后叠加平衡点、分界线等元素。否则平衡点被大量轨迹线盖住标注符号根本看不见。5.3 单位混乱弧度、角度、π/180相平面程序里单位问题几乎是每个人都会踩的坑。状态方程内部必须统一用弧度制但画图时β要用度r用rad/s。如果你在ode45的x0里直接写了beta0 8表示8弧度那画出来的横轴范围是正负几百度的奇怪图形。我个人的习惯是所有初值在代码里显式做转换比如beta0 8*pi/180绝不在脑子里做心算。侧偏刚度单位是N/rad垂直载荷单位是N轮胎力单位是N代入状态方程前先检查一遍省得后面出图数据整体偏移查半天。5.4 只画线不画方向图再好看也无效相轨迹是一条静态曲线如果不标注时间方向读者很难判断轨迹是“从外往里走”还是“从里往外走”。这在稳定/发散的判定上是致命问题。补充方向信息有三种常用方法等时间间隔取点用plot画点标记点与点之间间隔越大说明运动越慢沿轨迹每隔固定弧长用quiver画一个小箭头箭头指向局部速度方向用quiver直接画网格点上的向量场[beta_mesh, r_mesh] meshgrid(-0.3:0.03:0.3, -0.5:0.05:0.5); bdot zeros(size(beta_mesh)); rdot zeros(size(r_mesh)); for k 1:numel(beta_mesh) dxk vehicle_single_track_sat(0, [beta_mesh(k); r_mesh(k)], P); bdot(k) dxk(1); rdot(k) dxk(2); end quiver(beta_mesh*180/pi, r_mesh, bdot, rdot, 0.5, r);向量场加轨迹场双重显示既能看到全局面貌又能单条追踪具体路径是我最常用的组合。5.5 忘检查特征值高阶不稳定工况被忽略批量扫描不同车速或不同前轮转角时一定要每一步都检查平衡点稳定性。车辆在低附着路面或高速工况下原点可能已经不是稳定焦点而是变成了不稳定平衡点。这时候整个相平面图的结构会发生质变原来那个“鱼嘴”稳定域会消失所有轨迹都发散。我第一次做μ0.3低附着工况扫描时就遇到这个问题相平面图看起来乱七八糟一度以为代码写错了。后来把特征值打印出来才发现原点已经不稳定了气贯长虹的结果才是对的。所以批量扫描时先算特征值、再画相图顺序不要颠倒。6. 相平面程序的进阶玩法从“画图”到“评价指标”6.1 β-βdot平面的程序切换论文里更常见的β-βdot相平面怎么出其实不用改状态方程只需要把βdot作为纵轴数据存下来。在状态方程函数里beta_dot本来就计算了直接输出即可。批量轨迹不需要额外求导[t, x] ode45((t, x) vehicle_single_track_sat(t, x, P), [0 10], x0, options); % 计算每条轨迹对应的 beta_dot beta_dot_traj (P.Cf*(P.delta - x(:,1) - P.a*x(:,2)/P.u) ... P.Cr*(-x(:,1) P.b*x(:,2)/P.u)) / (P.m*P.u) - x(:,2); plot(x(:,1)*180/pi, beta_dot_traj, k-, LineWidth, 0.3);注意βdot的物理单位是rad/s和横摆角速度r一致画图时不用乘π/180。工程上βdot可以由侧向加速度ay按下式估计βdot ≈ ay / u - r这意味着你不需要高精度侧偏角观测器用IMU的侧向加速度和横摆角速度信号就能在实车上画出相平面。这也是β-βdot平面在工程标定中比β-r平面更受欢迎的原因。6.2 变车速扫描看稳定域怎么收缩固定δ0把车速从10 m/s扫到35 m/s你会看到稳定域的宽度随车速提升明显变窄。程序上做一个循环即可u_list [10 15 20 25 30 35]; for k 1:length(u_list) P.u u_list(k); % 重新计算平衡点、画相图、统计稳定域面积 end从控制角度看这个稳定域面积随车速的变化曲线比单一工况下的相平面图更有信息量。ESC介入时机、控制增益的调度很多底盘团队就是参考这条曲线定的。6.3 把稳定域面积程序化作为底盘调参指标多初值法得到stable矩阵后统计稳定网格点数占全部网格点数的比例或者换算成物理面积就得到了一个量化指标。这个指标可以作为底盘参数的优化目标比如调整前后轴侧偏刚度分配、质心位置看哪个参数组合让稳定域面积最大。我在实际项目中是把这套逻辑封装成了一个函数function area_ratio calc_stable_area(P, beta_lim, r_lim, grid_res) % 输入车辆参数、相平面范围、网格分辨率 % 输出稳定域面积占比 % 内部完成网格积分、收敛判定、面积统计 end这样每次调参数跑一遍就能看到稳定域是变大还是变小比肉眼看图判断客观得多。这一步做完相平面就从“一张图”变成了“一个可优化的评价函数”实用性立刻不同。最后再分享一个小技巧批量跑相图时记得用print导出矢量图格式选svg或eps线宽和标注在论文里清晰很多。print(gcf, phase_plane, -dsvg, -r600);我试过用默认png导出图上细线直接糊成一团改成矢量图之后导师一眼就看清了稳定域边界沟通效率高了不少。