系统辨识大作业2

班级：自动化1班张青义，获得系统的脉冲传递函数G(2)和传递函数G(S);应用最小二乘辨识，获得脉冲响应序列8(k);同图显示两种方法的辨识效果图；应用相关最小二乘法，拟合对象的差分方程模型；构建时变对象，用最小二乘法和带遗忘因子的最小二乘法，(可以用辨识工具箱)辨识模型的参数，比较两种方法的辨识效果差异；答：根据index2搭建结构框图：aulKKNOTaaa相关分析法：利用结构框图得到UY和toutX×1234567X6由于是二阶水箱系统，可以假设系统的传递函数已知G(s)的结构，用长除法求得G(s)的s展开式，其系数等于从g(r)求得的各阶矩，然后求G(s)的参数。得到结果：在commandwindow中输入下列指令得到传递函数：最小二乘一次算法相关参数Estimate=inv(H'*H)*H'*得到相关系数为：Estimate=h=[-z(i-1);-z(i-2);M(i);M(i-1);M(iTheta(:,i-1)=Theta(:,i-2)+K*(z(i)-h'plot(i,Theta(1,:),i,Theta(2,:),i,Theta(3,:),i,Theta(4,:),i,The)第二题设SISO系统差分方程为(40分)y(k)=-a,y(k-1)-a₂Y(k-2)+bu(k-D)+b₂μ(k-2)+声。1)用最小二乘及递推最小二乘法估计0;2)用辅助变量法及其递推算法估计0;答：1.基本最小二乘法：n=2;%根据题目，本系统为2阶系统Y=Data(n+1:n+N,2);%输出数据矩阵disp('使用最小二乘算法所得结果如下：)%辨识参数数据输出如下：lsa1=thita(1),Isa2=thita(22.递推最小二乘法：%首先将TXT中的数据装载在Matlab中：whileWbs<1disp('从下列文件中选择所需加载文件：uy1.txtuy2.txtuy3.txt)Data=fscanf(Wbs,%g%g',[2inf]);%生成数据矩阵%递推最小二乘法n=2;%根据题目已知该系统阶数为2%为辨识参数向量赋初值，这里均取为0.001%直接给出矩阵Pn的初始状态POthita=[thita0,zeros(n*2,L)];%需辨识参数矩阵的初值及大小fori=n+1:1:L%这里i从第n+1个数开始form=1:n*2%输入矩阵赋值X=FIA1*PO*FIA1+1;%引入中间变量X为K递推公式中的求逆项，1维PO=P1;%作为下次迭代初值thita1=thitaO+K1*D1;%估值补偿公式，thita1为求被辨识的参数向量thitaO=thita1;%所得的参数向量作为下次迭代初值disp('递推最小二乘法所得结果如下：)Dta1=thita1(1),Dta2=thita1(2),Dtb1=thita1(3),Dtb2=thita1(4)%显示被辨识参数%首先将TXT中的数据装载在Matlab中：disp('从下列文件中选择所需加载文件：uy1.txtuy2.txtuy3.txt)Data=fscanf(Wbs,'%g%g’[2inf]);%生成数据矩阵%使用递推辅助变量法对2阶系统进行辨识，取bO=0%分别取出输入、输出数据u,y%首先使用最小二乘法，粗略得到初始的被辨识参数向量Y1=Data(3:L,2);%构造输出向量，也用于辅助变量法的输出向量%获得初值，以下使用辅助变量法Yf=fzOL*thita;%辅助模型输出%构造辅助变量矩阵，改变新的矩阵中的输出量，输入量不变thita=inv(Zfz'*fzOL)*Zfz'*Y1;%求取被估计参数的辅助变量估值disp('使用辅助变量法，得到如下辨识结果：)4.递推辅助变量法：%首先将TXT中的数据装载在Matlab中：disp('从下列文件中选择所需加载文件：uy1.txtuy2.txtuy3.txt)Data=fscanf(Wbs,%g%g',[2inf]);%生成数据矩阵%使用递推辅助变量法对2阶系统进行辨识，取bO=0%分别取出输入、输出数据u,y%取前100个数据利用最小二乘辨识Y1=Data(3:n0,2);%构造输出向量，也用于辅助变量法的输出向量Z=[-Y(2:n0-1)-Y(1:n0-2)U(2:n0-1)U(1%后面的数据计算，使用递推辅助变量法Y11=[Y(1);Y(2);Z*thikesy=[-Y(n-1);-Y(n-2disp(递推辅助变量法辨识结果：);%首先将TXT中的数据装载在Matlab中：disp('从下列文件中选择所需加载文件：uy1.txtuy2.txtuy3.txt)[Wbs,Message]=fopen(filename,'r');%Data=fscanf(Wbs,%g%g',[2inf]);%生成数据矩阵%使用广义最小二乘法，辨识2阶系统的参数%分别取出输入、输出数据u,y%首先使用最小二乘法，求系统参数向量初值Zgy1=Data(3:L,2);%输出矩阵y%得到初值后，以下使用广义最小二乘法进行辨识，由残差代替噪声误差E(1)=Y(1);%k=1时，残差的值%fl表示由残差代替噪声项，而辨识得到的滤波器的估计参数fl=inv(omiga'*omiga)*omiga%得到fl后，计算经过滤波的u,y数据；其中，k=1,2时：Data(k,1)=U(k)+fl(1)*Udisp('由广义最小二乘法得到的辨识结果如下：)实验结果：Data=fscanf(WbsN=500-n:%用前50组数据利用1sm估计出初值phi(:,1)=-y1(n:n+n0-1phi(:,2)=-y1(n-1:n+nO-2phi(:,3)=ul(n:n+nO-1phi(:,4)=ul(n-1:n+n0-2%计算残差K2=p20*omega*inv(1+omega'*p20*f1=f0+K2*(e(n+i-2)-p2=p20-K2*omega'*p2p2O=p2;%滤波yy(1:i,2)=y1(n-1:n+iuu(1:i,1)=ul(n:(n+i-1uu(1:i,2)=ul((n-1):(n+i-2%yg(k)=y(k)+f(1)*y(k-1)+f(%ug(k)=u(k)+f(1)*u(k-1)+f(ksai=[-yg(n+i-2)-yg(n+i-3)ug(n+i-2)ug(n+iK1=p10*ksai*inv(1+ksai’*p10*ksai);%K(N+1)值phi(1:i+1,3)=ul(n:n+i-1+phi(1:i+1,4)=ul(n-1:n+i-2disp('递推广义最小二乘法辨识结果：’);al=seta0(1),a2=seta0(2),bl=seta0实验结果：7.夏氏偏差修正法：%首先将TXT中的数据装载在Matlab中：[Wbs,Message]=fopen(filename,'r');%打开文件；Data=fscanf(Wbs,%g%g,[2inf]);%生成数据矩阵%使用夏氏修正法，对2阶系统进行参数辨识%首先使用最小二乘法，求系统参数向量初值Xsthita=Sgl2*Sgl3;%计算参数矩阵%得到初值后，以下使用夏氏修正法进行参数辨识：thitab0=0;%设偏差项的偏差初值为0E1=Zgl1-glOL*Xsthita;%计算残差Edisp('夏氏修正法)Xsa1=Xsthita(1),Xsa2=Xsthita(2),Xsb1=Xsthita(3),Xsb2=Xsthita(4)实验结果：5Xsb2=8.夏氏改良法：%首先将TXT中的数据装载在Matlab中：disp('从下列文件中选择所需加载文件：uy1.txtuy2.txtuy3.txt)[Wbs,Message]=fopen(filename,'r);%打开文件；Data=fscanf(Wbs,%g%g',[2inf]);%生成数据矩阵%最小二乘法，取bO=0n=2;%根据题目，本系统为2阶系统L=length(Data);%辨识所需数据的总长度N=L-n;%构造测量矩阵的数据长度FIA=zeros(N,2*n);%构造测量矩阵fori=1:N%测量矩阵赋值Y=Data(n+1:n+N,2);%输出数据矩阵thitaLS=thita;%最小二乘估值disp('使用夏氏改良法辨识结果为：)Xga1=thita(1),Xga2=thita(2),Xgb1=thit%首先将TXT中的数据装载在Matlab中：disp('从下列文件中选择所需加载文件：uy1.txtuy2.txtuy3.txt')Data=fscanf(Wbs,’%g%g',[2inf]);%生成数据矩阵%这里选用450组数据进行递推迭代n=2;N=450-n;%初始赋值如下：=[000000]’;PO=10^10*ey%循环迭代FIA=[-Y(n+i)-Y(n+i-1)U(n+i)U(n+K=P0*FIA*inv(1+FIA'*POEQ\*jc3\*hps20\o\al(\s\up9(P0),be)EQ\*jc3\*hps20\o\al(\s\up9(=),l)EQ\*jc3\*hps20\o\al(\s\up9(1:),b)disp('取450组数据进行递推夏氏辨识的结果为：’);Dxal=beta(1),Dxa2=beta(2),Dxb1=beta(3),Dxb2Dxb2=%首先将TXT中的数据装载在Matlab中：disp('从下列文件中选择所需加载文件：uy1.txtuy2.txtuy3.txt')Data=fscanf(Wbs,’%g%g',[2inf]);%生成数据矩阵%使用增广矩阵法进行2阶系统的参数辨识fori=1:3*n%首先给辨识参数向量赋初值，这里取thita0为0P0=10^10*eye(3*n,3*n);%然后设定初始状态P0为充分大的对角阵thita=[thita0,zeros(3*n,L)];%设定参数向量的初值及大小zs=zeros(L,1);%构造初始噪声矩阵forj=n+1:1:Lhn=[-Y(j-1),-Y(j-2),U(j-1),U(j-2),zs(j-1),zs(j-K=P0*hn*inv(hn'*PO*hn+1);%计算修正矩阵Kthita0=thitaO+K*[Y(j)-hn'*tP0=PO-K*hn'*P0;%计算状态矩阵P(N)E=abs((thita(:,j-1)-thita0)ifE<=ERRObreak;%根据设定的误差限，判断何时跳出EQ\*jc3\*hps24\o\al(\s\up10(end),Zgm)Uy1.txtuy2.txtuy3.txt11.极大似然法%首先将TXT中的数据装载在Matlab中：[Wbs,Message]=fopen(filename,'r');%打开文件；Data=fscanf(Wbs,%g%g,[2inf]);%生成数据矩阵mthita=zeros(3*n,1);%给被辨识参数矩阵赋初始值%用最小二乘法求初始的被辨识参数矩阵OL=[-Y(2:L-1),-Y(1:L-2),U(2:L-thita=inv(OL'*OL)*OL'*Z%得到出之后，使用牛顿—拉布森极大似然法编程%给J的梯度和Hassian矩阵赋初始值Ja1=mthita(1);Ja2=mthita(2);Jby2(i)=-Ja1*Y(i-1)-Ja2*Y(i-2)+Jb1*U(i-1)+Jb2*U(i-2)+Jc1*E(i-1)+Jc2*forj=n+1:1:lEa1(j)=Y(j-1)-[Jc1,Jc2]*[Ea1(j-1),Ea1(j-2)]';Ea2(j)=Y(j-2)-[Jc1,Jc2]Eb1(j)-U(j-1)--U(j-2)-[Jc1,Jc2]*[Eb2(j-1),EEc1(j)=-E(j-1)-[Jc1,Jc2]*[Ec1(j-1),Ec1(j-2)]';Ec2(j)=-E(j-2)-[Jc1,Jc2]Ethita=[Ea1(j),Ea2(j),fd=fd+E(j)*Ethita';%迭代计算」的梯度矩阵Hassian=Hassian+Ethita'*Ethita;%迭代计算Hassian矩阵%用牛顿-拉卜森法计算被辨识的参数的估值%跳出迭代的判断disp('极大似然法(牛顿-拉卜森方法),得到的结果为：)可用LS法起步，否则没有可靠的收敛性；对数据u(k)在z(以上题的结果为例，进行：(10分)第二部分：在递推过程中，由于是变化的故而需要不断计算；(2)辨识精度高于最小LS估计法；(4)参数估计时需构造辅助变量矩阵。第四题.系统模型阶次的辨识：(10分)1.按残差方差定阶法(含估计值误差方差最小法和F检验法)和残差白色定阶法：disp(从下列文件中选择所需加载文件：uy1.txtuy2.txtuy3.txt')