系统辨识及自适应控制实验报告.pdf

《系统辨识及自适应控制实验报告.pdf》由会员分享，可在线阅读，更多相关《系统辨识及自适应控制实验报告.pdf（30页珍藏版）》请在淘文阁 - 分享文档赚钱的网站上搜索。

1、 中南大学 系统辨识及自适应控制实验 指导老师 贺建军 姓 名 专业班级 测控 1102 班 04 号 实验日期 2014 年 11 月￥实验一 递推二乘法参数辨识 设被辨识系统的数学模型由下式描述：2341231232.01.51()()()1 1.50.70.11 1.50.70.1zzzy ku kkzzzzzz 式中(k)为方差为的白噪声。要求：（1）当输入信号u(k)是方差为 1 的白噪声序列时，利用系统的输入输出值在线辨识上述模型的参数；（2）当输入信号u(k)是幅值为 1 的逆 M 序列时，利用系统的输入输出值在线辨识上述模型的参数；分析比较在不同输入信号作用下，对系统模型参数辨

2、识精度的影响。（1）clear all;close all;a=1 ;b=1 2;d=3;%对象参数 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*randn(L,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(

3、d:d+nb);%此处 phi 为列向量 y(k)=phi*theta+xi(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);end plot(1:L,thetae);%line(1:L,theta,the

4、ta);xlabel(k);ylabel(参数估计 a,b);legend(a_1,a_2,a_3,b_0,b_1,b_2);axis(0 L-2 2);（2）clear all;a=1 ;b=1 2;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,:)=-

5、yk;uk(d:d+nb);%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、型参数未知或慢时变，仿真实验时用下列模型：34112122.51.510.5()()()1 1.50.71 1.50.7zzzy ku kkzzzz 式中(k)为方差为的白噪声。要求：（1）当设定输入 yr(k)为幅值是 10 的阶跃信号时，设计最小方差直接自校正控制算法对上述对象进行闭环控制；（2）1）当设定输入 yr(k)为幅值是 10 的方波信号时，设计最小方差直接自校正控制算法对上述对象进行闭环控制；（3）如果被控对象模型改为：、34112120.51.510.5()()()1 1.50.71 1.50.7zzzy ku kkzzzz 重复上述（1）、（2）实验，控制结果如何分析原因。

7、（1）clear all;close all;a=1 ;b=;c=1;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);

8、ones(L/4,1);ones(L/4+d,1);%期望输出 xi=sqrt*randn(L,1);%方差为 的白噪声序列 thetaek=ones(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);：thetae(:,k)=thetaek(:,1)+K*(y(k)-phie*thetaek(

9、:,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)ce(2)=sign(ce(2)*;end if he(1)0 yek(1)=ye;yrk(1)=yr(k);xik(1)=xi(k);end end figure(1);subplot(2,1,1);plot(time,yr(1:L),r:,time,y);xlabel(

10、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);legend(g_0,g_1,c_1);axis(0 L-3 4);subplot(2,1,2);plot(1:L,thetae(

11、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);（2）clear all;close all;a=1 ;b=;c=1;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=ze

12、ros(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*randn(L,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+n

13、h);-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+3:ng+nh+nc+2,k);if abs(ce(2)ce(2)=sign(ce(2)*;end if he(1)0|yek(1)=ye;yrk(1)=yr(

14、k);xik(1)=xi(k);end end figure(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(

15、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-1）clear all;close all;a=1 ;b=5;c=1;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*randn(L,1);%方差为 的白噪声序列 thetaek=ones(na+nb+d+nc,d);P=106*eye(na+nb+d+nc);for k=1:L time(k)=k;y(

17、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);ce(2)=sign(ce(2)*;end）if he(1)0 yek(1)=ye;yrk(1)=yr(k);xik(1)=xi(k);end end figure(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

18、,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)

19、;axis(0 L 0 4);（3-2）clear all;close all;a=1 ;b=5;c=1;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/

20、4,1);-ones(L/4,1);ones(L/4,1);-ones(L/4+d,1);%期望输出 xi=sqrt*randn(L,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);thetae(:,k)=thetaek(:,1)+K

21、*(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)ce(2)=sign(ce(2)*;end?if he(1)0 yek(1)=ye;yrk(1)=yr(k);xik(1)=xi(k);end end figure(1);subplot(2,1,1);plot(time,yr(1:L

22、),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);legend(g_0,g_1,c_1);axis(0 L-3 4);subplot(2,1,

23、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);实验三 模型参考自适应控制实验 设被控对象模型参数未知或慢时变，但其状态变量完全可观测，仿真 时取状态方程为：010536ssXXu 选择参考模型：0101052mmrXXy 状态完全可观测的模型参考自适应控制系统如下图所示：控制器自适应规律为：10(,)()()(0)tTTmsF e tR B PeXdF，、20(,)()()(0)tTTmG e tR B PeudG 式中：*1122,T

24、TRGRG ，为mm矩阵（m为输入个数）。当参考输入为()4sin(0.2)ry tt时，要求选择三组合适的 P、R1.!图 4.2-2 利用状态反馈的 MRAC 系统框图 cG()()(ssssX tA X tB u()()()mmmmXtA XtB u t F X*11mGB P*12mGB P ()sXe()mXe()u eX()u t()e e ()ry()u t()u t()ry t(e(和 R2，实现对被控对象的控制，使被控对象的 2 个状态变量分别跟踪参考模型的 2 个状态变量，并分析 P、R1和 R2对控制系统性能的影响。clear all;close all;h=;L=100

25、/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);u0=zeros(m,1);e0=zeros(n,1);xs0=zeros(n,1);xm0=zeros(n,1);for k=1:L time(k)=k*h;yr(k)=4*sin*pi*time(k);%输入信号 x

26、s(:,k)=xs0+h*(As*xs0+Bs*u0);%对象状态 xm(:,k)=xm0+h*(Am*xm0+Bm*yr0);%参考模型状态 e(:,k)=xm(:,k)-xs(:,k);%状态误差 F=F0+h*(R1*Bm*P*e0*xs0);%自适应律 G=G0+h*(R2*Bm*P*e0*yr0);u(:,k)=G*yr(k)+F*xs(:,k);%控制量 parae(1,k)=norm(Am-As-Bs*F);%参数收敛于 Am 的偏差 parae(2,k)=norm(Bm-Bs*G);%参数收敛于 Bm 的偏差%更新数据 yr0=yr(:,k);u0=u(:,k);e0=e(:,

27、k);xs0=xs(:,k);xm0=xm(:,k);F0=F;G0=G;end figure(1)subplot(2,1,1);plot(time,xm(1,:),r,time,xs(1,:),:);xlabel(t);ylabel(x_m_1(t),x_s_1(t);legend(x_m_1(t),x_s_1(t);subplot(2,1,2);plot(time,xm(2,:),r,time,xs(2,:),:);xlabel(t);ylabel(x_m_2(t),x_s_2(t);legend(x_m_2(t),x_s_2(t);figure(2);plot(time,parae(1,:),:,time,parae(2,:),r);xlabel(t);ylabel(参数收敛偏差 E);legend(|A_m-A_s-B_s*F|_2,|B_m-B_s*G|_2);