简介资源包聚焦于车辆二自由度动态模型的搭建与分析面向汽车工程及自动驾驶方向的开发者、研究人员。模型从状态空间方程出发转换为传递函数形式重点考察横摆角速度与车辆侧偏角对操控稳定性的影响。包体共2个文件分别为MATLAB脚本与FIG图形文件脚本负责状态空间至传递函数的转换与仿真计算图形文件展示波特图、阶跃响应等可视化结果压缩包仅74KB结构精简便于快速部署。当前已有642人学习下载使用者可借助脚本自定义输入条件观察不同前轮转角或车速下的车辆动态响应评估急转向、紧急避障等工况的稳定性也可为横摆力矩控制器的参数整定提供参考。该文件包虽小但完整覆盖从建模、转换到响应分析的核心流程适合车辆控制入门者与工程师快速验证。1. 车辆二自由度模型的传递函数到底在解决什么问题做底盘电控或辅助驾驶横向控制时最容易碰到的场景是前轮转角给下去车身先建立侧偏角再慢慢形成横摆角速度两个响应耦合在一起而且随车速变化非常大。如果直接用整车的非线性动力学模型去看很难分清哪个环节在限制带宽哪个环节在贡献相位延迟。dof2_tf这类二自由度传递函数工具的用处就是先把车辆真正最核心的横向与横摆运动提炼成两个可解析的传递函数前轮转角到横摆角速度、前轮转角到车辆侧偏角。这篇文章不打算讲泛泛的“车辆动力学”而是直接把自行车模型的状态方程写出来推导出传递函数然后再用 Python 把增益、阻尼、零点和首一/尾一归一化这些工程上最容易踩坑的点全部算一遍。适合正在做车辆稳定性分析、EPS/ESP 控制标定或者模拟器整车模型简化的工程师阅读。2. 车辆二自由度模型的状态方程与横摆传递函数推导2.1 用自行车模型写线性动力学方程二自由度模型通常指“自行车模型”忽略左右轮差异把前后轴各自合并成一个等效轮胎假设纵向车速Vx恒定只保留侧向速度vy和横摆角速度r两个自由度。为了直接观察质心侧偏角这里用状态变量x [β, r]^T其中β vy / Vx。前轮转角δ_f是输入。在这个假设下侧向力平衡方程和横摆力矩方程可以写成$$ m V_x (\dot{\beta} r) C_f(\delta_f - \beta - \frac{l_f r}{V_x}) C_r(- \beta \frac{l_r r}{V_x}) $$$$ I_z \dot{r} l_f C_f(\delta_f - \beta - \frac{l_f r}{V_x}) - l_r C_r(- \beta \frac{l_r r}{V_x}) $$其中m是整车质量Iz是绕 Z 轴的横摆转动惯量lf、lr分别是质心到前轴和后轴的距离C_f、C_r是前、后轴的等效侧偏刚度都取正值。把两个方程整理成状态空间形式得到$$ \begin{bmatrix} \dot{\beta}\ \dot{r} \end{bmatrix}\begin{bmatrix} -\frac{C_fC_r}{mV_x} -1 - \frac{C_f l_f - C_r l_r}{mV_x^2}\ -\frac{C_f l_f - C_r l_r}{I_z} -\frac{C_f l_f^2 C_r l_r^2}{I_z V_x} \end{bmatrix} \begin{bmatrix} \beta\ r \end{bmatrix} \begin{bmatrix} \frac{C_f}{mV_x}\ \frac{l_f C_f}{I_z} \end{bmatrix} \delta_f $$注意矩阵第一行第二列里的-1来自侧向力方程中的离心项m Vx r这一项是最容易被忽略的。如果把它丢掉高速时横摆角速度直流增益会明显偏大。输出选择y [r, β]^T也就是第二行横摆角速度、第一行侧偏角所以输出矩阵为$$ C_{out} \begin{bmatrix} 0 1\ 1 0 \end{bmatrix}, \quad D 0 $$这里用C_out这个符号是为了和侧偏刚度C_f区分。实际代码中矩阵名不要直接用C否则很容易在调试时把两个含义弄混。2.2 传递函数分母与分子二阶系统加一个零点用G(s) C_out (sI-A)^{-1} B求传递函数。因为系统只有两个状态所以特征多项式一定是二阶$$ \Delta(s) s^2 2\zeta \omega_n s \omega_n^2 $$其中ωn是固有圆频率ζ是阻尼比。由于状态矩阵A中含有1/Vx项这两个指标都会随车速变化。横摆角速度对前轮转角的传递函数分子通常不是常数而是一个带正实零点的多项式。也就是说r/δ_f是“一个零点 两个极点”的系统零点会让高频相位提前这也是为什么实际车辆在快速转向时横摆响应比简单的二阶滞后模型更快一些。车辆侧偏角对前轮转角的传递函数分子同样含有一个零点但其常数项往往是负值这代表稳态时质心侧偏角与转向方向相反。直接用sympy做一次符号推导可以很清楚地看到这种结构import sympy as sp s sp.symbols(s) m, Iz, lf, lr, Cf, Cr, Vx sp.symbols(m Iz lf lr Cf Cr Vx) A sp.Matrix([ [-(Cf Cr) / (m * Vx), -1 - (Cf * lf - Cr * lr) / (m * Vx**2)], [-(Cf * lf - Cr * lr) / Iz, -(Cf * lf**2 Cr * lr**2) / (Iz * Vx)] ]) B sp.Matrix([[Cf / (m * Vx)], [lf * Cf / Iz]]) C_out sp.Matrix([[0, 1], [1, 0]]) G sp.simplify(C_out * (s * sp.eye(2) - A).inv() * B) print(sp.factor(G[0, 0])) # r/delta_f print(sp.factor(G[1, 0])) # beta/delta_f这段代码会输出两个传递函数的符号表达式。可以看到分子中确实都存在s的一次项而不是纯比例或纯积分形式。符号表达式的意义在于当你改变质心位置或前后侧偏刚度比例时零点会如何移动一眼就能看出来这比纯数值调参有用得多。2.3 从直流增益看不足转向与过度转向稳态横摆角速度增益指的就是s→0时的传递函数值工程上通常叫“横摆增益”或“转向灵敏度”。它可以从状态空间直接求直流增益$$ G_r(0) - C_{out, r} A^{-1} B $$也可以解稳态代数方程得到。对二自由度模型来说横摆增益随车速的变化会呈现出“先线性增大后增速放缓最终下降”的形状。峰值出现在什么车速取决于不足转向系数如果峰值一直不出现直到临界车速前都在上升那就是中性转向或过度转向性质。先给出一组真实的整车参数后面所有代码和表格都基于这组参数计算。这是一台典型 B 级轿车的等效参数参数符号含义数值单位m整车质量1500kgIz横摆转动惯量2500kg·m²lf质心到前轴距离1.2mlr质心到后轴距离1.4mCf前轴等效侧偏刚度45000N/radCr后轴等效侧偏刚度45000N/radVx纵向车速20m/s注意这里的C_f、C_r是整轴等效侧偏刚度不是单胎值通常来自轮胎试验数据换算后的结果。如果直接把单胎刚度乘 2 放进去往往会让车辆看起来比实际偏过度。3. 用 Python 把 dof2_tf 传递函数算出来3.1 最小状态空间实现与 ss2tf 转换现在用control库把上一章的状态空间矩阵变成传递函数。Python 环境需要安装numpy和control推荐再装一个matplotlib方便后面画阶跃响应。先实现最小可用版本import numpy as np import control as ct # 参数来自表 2-1 m 1500.0 Iz 2500.0 lf 1.2 lr 1.4 Cf 45000.0 Cr 45000.0 Vx 20.0 # 状态向量 [beta, r]输入 delta_f A np.array([ [-(Cf Cr) / (m * Vx), -1.0 - (Cf * lf - Cr * lr) / (m * Vx * Vx)], [-(Cf * lf - Cr * lr) / Iz, -(Cf * lf * lf Cr * lr * lr) / (Iz * Vx)] ]) B np.array([[Cf / (m * Vx)], [lf * Cf / Iz]]) # 输出顺序第一行 r第二行 beta C_out np.array([[0.0, 1.0], [1.0, 0.0]]) D_out np.zeros((2, 1)) sys_ss ct.ss(A, B, C_out, D_out) sys_tf ct.ss2tf(sys_ss) for idx, name in [(0, r/delta_f), (1, beta/delta_f)]: print(name, sys_tf[idx, 0])在 Vx 20 m/s 时输出大约为r/delta_f (21.6 s 70.2) / (s^2 6.06 s 12.7) beta/delta_f (1.5 s - 16.7) / (s^2 6.06 s 12.7)说明这段代码里的几个关键点矩阵A第一行第二列的-1后面还有一项-(Cf*lf - Cr*lr)/(mVx^2)因为这里前、后轴侧偏刚度相同但轴距载荷分配不同所以综合结果是0.015最终让A[0,1]变成约-0.985。如果你把这一项漏掉横摆角速度的直流增益会上偏约 3% 左右。C_out的行顺序直接决定ss2tf输出哪一行是r哪一行是β。很多人调试时发现两个传递函数对调就是因为这里的行顺序写反了。建议在注释里写清楚输出顺序。整理成尾一形式后r/δ_f的直流增益为$$ \frac{70.2}{12.7} 5.52 \ \text{rad/s per rad} $$这个数值意味着前轮转角给 0.1 rad约 5.7°时稳态横摆角速度大约为 0.55 rad/s。3.2 首一还是尾一开环传递函数增益藏在常数项比值里control库返回的传递函数默认是“首一”形式也就是分子分母都除以最高次项系数使分母最高次项系数为 1。比如上面的传递函数中分母已经写成s^2 6.06s 12.7最高次项系数就是 1。但工程上做增益调度时更习惯用“尾一”形式也就是分母常数项为 1。因为直流增益直接等于分子常数项与分母常数项的比值。两种归一化方式本身没有对错但换算错了会导致开环增益差好几个数量级。把首一形式转成尾一形式的代码是num np.asarray(sys_tf[0, 0].num[0][0], float) # 首一分子系数 den np.asarray(sys_tf[0, 0].den[0][0], float) # 首一分母系数 # 直流增益 分子常数项 / 分母常数项 dc_gain num[-1] / den[-1] print(直流增益 , dc_gain) # 约 5.516 # 转尾一分子分母同除以分母常数项 num_tail num / den[-1] den_tail den / den[-1] print(尾一 num:, num_tail) # 约 [1.697, 5.516] print(尾一 den:, den_tail) # 约 [0.0786, 0.476, 1.0]转换之后横摆角速度传递函数变成$$ \frac{r}{\delta_f}(s) \frac{1.697s 5.516}{0.0786s^2 0.476s 1} $$这时分子常数项就是直流增益5.516而分母常数项是 1。两种形式描述同一个系统但如果你拿着首一形式去找“开环传递函数增益”很可能会把70.2当成表调增益或者在 matlab 里用dcgain时忘记除以分母常数项。建议把这段转换逻辑封装成一个公共函数所有增益表都从尾一形式读取。3.3 批量计算不同速度下的传函参数实际工程中不会只算一个速度。把上面的逻辑包进函数循环计算 10、20、30、40 m/s 下的直流增益、固有频率和阻尼比def vehicle_tf(Vx, m1500.0, Iz2500.0, lf1.2, lr1.4, Cf45000.0, Cr45000.0): A np.array([ [-(Cf Cr) / (m * Vx), -1.0 - (Cf * lf - Cr * lr) / (m * Vx * Vx)], [-(Cf * lf - Cr * lr) / Iz, -(Cf * lf * lf Cr * lr * lr) / (Iz * Vx)] ]) B np.array([[Cf / (m * Vx)], [lf * Cf / Iz]]) sys ct.ss(A, B, np.array([[0.0, 1.0]]), np.zeros((1, 1))) return ct.ss2tf(sys) for V in [10, 20, 30, 40]: G vehicle_tf(V) num np.asarray(G[0, 0].num[0][0], float) den np.asarray(G[0, 0].den[0][0], float) dc num[-1] / den[-1] wn np.sqrt(den[-1]) zeta -den[1] / (2.0 * wn) print(f{V:2.0f} m/s: dc{dc:.3f}, wn{wn:.2f}, zeta{zeta:.2f})注意这里求阻尼比使用的是den[1]因为首一形式下numpy数组的索引 0 是s^2系数索引 1 是s系数索引 2 是常数项。这是信号处理库的常见排列顺序和手写多项式时的习惯不一样容易看错。4. 阶跃响应与转向特性分析横摆角速度增益怎么变4.1 速度对横摆增益与阻尼比的影响用上一节代码跑出来的结果整理成表格Vx (m/s)直流增益 r/δ (rad/s per rad)ωn (rad/s)ζ103.506.330.96205.523.570.85306.112.770.73405.972.430.62横向对比这组数据可以得出几个工程判断低速时横摆角速度直流增益随车速近似线性增大驾驶员会觉得转向“灵敏”。到了 30 m/s 以后增益不再继续上升这说明车辆已经有明显的不足转向趋势。如果增大后轴侧偏刚度Cr峰值速度会降低不足转向量变大减小Cr则相反车辆会变得更灵活但阻尼比下降速度会加快。阻尼比从 0.96 降到 0.62意味着高速情况下横摆振荡衰减变慢。你可以在阶跃响应里看到第一个超调量增大。对于线性模型而言只要阻尼比大于 1 就没有超调在 40 m/s 时虽然还有 0.62但实际车辆因为轮胎非线性衰减会更快所以不能只看这个线性阻尼比来标定转向手感。4.2 车辆侧偏角响应低速小、高速大的物理原因车辆侧偏角 β 的稳态值也是从传递函数直流增益里直接读。使用上一节的vehicle_tf函数把输出矩阵改为只输出 β 即可。从传递函数本身可以看出β 的分子常数项是负值所以稳态 β 与转向方向符号相反。前轮转向时车头往内转质心速度方向却对外偏β 为负说明质心“跟不上”车头方向。高速时这个负值会明显增大因为侧向力建立得更慢需要更大的侧偏角来产生足够的横摆力矩。如果根据参数表手动计算Vx 20 m/s 时 β/δ_f 的直流增益约为-1.31 rad/rad即前轮转角给 0.1 rad 时稳态侧偏角约为 -0.13 rad约 7.5°。这在真实车辆上已经算很大的侧偏角了说明该参数下车辆在高速紧急变道时稳定性余量较小。在 ESP/CMS 标定中常用质心侧偏角速度dβ/dt作为触发信号因为 β 本身变化较慢dβ/dt能更早暴露失稳趋势。用传递函数可以很容易求出β的导数的频域响应也就是在原传递函数前面乘一个s这比做差分滤波更干净。4.3 零点、极点与转向灵敏度的统一理解横摆角速度传递函数的零点随速度变化但它始终是一个负实数零点所以系统是“最小相位”的不会有初期反向响应。零点位置越接近虚轴相位提前作用越明显。以 Vx 20 m/s 为例零点在约-3.25 rad/s和固有频率3.57 rad/s很接近所以相位提前在中频段已经很显著。用以下代码可以同时打印零点和极点num np.asarray(sys_tf[0, 0].num[0][0], float) den np.asarray(sys_tf[0, 0].den[0][0], float) zeros np.roots(num) poles np.roots(den) print(零点:, zeros) print(极点:, poles) print(阻尼比:, -poles.real / np.abs(poles))如果某个参数组合下极点实部出现在右半平面直接说明该速度下车队是开环不稳定的。这个分析完全基于传递函数不需要做非线性仿真适合在概念设计阶段快速筛选参数范围。5. 用状态空间数值积分反查传递函数守住模型边界5.1 阶跃响应对照tf 和 lsim 必须一致ss2tf转换过程中涉及矩阵求逆和多项式化简代码里最容易出错的是输出矩阵行顺序。一个很有效的验证方法是用同一个状态空间模型分别做step_response和lsim再把结果和传递函数的step_response对比最大误差应接近机器精度。import numpy as np import control as ct # 复用前面的 sys_ss取第一输出 r sys_r_ss ct.ss(A, B, np.array([[0.0, 1.0]]), np.zeros((1, 1))) T np.linspace(0, 2.0, 1000) t_ss, y_ss ct.step_response(sys_r_ss, TT) t_tf, y_tf ct.step_response(sys_tf[0, 0], TT) err np.max(np.abs(y_ss[0, :] - y_tf[0, :])) print(最大误差:, err)这个误差如果超过 1e-8基本可以确定状态矩阵某个元素或输出矩阵写错了。也可以换成长度足够长的正弦输入用lsim对比频域响应能覆盖更多频率点。5.2 低速与高频边界二自由度传递函数在低速段和高速段都有自己的使用边界。车速低于 5 m/s 时侧偏角本身很小轮胎基本没有多少侧向力模型中的1/Vx项会迅速放大导致传递函数数值对参数极度敏感。此时更适合用运动学自行车模型。高频段方面线性轮胎模型没有考虑轮胎松弛长度也没有考虑悬架侧倾转向所以当输入频率超过悬架共振频率后真实车辆的横摆响应会出现额外相位滞后。用二自由度传递函数做控制设计时建议把有效带宽限制在 3 Hz 以内如果整车频响在这一段出现明显的 -60 dB/dec 斜率说明垂向和悬架自由度已经参与进来不能再简化成二自由度系统。把这段对照验证代码放进你的回归脚本里每次改一个参数就跑一遍确认直流增益和零极点都和表 4-1 对得上再谈下游的增益调度。本文还有配套的精品资源点击获取