第六章 离散时间和连续时间模型的仿真
§1 状态变量
6.1.1 状态变量的基本概念
1) 状态变量集
计算机仿真中必须搞清楚实体相互关系的规则。计算机记录描述变量的过去值，根据相互关系规则，可计算描述变量的未来值。
状态变量集是所有描述变量的一个子集，只要知道这些变量的现在值和输入变量值，就可计算模型的所有描述变量未来值。 2）模型完全描述
完全描述模型：
假设模型具有描述变量α1, α2, , αn ，如果在任一时间t ，变量α1的值为y 1，变量α2的值为y 2，„，若实体的相互关系规则对任一未来时
' '
间 t′（大于 t）确定了值y 1' , y 2的唯一集，那么该模型是完全, , y n
描述的。
模型完全描述的充要条件：
如果各描述变量的各个值只在任一时间t 唯一确定所有这些变量在任一未来时间t ′的值，就说描述变量集的某个子集是状态变量集。
如果模型是完全描述的，α1, α2, , αn 或它的真子集便是状态变量集。
模型是完全描述的充要条件是该模型的描述变量中存在状态变量集 例：二辆汽车面对而驶，V 1、V 2
6.1.2 状态变量的仿真性质
1） 程序预置
' '
假设程序给出计算t ′时的y 1' , y 2的任务。则仅需预置, , y n
（也即是初始化）那些与状态变量有关的存储单元。 2） 重复操作
假设给定t 时的y 1, y 2, , y n 值之后，因为丢失了第一次仿
' ' 真操作的记录，要重复计算t ′时的y 1' , y 2值，只要与状, , y n
态变量有关的单元，预置y 1, y 2, , y n 的相同值，则在不同计算机和不同时间作两次操作，结果仍然相同。 3） 程序中断和重新起动
' ' 设计算t ′时的y 1' , y 2值之后，安排中断程序。在, , y n
某时间之后可以重新起动。 4） 程序恢复
假设计算机在执行程序时发生事故，修复正常时，重新预置肯定将最终产生相同结果，但比从中断点重新起动要花费更多的时间。
4）离散时间仿真的定义
i
对于时间t 1，„ t i „，给出t i 的状态值y 1i y m ，由程序根据分量i +1相互关系能计算t i+1的唯一描述值y 1i +1 y m ，也就是对于任何时间对
偶（t 1，t i+1）均成立，称集{ t 1，t 2„}为计算时刻，若时间是步长h
的逐次倍数t i+1- t i =h，整个仿真称离散时间仿真。又假设分量相互关系规则不依赖于时间，仅仅与状态值y 1„y m 有关，模型就称时不变模型。
§6.2离散时间模型仿真
6.2.1时不变离散时间模型的仿真过程
1）仿真过程
给定t 时的状态值y 1„y m ，求t 1的描述变量值的问题。设t=tM 和t 1=tM+N，计算时刻集t M ，t M+1，„t N 在t 和t 1之间，仿真过程： 步1 预置状态变量α1 αm 的值分别为y 1，y 2，„y m 步2 预置时标为t M
步3 根据相互关系规则和状态变量现值，产生状态变量的新内容，并使其它剩余描述变量产生新的内容 步4 t=t+h，推进时标
步5 t ≥t m +Nh ，停止计算；否则，返回步3 2）仿真基本特性
（1）计算模型描述值的采样时刻，由t ，t+h，„„序列组成 （2）迭代次数(t 1-t ) /h ，h 越小，迭代次数越少 （3）步3体现模型的相互关系规则，是关键
6.2.2 离散时间模型的形式规范
1．范式
下面讨论离散时间仿真步3
假设共 n 个变量，其中m 个为状态变量，并假设无输入变量。
它的输入:输出:描述值
f ,
111
f (y 1, y 2, , y m ) =(y 1, , y 1m , y m +1, y n )
把f 分成二个部分
11δ(y 1, y m ) =(y 1, , y m )
1
1
1m
λ((y 1, y m ), (y , , y ) )=(y , , y , , y )
11
1m
1n
即
y 1,
11, „,y n 1
若λ不依赖y 1„,y m ，则
1
δ(u 1, y m ) =(y 1, y 1m ) 1111λ(y 1, , y m ) =(y 1, y n )
f (y 1, , y m ) =λ(δ(y 1, , y m ) )
映射δ称为状态转移函数—它取一列时间t i 的模型状态变量值，并产生一列时间t i+1的模型状态变量值。上面的简式称范式。
2．离散时间系统的组成和形式化描述
状态、输出、转移函数
假设DES.var 是描述变量集，state.var 和output.var 分别为状态变量和输出变量集，state.var 的范围集. 用RANGRE. βi 表示，这样状态集可用下式笛卡尔积表示
STATES ={(y β1, , y βi , ) y βi ∈RANGE . βi , βi ∈state . var }
=R A N G ⋅E β1⨯R A N G ⋅E β2 =⨯β∈s t a . v t a e R ⋅E β r A N G
同样，可用集合论来描述输出集OUTPUTS
对于时不变和离散时间步长为h 的离散时间仿真模型，就可给出它的形式描述：
⎧δh :STATES →STATES 和 λ:STATES →OUTPUTS
⎨
S (STATE ) ∈STATES λ(STATE ) ∈OUTPUTS ⎩h
若STATE 是t i 时刻的模型状态变量，δh (STATE)则t i +h时刻状态，λ(STATE)则是t i 的输出。
上述四元组是离散时间系统规范 例：线性同余发生器的离散时间系统规范
M=,Q为状态集[0，m]，Y 是输出集[0，1]
δ是转移函数，δ：Q →Q ，q ∈Q, δ(q)等于aq+c关于m 的模，λ是输出函数，λ：Q → Y, λ(q)=q/m
6.2.3 离散时间模型的结构与行为
假设仿真时间区间为[t, t'], 把它称为观察区间，计算时刻序列为t M , „, t M+Nh, 设序列q M , „, q M+N, 分别为计算时刻t M , „, t M+Nh的状态，同样输出变量序列λ(qM ), „, λ (qM+N) 把这二个序列分别称为状态轨迹和输出轨迹。 1）行为
系统的行为包括模型的状态行为和模型的轨迹行为 状态行为是所有状态轨迹的集
输出行为是所有输出轨迹的集
2) 结构
状态转移函数：q M +i =δh (q M +i -1)
i =1, , N
6.2.4 非自治离散时间模型
DES ⋅var =INPUT ⋅var non.INPUT. var
STATE. var ⊂non.INPUT. var OUTPVT. var ⊂non.INPUT. var
非自治时不变离散时间系统规范可用如下五元组形式表示，
其中INPUTS ，STATES ，OUTPUTS 分别是INPUT.V ARIABLES ， STATE.V ARIABLES 和OUTPUT.V ARIABLES 的范围集叉积的子集；δh 为状态转移函数
δh ：STATES ⨯ INPUTS → STATES
而λ是输出函数
λ：STATES ⨯ INPUTS → OUTPUTS
具有输入变量系统的行为
⎧状态行为:它是输入-状态轨迹
模型的行为⎨
⎩输入输出行为:它是输入-输出轨迹
§6.3 连续时间模型仿真
6.3.1 微分方程系统规范
特点：微分方程系统不能直接规定下次状态，而是根据所提供的状态变化信息来计算下次状态值。
微分方程系统的规范的五元组形式表示：
，
其中INPUTS ，STA TES ，OUTPUTS 分别是INPUT ·V ARIABLES ， STATE ·V ARIABLES 和OUTPUT ·V ARIABLES 的范围集叉积的子集；f 为导数函数
f ：STATES ⨯ INPUTS → STATES
而λ是输出函数
λ：STATES ⨯ INPUTS → OUTPUTS
导数函数f 是状态变量和输入变量的函数，可用下式表示
dq (t )
=f (q (t ), x (t )), q(0) =q dt
上式说明它们的解依赖于给定的初始值（ q(0) =q ），并在每个时刻必须满足微分方程。
微分方程系统的关键是导数函数。设模型有输入变量序列X 1, X 2..., X m 和输出变量序列Y 1, Y 2..., Y n ，导数函数可用以下一阶微分方程组表示
dY 1
=f 1(Y 1, Y 2..., Y n , X 1, X 2..., X m ) dt dY n
=f n (Y 1, Y 2..., Y n , X 1, X 2..., X m ) dt
可以把每个微分函数分解为积分器和用来计算积分器输入的函
数，把积分器的输出变量来构成模型的状态变量（详见下一节）。积分器是构成微分方程说明系统的基本环节，它可表示为Y=INTGRL（IY ，YRT ），它建立起它的输入YRT 和输出Y （用简化符号）所假定的值之间的关系
dY
=YRT (t ) dt
6.3.2 积分法
欧拉法的基本概念是
dY (t ) Y (t +h ) -Y (t )
=lim =X (t ) dt h h →0
对于足够小的h ，有
Y (t +h ) =Y (t ) +h ∙X (t )
设定较小的步长h ，对于给定运行长度，则会要求有较长的计算时间。但是，从计算精度的角度来说，h 越小计算精度越高。
积分法思想：
积分法利用输出和导数的估计过去值和（或）未来值来致力于更好地估计现在值。对于积分器Y=INTGRL（IY ；X ），计算时间t i 的Y 值的依据可能包括计算Y 和X 在先前计算时刻t i -1, t i -2, 和（或）后来时刻t i +1, t i +2, 的各个值。输出Y 和导数X 的各个值是相互依存的，即积分器本身只产生Y 对X 的依赖关系，而通过其它函数，可形成积分器输出反馈，这使X 对Y 产生依赖关系。 误差传播
（1）每步计算近似值, 因为即使X 和Y 在是正确的，由于在连
续区间（t ，t+h）X 和Y 的值是变化，程序不能在这区间使用这些在时间t 的值，来正确地计算Y 在时间(t＋h ）的值。
（2）系统传播前误差的影响而积累起来，Y 有误差，Y 与X 是依存关系，导致X 误差，再导致Y 误差。 因果法:
仅使用Y 和X 在前计算时刻t i-1,t i-2, ⋯, 的值和时间t i 的X 值，去计算t i 的Y 值，叫做因果法。
因果法的阶数：若Y i 的f (Xi-d , ⋯, Xi ; Yi-d , ⋯, Yi-1) 表明f 是自变量的线性组合，它的阶数为d 。
例: 亚当斯法 Y t+h = Yt +h (3Xt -X t-h ) 为二阶因果法
对于d 阶因果法，仿真时应保存d 个模型导数和状态（积分器的输入及输出）变量的过去值。
因果积分法对应的状态转移函数可用以下语句表示
Y=INT•METHOD （P1Y ，P2Y ，„，PdY ，P1X, „，PdX ，X ） 其相关的状态变量、状态转移过程和输出序列可见表6.1。
表6.1 因果积分法
非因果法:
非因果法是不但考虑导数和输出的过去值，而且使用Y 和X 在时间t i+1, ti+2, ⋯, 的估计值。
两次试探性估计现在值：（1）计算未来值的初始“预算”阶级
（2）计算最终值时的“校正”阶级
Y i =g (x i -d , , x i -1, , i +1, , i +e , Y i -d , , Y i -1, i , i +1, , i +e )
(x i -d , y i -d ), , (x i -1, y i -1) 是已计算过的过去值对偶Y i 是Y 在t i 的最终估计(i , i ), , (i +e , i +e ) 是预算现在和未来值，值。
简单的预算一校正法如下：
t =Y t -h +hX t -h
h
Y t =Y t -h +(t +X t -h )
2
非因果积分法的实现可描述如下：利用因果法f 函数，计算e+1个时刻t i , t i +1,..., t i +e 积分器的导数和输出对偶，并把获得的对偶
(, Y , X
i
i
i +1
, Y i +1, X i +e , Y i +e 存储起来；然后，把函数g 应用于这些值
))
和d 个过去值(X i -d , Y i -d ), (X i -d +1, Y i -d +1), , (X i -1, Y i -1) ，求得输出变量Y 在时间t i 的最终估计值Y i 。
§6.4离散时间和连续时间仿真模型的描述
6.4.1 污染模型
POL(t+h)=POL(t)+[POLG(t)-POLA(t)]h POLG(t)=P(t)*POLCM(t) POLCM(t)=POLCMT(CIR(t)) CIR(t)=
CI (t )
P (t )
POLA(t)=POL(t)⨯POLR(t) POLA(t)=POLRT(POL(t))
⎡⎤⎛CI (t ) ⎫
Pol (t +h ) =Pol (t ) +⎢P (t ) ⨯POLCMT -POL (t ) ⨯POLRT (POL (t )) ⎪⎥⨯h P (t ) ⎪
⎝⎭⎣⎦
对于状态变量POL
δPOL (pol , ci , p ) =pol +[p ⨯POLCMT () -pol ⨯POLRT (pol )]h
ci p
1．瞬时函数：根据时间t 的各输入值确定模型在时间t 的输出值的函数。
2．记忆或非瞬时函数：用微分方程说明，它可用连续时间和离散时间来进行仿真。
1． 已标识了瞬时函数及记忆函数的模型网络描述，每个模型网络图可容易地转换成某些仿真语言的语句。
2． 由模型网络描述，也可容易地明确状态转移函数和输出函数，离散时间模型的仿真就可按其仿真过程进行。
3． 仿真过程中，按某个固定顺序搜索状态变量，用有关的局部转移函数计算每个状态变量来获得下次状态值。由于某些瞬时函数可能会直接影响两个以上的转移函数或输出函数，会不可避免地出现不止一次地计算这个瞬时函数。如果能找出各个瞬时函数的计算次序，就可以这一问题。为此，给出一种非常类似于仿真语言的模型描述语言，它可描述编排模型描述语句的过程。
6.4.2 模型描述语言
模型描述语言是一种类似仿真语言的用来描述编排模型描述语句的过程。模型描述语言中的三种语句： 1) 瞬时函数
Y 1, Y2, ⋯, Ym = Instant. Func (X1,X 2, ⋯ Xn )
在t 时刻，能给定变量X 1, ⋯, X n 的值，输出Y 在时刻t 的值直接可由上式获得。 2) 输入时间函数
X 1, X2, ⋯, Xm = Time. Func
它是用来产生模型的输入轨迹，同瞬时函数不同，它仅与时间t 有关
例X 1, X2=sin, cos
X 1, X2的轨迹为 (sin(t), cos(t))
3) 记忆函数
Y 1, Y2, ⋯, Ym = Mem. Func (Q1, Q2, ⋯, Qn ; X1, ⋯, Xp )
它具有状态变量Q 1, ⋯, Qn , 其初值必须设定，本质上它包括输入轨迹、状态转移函数和输出函数三个部分。
对于记忆函数，用离散时间和连续时间仿真，其处理方式不同，它们分别可用一个延迟元件和积分器来描述
记忆函数的两种基本环节：
对于微分方程描述的系统，记忆函数的基本环节是积分器。积分器可用语句Y =INTGRL (IY , YRT ) 来表示，它的输出Y 是输入YRT 加上初始状态IY 的时间积分。它是构成微分方程描述的系统的基本环节。可把积分器转化成离散时间仿真形式，但其对应的转移函数只可能是近似的。
对于离散时间模型，记忆函数的基本环节是迟延元件。迟延元件可用语句 Y =DELAY (IY , YP ) 来表示。对于离散时间的一个DELAY(迟延) ，其输出轨迹Y （·）比输入轨迹YP(·) 滞后一步，它的初始值由IY 给定。DELAY 语句能够描述每一种离散时间模型。
模型描述语句与模型网络图存在一一对应关系，可从模型描述语句直接推导出模型网络图。例如，描述语句 U =PROD（Y1，Y2），表示瞬时函数PROD 有两条Y1和Y2输入线和一条U 输出
线；描述语句 Y1=DELAY（IY1,Y1P ），说明Y1也是DELAY 函数的输出线；这样，两个DELAY 方框的输出作为 PROD 的输入。 从上面的例子，发现描述序列中语句的语序与描述无关，也就是从语句序列S 1,S 2, „，S n 的任何排列中均可得到相同的网络。
考虑以下模型描述： Y2=DELAY（IY2,Y2P ）
U =PROD（Y1，Y2） Y1P=SUM（X1，U ） Y2P=SUM (X2，U)
X1，X2=SIN，COS Y1=DELAY（IY1,Y1P ）
U
图6.5 模型描述语句对应的网络图
6.4.3 模型描述语句序列分析
在描述序列中，语句的语序同描述无关，一个描述序列的任何排列均可得到相同的网络。但在模型仿真中，语句被翻译计算机的指令时，则必须正确次序执行这些动作。
描述语句序列必须经过分析后，才能确定它是否有效及正确执行次序。
1．无效语句的检查
检查是否有任一变量在两个不同语句的左边出现，若存在则是无效语句
2．无记忆和输入时间函数时的函数循环问题判别
（1）领先关系
对于Y 1, Y2, ⋯, Ym = Instant. Func (X1,X 2, ⋯ Xn )
若U 是X i 中的一个变量，V 是Y j 中的一个变量，则U 领先于V 。若U 领先V ，在计算t 时刻的V 值之前，必须要知道时间t 的U 值。 （2）领先关系的闭包传递
若有一变量序列W o , ⋯, Wn , Wo = U, Wn = V, 各个W i 领先W i+1, 则U 领先*V，
U 领先*V表示在网络图中存在一个沿箭头从U 到V 的路径 例：V = Prod (Y1, Y2)
Y 1=Sum (X1, V) Y 2=Sum (X2
,V)
结论:有存在V 领先*V，模型描述语句构成循环，存在无效描述
3．变量的次序排列
1）变量的级别
（a ）若对某个变量V ，变量U 领先于V ，且无变量W ，它使W 领先U ，那么变量U 的级别Lev （U ）=0
（b ）任一（其他）变量的级别是网络中从级别为0的变量到该变量最长路径的距离。
若U 领先*V，级别Lev （U ）=0，n 为最大数，U 领先W 1领先W 2⋯W n-1领先V ，级别Lev （V ）=n
举例
2）语句的级别（函数的级别）
对于语句S i : Y1, ⋯, Ym = Instant.func (X1, X2, ⋯, Xn )
语句S i 的级别是X j 的最大级别，即级别(Si )=max{级别(Xj ), j=1, n} 求Prod 和Sum 的级别。
根据这些排序层次，就可自动构造层次关系清楚的模型网络图。在构造模型网络图时，把相同排序层次的所有变量和语句置在一条线上。
对于例6.7的模型描述语句：
LEV[U=PROD（Y1，Y2）]=max｛LEV （Y1），LEV （Y2）｝=0
LEV [Y1P=SUM（X1，U ）]= LEV （U ）=1 LEV [Y2P=SUM（X2，U ）]= LEV（U ）=1
排序层次
图6.7 模型网络图
对于瞬时函数，排列次序按它的级别排列S 1, S2, ⋯, Sn , 级别(Si ) ≤级别(Si+1) 。一般来说，时间函数语句排在最前，记忆型函数则排列最后。最终序列形式为：
────
┆ TIME ·
FUNOT ────
┆ 排序层次为0 ────
┆ INSTANT ·FUNCT ────
┆ 排序层次为M ────
┆ MEM ·FUNCT ────
其中M 是语句的排序层次的最大值。 X 1, X2=Sin, Cos U = Prod (Y1, Y2) Y 1P = Sum (X1, U) Y 2P= Sum (X2, U) Y 1 = Delay (IY1, Y1P)
Y 2=Delay (IY2, Y2P),
6.4.4 记忆函数仿真
1、离散时间模型
当模型描述仅包含 DELAY 记忆函数语句时，可选择 DELAY 的输出变量来构造模型的状态变量。因而每个 Y=DELAY（IY ，X ）语句可解释如下：在预置阶段，设置 Y 等于IY ；在状态转移阶段，设置Y 在时间t i +h的值为X 在时间t i 的值。
对于例6.6，其模型描述语句对应的仿真的基本过程为： 1） 预置时标T 为要求的初始时间t ； 2） 置 Y1为IYl 、 Y2为IY2; 3） 置 X1，X2为SIN(T）,COS(T); 4） 置U 为PROD(Y1,Y2）；
5） 置Y1P 为 SUM （X1,U ）、Y2P 为SUM （X2,U ）； 6） 置Y1为 Y1P 、Y2为 Y2P ； 7） 推进仿真时标，置 T 为T+h； 8） 若T 小于终止时间，转向3）； 9） 停机。
常态形序列，则可以每个延迟元件为中心构造描述语句，可以起到对这个模型描述序列的分解作用。 设S 1，...,S m 是最终序列，并令
Y1=DELAY（IYl ；Y1P ）
┆
Yn=DELAY（IYn ，YnP ） 是记忆函数语句。对于每个迟延元件的输入变量 YiP ，设先前变量集PRIOR ·V ARIABLES i 是计算YiP 有关的所有变量集，也就是PRIOR ·V ARIABLES i 是这样的一个变量集，它的所有变量领先*于Y i P 。现设S i1,S i2，„，S im 是计算YiP 有关的最终序列的子序列，即每个S ij 包含YiP 和（或）包含至少一个PRIOR ·V ARIABLES i 中的变量。
这样可构造常态形序列，其基本形式如下：
──
┆ ── Y1的转移函数 Y1 = DELAY (IY1; Y1P)
──
┆ ── Y2的转移函数 Y2 = DELAY (IY2; Y2P) ┆
──
┆ ── Yn 的转移函数 Yn = DELAY (IYn; YnP)
── 输出函数 ┆ 剩余语句 ──2、微分方程说明的模型
在连续时间模型中，一般用微分方程来描述系统模型，其记忆函数语句就是积分器。积分器Y=INTGRL（IY ，YRT ），是建立起它的输入YRT 和输出Y 所假定的值之间的关系
dY
YRT (t ) dt
这样，就可用积分法根据输出变量去逐次计算描述变量值，对整个导
数函数可通过一定手段构造对应的状态转移函数，用迭代方式进行系统仿真。
以下是欧拉积分法的仿真过程： 1） 预置T=to, Y=IY；
2） 计算中间变量M=Y+h·YRT ； 3） 推进时标T=T+h, Y=M；
4） T 是否小于终止时间，否，转2）； 5） 结束。
作业：对于滞后变量的伪随机数发生器，其算法的步骤：
1）设S -5，S -4，S -3，S -2，S -1是区间[0，m]中的各个整数（逻辑种子） 2）置i=0
3）选取S i-3,S i-5, 计算：S i 为aS i-3+bSi-5关于m 的模 r i= Si /m.
4）i=i+1,转入3）。
1) 给出本伪随机数产生模型的离散时间系统规范； 2) 给出本模型网络图；
3) 给出本模型描述语句和去掉记忆环节后的各变量和语句的排序层次。
第六章 离散时间和连续时间模型的仿真
§1 状态变量
6.1.1 状态变量的基本概念
1) 状态变量集
计算机仿真中必须搞清楚实体相互关系的规则。计算机记录描述变量的过去值，根据相互关系规则，可计算描述变量的未来值。
状态变量集是所有描述变量的一个子集，只要知道这些变量的现在值和输入变量值，就可计算模型的所有描述变量未来值。 2）模型完全描述
完全描述模型：
假设模型具有描述变量α1, α2, , αn ，如果在任一时间t ，变量α1的值为y 1，变量α2的值为y 2，„，若实体的相互关系规则对任一未来时
' '
间 t′（大于 t）确定了值y 1' , y 2的唯一集，那么该模型是完全, , y n
描述的。
模型完全描述的充要条件：
如果各描述变量的各个值只在任一时间t 唯一确定所有这些变量在任一未来时间t ′的值，就说描述变量集的某个子集是状态变量集。
如果模型是完全描述的，α1, α2, , αn 或它的真子集便是状态变量集。
模型是完全描述的充要条件是该模型的描述变量中存在状态变量集 例：二辆汽车面对而驶，V 1、V 2
6.1.2 状态变量的仿真性质
1） 程序预置
' '
假设程序给出计算t ′时的y 1' , y 2的任务。则仅需预置, , y n
（也即是初始化）那些与状态变量有关的存储单元。 2） 重复操作
假设给定t 时的y 1, y 2, , y n 值之后，因为丢失了第一次仿
' ' 真操作的记录，要重复计算t ′时的y 1' , y 2值，只要与状, , y n
态变量有关的单元，预置y 1, y 2, , y n 的相同值，则在不同计算机和不同时间作两次操作，结果仍然相同。 3） 程序中断和重新起动
' ' 设计算t ′时的y 1' , y 2值之后，安排中断程序。在, , y n
某时间之后可以重新起动。 4） 程序恢复
假设计算机在执行程序时发生事故，修复正常时，重新预置肯定将最终产生相同结果，但比从中断点重新起动要花费更多的时间。
4）离散时间仿真的定义
i
对于时间t 1，„ t i „，给出t i 的状态值y 1i y m ，由程序根据分量i +1相互关系能计算t i+1的唯一描述值y 1i +1 y m ，也就是对于任何时间对
偶（t 1，t i+1）均成立，称集{ t 1，t 2„}为计算时刻，若时间是步长h
的逐次倍数t i+1- t i =h，整个仿真称离散时间仿真。又假设分量相互关系规则不依赖于时间，仅仅与状态值y 1„y m 有关，模型就称时不变模型。
§6.2离散时间模型仿真
6.2.1时不变离散时间模型的仿真过程
1）仿真过程
给定t 时的状态值y 1„y m ，求t 1的描述变量值的问题。设t=tM 和t 1=tM+N，计算时刻集t M ，t M+1，„t N 在t 和t 1之间，仿真过程： 步1 预置状态变量α1 αm 的值分别为y 1，y 2，„y m 步2 预置时标为t M
步3 根据相互关系规则和状态变量现值，产生状态变量的新内容，并使其它剩余描述变量产生新的内容 步4 t=t+h，推进时标
步5 t ≥t m +Nh ，停止计算；否则，返回步3 2）仿真基本特性
（1）计算模型描述值的采样时刻，由t ，t+h，„„序列组成 （2）迭代次数(t 1-t ) /h ，h 越小，迭代次数越少 （3）步3体现模型的相互关系规则，是关键
6.2.2 离散时间模型的形式规范
1．范式
下面讨论离散时间仿真步3
假设共 n 个变量，其中m 个为状态变量，并假设无输入变量。
它的输入:输出:描述值
f ,
111
f (y 1, y 2, , y m ) =(y 1, , y 1m , y m +1, y n )
把f 分成二个部分
11δ(y 1, y m ) =(y 1, , y m )
1
1
1m
λ((y 1, y m ), (y , , y ) )=(y , , y , , y )
11
1m
1n
即
y 1,
11, „,y n 1
若λ不依赖y 1„,y m ，则
1
δ(u 1, y m ) =(y 1, y 1m ) 1111λ(y 1, , y m ) =(y 1, y n )
f (y 1, , y m ) =λ(δ(y 1, , y m ) )
映射δ称为状态转移函数—它取一列时间t i 的模型状态变量值，并产生一列时间t i+1的模型状态变量值。上面的简式称范式。
2．离散时间系统的组成和形式化描述
状态、输出、转移函数
假设DES.var 是描述变量集，state.var 和output.var 分别为状态变量和输出变量集，state.var 的范围集. 用RANGRE. βi 表示，这样状态集可用下式笛卡尔积表示
STATES ={(y β1, , y βi , ) y βi ∈RANGE . βi , βi ∈state . var }
=R A N G ⋅E β1⨯R A N G ⋅E β2 =⨯β∈s t a . v t a e R ⋅E β r A N G
同样，可用集合论来描述输出集OUTPUTS
对于时不变和离散时间步长为h 的离散时间仿真模型，就可给出它的形式描述：
⎧δh :STATES →STATES 和 λ:STATES →OUTPUTS
⎨
S (STATE ) ∈STATES λ(STATE ) ∈OUTPUTS ⎩h
若STATE 是t i 时刻的模型状态变量，δh (STATE)则t i +h时刻状态，λ(STATE)则是t i 的输出。
上述四元组是离散时间系统规范 例：线性同余发生器的离散时间系统规范
M=,Q为状态集[0，m]，Y 是输出集[0，1]
δ是转移函数，δ：Q →Q ，q ∈Q, δ(q)等于aq+c关于m 的模，λ是输出函数，λ：Q → Y, λ(q)=q/m
6.2.3 离散时间模型的结构与行为
假设仿真时间区间为[t, t'], 把它称为观察区间，计算时刻序列为t M , „, t M+Nh, 设序列q M , „, q M+N, 分别为计算时刻t M , „, t M+Nh的状态，同样输出变量序列λ(qM ), „, λ (qM+N) 把这二个序列分别称为状态轨迹和输出轨迹。 1）行为
系统的行为包括模型的状态行为和模型的轨迹行为 状态行为是所有状态轨迹的集
输出行为是所有输出轨迹的集
2) 结构
状态转移函数：q M +i =δh (q M +i -1)
i =1, , N
6.2.4 非自治离散时间模型
DES ⋅var =INPUT ⋅var non.INPUT. var
STATE. var ⊂non.INPUT. var OUTPVT. var ⊂non.INPUT. var
非自治时不变离散时间系统规范可用如下五元组形式表示，
其中INPUTS ，STATES ，OUTPUTS 分别是INPUT.V ARIABLES ， STATE.V ARIABLES 和OUTPUT.V ARIABLES 的范围集叉积的子集；δh 为状态转移函数
δh ：STATES ⨯ INPUTS → STATES
而λ是输出函数
λ：STATES ⨯ INPUTS → OUTPUTS
具有输入变量系统的行为
⎧状态行为:它是输入-状态轨迹
模型的行为⎨
⎩输入输出行为:它是输入-输出轨迹
§6.3 连续时间模型仿真
6.3.1 微分方程系统规范
特点：微分方程系统不能直接规定下次状态，而是根据所提供的状态变化信息来计算下次状态值。
微分方程系统的规范的五元组形式表示：
，
其中INPUTS ，STA TES ，OUTPUTS 分别是INPUT ·V ARIABLES ， STATE ·V ARIABLES 和OUTPUT ·V ARIABLES 的范围集叉积的子集；f 为导数函数
f ：STATES ⨯ INPUTS → STATES
而λ是输出函数
λ：STATES ⨯ INPUTS → OUTPUTS
导数函数f 是状态变量和输入变量的函数，可用下式表示
dq (t )
=f (q (t ), x (t )), q(0) =q dt
上式说明它们的解依赖于给定的初始值（ q(0) =q ），并在每个时刻必须满足微分方程。
微分方程系统的关键是导数函数。设模型有输入变量序列X 1, X 2..., X m 和输出变量序列Y 1, Y 2..., Y n ，导数函数可用以下一阶微分方程组表示
dY 1
=f 1(Y 1, Y 2..., Y n , X 1, X 2..., X m ) dt dY n
=f n (Y 1, Y 2..., Y n , X 1, X 2..., X m ) dt
可以把每个微分函数分解为积分器和用来计算积分器输入的函
数，把积分器的输出变量来构成模型的状态变量（详见下一节）。积分器是构成微分方程说明系统的基本环节，它可表示为Y=INTGRL（IY ，YRT ），它建立起它的输入YRT 和输出Y （用简化符号）所假定的值之间的关系
dY
=YRT (t ) dt
6.3.2 积分法
欧拉法的基本概念是
dY (t ) Y (t +h ) -Y (t )
=lim =X (t ) dt h h →0
对于足够小的h ，有
Y (t +h ) =Y (t ) +h ∙X (t )
设定较小的步长h ，对于给定运行长度，则会要求有较长的计算时间。但是，从计算精度的角度来说，h 越小计算精度越高。
积分法思想：
积分法利用输出和导数的估计过去值和（或）未来值来致力于更好地估计现在值。对于积分器Y=INTGRL（IY ；X ），计算时间t i 的Y 值的依据可能包括计算Y 和X 在先前计算时刻t i -1, t i -2, 和（或）后来时刻t i +1, t i +2, 的各个值。输出Y 和导数X 的各个值是相互依存的，即积分器本身只产生Y 对X 的依赖关系，而通过其它函数，可形成积分器输出反馈，这使X 对Y 产生依赖关系。 误差传播
（1）每步计算近似值, 因为即使X 和Y 在是正确的，由于在连
续区间（t ，t+h）X 和Y 的值是变化，程序不能在这区间使用这些在时间t 的值，来正确地计算Y 在时间(t＋h ）的值。
（2）系统传播前误差的影响而积累起来，Y 有误差，Y 与X 是依存关系，导致X 误差，再导致Y 误差。 因果法:
仅使用Y 和X 在前计算时刻t i-1,t i-2, ⋯, 的值和时间t i 的X 值，去计算t i 的Y 值，叫做因果法。
因果法的阶数：若Y i 的f (Xi-d , ⋯, Xi ; Yi-d , ⋯, Yi-1) 表明f 是自变量的线性组合，它的阶数为d 。
例: 亚当斯法 Y t+h = Yt +h (3Xt -X t-h ) 为二阶因果法
对于d 阶因果法，仿真时应保存d 个模型导数和状态（积分器的输入及输出）变量的过去值。
因果积分法对应的状态转移函数可用以下语句表示
Y=INT•METHOD （P1Y ，P2Y ，„，PdY ，P1X, „，PdX ，X ） 其相关的状态变量、状态转移过程和输出序列可见表6.1。
表6.1 因果积分法
非因果法:
非因果法是不但考虑导数和输出的过去值，而且使用Y 和X 在时间t i+1, ti+2, ⋯, 的估计值。
两次试探性估计现在值：（1）计算未来值的初始“预算”阶级
（2）计算最终值时的“校正”阶级
Y i =g (x i -d , , x i -1, , i +1, , i +e , Y i -d , , Y i -1, i , i +1, , i +e )
(x i -d , y i -d ), , (x i -1, y i -1) 是已计算过的过去值对偶Y i 是Y 在t i 的最终估计(i , i ), , (i +e , i +e ) 是预算现在和未来值，值。
简单的预算一校正法如下：
t =Y t -h +hX t -h
h
Y t =Y t -h +(t +X t -h )
2
非因果积分法的实现可描述如下：利用因果法f 函数，计算e+1个时刻t i , t i +1,..., t i +e 积分器的导数和输出对偶，并把获得的对偶
(, Y , X
i
i
i +1
, Y i +1, X i +e , Y i +e 存储起来；然后，把函数g 应用于这些值
))
和d 个过去值(X i -d , Y i -d ), (X i -d +1, Y i -d +1), , (X i -1, Y i -1) ，求得输出变量Y 在时间t i 的最终估计值Y i 。
§6.4离散时间和连续时间仿真模型的描述
6.4.1 污染模型
POL(t+h)=POL(t)+[POLG(t)-POLA(t)]h POLG(t)=P(t)*POLCM(t) POLCM(t)=POLCMT(CIR(t)) CIR(t)=
CI (t )
P (t )
POLA(t)=POL(t)⨯POLR(t) POLA(t)=POLRT(POL(t))
⎡⎤⎛CI (t ) ⎫
Pol (t +h ) =Pol (t ) +⎢P (t ) ⨯POLCMT -POL (t ) ⨯POLRT (POL (t )) ⎪⎥⨯h P (t ) ⎪
⎝⎭⎣⎦
对于状态变量POL
δPOL (pol , ci , p ) =pol +[p ⨯POLCMT () -pol ⨯POLRT (pol )]h
ci p
1．瞬时函数：根据时间t 的各输入值确定模型在时间t 的输出值的函数。
2．记忆或非瞬时函数：用微分方程说明，它可用连续时间和离散时间来进行仿真。
1． 已标识了瞬时函数及记忆函数的模型网络描述，每个模型网络图可容易地转换成某些仿真语言的语句。
2． 由模型网络描述，也可容易地明确状态转移函数和输出函数，离散时间模型的仿真就可按其仿真过程进行。
3． 仿真过程中，按某个固定顺序搜索状态变量，用有关的局部转移函数计算每个状态变量来获得下次状态值。由于某些瞬时函数可能会直接影响两个以上的转移函数或输出函数，会不可避免地出现不止一次地计算这个瞬时函数。如果能找出各个瞬时函数的计算次序，就可以这一问题。为此，给出一种非常类似于仿真语言的模型描述语言，它可描述编排模型描述语句的过程。
6.4.2 模型描述语言
模型描述语言是一种类似仿真语言的用来描述编排模型描述语句的过程。模型描述语言中的三种语句： 1) 瞬时函数
Y 1, Y2, ⋯, Ym = Instant. Func (X1,X 2, ⋯ Xn )
在t 时刻，能给定变量X 1, ⋯, X n 的值，输出Y 在时刻t 的值直接可由上式获得。 2) 输入时间函数
X 1, X2, ⋯, Xm = Time. Func
它是用来产生模型的输入轨迹，同瞬时函数不同，它仅与时间t 有关
例X 1, X2=sin, cos
X 1, X2的轨迹为 (sin(t), cos(t))
3) 记忆函数
Y 1, Y2, ⋯, Ym = Mem. Func (Q1, Q2, ⋯, Qn ; X1, ⋯, Xp )
它具有状态变量Q 1, ⋯, Qn , 其初值必须设定，本质上它包括输入轨迹、状态转移函数和输出函数三个部分。
对于记忆函数，用离散时间和连续时间仿真，其处理方式不同，它们分别可用一个延迟元件和积分器来描述
记忆函数的两种基本环节：
对于微分方程描述的系统，记忆函数的基本环节是积分器。积分器可用语句Y =INTGRL (IY , YRT ) 来表示，它的输出Y 是输入YRT 加上初始状态IY 的时间积分。它是构成微分方程描述的系统的基本环节。可把积分器转化成离散时间仿真形式，但其对应的转移函数只可能是近似的。
对于离散时间模型，记忆函数的基本环节是迟延元件。迟延元件可用语句 Y =DELAY (IY , YP ) 来表示。对于离散时间的一个DELAY(迟延) ，其输出轨迹Y （·）比输入轨迹YP(·) 滞后一步，它的初始值由IY 给定。DELAY 语句能够描述每一种离散时间模型。
模型描述语句与模型网络图存在一一对应关系，可从模型描述语句直接推导出模型网络图。例如，描述语句 U =PROD（Y1，Y2），表示瞬时函数PROD 有两条Y1和Y2输入线和一条U 输出
线；描述语句 Y1=DELAY（IY1,Y1P ），说明Y1也是DELAY 函数的输出线；这样，两个DELAY 方框的输出作为 PROD 的输入。 从上面的例子，发现描述序列中语句的语序与描述无关，也就是从语句序列S 1,S 2, „，S n 的任何排列中均可得到相同的网络。
考虑以下模型描述： Y2=DELAY（IY2,Y2P ）
U =PROD（Y1，Y2） Y1P=SUM（X1，U ） Y2P=SUM (X2，U)
X1，X2=SIN，COS Y1=DELAY（IY1,Y1P ）
U
图6.5 模型描述语句对应的网络图
6.4.3 模型描述语句序列分析
在描述序列中，语句的语序同描述无关，一个描述序列的任何排列均可得到相同的网络。但在模型仿真中，语句被翻译计算机的指令时，则必须正确次序执行这些动作。
描述语句序列必须经过分析后，才能确定它是否有效及正确执行次序。
1．无效语句的检查
检查是否有任一变量在两个不同语句的左边出现，若存在则是无效语句
2．无记忆和输入时间函数时的函数循环问题判别
（1）领先关系
对于Y 1, Y2, ⋯, Ym = Instant. Func (X1,X 2, ⋯ Xn )
若U 是X i 中的一个变量，V 是Y j 中的一个变量，则U 领先于V 。若U 领先V ，在计算t 时刻的V 值之前，必须要知道时间t 的U 值。 （2）领先关系的闭包传递
若有一变量序列W o , ⋯, Wn , Wo = U, Wn = V, 各个W i 领先W i+1, 则U 领先*V，
U 领先*V表示在网络图中存在一个沿箭头从U 到V 的路径 例：V = Prod (Y1, Y2)
Y 1=Sum (X1, V) Y 2=Sum (X2
,V)
结论:有存在V 领先*V，模型描述语句构成循环，存在无效描述
3．变量的次序排列
1）变量的级别
（a ）若对某个变量V ，变量U 领先于V ，且无变量W ，它使W 领先U ，那么变量U 的级别Lev （U ）=0
（b ）任一（其他）变量的级别是网络中从级别为0的变量到该变量最长路径的距离。
若U 领先*V，级别Lev （U ）=0，n 为最大数，U 领先W 1领先W 2⋯W n-1领先V ，级别Lev （V ）=n
举例
2）语句的级别（函数的级别）
对于语句S i : Y1, ⋯, Ym = Instant.func (X1, X2, ⋯, Xn )
语句S i 的级别是X j 的最大级别，即级别(Si )=max{级别(Xj ), j=1, n} 求Prod 和Sum 的级别。
根据这些排序层次，就可自动构造层次关系清楚的模型网络图。在构造模型网络图时，把相同排序层次的所有变量和语句置在一条线上。
对于例6.7的模型描述语句：
LEV[U=PROD（Y1，Y2）]=max｛LEV （Y1），LEV （Y2）｝=0
LEV [Y1P=SUM（X1，U ）]= LEV （U ）=1 LEV [Y2P=SUM（X2，U ）]= LEV（U ）=1
排序层次
图6.7 模型网络图
对于瞬时函数，排列次序按它的级别排列S 1, S2, ⋯, Sn , 级别(Si ) ≤级别(Si+1) 。一般来说，时间函数语句排在最前，记忆型函数则排列最后。最终序列形式为：
────
┆ TIME ·
FUNOT ────
┆ 排序层次为0 ────
┆ INSTANT ·FUNCT ────
┆ 排序层次为M ────
┆ MEM ·FUNCT ────
其中M 是语句的排序层次的最大值。 X 1, X2=Sin, Cos U = Prod (Y1, Y2) Y 1P = Sum (X1, U) Y 2P= Sum (X2, U) Y 1 = Delay (IY1, Y1P)
Y 2=Delay (IY2, Y2P),
6.4.4 记忆函数仿真
1、离散时间模型
当模型描述仅包含 DELAY 记忆函数语句时，可选择 DELAY 的输出变量来构造模型的状态变量。因而每个 Y=DELAY（IY ，X ）语句可解释如下：在预置阶段，设置 Y 等于IY ；在状态转移阶段，设置Y 在时间t i +h的值为X 在时间t i 的值。
对于例6.6，其模型描述语句对应的仿真的基本过程为： 1） 预置时标T 为要求的初始时间t ； 2） 置 Y1为IYl 、 Y2为IY2; 3） 置 X1，X2为SIN(T）,COS(T); 4） 置U 为PROD(Y1,Y2）；
5） 置Y1P 为 SUM （X1,U ）、Y2P 为SUM （X2,U ）； 6） 置Y1为 Y1P 、Y2为 Y2P ； 7） 推进仿真时标，置 T 为T+h； 8） 若T 小于终止时间，转向3）； 9） 停机。
常态形序列，则可以每个延迟元件为中心构造描述语句，可以起到对这个模型描述序列的分解作用。 设S 1，...,S m 是最终序列，并令
Y1=DELAY（IYl ；Y1P ）
┆
Yn=DELAY（IYn ，YnP ） 是记忆函数语句。对于每个迟延元件的输入变量 YiP ，设先前变量集PRIOR ·V ARIABLES i 是计算YiP 有关的所有变量集，也就是PRIOR ·V ARIABLES i 是这样的一个变量集，它的所有变量领先*于Y i P 。现设S i1,S i2，„，S im 是计算YiP 有关的最终序列的子序列，即每个S ij 包含YiP 和（或）包含至少一个PRIOR ·V ARIABLES i 中的变量。
这样可构造常态形序列，其基本形式如下：
──
┆ ── Y1的转移函数 Y1 = DELAY (IY1; Y1P)
──
┆ ── Y2的转移函数 Y2 = DELAY (IY2; Y2P) ┆
──
┆ ── Yn 的转移函数 Yn = DELAY (IYn; YnP)
── 输出函数 ┆ 剩余语句 ──2、微分方程说明的模型
在连续时间模型中，一般用微分方程来描述系统模型，其记忆函数语句就是积分器。积分器Y=INTGRL（IY ，YRT ），是建立起它的输入YRT 和输出Y 所假定的值之间的关系
dY
YRT (t ) dt
这样，就可用积分法根据输出变量去逐次计算描述变量值，对整个导
数函数可通过一定手段构造对应的状态转移函数，用迭代方式进行系统仿真。
以下是欧拉积分法的仿真过程： 1） 预置T=to, Y=IY；
2） 计算中间变量M=Y+h·YRT ； 3） 推进时标T=T+h, Y=M；
4） T 是否小于终止时间，否，转2）； 5） 结束。
作业：对于滞后变量的伪随机数发生器，其算法的步骤：
1）设S -5，S -4，S -3，S -2，S -1是区间[0，m]中的各个整数（逻辑种子） 2）置i=0
3）选取S i-3,S i-5, 计算：S i 为aS i-3+bSi-5关于m 的模 r i= Si /m.
4）i=i+1,转入3）。
1) 给出本伪随机数产生模型的离散时间系统规范； 2) 给出本模型网络图；
3) 给出本模型描述语句和去掉记忆环节后的各变量和语句的排序层次。