系统辨识之经典辨识法.docx
《系统辨识之经典辨识法.docx》由会员分享，可在线阅读，更多相关《系统辨识之经典辨识法.docx（20页珍藏版）》请在冰豆网上搜索。
系统辨识之经典辨识法
系统辨识作业一
学院信息科学与工程学院
专业控制科学与工程
班级控制二班
姓名
学号
2018年11月
系统辨识
所谓辨识就是通过测取研究对象在认为输入作用的输出响应，或正常运行时
的输入输出数据记录，加以必要的数据处理和数学计算，估计出对象的数学模型。
辨识的内容主要包括四个方面：
①实验设计；
②模型结构辨识；
③模型参数辨识；
④模型检验。
辨识的一般步骤：
根据辨识目的，利用先验知识，初步确定模型结构；采集数据；然后进行模型参数和结构辨识；最终验证获得的最终模型。
根据辨识方法所涉及的模型形式来说，辨识方法可以分为两类：
一类是非参数模型辨识方法，另一类是参数模型辨识方法。
其中，非参数模型辨识方法又称为经典的辨识方法，它主要获得的是模型是非参数模型。
在假定过程是线性的前提下，不必事先确定模型的具体结构，广泛适用于一些复杂的过程。
经典辨识方法有很多，其中包括阶跃响应法、脉冲响应法、相关分析法和普分析法等等，本次实验所采用的辨识方法为阶跃响应法和脉冲响应法。
1．阶跃响应法
阶跃响应法是一种常用非参数模型辨识方法。
常用的方法有近似法、半对数法、切线法、两点法和面积法等。
本次作业采用面积法求传递函数。
1.1面积法
①当系统的传递函数无零点时，即系统传递函数如下：
G（S）=𝑎
𝑛𝑠𝑛+𝑎𝑛−1𝑠𝑛1−1+⋯+𝑎1𝑠+1（1-1）系统的传递函数与微分方程存在着一一对应的关系，因此，可以通过求取微分方程的系数来辨识系统的传递函数。
在求得系统的放大倍数K后，要得到无因次阶跃响应y（t）（设τ=0），其中y（t）用下式描述：
𝑑𝑛𝑦（𝑡）𝑛−1（𝑡）
𝑎𝑛
𝑑𝑡𝑑𝑡𝑑𝑡（1-2）面积法原则上可以求出n为任意阶的个系数。
以n为3为例。
有：
𝑑3𝑦（𝑡）𝑑2𝑦（𝑡）𝑑𝑦（𝑡）
{𝑑𝑡|𝑡→∞=𝑑𝑡|𝑡→∞=𝑑𝑡|𝑡→∞=0（1-3）
𝑦（𝑡）|𝑡→∞=1
将式
（1）中的y（t）移至右边，在[0，t]上积分，得
𝑑2𝑦（𝑡）
𝑎3
𝑑𝑡𝑑𝑡（1-4）定义：
𝐹1（𝑡）=∫0𝑡[1−𝑦（𝑡）]𝑑𝑡（1-5）由式（1-3）条件可知，当t→∞时，
𝑎
𝑑𝑡（1-6）
同理，定义
𝐹2
𝑑𝑡（1-7）
由式（1-,3）条件可知，当t→∞时，
𝑎
𝑑𝑡（1-8）
因此，可得
𝐹𝑛（𝑡）=∫0𝑡[𝐹𝑛−1（𝑡）−𝑎𝑛−1𝑦（𝑡）]dt（1-9）
𝑎𝑛=𝐹𝑛（∞）（1-10）
②当系统的传递函数存在零点时，传递函数如下：
G（s）=k
bsmmn+basmn-1-1smn-1-1++LL++asbs11+1+1，（nm）（1-11）
asn+
其中，Kh=（）/U0定义
1
G（s）=K
P（s）
其中，
P（s）=
bsasnmn++basmn-1-1smn-1-1++LL++asbs11+1+1=+1i=1Csii（1-12）m
根据[1−h*（t）]的Laplace变换，求出一阶面积A1，确定L[h（*1t]），并定义二
阶面积A2，以此类推，得到i阶面积Ai。
进一步利用e−st拉氏变换，得到
L[1−h*（t]）=Msii，进而得到Ai的值：
i=0
Ai=01−h*（t）（i1）!
−−t）i−1dt+tj−=20Ai−−j101−h*（t）
（−j!
t）jdt（1-13）（
根据ACi=i，可得：
𝑎𝑛𝑠𝑛+𝑎𝑛−1𝑠𝑛−1+⋯+𝑎1𝑠+1
=（𝑏𝑚𝑠𝑚+𝑏𝑚−1𝑠𝑚−1+⋯+𝑏1𝑠+1）（1+∑∞𝑖=1𝐴𝑖𝑠𝑖）。
比较上式两边s的各次幂，便可得到a,b,A之间的关系，如下：
b1AnAn−1LAnm−+1−1An+1
b2An+1AnLAnm−+2An+2=−
MLLLLM
bmAnm+−1Anm+−2LAnAnm+
b1
a1110LL0000bM2+AA12（1-14）
a2=A1
MLLLLLM
bmAn
anAn−1An−2LA110
由此可知，根据式（1-12）、（1-13）、（1-14）便可得到辨识传递函数的参数a,b。
1.2实验过程1.2.1无零点模型系统假设系统的传递函数模型为G（s）=
21,为无零点的模型，利用
10𝑠+6.5𝑠+1
Matlab编程，分别在没有噪声和有噪声两种情况下进行辨识，比较辨识结果。
1.没有噪声时，程序如下：
clear;
%==================获得原传递函数方程=======================%num=[1];den=[106.51];
%=====================产生阶跃采样序列======================%T=0.2;%采样周期t=0:
T:
30;%采样时间L=length（t）;%数据长度h=step（num,den,t）;%原传递函数的阶跃响应
K=h（L）%系统增益
%======================面积法求解参数======================%s1=0;fori=1:
Ls1=s1+（1-h（i））*T;F（i）=s1;enda1=s1;s2=0;
fori=1:
Ls2=s2+（F（i）-a1*h（i））*T;enda2=s2;num1=[1];den1=[a2a11];
disp（'原传递函数为:
'）
G1=tf（num,den）
disp（'通过辨识得到的传递函数为：
'）
G2=tf（num1,den1）
%=============原传递函数和辨识函数的阶跃响应对比图=============%step（G1,'b-',G2,'r-.'）
title（'原系统与辨识后所得到系统阶跃响应对比'）legend（'原响应曲线','辨识响应曲线'）
（1）当采样周期T=0.2秒，采样时间t=30s时，行程序后得到原传递函数G1和辨识得到的传递函数G2如图1.1：
图1.1
原系统和辨识后系统的阶跃响应对比图如下：
图1.2
（2）当采样周期T=0.2秒，采样时间t=50s时，行程序后得到原传递函数G1和辨识得到的传递函数G2如下：
图1.3
原系统和辨识后系统的阶跃响应对比图如下：
（3）当采样周期T=0.02秒，采样时间t=50s时，行程序后得到原传递函数G1和辨识得到的传递函数G2如下：
图1.5
原系统和辨识后系统的阶跃响应对比图如下：
2.有噪声的情况下，系统程序如下：
主程序还是用面积法，在程序中添加以下代码：
%产生期望为0，方差为0.01的噪声
figure
（1）w=randn（1,L）;%建立服从正态分布的随机矩阵。
w=w/std（w）;w=w-mean（w）;qw=0;fc=0.01;w=qw+sqrt（fc）*w;
%=====================阶跃采样序列中加入白噪声==================%h=h+w;plot（t,w）;
（1）加入的噪声如下图所示：
图1.7
（2）当采样周期T=0.02s,采样时间t=50s时，辨识结果如下：
图1.8
原系统与辨识系统阶跃响应如图所示：
结合上述无测量噪声和有测量噪声两种情况下的辨识结果，列出如下所示的表格：
表1-1
噪声情况
条件
增益
a1
a2
参数
采样时间
数据长度
1.0
10
6.5
真值
无测量噪声
0.2
30
0.9985
11.43
6.594
估计值
0.2
50
1.0000
11.31
6.6
0.02
50
1.0000
10.13
6.51
有测量噪
（方差为0.01）
0.02
50
1
7.784
6.709
分析：
通过对比不同的采样周期和不同的采样时间在无测量噪声情况的辨识结果可知，在相同的采样周期下，适当的增加采样时间，可以提高辨识精度，尤其是对增益的提高有很大影响；而在相同的采样时间下，适当的减小采样时间，对于系统参数的辨识精度有很大的提高。
因此，可以发现合理采样时间和数据长度，可以提高辨识的精度，令辨识后的传递函数系数与原传递函数系数更接近，差距小，从而得到满意的辨识结果。
通过对比无测量噪声和有测量噪声两种情况下的辨识结果，我们可以发现在白噪声的情况下，曲线拟合较无噪声情况下要差，说明白噪声对于面积法辨识系统存在较大的干扰，会对辨识结果产生一定的影响。
1.2.2有零点模型系统
17.5𝑆2+7.5𝑆+1假设系统的传递函数为G（s）=4
𝑆3+5𝑆2+8𝑆+1,为有零点的模型，其中n=3,m=2,用面积法需要求解𝐴1~𝐴5，利用Matlab编程，分别在没有噪声和有噪声两种情况下进行辨识，比较辨识结果。
（1）没有噪声时，程序如下：
clear;
%==================获得原传递函数方程=======================%num=[17.57.51];den=[4581];
%=====================产生阶跃采样序列======================%T=0.02;%采样周期t=0:
T:
100;%采样时间L=length（t）;%数据长度y=step（num,den,t）;k=y（L）%系统增益
%======================面积法求解参数======================%sum1=0;fori=1:
L-1;
sum1=sum1+（1-（y（i）+y（i+1））/2）*T;A（i）=sum1;endA1=sum1sum2=0;fori=1:
L-1;
sum2=sum2+（A（i）-A1*（y（i）+y（i+1））/2）*T;B（i）=sum2;endA2=sum2sum3=0;fori=1:
L-1;
sum3=sum3+（B（i）-A2*（y（i）+y（i+1））/2）*T;C（i）=sum3;endA3=sum3sum4=0;
fori=1:
L-1;
sum4=sum4+（C（i）-A3*（y（i）+y（i+1））/2）*T;D（i）=sum4;endA4=sum4sum5=0;fori=1:
L-1;sum5=sum5+（D（i）-A4*（y（i）+y（i+1））/2）*T;end
A5=sum5
%==============根据所得A（i），利用公式求取a、b=================%
M=（-1）*（inv（[A3,A2;A4,A3]））*[A4;A5];b1=M（1,1）;b2=M（2,1）;
N=[100;A110;A2A11]*[b1;b2;0]+[A1;A2;A3];a1=N（1,1）;a2=N（2,1）;a3=N（3,1）;
%================根据所求a、b，得到辨识后传递函数==============%num1=[b2b11];den1=[a3a2a11];
disp（'原传递函数为:
'）
G1=tf（num,den）
disp（'通过辨识得到的传递函数为：
'）
G2=tf（num1,den1）
%=============原传递函数和辨识函数的阶跃响应对比图=============%
step（G1,'b-',G2,'r-.'）
title（'原系统与辨识后所得到系统阶跃响应对比'）legend（'原响应曲线','辨识响应曲线'）
当采样时间取0.02秒，数据长度为100时，辨识结果如下：
图1.10
原系统与辨识后的系统阶跃响应对比图：
当采样时间为0.02，数据长度为400时，系统辨识结果如下：
图1-12
原系统与辨识后的系统阶跃响应对比图：
图1-13
当采样时间为0.2秒，数据长度为400时，系统辨识结果如下：
图1-14原系统与辨识后的系统阶跃响应对比图：
综上所述，结果如表
表1-2
条件
增益
a3
a2
a1
b2
b1
参数
采样时间
数据长度
1.0000
4
5
8
17.5
7.5
真值
0.02
100
1.0000
3.559
5.224
8.051
17.71
7.515
估计值
0.02
400
1.0000
4.016
4.914
7.979
17.42
7.479
0.2
400
1.0000
4.129
4.168
7.792
16.61
7.278
分析：
通过对比不同的采样周期和不同的采样时间在无测量噪声情况的辨识结果可知，对于存在有零点的系统来说，通过面积法辨识系统必须合理的选择分子分母的阶次，否则不能得出正确的辨识结果。
从表格中也可以发现，在相同的采样时间下，由于系统收敛过程较长从，增加数据长度对系统的辨识精度有了很大的提高；同样在相同的数据长度下，合理的减小采样时间，也可以提高系统的辨识精度。
2.脉冲响应法
脉冲响应法（impulseresponsemethod）是指线性系统在零初始条件输入脉冲信号，信号后输出的瞬态响应，即输出响应叫脉冲响应。
建立系统的非参数模来观测系统脉冲响应，以求得系统数学模型的待定参数，进而实现系统辨识的目的。
2.1Hankel矩阵法
一个n阶过程的脉冲传递函数为
G（z）=−1
1−−1+bzaz22−2−2++LL++bznazn−n−n（2-1）bz1
1+az1+
将传递函数转化成为状态方程后，进一步推导，可知
G（z）=g
（1）z−1−1+g
（2）z−2+g（3）z−3+L
根据上述公式，可得
bz1−1+bz2−2+L+bzn−n
=g
（1）z−1+g
（2）z−2+g（3）z−3+L1+az1−1+az2−2+L+azn−n
=g
（1）z−1+L+g
（2）+ag1
（1）z−2+L+g（n）+n−1g（i）ani−z−n（2-2）
i=1
n2n−1
+g（n++1）g（i）an+−iz−+（n1）+L+g（2n）+g（i）a2ni−z−2n
i=1i=1脉冲相应法参数计算公式如下：
g
（1）g
（2）Lgn（）angn（+1）
g
（2）g（3）Lgn（+1）an−1=−gn（+2）（2-3）
MMOMMM
gn（）gn（+!
）Lgn（2−1）a1gn
（2）
b110L00g
（1）
b2=a11L00g
（2）
MMMOMMM
bnan−1an−2La11gn（）（2-4）
Hankel矩阵的定义：
𝑔（𝑘）𝑔（𝑘+1）𝑔（𝑘+𝑙−1）
⋯
𝑔（𝑘+1）𝑔（𝑘+2）𝑔（𝑘+𝑙）
H（l,k）=（）（2-5）
⋮⋱⋮
𝑔（𝑘+𝑙−1）𝑔（𝑘+𝑙）⋯𝑔（𝑘+2𝑙−2）
公式（2-3）左边的矩阵是一种特定的Hankel矩阵。
记作H（n,l）,并且它是可逆的。
因此，只需要将获得的脉冲响应值g（k），k=1,2,3,…,2n填入式（3-3）、
（3-4）便可以求出脉冲传递函数的估计值，即传递函数系数：
𝑎1,𝑎2…,𝑎𝑛,𝑏1,𝑏2…,𝑏𝑛。
2.2.实验过程
假设系统模型的传递函数为G（s）=𝑠
3+155𝑠2𝑠++61𝑠+100,利用Matlab编程，分别在没有噪声和有噪声两种情况下进行辨识，比较辨识结果。
（1）无噪声的情况下，程序如下：
clear
%==================获得原传递函数方程=======================%
num=[51];den=[11550100];
%=====================产生脉冲采样序列======================%T=0.02;%采样周期t=0:
T:
30;%采样时间gk=impulse（num,den,t）;%原传递函数的脉冲响应
%=====================Hamkel矩阵求解参数=====================%H=[gk
（1）gk
（2）gk（3）;gk
（2）gk（3）gk（4）;gk（3）gk（4）gk（5）];%构造Hankel
矩阵求a、b系数矩阵
H1=H^-1;%H求逆G1=-[gk（4）;gk（5）;gk（6）];
a=H1*G1;%这里a=[a3;a2;a1]a=flipud（a）.';%a中的值倒转，a=[a1,a2,a3]
A=[100;a
（1）10;a
（2）a
（1）1];%构造下三角阵G2=[gk
（1）;gk
（2）;gk（3）];b=A*G2;a=[1a];
disp（'原传递函数及其Z函数'）
G1=tf（num,den）
dsys1=c2d（G1,T）
disp（'辨识得到的Z函数和传递函数'）
dsys2=filt（b.',a,0.02）%求得辨识得到的Z函数
G2=T*d2c（dsys2,'tustin',0.02）%双线性变换
%=============原传递函数和辨识函数的脉冲响应对比图=============%
impulse（G1,'b-',G2,'r-.'）%绘出原系统与辨识得到的系统的脉冲
响应曲线
title（'原模型与辨识后所得到系统脉冲响应对比'）legend（'原始响应曲线','辨识曲线'）程序运行结果如下：
原系统的传递函数和脉冲阶跃函数如下：
图2.1
辨识之后系统的传递函数和脉冲函数如下：
图2.2
在没有噪声的情况下，原系统与辨识系统脉冲响应对比图：
图2.3
（2）在噪声期望为0，方差为0.01的情况下，对系统进行参数辨识，在Mtlab中编程实现，程序如下：
clear
%==================获得原传递函数方程=======================%num=[51];den=[11550100];
%=====================产生脉冲采样序列======================%T=0.02;%采样周期t=0:
T:
30;%采样时间gk=impulse（num,den,t）;%原传递函数的脉冲响应
L=length（t）;
%产生均值为0，方差为0.01的白噪声
w=randn（1,L）;w=w/std（w）;w=w-mean（w）;qw=0;fc=0.01;w=qw+sqrt（fc）*w;
%=====================阶跃采样序列中加入白噪声==================%gk=gk+w;
%=====================Hamkel矩阵求解参数=====================%
H=[gk
（1）gk
（2）gk（3）;gk
（2）gk（3）gk（4）;gk（3）gk（4）gk（5）];%构造Hankel
H1=H^-1;%H求逆G1=-[gk（4）;gk（5）;gk（6）];a=H1*G1;%这里a=[a3;a2;a1]a=flipud（a）.';%a中的值倒转，a=[a1；a2；a3]
A=[100;a
（1）10;a
（2）a
（1）1];%构造下三角阵G2=[gk
（1）;gk
（2）;gk（3）];b=A*G2;a=[1a];
disp（'原传递函数及其Z函数'）
G1=tf（num,den）dsys1=c2d（G1,T）
disp（'辨识得到的Z函数和传递函数'）
dsys2=filt（b.',a,0.02）%求得辨识得到的Z函数
G2=T*d2c（dsys2,'tustin',0.02）%双线性变换
%=============原传递函数和辨识函数的脉冲响应对比图=============%impulse（G1,'b-',G2,'r-.'）%绘出原系统与辨识得到的系统的脉冲响应曲线title（'原模型与辨识后所得到系统脉冲响应对比'）legend（'原始响应曲线','辨识曲线'）程序运行结果如下：
辨识之后系统的传递函数和脉冲函数如下：
图2.4
原系统与辨识之后的系统脉冲响应图如下：
2.3结果分析
通过对比在没有噪声情况下原传递函数和辨识后系统的脉冲响应，我们发现两者的脉冲响应较为接近，但传递函数有一定差别。
而在加入白噪声的情况下，辨识后的系统脉冲响应与原传递函数的脉冲响应有较大的差别，存在很大
升级会员 