系统辨识与参数估计Fixed.docVIP

第一题 递推最小二乘估计参数 考虑如下图所示的仿真对象，图中是服从正态分布的不相关随机噪声；输入信号采用4阶逆M序列，特征多项式取，幅度为1，循环周期为；控制值，使数据的噪信比为73％。加权因子；数据长度L=500;初始条件取， 利用递推最小二乘算法在线估计参数， 利用模型阶次辨识方法（AIC准则），确定模型的阶次。 估计噪声的方差和模型静态增益 作出参数估计值随时间的变化图 答： 设过程的输入输出关系可以描述成 是输出量，是可观测的数据向量，是均值为0的随机噪声 选取的模型为结构是 加权最小二乘参数估计递推算法RWLS的公式如下， 为了把p(k)的对称性，可以把p(k)写成 如果把设成1的时候，加权最小二乘法就退化成最小二乘法。 用AIC准则定阶法来定阶，所用公式 其中模型参数和 噪声方差的极大似然估计值为 ， AIC 的定阶公式写成 取分别计算，找到使最小的那个作为模型的阶次。一般说来，这样得到的模型阶次都能比较接近实际过程的真实阶次。 源程序： clear close all; L=500; x1=1;x2=1;x3=1;x4=0;%移位寄存器初值 S=1;%方波初值 un=zeros(L,1); for i=1:L+4 IM=xor(S,x4);%进行异或运算，产生逆M序列 if IM==0 u(i)=-1; %fixed else u(i)=1; end un(i,1)=u(i); S=not(S);%产生方波 M(i)=xor(x1,x4);%进行异或运算，产生M序列 x4=x3;x3=x2;x2=x1;x1=M(i);%寄存器移位 end z(2)=0;z(1)=0;%取z的前二个初始值为零 for k=3:L+4; z(k)=1.5*z(k-1)-0.7*z(k-2)+1.0*u(k-1)+0.5*u(k-2)+normrnd(0,1,1,1);%给出理想的辨识输出采样信号 zn(k,1)=z(k); end c0=[0.001 0.001 0.001 0.001 ]; %a1 a2 b1 b2给出被辨识参数的初始值 p0=10^6*eye(4,4); E=0.000000005;%相对误差E=0.000000005 c=[c0,zeros(4,499)];%被辨识参数矩阵的初始值及大小 e=zeros(4,500);%相对误差的初始值及大小 for k=3:500; %开始求K h1=[-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;%求被辨识参数c e1=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;%给下次用 if abs(e2)=E break;%若参数收敛满足要求，终止计算 end%小循环结束 end%大循环结束 % c %显示被辨识参数 % e %显示辨识结果的收敛情况 % e2; a1=c(1,:); a2=c(2,:); b1=c(3,:); b2=c(4,:); figure(1);%第1个图形 i=1:500; plot(i,a1,r,i,a2,--,i,b1,g:,i,b2,k-.,LineWidth,2) %画出a1，a2，a3,b1，b2的各次辨识结果 title(各参数最小二乘法辨识结果)%图形标题 legend(a1,a2,b1,b2,0); grid on; %噪声方差的估计值和AIC准则定阶 % 模型阶次n=1 figure(2); for i=1:L H1(i,1)=z(i);; H1(i,2)=u(i); end estimate1=inv(H1*H1)*H1*z(2:501); D1=(z(2:501)-H1*estimate1)*(z(2:501)-H1*estimate1)/L; %噪声方差的估计值 AIC1=L*log(D1)+4*1; % 模型阶次n=2 for i=1:L H2(i,1)=z(i