连续系统的数字仿真

连续系统的数字仿真第1页，共134页，2023年，2月20日，星期二第3章连续系统的数字仿真-离散相似法

控制系统的数字仿真就是控制系统的数学模型在数字计算机上求解的过程。从前一章可以看到，控制系统的动态模型一般是用常微分方程、一阶常微分方程组(状态方程)或传递函数来描述的。因此，要在数字计算机上进行这类系统的仿真，必须先将连续模型变换为离散化的模型。这一章主要介绍连续模型离散化的方法及其差分方程的求取，最后介绍仿真程序的设计。第2页，共134页，2023年，2月20日，星期二3．1连续系统的离散化连续系统的响应是随时间连续变化的，但是连续系统的解析解无法用数字计算机求出，只能求出其数值解。也就是说，只能得到连续响应曲线上的有限个点。为此，必须把连续系统离散化，得到差分方程，再用数字计算机求解。这就是把微分运算转化为算术运算的过程。下面先讨论线性定常系统的离散过程。设一线性定常系统为

(3-1)

式中：X为n×1维状态向量；U为r×1维输入向量；A为n×n维状态矩阵；B为n×r维输入矩阵；Y为m×1维输出向量；C为m×n维输出矩阵；D为m×r维传递矩阵。此系统的方框，如图3．1所示。

第3页，共134页，2023年，2月20日，星期二此系统的方框，如图3．1所示。第4页，共134页，2023年，2月20日，星期二第5页，共134页，2023年，2月20日，星期二第6页，共134页，2023年，2月20日，星期二

由图3．3所示的离散结构即可导出连续系统离散后的离散数学模型，这个模型通常为差分方程的形式。由于这个过程使得离散系统与连续系统相似，因此可以认为差分方程式的解序列就是连续系统运动的采样值。这就是离散相似法。严格地讲，系统输出处的采样开关后面也应加上再现环节，才能与原系统相似。但是在仿真时，用计算机也只能得到离散序列的解，所以输出处的再现环节加与不加对于离散解序列都是一样的。实际上输出处的采样开关加与不加也无所谓，只要认为离散后的系统与原系统在采样点上的输出值近似相等就行了。

第7页，共134页，2023年，2月20日，星期二由香农定理可知，为使采样后的信号能无失真地再现，采样频率比信号最高频率要高两倍以上(或采样周期Ts≤0.5Tmin)，且加入的保持器应有如图3．4所示的频率特性。第8页，共134页，2023年，2月20日，星期二3．1．1零阶保持器

这是一种最常见的保持器，它的脉冲传递函数和频率特性如图3．5和图3．6所示。第9页，共134页，2023年，2月20日，星期二第10页，共134页，2023年，2月20日，星期二第11页，共134页，2023年，2月20日，星期二第12页，共134页，2023年，2月20日，星期二3．1．2一阶保持器

一阶保持器也称外推器，它的脉冲传递函数和频率特性如图3．8和图3．9所示。第13页，共134页，2023年，2月20日，星期二第14页，共134页，2023年，2月20日，星期二第15页，共134页，2023年，2月20日，星期二第16页，共134页，2023年，2月20日，星期二3．1．3三角保持器

三角保持器是一种理想保持器，物理上不能实现，数学上也是不能实现的，除非它所再现的信号为一已知信号。这一点从下面的定义中可以看出。它的脉冲传递函数曲线和频率特性曲线如图3．11和图3．12所示。第17页，共134页，2023年，2月20日，星期二第18页，共134页，2023年，2月20日，星期二第19页，共134页，2023年，2月20日，星期二

在实际中有时使用滞后一拍的三角保持器，它的脉冲传递函数曲线和频率特性曲线如图3．13和图3．14所示。第20页，共134页，2023年，2月20日，星期二

经过三角保持器及滞后三角保持器再现后的函数如图3．15和图3．16所示。第21页，共134页，2023年，2月20日，星期二

从上面的保持器特性可以看出，实际保持器与理想保持器的特性总是有差别的，所以要想使保持器引起的失真足够小，采样频率就要足够高。也就是说，在仿真计算时，为了使结果准确，计算步距就得足够小。这样势必要增加计算时间。为了使计算速度较快又不使误差过大，通常在保持器后加一补偿环节(也可加在保持器前)。保持器再现被采样的连续信号时，这些信号一般都有相位移，而且再现信号同其被采样的连续信号相比，都有所衰减。所以通常采用超前装置进行补偿。例如，当采用零阶保持器再现信号时，再现后的信号要比被采样的连续信号平均滞后T/2。故应采用超前半个周期的补偿(即取c＝eTs/2)去抵消零阶再现过程引入的滞后影响。在仿真中所采用的补偿器的数学表达式形式一般为

C=λerTs(3-11)第22页，共134页，2023年，2月20日，星期二3．2离散系统差分方程的求取在上一节中，已经叙述了怎样由一个连续系统求出它的离散相似系统。有了这个离散相似系统，就可以求出连续系统的离散化数学模型。离散化数学模型是用差分方程表示的，它的求解方法如下。设线性定常系统的状态方程描述如式(3-12)所示。

式中A、B均为常数阵。推导可得：第23页，共134页，2023年，2月20日，星期二第24页，共134页，2023年，2月20日，星期二第25页，共134页，2023年，2月20日，星期二第26页，共134页，2023年，2月20日，星期二系统的离散解：

(3-23)

在推导上式的过程中未作任何近似的假设，该式是一种精确的采样值计算公式。但是，当U(τ)是一个复杂的函数时，该式右端的积分是难以求得的。由于该积分的积分区间仅为T，当T较小时，一般来说U(τ)在这个积分区间的变化是不大的。因此，可以加入采样及再现环节，以使U(τ)在积分区间内为一个简单的特殊函数，从而使该积分计算容易进行。通常使用下面三种保持器作为再现环节。

第27页，共134页，2023年，2月20日，星期二1.零阶保持器第28页，共134页，2023年，2月20日，星期二2.一阶保持器第29页，共134页，2023年，2月20日，星期二3.三角保持器用与上述相同的方法，可以求出加三角保持器和滞后三角保持器时的系统差分方程。用各种保持器时的差分方程，详见表3．1。第30页，共134页，2023年，2月20日，星期二第31页，共134页，2023年，2月20日，星期二[例3．1]己知某控制系统框图如图3．17所示，求该系统的仿真模型，即差分方程。第32页，共134页，2023年，2月20日，星期二第33页，共134页，2023年，2月20日，星期二第34页，共134页，2023年，2月20日，星期二第35页，共134页，2023年，2月20日，星期二第36页，共134页，2023年，2月20日，星期二第37页，共134页，2023年，2月20日，星期二第38页，共134页，2023年，2月20日，星期二

如果把采样开关保持器加在系统入口R处，则得到的离散相似系统框图如图3．19所示。第39页，共134页，2023年，2月20日，星期二第40页，共134页，2023年，2月20日，星期二第41页，共134页，2023年，2月20日，星期二于是得到零阶保持器时的差分方程：第42页，共134页，2023年，2月20日，星期二

由此可见，采样器及保持器的位置不同，得到的差分方程也不相同。但应注意，无论离散-再现环节加到哪里，被离散再现的信号都应是状态方程中的输入量。

从保持器的定义式可以看出，零阶保持器能无失真地再现阶跃输入信号，即当输入信号为阶跃函数时，导出的差分方程是精确的。而三角保持器能无失真地再现斜坡输入信号。第43页，共134页，2023年，2月20日，星期二第44页，共134页，2023年，2月20日，星期二第45页，共134页，2023年，2月20日，星期二[例3．2]求如图3．21所示系统的差分方程。第46页，共134页，2023年，2月20日，星期二第47页，共134页，2023年，2月20日，星期二第48页，共134页，2023年，2月20日，星期二第49页，共134页，2023年，2月20日，星期二第50页，共134页，2023年，2月20日，星期二第51页，共134页，2023年，2月20日，星期二第52页，共134页，2023年，2月20日，星期二第53页，共134页，2023年，2月20日，星期二第54页，共134页，2023年，2月20日，星期二

对于许多控制系统，都可以把系统分解成由积分和惯性环节组成系统。因此，可以事先求出这两个环节的差分方程的通用式，以后就不需要每次求解差分方程了。对于积分环节，如图3.23所示。第55页，共134页，2023年，2月20日，星期二第56页，共134页，2023年，2月20日，星期二

对于惯性环节，如图3.24所示。第57页，共134页，2023年，2月20日，星期二第58页，共134页，2023年，2月20日，星期二第59页，共134页，2023年，2月20日，星期二3．3连续系统数字仿真程序设计前两节讲述了连续系统离散化及其差分方程的求取。有了系统的差分方程，就可以在计算机上编制仿真程序了。随系统性质的不同，对仿真程序要求也不同，一般的要求是计算速度快、精度高、使用方便、通用性强等。但这些要求往往是相互矛盾的，所以具体到某一问题时，应根据其特性突出某一要求而牺牲另外一些要求。

第60页，共134页，2023年，2月20日，星期二

在一般情况下，仿真程序由以下几个基本模块构成：

(1)初始化程序块；

(2)输入参数程序块；

(3)主运行程序块；

(4)输出仿真结果程序块。仿真程序流程图如图3．25所示。下面分别讨论这三个程序块。第61页，共134页，2023年，2月20日，星期二第62页，共134页，2023年，2月20日，星期二3．3．1初始化程序块这个功能块主要是对程序中所用到的变量、数组等进行定义，并赋以初值。在通用仿真程序里这个程序完成被仿真系统的结构组态。

这个程序块的内容随使用的程序设计语言的不同而不同，没有统一的格式。本书所有的程序均使用C，matlab和VisualC++6.0作为程序设计语言。第63页，共134页，2023年，2月20日，星期二

自20世纪90年代以来，随着多媒体技术和图像技术的蓬勃发展，可视化技术得到了广泛重视，越来越多的计算机专业人员和非专业人员都开始研究并应用可视化技术。一般讲，可视化技术包含两个方面的含义：一是软件开发阶段的可视化，即可视化编程，它使编程工作成为一件轻松愉快、饶有趣味的工作；二是通过可视化窗口将不改变的参数输入给计算机。第64页，共134页，2023年，2月20日，星期二

在通用仿真程序中，可视化输入参数程序块是非常复杂的，该程序块直接关系到人-机交互的方便性和程序的通用性。一般人们会花很大的精力来设计输入参数程序块。在后面的章节中将逐步介绍该程序块的设计方法。第65页，共134页，2023年，2月20日，星期二3．3．2输出仿真结果程序块

该程序块输出仿真结果，它们可以是状态变量、中间变量、输出变量的仿真结果。输出的形式总是数据表格或曲线的形式。下面介绍这两种输出程序的设计方法。在以后的程序中只写出调用这两个输出仿真结果函数的语句，而不再列写程序清单。第66页，共134页，2023年，2月20日，星期二1．打印数据函数PrintData()

一般要求打印数据函数能打印各变量的仿真数据及各数据对应的时间。由于仿真点数较多，输出数据量较大，用户不一定需要观察每一步的数据。因此，往往时隔一定的间隔打印一个数据。下面的函数可打印1～4个变量(变量个数用VN表示)的仿真数据及其对应的时间，每隔IN点打印一个数据。并要求在调用的程序里把需输出的变量值存放在数组Output[][]里，其中第一维下标表示第几个变量，第二维下标表示该变量的第几个仿真数据。LP、DT的意义表示计算步数和采样周期。第67页，共134页，2023年，2月20日，星期二第68页，共134页，2023年，2月20日，星期二

2．显示曲线函数Dispcurve()

该程序可以在显示器上同时显示1～8条响应曲线。横轴为响应时间，纵轴为响应数值。要求在调用的程序里把需要输出的变量值存放在数组Output[][]里，其中，第一维下标表示第几个变量，第二维下标表示该变量的第几个仿真数据。VN表示要显示变量的个数，LP、DT的意义表示计算步数和采样周期。第69页，共134页，2023年，2月20日，星期二第70页，共134页，2023年，2月20日，星期二第71页，共134页，2023年，2月20日，星期二第72页，共134页，2023年，2月20日，星期二第73页，共134页，2023年，2月20日，星期二第74页，共134页，2023年，2月20日，星期二3．3．3主运行程序块

这个程序块用来求解被仿真系统的差分方程。所选用的仿真算法不同，得到的差分方程也不同，仿真精度也不一样。不管怎样，这个程序要忠实于原差分方程，它不能改变原差分方程的意义，对于初编程序者来说在这方面是很容易出错的。从上一节可以看到，系统的仿真模型都是差分方程的形式。因为差分方程最容易用数字计算机求得数值解。差分方程的一般形式如式(2-8)所示。如果改变一下符号并写成递推计算的形式:第75页，共134页，2023年，2月20日，星期二

由式(3-67)可知，在求解此差分方程时，要用到计算时刻(k+1)T以前若干个采样时刻的输出值和输入值。这可以在内存中设置若干个存储单元，将这些数据存储起来，以便在计算时使用。

下面讨论输出变量的处理方法。对于式(3-67)所描述的系统，差分方程阶次为n。因此，需要在内存中设置n个单元，用以存放计算时刻(k+1)T以前的n个采样时刻的输出量。这些存储单元的安排如图3．26(a)所示。第76页，共134页，2023年，2月20日，星期二第77页，共134页，2023年，2月20日，星期二第78页，共134页，2023年，2月20日，星期二第79页，共134页，2023年，2月20日，星期二[例3．3]设计例3．2所述的控制系统仿真程序中的主运行程序块。取R为单位阶跃函数，状态变量初值为零。

第80页，共134页，2023年，2月20日，星期二第81页，共134页，2023年，2月20日，星期二

那么，需要多少步才能求出系统的响应呢?这取决于实际问题。一般计算到系统稳态以后就可以了，因为稳态后的响应我们己知道，再计算就没什么意义了。因此仿真之前应输入一个仿真时间，这个仿真时间一般是系统达到稳态的时间(在以后的程序中用变量ST来表示)。如果再给出计算步距(即采样周期，在以后的程序中用变量DT表示)，则可得到需要计算的步数

LP=ST/DT(3-68)因此，运行程序块要完成一个循环计算，其结构形式为

for(i=1；i<=LP；i++){/*差分方程表达式*/}第82页，共134页，2023年，2月20日，星期二

式(3-60)是一阶差分方程组，根据前面的讨论，每个变量用两个存储单元就够了。如果在程序中，用变量x11存放x1(k+1)，x10存放x1(k)，x21存放x2(k+1)，x20存放x2(k)，x31存放x3(k+1)，x30存放x3(k)，x41存放x4(k+1)，x40存放x4(k)，e(0)存放e(k)，用数组Outputy[][]存放输出y的各步值，则可得到加零阶保持器时的主运行程序块：第83页，共134页，2023年，2月20日，星期二第84页，共134页，2023年，2月20日，星期二第85页，共134页，2023年，2月20日，星期二第86页，共134页，2023年，2月20日，星期二第87页，共134页，2023年，2月20日，星期二第88页，共134页，2023年，2月20日，星期二第89页，共134页，2023年，2月20日，星期二第90页，共134页，2023年，2月20日，星期二

从上面的仿真程序设计过程中可以看出，由于串行计算的原因，三角保持器不是放在哪里都可以实现的。实际上，不只是三角保持器，有超前作用的其他保持器或加超前补偿的零阶、一阶保持器也不是放在哪里都可以实现的，因为在有超前作用的保持器下，其差分方程的输入项可能会需要k+1时刻的值。这样，在系统的反馈支路与主支路求和点后的第一个环节入口处，不能加有超前作用的离散-再现环节。第91页，共134页，2023年，2月20日，星期二第92页，共134页，2023年，2月20日，星期二第93页，共134页，2023年，2月20日，星期二3．4典型非线性环节的仿真程序设计

一般工业系统均存在非线性问题。例如，调节仪表的输出可能存在饱和非线性，执行机构可能存在齿轮间隙非线性，还有些系统采用非线性调节器。因此，应当考虑非线性系统的仿真问题。

实际上，前面介绍的按环节离散化的仿真方法，可以很容易地推广到具有典型非线性环节的非线性系统的仿真。只要事先编好这些典型非线性环节的仿真子程序，让信号通过它再输入到线性环节就行了。为了以后编程的方便，下面给出几种典型非线性环节的仿真子程序。

第94页，共134页，2023年，2月20日，星期二1．继电器继电器的非线性特性如图3．29所示。图中y0为继电器常数，它的特性可用如下的数学表主式描述：当u＜0时，y＝-y0

当u≥0时，y＝y0。第95页，共134页，2023年，2月20日，星期二2．有不灵敏区的继电器具有不灵敏区的继电器特性如图3．30所示，图中y0、c为常数。其数学表达式为当u≥c时，y＝y0

当|u|＜c时，y＝0；当u≤-c时，y＝-y0第96页，共134页，2023年，2月20日，星期二3．限幅器限幅器的特性如图3．31所示，其数学表达式为当u＞c时，y=c

当|u|≤c时，y=u

当u＜-c时，y=-c式中：c为常数。第97页，共134页，2023年，2月20日，星期二4．不灵敏区不灵敏区特性如图3．32所示，图中c为常数，其数学表达式为当u＞c时，y=u-c

当|u|≤c时，y=0

当u＜-c时，y=u+c第98页，共134页，2023年，2月20日，星期二5．齿轮间隙非线性第99页，共134页，2023年，2月20日，星期二第100页，共134页，2023年，2月20日，星期二6．具有死区和滞环的继电非线性第101页，共134页，2023年，2月20日，星期二[例3．4]某发电机励磁系统框图如图3．35所示，试求URFF为单位阶跃函数时系统的响应。第102页，共134页，2023年，2月20日，星期二根据图3.36可得到系统的差分方程如下：第103页，共134页，2023年，2月20日，星期二第104页，共134页，2023年，2月20日，星期二仿真程序如下：#include<stdio.h>/*函数头文件声明*/#include<stdlib.h>#include<math.h>#include<graphics.h>floatoutputy[1][1501],ST,DT;intLP,VN;voiddispcurve();/*显示曲线子程序*/第105页，共134页，2023年，2月20日，星期二main()/*主运行程序块*/{floatx11,x21,x31,x41,x10,x20,x30,x40;floate1,e2,e3;/*定义初始变量和参数*/floatR;floatURff,Ue,UR,UR1;/*定义输入*/inti;floata1,a2,a3,a4;/*中间变量参数*/x10=0;x20=0;x30=0;x40=0;/*参数赋初值*/ST=10;/*仿真时间*/DT=0.02;/*仿真步距*/LP=ST/DT;VN=1;R=1;URff=1;第106页，共134页，2023年，2月20日，星期二outputy[0][0]=0;/*存储输出量的初值*/a1=exp(-dt/0.1);a2=exp(dt*0.1);a3=exp(-dt);a4=exp(-dt/0.05);第107页，共134页，2023年，2月20日，星期二

for(i=0;i<LP;i++)/*计算各差分方程*/{Ue=URff-x40;x11=a1*x10+40*(1-a1)*Ue;if(x11>30)UR=30;if(abs(x11)<=30)UR=x11;if(x11<-30)UR=-30;UR1=UR-0.01*x20;x21=a2*x20+(a2-1)*UR1/0.05;x31=a3*x30+(1-a3)*x21;x41=a4*x40+0.05*(1-a4)*x31;outputy[0][i]=x31;x10=x11;x20=x21;x30=x31;x40=x41;}dispcurve();/*调用画图子程序*/}第108页，共134页，2023年，2月20日，星期二第109页，共134页，2023年，2月20日，星期二3．5关于采样周期(计算步距)和仿真时间的选择

计算步距和仿真时间是两个重要的仿真参数，特别是计算步距。如果选择得不恰当，就可能造成较大的计算误差，甚至可以使一个本来稳定的系统歪曲成一个不稳定的系统。计算步距不仅和被仿真的系统有关，还和仿真算法、精度要求等因素有关。因此，要在仿真计算之前准确地选好这两个参数是件不容易的事情。第110页，共134页，2023年，2月20日，星期二

在3．1节中曾提及过，为了使被采样的信号无失真的再现，必须满足式中:为最低采样频率；为被再现信号的频带限。

但是，在仿真中所遇到的大多数被再现信号是没有频带限的，所以一般取采样频率是再现信号主要频带中的最高频率的5～10倍。

第111页，共134页，2023年，2月20日，星期二

至于主要频带中的最高频率是没有确切的定义的，但对于像这样简单的低通滤波器(惯性环节)，其频谱的主要频带可以认为大约是截止角频率()的10倍左右所以可以选择即选择计算步距(采样周期)(3-70)

第112页，共134页，2023年，2月20日，星期二

对于复杂环节仍可取主要频带是开环频率特性的剪切频率的10倍，如仍用表示环节的剪切频率，则计算步距仍为式(3-70)。若系统中有几个小闭环，则应取最快的小闭环频率特性的剪切频率。

随着系统环节数目的增加，加入的离散-再现环节就可能增多。这时计算步距应相应地减小。

由于求系统开环频率特性比较麻烦，下面针对控制系统工程设计的需要，提供估计这两个仿真参数的简单方法。第113页，共134页，2023年，2月20日，星期二3．5．1计算步距的估计

计算步距可用式(3-71)估计

(3-71)式中：DT为计算步距；n为被控对象传递函数的阶次；T为被控对象传递函数的时间常数。第114页，共134页，2023年，2月20日，星期二

如果被控对象有若干个，则应以其中nT最小的为准。一般使用者按上述区间选择一个适当的计算步距，其仿真结果是令人满意的。一般来说，计算步距选择的越小，计算精度就越高(但步距太小，也可能由于计算机的舍入误差占了主导地位而降低了计算精度)，耗费的计算时间就越长；反之，就会导致计算精度差，计算时间短。但对于一般工程设计来说，最关心的是系统的稳定性，能够大致观察过渡过程曲线就可以了，对计算数据精确度的要求并不太高，所以只需采用试算的方法求计算步距即可。如果过渡过程曲线的第一个半波含有20到30个计算点，则说明计算步距选择是恰当的。否则，可适当增大或减小计算步距，再试算。第115页，共134页，2023年，2月20日，星期二[例3．5]某火力发电厂汽包水位控制系统采用单级三冲量调节方式，系统框图如图3．38所示。对此系统进行仿真，绘制蒸汽量D扰动时，水位H的变化曲线。第116页，共134页，2023年，2月20日，星期二第117页，共134页，2023年，2月20日，星期二第