系统辨识之经典辨识法Word文档格式.docx

1、1 （𝑡）𝑎𝑑𝑡 （1-2） 面积法原则上可以求出n为任意阶的个系数。以n为3为例。有：3𝑦） 𝑑2𝑦 |𝑡 = = 0 （1-3） ）|𝑡 = 1将式（1）中的y（t）移至右边，在0，t上积分，得 3 （1-4） 定义： 𝐹1（𝑡） = 0𝑡1 𝑦）𝑑 （1-5） 由式（1-3）条件可知，当t时， （1-6） 同理，定义 𝐹2 （1-7） 由式（1-,3）

2、条件可知，当t时， （1-8） 因此，可得 𝐹1（𝑡） 𝑎1𝑦） dt （1-9） 𝑎 = 𝐹（） （1-10） 当系统的传递函数存在零点时，传递函数如下： G（s）=k b sm mn +ba smn-1-1smn-1-1 +LL +a sbs11 +1+1，（n m） （1-11） a sn +其中，K h= （ ） / U0 定义 1 G（s）=K P（s）其中， P（s） = b sa sn mn +ba smn-1-1smn-1-1 +LL +a sbs11 +1+1 = +1 i =1

3、C si i （1-12） m根据1h*（t）的Laplace变换，求出一阶面积A1，确定Lh（*1 t ），并定义二阶面积A2 ，以此类推，得到i 阶面积Ai 。进一步利用est 拉氏变换，得到L1h*（t ）=M si i ，进而得到Ai 的值：i=0Ai = 0 1h*（t） （i 1）!t）i1dt + tj=20 Ai j 1 0 1h*（t） （j!t） j dt （1-13） （根据A Ci = i ，可得： + 𝑎1 + + 𝑎 + 1 = （𝑏𝑚 + 𝑏1 + + 𝑏 + 1）（1

4、 + 𝑖=1 𝐴𝑖）。比较上式两边s的各次幂，便可得到a, b, A之间的关系，如下： b1 An An1 L An m +1 1 An+1 b2 An+1 An L An m +2 An+2 = M L L L L M bm An m+ 1 An m+ 2 L An An m+ b1 a1 1 10 LL 0 00 0 bM2 + AA12 （1-14） a2 = A1 M L L L L L M bm An an An1 An2 L A1 1 0 由此可知，根据式（1-12）、（1-13）、（1-14）便可得到辨识传递函数的参数a, b。1.

5、2实验过程 1.2.1无零点模型系统 假设系统的传递函数模型为G（s） = 2 1 ,为无零点的模型，利用10𝑠 +6.5𝑠+1Matlab 编程，分别在没有噪声和有噪声两种情况下进行辨识，比较辨识结果。1.没有噪声时，程序如下：clear;%=获得原传递函数方程=% num=1; den=10 6.5 1;%=产生阶跃采样序列=% T=0.2; %采样周期 t=0:T:30; %采样时间 L=length（t）; %数据长度 h=step（num,den,t）; %原传递函数的阶跃响应 K=h（L） %系统增益 %=面积法求解参数=% s1=0; for i=

6、1:L s1=s1+（1-h（i）*T; F（i）=s1; end a1=s1; s2=0;for i=1:L s2=s2+（F（i）-a1*h（i）*T; end a2=s2; num1=1; den1=a2 a1 1;disp（原传递函数为:） G1=tf（num,den） 通过辨识得到的传递函数为：G2=tf（num1,den1） %=原传递函数和辨识函数的阶跃响应对比图=% step（G1,b-,G2,r-.title（原系统与辨识后所得到系统阶跃响应对比） legend（原响应曲线,辨识响应曲线（1） 当采样周期T=0.2秒，采样时间t=30s时，行程序后得到原传递函数G1和辨识得到

7、的传递函数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

8、; 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 估计值50 1.0000 11.31 6.6 0.02 10.13 6.51 有测量噪 （方差为 0.01）

9、 1 7.784 6.709 分析：通过对比不同的采样周期和不同的采样时间在无测量噪声情况的辨识结果可知，在相同的采样周期下，适当的增加采样时间，可以提高辨识精度，尤其是对增益的提高有很大影响；而在相同的采样时间下，适当的减小采样时间，对于系统参数的辨识精度有很大的提高。因此，可以发现合理采样时间和数据长度，可以提高辨识的精度，令辨识后的传递函数系数与原传递函数系数更接近，差距小，从而得到满意的辨识结果。通过对比无测量噪声和有测量噪声两种情况下的辨识结果，我们可以发现在白噪声的情况下，曲线拟合较无噪声情况下要差，说明白噪声对于面积法辨识系统存在较大的干扰，会对辨识结果产生一定的影响。1.2.2

10、有零点模型系统 17.5𝑆2+7.5𝑆+1 假设系统的传递函数为G（s） = 4𝑆3+5𝑆2+8𝑆+1,为有零点的模型，其中n=3,m=2, 用面积法需要求解𝐴1𝐴5，利用 Matlab 编程，分别在没有噪声和有噪声两种情况下进行辨识，比较辨识结果。（1） 没有噪声时，程序如下：%=获得原传递函数方程=% num=17.5 7.5 1; den=4 5 8 1;%=产生阶跃采样序列=% T=0.02;100; %数据长度 y=step（num,den,t）; k=y（L） %系统增益

11、 %=面积法求解参数=% sum1=0;L-1; sum1=sum1+（1-（y（i）+y（i+1）/2）*T; A（i）=sum1; end A1=sum1 sum2=0; sum2=sum2+（A（i）-A1*（y（i）+y（i+1）/2）*T; B（i）=sum2; end A2=sum2 sum3=0; sum3=sum3+（B（i）-A2*（y（i）+y（i+1）/2）*T; C（i）=sum3; end A3=sum3 sum4=0; sum4=sum4+（C（i）-A3*（y（i）+y（i+1）/2）*T; D（i）=sum4; end A4=sum4 sum5=0; sum5=

12、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=1 0 0;A1 1 0;A2 A1 1*b1;b2;0+A1;A2;A3; a1=N（1,1）; a2=N（2,1）; a3=N（3,1）;%=根据所求a、b，得到辨识后传递函数=% num1=b2 b1 1; den1=a3 a2 a1 1;%=原传递函数和辨识函数的阶跃响应对比图=% step（G1, 当采样时间取0.02秒，数据长度为100时，

13、辨识结果如下：图1.10 原系统与辨识后的系统阶跃响应对比图：当采样时间为0.02，数据长度为400时，系统辨识结果如下：图1-12 图1-13 当采样时间为0.2秒，数据长度为400时，系统辨识结果如下：图1-14 原系统与辨识后的系统阶跃响应对比图：综上所述，结果如表 表1-2 a 3 a 2 a 1 b 2 b 1 4 5 8 17.5 7.5 100 3.559 5.224 8.051 17.71 7.515 400 4.016 4.914 7.979 17.42 7.479 0. 2 4.129 4.168 7.792 16.61 7.278 通过对比不同的采样周期和不同的采样时间在

14、无测量噪声情况的辨识结果可知，对于存在有零点的系统来说，通过面积法辨识系统必须合理的选择分子分母的阶次，否则不能得出正确的辨识结果。从表格中也可以发现，在相同的采样时间下，由于系统收敛过程较长从，增加数据长度对系统的辨识精度有了很大的提高；同样在相同的数据长度下，合理的减小采样时间，也可以提高系统的辨识精度。2.脉冲响应法 脉冲响应法（ impulse response method ）是指线性系统在零初始条件输入脉冲信号，信号后输出的瞬态响应，即输出响应叫脉冲响应。建立系统的非参数模来观测系统脉冲响应，以求得系统数学模型的待定参数，进而实现系统辨识的目的。2.1 Hankel矩阵法 一个 n

15、 阶过程的脉冲传递函数为 G（z ）=1 1 1+b za z2222+LL+b zna znnn （2-1） b z 1 1+a z1 +将传递函数转化成为状态方程后，进一步推导，可知 G（z ）=g（1）z1 1+g（2）z2 +g（3）z3 +L 根据上述公式，可得 bz1 1 +b z2 2 +L +b zn n= g（1）z1 +g（2）z2 +g（3）z3 +L 1+a z1 1 +a z2 2 +L +a zn n = g（1）z1 +L + g（2）+ag1 （1） z2 +L + g（n）+ n1 g（i）an i zn （2-2） i=1 n 2n1+ g（n+ +1） g

16、（i）an+i z +（n 1） +L + g（2n）+ g（i）a2n i z2n i=1 i=1 脉冲相应法参数计算公式如下： g（1） g（2） L g n（ ） an g n（ +1） g（2） g（3） L g n（ +1） an1 = g n（ +2） （2-3） M M O M M M g n（ ） g n（ +!） L g n（2 1） a1 g n（2 ） b1 1 0 L 0 0 g（1） b2 = a1 1 L 0 0 g（2） M M M O M M M bn an1 an2 L a1 1 g n（ ） （2-4） Hankel矩阵的定义： 𝑔（

17、9896;） 𝑔 + 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.

18、实验过程 假设系统模型的传递函数为G（s） = 𝑠3+155𝑠2𝑠+61𝑠+100, 利用Matlab编程，分别在没有噪声和有噪声两种情况下进行辨识，比较辨识结果。（1）无噪声的情况下，程序如下： clear %=获得原传递函数方程=% num=5 1; den=1 15 50 100;%=产生脉冲采样序列=% T=0.02; %采样时间 gk=impulse（num,den,t）; %原传递函数的脉冲响应 %=Hamkel矩阵求解参数=% H=gk（1） gk（2） gk（3）;gk（2） gk（3） gk（4）;gk（3） g

19、k（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=1 0 0;a（1） 1 0;a（2） a（1） 1; %构造下三角阵 G2=gk（1）;gk（2）;gk（3）; b=A*G2; a=1 a;原传递函数及其Z函数dsys1=c2d（G1,T） 辨识得到的Z函数和传递函数dsys2=filt（b.,a,0.02） %求得辨识得到的Z函数 G2=T*d2c（dsys2,tustin,0.02） %双线

20、性变换 %=原传递函数和辨识函数的脉冲响应对比图=% impulse（G1,） %绘出原系统与辨识得到的系统的脉冲响应曲线 原模型与辨识后所得到系统脉冲响应对比原始响应曲线辨识曲线） 程序运行结果如下：原系统的传递函数和脉冲阶跃函数如下：图2.1 辨识之后系统的传递函数和脉冲函数如下：图2.2 在没有噪声的情况下，原系统与辨识系统脉冲响应对比图：图2.3 （2）在噪声期望为 0，方差为 0.01 的情况下，对系统进行参数辨识，在 Mtlab中编程实现，程序如下：clear %=获得原传递函数方程=% num=5 1;L=length（t）;%产生均值为0，方差为0.01的白噪声 w=randn

21、（1,L）; w=w/std（w）;%=阶跃采样序列中加入白噪声=% gk=gk+w;%=Hamkel矩阵求解参数=% H=gk（1） gk（2） gk（3）;%构造Hankel a=H1*G1; %a中的值倒转，a=a1；a2；a3 G1=tf（num,den） dsys1=c2d（G1,T） %=原传递函数和辨识函数的脉冲响应对比图=% impulse（G1,） %绘出原系统与辨识得到的系统的脉冲响应曲线 title（图2.4 原系统与辨识之后的系统脉冲响应图如下：2.3结果分析 通过对比在没有噪声情况下原传递函数和辨识后系统的脉冲响应，我们发现两者的脉冲响应较为接近，但传递函数有一定差别。而在加入白噪声的情况下，辨识后的系统脉冲响应与原传递函数的脉冲响应有较大的差别，存在很大的不确