过程辨识仿真作业.doc

例6.6仿真对象如图所示：

图中，为服从分布的不相关随机噪声；控制λ值，是数据的信噪比为23%。经过噪声模型后，迭加在输出y(k)上； 模型结构选用如下形式：

输入信号u(k)采用4阶的M序列，特征多项式取F(s)=1s，幅度为1, 循环周期为31bit，数据长度L=480, 初始条件=0.001, P(0)=I。利用辅助变量法辨识系统参数、。

一、 实验原理：
设SISO过程采用如下数学模型

描述，其中和表示过程的输入输出；是零均值的有色噪声；且

并且假定模型的阶次已知。由于是有色噪声，直接利用最小二乘法不能获得模型参数的无偏一致估计，这时可以利用辅助变量法，一起获得无偏一致估计。
1) 一次完成算法
将模型化为最小二乘格式

或
根据最小二乘知识，参数的最小估计值为

其中

由Frechet定理知

如果是白噪声，则，那么

如果不是白噪声，一般。为了使不是白噪声得情况下，仍有，定义一个辅助变量矩阵

=

使之满足以下两个条件：
①是非奇异阵；
②，即与独立。其中，称作辅助变量。
如果所用辅助变量能满足上述两个条件，则有

(*)
其中称作辅助变量参数估计值。可见只要适当选择辅助变量，使之满足以上两个条件，参数就可以是无偏一致估计。
2) 辅助变量的选择
选择辅助变量的基本原则就是上述两个条件必须满足。如果采用图a所示的作为辅助变量，置

图a 辅助变量的选择

当是持续激励信号时，必有是非奇异的，且因只与有关，即必与噪声无关；故有。这样辅助变量常有以下几种选择方法：
1． 自适应滤波法
把图a中辅助模型看成自适应滤波器时，辅助变量可按如下关系

或 (**)

取0.01~0.1；取0~10

求得。
2． 纯滞后法
当图a中的辅助模型为纯滞后环节时，辅助变量取作

其中，是多项式的阶次。此时，可写成
当是持续激励信号，且与噪声无关，则辅助变量可满足上述两个条件。
3． Tally原理
如果噪声可看成下列模型的输出

其中，是零均值的不相关噪声；且

则辅助变量取作

即
可满足上述两个条件，显然，只要输入信号是持续激励信号，条件①即可满足。另外，由于输入信号与噪声无关，故有

且

那么

则条件②也满足。
3) 递推算法
将(*)式一次完成算法写成

定义

则和推导最小二乘递推算法一样，辅助变量法的递推算式(简称RIV)可以写成

式中，为辅助向量，如果辅助变量选取(**)式，则需要用最小二乘法先递推若干步，以获取初步参数估计值，作为辅助变量法的递推初始状态。
二、 实验步骤
1） 输入信号的产生
输入信号采用4阶M序列，特征多项式取，幅度为1，循环周期为。
2） 噪声的产生

图b 仿真系统
如上图所示，由经过噪声模型后产生。要求为服从分布的不相关随机噪声；控制λ值，是数据的信噪比为23%。
3） 辅助变量的选择
使用Tally原理给出辅助变量，辅助变量取做，则辅助变量可写成：
其中，=2；=2；即

4） 利用辅助变量法进行参数估计
(1)初始化系统辨识初值：根据题目要求：=0.001, P(0)=I
(2)由于辅助最小二乘算法对初始值P(0)敏感，因此采用最小二乘算法（LS）进行辨识算法的起步，即用最小二乘辨识方法辨识前100步，然后用辅助最小二乘算法辨识。否则可能辨识没有可靠的收敛性；
(3)计算残差{ε(k)}的统计性质；
(5)计算阶跃响应。
三、 实验结果
根据实验要求产生的4阶M序列如下图所示：

图c 幅度为1的4阶M序列

图d RIV参数估计值的变化过程

图e ls参数估计值的变化过程

图f RIV阶跃响应比较

表1 RIV算法的辨识结果（噪信比η=23%）
参 数
a1
a2
b1
b2
静态增益
噪声均值
噪声方差
真 值
-1.5
0.7
1.0
0.5
7.5
0.0
0.23
估计值
-1.5143
0.7100
0.9064
0.4975
7.5332
2.24E-2
0.2738
表2 阶跃响应比较
k
模型阶跃响应
过程阶跃响应
k
模型阶跃响应
过程阶跃响应
0
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
0
0.0553 0.0543 1.2950 3.1685 5.4466 7.4350
8.9429 9.2900 9.4649 8.9030 8.5585 7.7741
7.0527 6.5679 6.4264
0
0.0553 0.0543 0.7868 2.6600 4.6122 6.7288
8.1230 8.9962 8.4569 8.4759 7.5047 7.3810
6.4403 5.9047 5.6819
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
6.4046 6.4181 6.5413
6.6312 6.5494 6.5866 7.3975 7.7134 7.7220
7.4495 7.3840 6.8860 6.6050 6.7591 6.2673
5.8103 5.8771 5.9129
6.0892 6.1381 5.9512 6.0651 7.2625 7.1674
6.9576 6.5403 6.6336 5.9283 5.8554 6.2861

表3 输出残差序列{ε(k)}的统计性质
均值{ε(k)}=-2.29E-2
ρ(0)
1.0

ρ(1)
0.2621
ρ(6)
-0.1803
ρ(11)
0.3579
ρ(16)
0.3454
ρ(2)
-0.4282
ρ(7)
-0.6902
ρ(12)
0.3255
ρ(17)
-0.4567
ρ(3)
-0.5079
ρ(8)
-0.7099
ρ(13)
0.7054
ρ(18)
-0.4462
ρ(4)
-0.3423
ρ(9)
-0.9590
ρ(14)
1.2514
ρ(19)
-0.3574
ρ(5)
-0.1340
ρ(10)
-0.3243
ρ(15)
0
ρ(20)
-0.0390
四、 实验结果分析
由于本题的输入是有色噪声，因此用最小二乘辨识方法一般得不到无偏一致估计。这里采用辅助变量法，递推计算480步后的辨识结果如表1所示。参数估计的变化过程如图d所示，模型与过程的阶跃响应比较吻合，如图f与表2所示，辨识结果是令人满意的。为进一步确认辨识结果，需要对所获得的模型进行检验。计算输出残差{ε(k)}的均值和自相关系数，结果如表3所示。表3表明，由于数值比较大，可以确信输出残差{ε(k)}序列是有色噪声，而输入v(k)就是有色噪声，因此所获得的模型是可靠的。仿真实验表明，辅助变量法的计算量与最小二乘相当，但辨识结果却比最小二乘法好的多。因此辅助变量法是一种很有价值的辨识方法，尤其当噪声是有色的，而噪声的模型又不好确定时，辅助变量法就显示出了它的优势。但辅助变量法不能像增广最小二乘法或广义最小二乘法那样可以同时获得噪声模型的参数估计。
附录 实验源程序
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);
end
end
r=a(:,1);
r1=[r;r;r;r;r;r;r;r];
r2=[r1;r1;r1;r1];
u=zeros(1,480);
for i=1:480
u(i)=1*(1-2*r2(i));% 生成4阶M序列
end
figure(1)
stairs(u);
axis([1 480 -1.5 1.5])
grid on;
p=10^6*eye(4);
theta=[0.001 0.001 0.001 0.001]';
I=eye(4);
a1=zeros(1,481);
a2=a1;
b1=a1;
b2=a1;
z=a1;
y=a1;
rou=zeros(1,20);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%¸辅助变量%%%%%%%%%%%%%%%%%%%
f=[0 0 0 0]';
A=[1 -1.5 0.7];
B=[0 1.0 0.5];
C=[1 -1 0.2];
E=zeros(480,1);
T=zeros(480,4);
T(1,:)=theta';
v=normrnd(0,1,480,1);%生成480均值，方差为1的正态分布噪声
V=0.23*v;
e=filter(C,A,V);%生成有色噪声
z(1)=v(1);
z(2)=v(2);
%计算输出
for k=3:1:480
z(k)=1.5*z(k-1)-0.7*z(k-2)+1.0*u(k-1)+0.5*u(k-2)+e(k);
end
%辅助变量法进行辨识
for i=2:1:40
theta=theta+p*h*((h'*p*h+1)^(-1))*(z(i)-h'*theta);
p=(I-p*h*(h'*p*h+1)^(-1)*h')*p;
h=[-z(i) -z(i-1) u(i) u(i-1)]';
T(i,:)=theta';
E(i)=z(i)-h'*T(i,:)';
end
for i=41:1:480
theta=theta+p*f*((h'*p*f+1)^(-1))*(z(i)-h'*theta);
p=(I-p*f*(h'*p*f+1)^(-1)*h')*p;
h=[-z(i) -z(i-1) u(i) u(i-1)]';
f=[-z(i-3) -z(i-4) u(i-1) u(i-2)]'; %Tally原理
T(i,:)=theta';
E(i)=z(i)-h'*T(i,:)';
end
k=1:1:480;
figure(2);
plot(k,T(:,1),'b');
axis([1 540 -1.6 1.6])
hold on;
grid on;
k=1:1:480;
plot(k,T(:,2),'r');
hold on;
grid on;
k=1:1:480;
plot(k,T(:,3),'g');
hold on;
grid on;
k=1:1:480;
plot(k,T(:,4),'k');
hold on;
grid on;
title('RIVµÄ±æʶ½á¹û');
text(490,1.0,'\leftarrow b1');
text(490,0.7,'\leftarrow a2');
text(490,0.5,'\leftarrow b2');
text(490,-1.5,'\leftarrow a1');
xlabel('k');
ylabel('theta');
%输出残差的计算
a=mean(E)
for i=1:1:480
A0=E(i)*E(i);
end
R0=mean(A0);
for i=1:1:479
A1=E(i)*E(i+1);
end
R(1)=mean(A1);
for i=1:1:478
A2=E(i)*E(i+2);
end
R(2)=mean(A2);
for i=1:1:477
A3=E(i)*E(i+3);
end
R(3)=mean(A3);
for i=1:1:476
A4=E(i)*E(i+4);
end
R(4)=mean(A4);
for i=1:1:475
A5=E(i)*E(i+5);
end
R(5)=mean(A5);
for i=1:1:474
A6=E(i)*E(i+6);
end
R(6)=mean(A6);
for i=1:1:473
A7=E(i)*E(i+7);
end
R(7)=mean(A7);
for i=1:1:472
A8=E(i)*E(i+8);
end
R(8)=mean(A8);
for i=1:1:471
A9=E(i)*E(i+9);
end
R(9)=mean(A9);
for i=1:1:470
A10=E(i)*E(i+10);
end
R(10)=mean(A10);
for i=1:1:469
A11=E(i)*E(i+11);
end
R(11)=mean(A11);
for i=1:1:468
A12=E(i)*E(i+12);
end
R(12)=mean(A12);
for i=1:1:467
A13=E(i)*E(i+13);
end
R(13)=mean(A13);
for i=1:1:466
A14=E(i)*E(i+14);
end
R(14)=mean(A14);
for i=1:1:465
A15=E(i)*E(i+15);
end
R(1)=mean(A1);
for i=1:1:464
A16=E(i)*E(i+16);
end
R(16)=mean(A16);
for i=1:1:463
A17=E(i)*E(i+17);
end
R(17)=mean(A17);
for i=1:1:462
A18=E(i)*E(i+18);
end
R(18)=mean(A18);
for i=1:1:461
A19=E(i)*E(i+19);
end
R(19)=mean(A19);
for i=1:1:460
A20=E(i)*E(i+20);
end
R(20)=mean(A20);
for i=1:1:20
rou(i)=R(i)/R0;
end
rou
%¸阶跃响应
Gain1=[];
Gain=[];
B=[1.0 -1 0.2];
A=[1 -1.5 0.7];
u=ones(1,30);
v=normrnd(0,1,30,1);
v=0.23*v;
e=filter(A,B,v);
junzhi=mean(e)%求噪声均值
fangcha=std(e)%求噪声方差
y=zeros(1,30);
z=zeros(1,30);
y(1)=v(1);
z(1)=v(1);
y(2)=v(2);
z(2)=v(2);
for i=3:1:30
y(i)=1.5*y(i-1)-0.7*y(i-2)+1.0*u(i-1)+0.5*u(i-2)+e(i);
end
Gain1=1.5*0.5*10;
for k=3:1:30
z(k)=-T(480,1)*y(k-1)-T(480,2)*y(k-2)+T(480,3)*u(k-1)+T(480,4)*u(k-2)+e(i);
Gain=-T(480,4)*T(480,1)*10;
end
figure(3);
plot(y,'*r');
hold on;
plot(z,'.b');
legend('¹À¼ÆÖµ','ÀíÂÛÖµ');
ylabel('h(k)');

《系统辨识》大作业

学院：电信学院
专业：检测技术与自动化装置
学号：102081102003
姓名：丑永新