#v系统识别matlab第8章.docx

#v系统识别matlab第8章.docx

《#v系统识别matlab第8章.docx》由会员分享，可在线阅读，更多相关《#v系统识别matlab第8章.docx（37页珍藏版）》请在冰豆网上搜索。

#v系统识别matlab第8章

第8章系统辨识MATLAB仿真实例

MATLAB是一套高性能数字计算和可视化软件,它集成概念设计,算法开发,建模仿真,实时实现于一体,构成了一个使用方便、界面友好的用户环境,其强大的扩展功能为各领域的应用提供了基础。

对于一个简单的系统,可以通过分析其过程的运动规律,应用一些已知的定理和原理,建立数学模型,即所谓的“白箱建模”。

但对于比较复杂的生产过程,该建模方法有很大的局限性。

因为过程的输入输出信号一般总是可以测量的,而且过程的动态特性必然表现在这些输入输出数据中,那么就可以利用输入输出数据所提供的信息来建立过程的数学模型。

这种建模方法就称为系统辨识。

把辨识建模称作“黑箱建模”。

本章列举了在系统辨识中常用的基本算法的MATLAB程序。

8.1白噪声的产生

如果在计算机上比较经济的产生统计上理想的各种不同分布的白噪声序列，则对系统辨识仿真研究提供了极大方便。

为了简单起见，常把各种不同分布的白噪声序列称为随机数。

从理论上讲，只要有了一种具有连续分布的随机数，就可以通过函数变换的方法产生其他任意分布的随机数。

显然，在具有连续分布的随机数中，<0,1）均匀分布的随机数是最简单、最基本的一种，有了<0,1）均匀分布的随机数，便可以产生其他任意分布的随机数和白噪声。

本文首先产生<0,1）均匀分布的随机序列，然后对程序稍加改动，将产生<0,1）均匀分布的随机数统统减去0.5，相当于将原随机序列图的横坐标向上平移0.5，原随机序列变成了<-0.5,0.5）的白噪声，然后乘以存储器f的系数，这里取法f=2，则可得到产生<-1,1）均匀分布的白噪声。

8.1.1乘同余法的介绍

用如下递推同余式产生正整数序列{

}

<8.1）

式中，mod为取M的余数，如：

M为2的方幂，即

，k为k>2是整数；若k=3，则

初值

取正奇数，如

=1。

再令

（8.2>

可以证明序列{

}是伪随机数序列；同时还可以证明伪随机序列{

}的循环周期达到最大值

。

将式<8.1）和<8.2）合并为

<8.3）

式中，初值为

。

8.1.2举例

利用乘同余法，选R=2,A=6,k=8,

递推100次，采用MATLAB的仿真语言编程，产生<-1,1）均匀分布的白噪声。

用乘同余法产生（见配套MATLAB程序algorithms1.m>

1编程如下：

A=6。

x0=1。

M=255。

f=2。

N=100；%初始化；

x0=1。

M=255。

fork=1:

N%乘同余法递推100次；

x2=A*x0。

%分别用x2和x0表示xi+1和xi-1；

x1=mod（x2,M>。

%取x2存储器的数除以M的余数放x1

v1=x1/256。

%将x1存储器中的数除以256得到小于1的随机数放v1中；

v（:

k>=（v1-0.5>*f。

%将v1中的数<

）减去0.5再乘以存储器f中的系数，存放在矩阵存储器v的第k列中，v（:

k>表示行不变、列随递推循环次数变化；

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=

-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