自适应控制程序.pdf

《自适应控制程序.pdf》由会员分享，可在线阅读，更多相关《自适应控制程序.pdf（23页珍藏版）》请在淘文阁 - 分享文档赚钱的网站上搜索。

1、% M% M 序列及其逆序列的产生序列及其逆序列的产生设 M 序列M(k)由如下 4 位移位寄存器产生：xi xi3 xi4S(k)为方波序列，逆 M 序列IM(k)= M(k)S(k)clear all; close all;L=60;%序列长度x1=1;x2=1;x3=1;x4=0;%移位寄存器初值S=1;%方波初值for k=1:LIM=xor(S,x4);%进行异或运算，产生逆M 序列if IM=0u(k)=-1;elseu(k)=1;endS=not(S);%产生方波M(k)=xor(x3,x4);%产生 M 序列x4=x3;x3=x2;x2=x1;x1=M(k);%寄存器移位end

2、subplot(2,1,1);stairs(M);grid;axis(0 L/2 -0.5 1.5);xlabel(k);ylabel(M序列幅值);title(M 序列);subplot(2,1,2);stairs(u);grid;axis(0 L -1.5 1.5);xlabel(k);ylabel(逆 M 序列幅值);title(逆 M 序列);%白噪声及有色噪声序列的产生白噪声及有色噪声序列的产生设(k) 为均值为 0，方差为 1 的高斯白噪声序列，e(k)为有色噪声序列：C(z1)10.5z10.2z2(k) (k)e(k) G(z )(k) 1123D(z )11.5z0.7z0.

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

4、);endek(1)=e(k);for i=nc:-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(有色噪声序列);%批处理最小二乘参数估计批处理最小二乘参数估计(LS)(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

5、.7;b=1 0.5;d=3;%对象参数na=length(a)-1;nb=length(b)-1;%计算阶次L=500;%数据长度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;%对象参数真值for k=1:Lphi(k,:)=-yk;uk(d:d+nb); %phi(k,:)为行向量，便于组成 phi 矩阵y(k)=phi(k,:)*theta+xi(k);%采集输出数据IM=xor(S,x4);if IM=0u(k)

6、=-1;elseu(k)=1;endS=not(S);M=xor(x3,x4);%产生 M 序列%更新数据x4=x3;x3=x2;x2=x1;x1=M;for i=nb+d:-1:2uk(i)=uk(i-1);enduk(1)=u(k);for i=na:-1:2yk(i)=yk(i-1);endyk(1)=y(k);endthetaevaluation=inv(phi*phi)*phi*y%计算参数估计值thetaevaluation =thetaevaluation =-1.5362-1.53620.68020.68021.00681.00680.48640.4864%遗忘因子递推最小二乘

7、参数估计遗忘因子递推最小二乘参数估计(FFRLS)(FFRLS)考虑如下系统：y(k)a1y(k 1)a2y(k 2) b0u(k 3)bu1(k 4)(k)式中(k)为均值为 0、方差为 0.1 的白噪声，T1.5,0.7,1,0.5 ,k 500T对象时变参数(k) a1,a2,b0,b1为：(k) T1,0.4,1.5, 0.2 ,k 500取遗忘因子=0.98，clear all; close all;a=1 -1.5 0.7;b=1 0.5;d=3;%对象参数na=length(a)-1;nb=length(b)-1;%计算阶次L=1000;%数据长度uk=zeros(d+nb,1)

8、;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=106*eye(na+nb+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:na+1);b;%对象参数真值phi=-yk;uk(d:d+nb);y(k)=phi

9、*theta(:,k)+xi(k);%采样输出数据%遗忘因子递推最小二乘公式K=P*phi/(lambda+phi*P*phi);thetae(:,k)=thetae_1+K*(y(k)-phi*thetae_1);P=(eye(na+nb+1)-K*phi)*P/lambda;%更新数据thetae_1=thetae(:,k);for i=d+nb:-1:2uk(i)=uk(i-1);enduk(1)=u(k);for i=na:-1:2yk(i)=yk(i-1);endyk(1)=y(k);endsubplot(2,1,1);plot(1:L,thetae(1:na,:);hold on;

10、plot(1:L,theta(1:na,:),k:);xlabel(k);ylabel(参数估计 a);legend(a_1,a_2);axis(0 L -2 2);subplot(2,1,2);plot(1:L,thetae(na+1:na+nb+1,:);hold on;plot(1:L,theta(na+1:na+nb+1,:),k:);xlabel(k);ylabel(参数估计 b);legend(b_0,b_1);axis(0 L -0.5 2);%增广递推最小二乘参数估计增广递推最小二乘参数估计(ERLS)(ERLS)考虑如下系统：y(k)1.5y(k 1)0.7y(k 2) u(

11、k 3)0.5u(k 4)(k)(k 1)0.2(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=length(a)-1;nb=length(b)-1;nc=length(c)-1;%计算阶次L=1000;%数据长度uk=zeros(d+nb,1);yk=zeros(na,1);%输入输出初值xik=zeros(nc,1);%噪声初值xiek=zeros(nc,1);%噪声估计初值u=randn(L,1);%输入采用方差为

12、 1 的白噪声序列xi=sqrt(0.1)*randn(L,1);% 方差为 0.1 的白噪声干扰序列theta=a(2:na+1);b;c(2:nc+1);%对象参数真值thetae_1=zeros(na+nb+1+nc,1);%参数初值,na+nb+1+nc 为辨识参数个数P=106*eye(na+nb+1+nc);lambda=0.98;%遗忘因子范围0.9 1for k=1:Lphi=-yk;uk(d:d+nb);xik;%测量向量y(k)=phi*theta+xi(k);%采样输出数据phie=-yk;uk(d:d+nb);xiek;%估计的测量向量%增广递推最小二乘公式K=P*ph

13、ie/(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=na:-1:2yk(i)=yk(i-1);endyk(1)=y(k);for i=nc:-1:2xik(i)=xik(i-1);xiek(i)=xiek(i-1);endxik(1)=xi(k);xi

14、ek(1)=xie;endfigure(1)plot(1:L,thetae(1:na,:);xlabel(k);ylabel(参数估计 a);legend(a_1,a_2);axis(0 L -2 2);figure(2)plot(1:L,thetae(na+1:na+nb+1,:);xlabel(k);ylabel(参数估计 b);legend(b_0,b_1);axis(0 L 0 1.5);figure(3)plot(1:L,thetae(na+nb+2:na+nb+nc+1,:);xlabel(k);ylabel(参数估计 c);legend(c_1,c_2);axis(0 L -2

15、2);递推最小二乘参数估计递推最小二乘参数估计(RLS)(RLS)考虑如下系统：y(k)1.5y(k 1)0.7y(k 2) u(k 3)0.5u(k 4)(k)式中(k)为方差为 0.1 的白噪声。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=sqr

16、t(0.1)*randn(L,1);%方差为 0.1 的白噪声干扰序列theta=a(2:na+1);b;%对象参数真值 theta=-1.5,0.7;1,0.5thetae_1=zeros(na+nb+1,1);%参数初值 为 4x1 的全零矩阵P=106*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+nb+1)

17、-K*phi)*P;%更新数据thetae_1=thetae(:,k);for i=d+nb:-1:2uk(i)=uk(i-1);enduk(1)=u(k);for i=na:-1:2yk(i)=yk(i-1);endyk(1)=y(k);endplot(1:L,thetae);%line(1:L,theta,theta);xlabel(k);ylabel(参数估计 a,b);legend(a_1,a_2,b_0,b_1);axis(0 L -2 2);%最小方差控制最小方差控制(MVC)(MVC)考虑如下系统：y(k)1.7y(k 1)0.7y(k 2) u(k 4)0.5u(k 5)(k)

18、0.2(k 1)式中(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=length(a)-1;nb=length(b)-1;nc=length(c)-1;%计算阶次nh=nb+d-1;%nh 为多项式 H 的阶次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);-on

19、es(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:na+1)*yk+b*uk(d:d+nb)+c*xi(k);xik;%采集输出数据u(k)=(-h(2:nh+1)*uk(1:nh)+c*yr(k+d:-1:k+d-min(d,nc);yrk(1:nc-d)-g*y(k);yk(1:na-1)/h(1);% 求控制量%更新数据for i=d+nb:-1:2uk(i)=uk

20、(i-1);enduk(1)=u(k);for i=na:-1:2yk(i)=yk(i-1);endyk(1)=y(k);for i=nc:-1:2yrk(i)=yrk(i-1);xik(i)=xik(i-1);endif nc0yrk(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);单步单步 Diopha

21、ntineDiophantine 方程求解方程求解求解下列系统的单步 Diophantine 方程：（1）（2）%单步 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) %调用函数 sindiophantinefunction e,f,g=singlediophantine(

22、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:ne+1e(i)=0;for j=2:ie(i)=e(i)+e(i+1-j)*ad(j);ende(i)=cd(i)-e(i);%计算 eiendfor i=1:ng+1y(k)1.5y(k 1)0.7y(k 2)u(k 3)0.

23、5u(k 4)(k)y(k)0.95y(k 1) u(k 2) 2u(k 3)(k)0.7(k 1)（3）y(k)1.7y(k 1)0.7y(k 2) 0.9u(k 4)u(k 5)(k)2(k 1)g(i)=0;for j=1:ne+1g(i)=g(i)+e(ne+2-1)*ad(i+j);endg(i)=cd(i+d)-g(i);%计算 giendf=conv(b,e);%计算 Fe =1.00001.50001.5500f =1.00002.00002.30000.7750g =1.2750-1.0850多步多步 DiophantineDiophantine 方程求解方程求解求解如下系统

24、的多步 Diophantine 方程，并去预测长度N=3y(k)2.7y(k 1) 2.4y(k 2)0.7y(k 3) 0.9u(k 1)u(k 2)(k) 2(k 1)%多步 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) %调用函数 multidiophantinefunction E,F,G=multidiophan

25、tine(a,b,c,N)%*%功能：多步 Diophanine 方程的求解%调用格式：E,F,G=sindiophantine(a,b,c,N)（注：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 的初值E=zeros(N); E(1,1)=1; F(1,:)=conv(b,E(1,:);if na=ncG(1,:)=c(2:nc+1) zeros(1,na-nc)-a(

26、2:na+1); %令 c(nc+2)=c(nc+3)=.=0elseG(1,:)=c(2:nc+1)-a(2:na+1) zeros(1,nc-na); %令 a(na+2)=a(na+3)=.=0end%求 E、G、Ffor 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,:)=conv(b,E(j+1,:);end%最小方差自校正控制最小方差自校正控制用最小方差自校正控制最小方差自校

27、正控制算法对以下系统进行闭环控制：y(k)1.7y(k 1)0.7y(k 2) u(k 4)0.5u(k 5)(k)0.2(k 1)式中(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=length(a)-1; nb=length(b)-1; nc=length(c)-1; %na、nb、nc 为多项式 A、B、C 阶次nf=nb+d-1; %nf 为多项式 F 的阶次L=400; %控制步数uk=zero

28、s(d+nb,1); %输入初值：uk(i)表示 u(k-i);yk=zeros(na,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 初值设置thetae_1=0.001*ones(na+nb+1+nc,1); %非常小的正数（这里不能为0）P=106*eye(na+nb+1+n

29、c);for k=1:Ltime(k)=k;y(k)=-a(2:na+1)*yk+b*uk(d:d+nb)+c*xi(k);xik; %采集输出数据%递推增广最小二乘法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); %白噪声的估计值%提取辨识参数ae=1 thetae(1:na,k); be=thetae(na+1:na+nb+1,k); ce=

30、1 thetae(na+nb+2:na+nb+1+nc,k);if abs(be(2)0.9be(2)=sign(ce(2)*0.9; %MVC算法要求 B 稳定endif abs(ce(2)0.9ce(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);fo

31、r i=d+nb:-1:2uk(i)=uk(i-1);enduk(1)=u(k);for i=na:-1:220yk(i)=yk(i-1);endyk(1)=y(k);for i=nc:-1:2yrk(i)=yrk(i-1);yr(k)、y(k)100-10-200yr(k)y(k)50100150200k25030035040040 xik(i)=xik(i-1);xiek(i)=xiek(i-1);u(k)200-20-40endif nc0yrk(1)=yr(k);xik(1)=xi(k);xiek(1)=xie;endendfigure(1);subplot(2,1,1);plot(t

32、ime,yr(1:L),r:,time,y);xlabel(k); ylabel(y_r(k)、y(k);02150100150200k250300350400a1a2参数估计a0-1-2-3050100150200k2503003504001.5b0b1c1参数估计b、clegend(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:na,:);xlabel(k); y

33、label(参数估计 a);legend(a_1,a_2); axis(0 L -3 2);subplot(212)plot(1:L,thetae(na+1:na+nb+1+nc,:);xlabel(k); ylabel(参数估计 b、c);legend(b_0,b_1,c_1); axis(0 L 0 1.5);10.50050100150200k250300350400%最小方差直接自校正控制最小方差直接自校正控制设被控对象为：y(k)1.7y(k 1)0.7y(k 2) u(k 4)0.5u(k 5)(k)0.2(k 1)式中(k)为方差为 0.1 的白噪声。%最小方差直接自校正控制cl

34、ear 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,1); %输出初值yek=zeros(nc,1); %最优输出预测估计初值yrk=zeros(nc,1); %期望输出初值xik

35、=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=106*eye(na+nb+d+nc);for k=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+nf);-yek(1:nc);K=P

36、*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)%提取辨识参数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.9ce(2)=sign(ce(2)*

37、0.9;endif fe(1)0yek(1)=ye;yrk(1)=yr(k);xik(1)=xi(k);endendfigure(1);subplot(2,1,1);20yr(k)10(ky(k)y、)0(kry-10-20050100150200250300350400k4020)u(k0-20-40050100150200250300350400k4gc0、2gg1计估0c1数参-2050100150200250300350400k4f30f计f1估2f数2参1f3f40050100150200250300350400kplot(time,yr(1:L),r:,time,y);xlabel

38、(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+nf+3:ng+2+nf+nc,:);xlabel(k); ylabel(参数估计 g、c);legend(g_0,g_1,c_1); axis(0 L -3 4);subplot(212)plot(1:L,theta

39、e(ng+2:ng+2+nf,:);xlabel(k); ylabel(参数估计 f);legend(f_0,f_1,f_2,f_3,f_4); axis(0 L 0 4);%广义最小方差控制（显示控制）广义最小方差控制（显示控制）设被控对象为如下开环不稳定非最小相位系统：y(k)1.7y(k 1)0.7y(k 2) u(k 4)0.5u(k 5)(k)0.2(k 1)式中(k)为方差为 0.1 的白噪声。%广义最小方差控制（显式控制律）clear all; close all;a=1 -1.7 0.7; b=1 2; c=1 0.2; d=4; %对象参数na=length(a)-1; nb

40、=length(b)-1; nc=length(c)-1; %na、nb、nc 为多项式 A、B、C 阶次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);yk=zeros(na,1); %输出初值yrk=zeros(nc,1); %期望输出初值xik=zeros(nc,1); %白噪声初值yr=10*ones(L/4,1);-on

41、es(L/4,1);ones(L/4,1);-ones(L/4+d,1); %期望输出xi=sqrt(0.1)*randn(L,1); %白噪声序列e,f,g=sindiophantine(a,b,c,d); %求解单步 Diophantine 方程CQ=conv(c,Q); FP=conv(f,P); CR=conv(c,R); GP=conv(g,P); %CQ=C*Qfor k=1:Ltime(k)=k;y(k)=-a(2:na+1)*yk+b*uk(d:d+nb)+c*xi(k);xik; %采集输出数据u(k)=(-Q(1)*CQ(2:nc+nq+1)*uk(1:nc+nq)/b(1

42、)-FP(2:np+nf+1)*uk(1:np+nf).+CR*yr(k+d:-1:k+d-min(d,nr+nc); yrk(1:nr+nc-d).-GP*y(k); yk(1:np+ng)/(Q(1)*CQ(1)/b(1)+FP(1); %求控制量%更新数据for i=d+nb:-1:2uk(i)=uk(i-1);enduk(1)=u(k);for i=na:-1:2yk(i)=yk(i-1);endyk(1)=y(k);u(k)yr(k)、y(k)100-10-20020yr(k)y(k)50100150200k2503003504005for i=nc:-1:2yrk(i)=yrk(i

43、-1);xik(i)=xik(i-1);endif nc0yrk(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);0-5050100150200k250300350400plot(time,u);xlabel(k); ylabel(u(k);%广义最小方差自校正控制（间接算法）广义最小方差自校正控制（间接算法）设被控对象为如下开环不稳定非最小相位系统：y(k)1.7y(k

44、1)0.7y(k 2) u(k 4)0.5u(k 5)(k)0.2(k 1)式中(k)为方差为 0.1 的白噪声。%广义最小方差自校正控制（间接算法）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; %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

45、(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(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 初值设置thetae_1=0.001*ones(na+nb+1+

46、nc,1);%非常小的正数（此处不能为0）P=106*eye(na+nb+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=-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(

47、1:na,k); be=thetae(na+1:na+nb+1,k); ce=1 thetae(na+nb+2:na+nb+1+nc,k);if abs(ce(2)0.9ce(2)=sign(ce(2)*0.9;ende,f,g=sindiophantine(ae,be,ce,d); %求解单步 Diophantine 方程CQ=conv(ce,Q); FP=conv(f,Pw); CR=conv(ce,R); GP=conv(g,Pw); %CQ=Ce*Qu(k)=(-Q(1)*CQ(2:nc+nq+1)*uk(1:nc+nq)/be(1)-FP(2:np+nf+1)*uk(1:np+nf)

48、.+CR*yr(k+d:-1:k+d-min(d,nr+nc); yrk(1:nr+nc-d).-GP*y(k); yk(1:np+ng)/(Q(1)*CQ(1)/be(1)+FP(1);%求控制量%更新数据thetae_1=thetae(:,k);yr(k)、y(k)20100-10-200yr(k)y(k)for i=d+nb:-1:2uk(i)=uk(i-1);enduk(1)=u(k);for i=na:-1:2yk(i)=yk(i-1);50100150200k250300350400105u(k)endyk(1)=y(k);for i=nc:-1:2yrk(i)=yrk(i-1);

49、xik(i)=xik(i-1);xiek(i)=xiek(i-1);endif nc00-5-10050100150200k2503003504002.5yrk(1)=yr(k);xik(1)=xi(k);xiek(1)=xie;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(0 L -20 20);subplot(2,1,2);plot(time,u);xlabel(k); ylabel(u(k); axis(0

50、 L -10 10);figure(2)plot(1:L,thetae);xlabel(k); ylabel(辨识参数 a、b、c);辨识参数a、b、c21.510.50-0.5-1-1.5-20a1a2b0b1c150100150200k250300350400legend(a_1,a_2,b_0,b_1,c_1); axis(0 L -2 2.5);%广义最小方差自校正控制（直接算法）广义最小方差自校正控制（直接算法）设被控对象为如下开环不稳定非最小相位系统：y(k)1.7y(k 1)0.7y(k 2) u(k 4)0.5u(k 5)(k)0.2(k 1)式中(k)为方差为 0.1 的白噪