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

12008 级硕士研究生系统建模理论试卷已知一个三阶线性离散系统的输入、输出数据，共有 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 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.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 动态过程的数学模型为（1.1.1）)()()(11 knuzBkzA其中， 为过程的输入量， 为过程的输出量， 是噪声，多项式 和)(ku )(1zA为：1zB nbzzbaaAa 211)(在本题中， = =3.即an3211)(zbzbBaA将此模型写成最小二乘格式 （1.1.2）)()(knkh2其中， 是过程的输出量； 是可观测的数据向量； 是均值为零的随机噪声。)(kz)(kh)(kn式中 , )3(,2)1(3)2(),1()( 21baukuzh对于 ，方程式（ 1.1.2）构成一个线性方程组，可以把它写成Lk,2（1.1.3 )3()2()1()3()2()1( 00)()0()(2)1( ,)2(1 lullulzlzlzzl lnnzzl l h）利用数据序列 和 ，极小化准则函数 （1.1.4）)(kz)(hLkkzJ12)()(h使 的 估计值记作 ，称为参数 的最小二乘估计值。通过极小化（1.1.4）式min)(J来计算 的方法称作最小二乘法，未知模型参数 最可能的值是在实际观测值与计算值之累次误差的平方和达到最小处所得到的，这种模型输出能最好地接近实际过程的输出。2、辨识原理考虑模型（1.1.2）式的辨识问题，其中 和 都是可观测的数据， 是待估计参数，)(kzh准则函数取（1.1.4）根据（1.1.3）的定义，准则函数 可写成二次型的形式)(J（1.2.1）)()()HzllllJ显然上式中的 代表模型的输出，或者说是过程的输出预报值。因此 可以看作来l )(J衡量模型输出与实际过程输出的接近情况。极小化 ，求得参数 的估计值)(J（1.2.2）llz1)(将使模型的输出最好的预报过程的输出。3、辨识结果 000.1.- 321321 ba3二、递推最小二乘法（RLS）1、数学模型在第一部分中建立了最小二乘法，并用一次完成算法进行了计算。但是由于具体使用时占用内存量大，而且不能用于在线辨识。所以引入了最小二乘参数估计的递推算法。递推算法的基本思想可以概括成：新的估计值 =老的估计值 +修正项 （2.1.1）)(k)1(k在此算法中，时不变 SISO 动态过程的数学模型仍与最小二乘法的一样。模型的最小二乘格式也相同，只是计算方法不同，具体计算方法，在递推最小二乘法的辨识原理中描述。2、辨识原理首先将 1.2.2 式一次完成算法写成（2.2.1）)()()()( 111 lililll izzz hhHPH 定义（2.2.2）111)()(ki kiiPh其中： ，)(2kkhH)1(2)(1kkh由（2.2.2）式可得：（2.2.3）)()()()()( 111 kiki hPP令： ,2,1kzzkz则： )()1()() 11kikk izhH于是有 （2.2.4）11 )()(iizhP令 ,2),kzkz利用（2.2.3）和（2.2.4）式，可得4（2.2.5）)1()()(1( )1)()()()1 kkzk zizkkikk hpPzH引进增益矩阵 ，定义为K（2.2.6）)()(khP则（2.2.5）式写成（2.2.7）)1()()1() kzkh进一步把（2.2.3）式写成（2.2.8 ）11)()()( kP为了避免矩阵求逆运算，利用矩阵反演公式可将（2.2.8）式演变成（2.2.9）)1()()()()( 11 kk PhKIh将（2.2.9）式代入（2.2.6）式，整理后有（2.2.10 ）1)()()() kkPPK综合（2.2.7） 、 （2.2.9） 、 （2.2.10）式便得到最小二乘参数估计递推算法。 )()()1()zkh )1(1kkPhKI1)()()() kPPK利用上述公式即可求得参数 的估计值3、辨识结果 0.1-.5-02-0.7-1.5- 321321 ba三、广义最小二乘法（GLS）1、数学模型设时不变 SISO 动态过程的数学模型为（3.1.1）)(1)()(11 kvzCkuzBkzA其中， 和 分别表示过程的输入和输出， 是均值为零的不相关随机噪声，多u )(kv5项式 、 和 为：)(1zA)(1B)(1zC cbannzzczCba 2112)(在本题中，在本题中， = = =3.即anbc（3.1.2 ） 3211 3211)( zczczCbBA广义最小二乘发的基本思想是基于对数据先进行一次滤波预处理，然后利用普通最小二乘法对滤波后的数据进行辨识。2、辨识原理由（3.1.1）式可得 （3.2.1）)()()()( 11 kvuzCBkzCzA令 ， （3.2.2）)()(1kukuzf及（3.2.3）, )3(),2(),1()3()2(),()( 211bakukuzzz fffffff h可将模型（3.1.1）式化成最小二乘格式 （3.2.4）)()(vff h由于上式 是白噪声，所以利用最小二乘法即可获得参数 的无偏估计。但是，数据向)(kv量 中的变量均需要按照（3.2.2）式计算，然而噪声模型 并不知道。为此需fh )(1zC要用迭代的方法来估计 。)(1zC令 （3.2.5）)(1)kvzke置 （3.2.6）,)()3()231ckkeeeh就把噪声模型（3.2.5）也化成最小二乘格式 )()(kveh由于上式的噪声已是白噪声，所以再次利用最小二乘法可获得噪声模型参数 的无偏估计。e但是，数据向量 依然包含着不可测得噪声量 ，它可用相应)(keh )3(,2),1(kek6的估计值代替，置 （3.2.7）)3(),2(),1() kekeeh其中 ，当 时，按照式子 计算，0,)(ke)(khz式子中 3,2(),1,(),(),(ukuzkz综上分析，广义最小二乘递推算法可归纳成： )()()1() kfff hK1kffff PIP 1)()()() kkffffff 1( keeee h)1()()kPKIP 1)( kkeeeeeh利用上述公式即可求得参数 的估计值 并能求出噪声模型的估计。3、估计结果 123123-1.05-.80.2-.0-.5-0.1ab噪声模型的估计值为 1.0e-003 * -0.0000 0.4115 0.0001 即 1()Cz7附录：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 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: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=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.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.8-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: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);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:40k(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: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.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;9u=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=106;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