广义最小方差控制

%广义最小方差控制(显示控制)考虑如下系统：()1.7(1)0.7(2)(4)0.5(5)()0.2(1)ykykykukukkk式中(k)为方差为0.1的白噪声。取111()1,()1,()2PzRzQz，期望输出yr(k)为幅值为10的方波信号。clearall;closeall;a=[1-1.70.7];b=[12];c=[10.2];d=4;na=length(a)-1;nb=length(b)-1;nc=length(c)-1;nf=nb+d-1;ng=na-1;P=1;R=1;Q=2;%加权多项式np=length(P)-1;nr=length(R)-1;nq=length(Q)-1;L=400;uk=zeros(d+nb,1);yk=zeros(na,1);yrk=zeros(nc,1);xik=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);[e,f,g]=singlediophantine(a,b,c,d);CQ=conv(c,Q);FP=conv(f,P);CR=conv(c,R);GP=conv(g,P);fork=1:Ltime(k)=k;y(k)=-a(2:na+1)*yk+b*uk(d:d+nb)+c*[xi(k);xik];u1=-Q(1)*CQ(2:nc+nq+1)*uk(1:nc+nq)/b(1)-FP(2:np+nf+1)*uk(1:np+nf);u2=CR*[yr(k+d:-1:k+d-min(d,nr+nc));yrk(1:nr+nc-d)];u(k)=(u1+u2-GP*[y(k);yk(1:np+ng)])/(Q(1)*CQ(1)/b(1)+FP(1));%更新数据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);fori=nc:-1:2yrk(i)=yrk(i-1);xik(i)=xik(i-1);endifnc0yrk(1)=yr(k);xik(1)=xi(k);endendsubplot(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)');subplot(2,1,2);plot(time,u);xlabel('k');ylabel('u(k)');