《系统辨识》大作业学号：********班级：自动化1班姓名：**信息与控制工程学院自动化系2015-07-11第一题模仿index2，搭建对象，由相关分析法，获得脉冲响应序列ˆ()g k，由ˆ()g k，参照讲义，获得系统的脉冲传递函数()G z和传递函数()G s;应用最小二乘辨识,获得脉冲响应序列ˆ()g k;同图显示两种方法的辨识效果图；应用相关最小二乘法，拟合对象的差分方程模型；构建时变对象，用最小二乘法和带遗忘因子的最小二乘法，(可以用辨识工具箱) 辨识模型的参数，比较两种方法的辨识效果差异;答：根据index2搭建结构框图：相关分析法：利用结构框图得到UY 和tout其中的U就是题目中要求得出的M序列，根据结构框图可知序列的周期是1512124=-=-=npN。

在command window中输入下列指令，既可以得到脉冲相应序列()g k：aa=5;NNPP=15;ts=2; RR=ones(15)+eye(15); for i=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) hold onstem(0:2:29,GG,'filled') Grid;title('脉冲序列g(τ)')最小二乘法建模的响应序列由于是二阶水箱系统，可以假设系统的传递函数为221101)(sa s a sb b s G +++=，已知)(τg ，求2110,,,a a b b已知G （s ）的结构，用长除法求得G(s)的s 展开式，其系数等于从 )( g 求得的各阶矩，然后求G(s)的参数。

得到结果： a1 =-1.1561 a2 =0.4283 b0 =-0.0028 b1=0.2961在command window 中输入下列指令得到传递函数：最小二乘一次算法相关参数%最小二乘法一次完成算法 M=UY(:,1); z=UY(:,2); H=zeros(100,4); for i=1:100 H(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.7866 0.1388 0.5707 0.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;%遗忘因数for i=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;while Wbs<1disp('从下列文件中选择所需加载文件：uy1.txt uy2.txt uy3.txt') filename = input('输入所需打开文件名:\','s');[Wbs,Message] = fopen(filename,'r');%打开文件；endData = fscanf(Wbs,'%g %g',[2 inf]);%生成数据矩阵Data = Data';%最小二乘法，取b0=0n = 2;%根据题目，本系统为2阶系统L = length(Data);%辨识所需数据的总长度N = L-n;%构造测量矩阵的数据长度FIA = zeros(N,2*n);%构造测量矩阵for i = 1:N %测量矩阵赋值for l = 1:n*2if(l>n)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.txt uy2.txt uy3.txt2.递推最小二乘法：%首先将TXT中的数据装载在Matlab中：clear;clc;Wbs = 0;while Wbs<1disp('从下列文件中选择所需加载文件：uy1.txt uy2.txt uy3.txt') filename = input('输入所需打开文件名:\','s');[Wbs,Message] = fopen(filename,'r');%打开文件；endData = fscanf(Wbs,'%g %g',[2 inf]);%生成数据矩阵Data = Data';%递推最小二乘法n = 2;%根据题目已知该系统阶数为2L = length(Data);%为辨识参数向量赋初值，这里均取为0.001for a = 1:1:n*2thita0(a,1) = 0.001;end%直接给出矩阵Pn的初始状态P0P0 = 10^9*eye(n*2,n*2);thita = [thita0,zeros(n*2,L)];%需辨识参数矩阵的初值及大小for i = n+1:1:L %这里i从第n+1个数开始for m = 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.txt uy2.txt uy3.txt3.辅助变量法：%首先将TXT中的数据装载在Matlab中：clear;clc;Wbs = 0;while Wbs<1disp('从下列文件中选择所需加载文件：uy1.txt uy2.txt uy3.txt')filename = input('输入所需打开文件名:\','s');[Wbs,Message] = fopen(filename,'r');%打开文件；endData = fscanf(Wbs,'%g %g',[2 inf]);%生成数据矩阵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;%得到初始参数向量%获得初值，以下使用辅助变量法for i=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.txt uy2.txt uy3.txt4.递推辅助变量法：%首先将TXT中的数据装载在Matlab中：clear;clc;Wbs = 0;while Wbs<1disp('从下列文件中选择所需加载文件：uy1.txt uy2.txt uy3.txt') filename = input('输入所需打开文件名:\','s');[Wbs,Message] = fopen(filename,'r');%打开文件；endData = fscanf(Wbs,'%g %g',[2 inf]);%生成数据矩阵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);%后面的数据计算，使用递推辅助变量法for n = 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.txt uy2.txt uy3.txt5.广义最小二乘法：%首先将TXT中的数据装载在Matlab中：clear;clc;Wbs = 0;while Wbs<1disp('从下列文件中选择所需加载文件：uy1.txt uy2.txt uy3.txt')filename = input('输入所需打开文件名:\','s');[Wbs,Message] = fopen(filename,'r');%打开文件；endData = fscanf(Wbs,'%g %g',[2 inf]);%生成数据矩阵Data = Data';%使用广义最小二乘法，辨识2阶系统的参数n = 2;L = length(Data);N = L-n;for m=1:1:300%分别取出输入、输出数据u,yU = Data(1:L,1); Y = Data(1:L,2);%首先使用最小二乘法，求系统参数向量初值gyOL = [-Y(2:L-1),-Y(1:L-2),U(2:L-1),U(1:L-2)];%LS测量矩阵Zgy1 = Data(3:L,2);%输出矩阵ythita = inv(gyOL'*gyOL)*gyOL'*Zgy1;%参数向量初值Gya1=thita(1);Gya2=thita(2);Gyb1=thita(3);Gyb2=thita(4);%得到初值后，以下使用广义最小二乘法进行辨识，由残差代替噪声误差E(1) = Y(1);%k=1时，残差的值E(2) = Y(2)+Gya1*Y(1)-Gyb1*U(1);%k=2时，残差的值for i = 3:N+n %残差公式如下：E(i)= Y(i)+Gya1*Y(i-1)+Gya2*Y(i-2)-Gyb1*U(i-1)-Gyb2*U(i-2);endE1 = E(n+1:n+N)';omiga(1:N,1) = -E(n:n+N-1)';omiga(1:N,2) = -E(n-1:n+N-2)';%fl表示由残差代替噪声项，而辨识得到的滤波器的估计参数fl=inv(omiga'*omiga)*omiga'*E1;%得到fl后，计算经过滤波的u，y数据；其中，k=1,2时：Data(1,1)=U(1);Data(1,2)=Y(1);Data(2,2)=Y(2)+fl(1)*Y(1);Data(2,1)=U(2)+fl(1)*U(1);%k>=3时：for k=3:N+nData(k,1) = U(k)+fl(1)*U(k-1)+fl(2)*U(k-2);Data(k,2) = Y(k)+fl(1)*Y(k-1)+fl(2)*Y(k-2);endenddisp('由广义最小二乘法得到的辨识结果如下：')Gya1=thita(1),Gya2=thita(2),Gyb1=thita(3),Gyb2=thita(4)实验结果：Uy1.txt uy2.txt uy3.txt6.递推广义最小二乘法：%首先将TXT中的数据装载在Matlab中：clear;clc;Wbs = 0;while Wbs<1disp('从下列文件中选择所需加载文件：uy1.txt uy2.txt uy3.txt') filename = input('输入所需打开文件名:\','s');[Wbs,Message] = fopen(filename,'r');%打开文件；endData = fscanf(Wbs,'%g %g',[2 inf]);%生成数据矩阵Data = Data';u1=Data(:,1);y1=Data(:,2);n=2;N=500-n;n0=50;%用前50组数据利用lsm估计出初值%构造矩阵phiphi(:,1)=-y1(n:n+n0-1);phi(:,2)=-y1(n-1:n+n0-2);phi(:,3)=u1(n:n+n0-1);phi(:,4)=u1(n-1:n+n0-2);y(:,1)=y1(n+1:n+n0);seta=(inv(phi'*phi))*phi'*y;seta0=seta;f0=[0 0]';p10=inv(phi'*phi);p20=10^6*eye(2,2);y(1:i,1)=y1(n+1:n+i);u(1:i,1)=u1(n+1:n+i);%计算残差e=y-phi*seta0;%f 估计omega=[-e(n+i-2) -e(i+1-1)]';K2=p20*omega*inv(1+omega'*p20*omega);f1=f0+K2*(e(n+i-2)-omega'*f0);p2=p20-K2*omega'*p20;p20=p2;f0=f1;%滤波yy(1:i,1)=y1(n:n+i-1);yy(1:i,2)=y1(n-1:n+i-2);uu(1:i,1)=u1(n:(n+i-1));uu(1:i,2)=u1((n-1):(n+i-2));yg=y+yy*f1; %yg(k)=y(k)+f(1)*y(k-1)+f(2)*y(k-2)ug=u+uu*f1; %ug(k)=u(k)+f(1)*u(k-1)+f(2)*u(k-2)% seta 估计ksai=[-yg(n+i-2) -yg(n+i-3) ug(n+i-2) ug(n+i-3)]';K1=p10*ksai*inv(1+ksai'*p10*ksai);% K(N+1)值p1=p10-K1*ksai'*p10; % P(N+1)值p10=p1;seta1=seta0+K1*(y1(n+i+1)-ksai'*seta0); % seta(N+1)值seta0=seta1;% 构造新phi 阵phi(1:i+1,1)=-y1(n:n+i-1+1);phi(1:i+1,2)=-y1(n-1:n+i-2+1);phi(1:i+1,3)=u1(n:n+i-1+1);phi(1:i+1,4)=u1(n-1:n+i-2+1);enddisp('递推广义最小二乘法辨识结果：');a1=seta0(1),a2=seta0(2),b1=seta0(3),b2=seta0(4)实验结果：Uy1.txt uy2.txt uy3.txt7.夏氏偏差修正法：%首先将TXT中的数据装载在Matlab中：clear;clc;Wbs = 0;while Wbs<1disp('从下列文件中选择所需加载文件：uy1.txt uy2.txt uy3.txt') filename = input('输入所需打开文件名:\','s');[Wbs,Message] = fopen(filename,'r');%打开文件；endData = fscanf(Wbs,'%g %g',[2 inf]);%生成数据矩阵Data = Data';%使用夏氏修正法，对2阶系统进行参数辨识n = 2;L = length(Data);N = L-n;U = Data(:,1);Y = Data(:,2);%首先使用最小二乘法，求系统参数向量初值glOL =[-Y(2:L-1),-Y(1:L-2),U(2:L-1),U(1:L-2)];Zgl1 = Data(3:L,2);Sgl1 = glOL'*glOL;Sgl2=inv(Sgl1);Sgl3=glOL'*Zgl1;Xsthita = Sgl2*Sgl3;%计算参数矩阵%得到初值后，以下使用夏氏修正法进行参数辨识：thitab0 = 0;%设偏差项的偏差初值为0Fa = Sgl2*glOL';M = eye(N)-glOL*Sgl2*glOL';F = eye(N);if(F>=10^-6*eye(N))E1 = Zgl1-glOL*Xsthita;%计算残差Eomiga(2:N,1) = -E1(1:N-1);omiga(3:N,2) = -E1(1:N-2);D = omiga'*M*omiga;Fx = inv(D)*omiga'*M*Zgl1;thitab = Fa*omiga*Fx;Xsthita = Xsthita - thitab;F = thitab - thitab0;thitab0 = thitab;enddisp('夏氏修正法')Xsa1 = Xsthita(1),Xsa2 = Xsthita(2),Xsb1 = Xsthita(3),Xsb2 = Xsthita(4)实验结果：Uy1.txt uy2.txt uy3.txt%首先将TXT中的数据装载在Matlab中：clear;clc;Wbs = 0;while Wbs<1disp('从下列文件中选择所需加载文件：uy1.txt uy2.txt uy3.txt') filename = input('输入所需打开文件名:\','s');[Wbs,Message] = fopen(filename,'r');%打开文件；endData = fscanf(Wbs,'%g %g',[2 inf]);%生成数据矩阵Data = Data';%最小二乘法，取b0=0n = 2;%根据题目，本系统为2阶系统L = length(Data);%辨识所需数据的总长度N = L-n;%构造测量矩阵的数据长度FIA = zeros(N,2*n);%构造测量矩阵for i = 1:N %测量矩阵赋值for l = 1:n*2if(l>n)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;%计算参数矩阵thitaLS = thita;%最小二乘估值m = inv(FIA'*FIA)*FIA';for i = 1:1:1000%构造OmegaErr = Y-FIA*thitaLS;%计算残差Eomiga(2:N,1) = -Err(1:N-1);omiga(3:N,2) = -Err(1:N-2);F = inv(omiga'*omiga)*omiga'*Err;thita = thitaLS-m*omiga*F;enddisp('使用夏氏改良法辨识结果为：')Xga1 = thita(1),Xga2 = thita(2),Xgb1 = thita(3),Xgb2 = thita(4)实验结果：Uy1.txt uy2.txt uy3.txt%首先将TXT中的数据装载在Matlab中：clear;clc;Wbs = 0;while Wbs<1disp('从下列文件中选择所需加载文件：uy1.txt uy2.txt uy3.txt')filename = input('输入所需打开文件名:\','s');[Wbs,Message] = fopen(filename,'r');%打开文件；endData = fscanf(Wbs,'%g %g',[2 inf]);%生成数据矩阵Data = Data';%这里选用450组数据进行递推迭代n = 2;N = 450-n;U = Data(:,1);Y = Data(:,2);%初始赋值如下：beta = [0 0 0 0 0 0]';P0= 10^10*eye(6,6);E1 = 0;E2 = 0;%循环迭代for i = 3:NE1 = Y(n+i)-beta(1)*Y(n+i-1)-beta(2)*Y(n+i-2)-beta(3)*Y(n+i-1)-beta(4)*U(n+i-2);E2 = Y(n+i-1)-beta(1)*Y(n+i-2)-beta(2)*Y(n+i-3)-beta(3)*U(n+i-2)-beta(4)*U(n+i-3); FIA = [-Y(n+i) -Y(n+i-1) U(n+i) U(n+i-1) -E1 -E2]';K = P0*FIA*inv(1+FIA'*P0*FIA);P1 = P0-K*FIA'*P0;P0 = P1;beta1 = beta+K*(Y(n+i+1)-FIA'*beta);beta = beta1;enddisp('取450组数据进行递推夏氏辨识的结果为:');Dxa1=beta(1),Dxa2=beta(2),Dxb1=beta(3),Dxb2=beta(4)实验结果：Uy1.txt uy2.txt uy3.txt10.增广矩阵法：%首先将TXT中的数据装载在Matlab中：clear;clc;Wbs = 0;while Wbs<1disp('从下列文件中选择所需加载文件：uy1.txt uy2.txt uy3.txt') filename = input('输入所需打开文件名:\','s');[Wbs,Message] = fopen(filename,'r');%打开文件；endData = fscanf(Wbs,'%g %g',[2 inf]);%生成数据矩阵Data = Data';%使用增广矩阵法进行2阶系统的参数辨识n = 2;L = length(Data);U = Data(1:L,1);Y = Data(1:L,2);for i = 1:3*n %首先给辨识参数向量赋初值，这里取thita0为0thita0(i,1) = 0;endP0 = 10^10*eye(3*n,3*n);%然后设定初始状态P0为充分大的对角阵ERR0 = 0.00000000001;%收敛情况指标，相对误差ERR0thita = [thita0,zeros(3*n,L)];%设定参数向量的初值及大小zs = zeros(L,1);%构造初始噪声矩阵hn = [0,0,0,0,0,0]';for j = n+1:1:Lzs(j-1) = Y(j-1)-hn'*thita0;hn = [-Y(j-1),-Y(j-2),U(j-1),U(j-2),zs(j-1),zs(j-2)]';K = P0*hn*inv(hn'*P0*hn+1);%计算修正矩阵Kthita0 = thita0+K*[Y(j)-hn'*thita0];P0 = P0-K*hn'*P0;%计算状态矩阵P(N)thita(:,j) = thita0;E = abs((thita(:,j-1)-thita0)./thita0);if E<=ERR0 break;%根据设定的误差限，判断何时跳出endendZgma1=thita0(1),Zgma2=thita0(2),Zgmb1=thita0(3),Zgmb2=thita0(4), Zgmc1=thita0(5),Zgmc2=thita0(6)Uy1.txt uy2.txt uy3.txt11.极大似然法：%首先将TXT中的数据装载在Matlab中：clear;clc;Wbs = 0;while Wbs<1disp('从下列文件中选择所需加载文件：uy1.txt uy2.txt uy3.txt')filename = input('输入所需打开文件名:\','s');[Wbs,Message] = fopen(filename,'r');%打开文件；endData = fscanf(Wbs,'%g %g',[2 inf]);%生成数据矩阵Data = Data';n = 2;L = length(Data);N = L-n;U = Data(:,1);Y = Data(:,2);mthita = zeros(3*n,1);%给被辨识参数矩阵赋初始值%用最小二乘法求初始的被辨识参数矩阵OL = [-Y(2:L-1),-Y(1:L-2),U(2:L-1),U(1:L-2)];Z1 = Data(3:L,2);thita = inv(OL'*OL)*OL'*Z1;mthita(1:2*n,1) = thita;%得到出之后，使用牛顿—拉布森极大似然法编程E = zeros(L,1);J = 0;sgma0 = 1;Ea1 = zeros(L,1);Ea2 = zeros(L,1);Eb1 = zeros(L,1);Eb2 = zeros(L,1);Ec1 = zeros(L,1);Ec2 = zeros(L,1);%给J的梯度和Hassian矩阵赋初始值Ethita = zeros(3*n,1);fd = zeros(3*n,1);Hassian = zeros(3*n,3*n);for p = 1:500Ja1 = mthita(1);Ja2 = mthita(2);Jb1 = mthita(3);Jb2 = mthita(4);Jc1 = mthita(5);Jc2 = mthita(6);for i = n+1:Ly2(i) = -Ja1*Y(i-1)-Ja2*Y(i-2)+Jb1*U(i-1)+Jb2*U(i-2)+Jc1*E(i-1)+Jc2*E(i-2);E(i) = Y(i)-y2(i);J = 0.5*(E(i)^2);sgma = (1/N)*(E(i)^2);endfor j = n+1:1:LEa1(j) = Y(j-1)-[Jc1,Jc2]*[Ea1(j-1),Ea1(j-2)]'; Ea2(j) = Y(j-2)-[Jc1,Jc2]*[Ea2(j-1),Ea2(j-2)]';Eb1(j) = -U(j-1)-[Jc1,Jc2]*[Eb1(j-1),Eb1(j-2)]';Eb2(j) = -U(j-2)-[Jc1,Jc2]*[Eb2(j-1),Eb2(j-2)]';Ec1(j) = -E(j-1)-[Jc1,Jc2]*[Ec1(j-1),Ec1(j-2)]';Ec2(j) = -E(j-2)-[Jc1,Jc2]*[Ec2(j-1),Ec2(j-2)]';Ethita = [Ea1(j),Ea2(j),Eb1(j),Eb2(j),Ec1(j),Ec2(j)];fd = fd+E(j)*Ethita';%迭代计算J的梯度矩阵Hassian = Hassian + Ethita'*Ethita;%迭代计算Hassian矩阵end%用牛顿-拉卜森法计算被辨识的参数的估值mthita = mthita - inv(Hassian)*fd;JFIA = (sgma-sgma0)/sgma0;sgma0 = sgma;%跳出迭代的判断if JFIA<0.000000001 breakendenddisp('极大似然法（牛顿-拉卜森方法），得到的结果为：')Ja1 = mthita(1),Ja2 = mthita(2),Jb1 = mthita(3),Jb2 = mthita(4),Jc1 = mthita(5),Jc2 = mthita(6)实验结果：Uy1.txt uy2.txt uy3.txt12.分析噪声)(k 特性:RLS 法：当噪声比较小时，辨识精度较高；当噪声比较大时，收敛点可能不唯一，参数估计值往往是有偏的；但是运用此方法时，数据要充分多，否则辨识精度明显下降；另外噪声模型阶次不宜过高，也会影响精度，所以误差的原因主要是数据量的大小以及噪声模型阶次。