系统辨识经典辨识方法

fprintf('原系统传递函数为:')
G=tf(num,den)
L=length(t);
cs=randn(1,L);%产生噪声
cs=cs/std(cs);
cs=cs-mean(cs);
cs=pj+fc*cs;
cs=cs';
y=step(num,den,t)+cs;%加入噪声
Length=length(y);
由上图可以看出，辨识所得结果比较准确。
1.2.1传递函数形式如式(1.9)的系统（无噪声）
取系统传递函数如下：
MATLAB程序如下：
clc%清空工作区
clear
dt=0.01; %设置采样时间
t=0:dt:20; %设置时间长度
num=[1 3 1]; %分子多项式系数
den=[2 5.5 3 1]; %分母多项式系数
legend('原传递函数','辨识所得传递函数');
title('原传递函数与辨识所得传递函数的阶跃响应对比')
grid;
xlabel('t/s');
ylabel('y(t)和h(t)');
fprintf('相关系数：');%求相关系数
r=corrcoef(y,h)
运行以上程序得到结果如下：
原系统传递函数为:
%绘制原传递函数阶跃响应曲线
fprintf('原系统传递函数为:')
G=tf(num,den)
y=step(num,den,t);
Length=length(y);
plot(t,y);
grid;
xlabel('t/s');
ylabel('y(t)');
%辨识程序设计，此系统m+n=5,故应计算A1-A5
a3=sum3
%绘制辨识后的传递函数
dt=0.01;
t=0:dt:50;
num2=1;
den2=[a3 a2 a1 1];
fprintf('系统辨识后的传递函数为:')
G=tf(num2,den2)
h=step(num2,den2,t);%辨识所得传递函数阶跃响应
plot(t,y,'black',t,h,'blue');
den1=[a3 a2 a1 1];
fprintf('辨识所得传递函数为:')
G=tf(num1,den1)
h=step(num1,den1,t);
plot(t,y,'black',t,h,'blue');
legend('原传递函数','辨识所得传递函数');
title('原传递函数与辨识所得传递函数的阶跃响应对比')
end
A5=sum5
%求分子系数b1,b2
M=(-1)*(inv([A3,A2;A4,A3]))*[A4;A5];
fprintf('分子多项式系数为：')
b1=M(1,1)
b2=M(2,1)
%求分母系数a1,a2,a3
N=[1 0 0;A1 1 0;A2 A1 1]*[b1;b2;0]+[A1;A2;A3];
legend('原传递函数','辨识所得传递函数');
title('原传递函数与辨识所得传递函数的阶跃响应对比')
grid;
xlabel('t/s');
ylabel('y(t)');
fprintf('相关系数：'); %求相关系数
r=corrcoef(y,h)
运行程序，结果如下：
原系统传递函数为:
G =
……………………………（1.2）
面积法原则上可以求出n为任意阶的各系数。以n=3为例，注意到
…………………………（1.3）
将式（2.1.2）的y(t)项移至右边，在[0，t]上积分，得
…………………………………（1.4）
定义
……………………………………………………………（1.5）
则由式（2.1.3）给出的条件可知，在t→∞
定义二阶面积为：
…(1.15)
同理，令
……………………………………(1.16)
定义 阶面积为 。由此可得：
…(1.17)
上式可写成如下形式：
………………………(1.18)
………………………(1.19)
通过该系数矩阵，即可求出传递函数分子分母系数的值。
1.2程序设计
1.2.1传递函数形式如式1.1的系统
Continuous-time transfer function.
A1-A5阶面积分别为：
A1 =
0.0076
A2 =
4.3364
A3 =
-9.6285
A4 =
15.7297
A5 =
7.1087
分子多项式系数为：
b1 =
7.4409
b2 =
12.8942
分母多项式系数为：
a1 =
7.4485
a2 =
……………………………………………………………（1.6）
将式a1y(t)移到等式右边，定义
…………………………………（1.7）
利用初始条件（2.1.3）当t→∞时
……………………………………………………………………（1.8）
同理有a3=F3(∞)
以此类推，若n≥2，有an=Fn(∞)
1.1.2分子、分母分别为m阶和n阶多项式的系统
经典辨识方法报告
1.面积法
1.1辨识原理
1.1.1分子多项式为1的系统
……………………………………………（1.1）
由于系统的传递函数与微分方程存在着一一对应的关系，因此，可以通过求取微分方程的系数来辨识系统的传递函数。在求得系统的放大倍数K后，要先得到无因次阶跃响应y(t)（设τ=0）。大多数自衡的工业过程对象的y(t)可以用下式描述来近似
b1=M(1,1)
b2=M(2,1)
%求分母系数a1,a2,a3
N=[1 0 0;A1 1 0;A2 A1 1]*[b1;b2;0]+[A1;A2;A3];
fprintf('分母多项式系数为：')
a1=N(1,1)
a2=N(2,1)
a3=N(3,1)
%求辨识所得传递函数
num1=[b2 b1 1];
D(i)=sum4;
end
A4=sum4
%求A5
sum5=0;
for(i=1:Length-1)
sum5=sum5+(D(i)-A4*(y(i)+y(i+1))/2)*dt;
end
A5=sum5
%求分子系数b1,b2
M=(-1)*(inv([A3,A2;A4,A3]))*[A4;A5];
fprintf('分子多项式系数为：')
G =
1
-------------------------
3 s^3 + 2 s^2 + 5.2 s + 1
Continuous-time transfer function.
辨识参数结果：
a1 =
5.2048
a2 =
2.0608
a3 =
2.8388
系统辨识后的传递函数为:
G =
1
-----------------------------------
F(i)=sum1;
end
a1=sum1
%求a2
sum2=0;
for(i=1:Length)
sum2=sum2+(F(i)-a1*y(i))*dt;
f(i)=sum2;
end
a2=sum2
%求a3
sum3=0;
for(i=1:Length)
sum3=sum3+(f(i)-a2*y(i))*dt;
end
17.2869
a3 =
22.7359
辨识所得传递函数为:
G =
12.89 s^2 + 7.441 s + 1
-----------------------------------
22.74 s^3 + 17.29 s^2 + 7.448 s + 1
Continuous-time transfer function.
相关系数：
r =
1.0000 0.9996
0.9996 1.0000
此时原传递函数和辨识所得传递函数的阶跃响应对比如下图：
图1.2原传递函数和辨识所得传递函数的阶跃响应对比
由上图可以看出，在未加入噪声之前，采用面积法辨识结果很精确，并且，分子可以为阶次低于分母的任意阶次。
1.2.2加入噪声
为便于对比，仍取系统传递函数如下：
MATLAB程序如下：
clc%清空工作区
clear
dt=0.01; %设置采样时间
t=0:dt:20; %设置时间长度
pj=0;
fc=sqrt(0.01);%加入方差为0.01的白噪声序列
num=[1 3 1]; %分子多项式系数
den=[2 5.5 3 1]; %分母多项式系数
%绘制原传递函数阶跃响应曲线
取系统传递函数如下：
MATLAB程序如下：
clc%清空工作区
clear
dt=0.01;%设置采样时间
t=0:dt:50;%设置时间长度
num=1;%此系统分子为1
den=[3 2 5.2 1];%分母多项式系数
%绘制原传递函数阶跃响应曲线
fprintf('原系统传递函数为:')
G=tf(num,den)
grid;
xlabel('t/s');
ylabel('y(t)');
fprintf('相关系数：'); %求相关系数
r=corrcoef(y,h)
运行以上程序结果如下：
原系统传递函数为:
G =
s^2 + 3 s + 1
-------------------------
2 s^3 + 5.5 s^2 + 3 s + 1
2.839 s^3 + 2.061 s^2 + 5.205 s + 1
Continuous-time transfer function.
相关系数：
r =
1.0000 0.9999
0.9999 1.0000
此时原传递函数和辨识所得传递函数的阶跃响应对比如下图：
图1.1原传递函数和辨识所得传递函数的阶跃响应对比
fprintf('A1-A5阶面积分别为：')
%求A1
sum1=0;
for(i=1:Length-1)
sum1=sum1+(1-(y(i)+y(i+1))/2)*dt;
A(i)=sum1;
end
A1=sum1
%求A2
sum2=0;
for(i=1:Length-1)
sum2=sum2+(A(i)-A1*(y(i)+y(i+1))/2)*dt;
当传递函数的形式如下所示时
…………………………………(1.9)
定义
………………………………(1.10)
由于 …………………………………………(1.11)
则 的Laplace变换为：
……………………………………(1.12)
定义一阶面积 为：
………(1.13)
令
……………………………………………………………(1.14)
A(i)=sum1;
end
A1=sum1
%求A2
sum2=0;
for(i=1:Length-1)
sum2=sum2+(A(i)-A1*(y(i)+y(i+1))/2)*dt;
B(i)=sum2;
end
A2=sum2
%求A3
sum3=0;
for(i=1:Length-1)
sum3=sum3+(B(i)-A2*(y(i)+y(i+1))/2)*dt;
C(i)=sum3;
end
A3=sum3
%求A4
sum4=0;
for(i=1:Length-1)
sum4=sum4+(C(i)-A3*(y(i)+y(i+1))/2)*dt;
D(i)=sum4;
end
A4=sum4
%求A5
sum5=0;
for(i=1:Length-1)
sum5=sum5+(D(i)-A4*(y(i)+y(i+1))/2)*dt;
fprintf('分母多项式系数为：')
a1=N(1,1)
a2=N(2,1)
a3=N(3,1)
%求辨识所得传递函数
num1=[b2 b1 1];
den1=[a3 a2 a1 1];
fprintf('辨识所得传递函数为:')
G=tf(num1,den1)
h=step(num1,den1,t);
plot(t,y,'black',t,h,'blue');
y=step(num,den,t);
Length=length(y);%数据长度
plot(t,y);
grid;
xlabel('t/s');
ylabel('y(t)');
%进行辨识设计
fprintf('辨识参数结果：');
%求a1
sum1=0;
for(i=1:Length)
sum1=sum1+(1-y(i))*dt;
B(i)=sum2;
end
A2=sum2
%求A3
sum3=0;
for(i=1:Length-1)
sum3=sum3+(B(i)-A2*(y(i)+y(i+1))/2)*dt;
C(i)=sum3;
end
A3=sum3
%求A4
sum4=0;
for(i=1:Length-1)
sum4=sum4+(C(i)-A3*(y(i)+y(i+1))/2)*dt;
s^2 + 3 s +ቤተ መጻሕፍቲ ባይዱ1
-------------------------
plot(t,y);
grid;
xlabel('t/s');
ylabel('y(t)');
%辨识程序设计，此系统m+n=5,故应计算A1-A5
fprintf('A1-A5阶面积分别为：')
%求A1
sum1=0;
for(i=1:Length-1)
sum1=sum1+(1-(y(i)+y(i+1))/2)*dt;