系统建模与仿真 第5次课 第三章课件

则式（2.87）可写为

（2.88）

N维输出向量2n+1维参数向量N维噪声向量N×(2n+1)维测量矩阵最小二乘法：回顾1则式（2.87）可写为（2.88）N维输出向量2n+1维可得的最小二乘估计

（2.98）

J为极小值的充分条件是

（2.99）

即矩阵为正定矩阵，或者说矩阵是非奇异的。

2可得的最小二乘估计（2.98）J为极小值的充分条件

最小二乘估计的概率性质

如果ξ(k)是不相关随机数序列，且均值为0。1）无偏性2）一致性3）渐进正态性性辅助变量法、广义最小二乘法如果ξ是均值为0且服从正态分布的白噪声向量，则最小二乘参数估计值服从正态分布。3最小二乘估计的概率性质如果ξ(k)是不相关随机数序列，且2.8辅助变量法

现在开始讨论如何克服最小二乘法的有偏估计问题。

对于原辨识方程

（2.157）当是不相关随机序列时，最小二乘法可以得到参数向量的一致性无偏估计。但是，在实际应用中往往是相关随机序列。

42.8辅助变量法现在开始讨论如何克服最小二乘法的式中Q是非奇异的。假定存在着一个的矩阵Z（与同阶数）,满足约束条件

（2.158）用乘以式（2.157）等号两边得（2.159）由上式可得

（2.160）5式中Q是非奇异的。假定存在着一个如果取

（2.161）作为的估值，则称估值为辅助变量估值，矩阵Z称为辅助变量矩阵，Z中的元素称为辅助变量。

从式（2.161）可以看到，与最小二乘法估值的计算公式（2.98）具有相同的形式，因此计算比较简单。

根据式（2.160）和式（2.161）可得

（2.162）6如果取（2.161）作为的估值，则称估值为辅（2.163）当N很大时，对上式等号两边取极限得

（2.164）根据式（2.157）所假定的约束条件，可得

因此辅助变量估计是无偏估计。

7（2.163）当N很大时，对上式等号两边取极限得（2.16剩下的问题是如何选择辅助变量，即如何确定辅助变量矩阵Z的各个元素。选择辅助变量的基本原则是式（2.158）所给出的两个条件必须得到满足。

这可以简单地理解为所选择的辅助变量应与不相关，但与和中的强烈相关。8剩下的问题是如何选择辅助变量，即如何确定辅助变量矩阵Z可以有各种选择方法，下面介绍几种常用的选择方法。1）递推辅助变量参数估计法

辅助变量取作,是辅助模型

（2.165）的输出向量的元素，辅助变量矩阵Z为

9Z可以有各种选择方法，下面介绍几种常用的选择方法。1）递推辅（2.166）（2.167）10（2.166）（2.167）10递推辅助变量参数估计法计算步骤：11递推辅助变量参数估计法计算步骤：112）自适应滤波法

这种方法所选择的辅助变量和辅助变量矩阵Z的形式与上一种方法完全相同，只是辅助模型中参数向量的估计方法与上一种方法有所不同。取

（2.168）式中：取；d取；为k时刻所得到的参数向量估计值。当是持续激励信号时，所选的辅助变量可以满足式（2.158）所给出的2个约束条件。

122）自适应滤波法这种方法所选择的辅助变量和2.9最小二乘类辨识算法的比较1.基本最小二乘法（LS）（1）对白噪声，参数估计是无偏的一致的最小方差估计。（2）对有色噪声，参数估计有偏，但具有收敛性。（3）对未知的直流分量敏感。（4）对高阶系统，优于其它算法，可作为其它辨识算法的较好的起始算法。（5）一次完成算法精度较高，但逆矩阵计算量大，不适于在线辨识。（6）数据较多时，占用存储量和机时较多。132.9最小二乘类辨识算法的比较1.基本最小二乘法（LS）（2.递推最小二乘法（RLS）（1）对白噪声，参数估计是无偏的一致的最小方差估计。（2）对有色噪声，参数估计有偏，但具有收敛性。（3）对未知的直流分量敏感。（4）对高阶系统，优于其它算法，可作为其它辨识算法的较好的起始算法。（5）计算量少，占用内存较少。（6）适用于在线辨识，采用自适应算法后，可用于时变系统的辨识。142.递推最小二乘法（RLS）（1）对白噪声，参数估计是无偏3.辅助变量法（IV）（1）估计值是弱一致估计的。（2）对初态敏感，初态选取得不合适，就可能不收敛。（3）算法简单，计算量少，但不能得到噪声模型。153.辅助变量法（IV）（1）估计值是弱一致估计的。（2）2.10系统结构辨识模型阶的确定

在一些实际问题中，模型的阶可以按理论推导获得，而在另一些实际问题中，模型的阶却无法用理论推导的方法确定，需要对模型的阶进行辨识。下面介绍几种常用的模型阶的确定方法。

2.10.1按残差方差定阶

一种简单而有效的方法就是选定模型阶数n的不同取值，按估计误差方差最小或F检验来确定模型的阶。1)按估计误差方差最小定阶

考虑系统模型

162.10系统结构辨识模型阶的确定在一些实际问题（2.178）式中：为输出；为输入。

设是均值为0、方差为的白噪声序列。用最小二乘法求出的估值。具体为

（2.179）（2.180）（2.181）17（2.178）式中：为输出；为输入。(2.182)(2.183)残差为

18(2.182)(2.183)残差为18如图2.19所示，对某一系统，当=1，2，…时，随着的增加而减小。如果为正确的阶，则在时，出现最后一次陡峭的下降，再增大，则保持不变或只有微小的变化。

图2.19所示的例子，=3。

19如图2.19所示，对某一系统，当=1，2，…时，图2.19曲线图

20图2.19曲线图202)确定模型阶的F检验法

由于随着的增加而减小，在阶数的增大过程中，我们对那个使显著减小的阶感兴趣。为此，引入准则

（2.187）式中表示具有N对输入和输出数据、有个模型参数的系统估计误差的平方和。212)确定模型阶的F检验法由于随着的增表2.1某一系统计算结果

0.993.159.439.6750.94416.56418.73426.40447.25469.64592.65654321计算时取，，。从表2.1可以看出：当时，t的减小是显著的；当时，t的减小是不显著的。所以该系统的阶数可选为3。

对某一系统的计算结果如表2.1所列。

22表2.1某一系统计算结果0.993.159.4392.10.2确定阶的Akaike信息准则

与上述2个准则不同，Akaike信息准则（AIC，AkaikeInformationCriterion）是一个考虑了模型复杂性的准则。这个准则定义为

（2.189）式中：L是模型的似然函数；是模型中的参数数目。当AIC为最小的那个模型就是最佳模型。这个准则是Akaike总结了时间序列统计建模的发展史，在企图对一个复杂系统寻找近似模型的概率论的大量探索启示下，借助信息论而提出的一个合理的确定阶的准则。在一组可供选择的随机模型中，AIC最小的那个模型是一个可取的模型。这个准则的优点就在于它是一个完全客观的准则，应用这个准则时，不要求建模人员主观地判断“陡峭的下降”。

232.10.2确定阶的Akaike信息准则与上述21)白噪声情况下的AIC定阶公式

假定是均值为0、方差为并且服从正态分布的不相关随机噪声。

考虑系统模型（2.190）式中

根据前面几节的讨论和定义，由式（2.190）可写出关系式

（2.191）241)白噪声情况下的AIC定阶公式假定是均值为0、方输出变量在条件下的似然函数为

（2.192）对上式取对数可得

（2.193）求使为最大的的估值。根据可得

与前述的最小二乘估计一致。按照可得

25输出变量在条件下的似然函数为（2.192）对上式取因此

即（2.194）式中c为一常数。

式（2.194）给出了AIC定义中的第一项。

为

待估的参数，，…，，，，…，及，共有个，即，因而有26因此即（2.194）式中c为一常数。式（2.19即

可去掉上式中的常数项，则（2.195）选取不同的阶数和，按式（2.195）计算AIC，可得最优阶数和。在式（2.195）中加进项，表示对不同，若相近时，则取较小的模型。

27即可去掉上式中的常数项，则（2.195）选取不同的例2.4系统模型为

式中：是均值为0、方差为1且服从正态分布的不相关随机噪声；输入信号采用伪随机数。辨识模型采用的形式为

数据长度取N=1024。为了避免非平稳过程的影响，去掉前300个数据，取=1，2，3，4，=1，2，3，4，分别计算AIC（，），即28例2.4系统模型为式中：是均值为0、方差为1且服计算结果如表2.2所列。显然，应取=3，=2，可见利用AIC确定的模型阶次与系统的真实阶次相同。

表2.2不同和所对应的AIC16.218

15.108

15.931

417.649

15.599

14.070

25.864

316.800

30.393

51.085

280.046

223.380

97.353

341.766

1022.94

1432129计算结果如表2.2所列。显然，应取=3，=2，可2)有色噪声情况下的AIC定阶公式

有色噪声情况下的系统模型可以表示为

（2.196）式中

是均值为0、方差为，且服从正态分布的不相关随机噪声。

与前面的讨论相类似，可得

（2.197）302)有色噪声情况下的AIC定阶公式有色噪声情况下的服从正态分布N（0，1），为二位式伪随机序列，数据长度N=300。假定模型阶，利用极大似然法估计模型参数，根据AIC的计算结果，可得：当时，AIC最小，其结果与所给系例2.5：设系统模型为

式中

统模型相符合。31服从正态分布N（0，1），为二位式伪随机序列，F检验法和AIC准则法，可以在参数估计的过程中同时辨识阶次，所以应用上比较方便。阶次辨识和参数辨识两者是互相依赖的，也就是说进行参数估计时需要已知阶次；而辨识阶次时又要利用参数估计值，二者是不可分离的。32F检验法和AIC准则法，可以在参数估计的过程中同时辨识阶次，第三章连续系统的数字仿真本章讨论经典的连续系统数字仿真的原理与方法，内容包括连续系统数字仿真的基本概念、经典的数值积分法、经典的线性多步法等。3.1节讨论离散化原理及要求，是连续系统仿真的基础；3.2节对经典的数值积分法——龙格—库塔法及其他典型的数值积分法仿真建模原理进行详细分析，并通过实例说明其应用要点；3.3节对经典的线性多步法进行了介绍。

33第三章连续系统的数字仿真本章讨论经典的连续系统3.1离散化原理及要求

在数字计算机上对连续系统进行仿真时，首先遇到的问题是，数字计算机的数值及时间均具有离散性，而被仿真系统的数值及时间均具有连续性，后者如何用前者来实现。

连续系统仿真，从本质上是从时间、数值两个方面对原系统进行离散化，并选择合适的数值计算方法来近似积分运算，由此得到离散模型来近似原连续模型。如何保证离散模型的计算结果从原理上确能代表原系统的行为，这是连续系统数字仿真首先必须解决的问题。

343.1离散化原理及要求在数字计算机上对连续系统图3.1相似原理示意图

35图3.1相似原理示意图35设系统模型(参见图3.1)为，其中为输入变量，为系统变量。令仿真时间间隔为h，离散化后的输入变量为，系统变量为，其中表示。如果，，即，

(对所有k＝0，1，2…)，则可认为两模型等价，这称为相似原理。

36设系统模型(参见图3.1)为，实际上，要完全保证，是很困难的。进一步分析离散化引入的误差，随着计算机技术的发展，由计算机字长引入的舍入误差可以忽略，关键是数值积分算法，也称为仿真建模方法。相似原理用于仿真时，对仿真建模方法有三个基本要求：

(1)稳定性

若原连续系统是稳定的，则离散化后得到的仿真模型也应是稳定的。(2)准确性

有不同的准确性评价准则，最基本的准则是：

绝对误差准则：

37实际上，要完全保证，是很困难的相对误差推则：

其中表示规定的误差量。

(3)快速性数字仿真是一步一步推进的，即由某一初始值出发，逐步计算，得到，，…。每一步计算所需时间决定了仿真速度。如果第k步计算对应的系统时间间隔为，计算机由计算需要的时间为，那么，若称为实时仿真，则称为超实时仿真，而大多数情况是，对应于离线仿真。

38相对误差推则：其中表示规定的误差量。(连续系统数字仿真中最基本的算法是数值积分算法。

对于形如的系统，已知系统变量y的初始条件，现在要计算y随时间变化的过程。39连续系统数字仿真中最基本的算法是数值积分算法。对于计算过程可以这样考虑(参见图3.2)：首先求出初始点的，对微分方程积分，可以写作

（3.1）图3.2数值积分原理

40计算过程可以这样考虑(参见图3.2)：首先求出初始点的，欧拉法用矩形面积近似表示积分结果，也就是当时，的近似值为：

重复上述做法，当时，

（3.2）所以，对任意时刻，有：

图3.2所示曲线下的面积就是，由于难以得到积分的数值表达式，人们对数值积分方法进行了长期探索，其中欧拉法是最经典的近似方法。

41欧拉法用矩形面积近似表示积分结果，也就是当时，为进一步提高计算精度，人们提出了“梯形法”。令，已知时的近似值，那么，梯形法近似积分形式如下式所示：

令称为第步的计算步距。若积分过程中步距不变，则可以证明，欧拉法的截断误差正比于。

42为进一步提高计算精度，人们提出了“梯形法”。令（3.3）可见，梯形法是隐函数形式。采用这种积分方法，其最简单的预报-校正方法是用欧拉法估计初值，用梯形法校正，即：

（3.4）（3.5）(3.4)式称为预报公式，(3.5)式称为校正公式。用欧拉法估计的值，代入校正公式得到的校正值。

43（3.3）可见，梯形法是隐函数形式。采用这种积分方法，其最简设是规定的足够小正数，称作允许误差。若，称作第一次校正；，称作第二次校正。通过反复迭代，直到满足，这时就是满足误差要求的校正值。44设是规定的足够小正数，称作允许误差。若上述方法是针对微分方程在已知初值情况下进行求解，因此也称为微分方程初值问题数值计算法，为统一起见，本课程中称为数值积分法。连续系统数字仿真的离散化方法有两类，即数值积分方法和离散相似方法，本课程讨论数值积分法。45上述方法是针对微分方程在已知初值情况下进行求数值积分方法采用递推方式进行运算，而采用不同的积分方法会引进不同的计算误差，为了提高计算精度，往往会增加运算量。就同一种积分算法而言，为提高计算精度，减小积分步距h，这样，将降低计算速度。

46数值积分方法采用递推方式进行运算，而采用不同的积分方因此，计算精度和速度是连续系统仿真中常见的一对矛盾，也是数字仿真中要求解决的问题之一。也就是说，选择合适的算法、合适的软硬件环境，在保证计算精度的前提下，考虑怎样提高仿真的速度。

经典的数值积分法可分为单步法与多步法，下面我们将分别来介绍这两类方法中最常用的算法。

47因此，计算精度和速度是连续系统仿真中常见的一对矛盾，也3.2龙格—库塔法

3.2.1龙格—库塔法基本原理

在连续系统仿真中，主要的数值计算工作是对一阶微分方程进行求解。

因为

(3.6)若令

(3.7)则有

(3.8)483.2龙格—库塔法3.2.1龙格—库塔法基本原理因此主要问题就是如何对进行数值求解，即如何对进行积分。这通常称作“右端函数”计算问题。

下面给出常见的几种龙格—库塔积分公式：其中，。

（3.9）这种递推公式的截断误差正比于的三次方。

1）二阶龙格—库塔法(简称RK—2)49因此主要问题就是如何对进行数值求解，即如何对2）四阶龙格—库塔法(简称RK—4)公式

（3.10）其中：由于这组计算公式有较高的精度，所以在数字仿真中应用较为普通。

502）四阶龙格—库塔法(简称RK—4)公式（3.10）其中：所有龙格—库塔公式都有以下特点：

(1)在计算时只用到，而不直接用，等项。换句话说，在后一步的计算中，仅仅利用前一步的计算结果，所以称为单步法。显然它不仅能使存储量减小，而且此法可以自启动，即已知初值后，不必用别的方法来帮助，就能由初值逐步计算得到后续各时间点上的仿真值。

(2)步长在整个计算中并不要求固定，可以根据精度要求改变，但是在一步中算若干个系数(俗称龙格—库塔系数)，则必须用同一个步长。

51所有龙格—库塔公式都有以下特点：(1)在计算时只用(3)龙格—库塔法的精度取决于步长的大小及方法的阶次。许多计算实例表明：为达到相同的精度，四阶方法的可以比二阶方法的大10倍，而四阶方法的每步计算量仅比二阶方法大1倍，所以总的计算量仍比二阶方法小。正是由于上述原因，一般系统进行数字仿真常用四阶龙格—库塔公式。值得指出的是：高于四阶的方法由于每步计算量将增加较多，而精度提高不快，因此使用得也比较少。

52(3)龙格—库塔法的精度取决于步长的大小及方法的52(4)若在展开台劳级数时，只取这一项，而将以及以上的项都略去，则可得这就是欧拉公式，因此欧拉公式也可看作是一阶龙格—库塔公式。它的截断误差正比于，是精度最低的一种数值积分公式。

53(4)若在展开台劳级数时，只取这一项，而53其中，。

（3.9）这种递推公式的截断误差正比于的三次方。

二阶龙格—库塔法(简称RK—2)54其中，其中：四阶龙格—库塔法(简称RK—4)公式

55其中：四阶龙格—库塔法(简称RK—4)公式55四阶龙格-库塔-墨森法误差估计公式56四阶龙格-库塔-墨森法误差估计公式56二阶龙格-库塔-费尔别格法误差估计公式57二阶龙格-库塔-费尔别格法误差估计公式57四阶龙格-库塔-夏普法误差估计公式58四阶龙格-库塔-夏普法误差估计公式583.2.2实时龙格—库塔法

仿真模型的运行速度往往与实际系统运行的速度不同。然而，当有实际的装置或被训练的人介入仿真过程时，就要求仿真模型的运行速度往往与实际系统运行的速度保持一致，这称为实时仿真。

一般的数值积分法难以满足实时仿真的要求，这不仅仅是因为由这些方法所得到的模型的执行速度较慢，而且这些方法的机理不符合实时仿真的特点。593.2.2实时龙格—库塔法仿真模型的运行速度往RK—2的公式如下：下面让我们以二阶龙格—库塔法为例来分析。假设对如下一般形式的系统进行仿真：

60RK—2的公式如下：下面让我们以二阶龙格—库塔法为例即在一个计算步长内分两子步：首先在时刻利用当前的，计算。假定在／2的时间内计算机正好计算一次右端函数，然后，在时刻计算，尽管此时已经得到，但无法得到。实时仿真除了要满足执行速度的要求外，还要求实时接收外部输入，并实时产生输出。为此，或对也进行预报(这将加大仿真误差)，或将仿真执行延迟／2。这就是说，输出要滞后半个计算步距。与此相类似，RK—4公式也不适用于实时仿真。

61即在一个计算步长内分两子步：首先在时刻利用当61（3.11）首先在时刻利用当前的，计算。然后在时刻计算，此时已经得到，由于计算一次右端函数需要／2的时间，也可得到，的计算就不会引入新的误差，并能实时输出。

为了克服这个缺陷，人们提出了如下形式的实时二阶龙格—库塔法：62（3.11）首先在时刻利用当前的，计四阶五级的实时RK计算公式如下：（3.12）该公式的特点是：将一个仿真步分为5个子步，每个子步均能对外部输入实时采样，并保证了实时输出。

63四阶五级的实时RK计算公式如下：（3.12）该公式的特点是：龙格-库塔法举例[例3.1]给定微分方程及初始条件：，，试用4阶龙格-库塔法计算两步，求得及，其中取。64龙格-库塔法举例[例3.1]给定微分方程及初始条件：[解]初始时刻，计算计算65[解]初始时刻，计算计算65计算66计算66计算67计算67计算因此68计算因此68现在，计算计算69现在，计算计计算70计算70计算计算因此71计算计算因此71[例3.2]

(1)给定微分方程组及其初始条件而，试用4阶龙格-库塔法计算一步，求得，其中。72[例3.2](1)给定微分方程组及其初始条件而，试用4阶[解]对于本题而言，除了和以外，其余的量都是2维列向量。令：，，而其中的73[解]对于本题而言，除了和以外，其余的首先，令计算74首先，令计算74即75即75计算76计算76即77即77计算78计算78即79即79计算80计算80即81即81即82即82(2)给定带有时滞的微分方程组及其初始条件试分析4阶龙格-库塔法的算法结构。83(2)给定带有时滞的微分方程组及其初始条件试分析4阶龙格-库这是一个带有状态时滞的微分方程组，通常的4阶龙格-库塔法是无法处理的。必须对4阶龙格-库塔法进行必要的改进。具体方法如下：84这是一个带有状态时滞的微分方程组，通常的4阶龙令85令85对本例而言，4阶龙格-库塔法公式具体为其中不妨称为改进的龙格-库塔法86对本例而言，4阶龙格-库塔法公式具体为其中不妨称为改进的龙格具体为87具体为87888889899090919192929393除了对上述过程进行改进以外，其余的都和常规的4阶龙格-库塔法一样。主要改进措施是在微分方程右端函数的自变量中增加一项带有时滞的状态向量，即：改进94除了对上述过程进行改进以外，其余的都和常规的4阶本次课内容小结连续系统数字仿真离散化离散相似法数值积分法欧拉法梯形法龙格-库塔法稳定性快速性准确性预报-校正法RK-2RK-4实时RK-2实时四阶五级RK95本次课内容小结连续系统数字仿真离散化离散相似法数值积分法欧拉则式（2.87）可写为

（2.88）

N维输出向量2n+1维参数向量N维噪声向量N×(2n+1)维测量矩阵最小二乘法：回顾96则式（2.87）可写为（2.88）N维输出向量2n+1维可得的最小二乘估计

（2.98）

J为极小值的充分条件是

（2.99）

即矩阵为正定矩阵，或者说矩阵是非奇异的。

97可得的最小二乘估计（2.98）J为极小值的充分条件

最小二乘估计的概率性质

如果ξ(k)是不相关随机数序列，且均值为0。1）无偏性2）一致性3）渐进正态性性辅助变量法、广义最小二乘法如果ξ是均值为0且服从正态分布的白噪声向量，则最小二乘参数估计值服从正态分布。98最小二乘估计的概率性质如果ξ(k)是不相关随机数序列，且2.8辅助变量法

现在开始讨论如何克服最小二乘法的有偏估计问题。

对于原辨识方程

（2.157）当是不相关随机序列时，最小二乘法可以得到参数向量的一致性无偏估计。但是，在实际应用中往往是相关随机序列。

992.8辅助变量法现在开始讨论如何克服最小二乘法的式中Q是非奇异的。假定存在着一个的矩阵Z（与同阶数）,满足约束条件

（2.158）用乘以式（2.157）等号两边得（2.159）由上式可得

（2.160）100式中Q是非奇异的。假定存在着一个如果取

（2.161）作为的估值，则称估值为辅助变量估值，矩阵Z称为辅助变量矩阵，Z中的元素称为辅助变量。

从式（2.161）可以看到，与最小二乘法估值的计算公式（2.98）具有相同的形式，因此计算比较简单。

根据式（2.160）和式（2.161）可得

（2.162）101如果取（2.161）作为的估值，则称估值为辅（2.163）当N很大时，对上式等号两边取极限得

（2.164）根据式（2.157）所假定的约束条件，可得

因此辅助变量估计是无偏估计。

102（2.163）当N很大时，对上式等号两边取极限得（2.16剩下的问题是如何选择辅助变量，即如何确定辅助变量矩阵Z的各个元素。选择辅助变量的基本原则是式（2.158）所给出的两个条件必须得到满足。

这可以简单地理解为所选择的辅助变量应与不相关，但与和中的强烈相关。103剩下的问题是如何选择辅助变量，即如何确定辅助变量矩阵Z可以有各种选择方法，下面介绍几种常用的选择方法。1）递推辅助变量参数估计法

辅助变量取作,是辅助模型

（2.165）的输出向量的元素，辅助变量矩阵Z为

104Z可以有各种选择方法，下面介绍几种常用的选择方法。1）递推辅（2.166）（2.167）105（2.166）（2.167）10递推辅助变量参数估计法计算步骤：106递推辅助变量参数估计法计算步骤：112）自适应滤波法

这种方法所选择的辅助变量和辅助变量矩阵Z的形式与上一种方法完全相同，只是辅助模型中参数向量的估计方法与上一种方法有所不同。取

（2.168）式中：取；d取；为k时刻所得到的参数向量估计值。当是持续激励信号时，所选的辅助变量可以满足式（2.158）所给出的2个约束条件。

1072）自适应滤波法这种方法所选择的辅助变量和2.9最小二乘类辨识算法的比较1.基本最小二乘法（LS）（1）对白噪声，参数估计是无偏的一致的最小方差估计。（2）对有色噪声，参数估计有偏，但具有收敛性。（3）对未知的直流分量敏感。（4）对高阶系统，优于其它算法，可作为其它辨识算法的较好的起始算法。（5）一次完成算法精度较高，但逆矩阵计算量大，不适于在线辨识。（6）数据较多时，占用存储量和机时较多。1082.9最小二乘类辨识算法的比较1.基本最小二乘法（LS）（2.递推最小二乘法（RLS）（1）对白噪声，参数估计是无偏的一致的最小方差估计。（2）对有色噪声，参数估计有偏，但具有收敛性。（3）对未知的直流分量敏感。（4）对高阶系统，优于其它算法，可作为其它辨识算法的较好的起始算法。（5）计算量少，占用内存较少。（6）适用于在线辨识，采用自适应算法后，可用于时变系统的辨识。1092.递推最小二乘法（RLS）（1）对白噪声，参数估计是无偏3.辅助变量法（IV）（1）估计值是弱一致估计的。（2）对初态敏感，初态选取得不合适，就可能不收敛。（3）算法简单，计算量少，但不能得到噪声模型。1103.辅助变量法（IV）（1）估计值是弱一致估计的。（2）2.10系统结构辨识模型阶的确定

在一些实际问题中，模型的阶可以按理论推导获得，而在另一些实际问题中，模型的阶却无法用理论推导的方法确定，需要对模型的阶进行辨识。下面介绍几种常用的模型阶的确定方法。

2.10.1按残差方差定阶

一种简单而有效的方法就是选定模型阶数n的不同取值，按估计误差方差最小或F检验来确定模型的阶。1)按估计误差方差最小定阶

考虑系统模型

1112.10系统结构辨识模型阶的确定在一些实际问题（2.178）式中：为输出；为输入。

设是均值为0、方差为的白噪声序列。用最小二乘法求出的估值。具体为

（2.179）（2.180）（2.181）112（2.178）式中：为输出；为输入。(2.182)(2.183)残差为

113(2.182)(2.183)残差为18如图2.19所示，对某一系统，当=1，2，…时，随着的增加而减小。如果为正确的阶，则在时，出现最后一次陡峭的下降，再增大，则保持不变或只有微小的变化。

图2.19所示的例子，=3。

114如图2.19所示，对某一系统，当=1，2，…时，图2.19曲线图

115图2.19曲线图202)确定模型阶的F检验法

由于随着的增加而减小，在阶数的增大过程中，我们对那个使显著减小的阶感兴趣。为此，引入准则

（2.187）式中表示具有N对输入和输出数据、有个模型参数的系统估计误差的平方和。1162)确定模型阶的F检验法由于随着的增表2.1某一系统计算结果

0.993.159.439.6750.94416.56418.73426.40447.25469.64592.65654321计算时取，，。从表2.1可以看出：当时，t的减小是显著的；当时，t的减小是不显著的。所以该系统的阶数可选为3。

对某一系统的计算结果如表2.1所列。

117表2.1某一系统计算结果0.993.159.4392.10.2确定阶的Akaike信息准则

与上述2个准则不同，Akaike信息准则（AIC，AkaikeInformationCriterion）是一个考虑了模型复杂性的准则。这个准则定义为

（2.189）式中：L是模型的似然函数；是模型中的参数数目。当AIC为最小的那个模型就是最佳模型。这个准则是Akaike总结了时间序列统计建模的发展史，在企图对一个复杂系统寻找近似模型的概率论的大量探索启示下，借助信息论而提出的一个合理的确定阶的准则。在一组可供选择的随机模型中，AIC最小的那个模型是一个可取的模型。这个准则的优点就在于它是一个完全客观的准则，应用这个准则时，不要求建模人员主观地判断“陡峭的下降”。

1182.10.2确定阶的Akaike信息准则与上述21)白噪声情况下的AIC定阶公式

假定是均值为0、方差为并且服从正态分布的不相关随机噪声。

考虑系统模型（2.190）式中

根据前面几节的讨论和定义，由式（2.190）可写出关系式

（2.191）1191)白噪声情况下的AIC定阶公式假定是均值为0、方输出变量在条件下的似然函数为

（2.192）对上式取对数可得

（2.193）求使为最大的的估值。根据可得

与前述的最小二乘估计一致。按照可得

120输出变量在条件下的似然函数为（2.192）对上式取因此

即（2.194）式中c为一常数。

式（2.194）给出了AIC定义中的第一项。

为

待估的参数，，…，，，，…，及，共有个，即，因而有121因此即（2.194）式中c为一常数。式（2.19即

可去掉上式中的常数项，则（2.195）选取不同的阶数和，按式（2.195）计算AIC，可得最优阶数和。在式（2.195）中加进项，表示对不同，若相近时，则取较小的模型。

122即可去掉上式中的常数项，则（2.195）选取不同的例2.4系统模型为

式中：是均值为0、方差为1且服从正态分布的不相关随机噪声；输入信号采用伪随机数。辨识模型采用的形式为

数据长度取N=1024。为了避免非平稳过程的影响，去掉前300个数据，取=1，2，3，4，=1，2，3，4，分别计算AIC（，），即123例2.4系统模型为式中：是均值为0、方差为1且服计算结果如表2.2所列。显然，应取=3，=2，可见利用AIC确定的模型阶次与系统的真实阶次相同。

表2.2不同和所对应的AIC16.218

15.108

15.931

417.649

15.599

14.070

25.864

316.800

30.393

51.085

280.046

223.380

97.353

341.766

1022.94

14321124计算结果如表2.2所列。显然，应取=3，=2，可2)有色噪声情况下的AIC定阶公式

有色噪声情况下的系统模型可以表示为

（2.196）式中

是均值为0、方差为，且服从正态分布的不相关随机噪声。

与前面的讨论相类似，可得

（2.197）1252)有色噪声情况下的AIC定阶公式有色噪声情况下的服从正态分布N（0，1），为二位式伪随机序列，数据长度N=300。假定模型阶，利用极大似然法估计模型参数，根据AIC的计算结果，可得：当时，AIC最小，其结果与所给系例2.5：设系统模型为

式中

统模型相符合。126服从正态分布N（0，1），为二位式伪随机序列，F检验法和AIC准则法，可以在参数估计的过程中同时辨识阶次，所以应用上比较方便。阶次辨识和参数辨识两者是互相依赖的，也就是说进行参数估计时需要已知阶次；而辨识阶次时又要利用参数估计值，二者是不可分离的。127F检验法和AIC准则法，可以在参数估计的过程中同时辨识阶次，第三章连续系统的数字仿真本章讨论经典的连续系统数字仿真的原理与方法，内容包括连续系统数字仿真的基本概念、经典的数值积分法、经典的线性多步法等。3.1节讨论离散化原理及要求，是连续系统仿真的基础；3.2节对经典的数值积分法——龙格—库塔法及其他典型的数值积分法仿真建模原理进行详细分析，并通过实例说明其应用要点；3.3节对经典的线性多步法进行了介绍。

128第三章连续系统的数字仿真本章讨论经典的连续系统3.1离散化原理及要求

在数字计算机上对连续系统进行仿真时，首先遇到的问题是，数字计算机的数值及时间均具有离散性，而被仿真系统的数值及时间均具有连续性，后者如何用前者来实现。

连续系统仿真，从本质上是从时间、数值两个方面对原系统进行离散化，并选择合适的数值计算方法来近似积分运算，由此得到离散模型来近似原连续模型。如何保证离散模型的计算结果从原理上确能代表原系统的行为，这是连续系统数字仿真首先必须解决的问题。

1293.1离散化原理及要求在数字计算机上对连续系统图3.1相似原理示意图

130图3.1相似原理示意图35设系统模型(参见图3.1)为，其中为输入变量，为系统变量。令仿真时间间隔为h，离散化后的输入变量为，系统变量为，其中表示。如果，，即，

(对所有k＝0，1，2…)，则可认为两模型等价，这称为相似原理。

131设系统模型(参见图3.1)为，实际上，要完全保证，是很困难的。进一步分析离散化引入的误差，随着计算机技术的发展，由计算机字长引入的舍入误差可以忽略，关键是数值积分算法，也称为仿真建模方法。相似原理用于仿真时，对仿真建模方法有三个基本要求：

(1)稳定性

若原连续系统是稳定的，则离散化后得到的仿真模型也应是稳定的。(2)准确性

有不同的准确性评价准则，最基本的准则是：

绝对误差准则：

132实际上，要完全保证，是很困难的相对误差推则：

其中表示规定的误差量。

(3)快速性数字仿真是一步一步推进的，即由某一初始值出发，逐步计算，得到，，…。每一步计算所需时间决定了仿真速度。如果第k步计算对应的系统时间间隔为，计算机由计算需要的时间为，那么，若称为实时仿真，则称为超实时仿真，而大多数情况是，对应于离线仿真。

133相对误差推则：其中表示规定的误差量。(连续系统数字仿真中最基本的算法是数值积分算法。

对于形如的系统，已知系统变量y的初始条件，现在要计算y随时间变化的过程。134连续系统数字仿真中最基本的算法是数值积分算法。对于计算过程可以这样考虑(参见图3.2)：首先求出初始点的，对微分方程积分，可以写作

（3.1）图3.2数值积分原理

135计算过程可以这样考虑(参见图3.2)：首先求出初始点的，欧拉法用矩形面积近似表示积分结果，也就是当时，的近似值为：

重复上述做法，当时，

（3.2）所以，对任意时刻，有：

图3.2所示曲线下的面积就是，由于难以得到积分的数值表达式，人们对数值积分方法进行了长期探索，其中欧拉法是最经典的近似方法。

136欧拉法用矩形面积近似表示积分结果，也就是当时，为进一步提高计算精度，人们提出了“梯形法”。令，已知时的近似值，那么，梯形法近似积分形式如下式所示：

令称为第步的计算步距。若积分过程中步距不变，则可以证明，欧拉法的截断误差正比于。

137为进一步提高计算精度，人们提出了“梯形法”。令（3.3）可见，梯形法是隐函数形式。采用这种积分方法，其最简单的预报-校正方法是用欧拉法估计初值，用梯形法校正，即：

（3.4）（3.5）(3.4)式称为预报公式，(3.5)式称为校正公式。用欧拉法估计的值，代入校正公式得到的校正值。

138（3.3）可见，梯形法是隐函数形式。采用这种积分方法，其最简设是规定的足够小正数，称作允许误差。若，称作第一次校正；，称作第二次校正。通过反复迭代，直到满足，这时就是满足误差要求的校正值。139设是规定的足够小正数，称作允许误差。若上述方法是针对微分方程在已知初值情况下进行求解，因此也称为微分方程初值问题数值计算法，为统一起见，本课程中称为数值积分法。连续系统数字仿真的离散化方法有两类，即数值积分方法和离散相似方法，本课程讨论数值积分法。140上述方法是针对微分方程在已知初值情况下进行求数值积分方法采用递推方式进行运算，而采用不同的积分方法会引进不同的计算误差，为了提高计算精度，往往会增加运算量。就同一种积分算法而言，为提高计算精度，减小积分步距h，这样，将降低计算速度。

141数值积分方法采用递推方式进行运算，而采用不同的积分方因此，计算精度和速度是连续系统仿真中常见的一对矛盾，也是数字仿真中要求解决的问题之一。也就是说，选择合适的算法、合适的软硬件环境，在保证计算精度的前提下，考虑怎样提高仿真的速度。

经典的数值积分法可分为单步法与多步法，下面我们将分别来介绍这两类方法中最常用的算法。

142因此，计算精度和速度是连续系统仿真中常见的一对矛盾，也3.2龙格—库塔法

3.2.1龙格—库塔法基本原理

在连续系统仿真中，主要的数值计算工作是对一阶微分方程进行求解。

因为

(3.6)若令

(3.7)则有

(3.8)1433.2龙格—库塔法3.2.1龙格—库塔法基本原理因此主要问题就是如何对进行数值求解，即如何对进行积分。这通常称作“右端函数”计算问题。

下面给出常见的几种龙格—库塔积分公式：其中，。

（3.9）这种递推公式的截断误差正比于的三次方。

1）二阶龙格—库塔法(简称RK—2)144因此主要问题就是如何对进行数值求解，即如何对2）四阶龙格—库塔法(简称RK—4)公式

（3.10）其中：由于这组计算公式有较高的精度，所以在数字仿真中应用较为普通。

1452）四阶龙格—库塔法(简称RK—4)公式（3.10）其中：所有龙格—库塔公式都有以下特点：

(1)在计算时只用到，而不直接用，等项。换句话说，在后一步的计算中，仅仅利用前一步的计算结果，所以称为单步法。显然它不仅能使存储量减小，而且此法可以自启动，即已知初值后，不必用别的方法来帮助，就能由初值逐步计算得到后续各时间点上的仿真值。

(2)步长在整个计算中并不要求固定，可以根据精度要求改变，但是在一步中算若干个系数(俗称龙格—库塔系数)，则必须用同一个步长。

146所有龙格—库塔公式都有以下特点：(1)在计算时只用(3)龙格—库塔法的精度取决于步长的大小及方法的阶次。许多计算实例表明：为达到相同的精度，四阶方法的可以比二阶方法的大10倍，而四阶方法的每步计算量仅比二阶方法大1倍，所以总的计算量仍比二阶方法小。正是由于上述原因，一般系统进行数字仿真常用四阶龙格—库塔公式。值得指出的是：高于四阶的方法由于每步计算量将增加较多，而精度提高不快，因此使用得也比较少。

147(3)龙格—库塔法的精度取决于步长的大小及方法的52(4)若在展开台劳级数时，只取这一项，而将以及以上的项都略去，则可得这就是欧拉公式，因此欧拉公式也可看作是一阶龙格—库塔公式。它的截断误差正比于，是精度最低的一种数值积分公式。

148(4)若在展开台劳级数时，只取这一项，而53其中，。

（3.9）这种递推公式的截断误差正比于的三次方。

二阶龙格—库塔法(简称RK—2)149其中，其中：四阶龙格—库塔法(简称RK—4)公式

150其中：四阶龙格—库塔法(简称RK—4)公式5