系统辨识经典辨识方法.doc

1、系统辨识 经典辨识方法经典辨识方法报告1. 面积法1.1 辨识原理1.1.1 分子多项式为 1 的系统（1.1）)( 11sasasGnn由于系统的传递函数与微分方程存在着一一对应的关系，因此，可以通过求取微分方程的系数来辨识系统的传递函数。在求得系统的放大倍数 K 后，要先得到无因次阶跃响应 y(t)（设 =0） 。大多数自衡的工业过程对象的 y(t)可以用下式描述来近似（1.2）1)()()(a11 tydadtydtynn面积法原则上可以求出 n 为任意阶的各系数。以 n=3 为例，注意到（1.3）1|)(,0|)(d|)(|)(d23 tttt yyy将式（2.1.2）的 y(t)项移

2、至右边，在 0，t 上积分，得（1.4）tdtyyadttya01223 )()()(定义（1.5）tdtyF01)()(则由式（2.1.3）给出的条件可知，在 t01)(adty（1.6）将式 a1y(t)移到等式右边，定义系统辨识 经典辨识方法（1.7）)()()(a 20123 tFdtyatFtydtt 利用初始条件（2.1.3）当 t时 )(a2（1.8）同理有 a3=F3()以此类推，若 n2，有 an=Fn()1.1.2 分子、分母分别为 m 阶和 n 阶多项式的系统当传递函数的形式如下所示时(1.9)011)( )(uhKsasabbsGnnm 定义(1.10)111)()(

3、iimnn scsbsbaasPG由于 0*)(1)(dtehthLs(1.11)则 的 Laplace 变换为：)(1*th(1.12)11*)()( iiiisCsPtL定义一阶面积 为：1A系统辨识 经典辨识方法(1.13)110001 lim)(*1li)(* csCthLdthAiiiiss 令(1.14)1()*1scthL定义二阶面积为：(1.15)2120*012 )()(limcscdthAiiiiist 同理，令 (1.16).1() 12*1 ii scscsthL定义 阶面积为 。由此可得：iiiA(1.17) 11111 )(.(. iimmnn sAsbsbsasa

4、上式可写成如下形式：(1.18) mnnnmnn mnnnn AAb 2121 21 1121系统辨识 经典辨识方法(1.19) nmnn AbAa 2121121121 00通过该系数矩阵，即可求出传递函数分子分母系数的值。1.2 程序设计1.2.1 传递函数形式如式 1.1 的系统取系统传递函数如下：()= 133+22+5.2+1MATLAB 程序如下：clc %清空工作区cleardt=0.01; %设置采样时间t=0:dt:50; %设置时间长度num=1; %此系统分子为 1den=3 2 5.2 1; %分母多项式系数%绘制原传递函数阶跃响应曲线fprintf(原系统传递函数为

5、:)G=tf(num,den)y=step(num,den,t);Length=length(y); %数据长度plot(t,y);grid;xlabel(t/s);ylabel(y(t);%进行辨识设计fprintf(辨识参数结果：);%求 a1sum1=0;for(i=1:Length)sum1=sum1+(1-y(i)*dt;F(i)=sum1;enda1=sum1%求 a2系统辨识 经典辨识方法sum2=0;for(i=1:Length)sum2=sum2+(F(i)-a1*y(i)*dt;f(i)=sum2;enda2=sum2%求 a3sum3=0;for(i=1:Length)s

6、um3=sum3+(f(i)-a2*y(i)*dt;enda3=sum3%绘制辨识后的传递函数dt=0.01;t=0:dt:50;num2=1;den2=a3 a2 a1 1;fprintf(系统辨识后的传递函数为 :)G=tf(num2,den2)h=step(num2,den2,t); %辨识所得传递函数阶跃响应plot(t,y,black,t,h,blue);legend(原传递函数 ,辨识所得传递函数);title(原传递函数与辨识所得传递函数的阶跃响应对比)grid;xlabel(t/s);ylabel(y(t)和 h(t);fprintf(相关系数：); %求相关系数r=corrc

7、oef(y,h)运行以上程序得到结果如下：原系统传递函数为:G =1-3 s3 + 2 s2 + 5.2 s + 1Continuous-time transfer function.辨识参数结果：a1 =系统辨识 经典辨识方法5.2048a2 =2.0608a3 =2.8388系统辨识后的传递函数为:G =1-2.839 s3 + 2.061 s2 + 5.205 s + 1Continuous-time transfer function.相关系数：r =1.0000 0.99990.9999 1.0000此时原传递函数和辨识所得传递函数的阶跃响应对比如下图：系统辨识 经典辨识方法0 5

8、10 15 20 25 30 35 40 45 5000.10.20.30.40.50.60.70.80.91t/sy(t)与h(t)与 与 与 与 与 与 与 与 与 与 与 与 与 与 与 与 与 与 与 与 与与 与 与 与 与与 与 与 与 与 与 与 与图 1.1 原传递函数和辨识所得传递函数的阶跃响应对比由上图可以看出，辨识所得结果比较准确。1.2.1 传递函数形式如式(1.9)的系统（无噪声）取系统传递函数如下：()= 2+3+123+5.52+3+1MATLAB 程序如下：clc %清空工作区cleardt=0.01; %设置采样时间t=0:dt:20; %设置时间长度num=

9、1 3 1; %分子多项式系数den=2 5.5 3 1; %分母多项式系数%绘制原传递函数阶跃响应曲线fprintf(原系统传递函数为 :)G=tf(num,den)y=step(num,den,t);Length=length(y); plot(t,y);grid;xlabel(t/s);ylabel(y(t);系统辨识 经典辨识方法%辨识程序设计，此系统 m+n=5,故应计算 A1-A5fprintf(A1-A5 阶面积分别为： )%求 A1sum1=0;for(i=1:Length-1)sum1=sum1+(1-(y(i)+y(i+1)/2)*dt;A(i)=sum1;endA1=su

10、m1%求 A2sum2=0;for(i=1:Length-1)sum2=sum2+(A(i)-A1*(y(i)+y(i+1)/2)*dt;B(i)=sum2;endA2=sum2%求 A3sum3=0;for(i=1:Length-1)sum3=sum3+(B(i)-A2*(y(i)+y(i+1)/2)*dt;C(i)=sum3;endA3=sum3%求 A4sum4=0;for(i=1:Length-1)sum4=sum4+(C(i)-A3*(y(i)+y(i+1)/2)*dt;D(i)=sum4;endA4=sum4%求 A5sum5=0;for(i=1:Length-1)sum5=sum

11、5+(D(i)-A4*(y(i)+y(i+1)/2)*dt;endA5=sum5%求分子系数 b1,b2M=(-1)*(inv(A3,A2;A4,A3)*A4;A5;fprintf(分子多项式系数为： )系统辨识 经典辨识方法b1=M(1,1)b2=M(2,1)%求分母系数 a1,a2,a3N=1 0 0;A1 1 0;A2 A1 1*b1;b2;0+A1;A2;A3;fprintf(分母多项式系数为： )a1=N(1,1)a2=N(2,1)a3=N(3,1)%求辨识所得传递函数num1=b2 b1 1;den1=a3 a2 a1 1;fprintf(辨识所得传递函数为 :)G=tf(num1

12、,den1)h=step(num1,den1,t);plot(t,y,black,t,h,blue);legend(原传递函数 ,辨识所得传递函数);title(原传递函数与辨识所得传递函数的阶跃响应对比)grid;xlabel(t/s);ylabel(y(t);fprintf(相关系数：); %求相关系数r=corrcoef(y,h)运行以上程序结果如下：原系统传递函数为:G =s2 + 3 s + 1-2 s3 + 5.5 s2 + 3 s + 1Continuous-time transfer function.A1-A5 阶面积分别为：A1 =0.0076A2 =4.3364系统辨识

13、经典辨识方法A3 =-9.6285A4 =15.7297A5 =7.1087分子多项式系数为：b1 =7.4409b2 =12.8942分母多项式系数为：a1 =7.4485a2 =17.2869a3 =22.7359辨识所得传递函数为:G =12.89 s2 + 7.441 s + 1-22.74 s3 + 17.29 s2 + 7.448 s + 1Continuous-time transfer function.相关系数：r =1.0000 0.99960.9996 1.0000此时原传递函数和辨识所得传递函数的阶跃响应对比如下图：系统辨识 经典辨识方法0 2 4 6 8 10 12

14、14 16 18 2000.20.40.60.811.21.4t/sy(t)与 与 与 与 与 与 与 与 与 与 与 与 与 与 与 与 与 与 与 与 与与 与 与 与 与与 与 与 与 与 与 与 与图 1.2 原传递函数和辨识所得传递函数的阶跃响应对比由上图可以看出，在未加入噪声之前，采用面积法辨识结果很精确，并且，分子可以为阶次低于分母的任意阶次。1.2.2 加入噪声为便于对比，仍取系统传递函数如下：()= 2+3+123+5.52+3+1MATLAB 程序如下：clc %清空工作区cleardt=0.01; %设置采样时间t=0:dt:20; %设置时间长度pj=0; fc=sqr

15、t(0.01); %加入方差为 0.01 的白噪声序列num=1 3 1; %分子多项式系数den=2 5.5 3 1; %分母多项式系数%绘制原传递函数阶跃响应曲线fprintf(原系统传递函数为 :)G=tf(num,den)L=length(t);cs=randn(1,L); %产生噪声cs=cs/std(cs); cs=cs-mean(cs); 系统辨识 经典辨识方法cs=pj+fc*cs;cs=cs;y=step(num,den,t)+cs; %加入噪声Length=length(y); plot(t,y);grid;xlabel(t/s);ylabel(y(t);%辨识程序设计，此

16、系统 m+n=5,故应计算 A1-A5fprintf(A1-A5 阶面积分别为： )%求 A1sum1=0;for(i=1:Length-1)sum1=sum1+(1-(y(i)+y(i+1)/2)*dt;A(i)=sum1;endA1=sum1%求 A2sum2=0;for(i=1:Length-1)sum2=sum2+(A(i)-A1*(y(i)+y(i+1)/2)*dt;B(i)=sum2;endA2=sum2%求 A3sum3=0;for(i=1:Length-1)sum3=sum3+(B(i)-A2*(y(i)+y(i+1)/2)*dt;C(i)=sum3;endA3=sum3%求

17、A4sum4=0;for(i=1:Length-1)sum4=sum4+(C(i)-A3*(y(i)+y(i+1)/2)*dt;D(i)=sum4;endA4=sum4%求 A5sum5=0;系统辨识 经典辨识方法for(i=1:Length-1)sum5=sum5+(D(i)-A4*(y(i)+y(i+1)/2)*dt;endA5=sum5%求分子系数 b1,b2M=(-1)*(inv(A3,A2;A4,A3)*A4;A5;fprintf(分子多项式系数为： )b1=M(1,1)b2=M(2,1)%求分母系数 a1,a2,a3N=1 0 0;A1 1 0;A2 A1 1*b1;b2;0+A1

18、;A2;A3;fprintf(分母多项式系数为： )a1=N(1,1)a2=N(2,1)a3=N(3,1)%求辨识所得传递函数num1=b2 b1 1;den1=a3 a2 a1 1;fprintf(辨识所得传递函数为 :)G=tf(num1,den1)h=step(num1,den1,t);plot(t,y,black,t,h,blue);legend(原传递函数 ,辨识所得传递函数);title(原传递函数与辨识所得传递函数的阶跃响应对比)grid;xlabel(t/s);ylabel(y(t);fprintf(相关系数：); %求相关系数r=corrcoef(y,h)运行程序，结果如下：

19、原系统传递函数为:G =s2 + 3 s + 1-2 s3 + 5.5 s2 + 3 s + 1Continuous-time transfer function.A1-A5 阶面积分别为：系统辨识 经典辨识方法A1 =0.0081A2 =4.7477A3 =-13.2541A4 =39.7921A5 =-121.2501分子多项式系数为：b1 =3.6415b2 =1.7846分母多项式系数为：a1 =3.6496a2 =6.5619a3 =4.0492辨识所得传递函数为:G =1.785 s2 + 3.641 s + 1-4.049 s3 + 6.562 s2 + 3.65 s + 1Co

20、ntinuous-time transfer function.相关系数：r =系统辨识 经典辨识方法1.0000 0.90760.9076 1.0000此时原传递函数和辨识所得传递函数的阶跃响应对比如下图：0 2 4 6 8 10 12 14 16 18 20-0.4-0.200.20.40.60.811.21.41.6t/sy(t)与 与 与 与 与 与 与 与 与 与 与 与 与 与 与 与 与 与 与 与 与与 与 与 与 与与 与 与 与 与 与 与 与图 1.3 加入噪声后原传递函数和辨识所得传递函数的阶跃响应对比由运行结果和上图可以看出，加入噪声后，系统的辨识结果受到了影响，准确

21、度下降，但基本上可以接受。2 差分方程法2.1 辨识原理设过程传递函数为(2.1)(,1)(10mnsasabbsGnnm则当特征方程有n个单根s 1,s2,sn时，传递函数可写成(2.2)nscsc21)(当特征方程有n-r个单根s 1,s2,sn-r，r阶重根s 0时，传递函数可写成(2.3)rnrnnrn sccsccG )()()( 0200121 系统辨识 经典辨识方法为了确定c i和 si，从获取的过程输出脉冲响应g(t)中，选取前n+1个坐标点，每个坐标点间隔相同的采样时间T 0，组成一个自回归模型(2.4)1()()()ngkgk其中 为待定系数。1， 2如果特征方程(2.5)

22、0xx0021nTT有一个单根，则其必是AR模型的解，它们的线性组合也是 AR模型的解。当其特征方程有n个单根时，自回归模型的解为(2.6)00021)(gkTnkTkTxx当其特征方程有n-r个单根r阶重根时，自回归模型的解为 (2.7)000000 12121)( kTrnkTrnkTrnkTrnkTk xxxxxg 根据分析，无论特征方程是单根还是重根，都有(2.8)iiiixTslog1c0显然，一旦求出x i和 i，便可的传递函数。2.1 程序设计选取系统传递函数如下：()= 2+1.053+5.52+3+1MATLAB 程序设计如下：clc %清空工作区cleardt=0.01;

23、%设置采样时间t=0:dt:50; %设置时间长度num=2 1.05; %分子多项式系数den=1 5.5 3 1; %分母多项式系数%绘制原传递函数阶跃响应曲线fprintf(原系统传递函数为 :)G=tf(num,den)系统辨识 经典辨识方法y=impulse(num,den,t); %脉冲响应曲线Length=length(y); plot(t,y);grid;xlabel(t/s);ylabel(y(t);%时刻为 0s,1s 2s 3s 4s 5s 时的脉冲响应,取 T0=1fprintf(时刻为 0s,1s 2s 3s 4s 5s 时的脉冲响应,取 T0=1);Y=y(1,:)

24、 y(101,:) y(201,:) y(301,:) y(401,:) y(501,:) A=Y(:,2:4);Y(:,3:5);Y(:,4:6); %构造系数矩阵B=Y(:,1:3); %构造方程右边的数值fprintf(系数 ：);a=(inv(A)*(-1)*B) %求得待定系数 ,在这里用 a 表示p=a(3,1) a(2,1) a(1,1) 1; %关于 x 的方程系数x=roots(p); %求得 x1,x2,x3fprintf(系统极点：);s=log(roots(p) %求得传递函数极点%构造系数矩阵，并求得 ，这里用 E 表示C=1 1 1;x(1,1) x(2,1) x(

25、3,1);(x(1,1)2) (x(2,1)2 (x(3,1)2);D=0 Y(:,2) Y(:,3);fprintf(系数 );E=(inv(C)*D % 绘制辨识所得传递函数的脉冲响应num1,den1=residue(E,s,);fprintf(辨识所得传递函数 );G1=tf(num1,den1) %传递函数h=impulse(num1,den1,t); %脉冲响应曲线plot(t,y,black,t,h,blue);legend(原传递函数 ,辨识所得传递函数);title(原系统与辨识所得系统脉冲响应对比);grid;xlabel(t/s);ylabel(y(t)和 h(t);运行

26、以上程序，得到结果如下：原系统传递函数为:G =2 s + 1.05-s3 + 5.5 s2 + 3 s + 1系统辨识 经典辨识方法Continuous-time transfer function.时刻为 0s,1s 2s 3s 4s 5s 时的脉冲响应 ,取 T0=1Y =0 0.3685 0.2956 0.2076 0.1261 0.0607系数 ：a =-141.2815348.0135-244.6919系统极点：s =-0.2835 + 0.3498i-0.2835 - 0.3498i-4.9329 系数 E =0.2028 - 0.1637i0.2028 + 0.1637i-0.

27、4055 辨识所得传递函数 Warning: The numerator or denominator of this transfer function has complex-valued coefficients. In warning at 26In tf.tftf.tf at 357In chafenfangchengfa at 36 G1 =(2-1.11e-16i) s + (1.05-2.22e-16i)-s3 + 5.5 s2 + 3 s + 1Continuous-time transfer function.Warning: The numerator or denomi

28、nator of this transfer function has complex-系统辨识 经典辨识方法valued coefficients. In warning at 26In tf.tftf.tf at 357In impulse at 117In chafenfangchengfa at 37 Warning: Imaginary parts of complex X and/or Y arguments ignored In chafenfangchengfa at 38 此时原传递函数和辨识所得传递函数的阶跃响应对比如下图：0 5 10 15 20 25 30 35 40

29、45 50-0.0500.050.10.150.20.250.30.350.4t/sy(t)与h(t)与 与 与 与 与 与 与 与 与 与 与 与 与 与 与 与与 与 与 与 与与 与 与 与 与 与 与 与图 2.1 原传递函数和辨识所得传递函数的阶跃响应对比由上图和运行结果可以看出，采用差分方程法所得到的辨识结果也是比较精确的。3 总结通过本次实验仿真，我掌握了面积法、差分方程法等经典辨识方法的一些基础编程实验方法，将所学理论知识运用在了实际操作中，对系统辨识的认识和掌握进一步加深。在加入噪声的仿真中，由于噪声时随机序列，所以每次仿真的结果都不相同，并且加入噪声会使辨识结果准确度降低。另外，采样时间选择不同，对辨识结果也有影响。