系统辨识经典辨识方法.docx
- 文档编号:29456182
- 上传时间:2023-07-23
- 格式:DOCX
- 页数:23
- 大小:641.30KB
系统辨识经典辨识方法.docx
《系统辨识经典辨识方法.docx》由会员分享,可在线阅读,更多相关《系统辨识经典辨识方法.docx(23页珍藏版)》请在冰豆网上搜索。
系统辨识经典辨识方法
经典辨识方法报告
1.面积法
1.1辨识原理
1.1.1分子多项式为1的系统
由于系统的传递函数与微分方程存在着一一对应的关系,因此,可以通过求取微分方程的系数来辨识系统的传递函数。
在求得系统的放大倍数K后,要先得到无因次阶跃响应y(t)(设T=0)o大多数自衡的工业过程对象的y(t)可以用下式描述來近似
面积法原则上可以求出n为任意阶的各系数。
以n=3为例,注意到
将式(2.1.2)的y(t)项移至右边,在[0,t]上积分,得
定义
F](t)=[[l-y(t)]dt(1.5)
则由式(2.1.3)给出的条件可知,在t-H»
(1.6)
将式aiy(t)移到等式右边,定义
(1.7)
a3畔+4y(t)=[[F】(t)_a】y(t)]dt=F2(t)
利用初始条件(2.13)当t->oo时
同理有a3=F3(oo)
以此类推,?
fn>2>有an=Fn(co)
1.1.2分子、分母分别为m阶和n阶多项式的系统
当传递函数的形式如下所示时
G(s)=K
bmSm+bwS10"+…+b(s+l
ansn++•••+ajS+1
(n>m)
K」(%
(19)
定义
G(s)=K
(1.10)
(111)
L[l-hXt)]=f[l-h*(t)]e^dt
则[1—h(t)]的Laplace变换为:
00
11工才
L[l-lf(t)]=—=
ssP(s)1+工阳
(112)
1=1
定义一阶而积A为:
A二『[1-h*(t)]dt=l^L[l-h*(t)]=lini
'1+工W
oo
EQs-1
二C]
巨'(1.13)
令
wC(t)i=
1
S(1+C]S)
(114)
定义一阶而积为,
A二『〔[】")-hSdzdt二hm
5->0
8
Zes"
i=2
(1+qsXl+EqsJ
1-1
同理,令
-(1.17)
上式可写成如下形式:
定义i阶面积为A=q。
由此可得:
gL++…+a】s+1=(bmsm+bm_1sm_1+...+bxs+1)(1+£举”)
通过该系数矩阵,即可求出传递函数分子分母系数的值。
1.2程序设计
1.2.1传递函数形式如式1.1的系统
取系统传递函数如下:
MATLAB程序如下:
clc%清空工作区
clear
dt=0.01;%设置采样时间
t=0:
dt:
50;%设置时间长度
num=l;%此系统分子为1
den=[325.21];%分母多项式系数
%绘制原传递函数阶跃响应曲线fprintff原系统传递函数为
G=tf(num,den)y=step(num/den,t);
Length=length(y);%数据长度
Plot(t,y);
grid;
xlabel('t/s');
ylabel('y(t)');
%进行辨识设计
fprintff辨识参数结果
%求al
suml=O;
for(i=l:
Length)
:
');
suml=suml+(l-y(i))*dt;F(i)=suml;
end
al=suml
%求a2
sum2=0;
for(i=l:
Length)sum2=sum2+(F(i)-al*y(i))*dt;f(i)=sum2;
end
a2=sum2
%求a3
sum3=0;
for(i=l:
Length)sum3=sum3+(f(i)-a2*y(i))*dt;
end
a3=sum3
%绘制辨识后的传递函数
dt=0.01;
t=0:
dt:
50;
num2=l;
den2=[a3a2al1];
fprintfC系统辨识后的传递函数为门
G=tf(num2,den2)
h=step(num2,den2/t);%辨识所得•传递函数阶跃响应
plot(t,y,'black'd,h,'blue');
legendC原传递函数T辨识所得传递函数J;
titler原传递函数与辨识所得传递函数的阶跃响应对比Jgrid;
xlabel('t/s');
ylabefy⑴和h(t)');
fprintffffl关系数:
J;%求相关系数
r=corrcoef(yzh)
运行以上程序得到结果如下:
原系统传递函数为:
G=
1
3sA3+2sA2+5.2s+l
Continuous-timetransferfunction.
辨识参数结果:
al=
5.2048
a2=
2.0608
a3=
2.8388
系统辨识后的传递函数为:
G=
1
2.839sA3+2.061sA2+5.205s+1Continuous-timetransferfunction.
相关系数:
1.00000.9999
0.99991.0000
此时原传递函数和辨识所得传递函数的阶跃响应对比如下图:
1.2.1传递函数形式如式(1・9)的系统(无噪声)
取系统传递函数如下:
s?
+3s+1G")2s3+5.5s2+3s+1
MATLAB程序如下:
clc%清空工作区
clear
dt=0.01;%设置采样时间
t=0:
dt:
20;%设置时间长度
num=[l31];%分子多项式系数
den=[25.531];%分母多项式系数
%绘制原传递函数阶跃响应曲线fprintfC原系统传递函数为门
G=tf(num,den)y=step(num,den,t);
Length=length(y);
Plot(t,y);
grid;
xlabelCt/s*);
ylabelCy(tr);
%辨识程仔设计,此系统m+n=5,故应计算A1-A5
fprintfCAl-A5阶面积分别为」)
%求A1
suml=O;
for(i=l:
Length-l)
suml=suml+(l-(y(i)+y(i+l))/2)*dt;A(i)=suml;
end
Al=suml
%求A2
sum2=0;
for(i=l:
Length-l)sum2=sum2+(A(i)-Al*(y(i)+y(i+l))/2)*dt;B(i)=sum2;
end
A2=sum2
%求A3
sum3=0;
for(i=l:
Length-l)sum3=sum3+(B(i)-A2*(y(i)+y(i+l))/2)*dt;C(i)=sum3;
end
A3=sum3
%求A4
sum4=0;
for(i=l:
Length-l)sum4=sum4+(C(i)-A3*(y(i)+y(i+l))/2)*dt;D(i)=sum4;
end
A4=sum4
%求A5
sum5=0;
for(i=l:
Length-l)sum5=sum5+(D(i)-A4*(y(i)+y(i+l))/2)*dt;
end
A5=sum5
%求分子系数bl,b2
M=(-l)*(inv([A3/A2;A4,A3]))*[A4;A5];
fprintf('分子多项式系数为:
J
bl=M(lzl)
b2=M(2/l)
%求分母系数31,32,33
N=[l0O;A110;A2All]*[bl;b2;0]+[Al;A2;A3];
fprintf('分母多项式系数为:
’)
al=N(l,l)
a2=N(2,l)
a3=N(3,l)
%求辨识所得传递函数
numl=[b2bl1];
denl=[a3a2al1);
fprintfC辨识所得传递函数为:
')
G=tf(numl,denl)
h=step(numl,denl/t);
plot(t,y,'black',t,h,'blue');
legendf原传递函数丁辨识所得传递函数J;
title/原传递函数与辨识所得传递函数的阶跃响应对比J
grid;
xlabel('t/s');
ylabelCy(t)');
fprintfC相关系数:
');%求相关系数
r=corrcoef(yzh)
运行以上程序结果如下:
原系统传递函数为:
G=
sA2+3s+1
2sA3+5.5sA2+3s+1Continuous-timetransferfunction.
A1-A5阶面积分别为:
Al=
0.0076
A2=
4.3364
A3=
-9.6285
A4=
15.7297
A5=
7.1087
分子多项式系数为:
bl=
7.4409
b2=
12.8942
分母多项式系数为:
al=
7.4485
a2=
17.2869
a3=
22.7359
辨识所得传递函数为:
G=
12.89sA2+7.441s+1
22.74sA3+17.29sA2+7.448s+1Continuous-timetransferfunction.相关系数:
r=
1.0000
0.9996
0.99961.0000
此时原传递函数和辨识所得传递函数的阶跃响应对比如下图:
MATLAB程序如下:
clc%清空工作区
clear
dt=0.01;%设置采样时间
t=0:
dt:
20;%设置时间长度
pj=o;
fc=sqrt(0.01);%加入方差为0.01的白噪再序列
num=[l31];%分子多项式系数
den=[25.531];%分母多项式系数
%绘制原传递函数阶跃响应曲线fprintff原系统传递函数为J
G=tf(num,den)
L=length(t);
cs=randn(lzL);%产生噪声
cs=cs/std(cs);
cs=cs-mean(cs);
cs=pj+fc*cs;
cs=csf;
y=step(numzden,t)+cs;%加入噪声
Length=length(y);
Plot(t,y);
grid;
xlabelCt/s');
ylabelCy(ty);
%辨识程序设计,此系统m+n=5,故应计算A1-A5fprintf('Al-A5阶面积分别为:
J
%求A1
suml=O;
for(i=l:
Length-l)
suml=suml+(l-(y(i)+y(i+l))/2)*dt;A(i)=suml;
end
Al=suml
%求A2
sum2=0;
for(i=l:
Length-l)sum2=sum2+(A(i)-Al*(y(i)+y(i+l))/2)*dt;B(i)=sum2;
end
A2=sum2
%求A3
sum3=0;
for(i=l:
Length-l)sum3=sum3+(B(i)-A2*(y(i)+y(i+l))/2)*dt;C(i)=sum3;
end
A3=sum3
%求A4
sum4=0;
for(i=l:
Length-l)
sum4=sum4+(C(i)-A3*(y(i)+y(i+l))/2)*dt;
D(i)=sum4;
end
A4=sum4
%求A5sum5=0;
for(i=l:
Length-l)
sum5=sum5+(D(i)-A4*(y(i)+y(i+l))/2)*dt;
end
A5=sum5
%求分子系数bl,b2
M=(-l)*(inv([A3/A2;A4/A3]))*[A4;A5];
fprintff分子多项式系数为:
J
b2=M(2zl)
%求分母系数alza2,a3
N=[l0O;A110;A2All]*[bl;b2;0]+[Al;A2;A3];
fprintff分母多项式系数为:
J
31=^1?
1)
a2=N(2,l)
a3=N(3,l)
%求辨识所得传递函数
numl=[b2bl1);
denl=[a3a2al1];
fprintfC辨识所得传递函数为
G=tf(numl#denl)
h=step(numl/denl,t);
plotl^y/black'.th/blue*);
legendf原传递函数T辨识所得传递函数J;
title/原传递函数与辨识所得传递函数的阶跃响应对比Jgrid;
xlabel('t/s');
ylabeKMt)*);
fprintfC相关系数:
');%求相关系数
r=corrcoef(y,h)
运行程序,结果如下:
原系统传递函数为:
G=
sA2+3s+l
2sA3+5.5sA2+3s+l
Continuous-timetransferfunction.
A1-A5阶面积分别为:
Al=
0.0081
A2=
4.7477
A3=
-13.2541
A4=
39.7921
A5=
-121.2501
分子多项式系数为:
bl=
3.6415
b2=
1.7846
分母多项式系数为:
al=
3.6496
a2=
6.5619
a3=
4.0492
辨识所得传递函数为:
G=
1.785sA2+3.641s+1
4.049sA3+6.562sA2+3.65s+1Continuous-timetransferfunction.相关系数:
r=
1.00000.9076
0.90761.0000
此时原传递函数和辨识所得传递函数的阶跃响应对比如下图:
则当特征方程有I】个单根S1.S2,••-Sn时,传递函数可写成
当特征方程有mr个单根r阶重根罚时,传递函数可写成
G(s)=-^―+—^―+•••+«+邑±1+7+2=+…+_5}_(23)
S—qs—SnS—S^rs—So(s-So)-(S_So)T
为了确定G和Sx,从获取的过程输出脉冲响应g(t)中,选取前11+1个坐标点,每个坐标点间隔相同的采样时间To,组成一个自回归模型
g(k)+qg(4W---+aDg弦羽)(2.4)
其中a2Qn为待定系数。
如果特征方程
1+qx%+a2x7J*+•••+c^xnli=0(2.5)
有一个单根,则其必是AR模型的解,它们的线性组合也是AR模型的解。
当其特征方程有n个单根时,自回归模型的解为
(26)
2(k)=人坷咲+0?
x严+…+°屁崛
当其特征方程有n-r个单根I阶重根时,自回归模型的解为
g(k)=人严+…严+02朋呵+0“冋呵+・.・+处「弋耳……(2.7)
根据分析,无论特征方程是单根还是重根,都有
(2.8)
显然,一口求出%和»便可的传递函数。
2.1程序设计
选取系统传递函数如下:
〜、2s+1.05
G(S)s3+5.5s2+3s+l
MATLAB程序设计如下:
clccleardt=0.01;
%清空工作区
%设置采样时间
t=0:
dt:
50;%设置时间长度
num=[21.05];%分子多项式系数
den=[l5.531];%分母多项式系数
%绘制原传递函数阶跃响应曲线fprintff原系统传递函数为G=tf(numzden)
y=impulse(num,den,t);%脉冲响应曲线Length=length(y);
Plot(t,y);
grid;
xlabelCt/s');ylabeKMt)*);
%时刻为Os,Is2s3s4s5s时的脉冲响应,取T0=lfprintfCW刻为Oszls2s3s4s5s时的脉冲响应,取T0=l');
Y=[y(l/:
)y(101,:
)y(201,:
)y(301,:
)y(401,:
)y(501z:
)]
A=[Y(:
2:
4);Y(:
3:
5);Y(:
4:
6)];
B=Y(:
l:
3y;
fprintff系数a:
');a=(inv(A))*((-l)*B)p=[a(3/l)a(2,l)a(lzl)l];
%构造系数矩阵
%构造方程右边的数值
%求得待定系数a,在这里用a表示
%关于X的方程系数
x=roots(p);
fprintff系统极点:
J;
%求得xl,x乙x3
s=log(roots(p))
%求得传递函数极点
%构造系数矩阵,并求得p,这里用E表示
C=[l1l;x(lzl)x(2,l)x(3,l);(x(m2)(x(2,1))A2(x(3,l)八2)];D=[0Y(:
2)Y(:
3)]:
fprintff系数P');
E=(inv(C))*D
%绘制辨识所得传递函数的脉冲响应
[numl,denl]=residue(E,s,[]);
fprintff辨识所得传递函数');
Gl=tf(numl,denl)%传递函数
h^mpulsefnum^denlj);%脉冲响应曲线plot(t,y,'black:
t,h,'blue‘);
legendf原传递函数丁辨识所得传递函数J;titleC®系统与辨识所得系统脉冲响应对比J;grid;
xlabelCt/s');
ylabel('y(t)fnh(t)');
运行以上程序,得到结果如下:
原系统传递函数为:
G=
2S+1.05
sA3+5.5sA2+3s+1
Continuous-timetransferfunction.
时刻为Os,Is2s3s4s5s时的脉冲响应,取T0=l
Y=
00.36850.29560.20760.12610.0607
系数a:
a=
-141.2815
348.0135
-244.6919
系统极点:
s=
-0.2835+0.3498i
-0.2835・0.3498i
-4.9329
系数P
E=
0.2028-0.1637i
0.2028+0.1637i
-0.4055
辨识所得传递函数Warning:
Thenumeratorordenominatorofthistransferfunctionhascomplex-valuedcoefficients.
>Inwarningat26
Inat357
Inchafenfangchengfaat36
G1=
(2-l.lle-16i)s+(1.05-2.22e-16i)
sA3+5.5sA2+3s+l
Continuous-timetransferfunction.
Warning:
Thenumeratorordenominatorofthistransferfunctionhascomplex-valuedcoefficients・
>Inwarningat26
Intftf>tftfat357
Inimpulseat117
Inchafenfangchengfaat37
Warning:
ImaginarypartsofcomplexXand/orYargumentsignored
>Inchafenfangchengfaat38
»
图2.1原传递函数和辨识所得传递函数的阶跃响应对比
由上图和运行结果町以看出,采用差分方程法所得到的辨识结果也是比较精确的。
3总结
通过本次实验仿真,我掌握了而积法、差分方程法等经典辨识方法的一些基础编程实验方法,将所学理论知识运用在了实际操作中,对系统辨识的认识和掌握进一步加深。
在加入噪声的仿真中,由于噪声时随机序列,所以每次仿真的结果都不相同,并且加入噪声会使辨识结果准确度降低。
另外,采样时间选择不同,对辨识结果也有彫响。
- 配套讲稿:
如PPT文件的首页显示word图标,表示该PPT已包含配套word讲稿。双击word图标可打开word文档。
- 特殊限制:
部分文档作品中含有的国旗、国徽等图片,仅作为作品整体效果示例展示,禁止商用。设计者仅对作品中独创性部分享有著作权。
- 关 键 词:
- 系统 辨识 经典 方法