Matlab系统辨识应用例子Word格式.docx
《Matlab系统辨识应用例子Word格式.docx》由会员分享，可在线阅读，更多相关《Matlab系统辨识应用例子Word格式.docx（20页珍藏版）》请在冰豆网上搜索。
-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）]%给样本矩阵赋值
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）从觥中分离出并显示◎、a2、b、
%End
程序运行结果：
>
u=[-1,1,-1,1,1,1，1,-1,-1,-1,1，-1,-1,1，1]
z=[00,,,,,,,,]
HL=
『00\
ZL=[,,,,,,,,,,,,,]c=[,,,]T
a1=
a2=b1=b2=
-1
10
图2最小二乘一次完成算法仿真实例中输入信号和输岀观测值
STEMPLO
TOFINF
UTP
MSERIE、
S
5
15
rm
从仿真结果表1可以看出，由于所用的输出观测值没有任何噪声成分，所以辨识结果也无任何误差。
例2根据热力学原理，对给定质量的气体，体积V与压力P之间的关系为PV，其中和为待定参数。
经实验获得如下一批数据，V的单位为立方英寸，P的单位为帕每平方英寸。
V
P
试用最小二乘一次完成算法确定参数和。
首先要写出系统的最小二乘表达式。
为此，把体积V与压力P之间的关系PV改为对数关系，即，logPlogVlog。
此式与式z（k）h（k）0e（k）,对比可得：
z（k）logP,h（k）[logV1],e[log]。
例2的程序如下。
%实际压力系统的最小二乘辨识程序，文件名：
clear%工作间清零
V=[,,,,,]'
，P=[,,,,,]'
%赋初值并显示V、P%logP=-alpha*logV+logbeita=[-logV,1][alpha,log（beita）]'
=HL*sita%注释P、V之间的关系
fori=1:
6;
Z（i）=log（P（i））;
%循环变量的取值为从1到6，系统的采样输出赋值
End%循环结束
ZL=Z'
%zL赋值HL=[-log（V
（1））,1;
-log（V
（2））,1;
-log（V（3））,1;
-log（V（4））,1;
-log（V（5））,1;
-log（V（6））,1]
%Hl赋值
%CalculatingParameters
c仁HL'
c2=inv（c1）;
c3二HL'
c4二c2*c3计算被辨识参数的值
%SeparationofParameters
alpha二c4
（1）%为c4的第一个元素
beita二exp（c4
（2））%为以自然数为底的c4的第二个元素的指数
V=[,,,,,]
P=[,,,,,]
ZL=[,,,,,]
HL=c4=
r，
alpha=
beita=+004
仿真结果表明，用最小二乘一次完成算法可以迅速辨识出系统参数,
即=，=+004。
例3考虑图3所示的仿真对象，图中，v（k）是服从N（O,1）分布
的不相关随机噪声。
且
11
G（z1）B（z），n（z1）D（z），
A（zJ）C（z1）
A（z1）11.5z10.7z2C（z1）
112
B（z）1.0z0.5z
1
D（z）1
经过计算，得到系统真实的模型：
选择图3所示的辨识模型。
仿真对象选择如下的模型结构
z（k）a1Z（k1）a2Z（k2）be（k1）bu（k2）v（k）
输入信号采用4位移位
u（k）
M/_1\
厶丿
»
G（z1）
-X
e（k）
z（k）
y（k）
图3最小二乘递推算法辨识实例结构图
寄存器产生的M序列，幅度为。
最小二乘递推算法辨识的程序流程如图4所示。
下面给出具体程序。
%最小二乘递推算法辨识程序，在光盘中的文件名:
clear%清理工作间变量
L=15;
%M序列的周期
y1=1;
y2=1;
y3=1;
y4=0;
%四个移位寄存器的输出初始值
L;
%开始循环，长度为L
x仁xor（y3,y4）;
%第一个移位寄存器的输入是第三个与第四个移位寄存器的输出的“或”
x2=y1;
%第二个移位寄存器的输入是第一个移位寄存器的输出
x3=y2;
%第三个移位寄存器的输入是第二个移位寄存器的输出
x4=y3;
%第四个移位寄存器的输入是第三个移位寄存器的输出
y（i）=y4;
%取出第四个移位寄存器的幅值为"
0"
和"
1"
的输出信号即M序列
ify（i）>
u（i）=;
%如果M序列的值为"
辨识的输入信号取“-0.03”
elseu（i）=;
辨识的输入信号取“0.03”
end%小、循环结束
y仁x1;
y2=x2;
y3=x3;
y4=x4;
%为下一次的输入信号做准备
end%大循环结束，产生输入信号u
figure
（1）;
%第一个图形
stem（u）,gridon%显示出输入信号径线图并给图形加上网格
z
（2）=0;
z
（1）=0;
%设z的前两个初始值为零
15;
%循环变量从3到15
z（k）二*z（k-1）*z（k-2）+u（k-1）+*u（k-2）;
瀚出采样信号
%RLS递推最小二乘辨识
cO=[]'
;
%直接给出被辨识参数的初始值，即一个充分小的实向量
p0=10"
6*eye（4,4）;
%直接给出初始状态P0,即一个充分大的实数单位矩阵
E=;
%取相对误差E=
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）
k仁pO*h1*x1;
%求出K的值
d1=z（k）-h1'
*c0;
c1=c0+k1*d1;
%求被辨识参数c
e1=c1-c0;
%求参数当前值与上一次的值的差值
e2=e1./c0;
%求参数的相对变化
e（:
k）=e2;
%把当前相对变化的列向量加入误差矩阵的第k列
c0=c1;
%新获得的参数作为下一次递推的旧参数
c（:
k）=c1;
%把辨识参数c列向量加入辨识参数矩阵的第k列
p仁p0-k1*h1'
*p0;
%求出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）;
%第二个图形
%横坐标从1到15
plot（i,a1,'
r'
i,a2,'
:
'
i,b1,'
g'
i,b2,'
）%画出a1，a2，b1，b2的各次辨识
结果
title（'
ParameterIdentificationwithRecursiveLeastSquaresMethod'
）%
图形标题
figure（3）;
%第三个图形
%横坐标从1至U15
plot（i,ea1,'
i,ea2,'
i,eb1,'
b'
i,eb2,'
r:
）%画出al,a2,bl,b2的各次辨
识结果的收敛情况
IdentificationPrecision'
）%图形标题
程序运行结果:
c=
0.03
0.02
0.01
-0.01
-0.02
-0.03
300
200
100
-100
-200
-300
-400
-500
r
■
-MConvergenceofb1
-
■Convergenceofa2
Convergence
-ofb2
Convergenceofa1
IdentificationPrecision
Fig.3IdentificationError
图5最小二乘递推算法的参数辨识仿真
e=
000
00
参数aia2bib2
真值
估计值
仿真结果表明，大约递推到第十步时，参数辨识的结果基本达到稳定状态，即ai=,a2=,bi=,b2=。
此时，参数的相对变化量。
从整个辨识过程来看，精度的要求直接影响辨识的速度。
虽然最终的精度可以达到很小，但开始阶段的相对误差会很大。
例四：
考虑图4所示的仿真对象，图中，v（k）是服从N（0,1）分布的不
相关随机噪声。
G（z1）B（z），N（z1）D（zI
A（z1）C（z1）
A（z1）11.5a1z10.7z2C（z1）
B（z1）1.0z10.5z2
D（z'
）1z10.2z2
工作间清零二＞
给M序列的长度L赋值
模型结构选用如下形式：
z（k）a1z（k1）a2（k2）b1u（k1）b2u（k2）v（k）d1v（k1）d2（v（k2）
增广最小二乘辨识程序流程上图所示。
程序如下。
%增广最小二乘辨识程序
clear
L=60;
%4位移位寄存器产生的M序列的周期y1=1;
%4个移位寄存器的输出初始值fori=1:
x1=xor（y3,y4）;
%第一个移位寄存器的输入信号
%第二个移位寄存器的输入信号
%第三个移位寄存器的输入信号
%第四个移位寄存器的输入信号
%第四个移位寄存器的输出信号，M序列,幅值"
，ify（i）>
u（i）=-1;
%M序列的值为"
时,辨识的输入信号取“-1”elseu（i）=1;
%M序列的值为"
时，辨识的输入信号取“1”end
y1=x1;
%为下一次的输入信号准备
%画第一个图形
subplot（2,1,1）;
%画第一个图形的第一个子图
stem（u）,gridon%画出M序列输入信号
v=randn（1,60）;
%产生一组60个正态分布的随机噪声subplot（2,1,2）;
%画第一个图形的第二个子图plot（v）,gridon;
%画出随机噪声信号
u，v%显示输入信号和噪声信号
z=zeros（7,60）;
zs二zeros（7,60）;
zm二zeros（7,60）;
zmd二zeros（7,60%输出采样矩阵、不考虑噪声时系统输出矩阵、不考虑噪声时模型输出矩阵、模型输出矩阵的大小
z
（1）=0;
zs
（2）=0;
zs
（1）=0;
zm
（2）=0;
zm
（1）=0;
zmd
（2）=0;
zmd
（1）=0;
%输出采样、不考虑噪声时系统输出、不考虑噪声时模型输出、模型输出的初值
c0=[]'
%直接给出被辨识参数的初始值,即一个充分小的实向量
p0=10A6*eye（7,7）;
%直接给出初始状态P0,即一个充分大的实数单位
矩阵
%相对误差E=
c=[c0,zeros（7,14）];
%被辨识参数矩阵的初始值及大小
e=zeros（7,15）;
%相对误差的初始值及大小
60;
%开始求Kz（k）=*z（k-1）*z（k-2）+u（k-1）+*u（k-2）+v（k）-v（k-1）+*v（k-2）;
%系统在M序列输入下的输出采样信号
h仁[-z（k-1）,-z（k-2）,u（k-1）,u（k-2）,v（k）,v（k-1）,v（k-2）]'
%为求K（k）作准
备
k1=p0*h1*x1;
%K
%辨识参数czs（k）=*z（k-1）+*z（k-2）+u（k-1）+*u（k-2）;
%系统在M序列的输入下的输出响应
zm（k）=[-z（k-1）,-z（k-2）,u（k-1）,u（k-2）]*[c1
（1）;
c1
（2）;
c1（3）;
c1（4）];
%模型
在M序列的输入下的输出响应
zmd（k）=h1'
*c1;
%模型在M序列的输入下的输出响应e1=c1-c0;
e2=e1./c0;
%求参数误差的相对变化e（:
%给下一次用
%把递推出的辨识参数c的列向量加入辨识参数矩阵p1=p0-k1*k1'
*[h1'
*p0*h1+1];
%findp（k）
%若收敛情况满足要求，终止计算
end%判断结束
end%循环结束
c,e,%显示被辨识参数及参数收敛情况
z,zs,zm%显示输出采样值、系统实际输出值、模型输出值
%分离变量
%分离出a1、a2、b1、b2d1=c（5,:
d2=c（6,:
d3=c（7,:
%分离出d1、d2、d3
ea1=e（1,:
%分离出a1、a2、b1、b2的收敛情况
ed1=e（5,:
ed2=e（6,:
ed3=e（7,:
%分离出d1、d2、d3的收敛情况figure
（2）;
%画第二个图形i=1:
b:
i,d1,'
i,d2,'
g:
i,d3,'
g+'
）%画出各个被辨识参数
标题
i=1:
%画出第三个图形
Plot（i,ea1,'
i,ed1,'
i,ed2,'
r+'
）%画
出各个参数收敛情况
）%标题
figure（4）;
subplot（4,1,1）;
%画出第四个图形，第一个子图
plot（i,zs（i）,'
）%画出被辨识系统在没有噪声情况下的实际输出响应
subplot（4,1,2）;
plot（i,z（i）,'
）%第二个子图,画出被辨识系统的采
样输出响应
subplot（4,1,3）;
plot（i,zm（i）,'
）%第三个子图,画出模型含有噪声
的输出响应
subplot（4,1,4）;
）%第四个子图，画出模型去除噪声
后的输出响应
程序执行结果：
u=
[1,-1,-1,-1,-1,1,1,1,-1,1,1,-1,-1,1,-1,1,-1,-1,1,-1,-1,1,1,-1,1,1,-1,-1,1,-1,1,-1,-1,-1,-1,1,1,1,-1,1,1,-1,-1,1,-1,1,-1,
-1,-1,-1,1,1,1,-1,1,1,-1,-1,1,-1]
v=
[,,,,,,,,,,,,,,
JJJJJJJJJJJJJJ
,,,,,,,,,,J,,,,,,,,,,,,,,,,,,]
SO0
0000
t00
z=
[0,0,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,，,,,,,,,,,]
zs=
[0,0,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,]
zm=
[0,0,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,]
-10
20
-20
；
System
ResponseW
ithoutNoises
l
30
40
50
60
ithNoises
■ModelResponseWithn
oises
i
ModelResponseWithoutN
oise
一
_—i—-
Fig.4Responses
图5增广最小二乘递推算法辨识仿真
表3增广最小二乘递推算法的辨识结果
参数a1a2b1b2d1d2d3
仿真结果表明，递推到第9步时，参数辨识的结果基本达到稳定状态,即a1=,a2=,b1=,b2=,d1=,d2=,d3=。
此时，辨识参数的相对变化量<
与最小二乘递推算法相比，增广最小二乘递推算法虽
clearall;
N=1000;
A=[1]；
B=[0];
C=[1]；
M仁idpoly（A,B,C）;
step（M1,[0100]）;
grid;
%
U=iddata（[],idinput（N,'
prbs'
））;
%E=iddata（[],idinput（N,'
rgs'
丫仁sim（M1,[U,E]）;
%outputdata
Z=iddata（Y1,U）;
%ARMAXorderguji
NN=struc（1:
2,1:
4,1:
4）;
Loss_fun=arxstruc（Z,Z,NN）;
order=selstruc（Loss_fun,'
aic'
order=[order
（1）,order
（2）,1,order（3）]%ARMAXprametersgujiModel_para=armax（Z,order）;
present（Model_para）;
compare（Z,Model_para）;
U=[0,,,-,-,,,]'
Z=[0,,,-,-,,,]'
%systeminputoutputinitalvalu
升级会员 