1、1、系統(tǒng)辨識部分1.1 列車制動系統(tǒng)介紹制動系統(tǒng)是列車操縱系統(tǒng)的重要組成部分，它是用來調(diào)節(jié)車速和進(jìn)行停車，是列車安全可靠運行的基本保障。現(xiàn)代列車通常用電空混合制動系統(tǒng)，它通過制動控制器協(xié)調(diào)制動和空氣制動的介入比例，基本能夠保證在車速變化情況下一致的制動性能20當(dāng)司機(jī)操縱控制手柄時，制動控制器收到司機(jī)請求的目標(biāo)加速度，它根據(jù)車重和車速進(jìn)行相應(yīng)補(bǔ)償，然后計算出所需制動力和分配方案，最后驅(qū)動執(zhí)行機(jī)構(gòu)完成制動力的輸出。這是一個反饋控制的過程，但是由于它并不補(bǔ)償由于坡道和彎道所帶來的附加影響。1.2 收集先驗知識ATO(自動駕駛系統(tǒng))的功能是代替司機(jī)駕駛，它根據(jù)目標(biāo)速度、行車許可及線路情況等自動產(chǎn)生牽引

2、或制動指令。ATO控制器和司機(jī)對車輛的控制都是通過車輛本身的牽引控制系統(tǒng)和制動控制系統(tǒng)完成的。對于ATO來說制動模型是包括列車制動控制系統(tǒng)在內(nèi)的車輛運行模型，它描述從輸入到制動控制器的制動指令到車輛運行狀態(tài)：加速度、速度間的動態(tài)關(guān)系。圖1是制動模型框圖12o圖1制動模型框圖1.3 數(shù)學(xué)模型分析3制動控制器通過反饋調(diào)節(jié)實現(xiàn)對目標(biāo)加速度的跟蹤，這是一個動態(tài)的過程。這個過程可以用一階動態(tài)系統(tǒng)近似，另外考慮到系統(tǒng)傳輸延時，則可以用下面方程進(jìn)行描述U11.a(t)=-a(t)+At。)(1)TTA式中，a(t)為控制加速度，它是由制動系統(tǒng)控制器的作用而是列車產(chǎn)生的加速度；A為目標(biāo)加速度；丁為系統(tǒng)響應(yīng)時間

3、常數(shù)；仃為傳輸延遲。車輛實際加速度a(t)由控制加速度a、t)和環(huán)境(如彎道、坡道)造成的附加加速度|_a：t)構(gòu)成，車輛速度v由實際速度決定。(2)a(t)=a(t)La(t)Lv(t)=a(t)(3)目標(biāo)加速度A(t)為列車制動控制器的輸入指令，它通過駕駛員或ATO控制(4)u(t)(制動指令)產(chǎn)生，它們之間的關(guān)心可以用式(4)靜態(tài)函數(shù)關(guān)系描述。At)=F(u(t)以上方程描述了從ATO控車角度看到的列車制動模型，它可以圖2框圖表示。a(t)圖2制動模型框圖其中輸入為駕駛員或ATO控制指令u(t),輸出為車輛實時速度v(t),式(1)中e和仃是需要辨識的參數(shù)。辨識式(1)中參數(shù)需要知道控制

4、加速度3t)，但是由于|_a(t)的干擾，無法從測量的實時速度中獲得。我們假設(shè)列車是在平直軌道上運行，環(huán)境影響可忽略不計的理想狀態(tài)，即：Ua(t)-0(5)AL那么a(t)定V(t)(6)1.4穩(wěn)態(tài)響應(yīng)參數(shù)估計6暫態(tài)響應(yīng)消失后a(t)=A(to)=F(u(t-cr)ta仃(7)可以看出，如果u(t)為一恒值，則A(t)也為恒值。根據(jù)式(5),可以利用測量的速度來計算在各制動級別上的穩(wěn)態(tài)加速度，得到如下關(guān)系：Lv(t)=A(t-0)=F(u(t-二)=F(u)t>仃i=1,2,3,，N(8)實驗測量數(shù)據(jù)為:U(t)12345A(t)-20.7-39.5-58.4-77.396.167-11

5、5.0-133.8用matlab對其進(jìn)行線性擬合，代碼為：clearu=1,2,3,4,5,6,7;A=-2O7,-39.5,-58.4,-77.3,-96.1,-115Q-133.8;fori=1:7;Z(i)=A(i)；endZL=Z'HL=u(1),1;u(2),1;u(3),1;u(4),1;u(5),1;u(6),1;u(7),1c1=HL'*HL;c2=inv(c1);c3=HL'*ZL;c4=c2*c3a=c4(1)b=c4(2)程序運行結(jié)果：ZL=-20.7000-39.5000-58.4000-77.3000-96.1000-115.0000-133.

6、8000HL=c4 =-18.8571-1.8286-18.8571-1.8286即可得1111A t ) = 1 8 . 8i6t ( )1(9)1.5動態(tài)響應(yīng)參數(shù)辨識與仿真A為了方便計算，我們將A(t)作為輸入，a(t)作為輸出。式(1)描述了從目標(biāo)加速度到控制加速度以7和仃為參數(shù)的傳遞函數(shù)，表示為：Gs；i,仃)，s為拉普拉斯變換變量。最優(yōu)的參數(shù)應(yīng)該是使輸入u得到響應(yīng)與實際響應(yīng)誤差e(t)最小。其中:(10)(11)e-s0s,Le(t)=g(t；,二)LF(u(t)-v(t)式中g(shù)(s；工,*是0；工,。)的反變換。誤差大小用式(11)進(jìn)行評估13。J(.,二)=e2(t )dt(12

7、)取邙=1.2s,0.4s,貝ij-1.2sGs,二)二(13)e-0.4s1對式(13)進(jìn)行Z變換，取采樣時間為T=1,得到離散系統(tǒng)的傳遞函數(shù):Gz)2 0.8647 z 0.05325z -z 一 0.08208(14)則有y(k)0.08208y(k1)=0.8647u(k-2)+0.05325Uk-3)+v(k)(15)式中，v(k)是方差為0.1的白噪聲干擾。制動控制系統(tǒng)為時滯系統(tǒng)，數(shù)據(jù)不能突變，將輸入信號定為階躍信號，用遞推最小二乘法進(jìn)行參數(shù)估計。算法步驟：1 .設(shè)置3(0)和P(0),輸入初始數(shù)據(jù)；2 .采樣當(dāng)前輸出y(k)和輸入u(k);3 .利用最小二乘法估計$遞推公式，計算

8、K(k),Q(k)和P(k);4 ,kTk+1返回第二步，繼續(xù)循環(huán)4。仿真程序如下：clearall;closeall;a=1-0.08208;b=0.86470.05325'd=2;%寸象參數(shù)na=length(a)-1;nb=length(b)-1;%na、nb為A、B的階次L=500;%仿真長度y1=1;y2=1;y3=1;y4=1;%生成M序列fori=1:L;%生成M序列x1=xor(y3,y4);x2=y1;x3=y2;x4=y3;y(i)=y4;ify(i)>0,u(i)=y(i);elseu(i)=0.5;endy1=x1;y2=x2;y3=x3;y4=x4;en

9、duk=zeros(d+nb,1);%輸入初值：uk（i）表示u（k-i）yk=zeros(na,1);%輸出初值xi=sqrt(0.1)*randn(L,1);%白噪聲干擾theta=a(2:na+1);b;%對象參數(shù)真值%thetae 初值thetae_1=zeros(na+nb+1,1);P=10A6*eye(na+nb+1);fork=1:Lphi=-yk;uk(d:d+nb);%此處phi 為列向量y(k)=phi'*theta+xi(k);%采集輸出數(shù)據(jù)%遞推最小二乘法K=P*phi/(1+phi'*P*phi);thetae(:,k)=thetae_1+K*(y(

10、k)-phi'*thetae_1);P=(eye(na+nb+1)-K*phi')*P;%更新數(shù)據(jù)thetae_1=thetae(:,k);fori=d+nb:-1:2uk(i)=uk(i-1);enduk(1)=u(k);fori=na:-1:2yk(i)=yk(i-1);endyk(1)=y(k);endplot(1:L,thetae);xlabel('k');ylabel('2?ey1a?a?0b');); axis(0 L -2 4);legend('a_1','b_0','b_1'對系統(tǒng)進(jìn)

11、行RIS辨識仿真,參數(shù)辨識如下:3210q盤較期O50O400025050圖3白色噪聲方差為0.01時參數(shù)辨識結(jié)果圖4白色噪聲方差為0.1時參數(shù)辨識結(jié)果q百工±=若卷2J,，1J.1.1.1.LL.050100150200250300350400450500k圖5白色噪聲方差為1時參數(shù)辨識結(jié)果以上參數(shù)曲線說明隨著擾動強(qiáng)度的增強(qiáng)，辨識的效果會有所下降。2、自適應(yīng)控制部分2.1基于最小方差間接自校正控制基于前節(jié)討論的列車制動模型，采用最小方差間接自校正控制。當(dāng)被控對象的參數(shù)未知時，可首先利用遞推增廣最小二乘法在線實時估計對象參數(shù)，然后再設(shè)計最小方差控制律，即將對象參數(shù)估計器和控制器的設(shè)計

12、分開進(jìn)行，就形成了最小方差自校正控制間接算法網(wǎng)4。2.3控制器設(shè)計根據(jù)前文，現(xiàn)對最小方差自適應(yīng)控制算法進(jìn)行推導(dǎo)。1)、單步輸出預(yù)測45：結(jié)合本系統(tǒng)，令輸出為y(k),輸入為u(k)。設(shè)系統(tǒng)數(shù)學(xué)模型為：-4,-d_1一，一、A(z)y(k)=zB(z)u(k)+C(z卜(k)(16)式中，C(z')為Hurwitz多項式，即其零點完全位于z平面的單位圓內(nèi)；u(k)和y(k)表示系統(tǒng)的輸入輸出，”k)為方差仃2的白噪聲，d>1是系統(tǒng)純延時，且有Ia(z,)=1,a1z,-a2z2,，anzzJaJb(z,)=b0+b1z，+b2z/+.+bnbz"b,b0#0(17)C(z

13、,)=1CiZJ-az，c%z"c對象(2.1)基于k時刻和之前時刻的輸入輸出數(shù)據(jù)記作Yk,Uk=y(k),y(k1),.,u(k),u(k1),.(18)基于Wk,Uk對k+d時刻的輸入記作?(k+d|k)(19)(輸出)預(yù)測誤差記作?(k+d|k)=y(k+d)-?(k+d|k)(20)那么，關(guān)于提前d步最小方差預(yù)測輸出可由如下定理給出。2)、定理(最優(yōu)d步預(yù)測輸出)45使如下性能指標(biāo)(預(yù)測誤差的方差)Ey2(k+d|k)為最小的d步最優(yōu)預(yù)測輸出y*(k+d|k)必滿足方程*1C(z)y(kd|k)=G(z)y(k)F(z)u(k)式中1_d1C(z)=A(z)E(z)z&quo

14、t;G(z)_1F(z)=B(z)E(z)且Rz,)=1+ez，+.+,zf,(ne=d1)G(z')=g0gz'.gnz,(ng=na-1)gF(z)=f°f1z.z,(%=踴d-1)此時，最有預(yù)測誤差的方差為d1Ey*(kd|k)2=(1八e2)-23)、最小方差控制律45(21)(22)(23)(24)i1設(shè)B(z,)也是Hurwitz多項式，則其零點完全位于z平面單位圓內(nèi)，即對象是最小相位或者逆穩(wěn)，可得如下定理：設(shè)控制目標(biāo)是使實際輸出y(k+d)跟蹤期望輸出yr(k+d),使性能指標(biāo)J=Ey(k+d)-y.(k+d)2(25)為最小。由上面定理可知：.J-*-

15、y(k+d)=E«k+d)+y(k+d|k)(26)將上式代入式(25)中，得2.J二或E(kd)y(kd|k)yr(kd)=EE(k+d)2+Ey*(k+d|k)-yr(k+d)2(27)上式右邊第一項不可控，若使J最小，需使上式右邊第二項為0,即*一_-,_.y(k+d|k)=yr(k+d)(28)將式(28)代入最優(yōu)預(yù)測輸出方程式(21)得Qz)y*(k+d|k)=Gz,)y(k)+F(z)u(k)(29)得Ctz)yr(k+d)=G(z,)y(k)+F(z)Uk)(30)上式經(jīng)過簡單的變形即得最小方差控制律為F(z')Uk)=az)y.(k+d)G(z,)y(k)(3

16、1)2.4仿真驗證用matlab對上面自適應(yīng)控制進(jìn)行仿真，其步驟如下5：已知：模型階次na，幾，幾及純延時do第一步：設(shè)置初值日(0)和R0),輸入初始數(shù)據(jù)；第二步：采樣當(dāng)前實際輸出y(k)和期望輸出yr(k+d)；第三步：利用遞推增廣最小二乘法在線估計被控對象參數(shù)即A,B和C;第四步：求解Sindiophantine方程，得到多項式E、F、G的系數(shù)；第五步：利用式(22)計算并實施u(k)；第六步：返回第二步(ktk+1),繼續(xù)循環(huán)。被控對象為y(k)一0.08208y(k-1)=0.8647u(k-2)0.05325Uk-3)(k)式中qk)為白噪聲。6取初值R0)=10I、40)=0.0

17、005,期望輸出yr(k)為幅值為5的萬波信號，采用最小方差自校正控制間接算法。由于參數(shù)估計初期，所獲得的估計參數(shù)不準(zhǔn)確，很可能偏離實際參數(shù)真值較遠(yuǎn)，使多項式B(z,)不穩(wěn)定，可以在實施自校正控制時，將B(z)的參數(shù)估計值限制在穩(wěn)定范圍內(nèi)6。仿真程序如下：clearall;closeall;clc;%模型參數(shù)%輸出階數(shù)a=1-0.08208;b=10.5;c=1;d=2;na=length(a)-1;nb=length(b)-1;nc=length(c)-1;nf=nb+d-1;%nf為多項式的階次L=400;%控制步數(shù)uk=0.91*ones(d+nb,1);臉入初值；u(k)表示u(k-1

18、)yk=4.8*ones(na,1);%出初值yrk=1*ones(nc,1);%11望輸出初值xik=zeros(nc,1);%白噪聲初值xiek=zeros(nc,1);yr=5*0.8*ones(L/2,1);-0.6*ones(L/2+d,1);%期望輸出xi=sqrt(0.1)*randn(L,1);硒噪聲序列%RELS初值thetae_1=0.0005*ones(na+nb+1+nc,1);%lt匕處不能為0P=10A6*eye(na+nb+1+nc)；fork=1:Ltime(k)=k;%采集輸出數(shù)據(jù)y(k)=-a(2:na+1)*yk+b*uk(d:d+nb)+c*xi(k);

19、%遞推增廣最小二乘法phie=-yk;uk(d:d+nb);K=P*phie/(1+phie'*P*phie);thetae(:,k)=thetae_1+K*(y(k)-phie'*thetae_1);P=(eye(na+nb+1+nc)-K*phie')*P;xie=y(k)-phie'*thetae(:,k);%白噪聲的估值%提取辨識參數(shù)ae=1thetae(1:na,k)'be=thetae(na+1:na+nb+1,k)'ce=1e,f,g=sindiophantine(ae,be,ce,d);%求解單步Diophantine方程u(k)

20、=(-f(2:nf+1)*uk(1:nf)+ce*yr(k+d:-1:k+d-min(d,nc);yrk(1:nc-d)-g*y(k);yk(1:na-1)/f(1);%更新數(shù)據(jù)thetae_1=thetae(:,k);fori=d+nb:-1:2uk(i)=uk(i-1);enduk(1)=u(k);fori=na:-1:2yk(i)=yk(i-1);endyk(1)=y(k);endfigure(1);subplot(2,1,1);plot(time,yr(1:L),'r:',time,y);xlabel('k');ylabel('y_r(k)?0y(k)');legend('y_r(k)','y(k)');axis(0L-55)