自動(dòng)化級(jí)系統(tǒng)辨識(shí)大作業(yè)王萬

1、系統(tǒng)辨識(shí)大作業(yè)班 級(jí)： 自動(dòng)化1101班 姓 名： 王萬秋 學(xué) 號(hào)： 11052204 第一題模仿index3，搭建如下的單輸入單輸出系統(tǒng)的差分方程取真值、和，輸入信號(hào)采用4階m序列，幅值為1。當(dāng)?shù)木禐?，方差分別為0.1和0.5的高斯噪聲時(shí)，分別用一般最小二乘法、遞推最小二乘法和增廣遞推最小二乘法估計(jì)參數(shù)。并通過對(duì)三種方法的辨識(shí)結(jié)果的分析和比較，說明上述三種參數(shù)辨識(shí)方法的優(yōu)缺點(diǎn)。（15分）利用simulink搭建的模型框圖如下：一般最小二乘法程序：u=uy(1:450,1)' %輸入矩陣z=uy(1:450,2)' %輸出矩陣h=zeros(400,4); for i=1:

2、400 h(i,1)=-z(i+1); h(i,2)=-z(i); h(i,3)=u(i+1); h(i,4)=u(i); end theta=inv(h'*h)*h'*(z(3:402)'遞推最小二乘法程序代碼：u=uy(1:450,1)' %輸入矩陣z=uy(1:450,2)' %輸出矩陣p=100*eye(4); %估計(jì)方差pstore=zeros(4,401); pstore(:,1)=p(1,1),p(2,2),p(3,3),p(4,4); theta=zeros(4,401); %參數(shù)的估計(jì)值，存放中間過程估值 theta(:,1)=3;3;

3、3;3; k=10;10;10;10; for i=3:402 h=-z(i-1);-z(i-2);u(i-1);u(i-2); k=p*h*inv(h'*p*h+1); theta(:,i-1)=theta(:,i-2)+k*(z(i)-h'*theta(:,i-2); p=(eye(4)-k*h')*p; pstore(:,i-1)=p(1,1),p(2,2),p(3,3),p(4,4); end i=1:401; theta = theta(:,length(theta(1,:)subplot(2,1,1)plot(i,theta(1,:),i,theta(2,:

4、),i,theta(3,:),i,theta(4,:) title('遞推最小二乘的估計(jì)值過度情況') legend('a1','a2','b1','b2')subplot(2,1,2)plot(i,pstore(1,:),i,pstore(2,:),i,pstore(3,:),i,pstore(4,:) title('遞推最小二乘估計(jì)方差過度情況') legend('a1','a2','b1','b2')增廣遞推最小二乘法程序代碼：u

5、=uy(1:450,1)' %輸入矩陣z=uy(1:450,2)' %輸出矩陣v=uy(1:450,3)' %噪聲矩陣p=100*eye(6); %估計(jì)方差 pstore=zeros(6,300); pstore(:,1)=p(1,1),p(2,2),p(3,3),p(4,4),p(5,5),p(6,6); theta=zeros(6,300); %參數(shù)的估計(jì)值，存放中間過程估值 theta(:,1)=3;3;3;3;3;3; k=10;10;10;10;10;10; for i=3:300 h=-z(i-1);-z(i-2);u(i-1);u(i-2);v(i-1);

6、v(i-2); k=p*h*inv(h'*p*h+1); theta(:,i-1)=theta(:,i-2)+k*(z(i)-h'*theta(:,i-2); p=(eye(6)-k*h')*p; pstore(:,i-1)=p(1,1),p(2,2),p(3,3),p(4,4),p(5,5),p(6,6); end i=1:300; theta=theta(:,length(theta(1,:)-10)subplot(2,1,1)plot(i,theta(1,:),i,theta(2,:),i,theta(3,:),i,theta(4,:),i,theta(5,:),

7、i,theta(6,:) title('增廣遞推最小二乘估計(jì)值過渡情況') legend('a1','a2','b1','b2')subplot(2,1,2)plot(i,pstore(1,:),i,pstore(2,:),i,pstore(3,:),i,pstore(4,:),i,pstore(5,:),i,pstore(6,:) title('增廣遞推最小二乘估計(jì)方差過度情況') legend('a1','a2','b1','b2'

8、)辨識(shí)結(jié)果：噪聲辨識(shí)方法一般最小二乘法遞推最小二乘法增廣遞推最小二乘法噪聲方差為0.1theta = 1.1141 0.5636 0.9813 0.2415theta = 1.1145 0.5640 0.9813 0.2420theta = 1.4021 0.52540.9853 0.2438噪聲方差為0.5theta = 0.8810 0.4127 0.9587 0.0558theta = 0.8812 0.4128 0.9587 0.0561theta = 1.41820.56240.9711 0.2647噪聲方差為0.1噪聲方差為0.5由圖表可得：一般最小二乘法和遞推最小二乘法的辨識(shí)結(jié)果

9、很接近，而增廣遞推最小二乘法的辨識(shí)結(jié)果和其他兩種方法差距較大，尤其a1和b2。此外，噪聲越大，對(duì)前兩種方法辨識(shí)效果影響較大，對(duì)于增廣遞推最小二乘法，抗噪聲能力較強(qiáng)。三種方法的優(yōu)缺點(diǎn)：一般最小二乘法優(yōu)點(diǎn)：白噪聲可得無偏漸進(jìn)無偏估計(jì)；算法簡(jiǎn)單、可靠；計(jì)算量最??；一次即可完成算法，適合于離線辨識(shí)。缺點(diǎn)：隨著數(shù)據(jù)的增長(zhǎng)，該方法將出現(xiàn)“數(shù)據(jù)飽和”現(xiàn)象。遞推最小二乘法優(yōu)點(diǎn)：可以減少數(shù)據(jù)存儲(chǔ)量，避免矩陣求逆，兼騷運(yùn)算量。缺點(diǎn)：為避免出現(xiàn)數(shù)據(jù)飽和現(xiàn)象，必須限制過去數(shù)據(jù)的影響。增廣遞推最小二乘法優(yōu)點(diǎn)：比遞推算法有著更高的準(zhǔn)確度。增廣最小二乘的遞推算法對(duì)應(yīng)的噪聲模型為滑動(dòng)平均噪聲，擴(kuò)充了參數(shù)向量和數(shù)據(jù)向量的維數(shù)，

10、把噪聲模型的辨識(shí)同時(shí)考慮進(jìn)去。以上兩種最小二乘法只能獲得過程模型的參數(shù)估計(jì)，而增廣最小二乘法同時(shí)又能獲得噪聲模型的參數(shù)估計(jì)。當(dāng)數(shù)據(jù)長(zhǎng)度較大時(shí)，辨識(shí)精度低于極大似然法。第二題模仿實(shí)驗(yàn)二，搭建對(duì)象，由相關(guān)分析法，獲得脈沖響應(yīng)序列，由，參照講義， 獲得系統(tǒng)的脈沖傳遞函數(shù)和傳遞函數(shù);應(yīng)用最小二乘辨識(shí),獲得脈沖響應(yīng)序列;應(yīng)用相關(guān)最小二乘法，擬合對(duì)象的差分方程模型；加階躍擾動(dòng)，用最小二乘法和帶遺忘因子的最小二乘法，(可以用辨識(shí)工具箱) 辨識(shí)二階差分方程模型的參數(shù)，比較兩種方法的辨識(shí)效果差異；采用不少于兩種定階方法，確定模型的階次。（40分）搭建模型如下：1. 相關(guān)分析法獲得脈沖響應(yīng)序列程序代碼：np =

11、24-1;deta = 2;a = 5; u = zeros(np,np);for i=1:1:np for j=1:1:np u(i,j) = uy(np+1-i+j,1); endendy = zeros(np,1);for i=1:1:np y(i,1)=uy(np+i,2);endr = ones(np,np);for i = 1:1:np r(i,i) = 2;endg = (1/(a*a*(np+1)*deta)*r*u*yt = 0:2:28plot(t',g,'-*')獲得的脈沖響應(yīng)序列如下：g = 0.0127 0.2105 0.3058 0.2690

12、 0.1888 0.1330 0.0796 0.0603 0.0595 0.0419 0.0236 0.0228 0.0146 0.0337 0.0069響應(yīng)序列的散點(diǎn)圖：2. 求解g(z)和g(s)求解g(z)程序代碼：clcn=3;ts=2;np=15;u=uy(31:200,1);y=uy(31:200,2);f=zeros(120,6);f(:,1)=-1*y(31:150);f(:,2)=-1*y(30:149);f(:,3)=-1*y(29:148);f(:,4)=-1*u(31:150);f(:,5)=-1*u(30:149);f(:,6)=-1*u(29:148);yy=y(3

13、2:151);q=inv(f'*f)*f'*yy;a=q(1:3);b=q(4:6);a=a',b=b'g=tf(b,1,a,2)則系統(tǒng)的離散傳遞函數(shù)：transfer function: -0.3946 z2 - 0.3493 z - 0.07706-z3 - 0.7322 z2 + 0.02371 z + 0.01708sampling time: 2因此，將離散傳函轉(zhuǎn)化為連續(xù)傳函：gs = d2c(g,'zoh')則有：transfer function: 0.1374 s3 - 0.4647 s2 - 0.3989 s - 1.571-s

14、4 + 3.063 s3 + 5.763 s2 + 3.892 s + 0.5904我覺得該連續(xù)傳函不合理（階次有點(diǎn)高）。3. 應(yīng)用最小二乘法獲得脈沖響應(yīng)序列yy=uy(31:54,2)/2;for i=24:-1:1 uu(i,:)=uy(i+30:-1:16+i,1)'endgg=inv(uu'*uu)*uu'*yy;plot(0:2:29,gg,'b')hold onstem(0:2:29,gg,'b','filled')title('最小二乘法獲得脈沖序列');獲得的脈沖響應(yīng)序列為：gg = 0.0

15、588 0.2617 0.3784 0.3243 0.2403 0.1786 0.1439 0.1048 0.0834 0.0811 0.0753 0.0853 0.0799 0.06830.0587響應(yīng)序列的散點(diǎn)圖：4. 相關(guān)最小二乘法擬合對(duì)象的差分方程u=uy(1:500,1)' %輸入矩陣z=uy(1:500,2)' %輸出矩陣h=zeros(400,4); for i=1:400 h(i,1)=-z(i+1); h(i,2)=-z(i); h(i,3)=u(i+1); h(i,4)=u(i); end theta=inv(h'*h)*h'*(z(3:40

16、2)'計(jì)算結(jié)果：theta = -0.9631 0.2145 0.38950.25925. 加擾動(dòng)后，用最小二乘法和帶遺忘因子的最小二乘法辨識(shí)差分方程（并比較兩種方法）一般最小二乘法辨識(shí)結(jié)果：theta = -1.2184 0.2441 0.4471 0.1591帶遺忘因子的最小二乘法程序代碼：na=2; nb=2;d=0; l=400;u=0.95; uk=uy(1:400,1)' yk=uy(1:400,2)' thetae_1=zeros(na+nb,1); p=106*eye(na+nb); for k=3:lphi=-yk(k-1):-1:(k-2) uk(k

17、-1):-1:(k-2); k=p*phi'/(u+phi*p*phi'); thetae(:,k)=thetae_1+k*(yk(k)-phi*thetae_1); p=(eye(na+nb)-k*phi)*p/u; thetae_1=thetae(:,k); for i=d+nb:-1:2 uk(i)=uk(i-1); end for i=na:-1:2 yk(i)=yk(i-1); endendplot(1:l,thetae);theta = thetae(:,400)legend('a1','a2','b1','b

18、2'); title('帶遺忘因子的最小二乘法便是參數(shù)')辨識(shí)結(jié)果：theta = -1.0907 0.1004 0.4598 0.2370 比較兩種辨識(shí)方法結(jié)果可得：a1和b1比較一致，而a2和b2差別較大，帶遺忘因子辨識(shí)過程波動(dòng)較大。6. 兩種定階方法確定模型的階次第一種方法aic定階方法：程序代碼：l = 200;u = uy(:,1);y = uy(:,2);for n = 1:1:5 n = 100-n; yd(1:n,1) = y(n+1:n+n); for i=1:n for l=1:1:2*n if(l<=n) fia(i,l) = -y(n+i-

19、l); else fia(i,l) = u(2*n+i-l); end end end thita = inv(fia'*fia)*fia'*yd; omiga = (yd-fia*thita)'*(yd-fia*thita)/n; aic(n) = n*log(omiga)+4*n;endplot(aic,'-');title('aic隨階次的變化曲線');xlabel('n');ylabel('aic');aic隨階次變化曲線：因此系統(tǒng)階次為2.第二種方法殘差定階方法：程序代碼：data = uy; l

20、 = length(data);%輸入輸出數(shù)據(jù)長(zhǎng)度jn = zeros(1,10);t = zeros(1,10);rm = zeros(10,10);for n=1:1:10; n = l-n; fia = zeros(n,2*n);%構(gòu)造測(cè)量矩陣 du = zeros(n,1); dy = zeros(n+1,1); r1 = 0;r0 = 0;for i = 1:n %測(cè)量矩陣賦值 for l = 1:n*2 if(l<=n) fia(i,l) = -data(n+i-l,2); elseif(n+i-l+2>0) fia(i,l) = data(n+i-l+2,1); en

21、d endendy = data(n+1:n+n,2);%輸出數(shù)據(jù)矩陣thita = inv(fia'*fia)*fia'*y;%計(jì)算參數(shù)矩陣 jn0 = 0; for k = n+1:n+n for j = 1:n du(j) = data(k-j,1); dy(j) = data(k+1-j,2); end dy(n+1) = data(1,2); e1(k) = 1,thita(1:n)'*dy-thita(n+1:2*n)'*du; jn1 = jn0+e1(k)2; %f檢驗(yàn)法 t(n) = abs(jn0-jn1)/jn1*(n-2*n-2)/2);

22、 jn0 = jn1; end jn(n) = jn0; for m = 0:1:9 for m2 = n+1:1:l-m r1 = r1+e1(m2)*e1(m2+m)/(l-m-n); r0 = r0+e1(m2)2/(l-m-n); end rm(n,m+1) = r1/r0; endendsubplot(2,1,1);plot(1:10,jn,'r');title('殘差平方和jn曲線圖');subplot(2,1,2);plot(1:1:10,t,'g');title('f檢驗(yàn)法結(jié)果圖');figure(2);plot(

23、0:9,rm(1,:),'g'),hold on;plot(0:9,rm(2,:),'b'),hold on;plot(0:9,rm(3,:),'r');title('殘差白色定階結(jié)果圖');hold off;定階結(jié)果：由殘差平方和jn曲線（殘差平方和越小，階次越合理）可定系統(tǒng)階次為2。綜上，系統(tǒng)階次定為2，也符合實(shí)際二階水箱系統(tǒng)的階次。第三題模仿實(shí)驗(yàn)三，搭建適合廣義最小二乘法的對(duì)象模型，實(shí)現(xiàn)廣義最小二乘法;把學(xué)過的最小二乘算法應(yīng)用于這個(gè)對(duì)象,比較辨識(shí)的效果差異（20分）搭建的模型：廣義最小二乘法程序代碼：u=uy(1:700,1

24、)' %輸入矩陣z=uy(1:700,2)' %輸出矩陣e=uy(1:700,3)'p=100*eye(4); %估計(jì)方差 pstore=zeros(4,400); pstore(:,2)=p(1,1),p(2,2),p(3,3),p(4,4); theta=zeros(4,400); %參數(shù)的估計(jì)值，存放中間過程估值 theta(:,2)=3;3;3;3; k=10;10;10;10; %增益 pe=10*eye(2); thetae=zeros(2,400); thetae(:,2)=0.5;0.3; ke=10;10; for i=3:400 h=-z(i-1);

25、-z(i-2);u(i-1);u(i-2); k=p*h*inv(h'*p*h+1); theta(:,i)=theta(:,i-1)+k*(z(i)-h'*theta(:,i-1); p=(eye(4)-k*h')*p; pstore(:,i-1)=p(1,1),p(2,2),p(3,3),p(4,4); he=-e(i-1);-e(i-2); ke=pe*he*inv(1+he'*pe*he); thetae(:,i)=thetae(:,i-1)+ke*(e(i)-he'*thetae(:,i-1); pe=(eye(2)-ke*he')*p

26、e; end disp('參數(shù)a1 a2 b1 b2的估計(jì)結(jié)果：') theta(:,400) disp('噪聲傳遞系數(shù)c1 c2的估計(jì)結(jié)果：') thetae(:,400) i=1:400; figure(1) plot(i,theta(1,:),i,theta(2,:),i,theta(3,:),i,theta(4,:) title('參數(shù)過渡過程') legend('a1','a2','b1','b2')figure(2) plot(i,pstore(1,:),i,pstore

27、(2,:),i,pstore(3,:),i,pstore(4,:) title('估計(jì)方差過渡過程') figure(3) plot(i,thetae(1,:),i,thetae(2,:); title('噪聲傳遞系數(shù)過渡過程') legend('e1','e2')運(yùn)行結(jié)果：參數(shù)的估計(jì)：theta = 1.3580 0.7768 0.9932 0.4771噪聲傳遞系數(shù)：thetae = -0.0089 0.0823各參數(shù)變化曲線：第四題對(duì)第三題的對(duì)象，用夏氏法和輔助變量法實(shí)現(xiàn)（15分）一、夏氏法：程序代碼：data = uy;%生

28、成數(shù)據(jù)矩陣%使用夏氏修正法，對(duì)2階系統(tǒng)進(jìn)行參數(shù)辨識(shí)n = 2;l = length(data);n = l-n;u = data(:,1);y = data(:,2);glol =-y(2:l-1),-y(1:l-2),u(2:l-1),u(1:l-2);zgl1 = data(3:l,2);sgl1 = glol'*glol;sgl2=inv(sgl1);sgl3=glol'*zgl1;xsthita = sgl2*sgl3;%計(jì)算參數(shù)矩陣thitab0 = 0;%設(shè)偏差項(xiàng)的偏差初值為0fa = sgl2*glol'm = eye(n)-glol*sgl2*glol&

29、#39;f = eye(n);if(f>=10-6*eye(n) e1 = zgl1-glol*xsthita;%計(jì)算殘差e omiga(2:n,1) = -e1(1:n-1); omiga(3:n,2) = -e1(1:n-2); d = omiga'*m*omiga; fx = inv(d)*omiga'*m*zgl1; thitab = fa*omiga*fx; xsthita = xsthita - thitab; f = thitab - thitab0; thitab0 = thitab;endtheta = xsthita(1) xsthita(2) xst

30、hita(3) xsthita(4)'結(jié)果：theta = 1.3587 0.7629 1.0017 0.4596二、輔助變量法：程序代碼：u=uy(1:700,1)' %輸入矩陣z=uy(1:700,2)' %輸出矩陣%遞推求解 p=100*eye(4); %估計(jì)方差 pstore=zeros(4,400); pstore(:,1)=p(1,1),p(2,2),p(3,3),p(4,4); theta=zeros(4,400); %參數(shù)的估計(jì)值，存放中間過程估值 theta(:,1)=3;3;3;3; theta(:,2)=3;3;3;3; theta(:,3)=0;

31、0;0;0; theta(:,4)=0;0;0;0; % k=zeros(4,400); %增益矩陣 k = 10;10;10;10; for i=5:400 h=-z(i-1);-z(i-2);u(i-1);u(i-2); hstar=-z(i-2-1);-z(i-2-2);u(i-1);u(i-2); %輔助變量 k=p*hstar*inv(h'*p*hstar+1); theta(:,i)=theta(:,i-1)+k*(z(i)-h'*theta(:,i-1); p=(eye(4)-k*h')*p; pstore(:,i-1)=p(1,1),p(2,2),p(3

32、,3),p(4,4); endtheta = theta(:,400) i=1:400; figure(1) plot(i,theta(1,:),i,theta(2,:),i,theta(3,:),i,theta(4,:) title('辨識(shí)參數(shù)過渡過程') figure(2) plot(i,pstore(1,:),i,pstore(2,:),i,pstore(3,:),i,pstore(4,:) title('估計(jì)方差過渡過程') 辨識(shí)結(jié)果：theta = 1.3813 0.8086 0.9787 0.4995第五題采用多變量系統(tǒng)的最小二乘辨識(shí)方法辨識(shí)如下mi

33、mo系統(tǒng)的參數(shù) 式中，和是同分布的隨機(jī)噪聲，且有；輸入信號(hào)采用4階序列，其幅值為1。模型的理想系數(shù)為（10分）不利用simulink中的模型對(duì)系統(tǒng)進(jìn)行仿真和測(cè)試，利用程序根據(jù)系統(tǒng)的差分方程組直接獲得系統(tǒng)的測(cè)試數(shù)據(jù)，并對(duì)系統(tǒng)參數(shù)進(jìn)行辨識(shí)。程序代碼：l=300; % m序列的長(zhǎng)度y1=1;y2=1;y3=1;y4=1;%4個(gè)移位積存器的輸出初始值for i=1:l; %開始循環(huán)，長(zhǎng)度為l x1=xor(y3,y4); x2=y1; x3=y2; x4=y3; y(i)=y4; if y(i)>0.5, u1(i)=-1,u2(i)=-1; %如果m序列的值為"1"時(shí),%辨識(shí)的輸入信號(hào)取“-1” else u1(i)=1,u2(i)=1; %當(dāng)m序列的值為"0"時(shí).辨識(shí)的輸入信號(hào)取“1” end y1=x1;y2=x2;y3=x3;y4=x4; %為下一次的輸入信號(hào)做準(zhǔn)備end %大循環(huán)結(jié)束，產(chǎn)生輸入信號(hào)u%高斯白噪聲，長(zhǎng)度為lv1=0.1*randn