用最小方差自校正控制算法对以下系统进行闭环控制： ( ) 1.7 ( y k y k  1) 0.7 ( y k    2)  ( u k  4) 0.5 ( u k   5)  ( ) 0.2 ( k k     1) 式中(k)为方差为 0.1 的白噪声, 取期望输出 yr(k)为幅值为 10 的方波 信号。 解：上式可以化为： (1 1.7  z 1   0.7 z 2  ) ( ) y k  z 4  (1 0.5  z 1  ) ( ) u k   (1 0.2 z 1  ) ( ) k 则有 ( A z ( B z ( C z d  1  1  1  4 1   0.7 z 2  ) 1 1.7 z   ) 1 0.5 z   ) 1 0.2 z   1  1  Diophantine 方程为： ( C z  1 )   ( ( A z F z ) 1  1 )  d (  z G z  1 ) 又有 取初值 6 P(0)=10 ,   (0) 0  ( H z  1 )   ( ( B z F z ) 1  1 ) 递推公式为： ( ) k ˆ         ˆ (  k  ( ) K k ( ) P k ˆ ) (  k  1)] 1)   ˆ T ( )[ ( ) ( k d K k y k    ˆ ) ( 1) ( k d P k    ˆ ˆ T 1) ( ) ( ( k d P k      ˆ T ( ) )] ( ( I K k k d P k    ) k d  1)   1 [  程序清单如下： clear all; close all; a=[1 -1.7 0.7]; b=[1 0.5]; c=[1 0.2]; d=4; %对象参数 

na=length(a)-1; nb=length(b)-1; nc=length(c)-1; %多项式 A、B、C 的阶次 nh=nb+d-1; ng=na-1; %nh、ng 为多项式 H、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); %白噪声初值 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); %白噪声序列 %递推估计初值 thetaek=zeros(na+nb+d+nc,d); P=10^6*eye(na+nb+d+nc); %P(k)的初始值 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(:,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+2+nc,k)']; 

if abs(ce(2))>0.9 ce(2)=sign(ce(2))*0.9; end if he(1)<0.1 %设 h0 的下界为 0.1 he(1)=0.1; end u(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); %控制量 %更新数据 for i=d:-1:2 thetaek(:,i)=thetaek(:,i-1); end thetaek(:,1)=thetae(:,k); for i=d+nh:-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 nc>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,2); plot(time,u); xlabel('k'); ylabel('u(k)'); axis([0 L -40 40]); figure(2) subplot(211) 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(212) 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]); 

) k ( y 、 ) k ( r y ) k ( u 20 10 0 -10 -20 0 40 20 0 -20 -40 0 c 、 g 计 估 数 参 h 计 估 数 参 4 2 0 -2 0 4 3 2 1 0 0 yr(k) y(k) 50 100 150 200 k 250 300 350 400 50 100 150 200 k 250 300 350 400 g0 g1 c1 50 100 150 200 k 250 300 350 400 h0 h1 h2 h3 h4 50 100 150 200 k 250 300 350 400