1、5.3 测试法建模测试法建模5.3.1 测试法建模的方法测试法建模的方法 机理法建模能够解决实际生产中一部分过程的建模问题。但是，还有很多生产过程由于工艺的复杂性、产品本身在加工中的变化性（如物理或化学变化），使得用机理法建模在技术上遇到了极大的困难，人们不得不考虑用其它的方法建模。从学科角度看，建立数学模型应该属于系统辨识（System Identification）与参数估计（Parametric Estimation）的范畴。简单地说，系统辨识主要是对被研究对象的结构进行判断，解决“是什么”的问题，例如一阶惯性环节，二阶系统；而参数估计则对支撑结构的参数进行估计，解决“是多少”的问题。事
2、实上，有很多比较复杂的过程，我们对其工作机理并不清楚，更难以用数学和物理的方法加以具体描述，此时用测试法建模是一个不得已而为之的方法。与前述的机理法建模相比，测试法建模不需要深入了解过程的工作机制，通常的做法是将其看作一个“黑箱”，通过从外部施加适当的输入信号，测得过程的输出信号，通过对这些输入和输出信号的处理和研究，获得其动态特性和数学模型。因此，问题主要归纳为：1）施加何种输入信号才能最大限度地激励被测过程，使得动态特性得以充分表现，并通过输出信号显露出来？2）对获得的数据或波形，通过什么方法和技术才能估算出适用于控制用的动态模型？一般说来，模型有非参数模型（Nonparametric M
3、odel）和参数模型（Parametric Model）之分。建立非参数模型的方法通常有：时域法（Time-domain Method）、频域法（Frequency-domain Method）和统计相关法（Statistical Correlation Method）等，这类建模不需要事先确定模型的结构，可用于广泛的被控过程；获得参数模型的方法主要有：最小二乘法（Least Square Method）、极大似然法（Maximum Likelihood Method）和梯度校正法（Gradient Correction Method）等，这类建模需要假设一定的模型结构，通过极小化模型与过程之
4、间的误差准则，来确定相应的模型参数。用时域法测定被控过程的数学模型：对过程施加阶跃信号，或者方波信号，测取响应曲线，并由此确定过程的传递函数。该方法具有测试简单、需用设备少的优点，但测试精度不高，其获得的模型可用于一般工业过程控制；用频域法测定被控过程的数学模型：对过程施加不同频率的正弦波输入信号，获得相应的输出幅值与相位，由此可得到该过程的频率特性，由频率特性获得传递函数。该方法需用专门的频率发生和测试设备，模型精度比用时域法高；用统计相关法测定被控过程的数学模型：对被控过程施加伪随机信号，采用统计相关法获得过程的动态特性。它的特点是，可在生产状态下施加随机信号，并测取相关数据，精度较高，但
5、需获得较多数据，并借助计算机协助处理。最小二乘法又称最小平方法，是估计离散时间数学模型参数的一常用种方法。随着计算机技术在控制中的应用，最小二乘法在过程辨识的实践中被越来越广泛地采用。本章主要讨论用时域响应法和最小二乘法获取数学模型。5.3.2 时域响应曲线法时域响应曲线法 1响应曲线的测取 时域响应曲线法是对被控过程施加阶跃信号，如果被控过程不允许长期施加阶跃信号，则改用矩形脉冲信号，然后测取响应曲线，并由此求取输入和输出之间的传递函数。（1）阶跃响应曲线的测取 当被控过程稳定之后，对调节阀施加一个幅值合适的阶跃信号，用记录仪或数据采集系统记录被控量的变化曲线，例如被控量为温度时，就记录温度
6、响应曲线，直到变化曲线进入稳定状态为止。下面几点是在这一过程中值得注意的：a)在施加阶跃信号之前，被控过程应处在较为稳定的工作状态，在下一轮施加输入信号时，应等前一过程结束、并恢复稳态一段时间之后进行。b)施加阶跃信号的幅度：通常为正常输入信号的 515，以不影响正在进行的生产为好。同时，幅度也不能太小，因为过小的输入容易被其它信号淹没，在响应曲线上难以表现出来。c)多次、全面测试，消除偶然性，获得真实结果：试验不仅应在相同条件下，重复几次，以获得两次及其以上的较为接近的响应曲线，而且也应选取不同负荷、不同输入值，测得相应的响应曲线，以获得全面的动态特性。（2）矩形脉冲相应曲线的测取 在有些情
7、况下，用阶跃信号输入时，可能会危及安全生产，或者影响产品的质量和数量。此时，输入可考虑用矩形脉冲信号代替阶跃信号，测得过程的矩形脉冲响应曲线。由于通过阶跃响应曲线求传递函数被人们所熟悉，所以，往往将矩形脉冲响应曲线再转化为阶跃响应曲线，进而按阶跃响应曲线法确定传递函数。图 5-12a)为矩形脉冲输入信号，它可以分解为图 b)所示的两个阶跃信号的叠加，即1211()()()()()u tu tu tu tu ta其中 a 为脉冲宽度，。如果被控过程是线性的，则其矩形脉冲响应曲线 可分解为阶跃响应曲线 和 ，即21()()u tu ta()y t1()y t2()y t图 5-12 矩形脉冲及其响
8、应分解图 1211()()()()()y ty ty ty ty ta见图5-12 c)所示。这里，和 分别为 和 的响应。于是，起源于零点的阶跃响应为：1()y t2()y t1()u t2()u t11()()()y ty ty ta（5-22）其中，为矩形脉冲响应，为起源于 点的阶跃响应（注意此时的符号为正）。式（5-22）为由矩形脉冲响应求阶跃响应的公式，可用作图法逐步求出：t 在 0a 时，当 时，此时的 前面已经有，为已知，所以可求得 ，如此类推，一直进行到进入稳态。()y t1()y tata1()()y ty t2ata 11()()()y ty ty ta1()y ta()y
9、 t1()y t 由无自平衡能力过程的矩形脉冲响应曲线转化为阶跃响应曲线，也可通过作图法实现，这可作为练习，留给读者来完成（见思考题与习题 5-5）。2由过程阶跃响应曲线确定传递函数 由阶跃响应曲线确定传递函数，通常需要确定传递函数的结构及其参数两部分。结构形式是指被控过程的传递函数形式，生产过程主要是：一阶惯性环节、二阶惯性环节、或等 n 阶惯性环节，并且这些环节时常含有纯滞后，其表达形式为()1KG sTs12()(1)(1)KG sT sT s()1nKG sTs(2,3,)n，()e1sKG sTs12()e(1)(1)sKG sT sT s()e1snKG sTs，对于无自平衡能力的
10、过程，也有类似的形式 1()G sTs121()(1)G sT s T s121()(1)nG sT s T s1()esG sTs121()e(1)sG sT s T s121()e(1)snG sT s T s，传递函数的参数是伴随结构形式出现的待定常数，如一阶惯性、具有纯时延、有自平衡能力的传递函数含有：K、T 和 三个需要确定的参数。关于传递函数结构形式的确定，主要有两方面的考虑：一是根据对被控过程的经验和知识（即通常所说的先验知识）来确定；二是根据控制的要求，尽量将一个原本较复杂的过程用低阶的传递函数来近似描述，因此产生的误差只要处在可接受的范围即可。下面的讨论，集中在参数的确定上。
11、（1）由阶跃响应曲线求一阶惯性加纯时延环节的参数 这里传递函数形式为e1sKTs（5-23）它有三个参数需要确定，即放大系数 K、时间常数 T 和纯时延 。原本就是一阶惯性加纯时延过程的阶跃响应曲线见图 5-5，其 T 和 从图中很容易确定，放大系数为 （其中 为阶跃输入信号的幅值）。下面讨论原本是二阶及其以上过程，且响应曲线呈“s”形，如何用式（5-23）来近似描述。000()/Khx0 x 现有阶跃响应曲线如图 5-13 所示，试图用式（5-23）来近似描述，需确定 K、T 和 三个参数。显然，放大系数可用下式求得0()(0)yyKx（5-24）图 5-13 由阶跃响应曲线确定 T 和 图
12、5-14 纵坐标的标么化其中 为输入幅值，为已知量。0 x 关于 T 和 的确定，有两种方法：一是作图法，二是计算法，下面分别介绍。用作图法求 T 和 ：首先找到响应曲线上凹和下凹的交接点-拐点 D，过 D 点作曲线的切线，切线与时间轴 t 相交于 A 点，与 相较于 C 点，该点在时间轴上的投影为 B，则OA为 ，AB 为 T。()y 作图法的问题是：曲线的拐点有时不容易找到，并且作切线时，有一定的随意性。所以，用作图法求 T 和 ，可能会因人而异，有一定的误差。计算法求 T 和 ：将阶跃响应曲线的纵坐标 标么化，即：实际值基准值，这里取基准值为 ，于是 的标么值（Unit Value）为(
13、)y t0()yKx()y t()()()y ty ty前面的图 5-13 则变为图 5-14。该标么化处理，并不改变响应曲线的横坐标和形状，仅方便求得参数。我们的目的是要通过图5-14所示的阶跃响应曲线，容易求出式（5-23）中的 T 和 。式（5-23）的阶跃响应标么化后，输出为 0()1etTty tt（5-25）为求 T 和 ，在图 5-14 曲线上取两点：E 和 F ，且 ，则有 11,()ty t22,()ty t21tt11()1etTy t 22()1etTy t 由此解出2112ln 1()ln 1()ttTy ty t211212ln 1()ln 1()ln 1()ln 1
14、()ty tty ty ty t当然，为了计算上的方便，也可取 ，代入上两式，有 1()0.39y t2()0.632y t212()Ttt122tt，算出 T 和 后，可检验一下用式（5-25）与实测曲线的误差大小，如果误差可接受，则所求的传递函数式（5-23）可用。否则，应考虑用其它型传递函数（例如高阶传递函数）来描述。具体方法为，另取三点：、和 ，具体为3t4t5t3t40.8tT52tT，由式（5-25）算得3()0y t4()0.55y t5()0.865y t，并分别与图5-13中 、对应的纵坐标比较即可。3tt4t5t （2）由阶跃响应曲线求二阶惯性及其以上环节的参数 当你用一阶
15、惯性环节近似被控过程传递函数，检验发现误差不能满足原定的精度时，可考虑二阶及其以上的惯性环节传递函数。设有阶跃响应曲线如图 5-15，现在，欲用二阶惯性环节的传递函数 12()(1)(1)KG sT sT s（5-26）来近似描述它，其中，k 、和 为待定的参数。1T2T图 5-15 阶跃响应曲线图 5-16 具有纯时延的阶跃响应曲线 当输入为 时，该传递函数的响应为 0()x tx12/1202121()1eet Tt TTTy tKxTTTT12/122121()1eet Tt TTTyTTTT（5-27）在图 5-15 所示的阶跃响应曲线上，找出两点：A(，)和 B(，),并将这两点分别
16、代入式（5-27），有 1t1()0.4()y ty2t2()0.8()y ty1112/122121ee0.6t Tt TTTTTTT 2122/122121ee0.2tTtTTTTTTT 其近似解为121212.16TTtt1 2122121.740.55TTttTT 研究表明，由式（5-27）表示的阶跃响应，应有 ，并且当 时，被控过程应为一阶惯性环节 ，且时间常数为 120.317/0.46tt12/0.317tt/(1)KTs 122.12ttT当 时，被控过程可为二阶等容惯性环节 ，且 12/0.46tt 2/(1)KTs 124.36ttT当 时，被控过程应为二阶以上惯性环节，可
17、用等 n 容惯性环节来描述 12/0.46tt(1)nKTs 其中 n 和 T 分别按下式计算 21211.0750.5tntt122.16ttTn如果算得 n 不为整数，应取最接近的整数。n 与 的关系也可用表 5-1 表示。12/tt表表5-1 多容过程的多容过程的 n 与与 之间的关系之间的关系 12/tt12ttn1234567891012140.3170.4600.5340.5840.6180.6400.6660.6840.6990.7120.7340.751 如果阶跃响应曲线有明显的纯时延，如图 5-16 所示，则应在式（5-26）右边乘上一个纯时延环节：，变为 es12()e(1
18、)(1)sKG sT sT s其中 见图5-16。具体用上面的公式求 、时，应在 和 中减去 时间段。1T2T1t2t 对于以上传递函数中的放大系数 K，仍可用式（5-24）求取。（3）由无自平衡过程的阶跃响应曲线求过程参数 当阶跃响应的曲线如图 5-17 所示时，该过程的传递函数则具有无自平衡特性，其特点是，随着 ，响应曲线的变化速率逐渐趋于某一常数。t 图 5-17 无自平衡过程阶跃响应曲线 根据该响应曲线，该过程可用1()esG sTs（5-28）来近似。当阶跃信号 作用于输入端时，其输出为0()1()x txt0()()xy ttTt()0y t 0t 下面将讨论 T 和 的确定。作阶
19、跃响应直线部分的延长线（见图中虚线段），与 t 轴相交于点 ，该线与时间轴 t 夹角为 ，于是 2t2t 由于 ，所以 0tan/yxT 001tan/xxTyt其中，和 见图 5-17。1yt 由图5-17可知，用式（5-28）近似原过程的最大误差发生在 这一段曲线上，即响应的起始段。为此可再加一个惯性环节来减小误差，即采用下列传递函数来描述过程 1Qt 11()e(1)sG sTs T s（5-29）其中 T 的确定如上所示，即不变，而 和 的确定如下：1T1t121Ttt 显然，用式（5-29）描述响应曲线表示的过程，比式（5-28）要精确些。5.4 基于最小二乘法的过程辨识基于最小二乘
20、法的过程辨识 应该说，用最小二乘法建模仍然是一种测试法建模，但含有较多的处理技巧与方法，这里将其单列为一节，主要是考虑内容稍多、篇幅较大。最小二乘法是系统辨识中的一种常用参数估计方法，它具有原理明了、算法简捷、收敛较快、相对容易理解的特点，因而被广泛用于参数估计之中。最小二乘法包括：最小二乘的批处理法、递推法、渐消记忆法和增广法等。5.4.1 离散时间系统模型离散时间系统模型 数学模型分为连续时间系统模型和离散时间系统模型。前面讨论的是连续时间系统模型，随着计算机的普及与应用，离散时间系统模型被越来越重视，最小二乘法采用的是离散时间系统模型。在这两类模型中按是否有随机扰动，每类又可分为确定性和
21、随机性两种形式。确定性离散系统（Deterministic Discrete Systems）的单输入/输出方程形式：11()()()()dA zy kzB zu k(1,2,)k（5-30）式中 aa11212()1nnA za za za z bb1120120()(0)nnB zbb zb zb zbd 为纯时延，且 ，和 分别为k时刻的输出和输入。1d()y k()u k 随机离散系统（Stochastic Discrete Systems）输入-输出差分模型一般有下面几种：自回归滑动平均（Auto-Regressive Moving Average，ARMA）模型 111()()()
22、()()()dA zy kzB zu kC zk（5-31）式中 、与式（5-30）中的相同，1()A z1()B zcc11212()1nnC zc zc zc z 为白噪声（White Noise）序列，且()kE()0k2()E()()0()ijijij，式（5-31）右边第 1 项称为滑动平均项，左边项称为自回归项。式（5-31）也可写为 1111()()()()()()()dzB zC zy ku kkA zA z这里11()()dzB zA z11()()C zA z，分别被称为过程模型和噪声模型，后者也被称为成形滤波器。而11()()()()C ze kkA z可以看作是白噪声经
23、线性环节的输出，它一般为有色噪声（Coloured Noise）。自回归积分滑动平均（Auto-Regressive Integrated Moving Average，ARIMA）模型111()()()(1)()()/A zy kB zu kC zk或者111()()()(1)()()A zy kB zu kC zk这里，、和 与式（5-31）中相同。与式（5-31）相比，这里假定 。当 时，多项式中前 项的系数为零。11z 1()A z1()B z1()C z1d 1d 1()B z1d 最小二乘模型（Least Square Model）111()()()()()()()()dA zy
24、kzB zu kkB zu kdk 或者111()1()()()()()dzB zy ku kkA zA z式中，、d 和 有与前面相同的含义。1()A z1()B z()k 滑动平均（Moving Average，MA）模型 11()()()()()dy kzB zu kC zk与式（5-31）相比，这里有：。1()1A z5.4.2 批处理最小二乘法批处理最小二乘法 考虑最小二乘模型 11()()()()()A zy kB zu kdk式中，为白噪声。aa11212()1nnA za za za z bb112012()nnB zbb zb zb z()k 设已知 、，现在的任务是根据可量
25、测的输入和输出，确定参数：anbna12naaa,b012nbb bb,，由输入-输出模型有 1201()(1)(2)()()(1)()()abnanby ka y ka y ka y knb u kdb u kdb u kdnk （5-32）令：Tab()(1),(2),(),(),(1),()ky ky ky knu kdu kdu kdn ab1(1)nnR为观测向量；为待估参数向量，则式（5-32）可写为另一形式abT1201,nnaaabbb 1(1)abnnRT()()()y kkk（5-33）现有 N 次观测数据组(),(),1,2,y iu iiN并且ab1Nnn 当 时，由式
26、（5-32）有：1,2,kNab12a0b(1)(0)(12)(1)(1)(1)(1)nnya ya ya ynb udb udn ab12a0b(2)(1)(0)(2)(2)(2)(2)nnya ya ya ynb udb udn ab12a0b()(1)(2)()()()()nny Na y Na y Na y Nnb u Ndb u NdnN 引入下列符号abTT1(1)T(1)(1)(2)(2),()()NNnnyyRRy NN Y，1(1)(2)NR由式（5-33），有下列矩阵形式：Y 由于真实的参数向量 并不知道，不妨用 来表示它的估计值，于是，基于 的输出估计为 T()()y k
27、k式中，abT1201,nnaaabbb 现在定义残差 （也是一随机变量）为实际输出与估计输出之差：()kTT()()()()()()ky ky kkkkT()()()kk（5-34）对于 ，则有1,2,kN()Y 其中 。T1(1)(2)()NNR 现在的任务是：求使目标函数TT()J Y()Y为最小的 （记为 ）。LS 展开上式，有TTTTT2()J Y YY 由求极值的方法，对 求一阶导，并令其为零 TT220J Y从而有T1TLS()Y （5-35）由于二阶导2T220J 所以由式（5-35）求得 的为极小值。式（5-35）即为批处理法的最小二乘估计。LS 最小二乘估计的统计特性讨论：
28、由于 为白噪声序列，故有()k E02EI，（1）无偏性：使 J 为最小的参数估计向量 的数学期望为参数真值向量，即 LS LSE 这是因为T1TT1TLSEE()()E Y它揭示最小二乘参数估计是围绕参数真值波动的统计性质。（2）估计误差（偏差）协方差 1T2LSLSE()()主对角线上各元表现参数估计的散度，非对角线上各元反映参数估计 分量相互影响程度或相关性大小。LS 上式的成立，是因为TTT1TT1TLSLSEE()()YY1TTTT1E()()()YY11TT2T I12T （3）最小方差估计：设 为 的任一其他线性无偏估计，则 TTLSLSEE()()即最小二乘估计是最小方差估计，
29、也就是参数估计 离参数真值 最近。LS 由于 为 的任一线性无偏差估计，所以 可表示为 PY式中，。由于，即 ab(1)nnNRP EEPYP PI从而TTTT2TEE E E()()E()()PYPPYPPYYPPP因为T1T1TTT10()()()PPPP所以2T2T1()PPTTLSLSEE （4）一致收敛性：若 T1lim()NN R存在且正定，则 是一致收敛的，即 LS LSlimN 由于T12LSLS0E 所以12TLSLS211lim ElimlimNNNNNN 0R 由 的无偏性知 LSE 0 LS 所以LSlimN 例例 5-1 现有最小二乘模型 11()()()()()dA
30、 zy kzB zu kk其中 、，为零均值白噪声，实验获得输入输出数据如表 5-2，试用批处理最小二乘法确定多项式 和 的参数。2an 1bn 1d()k1()A z1()B z 解解：首先组成矩阵 ，然后按式（5-31）求出。下面用Matlab 语言来做。u=-1;-1;-1;-1;-1;1;1;-1;-1;1;-1;1;1;-1;1;1;1;1;-1;1;-1;1;-1;-1;-1;1;-1;-1;1;1;1;%将u 输入到工作空间 y=0.0582;-0.6395;-1.8510;-2.0242;-1.4363;-0.9188;0.3053;2.2938;1.0769;-2.2943;
31、-1.9656;0.4587;1.3710;1.8783;0.2454;-1.1223;0.7848;2.4983;2.2147;-0.2424;-1.5523;-0.5707;0.5078;0.7394;-1.4378;-2.6328;-0.5359;1.4520;-0.4325;-1.2545;1.1510;%将y 输入到工作空间 表表5-2 实验获得的输入输出数据实验获得的输入输出数据 tuytuytuy0.0005.5-1-1.96561 1.01-0.57070.5-10.05826.010.45871 1.5-10.50781.0-1-0.63956.511.37101 2.0-1
32、0.73941.5-1-1.85107.0-11.87831 2.5-1-1.43782.0-1-2.02427.510.24541 3.01-2.63282.5-1-1.43638.01-1.12231 3.5-1-0.53593.01-0.91888.510.78481 4.0-11.45203.510.30539.012.49831 4.51-0.43254.0-12.29389.5-12.21471 5.01-1.25454.5-11.0769101-0.24241 5.511.15105.01-2.294310.5-1-1.5523-phi=0;-1*y(1:end-1),0;0;-
33、1*y(1:end-2),0;u(1:end-1),0;0;u(1:end-2);%组建 Theta=inv(phi*phi)*phi*y；%计算 1Y 运行结果为：T=-0.5076，0.6075，0.6854，0.7947 这与真值：-0.5,0.6,0.7,0.8 相差不大。另外，也可以用Matlab的辨识工具箱获得结果。T=iddata(y,u,0.5);%处理数据 G=arx(T,2,2,1);%调用辨识函数 运行结果为：Discrete-time IDPOLY model:A(q)y(t)=B(q)u(t)+e(t)A(q)=1-0.5077 q-1+0.6075 q-2 B(q)
34、=0.6863 q-1+0.7947 q-2Estimated using ARX from data set T Loss function 0.00146694 and FPE 0.00190159Sampling interval:0.5其中 即为 ，符号 q 就是本书中的符号 z。由此可见，两结果几乎相同。q-iiz(1,2)i 5.4.3 递推最小二乘法递推最小二乘法 前面介绍的批处理法需要大量的计算机内存（随着N增大），而且当 变化时，不能自动跟踪其变化，实时性不好。而递推最小二乘法（Recursive Least Square Method）则能解决这一问题。LS 引理：设A、C
35、 和 均为非奇异方阵，则有()ABCD1111111()()ABCDAA B CDA BDA 证明证明：用 左乘上式右端，若结果为单位矩阵，则说明其等式关系成立。()ABCD111111()()ABCD AA B CDA BDA1111111111()()IBCDAB CDA BDABCDA B CDA BDA111111()()IBCDAB ICDA B CDA BDA1111111()()IBCDABC CDA B CDA BDA I当 时，有 CI111111()()ABDAA B IDA BDA（5-36）递推最小二乘算法递推最小二乘算法 最小二乘模型描述的系统 11()()()()(
36、)dA zy kzB zu kk基于 及以前的输出 y、及以前的输入 u，未知参数向量(1)k(1)kd abT1201,nnaaabbb 的最小二乘估计 的递推公式为 T(1)()(1)(1)(1)()kkky kkkK T()(1)(1)1(1)()(1)kkkkkk PKPT(1)(1)(1)()kkkk PIKP（5-37）（5-38）（5-39）式中T1()()k PT(1),(2),()k，T()(1),(2),(),(),(1),()abky ky ky knu kdu kdu kdn 为观测向量，为未知参数向量 的估计值，为增益向量。(1)k K 证明：证明：设 是基于 和 的
37、估计，根据批处理最小二乘法有()k YT1T()()k Y对于 时的最小二乘估计为 1k T1T1111(1)kkkkkY式中 ab(1)(1)1T(1)knnkRk (1)11(1)kkRy k YY，于是TT1T(1)(1)(1)(1)(1)kkkky k Y（5-40）令TT1(1)(1)(1)kkkP 并视 、和 分别为式（5-36）中的 A、B、D，再考虑 在此为标量，由引理有 T (1)k T(1)k 1IDA BT1TT1T1TT1()(1)(1)()(1)()1(1)()(1)kkkkk P再令 ，则有 T1()()k PTT()(1)(1)()(1)()1(1)()(1)kk
38、kkkkkkkPPPPP 又令（5-41）T()(1)(1)1(1)()(1)kkkkkk PKP即为式（5-39），且式（5-41）可写为TT(1)()(1)(1)()(1)(1)()kkkkkkkk PPKPIKP它就是式（5-39）。将式（5-39）代入式（5-40），并注意交换标量和矩阵的位置，有TTTTTTT(1)(1)(1)()(1)(1)(1)(1)()(1)(1)()(1)(1)()(1)(1)()(1)(1)()(1)(1)1(1)()(1)(1)()kkkkky kkkkky kkkkkkkkky kkky kkP kkkkIKPYIKPIKIKPPIK T(1)1(1)(
39、)(1)kkkkP TT()(1)(1)()(1)(1)1(1)()(1)kkkkky kkkkPPPT()(1)(1)()(1)(1)kkkkky kKKT()(1)(1)(1)()kky kkkK即为式（5-37）。递推最小二乘算法的说明：（1）从式（5-37）看，新的参数向量估计值 为先前的参数向量估计值 加修正项(1)k ()k T(1)(1)(1)()ky kkkK 在不断更新过程中，、和 的行列数不变，但它们的旧数据不断被新数据替换；T(1)k (1)k K （2）为增益向量，为误差的协方差阵，一般 与 成正比，协方差越大，说明估计值与真值相差越大，增益向量也会越大，所产生的校正作
40、用也越大；(1)k K(1)k P(1)k P(1)k K （3）初值 和 的确定。(0)(0)P 方法1：若已有 组数据，则可批处理它们，并将结果作为初值，即 ab1NnnT1T(0)()YT1(0)()P，方法2：，其中 。(0)0(0)PI61010105.4.4 具有遗忘因子的递推最小二乘法具有遗忘因子的递推最小二乘法 递推最小二乘法有一个缺点：常常出现“数据饱和”。随着k的增加，和 变得越来越小，式（5-37）中的修正项对 的修正能力变得越来越弱，即新近加入的输入/输出数据对参数向量估计值的更新作用不大。这样导致的结果是：参数估计值难以接近真值；当参数真值时变时，该算法无法跟踪这种变
41、化，从而使实时参数辨识失败。(1)k P(1)k K(1)k 解决该问题的方法之一是用具有遗忘因子（Forgetting Factor）的递推最小二乘法。取性能指标函数：T()()J YW Y式中，为加权对角阵：W001 WN为观测数据组数，为遗忘因子：。01 按与前面相同的思路，可推出具有遗忘因子的递推最小二乘估计公式：T(1)()(1)(1)(1)()kkky kkkK T()(1)(1)(1)()(1)kkkkkk PKPT1(1)(1)(1)()kkkk PIKP式中 。T1()()kPW 几点说明：（1）当 时，该估计公式组即为递推最小二乘算法公式组；（2）的选取范围一般在：，参数变
42、化快时，取小点；变化慢时，取大些；10.950.99 （3）初值的选取与前面的递推最小二乘法相同。5.4.5 递推增广最小二乘法递推增广最小二乘法 以上用的是最小二乘模型，在很多情况下，所以需考虑更一般的情况，即ARMA模型：1()1C z111()()()()()()A zy kB zu kdC zk（5-42）其中，11212()1aannA za za za z 112012()bbnnB zbb zb zb z11212()1ccnnC zc zc zc z,式（5-42）可变为 111()1()()()()()1()()y kA zy kB zu kdC zkk将其表示为 T()()
43、()y kkk式中T()(1),(2),(),(),(1),aky ky ky knu kdu kd c,(),(1),(2),()bu kdnkkkncT120112,abnnnaaabbbccc 在 中，由于 不可测，所以只能用其估计值 来替换，设 T()k T()()()()()ky ky ky kk式中Ta()(1),(2),(),(),(1),ky ky ky knu kdu kd bc,(),(1),(2),()u kdnkkknabcT120112,nnnaaabbbccc 用与前面类似的方法，并考虑用 代替 中的 ，则有下列递推公式 T()k()y kT()k T()(1)()
44、()()(1)kkky kkkK T1()(1)()1()(1)()kkkkkkKPPT()()()(1)kkkkPIKP 的估计可用下列方法之一来进行：（1）预测形式：()kT()()(1)()()(1)ky ky k ky kkk （2）滤波形式：T()()()()()()ky ky k ky kkk 并且当 时，有 。0k()0k 前面介绍了最小二乘参数估计算法和特性，但是估计参数是否收敛到真实值尚未提及，由于这方面的内容涉及较多的定理及证明，这里只能给出结论，证明可参考有关资料。被估参数收敛定理被估参数收敛定理：若 是 N 阶持续激励（Persistently Exciting of
45、Order N），即充分丰富（Sufficient Rich），则式（5-33）式（5-35）所示的算法能保证指数收敛。()u k 例例 5.2 现有自回归滑动平均模型 111()()()()()()dA zy kzB zu kC zk已知 ，为白噪声：零均值，方差为 0.1，取采样周期 ，实验获得输入/输出数据如表 5-3，试用递推增广最小二乘法确定对象参数：、和 ，并将参数估计变化过程用图形表示出来。2,1,2,2abcnnnd()k1sT 1a2a0b1b1c2c表表5-3 输入输入/输出数据输出数据u1.00001.00000.00000.0000-1.0000-1.00001.000
46、01.0000y0.11650.08021.02512.34462.13011.5708-0.1742-1.6821u0.00000.00001.00001.00000.00000.0000-1.0000-1.0000y-0.90941.11281.32800.92811.69902.56812.41131.4639 解解：由递推增广最小二乘法公式，有参数估计算法：(1)初始化 、等，并给 u 和 y 赋值；(2)；(3)计算 ；(4)计算 ；(5)保存参数 、和 ；(6)计算 ；(7)求滤波形式的 ；(8)移位处理 u、y、和 ；(9)k 步终了吗?否，去第 2 步；是，去第 10 步；(1
47、0)打印结果。1kk(0)P(0)T()(1),(2),(2),(3),(1),(2)ky ky ku ku kkk T1()(1)()1()(1)()K kP kkk P kkT()(1)()()()(1)kkK ky kkk1()a k2()ak0()b k1()b k1()c k2()c kT()()()(1)P kIK kk P kT()()()()ky kkk 按算法编写程序，经调试，得出 、和 随拍数变化的规律，具体参见图5-18。具体程序为（用 Matlab 编写）：1a2a0b1b1c2cclearN=40;u=1,1,0,0,-1,-1,1,1,0,0,1,1,0,0,-1,
48、-1,1,1,0,0,1,1,0,0,-1,-1,1,1,0,0,1,1,.0,0,-1,-1,1,1,0,0;%录入输入数据 y=0.1165,0.0802,1.0251,2.3446,2.1301,1.5708,-0.1742,-1.6821,-0.9094,1.1128,1.3280,0.9281,1.6990,2.5681,2.4113,1.4639,-0.1749,-1.9893,-1.1059,1.0812,1.4497,0.9304,1.5019,2.6898,2.3573,1.5460,-0.1111,-1.7658,-0.9011,0.8914,1.4770,1.1226,1
49、.6800,2.4018,2.3039,1.4445,-0.1399,-1.9036,-0.9657,1.1687;%录入输出数据 I=eye(6);P=1e9*I;Q0=zeros(6,1);h0=1;t=0;%初始化 y1=0;y2=0;u1=0;u2=0;u3=0;x1=0;x2=0;for j=1:N f=-y1;-y2;u2;u3;x1;x2;%构成 T()(1),(2),(2),(3),(1),(2)ky ky ku ku kkk K=P*f/(1+f*P*f);Q=Q0+K*(y(j)-f*Q0);%T()(1)()()()(1)kkky kkkK a1(j)=Q(1);a2(j
50、)=Q(2);b0(j)=Q(3);b1(j)=Q(4);c1(j)=Q(5);c2(j)=Q(6);P=(eye(6)-K*f)*P;x(j)=y(j)-f*Q;%求 T()()()()ky kkk y2=y1;y1=y(j);u3=u2;u2=u1;u1=u(j);x2=x1;x1=x(j);Q0=Q;tt(j)=t+h0;t=tt(j);end plot(tt,a1,m,tt,a2,g,tt,b0,r,tt,b1,b,tt,c1,-.,tt,c2,-)%画参数图 从最后一拍看，有 ，和 ，与真值-0.8、0.15、1、0.5、-0.65 和 0.1 相差不大，估计成功。该题也可用 Mat
