行驶动力学建模仿真及主动悬架控制器设计.docx
《行驶动力学建模仿真及主动悬架控制器设计.docx》由会员分享，可在线阅读，更多相关《行驶动力学建模仿真及主动悬架控制器设计.docx（22页珍藏版）》请在冰豆网上搜索。
行驶动力学建模仿真及主动悬架控制器设计
以单轮车辆模型为例，介绍行驶动力学计算机建模、仿真分析以及利用线性二次最优控制理论进行主动悬架LQG控制器设计过程。
1.计算机仿真系统模型的建立
根据图7所示的主动悬架单轮车辆模型，运用牛顿运动定律，建立系统的运动方程，即：
（4）
（5）
这里，采用一个滤波白噪声作为路面输入模型，即：
（6）
式中，xg为路面垂向位移（m）；G0为路面不平度系数（m3/cycle）；u为车辆前进速度（m/s）；w为数字期望为零的高斯白噪声；f0为下截止频率（Hz）。
图7单轮车辆模型
结合式（4）、式（5）和式（6），将系统运动方程和路面输入方程写成矩阵形式，即得出系统的空间状态方程：
（7）
式中，
，为系统状态矢量；W=（w（t）），为高斯白噪声输入矩阵；U=（Ua（t）），为输入控制矩阵；
；
；
2.LOG控制器设计
车辆悬架设计中的主要指标包括：
①代表轮胎接地性的轮胎动载荷；②代表轮胎舒适性的车身垂向振动加速度；③影响车身姿态且与轮胎布置有关的悬架动行程。
因此，LQG控制器设计中的性能指标J即为轮胎动位移、悬架动行程和车身垂向振动加速度的加权平方和在时域T内的积分值，其表达式为：
（8）
式中，q1、q2和q3分别为轮胎动位移、悬架动行程和车身垂向振动加速度的加权系数。
加权系数的选取决定了设计者对悬架性能的倾向，如对车身垂向振动加速度项选择较大的权值，则考虑更多的是提高车辆操纵稳定性。
为方便起见，这里取车身垂向振动加速度的加权值q3=1。
将性能指标J的表达式（8）改写成矩阵形式，即：
（9）
式中，
；
；
当车辆参数值和加权系数值确定后，最优控制反馈增益矩阵可有黎卡提（Riccati）方程求出，其形式如下：
（10）
最优反馈控制增益矩阵
，由车辆参数和加权系数决定。
根据任意时刻的反馈状态变量X（t），就可得到t时刻作动器的最优控制力Ua，即：
（11）
3.计算实例
这里，以某轿车的后悬架为例，给出一个完整的计算实例，包括车辆模型参数、仿真路面输入参数、控制器的设计参数以及计算结果。
此例中车辆以20m/s的速度在某典型路面上行驶，仿真时间T=50s。
计算中输入的各参数及数值详见表2。
表2单轮车辆模型仿真输入参数值
车辆模型参数
符号
单位
数值
簧载质量
非簧载质量
悬架刚度
轮胎刚度
悬架工作空间
mb
mw
Ks
Kt
SWSc
Kg
Kg
N/m
N/m
mm
320
40
20000
200000
仿真路面输入参数
符号
单位
数值
路面不平度系数
车速
下截止频率
G0
U
f0
m3/cycle
m/s
Hz
5.0x10-6
20
0.1
性能指标加权参数
符号
单位
数值
轮胎动位移
悬架动行程
车身加速度
q1
q2
q3
—
—
—
80000
5
1
仿真计算中以式（6）所示的滤波白噪声作为路面输入模型。
白噪声的生成可直接调用MATLAB函数WGN（M，N，P）（此函数需要安装信号处理工具箱Communicationstoolbox），其中M为生成矩阵的行数，N为列数，P为白噪声的功率（单位为dB）。
本例中取M=10001，N=1，P=20。
这意味着仿真计算中去一条白噪声，共10001个采集点，噪声强度为20dB。
设定采样时间为0.005s、车速为20m/s时，相当于仿真路面长度为1000m，仿真时间为50s。
根据建立的系统状态方程式（7）及最优化性能指标函数式（9），利用已知的矩阵A、B、Q、R、N，调用MATLAB中的线性二次最优控制器设计函数[K，S，E]=LQR（A，B，Q，R，N），即可完成最优主动悬架控制器的设计。
输出的结果中，K为最优控制反馈增益矩阵，S为黎卡提方程的解，E为系统闭环特征根。
根据表2给出的仿真输入参数，本例中求得的最优反馈增益矩阵K为：
K=（711.88-1241.5-19284-2038.520864）
同时，还得到了黎卡提方程的解：
在相同的仿真条件下，可将所设计的主动悬架系统与一个被动系统进行对比分析。
在被动悬架系统中，取悬架刚度Ks=22000N/m，阻尼系数Cs=1000NS/m。
除此之外，其他输入参数值均与主动悬架系统完全相同。
4.MATLAB仿真过程
1）生成路面输入模型
代码如下：
a=wgn（10001,1,20）;
t=0:
0.005:
50;
road_file（:
1）=t';
road_file（:
2）=a;
saveroad_fileroad_file
2）参数输入
代码如下：
loadroad_file.mat%载入路面数据模型
Ks=22000;mb=320;Kt=200000;mw=40;f0=0.1;G0=0.000005;u=20;Kb=20000;Ks1=22000;Cs=1000;%输入仿真有关参数
A=[0,0,-Ks/mb,Ks/mb,0;%建立主动悬架的状态矩阵
0,0,Ks/mw,（-Kt-Ks）/mw,Kt/mw;
1,0,0,0,0;
0,1,0,0,0;
0,0,0,0,-2*pi*f0];
A1=[-Cs/mb,Cs/mb,-Ks1/mb,Ks1/mb,0;%建立被动悬架的状态矩阵
Cs/mw,-Cs/mw,Ks1/mw,（-Kt-Ks1）/mw,Kt/mw;
1,0,0,0,0;
0,1,0,0,0;
0,0,0,0,-2*pi*f0];
B=[1/mb,0;
-1/mw,0;
0,0;
0,0;
0,2*pi*sqrt（G0*u）];
B1=[0,0;
0,0;
0,0;
0,0;
0,2*pi*sqrt（G0*u）];
C=[1,0,0,0,0;
0,1,0,0,0;
0,0,1,0,0;
0,0,0,1,0;
0,0,0,0,1];
D=[0,0;
0,0;
0,0;
0,0;
0,0];
K=[711.88,-1241.5,-19284,-2038.5,20864];
K1=[0,0,0,0,0];
3）用Simulink创建仿真框图
●状态变量
●输入与系统模块，如下图：
●输出模块，如下图：
●整体程序框图如下：
4）结果分析
可以直接通过双击scope查看输出的波形图，为更好比较主动悬架与被动悬架的差别，下面通过输出到workspace的状态变量编程绘图并计算均方根值。
代码如下：
%%绘制车身加速度曲线，并计算均方根值
%ba-主动悬架车身加速度
%ba1-被动悬架车身加速度
ba=diff（X.data（:
1））./diff（X.time）;
ba1=diff（X1.data（:
1））./diff（X1.time）;
subplot（2,1,1）
plot（X.time（1:
end-1）,ba）
subplot（2,1,2）
plot（X1.time（1:
end-1）,ba1）
BA=norm（ba,2）./（length（ba）.^0.5）;
BA1=norm（ba1,2）./（length（ba1）.^0.5）;
%%绘制悬架动行程曲线，并计算其均方根值
%sws-主动悬架动行程
%sws1-被动悬架动行程
figure（）
sws=X.data（:
3）-X.data（:
4）;
sws1=X1.data（:
3）-X1.data（:
4）;
subplot（2,1,1）
plot（X.time,sws）
subplot（2,1,2）
plot（X.time,sws1）
SWS=norm（1000*sws,2）./（length（sws）.^0.5）;
SWS1=norm（1000*sws1,2）./（length（sws1）.^0.5）;
%%绘制轮胎动位移曲线，并计算其均方根值
%dtd-主动悬架动位移
%dtd1-被动悬架动位移
figure（）
dtd=X.data（:
4）-X.data（:
5）;
dtd1=X1.data（:
4）-X1.data（:
5）;
subplot（2,1,1）
plot（X.time,dtd）
subplot（2,1,2）
plot（X.time,dtd1）
DTD=norm（1000*dtd,2）./（length（dtd）.^0.5）;
DTD1=norm（1000*dtd1,2）./（length（dtd1）.^0.5）;
结果如下：
1 车身加速度曲线
2 悬架动行程曲线
3 轮胎动位移
●主动悬架与被动悬架性能指标均方根值比较
性能指标
主动悬架
被动悬架
车身加速度BA（m/s2）
1.51
4.60
悬架动行程SWS（mm）
34.42
67.10
轮胎动位移DTD（mm）
5.87
34.57
5.半车模型建模及仿真
●半车模型建模
悬架的半车模型有四个自由度，可选取前后轮垂直位移、前后悬架与车身连接处垂直位移四个自由度，写出状态空间方程：
以
作为系统状态变量，状态空间方程如下：
●LQG控制器
选取如下二次型性能指标：
写成矩阵形式：
由黎卡提方程求出反馈矩阵K，则
，因此状态空间方程可写成：
由此形式可方便求出其时域与频域响应。
5.1随机线性最优控制
●路面模型
本例仿真车速为20m/s，轴距为2.8m，因此滞后时间为0.14s。
若取仿真时间T为20s，采样时间间隔为0.005s，因而仿真点数4028。
程序代码如下：
%%生成路面模型
a=wgn（4029,1,20）;
t=0:
0.005:
20.14;
r=[a,t'];
saveroad_filer
%%车身模型参数输入
clear
clc
mb=690;
I=1222;
mwf=40;
mwr=45;
Ksf=17000;
Ksr=22000;
Ktf=200000;
Ktr=200000;
a=1.3;
b=1.5;
%%仿真路面参数输入
G0=5e-6;
u=20;
f0=0.1;
%%性能指标加权系数
q1=80000;
q2=100;
q3=80000;
q4=100;
%%半车悬架模型建模，计算A、B、F矩阵
a1=1/mb+b^2/I;
a2=1/mb-a*b/I;
a3=1/mb+a^2/I;
A=[0000a1*Ksr-a1*Ksra2*Ksf-a2*Ksf00;
0000-Ksr/mwr（Ksr-Ktr）/mwr00Ktr/mwr0;
0000a2*Ksr-a2*Ksra3*Ksf-a3*Ksf00;
000000-Ksf/mwf（Ksf-Ktf）/mwf0Ktf/mwf;
1000000000;
0100000000;
0010000000;
0001000000;
00000000-2*pi*f00;
000000000-2*pi*f0];
B=[a2a1;0-1/mwr;a3a2;-1/mwf0;zeros（6,2）];
F=[zeros（8,2）;2*pi*sqrt（G0*u）0;02*pi*sqrt（G0*u）];
%%LQG控制器设计，计算Q、R、N矩阵并求出反馈矩阵K
b1=a3^2+a2^2;
b2=a2*a3+a1*a2;
b3=a2^2+a1^2;
sq=[q4+Ksr^2*b3-q4-Ksr^2*b3Ksf*Ksr*b2-Ksf*Ksr*b200;
-q4-Ksr^2*b3q3+q4+Ksr^2*b3-Ksf*Ksr*b2Ksf*Ksr*b2-q30;
Ksf*Ksr*b2-Ksf*Ksr*b2q2+Ksf^2*b1-q2-Ksf^2*b100;
-Ksf*Ksr*b2Ksf*Ksr*b2-q2-Ksf^2*b1q1+q2+Ksf^2*b10-q1;
0-q300q30;
000-q10q1];
Q=[zeros（4,10）;zeros（6,4）sq];
R=[b1b2;b2b3];
N=[zeros（4,2）;Ksr*b2Ksr*b3;-Ksr*b2-Ksr*b3;Ksf*b1Ksf*b2;-Ksf*b1-Ksf*b2;zeros（2,2）];
K=lqr（A,B,Q,R,N）;
%%求主动悬架的时域响应
loadroad_file.mat
w1=road_file（:
2）;
w2=w1;
t=road_file（:
1）;
C=eye（10）;
D=zeros（10,2）;
[Y,X]=lsim（A-B*K,F,C,D,[w1';w2'],t）;
输出矩阵C为10X10单位矩阵，即将所有状态均输出，因此系统输出Y与系统状态变量X是完全相同的。
5.2预瞄控制
预瞄控制中，前后轮的路面输入不再相同，后轮的输入比前轮输入滞后一个预瞄时间，及前轮路面输入w1与后轮路面输入w2满足：
采用PaDa近似后，可转化为时域方程：
则预瞄控制系统空间状态方程可写为：
同样可由黎卡提方程求出其反馈矩阵
程序代码如下：
%%预瞄控制
tao=0.14;
a0=12/tao^2;
a1=6/tao;
a2=1;
An=[01;-a0-a1];
Bn=[-2*a1;6*a0];
Dn=[00;10];
En=[1;1];
A0=[AF*Dn;zeros（2,10）An];
B0=[B;zeros（2,2）];
C0=[C,zeros（10,2）;zeros（2,12）];
%C0=eye（12）;
D0=zeros（12,1）;
F0=[F*En;Bn];
Q0=[Q,zeros（10,2）;zeros（2,12）];
N0=[N;zeros（2,2）];
K0=lqr（A0,B0,Q0,R,N0）;
[Y0,X0]=lsim（A0-B0*K0,F0,C0,D0,w1,t）;
5.3结果比较
运行结果中的变量X、X0分别存储无预瞄和有预瞄的状态变量，预瞄控制主要改善后悬架的性能，因此选择后悬架的加速度、后悬架动行程及后轮动位移作为比较项目。
程序代码如下：
%%绘制车身加速度曲线，并计算均方根值
%ba-最优控制车身加速度
%ba0-预瞄控制车身加速度
ba=diff（X（:
1））/0.005;
ba0=diff（X0（:
1））/0.005;
subplot（2,1,1）
plot（t（1:
end-1）,ba）
subplot（2,1,2）
plot（t（1:
end-1）,ba0）
BA=norm（ba,2）./（length（ba）.^0.5）;
BA0=norm（ba0,2）./（length（ba0）.^0.5）;
%%绘制悬架动行程曲线，并计算其均方根值
%swsr-最优控制悬架动行程
%swsr0-预瞄控制后悬架动行程
figure（）
swsr=X（:
5）-X（:
6）;
swsr0=X0（:
5）-X0（:
6）;
subplot（2,1,1）
plot（t,swsr）
subplot（2,1,2）
plot（t,swsr0）
SWS=norm（1000*swsr,2）./（length（swsr）.^0.5）;
SWS0=norm（1000*swsr0,2）./（length（swsr0）.^0.5）;
%%绘制轮胎动位移曲线，并计算其均方根值
%dtdr-最优控制动轮胎位移
%dtdr0-预瞄控制后悬架轮胎动位移
figure（）
dtdr=X（:
6）-X（:
10）;
dtdr0=X0（:
6）-X0（:
10）;
subplot（2,1,1）
plot（t,dtdr）
subplot（2,1,2）
plot（t,dtdr0）
DTD=norm（1000*dtdr,2）./（length（dtdr）.^0.5）;
DTD0=norm（1000*dtdr0,2）./（length（dtdr0）.^0.5）;
①、后悬架加速度曲线
②、后悬架动行程曲线
③、后轮动位移比较
④、性能指标均方根值比较
后悬架加速度/（m/s2）
后悬架动行程（mm）
后轮动位移（mm）
无预瞄
1.5254
16.8165
6.0459
有预瞄
1.5254
16.8164
6.0427
可以看出，并没有得到理想的仿真结果，说明仿真过程有问题。
但能力有限，未找出问题所在。
升级会员 