过程控制技术-第5章课件2

1、5.3 测试法建模5.3.1 测试法建模的方法 机理法建模能够解决实际生产中一部分过程的建模问题。但是，还有很多生产过程由于工艺的复杂性、产品本身在加工中的变化性（如物理或化学变化），使得用机理法建模在技术上遇到了极大的困难，人们不得不考虑用其它的方法建模。 从学科角度看，建立数学模型应该属于系统辨识（System Identification）与参数估计（Parametric Estimation）的范畴。简单地说，系统辨识主要是对被研究对象的结构进行判断，解决“是什么”的问题，例如一阶惯性环节，二阶系统；而参数估计则对支撑结构的参数进行估计，解决“是多少”的问题。 事实上，有很多比较复杂的

2、过程，我们对其工作机理并不清楚，更难以用数学和物理的方法加以具体描述，此时用测试法建模是一个不得已而为之的方法。与前述的机理法建模相比，测试法建模不需要深入了解过程的工作机制，通常的做法是将其看作一个“黑箱”，通过从外部施加适当的输入信号，测得过程的输出信号，通过对这些输入和输出信号的处理和研究，获得其动态特性和数学模型。因此，问题主要归纳为： 1）施加何种输入信号才能最大限度地激励被测过程，使得动态特性得以充分表现，并通过输出信号显露出来？ 2）对获得的数据或波形，通过什么方法和技术才能估算出适用于控制用的动态模型？ 一般说来，模型有非参数模型（Nonparametric Model）和参数

3、模型（Parametric Model）之分。建立非参数模型的方法通常有：时域法（Time-domain Method）、频域法（Frequency- 5.3 测试法建模5.3.1 测试法建模的方法 -domain Method）和统计相关法（Statistical Correlation Method）等，这类建模不需要事先确定模型的结构，可用于广泛的被控过程；获得参数模型的方法主要有：最小二乘法（Least Square Method）、极大似然法（Maximum Likelihood Method）和梯度校正法（Gradient Correction Method）等，这类建模需要假设一

4、定的模型结构，通过极小化模型与过程之间的误差准则，来确定相应的模型参数。 用时域法测定被控过程的数学模型：对过程施加阶跃信号，或者方波信号，测取响应曲线，并由此确定过程的传递函数。该方法具有测试简单、需用设备少的优点，但测试精度不高，其获得的模型可用于一般工业过程控制； 用频域法测定被控过程的数学模型：对过程施加不同频率的正弦波输入信号，获得相应的输出幅值与相位，由此可得到该过程的频率特性，由频率特性获得传递函数。该方法需用专门的频率发生和测试设备，模型精度比用时域法高； 用统计相关法测定被控过程的数学模型：对被控过程施加伪随机信号，采用统计相关法获得过程的动态特性。它的特点是，可在生产状态下

5、施加随机信号，并测取相关数据，精度较高，但需获得较多数据，并借助计算机协助处理。 最小二乘法又称最小平方法，是估计离散时间数学模型参数的一常用种方法。随着计算机技术在控制中的应用，最小二乘法在过程辨识的实践中被越来越广泛地采用。 本章主要讨论用时域响应法和最小二乘法获取数学模型。5.3.2 时域响应曲线法 1响应曲线的测取 -domain Method）和统计相关法（Statisti时域响应曲线法是对被控过程施加阶跃信号，如果被控过程不允许长期施加阶跃信号，则改用矩形脉冲信号，然后测取响应曲线，并由此求取输入和输出之间的传递函数。 （1）阶跃响应曲线的测取 当被控过程稳定之后，对调节阀施加一个

6、幅值合适的阶跃信号，用记录仪或数据采集系统记录被控量的变化曲线，例如被控量为温度时，就记录温度响应曲线，直到变化曲线进入稳定状态为止。下面几点是在这一过程中值得注意的： a) 在施加阶跃信号之前，被控过程应处在较为稳定的工作状态，在下一轮施加输入信号时，应等前一过程结束、并恢复稳态一段时间之后进行。 b) 施加阶跃信号的幅度：通常为正常输入信号的 515，以不影响正在进行的生产为好。同时，幅度也不能太小，因为过小的输入容易被其它信号淹没，在响应曲线上难以表现出来。 c) 多次、全面测试，消除偶然性，获得真实结果：试验不仅应在相同条件下，重复几次，以获得两次及其以上的较为接近的响应曲线，而且也应

7、选取不同负荷、不同输入值，测得相应的响应曲线，以获得全面的动态特性。 （2）矩形脉冲相应曲线的测取 在有些情况下，用阶跃信号输入时，可能会危及安全生产，或者影响产品的质量和数量。此时，输入可考虑用矩形脉冲信号代替阶跃信号，测得过程的矩形脉冲响应曲线。由于通过阶跃响应曲线求传递函数被人们所熟悉，所以，往往将矩形脉冲响应曲线再转化为阶跃响应曲线，进而按阶跃响应曲线法确定传递函数。 图 5-12a) 为矩形脉冲输入信号，它可以分解为图 b) 所示的两个阶跃信号的叠加，即时域响应曲线法是对被控过程施加阶跃信号，如果被控过程不允许长其中 a 为脉冲宽度， 。如果被控过程是线性的，则其矩形脉冲响应曲线 可

8、分解为阶跃响应曲线 和 ，即图 5-12 矩形脉冲及其响应分解图 见图5-12 c)所示。这里， 和 分别为 和 的响应。于是，起源于零点的阶跃响应为： （5-22） 其中， 为矩形脉冲响应， 为起源于 点的阶跃响应（注意此时的符号为正）。式（5-22）为由矩形脉冲响应求阶跃响应的公式，可用作图法逐步求出： t 在 0a 时， ，当 时， ，此时的 前面已经有， 为已知，所以可求得 ，如此类推，一直进行到进入稳态。 其中 a 为脉冲宽度， 由无自平衡能力过程的矩形脉冲响应曲线转化为阶跃响应曲线，也可通过作图法实现，这可作为练习，留给读者来完成（见思考题与习题 5-5）。 2由过程阶跃响应曲线确

9、定传递函数 由阶跃响应曲线确定传递函数，通常需要确定传递函数的结构及其参数两部分。结构形式是指被控过程的传递函数形式，生产过程主要是：一阶惯性环节、二阶惯性环节、或等 n 阶惯性环节，并且这些环节时常含有纯滞后，其表达形式为 ， ， ， ， 对于无自平衡能力的过程，也有类似的形式 ， ， ， ， 传递函数的参数是伴随结构形式出现的待定常数，如一阶惯性、具有纯时延、有自平衡能力的传递函数含有：K、T 和 三个需要确定的参数。 由无自平衡能力过程的矩形脉冲响应曲线转化为阶 关于传递函数结构形式的确定，主要有两方面的考虑：一是根据对被控过程的经验和知识（即通常所说的先验知识）来确定；二是根据控制的要

10、求，尽量将一个原本较复杂的过程用低阶的传递函数来近似描述，因此产生的误差只要处在可接受的范围即可。下面的讨论，集中在参数的确定上。 （1）由阶跃响应曲线求一阶惯性加纯时延环节的参数 这里传递函数形式为（5-23） 它有三个参数需要确定，即放大系数 K、时间常数 T 和纯时延 。原本就是一阶惯性加纯时延过程的阶跃响应曲线见图 5-5，其 T 和 从图中很容易确定，放大系数为 （其中 为阶跃输入信号的幅值）。下面讨论原本是二阶及其以上过程，且响应曲线呈 “s” 形，如何用式（5-23）来近似描述。 现有阶跃响应曲线如图 5-13 所示，试图用式（5-23）来近似描述，需确定 K、T 和 三个参数。

11、 显然，放大系数可用下式求得（5-24） 关于传递函数结构形式的确定，主要有两方面的 图 5-13 由阶跃响应曲线确定 T 和 图5-14 纵坐标的标么化其中 为输入幅值，为已知量。 关于 T 和 的确定，有两种方法：一是作图法，二是计算法，下面分别介绍。 用作图法求 T 和 ：首先找到响应曲线上凹和下凹的交接点 - 拐点 D，过 D 点作曲线的切线，切线与时间轴 t 相交于 A 点，与 相较于 C 点，该点在时间轴上的投影为 B，则OA为 ，AB 为 T。 作图法的问题是：曲线的拐点有时不容易找到，并且作切线时，有一定的随意性。所以，用作图法求 T 和 ，可能会因人而异，有一定的误差。 计算

12、法求 T 和 ：将阶跃响应曲线的纵坐标 标么化，即：实际值基准值，这里取基准值为 ，于是 的标么值（Unit Value）为 图 5-13 由阶跃响应曲线确定 T 和 前面的图 5-13 则变为图 5-14。该标么化处理，并不改变响应曲线的横坐标和形状，仅方便求得参数。我们的目的是要通过图5-14所示的阶跃响应曲线，容易求出式（5-23）中的 T 和 。 式（5-23）的阶跃响应标么化后，输出为 （5-25） 为求 T 和 ，在图 5-14 曲线上取两点：E 和 F ，且 ，则有 由此解出前面的图 5-13 则变为图 5-14。该标么化处理，并不改当然，为了计算上的方便，也可取 ， ，代入上两

13、式，有 ， 算出 T 和 后，可检验一下用式（5-25）与实测曲线的误差大小，如果误差可接受，则所求的传递函数式（5-23）可用。否则，应考虑用其它型传递函数（例如高阶传递函数）来描述。 具体方法为，另取三点： 、 和 ，具体为， ， 由式（5-25）算得， ， 并分别与图5-13中 、 、 对应的纵坐标比较即可。 （2）由阶跃响应曲线求二阶惯性及其以上环节的参数 当你用一阶惯性环节近似被控过程传递函数，检验发现误差不能满足原定的精度时，可考虑二阶及其以上的惯性环节传递函数。 当然，为了计算上的方便，也可取 设有阶跃响应曲线如图 5-15，现在，欲用二阶惯性环节的传递函数 （5-26） 来近似

14、描述它，其中， k 、 和 为待定的参数。 图 5-15 阶跃响应曲线图 5-16 具有纯时延的阶跃响应曲线 当输入为 时，该传递函数的响应为 （5-27） 设有阶跃响应曲线如图 5-15，现在，欲 在图 5-15 所示的阶跃响应曲线上，找出两点：A( ， ) 和 B( ， ),并将这两点分别代入式（5-27），有 其近似解为 研究表明，由式（5-27）表示的阶跃响应，应有 ，并且当 时，被控过程应为一阶惯性环节 ，且时间常数为 当 时，被控过程可为二阶等容惯性环节 ，且 在图 5-15 所示的阶跃响应曲线上，找出两当 时，被控过程应为二阶以上惯性环节，可用等 n 容惯性环节来描述 其中 n

15、和 T 分别按下式计算 如果算得 n 不为整数，应取最接近的整数。n 与 的关系也可用表 5-1 表示。 表5-1 多容过程的 n 与 之间的关系 n1234567891012140.3170.4600.5340.5840.6180.6400.6660.6840.6990.7120.7340.751当 时，被控过程应为二 如果阶跃响应曲线有明显的纯时延，如图 5-16 所示，则应在式（5-26）右边乘上一个纯时延环节： ，变为 其中 见图5-16。具体用上面的公式求 、 时，应在 和 中减去 时间段。 对于以上传递函数中的放大系数 K，仍可用式（5-24）求取。 （3）由无自平衡过程的阶跃响应

16、曲线求过程参数 当阶跃响应的曲线如图 5-17 所示时，该过程的传递函数则具有无自平衡特性，其特点是，随着 ，响应曲线的变化速率逐渐趋于某一常数。 图 5-17 无自平衡过程阶跃响应曲线 如果阶跃响应曲线有明显的纯时延，如图 5- 根据该响应曲线，该过程可用（5-28） 来近似。 当阶跃信号 作用于输入端时，其输出为 下面将讨论 T 和 的确定。作阶跃响应直线部分的延长线（见图中虚线段），与 t 轴相交于点 ，该线与时间轴 t 夹角为 ，于是 由于 ，所以 其中， 和 见图 5-17。 根据该响应曲线，该过程可用（5-28） 来近似。 由图5-17可知，用式（5-28）近似原过程的最大误差发生

17、在 这一段曲线上，即响应的起始段。为此可再加一个惯性环节来减小误差，即采用下列传递函数来描述过程 （5-29） 其中 T 的确定如上所示，即不变，而 和 的确定如下： 显然，用式（5-29）描述响应曲线表示的过程，比式（5-28）要精确些。 5.4 基于最小二乘法的过程辨识 应该说，用最小二乘法建模仍然是一种测试法建模，但含有较多的处理技巧与方法，这里将其单列为一节，主要是考虑内容稍多、篇幅较大。最小二乘法是系统辨识中的一种常用参数估计方法，它具有原理明了、算法简捷、收敛较快、相对容易理解的特点，因而被广泛用于参数估计之中。最小二乘法包括：最小二乘的批处理法、递推法、渐消记忆法和增广法等。 由

18、图5-17可知，用式（5-28）近似原过程5.4.1 离散时间系统模型 数学模型分为连续时间系统模型和离散时间系统模型。前面讨论的是连续时间系统模型，随着计算机的普及与应用，离散时间系统模型被越来越重视，最小二乘法采用的是离散时间系统模型。在这两类模型中按是否有随机扰动，每类又可分为确定性和随机性两种形式。 确定性离散系统（Deterministic Discrete Systems）的单输入/输出方程形式： （5-30） 式中 d 为纯时延，且 ， 和 分别为k时刻的输出和输入。 随机离散系统（Stochastic Discrete Systems）输入-输出差分模型一般有下面几种： 自回归

19、滑动平均（Auto-Regressive Moving Average，ARMA）模型 （5-31）式中 、 与式（5-30）中的相同， 5.4.1 离散时间系统模型 数学模型分为连 为白噪声（White Noise）序列，且 ， 式（5-31）右边第 1 项称为滑动平均项，左边项称为自回归项。式（5-31）也可写为 这里， 分别被称为过程模型和噪声模型，后者也被称为成形滤波器。而可以看作是白噪声经线性环节的输出，它一般为有色噪声（Coloured Noise）。 为白噪声（White Noise）序列，且 自回归积分滑动平均（Auto-Regressive Integrated Moving

20、 Average，ARIMA）模型或者这里， 、 、 和 与式（5-31）中相同。与式（5-31）相比，这里假定 。当 时， 多项式中前 项的系数为零。 最小二乘模型（Least Square Model） 或者式中， 、 、d 和 有与前面相同的含义。 滑动平均（Moving Average，MA）模型 自回归积分滑动平均（Auto-Reg与式（5-31）相比，这里有： 。 5.4.2 批处理最小二乘法 考虑最小二乘模型 式中， ， ， 为白噪声。 设已知 、 ，现在的任务是根据可量测的输入和输出，确定参数： ， 由输入-输出模型有 （5-32） 令： 为观测向量； 为待估参数向量，则式（5

21、-32）可写为另一形式（5-33） 与式（5-31）相比，这里有： 现有 N 次观测数据组并且 当 时，由式（5-32）有： 引入下列符号， 现有 N 次观测数据组并且 当 由式（5-33），有下列矩阵形式： 由于真实的参数向量 并不知道，不妨用 来表示它的估计值，于是，基于 的输出估计为 式中， 现在定义残差 （也是一随机变量）为实际输出与估计输出之差： （5-34） 对于 ，则有其中 。 现在的任务是：求使目标函数为最小的 （记为 ）。 由式（5-33），有下列矩阵形式： 由于真实的参数向量 展开上式，有由求极值的方法，对 求一阶导，并令其为零 从而有（5-35） 由于二阶导所以由式（5-

22、35）求得 的为极小值。式（5-35）即为批处理法的最小二乘估计。 最小二乘估计的统计特性讨论： 由于 为白噪声序列，故有 ， 展开上式，有由求极值的方法，对 求 （1）无偏性：使 J 为最小的参数估计向量 的数学期望为参数真值向量，即 这是因为它揭示最小二乘参数估计是围绕参数真值波动的统计性质。 （2）估计误差（偏差）协方差 主对角线上各元表现参数估计的散度，非对角线上各元反映参数估计 分量相互影响程度或相关性大小。 上式的成立，是因为 （1）无偏性：使 J 为最小的参数估计向量 （3）最小方差估计：设 为 的任一其他线性无偏估计，则 即最小二乘估计是最小方差估计，也就是参数估计 离参数真值

23、 最近。 由于 为 的任一线性无偏差估计，所以 可表示为 式中， 。由于， ，即 从而因为 （3）最小方差估计：设 为 所以 （4）一致收敛性：若 存在且正定，则 是一致收敛的，即 由于所以所以 （4）一致收敛性：若 存在且正定，则 由 的无偏性知 所以 例 5-1 现有最小二乘模型 其中 、 、 ， 为零均值白噪声，实验获得输入输出数据如表 5-2，试用批处理最小二乘法确定多项式 和 的参数。 解：首先组成矩阵 ，然后按式（5-31）求出 。下面用Matlab 语言来做。 u = -1; -1; -1; -1; -1; 1; 1; -1; -1; 1; -1; 1; 1; -1; 1; 1;

24、 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; -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;

25、 1.4520; -0.4325; -1.2545; 1.1510 ; % 将y 输入到工作空间 由 的无偏性知 所以 表5-2 实验获得的输入输出数据 tuytuytuy0.0005.5-1-1.965611.01-0.57070.5-10.05826.010.458711.5-10.50781.0-1-0.63956.511.371012.0-10.73941.5-1-1.85107.0-11.878312.5-1-1.43782.0-1-2.02427.510.245413.01-2.63282.5-1-1.43638.01-1.122313.5-1-0.53593.01-0.91888

26、.510.784814.0-11.45203.510.30539.012.498314.51-0.43254.0-12.29389.5-12.214715.01-1.25454.5-11.0769101-0.242415.511.15105.01-2.294310.5-1-1.5523-表5-2 实验获得的输入输出数据 tuytuytuy0.00 phi = 0; -1*y(1:end-1), 0; 0; -1*y(1:end-2), 0; u(1:end-1), 0; 0; u(1:end-2); % 组建 Theta = inv (phi*phi) * phi*y；% 计算 运行结果为：

27、= -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) = 0.6863 q-1 + 0.7947 q-2Estimated using ARX from da

28、ta set T Loss function 0.00146694 and FPE 0.00190159 phi = 0;Sampling interval: 0.5其中 即为 ，符号 q 就是本书中的符号 z。由此可见，两结果几乎相同。 5.4.3 递推最小二乘法 前面介绍的批处理法需要大量的计算机内存（随着N增大），而且当 变化时， 不能自动跟踪其变化，实时性不好。而递推最小二乘法（Recursive Least Square Method）则能解决这一问题。 引理：设A、C 和 均为非奇异方阵，则有 证明：用 左乘上式右端，若结果为单位矩阵，则说明其等式关系成立。 当 时，有 Sampl

29、ing interval: 0.5其中 （5-36） 递推最小二乘算法 最小二乘模型描述的系统 基于 及以前的输出 y 、 及以前的输入 u ，未知参数向量 的最小二乘估计 的递推公式为 （5-37） （5-38） （5-39） （5-36） 递推最小二乘算法 式中， 为观测向量， 为未知参数向量 的估计值， 为增益向量。 证明： 设 是基于 和 的估计，根据批处理最小二乘法有 对于 时的最小二乘估计为 式中 ， 于是式中， 为观测向量， 为未知参数向量 的估计值， （5-40） 令并视 、 和 分别为式（5-36）中的 A、B、D，再考虑 在此为标量，由引理有 再令 ，则有 又令（5-41）

30、 （5-40） 令并视 、 即为式（5-39），且式（5-41）可写为它就是式（5-39）。 将式（5-39）代入式（5-40），并注意交换标量和矩阵的位置，有即为式（5-39），且式（5-41）可写为它就是式（5-39即为式（5-37）。 递推最小二乘算法的说明： （1）从式（5-37）看，新的参数向量估计值 为先前的参数向量估计值 加修正项 在不断更新过程中， 、 和 的行列数不变，但它们的旧数据不断被新数据替换； （2） 为增益向量， 为误差的协方差阵，一般 与 成正比，协方差越大，说明估计值与真值相差越大，增益向量也会越大，所产生的校正作用也越大； （3）初值 和 的确定。 方法1：若

31、已有 组数据，则可批处理它们，并将结果作为初值，即 ， 方法2： ， ，其中 。 5.4.4 具有遗忘因子的递推最小二乘法 递推最小二乘法有一个缺点：常常出现“数据饱和”。随着k的增加， 和 变得越来越小，式（5-37）中的修正项对 的修正能力变得越来越弱，即新近加入的输入/输出数据对参数向量估计值的更新作用不大。这样导致的结果是：参数估计值难以接近真值；当参数真值时变时，该算法无法跟踪这种变化，从而使实时参数辨识失败。 即为式（5-37）。 递推最小二乘算法的说明： 解决该问题的方法之一是用具有遗忘因子（Forgetting Factor）的递推最小二乘法。 取性能指标函数： 式中， 为加权

32、对角阵： N为观测数据组数， 为遗忘因子： 。 按与前面相同的思路，可推出具有遗忘因子的递推最小二乘估计公式： 解决该问题的方法之一是用具有遗忘因子（Fo式中 。 几点说明： （1）当 时，该估计公式组即为递推最小二乘算法公式组； （2） 的选取范围一般在： ，参数变化快时， 取小点；变化慢时，取大些； （3）初值的选取与前面的递推最小二乘法相同。 5.4.5 递推增广最小二乘法 以上用的是最小二乘模型，在很多情况下， ，所以需考虑更一般的情况，即ARMA模型： （5-42） 其中， , , 式（5-42）可变为 式中 将其表示为 式中在 中，由于 不可测，所以只能用其估计值 来替换，设 式中

33、 用与前面类似的方法，并考虑用 代替 中的 ，则有下列递推公式 将其表示为 式中在 中，由于 不可测 的估计可用下列方法之一来进行： （1）预测形式： （2）滤波形式：并且当 时，有 。 前面介绍了最小二乘参数估计算法和特性，但是估计参数是否收敛到真实值尚未提及，由于这方面的内容涉及较多的定理及证明，这里只能给出结论，证明可参考有关资料。 被估参数收敛定理：若 是 N 阶持续激励（Persistently Exciting of Order N），即充分丰富（Sufficient Rich），则式（5-33）式（5-35）所示的算法能保证指数收敛。 的估计可用下列方法之一来进行： 例 5.2

34、现有自回归滑动平均模型 已知 ， 为白噪声：零均值，方差为 0.1，取采样周期 ，实验获得输入 / 输出数据如表 5-3，试用递推增广最小二乘法确定对象参数： 、 、 、 、 和 ，并将参数估计变化过程用图形表示出来。 表5-3 输入/输出数据u1.00001.00000.00000.0000-1.0000-1.00001.00001.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.928

35、11.69902.56812.41131.4639u1.00001.00000.00000.00001.00001.00000.00000.0000y-0.1749-1.9893-1.10591.08121.44970.93041.50192.6898u-1.0000-1.00001.00001.00000.00000.00001.00001.0000y2.35731.5460-0.1111-1.7658-0.90110.89141.47701.1226u0.00000.0000-1.0000-1.00001.00001.00000.00000.0000y1.68002.40182.30391

36、.4445-0.1399-1.9036-0.96571.1687 例 5.2 现有自回归滑动平均模型 已 解：由递推增广最小二乘法公式，有参数估计算法： (1) 初始化 、 等，并给 u 和 y 赋值； (2) ； (3) 计算 ； (4) 计算 ； (5) 保存参数 、 、 、 、 和 ； (6) 计算 ； (7) 求滤波形式的 ； (8) 移位处理 u、y、 和 ； (9) k 步终了吗? 否， ，去第 2 步；是，去第 10 步； (10) 打印结果。 按算法编写程序，经调试，得出 、 、 、 、 和 随拍数变化的规律，具体参见图5-18。具体程序为（用 Matlab 编写）： clea

37、rN=40; 解：由递推增广最小二乘法公式，有参数估计算 u=1,1,0,0,-1,-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,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

38、,2.6898,2.3573,1.5460,-0.1111, -1.7658,-0.9011, 0.8914,1.4770,1.1226,1.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; %构成 K=P*f/(1+f*P*f); Q=Q0+K*(y(j)-f*Q0); % a1(j)=Q(1);a2(j)=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; %求 u=1,1,0,0,-1,-1,1, y2=y1;y1=y(j);u3=u2;u2=u1;u1=u(j);x2=x1