自适应控制程序.doc

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

1、如有侵权，请联系网站删除，仅供学习与交流自适应控制程序【精品文档】第 22 页% M 序列及其逆序列的产生设M序列M(k)由如下4位移位寄存器产生：S(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:L IM=xor(S,x4); %进行异或运算，产生逆M序列 if IM=0 u(k)=-1; else u(k)=1; end S=not(S); %产生方波 M(k)=xor(x3,x4); %产生M序列 x4=x3;x3=x2;x

2、2=x1;x1=M(k); %寄存器移位endsubplot(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序列);%白噪声及有色噪声序列的产生设x(k) 为均值为0，方差为1的高斯白噪声序列，e(k)为有色噪声序列： 高斯白噪声序列 x(k)在Matlab中由rand()函数产生，程序如下：clear all; close

3、 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:L e(k)=-d(2:nd+1)*ek+c*xi(k);xik; %产生有色噪声 %数据更新 for i=nd:-1:2 ek(i)=ek(i-1); end ek(1)=e(k); for i=nc:-1:2 xik(i)=xik(i-1); e

4、nd xik(1)=xi(k);endsubplot(2,1,1);plot(xi);xlabel(k);ylabel(噪声幅值);title(白噪声序列);subplot(2,1,2);plot(e);xlabel(k);ylabel(噪声幅值);title(有色噪声序列);%批处理最小二乘参数估计(LS)考虑如下系统：式中x(k)为方差为1的白噪声。clear all;a=1 -1.5 0.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);

5、%输入初值x1=1;x2=1;x3=1;x4=0;S=1;%移位寄存器初值，方波初值xi=rand(L,1);%白噪声序列theta=a(2:na+1);b; %对象参数真值for k=1:L phi(k,:)=-yk;uk(d:d+nb); %phi(k,:)为行向量，便于组成phi矩阵 y(k)=phi(k,:)*theta+xi(k); %采集输出数据 IM=xor(S,x4); if IM=0 u(k)=-1; else u(k)=1; end S=not(S);M=xor(x3,x4); %产生M序列 %更新数据 x4=x3;x3=x2;x2=x1;x1=M; for i=nb+d:

6、-1:2 uk(i)=uk(i-1); end uk(1)=u(k); for i=na:-1:2 yk(i)=yk(i-1); end yk(1)=y(k);endthetaevaluation=inv(phi*phi)*phi*y %计算参数估计值thetaevaluation = -1.5362 0.6802 1.00680.4864%遗忘因子递推最小二乘参数估计(FFRLS)考虑如下系统：式中x(k)为均值为0、方差为0.1的白噪声，对象时变参数为：取遗忘因子l=0.98，clear all; close all;a=1 -1.5 0.7;b=1 0.5;d=3; %对象参数na=le

7、ngth(a)-1;nb=length(b)-1; %计算阶次L=1000;%数据长度uk=zeros(d+nb,1);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:L if k=501 a=1 -1 0.4;b=1.5 0.2; %对象

8、参数突变 end theta(:,k)=a(2:na+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(na+nb+1)-K*phi)*P/lambda; %更新数据 thetae_1=thetae(:,k); for i=d+nb:-1:2 uk(i)=uk(i-1); enduk(1)=u(k); for i=na:-1

9、:2 yk(i)=yk(i-1); end yk(1)=y(k);endsubplot(2,1,1);plot(1:L,thetae(1:na,:);hold on;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(

10、0 L -0.5 2);%增广递推最小二乘参数估计(ERLS)考虑如下系统：式中x(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

11、,1); %输入采用方差为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:L phi=-yk;uk(d:d+nb);xik; %测量向量 y(k)=phi*theta+xi(k); %采样输出数据phie=-yk;uk(d:d+nb);xiek; %估计的测量向量

12、 %增广递推最小二乘公式 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:2 uk(i)=uk(i-1); end uk(1)=u(k); for i=na:-1:2 yk(i)=yk(i-1); end yk(1)=y(k); for i=nc:-1:2 xik(i)=xik(i-1);

13、xiek(i)=xiek(i-1); end xik(1)=xi(k); xiek(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);ylabe

14、l(参数估计c);legend(c_1,c_2);axis(0 L -2 2);递推最小二乘参数估计(RLS)考虑如下系统：式中x(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=sqrt(0.1)*randn(L,1);

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

16、据 thetae_1=thetae(:,k); for i=d+nb:-1:2 uk(i)=uk(i-1); end uk(1)=u(k); for i=na:-1:2 yk(i)=yk(i-1); end yk(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)考虑如下系统：式中x(k)为方差为0.1的白噪声。取期望输出yr(k)为幅值为10的方波信号。clear all;close a

17、ll;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);-ones(L/4+d,1);%期望输出xi=sqrt(0.1)*randn(L,1);%方差为0.1的白噪声序列h,f,g=singlediopha

18、ntine(a,b,c,d);%求解单步Diophantine方程for k=1:L time(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:2 uk(i)=uk(i-1); end uk(1)=u(k); for i=na:-1:2 yk(i)=yk(i-1); end yk(1)=y(k); for

19、 i=nc:-1:2 yrk(i)=yrk(i-1); xik(i)=xik(i-1); end if nc0 yrk(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); 单步Diophantine方程求解求解下列系统的单步Diophantine方程：（1）（2）（3）%单步Diophantine方程的求

20、解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(a,b,c,d) %单步Diophantine方程求解na=length(a)-1;nb=length(b)-1;nc=length(c)-1;%计算阶次ne=

21、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+1 e(i)=0; for j=2:i e(i)=e(i)+e(i+1-j)*ad(j); end e(i)=cd(i)-e(i);%计算eiendfor i=1:ng+1 g(i)=0; for j=1:ne+1 g(i)=g(i)+e(ne+2-1)*ad(i+j); end g(i)=cd(i+d)-g(i);%计算giendf=conv(b,e);%计算Fe = 1.0000 1

22、.5000 1.5500f = 1.0000 2.0000 2.3000 0.7750g =1.2750 -1.0850多步Diophantine方程求解求解如下系统的多步Diophantine方程，并去预测长度N=3%多步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

23、 E,F,G=multidiophantine(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=nc G(1,:)=c(2:nc+1) zeros(1

24、,na-nc)-a(2:na+1); %令c(nc+2)=c(nc+3)=.=0else G(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-1 for i=1:j E(j+1,i)=E(j,i); end E(j+1,j+1)=G(j,1); for i=2:na G(j+1,i-1)=G(j,i)-G(j,1)*a(i); end G(j+1,na)=-G(j,1)*a(na+1); F(j+1,:)=conv(b,E(j+1,:);end%最小方差自校正控制用最小方差自校正

25、控制算法对以下系统进行闭环控制：式中x(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=zeros(d+nb,1); %输入初值：uk(i)表示u(k-i);yk=zeros(na,1); %输出初值yrk=zero

26、s(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+nc);for k=1:L time(k)=k; y(k)=-a(2:na+1)*yk+b*uk(d:d+nb)+c*xi(k

27、);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=1 thetae(na+nb+2:na+nb+1+nc,k); if abs(be(2)0.9 be(2

28、)=sign(ce(2)*0.9; %MVC算法要求B稳定 end if abs(ce(2)0.9 ce(2)=sign(ce(2)*0.9; end e,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:2 uk(i)=uk(i-1); end uk(1)=u(k); f

29、or i=na:-1:2 yk(i)=yk(i-1); end yk(1)=y(k); for i=nc:-1:2 yrk(i)=yrk(i-1); xik(i)=xik(i-1); xiek(i)=xiek(i-1); end if nc0 yrk(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,

30、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); ylabel(参数估计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);%最小方差直接自校正控制设被控对象为：式中x(k)为方差

31、为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,1); %输出初值yek=zeros(nc,1); %最优输出预测估计初值yrk=zeros(n

32、c,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=106*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:

33、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) %提取辨识参数 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 ab

34、s(ce(2)0.9 ce(2)=sign(ce(2)*0.9; end if fe(1)0 yek(1)=ye; yrk(1)=yr(k); xik(1)=xi(k); 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 L -40 40);figure(2)subplot(2

35、11)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,thetae(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);%广义最小方差控制（显示控制）设被控对象为如下开环不稳定非最小相位系统：式中x(k)为方差为0.1的白噪声。%广义最小方差控制

36、（显式控制律）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的阶次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(n

37、a,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=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:L time(k)=k; y(k)=-a(2:na+1

38、)*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)-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:2 uk(i)=uk(i-1); end uk(1)=u(k); for i=na:-1:2 yk(i)=yk(i-1); end yk(1)

39、=y(k); for i=nc:-1:2 yrk(i)=yrk(i-1); xik(i)=xik(i-1); end if nc0 yrk(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); %广义最小方差自校正控制（间接算法）设被控对象为如下开环不稳定非最小相位系统：式中x(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-