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

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

2、w)這里，采用一個(gè)濾波白噪聲作為路面輸入模型，即:xg(t)2 foXg(t) 2G0uw(t)（6）（m3/cycle） ； u為車輛前進(jìn)速度（m/s）；Hz）。式中，xg為路面垂向位移（m）; GO為路面不平度系數(shù) w為數(shù)字期望為零的高斯白噪聲； fO為下截止頻率（結(jié)合式（4）、式（5）和式（6），將系統(tǒng)運(yùn)動(dòng)方程和路面輸入方程寫成矩陣形式，即得 出系統(tǒng)的空間狀態(tài)方程：(7)X AX BU FW式中，X （& & Xb Xw Xg），為系統(tǒng)狀態(tài)矢量；W=（w（t）,為高斯白噪聲輸入矩 陣；U=（ Ua（t）,為輸入控制矩陣；00KsmbKsmb01mb0 A0KsKt KsKa10mwmwm

3、wBmw0100000F0010000000002f；0；2 . Gu2. LOG控制器設(shè)計(jì)車輛懸架設(shè)計(jì)中的主要指標(biāo)包括：代表輪胎接地性的輪胎動(dòng)載荷；代表輪胎舒適性的車身垂向振動(dòng)加速度；影響車身姿態(tài)且與輪胎布置有關(guān)的懸架動(dòng)行程。因此，LQG控制器設(shè)計(jì)中的性能指標(biāo) J即為輪胎動(dòng)位移、懸架動(dòng)行程和車身垂向振動(dòng)加速度的加權(quán)平方和在時(shí)域T內(nèi)的積分值，其表達(dá)式為:TJ 何 1 5Xw(t) Xg(t)2玨區(qū)Xw(t)2 q3Xb2(t)dt(8)T 0式中，ql、q2和q3分別為輪胎動(dòng)位移、懸架動(dòng)行程和車身垂向振動(dòng)加速度的加權(quán)系數(shù)。加 權(quán)系數(shù)的選取決定了設(shè)計(jì)者對(duì)懸架性能的傾向，如對(duì)車身垂向振動(dòng)加速度項(xiàng)選

4、擇較大的權(quán)為方便起見，這里取車身垂向振動(dòng)加速度的加權(quán)值，則考慮更多的是提高車輛操縱穩(wěn)定性。 值 q3=1。將性能指標(biāo)J的表達(dá)式（8）改寫成矩陣形式，即：1 T tTT(9)J Jim （XtQX UtRU 2XTNU）dt T 0式中，000000K；00q22QmbK；00q22 qmb0000000K；q220mbK；Nq22q叫c1R 2qq .mb .Ks當(dāng)車輛參數(shù)值和加權(quán)系數(shù)值確定后，最優(yōu)控制反饋增益矩陣可有黎卡提（Riccati）方程求出，其形式如下：(10)PA AtP (PB N)R 1(BtP Nt) Q 0最優(yōu)反饋控制增益矩陣 K BTP NT，由車輛參數(shù)和加權(quán)系數(shù)決定。

5、根據(jù)任意時(shí)刻的反饋狀態(tài)變量X （ 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。表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

6、.0x10-6車速Um/s20下截止頻率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）， 其中M為生成矩陣的行數(shù)，N為列數(shù)，P為白噪聲的功率（單位為dB）。本例中取M=10001，N=1，P=20。這意味著仿真計(jì)算中去一條白噪聲，共10001個(gè)采集點(diǎn)，噪聲強(qiáng)度為 20dB。設(shè)定采樣時(shí)間為 0.005s、車速為20m/s時(shí)，相當(dāng)于仿真路面

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

8、987.52627.23642.47450.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 仿真過程1) 生成路面輸入模型代碼如下：a=wgn(10001,1,20);t=0:0.005:50;road_)=t;road_)=a;save road_2) 參

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

10、;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)變量 XXb,Xw,Xb,Xw,Xg輸入與系統(tǒng)模塊，如下圖：Gain輸

11、出模塊，如下圖:Drrvative車身加速虞O懸架動(dòng)行程 simoutTo Workspace輪胎動(dòng)位移踣面輸入d ?車身加速度一(Xb)dt懸架動(dòng)行程 xb x 輪胎動(dòng)位移 xw xg 路面輸入 Xg整體程序框圖如下:Gainroadjle)C -y- Cif+Du錯(cuò)出況二r呈FromWort(5p3ce du.iirDei VHlivt!O1x - Aa*0u y = SDuSte-Spa-2tJGa n1djd:D出ivaciwdO4 )結(jié)果分析可以直接通過雙擊scope查看輸出的波形圖，為更好比較主動(dòng)懸架與被動(dòng)懸架的差別，下面通過輸出到workspace的狀態(tài)變量編程繪圖并計(jì)算均方根值

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

13、主動(dòng)懸架動(dòng)行程%sws1-被動(dòng)懸架動(dòng)行程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)./(le ngth(sws).A0.5);SWS1=n orm(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

14、)-X.data(:,5);dtd 1=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)./(le ngth(dtd).A0.5);DTD仁 norm(1000*dtd1,2)./(le ngth(dtd1).A0.5);結(jié)果如下： 車身加速度曲線円度速加身車1050-5-10IJ I1 11/ Bl 1 r：-.佃 Mrj IE* i卜*|.Jl 鬧J.4 和 1 .i p 111 *!卜 T ；1k丁/1r dV tf忙1 i

15、f j主動(dòng)懸架車身加速度01234567時(shí)間被動(dòng)懸架車身加速度8910卩度速加身車20100-10-20FiJlII，山“MiHr-機(jī)丫21fl ft 卩1帀Hfl1l01234567時(shí)間8910 懸架動(dòng)行程曲線i/iili T皿.1: lih KJ1 1 h I j.J ,1卩11-叫N忖1 *，卩L j1-.:f1 I10.040.020-0.02-0.040 1 2主動(dòng)懸架輪胎動(dòng)位移34567時(shí)間被動(dòng)懸架輪胎動(dòng)位移8910卩移位b 1ft機(jī)1 * J1丹fiJk. Af旨 1L尸*e / jJf.飛1主動(dòng)懸架動(dòng)行程0.10.050-0.05-0.11 23456時(shí)間被動(dòng)懸架動(dòng)行程7891

16、0卩移位00.60.40.20-0.2ji; 1 ik ,a r 1J i:,丄嘰I,Ij.fl1-If .r4-A1 tI i1 23456時(shí)間789010輪胎動(dòng)位移3456時(shí)間789100.10.050-0.050 1 2-0.1主動(dòng)懸架與被動(dòng)懸架性能指標(biāo)均方根值比較性能指標(biāo)主動(dòng)懸架被動(dòng)懸架車身加速度BA(m/s 2)1.514.60懸架動(dòng)行程SWS(mm)34.4267.10輪胎動(dòng)位移DTD(mm)5.8734.575. 半車模型建模及仿真半車模型建模懸架的半車模型有四個(gè)自由度，可選取前后輪垂直位移、 前后懸架與車身連接處垂直位移四個(gè)自由度，寫出狀態(tài)空間方程： 以X Z4,Z3,Z2,Z

17、i,Z4,Z3,Z2,Zi,Xio, X20作為系統(tǒng)狀態(tài)變量，狀態(tài)空間方程如下:X(t) AX(t) BU(t) FW(t)UafUar,為控制輸入矩陣，即前后懸架作動(dòng)器力;W W1 ,為路面輸入模型的高斯 白噪聲輸入。 w2LQG控制器選取如下二次型性能指標(biāo)：TimToqi(ZiZ01)2q2(Z2 Z1)2? 21 Z2q3(Z3Z02 )2? 22Z4 dtq1 ,q3,前后懸架輪胎動(dòng)位移加 權(quán)系數(shù); q2,q4,前后懸架動(dòng)行程加權(quán)系 數(shù)；1, 2,車身加速度加權(quán)系數(shù)。寫成矩陣形式:utru2XtNU )dt1 T lim (XtQX T T1 0由黎卡提方程求出反饋矩陣 K，則U KX

18、，因此狀態(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。程序代碼如下：%生成路面模型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=

19、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/I;a2=1/mb-a*b/I;a3=1/mb+aA2/I;A=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 0

20、 0 0 0 0 -Ksf/mwf (Ksf-Ktf)/mwf 0 Ktf/mwf;1 0 0 0 0 0 0 0 0 0;0 1 0 0 0 0 0 0 0 0;0 0 1 0 0 0 0 0 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=a2 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 b

21、1=a3A2+a2A2;b2=a2*a3+a1*a2; b3=a2A2+a1A2;sq=q4+KsrA2*b3 -q4-KsrA2*b3 Ksf*Ksr*b2 -Ksf*Ksr*b2 0 0; -q4-KsrA2*b3 q3+q4+KsrA2*b3 -Ksf*Ksr*b2 Ksf*Ksr*b2 -q3 0;Ksf*Ksr*b2 -Ksf*Ksr*b2 q2+KsfA2*b1 -q2-KsfA2*b1 0 0;-Ksf*Ksr*b2 Ksf*Ksr*b2 -q2-KsfA2*b1 q1+q2+KsfA2*b1 0 -q1;0 -q3 0 0 q3 0;0 0 0 -q1 0 q1;Q=zeros

22、(4,10);zeros(6,4) sq;R=b1 b2;b2 b3;N=zeros(4,2);Ksr*b2Ksr*b3;-Ksr*b2-Ksr*b3;Ksf*b1-Ksf*b2;zeros(2,2);K=lqr(A,B,Q,R,N);%求主動(dòng)懸架的時(shí)域響應(yīng)load road_w1=road_);w2=w1;t=road_);C=eye(10);D=zeros(10,2);Y,X=lsim(A-B*K,F,C,D,w1;w2,t);輸出矩陣C為10X10單位矩陣，即將所有狀態(tài)均輸出，因此系統(tǒng)輸出X是完全相同的。Ksf*b2;-Ksf*b1Y與系統(tǒng)狀態(tài)變量5.2預(yù)瞄控制預(yù)瞄控制中，前后輪的路面輸

23、入不再相同，后輪的輸入比前輪輸入滯后個(gè)預(yù)瞄時(shí)間，及前輪路面輸入w1與后輪路面輸入w2滿足:w2采用PaDa近似后，可轉(zhuǎn)化為時(shí)域方程:(t) A (t) Bw/t)則預(yù)瞄控制系統(tǒng)空間狀態(tài)方程可寫為:A FD X BFEUw-i0 A0B同樣可由黎卡提方程求出其反饋矩陣程序代碼如下: %預(yù)瞄控制 tao=0.14;a0=12/taoA2; 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)

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

25、5;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)40.5);BAO=norm(ba0,2)./(le ngth(ba0)90.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)./(