系统的辨识大作业

1、实用标准文案系统辨识大作业学号:12051124班级:自动化1班精彩文档张青姓名:信息与控制工程学院自动化系2015-07-11第一题模仿index2 ,搭建对象，由相关分析法，获得脉冲响应序列g?(k),由於k),参照讲义，获得系统的脉冲传递函数 G(z)和传递函数G(s);应用最小二乘辨识，获得脉冲响应序 列(k)；同图显示两种方法的辨识效果图；应用相关最小二乘法，拟合对象的差分方程模型；构建时变对象，用最小二乘法和带遗忘因子的最小二乘法，(可以用辨识工具箱)辨识模型的参数，比较两种方法的辨识效果差异；答：根据index2搭建结构框图：相关分析法：利用结构框图得到 UY和tout其中的U就

2、是题目中要求得出的M序列，根据结构框图可知序列的周期是 _ n_ 4Np =2 -1 =2 1=15 。 p在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( t )脉冲序列gE 04ii-0

3、05-liJJ0510152630最小二乘法建模的响应序列一. b-bs由于是二阶水箱系统，可以假设系统的传递函数为G(s)= b0 b1s 2 ,已知g(T),求1als a2sb0,ha,a2已知G (s)的结构，用长除法求得 G(s)的s展开式，其系数等于从g(E)求得的各阶矩,然后求G(s)的参数。得到结果：al =-1.1561a2 =0.4283b0 =-0.0028 b1 =0.2961在command window中输入下列指令得到传递函数： num=0,2961 -0. 0028 ; den= 0. 4283 -1. 1561 1 ; sys=tf (num, den)Ira

4、nsfer funct ion:0.2961 s - 0* 002S0.4283 s 2 - L 156 s + 1最小二乘一次算法相关参数%最小二乘法一次完成算法M=UY(:,1);z=UY(:,2);H=zeros(100,4);for i=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(:

5、,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; %d1忘因数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

6、(2,:),i,Theta(3,:),i,Theta(4,:),i,Theta(5,:)title(带遗忘因子最小二乘法)grid%吉束带遗忘因子最小二乘 0.60.4020+ 2440.6-0BEstimate可由仿真图得出，可知两种方法参数确定十分接近。两种方法的辨识效果差异：采用相关分析法和最小二乘法辨识出的脉冲响应图形可看出，用最小二乘法辨识出的图形 更像脉冲响应，他在最后无偏差，衰减为零，其可实现无偏差估计。而相关分析法图形虽 与相关分析的相近，但它最后有偏差。带遗忘因子的递推最小二乘辨识的参数平均值可随 着实际参数变化而突变。但他对噪声比较敏感，只是参数围绕参数实际值上下波动，而辨

7、 识出参数的平均值接近实际值，而且比其他方法更加接近，并且可以防止数据饱和。第二题设SISO系统差分方程为(40分)y(k)=-a1y(k -1) -a2y(k -2) b1u(k -1) b2u(k -2)(k)辨识参数向量为日=aia2b b2 T,输入输出数据详见数据文件uy01.txt -uy03.txt 。 .k)为噪声方差各异的白噪声或有色噪声。试求解：1 )用最小二乘及递推最小二乘法估计9;2 )用辅助变量法及其递推算法估计 H ;3 )用广义最小二乘法及其递推算法估计9;4 )用夏氏偏差修正法、夏氏改良法及其递推算法估计9 ;5)用增广矩阵法估计9 ;6)用极大似然法估计9 ;

8、7)分析噪声:(k)特性；答：1.基本最小二乘法：喻先将TXT中的数据装载在 Matlab中： clear;clc;Wbs = 0;while Wbsn)FIA(i,l) = Data(n+2+i-l,1);elseFIA(i,l) = -Data(n+i-l,2);endendend丫 = 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

9、)实验结果：Uy1.txtuy2.txtuy3.txt sal =Isal =Isal =-h 4555-1. 4579-1.4693lsa2 -lsa2 =0. S5190. 85130.3639Isbl =Isbl =Isbl =0. 42500. 7344口.7451lsb2 工lsb2 =lsb2 =0. 56190. 57730.55432.递推最小二乘法：痂先将TXT中的数据装载在Matlab中：clear;clc;Wbs = 0;while Wbs1disp(从下列文件中选择所需加载文件：uy1.txt uy2.txt uy3.txt)filename = input(输入所需打

10、开文件名：,s);Wbs,Message = fopen(filename,r);%打开文件；endData = fscanf(Wbs,%g %g,2 inf);%生成数据矩阵Data = Data;%M推最小二乘法n = 2;%根据题目已知该系统阶数为2L = length(Data);私辨识参数向量赋初值，这里均取为0.001for a = 1:1:n*2thita0(a,1) = 0.001; end帆接给出矩阵Pn的初始状态P0P0 = 10A9*eye(n*2,n*2);thita = thita0,zeros(n*2,L);%需辨识参数矩阵的初值及大小for i = n+1:1:L

11、 % 这里i从第n+1个数开始for m = 1:n*2 % 输入矩阵赋值 if(mtb2 =Dtb2 =0. 56190. 57730. 55433.辅助变量法：砥先将TXT中的数据装载在 Matlab中：clear;clc;Wbs = 0;while Wbs别取出输入、输出数据 u,yU = Data(1:L,1);Y = Data(1:L,2);%T先使用最小二乘法，粗略得到初始的被辨识参数向量fzOL =卜丫(2:L-1),-Y(1:L-2),U(2:L-1),U(1:L-2);Y1 = Data(3:L,2);%构造输出向量，也用于辅助变量法的输出向量thita = inv(fzOL

12、*fzOL)*fzOL*Y1;%得到初始参数向量哪得初值，以下使用辅助变量法 fo门=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.

13、txt uy2.txt uy3.txtFzal =FM=-1.4755-1. 4479Fza2 =Fza2二0, S693。，8炎7Fzbl -Fzbl-43040. 7302Fzb2 =Fzb2=Fzal =-1.4662Fza2 =0.8620Fzbl =0. 7441Fzb2 =0.55540. 55890. 53214.递推辅助变量法：喻先将TXT中的数据装载在 Matlab中：clear;clc;Wbs = 0;while Wbs别取出输入、输出数据u,yU = Data(1:L,1);Y = Data(1:L,2);N = 498;n0 = 100;%取前100个数据利用最小二乘辨

14、识FIA =卜丫(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 =卜丫11(n-1) -Y11(n-2) U(n-1) U(n

15、-2);Z = Z; Zn;kesy =卜Y(n-1); -Y(n-2); U(n-1); U(n-2);K = P0*Zn7(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.txtQfal =-1.467 t=0.8624Dfbl =0.4291Ufb2 =0.5552Df al

16、=4657Dfa2 =0. 8625Dfbl =0. 7466Dfb2 =0.56SBIf al =-L4617Ifa2 =0.8559Dfb=Q,74MSfb2 =Q.E6385.广义最小二乘法：喻先将TXT中的数据装载在 Matlab中：clear;clc;Wbs = 0;while Wbs别取出输入、输出数据u,yU = Data(1:L,1); Y = Data(1:L,2);%T先使用最小二乘法，求系统参数向量初值gyOL = -Y(2:L-1),-Y(1:L-2),U(2:L-1),U(1:L-2);%LS测量矩阵Zgy1 = Data(3:L,2);%输出矩阵 ythita =

17、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) =

18、-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-

19、1)+fl(2)*Y(k-2);endenddisp(由广义最小二乘法得到的辨识结果如下：)Gya1=thita(1),Gya2=thita(2),Gyb1=thita(3),Gyb2=thita(4)实验结果：实验结果:Uy1.txt uy2.txtGyal =-L4694Gya2 =XSB47Gybl =0. 4288uy3.txtGy al =-I. 4669Gya2 =H8565Gyb L =0,7415Gy al =-1. 4630Gya2 二0.8607Gybl =0.7425。油2 =Gyb2 =Gyb2 =0.55490.5638Q“ 55536 .递推广义最小二乘法：喻先将T

20、XT中的数据装载在 Matlab中：clear;clc;Wbs = 0;while Wbs1disp( 从下列文件中选择所需加载文件： 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:

21、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=10A6*eye(2,2);for i=n0:N-2y(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+

22、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*fi;%yg(k尸y(k)+f*y(k-i)+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

23、-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(递推广义最

24、小二乘法辨识结果：)；a1=seta0,a2=seta0(2),b1=seta0(3),b2=seta0(4)实验结果：Uy1.txt uy2.txta1 =-0.7976a2 =-0.1941b1 =-0.0073b2 =0.3683uy3.txta1 =-0.7872 a2 =-0.2300 b1 =0.1577 b2 =0.5090a1 =-0.8052 a2 =-0.1857 b1 =0.1552 b2 =0.49127 .夏氏偏差修正法：喻先将TXT中的数据装载在Matlab中：clear;clc;Wbs = 0;while Wbs=10A-6*eye(N)E1 = Zgl1-glO

25、L*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 - thitabO;thitabO = thitab;enddisp(夏氏修正法)Xsa1 = Xsthita(1),Xsa2 = Xsthita(2),Xsb1 = Xsthita(3),Xsb2 = Xsthita(4)实验结果：Uy1.txt uy2.

26、txtXsa1 =-1.4789Xsa2 =0.8659Xsb1 =0.4259Xsb2 =0.5553uy3.txtXsa1 = -1.4680Xsa2 = 0.8548Xsb1 = 0.7392Xsb2 = 0.5712Xsa1 =-1.4651Xsa2 =0.8603Xsb1 =0.7423Xsb2 =0.55628 .夏氏改良法：喻先将TXT中的数据装载在Matlab中：clear;clc;Wbs = 0;while Wbsn)FIA(i,l) = Data(n+2+i-l,1); elseFIA(i,l) = -Data(n+i-l,2); endendendY = Data(n+1

27、: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

28、 = thita(2),Xgb1 = thita(3),Xgb2 = thita(4)Uy1.txtuy2.txtuy3.txtXga1 = -1.4779Xga2 = 0.8654Xgb1 = 0.4258Xgb2 = 0.5557Xga1 = -1.4674Xga2 = 0.8544Xgb1 = 0.7390Xgb2 =0.5716实验结果:Xga1 = -1.4651Xga2 = 0.8603Xgb1 = 0.7424Xgb2 =0.55629 .夏氏递推法：喻先将TXT中的数据装载在 Matlab中：clear;clc;Wbs = 0;while Wbs1disp(从下列文件中选择所需

29、加载文件：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);%U始赋值如下：beta = 0 0 0 0 0 0;P0= 10A10*eye(6,6);E1 = 0;E2 = 0;%1环迭代for i = 3:NE1Y(n+i)-beta

30、(1)*Y(n+i-1)-beta(2)*Y(n+i-2)-beta(3)*Y(n+i-1)-beta(4)*U(n+i-2);E2Y(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;betal = beta+K*(Y(n+i+1)-FIA*beta);beta = betal;enddisp(

31、取450组数据进行递推夏氏辨识的结果为：)；Dxa1=beta,Dxa2=beta(2),Dxb1=beta(3),Dxb2=beta(4)实验结果：Uy1.txt uy2.txtuy3.txtDxa1 =-1.4067Dxa2 =0.8049Dxb1 =0.4300Dxb2 =0.5650Dxa1 =-1.5289Dxa2 =0.9451Dxb1 =0.7493Dxb2 =0.5432Dxa1 =-1.4919Dxa2 =0.8810Dxb1 =0.7385Dxb2 =0.541010 .增广矩阵法：喻先将TXT中的数据装载在 Matlab中：clear;clc;Wbs = 0;while

32、Wbs1disp(从下列文件中选择所需加载文件：uy1.txt uy2.txt uy3.txt)filename = input(输入所需打开文件名：,s);Wbs,Message = fopen(filename,r);%打开文件；endData = fscanf(Wbs,%g %g,2 inf);%生成数据矩阵Data = Data;%使用增广矩B法进行2阶系统的参数辨识n = 2;L = length(Data);U = Data(1:L,1);Y = Data(1:L,2);for i = 1:3*n %首先给辨识参数向量赋初值，这里取 thitaO 为0thita0(i,1) = 0

33、;endP0 = 10A10*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);%计算修正矩阵 Kthit

34、a0 = 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.txtuy2.txtuy3.txtZgmal = -1.4683Zgma2 = 0.8626Zg

35、mbl =0.4246Zgmb2 = 0.5540Zgmcl = -0.9205Zgmc2 =0.0436Zgmal = -1.4629Zgma2 = 0.8520Zgmb1 = 0.7396Zgmb2 = 0.5751Zgmc1 = 0.0282Zgmc2 = -0.2689Zgma1 = -1.4573Zgma2 = 0.8543Zgmb1 = 0.7427Zgmb2 = 0.5563Zgmc1 = 0.5183Zgmc2 =0.280611 .极大似然法： 喻先将TXT中的数据装载在 Matlab中： clear;clc;Wbs = 0;while Wbs1disp(从下列文件中选择所需

36、加载文件：uy1.txt uy2.txt uy3.txt)filename = input(输入所需打开文件名：,s);Wbs,Message = fopen(filename,r);%打开文件；end Data = 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 =卜丫(2:L-1),-Y(1:L-2),U(2:L-1),U(1

37、: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 = z

38、eros(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)A2);sgma = (1/N)*(E(i)A2);endfor j = n+1:1:LEa1(j)=Y(j-1)-Jc1,J

39、c2*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的

40、梯度矩阵Hassian = Hassian + Ethita*Ethita;%迭代计算 Hassian 矩阵end 卿牛顿-拉卜森法计算被辨识的参数的估值mthita = mthita - inv(Hassian)*fd;JFIA = (sgma-sgma0)/sgma0;sgma0 = sgma;啕K出迭代的判断if JFIA0.000000001 breakend,得到的结果为：)=mthita(3),Jb2 = mthita(4),Jc1 = mthita(5),Jc2enddisp(极大似然法(牛顿-拉卜森方法)=mthita(2),Jb1Ja1 = mthita,Ja2 =mthit

41、a(6)Uy1.txtuy2.txtuy3.txtJa1-1.4789Ja1-1.4680Ja20.8659Ja20.8548Jb10.4259Jb10.7392Jb20.5553Jb20.5712Jc1-0.9669Jc1-0.0411Jc2-0.4682Jc2-0.3384实验结果:Ja1-1.4651Ja20.8603Jb10.7423Jb20.5562Jc10.5206Jc2-0.069412.分析噪声.k)特性：RLS法：当噪声比较小时，辨识精度较高；当噪声比较大时，收敛点可能不唯一，参数估计值往往是有偏的； 但是运用此方法时， 数据要充分多，否则辨识精度明显 下降；另外噪声模型阶次

42、不宜过高，也会影响精度，所以误差的原因主要是数据量的大小以 及噪声模型阶次。IV法：辨识效果较好，能适应较大的范围的噪声特性。对初值P(0) 敏感，可用LS法起步，否则没有可靠的收敛性；对数据 u(k) 在 z(k) 的直流分量敏感，对 z(k)的直流分量不敏感。所以误差的主要原因是初值P(0)以及数据计算过程产生的误差。ELS法：ELS是极大似然法的一种近似形式，当数据长度较小时，辨识精度可能优于极大似然法，但数据长度较大时，精度低于极大似然法。第三题以上题的结果为例，进行：(10分)1 .分析比较各种方法估计的精度；2 .分析其计算量；3 .分析噪声方差的影响；4 .比较白噪声和有色噪声对

43、辨识的影响。答：1 .分析比较各种方法估计的精度；最小二乘估计虽然不能满足量测方程中的每一个方程，使每个方程都有偏差， 但它使所有方程偏差的平方和达到最小，兼顾了所有方程的近似程度，使整体误差达到最小。辅助变量法、最小二乘法、极大似然法对输入输出模型有相同的精度增广最小二乘法所得结果中b2误差较大；辅助变量法辨识精度高于最小二乘法 夏氏法辨识精度较高；克服有偏估计。2 .分析其计算量；最小二乘估计计算量大、存储大；增广最小二乘递推算法，扩充了最小二乘法的参数向量和数据向量的维数，把噪声模型的辨识同时考虑进去，所以计算量较大； 递推广义LS中的计算分为二部分： 第一部分：按递推 LS法，随N增大

44、不断计算参数所以计算量较大， 第二部分：在递推过程中，由于是变化的故而需要不断计算；辅助变量法：(1)计算与基本最小二乘估计同样简单；(2) 辨识精度高于最小 LS估计法；(3) 是一种无偏估计方法；(4) 参数估计时需构造辅助变量矩阵。3 氏法：(1)无偏的估计方法；(2) 计算量较广义LS的要小许多；(3) 不需进行I/O数据的反复过滤，计算效率高；(4) 其递推算法可推广至 MIMO系统，而广义LS法则不行；(5) 是一种循环迭代方法，收敛速度较LS要慢；(6) 估计结果较好，可分为夏式修正法和夏式改良法两种极大似然；计算量较小。3 .分析噪声方差的影响；噪声越大，对最小二乘法辨识效果影响较大，对于增广递推最小二乘法， 抗噪声能力较强。通过编程计算，获得在噪声方差比较小的情况下，各种方法所获得的估值比较理想， 但随着噪声方差的增大，估值的偏差随之增大，横向比较看来夏式法与广义最小二乘法能够 更好地还原参数值，当观测值足够多时，各种方法都能很好地反映参数真实值。4 .比较白噪声和有色噪声对辨识的影响。由于所用的输出观测值有有色噪声成分和白噪声成分，