系统辨识作业2

系统辨识作业学院：专业：姓名：学号：日期：系统辨识作业：以下图为仿真对象图中，v(k)为服从N(0,1)正态分布的不相关随即噪声，输入信号采用循环周期Np>500的逆M 序列，幅值为1，选择辨识模型为：)()2()1()2()1()(2121k v k u b k u b k z a k z a k z +-+-=-+-+加权因子1)(=Λk ，数据长度L=500，初始条件取I P 610)0(= ，⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎣⎡=001.0001.0001.0)0(ˆ θ 要求：（1）采用一次完成最小二乘法对系统进行辨识，给出数据u(k)和z(k)，及L H，L Z 和θ 和)ˆ(θJ 的值。

（2）采用递推最小二乘法进行辨识，要给出参数收敛曲线以及新息)(~k Z，残差)(k ε，准则函数)(k J 随着递推次数K 的变化曲线。

（3）对仿真对象和辨识出的模型进行阶跃响应对比分析以检验辨识结果的实效。

1、一次完成法对系统进行辨识：估计LT LLT LLSZ H H H 1)(ˆ-=θ，其中 []2121,,,b b a a LS =θ⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎣⎡=L L Z Z Z Z 21⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎣⎡------------=⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎣⎡=)2()1()2()1()0()1()0()1()1()0()1()0()()2()1(L u L u L z L z u u z z u u z z L h h h H L 一次完成算法对系统辨识的Matlab 程序见附录:部分输入、输出数据如下，全部的输入输出数据用图1.1所示 输入数据u(k)=Columns 1 through 160 0 1 1 1 1 0 0 0 0 1 0 0 1 1 0输出数据z(k)= Columns 1 through 90 0 1.2372 3.9331 6.4987 7.9909 7.7132 6.5947 5.4523Columns 10 through 183.2212 0.8419 0.6243 1.7110 0.7126 1.0712 2.8037 3.80474.6372图1.1输入数据u(k)和输出数据z(k)I H的值为(部分)：HL=-5.4523 -6.5947 0 0 -3.2212 -5.4523 0 0 -0.8419 -3.2212 1.0000 0 … … … … -3.5706 -5.1944 0 0 -0.6787 -3.5706 0 02.3020 -0.6787 0 0 I Z的值为(部分)：ZL= 3.2212 0.8419 0.6243 1.71100.7126 … 0.6787 -2.3020 -3.8270一次完成辨识后的值为：theta= -1.4918 0.6932 1.0541 0.4691)ˆ(θ J 的值为：J(theta)=493.1837 2、递推最小二乘法辨识系统： )]1()(')()[()1()(--+-=k k h k z k K k k θθθ1])(1)()1()(')[()1()(-Λ+--=k k h k P k h k h k P k K )1()](')([)(--=k P k h k K I k P初始条件：I P 610)0(= ，⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎣⎡=001.0001.0001.0)0(ˆ θ 递推最小二乘法辨识系统的Matlab 程序见附录: 其参数收敛曲线如图2.1图2.1 参数收敛曲线新息)(~k Z ，残差)(k ε，准则函数)(k J 随着递推次数K 的变化曲线如图2.2中依次所示:图2.2 新息、残差、准则函数变化曲线3、仿真对象和辨识出的模型进行阶跃响应对比分析仿真对象和辨识模型阶跃响应对比Matlab程序见附录:图 3.1分别给出了一次完成算法辨识出来系统和辨识对象的阶跃响应对比图，递推算法辨识出来系统和辨识对象的阶跃响应对比图。

附录：一次完成和递推法系统辨识Matlab程序%%最小二乘法辨识系统;%叠加噪声为1/(1-1.5z^(-1)+0.7z^(-2));%化为差分方程形式为;%e(k)=v(k)+1.5e(k-1)-0.7e(k-2);%仿真对象为(1z^-1+0.5z^-2)/(1-1.5z^-1+0.7z^-2);%化为差分方程形式为;%y(k)-1.5y(k-1)+0.7y(k-2)=u(k-1)+0.5u(k-2);%辨识模型为;%z(k)=-a1z(k-1)-a2z(k-2)+b1u(k-1)+b2u(k-2)+v(k);%================================%主函数;function mainclose all;clc;clear;%================%产生逆M序列;X=[0,1,1,0,1,0,0,1]; %寄存器初始值;F=[0,1,1,1,0,0,0,1]; %特征多项式;N=1000; %生成长度;M=[]; %存放产生的M序列; %产生逆M序列函数调用;out=IMxulie(X,F,N);%阶梯图输出逆M序列;figure(1);M=out;subplot(2,1,1);stairs(M);xlabel('k')ylabel('逆M序列')title('移位寄存器产生的逆M序列');grid on;%=================%一次完成最小二乘法辨识系统;%e(k)=v(k)+1.5e(k-1)-0.7e(k-2);%y(k)=1.5y(k-1)-0.7y(k-2)+u(k-1)+0.5u(k-2);%z(k)=y(k)+e(k);%即z(k)=1.5y(k-1)-0.7y(k-2)+u(k-1)+0.5u(k-2)+v(k)+1.5e(k-1)-0.7e(k-2); %产生均值为0，方差为1的白噪声;v=[];v=randn(1,600);%产生系统辨识所需要数据;y(1)=0;y(2)=0;e(1)=0;e(2)=0;for i=3:600y(i)=1.5*y(i-1)-0.7*y(i-2)+M(i-1)+0.5*M(i-2);e(i)=v(i)+1.5*e(i-1)-0.7*e(i-2);z(i)=y(i)+e(i);endsubplot(2,1,2);stem(z);xlabel('k');ylabel('幅度值');title('输出数据');grid on;%辨识模型为;%z(k)=-a1z(k-1)-a2z(k-2)+b1u(k-1)+b2u(k-2)+v(k);for i=10:509H(i-9,:)=[-z(i-1),-z(i-2),M(i-1),M(i-2)];endZL=(z(10:509))';ES=inv(H'*H)*H'*ZL; %inv用来求矩阵的逆矩阵;%输出辨识的参数;disp('输入数据u(k)=');disp(M);disp('输出数据z(k)=');disp(z);disp('HL=');disp(H);disp('ZL=');disp(ZL);disp('theta=');disp(ES);%由估计值计算准则函数的值;J=(ZL-H*ES)'*(ZL-H*ES);disp('J(theta)=');disp(J);%================================%递推最小二乘法辨识系统;%一式s(k)=s(k-1)+k(k)(z(k)-h'(k)s(k-1));%二式k(k)=p(k-1)h(k)[h'(k)p(k-1)h(k)+1/w(k)]';%三式p(k)=[I-k(k)h'(k)]p(k-1);%s(0)=[0.01,0.01,...0.01]';%p(0)=10^6I;%新息ZX=z(k)-H'(k)S(k-1);%残差E=z(k)-H'(k)S(k);%准则函数J(k)%递推求解;z1=z(10:509);M1=M(10:509);P=10^6*eye(4);S(:,1)=[0.01;0.01;0.01;0.01]';S(:,2)=[0.01;0.01;0.01;0.01]';J1(1)=0;J1(2)=0;for i=3:500h=[-z1(i-1);-z1(i-2);M1(i-1);M1(i-2)];K=P*h*inv(h'*P*h+1);S(:,i)=S(:,i-1)+K*(z1(i)-h'*S(:,i-1)); %待估计参数变化; zx(i)=z1(i)-h'*S(:,i-1); %新息变化;E(i)=z1(i)-h'*S(:,i); %残差变化;J1(i)=J1(i-1)+(zx(i)*zx(i))/(h'*P*h+1); %准则函数变化; P=(eye(4,4)-K*h')*P;end%待估计参数过度曲线;figure(2);i=1:500;plot(i,S(1,:),i,S(2,:),i,S(3,:),i,S(4,:));title('待估计参数过度过程');legend('参数a1的过渡过程','参数a2的过渡过程','参数b1的过渡过程','参数b2的过渡过程');grid on;%disp(S(:,500));%新息变化曲线;figure(3);subplot(3,1,1);i=1:500;plot(i,zx(i));title('新息变化曲线');grid on;%残差变化曲线;subplot(3,1,2);i=1:500;plot(i,E(i));title('残差变化曲线');grid on;%准则函数变化曲线;subplot(3,1,3);i=1:500;plot(i,J1(i));title('准则函数变化曲线');grid on;%disp(J1(500));%================================%分别对一次完成，递推和辨识对象阶跃响应对比;%%一次完成法辨识阶跃响应对比;%一次辨识后对仿真对象和辨识模型进行阶跃响应对比;figure(4);subplot(2,1,1);%辨识模型阶跃响应;%z(k)+a1z(k-1)+a2z(k-2)=b1u(k-1)+b2u(k-2);a=[1,ES(1),ES(2)];b=[0,ES(3),ES(4)];%求离散系统阶跃响应;gn=dstep(b,a,30);stem(gn,'ro');grid on;hold on;%加阶跃响应后辨识对象的系统函数为：%(1+1z^-1+0.5z^-2)/(1-1.5z^-1+0.7z^-2);b1=[0,1.0,0.5];a1=[1,-1.5,0.7];gn1=dstep(b1,a1,30);stem(gn1,'bx');hold off;xlabel('k');ylabel('幅值h(k)');title('辨识对象和一次完成法辨识系统阶跃响应对比'); legend('估计值','理论值');%%递推法辨识阶跃响应对比;subplot(2,1,2);%辨识模型阶跃响应;%z(k)+a1z(k-1)+a2z(k-2)=b1u(k-1)+b2u(k-2);a=[1,S(1,500),S(2,500)];b=[0,S(3,500),S(4,500)];%求离散系统阶跃响应;gn=dstep(b,a,30);stem(gn,'ro');grid on;hold on;%加阶跃响应后辨识对象的系统函数为：%(1+1z^-1+0.5z^-2)/(1-1.5z^-1+0.7z^-2);b1=[0,1.0,0.5];a1=[1,-1.5,0.7];gn1=dstep(b1,a1,30);stem(gn1,'bx');hold off;电信工程学院xlabel('k');ylabel('幅值h(k)');title('辨识对象和递推法辨识系统阶跃响应对比'); legend('估计值','理论值');%================================end%================================%================================%子函数，产生周期大于500的逆M序列函数; function out=IMxulie(X,F,N)n=N;X1=X;F1=F;s=1;for i=1:nY=X1;t=F1(8)*Y(8);for j=7:-1:1t=xor(t,F1(j)*Y(j));endX1(1)=t;for k=2:8X1(k)=Y(k-1);endU(i)=Y(8);V(i)=xor(U(i),s);s=not(s);endout=V;end%================================- 10 -。