系统辨识经典辨识方法

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

2、7）利用初始条件（2.1.3）当t时 （1.8）同理有a3=F3()以此类推，若n2，有an=Fn()1.1.2 分子、分母分别为m阶和n阶多项式的系统当传递函数的形式如下所示时 (1.9)定义 (1.10)由于 (1.11)则的Laplace变换为：(1.12)定义一阶面积为：(1.13)令 (1.14)定义二阶面积为：(1.15)同理，令 (1.16)定义阶面积为。由此可得：(1.17)上式可写成如下形式：(1.18)(1.19)通过该系数矩阵，即可求出传递函数分子分母系数的值。1.2 程序设计1.2.1 传递函数形式如式1.1的系统取系统传递函数如下：Gs=13s3+2s2+5.2s+1

3、MATLAB程序如下：clc %清空工作区cleardt=0.01; %设置采样时间t=0:dt:50; %设置时间长度num=1; %此系统分子为1den=3 2 5.2 1; %分母多项式系数%绘制原传递函数阶跃响应曲线fprintf('原系统传递函数为:')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;f

4、or(i=1:Length) sum1=sum1+(1-y(i)*dt; F(i)=sum1;end a1=sum1 %求a2 sum2=0; for(i=1:Length) sum2=sum2+(F(i)-a1*y(i)*dt; f(i)=sum2; enda2=sum2 %求a3 sum3=0; for(i=1:Length) sum3=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,de

5、n2)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=corrcoef(y,h)运行以上程序得到结果如下：原系统传递函数为:G = 1 - 3 s3

6、+ 2 s2 + 5.2 s + 1 Continuous-time transfer function.辨识参数结果：a1 = 5.2048a2 = 2.0608a3 = 2.8388系统辨识后的传递函数为:G = 1 - 2.839 s3 + 2.061 s2 + 5.205 s + 1 Continuous-time transfer function.相关系数：r = 1.0000 0.99990.9999 1.0000此时原传递函数和辨识所得传递函数的阶跃响应对比如下图：图1.1 原传递函数和辨识所得传递函数的阶跃响应对比由上图可以看出，辨识所得结果比较准确。1.2.1 传递函数形式

7、如式(1.9)的系统（无噪声）取系统传递函数如下：Gs=s2+3s+12s3+5.5s2+3s+1MATLAB程序如下：clc %清空工作区cleardt=0.01; %设置采样时间t=0:dt:20; %设置时间长度num=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)'

8、);%辨识程序设计，此系统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

9、)=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=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;A

10、2 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,den1)h=step(num1,den1,t);plot(t,y,'black',t,h,'blue');legend('原传递函数','辨识所得传递函数');title('原传递函数与辨识所得传递函数的阶跃响应

11、对比')grid;xlabel('t/s');ylabel('y(t)');fprintf('相关系数：'); %求相关系数r=corrcoef(y,h)运行以上程序结果如下：原系统传递函数为:G = s2 + 3 s + 1 - 2 s3 + 5.5 s2 + 3 s + 1 Continuous-time transfer function.A1-A5阶面积分别为：A1 = 0.0076A2 = 4.3364A3 = -9.6285A4 = 15.7297A5 = 7.1087分子多项式系数为：b1 = 7.4409b2 = 12.8

12、942分母多项式系数为：a1 = 7.4485a2 = 17.2869a3 = 22.7359辨识所得传递函数为:G = 12.89 s2 + 7.441 s + 1 - 22.74 s3 + 17.29 s2 + 7.448 s + 1 Continuous-time transfer function.相关系数：r = 1.0000 0.9996 0.9996 1.0000此时原传递函数和辨识所得传递函数的阶跃响应对比如下图：图1.2 原传递函数和辨识所得传递函数的阶跃响应对比由上图可以看出，在未加入噪声之前，采用面积法辨识结果很精确，并且，分子可以为阶次低于分母的任意阶次。1.2.2 加

13、入噪声为便于对比，仍取系统传递函数如下：Gs=s2+3s+12s3+5.5s2+3s+1MATLAB程序如下：clc %清空工作区cleardt=0.01; %设置采样时间t=0:dt:20; %设置时间长度pj=0; fc=sqrt(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);

14、 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)');%辨识程序设计，此系统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

15、=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=sum5+(D(i)-A4*(y(i)+y(i+1)/2)*dt;en

16、dA5=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,den1)h=step(num1,

17、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 + 1 Continu

18、ous-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 + 1 Continuous-time transfer function.相关系数：r = 1.0000 0.907

19、6 0.9076 1.0000>>此时原传递函数和辨识所得传递函数的阶跃响应对比如下图：图1.3 加入噪声后原传递函数和辨识所得传递函数的阶跃响应对比由运行结果和上图可以看出，加入噪声后，系统的辨识结果受到了影响，准确度下降，但基本上可以接受。2 差分方程法2.1 辨识原理设过程传递函数为 (2.1)则当特征方程有n个单根s1,s2,···sn时，传递函数可写成 (2.2)当特征方程有n-r个单根s1,s2,···sn-r，r阶重根s0时，传递函数可写成 (2.3)为了确定ci和si，从获取的过程输出脉冲响应g(t)中，选

20、取前n+1个坐标点，每个坐标点间隔相同的采样时间T0，组成一个自回归模型 (2.4)其中1，2n为待定系数。如果特征方程 (2.5)有一个单根，则其必是AR模型的解，它们的线性组合也是AR模型的解。当其特征方程有n个单根时，自回归模型的解为 (2.6)当其特征方程有n-r个单根r阶重根时，自回归模型的解为 (2.7)根据分析，无论特征方程是单根还是重根，都有 (2.8)显然，一旦求出xi和i，便可的传递函数。2.1 程序设计选取系统传递函数如下：Gs=2s+1.05s3+5.5s2+3s+1MATLAB程序设计如下：clc %清空工作区cleardt=0.01; %设置采样时间t=0:dt:5

21、0; %设置时间长度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');

22、Y=y(1,:) 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=

23、1 1 1;x(1,1) x(2,1) x(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('原传递函数'

24、;,'辨识所得传递函数');title('原系统与辨识所得系统脉冲响应对比');grid;xlabel('t/s');ylabel('y(t)和h(t)');运行以上程序，得到结果如下：原系统传递函数为: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.2815 348.0135 -244.6919系统极点：s = -0.2835 + 0.3498i -0.2835 - 0.3498i -4.9329 系数E = 0.2028 - 0.1637i 0.2028 + 0.1637i -0.4055