现代控制理论极大似然法辨识PPT课件

1、因此，它是基于概率统计基础上的参数估计方法。本章主要讨论极大似然法和递推极大似然法辨识。最后简述模型阶的确定。第1页/共76页第一节 极大似然法辨识极大似然法是以观测值的出现概率为最大作为估计准则。 设有离散随机过程 与未知参数 有关，假设已知条件概率分布密度 。若得到 个独立的观测值 ，则可得分布密度 。要求根据这些观测值估计未知参数 ，估计的准则是观测值 的出现概率为最大。这此定义似然函数 kz/kf zn12nzzz， ， ，12/kf zf zf z， ， kz1211221/nnnniiiL zzzfzfzfzfz， ， ， ，(15-1)第2页/共76页 上式右边是 个概率密度函数

2、的边乘，似然函数L是 的函数。如果L达到极大值， 的出现概率为最大。因此，极大似然法的实质就是求解L达到极大值的 的估值 。为了便于求 ，对式(15-1)等号两边取对数，则把连乘变成连加，即n kz1ln/niiiLfz(15-2) 由于对数函数是单调增加函数，当L取极大时， 也同时取极大值。于是由lnL0L可得 的极大似然估计 。L(15-3)第3页/共76页 现在用极大似然法辨识系统差分方程中的参数。从式(13-6)和式(13-9)可得 11nniiiiy ka y kibu kik 11a zy kb zu kk或(15-4)以及 y(15-5)或NNN y式中下标N表示式(15-5)中

3、的方程个数。(15-6)第4页/共76页根据式(13-14)得残差 11ka zy kb zu k式中12TNe ne ne nNe， ， 假设是 均值为零的高斯分布的不相关随机序列，且与 不相关。 k u kNNNy(15-7)则NNNey(15-8)第5页/共76页 设 服从高斯分布， 具有相同的方差 ，则可得似然函数为Ne e k2 /222/1exp -22NNTNNNNNLLeyyy(15-9)对式(15-9)等号两边取对数 221ln/ln2ln2222TNNNNNNNL y yy(15-10)第6页/共76页 求 对未知参数 和 的偏导数，令偏导数等于零，可得ln/NL y22l

4、n10TTTNNNNL y 222ln022TNNNNLN yy(15-11)(15-12)解式(15-11)，得 的极大似然估计为1TTLNNNN y(15-13)第7页/共76页 从式(13-13)可看出，对于 为高斯白噪声序列这一特殊情况，极大似然估计与一般最小二乘法估计完全相同。 k 实际上， 往往不是白噪声序列，而相关噪声序列。下面讨论在残差相关情况下的极大似然辨识问题。 k解式(15-12)，得 22111n NTNNNNk nekNN yy(15-14)第8页/共76页 是均值为零的高斯分布白噪声序列。多项式 中的各系数 和序列 的均方差 都是未知参数。 k11a zb z，1c

5、 z120212nnnaaabbbccc， ， ， ； ， ， ， ； ， ， ， k将式(15-4)写成 111a zy kb zu kc zk(15-15) 1c zkk111121nnc zc zc zc z 式中(15-16)(15-17)第9页/共76页 式中 为预测误差， 和 分别为 和 的估值。预测误差可用下式表示：e kiiiab、iciiab、ic 11111111nnniiiiiinnnniiiie ky ky ky ka y kibu kic e kia za zy kbu kic e ki 设待估参数1021Tnnnaab bbcc(15-18)设 的预测估值为 y k

6、 111nnniiiiiiy ka y kibu kic e ki (15-19)第10页/共76页式中1111111011111()1()()1nnna za za zb zbb zb zc zc zc z 假设预测误差 服从高斯分布，并且序列 具有相同的方差 。因 与 有关，所以 是被估参数 的函数。 e k e k2 e k111(),c za zb z及2因此，预测误差 满足下式： e k 111c ze ka zy kb zu k(15-20)第11页/共76页式中10211212,TNTNTnnny ny ny nNe ne ne nNaabbbccye， ， ， ， ， ， ，

7、， ，现把式(15-20)写成 111nnniiiiiie ky ka y kibu kic e ki(15-21) 令 ，可得 的N个方程式。把这N个方程式写成向量矩阵形式：1,2,knnnN e kNNNey(15-22)第12页/共76页 11111222121Ny nyu nue ney nyu nue ney nNy Nu nNu Ne nNe N 因为已假设 是零均值的高斯随机序列。所以，极大似然函数 e k/22211/ ,exp -22TNNNNL ye e(15-23)对式(15-23)等号两边取对数，得221ln/ ,ln2ln222TNNNNNL ye e(15-24)第

8、13页/共76页即 22211ln/ ,ln2ln222n NNk nNNLek y求 的偏导数，令其等于零，得2lnL对 2222122211ln102212122n Nk nn Nn Nk nk nLNekekekNNN J式中 2112n Nk nek J(15-26)第14页/共76页 显然，方程(15-20)可理解为预测模型，而 可看作为预测误差。由式(15-26)可知，要使 最小，就要使预测误差的平方和为最小。即使对概率密度不作任何假设，这个准则也是有意义的。因此可按 最小这一准则求 的估值。 e kJJ120212nnnaaabbbccc， ， ， ； ， ， ， ； ， ， ，

9、在估计 时，总是希望 愈小愈好，则22min2NJ(15-27)第15页/共76页 由于 是参数 的线性函数，所以 是这些参数的二次函数。求使L为最大的 ，等价于在式(15-20)的约束条件下求 使 为最小。 e k120212nnnaaabbbccc， ， ， ； ， ， ， ； ， ， ，JJ 在用式(15-20)表示系统模型的情况下，若先算出 的梯度 和海赛矩阵 ，而后用牛顿拉卜森法，可使计算在为简化。下面用牛顿拉卜森法进行迭代计算，求出比较准确的 值。JJ22J0J022J 设 是 的初始值， 和 表示 在 处偏导数和海赛矩阵，则按牛顿拉卜森公式求 的新估值 为0J01012102JJ

10、(15-28)第16页/共76页整个迭代计算步骤如下： 选定初始值 。关于 中的 可按模型式(15-20)有001201,nna aab bb； 111c ze ka zy kb zu k 用最小二乘法求得。而 中的 可先假定某一组值，也可取一组全为零的值。012 ,nc cc第17页/共76页 计算预测误差 e ky ky k给出 2112n Nk nek J并计算 2211n Nk nekN 第18页/共76页 计算 和J22J 21112n Nk nn Nk neke ke k JJ(15-29)式中 111Tnnne ke ke ke ke ke ke kaabbcc 1njjie k

11、e kjy kicaa(15-30)第19页/共76页 1njjiie ke kju kicbb 1njjiie ke kje kiccc (15-31)(15-32)把(15-30)、式(15-31)和式(15-32)改写成 111iiie kc zy kiae kc zu kibe kc ze kic (15-33)(15-34)(15-35)第20页/共76页由以下三式分别得到下列方程组： 101111ijijije ke kije kiaaae ke kije kibbbe ke kije kiccc (15-36) 式(15-33) 、式(15-34) 和式(15-35)都是差分方程

12、，这些差分方程的初始条件都为零。可通过解这组方程求得 关于 的全部偏导数。 和 分别为 ， 和 的线性函数。 e k1naa， ， ，01nnbbcc， ， ，和 iie ke kab， ie kcy kiu kie ki第21页/共76页 利用式(15-34) 、式(15-34) 和式(15-35)，可很方便地求出 关于 的二阶混合偏导数。J 21202111ijiijiijje ke kie kjia caae ke kje kjib cbbe ke kie kjic ccc 下面再求 关于 的二阶偏导数。J 222211Tn Nn Nk nk ne ke ke ke ke k J(15-

13、37)第22页/共76页 的其余二阶偏导数都等于零。从上述三式可看出，二阶偏导数可用一阶偏导数来表示，因此计算比较简单。而且，二阶偏导数可用 ， 和 表示。 e k y k u k e k 当 接近于真值 时， 接近于零。在这种情况下， 接近于零，因此 可用下面的近似式表示： e k 221n Nk ne ke k 22J 221Tn Nk ne ke k J从而使计算更加简单。第23页/共76页 按牛顿拉卜森法计算 的新估值 为1012102 JJ 重复至的计算步骤。经过 为迭代计算之后，可得 进一步迭代可得rr01212rr JJ第24页/共76页 式(15-38)表明，当残差方差的计算误

14、差降到0.01%时，就停止计算。这一方法，即使在噪声比较大的情况下，也能得到较好的估值 。则可停止计算，否则继续迭代计算。如果2241210rrr(15-38)第25页/共76页第二节 递推极大似然法辨识 在线辨识需要用递推极大似然法辨识，下面只讨论用近似方法推导出递推极大似然法的计算公式。设系统的模型为 111a zy kb zu kc zk式中 为预测误差，即 k 1111kcza zy kb zu k(15-39)第26页/共76页 显然， 是模型参数 的函数，所以预测误差可表示为 k101nnnaabbcc， ， ， ， ， ， ， ， kk，由指标函数式(15-26)得21n NNk

15、 nk J，按 最小确定 的估值 。J(15-40)第27页/共76页 如果 是 的线性函数，则可用最小二乘法求 的递推公式。这里，用 的二次型函数逼近 ，从估值 的领域内展开，得k， NJ,Tkkk(15-41)式中，设,ke k则 1111,e kcza zy kb zu k(15-42)第28页/共76页而,TTkk(15-43)参照式(15-30)、式(15-31)和式 (15-32)，可得 的各分量：,e k 111njrjiinjrjiinjrjiie ke kjy kicykiaae ke kju kicukibbe ke kje kicekicc (15-44)(15-45)(

16、15-46)第29页/共76页设 分别表示为kkkyue，和121212TkrrrTkrrrTkrrrykykyknukukuknekekeknyue， ， ， ，则,TTTkkke kyue，从式(15-44)、式(15-45)和式(15-46)可知，只要将输出 ，输入 和误差 进行简单的移位和滤波，就可得到 。 y k u k e k e k第30页/共76页现在用二次型函数逼近 ，即假设存在 和余项 ，使 NJNNP，N 211n NNk nTNNNNk JP，(15-47)根据式(15-41)和式(15-47)可得 2112111n NNk nTTNNNNNNNke JP，(15-48

17、)第31页/共76页设 ，则式(15-48)可写成N 121111112TTTNNNNNNNNee JP式中1111NTNNee nNe(15-49)对上式配平方得 111111TNNNNNrr JP(15-50)第32页/共76页式中11111111121111111TNNNNNNNNTNNNNNNNNreeeePPPP(15-51) 在式(15-50)中， 为已知值。当 时，即当 时， 为极小。所以 的新估计 为1N10NrN 1Nr 1NJ1N11NNNr(15-52)第33页/共76页对式(15-51)应用矩阵求逆引理，得11111111111111TTNNNNNNNNNTNNNNNN

18、NrePPPPPPP(15-53)(15-54)把式(15-53)代入式(15-52)得111NNNNeK(15-55)式中11111111TNNNNNNTNNey nN KPP(15-56)(15-57)第34页/共76页而1 , ,11 , ,1, ,1Ty n Ny Nu n Nu Ne n Ne N， 最后，需要求 的递推关系式。根据 的定义，有1NN与1NN和111011111111NNnnTne nNee nNe nNe nNe nNaabbe nNe nNcc， ， ， ，第35页/共76页101NnnTne nNe nNe nNe nNaabbe nNe nNcc， ， ，由式

19、(15-30)得 的第一行为1N121111111ne Nne Nne Nny Nnccaaae Nca第36页/共76页根据式(15-36)可得12111ne Nne Nne Ne Nnaaaa， ，则12111111ne Nne Nne Ny Nnccaaae Nna第37页/共76页 同理，可得 的其他各行， 与 的递推关系以下式表示： 1N1NN1111100010011000100100100nnNNnccy Nnccu Nncce Nn(15-58)第38页/共76页 方程式(15-53)、(15-55)、(15-56)、和方程式(15-58)为递推极大似在法的一组计算公式。可以证

20、明，这个方法以概率1收敛到估计准则的一个局部极小值。这是一个比较好的方法。第39页/共76页第三节 导弹气动参数的极大似然法辨识建立了导弹数学模型并得到了观测数据后，需要解决模型中的参数估计问题，即所谓的参数辨识问题。由于我们采用的是Bryan的气动力系数的线性化数学模型，因此气动参数辨识只是在已知模型下的参数估计问题。目前在气动参数辨识中应用最广泛的是极大似然法。第40页/共76页一、极大似然法辨识的基本公式一、极大似然法辨识的基本公式 12ttttttxf xunzg xun，设非线性动力学系统为 式中， 为 维状态向量， 是 维输入向量， 是 维观向量， 是 维待估参数向量， 分别为过程

21、噪声和测量噪声。 txn tur tzmp12nn、(15-59)第41页/共76页定义似然函数如下： /2111/2exp2NNTmiiiiiLRttttzzzRzz(15-60) 式中， 为观测向量维数，N是采样点数， 是 在 时刻的估值。对式(15-60)数，并将等号右边“-”变“+”，即得式(15-61)。这样对式(15-60)取极大值，等价于对式(15-61)取极小值。 1111,ln2ln2212NTiiiiimNNttttJRRzzRzzm itz tzit(15-61)第42页/共76页 由上面的指标函数可知，求 的估值问题，归结为对式(15-62)求极小值问题。由于导弹的动力

22、学方程是由彼此间存在耦合的微分方程组描述的，故指标函一般不是气动系数的二次函数，必须通过优化算法来获得参数估计值。通常使用改进的牛顿拉卜森算法。 式(15-61)中的第一项与被估参数 无关，故指标函数可简化为 1111ln22NTiiiiiNttttJRzzRzz(15-62)关于测量噪声 的方差 ，可用下式计算：2nR 11NTiiiiittttN Rzzzz(15-63)第43页/共76页 设已知 的某一估值为 ，相应的指标函数是 ，要求 使指标函数 达到极小值。令kkJkJ22210kkkkkkkJ kJJJO (15-64)式(15-64)中略去 项，可得2k122kkkJJ 11kT

23、NkiJiii zRzz(15-65)(15-66)第44页/共76页二、模型分解后的气动参数辨识二、模型分解后的气动参数辨识 由于导弹飞行动力学模型很复杂，所需辨识的气动参数很多，计算量很大。为了减少计算量，节省CPU时间，必须对导弹的动力学模型进行分解。考虑到气动力系数主要反映在导弹质心运动力学方程上，而力矩系数则主要反映在绕质心旋转的运动分开，分别辨识力系数和力矩系数，也就是将一个六自由度复杂系统的辨识问题分解为两个三自由度简单系统的辨识问题，其中耦合项可作适当处理，这样可以大大减少计算量。第45页/共76页 采用模型分解辨识气动参数可节省CPU时间的另一个原因是，在实际的辨识过程中，力

24、系数比较容易辨识，一般只需迭代四、五次即可达到收敛精度的要求，而力矩系数的辨识则需要花费更多的时间。迭代运算公式如下： 211kTNkiJii zzR(15-67) 11111kTTNNkiiiiiii zzzRRzz(15-68)1kkk上面式中， 称为灵敏度。 iz(15-69)第46页/共76页 牛顿拉卜森算法的具体迭代过程为：根据导弹的风洞实验或理论计算结果给出待估参数的迭代初值 ，再根据导弹的动力学方程求出在 条件下的导弹状态值 、观测值 和灵敏度方程 ，然后由式(15-68)计算出 ，再以 代替原来的 ，重复上述过程。00 x zzk1000 每一次迭代过程都要计算 和 ，当J11

25、kkJJ11kkJJ时，则所得的参数估值 即为所求之值。 为事先给定的精度，通常取 。k0.01(15-70)第47页/共76页 若迭代过程出现发散现象，即 ，就采用变步长的方法进行搜索，即以 代替 ，其中 ，我们取 。因而就花费了许多不必要的时间，使其整个辨识时间要比模型分解法的辨识时间长得多。1kkJJ11kk1kk0112n1,2,3,4n 力系数辨识与力矩系数辨识方法完全相同，本节只介绍力矩系数辨识的过程和结果。下面进行力矩系数的辨识。取状态方程为012xxKRRxxxxxxxQS LLmmmVJV02yyzxRRRyxzyyxyyyyyyJJQS LLmmmmVJJV (15-71)

26、(15-72)第48页/共76页03zzxyRRRzxyzzzzzzzzJJQS LLmmmmVJJV (15-73)观测方程为1xxRW2yyRW3zzRW 式中，下标 表示真值， 和 分别为过程噪声和观测噪声， 由测量值进行估计得到。RiV1,2,3iW i 和(15-74)(15-75)(15-76)第49页/共76页由于导弹是轴对称的，故 。00yzyzzyyzyzyzmmmmmmmm，式(15-71)、式(15-72)和式(15-73)中待估参数共有8个，即000 xxyTxyyzxxymmmmmmm，(15-77)设 为测量值，而 为估计值，则可得到辨识力矩系数时的指标函数为xyz

27、，xyz， 111ln22NTxyzxyziNJiiiiiiRR 11TNxyzxyziiiiiiiN R(15-78)(15-79)第50页/共76页式中 xxxiii yyyiii zzziii(15-80)(15-81)(15-82)设11,2,3ijaijAR，(15-83)第51页/共76页 181112211888818xxxYzNkYYixYzzziiiiiJiiiiiii A 1111888xYzxNkYizxYziiiiJiiiii A(15-84)(15-85)由式(15-66)和(15-67)两式可得第52页/共76页计算中为了求得灵敏度 ，由状态方程两边对 求偏导数，得

28、xYzppp，11,2,32xxxRRRxpxpQS LLmUpptJV，24,5,6,7yyzxxzzxpyppyRRRyppJJtJQS L LmUppJV，35,6,7,8yxyyxzyxpzppRRRzyzpJJtJQS L LmUppJV，(15-86)(15-87)(15-88)第53页/共76页式中1,11,21,11,31,12,422,52,42,62,42,72,43,82,423,53,83,63,81,31,12RRxxRRRxyyRyzRzQS LUUUJLQS LUUUVJUUUULUUUUVUUUULUUV，其它 ,0U i j 第54页/共76页三、参数辨识中的

29、一些具体问题三、参数辨识中的一些具体问题1. 数据窗的选择 数据窗长度一般应尽量长一些，但在实际中数据窗的长度往往受到一定的限制，例如，当有外进入系统，或者系统的动态发生变化时，数据窗不得不缩短；另外记录装置对数据的个数也可能有些限制。这里所研究的是导弹时变，不同时刻导弹速度是不同的，因而各个时刻的气动参数也不同。因此在辨识过程中选取多少个样本点，即多长的数据窗是必须考虑的一个实际问题，若样本点太多，又必然会增长计算时间。通过仿真计算，最后选取了60个样本点，这样既能保证计算精度，又不会使计算时间太长。第55页/共76页关于数据窗的形状也是要考虑的问题，因为马赫数（导弹速度与所在高度音速之比）

30、是变化的，每一次迭代计算中，数据窗内的样本点处于不同的马赫数下。显然，离当前辨识时刻越近的数据中所包含的系统此时刻的信息量就越大；反之，离当前辨识越远的数据，其本身所包含的系统步时刻的信息量变越少。因此在实际计算过程中，按样本点靠近当前时刻的远近对其进行了指数加权处理。具体采用的数形状据窗形状如图15-1所示。第56页/共76页2. 辨识步长的选择一般对于气动参数辨识步长的要求是，相邻两个辨识时刻之间的马赫数变化不能大于0.1马赫，本例所研究的导弹，在其马赫数变化比较剧烈的飞行段内，采用0.1秒的的辨识步长，而其他时间内采用0.3秒的步长。本例采用了变步长和固定步长进行辨识，两种结果表明，采用

31、变步长在保证精度的前，大大缩短了计算时间。第57页/共76页3. 激励信号采用3211波形和方波两种激励信号。激励信号如图15-2所示。第58页/共76页四、辨识结果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.37xxmxxmymyymyym表15-1 “32111”激励，

32、变步长，秒，马赫数=0.73第59页/共76页表15-2 “32111”激励，定步长，秒，马赫数=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.27xxmxxmymyymyym第60页/共76页2表15-3为方波激励，采用变步长得到的一组气动参数辨识部 分结果。表15-1 “32111”激励，变步长，秒，马赫数=0.37参 数真 值估 值相对误差(%)-0.003

33、53123-0.003535540.12-0.3514840-0.35134700.04-0.0236169-0.02389391.17-0.0419191-0.04185220.16-1.769530-0.929184047.49xxmxxmymyymyym第61页/共76页 3结果分析由表15-1、表15-2和表15-3可知，采用变步长和固定步长所得到的气动参数辨识结果基本相同，均能满足精度要求，的辨识误差较大，主要原因是动导数不易辨识。对于正常布局导弹来说，这个气动参数辨识不准基本不影响整个导弹设计。第62页/共76页第四节 模型阶的确定前面所讨论的差分方程的参数辨识方法，都假定模型的结

34、构或差分方程的阶是已知的。实际上，模型的阶往往是求未知的，在这种情况下，存在模型阶的辨识问题。下面只介绍两种常用的确定模型阶的方法，即按残差方差定阶的AIC准则定阶。第63页/共76页一、按残差方差定阶一、按残差方差定阶1. 按估计误差方差最小定阶若线性系统模型为 11a zy kb zu kk式中， 为输出， 为输入，设 的均值为零、方差为 的白噪声序列。用最小二乘法求出 的估值 。 y k u k k2(15-89)第64页/共76页1TTyey 残差为 1121n Nak ne ka zy kb zu kek J(15-90)(15-91)若线性系统模型为 111a zy kb zu k

35、c zk(15-92)则残差为 111niiie ka zy kb zu kc z e k指标函数 仍为式(15-91)。nJ(15-93)第65页/共76页 如图15-3所示，对某一系统，当 随着的 增加而减小。如果 为正确的阶，则 时， 出现最后一次陡峭的下降，而以后 保持不变或只有微小的变化。图15-3所示的系统， 。1,2,nnJn0n01nnnJnJ03n 图15-3第66页/共76页2. F检验方法确定模型的阶 式中，J表示辨识系统的误差平方和。此系统有N对输入输出数据，有 个模型参数。对某一系统的计算结果如表15-4所示。12in由于 随着 的增大而减小，在阶数 的增大过程中，主要是找出使 显著减小的阶数，为此引入准则nJnnnJ111122iiiiiiJJNntJnn(15-94)第67页/共76页 从表15-4可知，当t3时， 的减小是显著的，当t3时， 的减小是不显著的，所以该系统的阶数可选为5。aJaJ第68页/共76页二、按二、按Akaike信息信息(AIC)准则确定系统的阶准则确定系统的阶 这个准则是日本学者Akaike（英译音）总结了时间序列统计建模的经验，借助于信息论所提出的一个合理定阶准则。在一组可供选择的随机模型中，AIC最