自适应控制程序

1、1 5D0*3I? £axis(0 L/2 -0.5 1.5);xlabel('k');ylabel('M序列幅值');title('M 序列');subplot(2,1,2);% M序列及其逆序列的产生设M序列 M(k)由如下4位移位存放器产生：x = Xj二 Xj /S(k)为方波序列，逆 M 序列IM(k)= M(k)二 S(k)clear all; close all;L=60;%序列长度x仁1;x2=1;x3=1;x4=0;%移位存放器初值S=1;%方波初值for k=1:LIM=xor(S,x4);%进行异或运算，产生逆M序

2、列if IM=0u(k)=-1;elseu(k)=1;endS=not(S);%产生方波M(k)=xor(x3,x4);%产生 M 序列x4=x3;x3=x2;x2=x1;x 仁 M(k);% 存放器移位end subplot(2,1,1); stairs(M);grid;stairs(u);grid;axis(0 L -1.5 1.5);xlabel('k');ylabel(' 逆M序列幅值');title('逆M序列');%白噪声及有色噪声序列的产生设(k)为均值为0,方差为1的高斯白噪声序列，e(k)为有色噪声序列：e(k)二 G(z

3、9;) (k)=Dll ("肿 A 時的高斯白噪声序列(k)在Matlab中由rand()函数产生，程序如下clear all; close all;L=500;%仿真长度%分子分母多项式系数d=1 -1.5 0.7 0.1; c=1 0.5 0.2;nd=le ngth(d)-1 ;n c=le n gth(c)-1;% 阶次xik=zeros( n c,1);% 白噪声初值ek=zeros (n d,1);xi=randn(L,1);%产生均值为0,方差为1的高斯白噪声序列for k=1:L%产生有色噪声%数据更新for i=n d:-1:2ek(i)=ek(i-1);endek

4、(1)=e(k);for i=n c:-1:2xik(i)=xik(i-1);endxik(1)=xi(k);endsubplot(2,1,1);plot(xi);xlabel('k');ylabel('噪声幅值');title('白噪声序列');subplot(2,1,2);plot(e);xlabel('k');ylabel('噪声幅值');title('有色噪声序列');有色厘声序列2:1:0-1J- k60 1DQ 16 : 2QD 2E0300 350A1 : 0 A5 : 5QDe(k)

5、=-d(2: nd+1)*ek+c*xi(k);xik;%批处理最小二乘参数估计(LS)考虑如下系统：y(k) -1.5y(k -1) 0.7y(k -2) = u(k -3) 0.5u(k -4)(k)式中(k)为方差为1的白噪声。clear all;a=1 -1.5 0.7'b=1 0.5'd=3;% 对象参数na=le n gth (a)-1; n b=le n gth(b)-1;% 计算阶次L=500;%数据长度uk=zeros(d+nb,1);yk=zeros (n a,1);% 输入初值x仁1;x2=1;x3=1;x4=0;S=1;%移位存放器初值，方波初值xi=r

6、and(L,1);% 白噪声序列theta=a(2: na+1);b;% 对象参数真值for k=1:Lphi(k,:)=-yk;uk(d:d+nb)： %phi(k,:)为行向量，便于组成phi 矩阵y(k)=phi(k,:)*theta+xi(k); %采集输出数据if IM=0 u(k)=-1;elseu(k)=1;endS=not(S);M=xor(x3,x4); %产生 M 序列%更新数据 x4=x3;x3=x2;x2=x1;x 仁 M;for i=n b+d:-1:2 uk(i)=uk(i-1);end uk(1)=u(k);for i=n a:-1:2 yk(i)=yk(i-1)

7、;endyk(1)=y(k);endthetaevaluatio n=i nv(phi'*phi)*phi'*y'%计算参数估计值 thetaevaluation =-1.53620.68021.00680.4864%遗忘因子递推最小二乘参数估计 (FFRLS )考虑如下系统：y(k) a,y(k -1) a 2y(k - 2) = b °u(k - 3) bu(k -4)(k)k 乞 500k 500T.5,对象时变参数-(kAai,a2,bo,biT为0.7,1,0.5T,式中(k)为均值为0、方差为0.1的白噪声，取遗忘因子 =0.98 , clear

8、all; close all;a=1 -1.5 0.7'b=1 0.5'd=3;% 对象参数na=le ngth (a)-1; nb=le ngth(b)-1;% 计算阶次L=1000;% 数据长度uk=zeros(d+nb,1);yk=zeros (n a,1);% 输入输出初值u=ra ndn (L,1);%输入采用方差为1的白噪声序列xi=sqrt(0.1)*randn(L,1);%方差为0.1的白噪声干扰序列%theta=a(2: na+1);b;% 对象参数真值thetae_ 仁 zeros( na+n b+1,1); % 参数初值P=10 A6*eye( na+nb

9、+1);lambda=0.98;%遗忘因子范围0.9 1for k=1:Lif k=501a=1 -1 0.4'b=1.5 0.2'% 对象参数突变endtheta(:,k)=a(2:n a+1);b;% 对象参数真值phi=-yk;uk(d:d+nb);y(k)=phi'*theta(:,k)+xi(k);% 采样输出数据%遗忘因子递推最小二乘公式K=P*phi/(lambda+phi'*P*phi);thetae(:,k)=thetae_1+K*(y(k)-phi'*thetae_1);P=(eye (n a+n b+1)-K*phi')*P

10、/lambda;%更新数据nrthetae_1=thetae(:,k);for i=d+n b:-1:2uk(i)=uk(i-1);end1 goo woouk(1)=u(k);for i=n a:-1:2yk(i)=yk(i-1);endyk(1)=y(k);endsubplot(2,1,1);plot(1:L,thetae(1:na,:);hold on;plot(1 丄,theta(1:na,:),'k:'); xlabel('k');ylabel('参数估计 a');lege nd('a_1','a_2')

11、;axis(0 L -2 2);8G0IISOO 7lil300400500ksubplot(2,1,2);参数估计 b');plot(1:L,thetae(na+1:na+nb+1,:);hold on;plot(1:L,theta(na+1:na+nb+1,:),'k:'); xlabel('k');ylabel(' legend('b 0','b 1');axis(0 L -0.5 2);%增广递推最小二乘参数估计 (ERLS)考虑如下系统：y(k) -1.5y(k -1) 0.7y(k - 2) = u(k

12、 - 3) 0.5u(k - 4) (k) -(k - 1) 02 (k - 2)式中(k)为方差为0.1的白噪声。选择方差为 1 的白噪声作为输入信号 u(k).clear all; close all;a=1 -1.5 0.7'b=1 0.5'c=1 -1 0.2'd=3; % 对象参数na=le ngth(a)-1; nb=le ngth(b)-1; nc=le ngth(c)-1; % 计算阶次L=1000;%数据长度uk=zeros(d+nb,1);yk=zeros (n a,1);% 输入输出初值xik=zeros( nc,1);% 噪声初值xiek=zer

13、os (n c,1); % 噪声估计初值u=randn (L,1);%输入采用方差为 1 的白噪声序列xi=sqrt(0.1)*randn(L,1); % 方差为 0.1 的白噪声干扰序列theta=a(2: na+1);b;c(2: nc+1); % 对象参数真值thetae_ 仁 zeros(na+nb+1+ nc,1);% 参数初值 ，na+nb+ 1 + nc 为辨识参数个数P=10 A6*eye( na+n b+1+ nc);lambda=0.98;%遗忘因子范围 0.9 1for k=1:Lphi=-yk;uk(d:d+nb);xik;%测量向量y(k)=phi'*thet

14、a+xi(k); % 采样输出数据phie=-yk;uk(d:d+nb);xiek; % 估计的测量向量%增广递推最小二乘公式K=P*phie/(1+phie'*P*phie); thetae(:,k)=thetae_1+K*(y(k)-phie'*thetae_1);P=(eye( na+nb+1+ nc)-K*phie')*P;xie=y(k)-phie'*thetae(:,k); % 白噪声估计值%更新数据thetae_1=thetae(:,k);for i=d+nb:-1:2uk(i)=uk(i-1);enduk(1)=u(k);for i=n a:-1

15、:2yk(i)=yk(i-1);4HD_ _G0DTOD80D oneilD1002D0302400500 fem ;kk1 2考虑如下系统:递推最小二乘参数估计RLSy(k) -1.5y(k -1) 0.7y(k -2) = u(k -3) 0.5u(k -4)式中(k)为方差为0.1的白噪声。(k)end yk(1)=y(k);for i=n c:-1:2xik(i)=xik(i-1);xiek(i)=xiek(i-1);endxik(1)=xi(k);xiek(1)=xie;endfigure(1)plot(1 丄,thetae(1: n a,:); xlabel('k'

16、);ylabel('参数估计a')；lege n d('a_1','a_2');axis(0 L -2 2); figure (2) plot(1:L,thetae( na+1: na+n b+1,:); xlabel('k');ylabel('参数估计 b'); legend('b_0','b_1');axis(0 L 0 1.5); figure(3) plot(1:L,thetae (na+n b+2: na+nb+n c+1,:); xlabel('k');y

17、label('参数估计 c');lege nd('c 1','c 2');axis(0 L -2 2);clear all; close all;a=1 -1.5 0.7'b=1 0.5' d=3;%对象参数na=length(a)-1;nb=length(b)-1;% 计算阶次,na=2,nb=1L=500;%数据长度(仿真长度)uk=zeros(d+nb,1);yk=zeros( na,1);% 输入输出初值 uk:4x1,ykx1u=randn (L,1);%输入采用方差为1的白噪声序列xi=sqrt(O.1)*randn(

18、L,1);%方差为0.1的白噪声干扰序列theta=a(2: na+1);b;% 对象参数真值 theta= -1.5,0.7;1,0.5thetae_仁zeros(na+nb+1,1);%参数初值 日为4x1的全零矩阵P=10 A6*eye( na+nb+1);for k=1:Lphi=-yk;uk(d:d+nb);%此处 phi 为列向量 4x1y(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

19、+n b+1)-K*phi')*P;%更新数据thetae_1=thetae(:,k);for i=d+nb:-1:2uk(i)=uk(i-1);enduk(1)=u(k);for i=n a:-1:2yk(i)=yk(i-1);endyk(1)=y(k);end plot(1:L,thetae); %li n e(1:L,theta,theta);xlabel('k');ylabel('参数估计 a,b');lege n d('a_1','a_2','b_0','b_1');axis(0

20、L -2 2);%最小方差控制 MVC考虑如下系统(k) 0.2 (k-1)y(k) -1.7y(k-1) 0.7y(k - 2) = u(k - 4) 0.5u(k-5) 式中(k)为方差为0.1的白噪声。取期望输出yr(k)为幅值为10的方波信号。clear all;close all;a=1 -1.7 0.7;b=1 0.5;c=1 0.2;d=4;%对象参数na=le ngth(a)-1; nb=le ngth(b)-1; nc=le ngth(c)-1;%计算阶次nh=nb+d-1;%nh为多项式H的阶次L=400;uk=zeros(d+nb,1);yk=zeros (n a,1);

21、yrk=zeros( n c,1);xik=zeros( n c,1);yr=10*o nes(L/4,1);-o nes(L/4,1);o nes(L/4,1);-o nes(L/4+d,1);%期望输出xi=sqrt(0.1)*randn(L,1);%方差为 0.1 的白噪声序列h,f,g=singlediophantine(a,b,c,d);%求解单步 Diophantine 方程for k=1:Ltime(k)=k;y(k)=-a(2: n a+1)*yk+b*uk(d:d+nb)+c*xi(k);xik;%采集输出数据u(k)=(-h(2: nh+1)*uk(1: nh)+c*yr(

22、k+d:-1:k+d-mi n(d, n c);yrk(1: nc-d)-g*y(k);yk(1: na- 1)/h(1);%求控制量%更新数据0for i=d+n b:-1:2 uk(i)=uk(i-1);enduk(1)=u(k);50100150200250300350400for i=n a:-1:2 yk(i)=yk(i-1);endyk(1)=y(k);for i=n c:-1:2yrk(i)=yrk(i-1); xik(i)=xik(i-1);endif n c>0yrk(1)=yr(k);endendsubplot(2,1,1);plot(time,yr(1:L),

23、9;r:',time,y); xlabel('k');ylabel('y_r(k)、 y(k)'); lege nd('y_r(k)','y(k)');subplot(2,1,2);plot(time,u);xlabel('k');ylabel('u(k)');单步 Diophantine 方程求解 求解以下系统的单步 Diophantine 方程：(1)0.7y(k - 2) = u(k - 3) 0.5u(k - 4)(2)2) 2u(k - 3)(3) y(k)-1.7y(k-1) 0

24、.7y(k-2) = 0.9u(k-4) u(k-5)y(k) - 1.5y(k -1)(k)y(k) - 0.95y(k-1) = u(k (k)-0.7 (k-1)(k) 2 (k-1)% 单步 Diophantine 方程的求解clear all;a=1 -1.5 0.7; b=1 0.5; c=1; d=3; % 例 4.1 (1)%a=1 -0.95; b=1 2; c=1 -0.7; d=2; % 例 4.1 (2) %a=1 -1.7 0.7; b=0.9 1; c=1 2; d=4; % 例 4.1 (3) e,f,g=sindiophantine(a,b,c,d) % 调用函

25、数 sindiophantine function e,f,g=singlediophantine(a,b,c,d) % 单步 Diophantine 方程求解 na=length(a)-1;nb=length(b)-1;nc=length(c)-1;% 计算阶次 ne=d-1;ng=na-1;%E,G 的阶次 ad=a,zeros(1,ng+ne+1-na);cd=c,zeros(1,ng+d-nc);% 令 a(na+2)=a(na+3)=.=0e(1)=1;for i=2:n e+1e(i)=0;for j=2:ie(i)=e(i)+e(i+1-j)*ad(j); ende(i)=cd(

26、i)-e(i);% 计算 eiendg(i)=o ；for j=1: n e+1g(i)=g(i)+e( ne+2-1)*ad(i+j);endg(i)=cd(i+d)-g(i);% 计算 gig =1.2750end-1.0850f=conv(b,e);% 计算 F1.00001.50001.55001.00002.00002.30000.7750多步Diophantine 方程求解求解如下系统的多步 Diophantine方程，并去预测长度N=3(k)2 (k -1)y(k) -2.7y(k -1) 2.4y(k - 2) - 0.7y(k - 3) = 0.9u(k -1) u(k -

27、2)%多步Diophantine方程的求解clear all;a=1 -2.7 2.4 -0.7; b=0.9 1; c=1 2;na=length(a)-1; nb=length(b)-1; nc=length(c)-1; %A 、B、C 的阶次N=3; %预测步数E,F,G=multidiophantine(a,b,c,N) %调用函数 multidiophantine function E,F,G=multidiopha ntin、(a bcN) °(%*%功能：多步Diophanine方程的求解%调用格式：E,F,G=si ndiopha ntin e(a,b,c,N)(注：

28、d=1)%输入参数:多项式A、B、C系数向量及预测步数(共 4个)%输出参数:Diophanine方程的解 E、F、G (共3个)(%* na=length(a)-1; nb=length(b)-1; nc=length(c)-1; %A 、B、C 的阶次%E、F、G的初值if na>=ncG(1,:)=c(2:nc+1) zeros(1,na-nc)-a(2:na+1); %令 c(nc+2)=c(nc+3)=.=0 elseG(1,:)=c(2:nc+1)-a(2:na+1) zeros(1,nc-na); %令 a(na+2)=a(na+3)=.=0 end% 求 E、 G、 Ff

29、or j=1:N-1for i=1:jE(j+1,i)=E(j,i);endE(j+1,j+1)=G(j,1);for i=2: naG(j+1,i-1)=G(j,i)-G(j,1)*a(i);endG(j+1, na)=-G(j,1)*a( na+1);F(j+1,:)=co nv(b,E(j+1,:);end%最小方差自校正控制用最小方差自校正控制 算法对以下系统进行闭环控制：y(k) -1.7y(k-1) 0.7y(k - 2) = u(k - 4) 0.5u(k-5)式中(k)为方差为0.1的白噪声。取期望输出yr(k)为幅值为10的方波信号。%最小方差间接自校正控制clear all

30、; 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; %na、 nb、nc 为多项式 A、B、 C 阶次nf=nb+d-1; %nf 为多项式 F 的阶次(k) 0.2 (k-1)L=400; % 控制步数uk=zeros(d+nb,1); % 输入初值： uk(i) 表示 u(k-i); yk=zeros(na,1); % 输出初值yrk=zeros(nc,1); % 期望输出初值xik=zeros(nc,1); % 白噪声初值xiek=ze

31、ros(nc,1); % 白噪声估计初值期望输出 xi=sqrt(0.1)*randn(L,1); % 白噪声序列求控制量yr=10*ones(L/4,1);-ones(L/4,1);ones(L/4,1);-ones(L/4+d,1); %RELS 初值设置thetae_1=0.001*ones(na+nb+1+nc,1); % 非常小的正数 ( 这里不能为 0)P=10 A6*eye( na+n b+1+ nc);for k=1:Ltime(k)=k; y(k)=-a(2:na+1)*yk+b*uk(d:d+nb)+c*xi(k);xik; %采集输出数据% 递推增广最小二乘法 phie=

32、-yk;uk(d:d+nb);xiek;K=P*phie/(1+phie'*P*phie);thetae(:,k)=thetae_1+K*(y(k)-phie'*thetae_1);P=(eye(na+nb+1+nc)-K*phie')*P;xie=y(k)-phie'*thetae(:,k); % 白噪声的估计值%提取辨识参数ae=1 thetae(1:na,k)' be=thetae(na+1:na+nb+1,k)' ce=1 thetae(na+nb+2:na+nb+1+nc,k)' if abs(be(2)>0.9 be(2

33、)=sign(ce(2)*0.9; %MVC 算法要求 B 稳定endif abs(ce(2)>0.9 ce(2)=sign(ce(2)*0.9;ende,f,g=sindiophantine(ae,be,ce,d); % 求解单步 Diophantine 方程 u(k)=(-f(2:nf+1)*uk(1:nf)+ce*yr(k+d:-1:k+d-min(d,nc);yrk(1:nc-d)-g*y(k);yk(1:na-1)/f(1); %更新数据thetae_1=thetae(:,k);for i=d+nb:-1:2uk(i)=uk(i-1);enduk(1)=u(k);yk(i)=y

34、k(i-1);endyk(1)=y(k);for i=n c:-1:2yrk(i)=yrk(i-1);xik(i)=xik(i-1); xiek(i)=xiek(i-1);endif nc>0yrk(1)=yr(k);xik(1)=xi(k);xiek(1)=xie;endendfigure(1);subplot(2,1,1);plot(time,yr(1 丄),'r:',time,y);xlabel('k'); ylabel('y_r(k) 、y(k)'); legend('y_r(k)','y(k)');

35、 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 丄,thetae(1: na,:); xlabel('k'); ylabel('参数估计 a');lege nd('a_1','a_2'); axis(0 L -3 2); subplot(212)plot(1:L,thetae( n a+1: na+n b+

36、1+ n c,:); xlabel('k'); ylabel('参数估计 b、c'); legend('b 0','b 1','c 1'); axis(0 L 0 1.5);4050100150200250设被控对象为:200-20-405b计估数 5y(k)，yr(k)300350400kT 5040D50225030O3502D-1 1-2-35050225030D35D40D.町 mc150T50225030D35D40D%最小方差直接自校正控制y(k) T.7y(kT) 0.7y(k - 2) = u(k

37、- 4) 0.5u(k-5)(k) 02 (kT)式中k为方差为0.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; %na、nb、nc 为多项式 A、 B、C 阶次 nf=nb+d-1;ng=na-1; %nf 、 ng 为多项式 F、 G 的阶次L=400; % 控制步数uk=zeros(d+nf,1); % 输入初值： uk(i) 表示 u(k-i); yk=zeros(d+ng,

38、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 A6*eye( na+n b+d+nc);for k=1:Ltime(k)=k; y(k)=-a(2:na+1)*yk(1:na)+b*

39、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); % 预测输出的估计值 ( 必须为 thetae(:,k-d) ) %ye=yr(k); % 预测输出的估计值可取 yr(k)%提取辨识参数g

40、e=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.9ce(2)=sign(ce(2)*0.9;endif fe(1)<0.1% 设 f0 的下界为 0.1fe(1)=0.1;endu(k)=(-fe(2:n f+1)*uk(1:n f)+ce*yr(k+d:-1:k+d-mina-1)/fe(1);量n(d, n c);yrk(1:n c-d)-ge*y(k);yk(1:%控制20%更新数据k 3for i=d:-1

41、:20、thetaek(:,i)=thetaek(:,i-1);k ry -10end-20 -thetaek(:,1)=thetae(:,k);040for i=d+n f:-1:2 uk(i)=uk(i-1);enduk(1)=u(k);for i=d+ng:-1:2yk(i)=yk(i-1);endyk(1)=y(k);yr(k)- y(k)1300350400200-20-40 (|x|G|pJuffrrI!50100150200250k50100150200250k300350400g0g1for i=n c:-1:2yek(i)=yek(i-1);yrk(i)=yrk(i-1);

42、xik(i)=xik(i-1);endif nc>0yek(1)=ye;yrk(1)=yr(k);xik(1)=xi(k);endend2- g计估数0 2 1 f计估数参50100150200k250300350f0f1f2f3f45023。0350400figure(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(

43、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+nf+3: ng+2+nf+n c,:);xlabel('k'); ylabel(' 参数估计 g、 c');lege nd('g_0','g_1','c_1'); axis(0 L -3 4);subplot(212)plot

44、(1:L,thetae( ng+2: ng+2+nf,:);xlabel('k'); ylabel(' 参数估计 f);lege nd('f_0','f_1','f_2','f_3','f_4'); axis(0 L 0 4);%广义最小方差控制 ( 显示控制 )设被控对象为如下开环不稳定非最小相位系统：(k) 0.2 (k-1)y(k) -1.7y(k -1) 0.7y(k _ 2) = u(k _ 4) 0.5u(k-5)式中(k)为方差为0.1的白噪声。%广义最小方差控制 ( 显式控制

45、律 )clear all; close all;a=1 -1.7 0.7; b=1 2; c=1 0.2; d=4; % 对象参数C 阶次na=length(a)-1; nb=length(b)-1; nc=length(c)-1; %na、 nb、nc 为多项式 A、B、nf=nb+d-1; ng=na-1; %nf 、ng 为多项式 F、G 的阶次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); % 输入初值： uk(i) 表示 u(k-i

46、);yk=zeros(na,1); % 输出初值yrk=zeros(nc,1); % 期望输出初值期望输出xik=zeros(nc,1); % 白噪声初值yr=10*o nes(L/4,1);-o nes(L/4,1);o nes(L/4,1);-o nes(L/4+d,1); %xi=sqrt(0.1)*randn(L,1); %白噪声序列e,f,g=sindiophantine(a,b,c,d); %求解单步 Diophantine 方程CQ=co n v(c,Q); FP=co nv(f,P); CR=co n v(c,R); GP=co n v(g,P); %CQ=C*Qfor k=1

47、:Ltime(k)=k;y(k)=-a(2: n a+1)*yk+b*uk(d:d+nb)+c*xi(k);xik; %采集输出数据u(k)=(-Q(1)*CQ(2: nc+n q+1)*uk(1: nc+ nq)/b(1)-FP(2: np+n f+1)*uk(1: np+nf)+CR*yr(k+d:-1:k+d-mi n(d, nr+n c); yrk(1: nr+n c-d)-GP*y(k); yk(1: np+n g)/(Q(1)*CQ(1)/b +FP(1); %求控制量%更新数据for i=d+n b:-1:2 uk(i)=uk(i-1);enduk(1)=u(k);for i=n

48、 a:-1:2 yk(i)=yk(i-1);endyk(1)=y(k);20100-10-2050100150200250300350yf(k) y(k)400for i=n c:-1:2yrk(i)=yrk(i-1);b o - -fxik(i)=xik(i-1); end if nc>0yrk(1)=yr(k); xik(1)=xi(k);end塔 JJ丄匚0501end subplot(2,1,1);kplot(time,yr(1:L),'r:',time,y); xlabel('k'); ylabel('y_r(k)、y(k)');

49、lege n d('y_r(k)','y(k)');plot(time,u);xlabel('k'); ylabel('u(k)');%广义最小方差自校正控制(间接算法)设被控对象为如下开环不稳定非最小相位系统：(k) 0.2 (k-1)y(k) -1.7y(k -1) 0.7y(k - 2) = u(k - 4) 0.5u(k-5)式中(k)为方差为0.1的白噪声。%广义最小方差自校正控制 间接算法 clear all; close all;a=1 -1.7 0.7; b=1 2; c=1 0.2; d=4; % 对象参数na=

50、length(a)-1; nb=length(b)-1; nc=length(c)-1; %na、 nb、nc 为多项式 A、B、C 阶次nf=nb+d-1; ng=na-1; %nf 、ng 为多项式 F、 G 的阶次Pw=1; R=1; Q=2; % 加权多项式 P、 R、 Qnp=length(Pw)-1; nr=length(R)-1; nq=length(Q)-1;L=400; % 控制步数uk=zeros(d+nb,1); % 输入初值： uk(i) 表示 u(k-i);yk=zeros(na,1); % 输出初值yrk=zeros(nc,1); % 期望输出初值xik=zeros

51、(nc,1); % 白噪声初值xiek=zeros(nc,1); % 白噪声估计初值yr=10*o nes(L/4,1);-o nes(L/4,1);o nes(L/4,1);-o nes(L/4+d,1); % 期望输出xi=sqrt(0.1)*randn(L,1); % 白噪声序列%RELS 初值设置thetae_ 仁 0.001*ones(na+nb+1+ nc,1);% 非常小的正数 ( 此处不能为 0)P=10 A6*eye( na+n b+1+ nc);for k=1:Ltime(k)=k;y(k)=-a(2:na+1)*yk+b*uk(d:d+nb)+c*xi(k);xik; %

52、采集输出数据%递推增广最小二乘法phie=-yk;uk(d:d+nb);xiek;K=P*phie/(1+phie'*P*phie);thetae(:,k)=thetae_1+K*(y(k)-phie'*thetae_1);P=(eye( na+nb+1+ n c)-K*phie')*P;xie=y(k)-phie'*thetae(:,k);%白噪声的估计值%提取辨识参数 ae=1 thetae(1: n a,k)' be=thetae( n a+1: na+n b+1,k)' ce=1 thetae( na+n b+2: na+n b+1+ n

53、 c,k)' if abs(ce (2) )>0.9ce(2)=sig n( ce(2)*0.9;ende,f,g=sindiophantine(ae,be,ce,d); %求解单步 Diophantine 方程CQ=co n v(ce,Q); FP=co nv(f,Pw); CR=co n v(ce,R); GP=co nv(g,Pw); %CQ=Ce*Q u(k)=(-Q(1)*CQ(2: nc+n q+1)*uk(1: nc+n q)/be(1)-FP(2: np+nf+1)*uk(1: np+nf)+CR*yr(k+d:-1:k+d-mi n(d, nr+n c); yr

54、k(1: nr+n c-d)-GP*y(k); yk(1: n p+n g)/(Q(1)*CQ(1)/be(1)+FP(1);%求控制量%更新数据thetae_1=thetae(:,k);for i=d+n b:-1:2 uk(i)=uk(i-1);enduk(1)=u(k);for i=n a:-1:2 yk(i)=yk(i-1);endyk(1)=y(k);for i=n c:-1:2 yrk(i)=yrk(i-1); xik(i)=xik(i-1); xiek(i)=xiek(i-1);end20100-10-201050-55050100150100150200250200k25030

55、0350yr(k)400300350400if nc>0a1yrk(1)=yr(k);xik(1)=xi(k);xiek(1)=xie;endendfigure(1);-0.5 - subplot(2,1,1);plot(time,yr(1:L),'r:',time,y);xlabel('k'); ylabel('y_r(k) 、 y(k)');b识 、a亠 Cl辨数参 b01.50.50legend('y_r(k)','y(k)'); axis(0 L -20 20);-2 csubplot(2,1,2);

56、 0 plot(time,u);xlabel('k'); ylabel('u(k)'); axis(0 L -10 10);figure(2)plot(1:L,thetae);xlabel('k'); ylabel(' 辨识参数 a、 b 、c'); legend('a 1','a 2','b 0','b 1','c 1'); axis(0 L -2 2.5);50b1c?I t100150200250300350400%广义最小方差自校正控制 (直接算法 )设被控对象为如下开环不稳定非最小相位系统：y(k) -1.7y(k-1) 0.7y(k - 2) = u(k - 4) 0.5u(k-5) (k) 0.2 (k-1)式中(k)为