行驶动力学建模仿真及主动悬架控制器设计.docx

行驶动力学建模仿真及主动悬架控制器设计.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

可以看出，并没有得到理想的仿真结果，说明仿真过程有问题。

但能力有限，未找出问题所在。