自适应控制程序

%M序列及其逆序列的产生设M序列{M(k)}由如下4位移位寄存器产生：34iiixxx{S(k)}为方波序列，逆M序列{IM(k)={M(k)S(k)}clearall;closeall;L=60;%序列长度x1=1;x2=1;x3=1;x4=0;%移位寄存器初值S=1;%方波初值fork=1:LIM=xor(S,x4);%进行异或运算，产生逆M序列ifIM==0u(k)=-1;elseu(k)=1;endS=not(S);%产生方波M(k)=xor(x3,x4);%产生M序列x4=x3;x3=x2;x2=x1;x1=M(k);%寄存器移位endsubplot(2,1,1);stairs(M);grid;axis([0L/2-0.51.5]);xlabel('k');ylabel('M序列幅值');title('M序列');subplot(2,1,2);stairs(u);grid;axis([0L-1.51.5]);xlabel('k');ylabel('逆M序列幅值');title('逆M序列');%白噪声及有色噪声序列的产生设(k)为均值为0，方差为1的高斯白噪声序列，e(k)为有色噪声序列：11211123()10.50.2()()()()()()11.50.70.1CzzzekGzkkkDzzzz高斯白噪声序列(k)在Matlab中由rand()函数产生，程序如下：clearall;closeall;L=500;%仿真长度d=[1-1.50.70.1];c=[10.50.2];%分子分母多项式系数nd=length(d)-1;nc=length(c)-1;%阶次xik=zeros(nc,1);%白噪声初值ek=zeros(nd,1);xi=randn(L,1);%产生均值为0，方差为1的高斯白噪声序列fork=1:Le(k)=-d(2:nd+1)*ek+c*[xi(k);xik];%产生有色噪声%数据更新fori=nd:-1:2ek(i)=ek(i-1);endek(1)=e(k);fori=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)考虑如下系统：()1.5(1)0.7(2)(3)0.5(4)()ykykykukukk式中(k)为方差为1的白噪声。clearall;a=[1-1.50.7]';b=[10.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];%对象参数真值fork=1:Lphi(k,:)=[-yk;uk(d:d+nb)]';%phi(k,:)为行向量，便于组成phi矩阵y(k)=phi(k,:)*theta+xi(k);%采集输出数据IM=xor(S,x4);ifIM==0u(k)=-1;elseu(k)=1;endS=not(S);M=xor(x3,x4);%产生M序列%更新数据x4=x3;x3=x2;x2=x1;x1=M;fori=nb+d:-1:2uk(i)=uk(i-1);enduk(1)=u(k);fori=na:-1:2yk(i)=yk(i-1);endyk(1)=y(k);endthetaevaluation=inv(phi'*phi)*phi'*y'%计算参数估计值thetaevaluation=-1.53620.68021.00680.4864%遗忘因子递推最小二乘参数估计(FFRLS)考虑如下系统：1201()(1)(2)(3)(4)()ykaykaykbukbukk式中(k)为均值为0、方差为0.1的白噪声，对象时变参数T1201()[,,,]kaabb为：[1.5,0.7,1,0.5],500()[1,0.4,1.5,0.2],500TTkkk取遗忘因子=0.98，clearall;closeall;a=[1-1.50.7]';b=[10.5]';d=3;%对象参数na=length(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=10^6*eye(na+nb+1);lambda=0.98;%遗忘因子范围[0.91]fork=1:Lifk==501a=[1-10.4]';b=[1.50.2]';%对象参数突变endtheta(:,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);fori=d+nb:-1:2uk(i)=uk(i-1);enduk(1)=u(k);fori=na:-1:2yk(i)=yk(i-1);endyk(1)=y(k);endsubplot(2,1,1);plot([1:L],thetae(1:na,:));holdon;plot([1:L],theta(1:na,:),'k:');xlabel('k');ylabel('参数估计a');legend('a_1','a_2');axis([0L-22]);subplot(2,1,2);plot([1:L],thetae(na+1:na+nb+1,:));holdon;plot([1:L],theta(na+1:na+nb+1,:),'k:');xlabel('k');ylabel('参数估计b');legend('b_0','b_1');axis([0L-0.52]);%增广递推最小二乘参数估计(ERLS)考虑如下系统：()1.5(1)0.7(2)(3)0.5(4)()(1)0.2(2)ykykykukukkkk式中(k)为方差为0.1的白噪声。选择方差为1的白噪声作为输入信号u(k).clearall;closeall;a=[1-1.50.7]';b=[10.5]';c=[1-10.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)*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=10^6*eye(na+nb+1+nc);lambda=0.98;%遗忘因子范围[0.91]fork=1:Lphi=[-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*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);fori=d+nb:-1:2uk(i)=uk(i-1);enduk(1)=u(k);fori=na:-1:2yk(i)=yk(i-1);endyk(1)=y(k);fori=nc:-1:2xik(i)=xik(i-1);xiek(i)=xiek(i-1);endxik(1)=xi(k);xiek(1)=xie;endfigure(1)plot([1:L],thetae(1:na,:));xlabel('k');ylabel('参数估计a');legend('a_1','a_2');axis([0L-22]);figure(2)plot([1:L],thetae(na+1:na+nb+1,:));xlabel('k');ylabel('参数估计b');legend('b_0','b_1');axis([0L01.5]);figure(3)plot([1:L],thetae(na+nb+2:na+nb+nc+1,:));xlabel('k');ylabel('参数估计c');legend('c_1','c_2');axis([0L-22]);递推最小二乘参数估计(RLS)考虑如下系统：()1.5(1)0.7(2)(3)0.5(4)()ykykykukukk式中(k)为方差为0.1的白噪声。clearall;closeall;a=[1-1.50.7]';b=[10.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);b];%对象参数真值theta=[-1.5,0.7;1,0.5]thetae_1=zeros(na+nb+1,1);%参数初值为4x1的全零矩阵P=10^6*eye(na+nb+1);fork=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)-K*phi')*P;%更新数据thetae_1=thetae(:,k);fori=d+nb:-1:2uk(i)=uk(i-1);enduk(1)=u(k);fori=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([0L