简介：资源包聚焦于车辆二自由度动态模型的搭建与分析，面向汽车工程及自动驾驶方向的开发者、研究人员。模型从状态空间方程出发，转换为传递函数形式，重点考察横摆角速度与车辆侧偏角对操控稳定性的影响。包体共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_f+C_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 | 整车质量 | 1500 | kg |
| Iz | 横摆转动惯量 | 2500 | kg·m² |
| lf | 质心到前轴距离 | 1.2 | m |
| lr | 质心到后轴距离 | 1.4 | m |
| Cf | 前轴等效侧偏刚度 | 45000 | N/rad |
| Cr | 后轴等效侧偏刚度 | 45000 | N/rad |
| Vx | 纵向车速 | 20 | m/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, m=1500.0, Iz=2500.0, lf=1.2, lr=1.4, Cf=45000.0, Cr=45000.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) | ζ |
|---|---|---|---|
| 10 | 3.50 | 6.33 | 0.96 |
| 20 | 5.52 | 3.57 | 0.85 |
| 30 | 6.11 | 2.77 | 0.73 |
| 40 | 5.97 | 2.43 | 0.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, T=T) t_tf, y_tf = ct.step_response(sys_tf[0, 0], T=T) 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 对得上，再谈下游的增益调度。
本文还有配套的精品资源，点击获取
