系统辨识与参数估计大作业

系统辨识与参数估计大作业第一题递推最小二乘估计参数考虑如上图所示的仿真对象，选择模型结构为：)()2()1()2()1()(2121kvkubkubkzakzakz，其中)(kv是服从)1,0(N正态分布的不相关随机噪声；输入信号)(ku采用4阶逆M序列，特征多项式取41)(sssF，幅度为1，循环周期为bitNp62；控制值，使数据的噪信比分别为10％，73％，100％三种情况。加权因子1)(k；数据长度L=500;初始条件取001.0)0(ˆ,10)0(6θIP，（1）利用递推最小二乘算法在线估计参数，（2）利用模型阶次辨识方法（AIC准则），确定模型的阶次。（3）估计噪声)(kv的方差和模型静态增益K（4）作出参数估计值随时间的变化图答：设过程的输入输出关系可以描述成()()()Tzkhknk()zk是输出量，()hk是可观测的数据向量，n(k)是均值为0的随机噪声)(kz)(ku)(kv21217.05.115.00.1zzzz217.05.111zz()(1),(2),(1),(2)Thkzkzkukuk1212,,,Taabb选取的模型为结构是1212()(1)(2)(1)(2)zkazkazkbukbuk12121.5,0.7,1.0,0.5aabb加权最小二乘参数估计递推算法RWLS的公式如下，11()(1)()()(1)()()()(1)()()()(1)()()()(1)TTTKkpkhkhkpkhkkkkKkzkhkkpkIKkhkpk为了把p(k)的对称性，可以把p(k)写成1()(1)()()()(1)()()TTpkpkKkKkhkpkhkk如果把()k设成1的时候，加权最小二乘法就退化成最小二乘法。用AIC准则定阶法来定阶，所用公式nnnnZHV(1),(2),(3),...,()TnZzzzzL1212,,...,,,...aTnnnaaabbb(0)(1)...(1)(0)(1)...(1)(1)(0)...(2)(1)(0)...(2).........................(1)(2)...()(1)(2)...()nzzznuuunzzznuuunHzLzLzLnuLuLuLn其中模型参数n和噪声()Vk方差的极大似然估计值为ML，2v12()1()()MLMLMLTTnnnnTvnnnnHHHZZHZHLAIC的定阶公式写成2()log4vAICnLn取1,2,3,4;n分别计算()AICn，找到使()AICn最小的那个n作为模型的阶次。一般说来，这样得到的模型阶次都能比较接近实际过程的真实阶次。信噪比为10%时：参数a1a2b1b2噪声方差静态增益模型阶次真值-1.50.710.51估计值-1.5190.722591.03140.509231.09517.56612信噪比为73%时：参数a1a2b1b2噪声方差静态增益模型阶次真值-1.50.710.51估计值-1.5190.722591.03140.509231.09517.56612信噪比为100%时：参数a1a2b1b2噪声方差静态增益模型阶次真值-1.50.710.51估计值-1.5190.722591.03140.50923估计值7.56612源程序：%function[a1a2b1b2nanbfangchaKk]=rwls(L,syn,Np)%na,nb为模型阶次，fangcha为噪声方差，Kk为静态增益a1=0;a2=0;b1=0;b2=0;na=0;nb=0;fangcha=0;Kk=0;L=500;Np=62;syn=1;x(1:4)=[1010];fori=1:Nptemp=xor(x(1),x(4));M(i)=x(4);forj=4:-1:2x(j)=x(j-1);endx(1)=temp;endS=ones(1,Np);%先产生一个全是1的序列Sifmod(Np,2)==0%判断Np是奇数还是偶数p=Np/2;elsep=(Np-1)/2;endforj=1:pS(2*j)=0;%将S序列的偶数位值均置为0，从而使S序列是0或1的方波序列endIM=xor(M,S);%使用M序列与方波序列S复合生成逆M序列IMu=IM*2-1;fori=(Np+1):Lu(i)=u(i-Np);endrandn('seed',2);v=randn(1,L);symsc;y(1)=0;y(2)=u(1);e(1)=c*v(1);e(2)=1.5*e(1)+c*v(2);fori=3:Ly(i)=1.5*y(i-1)-0.7*y(i-2)+u(i-1)+0.5*u(i-2);e(i)=1.5*e(i-1)-0.7*e(i-2)+c*v(i);endm=sum(e.^2);n=sum(y.^2);%c=solve('m/n-syn*syn','c');c=solve('m/n-syn*syn');c1=abs(double(c(1)));z(1)=c1*v(1);z(2)=u(1)+1.5*z(1)+c1*v(2);fori=3:Lz(i)=1.5*z(i-1)-0.7*z(i-2)+u(i-1)+0.5*u(i-2)+c1*v(i);endcta=zeros(4,L);cta(:,1)=[0.0010.0010.0010.001]';P=diag([10^610^610^610^6]);%k=2时h=[-z(1)0u(1)0]';K=P*h*inv(h'*P*h+1);P=P-K*K'*(h'*P*h+1);cta(:,2)=cta(:,1)+K*(z(2)-h'*cta(:,1));fork=3:Lh=[-z(k-1)-z(k-2)u(k-1)u(k-2)]';K=P*h*inv(h'*P*h+1);P=P-K*K'*(h'*P*h+1);cta(:,k)=cta(:,(k-1))+K*(z(k)-h'*cta(:,(k-1)));end%以上为参数估计值z=z';forna=1:4fornb=1:4A=zeros(L,na);B=zeros(L,nb);fori=1:Lforj=1:naifijA(i,j)=z(i-j);endendendfori=1:Lforj=1:nbifijB(i,j)=u(i-j);endendendH=[AB];cta1=inv(H'*H)*H'*z;%cta1为模型参数极大似然估计值cgm(na,nb)=(z-H*cta1)'*(z-H*cta1)/L;%cgm为噪声方差极大似然估计值AIC(na,nb)=L*log(cgm(na,nb))+2*(na+nb);endend[nanb]=find(AIC==min(min(AIC)));fangcha=cgm(na,nb)/(c1^2);a1=cta(1,500);a2=cta(2,500);b1=cta(3,500);b2=cta(4,500);Kk=(cta(3,L)+cta(4,L))/(1+cta(1,L)+cta(2,L));m=1:L;plot(m,cta(1,:),'b-',m,cta(2,:),'k-',m,cta(3,:),'y-',m,cta(4,:),'r-')第二大题卡尔曼滤波一个系统模型为)()()1(...1,0),()()()1(22211kwkxkxkkwkxkxkx同时有下列条件：（1）初始条件已知且有。T]0,0[)0(x（2）)(kw是一个标量零均值白高斯序列，且自相关函数已知为jkkwjwE)]()([，另外，我们有下列观测模型，即)1()1()1(...1,0),1()1()1(222111kvkxkzkkvkxkz，且有下列条件：（3）)1(1kv和)1(2kv是独立的零均值白高斯序列，且有jkkvjvE)]()([11，jkkvjvE2)]()([22，...2,1,0k（4）对于所有的j和k，)(kw与观测噪声过程)1(1kv和)1(2kv是不相关的，即0)]()([1kvjwE，0)]()([2kvjwE我们希望得到由观测矢量Tkzkzk)]1(),1([)1(21z估计状态矢量Tkxkxk)]1(),1([)1(21x的卡尔曼滤波的公式表示，并求解以下问题：（a）求出卡尔曼增益矩阵，并得出最优估计)1(kx和观测矢量)1(),...,2(),1(kzzz之间的递推关系。（b）用模拟数据确定状态矢量)(kx的估计值)|(ˆkkx，10...2,1,0k，并画出当10...2,1,0k时，)|(ˆ1kkx，)|(ˆ2kkx的图（c）通常，状态矢量的真实值是得不到的，但是为了用作图来说明问题，表1和表2给出了状态矢量元素的真实值。对于10...2,1,0k，在同一幅图中画出真实值和在（b）中确定的)(1kx的估计值。对)(2kx重复这一过程。当k从1变到10时，对每个元素2,1i，计算并画出各自的误差图，即)|(ˆ)(kkxkxii。（d）当k从1变到10时，通过用由卡尔曼滤波器决定的状态误差协方差矩阵画出)]|([21kkE和)]|([22kkE，而2,1),|(ˆ)()|(ikkxkxkkiii（e）讨论一下（C）中计算的误差与（d）中方差之间的关系。表1观测值时间下标k观测值)(1kz观测值)(2kz13．296919692．1013429423．387365150．4754079737．028306413．1768889849．712125212．49811140511．420183152．91992424615．979805836．17307616722．069342855．42519274828．302127813．05365741930．446838315．980511411038．758755954．51016361表2由模拟得到的实际状态值时间下标k实际状态值)(1kx实际状态值)(2kx00.000000000.0000000011.654287141.6542871423.503007021.8487198835.997852922.4755222249.150407403.17187816512.508739103.35833170616.921925944.41318684721.344833524.42290758825.893351444.54851792931.541353305.648001861036.936056705.39447034答：（a）卡尔曼增益矩阵：kkkTT1HPC(CPCR)’’估计值与观测值之间的递归关系为：kk-1k-1kkkkkXAXH(YCAX)（b）状态矢量估计值的计算框图：（c）)(1kkx和)(2kkx的图：（d）真实值与估计值的比较图：+1kH+1Z1kA1kC1ky1ˆkx1ˆkx1'ˆkx各自的误差图：（e）通过用卡尔曼滤波器的状态误差协方差矩阵画出的)]([21kkE和)]([22kkE：（f）分析：（e）中的方差是（d）中的误差平方后取均值，是均方误差。误差直接由真实值减去估计值，有正有负，而均方误差没有这个缺陷，更能综合的表示滤波的效果。源程序：%PKXZ均为2*2矩阵，用三维矩阵表示他们，记录每一个k值下的情况，如P(:,:,k)%k从1开始,k=1时取初值clearall;P(:,:,1)=[100;010];%P(0|0)初值，P表示P(k|k)X(:,:,1)=[0;0];%X(0)初值fy=[11;01]