《现代控制理论第十五章极大似然法辨识ppt课件》由会员分享，可在线阅读，更多相关《现代控制理论第十五章极大似然法辨识ppt课件（76页珍藏版）》请在金锄头文库上搜索。
1、第十五章 极大似然法辨识 极大似然法是一种得到广泛应用的辨识方法。这种方法要求引入有关随机变量的条件分布密度或似然函数，目的在于建立随机观测数据与未知参数之间的概率特性和统计关系，并通过它求出未知参数的估计值，例如导弹气动参数的估计。因此，它是基于概率统计基础上的参数估计方法。本章主要讨论极大似然法和递推极大似然法辨识。最后简述模型阶的确定。第一节 极大似然法辨识极大似然法是以观测值的出现概率为最大作为估计准则。 设有离散随机过程 与未知参数 有关，假设已知条件概率分布密度 。若得到 个独立的观测值 ，则可得分布密度 。要求根据这些观测值估计未知参数 ，估计的准则是观测值 的出现概率为最大。这
2、此定义似然函数(15-1) 上式右边是 个概率密度函数的边乘，似然函数L是 的函数。如果L达到极大值， 的出现概率为最大。因此，极大似然法的实质就是求解L达到极大值的 的估值 。为了便于求 ，对式(15-1)等号两边取对数，则把连乘变成连加，即(15-2) 由于对数函数是单调增加函数，当L取极大时， 也同时取极大值。于是由可得 的极大似然估计 。(15-3) 现在用极大似然法辨识系统差分方程中的参数。从式(13-6)和式(13-9)可得或(15-4)以及(15-5)或式中下标N表示式(15-5)中的方程个数。(15-6)根据式(13-14)得残差式中 假设是 均值为零的高斯分布的不相关随机序列
3、，且与 不相关。(15-7)则(15-8) 设 服从高斯分布， 具有相同的方差 ，则可得似然函数为(15-9)对式(15-9)等号两边取对数(15-10) 求 对未知参数 和 的偏导数，令偏导数等于零，可得(15-11)(15-12)解式(15-11)，得 的极大似然估计为(15-13) 从式(13-13)可看出，对于 为高斯白噪声序列这一特殊情况，极大似然估计与一般最小二乘法估计完全相同。 实际上， 往往不是白噪声序列，而相关噪声序列。下面讨论在残差相关情况下的极大似然辨识问题。解式(15-12)，得(15-14) 是均值为零的高斯分布白噪声序列。多项式 中的各系数 和序列 的均方差 都是未
4、知参数。将式(15-4)写成(15-15)式中(15-16)(15-17) 式中 为预测误差， 和 分别为 和 的估值。预测误差可用下式表示：设待估参数(15-18)设 的预测估值为(15-19)式中 假设预测误差 服从高斯分布，并且序列 具有相同的方差 。因 与 有关，所以 是被估参数 的函数。 因此，预测误差 满足下式：(15-20)式中现把式(15-20)写成(15-21) 令 ，可得 的N个方程式。把这N个方程式写成向量矩阵形式：(15-22) 因为已假设 是零均值的高斯随机序列。所以，极大似然函数(15-23)对式(15-23)等号两边取对数，得(15-24)即求 的偏导数，令其等于
5、零，得式中(15-26) 显然，方程(15-20)可理解为预测模型，而 可看作为预测误差。由式(15-26)可知，要使 最小，就要使预测误差的平方和为最小。即使对概率密度不作任何假设，这个准则也是有意义的。因此可按 最小这一准则求 的估值。在估计 时，总是希望 愈小愈好，则(15-27) 由于 是参数 的线性函数，所以 是这些参数的二次函数。求使L为最大的 ，等价于在式(15-20)的约束条件下求 使 为最小。 在用式(15-20)表示系统模型的情况下，若先算出 的梯度 和海赛矩阵 ，而后用牛顿拉卜森法，可使计算在为简化。下面用牛顿拉卜森法进行迭代计算，求出比较准确的 值。 设 是 的初始值，
6、 和 表示 在 处偏导数和海赛矩阵，则按牛顿拉卜森公式求 的新估值 为(15-28)整个迭代计算步骤如下： 选定初始值 。关于 中的 可按模型式(15-20)有 用最小二乘法求得。而 中的 可先假定某一组值，也可取一组全为零的值。 计算预测误差给出并计算 计算 和(15-29)式中(15-30)(15-31)(15-32)把(15-30)、式(15-31)和式(15-32)改写成(15-33)(15-34)(15-35)由以下三式分别得到下列方程组：(15-36) 式(15-33) 、式(15-34) 和式(15-35)都是差分方程，这些差分方程的初始条件都为零。可通过解这组方程求得 关于 的
7、全部偏导数。 和 分别为 ， 和 的线性函数。 利用式(15-34) 、式(15-34) 和式(15-35)，可很方便地求出 关于 的二阶混合偏导数。下面再求 关于 的二阶偏导数。(15-37) 的其余二阶偏导数都等于零。从上述三式可看出，二阶偏导数可用一阶偏导数来表示，因此计算比较简单。而且，二阶偏导数可用 ， 和 表示。 当 接近于真值 时， 接近于零。在这种情况下， 接近于零，因此 可用下面的近似式表示：从而使计算更加简单。 按牛顿拉卜森法计算 的新估值 为 重复至的计算步骤。经过 为迭代计算之后，可得 进一步迭代可得 式(15-38)表明，当残差方差的计算误差降到0.01%时，就停止计
8、算。这一方法，即使在噪声比较大的情况下，也能得到较好的估值 。则可停止计算，否则继续迭代计算。如果(15-38)第二节 递推极大似然法辨识 在线辨识需要用递推极大似然法辨识，下面只讨论用近似方法推导出递推极大似然法的计算公式。设系统的模型为式中 为预测误差，即(15-39) 显然， 是模型参数 的函数，所以预测误差可表示为由指标函数式(15-26)得按 最小确定 的估值 。(15-40) 如果 是 的线性函数，则可用最小二乘法求 的递推公式。这里，用 的二次型函数逼近 ，从估值 的领域内展开，得(15-41)式中，设则(15-42)而(15-43)参照式(15-30)、式(15-31)和式 (
9、15-32)，可得 的各分量：(15-44)(15-45)(15-46)设 分别表示为则从式(15-44)、式(15-45)和式(15-46)可知，只要将输出 ，输入 和误差 进行简单的移位和滤波，就可得到 。现在用二次型函数逼近 ，即假设存在 和余项 ，使(15-47)根据式(15-41)和式(15-47)可得(15-48)设 ，则式(15-48)可写成式中(15-49)对上式配平方得(15-50)式中(15-51) 在式(15-50)中， 为已知值。当 时，即当 时， 为极小。所以 的新估计 为(15-52)对式(15-51)应用矩阵求逆引理，得(15-53)(15-54)把式(15-53
10、)代入式(15-52)得(15-55)式中(15-56)(15-57)而 最后，需要求 的递推关系式。根据 的定义，有由式(15-30)得 的第一行为根据式(15-36)可得则 同理，可得 的其他各行， 与 的递推关系以下式表示： (15-58) 方程式(15-53)、(15-55)、(15-56)、和方程式(15-58)为递推极大似在法的一组计算公式。可以证明，这个方法以概率1收敛到估计准则的一个局部极小值。这是一个比较好的方法。第三节 导弹气动参数的极大似然法辨识建立了导弹数学模型并得到了观测数据后，需要解决模型中的参数估计问题，即所谓的参数辨识问题。由于我们采用的是Bryan的气动力系数
11、的线性化数学模型，因此气动参数辨识只是在已知模型下的参数估计问题。目前在气动参数辨识中应用最广泛的是极大似然法。一、极大似然法辨识的基本公式一、极大似然法辨识的基本公式设非线性动力学系统为 式中， 为 维状态向量， 是 维输入向量， 是 维观向量， 是 维待估参数向量， 分别为过程噪声和测量噪声。(15-59)定义似然函数如下：(15-60) 式中， 为观测向量维数，N是采样点数， 是 在 时刻的估值。对式(15-60)数，并将等号右边“-”变“+”，即得式(15-61)。这样对式(15-60)取极大值，等价于对式(15-61)取极小值。(15-61) 由上面的指标函数可知，求 的估值问题，归
12、结为对式(15-62)求极小值问题。由于导弹的动力学方程是由彼此间存在耦合的微分方程组描述的，故指标函一般不是气动系数的二次函数，必须通过优化算法来获得参数估计值。通常使用改进的牛顿拉卜森算法。 式(15-61)中的第一项与被估参数 无关，故指标函数可简化为(15-62)关于测量噪声 的方差 ，可用下式计算：(15-63) 设已知 的某一估值为 ，相应的指标函数是 ，要求 使指标函数 达到极小值。令(15-64)式(15-64)中略去 项，可得(15-65)(15-66)二、模型分解后的气动参数辨识二、模型分解后的气动参数辨识 由于导弹飞行动力学模型很复杂，所需辨识的气动参数很多，计算量很大。
13、为了减少计算量，节省CPU时间，必须对导弹的动力学模型进行分解。考虑到气动力系数主要反映在导弹质心运动力学方程上，而力矩系数则主要反映在绕质心旋转的运动分开，分别辨识力系数和力矩系数，也就是将一个六自由度复杂系统的辨识问题分解为两个三自由度简单系统的辨识问题，其中耦合项可作适当处理，这样可以大大减少计算量。 采用模型分解辨识气动参数可节省CPU时间的另一个原因是，在实际的辨识过程中，力系数比较容易辨识，一般只需迭代四、五次即可达到收敛精度的要求，而力矩系数的辨识则需要花费更多的时间。迭代运算公式如下：(15-67)(15-68)上面式中， 称为灵敏度。(15-69) 牛顿拉卜森算法的具体迭代过
14、程为：根据导弹的风洞实验或理论计算结果给出待估参数的迭代初值 ，再根据导弹的动力学方程求出在 条件下的导弹状态值 、观测值 和灵敏度方程 ，然后由式(15-68)计算出 ，再以 代替原来的 ，重复上述过程。 每一次迭代过程都要计算 和 ，当时，则所得的参数估值 即为所求之值。 为事先给定的精度，通常取 。(15-70) 若迭代过程出现发散现象，即 ，就采用变步长的方法进行搜索，即以 代替 ，其中 ，我们取 。因而就花费了许多不必要的时间，使其整个辨识时间要比模型分解法的辨识时间长得多。 力系数辨识与力矩系数辨识方法完全相同，本节只介绍力矩系数辨识的过程和结果。下面进行力矩系数的辨识。取状态方程
15、为(15-71)(15-72)(15-73)观测方程为 式中，下标 表示真值， 和 分别为过程噪声和观测噪声， 由测量值进行估计得到。(15-74)(15-75)(15-76)由于导弹是轴对称的，故 。式(15-71)、式(15-72)和式(15-73)中待估参数共有8个，即(15-77)设 为测量值，而 为估计值，则可得到辨识力矩系数时的指标函数为(15-78)(15-79)式中(15-80)(15-81)(15-82)设(15-83)(15-84)(15-85)由式(15-66)和(15-67)两式可得计算中为了求得灵敏度 ，由状态方程两边对 求偏导数，得(15-86)(15-87)(15
16、-88)式中其它 三、参数辨识中的一些具体问题三、参数辨识中的一些具体问题1. 数据窗的选择 数据窗长度一般应尽量长一些，但在实际中数据窗的长度往往受到一定的限制，例如，当有外进入系统，或者系统的动态发生变化时，数据窗不得不缩短；另外记录装置对数据的个数也可能有些限制。这里所研究的是导弹时变，不同时刻导弹速度是不同的，因而各个时刻的气动参数也不同。因此在辨识过程中选取多少个样本点，即多长的数据窗是必须考虑的一个实际问题，若样本点太多，又必然会增长计算时间。通过仿真计算，最后选取了60个样本点，这样既能保证计算精度，又不会使计算时间太长。关于数据窗的形状也是要考虑的问题，因为马赫数（导弹速度与所
17、在高度音速之比）是变化的，每一次迭代计算中，数据窗内的样本点处于不同的马赫数下。显然，离当前辨识时刻越近的数据中所包含的系统此时刻的信息量就越大；反之，离当前辨识越远的数据，其本身所包含的系统步时刻的信息量变越少。因此在实际计算过程中，按样本点靠近当前时刻的远近对其进行了指数加权处理。具体采用的数形状据窗形状如图15-1所示。2. 辨识步长的选择一般对于气动参数辨识步长的要求是，相邻两个辨识时刻之间的马赫数变化不能大于0.1马赫，本例所研究的导弹，在其马赫数变化比较剧烈的飞行段内，采用0.1秒的的辨识步长，而其他时间内采用0.3秒的步长。本例采用了变步长和固定步长进行辨识，两种结果表明，采用变
18、步长在保证精度的前，大大缩短了计算时间。3. 激励信号采用3211波形和方波两种激励信号。激励信号如图15-2所示。四、辨识结果辨识结果1. 表15-1和表15-2分别为3211激励，采用变步长和固定步长得到 的一组气动参数辨识部分结果参 数真 值估 值相对误差(%)-0.00353123-0.003530110.03-0.3514840-0.3518500.09-0.0236169-0.023677390.24-0.0419191-0.04207470.37-1.769530-1.12590036.37表15-1 “32111”激励，变步长，秒，马赫数=0.73表15-2 “32111”激励
19、，定步长，秒，马赫数=0.73参 数真 值估 值相对误差(%)-0.00353123-0.003555400.68-0.3514840-0.353778300.65-0.0236169-0.02385561.01-0.0419191-0.0420240.25-1.769530-0.915399048.272表15-3为方波激励，采用变步长得到的一组气动参数辨识部 分结果。表15-1 “32111”激励，变步长，秒，马赫数=0.37参 数真 值估 值相对误差(%)-0.00353123-0.003535540.12-0.3514840-0.35134700.04-0.0236169-0.0238
20、9391.17-0.0419191-0.04185220.16-1.769530-0.929184047.49 3结果分析由表15-1、表15-2和表15-3可知，采用变步长和固定步长所得到的气动参数辨识结果基本相同，均能满足精度要求，的辨识误差较大，主要原因是动导数不易辨识。对于正常布局导弹来说，这个气动参数辨识不准基本不影响整个导弹设计。第四节 模型阶的确定前面所讨论的差分方程的参数辨识方法，都假定模型的结构或差分方程的阶是已知的。实际上，模型的阶往往是求未知的，在这种情况下，存在模型阶的辨识问题。下面只介绍两种常用的确定模型阶的方法，即按残差方差定阶的AIC准则定阶。一、按残差方差定阶一
21、、按残差方差定阶1. 按估计误差方差最小定阶若线性系统模型为式中， 为输出， 为输入，设 的均值为零、方差为 的白噪声序列。用最小二乘法求出 的估值 。(15-89)残差为(15-90)(15-91)若线性系统模型为(15-92)则残差为指标函数 仍为式(15-91)。(15-93) 如图15-3所示，对某一系统，当 随着的 增加而减小。如果 为正确的阶，则 时， 出现最后一次陡峭的下降，而以后 保持不变或只有微小的变化。图15-3所示的系统， 。图15-32. F检验方法确定模型的阶 式中，J表示辨识系统的误差平方和。此系统有N对输入输出数据，有 个模型参数。对某一系统的计算结果如表15-4
22、所示。由于 随着 的增大而减小，在阶数 的增大过程中，主要是找出使 显著减小的阶数，为此引入准则(15-94) 从表15-4可知，当t3时， 的减小是显著的，当t3时， 的减小是不显著的，所以该系统的阶数可选为5。二、按二、按Akaike信息信息(AIC)准则确定系统的阶准则确定系统的阶 这个准则是日本学者Akaike（英译音）总结了时间序列统计建模的经验，借助于信息论所提出的一个合理定阶准则。在一组可供选择的随机模型中，AIC最小的那个模型是一个可取的模型。这个准则的优点就在于它是一个完全客观的准则，应用这个准则时，不需要建模人员作任何主观判断。AIC准则定义为 式中，L是模型的似然函数，
23、是模型中的参数数目，AIC为最小的那个模型就是最佳模型。(15-95)为了简单起见，先考虑下列模型：假设 服从正态布，似然函数为(15-96)取对数后得(15-97)因此求 为最大的 值，根据 可得(15-98)与前述的最小二乘法估计相一致。再根据 可得(15-99)即(15-100) 上式给出了AIC定义中的第一项，该式中 是常数， 。待估的参数为 ，共有 个，即式(15-95)中的 为 。因此(15-101)即(15-102) 选取不同的阶数 ，按式(15-102)计算AIC，可得最优阶数 。如果对不同的 ，当 相近时，应当取 较小的模型。 不管 取何值，上式中常数项 总是不变的。因此可以
24、去掉上式中的常数项，则(15-102)例15-1 设原系统模型为式中 系统的真实参数为： ， 。采样周期0.1秒， 为正态分布 ， 为二位式伪随机序列，记录了300对数据。试用AIC准则验该系统模型的阶数 。 假设模型阶数 ，分别代入式(15-102)进行计算，计算结果如表(15-5)所示。当 时，AIC为最小，因此系统模型的阶数 。AIC148723786-0.99562.101.00224561597-1.9790.9854.904.371.660.973278.6847-2.8512.7170.8651.060.811.050.720.200.034276.0850-2.2781.0800.697-0.4981.081.491.510.471.310.650.210.09表(15-5)
