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

中南大学系统辨识及自适应控制实验指导老师贺建军姓名专业班级测控1102班0909111814号实验日期2014年11月实验一递推二乘法参数辨识设被辨识系统的数学模型由下式描述：2341231232.01.51()()()11.50.70.111.50.70.1zzzykukkzzzzzz式中(k)为方差为0.1的白噪声。要求：（1）当输入信号u(k)是方差为1的白噪声序列时，利用系统的输入输出值在线辨识上述模型的参数；（2）当输入信号u(k)是幅值为1的逆M序列时，利用系统的输入输出值在线辨识上述模型的参数；分析比较在不同输入信号作用下，对系统模型参数辨识精度的影响。（1）clearall;closeall;a=[1-1.50.70.1]';b=[121.5]';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(0.1)*randn(L,1);%方差为0.1的白噪声干扰序列theta=[a(2:na+1);b];%对象参数真值thetae_1=zeros(na+nb+1,1);%参数初值P=10^6*eye(na+nb+1);fork=1:Lphi=[-yk;uk(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);fori=d+nb:-1:2uk(i)=uk(i-1);enduk(1)=u(k);fori=na:-1:2yk(i)=yk(i-1);endyk(1)=y(k);endplot([1:L],thetae);%line([1:L],[theta,theta]);xlabel('k');ylabel('参数估计a,b');legend('a_1','a_2','a_3','b_0','b_1','b_2');axis([0L-22]);（2）clearall;a=[1-1.50.70.1]';b=[121.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];%对象参数真值fork=1:Lphi(k,:)=[-yk;uk(d:d+nb)]';%phi(k,:)为行向量，便于组成phi矩阵y(k)=phi(k,:)*theta+xi(k);%采集输出数据IM=xor(S,x4);ifIM==0u(k)=-1;elseu(k)=1;endS=not(S);M=xor(x3,x4);%产生M序列%更新数据x4=x3;x3=x2;x2=x1;x1=M;fori=nb+d:-1:2uk(i)=uk(i-1);enduk(1)=u(k);fori=na:-1:2yk(i)=yk(i-1);endyk(1)=y(k);End实验二最小方差自校正控制实验设二阶纯滞后被控对象的数学模型参数未知或慢时变，仿真实验时用下列模型：34112122.51.510.5()()()11.50.711.50.7zzzykukkzzzz式中(k)为方差为0.1的白噪声。要求：（1）当设定输入yr(k)为幅值是10的阶跃信号时，设计最小方差直接自校正控制算法对上述对象进行闭环控制；（2）1）当设定输入yr(k)为幅值是10的方波信号时，设计最小方差直接自校正控制算法对上述对象进行闭环控制；（3）如果被控对象模型改为：34112120.51.510.5()()()11.50.711.50.7zzzykukkzzzz重复上述（1）、（2）实验，控制结果如何？分析原因。（1）clearall;closeall;a=[1-1.50.7];b=[2.51.5];c=[10.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)];%期望输出xi=sqrt(0.1)*randn(L,1);%方差为0.1的白噪声序列thetaek=ones(na+nb+d+nc,d);P=10^6*eye(na+nb+d+nc);fork=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+nc)-K*phie')*P;ye=phie'*thetaek(:,d);%预测输出估计值%提取辨识参数ge=thetae(1:ng+1,k)';he=thetae(ng+2:ng+nh+2,k)';ce=[1thetae(ng+nh+3:ng+nh+nc+2,k)'];ifabs(ce(2))0.9ce(2)=sign(ce(2))*0.9;endifhe(1)0.1he(1)=0.1;%设h0的下界为0.1endu(k)=(-he(2:nh+1)*uk(1:nh)+ce*[yr(k+d:-1:k+d-min(d,nc));yrk(1:nc-d)]-ge*[y(k);yk(1:na-1)])/he(1);%求控制量%更新数据fori=d:-1;2thetaek(:,i)=thetaek(:,i-1);endthetaek(:,1)=thetae(:,k);fori=d+nh:-1:2uk(i)=uk(i-1);enduk(1)=u(k);fori=d+ng:-1:2yk(i)=yk(i-1);endyk(1)=y(k);fori=nc:-1:2yek(i)=yek(i-1);yrk(i)=yrk(i-1);xik(i)=xik(i-1);endifnc0yek(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([0L-2020]);subplot(2,1,2);plot(time,u);xlabel('k');ylabel('u(k)');axis([0L-4040]);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([0L-34]);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([0L04]);（2）clearall;closeall;a=[1-1.50.7];b=[2.51.5];c=[10.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)];%期望输出xi=sqrt(0.1)*randn(L,1);%方差为0.1的白噪声序列thetaek=zeros(na+nb+d+nc,d);P=10^6*eye(na+nb+d+nc);fork=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+nc)-K*phie')*P;ye=phie'*thetaek(:,d);%预测输出估计值%提取辨识参数ge=thetae(1:ng+1,k)';he=thetae(ng+2:ng+nh+2,k)';ce=[1thetae(ng+nh+3:ng+nh+nc+2,k)'];ifabs(ce(2))0.9ce(2)=sign(ce(2))*0.9;endifhe(1)0.1he(1)=0.1;%设h0的下界为0.1endu(k)=(-he(2:nh+1)*uk(1:nh)+ce*[yr(k+d:-1:k+d-min(d,nc));yrk(1:nc-d)]-ge*[y(k);yk(1:na-1)])/he(1);%求控制量%更新数据fori=d:-1;2thetaek(:,i)=thetaek(:,i-1);endthetaek(:,1)=thetae(:,k);fori=d+nh:-1:2uk(i)=uk(i-1);enduk(1)=u(k);fori=d+ng:-1:2yk(i)=yk(i-1);endyk(1)=y(k);fori=nc:-1:2yek(i)=yek(i-1);yrk(i)=yrk(i-1);xik(i)=xik(i-1);endifnc0yek(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([0L-2020]);subplot(2,1,2);plot(time,u);xlabel('k');ylabel('u(k)');axis([0L-4040]);figure(2);subplot(2,1,1);p