行駛動力學建模、仿真及主動懸架控制器設計

1、目錄1. 計算機仿真系統(tǒng)模型的建立.- 1 -2. LOG 控制器設計.-2 -3. 計算實例.-3 -4. MATLAB 仿真過程.-4 -5. 半車模型建模及仿真.-8 -5.1 隨機線性最優(yōu)控制.-9 -5.2 預瞄控制 .-11-5.3 結果比較.-12 -Word 文檔以單輪車輛模型為例， 介紹行駛動力學計算機建模、 仿真分析以及利用線性二次最優(yōu)控制理論進行主動懸架 LQG 控制器設計過程。1. 計算機仿真系統(tǒng)模型的建立根據(jù)圖 7 所示的主動懸架單輪車輛模型，運用牛頓運動定律，建立系統(tǒng)的運動方程，即：&U aK s ( xbxw)mb xb（ 4）&U aK s ( xbxw ) K

2、t (xwX g )mw xw（ 5）這里，采用一個濾波白噪聲作為路面輸入模型，即：&2f0 xg (t )2 G0uw(t )xg (t )（ 6）式中，xg 為路面垂向位移 （ m）；G0 為路面不平度系數(shù) （ m3/cycle ）；u 為車輛前進速度 （ m/s）； w 為數(shù)字期望為零的高斯白噪聲； f0 為下截止頻率（ Hz）。圖 7單輪車輛模型結合式（ 4）、式（ 5）和式（ 6），將系統(tǒng)運動方程和路面輸入方程寫成矩陣形式，即得出系統(tǒng)的空間狀態(tài)方程：&AXBU FW（ 7）X&xbxwxg )TW= （ w（t），為高斯白噪聲輸入矩式中， X( xbxw，為系統(tǒng)狀態(tài)矢量；陣； U=（

3、 Ua （t），為輸入控制矩陣；Word 文檔00K sK s01mbmbmb000K sKt K sK a10mwmwmwBmwA100000F0010000000002 f0；0；2 G0u2. LOG 控制器設計車輛懸架設計中的主要指標包括：代表輪胎接地性的輪胎動載荷； 代表輪胎舒適性的車身垂向振動加速度；影響車身姿態(tài)且與輪胎布置有關的懸架動行程。因此，LQG控制器設計中的性能指標J 即為輪胎動位移、 懸架動行程和車身垂向振動加速度的加權平方和在時域 T 的積分值，其表達式為：T&212q2 xb (t) xw (t )2(t) dtJ lim q1 xw (t) xg (t)q3 xb

4、TT 0（ 8）式中， q1、q2 和 q3 分別為輪胎動位移、懸架動行程和車身垂向振動加速度的加權系數(shù)。加權系數(shù)的選取決定了設計者對懸架性能的傾向，如對車身垂向振動加速度項選擇較大的權值，則考慮更多的是提高車輛操縱穩(wěn)定性。為方便起見， 這里取車身垂向振動加速度的加權值 q3=1。將性能指標J 的表達式（ 8）改寫成矩陣形式，即：J lim1 T (XTQXUTRU2 X T NU )dtTT 0（ 9）式中，000000000000K s2q2K s200Qq2mb20mb200q2K s2K s2q1NK amb2q1 q21K smb2R000q1q1mb20；當車輛參數(shù)值和加權系數(shù)值確

5、定后，最優(yōu)控制反饋增益矩陣可有黎卡提（Riccati ）方程求出，其形式如下：PA ATP (PB N)R 1(BTP NT) Q 0（ 10）最優(yōu)反饋控制增益矩陣KBT PN T ，由車輛參數(shù)和加權系數(shù)決定。根據(jù)任意時刻的Word 文檔反 饋 狀 態(tài) 變 量 X （ t ）， 就 可 得 到 t時 刻 作 動 器 的 最 優(yōu) 控 制 力 Ua ， 即 ：U a (t )KX (t)（ 11）3. 計算實例這里， 以某轎車的后懸架為例，給出一個完整的計算實例，包括車輛模型參數(shù)、仿真路面輸入?yún)?shù)、 控制器的設計參數(shù)以及計算結果。此例中車輛以20m/s 的速度在某典型路面上行駛，仿真時間T=50s

6、。計算中輸入的各參數(shù)及數(shù)值詳見表2。表 2 單輪車輛模型仿真輸入?yún)?shù)值車輛模型參數(shù)符號單位數(shù)值簧載質(zhì)量mbKg320非簧載質(zhì)量mwKg40懸架剛度KsN/m20000輪胎剛度KtN/m200000懸架工作空間SWScmm100仿真路面輸入?yún)?shù)符號單位數(shù)值路面不平度系數(shù)G0m3/cycle5.0x10-6車速Um/s20下截止頻率f0Hz0.1性能指標加權參數(shù)符號單位數(shù)值輪胎動位移q180000懸架動行程q25車身加速度q31仿真計算中以式（6）所示的濾波白噪聲作為路面輸入模型。白噪聲的生成可直接調(diào)用MATLAB函數(shù) WGN（ M ，N ，P）（此函數(shù)需要安裝信號處理工具箱Communicati

7、ons toolbox ），其中 M 為生成矩陣的行數(shù)，N 為列數(shù)， P 為白噪聲的功率 （單位為dB）。本例中取M=10001 ，N=1 ，P=20 。這意味著仿真計算中去一條白噪聲，共10001 個采集點，噪聲強度為20dB 。設定采樣時間為0.005s、車速為 20m/s 時，相當于仿真路面長度為1000m ，仿真時間為50s。根據(jù)建立的系統(tǒng)狀態(tài)方程式（7）及最優(yōu)化性能指標函數(shù)式（9），利用已知的矩陣A 、B、 Q、 R、 N ，調(diào)用 MATLAB中的線性二次最優(yōu)控制器設計函數(shù)K ， S， E=LQR （ A ， B，Q， R， N），即可完成最優(yōu)主動懸架控制器的設計。輸出的結果中，K

8、為最優(yōu)控制反饋增益矩陣， S 為黎卡提方程的解，E 為系統(tǒng)閉環(huán)特征根。根據(jù)表 2 給出的仿真輸入?yún)?shù)，本例中求得的最優(yōu)反饋增益矩陣K 為：K= （ 711.88-1241.5-19284-2038.520864）同時，還得到了黎卡提方程的解：2.45590.02892.47458.66077.3090.02890.48860.02987.52627.2364S2.47450.02984.97448.67545.10338.66077.52628.67542710.12700.47.3097.23645.10332700.42693.7在相同的仿真條件下，可將所設計的主動懸架系統(tǒng)與一個被動系統(tǒng)進

9、行對比分析。在被動懸架系統(tǒng)中，取懸架剛度Ks=22000N/m ，阻尼系數(shù)Cs=1000NS/m 。除此之外，其他輸入?yún)?shù)Word 文檔值均與主動懸架系統(tǒng)完全相同。4. MATLAB仿真過程1）生成路面輸入模型代碼如下：a=wgn(10001,1,20);t=0:0.005:50;road_file(:,1)=t;road_file(:,2)=a;saveroad_fileroad_file2）參數(shù)輸入代碼如下：loadroad_file.mat%載入路面數(shù)據(jù)模型Ks=22000;mb=320;Kt=200000;mw=40;f0=0.1;G0=0.000005;u=20;Kb=20000;K

10、s1=22000;Cs=1000;%輸入仿真有關參數(shù)A=0,0,-Ks/mb,Ks/mb,0;%建立主動懸架的狀態(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; %建立被動懸架的狀態(tài)矩陣 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

11、;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;Word 文檔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 , xb , xw , xg 輸入與系統(tǒng)模塊，如下圖：輸出模塊，如下圖：車身加速度懸架動行程輪胎動位移路面輸入d?( xb )dtxbxwxwxgxgWord 文檔整體程序框圖如下：4）結果分析可以直接通過雙擊

12、scope 查看輸出的波形圖， 為更好比較主動懸架與被動懸架的差別， 下面通過輸出到 workspace 的狀態(tài)變量編程繪圖并計算均方根值。代碼如下：% 繪制車身加速度曲線，并計算均方根值%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

13、.5);BA1=norm(ba1,2)./(length(ba1).0.5);% 繪制懸架動行程曲線，并計算其均方根值%sws- 主動懸架動行程%sws1- 被動懸架動行程figure()sws=X.data(:,3)-X.data(:,4);Word 文檔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

14、.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);結果如下：車身加速度曲線主動懸架車身加速度10)2s/5m(度0速加身-5車-1

15、0012345678910時 間(s)被動懸架車身加速度20)2s/10m(度0速加身 -10車-20012345678910時 間 (s)懸架動行程曲線Word 文檔)m(主動懸架動行程0.10.05移0位-0.05-0.1)m(012345678910時 間(s)被動懸架動 行 程0.60.4移 0.2位0-0.2012345678910時 間 (s)輪胎動位移主動懸架輪胎動位移0.040.02)m(移0位-0.02-0.04012345678910時 間 (s)被動懸架輪胎動位移0.10.05)m(移0位-0.05-0.1012345678910時 間 (s)主動懸架與被動懸架性能指標均

16、方根值比較性能指標主動懸架被動懸架車身加速度 BA(m/s 2)1.514.60懸架動行程 SWS(mm)34.4267.10輪胎動位移 DTD(mm)5.8734.575. 半車模型建模及仿真半車模型建模Word 文檔懸架的半車模型有四個自由度， 可選取前后輪垂直位移、 前后懸架與車身連接處垂直位移四個自由度，寫出狀態(tài)空間方程：?以 X z4 , z3 , z2 , z1, z4 , z3 , z2 , z1, x10, x20 作為系統(tǒng)狀態(tài)變量，狀態(tài)空間方程如下：?X (t )AX (t )BU (t )FW (t)UU af ,為控制輸入矩陣，即前 后懸架作動器力；U arW w1 ,

17、為路面輸入模型的高斯 白噪聲輸入。 w2LQG 控制器選取如下二次型性能指標：J lim 1T? 2? 2 q1( z1 z01 ) 2q2 (z2 z1) 21 z2q3 ( z3 z02 )22 z4 dtTT 0q1 ,q3 , 前后懸架輪胎動位移加權系數(shù)；q2 ,q4, 前后懸架動行程加權系數(shù)；1 ,2 ,車身加速度加權系數(shù)。寫成矩陣形式：TJlim 1( X T QX U T RU 2X T NU ) dtTT 0由黎卡提方程求出反饋矩陣K，則 UKX ，因此狀態(tài)空間方程可寫成：?X (t ) ( ABK ) X (t)FW (t )( ABK ) - 狀態(tài)矩陣F 輸入矩陣由此形式可

18、方便求出其時域與頻域響應。5.1 隨機線性最優(yōu)控制路面模型本例仿真車速為20m/s，軸距為2.8m，因此滯后時間為0.14s。若取仿真時間T 為 20s，采樣時間間隔為0.005s，因而仿真點數(shù)4028。程序代碼如下：% 生成路面模型a=wgn(4029,1,20);t=0:0.005:20.14;r=a,t;saveroad_filerWord 文檔% 車身模型參數(shù)輸入clearclcmb=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;

19、f0=0.1;% 性能指標加權系數(shù)q1=80000;q2=100;q3=80000;q4=100;% 半車懸架模型建模，計算A、 B、 F矩陣a1=1/mb+b2/I;a2=1/mb-a*b/I;a3=1/mb+a2/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 0 0 0 0 -Ksf/mwf (Ksf-Ktf)/mwf 0 Ktf/mwf;1000000000;01

20、00000000;0010000000;0001000000;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控制器設計，計算 Q、 R、 N矩陣并求出反饋矩陣 K b1=a32+a22;b2=a2*a3+a1*a2;b3=a22+a12;sq=q4+Ksr2*b3 -q4-Ksr2*b3 Ksf*Ksr*b2 -Ksf*Ksr*b2 0 0;

21、 -q4-Ksr2*b3 q3+q4+Ksr2*b3 -Ksf*Ksr*b2 Ksf*Ksr*b2 -q3 0;Word 文檔Ksf*Ksr*b2 -Ksf*Ksr*b2 q2+Ksf2*b1 -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 Ksr*b3;-Ksr*b2 -Ksr*b3;Ksf*b1 Ksf*b2;-Ks

22、f*b1 -Ksf*b2;zeros(2,2);K=lqr(A,B,Q,R,N);% 求主動懸架的時域響應loadroad_file.matw1=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 單位矩陣，即將所有狀態(tài)均輸出，因此系統(tǒng)輸出Y 與系統(tǒng)狀態(tài)變量X 是完全相同的。5.2 預瞄控制預瞄控制中， 前后輪的路面輸入不再相同，后輪的輸入比前輪輸入滯后一個預瞄時間，及前輪路面輸入w1 與后輪路面輸入w2 滿足：w2se采用 PaDa

23、近似后，可轉(zhuǎn)化為時域方程：?(t )A(t )B w1(t )則預瞄控制系統(tǒng)空間狀態(tài)方程可寫為：?AFD XFEXB?0AUw10B同樣可由黎卡提方程求出其反饋矩陣K p程序代碼如下：% 預瞄控制tao=0.14;a0=12/tao2;a1=6/tao;a2=1;Word 文檔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,z

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