实验一 辩识离散线性系统脉冲响应特性的相关分析法
实验目的
(1) 掌握利用相关分析法辩识离散线性系统冲激响应特性的基本步骤和使用要点;
(2) 了解递推和批量两种算法如何具体编程;
(3) 熟悉估计算法中各种参数选择对辩识效果的影响。
实验原理
待辩识系统的基本情况和辩识工作的要求
待辩识系统的运行情况如图 1 所示,设其在静态工作点(U0,Y0)附近作局部线性化
所得动态模型为
y(k)+a1y(k-1)+a2y(k-2)=b1u(k-1-d)+b2u(k-2-d)
(1)
式中参数建议值为
a
1
,1
a
2
,5.0
b
1
,1
b
2
5.0
。
U0
U0+ u(k)
u(k)
待辩识系统
Y(k)=Y0+y(k)
η(k)
Y’(k)=Y0+y(k) +η(k)
图 1 待辩识系统的基本运行情况
为了进行系统冲激响应特性的辩识,在对象输入端施加电平幅值 V=1 的 M 序列信号(其
波形见图 2)作为动态测试信号 u(k)。本次实验的基本任务就是在输入、输出的实际观测值
Y’(k)=Y0+y(k) +η(k) 和 u(k)的基础上,用批量或递推算法估计出脉冲响应序列{h(k)}。
u(k)
0
1
2
3
4
5
6
7
8
9
10 11 12 13 14 15
图 2 M 序列信号的波形
M 序列测试信号的生成方法
M 序列又称最大长度伪随机二进制序列,它实际上是一种特殊的周期序列,其电平
值 u(k)只取+V 和-V (或者说其逻辑电平值
(ku 只取 1 或 0)。
)
u
(
k
)
u
(
kk
)
0
u
(
nk
)
M 序列的当前值可以由过去的值进行异或运算得到。选取恰当的 k0,便可以产生
最长周期为 Np=2n-1 的 M 序列。很容易由计算机程序来产生 M 序列。如取 n=4,
k0=3,
u
)0(
u
)1(
u
)2(
u
)3(
1
,便可按
u
(
k
)
u
(
k
)3
u
(
k
)4
依次推出
u
,
u
)4(
)5(
,
,最后得出图 2 所示的 M 序列。
M 序列的循环长度为 NP=2 N -1=15,由此得到 N=4.
M=[1 1 1 1 –1 –1 –1 1 –1 –1 1 1 –1 1 -1],寄存器初始状态为(1111)。
基本原理和基本关系式
在相关法中,采用伪随机信号 M 序列对系统模型进行在线辩识。由于 M 序列的自
相关函数很接近于周期性的δ函数,可近似地当作白噪声信号叠加在系统输入信号上,
按其与输出信号间的互相关函数,可得到系统的脉冲响应曲线。
设对象输入 u 为具有平稳遍历性的随机信号,其自相关函数为
)(R
E[u(t)u(t
u
)]
u(t)u(t
)(r
)
u
当系统初始过渡过程结束后,输入和输出之间的互相关函数为
)(R
E[u(t)y(t
uy
)]
u(t)y(t
)
r
uy
)(
式中 E[●]表示统计平均, ● 表示时间平均。由维纳 – 霍甫方程有
r
uy
)(
0
h(
)r
(
u
)
dv
1
T
p
pT
0
h(
)r
(
u
)
dv
当 u 为伪随机信号 M 序列时,自相关函数的离散值为
r
u
)(
l
r
u
|)
(
l
l
,0
N
2
V
,
V
N
2
p
,
2,
N
,
p
p
otherwise
(2)
(3)
(4)
(5)
式中,V、△和 Np 分别为 M 序列的电平值、基本电平时间和离散周期长度,ι为时差的
离散值。
相应的功率谱密度函数为
S
u
(
)
j
ˆ
S
0
i
2
sin(
2
2
)
0
(
i
)
(6)
式中
1
1
pN
ˆ
S
0
V
2
若将 M 序列的自相关函数近似为
r
uy
)(
ˆ
S
0
i
h
(
iT
)
p
并取的周期大于对象的整定时间,即
T
p
N
p
T
s
则有
r
uy
)(
ˆ
hS
0
),
(
0
T
p
这样就可以得到脉冲响应曲线
h
)(
r
uy
)(
ˆ
S
0
,
0
T
p
批量算法和递推算法
1) 批量算法(Batch algorithm):
r
uy
'
)(
l
r
uy
'
|)(
l
1
N
p
1
(
ku
N
p
k
0
)(')
kyl
h
)(
N
p
p
2
)1
V
(
N
r
uy
'
)(
r
uy
'
)0(
2) 递推算法(Recursive algorithm):
ˆ
h
)0(
ˆ
h
)1(
1
m
ˆ
h
m
ˆ
h
m
ˆ
h
(
)1
m
1
m
m
1
N
p
N
p
(
2
N
p
)1
V
u
)
u
(
m
u
)1
(
m
(
Nm
)1
p
y
'
(
m
)
ˆ
h
m
1
(7)
(8)
(9)
(10)
(11)
(12)
(13)
(14)
辩识对象运行情况的仿真
取 系 统 对 象 模 型
y
(
k
)
ya
1
(
k
)1
ya
2
(
k
)2
ub
(1
k
1
d
)
ub
2
(
k
2
d
)
中 初 值 为
y
)1(
y
)0(
0
。将 M 序列作为 )
(ku ,即可求
y k
(
,)
k
,3,2,1
。设工作点(U0,Y0)
分别为 0.1 和 0.2,则对象的实际输入、输出 U(k)和 Y’(k)分别为
U
(
k
)
U
0
u
(
k
)
Y
('
k
)
Y
0
y
(
k
)
(
k
)
(18)
其中正态白噪声序列η(k)根据中心极限定理采用 12 个在[-0.5,0,5]中均匀分布的
随机数叠加的方法生成,其方差可由遭受幅度因子调节。
在待辩识对象在大范围内是非线性的情况下,通过辩识得出的是对象在工作
点附近的局部线性化动态模型,因此还应去处工作点 Y0 ,以得出输出端动态分量
的 实 测 值 。 但 工 作 点 无 法 单 独 测 量 , 只 能 根 据 最 新 得 到 的
N
N
p
(
)4~3
个
Y
('
k
)
推算:
ˆ
Y
(0
k
)
1
N
或递推计算如下:
1
N
i
0
Y
'
(
k
)1
ˆ
Y
(0
k
)
ˆ
Y
(0
k
)1
1
N
Y
'
Y
'
(
k
)
(
Nk
)
式中 Y0(k)的初值 Y0(0)可设为 0。于是,输出端动态分量的实测值可近似估算为
实验步骤
批量算法如图 3:
y
'
(
k
)
Y
'
(
k
)
ˆ
Y
(
ko
)
初始化
(19)
(20)
(21)
生成数据 u’(k-d)
, Y’(k) , k=1-r-d,…,-1,0,…,Np-1
k1 ← Np + r + d
生成新的数据
u
(
k
)
,
u
(
k
)
,
Y
(
k
)
k ← k +1
k-k1>iTp?
r
uy
'
)(
l
r
uy
'
|)
(
l
1
N
p
h
(
)
N
p
p
)1
V
2
(
N
p
k
r
uy
0
N
1
(
ku
)
yl
('
k
)
(
)
'
r
uy
'
)0(
打印结果
图 3 批量算法流程图
第一步:产生 M 序列,如图所示
NP=15,N=4
第二步:求系统输出,编程结果如图所示
y(k)+a1y(k-1)+a2y(k-2)=b1u(k-1-d)+b2u(k-2-d)
a
,5.0
5.0
。
,1
,1
a
2
1
b
1
b
2
第三步:利用相关法求脉冲响应,如图所示
r
uy
'
)(
l
r
uy
'
|)(
l
1
N
p
1
(
ku
N
p
k
0
)(')
kyl
h
)(
N
p
p
2
)1
V
(
N
r
uy
'
)(
r
uy
'
)0(
递推算法如图 3:
初始化
生成数据 u’(k-d)
, Y’(k) , k=1-r-d,…,-1,0,…,Np-1
k1 ← Np + r + d
, 构造 h,
生成新的数据
向量,并送入 h1 向量,等待循环计算
k
)
)
(
k
)
u
(
k
,
Y
(
,
u
(ku
)
k ← k +1
ˆ
h
m
ˆ
h
(
m
m
ˆ
h
)0(
ˆ
h
)1(
N
p
)1
m
u
)
(
m
u
(
)1
m
u
(
Nm
)1
p
'
y
(
)
m
ˆ
h
1
m
ˆ
h
1
m
1
m
1
N
p
(
2
N
p
)1
V
k-k1>iTp?
打印结果
第一、 二步同前
第三步:脉冲响应如图所示
ˆ
h
m
m
m
ˆ
h
)0(
ˆ
h
)1(
ˆ
h
(
N
p
)1
m
ˆ
h
m
1
1
m
1
N
p
(
2
N
p
)1
V
u
)
u
(
m
u
)1
(
m
(
Nm
)1
p
y
'
(
m
)
ˆ
h
m
1
如下图所示,相关分析法很好地辨识了系统的脉冲响应
for t=1:15
Z(t,1)=Y(1,t);
end
g=H*R*M*Z
subplot(3,1,3)
k=1:1:i1;
plot(k,g,k,g,'rx')
xlabel('k')
ylabel('输出')
title('脉冲响应估计值')
for r=1:15
z(r)=m(r)*g(r);
end
%subplot(4,1,4)
%k=1:1:i1;
%plot(k,z,k,z,'rx')
%xlabel('k')
%ylabel('输出')
%title('辨识输出')
附录程序:
(1) 批量算法
x1=1;x2=1;x3=1;x4=1;
m=15;
for i=1:m
y4=x4;y3=x3;y2=x2;y1=x1;
x4=y3;x3=y2;x2=y1;x1=xor(y3,y4);
if y4==0
u(i)=-1;
u(i)=y4;
else
end
end
m=u
%grapher
i1=i;
k=1:1:i1;
subplot(3,1,1)
plot(k,u,k,u,'rx')
xlabel('k')
ylabel('m 序列')
title('移位寄存器产生的 m 序列')
a1=-1;a2=0.5;b1=1;b2=0.5;Y(1)=0;Y(2)=0;
for K=3:1:16
Y(K)=Y(K-1)-0.5*Y(K-2)+u(K-1)+0.5*u(K-2)
;
end
subplot(3,1,2)
k=1:1:i1+1;
plot(k,Y,k,Y,'rx')
xlabel('k')
ylabel('输出')
title('系统输出')
H=1/16;
p=ones(15);q=eye(15);
R=p+q;
for j=1:15
for i=1:16-j
M(i,i+j-1)=u(j);
end
end