系统辨识-最小二乘法

5.3 Matlab源程序：

%最小二乘估计

clear

u=[ 1.147 0.201 -0.787 -1.589 -1.052 0.866 1.152 1.573 0.626

0.433 -0.985 0.810 -0.044 0.947 -1.474 -0.719 -0.086 -1.099 1.450

1.151 0.485 1.633 0.043 1.326 1.706 -0.340 0.890 0.144 1.177

-0.390];

n=normrnd(0, sqrt(0.1), 1, 31); %方差为0.1

z=zeros(1,30);

for k=3:31

z(k)=-1.642*z(k-1)-0.715*z(k-2)+0.39*u(k-1)+0.35*u(k-2)+n(k)+1.642*n(k-1)+0.715 *n(k-2);

end

h0=[-z(2) -z(1) u(2) u(1)]';

HLT=[h0,zeros(4,28)];

for k=3:30

h1=[-z(k) -z(k-1) u(k) u(k-1)]';

HLT(:,k-1)=h1;

end

HL=HLT';

y=[z(3);z(4);z(5);z(6);z(7);z(8);z(9);z(10);z(11);z(12);z(13);z(14);z(15);z(16) ;z(17);z(18);z(19);z(20);z(21);z(22);z(23);z(24);z(25);z(26);z(27);z(28);z(29); z(30);z(31)];%求出FAI

c1=HL'*HL;

c2=inv(c1);

c3=HL'*y;

c=c2*c3;

%display('方差=0.1时，最小二乘法估计辨识参数θ如下：');

a1=c(1)

a2=c(2)

b1=c(3)

b2=c(4)

%递推最小二乘法估计

clear

u=[ 1.147 0.201 -0.787 -1.589 -1.052 0.866 1.152 1.573 0.626

0.433 -0.985 0.810 -0.044 0.947 -1.474 -0.719 -0.086 -1.099 1.450

1.151 0.485 1.633 0.043 1.326 1.706 -0.340 0.890 0.144 1.177

-0.390];

z(2)=0;

z(1)=0;

n=normrnd(0, sqrt(0.1), 1, 31); %方差为0.1

for k=3:31

z(k)=-1.642*z(k-1)-0.715*z(k-2)+0.39*u(k-1)+0.35*u(k-2)+n(k)+1.642*n(k-1)+0.715 *n(k-2);

end

c0=[0.001 0.001 0.001 0.001]'; %直接给出被辨识参数的初始值,即一个充分小的实向量

p0=10^6*eye(4,4); %直接给出初始状态P0，即一个充分大的实数单位矩阵

E=0.000000005; %取相对误差E=0.000000005

c=[c0,zeros(4,30)]; %被辨识参数矩阵的初始值及大小

e=zeros(4,30); %相对误差的初始值及大小

for k=3:30; %开始求K

h1=[-z(k-1),-z(k-2),u(k-1),u(k-2)]';

x=h1'*p0*h1+1; x1=inv(x); %开始求K(k)

k1=p0*h1*x1;%求出K的值

d1=z(k)-h1'*c0; c1=c0+k1*d1; %求被辨识参数c

e1=c1-c0; %求参数当前值与上一次的值的差值

e2=e1./c0; %求参数的相对变化

e(:,k)=e2; %把当前相对变化的列向量加入误差矩阵的最后一列

c0=c1; %新获得的参数作为下一次递推的旧参数

c(:,k)=c1; %把辨识参数c 列向量加入辨识参数矩阵的最后一列

p1=p0-k1*k1'*[h1'*p0*h1+1]; %求出 p(k)的值

p0=p1; %给下次用

if e2<=E

break; %如果参数收敛情况满足要求，终止计算

end

end

%display('方差为0. 1递推最小二乘法辨识后的结果是：');

a1=c(1,:);

a2=c(2,:);

b1=c(3,:);

b2=c(4,:);

%display('a1,a2,b1,b2经过递推最小二乘法辨识的结果是：');

for i=3:31;

if(c(1,i)==0)

q1=c(1,i-1);

break;

end

end

for i=3:31;

if(c(2,i)==0)

q2=c(2,i-1);

break;

end

end

for i=3:31;

if(c(3,i)==0)

q3=c(3,i-1);

break;

end

end

for i=3:31;

if(c(4,i)==0)

q4=c(4,i-1);

break;

end

end

a1=q1

a2=q2

b1=q3

b2=q4

%辅助变量递推最小二乘法估计

clear

na=2;

nb=2;

siitt=[1.642 0.715 0.39 0.35]';

siit0=0.001*eye(na+nb,1);

p=10^6*eye(na+nb);

siit(:,1)=siit0;

y(2)=0;y(1)=0;

x(1)=0;x(2)=0;

j=0;

u=[ 1.147 0.201 -0.787 -1.589 -1.052 0.866 1.152 1.573 0.626

0.433 -0.985 0.810 -0.044 0.947 -1.474 -0.719 -0.086 -1.099 1.450

1.151 0.485 1.633 0.043 1.326 1.706 -0.340 0.890 0.144 1.177

-0.390];

n=normrnd(0, sqrt(0.1), 1, 31); %方差为0.1

for k=3:31;

h=[-y(k-1),-y(k-2),u(k-1),u(k-2)]';

y(k)=h'*siitt+n(k)+1.642*n(k-1)+0.715*n(k-2);

hx=[-x(k-1),-x(k-2),u(k-1),u(k-2)]';

kk=p*hx/(h'*p*hx+1);

p=[eye(na+nb)-kk*h']*p;

siit(:,k-1)=siit0+kk*[y(k)-h'*siit0];

x(k)=hx'*siit(:,k-1);

j=j+(y(k)-h'*[1.642 0.715 0.39 0.35]')^2;

e=max(abs((siit(:,k-1)-siit0)./siit0));

ess(:,k-2)=siit(:,k-1)-siitt;