2动态过程数学模型参数估计的最小二乘方法

PAGE1PAGE8（2011教案）系统辨识第二章第二章参数估计的最小二乘方法LeastSquares§2—１静态线性模型参数的最小二乘估计（多元线性回归）什么是最小二乘估计系统辨识三要素：模型，数据，准则。例：y=ax+其中：y、x可测；—不可测的干扰项；a—未知参数。通过N次实验，得到测量数据yk和xkk=1、2、3…，确定未知参数a称“参数估计”。使准则J为最小：令：Ja=0,导出a=?称为“最小二乘估计”，即残差平方总和为最小的估计，Gauss于1792年提出。二、多元线性回归线性模型y=a0+a1x1++anxn+式(2引入参数向量:=[a0,a1,an]T(n+1)1进行N次试验，得出N个方程：yk=kT+k;k=1、2…、N式(2其中：k=[1,x1,x2,,xN]T(n+1)1方程组可用矩阵表示为y=+式（2-1-3）其中：y=[y1,y2,。。。,yN]T(N1)=[1,2,。。。,N]T(N1)N(n+1)估计准则有：=（y—）T(y—)(1N)(N1)J=yTy+TT-yT-TTy=yTy+TT-2TTy式（2-1-4）假设：（T）(n+1)(n+1)满秩，由利用线性代数的以下两个矩阵对向量求偏导数的公式:和(必须A为对称阵)由为对称阵有:和所以:令上式等于零，解出参数估计向量：Ls=（T）-1Ty式（2-1-5）令：P=（T）-1则参数估计向量Ls=PTy参数估计向量Ls被视为以下“正则方程”的解：（T）=Ty式（2-1-6）注：为了便于区别,我们用红体字符表示估计量或计算值，而用黑体表示为参数真值或实际测量值。又所以Ls使J为最小。三、关于参数最小二乘估计Ls性质的讨论以上求解参数最小二乘估计Ls时并为对{k}的统计特性做任何规定，这是最小二乘估计的优点。当{k}为平稳零均值白噪声时，则Ls有如下良好的估计性质：参数最小二乘估计Ls是y的线性估计Ls=PTy是y的线性表出；b)参数最小二乘估计Ls是无偏估计，即ELs=（参数真值）[证明]：ELs=E[PTy]=PTE(y)=PTE(+)=PT+E()=+0=最小二乘估计Ls的估计误差协方差阵是2P(n+1)(n+1)即：E[(Ls-)(Ls-)T]=2P[证明]：E[(Ls-)(Ls-)T]=E[PT(y-)(y-)TP]=E[PTTP]=PTE(T)P=

PT2INNP=PT2P=若{k}为正态分布零均值白噪声时，则Ls是线性无偏最小方差估计（证明从略）。如若{k}是有色噪声，则Ls不具有上述性质，即为有偏估计。四、最小二乘估计Ls的几何意义和计算问题1.最小二乘估计的几何意义

最小二乘估计的模型输出值为yk=kTLsk=1,2,…N输出实际测量值与模型输出值之差叫残差：k=yk–yk模型输出向量为y=Ls，而残差向量为：=y–y=y–LsT=Ty–T（T）-1Ty=Ty–Ty=0即残差向量与由测量数据矩阵的各个向量：1,2,…,N张成的超平面（估计空间）正交，而最小二乘模型输出向量y为实际输出向量y在估计空间上的正交投影，这就是最小二乘估计的几何意义。最小二乘估计的几何意义2.关于最小二乘估计计算中的病态问题估计参数向量Ls一般是求解正则方程：（T）=Ty式（2-1-6）得出。可以利用消元法等一系列求解多元线性一次方程组的方法，计算得出，其有解的条件是（T）=P–1矩阵非奇异（行列式数值大于零）。但有时在求解式（2-1-6）方程组是会出现矩阵接近于奇异（行列式数值接近于零），即所谓“病态”的情况。由此导致参数估计的结果不稳定，不可信。出现上述情况的原因可能是由于①被辨识的过程受到的外加激励不够；②采样间隔太密；③A/D转换的位数太短，计算舍入误差累计所致。为解决最小二乘计算中可能出现的病态问题，提出了不少改进算法，例如：Householder变换法、改进的平方根法和U—D分解算法。后者是Bierman1977提出的改善P阵计算性质（对称性、正定性和稳定性）而又不增加计算量的算法。正定P阵可以分解成一个上三角阵U（其对角线元素都为1）和一个对角阵DP=UDUT由此可解决最小二乘计算中可能出现的病态问题，具体可参阅关于《计算方法》的文献。总之，我们在使用最小二乘的辨识方法时，应该注意避免出现和克服病态问题。应用举例在建立生产过程的静态模型时，特别是在机理不清之时常用多元线性回归方法，例如：水泥凝固放热量与水泥成分的关系模型y=a0+a1x1+a2x2+a3x3+a4x4+y水泥凝固时的放热量（卡/克）;x1~x1水泥的几种成分。五、非线性最小二乘法（NonlinearLeastSquare）以“误差平方总和为最小”的估计准则，估计非线性模型参数的方法。假设非线性静态系统模型为y=f(x,)+非线性模型f的形式是已知的，参数未知。经过N次实验，取得N组数据（x1,y1）（xN,yN）。准则：需要用优化算法求解，常用的算法有两类——搜索法和迭代法。前者如单纯型法；后者如梯度法、高斯法、牛顿—拉夫森法、变尺度法等等。该类方法也还可应用于动态模型和时间连续模型的参数估计。单纯型法是先给出参数空间的几个猜想点，构成正多面体，计算各点的目标函数值，比较各值后舍去最差的点，按照反射、开拓、收缩等步骤确定新的估计点，直到预定的精度要求后停止搜索。迭代法是先猜想一个估计的初值，确定一个向量为可接受的方向和步长，进行迭代计算k+1–k=μ•v…。具体内容可阅有关计算方法的参考书籍。近年来发展出一系列基于生物进化论的优化新算法，如遗传算法（GeneticAlgorithms）(基于改进遗传算法的系统辨识方法.北方工业大学学报.第10卷,第1期,P62-67.)和免疫算法(ImmuneAlgorithms)，使得优化算法的性能得到了很大的改善。§2—2动态过程参数估计的线性最小二乘法一、模型类：考虑CAR模型式（2-2-1）其中{y(k)}和{u(k)}为可测的输出和输入，{(k)}为不可测的随机干扰。上式还可表示成：A(z-1)y(k)=B(z-1)u(k)+(k)式（2-2-2）其中：A(z-1)=1+a1z-1+a2z-2+。。。。+anz-nB(z-1)=b1z-1+b2z-2++bnz-n还可表示为

式（2-2-3）其中：当进行了k=1-n,2-n,..,0,1,2,…,N共计（N+n）次采样，得到N个方程：y(1)=-a1y(0)-…-any(1-n)+b1u(0)+…+bnu(1-n)+(1)y(2)=-a1y(1)-…-any(2-n)+b1u(1)+…+bnu(2-n)+(2)………………y(N)=-a1y(N-1)-…-any(N-n)+b1u(N-1)+…+bnu(N-n)+(N)用矩阵表示成yN=N+N式（2-2-4）其中：yN=[y(1),y(2),…,y(N)]TN=[(1),(2),…,(N)]T二、参数最小二乘估计Ls的导出估计准则为式（2-2-5）由解出：

N=（NTN）-1NTyN式（2-2-6）上式可视为以下正则方程的解（NTN）N=NTyN式（2-2-7）称为最小二乘的“一次完成算法”，是离线算法，有解的条件是（NTN）2n2n满秩。用消元法或平方根法解线性方程组，得出N。三、用配方法导出Ls（略）=(yN-N)T(yN-N)+yNTN（NTN）-1NTyN-yNTN（NTN）-1NTyN={-（NTN）-1NTyN}TNTN{-（NTN）-1NTyN}+yNTyN-yNTN（NTN）-1NTyN上式的后两项中均不含，能使得J=nim的条件是：=（NTN）-1NTyN即前式（2-2-6）用两种方法推证出相同结论。四、线性动态参数最小二乘估计Ls的性质在静态模型:yk=kT+k式(2中的k=[1,x1,,xN]T为确定性量，取值与yk统计性质无关；在动态模型:y(k)=kT+(k)式(2-2-的最小二乘估计Ls虽然形式上与静态的相同，但是式的k中包含y(k-1)、y(k-2)、…，导致有关估计的统计性质的证明要困难得多，不能简单地套用静态模型多元回归的结果。动态参数最小二乘估计Ls的估计性质的主要结果是：当N时:EN=0(渐进无偏估计)

a.s.

N(强一致性收敛)如若{(k)}为有色噪声，Ls是有偏估计。五、数字仿真y(k)-1.5y(k-1)+0.7y(k-2)=u(k-1)+0.5u(k-2)+(k)一般取N=100~200。仿真结果如下表：参数a1a2b1b2真值-1.50.710.5估计-1.50910.70741.00220.4656练习：P60例3.3，P62例3.4§2-3递推最小二乘方法(RLS法)一、递推算法的导出对式（2-2-1）的CAR过程{y(k)}和{u(k)}进行了k=1-n、…、N,共(n+N)次观测，组成了yN和N，可得出参数估计N，现在再进行一次新的采样，又得出N+1和新的估计N+1。如何由估计向量N经过递推直接得到由新的估计向量N+1，而不必反复做一次完成LS法的计算？先分析yN+1与yN以及由估计向量N+1与N式（2-3-1）式（2-3-2）由式（2-2-6）有由估计向量N+1=（N+1TN+1）-1N+1TyN+1=PN+1N+1TyN+1其中=(NTN+N+1N+1T)-1=(PN-1+N+1N+1T)-1根据以下矩阵求逆公式：若A是nn满秩矩阵，B和C是nm阵，且(A+BCT)nn和(I+CTA-1B)mm都满秩（I为m维单位阵），则有以下矩阵恒等式成立：(A+BCT)–1=A–1-A–1B(I+CTA-1B)–1CTA-1现令：P–1=A和N+1=B=Cm=1,考虑到前面得到的PN+1=(PN+N+1N+1T）-1得出式（2-3-3）其中：即相当于矩阵恒等式中的(I+CTA-1B)–1将式（2-3-3）、式（2-3-2）和式（2-3-1）代入由估计向量N+1=PN+1N+1TyN+1并令：式（2-3-4）将上式打开:N+1=PNNTyN+PNN+1y(N+1)–-KN+1N+1TPNNTyN–KN+1N+1TPNN+1y(N+1)N+1=N+(PNN+1-KN+1N+1TPNN+1)y(N+1)–KN+1N+1TN可以证明：(PNN+1-KN+1N+1TPNN+1)=KN+1式（2-3-4）N+1=N+KN+1(y(N+1)–N+1TN)式（2-3-5）二、递推算式的物理含义N+1=N+KN+1(y(N+1)–N+1TN)式（2-3-5）式（2-3-4）式（2-3-3）N+1T=[-y(N)，-y(N-n+1),u(N)，,u(N-n+1)]式（2-3-6（巧妙地设计出Ф使CAR能用LS法）模型：y(N+1)=N+1T+(N+1)N时刻对N+1时刻的预报y(N+1N)=N+1TN式（2-3-7（估计值N代入N+1T）预报误差(被称为新息)，用绿色表示ε(N+1)=y(N+1)-y(N+1N)=y(N+1)-N+1TN式（2-3-8则式（2-3-5）可表达成N+1=N+KN+1ε(N+1)式（2-3-5）物理意义：新的参数估计N+1是对上次老的估计N进行修正而得出的，修正是利用在N对新的输出y(N+1)预报的预报误差乘以一个修正系数向量。KN+1是修正系数向量，它需要递推计算得出，在递推计算KN+1时要用到估计误差的协方差阵PN，而后者也是递推得出的。三、初值选择和计算框图递推计算需要初值0和P0。初值的选择有两种方法：其一是，可以用初始的3n组数据用一次完成算法解出：P0=（T）-1和0=P0Ty,但是实际上并不用上述办法，而是用以下更为简便的方法。其二是令：0=0和P0=I2n2n,其中=103~106（即假设为足够大的正数）。可以证明：在递推2n步后，估计结果与前面介绍的精确初值相接近。递推计算框图如下图所示。四、LS法和RLS法的讨论（1）LS和RLS法数学等价均由使准则J=[y(k)-kT]2=min得出不要求对{（k）}的统计特性有任何验前知识。如果{（k）}为零均值白噪声，则可得渐进无偏估计，即当n时，E=0,且。若{（k）}为有色噪声，是有偏估计，但是因该算法简单，所以应用广泛。均可推广到多输入多输出系统。LS法和RLS法的比较LS法是一次完成算法,适于离线辩识,要记忆全部测量数据,算法复杂；RLS法是递推算法，适于在线辩识和时变过程，只需要记忆n+1步数据，程序简单；RLS法用粗糙初值时，如若N较小时，估计精度不如LS法。开始置初值0N由式（2-3-6）构成N+1进行第N+1次采样由式（2-3-4）NN+1计算KN+1由式（2-3-5）计算N+1由式（2-3-3）计算PN+1否终点判断是输出模型§2–4时变过程的参数估计（时变递推最小二乘法）参数的递推最小二乘估计(RLS)与一次完成的最小二乘估计(LS)是数学等价的，它们都仅适用于估计时间定常过程的参数，而不适用于估计时变过程的参数。时变过程的特点是过程的参数可能随着时间变化而改变。因此，它的数学模型参数具有“时间性”，在利用动态过程的输入—输出数据来辨识模型参数时，“老”的数据往往只能反映“老”的过程参数；而改变后的“新”参数，要靠用新的和比较新的实验数据来估计。因此时变过程参数估计的特点是，不同时段的实验数据的作用是有区别的。时变过程的参数估计有多种不同的算法，本节讲授最常用的“带遗忘因子的递推最小二乘估计算法”，以下简称为“遗忘因子法”。一、遗忘因子算法的思路算法的主要思路是“厚今薄古”，即对新--老数据给予不同的对待，逐渐遗忘老数据的影响。具体做法是每采得一个新的y(N+1)时，将以前的所有数据乘以小于1的加权因子（0<<1）,乘以加权因子后的估计准则为：即对不同时刻的残差予以不同加权的平方和，老数据的作用按照指数衰减，被遗忘。二、遗忘因子递推最小二乘算法的递推算式利用矩阵求逆公式得出：令：=2，遗忘因子（0<<1）;取值范围（0.95~0.995）,值愈小，“遗忘”愈快。带遗忘因子的RLS法由以下几个递推算式组成：N+1=N+KN+1(y(N+1)–N+1TN)式（2-3-5）式（2-4-1）式（2-4-2）三、遗忘因子的作用对于遗忘因子的不同值，得到了不同的遗忘效果：值较小时的估计跟踪参数时变的能力强（优点），但是噪声干扰影响造成的估计波动亦大（缺点），可通过以下的计算机仿真说明值的影响。模拟如下2阶过程：y(k)+a1y(k-1)+a2y(k-2)=b1u(k-1)+b2u(k-2)+(k)其中：a1=-1.5，a2=0.7,b1=1,b2=0.5,而给出值分别为1；0.995；0.99；0.95时的a1和b1：四、算法的改进1)变算法模型的一步预报误差为(k)=y(k)-kTk-1将2(k)与事先给定的m1值相比较(m1>0)当2(k)>m1时选择较小的值（？）当2(k)≤m1时选择较大的值（？）一般选择m1=（4~6）Var{(k)}。2)变P算法将2(k)与事先给定的m2值相比较(m2>0)当2(k)>m2时将P阵重新置为：106I比较：以2000步的预报误差平方总和Q做比较Q的理想值为2000RLS法变算法=0.999/0.9m1=4变P算法m2=25=0.995=1=0.99=1288962702502212221373）附加R阵法只要在P阵的递推公式上附加一个小小的正定附加方阵R，即可大大提高的跟踪能力。即：为事先给定的正值§2–5最小二乘法的应用例1．热交换器在线辨识（R.Isermann“Automatica“）

辩识以蒸汽流量调节阀为输入，热水出口温度为输出的动态模型，t=0.3s、N=31、n=3,用RLS法。结论：对采样的数据进行适当的滤波处理，对辩识的效果大有改善。例2．直流它激电动机参数在线辨识“自动化学报”直流电机模型的形式是已知的，主要参数如：机电时间常数Tm、电磁时间常数Ta、磁通量和负载转距Mc是需要辩识的，其中：假设Mc是常值，其方向总是与电机旋转方向相反，即：方案一、利用u和的采样值辩识Tm、Ta、Ce和Mc方案二、利用i和的采样值辩识J和Mc辨识结果如下表所示N81020304050J(km.m2)0.003650.004650.004480.004500.004510.00454MC(N.m)0.43800.46440.