过程辨识大作业

考虑如图所示的**对象。

图1 **对象。

图中，v(k)为服从n(0，1)分布的不相关随机噪声；经过噪声模型后，迭加在输出y(k)上；输出信号u(k)采用4阶m序列，特征多项式取f(s)=1⊕s⊕，幅度为1，循环周期np=31bit。选用如下的模型结构。

z(k)+a1z(k-1)+a2z(k-2)=b1u(k-1)+b2u(k-2)+e(k)

数据长度l=400，初始条件(0)=0.001，p(0)=106i。利用辅助变量法辨识结果。

1 递推算法。

辨识采用递推辅助变量算法，与推导最小二乘算法一样，辅助变量算法的递推算式（简称riv）可写成。

k)=(k-1)+k(k)[z(k)-(k)(k-1)]

k(k)=p(k-1)(k)[(k)p(k-1)(k)+1]-1

p(k)=[i-k(k)(k)]p(k-1)

2辅助变量的选择。

使用纯滞后的方法给出辅助变量，辅助变量取作x(k)=u(k-nb)，其中nb是多项式b(z-1)的阶次。此题nb=2。这时辅助变量(k)可写成。

k)=[u(k-nb-1)，…u(k-nb- na)，u(k-1)，…u(k- nb)]

[-u(k-3)，-u(k-4)，u(k-1)，u(k- 2)]

3 辨识步骤。

1）直接取辨识初值(0)=0.001，p(0)= 106i；

2）由于辅助最小二乘算法对初始值p(0)敏感，因此采用最小二乘算法（ls）进行辨识算法的起步，即用最小二乘辨识方法辨识前100步，然后用辅助最小二乘算法辨识。否则可能辨识没有可靠的收敛性；

3）计算残差的统计性质；

4）计算阶跃响应。

4 **结果。

由于本题的输入是有色噪声，因此用最小二乘辨识方法一般得不到无偏一致估计。这里采用辅助变量法，递推计算400步后的辨识结果如表1所示。参数估计的变化过程如图2所示，模型与过程的阶跃响应比较吻合，如图3与表2所示，辨识结果是令人满意的。

为进一步确认辨识结果，需要对所获得的模型进行检验。计算输出残差的均值和自相关系数，结果如表3所示。表3表明，由于数值比较大，可以确信输出残差序列是有色噪声，而输入v(k)就是有色噪声，因此所获得的模型是可靠的。

**实验表明，辅助变量法的计算量与最小二乘相当，但辨识结果却比最小二乘法好的多。因此辅助变量法是一种很有价值的辨识方法，尤其当噪声是有色的，而噪声的模型又不好确定时，辅助变量法就显示出了它的优势。但辅助变量法不能像增广最小二乘法或广义最小二乘法那样可以同时获得噪声模型的参数估计。

表1 riv算法的辨识结果（噪信比η=25%）

图2 参数估计值的变化过程图3 阶跃响应比较。

表2 阶跃响应比较。

表3 输出残差序列的统计性质。

附录1 辨识程序。

n=4;%生成4阶m序列。

a=zeros(2^n-1,n);

a(1,n)=1;

for i=2:1:2^n-1

a(i,1)=mod(a(i-1,n-1)+a(i-1,n),2);

for p=2:1:n

a(i,p)=a(i-1,p-1);

endend

r=a(:,1);

r1=[r;r;r;r;r;r];

r2=[r1;r1;r1;r1;r1];

u=zeros(1,400);

for i=1:400

u(i)=1*(1-2*r2(i));

endp=10^6*eye(4);

theta=[0.001 0.001 0.001 0.001]';

i=eye(4);

a1=zeros(1,401);

a2=a1;

b1=a1;

b2=a1;

z=a1;y=a1;

rou=zeros(1,20);

h=[-y(1) 0 u(1) 0]';

f=[0 0 0 0]';

a=[1 -1.5 0.7];

b=[1.0 0.5];

e=zeros(400,1);

t=zeros(400,4);

t(1,:)theta';

e=normrnd(0,0.25,400,1);%生成均值为0，方差为0.25的白噪声。

v=filter(2,1,e);%生成有色噪声。

y(1)=v(1);

y(2)=v(2);

for k=3:1:400

y(k)=1.5*y(k-1)-0.7*y(k-2)+1.0*u(k-1)+0.5*u(k-2)+v(k);

endtheta=theta+(p*h*((h'*p*h+1)^(1)))y(1)-h'*theta);

p=(i-p*h*(h'*p*h+1)^(1)*h')*p;

h=[-y(2) -y(1) u(2) u(1)]'

t(1,:)theta';

e(1)=y(1)-h'*t(1,:)

for i=2:1:100%最小二乘辨识方法辨识前100步。

theta=theta+p*h*((h'*p*h+1)^(1))*y(i)-h'*theta);

p=(i-p*h*(h'*p*h+1)^(1)*h')*p;

h=[-y(i) -y(i-1) u(i) u(i-1)]'

t(i,:)theta';

e(i)=y(i)-h'*t(i,:)

endfor i=101:1:400

theta=theta+p*f*((h'*p*f+1)^(1))*y(i)-h'*theta);%辅助变量法辨识后300步。

p=(i-p*f*(h'*p*f+1)^(1)*h')*p;

h=[-y(i) -y(i-1) u(i) u(i-1)]'

f=[-y(i-3) -y(i-4) u(i) u(i-1)]'

t(i,:)theta';

e(i)=y(i)-h'*t(i,:)

endk=1:1:400;

plot(k,t(:,1),'b');

hold on;

grid on;

k=1:1:400;

plot(k,t(:,2),'r');

hold on;

grid on;

k=1:1:400;

plot(k,t(:,3),'g');

hold on;

grid on;

k=1:1:400;

plot(k,t(:,4),'y');

hold on;

grid on;

title('riv的辨识结果');

text(380,1.2,'b1');

text(380,0.8,'a2');

text(380,0.4,'b2');

text(380,-1.4,'a1');

xlabel('k');

ylabel('theta');

附录2 输出残差序列的统计性质。

a=mean(e)

for i=1:1:400

a0=e(i)*e(i);

endr0=mean(a0);

for i=1:1:399

a1=e(i)*e(i+1);

endr(1)=mean(a1);

for i=1:1:398

a2=e(i)*e(i+2);

endr(2)=mean(a2);

for i=1:1:397

a3=e(i)*e(i+3);

endr(3)=mean(a3);

for i=1:1:396

a4=e(i)*e(i+4);

endr(4)=mean(a4);

for i=1:1:395

a5=e(i)*e(i+5);

endr(5)=mean(a5);

for i=1:1:394

a6=e(i)*e(i+6);

endr(6)=mean(a6);

for i=1:1:393

a7=e(i)*e(i+7);

endr(7)=mean(a7);

for i=1:1:392

a8=e(i)*e(i+8);

endr(8)=mean(a8);

for i=1:1:391

a9=e(i)*e(i+9);

endr(9)=mean(a9);

for i=1:1:390

a10=e(i)*e(i+10);

endr(10)=mean(a10);

for i=1:1:389

a11=e(i)*e(i+11);

endr(11)=mean(a11);

for i=1:1:388

a12=e(i)*e(i+12);

endr(12)=mean(a12);

过程辨识大作业

系统辨识大作业

系统辨识大作业

系统辨识大作业

其他用户还读了