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

5.3 测试法建模测试法建模5.3.1 测试法建模的方法测试法建模的方法 机理法建模能够解决实际生产中一部分过程的建模问题。但是，还有很多生产过程由于工艺的复杂性、产品本身在加工中的变化性（如物理或化学变化），使得用机理法建模在技术上遇到了极大的困难，人们不得不考虑用其它的方法建模。从学科角度看，建立数学模型应该属于系统辨识（System Identification）与参数估计（Parametric Estimation）的范畴。简单地说，系统辨识主要是对被研究对象的结构进行判断，解决“是什么”的问题，例如一阶惯性环节，二阶系统；而参数估计则对支撑结构的参数进行估计，解决“是多少”的问题。事实上，有很多比较复杂的过程，我们对其工作机理并不清楚，更难以用数学和物理的方法加以具体描述，此时用测试法建模是一个不得已而为之的方法。与前述的机理法建模相比，测试法建模不需要深入了解过程的工作机制，通常的做法是将其看作一个“黑箱”，通过从外部施加适当的输入信号，测得过程的输出信号，通过对这些输入和输出信号的处理和研究，获得其动态特性和数学模型。因此，问题主要归纳为：1）施加何种输入信号才能最大限度地激励被测过程，使得动态特性得以充分表现，并通过输出信号显露出来？2）对获得的数据或波形，通过什么方法和技术才能估算出适用于控制用的动态模型？一般说来，模型有非参数模型（Nonparametric Model）和参数模型（Parametric Model）之分。建立非参数模型的方法通常有：时域法（Time-domain Method）、频域法（Frequency-5.3 测试法建模5.3.1 测试法建模的方法 1-domain Method）和统计相关法（Statistical Correlation Method）等，这类建模不需要事先确定模型的结构，可用于广泛的被控过程；获得参数模型的方法主要有：最小二乘法（Least Square Method）、极大似然法（Maximum Likelihood Method）和梯度校正法（Gradient Correction Method）等，这类建模需要假设一定的模型结构，通过极小化模型与过程之间的误差准则，来确定相应的模型参数。用时域法测定被控过程的数学模型：对过程施加阶跃信号，或者方波信号，测取响应曲线，并由此确定过程的传递函数。该方法具有测试简单、需用设备少的优点，但测试精度不高，其获得的模型可用于一般工业过程控制；用频域法测定被控过程的数学模型：对过程施加不同频率的正弦波输入信号，获得相应的输出幅值与相位，由此可得到该过程的频率特性，由频率特性获得传递函数。该方法需用专门的频率发生和测试设备，模型精度比用时域法高；用统计相关法测定被控过程的数学模型：对被控过程施加伪随机信号，采用统计相关法获得过程的动态特性。它的特点是，可在生产状态下施加随机信号，并测取相关数据，精度较高，但需获得较多数据，并借助计算机协助处理。最小二乘法又称最小平方法，是估计离散时间数学模型参数的一常用种方法。随着计算机技术在控制中的应用，最小二乘法在过程辨识的实践中被越来越广泛地采用。本章主要讨论用时域响应法和最小二乘法获取数学模型。5.3.2 时域响应曲线法时域响应曲线法 1响应曲线的测取-domain Method）和统计相关法（Statisti2时域响应曲线法是对被控过程施加阶跃信号，如果被控过程不允许长期施加阶跃信号，则改用矩形脉冲信号，然后测取响应曲线，并由此求取输入和输出之间的传递函数。（1）阶跃响应曲线的测取 当被控过程稳定之后，对调节阀施加一个幅值合适的阶跃信号，用记录仪或数据采集系统记录被控量的变化曲线，例如被控量为温度时，就记录温度响应曲线，直到变化曲线进入稳定状态为止。下面几点是在这一过程中值得注意的：a)在施加阶跃信号之前，被控过程应处在较为稳定的工作状态，在下一轮施加输入信号时，应等前一过程结束、并恢复稳态一段时间之后进行。b)施加阶跃信号的幅度：通常为正常输入信号的 515，以不影响正在进行的生产为好。同时，幅度也不能太小，因为过小的输入容易被其它信号淹没，在响应曲线上难以表现出来。c)多次、全面测试，消除偶然性，获得真实结果：试验不仅应在相同条件下，重复几次，以获得两次及其以上的较为接近的响应曲线，而且也应选取不同负荷、不同输入值，测得相应的响应曲线，以获得全面的动态特性。（2）矩形脉冲相应曲线的测取 在有些情况下，用阶跃信号输入时，可能会危及安全生产，或者影响产品的质量和数量。此时，输入可考虑用矩形脉冲信号代替阶跃信号，测得过程的矩形脉冲响应曲线。由于通过阶跃响应曲线求传递函数被人们所熟悉，所以，往往将矩形脉冲响应曲线再转化为阶跃响应曲线，进而按阶跃响应曲线法确定传递函数。图 5-12a)为矩形脉冲输入信号，它可以分解为图 b)所示的两个阶跃信号的叠加，即时域响应曲线法是对被控过程施加阶跃信号，如果被控过程不允许长3其中 a 为脉冲宽度，。如果被控过程是线性的，则其矩形脉冲响应曲线 可分解为阶跃响应曲线 和 ，即图 5-12 矩形脉冲及其响应分解图 见图5-12 c)所示。这里，和 分别为 和 的响应。于是，起源于零点的阶跃响应为：（5-22）其中，为矩形脉冲响应，为起源于 点的阶跃响应（注意此时的符号为正）。式（5-22）为由矩形脉冲响应求阶跃响应的公式，可用作图法逐步求出：t 在 0a 时，当 时，此时的 前面已经有，为已知，所以可求得 ，如此类推，一直进行到进入稳态。其中 a 为脉冲宽度，4 由无自平衡能力过程的矩形脉冲响应曲线转化为阶跃响应曲线，也可通过作图法实现，这可作为练习，留给读者来完成（见思考题与习题 5-5）。2由过程阶跃响应曲线确定传递函数 由阶跃响应曲线确定传递函数，通常需要确定传递函数的结构及其参数两部分。结构形式是指被控过程的传递函数形式，生产过程主要是：一阶惯性环节、二阶惯性环节、或等 n 阶惯性环节，并且这些环节时常含有纯滞后，其表达形式为，对于无自平衡能力的过程，也有类似的形式，传递函数的参数是伴随结构形式出现的待定常数，如一阶惯性、具有纯时延、有自平衡能力的传递函数含有：K、T 和 三个需要确定的参数。由无自平衡能力过程的矩形脉冲响应曲线转化为阶5 关于传递函数结构形式的确定，主要有两方面的考虑：一是根据对被控过程的经验和知识（即通常所说的先验知识）来确定；二是根据控制的要求，尽量将一个原本较复杂的过程用低阶的传递函数来近似描述，因此产生的误差只要处在可接受的范围即可。下面的讨论，集中在参数的确定上。（1）由阶跃响应曲线求一阶惯性加纯时延环节的参数 这里传递函数形式为（5-23）它有三个参数需要确定，即放大系数 K、时间常数 T 和纯时延 。原本就是一阶惯性加纯时延过程的阶跃响应曲线见图 5-5，其 T 和 从图中很容易确定，放大系数为 （其中 为阶跃输入信号的幅值）。下面讨论原本是二阶及其以上过程，且响应曲线呈“s”形，如何用式（5-23）来近似描述。现有阶跃响应曲线如图 5-13 所示，试图用式（5-23）来近似描述，需确定 K、T 和 三个参数。显然，放大系数可用下式求得（5-24）关于传递函数结构形式的确定，主要有两方面的6 图 5-13 由阶跃响应曲线确定 T 和 图5-14 纵坐标的标么化其中 为输入幅值，为已知量。关于 T 和 的确定，有两种方法：一是作图法，二是计算法，下面分别介绍。用作图法求 T 和 ：首先找到响应曲线上凹和下凹的交接点-拐点 D，过 D 点作曲线的切线，切线与时间轴 t 相交于 A 点，与 相较于 C 点，该点在时间轴上的投影为 B，则OA为 ，AB 为 T。作图法的问题是：曲线的拐点有时不容易找到，并且作切线时，有一定的随意性。所以，用作图法求 T 和 ，可能会因人而异，有一定的误差。计算法求 T 和 ：将阶跃响应曲线的纵坐标 标么化，即：实际值基准值，这里取基准值为 ，于是 的标么值（Unit Value）为 图 5-13 由阶跃响应曲线确定 T 和 7前面的图 5-13 则变为图 5-14。该标么化处理，并不改变响应曲线的横坐标和形状，仅方便求得参数。我们的目的是要通过图5-14所示的阶跃响应曲线，容易求出式（5-23）中的 T 和 。式（5-23）的阶跃响应标么化后，输出为（5-25）为求 T 和 ，在图 5-14 曲线上取两点：E 和 F ，且 ，则有 由此解出前面的图 5-13 则变为图 5-14。该标么化处理，并不改8当然，为了计算上的方便，也可取 ，代入上两式，有，算出 T 和 后，可检验一下用式（5-25）与实测曲线的误差大小，如果误差可接受，则所求的传递函数式（5-23）可用。否则，应考虑用其它型传递函数（例如高阶传递函数）来描述。具体方法为，另取三点：、和 ，具体为，由式（5-25）算得，并分别与图5-13中 、对应的纵坐标比较即可。（2）由阶跃响应曲线求二阶惯性及其以上环节的参数 当你用一阶惯性环节近似被控过程传递函数，检验发现误差不能满足原定的精度时，可考虑二阶及其以上的惯性环节传递函数。当然，为了计算上的方便，也可取 9 设有阶跃响应曲线如图 5-15，现在，欲用二阶惯性环节的传递函数（5-26）来近似描述它，其中，k 、和 为待定的参数。图 5-15 阶跃响应曲线图 5-16 具有纯时延的阶跃响应曲线 当输入为 时，该传递函数的响应为（5-27）设有阶跃响应曲线如图 5-15，现在，欲10 在图 5-15 所示的阶跃响应曲线上，找出两点：A(，)和 B(，),并将这两点分别代入式（5-27），有 其近似解为 研究表明，由式（5-27）表示的阶跃响应，应有 ，并且当 时，被控过程应为一阶惯性环节 ，且时间常数为 当 时，被控过程可为二阶等容惯性环节 ，且 在图 5-15 所示的阶跃响应曲线上，找出两11当 时，被控过程应为二阶以上惯性环节，可用等 n 容惯性环节来描述 其中 n 和 T 分别按下式计算 如果算得 n 不为整数，应取最接近的整数。n 与 的关系也可用表 5-1 表示。表表5-1 多容多容过程的程的 n 与与 之间的关系之间的关系 n1234567891012140.3170.4600.5340.5840.6180.6400.6660.6840.6990.7120.7340.751当 时，被控过程应为二12 如果阶跃响应曲线有明显的纯时延，如图 5-16 所示，则应在式（5-26）右边乘上一个纯时延环节：，变为 其中 见图5-16。具体用上面的公式求 、时，应在 和 中减去 时间段。对于以上传递函数中的放大系数 K，仍可用式（5-24）求取。（3）由无自平衡过程的阶跃响应曲线求过程参数 当阶跃响应的曲线如图 5-17 所示时，该过程的传递函数则具有无自平衡特性，其特点是，随着 ，响应曲线的变化速率逐渐趋于某一常数。图 5-17 无自平衡过程阶跃响应曲线 如果阶跃响应曲线有明显的纯时延，如图 5-13 根据该响应曲线，该过程可用（5-28）来近似。当阶跃信号 作用于输入端时，其输出为 下面将讨论 T 和 的确定。作阶跃响应直线部分的延长线（见图中虚线段），与 t 轴相交于点 ，该线与时间轴 t 夹角为 ，于是 由于 ，所以 其中，和 见图 5-17。根据该响应曲线，该过程可用（5-28）来近似。14 由图5-17可知，用式（5-28）近似原过程的最大误差发生在 这一段曲线上，即响应的起始段。为此可再加一个惯性环节来减小误差，即采用下列传递函数来描述过程（5-29）其中 T 的确定如上所示，即不变，而 和 的确定如下：显然，用式（5-29）描述响应曲线表示的过程，比式（5-28）要精确些。5.4 基于最小二乘法的过程辨识基于最小二乘法的过程辨识 应该说，用最小二乘法建模仍然是一种测试法建模，但含有较多的处理技巧与方法，这里将其单列为一节，主要是考虑内容稍多、篇幅较大。最小二乘法是系统辨识中的一种常用参数估计方法，它具有原理明了、算法简捷、收敛较快、相对容易理解的特点，因而被广泛用于参数估计之中。最小二乘法包括：最小二乘的批处理法、递推法、渐消记忆法和增广法等。由图5-17可知，用式（5-28）近似原过程155.4.1 离散时间系统模型离散时间系统模型 数学模型分为连续时间系统模型和离散时间系统模型。前面讨论的是连续时间系统模型，随着计算机的普及与应用，离散时间系统模型被越来越重视，最小二乘法采用的是离散时间系统模型。在这两类模型中按是否有随机扰动，每类又可分为确定性和随机性两种形式。确定性离散系统（Deterministic Discrete Systems）的单输入/输出方程形式：（5-30）式中 d 为纯时延，且 ，和 分别为k时刻的输出和输入。随机离散系统（Stochastic Discrete Systems）输入-输出差分模型一般有下面几种：自回归滑动平均（Auto-Regressive Moving Average，ARMA）模型 （5-31）式中 、与式（5-30）中的相同，5.4.1 离散时间系统模型 数学模型分为连16 为白噪声（White Noise）序列，且，式（5-31）右边第 1 项称为滑动平均项，左边项称为自回归项。式（5-31）也可写为 这里，分别被称为过程模型和噪声模型，后者也被称为成形滤波器。而可以看作是白噪声经线性环节的输出，它一般为有色噪声（Coloured Noise）。为白噪声（White Noise）序列，且17 自回归积分滑动平均（Auto-Regressive Integrated Moving Average，ARIMA）模型或者这里，、和 与式（5-31）中相同。与式（5-31）相比，这里假定 。当 时，多项式中前 项的系数为零。最小二乘模型（Least Square Model）或者式中，、d 和 有与前面相同的含义。滑动平均（Moving Average，MA）模型 自回归积分滑动平均（Auto-Reg18与式（5-31）相比，这里有：。5.4.2 批处理最小二乘法批处理最小二乘法 考虑最小二乘模型 式中，为白噪声。设已知 、，现在的任务是根据可量测的输入和输出，确定参数：，由输入-输出模型有（5-32）令：为观测向量；为待估参数向量，则式（5-32）可写为另一形式（5-33）与式（5-31）相比，这里有：19 现有 N 次观测数据组并且 当 时，由式（5-32）有：引入下列符号，现有 N 次观测数据组并且 当 20由式（5-33），有下列矩阵形式：由于真实的参数向量 并不知道，不妨用 来表示它的估计值，于是，基于 的输出估计为 式中，现在定义残差 （也是一随机变量）为实际输出与估计输出之差：（5-34）对于 ，则有其中 。现在的任务是：求使目标函数为最小的 （记为 ）。由式（5-33），有下列矩阵形式：由于真实的参数向量 21 展开上式，有由求极值的方法，对 求一阶导，并令其为零 从而有（5-35）由于二阶导所以由式（5-35）求得 的为极小值。式（5-35）即为批处理法的最小二乘估计。最小二乘估计的统计特性讨论：由于 为白噪声序列，故有，展开上式，有由求极值的方法，对 求22 （1）无偏性：使 J 为最小的参数估计向量 的数学期望为参数真值向量，即 这是因为它揭示最小二乘参数估计是围绕参数真值波动的统计性质。（2）估计误差（偏差）协方差 主对角线上各元表现参数估计的散度，非对角线上各元反映参数估计 分量相互影响程度或相关性大小。上式的成立，是因为 （1）无偏性：使 J 为最小的参数估计向量 23 （3）最小方差估计：设 为 的任一其他线性无偏估计，则 即最小二乘估计是最小方差估计，也就是参数估计 离参数真值 最近。由于 为 的任一线性无偏差估计，所以 可表示为 式中，。由于，即 从而因为 （3）最小方差估计：设 为 24所以 （4）一致收敛性：若 存在且正定，则 是一致收敛的，即 由于所以所以 （4）一致收敛性：若 存在且正定，则 25 由 的无偏性知 所以 例例 5-1 现有最小二乘模型 其中 、，为零均值白噪声，实验获得输入输出数据如表 5-2，试用批处理最小二乘法确定多项式 和 的参数。解解：首先组成矩阵 ，然后按式（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;-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 输入到工作空间 由 的无偏性知 所以 26表表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-10.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-表5-2 实验获得的输入输出数据 tuytuytuy0.0027 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；%计算 运行结果为：=-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 data set T Loss function 0.00146694 and FPE 0.00190159 phi=0;28Sampling interval:0.5其中 即为 ，符号 q 就是本书中的符号 z。由此可见，两结果几乎相同。5.4.3 递推最小二乘法递推最小二乘法 前面介绍的批处理法需要大量的计算机内存（随着N增大），而且当 变化时，不能自动跟踪其变化，实时性不好。而递推最小二乘法（Recursive Least Square Method）则能解决这一问题。引理：设A、C 和 均为非奇异方阵，则有 证明证明：用 左乘上式右端，若结果为单位矩阵，则说明其等式关系成立。当 时，有 Sampling interval:0.5其中 29（5-36）递推最小二乘算法递推最小二乘算法 最小二乘模型描述的系统 基于 及以前的输出 y、及以前的输入 u，未知参数向量 的最小二乘估计 的递推公式为（5-37）（5-38）（5-39）（5-36）递推最小二乘算法 30式中，为观测向量，为未知参数向量 的估计值，为增益向量。证明：证明：设 是基于 和 的估计，根据批处理最小二乘法有 对于 时的最小二乘估计为 式中，于是式中，为观测向量，为未知参数向量 的估计值，31（5-40）令并视 、和 分别为式（5-36）中的 A、B、D，再考虑 在此为标量，由引理有 再令 ，则有 又令（5-41）（5-40）令并视 、32即为式（5-39），且式（5-41）可写为它就是式（5-39）。将式（5-39）代入式（5-40），并注意交换标量和矩阵的位置，有即为式（5-39），且式（5-41）可写为它就是式（5-3933即为式（5-37）。递推最小二乘算法的说明：（1）从式（5-37）看，新的参数向量估计值 为先前的参数向量估计值 加修正项 在不断更新过程中，、和 的行列数不变，但它们的旧数据不断被新数据替换；（2）为增益向量，为误差的协方差阵，一般 与 成正比，协方差越大，说明估计值与真值相差越大，增益向量也会越大，所产生的校正作用也越大；（3）初值 和 的确定。方法1：若已有 组数据，则可批处理它们，并将结果作为初值，即，方法2：，其中 。5.4.4 具有遗忘因子的递推最小二乘法具有遗忘因子的递推最小二乘法 递推最小二乘法有一个缺点：常常出现“数据饱和”。随着k的增加，和 变得越来越小，式（5-37）中的修正项对 的修正能力变得越来越弱，即新近加入的输入/输出数据对参数向量估计值的更新作用不大。这样导致的结果是：参数估计值难以接近真值；当参数真值时变时，该算法无法跟踪这种变化，从而使实时参数辨识失败。即为式（5-37）。递推最小二乘算法的说明：34 解决该问题的方法之一是用具有遗忘因子（Forgetting Factor）的递推最小二乘法。取性能指标函数：式中，为加权对角阵：N为观测数据组数，为遗忘因子：。按与前面相同的思路，可推出具有遗忘因子的递推最小二乘估计公式：解决该问题的方法之一是用具有遗忘因子（Fo35式中 。几点说明：（1）当 时，该估计公式组即为递推最小二乘算法公式组；（2）的选取范围一般在：，参数变化快时，取小点；变化慢时，取大些；（3）初值的选取与前面的递推最小二乘法相同。5.4.5 递推增广最小二乘法递推增广最小二乘法 以上用的是最小二乘模型，在很多情况下，所以需考虑更一般的情况，即ARMA模型：（5-42）其中，,式（5-42）可变为 式中 36将其表示为 式中在 中，由于 不可测，所以只能用其估计值 来替换，设 式中 用与前面类似的方法，并考虑用 代替 中的 ，则有下列递推公式 将其表示为 式中在 中，由于 不可测37 的估计可用下列方法之一来进行：（1）预测形式：（2）滤波形式：并且当 时，有 。前面介绍了最小二乘参数估计算法和特性，但是估计参数是否收敛到真实值尚未提及，由于这方面的内容涉及较多的定理及证明，这里只能给出结论，证明可参考有关资料。被估参数收敛定理被估参数收敛定理：若 是 N 阶持续激励（Persistently Exciting of Order N），即充分丰富（Sufficient Rich），则式（5-33）式（5-35）所示的算法能保证指数收敛。的估计可用下列方法之一来进行：38 例例 5.2 现有自回归滑动平均模型 已知 ，为白噪声：零均值，方差为 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.92811.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.4445-0.1399-1.9036-0.96571.1687 例 5.2 现有自回归滑动平均模型 已39 解解：由递推增广最小二乘法公式，有参数估计算法：(1)初始化 、等，并给 u 和 y 赋值；(2)；(3)计算 ；(4)计算 ；(5)保存参数 、和 ；(6)计算 ；(7)求滤波形式的 ；(8)移位处理 u、y、和 ；(9)k 步终了吗?否，去第 2 步；是，去第 10 步；(10)打印结果。按算法编写程序，经调试，得出 、和 随拍数变化的规律，具体参见图5-18。具体程序为（用 Matlab 编写）：clearN=40;解：由递推增广最小二乘法公式，有参数估计算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,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,41 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 相差不大，估计成功。该题也可用 Matlab 的辨识工具获得结果。在将输入输出数据键入Matlab工作空间后，T=armax(y,u,2,2,2,2)Discrete-time IDPOLY model:A(q)y(t)=B(q)u(t)+C(q)e(t)A(q)=1-0.7956 q-1+0.1438 q-2 B(q)=1.006 q-2+0.4969 q-3 C(q)=1-0.736 q-1+0.1019 q-2 Estimated using ARMAX Loss function 0.00841826 and FPE 0.0121597 Sampling interval:1 y2=y1;y1=42 这里符号 q 就是本书的符号 z，q-i 即为 。图 5-18 参数 、和 随拍数变化情况 5.4.6 模型阶次的确定模型阶次的确定 模型的结构在此表现为模型阶次，前面用最小二乘法估计参数时，我们假定差分方程的阶数，即模型阶次是已知的。事实上，现实情况往往并非如此，有一些实际问题的模型阶次可以用理论推导的办法获得，但更多的实际问题模型阶次是需要辨识的。辨识模型阶次的方法较多，如损失函数法（也称拟合度检验法）、Hankel矩阵法、行列式比法、最终预报误差法等。其中较为简单、用得较多的是损失函数法，因为它可以和参数估计 这里符号 q 就是本书的符号 z，q43一起求，具体方法如下：（1）依次设定模型的阶次 ，计算相应的最小二乘参数估计值 ，及其损失函数 （2）随着 n 的增加，将有明显的减小（在不同噪声水平下）。（3）当进行到 时，若 比 并无明显的减小，即 （为某一较小值），则 n 即为模型的阶次。例如，用最小二乘法估计参数时，在某噪声水平下，算得损失函数：，,，则有：，即阶次为 2。关于纯时延时间 ，离散系统中，通常取采样时间间隔 T 的整数倍 d，即 ，、2、3、。通常是可以事先知道的。如果由于工艺的复杂性，使一些变量无法测量，从而使纯时延未知时，可用阶跃响应曲线实验法确定，也可和阶次 n 一起进行估计，最佳的 d 应该是使目标函数（也称损失函数）J 趋于最小的值。具体算法如下：1设定模型阶次 ；2设定纯时延 ；3应用最小二乘法估计 ；4计算残差 E 和误差函数 J；一起求，具体方法如下：（1）依次设定模型的阶次 44 5.求最小误差函数 J；6d 为最佳值吗？是，继续下步，否则，返回第3步；7判断阶次 n 合适吗？是，继续下步，否则，返回第3步；8给出结果：d 和 n。5.求最小误差函数 J；45