Simulink车辆二自由度模型相平面分析：从建模到稳定性判断

做车辆动力学仿真的人一定绕不开一个问题怎么快速判断一辆车在某个车速、某个前轮转角、某条路面附着下到底是稳还是不稳我平时做底盘控制算法最常用的工具之一就是基于Simulink的车辆二自由度模型把质心侧偏角和横摆角速度放到相平面里看。这个标题里的组合就是为此准备的程序核心不复杂胜在可以任意改参数一套模型从线性研究一路用到非线性扩展看完这篇你也能照着搭一个出来。这篇文章适合三类人正在写毕业论文的车辆工程学生、刚接手底盘稳定性控制任务的工程师、以及任何想搞懂二自由度模型相平面到底怎么工程化落地的人。我会从模型推导开始把Simulink搭建、相平面绘制、参数扫描判断稳定性的完整流程都过一遍最后聊几个我不踩过一次根本意识不到的坑。1. 为什么偏偏看质心侧偏角横摆角速度这两个状态量先说结论车辆横向稳定性最直观的两个表现就是车身有没有甩出去侧偏和有没有转起来横摆支撑这两个表现的恰好就是质心侧偏角和横摆角速度。质心侧偏角β是车辆实际速度方向与车身纵轴之间的夹角。你可以把它理解成车身朝向和实际前进方向的偏差。正常直线行驶时β接近0高速过弯或者低附着路面上一脚重刹β会明显变大车身开始侧滑。横摆角速度r则是车辆绕垂直轴转动的角速度通俗说就是车头转得快不快。这两个量一个描述平移方向一个描述旋转快慢方向互相垂直、又通过轮胎力耦合在一起构成二自由度模型最自然的状态组合。为什么推荐放到相平面里看而不是只看时域响应曲线一个典型的场景给前轮一个阶跃转角仿真出来的β和r曲线可能慢慢稳定也可能发散。曲线能告诉你稳了没但很难告诉你能承受多大的初始扰动不失稳。相平面不一样它把β和r作为坐标轴每个点代表一个瞬时运动状态每个状态点按方程演化出一条轨迹。把所有轨迹画在一起一眼就能看出哪些初始状态会被吸引到稳定平衡点哪些会飞出去。这种全局视野是时域曲线给不了的。这套方法本质上是把非线性动力学的相图分析搬到车辆稳定性问题上。二自由度模型虽然简单但已经能解释车辆失稳的大部分机理尤其是相平面里的平衡点变化和稳定域收缩过程。工程上很多车身稳定控制系统的判稳逻辑最初就是从这类相平面分析里提炼出来的。2. 二自由度模型数学推导与脚本化参数准备2.1 模型假设与运动方程二自由度模型也叫自行车模型核心假设有三条忽略悬架运动车身只做平面运动忽略转向系统惯量前轮转角直接作为输入纵向车速恒定。这样整个车就简化成前轮、后轮两个等效车轮只剩横向运动和横摆运动两个自由度。在车身坐标系下牛顿第二定律和横摆力矩平衡可以写成m·V·(dβ/dt r) Fyf FyrIz·dr/dt a·Fyf - b·Fyr其中m是整车质量V是纵向车速Iz是绕垂直轴的转动惯量a和b分别是质心到前轴、后轴的距离Fyf和Fyr是前、后轴等效侧偏力。轮胎侧偏力在小侧偏角范围内近似线性Fyf -Cf·αfFyr -Cr·αrCf和Cr是前、后轴等效侧偏刚度αf和αr是前、后轮侧偏角推导几何关系和运动学后可以写成状态量表示αf β a·r/V - δαr β - b·r/Vδ是前轮转角。把这几式合并整理就能得到标准的线性二自由度状态空间方程dβ/dt -(CfCr)/(m·V)·β [-1 - (a·Cf - b·Cr)/(m·V²)]·r Cf/(m·V)·δdr/dt -(a·Cf - b·Cr)/Iz·β [-(a²·Cf b²·Cr)/(Iz·V)]·r a·Cf/Iz·δ这里有一个重要关系参数组合(a·Cf - b·Cr)的正负决定了车辆的稳态转向特性。小于0是不足转向大于0是过多转向。这个细节后面讲参数扫描时会专门展示。2.2 参数脚本化先建一个所有人都会改的参数文件既然是可以自己改参数的程序第一步就是把参数集中到一个脚本里而不是散落在Simulink各个模块的参数框中。我习惯建一个vehicle_params.m所有参数在这个文件里统一维护%% 车辆二自由度模型参数 m 1500; % 整车质量 kg Iz 2100; % 横摆转动惯量 kg·m^2 a 1.1; % 质心到前轴距离 m b 1.4; % 质心到后轴距离 m L a b; % 轴距 m Cf 80000; % 前轮等效侧偏刚度 N/rad Cr 90000; % 后轮等效侧偏刚度 N/rad Vx 20; % 纵向车速 m/s delta 0.02; % 前轮转角 rad %% 线性状态空间矩阵 A [-(CfCr)/(m*Vx), -1-(a*Cf-b*Cr)/(m*Vx^2); -(a*Cf-b*Cr)/Iz, -(a^2*Cfb^2*Cr)/(Iz*Vx)]; B [Cf/(m*Vx); a*Cf/Iz]; C eye(2); D [0; 0]; %% 计算特征值快速判断开环稳定性 poles eig(A); disp([当前车速: , num2str(Vx), m/s, 极点: , num2str(poles., 2)]);这里用A、B矩阵作为Simulink状态空间模块的直接输入。好处是改任何物理参数只需要编辑这个脚本再重新运行一次模型里的矩阵自动更新。特征值一行代码就能打印出来开环稳定与否立刻有个大致判断。2.3 模型要扩展非线性轮胎该如何提前考虑线性模型在横向加速度超过约0.4g时就不太可信了因为轮胎侧偏力会饱和。真正做失稳边界分析时建议至少加一个带饱和的非线性轮胎模型。在MATLAB Function块里写一个简化的饱和模型线性区用原来刚度饱和后限制峰值不超过附着极限function [Fyf, Fyr] tire_sat(alphaf, alphar, Fzf, Fzr, mu, Cf, Cr) % 简化轮胎饱和模型线性区 峰值限制 Fyf_lin Cf * alphaf; Fyr_lin Cr * alphar; Fyf max(min(Fyf_lin, mu*Fzf), -mu*Fzf); Fyr max(min(Fyr_lin, mu*Fzr), -mu*Fzr); end先用线性模型把相平面分析的流程跑通确认逻辑没问题再换成这个非线性版本对比一下稳定域的差异你会看到线性模型在低速小转角下非常够用但在低附着、大转角工况下会给出过于乐观的判断。3. Simulink模型搭建与改参数的三种实操路径3.1 方案一State-Space块最省事、最适合快速出结果打开Simulink新建模型从Continuous库拖一个State-Space模块从Sources库拖一个Constant作为前轮转角输入从Sinks库拖两个To Workspace分别命名beta_sim和r_sim。State-Space参数里直接填A:AB:BC:[1 0; 0 1]D:[0; 0]初始条件:[beta0; r0]beta0和r0可以在工作区里提前定义比如跑5组初始状态就预先赋值成向量。先运行vehicle_params.m生成A、B变量再CtrlD运行Simulink模型两个输出就在工作区里了。整个搭建过程不超过五分钟。关于初始条件我建议把工作区的beta0、r0设计成可批量赋值的形式比如beta0 -0.1; r0 0.15;这样每换一组初始条件改脚本重跑一次sim即可。3.2 方案二用积分器和Gain积木式搭建方便后续加非线性State-Space块虽然快但想在里面插入非线性环节就麻烦了。如果后续要换轮胎模型、加状态反馈控制器我更倾向用积分器做积木式搭建。模型结构不复杂核心是两组通路第一条通路算β把β、r、δ三路信号按状态方程第一行的系数缩放后求和送进第一个Integrator积分得到β再从β引出来反馈回路。第二条通路算r用状态方程第二行的系数组合β、r、δ求和后送进第二个Integrator积分得到r。具体操作按下面步骤来拖两个Integrator输出分别命名beta和r。以beta那路为例接一个Sum块输入分别是beta、r、delta三路信号对应的Gain系数分别是-(CfCr)/(m*Vx)、-1-(a*Cf-b*Cr)/(m*Vx^2)、Cf/(m*Vx)。信源用Gain把对应信号缩放后引入Sum。第二个Integrator同理Gain系数分别是-(a*Cf-b*Cr)/Iz、-(a^2*Cfb^2*Cr)/(Iz*Vx)、a*Cf/Iz。两路输出各接一个To Workspace。这个方案最大的优点是可以直接看到状态量的流动方向。想替换轮胎模型时把中间环节换成MATLAB Function块即可输入侧偏角输出侧偏力侧偏角由β、r、δ计算获得。3.3 三种改参数的方式我推荐你怎么选第一种改脚本变量法。适用于参数少、逻辑简单、一个人用的情况。所有参数集中在vehicle_params.m每次改完运行一下脚本再仿真成本最低。我已经习惯把可能调整的参数全部提成文件里的变量不在Simulink块里写死任何数字。第二种Mask封装法。适合把模型发给别人用或者自己频繁调参又不想来回翻脚本。右键选中整个子系统选择Create Mask在Parameters里定义Vx、delta、Cf等输入项用起来就像在用一个带参数对话框的独立组件。这样别人不需要知道模型内部结构也能快速改参数跑数据。第三种脚本批量扫描法。这是做参数研究最需要的模式后面参数扫描实验全靠它。核心思路是不改动模型本身在脚本里循环改工作区变量反复调用sim命令。需要注意State-Space模块里的A、B必须用变量名引用否则循环中改工作区变量不会反映到模型里。%% 批量扫描示例车速从15跑到35 m/s Vx_list [15, 20, 25, 30, 35]; colors lines(length(Vx_list)); figure; hold on; grid on; for i 1:length(Vx_list) Vx Vx_list(i); A [-(CfCr)/(m*Vx), -1-(a*Cf-b*Cr)/(m*Vx^2); -(a*Cf-b*Cr)/Iz, -(a^2*Cfb^2*Cr)/(Iz*Vx)]; B [Cf/(m*Vx); a*Cf/Iz]; assignin(base, A, A); assignin(base, B, B); beta0 -0.1; r0 0.1; % 固定的初始状态 assignin(base, beta0, beta0); assignin(base, r0, r0); sim(vehicle_model_lin); plot(beta_sim*180/pi, r_sim*180/pi, LineWidth, 1.3, Color, colors(i,:)); end legend(15m/s,20m/s,25m/s,30m/s,35m/s); xlabel(质心侧偏角 \beta (deg)); ylabel(横摆角速度 r (deg/s));用assignin把变量写到base工作区是批量仿真时保证Simulink模型能取到最新A、B的关键很多人卡在这一步。4. 相平面绘制的核心实现批量初始条件、方向场、零倾线4.1 方法一批量初始条件叠加轨迹最直观、最常用既然已经建好模型最简单的方式就是定义一组有代表性的初始状态分别仿真再画在同一张图上。初始条件怎么选我习惯按角度比例安排β取-15°到15°的范围r取-25°/s到25°/s的范围低速小转角工况下这些值足以覆盖常见的瞬态扰动范围。%% 批量初始条件相平面图 init_conds [ -0.20, -0.35; % 左前甩尾逆时针横摆 -0.10, -0.20; 0.05, 0.10; 0.00, 0.00; % 平衡点附近 0.10, 0.25; 0.20, 0.40; -0.05, 0.30; % 交叉扰动横摆与侧偏方向相反 0.15, -0.30; ]; figure; hold on; grid on; for i 1:size(init_conds, 1) beta0 init_conds(i, 1); r0 init_conds(i, 2); assignin(base, beta0, beta0); assignin(base, r0, r0); sim(vehicle_model_lin); plot(beta_sim*180/pi, r_sim*180/pi, LineWidth, 1.4); end xlabel(质心侧偏角 \beta (deg)); ylabel(横摆角速度 r (deg/s)); axis([-25 25 -40 40]);注意Simulink里状态单位全是rad、rad/s画图时统一乘180/pi转成角度。这个转换容易漏漏了以后图像比例会严重失衡横纵坐标数值完全对不上物理直觉。4.2 方法二方向场法用矢量箭头展示状态流动趋势轨迹只能看你选中的几条初始状态方向上还有没有盲区方向场补上这个空缺。用meshgrid生成网格点把每个网格点的dβ/dt、dr/dt算出来用quiver画箭头。箭头方向就是状态从该点出发的流动方向哪里汇聚、哪里发散一目了然。%% 方向场绘制 beta_grid deg2rad(-20:2:20); r_grid deg2rad(-30:3:30); [Bg, Rg] meshgrid(beta_grid, r_grid); dB A(1,1)*Bg A(1,2)*Rg B(1)*delta; dR A(2,1)*Bg A(2,2)*Rg B(2)*delta; % 归一化箭头长度避免箭头太长互相重叠 mag sqrt(dB.^2 dR.^2); dB dB ./ mag; dR dR ./ mag; figure; hold on; grid on; quiver(rad2deg(Bg), rad2deg(Rg), dB, dR, 0.6, Color, [0.7 0.7 0.7]); % 继续叠加轨迹 xlabel(质心侧偏角 \beta (deg)); ylabel(横摆角速度 r (deg/s));这个图的作用不只是好看你可以看到整个状态平面里的运动趋势平衡点附近箭头向里螺旋远场箭头如果朝外走说明初始偏差大到一定程度就不回来了稳定域边界大概就在箭头方向发生戏剧性转折的位置。4.3 方法三零倾线标注平衡点用特征值双重校验稳定性零倾线是让dβ/dt0和dr/dt0的两条曲线对于线性系统它们就是两条直线交点正是平衡点。把零倾线叠到相图里图的解读难度能降一大截。%% 零倾线 beta_range deg2rad(-20:0.1:20); r_null_beta -(A(1,1)*beta_range B(1)*delta) / A(1,2); beta_null_r -(A(1,2)*beta_range B(1)*delta) / A(1,1); plot(rad2deg(beta_range), rad2deg(r_null_beta), k--, LineWidth, 1.5); plot(rad2deg(beta_range), rad2deg(beta_null_r), k-., LineWidth, 1.5);线性系统的判断标准很简单特征值实部全为负平衡点是稳定焦点或结点所有初始状态都会收敛到它只要有一个特征值实部为正平衡点变为不稳定鞍点或结点发散不可避免。画相图之前先跑一下eig(A)心里有个预期再对着图看会清楚得多。对于非线性系统零倾线的交点可能不止一个每个交点的局部稳定性需要用雅可比矩阵线性化判断但相图上的全局趋势依然一目了然。5. 参数一改相图就变四组典型对照实验5.1 车速扫描稳定焦点没消失但吸力明显变弱把车速从15逐步升到35 m/s其他参数保持不变能看到的典型现象是相平面里所有轨迹仍朝平衡点收敛但收敛速度越来越慢轨迹螺旋圈数越来越多振荡幅度加大。物理解释其实不难横摆角速度方程里的阻尼项与车速成正比地变小系统等效阻尼下降所以它慢慢从稳定焦点滑向边界稳定。如果这辆车本身是过多转向特性车速超过临界车速后特征值实部会变正轨迹直接螺旋向外发散。这个临界车速在你的参数下是多少用脚本算一下特征值立刻就知道不需要肉眼硬猜。5.2 前轮转角阶跃平衡点绕着稳定域中心移动δ从0.02增加到0.06 rad最直观的变化是平衡点从坐标原点沿r轴正方向移动同时相轨迹也对原点不再对称。原因很直接转角产生一个持续的横摆力矩系统方程的常数输入项把稳定工作点推离了原点。这在工程上意味着转弯工况下车辆的稳定区整体平移了如果还用直线行驶时的β为0来判稳很可能会把正常弯道行驶误判为临界失稳。所以实际ESP逻辑里判稳边界往往是跟转角相关的动态包络而不是固定的一条线。5.3 低附着路面稳定域收缩得让人心里发凉路面附着系数从0.8降到0.3轮胎饱和提前等效侧偏刚度明显下降。这里我一般直接把Cf和Cr减小30%-50%来模拟。相平面上的表现非常扎眼轨迹收敛明显变慢同一初始状态在干沥青路面能轻松回到平衡点在冰雪路面上却跑出大半个圈才回来稍大一点的初始β直接螺旋发散。这个实验是我给学生演示时最常做的一步它直观解释了为什么雪天开车感觉方向没变车却不听话——稳定域缩小了同样的扰动幅度在低附着工况已经飞出相图边界。5.4 质心位置后移从不足转向滑向过多转向把质心往后挪让a变大、b变小导致a*Cf - b*Cr从负变正整车转向特性从不足转向变为过多转向。相图变化最剧烈低速时一切看起来还是正常的但车速一提高特征值实部很快跨过0平衡点失稳所有轨迹都会向外螺旋发散。这就是经典教科书里过多转向车辆超过临界车速后失去开环稳定性的相平面版本。你不需要真的改造一辆车只要改一下a和b参数就能在Simulink里把失稳全过程看透。参数调整相平面现象物理解释车速提高轨迹收敛变慢振荡加剧过多转向时发散系统阻尼随速度下降前轮转角增大平衡点沿r轴偏移相图不对称常量输入改变了工作点附着系数降低稳定域收缩初始偏差容易发散轮胎等效刚度降低饱和提前质心后移转向特性转为过多转向高速失稳aCf-bCr符号反转6. 从相平面到实际稳定性判断以及调试中的几个坑6.1 相平面结果怎么投射到工程判稳指标相平面图做出来不能只停在好看上。实际做车辆稳定性控制时最常参考的就是质心侧偏角门限一般干沥青路面β超过10°到12°就已经相当危险低附着路面这个门限还要更低。在相平面图上这个门限对应一条竖直的边界线轨迹一旦超过它就触发预警或介入控制。另外一个常用做法是画一个包络区域比如以平衡点为中心的椭圆或平行四边形轨迹完全落在包络内认为是稳定。这个包络边界可以通过反向仿真近似获得从一组处于临界状态的边界点倒着积分得到真正稳定域的边界曲线。反向轨迹逼近不稳定流形这个方法很有用尤其是非线性模型下稳定域形状不规则时比拍脑袋画椭圆靠谱得多。6.2 四个我实际踩过的坑第一个坑是单位混用。Simulink模型里我用rad和rad/s画图时转成deg和deg/s但初始条件设置时经常忘了乘pi/180导致轨迹方向完全错乱。检查方法和建议在vehicle_params.m开头写一行注释统一标注所有状态的单位或者直接在变量名里带_rad后缀别嫌麻烦。第二个坑是求解器设置。默认变步长ode45跑简单工况没问题但批量初始状态仿真时某些初始点会导致积分步长来回跳仿真时间拖得很长。我自己更习惯在批量扫描时固定步长用ode4步长设0.001秒每个工况跑3到5秒速度和稳定性都兼顾。注意改变步长后数值阻尼也会变不同题不要拿不同步长跑出来对比。第三个坑是线性模型别拿超大初始状态跑。β给到0.5 rad、r给到3 rad/s轨迹大幅发散看起来像失稳实际上早就超出线性轮胎的有效范围了模型本身已经失真。线性模型的初始扰动范围建议控制在β±10°、r±20°/s以内超出的工况改用非线性轮胎模型再分析。第四个坑是批量扫描时模型不刷新参数。用工作区变量填A、B后如果直接在循环里改工作区的值就调用sim个别版本模型会沿用上一次仿真时的旧矩阵。我的解决方式是显式用assignin(base, ...)写入并在sim调用前加一条set_param(vehicle_model_lin, SimulationCommand, update)强制模型同步一次实测稳定可靠。6.3 让程序更好用的小改进方向如果嫌手动改初始条件麻烦可以在vehicle_params.m里定义一组默认的工况场景比如低速干沥青高速低附着紧急双移线三种预置参数集用switch结构一键切换。我自己的版本就是这么组织的写论文做对照实验效率会高不少。另外可以把方向场、轨迹、零倾线、特征值四样东西全部塞进一个绘制函数输入参数只需要从工作区读一组当前参数输出就是一张完整的相平面图。这样后续做蒙地卡罗参数采样或者批量工况扫描时可以直接复用函数不需要每张图重新写一遍绘图代码。我自己调这套程序时最常干的事就是不停切换Vx、delta、附着条件反复刷新相图看到轨迹从稳定收敛变成螺旋发散再回头对一下特征值整车失稳机理的理解就会变得非常具象。如果你手里有一套自己的Simulink模型建议把相平面分析作为标准工具焊死在项目流程里准不会后悔。