别再死记硬背公式了！用Python/Simulink动手推导车辆动力学方程（附代码）

寒泊

101人浏览 · 2026-05-01 09:09:08

寒泊 · 2026-05-01 09:09:08 发布

用Python实战推导车辆动力学方程：从理论到代码的完整指南

在车辆工程和自动驾驶算法开发中，理解车辆动力学方程是核心基础。传统教材往往只给出最终公式，却很少展示这些方程是如何一步步推导出来的。本文将带你用Python的SymPy库，从零开始推导经典的二自由度车辆模型，并通过可视化理解轮胎侧偏特性。

1. 理解轮胎侧偏特性的物理本质

轮胎侧偏特性是车辆动力学中最关键也最容易被误解的概念之一。当车辆转弯时，轮胎实际运动方向与指向方向之间会产生一个夹角，这就是侧偏角(α)。侧偏力(Fy)则是地面作用于轮胎的横向力，与侧偏角密切相关。

关键物理现象：

刚性轮胎与弹性轮胎的侧偏行为完全不同
小侧偏角时，侧偏力与侧偏角近似线性关系
大侧偏角时，侧偏力趋于饱和，呈现非线性特性

用Python绘制典型的侧偏力-侧偏角曲线：

import numpy as np
import matplotlib.pyplot as plt

def magic_formula(alpha, B, C, D, E):
    return D * np.sin(C * np.arctan(B * alpha - E * (B * alpha - np.arctan(B * alpha))))

alpha = np.linspace(-0.3, 0.3, 100)  # 侧偏角范围(-0.3到0.3弧度)
Fy_linear = -50000 * alpha  # 线性区域
Fy_nonlinear = magic_formula(alpha, 10, 1.3, 8000, 0.5)  # 魔术公式模型

plt.figure(figsize=(10,6))
plt.plot(alpha, Fy_linear, '--', label='线性近似')
plt.plot(alpha, Fy_nonlinear, label='魔术公式模型')
plt.xlabel('侧偏角α (rad)')
plt.ylabel('侧偏力Fy (N)')
plt.title('轮胎侧偏特性曲线')
plt.grid()
plt.legend()
plt.show()

2. 建立二自由度车辆模型

二自由度车辆模型考虑车辆的横向运动和横摆运动，是分析车辆动态特性的基础模型。我们需要建立车身坐标系并定义关键参数：

模型参数表：

参数	描述	典型值
m	车辆质量	1500 kg
a	前轴到质心距离	1.2 m
b	后轴到质心距离	1.5 m
Cαf	前轮侧偏刚度	-50000 N/rad
Cαr	后轮侧偏刚度	-60000 N/rad
Iz	横摆转动惯量	3000 kg·m²
vx	纵向速度	20 m/s

使用SymPy进行符号推导：

from sympy import symbols, Matrix, cos, sin, Eq, solve, simplify

# 定义符号变量
m, a, b, C_af, C_ar, I_z, v_x = symbols('m a b C_af C_ar I_z v_x')
y_dot, phi_dot, y_ddot, phi_ddot, delta = symbols('y_dot phi_dot y_ddot phi_ddot delta')

# 计算前后轮侧偏角
alpha_f = (y_dot + a*phi_dot)/v_x - delta
alpha_r = (y_dot - b*phi_dot)/v_x

# 侧偏力计算
F_yf = C_af * alpha_f
F_yr = C_ar * alpha_r

# 动力学方程
eq1 = Eq(m*(y_ddot + v_x*phi_dot), F_yf + F_yr)
eq2 = Eq(I_z*phi_ddot, a*F_yf - b*F_yr)

# 解方程
solution = solve([eq1, eq2], [y_ddot, phi_ddot])
y_ddot_expr = simplify(solution[y_ddot])
phi_ddot_expr = simplify(solution[phi_ddot])

print(f"y_ddot = {y_ddot_expr}")
print(f"phi_ddot = {phi_ddot_expr}")

3. 推导状态空间方程

将二阶微分方程转化为状态空间形式是控制系统分析和设计的基础。我们定义状态向量X = [y_dot, phi_dot]^T，控制输入u = δ(前轮转角)。

推导得到的矩阵形式为：

Ẋ = A·X + B·u

其中系统矩阵A和控制矩阵B为：

# 状态空间矩阵
A = Matrix([
    [(C_af + C_ar)/(m*v_x), (a*C_af - b*C_ar)/(m*v_x) - v_x],
    [(a*C_af - b*C_ar)/(I_z*v_x), (a**2*C_af + b**2*C_ar)/(I_z*v_x)]
])

B = Matrix([
    [-C_af/m],
    [-a*C_af/I_z]
])

print("系统矩阵A:")
print(A)
print("\n控制矩阵B:")
print(B)

状态空间方程的物理意义：

A矩阵描述了车辆自身的动力学特性
B矩阵反映了前轮转角对系统的影响
该形式便于后续进行稳定性分析和控制器设计

4. 模型验证与动态特性分析

推导出方程后，我们需要验证模型的正确性并分析其动态特性。使用Python进行数值仿真：

from scipy.integrate import odeint

def vehicle_dynamics(X, t, A, B, delta):
    return A @ X + B * delta

# 参数值
params = {
    'm': 1500, 'a': 1.2, 'b': 1.5,
    'C_af': -50000, 'C_ar': -60000,
    'I_z': 3000, 'v_x': 20
}

# 计算数值化的A,B矩阵
A_num = np.array([
    [(params['C_af']+params['C_ar'])/(params['m']*params['v_x']), 
     (params['a']*params['C_af']-params['b']*params['C_ar'])/(params['m']*params['v_x'])-params['v_x']],
    [(params['a']*params['C_af']-params['b']*params['C_ar'])/(params['I_z']*params['v_x']),
     (params['a']**2*params['C_af']+params['b']**2*params['C_ar'])/(params['I_z']*params['v_x'])]
])

B_num = np.array([
    [-params['C_af']/params['m']],
    [-params['a']*params['C_af']/params['I_z']]
])

# 仿真步进转向输入
t = np.linspace(0, 5, 500)
delta_input = np.where(t < 1, 0, 0.05)  # 1秒后施加5度转向角

# 数值积分
X0 = np.zeros(2)  # 初始状态
X = np.zeros((len(t), 2))
for i in range(1, len(t)):
    dt = t[i] - t[i-1]
    X[i] = X[i-1] + vehicle_dynamics(X[i-1], t[i-1], A_num, B_num, delta_input[i-1]) * dt

# 绘制结果
plt.figure(figsize=(12,8))
plt.subplot(3,1,1)
plt.plot(t, delta_input*180/np.pi)
plt.ylabel('前轮转角(deg)')
plt.grid()

plt.subplot(3,1,2)
plt.plot(t, X[:,0])
plt.ylabel('横向速度(m/s)')
plt.grid()

plt.subplot(3,1,3)
plt.plot(t, X[:,1]*180/np.pi)
plt.ylabel('横摆角速度(deg/s)')
plt.xlabel('时间(s)')
plt.grid()
plt.tight_layout()
plt.show()

动态特性分析要点：