系统辨识大作业1201张青

《系统辨识》大作业

学号：12051124

班级：自动化1班

姓名：张青

信息与控制工程学院自动化系

2015-07-11

第一题

模仿index2，搭建对象，由相关分析法，获得脉冲响应序列?()g k ，由?

()g k ，参照讲义， 获得系统的脉冲传递函数()G z 和传递函数()G s ;应用最小二乘辨识,获得脉冲响应序列?

()g k ;同图显示两种方法的辨识效果图；应用相关最小二乘法，拟合对象的差分方程模型；构建时变对象，用最小二乘法和带遗忘因子的最小二乘法，(可以用辨识工具箱) 辨识模型的参数，比较两种方法的辨识效果差异;

答：根据index2搭建结构框图：

相关分析法：利用结构框图得到UY 和tout

其中的U 就是题目中要求得出的M 序列，根据结构框图可知序列的周期是

1512124=-=-=n p N 。

在command window 中输入下列指令，既可以得到脉冲相应序列()g k

：

aa=5;NNPP=15;ts=2; RR=ones(15)+eye(15); for i=15:-1:1

UU(16-i,:)=UY(16+i:30+i,1)'; end

YY=[UY(31:45,2)];

GG=RR*UU*YY/[aa*aa*(NNPP+1)*ts]; plot(0:2:29,GG) hold on

stem(0:2:29,GG,'filled') Grid;title('脉冲序列g(τ)')

最小二乘法建模的响应序列

由于是二阶水箱系统，可以假设系统的传递函数为2

21101)(s

a s a s

b 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); end

Estimate=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:201

h=[-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;

end

i=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

(

2

1

2

1

k

k

u

b

k

u

b

k

y

a

k

y

aξ

+

-

+

-

+

-

-

-

-

辨识参数向量为

θ＝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<1

disp('从下列文件中选择所需加载文件：uy1.txt uy2.txt uy3.txt') filename = input('输入所需打开文件名:\','s');

[Wbs,Message] = fopen(filename,'r');%打开文件；

end

Data = fscanf(Wbs,'%g %g',[2 inf]);%生成数据矩阵

Data = Data';

%最小二乘法，取b0=0

n = 2;%根据题目，本系统为2阶系统

L = length(Data);%辨识所需数据的总长度

N = L-n;%构造测量矩阵的数据长度

FIA = zeros(N,2*n);%构造测量矩阵

for i = 1:N %测量矩阵赋值

for l = 1:n*2

if(l>n)

FIA(i,l) = Data(n+2+i-l,1);

else

FIA(i,l) = -Data(n+i-l,2);

end

end

end

Y = 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.txt

2.递推最小二乘法：

%首先将TXT中的数据装载在Matlab中：

clear;clc;

Wbs = 0;

while Wbs<1

disp('从下列文件中选择所需加载文件：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;%根据题目已知该系统阶数为2

L = length(Data);

%为辨识参数向量赋初值，这里均取为0.001

for a = 1:1:n*2

thita0(a,1) = 0.001;

end

%直接给出矩阵Pn的初始状态P0

P0 = 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);

else

FIA1(m,1) = Data(i-m+2,1);

end

end

X = 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;%所得的参数向量作为下次迭代初值

end

disp('递推最小二乘法所得结果如下：')

Dta1=thita1(1),Dta2=thita1(2),Dtb1=thita1(3),Dtb2=thita1(4) %显示被辨识参数

实验结果：

Uy1.txt uy2.txt uy3.txt

3.辅助变量法：

%首先将TXT中的数据装载在Matlab中：

clear;clc;

Wbs = 0;

while Wbs<1

disp('从下列文件中选择所需加载文件：uy1.txt uy2.txt uy3.txt')

filename = input('输入所需打开文件名:\','s');

[Wbs,Message] = fopen(filename,'r');%打开文件；

end

Data = fscanf(Wbs,'%g %g',[2 inf]);%生成数据矩阵

Data = Data';

%使用递推辅助变量法对2阶系统进行辨识，取b0=0

n = 2;L = length(Data);N = L-n;

%分别取出输入、输出数据u,y

U = 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:400

Yf = 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;%求取被估计参数的辅助变量估值end

disp('使用辅助变量法，得到如下辨识结果：')

Fza1=thita(1),Fza2=thita(2),Fzb1=thita(3),Fzb2=thita(4)

实验结果：

Uy1.txt uy2.txt uy3.txt

4.递推辅助变量法：

%首先将TXT中的数据装载在Matlab中：

clear;clc;

Wbs = 0;

while Wbs<1

disp('从下列文件中选择所需加载文件：uy1.txt uy2.txt uy3.txt') filename = input('输入所需打开文件名:\','s');

[Wbs,Message] = fopen(filename,'r');%打开文件；

end

Data = fscanf(Wbs,'%g %g',[2 inf]);%生成数据矩阵

Data = Data';

%使用递推辅助变量法对2阶系统进行辨识，取b0=0

n = 2;L = length(Data);

%分别取出输入、输出数据u,y

U = 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:N

Y11 = [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);

end

disp('递推辅助变量法辨识结果：');

Dfa1=thita(1),Dfa2=thita(2),Dfb1=thita(3),Dfb2=thita(4)

实验结果：

Uy1.txt uy2.txt uy3.txt

5.广义最小二乘法：

%首先将TXT中的数据装载在Matlab中：

clear;clc;

Wbs = 0;

while Wbs<1

disp('从下列文件中选择所需加载文件：uy1.txt uy2.txt uy3.txt')

filename = input('输入所需打开文件名:\','s');

[Wbs,Message] = fopen(filename,'r');%打开文件；

end

Data = fscanf(Wbs,'%g %g',[2 inf]);%生成数据矩阵

Data = Data';

%使用广义最小二乘法，辨识2阶系统的参数

n = 2;

L = length(Data);

N = L-n;

for m=1:1:300

%分别取出输入、输出数据u,y

U = 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);%输出矩阵y

thita = 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);

end

E1 = 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+n

Data(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);

end

end

disp('由广义最小二乘法得到的辨识结果如下：')

Gya1=thita(1),Gya2=thita(2),Gyb1=thita(3),Gyb2=thita(4)

实验结果：

Uy1.txt uy2.txt uy3.txt

6.递推广义最小二乘法：

%首先将TXT中的数据装载在Matlab中：

clear;clc;

Wbs = 0;

while Wbs<1

disp('从下列文件中选择所需加载文件：uy1.txt uy2.txt uy3.txt') filename = input('输入所需打开文件名:\','s');

[Wbs,Message] = fopen(filename,'r');%打开文件；

end

Data = fscanf(Wbs,'%g %g',[2 inf]);%生成数据矩阵

Data = Data';

u1=Data(:,1);

y1=Data(:,2);

n=2;

N=500-n;

n0=50;

%用前50组数据利用lsm估计出初值

%构造矩阵phi

phi(:,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); end

disp('递推广义最小二乘法辨识结果：');

a1=seta0(1),a2=seta0(2),b1=seta0(3),b2=seta0(4) 实验结果：

Uy1.txt uy2.txt uy3.txt

7.夏氏偏差修正法：

%首先将TXT中的数据装载在Matlab中：

clear;clc;

Wbs = 0;

while Wbs<1

disp('从下列文件中选择所需加载文件：uy1.txt uy2.txt uy3.txt') filename = input('输入所需打开文件名:\','s');

[Wbs,Message] = fopen(filename,'r');%打开文件；

end

Data = 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;%设偏差项的偏差初值为0

Fa = Sgl2*glOL';

M = eye(N)-glOL*Sgl2*glOL';

F = eye(N);

if(F>=10^-6*eye(N))

E1 = Zgl1-glOL*Xsthita;%计算残差E

omiga(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;

end

disp('夏氏修正法')

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<1

disp('从下列文件中选择所需加载文件：uy1.txt uy2.txt uy3.txt') filename = input('输入所需打开文件名:\','s');

[Wbs,Message] = fopen(filename,'r');%打开文件；

end

Data = fscanf(Wbs,'%g %g',[2 inf]);%生成数据矩阵

Data = Data';

%最小二乘法，取b0=0

n = 2;%根据题目，本系统为2阶系统

L = length(Data);%辨识所需数据的总长度

N = L-n;%构造测量矩阵的数据长度

FIA = zeros(N,2*n);%构造测量矩阵

for i = 1:N %测量矩阵赋值

for l = 1:n*2

if(l>n)

FIA(i,l) = Data(n+2+i-l,1);

else

FIA(i,l) = -Data(n+i-l,2);

end

end

end

Y = Data(n+1:n+N,2);%输出数据矩阵

thita = inv(FIA'*FIA)*FIA'*Y;%计算参数矩阵

thitaLS = thita;%最小二乘估值

m = inv(FIA'*FIA)*FIA';

for i = 1:1:1000%构造Omega

Err = Y-FIA*thitaLS;%计算残差E

omiga(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;

end

disp('使用夏氏改良法辨识结果为：')

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<1

disp('从下列文件中选择所需加载文件：uy1.txt uy2.txt uy3.txt')

filename = input('输入所需打开文件名:\','s');

[Wbs,Message] = fopen(filename,'r');%打开文件；

end

Data = 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:N

E1 = 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;

end

disp('取450组数据进行递推夏氏辨识的结果为:');

Dxa1=beta(1),Dxa2=beta(2),Dxb1=beta(3),Dxb2=beta(4)

实验结果：

Uy1.txt uy2.txt uy3.txt

10.增广矩阵法：

%首先将TXT中的数据装载在Matlab中：

clear;clc;

Wbs = 0;

while Wbs<1

disp('从下列文件中选择所需加载文件：uy1.txt uy2.txt uy3.txt') filename = input('输入所需打开文件名:\','s');

[Wbs,Message] = fopen(filename,'r');%打开文件；

end

Data = 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为0

thita0(i,1) = 0;

end

P0 = 10^10*eye(3*n,3*n);%然后设定初始状态P0为充分大的对角阵

ERR0 = 0.00000000001;%收敛情况指标，相对误差ERR0

thita = [thita0,zeros(3*n,L)];%设定参数向量的初值及大小

zs = zeros(L,1);%构造初始噪声矩阵

hn = [0,0,0,0,0,0]';

for j = n+1:1:L

zs(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);%计算修正矩阵K

thita0 = 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;%根据设定的误差限，判断何时跳出

end

end

Zgma1=thita0(1),Zgma2=thita0(2),Zgmb1=thita0(3),Zgmb2=thita0(4), Zgmc1=thita0(5),Zgmc2=thita0(6)

Uy1.txt uy2.txt uy3.txt

11.极大似然法：

%首先将TXT中的数据装载在Matlab中：

clear;clc;

Wbs = 0;

while Wbs<1

disp('从下列文件中选择所需加载文件：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 = [-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:500

Ja1 = mthita(1);Ja2 = mthita(2);Jb1 = mthita(3);Jb2 = mthita(4);Jc1 = mthita(5);Jc2 = mthita(6);

for i = n+1:L

y2(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);

end

for j = n+1:1:L

Ea1(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 break

end

end

disp('极大似然法（牛顿-拉卜森方法），得到的结果为：')

Ja1 = mthita(1),Ja2 = mthita(2),Jb1 = mthita(3),Jb2 = mthita(4),Jc1 = mthita(5),Jc2 = mthita(6)

实验结果：

Uy1.txt uy2.txt uy3.txt

12.分析噪声)(k 特性:

RLS 法：当噪声比较小时，辨识精度较高；当噪声比较大时，收敛点可能

不唯一，参数估计值往往是有偏的；但是运用此方法时，数据要充分多，否则辨识精度明显下降；另外噪声模型阶次不宜过高，也会影响精度，所以误差的原因主要是数据量的大小以及噪声模型阶次。

IV 法：辨识效果较好，能适应较大的范围的噪声特性。对初值 P(0) 敏感，

可用LS 法起步，否则没有可靠的收敛性；对数据 u(k) 在 z(k) 的直流分量敏感，对z(k)的直流分量不敏感。所以误差的主要原因是初值P(0)以及数据计算过程产生的误差。 ELS 法：ELS 是极大似然法的一种近似形式，当数据长度较小时，辨识精 度可能优于极大似然法，但数据长度较大时，精度低于极大似然法。

第三题

以上题的结果为例，进行：（10分）

1. 分析比较各种方法估计的精度；

2. 分析其计算量；

3. 分析噪声方差的影响；

4. 比较白噪声和有色噪声对辨识的影响。

答：

1. 分析比较各种方法估计的精度；

最小二乘估计虽然不能满足量测方程中的每一个方程，使每个方程都有偏差，但它使所有方程偏差的平方和达到最小，兼顾了所有方程的近似程度，使整体误差达到最小。辅助变量法、最小二乘法、极大似然法对输入输出模型有相同的精度.

增广最小二乘法所得结果中b2误差较大；

辅助变量法辨识精度高于最小二乘法

夏氏法辨识精度较高；克服有偏估计。

2. 分析其计算量；

最小二乘估计计算量大、存储大；

增广最小二乘递推算法，扩充了最小二乘法的参数向量和数据向量的维数,把噪声模型的辨识同时考虑进去，所以计算量较大；

递推广义LS中的计算分为二部分：

第一部分：按递推LS法，随N增大不断计算参数所以计算量较大，

第二部分：在递推过程中，由于是变化的故而需要不断计算；

辅助变量法： (1) 计算与基本最小二乘估计同样简单；

(2) 辨识精度高于最小LS估计法；

(3) 是一种无偏估计方法；

(4) 参数估计时需构造辅助变量矩阵。

夏氏法：(1)无偏的估计方法；

(2)计算量较广义LS的要小许多；

(3)不需进行I/O数据的反复过滤，计算效率高；

(4)其递推算法可推广至MIMO系统，而广义LS 法则不行；

(5)是一种循环迭代方法，收敛速度较LS要慢；

(6)估计结果较好，可分为夏式修正法和夏式改良法两种

极大似然;计算量较小。

3. 分析噪声方差的影响；

噪声越大，对最小二乘法辨识效果影响较大，对于增广递推最小二乘法，抗噪声能力较强。

通过编程计算，获得在噪声方差比较小的情况下，各种方法所获得的估值比较理想，但随着噪声方差的增大，估值的偏差随之增大，横向比较看来夏式法与广义最小二乘法能够

更好地还原参数值，当观测值足够多时，各种方法都能很好地反映参数真实值。

4. 比较白噪声和有色噪声对辨识的影响。

由于所用的输出观测值有有色噪声成分和白噪声成分，所以辨识结果有误差，且用最小二乘法计算得到的参数估计误差较大。但是有色噪声误差明显大于白噪声。