系统建模灰箱方法.ppt

控制相关专业研究生选修课程系统建模方法马宏军东北大学 信息学院 控制理论与导航技术研究所2013年3月逸夫楼203第第9 9章章系统建模灰箱方法系统建模灰箱方法1、问题的提出、问题的提出1、问题的提出、问题的提出辨识目的：辨识目的：根据过程所提供的测量信息，在某种准则意 义下，估计模型的未知参数。ProcessInputOutputt()203251738895R()76582687394210101032a,b1、问题的提出、问题的提出辨识目的：辨识目的：根据过程所提供的测量信息，在某种准则意 义下，估计模型的未知参数。ProcessInputOutput工程实践 目 的模型结构参数辨识模型校验模型确定m次独立试验的数据1、问题的提出、问题的提出m次独立试验的数据1、问题的提出、问题的提出 未知量的最可能值是使各项实际观测值和计算值之间差的平方乘以其精确度的数值以后的和为最小。1795年，高斯提出的最小二乘的基本原理是1、问题的提出、问题的提出Gauss（1777-1855）使使 最小最小2、最小二乘辨识方法的基本概念、最小二乘辨识方法的基本概念通过试验确定热敏电阻阻值和温度间的关系 当测量没有任何误差时，仅需2个测量值。每次测量总是存在随机误差。2.1 利用最小二乘法求模型参数利用最小二乘法求模型参数根据最小二乘的准则有根据求极值的方法，对上式求导2.1 利用最小二乘法求模型参数利用最小二乘法求模型参数2.2 一般最小二乘法原理及算法一般最小二乘法原理及算法2.2 一般最小二乘法原理及算法一般最小二乘法原理及算法若考虑被辨识系统或观测信息中含有噪声 如果定义2.2 一般最小二乘法原理及算法一般最小二乘法原理及算法2.2 一般最小二乘法原理及算法一般最小二乘法原理及算法2.2 一般最小二乘法原理及算法一般最小二乘法原理及算法2.2 一般最小二乘法原理及算法一般最小二乘法原理及算法2.2 一般最小二乘法原理及算法一般最小二乘法原理及算法2.2 一般最小二乘法原理及算法一般最小二乘法原理及算法 最小二乘法的几何解释2.2 一般最小二乘法原理及算法一般最小二乘法原理及算法 最小二乘法的几何解释2.2 一般最小二乘法原理及算法一般最小二乘法原理及算法最小证明：2.2 一般最小二乘法原理及算法一般最小二乘法原理及算法如果由测量噪声及模型误差等引起的误差如果由测量噪声及模型误差等引起的误差V 的均值为的均值为0，且，且V与输入矢量与输入矢量Hm是统计独立，最小二乘的估计值是无偏的。是统计独立，最小二乘的估计值是无偏的。证明：根据第（1）式的证明，显然有2.2 一般最小二乘法原理及算法一般最小二乘法原理及算法证明：2.2 一般最小二乘法原理及算法一般最小二乘法原理及算法2.2 一般最小二乘法原理及算法一般最小二乘法原理及算法解：由题意得量测方程2.2 一般最小二乘法原理及算法一般最小二乘法原理及算法2.3 加权最小二乘法原理及算法加权最小二乘法原理及算法 一般最小二乘估计精度不高的原因之一是对测量数据同等对待 各次测量数据很难在相同的条件下获得的 有的测量值置信度高，有的测量值置信度低 对不同置信度的测量值采用加权的办法分别对待 置信度高的，权重取得大些；置信度低的，权重取的小些2.3 加权最小二乘法原理及算法加权最小二乘法原理及算法2.2 加权最小二乘法原理及算法加权最小二乘法原理及算法2.3 加权最小二乘法原理及算法加权最小二乘法原理及算法2.3 加权最小二乘法原理及算法加权最小二乘法原理及算法马尔可夫估计马尔可夫估计2.3 加权最小二乘法原理及算法加权最小二乘法原理及算法例例3.2 用用2台仪器对未知标量各直接测量一次，台仪器对未知标量各直接测量一次，量测量分别为量测量分别为z1和和z2，仪器的测量误差均值为，仪器的测量误差均值为0，方，方差分别为差分别为r和和4r的随机量，求其最小二乘估计，并的随机量，求其最小二乘估计，并计算估计的均方误差。计算估计的均方误差。2.3 加权最小二乘法原理及算法加权最小二乘法原理及算法解：由题意得量测方程例3.4 考虑仿真对象选择如下的辨识模型进行一般的最小二乘参数辨识。2.3 加权最小二乘法原理及算法加权最小二乘法原理及算法4阶M序列输出信号一般最小二乘参数辨识流程图3.3 递推最小二乘法原理及算法递推最小二乘法原理及算法 一般最小二乘或加权最小二乘为一次完成算法或批处理算法。计算量大、存储大、不适合在线辨识。采用参数递推估计递推最小二乘算法。3.3 递推最小二乘法原理及算法递推最小二乘法原理及算法3.3 递推最小二乘法原理及算法递推最小二乘法原理及算法3.3 递推最小二乘法原理及算法递推最小二乘法原理及算法如果设则有3.3 递推最小二乘法原理及算法递推最小二乘法原理及算法3.3 递推最小二乘法原理及算法递推最小二乘法原理及算法3.3 递推最小二乘法原理及算法递推最小二乘法原理及算法3.3 递推最小二乘法原理及算法递推最小二乘法原理及算法令3.3 递推最小二乘法原理及算法递推最小二乘法原理及算法3.3 递推最小二乘法原理及算法递推最小二乘法原理及算法3.3 递推最小二乘法原理及算法递推最小二乘法原理及算法3.3 递推最小二乘法原理及算法递推最小二乘法原理及算法3.3 递推最小二乘法原理及算法递推最小二乘法原理及算法3.3 递推最小二乘法原理及算法递推最小二乘法原理及算法例3.5 对3.4采用递推最小二乘估计辨识模型参数 选择如下的辨识模型进行递推最小二乘参数辨识。3.3 递推最小二乘法原理及算法递推最小二乘法原理及算法3.3 递推最小二乘法原理及算法递推最小二乘法原理及算法 数据饱和后，由于递推计算的舍入误差，不仅新的观测值数据饱和后，由于递推计算的舍入误差，不仅新的观测值对参数估计不起修正作用，反而使对参数估计不起修正作用，反而使 失去正定性，导致估计失去正定性，导致估计误差增加。误差增加。数据饱和数据饱和3.3 递推最小二乘法原理及算法递推最小二乘法原理及算法 当系统参数随时间变化时，因新数据被旧数据所当系统参数随时间变化时，因新数据被旧数据所淹没，递推算法无法直接使用。为适应时变参数的情淹没，递推算法无法直接使用。为适应时变参数的情况，修改算法时旧数据的权重况，修改算法时旧数据的权重(降低降低)，增加新数据的，增加新数据的作用。作用。3.3 递推最小二乘法原理及算法递推最小二乘法原理及算法 矩形窗矩形窗3.3 递推最小二乘法原理及算法递推最小二乘法原理及算法 矩形窗矩形窗3.3 递推最小二乘法原理及算法递推最小二乘法原理及算法 指数窗指数窗3.3 递推最小二乘法原理及算法递推最小二乘法原理及算法加加权权最最小小二二乘乘一一般般最最小小二二乘乘可适用时变参数系统只适用于时不变系统3.3 递推最小二乘法原理及算法递推最小二乘法原理及算法 指数窗指数窗参数参数快快时变时变 小小参数参数慢慢时变时变 大大 由上述最小二乘参数辨识的统计特性可知，当在由上述最小二乘参数辨识的统计特性可知，当在量测噪声的均值为量测噪声的均值为0时，才能保其估计值是无偏的。时，才能保其估计值是无偏的。在实际工程和社会系统的辨识中，量测噪声Vm 是 各种系统内外扰动和结构建模误差等因素的综合 反映；Vm不一定为统计独立的白噪声。3.4.1 处理有色噪声扰动的最小二乘类方法处理有色噪声扰动的最小二乘类方法3.4.1 处理有色噪声扰动的最小二乘类方法处理有色噪声扰动的最小二乘类方法 当量测噪声Vm不是统计独立的白噪声。量测噪声Vm是有色噪声，如何获得无偏估计？增广最小二乘法 广义最小二乘法 辅助变量法 多级最小二乘法 偏差补偿最小二乘法不同的有色噪声特性不同的有色噪声模型不同的辨识要求3.4.2 增广最小二乘法原理及算法增广最小二乘法原理及算法平稳相关序列 由关于有色噪声的结论和假设可知，平稳的相关扰动v(k)可被建模如下3.4.2 增广最小二乘法原理及算法增广最小二乘法原理及算法3.4.2 增广最小二乘法原理及算法增广最小二乘法原理及算法 当上述观测数据向量h(k)精确已知时，利用前面讨论的批处理最小二乘法可求得向量 的最小二乘估计值。向量h(k)中包含有不可测的噪声量v(k-1)，.，v(k-n)对自回归模型并不能直接用最小二乘方法.用递推参数估计在线估计在线估计噪声v(k)以实现模型参数在 线递推估计循环估计参数在递推估计过程中，假设当前或前一步的在线参数估计值已相当程度可用的前提下3.4.2 增广最小二乘法原理及算法增广最小二乘法原理及算法利利用用该该参参数数估估计计值值来来在在线线估估计计白白噪噪声声v(k)的的值值 以以替替代代数数据据向量向量h(k)中的白噪声中的白噪声v(k)噪噪声声v(k)的的具具体体的的估估计计算算法法是是如如下下的的事事后后估估计计或或事前估计算法事前估计算法:3.4.2 增广最小二乘法原理及算法增广最小二乘法原理及算法3.4.2 增广最小二乘法原理及算法增广最小二乘法原理及算法 渐消记忆递推增广最小二乘法渐消记忆递推增广最小二乘法3.4.2 增广最小二乘法原理及算法增广最小二乘法原理及算法增广最小二乘参数和噪声v(k)的估计可交替进行计算事事后后估计估计：事事前前估计：估计：3.4.2 增广最小二乘法原理及算法增广最小二乘法原理及算法例3.6 考虑理想数学模型为选择如下的辨识模型进行增广递推最小二乘参数辨识。3.4.2 增广最小二乘法原理及算法增广最小二乘法原理及算法3.4.2 增广最小二乘法原理及算法增广最小二乘法原理及算法下面给出随机线性离散系统在线辨识的伪代码/%第一步 初始化输入系统阶次na,nb和nc,以及加权因子输入系统模型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+1=0;wb=1:nc+1=0;1:na+nb+nc=0;P=106*I(na+nb+nc,na+nb+nc);3.4.2 增广最小二乘法原理及算法增广最小二乘法原理及算法/%第二步 辨识仿真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:nb+1+Cz1:nc+1*wf1:nc+1;3.4.2 增广最小二乘法原理及算法增广最小二乘法原理及算法/%输入输出数据检测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*h/(+h*P*h);=+K*yb(1)-h*;P=I-K*h*P/;修正矩阵P；输出在线递推参数估计值；3.4.2 增广最小二乘法原理及算法增广最小二乘法原理及算法vb1=yb(1)-h*;yb2:na+1=yb1:na;ub2:nb+1=ub1:nb;vb2:nc+1=vb1:nc;3.4.2 增广最小二乘法原理及算法增广最小二乘法原理及算法也可采用事也可采用事前估计前估计例例3.7 考虑如图下所示的仿真对象考虑如图下所示的仿真对象辨辨识识中中,选择选择如下模型如下模型结结构构y(k)+a1y(k-1)+a2y(k-2)=b1u(k-1)+b2u(k-2)+v(k)+c1v(k-1)+c2v(k-2)v(k)是服从均值为零是服从均值为零,方差为方差为1的的 正态分布的不相关随机噪声正态分布的不相关随机噪声;输入信号输入信号u(k)采用伪随机二进制序列采用伪随机二进制序列;通过控制通过控制 v值来改变数据的噪信比值来改变数据的噪信比.3.4.2 增广最小二乘法原理及算法增广最小二乘法原理及算法 计算机仿真结果(噪信比=23%,数据组数1000)参数a1a2b1b2c1c2真值-1.50.71.00.5-1.00.2估计值=1-1.4960.71170.99920.4982-0.90980.1193估计值=0.98-1.4650.69401.06600.5506-0.97780.33693.4.2 增广最小二乘法原理及算法增广最小二乘法原理及算法递推辨识过程的辨识值如下图所示递推辨识过程的辨识值如下图所示遗忘因子=1时递推辨识结果噪声估计误差3.4.2 增广最小二乘法原理及算法增广最小二乘法原理及算法遗忘因子=0.98时递推辨识结果噪声估计误差3.4.2 增广最小二乘法原理及算法增广最小二乘法原理及算法递推辨识过程的辨识值如下图所示递推辨识过程的辨识值如下图所示增广最小二乘法是最小二乘法的一种简单推广只是扩充了参数向量和数据向量h(k)的维数辨识过程模型参数的同时辨识噪声模型就这种意义上说，可称之为增广最小二乘法噪声模型参数估计的收敛过程比过程模型参数估计值的收敛速度慢从实用角度来说，噪声模型阶次不宜取太高3.4.2 增广最小二乘法原理及算法增广最小二乘法原理及算法增广最小二乘算法的特点增广最小二乘算法的特点3.4.2 增广最小二乘法原理及算法增广最小二乘法原理及算法 课后作业课后作业查阅相关系统辨识书籍根据递推最小二乘的工作流程图画出递推增广最小二乘的流程图3.4.2 增广最小二乘法原理及算法增广最小二乘法原理及算法对于有色噪声对于有色噪声v(k)，可通过对噪声建模的方式来使，可通过对噪声建模的方式来使其估计为无偏估计。其估计为无偏估计。噪声模型参数估计比过程模型参数估计的收敛速度慢 噪声模型的阶次不能太高 实际工程中存在一些系统，其噪声模型阶次很高 建模精度和应用比较困难广义最小二乘法广义最小二乘法 引入一个白色滤波器，将相关残差过滤成白引入一个白色滤波器，将相关残差过滤成白色残差。色残差。3.4.2 增广最小二乘法原理及算法增广最小二乘法原理及算法3.4.3 广义最小二乘法原理及算法广义最小二乘法原理及算法形成滤波器未知的、稳定未知的、稳定的、有限阶的的、有限阶的线性滤波器线性滤波器3.4.3 广义最小二乘法原理及算法广义最小二乘法原理及算法3.4.3 广义最小二乘法原理及算法广义最小二乘法原理及算法两边左乘线性滤波器N(z-1)记记为为用最小二乘法估计A(z-1)和B(z-1)利用计算用最小二乘法估计N(z-1)修正滤波器N(z-1)3.4.3 广义最小二乘法原理及算法广义最小二乘法原理及算法3.4.3 广义最小二乘法原理及算法广义最小二乘法原理及算法广义最小二乘算法的计算步骤如下:Step 1.确定模型的结构及A(z-1)、B(z-1)和N(z-1)的阶次；Step 2.选定稳定的稳定的初始滤波器N(z-1)；Step 3.采样获取新的观测数据y(k)和u(k)；3.4.3 广义最小二乘法原理及算法广义最小二乘法原理及算法Step 5.列成如下自回归方程:Step 4.基于滤波器N(z-1),进行如下滤波计算 y(k)=N(z-1)y(k)u(k)=N(z-1)u(k)3.4.3 广义最小二乘法原理及算法广义最小二乘法原理及算法Step 7.计算模型残差的估计值Step 6.用最小二乘法计算 Step 8.计算有色噪声v(k)和白噪声v(k)的自回归方程3.4.3 广义最小二乘法原理及算法广义最小二乘法原理及算法Step 10.修正滤波器N(z-1)Step 9.用最小二乘法计算 Step 11.如满足精度，辨识结束，否则转入Step3.GLS法法的的思思想想是是对对输输入入输输出出数数据据先先进进行行一一次次滤滤波波预预处处理理，然然后利用普通后利用普通LS法对滤波后的数据进行辨识，并反复迭代法对滤波后的数据进行辨识，并反复迭代受滤波模型好坏的影响较大受滤波模型好坏的影响较大滤波模型的好坏也直接与系统模型辨识结果有关系滤波模型的好坏也直接与系统模型辨识结果有关系.从优化理论的角度来说，从优化理论的角度来说，GLS法其实属于非线性优化方法法其实属于非线性优化方法难以避免出现非线性优化中的局部极值点情况难以避免出现非线性优化中的局部极值点情况该方法并不能保证得到的估计值是一致无偏的该方法并不能保证得到的估计值是一致无偏的这是这是GLS法的一个不太令人满意之处法的一个不太令人满意之处.3.4.3 广义最小二乘法原理及算法广义最小二乘法原理及算法思考题：思考题：量测噪声Vm是有色噪声，利用递推广义最小递推广义最小二乘法二乘法辨识参数的步骤分哪几步？能否画出流程图？3.4.3 广义最小二乘法原理及算法广义最小二乘法原理及算法递推广义最小二乘法递推广义最小二乘法递推递推GLSGLS法的基本思想是与成批型法的基本思想是与成批型LSLS法大致相同，不同的是法大致相同，不同的是:由于是递推估计，不能像成批型那样作反复迭代；由于是递推估计，不能像成批型那样作反复迭代；解解决决的的方方法法是是分分别别对对过过程程模模型型和和噪噪声声模模型型两两个个模模型型的的辨辨识识设设计计两两个个递递推推估估计计算算法法，并并在在每每一一个个递递推推步步中中，让让它它们依顺序递推一次；们依顺序递推一次；随随着着递递推推过过程程的的深深入入，将将不不断断改改进进噪噪声声模模型型N(z-1)的的辨辨识结果，同时亦得到较佳的识结果，同时亦得到较佳的A(z-1)和和B(z-1)的辨识结果的辨识结果。3.4.3 广义最小二乘法原理及算法广义最小二乘法原理及算法递推广义最小二乘算法的计算步骤如下:Step 1.确定模型的结构及A(z-1)、B(z-1)和N(z-1)的阶次；Step 2.初始化两个辨识过程，并选定稳定的稳定的初始 滤波器N(z-1)；Step 3.采样获取新的观测数据y(k)和u(k)；递推广义最小二乘法递推广义最小二乘法3.4.3 广义最小二乘法原理及算法广义最小二乘法原理及算法Step 5.列成如下自回归方程:Step 4.基于滤波器N(z-1),进行如下滤波计算 y(k)=N(z-1)y(k)u(k)=N(z-1)u(k)3.4.3 广义最小二乘法原理及算法广义最小二乘法原理及算法Step 7.计算模型残差的估计值Step 6.用最小二乘法计算 Step 8.计算有色噪声v(k)和白噪声v(k)的自回归方程3.4.3 广义最小二乘法原理及算法广义最小二乘法原理及算法Step 9.用最小二乘法计算 Step 10.修正滤波器N(z-1)Step 11.如满足精度，辨识结束，否则转入Step3.考虑如下有色噪声扰动的随机线性离散系统3.4.3 广义最小二乘法原理及算法广义最小二乘法原理及算法%第一步 初始化输入系统阶次na、nb和nc，以及加权因子输入系统模型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:nc+1=0;设定辨识变量初始值:yb1:na+1=0;ub1:nb+1=0;ybf1:na+1=0;ubf1:nb+1=0;vb=1:nc+1=0;_yu1:na+nb=0;P_yu=106*I(na+nb,na+nb);_v1:nc=0;P_v=106*I(nc,nc);输入初始白化滤波器Dze=1 de1 de2 3.4.3 广义最小二乘法原理及算法广义最小二乘法原理及算法%第二步 辨识仿真for k=1:最大仿真步数最大仿真步数/%被控对象模型仿真(产生系统输入输出信号,即数据)yf2:na+1=yf1:na;uf2:nb+1=uf1:nb;vf2:nd+1=vf1:nc;uf1=2*u*(rand()-0.5);wf=2*w*(rand()-0.5);vf1=-Dz2:nc+1*vf2:nc+1+wf;yf1=-2:na+1*yf2:na+1+Bz2:nb+1*uf2:nb+1+vf1;3.4.3 广义最小二乘法原理及算法广义最小二乘法原理及算法%输入输出数据检测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/(+_yu*P_yu*_yu);_yu=_yu+K_yu*ybf(1)-_yu*_yu;P_yu=I-K_yu*_yu*P_yu/;修正矩阵P_yu；3.4.3 广义最小二乘法原理及算法广义最小二乘法原理及算法%在线噪声估计vb1=y(1)+_yu1:na*yb(2:na+1)-_yuna+1:na+nb*ub(2:nb+1);%在线递推噪声模型 _v=-vb(2:nc+1);K_v=P_v*h_v/(+h_v*P_v*h_v);_v=_v+K_v*vf(1)-h_v*_v;P_v=I-K_v*h_v*P_v/;修正矩阵修正矩阵P_v；输出在线递推参数估计值输出在线递推参数估计值 _yu,_v；Dze(2:nd+1)=_v;3.4.3 广义最小二乘法原理及算法广义最小二乘法原理及算法%数据移位yb2:na+1=yb1:na;ub2:nb+1=ub1:nb;ybf2:na+1=ybf1:na;ubf2:nb+1=ubf1:nb;vf2:nc+1=vf1:nc;3.4.3 广义最小二乘法原理及算法广义最小二乘法原理及算法3.4.3 广义最小二乘法原理及算法广义最小二乘法原理及算法例3.8 考虑数学模型的结构为选择如下的辨识模型进行增广递推最小二乘参数辨识。计算机仿真结果计算机仿真结果参数a1a2b1b2c1c2真值-1.50.71.00.5-1.00.2噪信比=73%初始白化滤波器Df=1-1.49900.70170.99970.5003-1.00560.2002噪信比=73%初始白化滤波器Df=1+0.5z-1-0.5z-2-1.43840.60300.98700.5271-1.14260.4060噪信比=23%初始白化滤波器Df=1-1.46670.66130.98920.5055-1.16490.41183.4.3 广义最小二乘法原理及算法广义最小二乘法原理及算法递推辨识过程的辨识值如下图所示噪信比噪信比=73%,初始白化滤波器初始白化滤波器Df=1时递推辨识结果时递推辨识结果噪声估计误差3.4.3 广义最小二乘法原理及算法广义最小二乘法原理及算法噪信比噪信比=73%,初始白化滤波器初始白化滤波器Df=1+0.5z-1-0.5z-2时递推辨识结果时递推辨识结果噪声估计误差3.4.3 广义最小二乘法原理及算法广义最小二乘法原理及算法噪信比噪信比=23%,初始白化滤波器初始白化滤波器Df=1时递推辨识结果时递推辨识结果噪声估计误差3.4.3 广义最小二乘法原理及算法广义最小二乘法原理及算法3.4.3 广义最小二乘法原理及算法广义最小二乘法原理及算法广义最小二乘算法的特点广义最小二乘算法的特点 用于自回归输入模型，是一种迭代的算法 有色干扰下估计精度较高 基本思想是基于对数据先进行一次滤波处理，后利用 普通最小二乘法对滤波后的数据进行辨识 迭代收敛较快,但是收敛性未得到证明 能同时得到过程参数和噪声参数的估计 当过程的输出信噪比比较大或模型参数较多时，数据 白色化处理的可靠性就会下降，辨识结果可能 数据要充分多，否则辨识精度下降。模型阶次不宜过高，初始值对辨识结果有较大影响 计算量大,费机时 ELS和GLS可同时辨识系统模型和噪声模型 在一些实际系统中不需要知道噪声模型，即不 需要对噪声建模(辨识)系统量测噪声为有色噪声系统量测噪声为有色噪声 若采用ELS法和GLS法来辨识,需花费较多的计算 时间,而且辨识的参数越多则辨识的精度和效 果越差，广义最小二乘法尤其突出.问题：问题：系统噪声Vm是有色噪声，不需要对噪声 建模，还想达到无偏估计，如何做？3.4.3 广义最小二乘法原理及算法广义最小二乘法原理及算法3.4.4 辅助变量最小二乘法原理及算法辅助变量最小二乘法原理及算法辅助变量（辅助变量（Instrument Variable,IV)最小二乘算法最小二乘算法 系统噪声Vm是有色噪声，不需要对噪声建模，引入辅助系统，只要辅助系统选择恰当，可获得高精度的无偏估计。有偏估计3.4.4 辅助变量最小二乘法原理及算法辅助变量最小二乘法原理及算法（非奇异矩阵）（非奇异矩阵）（无偏估计）（无偏估计）零均值白噪声3.4.4 辅助变量最小二乘法原理及算法辅助变量最小二乘法原理及算法当当量测噪声Vm是有色噪声不一定成立不一定成立不一定是无偏估计不一定是无偏估计问题：问题：在有色噪声Vm下，如何获得无偏估计？3.4.4 辅助变量最小二乘法原理及算法辅助变量最小二乘法原理及算法定义如下辅助观测矩阵，并使下列极限成立定义如下辅助观测矩阵，并使下列极限成立（非奇异矩阵）（非奇异矩阵）问题：问题：如何构造或选择辅助变量？3.4.4 辅助变量最小二乘法原理及算法辅助变量最小二乘法原理及算法辅助变量的选择辅助变量的选择 3.4.4 辅助变量最小二乘法原理及算法辅助变量最小二乘法原理及算法辅助变量的选择辅助变量的选择。选取3.4.4 辅助变量最小二乘法原理及算法辅助变量最小二乘法原理及算法辅助变量的选择辅助变量的选择。（无偏估计）（无偏估计）辅助变量的选取还有更简单的方法辅助变量的选取还有更简单的方法辅助变量最小二乘法有递推形式，思路类同辅助变量最小二乘法有递推形式，思路类同3.4.4 辅助变量最小二乘法原理及算法辅助变量最小二乘法原理及算法递推辅助变量最小二乘法递推辅助变量最小二乘法 Step1：确定被辨识模型的结构及多项式确定被辨识模型的结构及多项式A(z-1)和和B(z-1)的阶次；的阶次；Step2：确定或设计所采用的辅助变量系统；确定或设计所采用的辅助变量系统；Step3：设定递推参数初值设定递推参数初值(0)，P(0)；Step4：采样获取新的观测数据采样获取新的观测数据y(k)和和u(k),并组成观测数据并组成观测数据h(k)；Step5：计算辅助变量计算辅助变量x(k),并组成辅助变量观测数据向量并组成辅助变量观测数据向量h*(k)；Step6：用递推辅助变量最小二乘法计算当前参数递推估计值；用递推辅助变量最小二乘法计算当前参数递推估计值；Step7：循环次数循环次数k加加1，然后转回到第，然后转回到第4步继续循环。步继续循环。3.4.4 辅助变量最小二乘法原理及算法辅助变量最小二乘法原理及算法递推辅助变量最小二乘法的步骤递推辅助变量最小二乘法的步骤 下面给出针对随机线性离散系统，给出辨识伪代码%第一步 初始化输入系统阶次na,nb和nc,以及辅助变量方法变量IV_method,辅助变量系统滞后IV-d，加权因子输入系统模型Az=1 a1 a2 和Bz=0 b1 b2;输入噪声模型Cz=1 c1 c2 输入系统输入信号u(k)的方差u和噪声v(k)的方差w设定系统变量初始值:yf1:na+1=0;uf1:nb+1=0;wf1:nc+1=0;设定辨识变量初始值:yb1:d+1=0;ub1:d+1=0;xi1:d+1=0;1:na+nb=0;P=106*I(na+nb,na+nb);3.4.4 辅助变量最小二乘法原理及算法辅助变量最小二乘法原理及算法/%第二步 辨识仿真for k=1:最大仿真步数最大仿真步数/%被控对象模型仿真(产生系统输入输出信号,即数据)yf2:na+1=yf1:na;uf2:nb+1=uf1:nb;wf2:nc+1=uf1:nc;uf1=2*u*(rand()-0.5);wf1=2*w*(rand()-0.5);yf1=-Az2:na+1*yf2:na+1+Bz2:nb+1*uf2:nb+1+Cz1:nc+1*wf1:nc+1;3.4.4 辅助变量最小二乘法原理及算法辅助变量最小二乘法原理及算法%输入输出数据检测ub1=uf1;yb1=yf1;%在线递推辨识过程仿真 _star=-xi(2:na+1)ub(2:nb+1);=-yb(2:na+1)ub(2:nb+1);K=P*h_star/(+h*P*h_star);=+K*yb(1)-h*;P=I-K*h*P/;输出在线递推参数估计值输出在线递推参数估计值；3.4.4 辅助变量最小二乘法原理及算法辅助变量最小二乘法原理及算法%辅助变量计算if IV_method=1%自适应滤波自适应滤波 xi(1)=h_star*;elseif IV_method=2%纯滞后纯滞后 xi(1)=ub(d+1);else%Tally原理原理 xi(1)=yb(d+1);endyb2:d+1=yb1:d;ub2:d+1=ub1:d;xi2:d+1=xi1:d;3.4.4 辅助变量最小二乘法原理及算法辅助变量最小二乘法原理及算法3.4.4 辅助变量最小二乘法原理及算法辅助变量最小二乘法原理及算法例3.9 考虑数学模型的结构为选择如下的辨识模型进行增广递推最小二乘参数辨识。参数a1a2 b1b2真值1.50.71.00.5纯滞后法IV估计值d=3-1.44470.64221.00690.5291d=4-1.57590.76551.02700.3981d=5-1.44030.67191.01270.5797Tally原理IV估计值d=3-1.47790.68681.01770.5291d=4-1.51320.70661.04710.4723d=5-1.49260.68620.98540.4910计算机仿真结果(噪信比=73%,C(z-1)=1-1.0z-1+0.2z-2)3.4.4 辅助变量最小二乘法原理及算法辅助变量最小二乘法原理及算法d=3时的纯滞后法辅助变量时的纯滞后法辅助变量递推辨识结果递推辨识结果参数估计误差的平方和3.4.4 辅助变量最小二乘法原理及算法辅助变量最小二乘法原理及算法d=4时的纯滞后法辅助变量时的纯滞后法辅助变量递推辨识结果递推辨识结果参数估计误差的平方和3.4.4 辅助变量最小二乘法原理及算法辅助变量最小二乘法原理及算法d=5时的纯滞后法辅助变量时的纯滞后法辅助变量递推辨识结果递推辨识结果参数估计误差的平方和3.4.4 辅助变量最小二乘法原理及算法辅助变量最小二乘法原理及算法d=3时的时的Tally原理辅助变量原理辅助变量递推辨识结果递推辨识结果参数估计误差的平方和3.4.4 辅助变量最小二乘法原理及算法辅助变量最小二乘法原理及算法d=4时的时的Tally原理辅助变量原理辅助变量递推辨识结果递推辨识结果参数估计误差的平方和3.4.4 辅助变量最小二乘法原理及算法辅助变量最小二乘法原理及算法d=5时的时的Tally原理辅助变量原理辅助变量递推辨识结果递推辨识结果参数估计误差的平方和3.4.4 辅助变量最小二乘法原理及算法辅助变量最小二乘法原理及算法3.4.4 辅助变量最小二乘法原理及算法辅助变量最小二乘法原理及算法辅助变量最小二乘（辅助变量最小二乘（IVLS）法的特点）法的特点 IVLS的计算量比LS估计增加不多，在相关噪声情 况下估计精度确有明显改善；IVLS法与递推IVLS法思路相似，但不等价；递推IVLS法的计算量与RLS法接近，在相关噪声 情况下估计精度优于RLS法；递推IVLS法对初值P(0)的选择非常敏感，为了提 高递推IV法的可靠性，在递推的前几十步最好用 RLS法过渡.3.5.1 多变量系统的最小二乘辨识的原理多变量系统的最小二乘辨识的原理 MIMO系统3.5.1 多变量系统的最小二乘辨识的原理多变量系统的最小二乘辨识的原理 MIMO系统的子系统 3.5.1 多变量系统的最小二乘辨识的原理多变量系统的最小二乘辨识的原理 3.5.1 多变量系统的最小二乘辨识的原理多变量系统的最小二乘辨识的原理 待辨识的参数：3.5.1 多变量系统的最小二乘辨识的原理多变量系统的最小二乘辨识的原理 3.5.2 多变量系统的最小二乘辨识的算法与设计多变量系统的最小二乘辨识的算法与设计)()()1()()()1()(101kVnkUBkUBkUBnkYAkYAkYnn+-+-+=-+-+LL3.5.2 多变量系统的最小二乘辨识的算法与设计多变量系统的最小二乘辨识的算法与设计)()()1()()()1()(101kVnkUBkUBkUBnkYAkYAkYnn+-+-+=-+-+LL3.5.2 多变量系统的最小二乘辨识的算法与设计多变量系统的最小二乘辨识的算法与设计 令，可得到个方程，并令3.5.2 多变量系统的最小二乘辨识的算法与设计多变量系统的最小二乘辨识的算法与设计 3.5.2 多变量系统的最小二乘辨识的算法与设计多变量系统的最小二乘辨识的算法与设计 令，可得到个方程，并令3.5.2 多变量系统的最小二乘辨识的算法与设计多变量系统的最小二乘辨识的算法与设计 3.5.2 多变量系统的最小二乘辨识的算法与设计多变量系统的最小二乘辨识的算法与设计 多变量系统的最小二乘辨识的算法的递推形式多变量系统的最小二乘辨识的算法的递推形式3.5.2 多变量系统的最小二乘辨识的算法与设计多变量系统的最小二乘辨识的算法与设计 例3.8 采用多变量系统的最小二乘辨识方法辨识如下MIMO 系统的参数3.5.2 多变量系统的最小二乘辨识的算法与设计多变量系统的最小二乘辨识的算法与设计 加速度计加速度计Z Z Z ZX X X XY Y Y Y陀螺仪陀螺仪陀螺仪陀螺仪x xy y0 0运动运动运动运动轨迹轨迹轨迹轨迹x x1 1y y1 14、最小二乘应用、最小二乘应用惯性器件标定惯性器件标定陀螺仪陀螺仪陀螺仪陀螺仪 加速度计加速度计加速度计加速度计 4、最小二乘应用、最小二乘应用惯性器件标定惯性器件标定零偏标度因数输出轴灵敏度误差系数摆轴灵敏度误差系数二阶非线性误差系数4、最小二乘应用、最小二乘应用惯性器件标定惯性器件标定4、最小二乘应用、最小二乘应用惯性器件标定惯性器件标定4、最小二乘应用、最小二乘应用景像匹配景像匹配4、最小二乘应用、最小二乘应用景像匹配景像匹配4、最小二乘应用、最小二乘应用景像匹配景像匹配摄摄像像机机坐坐标标系系ZXYOMmxypfy0Sx0o4、最小二乘应用、最小二乘应用摄像机标定摄像机标定4、最小二乘应用、最小二乘应用摄像机标定摄像机标定4、最小二乘应用、最小二乘应用摄像机标定摄像机标定标定结果：标定结果：内参数内参数标标定定结结果（果（单单位：像素）位：像素）偏心偏心x00+1.16245054039601910000e+001偏心偏心y00-4.54959287921496480000e+001焦距焦距fx+3.01024572283878210000e+003焦距焦距fy+3.00913305414917840000e+003畸畸变变系数系数k0+1.95850462242873280000e-008畸畸变变系数系数k1-2.32895263837667420000e-01516.59mm 16.56mm 16.4mm 4、最小二乘应用、最小二乘应用摄像机标定摄像机标定指指标标X方向方向Y方向方向均均值值（mm）0.0000860.000091标标准差（准差（mm）0.7934241.531048最大最大值值（mm）1.9562253.421633最小最小值值（mm）-2.36758-3.16439X方方向向误误差差Y方向误差方向误差4、最小二乘应用、最小二乘应用摄像机标定摄像机标定4、最小二乘应用、最小二乘应用摄像机标定摄像机标定标定原图校正后单级倒立摆示意图 4、最小二乘应用、最小二乘应用单级倒立摆单级倒立摆图中所示变量名的物理含义如表图中所示变量名的物理含义如表1所示。所示。4、最小二乘应用、最小二乘应用单级倒立摆单级倒立摆MruFPNS