系统辨识与参数估计大作业(共12页).doc
精选优质文档-倾情为你奉上系统辨识与参数估计大作业 第一题 递推最小二乘估计参数考虑如上图所示的仿真对象，选择模型结构为：，其中是服从正态分布的不相关随机噪声；输入信号采用4阶逆M序列，特征多项式取，幅度为1，循环周期为；控制值，使数据的噪信比分别为10，73，100三种情况。加权因子；数据长度L=500;初始条件取，（1） 利用递推最小二乘算法在线估计参数，（2） 利用模型阶次辨识方法（AIC准则），确定模型的阶次。（3） 估计噪声的方差和模型静态增益（4） 作出参数估计值随时间的变化图答：设过程的输入输出关系可以描述成是输出量，是可观测的数据向量，是均值为0的随机噪声选取的模型为结构是加权最小二乘参数估计递推算法RWLS的公式如下，为了把p(k)的对称性，可以把p(k)写成如果把设成1的时候，加权最小二乘法就退化成最小二乘法。用AIC准则定阶法来定阶，所用公式其中模型参数和 噪声方差的极大似然估计值为 ，AIC的定阶公式写成取分别计算，找到使最小的那个作为模型的阶次。一般说来，这样得到的模型阶次都能比较接近实际过程的真实阶次。信噪比为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 a1 a2 b1 b2 na nb fangcha Kk=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)=1 0 1 0;for i=1:Np temp=xor(x(1),x(4); M(i)=x(4); for j=4:-1:2 x(j)=x(j-1); end x(1)=temp;endS=ones(1,Np);%先产生一个全是1的序列Sif mod(Np,2)=0%判断Np是奇数还是偶数 p=Np/2;else p=(Np-1)/2;endfor j=1:p S(2*j)=0;%将S序列的偶数位值均置为0，从而使S序列是0或1的方波序列end IM=xor(M,S); %使用M序列与方波序列S复合生成逆M序列IMu=IM*2-1;for i=(Np+1):L u(i)=u(i-Np);endrandn('seed',2);v=randn(1,L);syms c;y(1)=0;y(2)=u(1);e(1)=c*v(1);e(2)=1.5*e(1)+c*v(2);for i=3:L y(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);for i=3:L z(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.001 0.001 0.001 0.001'P=diag(106 106 106 106);%k=2时h=-z(1) 0 u(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);for k=3:L h=-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' for na=1:4 for nb=1:4 A=zeros(L,na); B=zeros(L,nb);for i=1:L for j=1:na if i>j A(i,j)=z(i-j); end endendfor i=1:L for j=1:nb if i>j B(i,j)=u(i-j); end endendH=A B;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); endendna nb=find(AIC=min(min(AIC);fangcha=cgm(na,nb)/(c12);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） 初始条件已知且有。（2） 是一个标量零均值白高斯序列，且自相关函数已知为，另外，我们有下列观测模型，即，且有下列条件：（3） 和是独立的零均值白高斯序列，且有 ，（4） 对于所有的和，与观测噪声过程和是不相关的，即，我们希望得到由观测矢量估计状态矢量的卡尔曼滤波的公式表示，并求解以下问题：（a） 求出卡尔曼增益矩阵，并得出最优估计和观测矢量之间的递推关系。（b） 用模拟数据确定状态矢量的估计值，并画出当时，的图（c） 通常，状态矢量的真实值是得不到的，但是为了用作图来说明问题，表1和表2给出了状态矢量元素的真实值。对于，在同一幅图中画出真实值和在（b）中确定的的估计值。对重复这一过程。当从1变到10时，对每个元素，计算并画出各自的误差图，即。（d） 当从1变到10时，通过用由卡尔曼滤波器决定的状态误差协方差矩阵画出和，而（e） 讨论一下（C）中计算的误差与（d）中方差之间的关系。 表1观测值时间下标观测值观测值1322303734925112615672258283930510384表2由模拟得到的实际状态值时间下标实际状态值实际状态值00.0.11.1.23.1.35.2.49.3.512.3.616.4.721.4.825.4.931.5.1036.5.答：（a）卡尔曼增益矩阵：估计值与观测值之间的递归关系为：（b）状态矢量估计值的计算框图：+（c）和的图：（d）真实值与估计值的比较图：各自的误差图：（e）通过用卡尔曼滤波器的状态误差协方差矩阵画出的和：（f）分析：（e）中的方差是（d）中的误差平方后取均值，是均方误差。误差直接由真实值减去估计值，有正有负，而均方误差没有这个缺陷，更能综合的表示滤波的效果。源程序：%P K X Z均为2*2矩阵，用三维矩阵表示他们，记录每一个k值下的情况，如P(:,:,k)%k从1开始,k=1时取初值clear all;P(:,:,1)=10 0;0 10;%P(0|0)初值，P表示P(k|k)X(:,:,1)=0;0 ; %X(0)初值fy= 1 1;0 1; %tao=1;1;%H= 1 0;0 1;Q=1;R=1 0;0 2;Z(:,:,2)=3.;2.;Z(:,:,3)=3.;0.;Z(:,:,4)=7.;3.;Z(:,:,5)=9.;2.;Z(:,:,6)=11.;2.;Z(:,:,7)=15.;6.;Z(:,:,8)=22.;5.;Z(:,:,9)=28.;3.;Z(:,:,10)=30.;5.;Z(:,:,11)=38.;4.;for k=1:10 P1(:,:,k)=fy*P(:,:,k)*fy'+tao*Q*tao' %P1(:,:,k)表示P(k+1|k) K(:,:,k+1)=P1(:,:,k)*H'*inv(H*P1(:,:,k)*H'+R); X(:,:,k+1)=fy*X(:,:,k)+K(:,:,k+1)*(Z(:,:,k+1)-H*fy*X(:,:,k); P(:,:,k+1)=(eye(2,2)-K(:,:,k+1)*H)*P1(:,:,k);endgx1(1:11)=X(1,1,1:11); %gx1,gx2分别表示X的估计值gx2(1:11)=X(2,1,1:11);plot(0:10,gx1,0:10,gx2);legend('x1估计值','x2估计值','location','northwest');title('状态矢量X的估计值x1(k)和x2(k)');xlabel('k');ylabel('X');x1=0,1.,3.,5.,9.,12.,. 16.,21.,25.,31.,36.;x2=0,1.,1.,2.,3.,3.,. 4.,4.,4.,5.,5.;figure;plot(0:10,gx1,0:10,x1,0:10,gx2,0:10,x2);legend('x1估计值','x1真实值','x2估计值','x2真实值','location','northwest');title('状态矢量X的估计值与真实值之间的关系');xlabel('k');ylabel('X');delta1=x1-gx1;%delta1,delta2分别表示真值与估计值的差值delta2=x2-gx2;figure;plot(1:10,delta1(2:11),1:10,delta2(2:11);legend('x1的误差','x2的误差','location','northwest');title('状态矢量X的误差图');xlabel('k');ylabel('X');e12(1:11)=P(1,1,1:11);%e12,e22分别表示x1,x2方差,P287e22(1:11)=P(2,2,1:11);figure;plot(1:10,e12(2:11),1:10,e22(2:11);legend('x1的方差','x2的方差');title('状态矢量X每个单独分量估计的方差');专心-专注-专业
