闭环系统可辨识性实验研究

1、研1206 刘新菊 模式识别与智能系统 2012020176实验三 闭环系统可辨识性实验研究实验目的：通过实验掌握闭环系统可辨识性条件及其应用。实验内容：1. 模拟如下闭环系统前向通道方程其中v(k)是服从正态分布N（0，1）的高斯白噪声。反馈通道方程1） 二阶反馈且存在反馈通道噪声其中w(k)是服从正态分布N（0,1）的高斯白噪声。2） 二阶反馈但不加反馈噪声3） 一阶反馈并加反馈噪声4） 一阶反馈不加反馈噪声5） 比例反馈加反馈噪声6） 比例反馈不加反馈噪声2. 在以上各种不同的闭环状态下，利用最小二乘法直接辨识方法估计方法估计前向通道的模型参数。3. 根据实验结果分析闭环系统的可辨识性条

2、件。实验步骤：（1） 熟悉闭环系统的可辨识性概念及条件。（2） 设计实验方案，编制试验程序，上机调剂，记录结果。（3） 分析实验结果，完成书面实验报告，要求同时提供程序框图及其说明。实验三 闭环系统可辨识性实验报告一、熟悉闭环系统的可辨识性概念及条件系统辨识中的闭环问题作为一类特殊的辨识问题，近年来越来越受到人们的关注。Soderstrom 指出，由于输出信号的干扰噪声通过反馈环节与输入信号相关，直接采用频谱分析法，辨识结果将是对象传递函数与反馈传递函数倒数间的一个加权平均值；Ljung 的理论分析表明，如果采用预报误差法进行开环辨识时，只要对象模型集包含真实对象的动态特性

3、，即使噪声模型不足以描述噪声的真实动态特性，仍可获得对象的一致无偏估计，而将此方法直接用于闭环辨识时，只要噪声模型不能精确描述真实噪声，即使对象模型集包含了真实对象动态特性，得到的将是对象参数的有偏估计；Gustavsson 则举出了反馈环节是比例调节器时，直接采用预报误差法将导致对象不可辨识性的实例。  对上述问题的解决，存在三种经典的解决方法：直接法、间接法、联合输入输出法。 (1)直接辨识法。直接法是利用闭环条件下得到的数据，直接把对象当作在开环环境中进行辨识。但要保证系统的可辨识性及对象参数估计的一致无偏性，采用直接法时必须满足条件：存在足够的外

4、部激励信号；控制器的阶数足够高；控制器在不同的模式间切换；控制器时变、非线性。 (2)间接辨识法。间接法的思想是先根据闭环系统的输入输出信号用开环辨识方法得到闭环系统的传递函数，再由已知的控制器传递函数推导出对象的辨识模型。 (3)联合输入输出法。将闭环系统的输入输出统一看作为白噪声驱动系统产生的输出信号。先辨识出整个系统的传递函数矩阵，然后利用传递函数矩阵中的分块矩阵与系统各环节传递函数的关系，求出各环节的传递函数。 经典闭环辨识法在许多方面并不尽如人意：直接法虽然简洁,但必须满足较为苛刻的条件；间接法与联合输入输出法实验前要知道反馈控制器的传递函数，并且需要至

5、少三个可检测的信号，此外，闭环系统与开环对象模型间的参数转换导致了繁重的计算工作量。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) %*clearclc%*产生均值为0，方差为1 的正态分布噪声*v=random('Normal',0,1,1,402);w=random('Normal',0,1,1,402);%*产生观测序列z*z=zeros(402,1

6、);u=zeros(402,1);z(1)=-1'z(2)=0;u(1)=-1'u(2)=0;for i=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

7、;3;3;3;% K=zeros(4,400); %增益矩阵K=10;10;10;10;for i=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

8、(1)plot(i,Theta(1,:),i,Theta(2,:),i,Theta(3,:),i,Theta(4,:)title('待估参数过渡过程')gridfigure(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

9、.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) %*clearclc%*产生均值为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;for i=3:2002 z(i)=-1.4*z(i-1)-

10、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;for i=3:2002 h=-z(i-1);-z(i-2);u(i-1);u(i-

11、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'

12、,i,Theta(4,:),'b*')legend('a1','a2','b1','b2')title('待估参数过渡过程')gridfigure(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、 一阶反馈并加反馈噪声%*%闭环

13、系统%最小二乘辨识的递推算法%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) %*clearclc%*产生均值为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;for i=3:2002 z(i)=-

14、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;for i=3:2002 h=-z(i-1);-z(i-2);u(i-1);u(i-2)

15、; 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

16、,Theta(4,:),'b*')legend('a1','a2','b1','b2')title('待估参数过渡过程')gridfigure(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.01270.6569一阶反馈并加反馈噪声模型参数辨识过程4、 一阶反馈不加反馈噪声%*%闭环系统

17、%最小二乘辨识的递推算法%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)%*clearclc%*产生均值为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;for i=3:2002 z(i)=-1.4*z(i

18、-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;for i=3:2002 h=-z(i-1);-z(i-2);u(i-1);u(i-2); K=P*h*inv(

19、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,:

20、),'b*')legend('a1','a2','b1','b2')title('待估参数过渡过程')gridfigure(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.05780.6345一阶反馈并加反馈噪声模型参数辨识过程5、 比例反馈加反馈噪声%*%闭环系统%最小二乘辨识的递推算

21、法%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)%*clearclc%=产生均值为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;for i=3:2002 z(i)=-1.4*z(i-1)-0.45*z(i-2

22、)+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;for i=3:2002 h=-z(i-1);-z(i-2);u(i-1);u(i-2); K=P*h*inv(h'*P*h+1); T

23、heta(:,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*')l

24、egend('a1','a2','b1','b2')title('待估参数过渡过程')gridfigure(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

25、-1)-0.45*Z(k-2)+u(k-1)+0.7*u(k-2)+v(k)%u(k)=1.2*z(k)%*clearclc%*产生均值为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;for i=3:2002 z(i)=-1.4*z(i-1)-0.45*z(i-2)+u(i-1)+0.7*u(i-2)+

26、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;for i=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