#v系统识别matlab第8章Word格式.docx
《#v系统识别matlab第8章Word格式.docx》由会员分享，可在线阅读，更多相关《#v系统识别matlab第8章Word格式.docx（37页珍藏版）》请在冰豆网上搜索。
A=6。
x0=1。
M=255。
f=2。
N=100；
%初始化；
x0=1。
fork=1:
N%乘同余法递推100次；
x2=A*x0。
%分别用x2和x0表示xi+1和xi-1；
x1=mod（x2,M>
%取x2存储器的数除以M的余数放x1<
xi）中；
v1=x1/256。
%将x1存储器中的数除以256得到小于1的随机数放v1中；
v（:
k>
=（v1-0.5>
*f。
%将v1中的数<
）减去0.5再乘以存储器f中的系数，存放在矩阵存储器v的第k列中，v（:
表示行不变、列随递推循环次数变化；
x0=x1。
%xi-1=xi；
v0=v1。
end%递推100次结束；
v2=v%该语句后无‘；
’，实现矩阵存储器v中随机数放在v2中，且可直接显示在MATLAB的window中；
k1=k。
%grapher%以下是绘图程序；
k=1:
k1。
plot（k,v,k,v,'
r'
>
xlabel（'
k'
ylabel（'
v'
tktle（'
（-1,+1>
均匀分布的白噪声'
②程序运行结果如图8.1所示。
图8.1采用MATLAB产生的（-1,+1>
均匀分布的白噪声序列
③产生的（-1，1>
在程序运行结束后，产生的（-1，1>
均匀分布的白噪声序列，直接从MATLAB的window界面中copy出来如下<
v2中每行存6个随机数）：
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.2500-0.48440.1016-0.36720.8047-0.1328
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
=Y4。
end
end
M=U
%绘图
i1=i
1:
i1。
plot（k,U,k,U,'
rx'
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
Columns41through50
Columns51through60
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>
%用理想输出值作为观测值
subplot（3,1,1>
%画三行一列图形窗口中的第一个图形
stem（u>
%画输入信号u的径线图形
subplot（3,1,2>
%画三行一列图形窗口中的第二个图形
i=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>
-z（4>
-z（3>
u（4>
-z（5>
-z（4>
u（5>
-z（6>
-z（5>
u（6>
-z（7>
-z（6>
u（7>
-z（8>
-z（7>
u（8>
-z（9>
-z（8>
u（9>
-z（10>
-z（9>
u（10>
-z（11>
-z（10>
u（11>
-z（12>
-z（11>
u（12>
-z（13>
-z（12>
u（13>
-z（14>
-z（13>
u（14>
-z（15>
-z（14>
u（15>
]%给样本矩阵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。
%四个移位寄存器的输出初始值
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序列的值为"
辨识的输入信号取“-0.03”
elseu（i>
=0.03。
辨识的输入信号取“0.03”
end%小循环结束
y1=x1。
y2=x2。
y3=x3。
y4=x4。
%为下一次的输入信号做准备
end%大循环结束，产生输入信号u
figure（1>
%第一个图形
gridon%显示出输入信号径线图并给图形加上网格
z（2>
=0。
z（1>
%设z的前两个初始值为零
15。
%循环变量从3到15
z（k>
%输出采样信号
%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>
%相对误差的初始值及大小
%开始求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（:
=e2。
%把当前相对变化的列向量加入误差矩阵的最后一列
c0=c1。
%新获得的参数作为下一次递推的旧参数
c（:
=c1。
%把辨识参数c列向量加入辨识参数矩阵的最后一列
p1=p0-k1*k1'
*[h1'
*p0*h1+1]。
%求出p（k>
的值
p0=p1。
%给下次用
ife2<
=Ebreak。
%如果参数收敛情况满足要求，终止计算
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>
%第二个图形
%横坐标从1到15
plot（i,a1,'
i,a2,'
:
'
i,b1,'
g'
i,b2,'
%画出a1，a2，b1，b2的各次辨识结果
ParameterIdentificationwithRecursiveLeastSquaresMethod'
%图形标题
figure（3>
%第三个图形
plot（i,ea1,'
i,ea2,'
i,eb1,'
b'
i,eb2,'
r:
%画出a1，a2，b1，b2的各次辨识结果的收敛情况
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>
L=55。
辨识的输入信号取“-1”
辨识的输入信号取“1”
%设z的前四个初始值为零
45。
%循环变量从5到45
+1*u（k-1>
%RGLS广义最小二乘辨识
c0=[0.0010.0010.0010.0
升级会员 