系统辨识练习题方法一：%递推最小二乘参数估计（RLS）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、nb为A、B阶次L=480; %仿真长度uk=zeros(d+nb,1); %输入初值：uk(i)表示u(k-i)yk=zeros(na,1); %输出初值u=randn(L,1); %输入采用白噪声序列xi=sqrt(0.1)*randn(L,1); %白噪声序列theta=[a(2:na+1);b]; %对象参数真值thetae_1=zeros(na+nb+1,1); %thetae初值P=10^6*eye(na+nb+1);for k=1:Lphi=[-yk;uk(d:d+nb)]; %此处phi为列向量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+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]);方法三：%遗忘因子递推最小二乘参数估计（FFRLS）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、nb为A、B阶次L=1000; %仿真长度uk=zeros(d+nb,1); %输入初值：uk(i)表示u(k-i)yk=zeros(na,1); %输出初值u=randn(L,1); %输入采用白噪声序列xi=sqrt(0.1)*randn(L,1); %白噪声序列thetae_1=zeros(na+nb+1,1); %thetae初值P=10^6*eye(na+nb+1);lambda=0.98; %遗忘因子范围[0.9 1]for 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'*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(1,2,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(1,2,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]);方法四：%递推极大似然参数估计（RML）clear all; close all;a=[1 -1.5 0.7]'; b=[1 0.5]'; c=[1 -0.5]'; d=1; %对象参数na=length(a)-1; nb=length(b)-1; nc=length(c)-1; %na、nb、nc为A、B、C阶次nn=max(na,nc); %用于yf(k-i)、uf(k-i)更新L=480; %仿真长度uk=zeros(d+nb,1); %输入初值：uk(i)表示u(k-i)yk=zeros(na,1); %输出初值xik=zeros(nc,1); %白噪声初值xiek=zeros(nc,1); %白噪声估计初值yfk=zeros(nn,1); %yf(k-i)ufk=zeros(nn,1); %uf(k-i)xiefk=zeros(nc,1); %ξf(k-i)u=randn(L,1); %输入采用白噪声序列xi=randn(L,1); %白噪声序列thetae_1=zeros(na+nb+1+nc,1); %参数估计初值P=eye(na+nb+1+nc);for k=1:Ly(k)=-a(2:na+1)'*yk+b'*uk(d:d+nb)+c'*[xi(k);xik]; %采集输出数据%构造向量phi=[-yk;uk(d:d+nb);xiek];xie=y(k)-phi'*thetae_1;phif=[-yfk(1:na);ufk(d:d+nb);xiefk];%递推极大似然参数估计算法K=P*phif/(1+phif'*P*phif);thetae(:,k)=thetae_1+K*xie;P=(eye(na+nb+1+nc)-K*phif')*P;yf=y(k)-thetae(na+nb+2:na+nb+1+nc,k)'*yfk(1:nc); %yf(k)uf=u(k)-thetae(na+nb+2:na+nb+1+nc,k)'*ufk(1:nc); %uf(k)xief=xie-thetae(na+nb+2:na+nb+1+nc,k)'*xiefk(1:nc); %xief(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);xiefk(i)=xiefk(i-1);endxik(1)=xi(k);xiek(1)=xie;xiefk(1)=xief;for i=nn:-1:2yfk(i)=yfk(i-1);ufk(i)=ufk(i-1);endyfk(1)=yf;ufk(1)=uf;endfigure(1)plot([1:L],thetae(1:na,:),[1:L],thetae(na+nb+2:na+nb+1+nc,:)); xlabel('k'); ylabel('参数估计a、c');legend('a_1','a_2','c_1'); 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])自适应控制习题(1) %可调增益MIT-MRACclear all; close all;h=0.1; L=100/h; %数值积分步长、仿真步数num=[1]; den=[1 1 1]; n=length(den)-1; %对象参数kp=1; [Ap,Bp,Cp,Dp]=tf2ss(kp*num,den); %传递函数型转换为状态空间型km=1; [Am,Bm,Cm,Dm]=tf2ss(km*num,den); %参考模型参数gamma=0.1; %自适应增益yr0=0; u0=0; e0=0; ym0=0; %初值xp0=zeros(n,1); xm0=zeros(n,1); %状态向量初值kc0=0; %可调增益初值r=0.1; yr=r*[ones(1,L/4) -ones(1,L/4) ones(1,L/4) -ones(1,L/4)]; %输入信号for k=1:Ltime(k)=k*h;xp(:,k)=xp0+h*(Ap*xp0+Bp*u0);yp(k)=Cp*xp(:,k)+Dp*u0; %计算ypxm(:,k)=xm0+h*(Am*xm0+Bm*yr0);ym(k)=Cm*xm(:,k)+Dm*yr0; %计算yme(k)=ym(k)-yp(k); %e=ym-ypkc=kc0+h*gamma*e0*ym0; %MIT自适应律u(k)=kc*yr(k); %控制量%更新数据yr0=yr(k); u0=u(k); e0=e(k); ym0=ym(k);xp0=xp(:,k); xm0=xm(:,k);kc0=kc;endplot(time,ym,'r',time,yp,':');xlabel('t'); ylabel('y_m(t)、y_p(t)');%axis([0 L*h -10 10]);legend('y_m(t)','y_p(t)');(2)%可调增益MIT-MRACclear all; close all;h=0.1; L=100/h; %数值积分步长、仿真步数num=[1]; den=[1 1 1]; n=length(den)-1; %对象参数kp=1; [Ap,Bp,Cp,Dp]=tf2ss(kp*num,den); %传递函数型转换为状态空间型km=1; [Am,Bm,Cm,Dm]=tf2ss(km*num,den); %参考模型参数gamma=0.1; %自适应增益yr0=0; u0=0; e0=0; ym0=0; %初值xp0=zeros(n,1); xm0=zeros(n,1); %状态向量初值kc0=0; %可调增益初值r=1; yr=r*[ones(1,L/4) -ones(1,L/4) ones(1,L/4) -ones(1,L/4)]; %输入信号for k=1:Ltime(k)=k*h;xp(:,k)=xp0+h*(Ap*xp0+Bp*u0);yp(k)=Cp*xp(:,k)+Dp*u0; %计算ypxm(:,k)=xm0+h*(Am*xm0+Bm*yr0);ym(k)=Cm*xm(:,k)+Dm*yr0; %计算yme(k)=ym(k)-yp(k); %e=ym-ypkc=kc0+h*gamma*e0*ym0; %MIT自适应律u(k)=kc*yr(k); %控制量%更新数据yr0=yr(k); u0=u(k); e0=e(k); ym0=ym(k);xp0=xp(:,k); xm0=xm(:,k);kc0=kc;endplot(time,ym,'r',time,yp,':');xlabel('t'); ylabel('y_m(t)、y_p(t)');%axis([0 L*h -10 10]);legend('y_m(t)','y_p(t)');(3)(1)%可调增益MIT-MRAC归一化算法clear all; close all;h=0.1; L=100/h; %数值积分步长和仿真步数num=[1]; den=[1 1 1]; n=length(den)-1; %对象参数kp=1; [Ap,Bp,Cp,Dp]=tf2ss(kp*num,den); %传递函数型转换为状态空间型km=1; [Am,Bm,Cm,Dm]=tf2ss(km*num,den); %参考模型参数gamma=0.1; %自适应增益alpha=0.01; beta=2;yr0=0; u0=0; e0=0; ym0=0; %初值xp0=zeros(n,1); xm0=zeros(n,1); %状态向量初值kc0=0; %可调增益初值r=0.1; yr=r*[ones(1,L/4) -ones(1,L/4) ones(1,L/4) -ones(1,L/4)]; %输入信号for k=1:Ltime(k)=k*h;xp(:,k)=xp0+h*(Ap*xp0+Bp*u0);yp(k)=Cp*xp(:,k)+Dp*u0; %计算ypxm(:,k)=xm0+h*(Am*xm0+Bm*yr0);ym(k)=Cm*xm(:,k)+Dm*yr0; %计算yme(k)=ym(k)-yp(k); %e=ym-ypDD=e0*ym0/km/(alpha+(ym0/km)^2);if DD<-betaDD=-beta;endif DD>betaDD=beta;endkc=kc0+h*gamma*DD; %MIT自适应律u(k)=kc*yr(k); %控制量%更新数据yr0=yr(k); u0=u(k); e0=e(k); ym0=ym(k);xp0=xp(:,k); xm0=xm(:,k);kc0=kc;endplot(time,ym,'r',time,yp,':');xlabel('t'); ylabel('y_m(t)、y_p(t)');legend('y_m(t)','y_p(t)');(2)%可调增益MIT-MRAC归一化算法clear all; close all;h=0.1; L=100/h; %数值积分步长和仿真步数num=[1]; den=[1 1 1]; n=length(den)-1; %对象参数kp=1; [Ap,Bp,Cp,Dp]=tf2ss(kp*num,den); %传递函数型转换为状态空间型km=1; [Am,Bm,Cm,Dm]=tf2ss(km*num,den); %参考模型参数gamma=0.1; %自适应增益alpha=0.01; beta=2;yr0=0; u0=0; e0=0; ym0=0; %初值xp0=zeros(n,1); xm0=zeros(n,1); %状态向量初值kc0=0; %可调增益初值r=1; yr=r*[ones(1,L/4) -ones(1,L/4) ones(1,L/4) -ones(1,L/4)]; %输入信号for k=1:Ltime(k)=k*h;xp(:,k)=xp0+h*(Ap*xp0+Bp*u0);yp(k)=Cp*xp(:,k)+Dp*u0; %计算ypxm(:,k)=xm0+h*(Am*xm0+Bm*yr0);ym(k)=Cm*xm(:,k)+Dm*yr0; %计算yme(k)=ym(k)-yp(k); %e=ym-ypDD=e0*ym0/km/(alpha+(ym0/km)^2);if DD<-betaDD=-beta;endif DD>betaDD=beta;endkc=kc0+h*gamma*DD; %MIT自适应律u(k)=kc*yr(k); %控制量%更新数据yr0=yr(k); u0=u(k); e0=e(k); ym0=ym(k);xp0=xp(:,k); xm0=xm(:,k);kc0=kc;endplot(time,ym,'r',time,yp,':');xlabel('t'); ylabel('y_m(t)、y_p(t)');legend('y_m(t)','y_p(t)');。