系统辨识及自适应控制实验

1、中南大学系统辨识及自适应控制实验指导老师 贺建军 姓名 史伟东 专业班级 测控1102班0909111814号 实验日期 2014年11月 实验一 递推二乘法参数辨识设被辨识系统的数学模型由下式描述：式中x(k)为方差为0.1的白噪声。要求：（1） 当输入信号u(k)是方差为1的白噪声序列时，利用系统的输入输出值在线辨识上述模型的参数；（2） 当输入信号u(k)是幅值为1的逆m序列时，利用系统的输入输出值在线辨识上述模型的参数；分析比较在不同输入信号作用下，对系统模型参数辨识精度的影响。（1）clear all; close all;a=1 -1.5 0.7 0.1;b=1 2 1.5;d=3

2、; %对象参数na=length(a)-1;nb=length(b)-1; %计算阶次l=500; %数据长度uk=zeros(d+nb,1);yk=zeros(na,1); %输入输出初值u=randn(l,1); %输入采用方差为1的白噪声序列xi=sqrt(0.1)*randn(l,1); % 方差为0.1的白噪声干扰序列theta=a(2:na+1);b; %对象参数真值thetae_1=zeros(na+nb+1,1); %参数初值p=106*eye(na+nb+1);for k=1:l phi=-yk;uk(d:d+nb); %此处phi为列向量 y(k)=phi*theta+xi

3、(k); %采集输出数据 %递推公式 k=p*phi/(1+phi*p*phi); thetae(:,k)=thetae_1+k*(y(k)-phi*thetae_1); p=(eye(na+nb+1)-k*phi)*p; %更新数据 thetae_1=thetae(:,k); for i=d+nb:-1:2 uk(i)=uk(i-1); end uk(1)=u(k); for i=na:-1:2 yk(i)=yk(i-1); end yk(1)=y(k);endplot(1:l,thetae); %line(1:l,theta,theta);xlabel(k);ylabel(参数估计a,b)

4、;legend(a_1,a_2,a_3,b_0,b_1,b_2);axis(0 l -2 2);（2）clear all;a=1 -1.5 0.7 0.1;b=1 2 1.5;d=2; %对象参数na=length(a)-1;nb=length(b)-1; %计算阶次l=20; %数据长度uk=zeros(d+nb,1);yk=zeros(na,1); %输入初值x1=1;x2=1;x3=1;x4=0;s=1;%移位寄存器初值，方波初值xi=rand(l,1);%白噪声序列theta=a(2:na+1);b; %对象参数真值for k=1:l phi(k,:)=-yk;uk(d:d+nb);

5、% phi(k,:)为行向量，便于组成phi矩阵 y(k)=phi(k,:)*theta+xi(k); %采集输出数据 im=xor(s,x4); if im=0 u(k)=-1; else u(k)=1; end s=not(s);m=xor(x3,x4); %产生m序列 %更新数据 x4=x3;x3=x2;x2=x1;x1=m; for i=nb+d:-1:2 uk(i)=uk(i-1); end uk(1)=u(k); for i=na:-1:2 yk(i)=yk(i-1); end yk(1)=y(k);end实验二 最小方差自校正控制实验设二阶纯滞后被控对象的数学模型参数未知或慢时变

6、，仿真实验时用下列模型：式中x(k)为方差为0.1的白噪声。要求：（1） 当设定输入yr(k)为幅值是10的阶跃信号时，设计最小方差直接自校正控制算法对上述对象进行闭环控制；（2） 1）当设定输入yr(k)为幅值是10的方波信号时，设计最小方差直接自校正控制算法对上述对象进行闭环控制；（3） 如果被控对象模型改为： 重复上述（1）、（2）实验，控制结果如何？分析原因。（1）clear all;close all;a=1 -1.5 0.7; b=2.5 1.5; c=1 0.5; d=4; %对象参数na=length(a)-1; nb=length(b)-1; nc=length(c)-1;

7、%计算阶次nh=nb+d-1; ng=na-1; %nh 为多项式h 的阶次,ng 为多项式g 的阶次l=400;uk=zeros(d+nh,1);yk=zeros(d+ng,1);yek=zeros(nc,1); %最优输出预测估计初值yrk=zeros(nc,1);xik=zeros(nc,1);%xiek=zeros(nc,1); %白噪声估计值yr=10*ones(l/4,1);ones(l/4,1);ones(l/4,1);ones(l/4+d,1);% 期望输出xi=sqrt(0.1)*randn(l,1);%方差为0.1 的白噪声序列thetaek=ones(na+nb+d+nc

8、,d);p=106*eye(na+nb+d+nc);for k=1:l time(k)=k; y(k)=-a(2:na+1)*yk(1:na)+b*uk(d:d+nb)+c*xi(k);xik;%采集输出数据 phie=yk(d:d+ng);uk(d:d+nh);-yek(1:nc); k=p*phie/(1+phie*p*phie); thetae(:,k)=thetaek(:,1)+k*(y(k)-phie*thetaek(:,1); p=(eye(na+nb+d+nc)-k*phie)*p; ye=phie*thetaek(:,d);%预测输出估计值 %提取辨识参数 ge=thetae(

9、1:ng+1,k); he=thetae(ng+2:ng+nh+2,k); ce=1 thetae(ng+nh+3:ng+nh+nc+2,k); if abs(ce(2)0.9 ce(2)=sign(ce(2)*0.9; end if he(1)0 yek(1)=ye; yrk(1)=yr(k); xik(1)=xi(k); endendfigure(1);subplot(2,1,1);plot(time,yr(1:l),r:,time,y);xlabel(k);ylabel(y_r(k)、y(k);legend(y_r(k),y(k);axis(0 l -20 20);subplot(2,1

10、,2);plot(time,u);xlabel(k);ylabel(u(k);axis(0 l -40 40);figure(2);subplot(2,1,1);plot(1:l,thetae(1:ng+1,:),1:l,thetae(ng+nh+3:ng+2+nh+nc,:);xlabel(k);ylabel(参数估计g,c);legend(g_0,g_1,c_1);axis(0 l -3 4);subplot(2,1,2);plot(1:l,thetae(ng+2:ng+2+nh,:);xlabel(k);ylabel(参数估计h);legend(h_0,h_1,h_2,h_3,h_4);

11、axis(0 l 0 4);（2）clear all;close all;a=1 -1.5 0.7; b=2.5 1.5; c=1 0.5; d=4; %对象参数na=length(a)-1; nb=length(b)-1; nc=length(c)-1; %计算阶次nh=nb+d-1; ng=na-1; %nh 为多项式h 的阶次,ng 为多项式g 的阶次l=400;uk=zeros(d+nh,1);yk=zeros(d+ng,1);yek=zeros(nc,1); %最优输出预测估计初值yrk=zeros(nc,1);xik=zeros(nc,1);%xiek=zeros(nc,1); %

12、白噪声估计值yr=10*ones(l/4,1);-ones(l/4,1);ones(l/4,1);-ones(l/4+d,1);% 期望输出xi=sqrt(0.1)*randn(l,1);%方差为0.1 的白噪声序列thetaek=zeros(na+nb+d+nc,d);p=106*eye(na+nb+d+nc);for k=1:l time(k)=k; y(k)=-a(2:na+1)*yk(1:na)+b*uk(d:d+nb)+c*xi(k);xik;%采集输出数据 phie=yk(d:d+ng);uk(d:d+nh);-yek(1:nc); k=p*phie/(1+phie*p*phie)

13、; thetae(:,k)=thetaek(:,1)+k*(y(k)-phie*thetaek(:,1); p=(eye(na+nb+d+nc)-k*phie)*p; ye=phie*thetaek(:,d);%预测输出估计值 %提取辨识参数 ge=thetae(1:ng+1,k); he=thetae(ng+2:ng+nh+2,k); ce=1 thetae(ng+nh+3:ng+nh+nc+2,k); if abs(ce(2)0.9 ce(2)=sign(ce(2)*0.9; end if he(1)0 yek(1)=ye; yrk(1)=yr(k); xik(1)=xi(k); ende

14、ndfigure(1);subplot(2,1,1);plot(time,yr(1:l),r:,time,y);xlabel(k);ylabel(y_r(k)、y(k);legend(y_r(k),y(k);axis(0 l -20 20);subplot(2,1,2);plot(time,u);xlabel(k);ylabel(u(k);axis(0 l -40 40);figure(2);subplot(2,1,1);plot(1:l,thetae(1:ng+1,:),1:l,thetae(ng+nh+3:ng+2+nh+nc,:);xlabel(k);ylabel(参数估计g,c);le

15、gend(g_0,g_1,c_1);axis(0 l -3 4);subplot(2,1,2);plot(1:l,thetae(ng+2:ng+2+nh,:);xlabel(k);ylabel(参数估计h);legend(h_0,h_1,h_2,h_3,h_4);axis(0 l 0 4);（3-1）clear all;close all;a=1 -1.5 0.7; b=5 1.5; c=1 0.5; d=4; %对象参数na=length(a)-1; nb=length(b)-1; nc=length(c)-1; %计算阶次nh=nb+d-1; ng=na-1; %nh 为多项式h 的阶次,

16、ng 为多项式g 的阶次l=400;uk=ones(d+nh,1);yk=ones(d+ng,1);yek=ones(nc,1); %最优输出预测估计初值yrk=ones(nc,1);xik=ones(nc,1);%xiek=zeros(nc,1); %白噪声估计值yr=10*ones(l/4,1);ones(l/4,1);ones(l/4,1);ones(l/4+d,1);% 期望输出xi=sqrt(0.1)*randn(l,1);%方差为0.1 的白噪声序列thetaek=ones(na+nb+d+nc,d);p=106*eye(na+nb+d+nc);for k=1:ltime(k)=k

17、;y(k)=-a(2:na+1)*yk(1:na)+b*uk(d:d+nb)+c*xi(k);xik;%采集输出数据phie=yk(d:d+ng);uk(d:d+nh);-yek(1:nc);k=p*phie/(1+phie*p*phie);thetae(:,k)=thetaek(:,1)+k*(y(k)-phie*thetaek(:,1);p=(eye(na+nb+d+nc)-k*phie)*p;ye=phie*thetaek(:,d);%预测输出估计值%提取辨识参数ge=thetae(1:ng+1,k);he=thetae(ng+2:ng+nh+2,k);ce=1 thetae(ng+nh

18、+3:ng+nh+nc+2,k);if abs(ce(2)0.9ce(2)=sign(ce(2)*0.9;endif he(1)0yek(1)=ye;yrk(1)=yr(k);xik(1)=xi(k);endendfigure(1);subplot(2,1,1);plot(time,yr(1:l),r:,time,y);xlabel(k);ylabel(y_r(k)、y(k);legend(y_r(k),y(k);axis(0 l -20 20);subplot(2,1,2);plot(time,u);xlabel(k);ylabel(u(k);axis(0 l -40 40);figure(

19、2);subplot(2,1,1);plot(1:l,thetae(1:ng+1,:),1:l,thetae(ng+nh+3:ng+2+nh+nc,:);xlabel(k);ylabel(参数估计g,c);legend(g_0,g_1,c_1);axis(0 l -3 4);subplot(2,1,2);plot(1:l,thetae(ng+2:ng+2+nh,:);xlabel(k);ylabel(参数估计h);legend(h_0,h_1,h_2,h_3,h_4);axis(0 l 0 4);（3-2）clear all;close all;a=1 -1.5 0.7; b=5 1.5; c

20、=1 0.5; d=4; %对象参数na=length(a)-1; nb=length(b)-1; nc=length(c)-1; %计算阶次nh=nb+d-1; ng=na-1; %nh 为多项式h 的阶次,ng 为多项式g 的阶次l=400;uk=zeros(d+nh,1);yk=zeros(d+ng,1);yek=zeros(nc,1); %最优输出预测估计初值yrk=zeros(nc,1);xik=zeros(nc,1);%xiek=zeros(nc,1); %白噪声估计值yr=10*ones(l/4,1);-ones(l/4,1);ones(l/4,1);-ones(l/4+d,1)

21、;% 期望输出xi=sqrt(0.1)*randn(l,1);%方差为0.1 的白噪声序列thetaek=zeros(na+nb+d+nc,d);p=106*eye(na+nb+d+nc);for k=1:ltime(k)=k;y(k)=-a(2:na+1)*yk(1:na)+b*uk(d:d+nb)+c*xi(k);xik;%采集输出数据phie=yk(d:d+ng);uk(d:d+nh);-yek(1:nc);k=p*phie/(1+phie*p*phie);thetae(:,k)=thetaek(:,1)+k*(y(k)-phie*thetaek(:,1);p=(eye(na+nb+d+

22、nc)-k*phie)*p;ye=phie*thetaek(:,d);%预测输出估计值%提取辨识参数ge=thetae(1:ng+1,k);he=thetae(ng+2:ng+nh+2,k);ce=1 thetae(ng+nh+3:ng+nh+nc+2,k);if abs(ce(2)0.9ce(2)=sign(ce(2)*0.9;endif he(1)0yek(1)=ye;yrk(1)=yr(k);xik(1)=xi(k);endendfigure(1);subplot(2,1,1);plot(time,yr(1:l),r:,time,y);xlabel(k);ylabel(y_r(k)、y(

23、k);legend(y_r(k),y(k);axis(0 l -20 20);subplot(2,1,2);plot(time,u);xlabel(k);ylabel(u(k);axis(0 l -40 40);figure(2);subplot(2,1,1);plot(1:l,thetae(1:ng+1,:),1:l,thetae(ng+nh+3:ng+2+nh+nc,:);xlabel(k);ylabel(参数估计g,c);legend(g_0,g_1,c_1);axis(0 l -3 4);subplot(2,1,2);plot(1:l,thetae(ng+2:ng+2+nh,:);xl

24、abel(k);ylabel(参数估计h);legend(h_0,h_1,h_2,h_3,h_4);axis(0 l 0 4);实验三 模型参考自适应控制实验设被控对象模型参数未知或慢时变，但其状态变量完全可观测，仿真时取状态方程为：选择参考模型：状态完全可观测的模型参考自适应控制系统如下图所示：. (t)e(t)(t)控制器自适应规律为：，式中： ， 为mm矩阵（m为输入个数）。当参考输入为 时，要求选择三组合适的p、r1和r2，实现对被控对象的控制，使被控对象的2个状态变量分别跟踪参考模型的2个状态变量，并分析p、r1和r2对控制系统性能的影响。clear all; close all;h=0.01;l=100/h;as=0 1;-5 -3;bs=0;6;%对象参数am=0 1;-10 -5;bm=0;2;%参考模型参数sz=size(bs);n=sz(1);m=sz(2);%状态向量、输入维数p=3 1;1 1;%正定矩阵r1=7*eye(m);r2=5*eye(m);%自适应律参数矩阵f0=zeros(m,n);g0=zeros(m);%初值yr0=zeros(m,1)