闭环系统可辨识性实验研究.docx
《闭环系统可辨识性实验研究.docx》由会员分享，可在线阅读，更多相关《闭环系统可辨识性实验研究.docx（23页珍藏版）》请在冰豆网上搜索。
闭环系统可辨识性实验研究
实验三闭环系统可辨识性实验研究
实验目的：
通过实验掌握闭环系统可辨识性条件及其应用。
实验内容：
1.模拟如下闭环系统
前向通道方程
其中v（k）是服从正态分布N（0，1）的高斯白噪声。
反馈通道方程
1）二阶反馈且存在反馈通道噪声
其中w（k）是服从正态分布N（0,1）的高斯白噪声。
2）二阶反馈但不加反馈噪声
3）一阶反馈并加反馈噪声
4）一阶反馈不加反馈噪声
5）比例反馈加反馈噪声
6）比例反馈不加反馈噪声
2.在以上各种不同的闭环状态下，利用最小二乘法直接辨识方法估计方法估计前向通道的模型参数。
3.根据实验结果分析闭环系统的可辨识性条件。
实验步骤：
（1）熟悉闭环系统的可辨识性概念及条件。
（2）设计实验方案，编制试验程序，上机调剂，记录结果。
（3）分析实验结果，完成书面实验报告，要求同时提供程序框图及其说明。
实验三闭环系统可辨识性实验报告
一、熟悉闭环系统的可辨识性概念及条件
系统辨识中的闭环问题作为一类特殊的辨识问题，近年来越来越受到人们的关注。
Soderstrom 指出，由于输出信号的干扰噪声通过反馈环节与输入信号相关，直接采用频谱分析法，辨识结果将是对象传递函数与反馈传递函数倒数间的一个加权平均值；Ljung 的理论分析表明，如果采用预报误差法进行开环辨识时，只要对象模型集包含真实对象的动态特性，即使噪声模型不足以描述噪声的真实动态特性，仍可获得对象的一致无偏估计，而将此方法直接用于闭环辨识时，只要噪声模型不能精确描述真实噪声，即使对象模型集包含了真实对象动态特性，得到的将是对象参数的有偏估计；Gustavsson 则举出了反馈环节是比例调节器时，直接采用预报误差法将导致对象不可辨识性的实例。
对上述问题的解决，存在三种经典的解决方法：
直接法、间接法、联合输入输出法。
（1）直接辨识法。
直接法是利用闭环条件下得到的数据，直接把对象当作在开环环境中进行辨识。
但要保证系统的可辨识性及对象参数估计的一致无偏性，采用直接法时必须满足条件：
①存在足够的外部激励信号；②控制器的阶数足够高；③控制器在不同的模式间切换；④控制器时变、非线性。
（2）间接辨识法。
间接法的思想是先根据闭环系统的输入输出信号用开环辨识方法得到闭环系统的传递函数，再由已知的控制器传递函数推导出对象的辨识模型。
（3）联合输入输出法。
将闭环系统的输入输出统一看作为白噪声驱动系统产生的输出信号。
先辨识出整个系统的传递函数矩阵，然后利用传递函数矩阵中的分块矩阵与系统各环节传递函数的关系，求出各环节的传递函数。
经典闭环辨识法在许多方面并不尽如人意：
直接法虽然简洁,但必须满足较为苛刻的条件；间接法与联合输入输出法实验前要知道反馈控制器的传递函数，并且需要至少三个可检测的信号，此外，闭环系统与开环对象模型间的参数转换导致了繁重的计算工作量。
1、二阶反馈且存在反馈通道噪声
%*******************************
%闭环系统
%最小二乘辨识的递推算法
%Z（k）=-1.4*Z（k-1）-0.45*Z（k-2）+u（k-1）+0.7*u（k-2）+v（k）
%u（k）=0.33*z（k）+0.033*z（k-1）-0.4*z（k）+w（k）
%*******************************
clear
clc
%**************产生均值为0，方差为1的正态分布噪声**************
v=random（'Normal',0,1,1,402）;
w=random（'Normal',0,1,1,402）;
%************产生观测序列z*******************
z=zeros（402,1）;
u=zeros（402,1）;
z
（1）=-1';
z
（2）=0;
u
（1）=-1';
u
（2）=0;
fori=3:
402
z（i）=-1.4*z（i-1）-0.45*z（i-2）+u（i-1）+0.7*u（i-2）+v（i）;
u（i）=0.33*z（i）+0.033*z（i-1）-0.4*z（i-2）+w（i）;
end
%************递推求解************
P=100*eye（4）;%估计方差
Pstore=zeros（4,401）;
Pstore（:
1）=[P（1,1）,P（2,2）,P（3,3）,P（4,4）];
Theta=zeros（4,401）;%参数的估计值，存放中间过程估值
Theta（:
1）=[3;3;3;3];
%K=zeros（4,400）;%增益矩阵
K=[10;10;10;10];
fori=3:
402
h=[-z（i-1）;-z（i-2）;u（i-1）;u（i-2）];
K=P*h*inv（h'*P*h+1）;
Theta（:
i-1）=Theta（:
i-2）+K*（z（i）-h'*Theta（:
i-2））;
P=（eye（4）-K*h'）*P;
Pstore（:
i-1）=[P（1,1）,P（2,2）,P（3,3）,P（4,4）];
end
%***************************************************
disp（'参数a1、a2、b1、b2估计结果：
'）
Theta（:
401）
i=1:
401;
figure
（1）
plot（i,Theta（1,:
）,i,Theta（2,:
）,i,Theta（3,:
）,i,Theta（4,:
））
title（'待估参数过渡过程'）
grid
figure
（2）
plot（i,Pstore（1,:
）,i,Pstore（2,:
）,i,Pstore（3,:
）,i,Pstore（4,:
））
title（'估计方差变化过程'）
Grid
参数a1、a2、b1、b2估计结果：
ans=
1.3842
0.4367
1.0082
0.6828
二阶反馈且存在反馈通道噪声模型参数辨识过程
2、二阶反馈但不加反馈噪声
%*******************************
%闭环系统
%最小二乘辨识的递推算法
%Z（k）=-1.4*Z（k-1）-0.45*Z（k-2）+u（k-1）+0.7*u（k-2）+v（k）
%u（k）=0.33*z（k）+0.033*z（k-1）-0.4*z（k）
%*******************************
clear
clc
%**************产生均值为0，方差为1的正态分布噪声**************
v=random（'Normal',0,1,1,2002）;
%w=random（'Normal',0,1,1,402）;
%************产生观测序列z*******************
z=zeros（2002,1）;
u=zeros（2002,1）;
z
（1）=-1';
z
（2）=0;
u
（1）=-1';
u
（2）=0;
fori=3:
2002
z（i）=-1.4*z（i-1）-0.45*z（i-2）+u（i-1）+0.7*u（i-2）+v（i）;
u（i）=0.33*z（i）+0.033*z（i-1）-0.4*z（i-2）;
end
%************递推求解************
P=100*eye（4）;%估计方差
Pstore=zeros（4,2001）;
Pstore（:
1）=[P（1,1）,P（2,2）,P（3,3）,P（4,4）];
Theta=zeros（4,2001）;%参数的估计值，存放中间过程估值
Theta（:
1）=[3;3;3;3];
%K=zeros（4,400）;%增益矩阵
K=[10;10;10;10];
fori=3:
2002
h=[-z（i-1）;-z（i-2）;u（i-1）;u（i-2）];
K=P*h*inv（h'*P*h+1）;
Theta（:
i-1）=Theta（:
i-2）+K*（z（i）-h'*Theta（:
i-2））;
P=（eye（4）-K*h'）*P;
Pstore（:
i-1）=[P（1,1）,P（2,2）,P（3,3）,P（4,4）];
end
%***************************************************
disp（'参数a1、a2、b1、b2估计结果：
'）
Theta（:
2001）
i=1:
2001;
figure
（1）
plot（i,Theta（1,:
）,'r',i,Theta（2,:
）,'r*',i,Theta（3,:
）,'b',i,Theta（4,:
）,'b*'）
legend（'a1','a2','b1','b2'）
title（'待估参数过渡过程'）
grid
figure
（2）
plot（i,Pstore（1,:
）,i,Pstore（2,:
）,i,Pstore（3,:
）,i,Pstore（4,:
））
title（'估计方差变化过程'）
grid
参数a1、a2、b1、b2估计结果：
ans=
1.3812
0.4188
0.9284
0.7231
二阶反馈但不加反馈噪声模型辨识过程
3、一阶反馈并加反馈噪声
%*******************************
%闭环系统
%最小二乘辨识的递推算法
%Z（k）=-1.4*Z（k-1）-0.45*Z（k-2）+u（k-1）+0.7*u（k-2）+v（k）
%u（k）=z（k）+0.2*z（k-1）+w（k）
%*******************************
clear
clc
%*************产生均值为0，方差为1的正态分布噪声*************
v=random（'Normal',0,1,1,2002）;
w=random（'Normal',0,1,1,2002）;
%************产生观测序列z*******************
z=zeros（2002,1）;
u=zeros（2002,1）;
z
（1）=-1';
z
（2）=0;
u
（1）=-1';
u
（2）=0;
fori=3:
2002
z（i）=-1.4*z（i-1）-0.45*z（i-2）+u（i-1）+0.7*u（i-2）+v（i）;
u（i）=z（i）+0.2*z（i-1）+w（i）;
end
%************递推求解************
P=100*eye（4）;%估计方差
Pstore=zeros（4,2001）;
Pstore（:
1）=[P（1,1）,P（2,2）,P（3,3）,P（4,4）];
Theta=zeros（4,2001）;%参数的估计值，存放中间过程估值
Theta（:
1）=[3;3;3;3];
%K=zeros（4,400）;%增益矩阵
K=[10;10;10;10];
fori=3:
2002
h=[-z（i-1）;-z（i-2）;u（i-1）;u（i-2）];
K=P*h*inv（h'*P*h+1）;
Theta（:
i-1）=Theta（:
i-2）+K*（z（i）-h'*Theta（:
i-2））;
P=（eye（4）-K*h'）*P;
Pstore（:
i-1）=[P（1,1）,P（2,2）,P（3,3）,P（4,4）];
end
%***************************************************
disp（'参数a1、a2、b1、b2估计结果：
'）
Theta（:
2001）
i=1:
2001;
figure
（1）
plot（i,Theta（1,:
）,'r',i,Theta（2,:
）,'r*',i,Theta（3,:
）,'b',i,Theta（4,:
）,'b*'）
legend（'a1','a2','b1','b2'）
title（'待估参数过渡过程'）
grid
figure
（2）
plot（i,Pstore（1,:
）,i,Pstore（2,:
）,i,Pstore（3,:
）,i,Pstore（4,:
））
title（'估计方差变化过程'）
grid
参数a1、a2、b1、b2估计结果：
ans=
1.4047
0.4392
1.0127
0.6569
一阶反馈并加反馈噪声模型参数辨识过程
4、一阶反馈不加反馈噪声
%*******************************
%闭环系统
%最小二乘辨识的递推算法
%Z（k）=-1.4*Z（k-1）-0.45*Z（k-2）+u（k-1）+0.7*u（k-2）+v（k）
%u（k）=z（k）+0.2*z（k-1）
%*******************************
clear
clc
%**************产生均值为0，方差为1的正态分布噪声**************
v=random（'Normal',0,1,1,2002）;
%w=random（'Normal',0,1,1,2002）;
%************产生观测序列z*******************
z=zeros（2002,1）;
u=zeros（2002,1）;
z
（1）=-1';
z
（2）=0;
u
（1）=-1';
u
（2）=0;
fori=3:
2002
z（i）=-1.4*z（i-1）-0.45*z（i-2）+u（i-1）+0.7*u（i-2）+v（i）;
u（i）=z（i）+0.2*z（i-1）;
end
%************递推求解************
P=100*eye（4）;%估计方差
Pstore=zeros（4,2001）;
Pstore（:
1）=[P（1,1）,P（2,2）,P（3,3）,P（4,4）];
Theta=zeros（4,2001）;%参数的估计值，存放中间过程估值
Theta（:
1）=[3;3;3;3];
%K=zeros（4,400）;%增益矩阵
K=[10;10;10;10];
fori=3:
2002
h=[-z（i-1）;-z（i-2）;u（i-1）;u（i-2）];
K=P*h*inv（h'*P*h+1）;
Theta（:
i-1）=Theta（:
i-2）+K*（z（i）-h'*Theta（:
i-2））;
P=（eye（4）-K*h'）*P;
Pstore（:
i-1）=[P（1,1）,P（2,2）,P（3,3）,P（4,4）];
end
%***************************************************
%***************************************************
disp（'参数a1、a2、b1、b2估计结果：
'）
Theta（:
2001）
i=1:
2001;
figure
（1）
plot（i,Theta（1,:
）,'r',i,Theta（2,:
）,'r*',i,Theta（3,:
）,'b',i,Theta（4,:
）,'b*'）
legend（'a1','a2','b1','b2'）
title（'待估参数过渡过程'）
grid
figure
（2）
plot（i,Pstore（1,:
）,i,Pstore（2,:
）,i,Pstore（3,:
）,i,Pstore（4,:
））
title（'估计方差变化过程'）
grid
参数a1、a2、b1、b2估计结果：
ans=
8.4624
1.7891
8.0578
0.6345
一阶反馈并加反馈噪声模型参数辨识过程
5、比例反馈加反馈噪声
%*******************************
%闭环系统
%最小二乘辨识的递推算法
%Z（k）=-1.4*Z（k-1）-0.45*Z（k-2）+u（k-1）+0.7*u（k-2）+v（k）
%u（k）=1.2*z（k）+w（k）
%*******************************
clear
clc
%===========产生均值为0，方差为1的正态分布噪声=========
v=random（'Normal',0,1,1,2002）;
w=random（'Normal',0,1,1,2002）;
%************产生观测序列z*******************
z=zeros（2002,1）;
u=zeros（2002,1）;
z
（1）=-1';
z
（2）=0;
u
（1）=-1';
u
（2）=0;
fori=3:
2002
z（i）=-1.4*z（i-1）-0.45*z（i-2）+u（i-1）+0.7*u（i-2）+v（i）;
u（i）=1.2*z（i）+w（i）;
end
%************递推求解************
P=100*eye（4）;%估计方差
Pstore=zeros（4,2001）;
Pstore（:
1）=[P（1,1）,P（2,2）,P（3,3）,P（4,4）];
Theta=zeros（4,2001）;%参数的估计值，存放中间过程估值
Theta（:
1）=[3;3;3;3];
%K=zeros（4,400）;%增益矩阵
K=[10;10;10;10];
fori=3:
2002
h=[-z（i-1）;-z（i-2）;u（i-1）;u（i-2）];
K=P*h*inv（h'*P*h+1）;
Theta（:
i-1）=Theta（:
i-2）+K*（z（i）-h'*Theta（:
i-2））;
P=（eye（4）-K*h'）*P;
Pstore（:
i-1）=[P（1,1）,P（2,2）,P（3,3）,P（4,4）];
end
%***************************************************
%***************************************************
disp（'参数a1、a2、b1、b2估计结果：
'）
Theta（:
2001）
i=1:
2001;
figure
（1）
plot（i,Theta（1,:
）,'r',i,Theta（2,:
）,'r*',i,Theta（3,:
）,'b',i,Theta（4,:
）,'b*'）
legend（'a1','a2','b1','b2'）
title（'待估参数过渡过程'）
grid
figure
（2）
plot（i,Pstore（1,:
）,i,Pstore（2,:
）,i,Pstore（3,:
）,i,Pstore（4,:
））
title（'估计方差变化过程'）
grid
参数a1、a2、b1、b2估计结果：
ans=
1.4508
0.4730
1.0414
0.7135
比例反馈加反馈噪声模型参数辨识过程
6、比例反馈不加反馈噪声
%*******************************
%闭环系统
%最小二乘辨识的递推算法
%Z（k）=-1.4*Z（k-1）-0.45*Z（k-2）+u（k-1）+0.7*u（k-2）+v（k）
%u（k）=1.2*z（k）
%*******************************
clear
clc
%**************产生均值为0，方差为1的正态分布噪声**************
v=random（'Normal',0,1,1,2002）;
w=random（'Normal',0,1,1,2002）;
%************产生观测序列z*******************
z=zeros（2002,1）;
u=zeros（2002,1）;
z
（1）=-1';
z
（2）=0;
u
（1）=-1';
u
（2）=0;
fori=3:
2002
z（i）=-1.4*z（i-1）-0.45*z（i-2）+u（i-1）+0.7*u（i-2）+v（i）;
u（i）=1.2*z（i）;
end
%************递推求解************
P=100*eye（4）;%估计方差
Pstore=zeros（4,2001）;
Pstore（:
1）=[P（1,1）,P（2,2）,P（3,3）,P（4,4）];
Theta=zeros（4,2001）;%参数的估计值，存放中间过程估值
Theta（:
1）=[3;3;3;3];
%K=zeros（4,400）;%增益矩阵
K=[10;10;10;10];
fori=3:
2002
h=[-z（i-1）;-z（i-2）;u（i-1）;u（i-2）];
K=P*h*inv（h'*P*h+1）;
Theta（:
i-1）=Theta（:
i-2）+K*（z（i）-h'*Theta（:
i-2））;
P=（eye（4）-K*h'）*P;
Pstore（:
i-1）=[P（1,1）,P（2,2）,P（3,3）,P（4,4）];
end
%***************************************************
disp（'参数a1、a2、b1、b2估计结果：
'）
Theta（:
2001）
i=1:
2001;
figure
（1）
plot（i,Theta（1,:
）,'r',i,Theta（2,:
）,'r*',i,Theta（3,:
）,'b',i,Theta（4,:
）,'b*'）
legend（'a1','a2','b1','b2'）
title（'待估参数过渡过程'）
grid
figure
（2）
plot（i,Pstore（1,:
）,i,Pstore（2,:
）,i,Pstore（3,:
）,i,Pstore（4,:
））
title（'估计方差变化过程'）
grid
参数a1、a2、b1、b2估计结果：
ans=
3.3309
3.5410
2.6029
3.2713
比例反馈不加反馈噪声模型参数辨识过程
闭环辨识结果
反馈律
反馈噪声
参数估计值
可辨识性
a1（1.4）
a2（0.45）
b1（1.0）
b2（0.7）
2阶反馈：
u（k）=0.33*z（k）+0.033*z（k-1）-0.4*z（k）
w（k）
0
1.3842
0.4367
1.0082
0.6828
唯一，一致
w（k）=0
1.3812
0.4188
0.9284
0.7231
唯一，一致
1阶反馈：
u（k）=z（k）+0.2*z（k-1）
w（k）
0
1.4047
0.4392
1.0127
0.6569
唯一，一致
w（k）=0
8.4624
1.7891
8.0578
0.6345
不可辨识
升级会员 