阶跃响应法——面积辨识法

转自：http://ljd19900603.blog.163.com/blog/static/13803759420139179925635/

作业要求：

1.   假设模型求阶跃响应与脉冲响应

2.   软件计算求出响应曲线

3.   然后用所讲的一种方法来辨识求出传递函数以及误差

4.   在输入端加误差，由测量值得到传递函数与原来模型作比较

 

 

原模型分析

定义函数模型：

G(s) = 1/(s + 1)/(2*s + 1)/(3*s + 1)=1/(6*s3 + 11*s2 + 6*s + 1)

稳定性分析：

由于对模型进行分析时必需要验证系统是稳定的，这是进行控制辨识的基础。在定义的函数模型中可以看出，系统的极点全部在S的左半平面，所以是稳定的。

软件判定：stable.m

clc                    %清屏    

clear                  %清除变量

sym  s                %定义S

num = [1];              %分子

den = [6 11 6 1];         %分母

G = tf(num,den);         %传递函数

G1 = Zpk(G);             %把传递函数转换为零极点形式

Pole = G1.p{:};             %求极点

Ii = find(real(pole) >= 0)      %判断极点的实部是否在S右平面

If ii > 0 disp{‘模型不稳定’}    %若S右平面极点个数大于零不稳定

  else   disp{‘模型稳定’}    %没有S右平面极点则稳定

end

实验结果：模型稳定

 

脉冲响应

软件计算：

impulse.m

clc                      %清屏

clear                    %清除变量

sym  s                  %定义S

num = [1];                %分子

den = [6 11 6 1];           %分母

G = tf(num,den);           %传递函数

t0 = 0.01;                 %步长

t = 0 : t0 : 300;             %时间范围

Length = 30000;           %采样次数

y = impulse(G,t)            %脉冲响应

plot(t,y),grid;         %画出对应时间下的脉冲响应曲线，并加网格

title(‘原函数的脉冲响应曲线’)    %显示的标题

阶跃响应：

软件计算：

step.m

clc                      %清屏

clear                    %清除变量

sym  s                  %定义S

num = [1];                %分子

den = [6 11 6 1];           %分母

G = tf(num,den);           %传递函数

t0 = 0.01;                 %步长

t = 0 : t0 : 300;             %时间范围

y = step(G,t)               %阶跃响应

plot(t,y),grid;          %画出对应时间下的脉冲响应曲线，并加网格

title(‘原函数的阶跃响应曲线’)   %显示的标题

e(ss) = dcgain(G)             %求稳态误差

e(ss) = 2%

 

  

二、无噪声时的辨识分析

 

采用阶跃响应辨识的步骤：

第一步：实际测取过程的阶跃响应

第二步：由阶跃响应求过程的传递函数

常用方法：近似法、切线法、两点法、面积法，在这里我选用面积法来辨识。

面积法的原理：

考虑传递函数如下G（s）= 1/(an* sn + an-1* sn-1 +…+ a1*s+1),系统的传递函数与微分方程存在着一一对应关系，因此，可以通过求取微分方程的系数来辨识系统的传递函数。在求得系统的放大倍数以后，我们要得到无因次阶跃响应y(t). 大多数自平衡的工业控制对象的y(t)可近似为an*dny(t)/dtn+an-1*dn-1y(t)/dtn-1+······+a1*dy(t)/dt+y(t)=u(t)面积法原则上可求出民为任意阶的各系数，以n = 3为例，

将y(t)移至右边，在[0,t]上积分

a3*d2y(t)/dt3+a2*dy(t)/dt+a1y(t)=∫t0[1-y(t)]dt

由初始条件可知，当t趋于无穷时，

d3y(t)/dt3 =d2y(t)/dt2=dy(t)/dt=0，y(t)=1

所以当t趋于无穷时,a1=∫t0[1-y(t)]dt;

将式a1y(t)移到右边，定义a3*dy(t)/dt+a2*dy(t)=∫t0[a1 - a1y(t)]dt

                        a2 = ∫t0[a1 - a1y(t)]dt

同理可求:                a3 = ∫t0[a2 – a2y(t)]dt

软件计算：

Identification1.m

clc                      %清屏

clear                    %清除变量

sym  s                  %定义S

num = [1];                %分子

den = [6 11 6 1];           %分母

G = tf(num,den);           %传递函数

y = ilaplace(G)             %反拉氏变换求得微分方程

t0 = 0.05;                  %步长

t1 = 0:t0:30000;            %采样范围

y1 = subs(y,t,t1);            %求取每个采样点下的函数值

 

k1 = 0; k2 = 0; k3 = 0;        %系数初值

for i = 1：29999             %求系数a1，求积分变成求和

k1 = k1 + (1 – (y1(i) + y1(i + 1))/2) * t0   %求连续两个值的均值

end                                 %作为一次采集数

a1 = k1; k1 =0;                 %求得第一个系数a1

 

for i = 1; 29999                 %求系数a2

k1 = k1 + (1 – (y1(i) + y1(i + 1))/2) * t0

k2 = k2 + (k1 – a1*(y1(i) + y1(i + 1))/2) * t0

end

a2 = k2; k1 = 0;k2 = 0;

 

for i = 1; 29999                  %求系数a3

k1 = k1 + (1 – (y1(i) + y1(i + 1))/2) * t0

k2 = k2 + (k1 – a1 * (y1(i) + y1(i + 1))/2) * t0

k3 = k3 + (k3 – a2 * (y1(i) + y1(i + 1))/2) * t0

end

a3 = k3;

 

num1 = [1]                       %辨识的传递函数分子不变

den1 = [a3 a2 a1 1]              %辨识的传递函数母子

disp(‘辨识系统传递函数：’)

G1 = tf(num1,den1)                %求得辨识的传递函数

step(G);hold on;                 %求两个传递函数的阶跃响应

step(G1);

grid

legend(‘给定系统阶跃响应曲线’，‘辨识系统阶跃响应曲线’)

%添加图例的标注

a = step(G);

b = step(G1);

for i = 1: 100

 c(i) = a(i);

 d(i) = b(i);

end

R = corrcoef(c,d)   %相关系数，表征两个系统的相似度 

辨识结果：

G1(s) = 1/(5.9862*s^3 + 10.9886*s^2 + 5.9896*s + 1 )

R = 0.9988

从结果可以看出，用面积法辨识三阶定常经性系统的传递函数时，误差不大，但其中一个很重要的影响因素是采集间隔t0。

 

三、有噪声时的辨识分析

 

在有噪声时的辨识相对与没有噪声的辨识，对辨识方法的精确度有了更大的要求。在具有的辨识实现方法上，在输入端加一个Gn(s) = 1 / s的扰动即可。软件方面只需要加：G(s) = G(s) * Gn(s)

 

辨识结果：

G1(s) = 1/(5.1012*s^3 + 9.9024*s^2 + 5.123*s + 1 )

R = 0.8502

由以上结果可以看出，加噪声后辨识效果不及无噪声的效果。