《系统建模与仿真-第5次课-第三章课件.ppt》由会员分享，可在线阅读，更多相关《系统建模与仿真-第5次课-第三章课件.ppt（95页珍藏版）》请在taowenge.com淘文阁网|工程机械CAD图纸|机械工程制图|CAD装配图下载|SolidWorks_CaTia_CAD_UG_PROE_设计图分享下载上搜索。
1、则式（2.87）可写为（2.88）N维输出向量2n+1维参数向量N维噪声向量N(2n+1)维测量矩阵最小二乘法：最小二乘法：回顾回顾1可得 的最小二乘估计（2.98）J为极小值的充分条件是（2.99）即矩阵 为正定矩阵正定矩阵，或者说矩阵 是非奇异的。2 最小二乘估计的概率性质最小二乘估计的概率性质 如果如果(k)是不相关随机数序列，且均值为是不相关随机数序列，且均值为0。1）无偏性无偏性 2）一致性）一致性3）渐进正态性性渐进正态性性辅助变量法、广义最小二乘法辅助变量法、广义最小二乘法 如果如果是均值为是均值为0且服从正态分布的白噪声向量，则最小且服从正态分布的白噪声向量，则最小二乘参数估计
2、值服从正态分布。二乘参数估计值服从正态分布。32.8 辅助变量法辅助变量法 现在开始讨论如何克服最小二乘法的有偏估计问题。对于原辨识方程（2.157）当 是不相关随机序列时，最小二乘法可以得到参数向量 的一致性无偏估计。但是，在实际应用中 往往是相关随机序列。4式中Q是非奇异的。假定存在着一个 的矩阵Z（与 同阶数）,满足约束条件（2.158）用 乘以式（2.157）等号两边得（2.159）由上式可得（2.160）5如果取（2.161）作为 的估值，则称估值 为辅助变量估值，矩阵Z称为辅助变量矩阵，Z中的元素称为辅助变量。从式（2.161）可以看到，与最小二乘法估值 的计算公式（2.98）具有
3、相同的形式，因此计算比较简单。根据式（2.160）和式（2.161）可得（2.162）6（2.163）当N很大时，对上式等号两边取极限得（2.164）根据式（2.157）所假定的约束条件，可得 因此辅助变量估计是无偏估计。7 剩下的问题是如何选择辅助变量，即如何确定辅助变量矩阵Z的各个元素。选择辅助变量的基本原则是式（2.158）所给出的两个条件必须得到满足。这可以简单地理解为所选择的辅助变量应与 不相关，但与 和 中的 强烈相关。8Z可以有各种选择方法，下面介绍几种常用的选择方法。1）递推辅助变量参数估计法 辅助变量取作 ,是辅助模型（2.165）的输出向量 的元素，辅助变量矩阵Z为 9（2
4、.166）（2.167）10递推辅助变量参数估计法递推辅助变量参数估计法 计算步骤：计算步骤：112）自适应滤波法 这种方法所选择的辅助变量 和辅助变量矩阵Z的形式与上一种方法完全相同，只是辅助模型中参数向量 的估计方法与上一种方法有所不同。取（2.168）式中：取 ；d 取 ；为k时刻所得到的参数向量估计值。当 是持续激励信号时，所选的辅助变量可以满足式（2.158）所给出的2个约束条件。122.9最小二乘类辨识算法的比较最小二乘类辨识算法的比较1.基本最小二乘法（基本最小二乘法（LS）（1）对白噪声，参数估计是无偏的一致的最小方差估计。（2）对有色噪声，参数估计有偏，但具有收敛性。（3）对
5、未知的直流分量敏感。（4）对高阶系统，优于其它算法，可作为其它辨识算法的较好的起始算法。（5）一次完成算法精度较高，但逆矩阵计算量大，不适于在线辨识。（6）数据较多时，占用存储量和机时较多。132.递推最小二乘法（递推最小二乘法（RLS）（1）对白噪声，参数估计是无偏的一致的最小方差估计。（2）对有色噪声，参数估计有偏，但具有收敛性。（3）对未知的直流分量敏感。（4）对高阶系统，优于其它算法，可作为其它辨识算法的较好的起始算法。（5）计算量少，占用内存较少。（6）适用于在线辨识，采用自适应算法后，可用于时变系统的辨识。143.辅助变量法（辅助变量法（IV）（1）估计值是弱一致估计的。（2）对初
6、态敏感，初态选取得不合适，就可能不收敛。（3）算法简单，计算量少，但不能得到噪声模型。152.10 系统结构辨识系统结构辨识 模型阶的确定模型阶的确定 在一些实际问题中，模型的阶可以按理论推导获得，而在另一些实际问题中，模型的阶却无法用理论推导的方法确定，需要对模型的阶进行辨识。下面介绍几种常用的模型阶的确定方法。2.10.1 按残差方差定阶按残差方差定阶 一种简单而有效的方法就是选定模型阶数n的不同取值，按估计误差方差最小或F检验来确定模型的阶。1)按估计误差方差最小定阶 考虑系统模型 16（2.178）式中：为输出；为输入。设 是均值为0、方差为 的白噪声序列。用最小二乘法求出 的估值。具
7、体为（2.179）（2.180）（2.181）17(2.182)(2.183)残差为 18如图2.19所示，对某一系统，当 =1，2，时，随着 的增加而减小。如果 为正确的阶，则在 时，出现最后一次陡峭的下降，再增大，则 保持不变或只有微小的变化。图2.19所示的例子，=3。19图2.19 曲线图 202)确定模型阶的F检验法 由于 随着 的增加而减小，在阶数 的增大过程中，我们对那个使 显著减小的阶 感兴趣。为此，引入准则（2.187）式中 表示具有N对输入和输出数据、有 个模型参数的系统估计误差的平方和。21表2.1 某一系统计算结果 0.99 3.15 9.43 9.67 50.94 4
8、16.56 418.73 426.40 447.25 469.64 592.65 6 5 4 3 2 1计算时取 ，。从表2.1可以看出：当 时，t 的减小是显著的；当 时，t 的减小是不显著的。所以该系统的阶数可选为3。对某一系统的计算结果如表2.1所列。222.10.2 确定阶的确定阶的Akaike信息准则信息准则 与上述2个准则不同，Akaike信息准则（AIC，Akaike Information Criterion）是一个考虑了模型复杂性的准则。这个准则定义为（2.189）式中：L是模型的似然函数；是模型中的参数数目。当AIC为最小的那个模型就是最佳模型。这个准则是Akaike总结了
9、时间序列统计建模的发展史，在企图对一个复杂系统寻找近似模型的概率论的大量探索启示下，借助信息论而提出的一个合理的确定阶的准则。在一组可供选择的随机模型中，AIC最小的那个模型是一个可取的模型。这个准则的优点就在于它是一个完全客观的这个准则的优点就在于它是一个完全客观的准则，应用这个准则时，不要求建模人员主观地判断准则，应用这个准则时，不要求建模人员主观地判断“陡峭的下降陡峭的下降”。231)1)白噪声情况下的白噪声情况下的AIC定阶公式定阶公式 假定 是均值为0、方差为 并且服从正态分布的不相关随机噪声。考虑系统模型（2.190）式中 根据前面几节的讨论和定义，由式（2.190）可写出关系式（
10、2.191）24输出变量 在条件 下的似然函数为（2.192）对上式取对数可得（2.193）求使 为最大的 的估值 。根据 可得 与前述的最小二乘估计一致。按照 可得 25因此 即（2.194）式中c为一常数。式（2.194）给出了AIC定义中的第一项。为 待估的参数，及 ，共有 个，即 ，因而有26即 可去掉上式中的常数项，则（2.195）选取不同的阶数 和 ，按式（2.195）计算AIC，可得最优阶数 和 。在式（2.195）中加进 项，表示对不同 ，若 相近时，则取 较小的模型。27例2.4 系统模型为 式中：是均值为0、方差为1且服从正态分布的不相关随机噪声；输入信号 采用伪随机数。辨
11、识模型采用的形式为 数据长度取N=1024。为了避免非平稳过程的影响，去掉前300个数据，取 =1，2，3，4，=1，2，3，4，分别计算AIC（，），即28计算结果如表2.2 所列。显然，应取 =3，=2，可见利用AIC确定的模型阶次与系统的真实阶次相同。表2.2 不同 和 所对应的AIC 16.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 14321292)2)有色噪声情况下的有色噪声情况下的AIC定阶公式定阶公式 有色噪
12、声情况下的系统模型可以表示为（2.196）式中 是均值为0、方差为 ，且服从正态分布的不相关随机噪声。与前面的讨论相类似，可得（2.197）30 服从正态分布N（0，1），为二位式伪随机序列，数据长度N=300。假定模型阶 ，利用极大似然法估计模型参数，根据AIC的计算结果，可得：当 时，AIC最小，其结果与所给系例2.5：设系统模型为 式中 统模型相符合。31F检验法和AIC准则法，可以在参数估计的过程中同时辨识阶次，所以应用上比较方便。阶次辨识和参数辨识两者是互相依赖的，也就是说进行参数估计时需要已知阶次；而辨识阶次时又要利用参数估计值，二者是不可分离的。32第三章第三章 连续系统的数字仿
13、真连续系统的数字仿真 本章讨论经典的连续系统数字仿真的原理与方法，本章讨论经典的连续系统数字仿真的原理与方法，内容包括连续系统数字仿真的基本概念、经典的数值积内容包括连续系统数字仿真的基本概念、经典的数值积分法、经典的线性多步法等。分法、经典的线性多步法等。3.1节讨论离散化原理及要求，是连续系统仿真的基础；节讨论离散化原理及要求，是连续系统仿真的基础；3.2节对经典的数值积分法节对经典的数值积分法龙格龙格库塔法及其他典库塔法及其他典型的数值积分法仿真建模原理进行详细分析，并通过实型的数值积分法仿真建模原理进行详细分析，并通过实例说明其应用要点；例说明其应用要点；3.3节对经典的线性多步法进行
14、了介绍。节对经典的线性多步法进行了介绍。333.1 离散化原理及要求离散化原理及要求 在数字计算机上对连续系统进行仿真时，首先遇到的问题是，数字计算机的数值及时间均具有离散性，而被仿真系统的数值及时间均具有连续性，后者如何用前者来实现。连续系统仿真，从本质上是从时间、数值时间、数值两个方面对原系统进行离散化离散化，并选择合适的数值计算方法来近似积分运算，由此得到离散模型来近似近似原连续模型。如何保证离散模型的计算结果从原理上确能代表原系统的行为，这是连续系统数字仿真首先必须解决的问题。34图3.1 相似原理示意图 35 设系统模型(参见图3.1)为 ，其中 为输入变量，为系统变量。令仿真时间间
15、隔为h，离散化后的输入变量为 ，系统变量为 ，其中 表示 。如果 ，即 ，(对所有k0，1，2)，则可认为两模型等价，这称为相似原理相似原理。36实际上，要完全保证 ，是很困难的。进一步分析离散化引入的误差，随着计算机技术的发展，由计算机字长引入的舍入误差可以忽略，关键是数值积分算法，也称为仿真建模方法。相似原理用于仿真时，对仿真建模方法有三个基本要求：(1)稳定性 若原连续系统是稳定的，则离散化后得到的仿真模型也应是稳定的。(2)准确性 有不同的准确性评价准则，最基本的准则是：绝对误差准则：37相对误差推则：其中 表示规定的误差量。(3)快速性 数字仿真是一步一步推进的，即由某一初始值 出发
16、，逐步计算，得到 ，。每一步计算所需时间决定了仿真速度。如果第k步计算对应的系统时间间隔为 ，计算机由计算 需要的时间为 ，那么，若 称为实时仿真，则 称为超实时仿真，而大多数情况是 ，对应于离线仿真。38 连续系统数字仿真中最基本的算法是数值积分算法。对于形如 的系统，已知系统变量y的初始条件 ，现在要计算y随时间变化的过程 。39计算过程可以这样考虑(参见图3.2)：首先求出初始点 的，对微分方程积分，可以写作（3.1）图3.2 数值积分原理 40欧拉法用矩形面积近似表示积分结果，也就是当 时，的近似值为 ：重复上述做法，当 时，（3.2）所以，对任意时刻 ，有：图3.2所示曲线下的面积就
17、是 ，由于难以得到 积分的数值表达式，人们对数值积分方法进行了长期探索，其中欧拉法欧拉法是最经典的近似方法。41 为进一步提高计算精度，人们提出了“梯形法梯形法”。令 ，已知 时 的近似值 ，那么，梯形法近似积分形式如下式所示：令 称为第 步的计算步距。若积分过程中步距 不变，则可以证明，欧拉法的截断误差正比于 。42（3.3）可见，梯形法是隐函数形式。采用这种积分方法，其最简单的预报-校正方法是用欧拉法欧拉法估计初值，用梯形法梯形法校正，即：（3.4）（3.5）(3.4)式称为预报公式预报公式，(3.5)式称为校正公式校正公式。用欧拉法估计 的值，代入校正公式得到 的校正值 。43 设 是规
18、定的足够小正数，称作允许误差。若 ，称作第一次校正；，称作第二次校正。通过反复迭代，直到满足 ，这时 就是满足误差要求的校正值。44 上述方法是针对微分方程在已知初值情况下进行求解，因此也称为微分方程初值问题数值计算法，为统一起见，本课程中称为数值积分法。连续系统数字仿真的离散化方法有两类，即数值积数值积分方法分方法和离散相似方法离散相似方法，本课程讨论数值积分法。45 数值积分方法采用递推方式进行运算，而采用不同的积分方法会引进不同的计算误差，为了提高计算精度，往往会增加运算量。就同一种积分算法而言，为提高计算精度，减小积分步距h，这样，将降低计算速度。46 因此，计算精度和速度是连续系统仿
19、真中常见的一对矛盾，也是数字仿真中要求解决的问题之一。也就是说，选择合适的算法、合适的软硬件环境，在保证计算精度的前提下，考虑怎样提高仿真的速度。经典的数值积分法可分为单步法与多步法，下面我们将分别来介绍这两类方法中最常用的算法。473.2 龙格龙格库塔法库塔法 3.2.1 龙格龙格库塔法基本原理库塔法基本原理 在连续系统仿真中，主要的数值计算工作是对一阶微分方程 进行求解。因为(3.6)若令(3.7)则有(3.8)48因此主要问题就是如何对 进行数值求解，即如何对 进行积分。这通常称作“右端函数”计算问题。下面给出常见的几种龙格库塔积分公式：其中 ，。（3.9）这种递推公式的截断误差正比于
20、的三次方。1）二阶龙格）二阶龙格库塔法库塔法(简称简称RK2)492）四阶龙格四阶龙格库塔法库塔法(简称简称RK4)公式公式（3.10）其中：由于这组计算公式有较高的精度，所以在数字仿真中应用较为普通。50所有龙格所有龙格库塔公式都有以下特点：库塔公式都有以下特点：(1)在计算 时只用到 ，而不直接用 ，等项。换句话说，在后一步的计算中，仅仅利用前一步的计算结果，所以称为单步法单步法。显然它不仅能使存储量减小，而且此法可以自启动，即已知初值后，不必用别的方法来帮助，就能由初值逐步计算得到后续各时间点上的仿真值。(2)步长 在整个计算中并不要求固定，可以根据精度要求改变，但是在一步中算若干个系数
21、 (俗称龙格库塔系数)，则必须用同一个步长 。51(3)龙格库塔法的精度取决于步长 的大小及方法的阶次。许多计算实例表明：为达到相同的精度，四阶方法的 可以比二阶方法的 大10倍，而四阶方法的每步计算量仅比二阶方法大1倍，所以总的计算量仍比二阶方法小。正是由于上述原因，一般系统进行数字仿真常用四阶龙格库塔公式。值得指出的是：高于四阶的方法由于每步计算量将增加较多，而精度提高不快，因此使用得也比较少。52(4)若在展开台劳级数时，只取 这一项，而将 以及 以上的项都略去，则可得这就是欧拉公式，因此欧拉公式也可看作是一阶龙格库塔公式。它的截断误差正比于 ，是精度最低的一种数值积分公式。53其中 ，
22、。（3.9）这种递推公式的截断误差正比于 的三次方。二阶龙格二阶龙格库塔法库塔法(简称简称RK2)54其中：四阶龙格四阶龙格库塔法库塔法(简称简称RK4)公式公式 55四阶龙格四阶龙格-库塔库塔-墨森法墨森法误差估计公式误差估计公式56二阶龙格二阶龙格-库塔库塔-费尔别格法费尔别格法误差估计公式误差估计公式57四阶龙格四阶龙格-库塔库塔-夏普法夏普法误差估计公式误差估计公式583.2.2 实时龙格实时龙格库塔法库塔法 仿真模型的运行速度往往与实际系统运行的速度不同。然而，当有实际的装置或被训练的人介入仿真过程时，就要求仿真模型的运行速度往往与实际系统运行的速度保持一致，这称为实时仿真。一般的数
23、值积分法难以满足实时仿真的要求，这一般的数值积分法难以满足实时仿真的要求，这不仅仅是因为由这些方法所得到的模型的执行速度较不仅仅是因为由这些方法所得到的模型的执行速度较慢，而且这些方法的慢，而且这些方法的机理机理不符合实时仿真的特点。不符合实时仿真的特点。59RK2的公式如下：下面让我们以二阶龙格库塔法为例来分析。假设对如下一般形式的系统进行仿真：60即在一个计算步长内分两子步：首先在 时刻利用当前的 ，计算 。假定在 2的时间内计算机正好计算一次右端函数 ，然后，在 时刻计算 ，尽管此时已经得到 ，但 无法得到。实时仿真除了要满足执行速度的要求外，还要求实时接收外部输入，并实时产生输出。为此
24、，或对 也进行预报(这将加大仿真误差)，或将仿真执行延迟 2。这就是说，输出要滞后半个计算步距。与此相类似，RK4公式也不适用于实时仿真。61（3.11）首先在 时刻利用当前的 ，计算 。然后在 时刻计算 ，此时 已经得到，由于计算一次右端函数 需要 2的时间，也可得到，的计算就不会引入新的误差，并能实时输出 。为了克服这个缺陷，人们提出了如下形式的实时二阶龙格实时二阶龙格库塔法库塔法：62四阶五级的实时四阶五级的实时RK计算公式计算公式如下：（3.12）该公式的特点是：将一个仿真步分为5个子步，每个子步均能对外部输入实时采样，并保证了实时输出。63龙格龙格-库塔法举例库塔法举例例例 3.1
25、给定微分方程及初始条件：，试用4阶龙格-库塔法计算两步，求得 及 ，其中取 。64解解 初始时刻 ，计算计算65计算66计算67计算因此68现在 ，计算计算69计算70计算计算因此71例例3.2 (1)给定微分方程组及其初始条件而，试用4阶龙格-库塔法计算一步，求得，其中 。72解解对于本题而言，除了 和 以外，其余的量都是2维列向量。令：，而其中的73首先，令计算74即75计算76即77计算78即79计算80即81即82(2)给定带有时滞的微分方程组及其初始条件试分析4阶龙格-库塔法的算法结构。83 这是一个带有状态时滞的微分方程组，通常的4阶龙格-库塔法是无法处理的。必须对4阶龙格-库塔法
26、进行必要的改进。具体方法如下：84令85对本例而言，4阶龙格-库塔法公式具体为其中不妨称为改进的龙格不妨称为改进的龙格-库塔法库塔法86具体为87888990919293 除了对上述过程进行改进以外，其余的都和常规的4阶龙格-库塔法一样。主要改进措施是在微分方程右端函数的自变量中增加一项带有时滞的状态向量 ，即：改进94本次课内容小结本次课内容小结连续系统数字仿真连续系统数字仿真离散化离散化离散相似法离散相似法数值积分法数值积分法欧拉法欧拉法梯形法梯形法龙格龙格-库塔法库塔法稳定性稳定性快速性快速性准确性准确性预报预报-校正法校正法RK-2RK-4实时实时RK-2实时四阶实时四阶 五级五级RK95