系统辨识与参数估计大作业教案_第1页
系统辨识与参数估计大作业教案_第2页
系统辨识与参数估计大作业教案_第3页
系统辨识与参数估计大作业教案_第4页
系统辨识与参数估计大作业教案_第5页
已阅读5页,还剩8页未读 继续免费阅读

下载本文档

版权说明:本文档由用户提供并上传,收益归属内容提供方,若内容存在侵权,请进行举报或认领

文档简介

1、系统辨识与参数估计大作业第一题递推最小二乘估计参数考虑如上图所示的仿真对象,选择模型结构为: z(k) a1z(k -1) a2z(k -2) =bu(k -1) b2u(k -2) v(k),其中v(k)是服从N(0,1)正态分布的不相关随机噪声; 输入信号u(k)采用4阶逆M序列, 特征多项式取F(s) =13s$ s4,幅度为1,循环周期为Np =62bit ;控制九值,使数据 的噪信比分别为 10%, 73%, 100%三种情况。加权因子 A(k)=1;数据长度L=500;初始 条件取 P (0) =1061 , ?(0) =0.001 ,利用递推最小二乘算法在线估计参数,利用模型阶次

2、辨识方法(AIC准则),确定模型的阶次。估计噪声v(k)的方差和模型静态增益 K作出参数估计值随时间的变化图答:设过程的输入输出关系可以描述成z(k) = hT (k)H + n(k)z(k)是输出量,h(k)是可观测的数据向量,n(k)是均值为0的随机噪声,Jh(k) - I -z(k -1), -z(k -2), u(k -1),u(k -2)3 - 1al ,a2,b1, b2 T选取的模型为结构是z(k) = -a1Z(k 一1) 一a2z(k 2) bu1 (k -1) bu2(k -2) ai 二1.5,a2 =0.7,“ =1.0,b2 =0.5加权最小二乘参数估计递推算法RWL

3、S的公式如下,1K(k) = p(k -1)h(k) hT(k)p(k-1)h(k) TOC o 1-5 h z IL(k)9(k) =0(k 1) + K(k) z(k) hT(k声(k 1)p(k)-|I -K(k)hT(k) p(k -1)为了把p(k)的对称性,可以把 p(k)写成T T1p(k) = p(k-1)-K(k)KT(k) hT(k)p(k -1)h(k)IL”(k)如果把A(k)设成1的时候,加权最小二乘法就退化成最小二乘法。用AIC准则定阶法来定阶,所用公式乙VnZn = lz(1),z(2), z(3),., z(L)T-jT*n =回,-a2,-ana回以。z(0)

4、z(-1).z(1 - n)u(0)u(-1).u(1-n)z(1)z(0).z(2 -n)u(1)u(0).u(2 - n)Hn= ._z(L-1) z(L-2) . z(L-n) u(L-1) u(L-2) . u(L-n)_2其中模型参数 我和 噪声V(k)方差的极大似然估计值为19 ML ,仃v=(HnTHn)HnTZn21T(Z = (7 - H ) (Z - H )v l (乙n n n ml )(乙 n n ml )AIC的定阶公式写成2AIC (n) = L log 4n取n=1,2,3,4;分别计算AIC(n),找到使AIC(n)最小的那个n作为模型的阶次。一般说来,这样得到

5、的模型阶次都能比较接近实际过程的真实阶次。 信噪比为10%时:060100150 200250300350400 450500参数ala2b1b2噪声力差静态增益模型阶次真值-1.50.710.51估计值-1.5190.722591.03140.509231.09517.56612信噪比为 参数73% 时:a1a2b1b2噪声力差静态增益模型阶次真值-1.50.710.51传计值-1.5190.722591.03140.509231.09517.566124111fli111nl050100150200250300350400450500信噪比为 参数100% 时:ala2b1b2噪声力差静态

6、增益模型阶次真值-1.50.710.51传计值-1.5190.722591.03140.50923情计值7.56612源程序:%function a1 a2 bl 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:Nptemp=xor(x(1),x(4);M(i)=x(4);for j=4:-1:2x(j)=x(j-1);endx(1)

7、=temp;endS=ones(1,Np);%先产生一个全是 1的序列Sif mod(Np,2)=0% 判断Np是奇数还是偶数p=Np/2;else p=(Np-1)/2;endfor j=1:pS(2*j)=0;%将S序列的偶数位值均置为0,从而使S序列是0或1的方波序列endIM=xor(M,S); %使用M序列与方波序列 S复合生成逆 M序列IMu=IM*2-1;for i=(Np+1):Lu(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);

8、for i=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.A2);n=sum(y.A2);%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: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);c

9、ta(:,1)=0.001 0.001 0.001 0.001;P=diag(10A6 10A6 10A6 10A6);%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: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;for

10、na=1:4for nb=1:4A=zeros(L,na);B=zeros(L,nb);for i=1:Lfor j=1:naif ijA(i,j)=z(i-j);endendendfor i=1:Lfor j=1:nbif ijB(i,j)=u(i-j);endendendH=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

11、);fangcha=cgm(na,nb)/(c1A2);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-) 第二大题卡尔曼滤波一个系统模型为x1(k 1) = x1(k) x2(k) w(k), k = 0,1 x2 (k 1) = x2 (k) w(k)同时有下列条件:初始条件已知且有。x(0)=0,0Tw(k)是

12、一个标量零均值白高斯序列,且自相关函数已知为 Ew(j)w(k) =3jk ,乙(k 1) = x1(k 1) v1(k 1), k =0,1另外,我们有下列观测模型,即,,且z2 (k 1) = x2(k 1) v2(k 1)有下列条件:v1(k+1)和v2(k+1)是独立的零均值白高斯序列,且有Ev1(jW(k)=染,Ev2(j)v2(k)=25jk, k = 0,1,2.对于所有的j和k , w(k)与观测噪声过程 *+1)和丫2%+1)是不相关的,即Ew(j)v1(k) =0, Ew(jM(k)=0我们希望得到由观测矢量z(k+1)=乙(k+1),z2(k+1)T估计状态矢量x(k +

13、1) =x1(k+1),x2(k +1)T的卡尔曼滤波的公式表示,并求解以下问题:求出卡尔曼增益矩阵,并得出最优估计x(k +1)和观测矢量z(1), z(2),., z(k +1)之间的递推关系。用模拟数据确定状态矢量x(k)的估计值 父(k|k), k =0,1,2.10 ,并画出当k=0,1,2.10 时,R(k|k), x2(k|k)的图通常,状态矢量的真实值是得不到的,但是为了用作图来说明问题,表 1和表2给出了状态矢量元素的真实值。对于k = 0,1,2.10 ,在同一幅图中画出真实值和在(b)中确定的x1(k)的估计值。对x2(k)重复这一过程。当 k从1变到10时,对每个元素i

14、 =1,2 ,计算并画出各自的误差图,即xi(k)-?i(k|k)o当k从i变到10时,通过用由卡尔曼滤波器决定的状态误差协方差矩阵画出E 2 (k | k)和 Ezl (k | k),而备(k | k)=为(k) ? (k | k),i = 1,2讨论一下(C)中计算的误差与(d)中方差之间的关系。表1观测值时间卜标k观测值z1 (k)观测值Z2(k)13. 296919692. 1013429423. 387365150. 4754079737. 028306413. 1768889849. 712125212. 49811140511, 420183152. 91992424615. 9

15、79805836. 069342855. 42519274828, 302127813. 05365741930. 446838315. 980511411038. 758755954. 51016361表2由模拟得到的实际状态值时间卜标k实际状态值x1(k)实际状态值x2(k)00.000000000.0000000011.654287141.6542871423.503007021.8487198835.997852922.4755222249.150407403508739103.35833170616.921925944.413186

16、84721.344833524.42290758825.893351444.54851792931.541353305.648001861036.936056705.39447034答:(a)卡尔曼增益矩阵:Hk=Pk*CTM(CMPkMCT+R),估计值与观测值之间的递归关系为:Xk=Ak Xk-i Hk (Yk-Ck Ak Xk-i)(b)状态矢量估计值的计算框图:(c) xi (k/k)和 x2 (k/k)的图:(d)真实值与估计值的比较图:(e)通过用卡尔曼滤波器的状态误差协方差矩阵画出的E与2(k/k)和E42(k/k):状态矢量X每个单独分量估计的方差(f)分析:(e)中的方差是(

17、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;% rH= 1 0;0 1;Q=1;R=1 0;0 2;Z(:,:,2)=3.29691969;2.10134294;Z(:,:,3)=3.387

18、36515;0.47540797;Z(:,:,4)=7.02830641;3.17688898;Z(:,:,5)=9.71212521;2.49811140;Z(:,:,6)=11.42018315;2.91992424;Z(:,:,7)=15.97870583;6.17307616;Z(:,:,8)=22.06934285;5.42519274;Z(:,:,9)=28.30212781;3.05365741;Z(:,:,10)=30.44683831;5.98051141;Z(:,:,11)=38.75875595;4.51016361;for k=1:10P1(:,:,k尸fy*P(:,:

19、,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.65428714,3.50300702,5.997852924,9.15040740,12.50873910,.16.92192594,21.344

温馨提示

  • 1. 本站所有资源如无特殊说明,都需要本地电脑安装OFFICE2007和PDF阅读器。图纸软件为CAD,CAXA,PROE,UG,SolidWorks等.压缩文件请下载最新的WinRAR软件解压。
  • 2. 本站的文档不包含任何第三方提供的附件图纸等,如果需要附件,请联系上传者。文件的所有权益归上传用户所有。
  • 3. 本站RAR压缩包中若带图纸,网页内容里面会有图纸预览,若没有图纸预览就没有图纸。
  • 4. 未经权益所有人同意不得将文件中的内容挪作商业或盈利用途。
  • 5. 人人文库网仅提供信息存储空间,仅对用户上传内容的表现方式做保护处理,对用户上传分享的文档内容本身不做任何修改或编辑,并不能对任何下载内容负责。
  • 6. 下载文件中如有侵权或不适当内容,请与我们联系,我们立即纠正。
  • 7. 本站不保证下载资源的准确性、安全性和完整性, 同时也不承担用户因使用这些下载资源对自己和他人造成任何形式的伤害或损失。

评论

0/150

提交评论