v2=
-0.9531-0.71880.6875-0.8359-0.01560.9219
0.57030.4531-0.2500-0.48440.1016-0.3672
0.8047-0.13280.21880.3359-0.9531-0.7188
0.6875-0.8359-0.01560.92190.57030.4531
-0.2500-0.48440.1016-0.36720.8047-0.1328
0.21880.3359-0.9531-0.71880.6875-0.8359
-0.01560.92190.57030.4531-0.2500-0.4844
0.1016-0.36720.8047-0.13280.21880.3359
-0.9531-0.71880.6875-0.8359-0.01560.9219
0.57030.4531-0.2500-0.48440.1016-0.3672
0.8047-0.13280.21880.3359-0.9531-0.7188
0.6875-0.8359-0.01560.92190.57030.4531
-0.2500-0.48440.1016-0.36720.8047-0.1328
0.21880.3359-0.9531-0.71880.6875-0.8359
-0.01560.92190.57030.4531-0.2500-0.4844
0.1016-0.36720.8047-0.13280.21880.3359
-0.9531-0.71880.6875-0.8359
8.2M序列的产生
在进行系统辨识时,选用白噪声作为辨识输入信号可以保证获得较好的辨识效果,但在项目上难以实现。
M序列是一种很好的辨识输入信号,它具有近似白噪声的性质,不仅可以保证有较好的辨识效果,而且项目上又易于实现。
8.2.1M序列的产生方式
M序列是一种离散二位式随机序列,所谓“二位式”是指每个随机变量只有两种状态。
离散二位式随机序列是按照确定的方式产生的,实际上是一种确定序列。
可用多级线性反馈移位寄存器产生M序列。
每级移位寄存器由双稳态触发器和门电路组成,称为1位,分别以0和1来表示2中状态。
当移位脉冲来到时,每位的内容<0或1)移到下一位,最后1位<即第n位)移出的内容即为输出。
为了保持连续工作,将最后2级寄存器的内容经过适当的逻辑运算后反馈到第1级寄存器作为输入。
8.2.2用四级移位寄存器产生M序列
用移位寄存器产生M序列的MATLAB软件实现(见配套MATLAB程序algorithms2.m>
①编程如下:
X1=1。
X2=0。
X3=1。
X4=0。
%移位寄存器输入Xi初T态<0101),Yi为移位寄存器各级输出
m=60。
%置M序列总长度
fori=1:
m%1#
Y4=X4。
Y3=X3。
Y2=X2。
Y1=X1。
X4=Y3。
X3=Y2。
X2=Y1。
X1=xor(Y3,Y4>。
%异或运算
ifY4==0
U(i>=-1。
else
U(i>=Y4。
end
end
M=U
%绘图
i1=i
k=1:
1:
i1。
plot(k,U,k,U,'rx'>
xlabel('k'>
ylabel('M序列'>
title('移位寄存器产生的M序列'>
②程序运行结果如图8.2所示。
图8.2软件实现的移位寄存器产生的M序列图
③四级移位寄存器产生的M序列
M=
Columns1through10
-11-11111-1-1-1
Columns11through20
1-1-111-11-111
Columns21through30
11-1-1-11-1-111
Columns31through40
-11-11111-1-1-1
Columns41through50
1-1-111-11-111
Columns51through60
11-1-1-11-1-111
i1=
60
8.3最小二乘一次完成算法的产生
考虑仿真对象
(8.4>其中,
是服从正态分布的白噪声N(0,1>。
输入信号采用4阶M序列,幅度为1。
选择如下形式的辨识模型
(8.5>
设输入信号的取值是从k=1到k=16的M序列,则待辨识参数
为
=
。
其中,被辨识参数
、观测矩阵zL、HL的表达式为
,
(8.6>
程序框图如图8.3所示。
Matlab仿真程序如下:
%二阶系统的最小二乘一次完成算法辨识程序,<见配套MATLAB程序algorithms3.m>
u=[-1,1,-1,1,1,1,1,-1,-1,-1,1,-1,-1,1,1]。
%系统辨识的输入信号为一个周期的M序列
z=zeros(1,16>。
%定义输出观测值的长度
fork=3:
16
z(k>=1.5*z(k-1>-0.7*z(k-2>+u(k-1>+0.5*u(k-2>。
%用理想输出值作为观测值
end
subplot(3,1,1>%画三行一列图形窗口中的第一个图形
stem(u>%画输入信号u的径线图形
subplot(3,1,2>%画三行一列图形窗口中的第二个图形
i=1:
1:
16。
%横坐标范围是1到16,步长为1
plot(i,z>%图形的横坐标是采样时刻i,纵坐标是输出观测值z,图形格式为连续曲线
subplot(3,1,3>%画三行一列图形窗口中的第三个图形
stem(z>,gridon%画出输出观测值z的径线图形,并显示坐标网格
u,z%显示输入信号和输出观测信号
%L=14%数据长度
HL=[-z(2>-z(1>u(2>u(1>。
-z(3>-z(2>u(3>u(2>。
-z(4>-z(3>u(4>u(3>。
-z(5>-z(4>u(5>u(4>。
-z(6>-z(5>u(6>u(5>。
-z(7>-z(6>u(7>u(6>。
-z(8>-z(7>u(8>u(7>。
-z(9>-z(8>u(9>u(8>。
-z(10>-z(9>u(10>u(9>。
-z(11>-z(10>u(11>u(10>。
-z(12>-z(11>u(12>u(11>。
-z(13>-z(12>u(13>u(12>。
-z(14>-z(13>u(14>u(13>。
-z(15>-z(14>u(15>u(14>]%给样本矩阵HL赋值
ZL=[z(3>。
z(4>。
z(5>。
z(6>。
z(7>。
z(8>。
z(9>。
z(10>。
z(11>。
z(12>。
z(13>。
z(14>。
z(15>。
z(16>]%给样本矩阵zL赋值
%CalculatingParameters
c1=HL'*HL。
c2=inv(c1>。
c3=HL'*ZL。
c=c2*c3%计算并显示
%DisplayParameters
a1=c(1>,a2=c(2>,b1=c(3>,b2=c(4>%从
中分离出并显示a1、a2、b1、b2
%End
程序运行结果:
>>
u=[-1,1,-1,1,1,1,1,-1,-1,-1,1,-1,-1,1,1]
z=[0,0,0.5000,0.2500,0.5250,2.1125,4.3012,6.4731,6.1988,3.2670,-0.9386,-3.1949,-4.6352,6.2165,-5.5800,-2.5185]
HL=
001.0000-1.0000
-0.50000-1.00001.0000
-0.2500-0.50001.0000-1.0000
-0.5250-0.25001.00001.0000
-2.1125-0.52501.00001.0000
-4.3012-2.11251.00001.0000
-6.4731-4.3012-1.00001.0000
-6.1988-6.4731-1.0000-1.0000
-3.2670-6.1988-1.0000-1.0000
0.9386-3.26701.0000-1.0000
3.19490.9386-1.00001.0000
4.63523.1949-1.0000-1.0000
6.21654.63521.0000-1.0000
5.58006.21651.00001.0000
ZL=[0.5000,0.2500,0.5250,2.1125,4.3012,6.4731,6.1988,3.2670,-0.9386,-3.1949,-4.6352,-6.2165,-5.5800,-2.5185]T
c=[-1.5000,0.7000,1.0000,0.5000]T
a1=-1.5000
a2=0.7000
b1=1.0000
b2=0.5000
>>
从仿真结果表8.4可以看出,因为所用的输出观测值没有任何噪声成分,所以辨识结果也无任何误差。
8.4最小二乘递推算法的产生
考虑图8.5所示的仿真对象,
图中,
是服从N<0,1)分布的不相关随机噪声。
且
,
,
选择图8.5所示的辨识模型。
仿真对象选择如下的模型结构
(8.7>
其中,
是服从正态分布的白噪声N<0,1)。
输入信号采用4位移位寄存器产生的M序列,幅度为0.03。
按式
(8.8>
构造h(k>;加权阵取单位阵
;计算K(k>、和P(k>,计算各次参数辨识的相对误差,精度满足要求后停机。
最小二乘递推算法辨识的Malab程序流程如图8.6所示。
图8.6最小二乘递推算法辨识的Malab6.0程序流程图
下面给出具体程序。
%最小二乘递推算法辨识程序,<见配套MATLAB程序algorithms4.m>
clear%清理工作间变量
L=15。
%M序列的周期
y1=1。
y2=1。
y3=1。
y4=0。
%四个移位寄存器的输出初始值
fori=1:
L。
%开始循环,长度为L
x1=xor(y3,y4>。
%第一个移位寄存器的输入是第三个与第四个移位寄存器的输出的“或”
x2=y1。
%第二个移位寄存器的输入是第一个移位寄存器的输出
x3=y2。
%第三个移位寄存器的输入是第二个移位寄存器的输出
x4=y3。
%第四个移位寄存器的输入是第三个移位寄存器的输出
y(i>=y4。
%取出第四个移位寄存器的幅值为"0"和"1"的输出信号,即M序列
ify(i>>0.5,u(i>=-0.03。
%如果M序列的值为"1",辨识的输入信号取“-0.03”
elseu(i>=0.03。
%如果M序列的值为"0",辨识的输入信号取“0.03”
end%小循环结束
y1=x1。
y2=x2。
y3=x3。
y4=x4。
%为下一次的输入信号做准备
end%大循环结束,产生输入信号u
figure(1>。
%第一个图形
stem(u>,gridon%显示出输入信号径线图并给图形加上网格
z(2>=0。
z(1>=0。
%设z的前两个初始值为零
fork=3:
15。
%循环变量从3到15
z(k>=1.5*z(k-1>-0.7*z(k-2>+u(k-1>+0.5*u(k-2>。
%输出采样信号
end
%RLS递推最小二乘辨识
c0=[0.0010.0010.0010.001]'。
%直接给出被辨识参数的初始值,即一个充分小的实向量
p0=10^6*eye(4,4>。
%直接给出初始状态P0,即一个充分大的实数单位矩阵
E=0.000000005。
%取相对误差E=0.000000005
c=[c0,zeros(4,14>]。
%被辨识参数矩阵的初始值及大小
e=zeros(4,15>。
%相对误差的初始值及大小
fork=3:
15。
%开始求K
h1=[-z(k-1>,-z(k-2>,u(k-1>,u(k-2>]'。
x=h1'*p0*h1+1。
x1=inv(x>。
%开始求K(k>
k1=p0*h1*x1。
%求出K的值
d1=z(k>-h1'*c0。
c1=c0+k1*d1。
%求被辨识参数c
e1=c1-c0。
%求参数当前值与上一次的值的差值
e2=e1./c0。
%求参数的相对变化
e(:
k>=e2。
%把当前相对变化的列向量加入误差矩阵的最后一列
c0=c1。
%新获得的参数作为下一次递推的旧参数
c(:
k>=c1。
%把辨识参数c列向量加入辨识参数矩阵的最后一列
p1=p0-k1*k1'*[h1'*p0*h1+1]。
%求出p(k>的值
p0=p1。
%给下次用
ife2<=Ebreak。
%如果参数收敛情况满足要求,终止计算
end%小循环结束
end%大循环结束
c,e%显示被辨识参数及其误差(收敛>情况
%分离参数
a1=c(1,:
>。
a2=c(2,:
>。
b1=c(3,:
>。
b2=c(4,:
>。
ea1=e(1,:
>。
ea2=e(2,:
>。
eb1=e(3,:
>。
eb2=e(4,:
>。
figure(2>。
%第二个图形
i=1:
15。
%横坐标从1到15
plot(i,a1,'r',i,a2,':
',i,b1,'g',i,b2,':
'>%画出a1,a2,b1,b2的各次辨识结果
title('ParameterIdentificationwithRecursiveLeastSquaresMethod'>
%图形标题
figure(3>。
%第三个图形
i=1:
15。
%横坐标从1到15
plot(i,ea1,'r',i,ea2,'g',i,eb1,'b',i,eb2,'r:
'>%画出a1,a2,b1,b2的各次辨识结果的收敛情况
title('IdentificationPrecision'>%图形标题
程序运行结果:
c=
0.001000.0010-0.4984-1.2328-1.4951-1.4962-1.4991-1.4998-1.4999
0.001000.00100.0010-0.23500.69130.69410.69900.69980.6999
0.001000.25091.24971.06651.00171.00201.00020.99990.9998
0.00100-0.24890.75000.56680.50200.50160.50080.50020.5002
-1.5000-1.5000-1.5000-1.4999-1.4999
0.70000.70000.70000.7000-0.7000
0.99990.99990.99990.99990.9999
0.50000.50000.50000.50000.5000
e=
000-499.42001.47340.21280.00070.00200.00040.0000
0000-235.9916-3.94160.00420.00700.00120.0001
00249.86123.9816-0.1466-0.06070.0003-0.0018-0.0003-0.0001
00-249.8612-4.0136-0.2443-0.1143-0.0007-0.0016-0.0012-0.0001
0.00010.0000-0.0000-0.00000.0000
0.0001-0.00000.00000.00000.0000
0.00010.00000.00000.0000-0.0000
-0.00040.0000-0.00000.0000-0.0000
图8.7最小二乘递推算法的参数辨识数
图8.8最小二乘递推算法的参数辨识结果
图8.9最小二乘递推算法的参数辨识结果收敛情况
表8.2最小二乘递推算法的辨识结果
参数a1a2b1b2
真值-1.50.71.00.5
估计值-1.49990.70.99990.5000
仿真结果表明,大约递推到第十步时,参数辨识的结果基本达到稳定状态,即a1=-1.4999,a2=0.7000,b1=0.9999,b2=0.5000。
此时,参数的相对变化量E0.000000005。
从整个辨识过程来看,精度的要求直接影响辨识的速度。
虽然最终的精度可以达到很小,但开始阶段的相对误差会很大,从图8.7(3>图形中可看出,参数的最大相对误差会达到三为数
8.5广义最小二乘算法的产生
广义最小二乘法所用的模型是:
<8.9)
<8.10)
其中我们选定模型参数:
a1=1.5,a2=-0.7,b1=1,b2=0.5,c1=0,c2=0
广义最小二乘递推算法的计算步骤:
1.给定初始条件
2.利用式
,计算
和
;
3.利用式
,构造
;
4.利用式
递推计算
;
5.利用
和
,计算
;
6.根据
来构造
;
7.利用
广义最小二乘法程序代码<见配套MATLAB程序algorithms5.m>
clear%清理工作间变量
L=55。
%M序列的周期
y1=1。
y2=1。
y3=1。
y4=0。
%四个移位寄存器的输出初始值
fori=1:
L。
%开始循环,长度为L
x1=xor(y3,y4>。
%第一个移位寄存器的输入是第三个与第四个移位寄存器的输出的“或”
x2=y1。
%第二个移位寄存器的输入是第一个移位寄存器的输出
x3=y2。
%第三个移位寄存器的输入是第二个移位寄存器的输出
x4=y3。
%第四个移位寄存器的输入是第三个移位寄存器的输出
y(i>=y4。
%取出第四个移位寄存器的幅值为"0"和"1"的输出信号,即M序列
ify(i>>0.5,u(i>=-1。
%如果M序列的值为"1",辨识的输入信号取“-1”
elseu(i>=1。
%如果M序列的值为"0",辨识的输入信号取“1”
end%小循环结束
y1=x1。
y2=x2。
y3=x3。
y4=x4。
%为下一次的输入信号做准备
end%大循环结束,产生输入信号u
z(2>=0。
z(1>=0。
%设z的前四个初始值为零
fork=3:
45。
%循环变量从5到45
z(k>=1.5*z(k-1>-0.7*z(k-2>+1*u(k-1>+0.5*u(k-2>。
%输出采样信号
end
%RGLS广义最小二乘辨识
c0=[0.0010.0010.0010.0