系统辨识-最小二乘法

5.3Matlab源程序：%最小二乘估计clearu=[1.1470.201-0.787-1.589-1.0520.8661.1521.5730.6260.433-0.9850.810-0.0440.947-1.474-0.719-0.086-1.0991.4501.1510.4851.6330.0431.3261.706-0.3400.8900.1441.177-0.390];n=normrnd(0,sqrt(0.1),1,31);%方差为0.1z=zeros(1,30);fork=3:31z(k)=-1.642*z(k-1)-0.715*z(k-2)+0.39*u(k-1)+0.35*u(k-2)+n(k)+1.642*n(k-1)+0.715*n(k-2);endh0=[-z(2)-z(1)u(2)u(1)]';HLT=[h0,zeros(4,28)];fork=3:30h1=[-z(k)-z(k-1)u(k)u(k-1)]';HLT(:,k-1)=h1;endHL=HLT';y=[z(3);z(4);z(5);z(6);z(7);z(8);z(9);z(10);z(11);z(12);z(13);z(14);z(15);z(16);z(17);z(18);z(19);z(20);z(21);z(22);z(23);z(24);z(25);z(26);z(27);z(28);z(29);z(30);z(31)];%求出FAIc1=HL'*HL;c2=inv(c1);c3=HL'*y;c=c2*c3;%display('方差=0.1时，最小二乘法估计辨识参数θ如下：');a1=c(1)a2=c(2)b1=c(3)b2=c(4)%递推最小二乘法估计clearu=[1.1470.201-0.787-1.589-1.0520.8661.1521.5730.6260.433-0.9850.810-0.0440.947-1.474-0.719-0.086-1.0991.4501.1510.4851.6330.0431.3261.706-0.3400.8900.1441.177-0.390];z(2)=0;z(1)=0;n=normrnd(0,sqrt(0.1),1,31);%方差为0.1fork=3:31z(k)=-1.642*z(k-1)-0.715*z(k-2)+0.39*u(k-1)+0.35*u(k-2)+n(k)+1.642*n(k-1)+0.715*n(k-2);endc0=[0.0010.0010.0010.001]';%直接给出被辨识参数的初始值,即一个充分小的实向量p0=10^6*eye(4,4);%直接给出初始状态P0，即一个充分大的实数单位矩阵E=0.000000005;%取相对误差E=0.000000005c=[c0,zeros(4,30)];%被辨识参数矩阵的初始值及大小e=zeros(4,30);%相对误差的初始值及大小fork=3:30;%开始求Kh1=[-z(k-1),-z(k-2),u(k-1),u(k-2)]';x=h1'*p0*h1+1;x1=inv(x);%开始求K(k)k1=p0*h1*x1;%求出K的值d1=z(k)-h1'*c0;c1=c0+k1*d1;%求被辨识参数ce1=c1-c0;%求参数当前值与上一次的值的差值e2=e1./c0;%求参数的相对变化e(:,k)=e2;%把当前相对变化的列向量加入误差矩阵的最后一列c0=c1;%新获得的参数作为下一次递推的旧参数c(:,k)=c1;%把辨识参数c列向量加入辨识参数矩阵的最后一列p1=p0-k1*k1'*[h1'*p0*h1+1];%求出p(k)的值p0=p1;%给下次用ife2=Ebreak;%如果参数收敛情况满足要求，终止计算endend%display('方差为0.1递推最小二乘法辨识后的结果是：');a1=c(1,:);a2=c(2,:);b1=c(3,:);b2=c(4,:);%display('a1,a2,b1,b2经过递推最小二乘法辨识的结果是：');fori=3:31;if(c(1,i)==0)q1=c(1,i-1);break;endendfori=3:31;if(c(2,i)==0)q2=c(2,i-1);break;endendfori=3:31;if(c(3,i)==0)q3=c(3,i-1);break;endendfori=3:31;if(c(4,i)==0)q4=c(4,i-1);break;endenda1=q1a2=q2b1=q3b2=q4%辅助变量递推最小二乘法估计clearna=2;nb=2;siitt=[1.6420.7150.390.35]';siit0=0.001*eye(na+nb,1);p=10^6*eye(na+nb);siit(:,1)=siit0;y(2)=0;y(1)=0;x(1)=0;x(2)=0;j=0;u=[1.1470.201-0.787-1.589-1.0520.8661.1521.5730.6260.433-0.9850.810-0.0440.947-1.474-0.719-0.086-1.0991.4501.1510.4851.6330.0431.3261.706-0.3400.8900.1441.177-0.390];n=normrnd(0,sqrt(0.1),1,31);%方差为0.1fork=3:31;h=[-y(k-1),-y(k-2),u(k-1),u(k-2)]';y(k)=h'*siitt+n(k)+1.642*n(k-1)+0.715*n(k-2);hx=[-x(k-1),-x(k-2),u(k-1),u(k-2)]';kk=p*hx/(h'*p*hx+1);p=[eye(na+nb)-kk*h']*p;siit(:,k-1)=siit0+kk*[y(k)-h'*siit0];x(k)=hx'*siit(:,k-1);j=j+(y(k)-h'*[1.6420.7150.390.35]')^2;e=max(abs((siit(:,k-1)-siit0)./siit0));ess(:,k-2)=siit(:,k-1)-siitt;siit0=siit(:,k-1);enda1=siit0(1)a2=siit0(2)b1=siit0(3)b2=siit0(4)clear%广义最小二乘估计clear;nn=normrnd(0,sqrt(0.1),1,31)';%方差为0.1uk=[1.1470.201-0.787-1.589-1.0520.8661.1521.5730.6260.433-0.9580.810-0.0440.947-1.474-0.719-0.086-1.0991.4501.1510.4851.6330.0431.3261.706-0.3400.8901.1441.177-0.390];yk(1)=0;yk(2)=0;fori=1:29;yk(i+2)=-1.642*yk(i+1)-0.715*yk(i)+0.39*uk(i+1)+0.35*uk(i)+nn(i+2)+1.642*nn(i+1)+0.715*nn(i);end;fori=1:29;A(i,:)=[-yk(i+1)-yk(i)uk(i+1)uk(i)];endsiit=inv(A'*A)*A'*(yk(3:31)+nn(2:30)')';e(1)=yk(1);e(2)=yk(2)+siit(1)*yk(1)-siit(3)*uk(1);fori=3:31;e(i)=yk(i)+siit(1)*yk(i-1)+siit(2)*yk(i-2)-siit(3)*uk(i-1)-siit(4)*uk(i-2);endfori=1:29;fai(i,:)=[-e(i+1)-e(i)];endf=inv(fai'*fai)*fai'*e(3:31)';fori=3:31;yk(i)=yk(i)+f(1)*yk(i-1)+f(2)*yk(i-2);endyk(2)=yk(2)+f(1)*yk(1);fori=3:30;uk(i)=uk(i)+f(1)*uk(i-1)+f(2)*uk(i-2);enduk(2)=uk(2)+f(1)*uk(1);forj=1:30fori=1:29;A(i,:)=[-yk(i+1)-yk(i)uk(i+1)uk(i)];endsiit=inv(A'*A)*A'*yk(3:31)';e(1)=yk(1);e(2)=yk(2)+siit(1)*(yk(1))-siit(3)*uk(1);fori=3:31;e(i)=yk(i)+siit(1)*(yk(i-1))+siit(2)*(yk(i-2))-siit(3)*uk(i-1)-siit(4)*uk(i-2);endfori=1:29;fai(i,:)=[-e(i+1)-e(i)];endf=inv(fai'*fai)*fai'*e(3:31)';k1(j)=f(1);k2(j)=f(2);fori=3:31;yk(i)=yk(i)+f(1)*(yk(i-1))+f(2)*(yk(i-2));endyk(2)=yk(2)+f(1)*yk(1);fori=3:30uk(i)=uk(i)+f(1)*uk(i-1)+f(2)*uk(i-2);enduk(2)=uk(2)+f(1)*uk(1);endsiit'%增广矩阵估计参数clearu=[1.1470.201-0.787-1.589-1.0520.8661.1521.5730.6260.433-0.9850.810-0.0440.947-1.474-0.719-0.086-1.0991.4501.1510.4851.6330.0431.3261.706-0.3400.8900.1441.177-0.390];n=normrnd(0,sqrt(0.1),1,31);%方差为0.1z=zeros(7,30);zs=zeros(7,30);zm=zeros(7,30);zmd=zeros(7,30);%输出采样、不考虑噪声时系统输出、不考虑噪声时模型输出、模型输出矩阵的大小z(2)=0;z(1)=0;zs(2)=0;zs(1)=0;zm(2)=0;zm(1)=0;zmd(2)=0;zmd(1)=0;%给输出采样、不考虑噪声时系统输出、不考虑噪声时模型输出、模型输出赋初值fork=3:31z(k)=-1.642*z(k-1)-0.715*z(k-2)+0.39*u(k-1)+0.35*u(k-2)+n(k)+1.642*n(k-1)+0.715*n(k-2);endc0=[0.0010.0010.0010.001]';%直接给出被辨识参数的初始值,即一个充分小的实向量p0=10^6*eye(4,4);%直接给出初始状态P0，即一个充分大的实数单位矩阵E=0.000000005;%取相对误差E=0.000000005c=[c0,zeros(4,30)];%被辨识参数矩阵的初始值及大小e=zeros(4,30);%相对误差的初始值及大小fork=3:30;%开始求Kh1=[-z(k-1),-z(k-2),u(k-1),u(k-2)]';x=h1'*p0*h1+1;x1=inv(x);%开始求K(k)k1=p0*h1*x1;%求出K的值d1=z(k)-h1'*c0;c1=c0+k1*d1;%求被辨识参数czs(k)=-1.642*z(k-1)-0.715*z(k-2)+0.39*u(k-1)+0.35*u(k-2);%系统在M序列的输入下不考虑扰动时的输出响应zm(k)=[-z(k-1),-z(k-2),u(k-1),u(k-2)]*[c1(1);c1(2);c1(3);c1(4)];%模型在输入下不考虑扰动时的的输出响应zmd(k)=h1'*c1;%模型在输入下的的输出响应e1=c1-c0;e2=e1./c0;%求参数的相对变化e(:,k)=e2;%把当前相对