系统辨识算法

1、传统系统辨识算法1. 引言迄今为止，已经有许多不同的辨识方法。这些辨识方法就其所涉及的模型的形式来说可以分为两类。一类是非参数模型的辨识方法，一类是参数模型的辨识方法。非参数模型的辨识方法（亦称经典的辨识放法）获得的模型是非参数模型。它在假定过程是线性的前提下，不必事先确定模型的具体结构，因而这类方法可适用于任意复杂的过程，工程上至今仍经常采用。参数模型的辨识方法（亦称现代的辨识方法）必须假定一种模型结构，通过极小化模型与过程之间的误差准则函数来确定模型的参数。如果模型的结构无法实现确定，则必须利用结构辨识方法先确定模型的结构参数（比如阶次、纯延迟等），再进一步确定模型参数。参数模型的辨识方法

2、又可以分为：最小二乘法辨识、梯度校正法辨识以及极大似然法辨识。根据计算机与过程之间的不同联接方式，辨识又可以分为离线辨识和在线辨识。离线辨识首先将采集到的数据储存在磁盘或磁带中，然后将数据成批输入计算机进行辨识计算。这种辨识方式多采用成批处理的算法，或称一次完成算法，其缺点是占用内存较大。在线辨识通常要在正常运行工况下进行，它一般采用实际处理算法，即每采样一组数据就进行一次辨识计算。这种辨识方式占用内存量比较小；尤其对时变过程的辨识或自适应控制问题来说，它比离线辨识方式具有更多的优势。本次作业使用经典的辨识中的一些方法对系统进行辨识。在经典的控制理论中，线性过程的动态特性通常用： 传递函数 G

3、(s) 频率响应 G(jw) 脉冲响应 g(t) 阶跃响应 h(t)来表示。后三种为非参数模型，其表现形式是以时间或频率为自变量的实验曲线。对过程施加特定的实验信号，同时测定过程的输出，可以求得这些非参数模型。经过适当的数学处理，它们又可以转变成参数模型传递函数的形式。获取上述非参数模型，并把它们转化成传递函数的主要方法有： 阶跃响应法 脉冲响应法 频率响应法 相关分析法 谱分析法这些辨识方法在工程上有广泛的应用，至今仍受到普遍重视。在本次作业中，我主要使用了阶跃响应法中的面积法和脉冲响应法来对系统进行辨识，并且对系统施加一定的噪声干扰，比较在有误噪声情况下辨识结果的不同。下面将对其分别进行介

4、绍。2. 面积法2.1 面积法的基本原理面积法为阶跃响应法之中的一种，阶跃响应法还包括两点法、切线法、近似法等。在系统的阶跃响应曲线呈现不规则的形状时，使用面积法比较方便。在使用面积法对系统进行辨识，首先就是要测得过程的阶跃响应，测试路线如图2-1所示。图2-1 过程阶跃响应的测取路线经过试验获取的阶跃响应要将其转化成无因此的形式，如图2-2所示：图2-2 把阶跃响应转化成无因次的形式设过程的传递函数为： （2.1.1）其中阶次m，n需要根据先验知识和精度考虑实现选定。显然有 （2.1.2）定义： （2.1.3）其中： （2.1.4）那么的Laplace变换为： （2.1.5）则一阶面积A1（

5、图2-2中阴影部分）为： （2.1.6）再令 （2.1.7）定义二阶面积为： （2.1.8）同理，令 （2.1.9）可得三阶面积为： （2.1.10）以此类推i阶面积Ai为： （2.1.11）其中： （2.1.12）进一步利用： （2.1.13）可得： （2.1.14）其中： （2.1.15）由式（2.1.5）、（2.1.14）以及可得： （2.1.16）显然，式（2.1.16）左边s各次幂项的系数均为零，固有： （2.1.17）即： （2.1.18）又根据式（2.1.14）以及可得： （2.1.19）比较上式两边s各次幂的系数，有： （2.1.20）有上式便可以求出传递函数的系数。应当指出，

6、为了求（n+m)个未知数，需要（n+m）个方程，因此需要利用式（2.1.18）求阶的面积。注意到以及，则式（2.1.20）可以写成矩阵形式： （2.1.21）由此，便可以利用面积法对系统进行辨识。但是，我们可以看出，利用面积法进行辨识计算量是比较大的，但是，计算机有着强大的计算功能，所以，将面积法（主要是式（2.1.18）以及（2.1.21）编制成计算机程序之后，再进行仿真计算，就会变得非常简单。在下一小节我们将进行详细的仿真介绍。2.2 MATLAB仿真对于面积法，应用MATLAB进行仿真时最主要的几个公式在上一节已经提到，就是公式（2.1.8）以及（2.1.21），公式（2.1.8）使用来

7、求取对应的面积，公式（2.1.21）使用来求取函数分子分母多项式的系数。下面我们对课本中的例题4.1进行仿真，运行程序如下所示：clear allclcdt=0.01;dt1=0.01;tmax=40;t=0:dt:tmax;M=1+tmax/dt;num=1;den=10,6.5,1;sys=tf(num,den)y=step(sys,t);yn=y/y(M);%把阶跃响应转化成无因次型A=0,0;%计算面积A1；A(1)=0;for i=1:M A(1)=A(1)+(1-yn(i)*dt;endA(1);%计算面积A2A21=0;for i=1:M A21=A21+(1-yn(i)*(-1

8、)*t(i)*dt;endA22=0;for i=1:M A22=A22+(1-yn(i)*dt;endA22=A22*A(1);A(2)=A21+A22;A%bnum1=1;%求aC2=1,0; A(1),1;D2=0;0;den1=;AA=A(1);A(2);den11=(C2*D2+AA)'den1=den11(2),den11(1),1;sysc=tf(num1,den1)yc=step(sysc,t);%没有噪声时的辨识系统yc=yc/yc(M); %加入噪声之后再进行辨识yz=yn+0.5*yn.*rand(size(yn);%yz=yz/yz(M);figure(2)pl

9、ot(t,yz)xlabel('Time Second');ylabel('Step Response');title('加入噪声之后系统的单位阶跃响应');Az=0,0;%计算面积Az1；Az(1)=0;for i=1:M Az(1)=Az(1)+(1-yz(i)*dt;endAz(1);%计算面积Az2Az21=0;for i=1:M Az21=Az21+(1-yz(i)*(-1)*t(i)*dt;endAz22=0;for i=1:M Az22=Az22+(1-yz(i)*dt;endAz22=Az22*Az(1);Az(2)=Az21+A

10、z22;Az%bznumz1=1;%azCz2=1,0; Az(1),1;Dz2=0;0;denz1=;AAz=Az(1);Az(2);denz11=(Cz2*Dz2+AAz)'denz1=denz11(2),denz11(1),1;syszc=tf(numz1,denz1)yzc=step(syszc,t);%加入噪声之后的辨识系统yzc=yzc/yzc(M);figure(3)plot(t,y,'r-',t,yc,'b-.',t,yzc,'k')xlabel('Time Second');ylabel('Ste

11、p Response');title('有无噪声辨识的结果');legend('原系统阶跃响应','无噪声情况下的辨识','加入噪声之后的辨识')仿真结果：原系统的传递函数：sys = 1 - 10 s2 + 6.5 s + 1 没有噪声情况下辨识系统的传递函数：sysc = 1 - 10.12 s2 + 6.5 s + 1有噪声情况下辨识系统的传递函数：syszc = 1 - 19.17 s2 + 5.659 s + 1 系统的单位阶跃响应曲线如图2-2-1所示：图2-2-1各系统阶跃响应曲线仿真分析：我们可以看到，在没

12、有噪声的情况下辨识系统和原系统的响应曲线完全重合，有噪声的情况下辨识的结果就会出先一定的偏差，反复检查修改程序之后结果还是没有变化，而且这个程序存在一个很大的缺点，那就是只能辨识分母为二阶，分子为零阶的系统，因为面积A1，A2是分别编写程序求解的，求解面积的程序并不是一个通用的程序，想要求解A3，A4或者更多的面积就得重新再加入求解A3，A4的程序，求解系数矩阵a,b时所用到的矩阵也是根据公式人为直接输入的，并不是通过程序对n,m的判断来自动构建的。所以，鉴于以上的种种缺点，对程序进行了改进，改进后的程序如下所示：clearclcAs=0;As1=0;Azs=0;Azs1=0;dt=0.005

13、;tmax=40;t=0:dt:tmax;M=1+tmax/dt;num=1;den=10,6.5,1;n=2;%可以更改m=0;%可以更改sys=tf(num,den)%原系统的传递函数y=step(sys,t);yn=y/y(M);%把阶跃响应转化成无因次型A=;%无噪声情况下的面积%无噪声%对求解面积进行了改进for i=1:n+m Af=0; for j1=1:M Af=Af+(1-yn(j1)*(-t(j1)(i-1)*dt/(factorial(i-1); end for j2=0:i-2 As1=0; As=0; for j3=1:M As1=As1+(1-yn(j3)*(-t(

14、j3)j2)*dt/(factorial(j2); end As=A(i-j2-1)*As1+As; end A(i)=As+Af;endA;%求解b%对求解a,b时用到的矩阵的构建方法进行了改进Ab=zeros(m,m);for i=1:m; Ab(i,:)=A(n+i-1):-1:(n-m+i);endAb;Anb=A(n+1:n+m);b=-inv(Ab)*Anb'b=b'num1=b(m:-1:1),1;%求解aAa=ones(n,n);for i=1:n for j=1:n if(j>i) Aa(i,j)=0; end endendfor i=1:n for j

15、=1:i-1 Aa(i,j)=A(i-j); endendAa; for i=1:n if(i<=m) Aaa(i,1)=b(i); else Aaa(i,1)=0; end endAaa;Aaaa=A(1:n);a=Aa*Aaa+Aaaa'a=a'den1=a(n:-1:1),1;sysc=tf(num1,den1)%没有噪声时辨识系统的传递函数yc=step(sysc,t);yc=yc/yc(M);%加入噪声之后再进行辨识yz=yn+0.5*yn.*rand(size(yn);figure(1)plot(t,yz)xlabel('Time Second'

16、;);ylabel('Step Response');title('加入噪声之后系统的单位阶跃响应');Az=;%有噪声情况下的面积for i=1:n+m Azf=0; for j1=1:M Azf=Azf+(1-yz(j1)*(-t(j1)(i-1)*dt/(factorial(i-1); end for j2=0:i-2 Azs1=0; Azs=0; for j3=1:M Azs1=As1+(1-yz(j3)*(-t(j3)j2)*dt/(factorial(j2); end Azs=Az(i-j2-1)*Azs1+Azs; end Az(i)=Azs+Az

17、f;endAz;%求解bzAzb=zeros(m,m);for i=1:m; Azb(i,:)=Az(n+i-1):-1:(n-m+i);endAzb;Aznb=Az(n+1:n+m);bz=-inv(Azb)*Aznb'bz=bz'numz1=bz(m:-1:1),1;%求解azAza=ones(n,n);for i=1:n for j=1:n if(j>i) Aza(i,j)=0; end endendfor i=1:n for j=1:i-1 Aza(i,j)=Az(i-j); endendAza; for i=1:n if(i<=m) Azaa(i,1)=b

18、z(i); else Azaa(i,1)=0; end endAzaa;Azaaa=Az(1:n);az=Aza*Azaa+Azaaa'az=az'denz1=az(n:-1:1),1;syszc=tf(numz1,denz1)%加入噪声之后辨识系统的传递函数yzc=step(syszc,t);yzc=yzc/yzc(M);figure(2)plot(t,y,'r-',t,yc,'b-.',t,yzc,'k')xlabel('Time Second');ylabel('Step Response')

19、;title('有无噪声辨识的结果');legend('原系统阶跃响应','无噪声情况下的辨识','加入噪声之后的辨识')改进之后程序的仿真结果为：原系统的传递函数：sys = 1 - 10 s2 + 6.5 s + 1 无噪声时辨识系统的传递函数：sysc = 1 - 10.09 s2 + 6.498 s + 1 有噪声时辨识系统的传递函数：syszc = 1 - 23.89 s2 + 5.658 s + 1原系统以及各辨识系统单位阶跃响应曲线如图2-2-2所示：图2-2-2各系统单位阶跃响应曲线仿真分析：程序改进之后对原来的系

20、统进行仿真辨识，发现结果和原系统的几乎没有区别，说明辨识结果的可信度还是有的。因为改进的程序可以选择传递函数的阶次，所以我们改变原来的系统再次进行仿真，将系统的传递函数改成 此时，令n=4,m=2,得出仿真结果：原系统的传递函数：sys = 2 s2 + 3 s + 1 - 2 s4 + 12 s3 + 10 s2 + 6.5 s + 1 无噪声时辨识系统的传递函数：sysc = 10.37 s2 + 6.049 s + 1 - 3.857 s4 + 41.07 s3 + 29.17 s2 + 9.544 s + 1有噪声时辨识系统的传递函数：syszc = -6.301 s2 + 6.308

21、 s + 1 - 4924 s4 + 2.788e04 s3 + 3854 s2 + 8.893 s + 1 图2-2-3各系统单位阶跃响应曲线我们可以看到在没有噪声的情况下辨识结果还是不错的，但是在噪声的影响下辨识效果就变得很差，说明面积法不能很好的消除噪声的影响。并且，n,m的选择也很关键，上述仿真是在令n=4,m=2的情况下进行的仿真，这是最理想的结果，但是事实往往是我们要对m,n的值进行推测，如果我们令n=4,m=1，效果便非常的不理想，辨识之后系统的阶跃响应曲线如图2-2-4所示，我们可以看出，辨识系统和原系统之间差距非常大，所以在利用面积法求传递函数时必须正确选择传递函数的阶次。图

22、2-2-4各系统单位阶跃响应曲线其次在仿真的过程之中还发现了一个问题，当原系统系统分子分母之间的阶次之差大于2时，即使m,n选择正确，也无法正确的对系统进行仿真，例如设系统传递函数为：，我们令n=3,m=0,仿真结果如图2-2-5所示，造成这结果的原因现在还不清楚。图2-2-5各系统单位阶跃响应曲线3. 脉冲响应法3.1 脉冲响应法的基本原理脉冲响应是在理想脉冲输入作用下过程的输出响应。考虑到工程上实际输入理想脉冲是不可能的，因此通常采用矩形脉冲输入，如图3-1-11所示：图3-1-1 过程的脉冲响应当矩形脉冲的宽度比过程的过渡时间小得多、且矩形脉冲的面积等于1时，过程的输出可以近似为脉冲输出

23、。由脉冲响应确定传递函数的方法有很多，下面介绍几种常用的方法。1.一阶过程如果过程能用一阶传递函数描述 ，则传递函数的参数T和K可以直接在脉冲响应曲线上确定，如图3-1-2所示： （3.1.1）图3-1-2一阶惯性环节脉冲响应与传递函数参数的关系2. 二阶过程如果过程能用二阶传递函数描述 ，则传递函数的参数也可以直接由脉冲响应曲线确定，如图3-1-3所示： （3.1.2）图3-1-3二阶惯性环节脉冲响应与传递函数参数的关系3. 差分方程法设过程的传递函数为： （3.1.3）当特征方程具有n个单根s1,s2,sn时，传递函数可以写成 （3.1.4）对应的脉冲响应为： （3.1.5） 当特征方程有

24、n-r个单根s1,s2,sn-r，r 阶重根s0时，传递函数可写成 （3.1.6）对应的脉冲响应为：（3.1.7）为了确定ci和si，从获取的过程输出脉冲响应g(t)中，选取前n+1个坐标点，每个坐标点间隔相同的采样时间T0，如图所示。各坐标点上的脉冲响应记为g(k)，g(k+1)，g(k+n)，组成一个自回归模型（AR）： （3.1.8）其中为待定系数。如果特征方程 （3.1.9）有一个单根 ，则必是AR模型的解，它们的线性组合也是AR模型的解.当其特正方程有n个单根时，自回归模型的解为： （3.1.10）当其特征方程有n-r个单根，r阶重根时，自回归模型的解为： (3.1.11)对照上述公

25、式，不论单根还是重根，都有： （3.1.12）3.2 MATLAB仿真对于脉冲响应法一阶和二阶过程比较简单，在这里就针对第三种差分方程法进行仿真。在第一次进行仿真时和面积法的情况一样，编写的程序只能对固定阶数的系统进行仿真，在这实际情况中是不允许的，所以经过改进的程序可以对任意阶的系统进行仿真，但是，后来发现程序有一个问题，就是当特征方程具有重根时，虽然能够求得和，但是在求解传递函数时还没有编写出一个通用的程序，只能是根据已得到的参数，然后根据式（3.1.6）进行认为的求解传递函数，这是本程序的一个不足之处。本次的MATLAB仿真是根据式（3.1.4）、（3.1.8）、（3.1.9）、（3.1

26、.10）、（3.1.12）进行编写的，程序如下所示：clear allclcnum=1.92;den=1,3.8,3.04,1.92;n=3;%可以手动修改阶次进行辨识直到最佳。disp('原系统传递函数')sys=tf(num,den)t=0:0.01:15;y=impulse(sys,t);%脉冲响应；%在T0=1的情况下提取一组数据,最少需要提取2n个数据a=y(1:100:(2*n-1)*100+1);%取k=0,1，2，.，n-1；A=;for i=1:n A(i,:)=a(i+1:1:i+n);endA;for i=1:n B(i,:)=-a(i);endB;X=A

27、B;%求解方程组的解X=X'den1=;den1=X(n:-1:1),1;P=roots(den1);%求解对应的特特征方程的解%求解系统的极点；for i=1:n s(i)=log(P(i);ends;%求取留数for i=1:n A1(i,:)=(P(1:n).(i-1);endA1;for i=1:n B1(i,:)=a(i);endB1;Y=;Y=A1B1;%求取系数之后得到传递函数den1=;num1=;for i=1:n zj=1,-s(i); den1(i,:)=zj; num1(i)=Y(i);endsyscom=tf(num1(1),den1(1,:);for i=2

28、:n syscom=parallel(syscom,tf(num1(i),den1(i,:);enddisp('无噪声时辨识系统的传递函数')syscomy1=impulse(syscom,t); %加入噪声之后再辨识zs=y+y.*rand(size(y)*0.5;%加入噪声figure(1)plot(t,zs)xlabel('Time Second')ylabel('Impulse Response')title('加入噪声后系统的脉冲响应图')%在T0=1的情况下提取一组数据az=zs(1:100:(2*n-1)*100+1);Az=;for i=1:n Az(i,:)=az(i+1:1:i+n);endfor i=1:n Bz(i,:)=-az(i);endBz;Xz=AzBz;%求解方程组的解Xz=Xz'den1z=Xz(n:-1:1),1;Pz=roots(den1z);%求解对应的特特征方程的解%求解系统的极点；for i=1:n sz(i)=log(Pz(i);end%求取留数for i=1:n Az(