用Python和SymPy搞定汽车二自由度模型：从理论方程到代码仿真（保姆级教程）-编程阁

用Python和SymPy搞定汽车二自由度模型：从理论方程到代码仿真（保姆级教程）

当你在车辆动力学教材中看到那些复杂的微分方程时，是否曾想过如何将它们变成可运行的代码？本文将带你从零开始，用Python实现经典的汽车二自由度模型仿真。不同于纯理论推导，我们会重点关注如何将数学公式转化为可执行的计算机程序，并直观地观察车辆动态响应。

1. 环境准备与基础概念

在开始编码之前，我们需要明确几个关键概念。二自由度模型主要研究车辆的侧向运动和横摆运动，这是分析车辆操纵稳定性的基础。为了简化问题，模型做了以下假设：

忽略悬架作用（无垂直、侧倾和俯仰运动）
恒定前进速度（无纵向加速度）
轮胎侧偏特性处于线性范围

安装必要的Python库：

pip install sympy numpy matplotlib scipy

核心工具介绍：

SymPy：符号计算库，用于处理微分方程
NumPy：数值计算基础
SciPy：提供ODE求解器
Matplotlib：结果可视化

提示：建议使用Jupyter Notebook进行交互式开发，方便实时查看计算结果和图形。

2. 建立数学模型

2.1 运动学方程推导

考虑车辆在平面内的运动，我们定义以下变量：

u：纵向速度（假设恒定）
v：侧向速度
r：横摆角速度

车辆坐标系下的加速度分量可表示为：

a_y = dv/dt + u*r

2.2 动力学方程建立

根据牛顿第二定律，建立力和力矩平衡方程：

参数	符号	说明
前轮侧偏角	α_f	前轮速度方向与轮胎指向的夹角
后轮侧偏角	α_r	后轮速度方向与轮胎指向的夹角
前轮侧偏刚度	C_f	前轮产生单位侧偏角所需的力
后轮侧偏刚度	C_r	后轮产生单位侧偏角所需的力

侧向力方程：

m*(dv/dt + u*r) = F_yf + F_yr

横摆力矩方程：

I_z*dr/dt = a*F_yf - b*F_yr

其中：

F_yf = -C_f * α_f
F_yr = -C_r * α_r
α_f ≈ δ - (v + a*r)/u
α_r ≈ -(v - b*r)/u

2.3 状态空间表示

将方程整理为矩阵形式：

[ dv/dt ] = A * [ v ] + B * [ δ ] [ dr/dt ] [ r ]

其中A和B是系统矩阵，δ是前轮转向角输入。

3. Python实现

3.1 符号推导

使用SymPy进行符号推导：

from sympy import symbols, Matrix, Eq, solve # 定义符号变量 v, r, δ = symbols('v r delta') u, m, Iz, Cf, Cr, a, b = symbols('u m Iz Cf Cr a b', positive=True) # 建立方程 αf = δ - (v + a*r)/u αr = -(v - b*r)/u Fyf = -Cf * αf Fyr = -Cr * αr # 动力学方程 eq1 = Eq(m*(v.diff() + u*r), Fyf + Fyr) eq2 = Eq(Iz*r.diff(), a*Fyf - b*Fyr) # 解耦微分方程 solution = solve((eq1, eq2), (v.diff(), r.diff())) dvdt = solution[v.diff()] drdt = solution[r.diff()]

3.2 数值求解

将符号表达式转换为数值函数：

import numpy as np from scipy.integrate import odeint def vehicle_model(state, t, u_val, m_val, Iz_val, Cf_val, Cr_val, a_val, b_val, delta_func): v, r = state δ = delta_func(t) # 参数替换字典 subs_dict = { 'u': u_val, 'm': m_val, 'Iz': Iz_val, 'Cf': Cf_val, 'Cr': Cr_val, 'a': a_val, 'b': b_val, 'delta': δ, 'v': v, 'r': r } dvdt_val = float(dvdt.subs(subs_dict)) drdt_val = float(drdt.subs(subs_dict)) return [dvdt_val, drdt_val]

3.3 仿真参数设置

典型轿车参数示例：

# 车辆参数 params = { 'u_val': 20.0, # 车速 [m/s] 'm_val': 1500.0, # 质量 [kg] 'Iz_val': 2500.0, # 横摆转动惯量 [kg·m²] 'Cf_val': 80000.0, # 前轮侧偏刚度 [N/rad] 'Cr_val': 100000.0, # 后轮侧偏刚度 [N/rad] 'a_val': 1.2, # 前轴到质心距离 [m] 'b_val': 1.5 # 后轴到质心距离 [m] } # 转向输入函数 def step_steering(t): return 0.1 if t >= 1.0 else 0.0 # 1秒后施加10°转向角 # 时间点 t = np.linspace(0, 5, 500) # 0到5秒，500个点

4. 仿真与结果分析

4.1 运行仿真

# 初始状态 state0 = [0.0, 0.0] # 初始侧向速度和横摆角速度为0 # 求解ODE states = odeint( vehicle_model, state0, t, args=( params['u_val'], params['m_val'], params['Iz_val'], params['Cf_val'], params['Cr_val'], params['a_val'], params['b_val'], step_steering ) ) v_sim = states[:, 0] # 侧向速度 r_sim = states[:, 1] # 横摆角速度

4.2 结果可视化

绘制响应曲线：

import matplotlib.pyplot as plt plt.figure(figsize=(12, 8)) # 侧向速度 plt.subplot(2, 1, 1) plt.plot(t, v_sim, label='侧向速度 [m/s]') plt.ylabel('v [m/s]') plt.grid(True) plt.legend() # 横摆角速度 plt.subplot(2, 1, 2) plt.plot(t, r_sim, label='横摆角速度 [rad/s]') plt.xlabel('时间 [s]') plt.ylabel('r [rad/s]') plt.grid(True) plt.legend() plt.tight_layout() plt.show()

4.3 轨迹计算

# 计算车辆轨迹 psi = np.cumsum(r_sim) * (t[1]-t[0]) # 横摆角积分 X = np.cumsum(params['u_val'] * np.cos(psi) - v_sim * np.sin(psi)) * (t[1]-t[0]) Y = np.cumsum(params['u_val'] * np.sin(psi) + v_sim * np.cos(psi)) * (t[1]-t[0]) plt.figure(figsize=(8, 6)) plt.plot(X, Y) plt.xlabel('X位置 [m]') plt.ylabel('Y位置 [m]') plt.title('车辆轨迹') plt.grid(True) plt.axis('equal') plt.show()

5. 参数影响分析

5.1 侧偏刚度影响

Cf_values = [40000, 80000, 120000] # 不同的前轮侧偏刚度 plt.figure(figsize=(12, 6)) for Cf in Cf_values: states = odeint( vehicle_model, state0, t, args=( params['u_val'], params['m_val'], params['Iz_val'], Cf, params['Cr_val'], params['a_val'], params['b_val'], step_steering ) ) plt.plot(t, states[:, 1], label=f'Cf={Cf/1000:.0f} kN/rad') plt.xlabel('时间 [s]') plt.ylabel('横摆角速度 [rad/s]') plt.title('不同前轮侧偏刚度下的响应') plt.legend() plt.grid(True) plt.show()

5.2 质心位置影响

a_values = [0.8, 1.2, 1.6] # 不同的前轴到质心距离 plt.figure(figsize=(12, 6)) for a in a_values: states = odeint( vehicle_model, state0, t, args=( params['u_val'], params['m_val'], params['Iz_val'], params['Cf_val'], params['Cr_val'], a, params['a_val'] + params['b_val'] - a, # 保持轴距不变 step_steering ) ) plt.plot(t, states[:, 1], label=f'a={a:.1f} m') plt.xlabel('时间 [s]') plt.ylabel('横摆角速度 [rad/s]') plt.title('不同质心位置下的响应') plt.legend() plt.grid(True) plt.show()

6. 进阶应用与扩展

6.1 正弦扫频测试

分析车辆在不同频率转向输入下的响应：

def sine_steering(t, freq=1.0, amplitude=0.1): return amplitude * np.sin(2*np.pi*freq*t) # 频率响应分析 frequencies = [0.1, 0.5, 1.0, 2.0] plt.figure(figsize=(12, 8)) for freq in frequencies: states = odeint( vehicle_model, state0, t, args=( params['u_val'], params['m_val'], params['Iz_val'], params['Cf_val'], params['Cr_val'], params['a_val'], params['b_val'], lambda t: sine_steering(t, freq) ) ) plt.plot(t, states[:, 1], label=f'{freq} Hz') plt.xlabel('时间 [s]') plt.ylabel('横摆角速度 [rad/s]') plt.title('不同频率转向输入下的响应') plt.legend() plt.grid(True) plt.show()

6.2 稳定性分析

计算特征值判断系统稳定性：

# 构造状态矩阵A A = np.array([ [-(params['Cf_val'] + params['Cr_val'])/(params['m_val']*params['u_val']), (params['b_val']*params['Cr_val'] - params['a_val']*params['Cf_val'])/(params['m_val']*params['u_val']) - params['u_val']], [(params['b_val']*params['Cr_val'] - params['a_val']*params['Cf_val'])/(params['Iz_val']*params['u_val']), -(params['a_val']**2*params['Cf_val'] + params['b_val']**2*params['Cr_val'])/(params['Iz_val']*params['u_val'])] ]) # 计算特征值 eigenvalues = np.linalg.eig(A)[0] print("系统特征值:", eigenvalues)