第二章 经典的连续系统仿真建模方法学
本章讨论经典的连续系统数字仿真的原理与方法，内容包括连续系统数字仿真的基本概念、经典的数值积分法、经典的线性多步法等。在数字计算机上进行连续系统仿真，首先要将连续模型离散化，因此，2.1节首先讨论离散化原理及要求，这是连续系统仿真的基础。然后，2.2节对经典的数值积分法----龙格-库塔法及其它典型的数值积分法仿真建模原理进行详细分析，并通过实例说明其应用要点；而2.3节对经典的线性多步法进行了介绍．
2.1 离散化原理及要求
在数字计算机上对连续系统进行仿真时，首先遇到的问题是如何解决数字计算机在数值及时间上的离散性与被仿真系统数值及时间上的连续性这一基本问题。
从根本意义上讲，数字计算机所进行的数值计算仅仅是“数字”计算，它表示数值的精度受限于字长，这将引入舍入误差；另一方面，这种计算是按指令一步一步进行的，因而，还必须将时间离散化，这样就只能得到离散时间点上系统性能。用数字仿真的方法对微分方程的数值积分是通过某种数值计算方法来实现的。任何一种计算方法都只能是原积分的一种近似。因此，连续系统仿真，从本质上是对原连续系统从时间、数值两个方面对原系统进行离散化，并选择合适的数值计算方法来近似积分运算，由此得到的离散模型来近似原连续模型。如何保证离散模型的计算结果从原理上确能代表原系统的行为，这是连续系统数字仿真首先必须解决的问题。
设系统模型为：),,(t u y f y
= ，其中u (t )为输入变量，y (t )为系统变量；令仿真时间间隔为h ，离散化后的输入变量为)(?k t u
，系统变量为)(?k t y ，其中k t 表示t=kh 。如果)()(?k k t u t u ≈，)()(?k k t y t y ≈，即0)()(?)(≈-=k k k u t u t u
t e ，0)()(?)(≈-=k k k y t y t y
t e （对所有k=0,1,2,…），则可认为两模型等价，这称为相似原理（参见图2.1）。
实际上，要完全保证0)(,0)(==k y k u t e t e 是很困难的。进一步分析离散化引入的误差，随着计算机技术的发展，由计算机字长引入的舍入误差可以忽略，关键是数值积分算法，也称为仿真建模方法。相似原理用于仿真时，对仿真建模方法有三个基本要求：
（1）稳定性：若原连续系统是稳定的，则离散化后得到的仿真模型也应是稳定的。关于稳定
性的详细讨论将在2.4节中进行。
（2）准确性：有不同的准确性评价准则，最基本的准则是：
图2.1 相 似 原 理
绝对误差准则：δ≤-=)()(?)(k k k y t y t y
t e 相对误差准则：δ≤-=)(?)()(?)(k k k k y t y
t y t y
t e
其中δ 规定精度的误差量。
（3）快速性：如前所述，数字仿真是一步一步推进的，即由某一初始值)(0t y 出发，逐步计
算，得到)(,),(),(21k t y t y t y ，每一步计算所需时间决定了仿真速度。若第k 步计算对应的系统时间间隔为,1k k k t t h -=+计算机由)(k t y 计算)(1+k t y 需要的时间为T k ，则，若T k =h k 称为实时仿真，T k
连续系统数字仿真中离散化最基本的算法是数值积分算法。对于形如),,(t u y f y
= 的系统，已知系统变量y 的初始条件y t y ()00=，现在要求y 随时间变化的过程y t ()。计算过程可以这样考虑（参见图2.2）：首先求出初始点y t y ()00=的f t y ()00，，微分方程可以写作：
y t y f t y dt t t
()(,)=+?00
（2.1）
图2.2所示曲线下的面积就是y t ()，由于难以得到f(y,u,t)积分的数值表达式，人们对数值积分方法进行了长期探索，其中欧拉法是最经典的近似方法。 欧拉法用矩形面积近似表示积分结果，也就是当
t=t 1时，y t ()1的近似值为y 1 ：
y y t y t f t y 11000=?+?()()?， (2.2) 重复上述作法，当t t =2时
y y t y t t f t y 2212111=?+-?()()()， 所以，对任意时刻t k+1，有：
y y t y t t f t y κκκκκκκ+++=?+-?111()()()， （2.3） 令t t h κκκ+-=1称为第κ步的计算步距。若积分过程中步距不变h h κ=，可以证明，欧拉法的截断误差正比于h 2
。
为进一步提高计算精度，人们提出了“梯形法”。 梯形法近似积分形式如式(2.4)所示，令：t t h h κκκ+-=1＝已知：t t =κ时y t ()κ的近似值y κ,那么：
y y t y h f t y f t y κκκκκκκ++++=?+
+11111
2
()[(,)(,)] (2.4) 可见，梯形法是隐函数形式。采用这种积分方法最简单的预报?校正方法是用欧拉法估计初值，用梯形法校正，即：
y y h f t y f t y i i κκκκκκ++++?+
+111112
()
()[(,)(,)] (2.5) y y h f t y i κκκκ+?+?1()(,) (2.6)
式(2.6)称作预报公式，采用欧拉法，式(2.5)为校正公式，采用梯形法。用欧拉法估计一次
y i κ+1()的值，代入校正公式得到y κ+1的校正值y i κ++11()。设ε是规定的足够小正整数，称作允许误差，若i =0, i +1=1称作第一次校正；i i =+=112，称作第二次校正；通过反复迭代，直到满足y y i i
κκε+++-≤111，这时y i κ++11()是满足误差要求的校正值。
上述方法是针对(2.3)式所示的微分方程在已知初值情况下进行求解，因此也称为微分方程
初值问题数值计算法，为统一起见，本书中称为数值积分法。
连续系统数字仿真的离散化方法有两类，它们是数值积分方法和离散相似方法，本章讨论数值积分法。
数值积分方法采用递推方式进行运算，而采用不同的积分方法会引进不同的计算误差，为了提高计算精度，往往会增加运算量。就同一种积分算法而言，为提高计算精度，减小积分步距h ,计算量增大，影响系统运算速度。因此，计算精度和速度是连续系统仿真中常迂到的一对矛盾，也是数字仿真中要求解决的问题之一。也就是说，选择合适的算法、合适的软、硬件环境，在保证计算精度的前提下，考虑怎样提高仿真的速度。
经典的数值积分法可分为两类：单步法与多步法，下面我们将分别来介绍这两类方法中的最常用的算法。
2.2 龙格库塔法
2.2.1龙格-库塔法基本原理
在上一节中我们已经讲过，在连续系统的仿真中，主要的数值计算工作是对
dy
dt
f t y =(,)的一阶微分方程进行求解。因为 y t y t f t y dt t t ()()(,)κκκ
κ+=+
+?
11
若令： y y t κκ?(),Q κκ
κ?
+?
f t y dt t t k (,)1
则有 y t y y ()κκκκ++?=+11Q (2.7)
因此主要问题就是如何对Q κ进行数值求解，即如何对f t y (,)进行积分。通常称作“右端函数”计算问题。已知y t y ()00=，假设我们从t 0跨出一步，t t h 10=+，t 1时刻为
y y t h 10=+()，可以在t 0附近展开台劳级数，只保留h 2项，则有：
20
0010)(21),(h t
f
dt dy y f h y t f y y t ????+++= (2.8) 在式(2.8)中括号后的下标0表示括号中的函数将用t t y y ==00，均同。我们假设这个解可以写成如下形式：
y y a k a k h 101122=++() (2.9) 其中， k f t y 100=(,)
k f t b h y b k h 201021=++()，
对k 2式右端的函数在t t y y ==00，处展成台劳级数，保留h 项，可得：
h k b t
f
b y t f k t 0
)(),(121002++??? 将k k 12，代入(2.9)式，则有：
])(),([),(0
12100200101h y
f
k b t f b y t f h a y t hf a y y t ????++++= 将上式与(2.8)式进行比较，可得：
可见有四个未知数a a b b 1212，，，，但只有三个方程，因此有无穷多个解，若限定
a a 12=，则可得其中一个解：
a a
b b 12121
2
1====，
将它们代入(2.9)式可得一组计算公式：
y y h
101
22=+
+()K K (2.10) 其中，K K K 1002001==++f t y f t h y h (,)()，， 若写成一般递推形式，即为：
y t y y h
k k k ()()++?=+
+111
22K K (2.11) 其中， K K K 121==++f t y f t h y h k k k k (,)(,)，
由于(2.8)式只取了h h ，2两项，而将h 3
以上的高阶项略去了，所以这种递推公式的截断误差正比于h 的三次方，又由于计算时只取了h 及h 2
项，故这种方法被称为二阶龙格?库塔法（简称RK-2）。
根据上述原理，若在展成台劳级数时保留h h h ，，23及h 4
项，那么可得一套截断误差正比于h 5
的四阶龙格--库塔法（简称RK-4）公式：
y t y y k k k
()()++?=++++111
2341
622K K K K (2.12) 其中： K 1=f t y (,)κκ
K K 2122=++f t h y h
()κκ， K K 3222
=+
+f t h y h
()κκ， K K 43=++f t h y h ()κκ，
由于这组计算公式有较高的精度，所以在数字仿真中应用较为普遍。
为了帮助读者更好地掌握这种方法的使用，下面我们对龙格--库塔法的特点作一些介绍。 1. 由于在解a a b b 1212，，，时，可以得到许多种龙格--库塔公式，经常使用的除(2.11)及(2.12)式给出的两组公式外，还有：
y y h κκ+=+12K (2.13) a a a b a b 12212211212
+===,/,/
其中 K 1=f t y (,)κκ K K 2122
=++f t h y h
()κκ， y y h
κκ+=+
+++11
234833()K K K K (2.14) 其中 K 1=f t y (,)κκ
K K 2133
=++f t h y h
()κκ， K K K 312233
=+
-+f t h y h
h ()κκ， K K K K 4123=++-+f t h y h h h ()κκ，
(2. 13)式为二阶龙格--库塔公式，(2.14)式是Shampine 提出的，称为RKS3-4公式。 对上述公式可以这样解释：因为y t ()满足
dy
dt f t y =(,)这个微分方程，已知(,)
t y κκ处y t ()的导数为f t y (,)κκ=K 1。根据外推原理，假设在t κ--t h κ+这个区间中y t ()的导数不变，则可得y t ()在t t h =κ＋的估计值y h i κ+K 0。令K i 为这个估计值处的导数，因此(2.11)式就是按
()
K K 122
+这个导数线性外推而得的y κ+1。而(2.12)式则是4点导数的加权平均和，由
于第2及第3点是以h 2为步长，精度较高，所以加权平均时各取2份，第1、4点各取1份。不论几阶龙格--库塔法，它们的计算公式总是由两部分组成的，即是上一步的结果y κ及步长
h 乘以t κ--t κ+1中各点导数的加权平均和。因此各种龙格---库塔法可以写成如下一般形式：
y y h
C i i
i s
κκ+==+∑11
K
(2.15)
其中 K K
i i ij j
j i f t a h y h b =++=-∑()κκ，1
1
i s =12，，，
式中各系数满足以下关系
a a
b i s
C i ij j i i i s
11
11
231====??
?
??????=-=∑∑，，，
(2.16)
一般来讲，为了减少计算量总是希望减少每步计算右端函数f 的次数（即减少计算K 的次数）。可以证明，1阶龙格--库塔公式至少要计算一次右端函数，即(2.11)式中的s(称为级数)至少等于1，而2阶公式s min =2；3阶公式s min =3；4阶公式s min =4；依此类推。有时为了某种特殊需要，可以选择s >s min 的计算公式。
表2.1为常见的各阶龙格---库塔公式，其中第二项是4阶4级RK-4公式1，而第三项是
4阶5级RK-4计算公式(简称RKM45)。其中，后三项有误差估计公式，所以可以用误差来控制积分步长，即每积分一步，利用误差估计公式估计出这一步的截断误差。若它大于允许的误差，则缩小步长，反之可以放大步长。
表2.1 常见的龙格--库塔公式
2. 所有龙格---库塔公式都有以下特点：
(1) 在计算y κ+1时只用到y κ，而不直接用y y κκ--12，等项。换句话说，在后一步的计算中，仅仅利用前一步的计算结果，所以称为单步法。显然它不仅能使存储量减小，而且此法可以自启动，即已知初值后，不必用别的方法来帮助，就能由初值逐步计算得到后续各时间点上的仿真值。
(2) 步长h 在整个计算中并不要求固定，可以根据精度要求改变，但是在一步中，为计算若干个系数K i (俗称：龙格---库塔系数)，则必须用同一个步长h 。
3. 龙格--库塔法的精度取决于步长h 的大小及方法的阶次。许多计算实例表明：为达到相同的精度，四阶方法的h 可以比二阶方法的h 大10倍，而四阶方法的每步计算量仅比二阶方法大一倍，所以总的计算量仍比二阶方法小。正是由于上述原因，一般系统进行数字仿真常用四阶龙格--库塔公式。值得指出的是：高于四阶的方法由于每步计算量将增加较多，而精度提高不快，因此使用得也比较少。
4. 若在展开台劳级数时，只取h 这一项，而将h 2
以及h 2
以上的项都略去，则可得： y y hf t y κκκκ+=+1(,)
这就是欧拉公式，因此欧拉公式也可看作是一阶龙格--库塔公式。它的截断误差正比于h 2
，是精度最低的一种数值积分公式。
读者也许会问：如果我们将步长取得十分小，那么欧拉公式也可以获得很高的计算精度吗？从理论上讲是这样，但是实际上由于计算机字长有限，在计算中存在着舍入误差，它与计算次数成正比。计算步长很小，则一个系统的过渡过程的计算就会增加许多，舍入误差就会十分明显地表现出来，这样一来，很难保证有较高的精度。 2.2.2 龙格--库塔法的误差估计及步长控制
一个高精度的仿真方法必须将步长的控制作为必要的手段。实现步长控制涉及局部误差估计和步长控制策略两方面的问题。本节将对数值积分中单步法的误差估计和步长控制作一介绍。
龙格--库塔数值积分方法的“误差估计”和“步长控制”的基本想法如图2.2所示,也就是，每积分一步都设法估计出本步的计算误差εk ，然后判断是否满足允许误差E ，据此，选择相应的步长控制策略，调整步长，再作下一步积分运算。
(1) 龙格--库塔法的误差估计
龙格--库塔法的误差估计，通常采用的方法是设法找到另一个低阶(一般是低一阶)的 龙格--库塔公式，则两个公式计算结果之差可以被看作是误差。例如，龙格--库塔--默森法(Runge-Kutta-Merson)法[6]计算公式为： )4(6
5411K K K +++
=+h
y y k k (2.17) 其中 ))
34(2
())
3(8
2())(63()
33()
,(3415314213121K K K K K K K K K K K K K -+++=+++=+++=++==h
y h t f h
y h t f h
y h t f h
y h t f y t f k k k k k k k k k k ，，，，
它是一个4阶5级公式。另外还可推出一个3阶4级公式
)1293(6
?4311K K K +-+=+h
y y
k k (2.18) 用11?++-k k y y
作为误差估计，即 )892(6
5431K K K K -+-=
h
E k (2.19) 式(2.17)、(2.18)、及(2.19)简称为RKM3-4法。又比如，E.Fehlberg 推导出的5阶龙格--库塔--费尔别格法[4]，它的计算公式为一个5阶6级方法，另外用了一个四阶5级方法求
1?+k y
，也是用11?++-k k y y 来估计误差，这一套计算公式被公认是对非病态系统进行仿真最为有效的方法之一。由于它是5阶精度，4阶的误差估计，因此被称为RKF4-5阶公式对，简称RK4-5法。具体系数如表2.2。
表2.2 Fehlberg 公式中的各系数
图2.2龙格?库塔法步长控制示意
???????
++=+=∑∑-==+1
16
1
1)
(i j i
ij k i k i i i i k k b h y h a t f c h y y K K K ， (i =126，， ) (2.20) 表中c i *表示4阶公式中的系数，因此误差估计为
∑
++-=-=b
i i i
k k k c c
h y y
E *
11)(?K (2.21)
当精度要求不太高时，RKF1-2法被广采用。它的计算公式为
)510(512
3211K K K +++
=+h
y y k k (2.22) 其中 ??
?
?
?
?
???
+++=++==))255(256()22(),(213121K K K K K K h y h t f h y h t f y t f k k k k k ，，
用另一个一阶公式
)255(256
?211K K ++=+h
y y
k k (2.23) 来估计误差
)(512
?3111K K -=-=++h
y y
E k k k (2.24) RKM3-4法与RKF4-5法是两种应用得极为广泛的数值积分算法。它们的共同特点是计算量较大，比如RKM3-4法，每次要计算5次f ,而只获得四阶精度及三阶误差估计，RKF4-5法每步要计算6次f ,能获得5阶精度及4阶误差估计。对于一般仿真问题，中等精度即可满足要求，因此用RKM3-4公式对就够了。1978年，夏普勒(Shampine)提出了RKS3-4公式对[7]，它每次只计算4次f 却能获得4阶精度与3阶误差估计，如(2.25)式所示。
)33(8
43211K K K K ++++
=+h
y y k k (2.25)
其中 ?????
??
?
?+-++=+-++=++==κκ))
(())
3(3
32()
3
3()
,(3214
213121K K K K K K K K K K h y h t f h
y h t f h y h t f y t f k k k k k k ，，，
另外引入了一个3阶公式
)49153(32
?543211K K K K K +++++=+h
y y
k k (2.26) 其中 ))33(8(43215K K K K K ++++
+=h
y h t f k k ，
(2.27) 它正好是下一次计算1+k y 时的K 1，因此只是在第一步要多算一次f ，以后仍每步计算四次
f 。RKS3-4的误差估计为
)4333(32
?5432111K K K K K +--+-=-=++h
y y
E k k k (2.28) (2). 步长控制 有以下几种步长控制的方法：
1. 加倍－减半法
设定一个最小误差限和最大误差限，当估计的局部误差大于最大误差时将步长减半，并重新计算这一步；当误差在最大误差与最小误差之间时，步长不变；当误差小于最小误差时将步长加倍。
每步的局部误差通常取以下形式：
)1(+=k k k y E e (2.29)
其中k E 为利用本节的公式计算出来的误差估计。由(2.29)式可知，当k y 较大时，k e 是相对误差，而当k y 的绝对值很小时，k e 就成了绝对误差。这样做的目的是避免当y 的值很小时，
k e 变得过大。
上述步长控制的策略可以表示为下式：
如果???
??≥εε<<εε≥k k k e
e e min
max
min max
则选择???
?
???
===+++k
k k k k k h
h h h h h 221111 (2.30)
其中εmax 和εmin 为最大、最小误差限。这种步长控制的方法简便易行，每步附加计算量小，但是这种方法不能实现每步都是最优步长。
2. 最优步长法
为了使每个积分步在保证精度的前提下能取最大步长（或称最优步长），可以设法根据本步误差的估计，近似确定下一步可能的最大步长。。这种方法可以做到在规定的精度下取得最大步长，因此减少了计算量。具体策略如下：
给定相对误差限ε0，设本步步长为k h ,本步相对误差估计值为k e (如2.29式所示)。假定
所采用的积分算法为K 阶，则(2.29)式中的k E 可表示为
K h E m ?ζ?=)( (2.31)
其中()
?ζ是f t y (,)在积分区间),(h t t k k + 中某处一些偏导数的组合，通常可取ζκ=t ，因此有
)1/()(+??=k k k y h t e K
(2.32) 若0ε≤k e 则本步积分成功。现在来确定下一步的最大步长1+k h ，假定1+k h 足够小，即可认为)()(1k k k t h t ?≈+?+，故下一步的误差可能为
)1/()()1/()(111
1+??≈+??=++++k k k k k y h t y h
t e K K
为使01ε≤＋k e ，则有 K
/101))(/)1((k k k t y h ?+ε≈+
将(2.32)代入上式，得
K
K
/101)
/(k k e h h ?ε≈+h e k ?ε=K /10)/( (2.33)
若0ε>k e ，则本步失败。此时仍可采用(2.33)式，但它表示重新积分的本步步长。由于我们假定了1+k h 足够小，因此)(k t ?基本不变，故必须限制步长的无限放大和缩小，一般可限制h 的最大放、缩系数为10，即要求
k h h h ?<<+κκ101.01 (2.34) 当f 函数中含有间断特性时，采用上述两这控制方法在间断点附近会出现步长频繁放大、缩小的振荡现象。由于最优步长控制方法是以本步误差外推下一步步长，因此振荡现象更为严重。解决这一问题的方法将在第7章讨论。 2.2.3 实时龙格－库塔法
仿真模型的运行速度往往与实际系统运行的速度不同。然而，当有实际的装置或被训练的人介入仿真过程时，就要求仿真模型的运行速度往往与实际系统运行的速度保持一致，这称为实时仿真。
一般的数值积分法难以满足实时仿真的要求，这不仅仅是因为由这些方法所得到的模型的执行速度较慢，而且这些方法的机理不符合实时仿真的特点。下面让我们以二阶龙格-库塔为例来分析。假设对如下一般形式的系统进行仿真：
))((t u y f dt
dy
，= RK-2公式如下：
)
,,()
,,()(2
11121211h y u t f y u t f h
y y k k k k k k k k K +=K =K K +K +
=+++
即在一个计算步内分两子步：首先在t k 时刻利用当前的u k ，y k 计算K 1。假定在h /2的时间内计算机正好计算一次右端函数f ，然后，在t k +h/2时刻计算K 2，尽管此时y k +1/2已经得到，
但u k +1则无法得到。实时仿真除了对模型执行速度满足要求外，还要求实时接收外部输入，并实时产生输出。为此，要么对u k +1也进行预报（这将加大仿真误差），要么将仿真执行延迟h/2。后者的RK-2公式的计算流程可表示成图2.4。这就是说，输出要迟后半个计算步距。
与此相类似，(2.12)式所示的RK-4公式也不适合于实时仿真，对此，读者可自行分析。
为了克服这个缺陷，人们提出了如下形式的实时2阶龙格－库塔法：
)2
,,()
,,(12/12/1212
1K K K K h y u t f y u t f h y y k k k k k k k k +
==+=+++
它的计算流程如图2.5所示。由图不难看出，首先在t k 时刻利用当前的u k ，y k 计算K 1。然后，在t k +h/2时刻计算K 2，此时y k +1/2已经得到，由于计算一次右端函数f 需要h /2的时间，u k +1/2也可得到，
文献[1] ]105515[24
543211K K K K K ++-+-+
=+h
y y k k 其中
)
2
10354()
5253()
52()
5
5()
,(41521413121K K K K K K K K K K K h
h y h t f h h
y h t f h
y t f h
y h t f y t f k k k k k k k k k k +++=+-+=++=++==，，， 它的计算流程如图2.6所示： 图2.4 RK-2的计算流程
该公式的特点是：将一个仿真步分为5个子步，每个子步均能对外部输入实时采样，并
保证了实时输出。
2.2.4 面向方程的龙格－库塔法仿真举例
为了加深对数值积分法的理解，以便能使用仿真程序来解决一两个实际问题，本小节将给出面向方程的龙格－库塔法仿真实例。“面向方程”是指仿真对象是以一阶微分方程组形式描述的系统，如(2.35)式所示：
)
0()0()
,,,(Y t Y t U P Y G dt dY
=== (2.35)
其中，Y y y y n =(,,)12 是n 维状态向量；),,(21k p p p P =是k 维参数向量；
U u u u m =(,),12 是m 维控制向量；G g g g n =(,,)12 是n 维函数向量；Y y y y n ()[(),(),()]000012= 是给定的初始状态。
本例所用仿真程序CSS1.C 采用定步长RK-4公式作数值积分。CSS1.C 是由C 语言编制的面向方程的仿真程序(原程序见附录)，其程序结构如图2.7所示。
图2.7 C S S 1.C 程序结构
程序中主要变量说明：
T 1 仿真总时间 T 2 积分步长
T 仿真时间 Y(1) ?Y(20) 系统状态变量 Y(21) ?Y(30) 用户自定义变量 G(1) ?G(20) 系统状态变量的导数 A(1) ?A(5) 用户输出变量 P(1) ?P(20) 系统的参数变量 N 1 方程阶次 N 2 参数个数
N 3 输出点数（在0≤T ≤T 1区间内，设输出时间间隔为T 3，则N 3＝T 1/ T 3+1） 下面对仿真对象作一说明。设卫星在空中运行的运动方程为：
d r dt k r
r d dt d dt r dr dt d dt 222
2
2
2
2=-+=-???????()θθθ (2.36)
其中k 是重力系数(k =401408km 3/s)。卫星轨道采用极坐标表示，通过仿真，研究发射速度v 对卫星轨道的影响。
这是一个二阶微分方程组，为此，首先要将其转换成一阶微分方程组。若设
y r y (),()12==θ，则有：
dy dt y dy dt
y dy dt k
y y y y y dy dt y y y ()
()()()()()()()()()()()()/()132431114442341===-*+**=-**????
?????
??
(2.37) 这是四个一阶微分方程，有4个状态变量y y y (),(),()123及y ()4。今希望用直角坐标输出，故还引入两个定义变量x y y y ==()()2122，，这样又得到两个代数方程：
y y y y y y ()()cos[()]
()()sin[()]21122212=*=*???
(2.38)
根据卫星的发射速度，可以建立起方程组(2.37)的初始条件：y ()164000=km(卫星到地心的距离，即为地球之半径)，y y y v y ()()()/()230410000===，，仿真研究的目的是要计算出当发射速度为8km/s ，10km/s 及12km/s 时的卫星轨道。
已知y(4)=y/y(1)，由y 及y ()10可求出y ()40，它们分别是：0.00125，0.0015625，0.001875(1/s)，按CSS1.C 的要求，为了仿真该系统，用户要输入系统的运动方程，即要将(2.37)及(2.38)式写到程序中去。CSS01.C 规定这部分写在系统模型输入(modsub)程序块中；比如例中要求输出y y y y ()()()()212212，，，则A(1)，A(2)，A(3)，A(4)四个数组都要用到，CSS1.C 规定这部分写在输出转换(repsub)程序块中；状态变量的初始数据及要改变的参量数值，CSS1.C 规定它放在主程序中。此系统没有参数变量。
下面就是为仿真卫星轨道系统由用户写的一段仿真源程序。在系统模型输入程序块(modsub)中写入以下程序：
Y(21)=Y(1) *cos(Y(2))
Y(22)=Y(1) *sin(Y(2))
G(1)=Y(3)
G(2)=Y(4)
G(3)=-401408/(Y(1) *Y(1))+Y(1) *Y(4) *Y(4))
G(4)=-2*Y(3) *Y(4)/Y(1)
在输出转换(repsub)程序块中写入以下程序：
A(1)=Y(21)
A(2)=Y(22)
A(3)=Y(1)
A(4)=Y(2)
如果需要，可以在输出打印(outsut)程序块中加入Y(4)0=Y(4)的打印语句。
程序写好后则可按以下步骤进行操作：
(1)将CSS1.C程序及用户写好的仿真源程序输入到计算机中。
(2)运行CSS1.C按计算机提问输入必要的数据。
(3)计算机执行CSS1.C并在终端上或打印机上输出运行结果。
下面给出用CSS1对发射速度为8Km/s时卫星轨道仿真的运行过程和结果。C:/>CSS1.C
Input total simulation time T1, step t2=10000, 200 ?
Input the order of system N1=4 ?
Input the number of unknown parameters N2=0 ?}
Input the number ofoutput points N3=51symbol 191 \f "Symbol" \s 10.5?
Input the number of runs J1=3 ?
Input initial values of staste variables=6400，0，0，00125 ?
Input the values of parameters=0 ?
OUTPUT DEVICE ---2(PRINTER),5---CRT?
然后机器开始运算，下面是运算一次后的结果：
Y(I) Y(I)max Y(I)min
1 6666.64600 .00000
2 11.9750
3 .00000
3 .15935 -.16005
4 .0012
5 .00000
21 6400.00000 -6665.12300
22 6484.59500 -6475.95700
SIMULA TION STEP T2=200.00000
T A(1) A(2) A(3) A(4) 200.00000 6205.06800 -2625.19900 6403.97800 .24990 400.00000 5632.80400 460.99630 6415.65600 .49918 600.00000 4719.81600 5857.74100 6434.28800 .74728 800.00000 3523.71300 -5904.65500 6458.69700 .99372 1000.00000 2118.52300 5751.83300 6487.37500 1.23814 1200.00000 589.16300 -6475.95700 6518.58600 1.48029 1400.00000 -974.40360 -3196.89000 6550.49800 1.72010 1600.00000 -2482.78100 -5245.30500 6581.28300 1.95762 1800.00000 -3852.12900 -3377.11100 6609.23000 2.19302
2000.00000 -5008.32900 -3907.27800 6632.82900 2.42657
2200.00000 -5890.24000 -1619.59200 6650.82900 2.65867
2400.00000 -6452.09600 4470.64700 6662.29700 2.88973
2600.00000 -6665.12300 6484.59500 6666.64600 3.12024
2800.00000 -6518.49300 -2042.90700 6663.65500 3.35070
3000.00000 -6019.65600 -2370.21500 6653.47500 3.58162
3200.00000 -5194.11300 5789.94800 6636.62700 3.81348
3400.00000 -4084.58200 -3236.51900 6613.98200 4.04674
3600.00000 -2749.50700 3792.25000 6586.73400 4.28177
3800.00000 -1260.82600 5677.02900 6556.35000 4.51888
4000.00000 299.14620 -4177.59700 6524.51100 4.75825
4200.00000 1841.68600 4245.86300 6493.02600 4.99998
4400.00000 3276.52700 990.17270 6463.74000 5.24398
4600.00000 4517.38300 -1451.66000 6438.41200 5.49007
4800.00000 5487.76900 3588.77100 6418.60400 5.73790
5000.00000 6126.65200 3327.23000 6405.56100 5.98700
5200.00000 6393.24700 -674.16460 6400.11700 6.23683
5400.00000 6270.39300 -1436.94300 6402.62200 6.48678
5600.00000 5766.00400 -5932.79600 6412.91500 6.73622
5800.00000 4912.40200 -5574.78600 6430.33500 6.98459
6000.00000 3763.61100 -109.16040 6453.78200 7.23137
6200.00000 2391.00900 -1646.15900 6481.79800 7.47618
6400.00000 877.94370 -6456.50100 6512.68200 7.71877
6600.00000 -686.12060 -6217.70800 6544.60700 7.95901
6800.00000 -2211.28600 2530.97000 6575.74000 8.19695
7000.00000 -3612.34500 3089.64900 6604.34000 8.43272
7200.00000 -4813.09400 -1149.51900 6628.85200 8.66659
7400.00000 -5749.79000 -4156.80000 6647.97500 8.89892
7600.00000 -6373.67300 -3862.74100 6660.71200 9.13013
7800.00000 -6652.72900 6121.71700 6666.41100 9.36071
8000.00000 -6572.76000 -3471.85100 6664.78000 9.59114
8200.00000 -6137.82900 4969.33100 6655.90400 9.82193
8400.00000 -5370.12700 6048.979000 6640.23500 10.05359
8600.00000 -4309.28600 5492.98900 6618.58100 10.28655
8800.00000 -3010.96300 -6381.13500 6592.07900 10.52123
9000.00000 -1544.80700 4957.36200 6562.14900 10.75793
9200.00000 8.51168 5168.19600 6530.44000 10.99688
9400.00000 1561.12700 1571.45800 6498.74500 11.23817
9600.00000 3022.65100 2749.46000 6468.90800 11.48176
9800.00000 4305.62200 6427.46800 6442.71400 11.72749
10000.00000 5331.29900 -106.48930 6421.77300 11.97503
由运行结果可知，输出结果共分三部分：第一部分是说明本次仿真的参数(在此例中Y(4)=0.00125)，第二部分是各状态变量及定义变量的最大值与最小值，第三部分是要求输出的四个变量随时间变化的表格。
第一章控制系统及仿真概述 控制系统的计算机仿真是一门涉及到控制理论、计算数学与计算机技术的综合性新型学科。这门学科的产生及发展差不多是与计算机的发明及发展同步进行的。它包含控制系统分析、综合、设计、检验等多方面的计算机处理。计算机仿真基于计算机的高速而精确的运算，以实现各种功能。 第一节控制系统仿真的基本概念 1．系统： 系统是物质世界中相互制约又相互联系着的、以期实现某种目的的一个运动整体，这个整体叫做系统。 “系统”是一个很大的概念，通常研究的系统有工程系统和非工程系统。 工程系统有：电力拖动自动控制系统、机械系统、水力、冶金、化工、热力学系统等。 非工程系统：宇宙、自然界、人类社会、经济系统、交通系统、管理系统、生态系统、人口系统等。 2．模型： 模型是对所要研究的系统在某些特定方面的抽象。通过模型对原型系统进行研究，将具有更深刻、更集中的特点。 模型分为物理模型和数学模型两种。数学模型可分为机理模型、统计模型与混合模型。 3．系统仿真： 系统仿真，就是通过对系统模型的实验，研究一个存在的或设计中的系统。更多的情况是指以系统数学模型为基础，以计算机为工具对系统进行实验研究的一种方法。 要对系统进行研究，首先要建立系统的数学模型。对于一个简单的数学模型，可以采用分析法或数学解析法进行研究，但对于复杂的系统，则需要借助于仿真的方法来研究。 那么，什么是系统仿真呢？顾名思义，系统仿真就是模仿真实的事物，也就是用一个模型（包括物理模型和数学模型）来模仿真实的系统，对其进行实验研究。用物理模型来进行仿真一般称为物理仿真，它主要是应用几何相似及环境条件相似来进行。而由数学模型在计算机上进行实验研究的仿真一般则称为数字仿真。我们这里讲的是后一种仿真。 数字仿真是指把系统的数学模型转化为仿真模型，并编成程序在计算机上投入运行、实验的全过程。通常把在计算机上进行的仿真实验称为数字仿真，又称计算机仿真。
航空工业管理学院 《电子信息系统仿真》课程设计 级电子信息工程专业班级 题目连续时间LTI系统仿真和时域分析 姓名学号 指导教师 二О一年月日
MATLAB软件简介 MATLAB 是MathWork 公司于1984 年推出的一套面向工程和科学运算的高性能软件，它具有强大的图形处理功能及符号运算功能，为我们实现信号的可视化及系统分析提供了强有力的工具。MATLAB 强大的工具箱函数可以分析连续信号、连续系统，同样也可以分析离散信号、离散系统，并可以对信号进行各种分析域计算，如相加、相乘、移位、反折、傅里叶变换、拉氏变换、Z 变换等等多种计算MATLAB 用于算法开发、数据可视化、数据分析以数值计算的高级技术计算语言和交互式环境，主要包括MATLAB和Simulink两大部分。MATLAB 可以进行矩阵运算、绘制函数和数据、实现算法、创建用户界面、连Matlab开发工作界面接其他编程语言的程序等，主要应用于工程计算、控制设计、信号处理与通讯、图像处理、信号检测、金融建模设计与分析等领域。Simulink是MATLAB最重要的组件之一，它提供一个动态系统建模、仿真和综合分析的集成环境。Simulink是MATLAB最重要的组件之一，它提供一个动态系统建模、仿真和综合分析的集成环境。在该环境中，无需大量书写程序，而只需要通过简单直观的鼠标操作，就可构造出复杂的系统。Simulink具有适应面广、结构和流程清晰及仿真精细、贴近实际、效率高、灵活等优点，并基于以上优点Simulink已被广泛应用于控制理论和数字信号处理的复杂仿真和设计。
设计目的 掌握信号经过LTI 系统的时域分析方法。根据连续时不变信号处理的基本概念、理论和方法对信号进行分析和处理，实现卷积积分或卷积和，零输入响应和零状态响应，学会应用MATLAB 对实际问题进行仿真，并对仿真结果进行分析。 在本次课程设计中，利用MATLAB 软件对LTI 连续系统时域进行仿真与分析。根据连续时不变信号处理的基本概念、理论和方法对信号进行分析和处理，实现卷积积分或卷积和，零输入响应和零状态响应，熟悉卷积和conv 函数，并会利用卷积求零状态响应，并对输出的波形和仿真结果进行分析。 理论分析 连续时间系统卷积分原理 连续时间信号1()f t 和2()f t 的卷积运算可用信号的分段求和来实现，即： 1212120 ()()*()()()lim ()()k f t f t f t f t f t d f k f t k ττ∞ ∞ -∞ ?→=-∞ ==-=?-???∑ ? 如果只求当t （n ）（n 为整数）时f (t )的值f (n ) ，则上式可得： 1212()()()()[()]k k f n f k f t k f k f n k ∞ ∞ =-∞ =-∞ ?= ?-???=??-?∑ ∑ 式中的12()[()]k f k f n k ∞ =-∞ ??-?∑ 实际上就是连续时间信号1()f t 和2()f t 经等 时间间隔均匀抽样的离散序列1()f k ?和2()f k ?的-。当 足够小时， ()f n ?就是卷积积分的结果——连续时间信号f (t )的较好数值近似。
第3章 连续系统仿真的方法 3.1 数值积分法 连续系统数值积分法，就是利用数值积分方法对广微分方程建立离散化形式的数学模型——差分方程，并求其数值解。可以想象在数学计算机上构造若干个数字积分器，利用这些数字积分器进行积分运算。在数字计算机上构造数字积分器的方法就是数值积分法，因而数字机的硬件特点决定了这种积分运算必须是离散和串行的。 把被仿真系统表示成一阶微分方程组或状态方程的形式。一阶向量微分方程及初值为 () (),00t Y Y t Y ???? ?????? Y =F = （3-1） 其中，Y 为n 维状态向量，F （t ，Y ）为n 维向量函数。 设方程（3-1）在011,,,,n n t t t t t +=…处的形式上的连续解为 ()()()()n+1n+1 t t n+10t t t =Y t +,(),n Y F t Y dt Y t F t Y dt =+ ?? （3-2） 设 n =() n Y Y t ，令 1n n n Y Y Q +=+ （3-3） 则有： ()1n+1t n Y Y += 也就是说， 1 (,)n n t n t Q F t Y dt +≈ ? （3-4） 如果n Y 准确解()n Y t 为近似值，n Q 是准确积分值的近似值，则式（3-4）
就是式（3-2）的近似公式。换句话说，连续系统的数值解就转化为相邻两个时间点上的数值积分问题。 因此，所谓数值解法，就是寻求初值问题（3-1）的真解在一系列离散点12n t t t <…<…上的近似解12,,,n Y Y Y ……，相邻两个时间离散点的间隔 1n n n t t +=-h ，称为计算步距或步长，通常取n =h h 为定值。可见，数值积分法的主要问题归结为对函数(,)F t y 的数值积分问题，即如何求出该函数定积分的近似解。为此，首先要把连续变量问题用数值积分方法转化成离散的差分方程的初值问题，然后根据已知的初值条件0y ，逐步地递推计算后续时刻的数值解(1,2,)i y i =…。所以，解初值问题的数值方法的共同特点是步进式的，采用不同的递推算法，就出现各种不同的数值积分方法。 3.2 替换法 基于数值积分的连续系统仿真方法具有成熟、计算精度比较高的优点，但算法公式比较复杂、计算量比较大，通常只有在对速度要求不高的纯数字仿真时使用。当进行实时仿真或在计算机控制系统中实现数字控制器的算法时，要求计算速度快，以便能在一个采样周期内完成全部计算任务，这就需要一些快速计算方法。 用数值积分方法在数字机上对一个连续系统进行仿真时，实际上已经进行了离散化处理，只不过在离散化过程中每一步都用到连续系统的模型，离散一步计算一步。那么，能否先对连续的模型进行离散化处理，得到一个“等效”的离散化模型，以后的每一步计算都直接在这个离散化模型基础上进行，而原来的连续数学模型不再参与计算呢？回答是肯定的。这些结构上比较简单的离散化模型，便于在计算机上求解，不仅用于连续系统数字仿真，而且也可用于数字控制器在计算机上实现。 替换法的基本思想是：对于给定的函数G （s ），设法找到s 域到z 域的的某种映射关系，它将S 域的变量s 映射到z 平面上，由此得到与连续系统传递函数G （s ）相对应的离散传函G （z ）。进而再根据G （z ）由z 反变换求的系统的时域离散模型——差分方程，据此便可以进行快速求解。
第二章 经典的连续系统仿真建模方法学 本章讨论经典的连续系统数字仿真的原理与方法，内容包括连续系统数字仿真的基本概念、经典的数值积分法、经典的线性多步法等。在数字计算机上进行连续系统仿真，首先要将连续模型离散化，因此，2.1节首先讨论离散化原理及要求，这是连续系统仿真的基础。然后，2.2节对经典的数值积分法----龙格-库塔法及其它典型的数值积分法仿真建模原理进行详细分析，并通过实例说明其应用要点；而2.3节对经典的线性多步法进行了介绍． 2.1 离散化原理及要求 在数字计算机上对连续系统进行仿真时，首先遇到的问题是如何解决数字计算机在数值及时间上的离散性与被仿真系统数值及时间上的连续性这一基本问题。 从根本意义上讲，数字计算机所进行的数值计算仅仅是“数字”计算，它表示数值的精度受限于字长，这将引入舍入误差；另一方面，这种计算是按指令一步一步进行的，因而，还必须将时间离散化，这样就只能得到离散时间点上系统性能。用数字仿真的方法对微分方程的数值积分是通过某种数值计算方法来实现的。任何一种计算方法都只能是原积分的一种近似。因此，连续系统仿真，从本质上是对原连续系统从时间、数值两个方面对原系统进行离散化，并选择合适的数值计算方法来近似积分运算，由此得到的离散模型来近似原连续模型。如何保证离散模型的计算结果从原理上确能代表原系统的行为，这是连续系统数字仿真首先必须解决的问题。 设系统模型为：),,(t u y f y = ，其中u (t )为输入变量，y (t )为系统变量；令仿真时间间隔为h ，离散化后的输入变量为)(?k t u ，系统变量为)(?k t y ，其中k t 表示t=kh 。如果)()(?k k t u t u ≈，)()(?k k t y t y ≈，即0)()(?)(≈-=k k k u t u t u t e ，0)()(?)(≈-=k k k y t y t y t e （对所有k=0,1,2,…），则可认为两模型等价，这称为相似原理（参见图2.1）。 实际上，要完全保证0)(,0)(==k y k u t e t e 是很困难的。进一步分析离散化引入的误差，随着计算机技术的发展，由计算机字长引入的舍入误差可以忽略，关键是数值积分算法，也称为仿真建模方法。相似原理用于仿真时，对仿真建模方法有三个基本要求： （1）稳定性：若原连续系统是稳定的，则离散化后得到的仿真模型也应是稳定的。关于稳定 性的详细讨论将在2.4节中进行。 （2）准确性：有不同的准确性评价准则，最基本的准则是： 图2.1 相 似 原 理
系统建模与仿真 开课对象：工业工程开课学期：6 学分：2学分；总学时：48学时；理论课学时：40学时； 实验学时：0 学时；上机学时：8学时 先修课程：概率论与数理统计 教材：系统建模与发展，齐欢，王小平编著，清华大学出版社，2004.7 参考书： 【1】离散事件系统建模与仿真，顾启泰，清华大学出版社 【2】现代系统建模与仿真技术，刘兴堂，西北工业大学出版社 【3】离散事件系统建模与仿真，王维平，国防科技大学出版社 【4】系统仿真导论，肖田元，清华大学出版社 【5】建模与仿真，王卫红，科学出版社 【6】仿真建模与分析(Simulaton Modeling and Analysis)(3rd eds.)，Averill M. Law, W.David Kelton，清华大学出版社/McGraw-Hill 一、课程的性质、目的和任务 建模与仿真是当代现代科学技术的主要内容，其技术已渗透到各学科和工程技术领域。本课程以一般系统理论为基础，让学生掌握适用于任何领域的建模与仿真的一般理论框架和基本方法。 本课程的目的和任务是使学生： 1．掌握建模基本理论； 2．掌握仿真的基本方法； 3．掌握一种仿真语言及仿真软件； 4．能够运用建模与仿真方法分析、解决工业工程领域的各种常见问题。 二、课程的基本要求 1．了解建模与仿真的作用和发展，理解组成要素。 2．掌握建模的几种基本方法，及模型简化的技术手段。 3．掌握建模的一般系统理论，认识随机数的产生的原因及统计控制方式。 4．能对离散事件进行仿真，并能分析运行结果。 三、课程的基本内容及学时分配 第一章绪论（3学时） 1．系统、模型、仿真的基本概念
机械系统建模与仿真 机理建模法 所谓机理模型，实际上就是采用由一般到特殊的推理理演绎方法，对巳知结构、参数的物理系统运用相应的物理定律或定理，经过合理分析简化建立起来描述系统各物理员动、静态变化性能的数学模型。 因此，机理建模法主要是通过理论分析推导方法建立系统模型。根据确定元件或系统行为所遵循的自然机理，如常用的物质不灭定律(用于液位、压力调节等)、能量守恒定律(用于温度调节等)、牛顿第二定律(用于速度、加速度调节等)、基尔霍夫定律(用于电气网络)等等，对系统各种运功规律的本质进行描述，包括质量、能量的变换和传递等过程，从而建立起变量间相互制约又相互依存的精确的数学关系。通常情况下，是给出微分方程形式或其派生形式——状态方程、传
递函数等。 实验建模法 所谓实验建模法，就是采用出特殊到一般的逻辑归纳方法，根据一定数量的在系统运行过程中的实测、观察的物理量数据，运用统计规律、系统辨识等理论合理估计出反映系统各物理量相互制约关系的数学模型。其主要依据是来自系统的大量实测数据．出此义称之为实验测定法。 当对所研究系统的内部结构和特性尚不清楚、甚至无法了解时,系统内部的机理变化规律就不能确定,通常称之为·黒箱”或“灰箱”问题,机理建模法也就无法应用。而根据所测到的系统输入输出数据,采用一定方法进行分析及处理来获得数学模型的统计模型法正好适应这种情况。通过对系统施加激励,观察和测取其响应,了解其内部変量的特性,并建立能近似反映同样变化的模拟系统的数学模型,就相当于建立起实际系统的数学描述(方程、曲线或图表等)。 （1）频率特性法 频率特性法是研究控制系统的一种应用广泛的工程实用方法。其特点在是通过建立系统频率响应与正弦输入信号之间的稳态特性关系，不仅可以反映系统的稳态性能，而且可以用来研究系统的稳定性和暂态性能;可以根据系统的开环频率特性，判别系统闭环后的各种性能；可以较方便地分析系统参数对动态性能的影响，并能大致指出改善系统性能的途径。 （2）系统辨识法 系统辨识法是现代控制理论与系统建模中常用的方法，它是依据测量到的输人与输出数据来建方静态与动态系统的数学模型，但其输出响应不局限于频率响应，阶压响应或脉冲响应等时间响应都可作为反映系统模型静态与动态特性的重要信息；而且，确定模型的过程更依赖于各种高效率的最优算法以及如何保证所测取数据的可靠性。因
课程设计 面向结构图的连续系统数题目 字仿真 学院计算机科学与信息工程学院 专业自动化 班级2010级2班 学生姓名小 指导教师吴诗贤 2013 年12 月20 日
面向结构图的连续系统数字仿真 姓名：陶园班级：10自动化3班学号：2010133330 摘要 根据自动控制系统中面向结构图的数字仿真的基本思想，探讨了仿真过程中典型环节的规范性、系统的连接矩阵、仿真求解、程序框图问题，并应用到实际的范例当中，并分析了结果总结了相关特点和相关结论。 自动控制系统常常是由许多环节组成的，要应用数字仿真方法对系统进行分析和研究，首先需要求出总的传递函数，再转化为状态空间表达式的形式，然后对其求解。当改变系统某一环节的参数时，尤其是要改变小闭环中某一环节的参数时，以上整个过程又需要重新计算，这对研究对象参数变化对整个控制系统的影响是十分不便的，为了克服这些缺点，同时大多数从事自动化工作的科技人员更习惯于用结构图的形式来分析和研究控制系统，为此产生了面向结构图的仿真方法。该方法只需将各个环节的参数及各环节间的连接方式输入计算机，仿真程序就能自动求出闭环系统的状态空间表达式。本课程设计主要介绍典型环节参数和连接关系构成闭环系统的状态方程的方法，而动态响应的计算，仍采用四阶龙格-库塔法。这种方法具有便于研究各个环节参数对系统的影响，并可以得到每个环节的动态响应，以及对多输入输出系统的进行仿真的有点。 关键字：结构图；典型环节；连接矩阵；数字仿真；
1、设计任务 已知某一系统结构如下图所示，编写matlab程序求a分别为2,4,6,8,10,12时输出量y的动态响应。 图1 2、需求分析及概要设计 2.1 需求分析 根据上述设计任务我们可以基本明确在我们课程设计当中应该明确以下几个方面： ?熟悉在数字计算机仿真技术中常用的四阶龙格-库塔算法。 ?明确在面向结构图的连续系统数字仿真，典型环节及其系数矩阵确定。 ?明确各连接矩阵的确定。 ?能够熟练运用MATLAB仿真软件。 2.2 设计思路 自动控制系统常常是由许多环节组成的，要应用数字仿真方法对系统进行分析和研究，首先需要求出总的传递函数，再转化为状态空间表达式的形式，然后对其求解。当改变系统某一环节的参数时，尤其是要改变小闭环中某一环节的参数时，以上整个过程又需要重新计算，这对研究对象参数变化对整个控制系统的影响是十分不便的，为了克服这些缺点，同时大多数从事自动化工作的科技人员更习惯于用结构图的形式来分析和研究控制系统，为此产生了面向结构图的仿真方法。该方法只需将各个环节的参数及各环节间的连接方式输入计算机，仿真程序就能自动求出闭环系统的状态空间表达式。以下是我们课程设计的主要设计思 图2
《建模与仿真》课程教学大纲 (Modeling and Simulation) 课程编码： 学分：2.5 总学时：40 适用专业：工业工程 先修课程：生产计划与控制、工程统计学、工程数学、运筹学、计算机编程技术 一、课程的性质、目的和任务 《建模与仿真》是面向工程实际的应用型课程，是工业工程系的主导课程之一。学生通过本课程的学习能够初步运用仿真技术来发现生产系统中的关键问题，并通过改进措施的实现，提高生产能力和生产效率。本课程的目的是要求学生通过学习、课堂教育和上机训练，能了解如何运用计算机仿真技术模拟生产系统的布置和调度管理。并熟悉和掌握计算机仿真软件的基本操作和能够实现的功能。使学生了解计算机仿真的基本步骤。结合本课程的特点，使学生掌握或提高系统化分析问题和解决问题的能力,为系统化管理生产打下基础。二、教学基本要求 具体在教学过程中要求学生应该达到： 1.全面了解本课程的性质与任务、框架内容以及理论和方法； 2.掌握仿真的概率统计基础知识。 3.掌握供理论模型建模方法。 4.掌握仿真模型的设计与实现方法。 5.熟练应用建模理论,对排队系统、库存系统、加工制造系统进行建模仿真。 三、教学内容与学时分配 离散事件系统仿真是仿真技术的重要领域，在规划论证、方案评估、计划调度、 加工制造、产品试验、生产培训、训练模拟、管理决策等方面得到广泛应用。本课程 深入地介绍了离散事件系统建模仿真的理论、方法和技术，突出对理论建模方法和计 算机实现技术的讲解，对离散事件系统建模仿真的发展和应用情况做了比较详尽的介 绍。 具体教学内容如下： 第一章绪论 4学时
本章分析了系统和制造系统定义、组成与特点，介绍了系统建模与仿真的基本概念和使用步骤，并给出应用案例。 本章教学目标： 本章教学基本要求： 了解常用术语及常用的仿真软件，了解仿真技术的的发展状况及应用。 理解系统与制造系统的定义及系统建模与仿真的概念及系统、模型与仿真之间的关系。 掌握制造系统建模与仿真的基本概念及基本步骤。 本章教学重点：制造系统建模与仿真的原则及基本步骤。 本章教学难点：制造系统建模与仿真的原则及基本步骤 第一节系统与制造系统 0.3学时 （一）什么是系统 （二）制造系统的组成与特点 第二节系统建模与仿真的基本概念。 0.3学时 (一)系统、模型与仿真的关系 (二)系统建模与仿真技术的特点 第三节制造系统建模与仿真的基本概念。 0.3学时 （一）制造系统建模与仿真的特点分析 （二）制造系统类型及建模元素 （三）制造系统仿真的功能分析 第四节系统建模与仿真的基本步骤 0.4学时 第五节系统建模与仿真的案例分析 0.5学时 （一）连杆生产线的组成与功能分析 （二）连杆生产线仿真模型的构建 （三）仿真逻辑的分析与定义 （四）仿真结果分析及系统优化 第二章系统建模与仿真的基本原理 2学时 本章在分析离散事件系统模型的分类和元素组成的基础上，介绍了建立系统模型的常用方法。 本章教学目标：使学生掌握常用的系统建模方法 本章教学基本要求：
制造系统建模与仿真学习心得 一、制造系统建模与仿真的含义 1．制造系统制造系统是制造过程及其所涉及的硬件、软件和人员所组成的一个将制造资源转变为产品或半成品的输入／输出系统，它涉及产品生命周期(包括市场分析、产品设计、工艺规划、加工过程、装配、运输、产品销售、售后服务及回收处理等)的全过程或部分环节。其中，硬件包括厂房、生产设备、工具、刀具、计算机及网络等；软件包括制造理论、制造技术(制造工艺和制造方法等)、管理方法、制造信息及其有关的软件系统等；制造资源包括狭义制造资源和广义制造资源；狭义制造资源主要指物能资源，包括原材料、坯件、半成品、能源等；广义制造资源还包括硬件、软件、人员等。随着科技的进步，制造系统的发展也经历了传统手工生产、机械化、自动化孤岛、集成制造、并行工程和敏捷制造等几个阶段。 2．模型与仿真模型是对真实对象和真实关系中那些有用的和让人感兴趣的特性的抽象，是对系统某些本质方面的描述。它以各种可用的形式描述被研究系统的信息。系统模型并不是对真实系统的完全复现，而是对系统的抽象，而仿真是通过对模型的实验以达到研究系统的目的，当制造系统尚未建立或者研究时间长成本高以及从安全性考虑我们有必要对制造系统预先进行建模并仿真以确定系统的最佳结构和配置方案、防止较大的经济损失、确定合理高效的作业计划，从而提高经济效益。 制造系统建模与仿真技术是以相似原理、模型理论、系统技术、信息技术以及建模与仿真应用领域的有关专业技术为基础，以计算机系统、与应用相关的物理效应设备及仿真器为工具，利用模型参与已有或设想的制造系统进行研究、分析、设计、加工生产、试验、运行、评估、维护、和报废(全生命周期)活动的一门多学科的综合性技术。 二、系统建模与仿真的发展及类型 1．系统建模与仿真的发展大致经历了这么几个阶段：1600—1940年左右，这一时期的建模仿真主要是在物理科学基础上的建模；20世纪40年代，由于电子计算机的出现，建模仿真技术开始飞速发展；20世纪50年代中期，建模仿真开始应用与航空领域；20世纪60年代，这一阶段主要是工业控制过程中的仿真；20世纪70年代，开始出现了包括经济、社会和环境因素的大系统仿真。到70年代中期，出现了系统与仿真的结合，如用于随机网络建模的SLAM仿真系统。在这一时期，系统仿真开始与更高级的决策结合，出现了决策支
电子科技大学中山学院学生实验报告 院别：电子信息学院 课程名称：信号与系统实验 一、实验目的 1.掌握连续系统Simulink 的建模方法。 2.掌握连续系统时域响应、频域响应的Simulink 仿真方法。 二、实验原理 连续系统的Simulink 仿真分析包括系统模型的创建和仿真分析两个过程。 利用Simulink 模块库中的有关功能模块创建的系统模型，主要有S 域模型、传输函数模型和状态空间模型等形式。 若将信号源子模块库（Sources ）中某种波形的信号源（如正弦或阶跃信号源）加于系统模型的输入端，则在系统模型的输出端用示波器观察零状态响应的 图1 系统时域响应Simulink 仿真的模型 以Sources 子模块库中的“lnl ”、Sinks 子模块中的“Outl ”分别作为系统模型的输入端和输出端，如图2所示。 ln1 out1 图2 系统响应Simulink 仿真的综合模型 建立图2形式系统模型并保存之后，利用如下响应的命令，可得到系统的 状态空间变量、频率响应曲线、单位阶跃响应和单位冲激响应的波形。 [A,B,C,D]=linmod(‘模型文件名’) %求状态空间矩阵，注意：‘模型文件名’不含扩展名 bode(A,B,C,D)；%绘制系统的频率特性曲线 bode(A,B,C,D,l u :ω :ωω:?1)； %绘制系统在10~ωω频率范围内，歩长为ω?的频率特性曲线；u i 为输入
端口编号，一般取1 Impulse(A,B,C,D) %绘制系统冲击响应的波形 Impulse(A,B,C,D,i u ,t : 1 :t t?) %绘制系统在时间范围内、歩长为的冲击 响应的波形 Step(A,B,C,D) %绘制系统阶跃响应的波形 Step(A,B,C,D,i u ,t : 1 :t t?) %绘制系统在 1 ~t t时间范围内、歩长为t?的 阶跃响应的波形 以上命令，可以逐条在命令窗口输入、执行，也可编写成M文件并运行，获得所需结果。 三、实验内容 1 、已知连续系统的系统函数为。用系统函数的形式建立系统 模型，进行Simulink仿真，（1）绘出阶跃响应波形（2）绘出系统的频率特性图。 2、已知连续系统的微分方程为 建立系统模型，进行Simulink仿真。（1）若f(t)=,绘出系统零状 态响应波形（2）分析系统的频率特性 3、线性系统如图17-13所示。要求：建立系统的S域模型，编写执行Simullink 仿真命令的M文件，求系统的状态空间变量，绘出系统的冲击响应波形和频率响应特性曲线。
《系统辨识》 实验手册 哈尔滨工业大学控制与仿真中心 2018年5月
目录 实验1 白噪声和M序列的产生---------------------------------------------------------- 2 实验2 脉冲响应法的实现---------------------------------------------------------------- 5 实验3 递推最小二乘法的实现---------------------------------------------------------- 9 附录实验报告模板---------------------------------------------------------------------- 13
实验1 白噪声、M 序列的产生 一、实验目的 1、熟悉并掌握产生均匀分布随机序列方法以及进而产生高斯白噪声方法 2、熟悉并掌握M 序列生成原理及仿真生成方法 二、实验原理 1、混合同余法 混合同余法是加同余法和乘同余法的混合形式，其迭代式如下： 11 1(*)mod /n n n n x a x b M R x M +++=+?? =? 式中a 为乘子，0x 为种子，b 为常数，M 为模。混合同余法是一种递归算法，即先提供一个种子0x ，逐次递归即得到一个不超过模M 的整数数列。 2、正态分布随机数产生方法 由独立同分布中心极限定理有：设随机变量12,,....,,...n X X X 相互独立，服从同一分布，且具有数学期望和方差： 2(),()0,(1,2,...)k k E X D X k μσ==>= 则随机变量之和1 n k i X =∑的标准化变量 : () n n n k k k X E X X n Y μ --= = ∑∑∑近似服从(0,1)N 分布。 如果n X 服从[0, 1]均匀分布，则上式中0.5μ=，2 1 12 σ= 。即 0.5n k X n Y -= ∑近似服从(0,1)N 分布。
系统建模仿真技术的历史现状和发展趋势分析 工程133 胡浩3130212026 【摘要】：经过半个多世纪的发展,仿真技术已经成为对人类社会发展进步具有重要影响的一门综合性技术学科。本文对建模与仿真技术发展趋势作了较全面分析。仿真建模方法更加丰富,更加需要仿真模型具有互操作性和可重用性,仿真建模VVA与可信度评估成为仿真建模发展的重要支柱;仿真体系结构逐渐形成标准,仿真系统层次化、网络化已成为现实,仿真网格将是下一个重要发展方向;仿真应用领域 更加丰富,向复杂系统科学领域发展,并将更加贴近人们的生活。 工程系统的仿真，起源于自动控制技术领域。从最初的简单电子、机械系统，逐步发展到今天涵盖机、电、液、热、气、电、磁等各个专业领域，并且在控制器和执行机构两个方向上飞速发展。 控制器的仿真软件，在研究控制策略、控制算法、控制系统的品质方面提供了强大的支持。随着执行机构技术的发展，机、电、液、热、气、磁等驱动技术的进步，以高可靠性、高精度、高反应速度和稳定性为代表的先进特征，将工程系统的执行品质提升到了前所未有的水平。相对控制器本身的发展，凭借新的加工制造技术的支持，执行机构技术的发展更加富于创新和挑战，而对于设计、制造和维护高性能执行机构，以及构建一个包括控制器和执行机构的完整的自动化系统也提出了更高的要求。 AMESIM软件正是能够提供平台级仿真技术的工具。从根据用户需求，提供液压、机械、气动等设计分析到复杂系统的全系统分析，
到引领协同仿真技术的发展方向，AMESIM的发展轨迹和方向代表了工程系统仿真技术的发展历程和趋势。 一、系统仿真技术发展的现状 工程系统仿真作为虚拟设计技术的一部分，与控制仿真、视景仿真、结构和流体计算仿真、多物理场以及虚拟布置和装配维修等技术一起，在贯穿产品的设计、制造和运行维护改进乃至退役的全寿命周期技术活动中，发挥着重要的作用，同时也在满足越来越高和越来越复杂的要求。因此，工程系统仿真技术也就迅速地发展到了协同仿真阶段。其主要特征表现为： 1、控制器和被控对象的联合仿真：MATLAB＋AMESIM，可以覆盖整个自动控制系统的全部要求。 2、被控对象的多学科、跨专业的联合仿真：AMESIM＋机构动力学＋CFD＋THERMAL＋电磁分析 3、实时仿真技术 实时仿真技术是由仿真软件与仿真机等半实物仿真系统联合实现的，通过物理系统的实时模型来测试成型或者硬件控制器。 4、集成进设计平台 现代研发制造单位，尤其是设计研发和制造一体化的大型单位，引进PDM/PLM系统已经成为信息化建设的潮流。在复杂的数据管理流程中，系统仿真作为CAE工作的一部分，被要求嵌入流程，与上下游工具配合。
实验一经典的连续系统仿真建模方法 一实验目的: 1 了解和掌握利用仿真技术对控制系统进行分析的原理和步骤。 2 掌握机理分析建模方法。 3 深入理解阶常微分方程组数值积分解法的原理和程序结构，学习用Matlab编写 数值积分法仿真程序。 4 掌握和理解四阶Runge-Kutta法，加深理解仿真步长与算法稳定性的关系。 二实验原理: 1非线性模型仿真 三实验内容: 1. 编写四阶 Runge_Kutta 公式的计算程序，对非线性模型（3）式进行仿真。 （1）将阀位u 增大10％和减小10％，观察响应曲线的形状; （2）研究仿真步长对稳定性的影响，仿真步长取多大时RK4 算法变得不稳定？ （3）利用 MATLAB 中的ode45()函数进行求解，比较与（1）中的仿真结果有何区别。 2. 编写四阶 Runge_Kutta 公式的计算程序，对线性状态方程（18）式进行仿真 （1）将阀位增大10％和减小10％，观察响应曲线的形状; （2）研究仿真步长对稳定性的影响，仿真步长取多大时RK4 算法变得不稳定？ （4）阀位增大10％和减小10％，利用MATLAB 中的ode45()函数进行求解阶跃响 应，比较与（1）中的仿真结果有何区别。 四程序代码: 龙格库塔: %RK4文件 clc close H=[1.2,1.4]';u=0.55; h=1; TT=[]; XX=[]; for i=1:h:200
k1=f(H,u); k2=f(H+h*k1/2,u); k3=f(H+h*k2/2,u); k4=f(H+h*k3,u); H=H+h*(k1+2*k2+2*k3+k4)/6; TT=[TT i]; XX=[XX H]; end; hold on plot(TT,XX(1,:),'--',TT,XX(2,:)); xlabel('time') ylabel('H') gtext('H1') gtext('H2') hold on 水箱模型: function dH=f(H,u) k=0.2; u=0.5; Qd=0.15; A=2; a1=0.20412; a2=0.21129; dH=zeros(2,1); dH(1)=1/A*(k*u+Qd-a1*sqrt(H(1))); dH(2)=1/A*(a1*sqrt(H(1))-a2*sqrt(H(2))); 三实验结果: 2编写四阶 Runge_Kutta 公式的计算程序，对线性状态方程（18）式进行仿真： 1 阀值u对仿真结果的影响 U=0.45;h=1; U=0.5;h=1;
新螺丝湾公交枢纽站 公交系统建模与仿真专业 姓名 学校 起止日期 报告提交日期
摘要 自新螺丝湾国际商贸城和新南部客运站建成并运营以来，该地区的人流量突增，同时给新螺丝湾公交枢纽站带来了新的挑战。因此，怎样能更好的做好该地区的公交车的调度和运营，满足乘客需求成了至关重要的问题。本课题在以满足乘客的需求提高该系统的效率为目标的前提下，采用了面向对象的离散虚拟仿真软件3D max和Flexsim,对昆明新螺丝湾公交枢纽站站内公交调度时间及库存时间进行了仿真系统研究;利用计算机仿真搭建了一个公交运输的仿真平台, 建立了昆明新螺丝湾公交枢纽站站内仿真模型，针对该系统在公交车调度和库存中存在的问题提出相应的改进方案，效率高,乘客满意度高,为公交公司在做公交调度方面的决策时提供有用方案。 关键词：新螺丝湾公交枢纽站;面向对象;flexsim;仿真模型
目录 摘要 1．绪言 (1) 2.理论分析 (1) 2.1数据采集方法： (1) 2.2数据样本分析： (1) 2.3数据分析方法： (1) 2.4建模方法 (2) 3.实验研究 (2) 3.1研究地点 (2) 3.2研究时间 (2) 3.3研究目的 (2) 3.4研究方法 (2) 4.数据分析 (2) 4.1新螺蛳湾公交调度分析 (2) 4.1.1 154路公交车调度分析 (2) 4.1.2 186路公交车调度分析 (5) 4.1.3 C7路公交车调度分析 (8) 4.2新螺蛳湾公交库存分析 (11) 4.2.1 154路公交车库存分析 (11) 4.2.2 186路公交车库存分析 (12)
4.2.3 C7路公交车库存分析 (13) 5.建模与仿真 (14) 5.1新螺蛳湾公交枢纽站布局图 (14) 5.2新螺蛳湾公交枢纽站公交车流程 (14) 5.3新螺蛳湾公交枢纽站设施的F LEXSIM实体映射对照 (15) 6.结果分析与结论 (15) 6.1154路公交 (15) 6.2186路公交 (16) 6.3C7路公交 (16) 7.感想 (16) 参考文献 (18)
电子科技大学中山学院学生实验报告 院别：电子信息学院 课程名称：信号与系统实验 一、实验目的 1.掌握连续系统的Simulnk 建模方法； 2.掌握连续系统时域响应、频域响应的Simulink 仿真方法。 二、实验原理 连续系统的Simulink 仿真分析包括系统模型的创建和仿真分析两个过程。 利用Simulink 模块库中的有关功能模块创建的系统模型，主要有S 域模型、传输函数模型和状态空间模型等形式。 若将信号源子模块库（Sources ）中某种波形的信号源（如正弦或阶跃信号源）加于系统模型的输入端， 图1 系统时域响应Simulink 仿真的模型 以Sources 子模块库中的”lnl ”、Sinks 子模块中的”Outl ”分别作为系统模型的输入端和输出端，如图2所示。 ln1 out1 图2 系统响应Simulink 仿真的综合模型 建立图2形式系统模型并保存之后，利用如下响应的命令，可得到系统的状态空间变量、频率响应曲线、单位冲击响应的波形。 [A,B,C,D]=linmod(‘模型文件名’) %求状态空间矩阵。注意：‘模型文件名’不含扩展名 bode(A,B,C,D)； %绘制系统的频率特性曲线 bode(A,B,C,D,l u :ω0:ωω:?1)； %绘制系统在10~ωω频率范围内，歩长为ω?的频率特性曲线；u i 为输入端口编号，一般取1 Impulse(A,B,C,D) %绘制系统冲击响应的波形 Impulse(A,B,C,D,i u ,t 0:1:t t ?) %绘制系统在时间范围内、歩长为的冲击响应的波形
Step(A,B,C,D) %绘制系统阶跃响应的波形 Step(A,B,C,D,i u ,t : 1 :t t?) %绘制系统在 1 ~t t时间范围内、歩长为t?的阶跃响应的波形 以上命令，可以逐条在命令窗口输入、执行，也可编写成M文件并运行，获得所需结果。 三、实验内容 （题目）3、线性系统如图17-13所示。要求：建立系统的S域模型，编写执行Simullink仿真命令的M文件，求系统的状态空间变量，绘出系统的冲击响应波形和频率响应特性曲线。 解答：（模型图） （m文件） [A,B,C,D]=linmod('OuDi') figure(1) impulse(A,B,C,D,1,0:0.01:100);grid,ylabel('12无线,欧迪,33') figure(2) step(A,B,C,D,1,0:0.01:100);grid,ylabel('12无线,欧迪,33') figure(3) bode(A,B,C,D,1,1:0.01:10);grid,ylabel('12无线,欧迪,33')
系统建模与仿真教学大纲 课程名称：系统建模与仿真课程编号： 英文名称：System Modeling and Simulation 学时：64 学分：3.5 适用专业：工业工程课程类别：必修 课程性质：学科基础课 先修课程：工程数学、运筹学、统计学、计算机编程技术 教材：《离散事件系统仿真》，Jerry Banks等著，肖田元等译，机械工业出版社，2007.7 一、本课程的性质与任务 《系统建模与仿真》是面向工程实际的应用型课程，是工业工程系的主导课程之一。学生通过本课程的学习能够初步运用仿真技术来发现生产系统中的关键问题，并通过改进措施的实现，提高生产能力和生产效率。 二、课程教学的基本要求： 本课程以制造型生产企业为核心，阐述了离散事件系统建模与仿真技术在生产企业分析中的基本原理和方法。其内容涉及计算机仿真技术在生产系统分析中的作用和原理、仿真软件的介绍，重点介绍排队系统、库存系统、加工系统以及输入、输出数据分析。本课程的目的是要求学生通过学习、课堂教育和上机训练，能了解如何运用计算机仿真技术模拟生产系统的布置和调度管理；并熟悉和掌握计算机仿真软件的基本操作和能够实现的功能；使学生了解计算机仿真的基本步骤。 三、课程内容及教学要求 第一章绪论
教学基本内容: 生产系统的基本特征、生产系统仿真的基本概念、生产系统仿真模型的建立思路、以及生产系统仿真研究的步骤。 重点： 系统、系统模型、系统仿真等建模与仿真相关的基本概念；离散系统与连续系统的区别；生产系统建模的方法与仿真研究的步骤。 难点： 系统、系统模型、系统仿真等建模与仿真相关的基本概念；离散系统与连续系统的区别。 教学基本要求：了解生产系统的基本特征；理解掌握系统、系统模型、系统仿真等建模与仿真相关的基本概念；了解系统仿真的类型；理解离散系统与连续系统的区别；熟悉生产系统建模的方法与仿真研究的步骤。 第二章生产仿真用概率统计 教学基本内容: 介绍随机变量、概率函数、随机数；均匀的连续分布随机数及其生成；各种离散分布随机数的产生；非均匀的连续分布随机数及其产生。 重点： 随机变量的数字特征；产生连续均匀分布随机数的几种方法；计算机产生随机数的方法；随机数的统计检验；各种离散分布随机数的产生方法；非均匀连续分布随机数产生方法。 难点： 各种随机数的产生方法。 教学基本要求：理解掌握随机变量（离散、连续），以及连续随机变量的密度函数的概念；掌握随机变量的数字特征；理解掌握随机数的概念，熟悉产生连续均匀分布随机数的几种方法，掌握计算机产生随机数的方法；熟悉随机数的统计检验；掌握各种离散分布随机数的产生方法；熟悉非均匀连续分布随机数及其产生方法。 第三章排队系统 教学基本内容: 介绍排队论的基本概念；到达时间间隔和服务时间的分布；排队系统的特征参数与分析方法；排队系统的仿真及分析。
第一篇 连续系统仿真篇 第一章 连续系统仿真概论 按系统模型的特征分类，可以有连续系统仿真及离散事件系统仿真两大类。本篇讨论连续系统仿真问题。过程控制系统、调速系统、随动系统等这类系统称作连续系统，它们共同之处是系统状态变化在时间上是连续的，可以用方程式（常微分方程、偏微分方程、差分方程）描述系统模型。 1.1 连续系统模型描述 连续系统仿真中的数学模型有很多种，但基本上可分为三类：连续时间模型、离散时间模型及连续－离散混合模型。本节将对它们的线性定常形式作一介绍，1.2节将介绍几种模型结构变换的方法。 1.1.1 连续时间模型 如果一个系统的输入量u(t)，输出量y(t)，系统的内部状态变量x(t)都是时间的连续函数，那么我们可以用连续时间模型来描述它。系统的连续时间模型通常可以有以下几种表示方式：常微分方程，传递函数，权函数，状态空间描述．本节仅对其一般描述形式作一简要介绍，偏微分方程将在第8章讨论。 1. 常微分方程 常微分方程可用(1.1)式表示： u c dt u d c dt u d c y a dt dy a dt y d a dt y d n n n n n n n n n n n +++=++++------- 1221111111 (1.1) 其中n 为系统的阶次，a i n i (,,,,)=012 为系统的结构参数，),,2,1(n j c j =为输入函数的结构参数，它们均为实常数。 2. 传递函数 若系统的初始条件为零，即系统在 t=0 时已处于一个稳定状态，那么对（1.1）式两边 取拉氏变换后可得： ) ()()() ()()()(2 21 1111s U c s U s c s U s c s Y a s sY a s Y s a s Y s n n n n n n n +++=++++---- 稍加整理，并记： j n j j n n j j j n s a s c s U s Y s G ∑∑=--=-== 1 0) () ()( (1.2) (1.2)式称为系统的传递函数。 3. 权函数