系统辨识算法

传统系统辨识算法1.引言迄今为止，已经有许多不同的辨识方法。这些辨识方法就其所涉及的模型的形式来说可以分为两类。一类是非参数模型的辨识方法，一类是参数模型的辨识方法。非参数模型的辨识方法（亦称经典的辨识放法）获得的模型是非参数模型。它在假定过程是线性的前提下，不必事先确定模型的具体结构，因而这类方法可适用于任意复杂的过程，工程上至今仍经常采用。参数模型的辨识方法（亦称现代的辨识方法）必须假定一种模型结构，通过极小化模型与过程之间的误差准则函数来确定模型的参数。如果模型的结构无法实现确定，则必须利用结构辨识方法先确定模型的结构参数（比如阶次、纯延迟等），再进一步确定模型参数。参数模型的辨识方法又可以分为：最小二乘法辨识、梯度校正法辨识以及极大似然法辨识。根据计算机与过程之间的不同联接方式，辨识又可以分为离线辨识和在线辨识。离线辨识首先将采集到的数据储存在磁盘或磁带中，然后将数据成批输入计算机进行辨识计算。这种辨识方式多采用成批处理的算法，或称一次完成算法，其缺点是占用内存较大。在线辨识通常要在正常运行工况下进行，它一般采用实际处理算法，即每采样一组数据就进行一次辨识计算。这种辨识方式占用内存量比较小；尤其对时变过程的辨识或自适应控制问题来说，它比离线辨识方式具有更多的优势。本次作业使用经典的辨识中的一些方法对系统进行辨识。在经典的控制理论中，线性过程的动态特性通常用：传递函数G(s)频率响应G(jw)脉冲响应g(t)阶跃响应h(t)来表示。后三种为非参数模型，其表现形式是以时间或频率为自变量的实验曲线。对过程施加特定的实验信号，同时测定过程的输出，可以求得这些非参数模型。经过适当的数学处理，它们又可以转变成参数模型——传递函数的形式。获取上述非参数模型，并把它们转化成传递函数的主要方法有：阶跃响应法脉冲响应法频率响应法相关分析法谱分析法这些辨识方法在工程上有广泛的应用，至今仍受到普遍重视。在本次作业中，我主要使用了阶跃响应法中的面积法和脉冲响应法来对系统进行辨识，并且对系统施加一定的噪声干扰，比较在有误噪声情况下辨识结果的不同。下面将对其分别进行介绍。2.面积法2.1面积法的基本原理面积法为阶跃响应法之中的一种，阶跃响应法还包括两点法、切线法、近似法等。在系统的阶跃响应曲线呈现不规则的形状时，使用面积法比较方便。在使用面积法对系统进行辨识，首先就是要测得过程的阶跃响应，测试路线如图2-1所示。图2-1过程阶跃响应的测取路线经过试验获取的阶跃响应要将其转化成无因此的形式，如图2-2所示：图2-2把阶跃响应转化成无因次的形式设过程的传递函数为：1111111()()1mmmmnnnnbsbsbsGsKnmasasas（2.1.1）其中阶次m，n需要根据先验知识和精度考虑实现选定。显然有0()hKu（2.1.2）定义：1()()GsKPs（2.1.3）其中：1111111()1nnnnmmmmasasasPsbsbsbs11iiics（2.1.4）那么[1*()]ht的Laplace变换为：11111[1*()]()1iiiiiicsLhtssPscs（2.1.5）则一阶面积A1（图2-2中阴影部分）为：10011101[1*()][1*()]1limlimsiiiisiiAhtdtLhtcsccs（2.1.6）再令111[*()](1)Lhtscs（2.1.7）定义二阶面积为：211000010222011[*()*()]{[*()*()]}[*()][*()](1)(1)limlimlimttssiiiisiiAhhddtLhthtdtLhtLhtscsccscs（2.1.8）同理，令22121[*()](1)Lhtscscs（2.1.9）可得三阶面积为：2323000[*()*()]tAhhddtc（2.1.10）以此类推i阶面积Ai为：11000[*()*()]tiiiiAhhddtc（2.1.11）其中：1211211[*()](1)iiiLhtscscscs（2.1.12）进一步利用：221()()()1!2!!istisssettti（2.1.13）可得：00[1*()][1*()]stiiiLhthtedtMs（2.1.14）其中：0()[1*()]!iitMhtdti（2.1.15）由式（2.1.5）、（2.1.14）以及iicA可得：111111iiiiiiAsMs（2.1.16）显然，式（2.1.16）左边s各次幂项的系数均为零，固有：1021103211202110iiiijjjAMAMAMAMAMAMAMAM（2.1.17）即：121000()()[1*()][1*()](1)!!ijtiijjttAhtdtAhtdtij（2.1.18）又根据式（2.1.14）以及iicA可得：11111111(1)(1)nnnnmmimmiiasasasbsbsbsAs（2.1.19）比较上式两边s各次幂的系数，有：111222113332112111,2,,iiiijijjaAbaAbbAaAbbAbAaAbbAinm（2.1.20）有上式便可以求出传递函数的系数1212,,,;,,,nmaaabbb。应当指出，为了求（n+m)个未知数，需要（n+m）个方程，因此需要利用式（2.1.18）求1,2,......,inm阶的面积。注意到0()iain以及0()ibim，则式（2.1.20）可以写成矩阵形式：111112122121112212121100010010nnnmnnnnmnmnmnmnnmmnnnnbAAAAbAAAAbAAAAbaAbaAAbaAAAA（2.1.21）由此，便可以利用面积法对系统进行辨识。但是，我们可以看出，利用面积法进行辨识计算量是比较大的，但是，计算机有着强大的计算功能，所以，将面积法（主要是式（2.1.18）以及（2.1.21））编制成计算机程序之后，再进行仿真计算，就会变得非常简单。在下一小节我们将进行详细的仿真介绍。2.2MATLAB仿真对于面积法，应用MATLAB进行仿真时最主要的几个公式在上一节已经提到，就是公式（2.1.8）以及（2.1.21），公式（2.1.8）使用来求取对应的面积，公式（2.1.21）使用来求取函数分子分母多项式的系数。下面我们对课本中的例题4.1进行仿真，运行程序如下所示：clearallclcdt=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;fori=1:MA(1)=A(1)+(1-yn(i))*dt;endA(1);%计算面积A2A21=0;fori=1:MA21=A21+(1-yn(i))*(-1)*t(i)*dt;endA22=0;fori=1:MA22=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)plot(t,yz)xlabel('TimeSecond');ylabel('StepResponse');title('加入噪声之后系统的单位阶跃响应');Az=[0,0];%计算面积Az1；Az(1)=0;fori=1:MAz(1)=Az(1)+(1-yz(i))*dt;endAz(1);%计算面积Az2Az21=0;fori=1:MAz21=Az21+(1-yz(i))*(-1)*t(i)*dt;endAz22=0;fori=1:MAz22=Az22+(1-yz(i))*dt;endAz22=Az22*Az(1);Az(2)=Az21+Az22;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('TimeSecond');ylabel('StepResponse');title('有无噪声辨识的结果');legend('原系统阶跃响应','无噪声情况下的辨识','加入噪声之后的辨识')仿真结果：原系统的传递函数：sys=1------------------10s^2+6.5s+1没有噪声情况下辨识系统的传递函数：sysc=1---------------------10.12s^2+6.5s+1有噪声情况下辨识系统的传递函数：syszc=1-----------------------19.17s^2+5.659s+1系统的单位阶跃响应曲线如图2-2-1所示：图2-2-1各系统阶跃响应曲线仿真分析：我们可以看到，在没有噪声的情况下辨识系统和原系统的响应曲线完全重合，有噪声的情况下辨识的结果就会出先一定的偏差，反复检查修改程序之后结果还是没有变化，而且这个程序存在一个很大的缺点，那就是只能辨识分母为二阶，分子为零阶的系统，因为面积A1，A2是分别编写程序求解的，求解面积的程序并不是一个通用的程序，想要求解A3，A4或者更多的面积就得重新再加入求解A3，A4的程序，求解系数矩阵a,b时所用到的矩阵也是根据公式人为直接输入的，并051015202530354000.20.40.60.811.21.4TimeSecondStepResponse有无噪声辨识的结果原系统阶跃响应无噪声情况下的辨识加入噪声之后的辨识不是通过程序对n,m的判断来自动构建的。所以，鉴于以上的种种缺点，对程序进行了改进，改进后的程序如下所示：clearclcAs=0;As1=0;Azs=0;Azs1=0;dt=0.005;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=[];%无噪声情况下的面积%无噪声%对求解面积进行了改进fori=1:n+mAf=0;forj1=1:MAf=Af+(1-yn(j1))*((-t(j1))^(i-1))*dt/(factori