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

本文格式为Word版，下载可任意编辑——系统辨识及自适应控制试验报告中南大学

系统辨识及自适应控

制试验

指导老师贺建军姓名

专业班级测控1102班0909111814号试验日期2023年11月

试验一递推二乘法参数辨识

设被辨识系统的数学模型由下式描述：

z?2?2.0z?3?1.5z?41y(k)?u(k)??(k)?1?2?3?1?2?31?1.5z?0.7z?0.1z1?1.5z?0.7z?0.1z式中?(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:L

phi=[-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);end

plot([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:L

phi(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;end

S=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

试验二最小方差自校正控制试验

设二阶纯滞后被控对象的数学模型参数未知或慢时变，仿真试验时用以下模型：

2.5z?3?1.5z?41?0.5z?1y(k)?u(k)??(k)?1?2?1?21?1.5z?0.7z1?1.5z?0.7z式中?(k)为方差为0.1的白噪声。要求：

（1）当设定输入yr(k)为幅值是10的阶跃信号时，设计最小方差直接

自校正控制算法对上述对象进行闭环控制；

（2）1）当设定输入yr(k)为幅值是10的方波信号时，设计最小方差

直接自校正控制算法对上述对象进行闭环控制；（3）假使被控对象模型改为：

0.5z?3?1.5z?41?0.5z?1u(k)??(k)y(k)??1?2?1?21?1.5z?0.7z1?1.5z?0.7z重复上述（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.9

ce(2)=sign(ce(2))*0.9;endifhe(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([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.9

ce(2)=sign(ce(2))*0.9;endifhe(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([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]);

（3-1）clearall;closeall;

a=[1-1.50.7];b=[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=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=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)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([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]);

（3-2）clearall;closeall;

a=[1-1.50.7];b=[51.5];c=[10.5];d=4;%对象参数na=length(a)-1;nb=length(b)-1;nc=length(c)-1;%计算阶次

.(t)?AX(t)?Buyrt(tX))mmmm(y((t)u)r)Xm((et)?ee((et))Gc??u))u('(te.(t)?AX(t)?Bu(t)Xssss'(t)F?t))Xs((e??

??1G*?1BmPX?X???2G*?1BmP图4.2-2利用状态反馈的MRAC系统框图

控制器自适应规律为：

tTF(e,t)???R1BmPe(?)XsT(?)d??F(0)，

0TG(e,t)??R2BmPe(?)uT(?)d??G(0)

0t式中：R1??1G*,R2??2G*?T?T，为m?m矩阵（m为输入个数）。

当参考输入为yr(t)?4sin(0.2?t)时，要求选择三组适合的P、R1和R2，实现对被控对象的控制，使被控对象的2个状态变量分别跟踪参考模型的2个状态变量，并分析P、R1和R2对控制系统性能的影响。

clearall;closeall;h=0.01;L=100/h;

As=[01;-5-3];Bs=[0;6];%对象参数Am=[01;-10-5];Bm=[0;2];%参考模型参数

Sz=size(Bs);n=Sz(1);m=Sz(2);%状态向量、输入维数P=[31;11];%正定矩阵

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);fork=1:Ltime(k)=k*h;

yr(k)=4*sin(0.2*pi*time(k));%输入信号

xs(:,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(:,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');

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+p