系统辨识大作业1201张青

《系统辨识》大作业信息与控制工程学院自动化系2015-07-11学号：12051124班级：自动化1班姓名：张青第一题模仿index2，搭建对象，由相关分析法，获得脉冲响应序列ˆ()gk，由ˆ()gk，参照讲义，获得系统的脉冲传递函数()Gz和传递函数()Gs;应用最小二乘辨识,获得脉冲响应序列ˆ()gk;同图显示两种方法的辨识效果图；应用相关最小二乘法，拟合对象的差分方程模型；构建时变对象，用最小二乘法和带遗忘因子的最小二乘法，(可以用辨识工具箱)辨识模型的参数，比较两种方法的辨识效果差异;答：根据index2搭建结构框图：相关分析法：利用结构框图得到UY和tout其中的U就是题目中要求得出的M序列，根据结构框图可知序列的周期是1512124npN。在commandwindow中输入下列指令，既可以得到脉冲相应序列()gk：aa=5;NNPP=15;ts=2;RR=ones(15)+eye(15);fori=15:-1:1UU(16-i,:)=UY(16+i:30+i,1)';endYY=[UY(31:45,2)];GG=RR*UU*YY/[aa*aa*(NNPP+1)*ts];plot(0:2:29,GG)holdonstem(0:2:29,GG,'filled')Grid;title('脉冲序列g(τ)')最小二乘法建模的响应序列由于是二阶水箱系统，可以假设系统的传递函数为221101)(sasasbbsG，已知)(g，求2110,,,aabb已知G（s）的结构，用长除法求得G(s)的s展开式，其系数等于从)(g求得的各阶矩，然后求G(s)的参数。得到结果：a1=-1.1561a2=0.4283b0=-0.0028b1=0.2961在commandwindow中输入下列指令得到传递函数：最小二乘一次算法相关参数%最小二乘法一次完成算法M=UY(:,1);z=UY(:,2);H=zeros(100,4);fori=1:100H(i,1)=-z(i+1);H(i,2)=-z(i);H(i,3)=M(i+1);H(i,4)=M(i);endEstimate=inv(H'*H)*H'*(z(3:102))%结束得到相关系数为：Estimate=-0.78660.13880.57070.3115带遗忘因子最小二乘法：%带遗忘因子最小二乘法程序M=UY(:,1);z=UY(:,2);P=1000*eye(5);%Theta=zeros(5,200);%Theta(:,1)=[0;0;0;0;0];K=zeros(4,400);%K=[10;10;10;10;10];lamda=0.99;%遗忘因数fori=3:201h=[-z(i-1);-z(i-2);M(i);M(i-1);M(i-2)];K=P*h*inv(h'*P*h+lamda);Theta(:,i-1)=Theta(:,i-2)+K*(z(i)-h'*Theta(:,i-2));P=(eye(5)-K*h')*P/lamda;endi=1:200;figure(1)plot(i,Theta(1,:),i,Theta(2,:),i,Theta(3,:),i,Theta(4,:),i,Theta(5,:))title('带遗忘因子最小二乘法')grid%结束Estimate可由仿真图得出，可知两种方法参数确定十分接近。两种方法的辨识效果差异:采用相关分析法和最小二乘法辨识出的脉冲响应图形可看出，用最小二乘法辨识出的图形更像脉冲响应，他在最后无偏差，衰减为零，其可实现无偏差估计。而相关分析法图形虽与相关分析的相近，但它最后有偏差。带遗忘因子的递推最小二乘辨识的参数平均值可随着实际参数变化而突变。但他对噪声比较敏感，只是参数围绕参数实际值上下波动，而辨识出参数的平均值接近实际值，而且比其他方法更加接近，并且可以防止数据饱和。第二题设SISO系统差分方程为（40分）y(k)=)()2()1()2()1(2121kkubkubkyakya辨识参数向量为＝1[a2a1b2b]T，输入输出数据详见数据文件uy01.txt—uy03.txt。)(k为噪声方差各异的白噪声或有色噪声。试求解：1）用最小二乘及递推最小二乘法估计；2）用辅助变量法及其递推算法估计；3）用广义最小二乘法及其递推算法估计；4）用夏氏偏差修正法、夏氏改良法及其递推算法估计；5）用增广矩阵法估计；6)用极大似然法估计；7）分析噪声)(k特性；答：1.基本最小二乘法：%首先将TXT中的数据装载在Matlab中：clear;clc;Wbs=0;whileWbs1disp('从下列文件中选择所需加载文件：uy1.txtuy2.txtuy3.txt')filename=input('输入所需打开文件名:\','s');[Wbs,Message]=fopen(filename,'r');%打开文件；endData=fscanf(Wbs,'%g%g',[2inf]);%生成数据矩阵Data=Data';%最小二乘法，取b0=0n=2;%根据题目，本系统为2阶系统L=length(Data);%辨识所需数据的总长度N=L-n;%构造测量矩阵的数据长度FIA=zeros(N,2*n);%构造测量矩阵fori=1:N%测量矩阵赋值forl=1:n*2if(ln)FIA(i,l)=Data(n+2+i-l,1);elseFIA(i,l)=-Data(n+i-l,2);endendendY=Data(n+1:n+N,2);%输出数据矩阵thita=inv(FIA'*FIA)*FIA'*Y;%计算参数矩阵disp('使用最小二乘算法所得结果如下：')%辨识参数数据输出如下：lsa1=thita(1),lsa2=thita(2),lsb1=thita(3),lsb2=thita(4)实验结果：Uy1.txtuy2.txtuy3.txt2.递推最小二乘法：%首先将TXT中的数据装载在Matlab中：clear;clc;Wbs=0;whileWbs1disp('从下列文件中选择所需加载文件：uy1.txtuy2.txtuy3.txt')filename=input('输入所需打开文件名:\','s');[Wbs,Message]=fopen(filename,'r');%打开文件；endData=fscanf(Wbs,'%g%g',[2inf]);%生成数据矩阵Data=Data';%递推最小二乘法n=2;%根据题目已知该系统阶数为2L=length(Data);%为辨识参数向量赋初值，这里均取为0.001fora=1:1:n*2thita0(a,1)=0.001;end%直接给出矩阵Pn的初始状态P0P0=10^9*eye(n*2,n*2);thita=[thita0,zeros(n*2,L)];%需辨识参数矩阵的初值及大小fori=n+1:1:L%这里i从第n+1个数开始form=1:n*2%输入矩阵赋值if(m=n)FIA1(m,1)=-Data(i-m,2);elseFIA1(m,1)=Data(i-m+2,1);endendX=FIA1'*P0*FIA1+1;%引入中间变量X为K递推公式中的求逆项,1维K1=P0*FIA1/X;D1=Data(i,2)-FIA1'*thita0;P1=P0-K1*K1'*X;%求出P(K)矩阵P0=P1;%作为下次迭代初值thita1=thita0+K1*D1;%估值补偿公式，thita1为求被辨识的参数向量thita0=thita1;%所得的参数向量作为下次迭代初值enddisp('递推最小二乘法所得结果如下：')Dta1=thita1(1),Dta2=thita1(2),Dtb1=thita1(3),Dtb2=thita1(4)%显示被辨识参数实验结果：Uy1.txtuy2.txtuy3.txt3.辅助变量法：%首先将TXT中的数据装载在Matlab中：clear;clc;Wbs=0;whileWbs1disp('从下列文件中选择所需加载文件：uy1.txtuy2.txtuy3.txt')filename=input('输入所需打开文件名:\','s');[Wbs,Message]=fopen(filename,'r');%打开文件；endData=fscanf(Wbs,'%g%g',[2inf]);%生成数据矩阵Data=Data';%使用递推辅助变量法对2阶系统进行辨识，取b0=0n=2;L=length(Data);N=L-n;%分别取出输入、输出数据u,yU=Data(1:L,1);Y=Data(1:L,2);%首先使用最小二乘法，粗略得到初始的被辨识参数向量fzOL=[-Y(2:L-1),-Y(1:L-2),U(2:L-1),U(1:L-2)];Y1=Data(3:L,2);%构造输出向量，也用于辅助变量法的输出向量thita=inv(fzOL'*fzOL)*fzOL'*Y1;%得到初始参数向量%获得初值，以下使用辅助变量法fori=1:400Yf=fzOL*thita;%辅助模型输出%构造辅助变量矩阵，改变新的矩阵中的输出量，输入量不变Zfz=fzOL;Zfz(2:N,1)=-Yf(1:(N-1),1);Zfz(3:N,2)=-Yf(1:(N-2),1);thita=inv(Zfz'*fzOL)*Zfz'*Y1;%求取被估计参数的辅助变量估值enddisp('使用辅助变量法，得到如下辨识结果：')Fza1=thita(1),Fza2=thita(2),Fzb1=thita(3),Fzb2=thita(4)实验结果：Uy1.txtuy2.txtuy3.txt4.递推辅助变量法：%首先将TXT中的数据装载在Matlab中：clear;clc;Wbs=0;whileWbs1disp('从下列文件中选择所需加载文件：uy1.txtuy2.txtuy3.txt')filename=input('输入所需打开文件名:\','s');[Wbs,Message]=fopen(filename,'r');%打开文件；endData=fscanf(Wbs,'%g%g',[2inf]);%生成数据矩阵Data=Data';%使用递推辅助变量法对2阶系统进行辨识，取b0=0n=2;L=length(Data);%分别取出输入、输出数据u,yU=Data(1:L,1);Y=Data(1:L,2);N=498;n0=100;%取前100个数据利用最小二乘辨识FIA=[-Y(2:n0-1),-Y(1:n0-2),U(2:n0-1),U(1:n0-2)];Y1=Data(3:n0,2);%构造输出向量，也用于辅助变量法的输出向量thita=inv(FIA'*FIA)*FIA'*Y1;%得到初始参数向量Z=[-Y(2:n0-1)-Y(1:n0-2)U(2:n0-1)U(1:n0-2)];thita=inv(Z'*FIA)*Z'*Y1;P0=inv(Z'*FIA);%后面的数据计算，使用递推辅助变量法forn=n0+1:NY11=[Y(1);Y(2);Z*thita];Zn=[-Y11(n-1)-Y11(n-2)U(n-1)U(n-2)];Z=[Z;Zn];kesy=[-Y(n-1);-Y(n-2);U(n-1);U(n-2)];K=P0*Zn'/(1+kesy'*P0*Zn');P1=P0-K*kesy'*P0;P0=P1;thita=thita+K*(Y(n)-kesy'*thita);enddisp('递推辅助变量法辨识结果：');Dfa1=thita(1),Dfa2=thita(2),Dfb1=thita(3),Dfb2=thita(4)实验结果：Uy1.txtuy2.txtuy3.txt5.广义最小二乘法：%首先将TXT中的数据装载在Matlab中：clear;clc;Wbs=0;whileWbs1disp('从下列文件中选择