行駛動(dòng)力學(xué)建模,仿真及主動(dòng)懸架控制器設(shè)計(jì)

1、1. 計(jì)算機(jī)仿真系統(tǒng)模型的建立2. LOG 控制器設(shè)計(jì)3. 計(jì)算實(shí)例4. MATLAB 仿真過(guò)程5. 半車模型建模及仿真5.1 隨機(jī)線性最優(yōu)控制5.2 預(yù)瞄控制5.3 結(jié)果比較目錄-.1 -.-.2 -.3 -.-.4 -.-.8 -.9 - 11 -.12 -以單輪車輛模型為例，介紹行駛動(dòng)力學(xué)計(jì)算機(jī)建模、仿真分析以及利用線性二次最優(yōu)控制 理論進(jìn)行主動(dòng)懸架 LQG控制器設(shè)計(jì)過(guò)程。1.計(jì)算機(jī)仿真系統(tǒng)模型的建立根據(jù)圖7所示的主動(dòng)懸架單輪車輛模型，運(yùn)用牛頓運(yùn)動(dòng)定律，建立系統(tǒng)的運(yùn)動(dòng)方程，即:(4)Xg)mbXb Ua Ks(Xb Xw)mwXwU a Ks(xb Xw) Kt (xw這里，采用一個(gè)濾波

2、白噪聲作為路面輸入模型，即:式中，xg為路面垂向位移（m）; GO為路面不平度系數(shù) w為數(shù)字期望為零的高斯白噪聲； fO為下截止頻率（xg(t)2 foXg(t)2 y/GUw©)（6）（m3/cycle）； u為車輛前進(jìn)速度（m/s）；Hz）。JC方簧質(zhì)量位移q控制力X屮非簧載質(zhì)量1立移«講非黃載質(zhì)童產(chǎn)度叫路面（立移結(jié)合式（4）、式（5） 出系統(tǒng)的空間狀態(tài)方程：圖7單輪車輛模型和式（6），將系統(tǒng)運(yùn)動(dòng)方程和路面輸入方程寫成矩陣形式，即得(7)AX BU FW式中，X （Xb & Xb Xw Xg）T，為系統(tǒng)狀態(tài)矢量；W=（w（t）,為高斯白噪聲輸入矩 陣；U=（ U

3、a（t）,為輸入控制矩陣；Ksmb0mw00mbKt Kmw0Kamw0mb12. LOG控制器設(shè)計(jì)mw002 JGoU車輛懸架設(shè)計(jì)中的主要指標(biāo)包括： 的車身垂向振動(dòng)加速度；影響車身姿態(tài)且與輪胎布置有關(guān)的懸架動(dòng)行程。因此， 制器設(shè)計(jì)中的性能指標(biāo) J即為輪胎動(dòng)位移、懸架動(dòng)行程和車身垂向振動(dòng)加速度的加權(quán)平方和 在時(shí)域T內(nèi)的積分值，其表達(dá)式為：1 TJ Tim - q1Xw(t) Xg(t)2 區(qū) Xw(t)2 q3X&2(t)dtT 0代表輪胎接地性的輪胎動(dòng)載荷;代表輪胎舒適性LQG控(8) 式中，q1、q2和q3分別為輪胎動(dòng)位移、懸架動(dòng)行程和車身垂向振動(dòng)加速度的加權(quán)系數(shù)。加 權(quán)系數(shù)的選取

4、決定了設(shè)計(jì)者對(duì)懸架性能的傾向，如對(duì)車身垂向振動(dòng)加速度項(xiàng)選擇較大的權(quán) 值，則考慮更多的是提高車輛操縱穩(wěn)定性。為方便起見，值 q3=1。將性能指標(biāo)J的表達(dá)式(8)改寫成矩陣形式，即：T(XtQX UtRU 2XTNU)dt這里取車身垂向振動(dòng)加速度的加權(quán)式中,JHrn(9)q2q202mbKf2mbq202mbqiq22mbqiqiKaKs0(10)最優(yōu)控制反饋增益矩陣可有黎卡提(Riccati)方程當(dāng)車輛參數(shù)值和加權(quán)系數(shù)值確定后, 求出，其形式如下：PA AtP (PB N)R1(Bt P Nt) Q 0根據(jù)任意時(shí)刻的最優(yōu)反饋控制增益矩陣 K BTP NT，由車輛參數(shù)和加權(quán)系數(shù)決定。反饋狀態(tài)變量X

5、 （ t ）,就可得到t時(shí)刻作動(dòng)器的最優(yōu)控制力Ua ,即：(11)Ua(t) KX(t)3.計(jì)算實(shí)例這里，以某轎車的后懸架為例， 給出一個(gè)完整的計(jì)算實(shí)例，包括車輛模型參數(shù)、仿真路 面輸入?yún)?shù)、控制器的設(shè)計(jì)參數(shù)以及計(jì)算結(jié)果。此例中車輛以20m/s的速度在某典型路面上行駛，仿真時(shí)間T=50s。計(jì)算中輸入的各參數(shù)及數(shù)值詳見表 表2單輪車輛模型仿真輸入?yún)?shù)值車輛模型參數(shù)符號(hào)單位數(shù)值黃載質(zhì)量mbKg320非簧載質(zhì)量mwKg40懸架剛度KsN/m20000輪胎剛度KtN/m200000懸架工作空間SWScmm100仿真路面輸入?yún)?shù)符號(hào)單位數(shù)值路面不平度系數(shù)G0m3/cycle5.0x10-6車速Um/s2

6、0下截止頻率f0Hz0.1性能指標(biāo)加權(quán)參數(shù)符號(hào)單位數(shù)值輪胎動(dòng)位移q1一80000懸架動(dòng)行程q2一5車身加速度q3一1仿真計(jì)算中以式（6）所示的濾波白噪聲作為路面輸入模型。白噪聲的生成可直接調(diào)用MATLAB 函數(shù) WGN（ M , N , P）（此函數(shù)需要安裝信號(hào)處理工具箱Communications toolbox），20dB。1000m，仿真時(shí)間為50s。9），利用已知的矩陣 A、K，S，E=LQR（ A，B，K為最優(yōu)控制反饋增益其中M為生成矩陣的行數(shù)，N為列數(shù)，P為白噪聲的功率（單位為dB）。本例中取M=10001， N=1，P=20。這意味著仿真計(jì)算中去一條白噪聲，共10001個(gè)采集點(diǎn)，

7、噪聲強(qiáng)度為設(shè)定采樣時(shí)間為0.005s、車速為20m/s時(shí)，相當(dāng)于仿真路面長(zhǎng)度為根據(jù)建立的系統(tǒng)狀態(tài)方程式（7）及最優(yōu)化性能指標(biāo)函數(shù)式（ B、Q、R、N，調(diào)用MATLAB中的線性二次最優(yōu)控制器設(shè)計(jì)函數(shù) Q, R，N），即可完成最優(yōu)主動(dòng)懸架控制器的設(shè)計(jì)。輸出的結(jié)果中， 矩陣，S為黎卡提方程的解，E為系統(tǒng)閉環(huán)特征根。K為:根據(jù)表2給出的仿真輸入?yún)?shù)，本例中求得的最優(yōu)反饋增益矩陣K=（ 711.88-1241.5-19284-2038.520864）同時(shí)，還得到了黎卡提方程的解：2.45590.02892.47458.66077.3090.02890.48860.02987.52627.23642.47

8、450.02984.97448.67545.10338.66077.52628.67542710.12700.47.3097.23645.10332700.42693.7在相同的仿真條件下，可將所設(shè)計(jì)的主動(dòng)懸架系統(tǒng)與一個(gè)被動(dòng)系統(tǒng)進(jìn)行對(duì)比分析。在被動(dòng)懸架系統(tǒng)中，取懸架剛度Ks=22000N/m，阻尼系數(shù) Cs=1000NS/m。除此之外，其他輸入?yún)?shù)值均與主動(dòng)懸架系統(tǒng)完全相同。4. MATLAB 仿真過(guò)程1）生成路面輸入模型 代碼如下：road_filea=wgn(10001,1,20); t=0:0.005:50; road_file(:,1)=t' road_file(:,2)=a;

9、 save road_file2) 參數(shù)輸入%載入路面數(shù)據(jù)模型Ks=22000;mb=320;Kt=200000;mw=40;f0=0.1;G0=0.000005;u=20;%輸入仿真有關(guān)參數(shù)%建立主動(dòng)懸架的狀態(tài)矩陣load road_file.mat代碼如下：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;%建立被動(dòng)懸架的狀態(tài)矩陣Cs/mw,-Cs/m

10、w,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創(chuàng)建仿真框圖狀態(tài)變量 X Xb,Xw,

11、Xb,Xw,Xg輸入與系統(tǒng)模塊，如下圖:Gain輸出模塊，如下圖:踣面輸入車身加速度懸架動(dòng)行程弘）XbXw輪胎動(dòng)位移XwXg路面輸入Xg整體程序框圖如下:Gainy- Cif+DuSUie-SpaceI Toadjiile IFromW0rt(5P3C*盲出媲痙量Garni路面輸人- Ai+Buk y = Cx4Du>j dLitit|Slate-Space1X14 ）結(jié)果分析可以直接通過(guò)雙擊 scope查看輸出的波形圖，為更好比較主動(dòng)懸架與被動(dòng)懸架的差別,通過(guò)輸出到works pace的狀態(tài)變量編程繪圖并計(jì)算均方根值。代碼如下: %繪制車身加速度曲線，并計(jì)算均方根值 %ba-主動(dòng)懸架車

12、身加速度%ba1-被動(dòng)懸架車身加速度ba=dif(X.data(:,1)./diff(X.time);ba仁diff(X1.data(:,1)./dif(X1.time);sub plot(2,1,1)plot(X.time(1:e nd-1),ba)sub plot(2,1,2)p lot(X1.time(1:e nd-1),ba1)BA= norm(ba,2)./(le ngth(ba)A0.5);BA仁 norm(ba1,2)./(le ngth(ba1)A0.5);%繪制懸架動(dòng)行程曲線，并計(jì)算其均方根值%sws-主動(dòng)懸架動(dòng)行程%sws1-被動(dòng)懸架動(dòng)行程figure。sws=X.data

13、(:,3)-X.data(:,4);sws1=X1.data(:,3)-X1.data(:,4);sub plot(2,1,1)plot(X.time,sws)sub plot(2,1,2)plot(X.time,sws1)SWS=norm(1000*sws,2)./(le ngth(sws).A0.5);SWS仁 norm(1000*sws1,2)./(le ngth(sws1).A0.5);%繪制輪胎動(dòng)位移曲線，并計(jì)算其均方根值%dtd-主動(dòng)懸架動(dòng)位移%dtd1-被動(dòng)懸架動(dòng)位移figure。dtd=X.data(:,4)-X.data(:,5);dtd1=X1.data(:,4)-X1.d

14、ata(:,5);sub plot(2,1,1)plot(X.time,dtd)sub plot(2,1,2)plot(X.time,dtd1)DTD=norm(1000*dtd,2)./(le ngth(dtd).A0.5);DTD1= norm(1000*dtd1,2)./(le ngth(dtd1).A0.5);結(jié)果如下:車身加速度曲線10-10020XZS/卩度速加身車r|H IJ JN,' "! i'lll-ii'；卜i1 LF / cl Fl!./ :耐【J?J I h 11 J F 川！1 TI'b'j' lyi. 1丫冃

15、H 1 I' f1 f M主動(dòng)懸架車身加速度50-5123789456時(shí)間被動(dòng)懸架車身加速度10-200 II r.川門ih.;川t八，丹'1-JIV, 11.艸：; ”->a =11L *' |100-10123789456時(shí)間10懸架動(dòng)行程曲線0.1卩移位0.6卩移位hiX嚴(yán)1 .11 11I-II 1I-'r-H1 H11.1.11- rJJ-叮廠1 -I'JJ=h 1 11主動(dòng)懸架動(dòng)行程0.050-0.05-0.1123467895時(shí)間 被動(dòng)懸架動(dòng)行程10-0.20Wrrf-.1 e-.J.M ；"：1 TIL匸：'卜 ，

16、-(.ifr- N. J' lll . I叫和L'-V川T kql11 h41 ,0.40.20123467895時(shí)間10輪胎動(dòng)位移0.04卩移位，占 l；p1lF 4.'« .j!：v h/1:摘IiLJ. Il-. 1 - 1 小訐沁J1川冷Ji :匸川£.悩：| FL；.if fl，巾111主動(dòng)懸架輪胎動(dòng)位移0.020-0.02-0.04123467895時(shí)間被動(dòng)懸架輪胎動(dòng)位移100.1-0.10.050-0.051236789105時(shí)間04主動(dòng)懸架與被動(dòng)懸架性能指標(biāo)均方根值比較性能指標(biāo)主動(dòng)懸架被動(dòng)懸架車身加速度BA(m/s2)1.514.60

17、懸架動(dòng)行程SWS(mm)34.4267.10輪胎動(dòng)位移DTD(mm)5.8734.575.半車模型建模及仿真半車模型建模懸架的半車模型有四個(gè)自由度，可選取前后輪垂直位移、 前后懸架與車身連接處垂直位移四個(gè)自由度，寫出狀態(tài)空間方程： 以X Z4,Z3,Z2,Zi,Z4,Z3,Z2,Zi,Xi0, X20作為系統(tǒng)狀態(tài)變量，狀態(tài)空間方程如下:X(t) AX(t) BU(t) FW(t)后懸架作動(dòng)器力;u Uaf ,為控制輸入矩陣，即前 UarWi,為路面輸入模型的高斯W2白噪聲輸入。LQG控制器選取如下二次型性能指標(biāo):1 T 2J Tim - q1(Z1 Z01)qzZ zjT 0權(quán)系數(shù);數(shù)；? 2

18、1 Z2q3( Z3202 )? 22Z4 dtq1 ,q3,前后懸架輪胎動(dòng)位移加q2,q4,前后懸架動(dòng)行程加權(quán)系2,車身加速度加權(quán)系數(shù)。寫成矩陣形式:1 TJ lim (XtQX UtRUT T 02XtNU )dt由黎卡提方程求出反饋矩陣K，則U KX，因此狀態(tài)空間方程可寫成:X(t) (A BK)X(t) FW(t) (A BK)-狀態(tài)矩陣F輸入矩陣由此形式可方便求出其時(shí)域與頻域響應(yīng)。5.1隨機(jī)線性最優(yōu)控制路面模型本例仿真車速為 20m/s,軸距為2.8m，因此滯后時(shí)間為 0.14s。若取仿真時(shí)間 T為20s, 采樣時(shí)間間隔為0.005s，因而仿真點(diǎn)數(shù)4028。程序代碼如下： %生成路面

19、模型a=wg n(4029,1,20);t=0:0.005:20.14; r=a,t'save road_file r% 車身模型參數(shù)輸入 clear clc mb=690;I=1222;mwf=40;mwr=45;Ksf=17000;Ksr=22000;Ktf=200000;Ktr=200000;a=1.3;b=1.5;% 仿真路面參數(shù)輸入G0=5e-6;u=20;f0=0.1;% 性能指標(biāo)加權(quán)系數(shù)q1=80000;q2=100;q3=80000;q4=100;%半車懸架模型建模，計(jì)算A、B、F矩陣a1=1/mb+bA2/l;a2=1/mb-a*b/I;a3=1/mb+a2/l;A=

20、0 0 0 0 a1*Ksr -a1*Ksr a2*Ksf -a2*Ksf 0 0;0 0 0 0 -Ksr/mwr (Ksr-Ktr)/mwr 0 0 Ktr/mwr 0;0 0 0 0 a2*Ksr -a2*Ksr a3*Ksf -a3*Ksf 0 0;0 00 0 0 0-Ksf/mwf (Ksf-Ktf)/mwf 0 Ktf/mwf;1 00 0 0 000 0 0;0 10 0 0 000 0 0;0 01 0 0 000 0 0;0 0 0 1 0 0 0 0 0 0;0 0 0 0 0 0 0 0 -2*pi*f0 0;0 0 0 0 0 0 0 0 0 -2*pi*f0;B=a

21、2 a1;0 -1/mwr;a3 a2;-1/mwf 0;zeros(6,2);F=zeros(8,2);2*pi*sqrt(G0*u) 0;0 2*pi*sqrt(G0*u);% LQG控制器設(shè)計(jì)，計(jì)算 Q、R、N矩陣并求出反饋矩陣 K b1=a32+a22;b2=a2*a3+a1*a2;b3=a22+a12;sq=q4+KsrA2*b3 -q4-Ksr2*b3 Ksf*Ksr*b2 -Ksf*Ksr*b2 0 0;-q4-KsrA2*b3 q3+q4+Ksr2*b3 -Ksf*Ksr*b2 Ksf*Ksr*b2 -q3 0;Ksf*Ksr*b2 -Ksf*Ksr*b2 q2+Ksf2*b1

22、 -q2-Ksf2*b1 0 0;-Ksf*Ksr*b2 Ksf*Ksr*b2 -q2-Ksf2*b1 q1+q2+Ksf2*b1 0 -q1; 0 -q3 0 0 q3 0;0 0 0 -q1 0 q1;Q=zeros(4,10);zeros(6,4) sq;R=b1 b2;b2 b3;N=zeros(4,2);Ksr*b2-Ksf*b2;zeros(2,2);K=lqr(A,B,Q,R,N);%求主動(dòng)懸架的時(shí)域響應(yīng)load road_file.mat w1=road_file(:,2);w2=w1;t=road_file(:,1);C=eye(10);D=zeros(10,2);Y,X=l

23、sim(A-B*K,F,C,D,w1'w2',t);Ksr*b3;-Ksr*b2-Ksr*b3;Ksf*b1Ksf*b2;-Ksf*b1輸出矩陣C為10X10單位矩陣，即將所有狀態(tài)均輸出，因此系統(tǒng)輸出 X是完全相同的。Y與系統(tǒng)狀態(tài)變量5.2預(yù)瞄控制預(yù)瞄控制中，前后輪的路面輸入不再相同，后輪的輸入比前輪輸入滯后一個(gè)預(yù)瞄時(shí)間,輪路面輸入w1與后輪路面輸入w2滿足：及前w2sew,采用PaDa近似后，可轉(zhuǎn)化為時(shí)域方程:則預(yù)瞄控制系統(tǒng)空間狀態(tài)方程可寫為:X A FD X B U 0FEw,同樣可由黎卡提方程求出其反饋矩陣Kp程序代碼如下: %預(yù)瞄控制 tao=0.14;aO=12/ta

24、oA2; a1=6/tao; a2=1;An=0 1;-a0 -a1;Bn=-2*a1;6*a0;Dn=0 0;1 0;En=1;1;A0=A F*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 結(jié)果比較運(yùn)行結(jié)果中的變量 X 、X0

25、分別存儲(chǔ)無(wú)預(yù)瞄和有預(yù)瞄的狀態(tài)變量，預(yù)瞄控制主要改善后懸 架的性能，因此選擇后懸架的加速度、后懸架動(dòng)行程及后輪動(dòng)位移作為比較項(xiàng)目。程序代碼如下：% 繪制車身加速度曲線，并計(jì)算均方根值 %ba- 最優(yōu)控制車身加速度 %ba0- 預(yù)瞄控制車身加速度 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)./(le ngth(ba).A0.5);BAO=norm(ba0,2)./(le ngth(ba0).A0.5);% 繪制懸架動(dòng)行程曲線，并計(jì)算其均方根值 %swsr- 最優(yōu)控制懸架動(dòng)行程 %swsr0- 預(yù)瞄控制后懸架動(dòng)行程 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)./(le ngth(swsr).A0.5);SWS