辩识最小二乘法

1、1最小二乘法(Least Squares LS)是一种经典有效的数据处理方法。它是1795年高斯在预测行星和彗星运动的轨道时提出并实际使用的。在系统辩识和参数估计领域中，最小二乘法是一种最基本的估计方法。它可用于动态系统，也可用于静态系统；可用于线性系统，也可以用于非线性系统；可用于离线估计，也可以用于在线估计。在随机的环境下利用最小二乘法时，并不要求知道观测数据的概率统计信息，而用这种方法所得到的估计结果，却有相当好的统计性质。第五章辩识的最小二乘法第五章辩识的最小二乘法2这时的问题是已知系统的输入和输出，求参数ai和bi的估计值。模型可以改写为：5.1 最小二乘估计最小二乘估计设时不变SI

2、SO动态系统的数学模型为：11() ( )() ( )( )A zy kB zu ke k ( ) ,( )u ky k10( )(1)()( )()( )nny ka y ka y knbu kbu kne k1111101()1()nnnnA za za zB zbb zb z 其中：3 ( ) (1) (1) (1) (1) (2) (2) (2) (1) () () ()y nyu nuy nyu nuy nNy Nu nNu N (1) (2) ()TYy ny ny nNYe 将上述模型改写为以下最小二乘格式：1201( )(1) () ( ) (1) () TTnnky ky k

3、n u k u ku kna aa b bb 令K=n+1n+N，共N次观测。记(1) (2) ()Tee ne ne nN( )( )( )Ty kke k4 (2) (1) (2) (2) (2) (1)( ) (1) ()(1) (1) (1)y ny nyuy ny nyu nnuy n Nuy n N1201(2) ( ) () ( ) (1) ()nneny N u n Nu Nen NaaenabbbYe 5可见，残差包括两个误差因素。一是参数估计误差带来的拟合误差；二是随机噪声带来的误差。对于上述模型的辩识问题，其中都是可以观测的数据，是待估计的参数。引入最小二乘准则:( ),

4、 ( )y kk21 ( )Nk nJe k ( )( )( )( )( )( ) ( )()( )TTTTe ky kkke kkke k10( )( )(1)()( )()nne ky ka y ka y k nbu kbu k n( )e k其中称为残差或方程误差， 是参数估计值。进一步可以得到：( )e k10 Tnnaa bb6若非奇异，可以得到：来确定估计值。求J对的偏导数并令其等于0，可得：可以看出，指数函数即残差的平方和。最小二乘估计是在残差二乘准则函数极小意义下的最优估计，即按照准则函数：21 ( )Nk nJe k 1LSTTY 0TTTJYYYY minTTJe eYYT

5、TY T 称 为最小二乘估计值，对应方法称为最小二乘法LS7其中其中s st t表示钢产量，表示钢产量，t t表示年代。试用最小二乘问表示年代。试用最小二乘问题的一次完成算法确定参数题的一次完成算法确定参数 。并以此来。并以此来预测该企业预测该企业19501950年的钢产量。年的钢产量。若钢产量用如下模型描述：若钢产量用如下模型描述：例例: : 已知某企业已知某企业19461946年年-1949-1949年钢产量如表年钢产量如表年份年份 1946 1946 1947 1947 1948 1948 1949 1949产量产量 10.0 10.0 12.2 12.2 13.5 13.5 16.1

6、16.12012tstt210,8在推导最小二乘法的结果时，并没有考虑噪声e(k)的统计特性。但在评价最小二乘估计的性质时，则必须假设噪声e(k)是不相关的，而且是同分布的随机变量，也即假设e(k)是白噪声序列，即20, NNeNNE eCov eI2 22 (1)(1) (2) (1) ()(1) (2) (2) (2) () (1) () (2) () ()TNNNCov eE e eE eE eeE ee NE eeE eE ee NE ee NE ee NE eN9其中w(k)称为加权因子，对所有的k，w(k)都必须是正数。引进加权因子是为了考虑观测数据的可信度。如果有理由认为现时刻的

7、数据比过去时刻的数据可靠，那么现时刻的加权值就要大于过去时刻的加权值。比如，可选当 ，当 ，这就体现了对不同时刻的数据给予了不同程度的信任。一般来说，w(k)的选择多少取决于人的主观因素，并无规律可循。如果准则函数取为加权函数，即为21( )( )NTknJw k eke We ( ), 01Nkw k11, (1)1Nkw, ( )1kN w k10通过极小化 计算 的方法称为加权最小二乘法，对应的 称为加权最小二乘估计值。加权最小二乘估计的解为：TJe WeWLSLSWLS1TTWLSWWY WLS其中W是一对称正定阵若取W=I,则。所以，最小二乘法是加权最小二乘法的一种特例。当获得一批数

8、据之后，利用最小二乘法或者加权最小二乘法可一次求得相应的参数估计值，这样处理的方法称为一次完成算法或者批处理算法。这在理论研究方面有很多方便之处，但在计算方面要碰到矩阵求逆的困难。但当维数增加时，矩阵求逆的运算量将急剧增加，会给计算速度和存储带来负担。11可以用高斯消去法进行求解方程式，以便更快地求得参数的估计值。但是，更实用的方法还是设法化为递推计算的形式，以便在线辩识，大大减少数据的存储。高斯“未知量的最适合值(最可能值)是使各项实际观测值和计算值之间差的平方乘以度量其精确度的数值以后的和为最小”。12(3) 如果出现列相关，即不满秩的情况，为病态矩阵，则不能得到最小二乘估计值。(2) 每

9、增加一次观测量，必须重新计算一次 上一节给出了最小二乘一次完成算法，但具体使用时不仅占用内存量较大，而且不能用于在线辩识。进一步，一次完成算法还有如下缺陷：5.3 最小二乘递推算法最小二乘递推算法(RLS)(1) 数据量越多，系统参数估计的精度就越高。为了获得满意的辩识结果，矩阵的阶数常常取得非常大。这样，矩阵求逆的计算量很大，存储量也很大。1,T T T13系统用线性差分方程来描述： 解决这个问题的办法是把它化成递推算法。依观测次序的递推算法就是每获得一次新的观测数据就修正一次参数的估计值，随着时间的推移，便能够获得满意的辩识结果。递推辩识算法具有无矩阵求逆，以及跟踪时变系统等特点，这样不仅

10、可以减少计算量和存储量，而且能实现在线辩识。( )( )( )( )(1) () ( ) ()TTy kke kky ky knu ku kn 14上式的最小二乘解为：1TTNNNNNNTNNPYY (1)(1)(1)(2)(2)(2), , ()()()TTNNNTny ne ny ne nnYey nNe nNnN 10(1) ( ) (1) (1) (1)(2) (1) (2) (2) (2) ()(1) ( ) () ( )nnay ny nyu nuay ny nyu nuby nNy nNy N u nNu Nb(1)(2) ()e ne ne nNNNNYe 令：可得：151TT

11、NNNNNNNTNPYY NNNYe 111TNNNYe 1111 NNNNTNNYYy111111111TTTNNNNNNNNYYP如果再增加一组新的观测值u(n+N+1),y(n+N+1)，记作uN+1,yN+1，则又增加一个方程：消除矩阵求逆过程，用PN来表示PN+1。1611111111111111111 NTTNNNNNTNTTTNNNNNNNNTNNNNIPPPIPP 17令，则得到：111111ABCAA B ICA BCA1111,TM MMNNMNAI BPC矩阵求逆引理：设A为nn矩阵， B为nm矩阵，C为mn矩阵，并且A，A+BC和I+CA-1B都是非奇异阵，则有恒等式：

12、1111111111111 1 TNNNNNTTNNNNNNNNTNNNNNTNNNPIPPPPIPPPPPP18111111111111111111111111111 111TTNTTNNNNNNNNTNNNNNTTNNNNNNNNNNNNTTNNNNNNTNNNNNNNNNNNTNNNNPYYPYPyPYyPPYPPyPyPPPyP式中为增益矩阵，记为GN+1，而为预报误差。11TNNNy1111NNTNNNPP19递推过程如下：1111TNNNNNNGy综上所述，得到最小二乘估计递推算法如下：111TNNNNNPPGP00111222,PGPGP11111NNNTNNNPGP20方法1:

13、m组数据，用LS一次算法，得到，再从m+1开始递推11( )( )max( )NNiNiii对于初值的选取：00,P1( )Ni20PI,mmP11TTmmmm mY 1TmmmP 方法2: ，任取00,另外，可以用下式作为递推算法的停机准则：式中为参数向量的第i个元素在N+1次递推计算结果， 为给定的表示精度要求的某一正数。21所谓数据饱和，是指随着时间的推移，采集的数据越来越多，新数据提供的信息被旧数据淹没。如果辩识算法对新、旧数据给予相同的信度，那么随着从新数据中获得的信息量相对下降，算法就会慢慢失去修正能力。在实际应用中，这时参数估计值可能偏离真值较远而无法更新。(1) 数据饱和现象

14、随着数据量的增长，递推的最小二乘法将出现所谓的“数据饱和”现象。这是由于增益矩阵随着数据的增加将逐渐趋于零，以致递推算法失去修正能力的缘故。下面针对数据的饱和现象，讨论渐消记忆法、限定记忆法和振荡记忆法等适应性算法。5.4 数据递推的饱和及解决办法数据递推的饱和及解决办法22由此可知，PN是递减的正定阵。随着递推次数的增加，这会导致PN0。所以增益矩阵GN+1也随着N的增加而逐渐趋于零向量，从而使得算法失去修正能力。因为：1111101TNNNNNNTNNNPPPPP1NNPP另外，由于递推在有穷字长上的计算机上实现时，每步都存在舍入误差。因此数据饱和后，由于这些原因致使新的采样值不仅对参数估

15、计不起改进作用，反而使得所计算的PN失去正定性，甚至失去对称性，造成参数的估计量与真实参数之间的偏差越来越大。23为了克服数据饱和现象，可以用降低旧数据影响的办法来修改算法。对于时变系统，参数随着时间变化，在辩识算法中必须充分利用新数据所包含的信息，尽可能降低旧数据的影响，以便获得跟踪参数变化的实时估计。24加衰减因子后的数据阵为：渐消记忆法又称为遗忘因子法，这种方法的思想是对旧数据加上遗忘因子，按指数加权来使得旧数据的作用衰减。(2) 渐消记忆法(RFF)最小二乘估计值为：(01)1111, NNNNTNNYYy1TTNNNNNY 指数函数为：21( )NTNNNkJe ke e新增观测yN

16、+1之后的数据阵为：1111, NNNNTNNYYy2511211NNNTNNNPGP得到如下递推公式：11122111TNNNNNTNNNPPIPP211111111TNNNNNNNNNNNPPyGy112111111TTNNNNNNPP 261111NNNTNNNPGP令 为遗忘因子。综上分析，渐消记忆法的递推算法(RFF)可以归纳为： 1111TNNNNPIGP1111TNNNNNNGy2,01RFF算法的结构和计算流程与RLS算法基本是一致的。初始状态的选取也可以用上一节的方法。但是遗忘因子必须选择接近于1的正数，通常不小于0.9。27限定记忆法每次估值只依据最新的N个数据，在此之前的

17、数据则全部删除。如考虑一个固定长度的矩形窗，每一个时刻有一个新数据点增加进来，一个旧数据点剔除出去，这样保持了每次都只取最新的N个数据进行计算。具体算法分为两个部分：(3) 限定记忆法,1,1,Ti N ii Nii N ii Ni Ni NiGy1.先进一个观测数据yi+N，即在i+N时刻，进一个数据递推公式为：,1,1,1iNiN iiNiTiNiNiiNGPP,1,Ti N ii N ii Ni NiPIGP28,1,1,TiN iiN iiN iiiiN iGy 2.再出一个数据yi，即剔除一个旧的观测数据递推公式为：,1,1iiN iiN iTiiN iiGPP,1,1,Ti N i

18、i N iii N iPIGP(4) 振荡记忆法振荡记忆法采用整段删除N组数据的方法，即当数据长度以达2N时，可以删除开始的N个数据，使得数据在N到2N-1个之间变化。29 我们在讨论最小二乘估计的统计性质时，发现当系统的噪声满足白噪声性质时，参数估计是无偏一致最小方差估计。但在一般情况下，系统的噪声都不是白噪声。为了获得好的估计效果，我们考虑广义最小二乘法、扩展最小二乘法、辅助变量法和相关二步法等参数辩识方法，每种方法都对应着一种特定的噪声模型结构。 广义最小二乘法(GLS)的基本思想在于对数据进行一次白化滤波处理，然后利用基本的最小二乘法对滤波后的数据进行辩识。如滤波模型选得合适，为数据进

19、行了较好的白化处理，则利用基本的最小二乘法就能得到无偏一致估计。5.5 广义最小二乘法广义最小二乘法30这类问题的辩识可用广义最小二乘法，以便获得无偏一致估计。令设SISO动态系统的数学模型为：11() ( )() ( )( )A zy kB zu ke k11( )() ( )( )() ( )ffy kC zy ku kC zu k( )k1() ( )( )C ze kk其中e(k)可以表示为一个以白噪声序列 为输入的线性系统的输出，即它满足自回归模型：111()1mmC zc zc z 分别为白化处理后的输出和白化后的输入。31上式与基本最小二乘法的模型是一样的，若已知这样可以估计出；

20、若 未知，通常可以用松弛算法来估计参数。则数学模型可以表示为：11()( )()( )( )ffA zykB zukk1()C z ,iia b1()C z11( )() ( )( )() ( )ffy kC zy ku kC zu k323. 计算C(z-1)，即通过得到参数估计值1. 先猜一个C(z-1)的值，即设C(z-1)=1，利用基本的最小二乘法对方程中的A(z-1)和B(z-1)进行估计，得到。(1) 广义最小二乘法(GLS)的松弛算法11( )() ( )() ( )e kA zy kB zu k11()( )()( )( )ffA zy kB zu kk1() ( )( )C

21、ze kk11(),()A zB z11(), ()A zB z1,mcc2. 进一步计算e(k)，即4. 利用C(z-1)可以计算yf(k)和uf(k)，通过等式计算新的估计值5. 重复步骤(2)，直到估计的精度达到要求为止。33所谓精度要求，可以看作下列不等式是否满足：(1)ik(1)( )max(1)iiiikkk式中为参数向量的第i个元素在k+1次迭代的计算结果广义最小二乘法就是根据上述过程，反复估计噪声参数C(z-1)与系统模型参数A(z-1),B(z-1)的一种迭代算法广义最小二乘法的收敛速度是比较慢的，需要经过多次迭代计算，才能得到较准确的参数估计值。一般情况下，估计值能够收敛到

22、稳态值。34(2) 递推算法(RGLS)假设前一时刻算出的系统模型参数估计为，噪声模型参数为 ，相应的算子多项式为，现时刻采集数据为yk+1和uk+1，则一次递推过程的步骤如下：(1),()Tkkkcccmk()()()()111 fffffkkk nkk nyyuu 1. 对新数据进行滤波并构成滤波后的输入输出数据向量1()kCz( )111111()()fkkkfkkkyC zyuC zu35()()()11()()()111fffkkkf TffkkkPGP111 Tkkk nkk nyyuu 1k2.对滤波后的输入输出数据作RLS估计，修正系统模型参数，计算公式为：()()()()111fff TfkkkkPIGP3.由此得到新的计算出新的残差估计值1111Tkkkkey()()()1111fffkkkkkkGy( )1-1-1 Tekkkk meee 36每一次递推过程都包括两次RLS的计算，通过滤波计算和残差计算将它们联系起来。( )( )( )11( )( )( )111eeekkke TeekkkPGP4.对新残差数据作RLS估计，修正系统模型参数，计算公式为：( )( )( )( )111eee TekkkkPIGP( )( )1111ee TkkkkkkccGe