1、第五讲第五讲 处理有色噪声扰动的处理有色噪声扰动的LS类方法类方法(1/3)第五讲第五讲 有色噪声扰动的最小二乘类方法有色噪声扰动的最小二乘类方法q 由上一讲由上一讲中所讨论的中所讨论的LS估计的统计特性可知估计的统计特性可知,在动态系统的在动态系统的参数估计中参数估计中,LS法仅当随机扰动为零均值的白噪声时法仅当随机扰动为零均值的白噪声时,才能保才能保证得到无偏一致的估计值证得到无偏一致的估计值.在实际工程和社会系统的辨识中在实际工程和社会系统的辨识中,随机扰动随机扰动w(k)只是各种只是各种系统内外扰动和结构建模误差等因素的综合反映系统内外扰动和结构建模误差等因素的综合反映.因此因此,并不
2、能保证随机扰动为统计独立的白噪声并不能保证随机扰动为统计独立的白噪声.故故,实际系统辨识时实际系统辨识时,并不能保证并不能保证LS估计值为无偏一致估估计值为无偏一致估计值计值.第五讲第五讲 处理有色噪声扰动的处理有色噪声扰动的LS类方法类方法(2/3)q 从本讲开始从本讲开始,我们将讨论具有随机扰动为我们将讨论具有随机扰动为相关性扰动相关性扰动(有色噪有色噪声声)的动态系统的无偏一致的参数估计问题的动态系统的无偏一致的参数估计问题,具体的具体的LS类方法类方法有有 增广最小二乘法增广最小二乘法(Extended Least-square Method,ELS)广义最小二乘法广义最小二乘法(Ge
3、neralized Least-square Method,GLS)辅助变量法辅助变量法(Instrumental Variables Method,IV)等等.其它无偏参数估计算法还有其它无偏参数估计算法还有 多步最小二乘法多步最小二乘法 偏差补偿最小二乘法偏差补偿最小二乘法第五讲第五讲 处理有色噪声扰动的处理有色噪声扰动的LS类方法类方法(3/3)q 这些不同的处理有色噪声的辨识方法主要是针对不同的这些不同的处理有色噪声的辨识方法主要是针对不同的 有色噪声的特性、有色噪声的特性、有色噪声的不同模型表达、以及有色噪声的不同模型表达、以及 不同的辨识要求不同的辨识要求提出的提出的.下面就分别讨
4、论下面就分别讨论:增广最小二乘法增广最小二乘法(ELS)广义最小二乘法广义最小二乘法(GLS)辅助变量法辅助变量法(IV)其它方法可参阅其它教材和文献其它方法可参阅其它教材和文献.1.ELS法法(1/19)1.增广最小二乘法增广最小二乘法q 考虑如下考虑如下SISO随机离散系统随机离散系统A(z-1)y(k)=B(z-1)u(k)+v(k)(1)其中其中y(k),u(k)和和v(k)分别为系统的输出、输入和相关的随机扰动分别为系统的输出、输入和相关的随机扰动;对所考虑的模型的相关随机扰动对所考虑的模型的相关随机扰动v(k),一般假定其一般假定其为平稳相关序列为平稳相关序列.baniiiniii
5、zbzzaz1111)(1)(BA1.ELS法法(2/19)由由第二讲第二讲中的关于有色噪声的结论中的关于有色噪声的结论1和假设和假设2可可知知,平稳的相关扰动平稳的相关扰动v(k)可被建模如下可被建模如下v(k)=C(z-1)w(k)(2)其中其中w(k)为白噪声序列为白噪声序列;C(z-1)为为未知的未知的、稳定的稳定的、有限阶的有限阶的线性滤线性滤波器波器,并可表示为如下首一多项式并可表示为如下首一多项式cniiizcz111)(C 因此因此,由系统模型由系统模型(1)和噪声模型和噪声模型(2),可得如下表示可得如下表示A(z-1)y(k)=B(z-1)u(k)+C(z-1)w(k)1.
6、ELS法法(3/19)当假定噪声当假定噪声w(k)可测已知时可测已知时,上式又可表示为如上式又可表示为如下自回归方程下自回归方程y(k)=t t(k-1)q)qw(k)(3)其中其中111.,.,-.,-)-(.,),1-()-(.,),1-()-(.,),1-()1-(cbannncbaccbbaankwkwnkukunkykykA(z-1)y(k)=B(z-1)u(k)+C(z-1)w(k)1.ELS法法(4/19)q当当上述上述观测数据向量观测数据向量(k-1)精确已知时精确已知时,利用前面讨利用前面讨论的论的成批或成批或RLS法可求得回归参数向量法可求得回归参数向量q q的的LS估计估
7、计值值.但是但是,实际上实际上上述上述数据向量数据向量(k-1)中包含有不可测中包含有不可测的噪声量的噪声量w(k-1),.,w(k-nc),因此对自回归式因此对自回归式(3)并并不能直接套用不能直接套用LS估计方法估计方法.为此为此,引入通过在递推参数估计过程中引入通过在递推参数估计过程中在线估计在线估计噪声噪声w(k)以实现模型参数在线递推估计的以实现模型参数在线递推估计的ELS法法.1.ELS法法(5/19)ELS法的思想就是法的思想就是:利用该参数估计值来在线估计白噪声利用该参数估计值来在线估计白噪声w(k-i)的值的值 (k-i)以替代数据向量以替代数据向量(k-1)中的白噪声中的白
8、噪声w(k-i),w 然后进行下一步的参数估计然后进行下一步的参数估计.在递推估计过程中在递推估计过程中,假设当前或前一步的在线参数估假设当前或前一步的在线参数估计值已相当程度可用的前提下计值已相当程度可用的前提下,选择噪声w(0)的估计值 对输入u(k)和输出y(k)采样 构成观测数据向量)1(k 下一个采样周期 根据辨识模型估计噪声w(k)辨识A(z-1)、B(z-1)和C(z-1)1.ELS法法(6/19)q噪声噪声w(k)的具体的估计算法是如下的事后估计或事的具体的估计算法是如下的事后估计或事前估计算法前估计算法:其中其中(k-1)为数据向量为数据向量(k-1)的如下在线估计值的如下在
9、线估计值)4()1-()1-(-)()()1-(-)()(事前估计事后估计kkkykkkykw)-(.,),1-()-(.,),1-()-(.,),1-()1-(cbankwkwnkukunkykyk1.ELS法法(7/19)q 因此因此,基于渐消记忆的基于渐消记忆的RLS估计算法可得如下递推的估计算法可得如下递推的ELS法法11)1()()1(kkkwkwkwk事后估计其中加权因子其中加权因子l l为为1时就为普通的时就为普通的ELS法法.上述上述ELS辨识实际上是辨识实际上是q q的的参数估计和噪声参数估计和噪声wk的估计交的估计交替进行替进行,计算顺序为计算顺序为:)7()1-(2)-(
10、)1-()1-(2)-()1-()6(-1)2),-(k)1-(1)-(-11)-()5()1-()1-(-)()1-()1-()(2kkkkkkkkkkkkykkkPPKIPPKIPK)1()()1(11kwkwkwkkk事前估计1.ELS法法(8/19)q上述分析过程表明上述分析过程表明,ELS法是法是LS法的一种简单推广法的一种简单推广.它只是扩充了参数向量它只是扩充了参数向量q q和数据向量和数据向量(k)的维数的维数,也同时辨识噪声模型也同时辨识噪声模型.就这种意义上说就这种意义上说,可称之为可称之为ELS法法.值得指出的是值得指出的是,ELS法虽然可同时得到噪声模型法虽然可同时得到
11、噪声模型的参数估计的参数估计,但其收敛过程却比过程模型但其收敛过程却比过程模型A(z-1)和和B(z-1)的估计值的收敛慢许多的估计值的收敛慢许多.从实用角度来说从实用角度来说,噪声模型噪声模型C(z-1)的阶次不宜的阶次不宜取得太高取得太高.1.ELS法法(9/19)-ELS法计算步骤法计算步骤q综上所述综上所述,ELS法的基本计算步骤可总结如下法的基本计算步骤可总结如下 确定被辨识系统模型的结构确定被辨识系统模型的结构,以及多项式以及多项式A(z-1)、B(z-1)和和C(z-1)的阶次的阶次;设定递推参数初值设定递推参数初值q q(0),P(-1),w(0);采样获取新的观测数据采样获取
12、新的观测数据y(k)和和u(k),并组成观测数据并组成观测数据向量向量(k-1);用式用式(5)(7)所示的所示的ELS法计算当前参数递推估计值法计算当前参数递推估计值;用用(4)式计算白噪声式计算白噪声w(k)的事后或事前在线估计值的事后或事前在线估计值w(k);1.循环次数循环次数k加加1,然后转回到第然后转回到第3步继续循环步继续循环.1.ELS法法(10/19)-ELS法仿真程序法仿真程序q 下面给出针对随机线性离散系统下面给出针对随机线性离散系统(XARMA)()1()()()1()()1()(111cnbnank-nwck-wckwk-nubk-ubk-nya.k-yakycba的
13、的ELS在线辨识仿真程序伪代码在线辨识仿真程序伪代码./%第一步第一步 初始化初始化输入系统阶次输入系统阶次na,nb和和nc,以及加权因子以及加权因子l l输入系统模型输入系统模型Az=1 a1 a2 和和Bz=0 b1 b2;输入噪声模型输入噪声模型Cz=1 c1 c2 输入系统输入信号输入系统输入信号u(k)的方差的方差 u、过程噪声、过程噪声w(k)的方差的方差 w和输入输出测和输入输出测量噪声量噪声 uw、yw设定系统变量初始值设定系统变量初始值:yf1:na+1=0;uf1:nb+1=0;wf1:nc+1=0;设定辨识变量初始值设定辨识变量初始值:yb1:na+1=0;ub1:nb
14、+1=0;wb=1:nc+1=0;q q1:na+nb+nc=0;P=106*I(na+nb+nc,na+nb+nc);1.ELS法法(11/19)-ELS法仿真程序法仿真程序)/%第二步第二步 辨识仿真辨识仿真for k=1:最大仿真步数最大仿真步数/%被控对象模型仿真被控对象模型仿真(产生系统输入输出信号产生系统输入输出信号,即数据即数据)yf2:na+1=yf1:na;uf2:nb+1=uf1:nb;wf2:nc+1=wf1:nc;uf1=2*u*(rand()-0.5);wf1=2*w*(rand()-0.5);yf1=-Az2:na+1*yf2:na+1+Bz2:nb+1*uf2:n
15、b+1+Cz1:nc+1*wf1:nc+1;1.ELS法法(12/19)-ELS法仿真程序法仿真程序/%输入输出数据检测输入输出数据检测ub1=uf1;yb1=yf1;/%或模拟检测噪声或模拟检测噪声%ub1=uf1+2*uw*(rand()-0.5);%yb1=yf1+2*yw*(rand()-0.5);/%在线递推辨识过程仿真在线递推辨识过程仿真=-yb(2:na+1)ub(2:nb+1)wb(2:nc+1);K=P*/(l l+*P*);q q=q q+K*yb(1)-*q q;P=I-K*P/l l;修正矩阵修正矩阵P；输出在线递推参数估计值输出在线递推参数估计值q q；1.ELS法法
16、(13/19)-ELS法仿真程序法仿真程序wb1=yb(1)-*q q;yb2:na+1=yb1:na;ub2:nb+1=ub1:nb;wb2:nc+1=wb1:nc;也可采用事前估计也可采用事前估计1.ELS法法(14/19)例例1q 例例1 考虑如图下所示的仿真对考虑如图下所示的仿真对象象,图中图中 1-1.5z-1+0.7z-2 1-1.0z-1+0.2z-2 1-1.5z-1+0.7z-2 1.0z-1+0.5z-2 ww(k)u(k)y(k)+图1 仿真对象 通过控制通过控制 w值来改变数据的噪信比值来改变数据的噪信比.q 辨识中辨识中,选择如下模型结构选择如下模型结构y(k)+a1
17、y(k-1)+a2y(k-2)=b1u(k-1)+b2u(k-2)+w(k)+c1w(k-1)+c2w(k-2)(8)其它条件同第四讲中的例其它条件同第四讲中的例1.结果如下表所示结果如下表所示.w(k)是服从均值为零是服从均值为零,方差方差为为1的正态分布的不相关的正态分布的不相关随机噪声随机噪声;输入信号输入信号u(k)采用伪随机采用伪随机二进制序列二进制序列;1.ELS法法(15/19)表表1 计算机仿真结果计算机仿真结果(噪信比噪信比=23%,数据组数数据组数1000)参数参数a1a2b1b2c1c2真值真值-1.50.71.00.5-1.00.2估计值估计值l l=1-1.4960.
18、71170.99920.4982-0.90980.1193估计值估计值l l=0.98-1.4650.69401.06600.5506-0.97780.33691.ELS法法(16/19)q 递推辨识过程的辨识值如下图所示遗忘因子l=1时递推辨识结果(1)参数估计误差的平方和噪声估计误差1.ELS法法(17/19)遗忘因子l=1时递推辨识结果(2)参数估计误差的平方和噪声估计误差1.ELS法法(18/19)遗忘因子l=0.98时递推辨识结果(1)参数估计误差的平方和噪声估计误差1.ELS法法(19/19)遗忘因子l=0.98时递推辨识结果(2)参数估计误差的平方和噪声估计误差2.GLS2.GL
19、S法法(1/3)(1/3)2.广义最小二乘法广义最小二乘法q 上一讲上一讲讨论了相关随机扰动讨论了相关随机扰动v(k)可用可用C(z-1)w(k)建模情况下的建模情况下的参数估计问题参数估计问题.但有些相关扰动用但有些相关扰动用C(z-1)w(k)来建模的话来建模的话,线性滤波器线性滤波器C(z-1)的阶次相当高的阶次相当高,这加大了参数估计的工作量这加大了参数估计的工作量,也极大地影也极大地影响了建模的精度和使用上的困难性响了建模的精度和使用上的困难性.在在此此情况情况下下,相关扰动相关扰动v(k)可用如下方式建模可用如下方式建模v(k)=w(k)/D(z-1)(1)其中其中w(k)为零均值
20、的白噪声为零均值的白噪声;2.GLS2.GLS法法(2/3)(2/3)1/D(z-1)为为未知的未知的、稳定的稳定的、有限阶的有限阶的线性滤线性滤波器波器,且且D(z-1)可表示为如下稳定的可表示为如下稳定的、阶次较、阶次较低的低的首一多项式首一多项式dniiizdz111)(D 因此因此,具有该类相关随机扰动的随机离散系统的具有该类相关随机扰动的随机离散系统的数学模型为数学模型为:A(z-1)y(k)=B(z-1)u(k)+v(k)=B(z-1)u(k)+w(k)/D(z-1)(2)2.GLS2.GLS法法(3/3)(3/3)q 广义最小二乘广义最小二乘(Generalized Least-
21、squares,GLS)法讨论的法讨论的是可用模型是可用模型(2)表示的随机离散系统的参数估计问题表示的随机离散系统的参数估计问题.下面下面,将分别讨论将分别讨论 成批型和成批型和 递推型递推型GLS法法的思想和算法的思想和算法,以及以及 GLS仿真仿真 GLS评述评述一、成批型一、成批型GLS法法(1/9)(1/9)一、成批型一、成批型GLS法法q 成批型成批型GLS法的基本思想是法的基本思想是:先预选定一个线性滤波器先预选定一个线性滤波器Df(z-1)将模型中的输出将模型中的输出y(k)、输、输入入u(k)和有色噪声和有色噪声w(k)/D(z-1)白色化白色化,即将模型即将模型A(z-1)
22、y(k)=B(z-1)u(k)+w(k)/D(z-1)两边左乘两边左乘线性滤波器线性滤波器Df(z-1)后后,有有A(z-1)Df(z-1)y(k)B(z-1)Df(z-1)u(k)+w(k)并记为并记为A(z-1)yf(k)B(z-1)uf(k)+w(k)再对再对上述上述白色化后的模型中多项式白色化后的模型中多项式A(z-1)和和B(z-1)的系数的系数进行进行LS估计估计;一、成批型一、成批型GLS法法(2/9)(2/9)然后基于对然后基于对A(z-1)和和B(z-1)的估计的估计,利用模型利用模型A(z-1)y(k)=B(z-1)u(k)+v(k)来计算模型残差来计算模型残差(相关噪声相
23、关噪声)v(k)的估计值的估计值 再基于噪声模型再基于噪声模型D(z-1)v(k)=w(k)和相关噪声和相关噪声v(k)的估计值的估计值,利用利用LS法法辨识辨识D(z-1),并基于并基于此此修正选定的滤波器修正选定的滤波器Df(z-1);如此循环往复地进行如此循环往复地进行Df(z-1)的修正及的修正及A(z-1)和和B(z-1)的迭代估计的迭代估计,直到模型的估计残差满足给定的精度直到模型的估计残差满足给定的精度或者或者Df(z-1)、A(z-1)和和B(z-1)的变化量相对较小为止的变化量相对较小为止.一、成批型一、成批型GLS法法(3/9)(3/9)计算步骤计算步骤q 上述上述成批型成
24、批型GLS法基本思想法基本思想可用如下流程图图示可用如下流程图图示.预选 白化 滤波器 输入输出数据 采样 输入输出数据滤波(白化)处理 辨识A(z-1)和B(z-1)计算模 型辨识 残差 辨识 D(z-1)下一次迭代 q 综上所述综上所述,成批型的成批型的GLS法的算法和计算步骤如下法的算法和计算步骤如下:Step 1.确定被辨识系统模型的结构确定被辨识系统模型的结构,以及多项式以及多项式A(z-1)、B(z-1)和和D(z-1)的阶次的阶次;Step 2.选定选定稳定的稳定的初始滤波器初始滤波器Df(z-1);一、成批型一、成批型GLS法法(4/9)(4/9)计算步骤计算步骤Step 3.
25、采样获取新的观测数据采样获取新的观测数据y(k)和和u(k);Step 4.基于滤波器基于滤波器Df(z-1),进行如下滤波计算进行如下滤波计算yf(k)=Df(z-1)y(k)(3)uf(k)=Df(z-1)u(k)(4)Step 5.由由(2)(4)式式,列成如下自回归方程列成如下自回归方程:)5()()1()(fffkwkkyt其中11ffffff,.,-,.,-)-(.,),1-()-(.,),1-()1-(bannbabbaankukunkykyk一、成批型一、成批型GLS法法(5/9)(5/9)计算步骤计算步骤 并由此得如下向量式的自回归方程组并由此得如下向量式的自回归方程组Yf=
26、fq qf+W (6)其中其中Yf=yf(1),yf(2),.,yf(L)t tW=w(1),w(2),.,w(L)t t f=f(0),f(1),.,f(L-1)t tStep 6.用如下用如下LS法计算法计算f参数估计值参数估计值)7()(ff1fffY一、成批型一、成批型GLS法法(6/9)(6/9)计算步骤计算步骤Step 7.由由(2)式有式有y(k)=t t(k-1)q qf+v(k)(8)其中其中(k-1)=y(k-1),y(k-na)u(k-1),u(k-nb)t t 因此因此,基于基于(7)式的参数估计值式的参数估计值 f,计算相关扰动计算相关扰动v(k)的估计的估计值值,即
27、回归式即回归式(8)的如下估计残差的如下估计残差)9(,.,2,1)1-(-)()(fLkkkykv一、成批型一、成批型GLS法法(7/9)(7/9)计算步骤计算步骤Step 8.由由(1)式式,可得如下关于相关扰动可得如下关于相关扰动v(k)和白噪声和白噪声w(k)的的自回归方程自回归方程)10()()1-()(vvkwkkv其中1fv-.,-)-(.,),1-()1-(dndddnkvkvk 并由此得如下向量式的自回归方程组并由此得如下向量式的自回归方程组V=vv+W (11)其中其中V=v(1),v(2),.,v(L)v=v(0),v(1),.,v(L-1)一、成批型一、成批型GLS法法
28、(8/9)(8/9)计算步骤计算步骤Vv1vvv)(Step 9.选取上一步估计得的选取上一步估计得的D(z-1)的估计值作进行的估计值作进行系统白色化处理的线性滤波器系统白色化处理的线性滤波器Df(z-1).若模型的估计残差满足给定的精度或者若模型的估计残差满足给定的精度或者Df(z-1)、A(z-1)和和B(z-1)的变化量相对较小则循环结束的变化量相对较小则循环结束,否则返回第否则返回第4步继续循环步继续循环.因此因此,将用式将用式(9)计算好的相关扰动计算好的相关扰动v(k)的估计值的估计值 (k)来替代自回归方程组来替代自回归方程组(11)的相关扰动的相关扰动v(k),则可则可得如下
29、回归参数向量得如下回归参数向量v,即噪声模型即噪声模型D(z-1)的系数的系数的的LS估计值估计值v 一、成批型一、成批型GLS法法(9/9)(9/9)q 以上辨识过程表明以上辨识过程表明,GLS法的思想是对输入输出数据先进行一法的思想是对输入输出数据先进行一次滤波预处理次滤波预处理,然后利用普通然后利用普通LS法对滤波后的数据进行辨识法对滤波后的数据进行辨识,并反复迭代此过程并反复迭代此过程.可想而知可想而知,这种方法受滤波模型的好坏的影响较大这种方法受滤波模型的好坏的影响较大,而在而在迭代过程中滤波模型的好坏也直接与系统模型辨识结果迭代过程中滤波模型的好坏也直接与系统模型辨识结果有直接的关
30、系有直接的关系.q 从优化理论的角度来说从优化理论的角度来说,GLS法其实属于非线性优化方法法其实属于非线性优化方法.因此因此,难以避免出现非线性优化中的局部极值点情况难以避免出现非线性优化中的局部极值点情况,即即该方法并不能保证得到的估计值是一致无偏的该方法并不能保证得到的估计值是一致无偏的.这是这是GLS法的一个不太令人满意之处法的一个不太令人满意之处.二、递推型二、递推型GLS法法(1/6)(1/6)二、递推二、递推GLS法法q 递推递推GLS法的基本思想是与成批型法的基本思想是与成批型LS法大致相同法大致相同,所不同的所不同的是是:由于所考虑的是递推估计由于所考虑的是递推估计,不能像成
31、批型那样作反复迭代不能像成批型那样作反复迭代.解决的方法是分别对解决的方法是分别对(5)和和(10)所示的两个模型的辨识设所示的两个模型的辨识设计两个递推估计算法计两个递推估计算法,并在每一个递推步中并在每一个递推步中,让它们依顺让它们依顺序递推一次序递推一次.随着递推过程的深入随着递推过程的深入,将不断改进噪声模型将不断改进噪声模型D(z-1)的辨识的辨识结果结果,同时亦得到较佳的同时亦得到较佳的A(z-1)和和B(z-1)的辨识结果的辨识结果.二、递推型二、递推型GLS法法(2/6)(2/6)q根据上述基本思想根据上述基本思想,有如下递推有如下递推GLS估计算法估计算法:Step 1.确定
32、被辨识系统模型的结构确定被辨识系统模型的结构,以及多项式以及多项式A(z-1)、B(z-1)和和D(z-1)的阶次的阶次;Step 2.选定递推参数初值选定递推参数初值 f(0)和和Pf(-1),v(0)和和Pv(-1),以及以及稳定的稳定的初始滤波器初始滤波器Df(z-1);Step 3.采样获取新的观测数据采样获取新的观测数据y(k)和和u(k);Step 4.基于滤波器基于滤波器Df(z-1),分别由分别由下下式计算式计算y(k),u(k)的滤波值的滤波值yf(k)和和uf(k);yf(k)=Df(z-1)y(k)uf(k)=Df(z-1)u(k)二、递推型二、递推型GLS法法(3/6)
33、(3/6)Step 5.构造观测数据向量构造观测数据向量f(k-1)15()1-(2)-()1-()1-(2)-()1-()14(0)2),-()1-(1)-(-11)-()13()1-()1-(-)()1-()1-()(ffffff2fffffffffffkkkkkkkkkkkkkykkkPPKIPPKIPKfffff)-(.,),1-()-(.,),1-()1-(bankukunkykyk并依如下并依如下RLS法估计多项式法估计多项式A(z-1)和和B(z-1)的系数的系数,即回归参数向量即回归参数向量f二、递推型二、递推型GLS法法(4/6)(4/6)Step 6.基于式基于式计算相关扰
34、动计算相关扰动v(k)的估计值的估计值,并构造观测数据向量并构造观测数据向量v(k-1).f)1-(-)()(kkykvv)-(.,),1-()1-(dnkvkvk二、递推型二、递推型GLS法法(5/6)(5/6)Step 7.依如下依如下RLS法估计多项式法估计多项式D(z-1)的系数的系数,即回归参数向即回归参数向量量vStep 8.选取上一步估计得的选取上一步估计得的D(z-1)的估计值作进行系统白色的估计值作进行系统白色化处理的线性滤波器化处理的线性滤波器Df(z-1).将循环次数加将循环次数加1,返回第返回第3步继续循环递推辨识步继续循环递推辨识.)18()1-(2)-()1-()1
35、-(2)-()1-()17(0)2),-()1-(1)-(-11)-()16()1-()1-(-)()1-()1-()(vvvvvv2vvvvvvvvvvkkkkkkkkkkkkkvkkkPPKIPPKIPK三、三、GLS仿真仿真(1/14)三、三、广义最小二乘法仿真广义最小二乘法仿真q考虑如下有色噪声扰动的随机线性离散系统考虑如下有色噪声扰动的随机线性离散系统)()()1()()1()(11kvk-nubk-ubk-nya.k-yakybnanba的的GLS在线辨识仿真程序伪代码在线辨识仿真程序伪代码.)()()1()(1kwk-nvdk-vdkvcnd三、三、GLS仿真仿真(2/14)%第
36、一步第一步 初始化初始化输入系统阶次输入系统阶次na,nb和和nd,以及加权因子以及加权因子l l输入系统模型输入系统模型Az=1 a1 a2 和和Bz=0 b1 b2;输入噪声模型输入噪声模型Dz=1 d1 d2 输入系统输入信号输入系统输入信号u(k)的方差的方差 u和噪声和噪声w(k)的方差的方差 w设定系统变量初始值设定系统变量初始值:yf1:na+1=0;uf1:nb+1=0;vf1:nd+1=0;设定辨识变量初始值设定辨识变量初始值:yb1:na+1=0;ub1:nb+1=0;ybf1:na+1=0;ubf1:nb+1=0;vb=1:nc+1=0;q q_yu1:na+nb=0;P
37、_yu=106*I(na+nb,na+nb);q q_v1:nd=0;P_v=106*I(nd,nd);输入初始白化滤波器输入初始白化滤波器Dze=1 de1 de2 三、三、GLS仿真仿真(3/14)%第二步第二步 辨识仿真辨识仿真for k=1:最大仿真步数/%被控对象模型仿真被控对象模型仿真(产生系统输入输出信号产生系统输入输出信号,即数据即数据)yf2:na+1=yf1:na;uf2:nb+1=uf1:nb;vf2:nd+1=vf1:nd;uf1=2*u*(rand()-0.5);wf=2*w*(rand()-0.5);vf1=-Dz2:nd+1*vf2:nd+1+wf;yf1=-2:
38、na+1*yf2:na+1+Bz2:nb+1*uf2:nb+1+vf1;三、三、GLS仿真仿真(4/14)%输入输出数据检测输入输出数据检测ub1=uf1;yb1=yf1;%输入输出滤波输入输出滤波ubf1=Dze(1:nd+1)*ub(1:nd+1);ybf1=Dze(1:nd+1)*yb(1:nd+1);%在线递推辨识系统模型在线递推辨识系统模型 _yu=-ybf(2:na+1)ubf(2:nb+1);K_yu=P_yu*_yu/(l l+_yu*P_yu*_yu);q q_yu=q q_yu+K_yu*ybf(1)-_yu*q q_yu;P_yu=I-K_yu*_yu*P_yu/l l;
39、修正矩阵修正矩阵P_yu；三、三、GLS仿真仿真(5/14)%在线噪声估计在线噪声估计vb1=y(1)+q_yu1:na*yb(2:na+1)-q_yuna+1:na+nb*ub(2:nb+1);%在线递推噪声模型在线递推噪声模型_v=-vb(2:nd+1);K_v=P_v*_v/(l+_v*P_v*_v);q_v=q_v+K_v*vf(1)-_v*q_v;P_v=I-K_v*_v*P_v/l;修正矩阵P_v；输出在线递推参数估计值q_yu,q_v；Dze(2:nd+1)=q_v;三、三、GLS仿真仿真(6/14)%数据移位数据移位yb2:na+1=yb1:na;ub2:nb+1=ub1:nb
40、;ybf2:na+1=ybf1:na;ubf2:nb+1=ubf1:nb;vf2:nd+1=vf1:nd;三、三、GLS仿真仿真(7/14)q 例例1 本例中的仿真对象、实验条件与本例中的仿真对象、实验条件与第四讲中第四讲中例例1相同相同,模型模型结构选择如下结构选择如下y(k)+a1y(k-1)+a2y(k-2)=b1u(k-1)+b2u(k-2)+v(k)(19)v(k)+d1v(k-1)+d2v(k-2)=w(k)(20)q 仿真结果如表仿真结果如表1所示所示.三、三、GLS仿真仿真(8/14)表表1 计算机仿真结果计算机仿真结果参数参数a1a2b1b2d1d2真值真值-1.50.71.
41、00.5-1.00.2噪信比噪信比=73%初始白化滤波器初始白化滤波器Df=1仿真仿真1-1.4990 0.70170.9997 0.5003-1.0056 0.2019仿真仿真2-1.5146 0.71091.0030 0.4869-0.9943 0.2019噪信比噪信比=73%初始白化滤波器初始白化滤波器Df=1+0.5z-1-0.5z-2仿真仿真1-1.4384 0.60300.9870 0.5271-1.1426 0.4060仿真仿真2-1.5740 0.73711.0013 0.4730-0.9114 0.2088噪信比噪信比=23%初始白化滤波器初始白化滤波器Df=1仿真仿真1-1
42、.4667 0.66130.9892 0.5055-1.1649 0.4118仿真仿真2 -1.4628 0.67170.9841 0.5010-1.2130 0.3902三、三、GLS仿真仿真(9/14)q 递推辨识过程的辨识值如下图所示递推辨识过程的辨识值如下图所示噪信比噪信比=73%,初始白化滤波器初始白化滤波器Df=1时递推辨识结果时递推辨识结果(1)参数估计误差的平方和噪声估计误差三、三、GLS仿真仿真(10/14)噪信比噪信比=73%,初始白化滤波器初始白化滤波器Df=1时递推辨识结果时递推辨识结果(2)参数估计误差的平方和噪声估计误差三、三、GLS仿真仿真(11/14)噪信比噪信
43、比=73%,初始白化滤波器初始白化滤波器Df=1+0.5z-1-0.5z-2时递推辨识结果时递推辨识结果(1)参数估计误差的平方和噪声估计误差三、三、GLS仿真仿真(12/14)噪信比噪信比=73%,初始白化滤波器初始白化滤波器Df=1+0.5z-1-0.5z-2时递推辨识结果时递推辨识结果(2)参数估计误差的平方和噪声估计误差三、三、GLS仿真仿真(13/14)噪信比噪信比=23%,初始白化滤波器初始白化滤波器Df=1时递推辨识结果时递推辨识结果(1)参数估计误差的平方和噪声估计误差三、三、GLS仿真仿真(14/14)噪信比噪信比=23%,初始白化滤波器初始白化滤波器Df=1时递推辨识结果时
44、递推辨识结果(2)参数估计误差的平方和噪声估计误差三、三、GLS评述评述(1/1)(1/1)三、三、GLS评述评述q 对对GLSGLS算法算法,有如下评述有如下评述:有色干扰下估计精度较高有色干扰下估计精度较高 迭代收敛较快迭代收敛较快,但是收敛性未得到证明但是收敛性未得到证明 能同时得到过程参数和噪声参数的估计能同时得到过程参数和噪声参数的估计 计算量大计算量大,费机时费机时3.IV法(1/2)3.辅助变量法q 上两讲上两讲分别讨论了同时辨识系统模型和噪声模型的分别讨论了同时辨识系统模型和噪声模型的ELS法和法和GLS法法.在有些建模及其模型的运用在有些建模及其模型的运用问题问题中中,并不需
45、要知道噪声模并不需要知道噪声模型型,即不需要对噪声建模即不需要对噪声建模(辨识辨识).此时此时,若采用若采用ELS法和法和GLS法来辨识法来辨识,需花费较多的计算时需花费较多的计算时间间,而且辨识的参数越多则辨识的精度和效果越差而且辨识的参数越多则辨识的精度和效果越差.这一点这一点,GLS法尤其突出法尤其突出.3.IV法(2/2)q 本讲讨论有相关随机扰动时的随机离散系统参数估计的辅本讲讨论有相关随机扰动时的随机离散系统参数估计的辅助变量助变量(Instrument Variable,IV)法法.该方法不需要辨识系统的噪声模型该方法不需要辨识系统的噪声模型,只要辅助系统选择只要辅助系统选择得恰
46、当得恰当,便可获得较高精度的无偏一致估计便可获得较高精度的无偏一致估计.本讲介绍的主要内容为本讲介绍的主要内容为:成批型算法成批型算法 辅助变量的选择辅助变量的选择 递推型算法递推型算法 IV法的评价法的评价一、成批型算法(1/6)一、成批型算法一、成批型算法q 考虑如下自回归方程考虑如下自回归方程的一致辨识问题的一致辨识问题y(k)=t t(k-1)q q+w(k),k=1,2,.,L (1)或或YL=Lq q+WL (2)其中对扰动量其中对扰动量w(k)未加白噪声的假定未加白噪声的假定,其它符号与其它符号与第三讲第三讲一致一致.下面先考察一般下面先考察一般LS法的参数辨识值的一致收敛性法的
47、参数辨识值的一致收敛性.一、成批型算法(2/6)q 由第三讲的提出由第三讲的提出的的一般一般LS估计式估计式有有LLLLY1LS)()3(11)()()(101001LSLLLLLLLLLLLLLLLWWW一、成批型算法(3/6)由式由式(3)可知可知,欲获得无偏一致的参数估计值欲获得无偏一致的参数估计值,若若)4(1Lim可逆矩阵LLLL)5(1Lim零向量LLLLW 依据依据Frechek定理定理:u若随机序列若随机序列xk,则则f(xk)f()则由条件则由条件(4)与与(5),有有q qLS q q0即即LS估计为无偏估计估计为无偏估计.一、成批型算法(4/6)但但,当当w(k)不为白噪
48、声时不为白噪声时,条件条件零向量LLLLW1Lim并不能保证成立并不能保证成立,相应的相应的LS估计值估计值LLLLLLW10LS11亦不能保证为无偏一致的亦不能保证为无偏一致的.因此下面的问题是因此下面的问题是:如何在有色噪声干扰条件下得到无偏估计？如何在有色噪声干扰条件下得到无偏估计？一、成批型算法(4/6)q 为了对为了对w(k)不为白噪声时不为白噪声时,可获得可获得q q的无偏一致估计值的无偏一致估计值,定义一定义一个如下辅助观测数据矩阵个如下辅助观测数据矩阵*L=*(0),*(1),*(L-1)t t (6)其中辅助变量矩阵须满足其中辅助变量矩阵须满足)7(1Lim*可逆矩阵LLLL
49、)8(1Lim*零向量LLLLW一、成批型算法(5/6)因此因此,使得对不定使得对不定(不相容不相容)方程组方程组YL=Lq q (9)和回归方程和回归方程YL=Lq q+WL可获得如下可获得如下IV解解)10()(*1*IVLLLLY一、成批型算法(6/6)当辅助变量矩阵满足条件当辅助变量矩阵满足条件(7)和和(8)时时,则有则有LLLLLLLLLLLLLLLLLLLWWWY*1*0*1*00*1*1*IV11)()()()(可见可见,只要选择适当的只要选择适当的辅助变量辅助变量,使之满足上述的两个条使之满足上述的两个条件件,IV估计值就可以是无偏一致的估计值就可以是无偏一致的.IV估计法的
50、关键在于如何具体构造出辅助变量矩阵估计法的关键在于如何具体构造出辅助变量矩阵*L.)11(0依据Frechek定理与条件(7)和(8)二、辅助变量的选择(1/9)二、辅助变量的选择二、辅助变量的选择q IV法的关键是选择适当的辅助变量使得条件法的关键是选择适当的辅助变量使得条件(7)和和(8)得到满得到满足足.如果采用如果采用下下图所示的输出图所示的输出x(k)作为辅助变量作为辅助变量,置置*(k-1)=x(k-1),.,x(k-na),u(k-1),.,u(k-nb)(12)那么那么,当当u(k)为高阶持续激励为高阶持续激励信号时信号时,条件条件(7)必定满足必定满足.又因又因x(k)只与只
