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