广义最小方差自校正控制（直接算法）.docVIP

%广义最小方差自校正控制(直接算法) 式中((k)为方差为0.1的白噪声。 取，期望输出yr(k)为幅值为10的方波信号。 clear all;close all; a=[1 -1.7 0.7];b=[1 2];c=[1 0.2];d=4; na=length(a)-1;nb=length(b)-1;nc=length(c)-1; nf=nb+d-1;ng=na-1; Pw=1;R=1;Q=2; %加权多项式 np=length(Pw)-1;nr=length(R)-1;nq=length(Q)-1; L=400; uk=zeros(d+nf,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); %RELS初值设定 thetaek=zeros(na+nb+d+nc,d); P=10^6*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+nf);-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); %最优预测输出估计值 % xie=y(k)-phie*thetae(:,k);%白噪声估计值 %提取辨识参数 ge=thetae(1:ng+1,k); fe=thetae(ng+2:ng+nf+2,k); ce=[1 thetae(ng+nf+3:ng+nf+2+nc,k)]; if abs(ce(2))0.9 ce(2)=sign(ce(2))*0.9; end if fe(1)0.1 fe(1)=0.1; end %[e,f,g]=singlediophantine(ae,be,ce,d); CQ=conv(ce,Q);FP=conv(fe,Pw);CR=conv(ce,R);GP=conv(ge,Pw); u1=-Q(1)*CQ(2:nc+nq+1)*uk(1:nc+nq)/fe(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)*Q(1)/fe(1)+fe(1)); %更新数据 for i=d:-1:2 thetaek(:,i)=thetaek(:,i-1); end thetaek(:,1)=thetae(:,k); for i=d+nf:-1:2 uk(i)=uk(i-1); end uk(1)=u(k); for i=d+ng:-1:2 yk(i)=yk(i-1); end yk(1)=y(k); for i=nc:-1:2 yek(i)=yek(i-1); yrk(i)=yrk(i-1); xik(i)=xik(i-1); end if nc0 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,2); plot(time,u); xlabel(k);ylabel(u(k));axis([0 L -10 10]); figure(2)