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

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

第一题 递推最小二乘估计参数

?v(k) 1?1?21?1.5z?0.7z u(k) 1.0z?0.5z?1?21?1.5z?0.7z?1?2z(k)

考虑如上图所示的仿真对象，选择模型结构为：

z(k)?a1z(k?1)?a2z(k?2)?b1u(k?1)?b2u(k?2)?v(k)

，其中v(k)是服从N(0,1)正态分布的不相关随机噪声；输入信号u(k)采用4阶逆M序列，特征多项式取F(s)?1?s?s4，幅度为1，循环周期为Np?62bit；控制?值，使数据的噪信比分别为10％，73％，100％三种情况。加权因子?(k)?1；数据长度L=500;初始

?(0)?0.001， 条件取P(0)?106I,θ（1） 利用递推最小二乘算法在线估计参数，

（2） 利用模型阶次辨识方法（AIC准则），确定模型的阶次。 （3） 估计噪声v(k)的方差和模型静态增益K （4） 作出参数估计值随时间的变化图 答：

设过程的输入输出关系可以描述成z(k)?h(k)??n(k)

Tz(k)是输出量，h(k)是可观测的数据向量，n(k)是均值为0的随机噪声

h(k)???z(k?1),?z(k?2),u(k?1),u(k?2)?

T???a1,a2,b1,b2?

选取的模型为结构是

Tz(k)??a1z(k?1)?a2z(k?2)?bu1(k?1)?bu2(k?2) a1??1.5,a2?0.7,b1?1.0,b2?0.5

加权最小二乘参数估计递推算法RWLS的公式如下，

?1?1?K(k)?p(k?1)h(k)?hT(k)p(k?1)h(k)??(k)????(k)??(k?1)?K(k)?z(k)?hT(k)?(k?1)???Tp(k)??I?K(k)h(k)???p(k?1)

为了把p(k)的对称性，可以把p(k)写成

?1?p(k)?p(k?1)?K(k)KT(k)?hT(k)p(k?1)h(k)??

?(k)??如果把?(k)设成1的时候，加权最小二乘法就退化成最小二乘法。 用AIC准则定阶法来定阶，所用公式

Zn?Hn?n?Vn

Zn??z(1),z(2),z(3),...,z(L)?

T?n????a1,?a2,...?an,b1,b2,...bn??

aTz(?1)?z(0)?z(1)z(0)Hn???......??z(L?1)z(L?2)...z(1?n)u(0)u(?1)...z(2?n)u(1)u(0)...............z(L?n)u(L?1)u(L?2)...u(1?n)?...u(2?n)?.? ?......?...u(L?n)?2其中模型参数?n和 噪声V(k)方差的极大似然估计值为?ML ，?v

?ML?(HnTHn)?1HnTZn 1T??(Zn?Hn?ML)(Zn?Hn?ML)L2vAIC 的定阶公式写成

AIC(n)?Llog?v?4n

取n?1,2,3,4;分别计算AIC(n)，找到使AIC(n)最小的那个n作为模型的阶次。一般说来，

2这样得到的模型阶次都能比较接近实际过程的真实阶次。 信噪比为10%时：

参数 a1 a2 b1 b2 噪声方差 静态增益

真值 -1.5 0.7 1 0.5 1 估计值

-1.519

0.72259

1.0314

0.50923

1.0951

7.5661

信噪比为73%时：

参数 a1 a2 b1 b2 噪声方差 静态增益

真值 -1.5 0.7 1 0.5 1 估计值

-1.519

0.72259

1.0314

0.50923

1.0951

7.5661

模型阶次

2

模型阶次

2

信噪比为100%时：

参数 a1 a2 b1 b2 噪声方差 静态增益真值 -1.5 0.7 1 0.5 1 估计值

-1.519

0.72259

1.0314

0.50923

估计值

7.5661

源程序：

模型阶次

2

%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; end

S=ones(1,Np);%先产生一个全是1的序列S if mod(Np,2)==0%判断Np是奇数还是偶数 p=Np/2; else p=(Np-1)/2; end

for j=1:p

S(2*j)=0;%将S序列的偶数位值均置为0，从而使S序列是0或1的方波序列 end

IM=xor(M,S); %使用M序列与方波序列S复合生成逆M序列IM u=IM*2-1; for i=(Np+1):L u(i)=u(i-Np); end

randn('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); end

m=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); end

cta=zeros(4,L);

cta(:,1)=[0.001 0.001 0.001 0.001]'; P=diag([10^6 10^6 10^6 10^6]); %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 end end

for i=1:L

for j=1:nb if i>j

B(i,j)=u(i-j); end end end

H=[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); end end

[na nb]=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-') 第二大题 卡尔曼滤波 一个系统模型为

x1(k?1)?x1(k)?x2(k)?w(k),k?0,1...x2(k?1)?x2(k)?w(k)

同时有下列条件：

（1） 初始条件已知且有。x(0)?[0,0]T

（2） w(k)是一个标量零均值白高斯序列，且自相关函数已知为E[w(j)w(k)]??jk，

另外，我们有下列观测模型，即有下列条件：

（3） v1(k?1)和v2(k?1)是独立的零均值白高斯序列，且有

z1(k?1)?x1(k?1)?v1(k?1),k?0,1...z2(k?1)?x2(k?1)?v2(k?1)，且

E[v1(j)v1(k)]??jk，E[v2(j)v2(k)]?2?jk，k?0,1,2...

（4） 对于所有的j和k，w(k)与观测噪声过程v1(k?1)和v2(k?1)是不相关的，即

E[w(j)v1(k)]?0，E[w(j)v2(k)]?0

我们希望得到由观测矢量z(k?1)?[z1(k?1),z2(k?1)]T估计状态矢量

x(k?1)?[x1(k?1),x2(k?1)]T的卡尔曼滤波的公式表示，并求解以下问题：

（a） 求出卡尔曼增益矩阵，并得出最优估计x(k?1)和观测矢量

z(1),z(2),...,z(k?1)之间的递推关系。

?(k|k)，k?0,1,2...10，并画出当（b） 用模拟数据确定状态矢量x(k)的估计值x?1(k|k)，x?2(k|k)的图 k?0,1,2...10时，x（c） 通常，状态矢量的真实值是得不到的，但是为了用作图来说明问题，表1和表2

给出了状态矢量元素的真实值。对于k?0,1,2...10，在同一幅图中画出真实值和在（b）中确定的x1(k)的估计值。对x2(k)重复这一过程。当k从1变到10

?i(k|k)。 时，对每个元素i?1,2，计算并画出各自的误差图，即xi(k)?x（d） 当k从1变到10时，通过用由卡尔曼滤波器决定的状态误差协方差矩阵画出

2?i(k|k),i?1,2 E[?12(k|k)]和E[?2(k|k)]，而?i(k|k)?xi(k)?x（e） 讨论一下（C）中计算的误差与（d）中方差之间的关系。

表1观测值 时间下标k 1 2 3 4 5 6 7 8 9 10 表2由模拟得到的实际状态值 时间下标k 0 1 2 3 4 5 6 7 8 9 10 答： 实际状态值x1(k) 0.00000000 1.65428714 3.50300702 5.99785292 9.15040740 12.50873910 16.92192594 21.34483352 25.89335144 31.54135330 36.93605670 实际状态值x2(k) 0.00000000 1.65428714 1.84871988 2.47552222 3.17187816 3.35833170 4.41318684 4.42290758 4.54851792 5.64800186 5.39447034 观测值z1(k) 3．29691969 3．38736515 7．02830641 9．71212521 11．42018315 15．97980583 22．06934285 28．30212781 30．44683831 38．75875595 观测值z2(k) 2．10134294 0．47540797 3．17688898 2．49811140 2．91992424 6．17307616 5．42519274 3．05365741 5．98051141 4．51016361 ?CT?(C?Pk’?CT?R)?1 （a）卡尔曼增益矩阵：Hk?Pk’估计值与观测值之间的递归关系为：Xk?Ak?Xk-1?Hk?(Yk?Ck?Ak?Xk-1) （b）状态矢量估计值的计算框图：

????k?1 xyk?1 + ?k?1 xHk?1 + Z?1 Ck?1 Ak?1 ?'k?1 x

（c）x1(kk)和x2(kk)的图：

??

（d）真实值与估计值的比较图：

各自的误差图：

（e）通过用卡尔曼滤波器的状态误差协方差矩阵画出的E[?1(kk)]和

2E[?2(kk)]：

2

（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.29691969;2.10134294]; Z(:,:,3)=[3.38736515;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: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); end

gx1(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.34483352,25.89335144,31.54135330,36.93605670]; x2=[0,1.65428714,1.84871988,2.47552222,3.17187816,3.35833170,... 4.41318684,4.42290758,4.54851792,5.64800186,5.39447034]; 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;Tlta1,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);á2,e22分别表示x1,x2方差,P287 e22(1:11)=P(2,2,1:11); figure;

plot(1:10,e12(2:11),1:10,e22(2:11)); legend('x1的方差','x2的方差');

title('状态矢量X每个单独分量估计的方差');

本文来源：https://www.bwwdw.com/article/vvvx.html