最小二乘法在系统辨识中的应用包含相关的三种算法

2008级硕士研究生《系统建模理论》试卷

已知一个三阶线性离散系统的输入、输出数据，共有40个采样值，试分别用：最小二乘法（LS）、递推最小二乘法（RLS）、广义最小二乘法（GLS）进行参数估计，给出相应的数学模型，并阐述相应的辨识原理。

k 1 2 3 4 5 6 7 8 9 10

u(k) 0.8251 0.0988 0.4628 -0.9168 2.2325 0.0777 2.3654 0.3476 1.1473 -1.9035

y(k) 1.5333 -1.0680 1.0666 -0.5284 -0.5835 3.1471 -3.7185 6.2149 -6.3026 7.2705

k 11 12 13 14 15 16 17 18 19 20

u(k) -0.9229 1.6400 -0.8410 0.7599 -0.4739 -0.1784 -1.7760 -1.6722 1.2959 -0.0591

y(k) -9.0552 8.1735 -5.9004 3.9870 -2.2486 0.9525 -0.5325 -1.5227 0.4200 1.0786

k 21 22 23 24 25 26 27 28 29 30

u(k) -1.0576 -1.0071 1.1342 -0.0740 0.6759 0.5221 0.9954 0.5271 -1.7656 0.4936

y(k) -1.5579 0.6640 -1.4222 2.6444 -2.9572 3.6340 -3.1281 3.8334 -3.2542 1.1568

k 31 32 33 34 35 36 37 38 39 40

u(k) 1.4810 0.9591 -3.1293 -0.3604 -0.4251 0.4185 -0.6728 -0.0027 2.1145 1.1157

y(k) 0.0615 0.9120 -0.0692 -3.2731 3.7486 -4.3194 4.7230 -5.2781 5.1507 -2.7235

一、最小二乘法（LS）

1、数学模型

设时不变SISO动态过程的数学模型为

（1.1.1）

其中，为过程的输入量，为过程的输出量，是噪声，多项式和为：

在本题中，==3.即

将此模型写成最小二乘格式 （1.1.2）

其中，是过程的输出量；是可观测的数据向量；是均值为零的随机噪声。

式中

对于，方程式（1.1.2）构成一个线性方程组，可以把它写成

（1.1.3）

利用数据序列和，极小化准则函数 （1.1.4）

使的估计值记作，称为参数的最小二乘估计值。通过极小化（1.1.4）式来计算的方法称作最小二乘法，

未知模型参数最可能的值是在实际观测值与计算值之累次误差的平方和达到最小处所得到的，这种模型输出能最好地接近实际过程的输出。

2、辨识原理

考虑模型（1.1.2）式的辨识问题，其中和都是可观测的数据，是待估计参数，准则函数取（1.1.4）根据（1.1.3）的定义，准则函数可写成二次型的形式

（1.2.1）

显然上式中的代表模型的输出，或者说是过程的输出预报值。因此可以看作来衡量模型输出与实际过程输出的接近情况。极小化，求得参数的估计值

（1.2.2）

将使模型的输出最好的预报过程的输出。

3、辨识结果

二、递推最小二乘法（RLS）

1、数学模型

在第一部分中建立了最小二乘法，并用一次完成算法进行了计算。但是由于具体使用时占用内存量大，而且不能用于在线辨识。所以引入了最小二乘参数估计的递推算法。

递推算法的基本思想可以概括成：

新的估计值=老的估计值+修正项 （2.1.1）

在此算法中，时不变SISO动态过程的数学模型仍与最小二乘法的一样。模型的最小二乘格式也相同，只是计算方法不同，具体计算方法，在递推最小二乘法的辨识原理中描述。

2、辨识原理

首先将1.2.2式一次完成算法写成

（2.2.1）

定义

（2.2.2）

其中：，

由（2.2.2）式可得：

（2.2.3）

令：

则：

于是有 （2.2.4）

令

利用（2.2.3）和（2.2.4）式，可得

（2.2.5）

引进增益矩阵，定义为

（2.2.6）

则（2.2.5）式写成

（2.2.7）

进一步把（2.2.3）式写成

（2.2.8）

为了避免矩阵求逆运算，利用矩阵反演公式可将（2.2.8）式演变成

（2.2.9）

将（2.2.9）式代入（2.2.6）式，整理后有

（2.2.10）

综合（2.2.7）、（2.2.9）、（2.2.10）式便得到最小二乘参数估计递推算法。

利用上述公式即可求得参数的估计值

3、辨识结果

三、广义最小二乘法（GLS）

1、数学模型

设时不变SISO动态过程的数学模型为

（3.1.1）

其中，和 分别表示过程的输入和输出，是均值为零的不相关随机噪声，多项式、和为：

在本题中，在本题中，===3.即

（3.1.2）

广义最小二乘发的基本思想是基于对数据先进行一次滤波预处理，然后利用普通最小二乘法对滤波后的数据进行辨识。

2、辨识原理

由（3.1.1）式可得 （3.2.1）

令， （3.2.2）

及（3.2.3）

可将模型（3.1.1）式化成最小二乘格式 （3.2.4）

由于上式是白噪声，所以利用最小二乘法即可获得参数的无偏估计。但是，数据向量中的变量均需要按照（3.2.2）式计算，然而噪声模型并不知道。为此需要用迭代的方法来估计。

令 （3.2.5）

置 （3.2.6）

就把噪声模型（3.2.5）也化成最小二乘格式

由于上式的噪声已是白噪声，所以再次利用最小二乘法可获得噪声模型参数的无偏估计。但是，数据向量依然包含着不可测得噪声量，它可用相应的估计值代替，置 （3.2.7）

其中，当时，按照式子计算，

式子中

综上分析，广义最小二乘递推算法可归纳成：

利用上述公式即可求得参数的估计值并能求出噪声模型的估计。

3、估计结果

噪声模型的估计值为1.0e-003 *[ -0.0000 0.4115 0.0001 ]即

附录：MATLAB程序

1、一次完成的最小二乘法

%%一次完成的最小二乘法

clear

s=(1:40);

z=[1.5333 -1.0680 1.0666 -0.5284 -0.5835 3.1471 -3.7185 6.2149 -6.3026 7.2705...

-9.0552 8.1735 -5.9004 3.9870 -2.2486 0.9525 -0.5325 -1.5227 0.4200 1.0786...

-1.5579 0.6640 -1.4222 2.6444 -2.9572 3.6340 -3.1281 3.8334 -3.2542 1.1568...

0.0615 0.9120 -0.0692 -3.2731 3.7486 -4.3194 4.7230 -5.2781 5.1507 -2.7235];

u=[0.8251 0.0988 0.4628 -0.9168 2.2325 0.0777 2.3654 0.3476 1.1473 -1.9035...

-0.9229 1.6400 -0.8410 0.7599 -0.4739 -0.1784 -1.7760 -1.6722 1.2959 -0.0591...

-1.0576 -1.0071 1.1342 -0.0740 0.6759 0.5221 0.9954 0.5271 -1.7656 0.4936...

1.4810 0.9591 -3.1293 -0.3604 -0.4251 0.4185 -0.6728 -0.0027 2.1145 1.1157];

h=zeros(40,6);

h(1,:)=[-z(1) 0 0 u(1) 0 0];

h(2,:)=[-z(2) -z(1) 0 u(2) u(1) 0];

for i=3:1:40

h(i,:)=[-z(i) -z(i-1) -z(i-2) u(i) u(i-1) u(i-2)];

end

o=inv(h'*h)*h'*z'

%误差分析

J=0;

for i=1:1:40

e(i)=[z(i)-h(i,:)*o];

J=J+e(i)^2;

end

J

plot(s,e);

2、递推的最小二乘法

%%递推的最小二乘法

clear

s=(1:40);

z=[1.5333 -1.0680 1.0666 -0.5284 -0.5835 3.1471 -3.7185 6.2149 -6.3026 7.2705...

-9.0552 8.1735 -5.9004 3.9870 -2.2486 0.9525 -0.5325 -1.5227 0.4200 1.0786...

-1.5579 0.6640 -1.4222 2.6444 -2.9572 3.6340 -3.1281 3.8334 -3.2542 1.1568...

0.0615 0.9120 -0.0692 -3.2731 3.7486 -4.3194 4.7230 -5.2781 5.1507 -2.7235];

u=[0.8251 0.0988 0.4628 -0.9168 2.2325 0.0777 2.3654 0.3476 1.1473 -1.9035...

-0.9229 1.6400 -0.8410 0.7599 -0.4739 -0.1784 -1.7760 -1.6722 1.2959 -0.0591...

-1.0576 -1.0071 1.1342 -0.0740 0.6759 0.5221 0.9954 0.5271 -1.7656 0.4936...

1.4810 0.9591 -3.1293 -0.3604 -0.4251 0.4185 -0.6728 -0.0027 2.1145 1.1157];

h=zeros(40,6);

h(1,:)=[-z(1) 0 0 u(1) 0 0];

h(2,:)=[-z(2) -z(1) 0 u(2) u(1) 0];

for i=3:1:40

h(i,:)=[-z(i) -z(i-1) -z(i-2) u(i) u(i-1) u(i-2)];

end

p0=10^6;

o0=0.001;

k(1,:)=[p0*h(1,:)'*inv(h(1,:)*p0*h(1,:)'+1)]';

K=zeros(6,40);

K(:,1)=k(1,:)';

p=zeros(6,6,40);

p(:,:,1)=[eye(6)-K(:,1)*h(1,:)]*p0;

o=zeros(6,40)';

o(1,:)=o0+K(:,1)*[z(1)-h(1)*o0];

for i=2:1:40

k(i,:)=[p(:,:,i-1)*h(i,:)'*inv(h(i,:)*p(:,:,i-1)*h(i,:)'+1)]';

p(:,:,i)=[eye(6)-k(i,:)'*h(i,:)]*p(:,:,i-1);

o(i,:)=[o(i-1,:)'+k(i,:)'*[z(i)-h(i,:)*o(i-1,:)']]';

end

%辨识结果

o(40,:)

%误差分析

J=0;

for i=2:1:40

e(i)=[z(i)-h(i,:)*o(i-1,:)'];

J=J+e(i)^2;

end

plot(s,e);

J

3、广义最小二乘法

%%广义预测的最小二乘法

clear

s=(1:40);

z=[1.5333 -1.0680 1.0666 -0.5284 -0.5835 3.1471 -3.7185 6.2149 -6.3026 7.2705...

-9.0552 8.1735 -5.9004 3.9870 -2.2486 0.9525 -0.5325 -1.5227 0.4200 1.0786...

-1.5579 0.6640 -1.4222 2.6444 -2.9572 3.6340 -3.1281 3.8334 -3.2542 1.1568...

0.0615 0.9120 -0.0692 -3.2731 3.7486 -4.3194 4.7230 -5.2781 5.1507 -2.7235];

u=[0.8251 0.0988 0.4628 -0.9168 2.2325 0.0777 2.3654 0.3476 1.1473 -1.9035...

-0.9229 1.6400 -0.8410 0.7599 -0.4739 -0.1784 -1.7760 -1.6722 1.2959 -0.0591...

-1.0576 -1.0071 1.1342 -0.0740 0.6759 0.5221 0.9954 0.5271 -1.7656 0.4936...

1.4810 0.9591 -3.1293 -0.3604 -0.4251 0.4185 -0.6728 -0.0027 2.1145 1.1157];

zf=z;

uf=u;

hf=zeros(40,6);

hf(1,:)=[-zf(1) 0 0 uf(1) 0 0];

hf(2,:)=[-zf(2) -zf(1) 0 uf(2) uf(1) 0];

o0=0.001;

pf0=10^6;

oe0=0;

pe0=1;

%k=1时

kf(1,:)=[pf0*hf(1,:)'*inv(hf(1,:)*pf0*hf(1,:)'+1)]';

pf=zeros(6,6,40);

pf(:,:,1)=[eye(6)-kf(1,:)'*hf(1,:)]*pf0;

of=zeros(6,40)';

of(1,:)=o0+kf(1,:)'*[z(1)-hf(1)*o0];

ke=zeros(40,3);

e(1)=z(1)-hf(1,:)*of(1,:)';

he(1,:)=[0 0 0];

ke(1,:)=[pe0*he(1,:)'*inv(he(1,:)*pe0*he(1,:)'+1)]';

pe=zeros(3,3,40);

pe(:,:,1)=[eye(3)-ke(1,:)'*he(1,:)]*pe0;

oe=zeros(3,40)';

oe(1,:)=oe0+ke(1,:)'*[e(1)-he(1)*oe0];

%k=2时

kf(2,:)=[pf(:,:,2-1)*hf(2,:)'*inv(hf(2,:)*pf(:,:,2-1)*hf(2,:)'+1)]';

pf(:,:,2)=[eye(6)-kf(2,:)'*hf(2,:)]*pf(:,:,2-1);

of(2,:)=[of(2-1,:)'+kf(2,:)'*[z(2)-hf(2,:)*of(2-1,:)']]';

e(2)=z(2)-hf(2,:)*of(2,:)';

he(2,:)=[-e(1) 0 0];

ke(2,:)=[pe(:,:,2-1)*he(2,:)'*inv(he(2,:)*pe(:,:,2-1)*he(2,:)'+1)]';

pe(:,:,2)=[eye(3)-ke(2,:)'*he(2,:)]*pe(:,:,2-1);

oe(2,:)=[oe(2-1,:)'+ke(2,:)'*[e(2)-he(2,:)*oe(2-1,:)']]';

%k=3时

kf(3,:)=[pf(:,:,3-1)*hf(3,:)'*inv(hf(3,:)*pf(:,:,3-1)*hf(3,:)'+1)]';

pf(:,:,3)=[eye(6)-kf(3,:)'*hf(3,:)]*pf(:,:,3-1);

of(3,:)=[of(3-1,:)'+kf(3,:)'*[z(3)-hf(3,:)*of(3-1,:)']]';

e(3)=z(3)-hf(3,:)*of(3,:)';

he(3,:)=[-e(2) -e(1) 0];

ke(3,:)=[pe(:,:,3-1)*he(3,:)'*inv(he(3,:)*pe(:,:,3-1)*he(3,:)'+1)]';

pe(:,:,3)=[eye(3)-ke(3,:)'*he(3,:)]*pe(:,:,3-1);

oe(3,:)=[oe(3-1,:)'+ke(3,:)'*[e(3)-he(3,:)*oe(3-1,:)']]';

%k>4时

for i=4:1:40

zf(i)=z(i)+z(i-1)*oe(i-1,1)+z(i-2)*oe(i-1,2)+z(i-3)*oe(i-1,3);

uf(i)=u(i)+u(i-1)*oe(i-1,1)+u(i-2)*oe(i-1,2)+u(i-3)*oe(i-1,3);

hf(i,:)=[-zf(i) -zf(i-1) -zf(i-2) uf(i) uf(i-1) uf(i-2)];

kf(i,:)=[pf(:,:,i-1)*hf(i,:)'*inv(hf(i,:)*pf(:,:,i-1)*hf(i,:)'+1)]';

pf(:,:,i)=[eye(6)-kf(i,:)'*hf(i,:)]*pf(:,:,i-1);

of(i,:)=[of(i-1,:)'+kf(i,:)'*[z(i)-hf(i,:)*of(i-1,:)']]';

e(i)=z(i)-hf(i,:)*of(i,:)';

he(i,:)=[-e(i-1) -e(i-2) -e(i-3)];

ke(i,:)=[pe(:,:,i-1)*he(i,:)'*inv(he(i,:)*pe(:,:,i-1)*he(i,:)'+1)]';

pe(:,:,i)=[eye(3)-ke(i,:)'*he(i,:)]*pe(:,:,i-1);

oe(i,:)=[oe(i-1,:)'+ke(i,:)'*[e(i)-he(i,:)*oe(i-1,:)']]';

end

%辨识结果

of(40,:)

oe(40,:)

%误差分析

plot(s,e);

本文来源：https://www.2haoxitong.net/k/doc/90ce5cfb54270722192e453610661ed9ac5155e6.html