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

1、2008级硕士研究生系统建模理论试卷已知一个三阶线性离散系统的输入、输出数据，共有40个采样值，试分别用：最小二乘法（LS）、递推最小二乘法（RLS）、广义最小二乘法（GLS）进行参数估计，给出相应的数学模型，并阐述相应的辨识原理。k 1 2 3 4 5 6 7 8 9 10u(k) 0.8251 0.0988 0.4628 -0.9168 2.2325 0.0777 2.3654 0.3476 1.1473 -1.9035y(k) 1.5333 -1.0680 1.0666 -0.5284 -0.5835 3.1471 -3.7185 6.2149 -6.3026 7.2705k 11 12

2、 13 14 15 16 17 18 19 20u(k) -0.9229 1.6400 -0.8410 0.7599 -0.4739 -0.1784 -1.7760 -1.6722 1.2959 -0.0591y(k) -9.0552 8.1735 -5.9004 3.9870 -2.2486 0.9525 -0.5325 -1.5227 0.4200 1.0786k 21 22 23 24 25 26 27 28 29 30u(k) -1.0576 -1.0071 1.1342 -0.0740 0.6759 0.5221 0.9954 0.5271 -1.7656 0.4936y(k) -1

3、.5579 0.6640 -1.4222 2.6444 -2.9572 3.6340 -3.1281 3.8334 -3.2542 1.1568k 31 32 33 34 35 36 37 38 39 40u(k) 1.4810 0.9591 -3.1293 -0.3604 -0.4251 0.4185 -0.6728 -0.0027 2.1145 1.1157y(k) 0.0615 0.9120 -0.0692 -3.2731 3.7486 -4.3194 4.7230 -5.2781 5.1507 -2.7235一、最小二乘法（LS）1、数学模型设时不变SISO动态过程的数学模型为 （）其

4、中，为过程的输入量，为过程的输出量，是噪声，多项式和为：在本题中，=3.即将此模型写成最小二乘格式 （）其中，是过程的输出量；是可观测的数据向量；是均值为零的随机噪声。式中对于，方程式（）构成一个线性方程组，可以把它写成（）利用数据序列和，极小化准则函数 （）使的估计值记作，称为参数的最小二乘估计值。通过极小化（）式来计算的方法称作最小二乘法，未知模型参数最可能的值是在实际观测值与计算值之累次误差的平方和达到最小处所得到的，这种模型输出能最好地接近实际过程的输出。2、辨识原理考虑模型（）式的辨识问题，其中和都是可观测的数据，是待估计参数，准则函数取（）根据（）的定义，准则函数可写成二次型的形式

5、 （）显然上式中的代表模型的输出，或者说是过程的输出预报值。因此可以看作来衡量模型输出与实际过程输出的接近情况。极小化，求得参数的估计值 （）将使模型的输出最好的预报过程的输出。3、 辨识结果二、递推最小二乘法（RLS）1、数学模型在第一部分中建立了最小二乘法，并用一次完成算法进行了计算。但是由于具体使用时占用内存量大，而且不能用于在线辨识。所以引入了最小二乘参数估计的递推算法。递推算法的基本思想可以概括成：新的估计值=老的估计值+修正项 （）在此算法中，时不变SISO动态过程的数学模型仍与最小二乘法的一样。模型的最小二乘格式也相同，只是计算方法不同，具体计算方法，在递推最小二乘法的辨识原理中

6、描述。2、辨识原理首先将式一次完成算法写成 （）定义 （）其中：，由（）式可得：（）令：则：于是有 （）令利用（）和（）式，可得 （）引进增益矩阵，定义为 （）则（）式写成 （）进一步把（）式写成 （）为了避免矩阵求逆运算，利用矩阵反演公式可将（）式演变成 （）将（）式代入（）式，整理后有 （）综合（）、（）、（）式便得到最小二乘参数估计递推算法。利用上述公式即可求得参数的估计值3、辨识结果三、广义最小二乘法（GLS）1、数学模型设时不变SISO动态过程的数学模型为 （）其中，和 分别表示过程的输入和输出，是均值为零的不相关随机噪声，多项式、和为：在本题中，在本题中，=3.即 （）广义最小二乘

7、发的基本思想是基于对数据先进行一次滤波预处理，然后利用普通最小二乘法对滤波后的数据进行辨识。2、辨识原理由（）式可得 （）令， （）及（）可将模型（）式化成最小二乘格式 （）由于上式是白噪声，所以利用最小二乘法即可获得参数的无偏估计。但是，数据向量中的变量均需要按照（）式计算，然而噪声模型并不知道。为此需要用迭代的方法来估计。令 （）置 （）就把噪声模型（）也化成最小二乘格式由于上式的噪声已是白噪声，所以再次利用最小二乘法可获得噪声模型参数的无偏估计。但是，数据向量依然包含着不可测得噪声量，它可用相应的估计值代替，置 （）其中，当时，按照式子计算，式子中综上分析，广义最小二乘递推算法可归纳成：

8、利用上述公式即可求得参数的估计值并能求出噪声模型的估计。3、估计结果噪声模型的估计值为1.0e-003 * -0.0000 0.4115 0.0001 即附录：MATLAB程序1、一次完成的最小二乘法%一次完成的最小二乘法clears=(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

9、 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.13

10、42 -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:40h(i,:)=-z(i) -z(i-1) -z(i-2) u(i) u(i-1) u(i-2);endo=inv(h'*h)*h'*z'%误差分析J

11、=0;for i=1:1:40e(i)=z(i)-h(i,:)*o;J=J+e(i)2;endJplot(s,e);2、递推的最小二乘法%递推的最小二乘法clears=(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.8

12、334 -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.

13、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:40h(i,:)=-z(i) -z(i-1) -z(i-2) u(i) u(i-1) u(i-2);endp0=106;o0=0.001;k(1,:)=p0*h(1,:)'*inv(h(1,:)*p0*h(1,:)'+1)

14、9;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)

15、-h(i,:)*o(i-1,:)''end%辨识结果o(40,:)%误差分析J=0;for i=2:1:40e(i)=z(i)-h(i,:)*o(i-1,:)'J=J+e(i)2;endplot(s,e);J3、 广义最小二乘法%广义预测的最小二乘法clears=(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.0

16、786.-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

17、.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=106;oe0=0;pe0=1;%k=1时kf(1,:)=pf0*hf(1,:)

18、'*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(:,

19、:,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)

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

21、*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,:)