系统辨识经典辨识方法

系统辨识经典辨识方法经典辨识方法报告1.面积法1.1辨识原理1.1.1分子多项式为1的系统11)(111sasasasGnnnn……………………………………………（1.1）由于系统的传递函数与微分方程存在着一一对应的关系，因此，可以通过求取微分方程的系数来辨识系统的传递函数。在求得系统的放大倍数K后，要先得到无因次阶跃响应y(t)（设τ=0）。大多数自衡的工业过程对象的y(t)可以用下式描述来近似1)()()()(a111tydttdyadttydadttydnnnn……………………………（1.2）面积法原则上可以求出n为任意阶的各系数。以n=3为例，注意到1|)(,0|)(d|)(d|)(d23tttttydttydttydtty…………………………（1.3）将式（2.1.2）的y(t)项移至右边，在[0，t]上积分，得tdttytyadttdyadttyda01223)](1[)()()(…………………………………（1.4）定义tdttytF01)](1[)(……………………………………………………………（1.5）则由式（2.1.3）给出的条件可知，在t→∞01)](1[adtty……………………………………………………………（1.6）将式a1y(t)移到等式右边，定义)()]()([)()(a201123tFdttyatFtyadttdyt…………………………………（1.7）利用初始条件（2.1.3）当t→∞时)(a22F……………………………………………………………………（1.8）同理有a3=F3(∞)以此类推，若n≥2，有an=Fn(∞)1.1.2分子、分母分别为m阶和n阶多项式的系统系统辨识经典辨识方法当传递函数的形式如下所示时0111111)()(11)(uhKmnsasasasbsbsbKsGnnnnmmmm…………………………………(1.9)定义1111111111)()(1)(iiimmmmnnnnscsbsbsbsasasasPsPKsG………………………………(1.10)由于0**)](1[)](1[dteththLst…………………………………………(1.11)则)](1[*th的Laplace变换为：111*1)(11)](1[iiiiiisCsCssPsthL……………………………………(1.12)定义一阶面积1A为：111100011lim)](*1[lim)](*1[csCsCthLdtthAiiiiiiss………(1.13)令)1(1)]([1*1scsthL……………………………………………………………(1.14)定义二阶面积为：系统辨识经典辨识方法21220**0012)1)(1()]()([limcscscscdtdhhAiiiiiiist…(1.15)同理，令)...1(1)]([11221*1iiiscscscsthL……………………………………(1.16)定义i阶面积为iicA。由此可得：1111111)1)(1...(1...iiimmmmnnnnsAsbsbsbsasasa…(1.17)上式可写成如下形式：mnnnnmnmnmnnnmnnnnAAAAAAAAAAAAbbb21121211121………………………(1.18)nmnnnAAAbbbAAAAaaa2121121121010010001………………………(1.19)通过该系数矩阵，即可求出传递函数分子分母系数的值。1.2程序设计1.2.1传递函数形式如式1.1的系统取系统传递函数如下：G(s)=13𝑠3+2𝑠2+5.2𝑠+1MATLAB程序如下：clc%清空工作区clear系统辨识经典辨识方法dt=0.01;%设置采样时间t=0:dt:50;%设置时间长度num=1;%此系统分子为1den=[325.21];%分母多项式系数%绘制原传递函数阶跃响应曲线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;for(i=1:Length)sum1=sum1+(1-y(i))*dt;F(i)=sum1;enda1=sum1%求a2sum2=0;for(i=1:Length)sum2=sum2+(F(i)-a1*y(i))*dt;f(i)=sum2;enda2=sum2%求a3sum3=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=[a3a2a11];系统辨识经典辨识方法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=corrcoef(y,h)运行以上程序得到结果如下：原系统传递函数为:G=1-------------------------3s^3+2s^2+5.2s+1Continuous-timetransferfunction.辨识参数结果：a1=5.2048a2=2.0608a3=2.8388系统辨识后的传递函数为:G=1-----------------------------------2.839s^3+2.061s^2+5.205s+1Continuous-timetransferfunction.相关系数：系统辨识经典辨识方法r=1.00000.99990.99991.0000此时原传递函数和辨识所得传递函数的阶跃响应对比如下图：图1.1原传递函数和辨识所得传递函数的阶跃响应对比由上图可以看出，辨识所得结果比较准确。1.2.1传递函数形式如式(1.9)的系统（无噪声）取系统传递函数如下：G(s)=𝑠2+3𝑠+12𝑠3+5.5𝑠2+3𝑠+1MATLAB程序如下：clc%清空工作区cleardt=0.01;%设置采样时间t=0:dt:20;%设置时间长度num=[131];%分子多项式系数den=[25.531];%分母多项式系数%绘制原传递函数阶跃响应曲线fprintf('原系统传递函数为:')G=tf(num,den)y=step(num,den,t);0510152025303540455000.10.20.30.40.50.60.70.80.91t/sy(t)和h(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=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%求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;end系统辨识经典辨识方法A5=sum5%求分子系数b1,b2M=(-1)*(inv([A3,A2;A4,A3]))*[A4;A5];fprintf('分子多项式系数为：')b1=M(1,1)b2=M(2,1)%求分母系数a1,a2,a3N=[100;A110;A2A11]*[b1;b2;0]+[A1;A2;A3];fprintf('分母多项式系数为：')a1=N(1,1)a2=N(2,1)a3=N(3,1)%求辨识所得传递函数num1=[b2b11];den1=[a3a2a11];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)运行以上程序结果如下：原系统传递函数为:G=s^2+3s+1-------------------------2s^3+5.5s^2+3s+1Continuous-timetransferfunction.A1-A5阶面积分别为：A1=0.0076系统辨识经典辨识方法A2=4.3364A3=-9.6285A4=15.7297A5=7.1087分子多项式系数为：b1=7.4409b2=12.8942分母多项式系数为：a1=7.4485a2=17.2869a3=22.7359辨识所得传递函数为:G=12.89s^2+7.441s+1-----------------------------------22.74s^3+17.29s^2+7.448s+1Continuous-timetransferfunction.相关系数：r=1.00000.9996系统辨识经典辨识方法0.99961.0000此时原传递函数和辨识所得传递函数的阶跃响应对比如下图：图1.2原传递函数和辨识所得传递函数的阶跃响应对比由上图可以看出，在未加入噪声之前，采用面积法辨识结果很精确，并且，分子可以为阶次低于分母的任意阶次。1.2.2加入噪声为便于对比，仍取系统传递函数如下：G(s)=𝑠2+3𝑠+12𝑠3+5.5𝑠2+3𝑠+1MATLAB程序如下：clc%清空工作区cleardt=0.01;%设置采样时间t=0:dt:20;%设置时间长度pj=0;fc=sqrt(0.01);%加入方差为0.01的白噪声序列num=[131];%分子多项式系数den=[25.531];%分母多项式系数%绘制原传递函数阶跃响应曲线fprintf('原系统传递函数为:')G=tf(num,den)L=length(t);cs=randn(1,L);%产生噪声cs=cs/std(cs)