自动化级系统辨识大作业 王万.doc

系统辨识大作业班 级： 自动化1101班 姓 名： 王万秋 学 号： 11052204 第一题模仿index3，搭建如下的单输入单输出系统的差分方程取真值、和，输入信号采用4阶M序列，幅值为1。当的均值为0，方差分别为0.1和0.5的高斯噪声时，分别用一般最小二乘法、递推最小二乘法和增广递推最小二乘法估计参数。并通过对三种方法的辨识结果的分析和比较，说明上述三种参数辨识方法的优缺点。（15分）利用simulink搭建的模型框图如下：一般最小二乘法程序：u=UY(1:450,1)' %输入矩阵z=UY(1:450,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:402)'递推最小二乘法程序代码：u=UY(1:450,1)' %输入矩阵z=UY(1:450,2)' %输出矩阵P=100*eye(4); %估计方差Pstore=zeros(4,401); Pstore(:,1)=P(1,1),P(2,2),P(3,3),P(4,4); Theta=zeros(4,401); %参数的估计值，存放中间过程估值 Theta(:,1)=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,:),i,Theta(3,:),i,Theta(4,:) title('递推最小二乘的估计值过度情况') legend('a1','a2','b1','b2')subplot(2,1,2)plot(i,Pstore(1,:),i,Pstore(2,:),i,Pstore(3,:),i,Pstore(4,:) title('递推最小二乘估计方差过度情况') legend('a1','a2','b1','b2')增广递推最小二乘法程序代码：u=UY(1:450,1)' %输入矩阵z=UY(1:450,2)' %输出矩阵v=UY(1:450,3)' %噪声矩阵P=100*eye(6); %估计方差 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); %参数的估计值，存放中间过程估值 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);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,:),i,Theta(6,:) title('增广递推最小二乘估计值过渡情况') 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('增广递推最小二乘估计方差过度情况') legend('a1','a2','b1','b2')辨识结果：噪声辨识方法一般最小二乘法递推最小二乘法增广递推最小二乘法噪声方差为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由图表可得：一般最小二乘法和递推最小二乘法的辨识结果很接近，而增广递推最小二乘法的辨识结果和其他两种方法差距较大，尤其a1和b2。此外，噪声越大，对前两种方法辨识效果影响较大，对于增广递推最小二乘法，抗噪声能力较强。三种方法的优缺点：一般最小二乘法优点：白噪声可得无偏渐进无偏估计；算法简单、可靠；计算量最小；一次即可完成算法，适合于离线辨识。缺点：随着数据的增长，该方法将出现“数据饱和”现象。递推最小二乘法优点：可以减少数据存储量，避免矩阵求逆，兼骚运算量。缺点：为避免出现数据饱和现象，必须限制过去数据的影响。增广递推最小二乘法优点：比递推算法有着更高的准确度。增广最小二乘的递推算法对应的噪声模型为滑动平均噪声，扩充了参数向量和数据向量的维数，把噪声模型的辨识同时考虑进去。以上两种最小二乘法只能获得过程模型的参数估计，而增广最小二乘法同时又能获得噪声模型的参数估计。当数据长度较大时，辨识精度低于极大似然法。第二题模仿实验二，搭建对象，由相关分析法，获得脉冲响应序列，由，参照讲义， 获得系统的脉冲传递函数和传递函数;应用最小二乘辨识,获得脉冲响应序列;应用相关最小二乘法，拟合对象的差分方程模型；加阶跃扰动，用最小二乘法和带遗忘因子的最小二乘法，(可以用辨识工具箱) 辨识二阶差分方程模型的参数，比较两种方法的辨识效果差异；采用不少于两种定阶方法，确定模型的阶次。（40分）搭建模型如下：1. 相关分析法获得脉冲响应序列程序代码：Np = 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,'-*')获得的脉冲响应序列如下：G = 0.0127 0.2105 0.3058 0.2690 0.1888 0.1330 0.0796 0.0603 0.0595 0.0419 0.0236 0.0228 0.0146 0.0337 0.0069响应序列的散点图：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(32: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)则系统的离散传递函数：Transfer function: -0.3946 z2 - 0.3493 z - 0.07706-z3 - 0.7322 z2 + 0.02371 z + 0.01708Sampling time: 2因此，将离散传函转化为连续传函：Gs = d2c(G,'zoh')则有：Transfer function: 0.1374 s3 - 0.4647 s2 - 0.3989 s - 1.571-s4 + 3.063 s3 + 5.763 s2 + 3.892 s + 0.5904我觉得该连续传函不合理（阶次有点高）。3. 应用最小二乘法获得脉冲响应序列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('最小二乘法获得脉冲序列');获得的脉冲响应序列为：GG = 0.0588 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响应序列的散点图：4. 相关最小二乘法拟合对象的差分方程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:402)'计算结果：Theta = -0.9631 0.2145 0.38950.25925. 加扰动后，用最小二乘法和带遗忘因子的最小二乘法辨识差分方程（并比较两种方法）一般最小二乘法辨识结果：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-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','b2'); title('带遗忘因子的最小二乘法便是参数')辨识结果：theta = -1.0907 0.1004 0.4598 0.2370 比较两种辨识方法结果可得：a1和b1比较一致，而a2和b2差别较大，带遗忘因子辨识过程波动较大。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-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随阶次变化曲线：因此系统阶次为2.第二种方法残差定阶方法：程序代码：Data = UY; L = length(Data);%输入输出数据长度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);%构造测量矩阵 du = zeros(n,1); dy = zeros(n+1,1); r1 = 0;r0 = 0;for i = 1:N %测量矩阵赋值 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); end endendY = Data(n+1:n+N,2);%输出数据矩阵thita = inv(FIA'*FIA)*FIA'*Y;%计算参数矩阵 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检验法 t(n) = abs(Jn0-Jn1)/Jn1*(N-2*n-2)/2); 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检验法结果图');figure(2);plot(0:9,rm(1,:),'g'),hold on;plot(0:9,rm(2,:),'b'),hold on;plot(0:9,rm(3,:),'r');title('残差白色定阶结果图');hold off;定阶结果：由残差平方和Jn曲线（残差平方和越小，阶次越合理）可定系统阶次为2。综上，系统阶次定为2，也符合实际二阶水箱系统的阶次。第三题模仿实验三，搭建适合广义最小二乘法的对象模型，实现广义最小二乘法;把学过的最小二乘算法应用于这个对象,比较辨识的效果差异（20分）搭建的模型：广义最小二乘法程序代码：u=UY(1:700,1)' %输入矩阵z=UY(1:700,2)' %输出矩阵e=UY(1:700,3)'P=100*eye(4); %估计方差 Pstore=zeros(4,400); Pstore(:,2)=P(1,1),P(2,2),P(3,3),P(4,4); Theta=zeros(4,400); %参数的估计值，存放中间过程估值 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);-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')*PE; end disp('参数a1 a2 b1 b2的估计结果：') Theta(:,400) disp('噪声传递系数c1 c2的估计结果：') ThetaE(:,400) i=1:400; figure(1) plot(i,Theta(1,:),i,Theta(2,:),i,Theta(3,:),i,Theta(4,:) title('参数过渡过程') legend('a1','a2','b1','b2')figure(2) plot(i,Pstore(1,:),i,Pstore(2,:),i,Pstore(3,:),i,Pstore(4,:) title('估计方差过渡过程') figure(3) plot(i,ThetaE(1,:),i,ThetaE(2,:); title('噪声传递系数过渡过程') legend('e1','e2')运行结果：参数的估计：theta = 1.3580 0.7768 0.9932 0.4771噪声传递系数：thetae = -0.0089 0.0823各参数变化曲线：第四题对第三题的对象，用夏氏法和辅助变量法实现（15分）一、夏氏法：程序代码：Data = UY;%生成数据矩阵%使用夏氏修正法，对2阶系统进行参数辨识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;%计算参数矩阵thitab0 = 0;%设偏差项的偏差初值为0Fa = Sgl2*glOL'M = eye(N)-glOL*Sgl2*glOL'F = eye(N);if(F>=10-6*eye(N) E1 = Zgl1-glOL*Xsthita;%计算残差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) Xsthita(3) Xsthita(4)'结果：theta = 1.3587 0.7629 1.0017 0.4596二、辅助变量法：程序代码：u=UY(1:700,1)' %输入矩阵z=UY(1:700,2)' %输出矩阵%递推求解 P=100*eye(4); %估计方差 Pstore=zeros(4,400); Pstore(:,1)=P(1,1),P(2,2),P(3,3),P(4,4); Theta=zeros(4,400); %参数的估计值，存放中间过程估值 Theta(:,1)=3;3;3;3; Theta(:,2)=3;3;3;3; Theta(:,3)=0;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,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('辨识参数过渡过程') figure(2) plot(i,Pstore(1,:),i,Pstore(2,:),i,Pstore(3,:),i,Pstore(4,:) title('估计方差过渡过程') 辨识结果：theta = 1.3813 0.8086 0.9787 0.4995第五题采用多变量系统的最小二乘辨识方法辨识如下MIMO系统的参数 式中，和是同分布的随机噪声，且有；输入信号采用4阶序列，其幅值为1。模型的理想系数为（10分）不利用simulink中的模型对系统进行仿真和测试，利用程序根据系统的差分方程组直接获得系统的测试数据，并对系统参数进行辨识。程序代码：L=300; % M序列的长度y1=1;y2=1;y3=1;y4=1;%4个移位积存器的输出初始值for i=1:L; %开始循环，长度为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"时,%辨识的输入信号取“-1” else u1(i)=1,u2(i)=1; %当M序列的值为"0"时.辨识的输入信号取“1” end y1=x1;y2=x2;y3=x3;y4=x4; %为下一次的输入信号做准备end %大循环结束，产生输入信号u%高斯白噪声，长度为Lv1=0.1*randn(1,L);v2=0.1*randn(1,L);Y1(2)=0;Y1(1)=0; %取Y的前两个初始值为零Y2(2)=0;Y2(1)=0; %取Y的前两个初始值为零for k=3:L; Y1(k)=-0.5*Y1(k-1)+0.2*Y2(k-1)-1.2*Y1(k-2)+0.6*Y2(k-2)+u1(k)+0.5*u1(k-1)-0.4*u2(k-1)+0.4*u1(k-2)-0.3*u2(k-2)+v1(k);Y2(k)=0.3*Y1(k-1)-0.6*Y2(k-1)-0.1*Y1(k-2)+0.6*Y2(k-2)+u2(k)+0.2*u1(k-1)-0.3*u2(k-1)-0.2*u1(k-2)+0.1*u2(k-2)+v2(k);endN=150; n=2;A1=0.5 -0.2;-0.3 0.6;A2=1.2 -0.6;0.1 -0.6;B0=1.0 0.0;0.0 1.0;B1=0.5 -0.4;0.2 -0.3;B2=0.4 -0.3;-0.2 0.1;c0=zeros(5,1); %被辨识参数矩阵的初始值c0 = c0 + 3;3;3;3;3H=zeros(N,2*n+1);for j=1:2 for k=0:N-1 for i=2*n+1:-1:1 if i>n+1 H(k+1,2*n+2-i) = -Y(j,i-(n+1)+k)' else H(k+1,2*n+2-i) = -U(j,i+k)' end end end I = eye(N); h=H' cc = ; c = c0; P=inv(H'*H); for k=1:N P=P-P*h*inv(I+h'*P*h)*h'*P; K=P*h*inv(1+h'*P*h); c=c+K*(y(1:N)'-h'*c); cc = cc c; end a1=cc(1,:); a2=cc(2,:); b1=cc(3,:); b2=cc(4,:);b3=cc(5,:); theta(:,j) = cc(:,N); figure(j+1);i=1:N; plot(i,a1,'r',i,a2,':',i,b1,'g',i,b2,':') title('MIMO最小二乘各参数变化情况')endTheta辨识结果：theta = 0.0000 0.0000 0.0000 0.0000 0.0016 -0.0045 -0.0022 -0.0013 0.4998 0.4970各参数变化情况：