第四章线性系统参数估计的最小二乘法.pdf
- 文档编号:3439522
- 上传时间:2023-05-05
- 格式:PDF
- 页数:16
- 大小:397.52KB
第四章线性系统参数估计的最小二乘法.pdf
《第四章线性系统参数估计的最小二乘法.pdf》由会员分享,可在线阅读,更多相关《第四章线性系统参数估计的最小二乘法.pdf(16页珍藏版)》请在冰点文库上搜索。
第四章第四章线性系统参数估计的最小二乘法线性系统参数估计的最小二乘法4.1引言引言最小二乘(LeastSquares)法是用于参数估计的数学方法,它使数学模型在误差平方和最小的意义上拟合实验数据。
1801年,意大利天文学家朱赛普皮亚齐发现了第一颗小行星谷神星。
经过40天的跟踪观测后,由于谷神星运行至太阳背后,使得皮亚齐失去了谷神星的位置。
随后全世界的科学家利用皮亚齐的观测数据开始寻找谷神星,但是根据大多数人计算的结果来寻找谷神星都没有结果。
时年24岁的高斯也计算了谷神星的轨道。
奥地利天文学家海因里希奥尔伯斯根据高斯计算出来的轨道重新发现了谷神星。
高斯使用的最小二乘法的方法发表于1809年他的著作天体运动论中,而法国科学家勒让德于1806年独立发现“最小二乘法”,但因不为时人所知而默默无闻。
勒让德曾与高斯为谁最早创立最小二乘法原理发生争执。
1829年,高斯提供了最小二乘法的优化效果强于其他方法的证明,因此被称为高斯-马尔可夫定理。
此后LS法成了处理观测所得实验数据的有力工具。
现在成为辨识的主要算法。
最小二乘辨识是一种一致的、无偏的、有效的方法4.2基本最小二乘方法基本最小二乘方法4.2.1最小二乘曲线拟合4.2.1最小二乘曲线拟合最小二乘技术提供给我们一个数学程序,通过它能获得一个在最小方差意义上与实验数据拟合最好的模型。
假定有一个变量Y,它与一个n维变量X=(x1,x2,xn)是线性关系,即Y=1x1+2x2+nxn(4.1)其中,=(1,2,n)是一个参数集。
在系统辨识中它们是未知的。
我们希望通过不同时刻对Y及X的观测值来估计出它们的数值。
例如,在研究两个变量(x,y)之间的关系时,通常的做法是取一个变量作为自变量,另一个作为因变量。
改变自变量可得到相应的因变量。
将所得到的一系列数据对描绘在直角坐标系中,得到一系列的点。
若这些点落在一条直线附近,就可以用一个直线方程来描述这个关系。
设已得到的一系列点分别为(1,1.8),(2,2.2),(3,3),(4,2.8),(5,3.211.522.533.544.555.561.522.533.54(1,1.8)(2,2.2)(3,3.3)(4,2.8)(5,3.2)(6,3.3),(6,3.3),直线方程为y=a0+a1x使用(1,1.8),(2,2.2)两个点得到的方程为y=1.4+0.4x;使用(1,1.8),(6,3.3)两个点得到的方程为y=1.5+0.3x,而使用(3,3)和(6,3.3)两个点得到的方程是y=2.7+0.1x。
从图中可看到,前两条线都仅能满足两个点的要求,而对其它点的误差都很大,其6个点的误差平方累计分别为0.49和0.42。
第三条线能满足三个点的要求,但误差平方累计更大,为1.58。
显然我们需要找到一条更为理想的直线来取得较小的误差。
例如图中的红色短划线,它的方程为y=1.697+0.294x,误差平方累计为0.25。
这条线是怎样得到的呢?
它是用最小二乘法得到的。
下面讨论更为一般的情况。
假设在t1,t2,tm时刻对Y及X的观测值序列已经被我们获得,并且用y(i),x1(i),x2(i),x3(i),i=1,2,m来表示这些观测数据。
显然,可以用m个方程组来表示量测数据与估计值之间的关系(4.2)+=+=+=)()()()()2()2()2()2()1()1()1()1(221122112211mxmxmxmyxxxyxxxynnnnnnLMLL写成矩阵形式,为Y=X(4.3)其中1)()1(=mmyyYMnmnnmxmxxxX=)()()1()1(11LMML11=nnM根据m相对n的大小不同,可能有三种情况:
1)mn如果无任何测量误差和噪声,也应有唯一解。
但它们是不可避免的,这样就可能造成无解的情况。
例:
对于系统2211xxy+=,采样得任取其中2个方程都可解得212121373241221;若有测量噪声或误差,则无解。
212121301.7324995.0005.2;12;-2-因为此时实际方程变为+=XY,在这种情况下只能求近似解。
通常使用的方法为最小平方误差法4.2.2最小二乘估计准则和正则方程4.2.2最小二乘估计准则和正则方程已知+=XY(4.4)则误差应为XY(4.5)设目标函数或估计的判据为误差的平方=miiTJ12(4.6)取J最小为优化指标XXXYYXYYXYXYJTTTTTTT+=)()(2121998.0004.221;将J对微分并令其为零(关于矩阵微分的方法见附录),有022=XXYXJTTYXXXTT=(4.7)得(4.8)YXXXTT1)(这样求得的就称为的最小二乘估计(LSE),在统计学上,方程(4.7)称为正则方程,称为残差。
在前面讨论的例子中,把6个数据对分别代入直线方程y=a0+a1x中可得到1个由6个直线方程构成的方程组XY=,其中各参数分别为;=10aaTY3.32.38.232.28.1=;。
将这些参数分别代入式(4.8)解出TX=654321111111294.0;687.110=aa,并可求出误差向量XY=。
误差的平方和为253.0*=例例4.1测得铜导线在温度时的电阻)(CTio)(iR如表6-1,求电阻R与温度T的近似函数关系。
i1234567)(CTio19.125.030.136.040.045.150.0)(iR76.3077.8079.2580.8082.3583.9085.10解:
解:
画出散点图,可见测得的数据接近一条直线,故取n=1,拟合函数为TaaR10+=于是可列出方程组-3-+=+=+=)7()7()2()2()1()1(101010TaaRTaaRTaaRM写成最小二乘标准格式为+=XY,其中TRY10.8590.8335.8280.8025.7980.7730.76=TX=0.501.450.400.361.300.251.191111111=10aa解方程组得YXXXTT1)(=572.700=a292.01=a故得R与T的拟合直线为TR292.0572.70+=利用上述关系式,可以预测不同温度时铜导线的电阻值。
例如,由R=0得T=-242.5,即预测温度T=-242.5时,铜导线无电阻。
4.2.3最小二乘估计的几何解释4.2.3最小二乘估计的几何解释Yx1xmY已知最小二乘输出估计向量;输出估计误差向量如图(输出残差向量)XY=YY=;输入测量矩阵=TmTTxxxXM21输出估计向量Y是输出测量Y在由x1,x2,xm所张成的空间上的正交投影,或者说,输出残差向量垂直于由x1,x2,xm所张成的空间。
这个几何解释说明,输出向量Y可以分解成属于估计空间的估计向量Y和垂直于估计空间的输出残差向量。
最小二乘估计的解的形式为最小二乘估计的解的形式为(4.9)YXXXTT1)(=很显然,只有正则方程的逆存在才能有解。
而正则方程的逆存在(或者说最小二乘估计的解存在)的充分必要条件为:
过程的输入信号必须是2n阶持续激励信号。
这个条件又被称为开环可辩识条件。
它意味着辩识所用的输入信号不能任意选择,否则就可能造成系统不可辩识。
通常所用的信号有
(1)随机序列(如白噪声);
(2)伪随机序列(如M序列);(3)离散序列(对含有n种其频率不是倍频关系的正弦信号的组合信号进行采样所得到的离散序列);4.3最小二乘在线性系统参数估计中的应用最小二乘在线性系统参数估计中的应用(差分方程辩识)(差分方程辩识)4.3.1线性系统模型的参数估计4.3.1线性系统模型的参数估计在连续系统中,系统的动态特性可用微分方程来描述,在离散系统中,系统的动态特性则用差分方程来描述。
对于n阶系统,它的差分方程可写为(4.10)()()()(11kuqBkyqA=传递函数可表示为)()()()()(111=qAqBkukyqH(4.11)其中A,B分别是引入向后移位算子q-1的多项式(4.12)nnqaqaqA+=L1111)(nnqbqbbqB+=L1101)(式(4.10)可写为)()1()()()1()(101nkubkubkubnkyakyakynn+=+LL-4-对有测量误差的系统,则有(4.13)()()()()(11kkuqBkyqA+=通过移项,还可写为)()()()(01kikubikyakyniinii+=)()(kkXT+=(4.14)当观测得到的输出输入数据为2n+1个时(也就是n+1次观测后),可写出)(,),(),(,),1(nkukunkykyXT=LL(4.15)这时被辩识参数。
当对输入输出观测了N+n次后,就可建立起有N个方程的方程组(N2n),即,101nnTbbbaaLL=+=XY(4.16);+=)()2()1(NnynynyYM+=)()2()1(NnnnM;+=)()2()1(NnXnXnXXTTTM+=)()()()1()2()2()2()1()1()1()1()(NuNnuNyNnyunuynyunuynyLLMMMMLLLL(4.17)这样我们就可以用前面介绍的最小二乘法求得估计量YXXXTT1)(=例例4.2已知一个二阶系统的传递函数为22111101)(+=zazazbbzH,在其输入端加入M序列输入后所得到的输出输入数据见下表,请利用这些数据辨识出系统的传递函数的系数。
k12345678910输入u1011001110输出y-0.45-0.011.152.561.92-0.30-0.800.912.922.40解:
解:
已知系统阶数n=2,有4个未知数。
将式(4.4)展开)1()()2()1()(1021+=kubkubkyakyaky根据要求,观测次数N2n+1,取N为6,k=3-5-=91.08.03.092.156.215.1)8()7()6()5()4()3(yyyyyyY=113.08.00192.13.00056.292.11015.156.21101.015.10145.001.0)7()8()6()7()6()7()5()6()5()6()4()5()4()5()3()4()3()4()2()3()2()3()1()2(uuyyuuyyuuyyuuyyuuyyuuyyX=57.089.076.088.0)(11021YXXXbbaaTT21176.088.0157.089.0)(+=zzzzH当N取8时有21177.088.0156.090.0)(+=zzzzH;当N取80时有21178.088.0157.094.0)(+=zzzzH而其真值为2118.09.015.01)(+=zzzzH例例4.3如图所示的一个差分方程系统,已知A为二阶,B为零阶,结构如下:
B(z-1)u(t)y(t)w(t)v(t)+A(z-1)=+=1)
(1)(122111zBzazazA其200次测量的最后10次观测数据为观测次数k191192193194195196197198199200观测输出Y0.580.530.550.610.550.530.550.550.590.54系统输入U0.790.770.760.840.780.760.760.780.810.77请用最小二乘辨识一次完成算法辨识多项式A的系数。
解:
解:
已知系统的传递函数为221111111)()()()()(+=zazazAzBkukyzH将其展开后为)()2()1()(21kukyakyaky=+;写成标准最小二乘格式后成为)()2()1()(21kukyakyaky+=由此可知待辨识参数是这里b021baaT=0=1。
由于b0已知,所以可将待辨识参数简化为,这样就将待辨识系统由3阶降为2阶,这时标准最小二乘格式变为21aaT=)2()1()()(21=kyakyakuky。
取6组观测数据-6-)196()197()198()198()195()196()197()197()194()195()196()196()192()193()195()195()192()193()194()194()191()192()193()193(212121212121yayauyyayauyyayauyyayauyyayauyyayauy=那么Y和X矩阵则分别为=23.021.023.023.023.021.078.055.076.055.076.053.078.055.084.061.076.055.0)198()198()197()197()196()196()195()195()194()194()193()193(uyuyuyuyuyuyY;=53.055.055.053.061.055.055.061.053.055.058.053.0)196()197()195()196()194()195()193()194()192()193()191()192(yyyyyyyyyyyyX=1.85731.85291.85291.8414XXT=92.3592.93-92.93-94.05.)(1XXT=0.7479-0.7424-YXT=0.0797-0.3229-)(211aaYXXXTT4.3.2稳态系统4.3.2稳态系统对于稳态系统,应有z1或q1。
那么对于单输入单输出系统其方程简化为)()()1(0kkubky+=+其实这里只是要求辨识出系统的比例系数,无较大的实用价值。
因而这里将讨论多输入单输出(MISO)及多输入多输出(MIMO)系统的稳态模型的辨识4.3.2.1最小二乘解最小二乘解稳态MISO系统表达式为)()()()(110kkubkubbkynn+=L(4.18)进行了N(Nn)次观测后,可建立起观测方程组TNyyY)()1(L=;TnbbL0=TN)()1(L=)()
(1)2()2
(1)1()1(1111NuNuuuuuXnnnLMMLL(4.19)于是可以形成最小二乘标准方程式+=XY可以解得YXXXTT1)(=-7-4.3.2.2从从MISO系统至系统至MIMO系统的推广系统的推广如,则有mrRURY;+=)()()()()
(1)()()(21110221201111021kkkkukubbbbbbbbbkykykyrmmrrrmmrMMLMMLLM(4.20)可先分解成r个形同式(4.18)的具有m个输入的单输出系统,然后用上述同样的方法求最小二乘解。
4.3.2.3至非线性系统的推广至非线性系统的推广只要能使输出与被辨识参数关系描述为线性方程,非线性系统参数估计也可用LS方法,例如)()()()(2210kkuakuaaky+=(4.21)可令x1(k)=u(k);x2(k)=u2(k),即可得到)()()()(22110kkxakxaaky+=(4.22)这样就可用LS方法进行估计。
4.3.3自回归(AR)模型的参数估计4.3.3自回归(AR)模型的参数估计考虑AR模型)()()(1kikyakynii=+=(4.23)可写成)()()(1kikyakynii+=(4.24)对于Nn次观测可建立方程+=XY,其中+=NnynyY()1(M;=naaM1+=)()1()2()1()1()()()2()1(NyNnyynyynyNXXXXTTTLMMLLM这样就可利用前面所学过的求解正则方程的方法辨识出AR模型的参数。
-8-9-4.3.4系统阶数与辨识参数个数的关系4.3.4系统阶数与辨识参数个数的关系前面的公式推导,是建立在系统传递函数的分母阶数na与分子阶数nb相同的基础上,即na=nb,但在实际系统中往往两者不相等,因而在实际的辨识中,如果系统的结构(或未知参数的个数)已知,则按已知的参数的位置组织数据进行辨识,如果仅知系统的阶数na,则应假设na=nbn,进行全部参数的辨识。
4.4递推最小二乘递推最小二乘(RLS)4.4.1算法推导4.4.1算法推导设在测量了m次后的观测数据为Ym,Xm。
所估计出的参数为,则有)(mmmXY=;mTmmTmYXXXm1)()(=现在,又进行了一次观测,取得了m+1时刻的数据y(m+1),x1(m+1),x2(m+1),xn(m+1),这时就要根据新的m+1个方程来对参数进行估计11mmXY=(4.25)新的LS估计应为(4.26)11111)()1(+=+mTmmTmYXXXm由此可见,在辨识中,当一组新的测量数据被引入后,为了获得新的辨识结果,需要将新的数据与过去的旧数据一起代入公式,重新估计。
更为重要的是,由于所有的数据都要参加运算,随着测量次数的增加,计算机所要存储的数据量也在不断的增加。
而且在最小二乘估计的运算中有一个矩阵求逆运算,当观测次数很大时,耗时巨大。
考虑到在前此次估计时把所有输出输入数据的信息都包含在了运算结果中,这次不过是仅增加了一组数据,完全可以利用前次估计的结果,添加上最新的数据信息,形成新的估计值。
由此将新的输出输入矩阵按新老数据进行分解:
+=+=+)1()1()()1(1myYmymyyYmmM;+=+=+)1()1()1()()()1()1(1111mXXmxmxmxmxxxXTmnnnmLLMML这样就找出了Ym+1与Ym以及Xm+1与Xm之间的关系,现在可以进一步寻找与之间的关系。
在寻找这个关系时,要用到下面这个矩阵运算恒等式:
)1(+m)(m如果如果A、C和和A+BCD为非奇异方阵,则有为非奇异方阵,则有(4.27)1111111)()(+=+DABDACBAABCDA为了简化推导运算,定义一个辅助矩阵P(m)(4.28)1)()(=mTmXXmP则(4.29)111)()1(+=mTmXXmP1)1()1(+=mXXmXXTmTm()1)1()1(+=mXmXXXTmTm()11)1()1()(+=mXmXmPT)()1()1()()1
(1)(1()()(1mPmXmXmPmXmXmPmPTT+=-10-最后一个等式用到了矩阵恒等式,其中A=PP-1(m),BX(m+1),CI,DXT(m+1)。
记,应注意,是一个标量!
则:
=+1)1()()1(1(mXmPmXT)()1()1()()()1(mPmXmXmPmPmPT+=+(4.30)11)1()1(+=+mTmYXmPm)1()1()1(+=mymXYXmPmTm+=mTmTmTmYXmPmXmXmPYXmP)()1()1()()()1()1()()1()1()()1()1()(+mymXmPmXmXmPmymXmPT上式最后一个等式右边最后两项可写为)1()1()()1()1()()1()1()(+mymXmPmXmXmPmymXmPT)1()1()()1
(1)1()(+mymXmPmXmXmPT)1()1()()1()1()(1+mymXmPmXmXmPT)1()1()()1()1()()1
(1)1()(+mymXmPmXmXmPmXmXmPTT)1()1()(+mymXmP而最后一个等式右边第一项为,所以)()(mYXmPmTm)1()1()()()1()1()()()1(+=+mymXmPmmXmXmPmmT)()1()1()1()()(mmXmymXmPmT+=(4.31)()1()1()1()(mmXmymKmT+=其中)1()()1(+=+mXmPmK(4.32)1()()1(1/)1()(+=mXmPmXmXmPT由式(4.31)可以看出,本次的估计值是由上次的估计值加上本次修正项构成。
而本次修正项中,表示由前m次测量数据的估计值来拟合本次采样所得数据所产生的误差,称为预测。
而K(m+1)可认为是对预测误差进行修正时的加权因子。
)()1()1(mmXmyT+)(m4.4.2递推公式4.4.2递推公式由式(4.30),(4.31)和(4.32)可归纳出递推最小二乘估计的递推公式(4.33a)1()()1(1/)1()()1(+=+mXmPmXmXmPmKT(4.33b)()1()1()()1(mPmXmKmPmPT+=+-11-(4.33c)()1()1()1()()1(mmXmymKmmT+=+初始值可选0或1,P(0)选为)0(IP=)0(,其中是一个很大的正数。
这样我们就可以用(4.33)这组公式开始迭代,直到求出最后需要的结果为止。
4.4.3遗忘因子法(适应算法)4.4.3遗忘因子法(适应算法)4.4.3.1“数据饱和”现象由式(4.33)可知,第m+1次的估值是在的基础上,依据新息)1(+m)(m1(预报误差)与增益矩阵)()1()1()1(mmXmymyT+=+)1(+mK的乘积来修正的。
K值越大,修正能力越强。
如果K趋近于0,算法就失去了修正能力。
由P的定义可知,PP-1是正定的,P也是正定的,那么必有(4.35)0)1()()1(1+mXmPmXT(4.36)0)1()()1()(+TmXmPmXmP由此可得0)1()()1
(1)1()()1()(+mXmPmXmXmPmXmPTT(4.37)由P的递推式(4.33b)可得(4.38)1()()()1()1(+=+mPmPmPmXmKT式(4.37)表明了式(4.38)的左边是大于零的,故可知式(4.38)的右边也是大于零的0)1()(+mPmP(4.39)这意味着P(m)是递减的正定矩阵。
当m时,P(m)0。
就是说P(m)将逐渐变成奇异的。
但因P(m)的秩不会变化,只能是P的各元素0,所以K也将随着m增加趋近于零,使RLS算法的修正能力失效。
适应算法就是针对“数据饱和”现象发展起来的辩识方法,它包括有遗忘因子法和限定记忆法。
我们在这里只介绍遗忘因子法。
4.4.3.2加权最小二乘估计(WLSE)4.4.3.2加权最小二乘估计(WLSE)加权的必要性对各次测量值,因条件不同,对各次误差的重视程度不同,同时各量测量的精度也不同,因而需进行相应的处理。
性能指标。
写成矩阵形式为2222211mmwwwwJ+=L)()(TTXYWXYWJw=W为对角矩阵,称为权矩阵。
对Jw求的极小化,得到的加权最小二乘估计wWYXWXXwT1T)(=1新息:
(预报误差))()1()1()1(mmXmymyT+=+残差:
(估计误差))1()1()1()1(+=+mmXmymeT-12-当W为I矩阵时,上式蜕变为普通最小二乘。
对W的限定:
对称正定阵对称正定阵一个一个n阶对称矩阵,若对于任何阶对称矩阵,若对于任何n维向量维向量A,(A0),都有,都有ATWA0,则称,则称W为正定的。
如果为正定的。
如果W是正定的,是正定的,X是满秩的,则是满秩的,则XTX的逆存在,即的逆存在,即(XTWX)-1存在。
存在。
若W各元素按时间顺序m-i01i=1,2,m进行选择,就称为遗忘因子法。
4.4.3.3遗忘因子法遗忘因子法遗忘因子法的目标函数为)()1()(202112miJmmTmmiimmL+=(4.40)这里=)()1(1mmmmmM。
已知估计误差为mmmXY=(参见式(4.5))。
参数的第m次最小二乘估计为(4.42)mTmmTmmYXXX1)(=在第m+1估计时,由于+=+=1121)(miimmiJ)1()()1(2212+=mmmL)1()(212+=mimiim)1(2+=mmTm+=)1()1(mmmTm11+=mTm(4.43)将式(4.41)两边同乘得mmmXY=(4.44)并且(4.45)1()1()1(+=+mXmymT由式(4.43)最后一个等式可知
- 配套讲稿:
如PPT文件的首页显示word图标,表示该PPT已包含配套word讲稿。双击word图标可打开word文档。
- 特殊限制:
部分文档作品中含有的国旗、国徽等图片,仅作为作品整体效果示例展示,禁止商用。设计者仅对作品中独创性部分享有著作权。
- 关 键 词:
- 第四 线性 系统 参数估计 最小二乘法