系统辨识大作业.docx

1、系统辨识大作业系统辨识大作业第一题x=1.01 2.03 3.02 4.01 5 6.02 7.03 8.04 9.03 10;y=9.6 4.1 1.3 0.4 0.05 0.1 0.7 1.7 3.8 9.1;ma,na=size(x);sumx=0;sumy=0;sumxy=0;sumx2=0;sumx2y=0;sumx3=0;sumx4=0;for k=1:na sumx=sumx+x(k); sumy=sumy+y(k) sumxy=sumxy+x(k)*y(k); sumx2=sumx2+x(k)*x(k); sumx2y=sumx2y+x(k)*x(k)*y(k); sumx3=

2、sumx3+x(k)3; sumx4=sumx4+x(k)4;endA=na,sumx,sumx2;sumx,sumx2,sumx3;sumx2,sumx3,sumx4;B=sumy;sumxy;sumx2y;C=AB；利用最小二乘的求得作图验证如下：第二题，搭建的模型如图（m序列发生器采用三个移位寄存器，其他参数同实验二）aa=5;NNPP=7;ts=2;RR=ones(7)+eye(7);UU=UY(15:21,1);UY(14:20,1);UY(13:19,1);UY(12:18,1);UY(11:17,1); UY(10:16,1);UY(9:15,1);YY=UY(8:14,2);G

3、G=(RR*UU*YY+4.4474)/aa*aa*(NNPP+1)*ts; plot(0:2:13,GG);hold onstem(0:2:13,GG,filled);grid on;title(脉冲响应序列)求系统的脉冲传递函数 g= -0.0036 0.2922 0.3683 0.2673 0.1705 0.0820 0.0489; a=;for k=1:1:6a=a;g(k),g(k+1);endG=a;g(7),g(1);G1=g(3:7,1);g(1);g(2);A=G(-G1);a1=A(2)a2=A(1)D=1 0;a1 1;C=g(1);g(2);B=D*C;b1=B(1)b

4、2=B(2)zuixiao 相关二步法z=UY(:,2);M=UY(:,1);P=100*eye(4); %Theta=zeros(4,196); %Theta(:,1)=3;3;3;3;Theta(:,2)=3;3;3;3;Theta(:,3)=3;3;3;3;Theta(:,4)=3;3;3;3;% K=zeros(4,400); %K=10;10;10;10;for i=5:196h=-z(i-1);-z(i-2);M(i-1);M(i-2);hstar=M(i-1);M(i-2);M(i-3);M(i-4); %K=P*hstar*inv(h*P*hstar+1);Theta(:,i)

5、=Theta(:,i-1)+K*(z(i)-h*Theta(:,i-1);P=(eye(4)-K*h)*P;end%=i=1:196;plot(i,Theta(1,:),i,Theta(2,:),i,Theta(3,:),i,Theta(4,:)title();grid on得到相关系数为ans = -0.7982 0.1377 0.5631 0.3149最小二乘递推算法相关参数clc;lamt=1;z=UY(:,2);u=UY(:,1);c0=0.001 0.001 0.001 0.001;p0=104*eye(4,4); c=c0,zeros(4,198); for k=3:200 h1=

6、-z(k-1),-z(k-2),u(k-1),u(k-2); x=h1*p0*h1+1*lamt; x1=inv(x); k1=p0*h1*x1; d1=z(k)-h1*c0; c1=c0+k1*d1; p1=1/lamt*(eye(4)-k1*h1)*p0; c(:,k-1)=c1; c0=c1; p0=p1;%C1enda1=c(1,:); a2=c(2,:); b1=c(3,:); b2=c(4,:); i=1:199;plot(i,a1,r,i,a2,b,i,b1,y,i,b2,black) ;title( );grid on得到相关系数为：ans = -0.7971 0.1423 0

7、.5616 0.3127带遗忘因子的最小二乘法将上述程序中的lamt=0.9；辨识系统阶次方法一、利用残差最小辨识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(l0) FIA(i,l) = Data(n+i-l+2,1); end endendY =

8、 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; t(n) = abs(Jn0-Jn1)/Jn1*(N-2*n-2)/2); Jn0 = Jn1; end Jn(n) = Jn0; for m = 0:1:9 for m2 =

9、 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(残差白色丁阶结果图);方法二

10、、AIC方法订阶次Data = UY;%AkkaikeL = length(Data);U = Data(:,1);Y = Data(:,2);for n = 1:1:5 N = 201-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(o

11、miga)+4*n;endplot(AIC,-);title(AI订阶法);xlabel(n);ylabel(AIC);grid on第三题，搭建的对象如图：广义最小二乘法n=800; %n M=UY(1:800,1);%=01 =v=randn(1,800);e=;e(1)=v(1);e(2)=v(2);for i=3:800e(i)=0*e(i-1)+0*e(i-2)+v(i);end%=z=z=UY(1:800,2);%zf=;zf(1)=-1;zf(2)=0;for i=3:800zf(i)=z(i)-0*z(i-1)-0*z(i-2);end%uf=;uf(1)=M(1);uf(2)

12、=M(2);for i=3:800uf(i)=M(i)-0*M(i-1)-0*M(i-2);end%P=100*eye(4); %Theta=zeros(4,800); %Theta(:,2)=3;3;3;3;K=10;10;10;10; %PE=10*eye(2);ThetaE=zeros(2,800);ThetaE(:,2)=0.5;0.3;KE=10;10;for i=3:800h=-zf(i-1);-zf(i-2);uf(i-1);uf(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;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%=i=1:800;figure(1);plot(i,Theta(1,:),i,Theta(2,:),i,Theta(3,:),i,Theta(4,:)title()figure(3);plot(i,ThetaE(1,:),