汽车半主动悬架最优控制仿真

汽车半主动悬架最优控制仿真
- 简介
- 使用说明

简介

半主动悬架介于被动悬架和全主动悬架之间，它通过改变减振器的阻尼特性适应不同的道路和行驶状况的需要，改善乘坐舒适性和操纵稳定性。由于半主动悬架在控制品质上接近全主动悬架，且结构简单，不需要力源，能量损耗小，因此被广泛使用。

使用说明

一、实验目的

1.建立汽车半主动悬架数学模型

2.建立汽车半主动悬架最优控制数学模型

3.求汽车半主动悬架最优控制参数

4.绘制汽车半主动悬架输出变量的时域特性曲线

二、仿真数据

汽车半主动悬架最优控制仿真所需参数见表6-11-1。

表6-11-1 汽车半主动悬架最优控制仿真所需参数

悬挂质量/kg	非悬挂质量/kg	悬架刚度/(N/m)	悬架不变阻尼系数/(N·s/m)
280	50	19000	1400
轮胎刚度/(N/m)	下截止频率/Hz	路面不平度系数	仿真时间/s
180000	0.07	5×10⁻⁶	10

三、实验步骤

1.建立汽车半主动悬架数学模型

原理参考教材第六章实例11

2.建立汽车半主动悬架最优控制数学模型

原理参考教材第六章实例11

3.求汽车半主动悬架最优控制参数

根据汽车半主动悬架最优控制数学模型，编写求汽车半主动悬架最优控制参数的MWORKS程序如下。

ms = 280 # 车身质量
mw = 50 # 轮胎和车轮组质量
ks = 19000 # 车身弹簧刚度
Cs = 1400 # 车身阻尼系数
kw = 180000 # 轮胎和车轮组刚度
q1 = 30000000
q2 = 4000000
q3 = 200
q4 = 1
A = [0 1 0 -1
-ks/ms -Cs/ms 0 Cs/ms
0 0 0 -1
ks/mw Cs/mw kw/mw -Cs/mw] # 状态矩阵
B = [0 1.0 / ms 0 -1 / mw]' # 输入矩阵
E = [0 0 1 0]' # 扰动输入矩阵
Qd = [q1 * ks^2 / (ms^2) + q2 q1 * ks * Cs / (ms^2) 0 -q1 * ks * Cs / (ms^2)
q1 * ks * Cs / (ms^2) q1 * Cs^2 / (ms^2) 0 -q1 * Cs^2 / (ms^2)
0 0 q3 * ks^2 0
-q1 * ks * Cs / (ms^2) -q1 * Cs^2 / (ms^2) -0 q1 * Cs^2 / ms^2] # 状态权重矩阵
Rd = q4 + q1 / ms^2 # 控制输入权重
Nd = 1 / ms^2 * [-q1 * ks -q1 * Cs -ks q1 * Cs]' # 扰动输入权重
K, S, E = lqr(A, B, Qd, Rd, Nd) # 计算LQR增益矩阵
println("K = ", K)

可得汽车半主动悬架最优控制参数为

k1= - 17548.95,k2= - 366.26,k3 = - 91.11,k4 = 1312.93

4.绘制汽车半主动悬架输出变量的时域特性曲线

根据汽车半主动悬架特性数学模型，编写绘制路面位移、车身垂直加速度、悬架动挠度和轮胎动载荷时域特性曲线的MWORKS程序如下。

ms = 280 # 车身质量
mw = 50 # 轮胎和车轮组质量
Ks = 19000 # 车身弹簧刚度
Cs = 1400 # 车身阻尼系数
Kw = 180000 # 轮胎和车轮组刚度
k1 = -18024.61
k2 = -675.99
k3 = 196.44
k4 = 1160.70
u = 16.67 # 系统输入
f0 = 0.07 # 系统频率
Sq = 0.000005 # 噪声方差
A = [0 1 0 -1
-Ks/ms-k1/ms -Cs/ms-k2/ms -k3/ms Cs/ms-k4/ms
0 0 0 -1
Ks/mw+k1/mw Cs/mw+k2/mw Kw/mw+k3/mw -Cs/mw+k4/mw] # 状态矩阵
B = [0 0 1 0]' # 输入矩阵
C = [-Ks/ms-k1/ms -Cs/ms-k2/ms -k3/ms Cs/ms-k4/ms
1 0 0 0
0 0 Kw 0] # 输出矩阵
white_noise = 1 .- 2 * 1 * rand(1, 1000) # 白噪声输入
q = Array{Float64}(undef, 1001)
dq = Array{Float64}(undef, 1000)
for i = 1:1000
    q[1] = 0
    q[i] = 0.01 * (-2 * pi * f0 * q[i] + 2 * pi * (sqrt(Sq * u)) * white_noise[i]) + q[i]
    global dq = 100 * diff(q)
end

ms = 280 # 车身质量
mw = 50 # 轮胎和车轮组质量
Ks = 19000 # 车身弹簧刚度
Cs = 1400 # 车身阻尼系数
Kw = 180000 # 轮胎和车轮组刚度
k1 = -18024.61
k2 = -675.99
k3 = 196.44
k4 = 1160.70
u = 16.67 # 系统输入
f0 = 0.07 # 系统频率
Sq = 0.000005 # 噪声方差
A = [0 1 0 -1
-Ks/ms-k1/ms -Cs/ms-k2/ms -k3/ms Cs/ms-k4/ms
0 0 0 -1
Ks/mw+k1/mw Cs/mw+k2/mw Kw/mw+k3/mw -Cs/mw+k4/mw] # 状态矩阵
B = [0 0 1 0]' # 输入矩阵
C = [-Ks/ms-k1/ms -Cs/ms-k2/ms -k3/ms Cs/ms-k4/ms
1 0 0 0
0 0 Kw 0] # 输出矩阵
white_noise = 1 .- 2 * 1 * rand(1, 1000) # 白噪声输入
q = Array{Float64}(undef, 1001)
dq = Array{Float64}(undef, 1000)
for i = 1:1000
    q[1] = 0
    q[i] = 0.01 * (-2 * pi * f0 * q[i] + 2 * pi * (sqrt(Sq * u)) * white_noise[i]) + q[i]
    global dq = 100 * diff(q)
end

x = Array{Float64}(undef, 4, 1001) # 存储系统状态的数组
y = Array{Float64}(undef, 3, 1000) # 存储系统输出的数组
for i = 1:1000
    x[:, 1] = [0; 0; 0; 0] # 初始化状态向量为零
    sysc = ss(A, B, zeros(size(A)), zeros(size(B))) # 创建连续时间系统
    sysd = c2d(sysc, 0.01) # 将连续时间系统离散化为0.01秒的采样时间
    G = sysd.A # 离散化后的状态转移矩阵
    H = sysd.B # 离散化后的输入矩阵
    x[:, i+1] = G * x[:, i] + H * dq[i] # 计算下一个时刻的状态向量
    y[:, i] = C * x[:, i] # 计算当前时刻的输出向量
end

figure(1)
t = 0:0.01:10 # 时间向量
plot(t, q) # 绘制q随时间的变化图
xlabel("时间/s")
ylabel("路面位移/m")
figure(2)
t1 = 0.01:0.01:10 # 时间向量
plot(t1, y[1, :]) # 绘制输出y1随时间的变化图
xlabel("时间/s")
ylabel("车身垂直加速度/(m/s²)")
figure(3)
t1 = 0.01:0.01:10 # 时间向量
plot(t1, y[2, :]) # 绘制输出y2随时间的变化图
xlabel("时间/s")
ylabel("悬架动挠度/m")
figure(4)
t1 = 0.01:0.01:10 # 时间向量
plot(t1, y[3, :]) # 绘制输出y3随时间的变化图
xlabel("时间/s")
ylabel("轮胎动载荷/N")

在MWORKS编辑器中输入这些程序，点击运行按钮，就会得到路面位移时域特性曲线(图6-11-2)、车身垂直加速度时域特性曲线(图6-11-3)、悬架动挠度时域特性曲线(图6-11-4)、轮胎动载荷时域特性曲线(图6-11-5)。

图6-11-2 轮胎动载荷时域特性曲线

图6-11-3 车身垂直加速度时域特性曲线

图6-11-4 悬架动挠度时域特性曲线

图6-11-5 轮胎动载荷时域特性曲线