Matlab系统辨识应用例子.docx
- 文档编号:9382199
- 上传时间:2023-02-04
- 格式:DOCX
- 页数:20
- 大小:172KB
Matlab系统辨识应用例子.docx
《Matlab系统辨识应用例子.docx》由会员分享,可在线阅读,更多相关《Matlab系统辨识应用例子.docx(20页珍藏版)》请在冰豆网上搜索。
Matlab系统辨识应用例子
例1、考虑仿真对象
z(k)1.5z(k1)0.7z(k2)u(k1)0.5u(k2)v(k)
其中,v(k)是服从正态分布的白噪声N(0,1)。
输入信号采用4阶M序
列,幅度为1。
选择如下形式的辨识模型
z(k)a1z(k1)a2z(k2)b1u(k1)pu(k2)v(k)
设输入信号的取值是从k=1到k=16的M序列,则待辨识参数?
ls
为?
LS=(T)1Tz。
其中,被辨识参数?
LS、观测矩阵的表达式为:
a1
z(3)
?
a2
址s.
z(4)
z
L
b1
b2
z(16)
z
(2)z
(1)u
(2)u
(1)
z(3)z
(2)u(3)u
(2)
LL
z(15)z(14)u(15)u(14)
程序框图如图1所示。
Matlab仿真程序如下:
图1最小二乘一次完成算法程序框图
%二阶系统的最小二乘一次完成算法辨识程序,文件名:
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)=*z(k-1)*z(k-2)+u(k-1)+*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)]%给样本矩阵赋值
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计算并显示?
ls
%DisplayParameters
a1=c
(1),a2=c
(2),b1=c(3),b2=c(4)从觥中分离出并显示◎、a2、b、
b2
%End
程序运行结果:
>>
u=[-1,1,-1,1,1,1,1,-1,-1,-1,1,-1,-1,1,1]
z=[00,,,,,,,,]
HL=
『00\
0
ZL=[,,,,,,,,,,,,,]c=[,,,]T
a1=
a2=b1=b2=
-1
10
图2最小二乘一次完成算法仿真实例中输入信号和输岀观测值
STEMPLO
TOFINF
UTP
MSERIE、
S
0
5
10
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'*HL;c2=inv(c1);c3二HL'*ZL;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
经过计算,得到系统真实的模型:
z(k)1.5z(k1)0.7z(k2)u(k1)0.5u(k2)v(k)
选择图3所示的辨识模型。
仿真对象选择如下的模型结构
z(k)a1Z(k1)a2Z(k2)be(k1)bu(k2)v(k)
其中,v(k)是服从正态分布的白噪声N(0,1)。
输入信号采用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;%四个移位寄存器的输出初始值
fori=1:
L;%开始循环,长度为L
x仁xor(y3,y4);%第一个移位寄存器的输入是第三个与第四个移位寄存器的输出的“或”
x2=y1;%第二个移位寄存器的输入是第一个移位寄存器的输出
x3=y2;%第三个移位寄存器的输入是第二个移位寄存器的输出
x4=y3;%第四个移位寄存器的输入是第三个移位寄存器的输出
y(i)=y4;%取出第四个移位寄存器的幅值为"0"和"1"的输出信号即M序列
ify(i)>,u(i)=;%如果M序列的值为"1",辨识的输入信号取“-0.03”
elseu(i)=;%如果M序列的值为"0",辨识的输入信号取“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的前两个初始值为零
fork=3:
15;%循环变量从3到15
z(k)二*z(k-1)*z(k-2)+u(k-1)+*u(k-2);瀚出采样信号
end
%RLS递推最小二乘辨识
cO=[]';%直接给出被辨识参数的初始值,即一个充分小的实向量
p0=10"6*eye(4,4);%直接给出初始状态P0,即一个充分大的实数单位矩阵
E=;%取相对误差E=
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)
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);%第二个图形
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至U15
plot(i,ea1,'r',i,ea2,'g',i,eb1,'b',i,eb2,'r:
')%画出al,a2,bl,b2的各次辨
识结果的收敛情况
title('IdentificationPrecision')%图形标题
程序运行结果:
>>
c=
0
0
0.03
0.02
0.01
-0.01
-0.02
-0.03
300
200
100
0
-100
-200
-300
-400
-500
1
r
■
-MConvergenceofb1
-
-
-
■Convergenceofa2
-
Convergence
-ofb2
Convergenceofa1
-
■
IdentificationPrecision
5
10
15
Fig.3IdentificationError
图5最小二乘递推算法的参数辨识仿真
e=
000
000
00
00
参数aia2bib2
真值
估计值
仿真结果表明,大约递推到第十步时,参数辨识的结果基本达到稳定状态,即ai=,a2=,bi=,b2=。
此时,参数的相对变化量。
从整个辨识过程来看,精度的要求直接影响辨识的速度。
虽然最终的精度可以达到很小,但开始阶段的相对误差会很大。
例四:
考虑图4所示的仿真对象,图中,v(k)是服从N(0,1)分布的不
相关随机噪声。
11
G(z1)B(z),N(z1)D(zI
A(z1)C(z1)
A(z1)11.5a1z10.7z2C(z1)
B(z1)1.0z10.5z2
112
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;y2=1;y3=1;y4=0;%4个移位寄存器的输出初始值fori=1:
L;
x1=xor(y3,y4);%第一个移位寄存器的输入信号
x2=y1;%第二个移位寄存器的输入信号
x3=y2;%第三个移位寄存器的输入信号
x4=y3;%第四个移位寄存器的输入信号
y(i)=y4;%第四个移位寄存器的输出信号,M序列,幅值"0"和"1",ify(i)>,u(i)=-1;%M序列的值为"1"时,辨识的输入信号取“-1”elseu(i)=1;%M序列的值为"0"时,辨识的输入信号取“1”end
y1=x1;y2=x2;y3=x3;y4=x4;%为下一次的输入信号准备
end
figure
(1);%画第一个图形
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
(2)=0;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=;%相对误差E=
c=[c0,zeros(7,14)];%被辨识参数矩阵的初始值及大小
e=zeros(7,15);%相对误差的初始值及大小
fork=3:
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)作准
备
x=h1'*p0*h1+1;x1=inv(x);k1=p0*h1*x1;%K
d1=z(k)-h1'*c0;c1=c0+k1*d1;%辨识参数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(:
k)=e2;
c0=c1;%给下一次用
c(:
k)=c1;%把递推出的辨识参数c的列向量加入辨识参数矩阵p1=p0-k1*k1'*[h1'*p0*h1+1];%findp(k)
p0=p1;%给下次用
ife2<=Ebreak;%若收敛情况满足要求,终止计算
end%判断结束
end%循环结束
c,e,%显示被辨识参数及参数收敛情况
z,zs,zm%显示输出采样值、系统实际输出值、模型输出值
%分离变量
a1=c(1,:
);a2=c(2,:
);b1=c(3,:
);b2=c(4,:
);%分离出a1、a2、b1、b2d1=c(5,:
);d2=c(6,:
);d3=c(7,:
);%分离出d1、d2、d3
ea1=e(1,:
);ea2=e(2,:
);eb1=e(3,:
);eb2=e(4,:
);%分离出a1、a2、b1、b2的收敛情况
ed1=e(5,:
);ed2=e(6,:
);ed3=e(7,:
);%分离出d1、d2、d3的收敛情况figure
(2);%画第二个图形i=1:
60;plot(i,a1,'r',i,a2,'r:
',i,b1,'b',i,b2,'b:
',i,d1,'g',i,d2,'g:
',i,d3,'g+')%画出各个被辨识参数
title('ParameterIdentificationwithRecursiveLeastSquaresMethod')%
标题
figure(3);i=1:
60;%画出第三个图形
Plot(i,ea1,'r',i,ea2,'r:
',i,eb1,'b',i,eb2,'b:
',i,ed1,'g',i,ed2,'g:
',i,ed2,'r+')%画
出各个参数收敛情况
title('IdentificationPrecision')%标题
figure(4);subplot(4,1,1);%画出第四个图形,第一个子图
i=1:
60;plot(i,zs(i),'r')%画出被辨识系统在没有噪声情况下的实际输出响应
subplot(4,1,2);i=1:
60;plot(i,z(i),'g')%第二个子图,画出被辨识系统的采
样输出响应
subplot(4,1,3);i=1:
60;plot(i,zm(i),'b')%第三个子图,画出模型含有噪声
的输出响应
subplot(4,1,4);i=1:
60;plot(i,zs(i),'b')%第四个子图,画出模型去除噪声
后的输出响应
程序执行结果:
>>
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,,,,,,,,,,,,,,,,,,]
c=
0
0
0
0
e=
SO0
0000
00
00
00
00
t00
z=
[0,0,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,]
zs=
[0,0,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,]
zm=
[0,0,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,]
10
0
-10
10
0
-10
10
0
-10
20
0
-20
;System
ResponseW
ithoutNoises
l
L
0
10
20
30
40
50
60
'System
ResponseW
ithNoises
0
10
20
30
40
50
60
■ModelResponseWithn
'oises
i
r
r
0
10
20
30
40
50
60
'ModelResponseWithoutN
oise
i
L
一
_—i—-
0
10
20
30
40
50
60
Fig.4Responses
Fig.3IdentificationError
图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);
clear
U=[0,,,-,-,,,]';
Z=[0,,,-,-,,,]';
%systeminputoutputinitalvalu
- 配套讲稿:
如PPT文件的首页显示word图标,表示该PPT已包含配套word讲稿。双击word图标可打开word文档。
- 特殊限制:
部分文档作品中含有的国旗、国徽等图片,仅作为作品整体效果示例展示,禁止商用。设计者仅对作品中独创性部分享有著作权。
- 关 键 词:
- Matlab 系统 辨识 应用 例子