自适应控制程序.docx

1、自适应控制程序% 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;x2=x1;x1=M(k); %寄存器移位endsubplot(

2、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)为有色噪声序列： 高斯白噪声序列 (k)在Matlab中由rand()函数产生，程序如下：clear all; close all;L=500; %仿真长度d=1 -1.5 0.7 0.

3、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); end xik(1)=xi(k);endsubplot(2,1,1

4、);plot(xi);xlabel(k);ylabel(噪声幅值);title(白噪声序列);subplot(2,1,2);plot(e);xlabel(k);ylabel(噪声幅值);title(有色噪声序列);%批处理最小二乘参数估计(LS)考虑如下系统：式中 (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); %输入初值x1=1;x2=1;x3=1;x4=0;S=1;%移位

5、寄存器初值，方波初值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:-1:2 uk(i)=uk(i-1); end uk(1)=u(

6、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)考虑如下系统：式中 (k)为均值为0、方差为0.1的白噪声，对象时变参数为：取遗忘因子 =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

7、=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; %对象参数突变 end theta(:,k)=a(2:na+1);b;

8、 %对象参数真值 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:2 yk(i)=yk(i-1); end yk(1)=y(k)

9、;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(0 L -0.5 2); %增广递推最小二乘参数估计(ERLS)

10、考虑如下系统：式中 (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); %输入采用方差为1的白噪声序列xi=sqrt(0.1

11、)*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; %估计的测量向量 %增广递推最小二乘公式 K=P*phie/(1+phie*P

12、*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); xiek(i)=xiek(i-1); end xik(1)=x

13、i(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);ylabel(参数估计c);legend(c_1,c_2);axis(0

14、 L -2 2);递推最小二乘参数估计(RLS)考虑如下系统：式中 (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); %方差为0.1的白噪声干扰序列theta=a(2:na+1)

15、;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: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; %更新数据 thetae_1=thetae(:,k); for i=d

16、+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)考虑如下系统：式中 (k)为方差为0.1的白噪声。取期望输出yr(k)为幅值为10的方波信号。clear all;close all;a=1 -1.7 0.7;b=1 0.5;c=1 0.2

17、;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=singlediophantine(a,b,c,d);%求解单步Diophantine

18、方程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 i=nc:-1:2 yrk(i)=yrk(i-1); xik

19、(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方程的求解clear all;a=1 -1.5 0.7; b=1 0.

20、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=d-1;ng=na-1;%E,G的阶次ad=a,zeros(1

21、,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.5000 1.5500f = 1.0000 2.0000 2

22、.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 E,F,G=multidiophantine(a,b,c,

23、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,na-nc)-a(2:na+1); %令c(nc+

24、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%最小方差自校正控制用最小方差自校正控制算法对以下系统进行闭环控制：式中 (k)为方差为

25、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=zeros(nc,1); %期望输出初值xik=zeros(

26、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);xik; %采集输出数据 %递推增广最小二乘法

27、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)=sign(ce(2)*0.9; %MVC算法要求

28、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); for 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(