系统辨识与自适应控制作业

系统辨识与自适应控制学院：专业：学号：姓名：系统辨识与自适应控制作业一、对时变系统进行参数估计。系统方程为：y(k)+a(k)y(k-1)=b(k)u(k-1)+e(k)其中：e(k)为零均值噪声，a(k)=b(k)=要求：1对定常系统（a=0.8,b=0.5）进行结构（阶数）确定和参数估计；2对时变系统，λ取不同值（0.9——0.99）时对系统辨识结果和过程进行比较、讨论3对辨识结果必须进行残差检验解：一（1）：分析：采用最小二乘法（LS）：最小二乘的思想就是寻找一个的估计值ˆ，使得各次测量的),1(miZi与由估计ˆ确定的量测估计ˆˆiiHZ之差的平方和最小，由于此方法兼顾了所有方程的近似程度，使整体误差达到最小，因而对抑制误差是有利的。在此，我应用批处理最小二乘法，收敛较快，易于理解，在系统参数估计应用中十分广泛。作业程序：clearall;a=[10.8]';b=[0.5]';d=3;%对象参数na=length(a)-1;nb=length(b)-1;%na、nb为A、B阶次L=500;%数据长度uk=zeros(d+nb,1);%输入初值：uk(i)表示u(k-i)yk=zeros(na,1);%输出初值x1=1;x2=1;x3=1;x4=0;S=1;%移位寄存器初值、方波初值xi=randn(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);%产生逆M序列ifIM==0u(k)=-1;elseu(k)=1;endS=not(S);M=xor(x3,x4);%产生M序列%更新数据x4=x3;x3=x2;x2=x1;x1=M;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);endthetae=inv(phi'*phi)*phi'*y'%计算参数估计值thetae结果：thetae=0.7787，0.5103真值=0.8，0.5解：一（2）：采用遗忘因子递推最小二乘参数估计；其仿真算法如下：Step1：设置初值、，及遗忘因子，输入初始数据；Step2：采样当前输入和输出数据；Step3：利用含有遗忘因子的递推公式计算、和；Step4：k=k+1，返回Step2继续循环。其中：对时变系统，λ取不同值（0.9——0.99）时对系统辨识结果和过程进行比较、讨论作业程序：clearall;closeall;a=[10.8]’;b=[0.5]’;d=2;%对象参数na=length(a)-1;nb=length(b)-1;%na、nb为A、B阶数L=500;%仿真长度uk=zeros(d+nb,1);%输入初值yk=zeros(na,1);%输出初值u=randn(L,1);%输入采用白噪声序列xi=sqrt(0.1)*randn(L,1);%白噪声序列thetae_1=zeros(na+nb+1,1);p=10^6*eye(na+nb+1);lambda=0.95;%遗忘因子fork=1:L;ifk300;a=[10.6]';b=[0.3]';endthetae(:,k)=[a(2:na+1);b];phi=[-yk;uk(d:d+nb)];y(k)=phi'*thetae(:,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(1,2,1)plot([1:L],thetae(1:na,:));holdon;plot([1:L],thetae(1:na,:),'k:');xlabel('k');ylabel('参数估计a');axis([0L-0.52]);subplot(1,2,2)plot([1:L],thetae(na+1:na+nb+1,:));holdon;plot([1:L],thetae(na+1:na+nb+1,:),'k:');xlabel('k');ylabel('参数估计b');axis([0L-0.52]);仿真结果图1-1遗忘因子递推最小二乘法的参数估计结果（=0.95）图1-2遗忘因子递推最小二乘法的参数估计结果（λ=0.99）综上所述，可知从仿真图可以看出，当0k300时，a的值大约在0.8左右波动，b的值大约在0.5左右波动。当k300时，a、b的值发生变化，此时a的值大约在0.6左右波动，b的值大约在0.3左右波动。所以在此种情况下，可以达到辨识目的，但是由于波动较大，辨识效果不理想。同时，当λ=0.99时比前面所述λ=0.96时的辨识结果更为明显。这种情况下当估计值变化稳定时，估计值的变化率比较小，辨识效果较为理想。一（3）：根据残差平方和估计模型的阶次SISO过程的差分方程模型的输出残差为)(~kz，数据长度L，nHˆ为nˆ阶时的数据矩阵，nˆˆ为nˆ阶时的参数的估计量，nˆ为模型阶次估计值，0n为真实阶次，则残差平方和函数J：)(~1)ˆ()ˆ(1~~1)ˆ(12ˆˆˆˆˆˆ00kzLHzHzLzzLnJLnnknnnTnnnnTn残差平方和有这样的性质：当L足够大时，随着nˆ增加)ˆ(nJ先是显著地下降，当nˆ0n时，)ˆ(nJ值显著下降的现象就终止。这就是损失函数法来定阶的原理。程序：N=15;%4位移位寄存器产生的M序列的周期y1=1;y2=1;y3=1;y4=0;fori=1:N;x1=xor(y3,y4);x2=y1;x3=y2;x4=y3;y(i)=y4;ify(i)0.5,u(i)=-1;elseu(i)=1;endy1=x1;y2=x2;y3=x3;y4=x4;end%白噪声的产生A=19;x0=12;M=1024;fork=1:Nx=A*x0;x1=mod(x,M);v(k)=x1/512;x0=x1;end%辨识主程序z=zeros(7,N);zs=zeros(7,N);zm=zeros(7,N);zmd=zeros(7,N);z(1)=0;z(2)=0;zs(1)=0;zs(2)=0;zm(1)=0;zm(2)=0;zmd(1)=0;zmd(2)=0;theta=zeros(7,1);p0=10^6*eye(7,7);the=[theta,zeros(7,14)];e=zeros(7,15);fork=3:N;h=[-z(k-1),-z(k-2),u(k-1),u(k-2),v(k),v(k-1),v(k-2)]';x=inv(h'*p0*h+1);q=p0*h*x;d1=z(k)-h'*theta;theta1=theta+q*d1;zs(k)=-1.5*z(k-1)+0.7*z(k-2)+u(k-1)+0.5*u(k-2);zm(k)=[-z(k-1),-z(k-2),u(k-1),u(k-2)]*[theta1(1);theta1(2);theta1(3);theta1(4)];zmd(k)=h'*theta1;e(:,k)=theta1-theta;theta=theta1;the(:,k)=theta1;p1=p0-q*q'*[h'*p0*h+1];p0=p1;endfigure(1);i=1:N;plot(i,the(1,:),'r',i,the(2,:),'r:',i,the(3,:),'b',i,the(4,:),'b:',i,the(5,:),'g',i,the(6,:),'g:',i,the(7,:),'g+')title('辨识曲线')figure(2);i=1:N;plot(i,e(1,:),'r',i,e(2,:),'r:',i,e(3,:),'b',i,e(4,:),'b:',i,e(5,:),'g',i,e(6,:),'g:',i,e(7,:),'r+')title('辨识误差曲线')图1-3-1辨识曲线图1-3-2辨识误差曲线由以上对比曲线可以看出，经过最小二乘法估计得到的数据与实际数据之间虽然存在区别，但是基本符合要求,说明辨识效果较好。二、已知系统方程：y(k)+1.2y(k-1)+0.35y(k-2)=u(k-2)+0.4u(k-3)+e(k)+1.1e(k-1)+0.28e(k-2)其中e(k)为白噪声序列。要求：1、对系统进行参数递推估计；2、进行最小方差控制器设计；3、进行基于最小方差控制的自适应控制系统仿真分析（设输入信号为方波信号）。二（1）：程序：clearall;closeall;a=[11.20.35]';b=[10.4]';c=[11.10.28]';d=3;%对象参数na=length(a)-1;nb=length(b)-1;nc=length(c)-1;%na、nb、nc为A、B、C阶次L=1000;%仿真长度uk=zeros(d+nb,1);%输入初值：uk(i)表示u(k-i)yk=zeros(na,1);%输出初值xik=zeros(nc,1);%噪声初值xiek=zeros(nc,1);%噪声估计初值u=randn(L,1);%输入采用白噪声序列xi=sqrt(0.1)*randn(L,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);fork=1:Lphi=[-yk;uk(d:d+nb);xik];y(k)=phi'*theta+xi(k);%采集输出数据phie=[-yk;uk(d:d+nb);xiek];%组建phie%递推增广最小二乘法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([0L-22]);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]);对系统进行最广最小二乘的递推估计（RELS），仿真结果如下图2.1所示。当仿真步数K