不需矩阵求逆的最小二乘法

1、一种不需矩阵求逆的最小二乘法本文着重介绍了一种新的不需要矩阵求逆的最小二乘法。在 LSLS 估计中，需要计算矩阵( (T) )的逆，由于矩阵，求逆比较影响辨识速度。能否在 LSLS 中避免( (T) )的求逆运算？从而提高 LSLS 的辨识速度。关键字：LSLS 估计求逆运算辨识速度最小二乘法1 问题提出在LS估计中，需要计算矩阵()的逆，由于矩阵，求逆比较影响辨识速度。能否在LS中避免()的求逆运算？从而提高LS的辨识速度。(1)不需矩阵求逆的LS算法特点不需计算矩阵的逆；辨识精度与基本LS相同；辨识速度比基本LS有较大提高；适合模型阶次n未知的情况下应用；是一种按模型阶次n递推的算法。(2

2、)算法推导系统的差分方程为：y(k)a1y(k-1)any(k-n)=t0u(k)tu(k-1)0u(k-n)(k)二y(k尸u(k)-31y(k-1)tu(k-1)-any(k-n)bnu(k-n)(k)上式写成矩阵形式有：y(k)=(k)M(k)式中：(k)=u(k)y(k-1)u(k-1)y(k-n)u(k-n)】b=bo-31b1-anbnT令k分别等于1,2n,有n个方程，则有：Y=8+E式中：丫=勺(1)y(2)y(N)1T底1-(1)(2)(N)1T珑(1)1u(1)|6T(2)u(2)：一：.3T(N)jMN)y(0)u(0)y(i)u(i)a+y(N-1)u(N-1)y(1-

3、n)u(1-n)y(2-n)u(2-n)aay(N-n)u(N-n)-0。1。22nJ62n2 任务提出由u(k)及y(k)辨识ai、bi及n3 采用方法采用模型阶次n的递推方法。即从n=0开始辨识一n=1一n=2一由Y=什己可得其LS估计值为：鼻文飞尸Y下面推导按模型阶次n的递推算法。模型阶次n=0的辨识o=L(1)u(2)u(N)1T=乳N-N=?o=gT0)T丫=Eu(i)y(i)/Zu2(i)i=1.i工引入中间变量X:T1X=(犷)记n=0时的X为X:*0=(中人0)=1u2(i)(2)求X(:J，假设模型阶次为n-1时的参数辨识结果已知，也就是？n，Xn已知，现欲求模型阶次为n时的

4、辨识结果。即求：小：1：Y这里，我们先求Xn。由可知:(=卬0巾1巾2巾2nd.巾2nn-n1巾2n1巾2nl式中：巾2n-y(1-n)y(2-n)y(N-n)1T巾2n=U(1-n)u(2-n)u(N-n)T令n=Ln_L巾2n_J,则：式中：A22=1/(巾2n巾2n一巾2nnXnn巾2n)A12=A。=xnnT巾2nA22A11=X；-A124n:XnT至此，Xn=(nTn)(nig=j/匕n-2nJ.-一二d士n 士niiT不2nd.T2nJ2nd根据分块矩阵求逆可得:Xn=(；)尸=B1121B12B22式中:B22=1/（巾2nJ巾2n一巾2nn1Xnn巾2n）B12=B21=Xn

5、J.n巾2n4B22B11-B12巾同理，令n=1n巾2nL则可推得:已求得。求解步骤：；=In,LX；=(；)T；产 Tn=；W2nLX0Q；口充求丫：丫=jnJ“Tn丫求?n2n=(Tnn)TnY=XnTnY4 总结。,X,2T1,XI，T2,X2，221通过与矩P$求逆的LS算法对比可得出： 不需矩阵求逆的LS算法优势是不需计算矩阵的逆，简化计算，减少计算量，并且辨识精度与基本LS相同，且辨识速度比基本LS有较大提高和适合模型阶次n未知的情况下应用， 是一种按模型阶次n递推的算法。:丫N工y(i-n)y(i)i=1N工u(i-n)y(i)i4动态模型参数极大似然辨识及其MTLAER现极大

6、似然法辨识方法特点：(1)无偏估计方法；(2)适用于己(k)相关情况；(3)当系统信噪比较小时有较好的估计效果；(4)算法稳定度好；(5)是一种递推算法；(6)实际工程中广泛使用。设动态系统的模型表示为11A(z)z(k)=B(z)u(k)+e(k)-e(k)=D(z)v(k)式中，v(k)是均值为0,方差为仃2,服从正态分布的不相关随机噪声；u(k)和z(k)表示系统的输入输出变量。现给出一系统模型为z(k)-1.2z(k-1)+0.6z(k-2)=u(k-1)+0.5(k-2)+e(k)e(k)=v(k)-v(k-1)+0.2v(k-2)其中v(k)为随机信号，输入信号是幅值为1的M系列或

7、随机信号，试用递推的极大似然法求系统辨识的参数仃0程序如下：cleara(1)=1;b(1)=0;d(1)=0;u(1)=d(1);z(1)=0;z(2)=0;%初始化fori=2:1200%产生m序歹Uu(i)a(i)=xor(c(i-1),d(i-1);b(i)=a(i-i);c(i)=b(i-1);d(i)=c(i-1);u(i)=d(i);Endu；v=randn(1200,1);%产生正态分布随机数V=0;%计算噪声方差fori=1:1200V=V+v(i)*v(i);endV1=V/1200;fork=3:1200%根据v和u计算zz(k)=1.2*z(k-1)-0.6*z(k-2

8、)+u(k-1)+0.5*u(k-2)+v(k)-v(k-1)+0.2*v(k-2);endo1=0.001*ones(6,1);p0=eye(6,6);%幅初值zf(1)=0.1;zf(2)=0.1;vf(2)=0.1;vf(1)=0.1;uf(2)=0.1;uf(1)=0.1;%迭代计算参数值和误差值fork=3:1200h=-z(k-1);-z(k-2);u(k-1);u(k-2);v(k-1);v(k-2);hf=h;K=p0*hf*inv(hf*p0*hf+1);p=eye(6,6)-K*hf*p0;v(k)=z(k)-h*o1;o=o1+K*v(k);p0=p;o1=o;a2(k)

9、=o(2);b1(k)=o(3);b2(k)=o(4);d1(k)=o(5);d2(k)=o(6);e1(k)=abs(a1(k)+1.2);e2(k)=abs(a2(k)-0.6);e3(k)=abs(b1(k)-1.0);e4(k)=abs(b2(k)-0.5);e5(k)=abs(d1(k)+1.0);e6(k)=abs(d2(k)-0.2);zf(k)=z(k)-d1(k)*zf(k-1)-d2(k)*zf(k-2);uf(k)=u(k)-d1(k)*uf(k-1)-d2(k)*uf(k-2);vf(k)=v(k)-d1(k)*vf(k-1)-d2(k)*vf(k-2);hf=-zf(k-1);-zf(k-2);uf(k-1);uf(k-2);vf(k-1);vf(k-2);endo1V1运行几次选择了一组比较接近真值的结果，如下所示。o1=-1.18920.56390.95480.4505-0.95000.1873V1=0.9849可知，估计值a1=