系统辨识之经典辨识法

系统辨识作业一学院信息科学与工程学院专业控制科学与工程班级控制二班姓名学号2018 年 11 月系统辨识所谓辨识就是通过测取研究对象在认为输入作用的输出响应，或正常运行时的输入输出数据记录，加以必要的数据处理和数学计算，估计出对象的数学模型。

辨识的内容主要包括四个方面：①实验设计；②模型结构辨识；③模型参数辨识；④模型检验。

辨识的一般步骤：根据辨识目的，利用先验知识，初步确定模型结构；采集数据；然后进行模型参数和结构辨识；最终验证获得的最终模型。

根据辨识方法所涉及的模型形式来说，辨识方法可以分为两类：一类是非参数模型辨识方法，另一类是参数模型辨识方法。

其中，非参数模型辨识方法又称为经典的辨识方法，它主要获得的是模型是非参数模型。

在假定过程是线性的前提下，不必事先确定模型的具体结构，广泛适用于一些复杂的过程。

经典辨识方法有很多，其中包括阶跃响应法、脉冲响应法、相关分析法和普分析法等等，本次实验所采用的辨识方法为阶跃响应法和脉冲响应法。

1．阶跃响应法阶跃响应法是一种常用非参数模型辨识方法。

常用的方法有近似法、半对数法、切线法、两点法和面积法等。

本次作业采用面积法求传递函数。

1.1面积法① 当系统的传递函数无零点时，即系统传递函数如下：G(S) = + −11−1+⋯+ 1+1(1-1) 系统的传递函数与微分方程存在着一一对应的关系，因此，可以通过求取微分方程的系数来辨识系统的传递函数。

在求得系统的放大倍数K后，要得到无因次阶跃响应y(t)(设τ=0)，其中y(t)用下式描述：() −1()(1-2) 面积法原则上可以求出n为任意阶的个系数。

以n为3为例。

有：3() 2() (){| →∞ =| →∞ =| →∞ = 0 (1-3)()| →∞ = 1将式（1）中的y(t)移至右边，在[0，t]上积分，得2()3(1-4) 定义：1( ) = ∫0[1 − ()] (1-5) 由式（1-3）条件可知，当t→∞时，(1-6)同理，定义2(1-7)由式（1-,3）条件可知，当t→∞时，(1-8)因此，可得( ) = ∫0[ −1( ) − −1()] dt (1-9)= (∞) (1-10)② 当系统的传递函数存在零点时，传递函数如下：=kG（s）b s mmn +ba s mn-1-1s mn-1-1 ++LL ++a sbs11 +1+1，（n m）（1-11）a s n +其中，K h= ( )/ U0定义1G(s)=KP(s)其中，P(s) = n mn ++ba s mn-1-1s mn-1-1++LL ++a sbs11 +1+1 = +1 i 1 C s i i（1-12）m根据[1−h*(t)]的Laplace变换，求出一阶面积A1，确定L[h（*1 t ]），并定义二阶面积A2 ，以此类推，得到i 阶面积A i 。

进一步利用e−st 拉氏变换，得到L[1−h*（t ]）=M s i i ，进而得到A i 的值：i=0A= 01−h*(t)(i 1)!−−t)i−1dt +tj−=20 A i−−j101−h*(t)−j!t) j dt（1-13） (根据A C i = i ，可得：+ −1 −1 + ⋯ + 1 + 1= (+ −1 −1 + ⋯ + 1+ 1)(1 + ∑∞=1 )。

比较上式两边s的各次幂，便可得到a, b, A之间的关系，如下：b1 A n A n−1 L A n m− +1−1 A n+1b2 A n+1 A n L A n m− +2A n+2 =−M L L L L Mb m A n m+ −1 A n m+ −2 L A n A n m+b1a1 110 LL 0 00 0b M2 AA12（1-14）a2 = A1ML L L L L Mb m A na n A n−1 A n−2 L A1 10由此可知，根据式(1-12)、(1-13)、(1-14)便可得到辨识传递函数的参数a, b。

1.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=1:L s1=s1+(1-h(i))*T; F(i)=s1; enda1=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)disp('通过辨识得到的传递函数为：')G2=tf(num1,den1)%=============原传递函数和辨识函数的阶跃响应对比图=============% step(G1,'b-',G2,'r-.')title('原系统与辨识后所得到系统阶跃响应对比') legend('原响应曲线','辨识响应曲线') （1）当采样周期T=0.2秒，采样时间t=30s时，行程序后得到原传递函数G1和辨识得到的传递函数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;w=qw+sqrt(fc)*w;%=====================阶跃采样序列中加入白噪声==================%h=h+w; plot(t,w);（1）加入的噪声如下图所示：图1.7（2）当采样周期T=0.02s,采样时间t=50s时，辨识结果如下：图1.8原系统与辨识系统阶跃响应如图所示：结合上述无测量噪声和有测量噪声两种情况下的辨识结果，列出如下所示的表格：表1-1结果可知，在相同的采样周期下，适当的增加采样时间，可以提高辨识精度，尤其是对增益的提高有很大影响；而在相同的采样时间下，适当的减小采样时间，对于系统参数的辨识精度有很大的提高。

因此，可以发现合理采样时间和数据长度，可以提高辨识的精度，令辨识后的传递函数系数与原传递函数系数更接近，差距小，从而得到满意的辨识结果。

通过对比无测量噪声和有测量噪声两种情况下的辨识结果，我们可以发现在白噪声的情况下，曲线拟合较无噪声情况下要差，说明白噪声对于面积法辨识系统存在较大的干扰，会对辨识结果产生一定的影响。

1.2.2有零点模型系统17.52+7.5+1 假设系统的传递函数为G(s) =43+52+8+1,为有零点的模型，其中n=3,m=2, 用面积法需要求解~5，利用 Matlab 编程，分别在没有噪声和有噪声两种情况下进行辨识，比1较辨识结果。

（1）没有噪声时，程序如下：clear;%==================获得原传递函数方程=======================%num=[17.5 7.5 1]; den=[4 5 8 1];%=====================产生阶跃采样序列======================%T=0.02; %采样周期 t=0:T:100; %采样时间 L=length(t); %数据长度 y=step(num,den,t);k=y(L) %系统增益%======================面积法求解参数======================%sum1=0; for i=1:L-1;sum1=sum1+(1-(y(i)+y(i+1))/2)*T;A(i)=sum1; end A1=sum1 sum2=0; fori=1:L-1;sum2=sum2+(A(i)-A1*(y(i)+y(i+1))/2)*T;B(i)=sum2; end A2=sum2 sum3=0; for i=1:L-1;sum3=sum3+(B(i)-A2*(y(i)+y(i+1))/2)*T;C(i)=sum3; end A3=sum3 sum4=0;for i=1:L-1;sum4=sum4+(C(i)-A3*(y(i)+y(i+1))/2)*T;D(i)=sum4; end A4=sum4 sum5=0; fori=1:L-1; sum5=sum5+(D(i)-A4*(y(i)+y(i+1))/2)*T; endA5=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];disp('原传递函数为:')G1=tf(num,den)disp('通过辨识得到的传递函数为：')G2=tf(num1,den1)%=============原传递函数和辨识函数的阶跃响应对比图=============%step(G1,'b-',G2,'r-.')title('原系统与辨识后所得到系统阶跃响应对比') legend('原响应曲线','辨识响应曲线') 当采样时间取0.02秒，数据长度为100时，辨识结果如下：图1.10原系统与辨识后的系统阶跃响应对比图：当采样时间为0.02，数据长度为400时，系统辨识结果如下：原系统与辨识后的系统阶跃响应对比图：图1-13当采样时间为0.2秒，数据长度为400时，系统辨识结果如下：图1-14 原系统与辨识后的系统阶跃响应对比图：综上所述，结果如表结果可知，对于存在有零点的系统来说，通过面积法辨识系统必须合理的选择分子分母的阶次，否则不能得出正确的辨识结果。