自適應(yīng)控制結(jié)課作業(yè)

1、自適應(yīng)控制結(jié)課作業(yè)2016年1月目錄1 遺忘因子遞推最小二乘法11.1最小二乘理論11.2帶遺忘因子的遞推最小二乘法1白噪聲與白噪聲序列1遺忘因子遞推最小二乘法21.3仿真實例4程序1 帶遺忘因子的遞推最小二乘法：62 廣義最小方差自校正控制82.1廣義最小方差自校正控制82.1.1 逆M序列82.1.2 廣義最小方差自校正控制82.2仿真實例10程序2 廣義最小方差自校正控制：133 參考模型自適應(yīng)控制163.1參考模型自適應(yīng)控制163.2仿真實例18數(shù)值積分18仿真結(jié)果18程序3 參考模型自適應(yīng)控制：231 遺忘因子遞推最小二乘法1.1最小二乘理論最小二乘最早的想法是高斯在1795年預(yù)測行

2、星和彗星運動軌道時提出來的，“未知量的最大可能的值是這樣一個數(shù)值，它使各次實際觀測和計算值之間的差值的平方乘以度量其精確度的數(shù)值以后的和為最小”1。這一估計方法原理簡單，不需要隨機變量的任何統(tǒng)計特性，目前已經(jīng)成為動態(tài)系統(tǒng)辨識的主要手段。最小二乘辨識方法使其能得到一個在最小方差意義上與實驗數(shù)據(jù)最好擬合的數(shù)學(xué)模型。由最小二乘法獲得的估計在一定條件下有最佳的統(tǒng)計特性，即統(tǒng)計結(jié)果是無偏的、一致的和有效的。1.2帶遺忘因子的遞推最小二乘法1.2.1白噪聲與白噪聲序列系統(tǒng)辨識中所用到的數(shù)據(jù)通常含有噪聲。從工程實際出發(fā)，這種噪聲往往可以視為具有理想譜密度的平穩(wěn)隨機過程。白噪聲是一種最簡單的隨機過程，是由一系

3、列不相關(guān)的隨機變量組成的理想化隨機過程。白噪聲的數(shù)學(xué)描述如下：如果隨機過程均值為0，自相關(guān)函數(shù)為，即式中，為單位脈沖函數(shù)(亦稱為Dirac函數(shù))，即，且則稱該隨機過程為白噪聲，其離散形式是白噪聲序列。如果隨機序列均值為零，且兩兩互不相關(guān)，即對應(yīng)的相關(guān)函數(shù)為：則這種隨機序列稱為白噪聲序列。其譜密度函數(shù)為常數(shù)。白噪聲序列的功率在到的全頻段內(nèi)均勻分布2。建立系統(tǒng)的數(shù)學(xué)模型時，如果模型結(jié)構(gòu)正確，則模型參數(shù)辨識的精度將直接依賴于輸入信號，因此合理選用辨識輸入信號是保證能否獲得理想的辨識結(jié)果的關(guān)鍵之一。理論分析表明，白噪聲作為被辨識系統(tǒng)的輸入時，可以激發(fā)系統(tǒng)的所有模態(tài)，可對系統(tǒng)充分激勵，可防止數(shù)據(jù)病態(tài)，保

4、證辨識精度，可以保證獲得較好的辨識效果。圖1-1 白噪聲序列1.2.2遺忘因子遞推最小二乘法假設(shè)被辨識的系統(tǒng)為一單入單出的離散時間系統(tǒng)，且已知為CAR模型，如圖2所示：(k)u(k)z-dB(z-1)A(z-1)y(k)圖1-2 辨識系統(tǒng)模型即式中，為輸入變量，為輸出變量，為白噪聲，且則上式可轉(zhuǎn)化為如下最小二乘格式：式中，為數(shù)據(jù)向量，待估參數(shù)向量，且取算法的性能指標(biāo)為式中，為遺忘因子（）。帶遺忘因子的遞推最小二乘估計的算法公式為：公式表明，新的參數(shù)估計是用新的實際測量值與基于老模型進行預(yù)測得到的量之偏差，對前面的參數(shù)估計加以修正得到的，修正系數(shù)陣為。的物理意義是參數(shù)估計誤差的方差，作為參數(shù)估計

5、精度的一種度量。遺忘因子的作用是削弱過去數(shù)據(jù)的作用，通常選擇0.95到0.998之間的數(shù)。帶遺忘因子的遞推最小二乘估計算法屬于在線辨識所用方法的一種，它既能克服離線辨識的缺點，也能克服遞推最小二乘估計中的“數(shù)據(jù)飽和”現(xiàn)象。遺忘因子最小二乘法的算法：已知式階次、式階次以及延遲。步驟1：設(shè)置初值和及遺忘因子，輸入初始數(shù)據(jù)；步驟2：采樣當(dāng)前輸出和輸入；步驟3：利用遞推公式，計算、和；步驟4：，返回步驟2，繼續(xù)循環(huán)。1.3仿真實例系統(tǒng)模型如下：其中，輸入為方差為1的白噪聲，為方差為0.1的白噪聲。由于和的選擇可按如下方法：式中，為充分大的正實數(shù)，為零向量或充分小的正的實向量。因此，取初值、。仿真結(jié)果如

6、下：圖 1-3 參數(shù)估計結(jié)果圖 1-4實際輸出與辨識輸出對比圖 1-5 實際輸出與辨識輸出對比（局部放大）附錄程序1 帶遺忘因子的遞推最小二乘法：clear all; close all;a=1 1.2 0.8 0.1' b=1 0.5 2 1.5' d=3;na=length(a)-1; nb=length(b)-1;L=1000;uk=zeros(d+nb,1); %輸入初值：uk(i)表示u(k-i)yk=zeros(na,1);%輸出初值ymk=zeros(na,1); %辨識輸出的初值u=randn(L,1); %輸入采用白噪聲序列xi=sqrt(0.1)*randn

7、(L,1); %白噪聲序列thetaek=zeros(na+nb+1,1); %thetae初值P=106*eye(na+nb+1);lambda=0.998;for k=1:L theta(:,k)=a(2:na+1);b;%對象參數(shù)真實值 phi=-yk;uk(d:d+nb); phim=-ymk;uk(d:d+nb); y(k)=phi'*theta(:,k)+xi(k); ym(k)=phim'*thetaek; %遞推公式： K=P*phi/(lambda+phi'*P*phi); thetae(:,k)=thetaek+K*(y(k)-phi'*th

8、etaek); P=(eye(na+nb+1)-K*phi')*P/lambda; thetaek=thetae(:,k); %產(chǎn)生新的phi和phim for i=d+nb:-1:2 uk(i)=uk(i-1); end uk(1)=u(k); for i=na:-1:2 yk(i)=yk(i-1); ymk(i)=ymk(i-1); end ymk(1)=ym(k); yk(1)=y(k);endfigure(1)plot(1:L,thetae(1:na+nb+1,:); hold on; plot(1:L,theta(1:na+nb+1,:),'k:');lege

9、nd('a_1','a_2','a_3','b_0','b_1','b_2','b_3'); axis(0 L -0.5 2.5);%x軸范圍，y軸范圍figure(2)plot(1:L,y);hold onplot(1:L,ym,'-r');legend('y','ym');axis(0 L -30 30);2 帶有輔助變量的最小方差控制器2.1帶有輔助變量的最小方差控制器2.1.1 逆M序列設(shè)M(k)是周期為Np bit、元素取值為

10、0或1的M序列，S(k)是周期為2bit、元素依次取值為0或1的方波序列，將這兩個序列按位進行異或運算，得到的復(fù)合序列就是周期為2Np bit、元素取值為0或1的逆M序列，記作IM(k)，即有IM(k)=M(k)S(k)將上述逆M序列的邏輯值“0”或“1”分別變換為1或1，此時逆M序列均值為0.雖然逆M序列只是M序列與方波序列簡單復(fù)合的結(jié)果，但其性質(zhì)卻優(yōu)于M序列，使其在辨識領(lǐng)域有著更為廣泛的應(yīng)用。2.1.2 廣義最小方差自校正控制器當(dāng)考慮干擾對系統(tǒng)的作用時，控制器的設(shè)計就是要最大限度的減小干擾對系統(tǒng)的影響。鑒于一般被控對象或過程都存在不同程度的純遲延，控制對系統(tǒng)的作用要到時刻才有響應(yīng)。在這段純

11、遲延的時間內(nèi)，干擾仍會作用于系統(tǒng)，所以在時刻預(yù)測時刻的輸出，并按照預(yù)測誤差的方差最小的原則，設(shè)計現(xiàn)時控制，并加以實施。當(dāng)過程參數(shù)未知，或者時變時，用遞推最小二乘法估計，或者直接估計控制器參數(shù)，然后算出控制量來，這就是最小方差自校正控制的基本思想。 而廣義最小方差自校正控制就是在性能指標(biāo)函數(shù)中引入加權(quán)多項式。被控對象為：其中：為延遲因子，為輸入變量，為輸出變量，為白噪聲。為Hurwitz多項式。選擇性能指標(biāo)函數(shù)為：式中，為期望輸出；為第拍的輸出；為第拍的控制；、和分別為輸出、期望輸出和控制的加權(quán)多項式，它們分別具有改善閉環(huán)系統(tǒng)性能，軟化輸入和約束控制量的作用。并且上述多項式的階次及參數(shù)根據(jù)實際需

12、要確定。由此，據(jù)文獻1知，廣義最小方差控制律為：uk=Cz-1Rz-1yrk+d-G(z-1)P(z-1)y(k)q0b0Cz-1Qz-1+F(z-1)P(z-1) (2-1)在進行控制系統(tǒng)設(shè)計時，一般可以取加權(quán)多項式、和，而大小的選取需要在快速性和穩(wěn)定性方面進行權(quán)衡。當(dāng)被控對象的參數(shù)未知時，也可以利用遞推算法直接估計廣義最小方差控制的參數(shù)，即廣義最小方差自校正控制直接算法。由式21知，在選定加權(quán)多項式P(z-1)、Rz-1、Qz-1后，只需再確定多項式G(z-1)、F(z-1)、Cz-1，即可確定廣義最小方差控制律u(k)?？刂破鲄?shù)的估計模型及遞推公式如下：k=k-1+Kkyk-Tk-dk

13、-1Kk=PK-1k-d+Tk-dPk-1k-dPk=1I-K(k)Tk-dP(k-1) (2-2)式中k-dyk-d,yk-d-ng,uk-d,uk-d-nf,-y*k-1,-y*k-ncTRng+nf+nc+2×1=g0,gng,f0,fnf,c1,cncTRng+nf+nc+2×1y*k=Tk-dk-d或y*k=yrk（2-3）由式21得廣義最小方差控制律為uk=1q02f0+f0q0f0q0-CQ+f0-FPuk+CRyrk+d-GPy(k) (2-4)由上式可知，在控制算法實施過程中，如果f0趨于零，則會出現(xiàn)零除現(xiàn)象。為此，應(yīng)對f0的最小值加以約束，這就意味著應(yīng)該

14、事先知道f0的符號和下界，或事先確定f0的值。2.1.3 帶有輔助變量的最小方差控制器廣義最小方差尚需解決偏差問題，因為采用廣義最小方差性能指標(biāo)中含有u2(k)，穩(wěn)態(tài)時必將使y2(k+d)產(chǎn)生某一數(shù)值，來補償u2(k)的值，而在穩(wěn)態(tài)時y(k+d)0產(chǎn)生了穩(wěn)態(tài)偏差，而且不便于處理設(shè)定值問題。定義一種輔助系統(tǒng)的廣義輸出為：xk=yk-yrk-d+uk-d (2-5)求使性能指標(biāo)函數(shù)：J=Ez2(k+d)為最小的控制u(k)為：uk=-1g0Hz-1yrk-Mz-1yk+i=1nggiuk-i (2-6)采用帶遺k=k-1+Kkzk-Tk-dk-1Kk=PK-1k-d+Tk-dPk-1k-dPk=1

15、I-K(k)Tk-dP(k-1) (2-7)式中：k-dyk-d,yk-d-nm,uk-d,uk-d-ng,-yrk-d,-yrk-d-nhTRnm+ng+nh+2×1=m0,mnm,g0,gng,h0,hnhTRnm+ng+nh+2×1（2-8）帶有輔助變量的最小方差控制器的算法：已知：模型階次、以及延遲。步驟1：設(shè)置初值和，輸入初始數(shù)據(jù)，并設(shè)置加權(quán)多項式、；步驟2：采樣當(dāng)前實際輸出和期望輸出yr(k)；步驟3：根據(jù)式28構(gòu)造觀測數(shù)據(jù)向量k-d，利用遺忘因子遞推最小二乘法遞推公式27在線實時估計控制器參數(shù)，即m、G和h;步驟4：利用式26計算并實施u(k);步驟5：返回步

16、驟2（），繼續(xù)循環(huán)。2.2仿真實例設(shè)系統(tǒng)模型如下：yk-1.8yk-1+0.7yk-2+0.1yk-3=uk-4+0.2uk-5+k+0.1(k-1)其中，為取值為0.1的逆M序列，采用帶有輔助變量的最小方差控制器。取初值、；設(shè)置加權(quán)多項式、Rz-1=1、2。期望輸出采用幅值為10的方波型號，其控制結(jié)果如下：圖21期望輸出與實際輸出對比圖22控制量u(k)圖23 參數(shù)辨識附錄程序2 帶有輔助變量的最小方差控制器：clear all;close all;clca=1 -1.8 0.7 0.1;b=1 1.8;c=1 0.2;d=4;na=length(a)-1; nb=length(b)-1;

17、nc=length(c)-1;nf=nb+d-1; ng=na-1; Pw=1;R=1;Q=2;np=length(Pw)-1; nh=length(R)-1+nc; nq=length(Q)-1;nm=np+ng;nG=max(np+nf,nc+nq);L=400;uk=zeros(d+nG,1);yk=zeros(d+nm,1);yrk=zeros(d+nh,1);xmk=zeros(nc,1);yr=10*ones(L/4,1);-ones(L/4,1);ones(L/4,1);-ones(L/4+d,1);x1=1;x2=1;x3=1;x4=0;S=1;for k=1:L IM=xor

18、(S,x4); if IM=0 xm(k)=-0.1; else xm(k)=0.1; end S=not(S); M(k)=xor(x3,x4); x4=x3;x3=x2;x2=x1;x1=M(k);endthetae_1=zeros(nm+nG+nh+3,d);P=106*eye(nm+nG+nh+3);lambda=0.998;for k=1:L time(k)=k; y(k)=-a(2:na+1)*yk(1:na)+b*uk(d:d+nb)+c*xm(k);xmk; phi=yk(d:d+nm);uk(d:d+nG);-yrk(d:nh+d); z(k)=Pw*y(k)-R*yrk(d

19、)+Q/b(1)*Q*uk(d); K=P*phi/(lambda+phi'*P*phi); thetae(:,k)=thetae_1(:,1)+K*(z(k)-phi'*thetae_1(:,1); P=(eye(nm+nG+nh+3)-K*phi')*P/lambda; me=thetae(1:nm+1,k)'ge=thetae(nm+2:nm+nG+2,k)'he=thetae(nm+nG+3:nm+nG+3+nh,k)' if ge(1)<0.1 ge(1)=0.1; end u(k)=(he(1)*yr(k)+he(2:nh+1)

20、*yrk(1:nh)-me(1)*y(k)-me(2:nm+1)*yk(1:nm)-ge(2:nG+1)*uk(1:nG)/ge(1); thetae; for i=d:-1:2 thetae_1(i)=thetae_1(i-1); end thetae_1(:,1)=thetae(:,k); for i=nG+d:-1:2 uk(i)=uk(i-1); end uk(1)=u(k); for i=nm+d:-1:2 yk(i)=yk(i-1); end yk(1)=y(k); for i=d+nh:-1:2 yrk(i)=yrk(i-1); end yrk(1)=yr(k); for i=n

21、c:-1:2 xmk(i)=xmk(i-1); end xmk(1)=xm(k);endfigure(1);plot(time,yr(1:L),'r:',time,y);xlabel('k');ylabel('y_r(k)¡¢y(k)');legend('y_r(k)','y(k)');axis(0 L -20 20); figure(2)plot(time,u);xlabel('k');ylabel('u(k)');axis(0 L -10 10); figur

22、e(3)plot(1:L,thetae(1:nm+1,:),1:L,thetae(nm+2:nm+nG+2,:),1:L,thetae(nm+nG+3:nm+nG+3+nh,:);xlabel('k');ylabel('參數(shù)估計m,G,h');legend('m_0','m_1','m_2','G_0','G_1','G_2','G_3','G_4','h_0','h_1')axis(0 L -3 10)

23、;3 參考模型自適應(yīng)控制3.1參考模型自適應(yīng)控制模型參考自適應(yīng)控制器（MRAC，model reference adaptive control），即為利用可調(diào)系統(tǒng)（包含被控對象）的各種信息，度量或測出某種性能指標(biāo)，把它與參考模型期望的性能指標(biāo)相比較；用性能指標(biāo)偏差（廣義誤差）通過非線性反饋的自適應(yīng)律來調(diào)節(jié)可調(diào)系統(tǒng)，以削弱可調(diào)系統(tǒng)因“不確定性”所造成的性能指標(biāo)的偏差，最后達到使被控的可調(diào)系統(tǒng)獲得較好的性能指標(biāo)的目的。MARC中的可調(diào)系統(tǒng)一般包括被控對象和調(diào)節(jié)器，它們形成一常規(guī)的反饋控制系統(tǒng)。 這個系統(tǒng)相對于MARC 系統(tǒng)來說是一個子系統(tǒng)或稱“內(nèi)回路”。另外， MARC系統(tǒng)還有一個自適應(yīng)反饋回路

24、，稱為外回路，它用來調(diào)節(jié)可調(diào)系統(tǒng)。有內(nèi)外回路組成雙回路系統(tǒng)是MARC 系統(tǒng)的結(jié)構(gòu)特點5。模型參考自適應(yīng)控制可以處理緩慢變化的不確定性對象的控制問題。 它由于可以不必經(jīng)過系統(tǒng)辨識而度量性能指標(biāo)，因而有可能或得快速跟蹤控制。由于被控對象的全部狀態(tài)要準(zhǔn)確得到很困難， 按被控對象輸入和輸出直接設(shè)計自適應(yīng)控制系統(tǒng)更有價值，一般有直接法和間接法。所謂直接和間接，指的是對未知的被控對象進行直接控制和間接控制。 間接控制的基本思想是用未知的被控對象的輸入輸出數(shù)據(jù)來估計被控對象的參數(shù)， 并用這些參數(shù)估計值產(chǎn)生一個反饋函數(shù)去調(diào)整調(diào)節(jié)器參數(shù)。直接控制和間接控制不同，在產(chǎn)生反饋控制信號之前沒有明顯的被控對象的辨識。所

25、以二者之間主要不同在于：在直接控制中要有一個顯式的理想特性的參考模型， 而間接控制則需要被控對象模型進行在線辨識并用隱式方法去產(chǎn)生自適應(yīng)律。此外，間接控制是用辨識誤差，而直接控制則用控制誤差去修改調(diào)節(jié)器參數(shù)6?，F(xiàn)用一種直接法，即 K.S.Narendra 提出的穩(wěn)定性自適應(yīng)控制方案。原理框圖如圖3-1所示：圖3-1 模型參考自適應(yīng)控制原理框圖由于系統(tǒng)的對象的分母階次比分子大兩階，自適應(yīng)控制器可設(shè)計如下。(1)選擇，并構(gòu)造輔助信號狀態(tài)方程(2)可調(diào)參數(shù)自適應(yīng)律為：式中，為正定矩陣。(3)自適應(yīng)控制律為：參考模型自適應(yīng)控制的算法：已知：被控對象的階數(shù)、。步驟1：選擇參考模型為嚴(yán)格正實、穩(wěn)定最小相位

26、系統(tǒng)，與階數(shù)及相對階相同，并具有理想的動態(tài)性能；并利用構(gòu)造輔助信號發(fā)生器狀態(tài)矩陣；步驟2：設(shè)置初值，選擇自適應(yīng)增益矩陣和輸入信號，并初始化數(shù)據(jù)；步驟3：采樣當(dāng)前參考模型輸出和系統(tǒng)實際輸出，并計算；步驟4：輔助信號狀態(tài)方程計算和；步驟5：利用自適應(yīng)律計算；步驟6：組件，并有自適應(yīng)控制律計算；步驟7：，返回步驟3，繼續(xù)循環(huán)。3.2仿真實例3.2.1數(shù)值積分對于如下連續(xù)系統(tǒng)的一階微分方程：在離散系統(tǒng)中可采用歐拉法求解，其遞推公式如下式中，稱為計算步長或者步距。該方法簡單、計算量小，由前一點即可推出后一點的值，屬于單步法。適當(dāng)減小計算步長h有助于提高計算精度。3.2.2仿真結(jié)果采用的對象模型為模型靜態(tài)

27、增益為1，故采用的參考模型為改模型是嚴(yán)格正實的最小相位系統(tǒng)，并且具有良好的動態(tài)性能。取自適應(yīng)增益矩陣，輸入信號為方波信號，使用模型參考自適應(yīng)控制后的仿真結(jié)果如下：圖 3-1 實際輸出與參考輸出對比(左圖無噪聲，右圖有噪聲) 圖 3-2 偏差值e(t) 圖 3-3 控制量u(t)圖 3-4 可調(diào)參數(shù)kc、c1、c2、c3曲線 圖 3-5 可調(diào)參數(shù)d0、d1、d2、d3曲線附錄程序3 參考模型自適應(yīng)控制：clear all; close all;h=0.1; L=800/h; F=1 2;nump=1 4 4; denp=1 8 22 24 9; Ap,Bp,Cp,Dp=tf2ss(nump,de

28、np); n=length(denp)-1; %對象參數(shù) xi=sqrt(0.01)*randn(L,1); %白噪聲 numm=3.375 10.125 6.75; denm=1 8 21.75 22.5 6.75; Am,Bm,Cm,Dm=tf2ss(numm,denm); Df=conv(F,numm); %輔助信號發(fā)生器傳遞函數(shù)分母Af=zeros(n-2,1),eye(n-2);-Df(n:-1:2); %輔助信號發(fā)生器狀態(tài)矩陣Bf=zeros(n-2,1);1; %輔助信號發(fā)生器輸入矩陣 yr0=0; yp0=0; u0=0; e0=0; %初值v10=zeros(n-1,1);

29、v20=zeros(n-1,1); %輔助信號發(fā)生器狀態(tài)初值xp0=zeros(n,1); xm0=zeros(n,1); %狀態(tài)向量初值theta0=zeros(2*n,1); %可調(diào)參數(shù)向量初值zeta0=zeros(2*n,1); yr=ones(1,L/80) -ones(1,L/80);for i=1:39 yr=yr ones(1,L/80) -ones(1,L/80);end Gamma=diag(0.8 1 0.2 0.2 0.3 0.1 0.2 10); %自適應(yīng)增益矩陣 for k=1:L time(k)=k*h; xp(:,k)=xp0+h*(Ap*xp0+Bp*u0); yp(k)=Cp*xp(:,k)+Dp*u0; yp(k)=yp(k)+xi(k); %加入白噪聲 xm(:,k)=xm0+h*(Am*xm0+Bm*yr0); ym(k)=Cm*xm(:,k)+D