排队论模型的计算机模拟数值模拟是依据被模拟对象的数学或逻辑模型,
利用计算机进行实验的一种技术,已成为与理论分析,实验室实验并列的重要研究方法,如飞机船舶的设计等大型项目,美国军队的计算机推演,排队论模型的计算机模拟,目的是研究系统性状随时间的发展,特别是它们能否具有某种稳定的时间平均性质,研究系统对参数的敏感性以及系统的优化与设计,由于它涉及从一定的概率分布抽取若干随机变量的值,因而是随机模拟,属 Monte
Carlo模拟的范畴,
一,实例考虑一小型机械加工车间,该车间加工四种不同类型的零件,零件相继到达的时间间隔 T是随机的,无妨设其按照 P(T<=x)=F(x)分布 ;每次来到的零件类型也是随机的,假设第 i种类型的零件出现的概率为 Pi.车间里有三台机床,每台机床可以加工任何一种零件,零件的加工时间由类型决定,是确定的,规则,如果零件到达车间时有空闲机床,则该零件被立即加工,否则需排队等待 ;先来到的零件排在前面,当机床加工完零件时,如果有零件在等待加工,则该机床立即加工排在队列前面的零件,加工完的零件分送不同部门,不再跟踪,但记录下来,
1.系统图像首先,我们用一组数来描述系统在任何时刻所处的状态,这组数称为 系统图像,随着时间的发展,
系统的状态不断变化,相应地我们只有修改系统图像即可,
我们以下表为例,设目前的时间是 2000(分 ).
零件类型 加工时间 到达时间 下一事件时刻下一零件 3 75 2002 2002
- - - -
排队零件 1 52 1992 -
4 43 1976 -
4 43 1972 2040
加工零件 2 21 1936 2017
3 75 1896 2003
现在时间 2000
已加工数 (1)33 (2)14 (3)24 (4)22
下一零件的信息可由 F(x) 及 Pi(x)
(i=1,2,3,4)对进行随机抽样得到,
表 1
零件类型 加工时间 到达时间 下一事件时刻下一零件 3 75 2002 2002
- - - -
排队零件 1 52 1992 -
4 43 1976 -
4 43 1972 2040
加工零件 2 21 1936 2017
3 75 1896 2003
现在时间 2001
已加工数 (1)33 (2)14 (3)24 (4)22
时间过了 1分钟,则图像基本不变,就是将现在时间调为 2001.
2.过程模拟从 2000分的初始表格,即系统图像开始,我们来进行过程模拟,先查看所以下一步可能事件发生的时间,由于加工零件的结束时间在表中是按顺序排列的,因此我们只需要考虑表中被加工零件的最后一行,并与第一行中的下一零件的达到时刻比较,看看那个更早,
零件类型 加工时间 到达时间 下一事件时刻下一零件 2 21 2018 2018
3 75 2002 -
排队零件 1 52 1992 -
4 43 1976 -
4 43 1972 2040
加工零件 2 21 1936 2017
3 75 1896 2003
现在时间 2002
已加工数 (1)33 (2)14 (3)24 (4)22
下一零件的信息可由 F(x) 及 Pi(x)
(i=1,2,3,4)对进行随机抽样得到,
表 2
注,模拟程序交替地在处理系统图像和计算抽样值的子程序间运行,
从表 2可知,下一事件是加工完毕一个 3型零件,
将时钟拨到 2003分,将这个 3型零件移出系统,并在表的最后一行的统计数据中将 3型零件的数量加 1,得到表 3.
零件类型 加工时间 到达时间 下一事件时刻下一零件 2 21 2018 2018
- - - -
排队零件 3 75 2002 -
1 52 1992 -
4 43 1976 2046
加工零件 4 43 1972 2040
2 21 1936 2017
现在时间 2003
已加工数 (1)33 (2)14 (3)25 (4)22
表 3
由表可知,下一事件的时间是 2017.
零件类型 加工时间 到达时间 下一事件时刻下一零件 2 21 2018 2018
- - - -
排队零件 - - - -
3 75 2002 -
1 52 1992 2069
加工零件 4 43 1976 2046
4 43 1972 2040
现在时间 2017
已加工数 (1)33 (2)15 (3)25 (4)22
表 4
由表可知,下一事件的时间是 2018.我们再得到新的表格,即系统图像,
3.统计量的计算模拟的目的是为了得到系统的统计性质,本例中我们只按类型统计了加工完毕的零件数,在不同的问题中,根据模拟的实际目的,可以得到各种有关参数的点估计或区间估计,如一个顾客的平均服务时间,平均等待时间,平均队列的长度等等,
初值的选择,许多时候,模拟的目的是为了得到系统平衡是的统计参数,在初始阶段往往受到初始值的影响,因此我们舍弃一适当长度的初始模拟数据,当然,若我们选择了恰当的初值,便不必舍弃了,
4.表处理技术在处理排队论模拟中的表格时,形成了一定的技巧,主要是使用 指针 来处理链接表格中的每一条记录,表格中的一行表示一个零件的相关信息,
我们称之为一条记录,在计算机内存中,每条记录用固定长度的若干连续单元存放,而系统图像中的连续两条记录则不必连续放置,因为在每条记录所使用的单元中其实存放了两部分的信息,
一是模拟信息,二是加入一个指针变量的值,用以表示此记录的下一记录的首地址,当然在表头
（第一条记录）、表末 (最后一条记录 )都要引入特殊标志,这样就形成一个链表,对记录的处理时,我们就只要对指针进行操作即可,
二,随机数的生成在排队模型的数值模拟中,必须按照各种给定的概率分布,得到所需要的随机数,在上面的实例中,按照给定的分布,决定下一零件的到来时间与类型,但是在计算机上,所有的这些数据都是按照一定的算法生成的,因此不是真正随机的,但是,只要其统计特性与真正的随机数很相似,就可以达到我们的目的,我们称它们为 伪随机数,下面我们来介绍 (0,1)区间上的均匀分布,
一般的伪随机数如不特别指明分布时,就是专指这种分布,
均匀分布伪随机数的生成
1.混同余法 适当取定非负整数 a,c,m及 x0,称
x0为种子,按照以下公式生成一个整数序列,
./
.,,2,1,0)(m od1
mxu
nimcaxx
ii
ii
令
则序列 ui即给出了 (0,1)区间上的伪随机数,容易看出,这样产生的伪随机数序列有一个周期,设其长度伪 p,最大可能的周期是 p=m.周期为 m的伪随机数序列称为具有 完全周期,利用数论中的知识来研究它,
现状,当 m=235时,取 a=27+1=129,c=1.
2,乘同余法 混同余法中包含乘法和加法两种运算,下面介绍的算法仅包含乘法,故名,
取定非负整数 a,m及 x0,令
,)( m o d1 maxx ii
仍以序列 ui为 (0,1)区间上的伪随机数,一般来说,乘同余法给出的序列不能达到完全周期 ;
但是,当 m与 x0互素,a满足一定的条件时,可达到一个最大的周期,非常接近 m,
m=235时,取 a=217+3,x0为奇数,
3,伪随机数在微机上的生成算法上面介绍的两种算法中的运算都是整数运算,
不能进行舍入,而现在的微机的整数不能超过 15个二进位 (-32767~32768),因此上面的算法不能在微机上实现,许多人想方设法来克服它,下面介绍 Wichmann和 Hill的公式,
),30323(m o d170
),30307(m o d172),30269(m o d171
1
11
ii
iiii
zz
yyxx
,30323/30307/30269/ iiii zyxT取则伪随机数为 ).(
iii TTR T R U N
从均匀分布随机数生成其他分布的随机数有各种方法,如变换法,舍选法等待,我们不再介绍,
三,过程模拟的更新方法我们在实例中介绍的方法可以视为一次随机实验,目的是得到一些统计推断,但是,由于相继时刻的系统状态是高度相关的,这样就不利于统计分析,还有,传统的方法不能回答应该模拟多长时间以及得到的结果精度如何等问题,
克服这一困难的方法是 更新模拟,
.,
.
,,
),,2,1(
,0
},0);({,
1
210
服从相同的分布度个循环的长为第称即构成一次次的循环一的概率规律所支配而且过去与未来都被同的历史立于过去过程未来的统计性质独对任何随机时间系列由过程本身决定的如果存在一考虑一个随机过程粗略地说
i
iii
i
l
iTTl
iT
TTT
ttx
本节前面介绍的例子中,我们假设输入过程是时齐的,即转移矩阵 Q=(pij)与时间无关,设初始时刻为三台机床都闲置,没有零件等待,当恰有一个零件到来的瞬间我们记为 T0,那么今后所以的这种时刻都是更新时刻,因此模拟更新法的思想是,我们将实验看出是重复实验,
)),(()(,,
)),(()(,
1
1
xfExfN
xfExfxjx
N
j
j
N
j
jj
时一般我们有当然取的取值为在时刻设设 x(t)为一随时间变化的描述系统状态的随机过程,如队列长度,我们关心的是函数 f(x)
的均值,则一般的处理方法是,
对一个具有更新性质的排队系统,我们将时间取离散的整数值,Ti为更新时刻,
.
)(
)(
))((
.}1),,{(,
.,)(
1
1
1
1
1
E
yE
xfEr
iy
TTxfy
ii
iii
T
Tj
ji
i
i
且量是独立同分布的随机向可知由更新过程的性质令此式给出了 r的点估计的令一种方式,但是现在我们可以估计概率,
.),cov (),va r(
),va r(),cov (2)va r(
,0)(,
22
分别表示方差和协方差这里的方差是且则令
iiii
iiiii
ryry
zzEryz
.),1,0(/)(
,,
1
,
1
11
nNryn
n
y
n
y
n
i
i
n
i
i
则由中心极限定理知令由于 z的方差实际上不知道,因此我们用样本估计值来代替,这样
.))((
1
1
,)(
1
1
,)(
1
1
,2,/?
.),1,0(
/
)?(
1
12
1
2
22
1
2
11
22
2
1211
2
n
i
ii
n
i
i
n
i
i
yy
n
s
n
syy
n
s
srsrssyr
nN
s
rrn
其中
).?,?(
1
n
sZ
r
n
sZ
r
r
的置信区间为时如度为从而我们可以求得置信
利用计算机进行实验的一种技术,已成为与理论分析,实验室实验并列的重要研究方法,如飞机船舶的设计等大型项目,美国军队的计算机推演,排队论模型的计算机模拟,目的是研究系统性状随时间的发展,特别是它们能否具有某种稳定的时间平均性质,研究系统对参数的敏感性以及系统的优化与设计,由于它涉及从一定的概率分布抽取若干随机变量的值,因而是随机模拟,属 Monte
Carlo模拟的范畴,
一,实例考虑一小型机械加工车间,该车间加工四种不同类型的零件,零件相继到达的时间间隔 T是随机的,无妨设其按照 P(T<=x)=F(x)分布 ;每次来到的零件类型也是随机的,假设第 i种类型的零件出现的概率为 Pi.车间里有三台机床,每台机床可以加工任何一种零件,零件的加工时间由类型决定,是确定的,规则,如果零件到达车间时有空闲机床,则该零件被立即加工,否则需排队等待 ;先来到的零件排在前面,当机床加工完零件时,如果有零件在等待加工,则该机床立即加工排在队列前面的零件,加工完的零件分送不同部门,不再跟踪,但记录下来,
1.系统图像首先,我们用一组数来描述系统在任何时刻所处的状态,这组数称为 系统图像,随着时间的发展,
系统的状态不断变化,相应地我们只有修改系统图像即可,
我们以下表为例,设目前的时间是 2000(分 ).
零件类型 加工时间 到达时间 下一事件时刻下一零件 3 75 2002 2002
- - - -
排队零件 1 52 1992 -
4 43 1976 -
4 43 1972 2040
加工零件 2 21 1936 2017
3 75 1896 2003
现在时间 2000
已加工数 (1)33 (2)14 (3)24 (4)22
下一零件的信息可由 F(x) 及 Pi(x)
(i=1,2,3,4)对进行随机抽样得到,
表 1
零件类型 加工时间 到达时间 下一事件时刻下一零件 3 75 2002 2002
- - - -
排队零件 1 52 1992 -
4 43 1976 -
4 43 1972 2040
加工零件 2 21 1936 2017
3 75 1896 2003
现在时间 2001
已加工数 (1)33 (2)14 (3)24 (4)22
时间过了 1分钟,则图像基本不变,就是将现在时间调为 2001.
2.过程模拟从 2000分的初始表格,即系统图像开始,我们来进行过程模拟,先查看所以下一步可能事件发生的时间,由于加工零件的结束时间在表中是按顺序排列的,因此我们只需要考虑表中被加工零件的最后一行,并与第一行中的下一零件的达到时刻比较,看看那个更早,
零件类型 加工时间 到达时间 下一事件时刻下一零件 2 21 2018 2018
3 75 2002 -
排队零件 1 52 1992 -
4 43 1976 -
4 43 1972 2040
加工零件 2 21 1936 2017
3 75 1896 2003
现在时间 2002
已加工数 (1)33 (2)14 (3)24 (4)22
下一零件的信息可由 F(x) 及 Pi(x)
(i=1,2,3,4)对进行随机抽样得到,
表 2
注,模拟程序交替地在处理系统图像和计算抽样值的子程序间运行,
从表 2可知,下一事件是加工完毕一个 3型零件,
将时钟拨到 2003分,将这个 3型零件移出系统,并在表的最后一行的统计数据中将 3型零件的数量加 1,得到表 3.
零件类型 加工时间 到达时间 下一事件时刻下一零件 2 21 2018 2018
- - - -
排队零件 3 75 2002 -
1 52 1992 -
4 43 1976 2046
加工零件 4 43 1972 2040
2 21 1936 2017
现在时间 2003
已加工数 (1)33 (2)14 (3)25 (4)22
表 3
由表可知,下一事件的时间是 2017.
零件类型 加工时间 到达时间 下一事件时刻下一零件 2 21 2018 2018
- - - -
排队零件 - - - -
3 75 2002 -
1 52 1992 2069
加工零件 4 43 1976 2046
4 43 1972 2040
现在时间 2017
已加工数 (1)33 (2)15 (3)25 (4)22
表 4
由表可知,下一事件的时间是 2018.我们再得到新的表格,即系统图像,
3.统计量的计算模拟的目的是为了得到系统的统计性质,本例中我们只按类型统计了加工完毕的零件数,在不同的问题中,根据模拟的实际目的,可以得到各种有关参数的点估计或区间估计,如一个顾客的平均服务时间,平均等待时间,平均队列的长度等等,
初值的选择,许多时候,模拟的目的是为了得到系统平衡是的统计参数,在初始阶段往往受到初始值的影响,因此我们舍弃一适当长度的初始模拟数据,当然,若我们选择了恰当的初值,便不必舍弃了,
4.表处理技术在处理排队论模拟中的表格时,形成了一定的技巧,主要是使用 指针 来处理链接表格中的每一条记录,表格中的一行表示一个零件的相关信息,
我们称之为一条记录,在计算机内存中,每条记录用固定长度的若干连续单元存放,而系统图像中的连续两条记录则不必连续放置,因为在每条记录所使用的单元中其实存放了两部分的信息,
一是模拟信息,二是加入一个指针变量的值,用以表示此记录的下一记录的首地址,当然在表头
（第一条记录）、表末 (最后一条记录 )都要引入特殊标志,这样就形成一个链表,对记录的处理时,我们就只要对指针进行操作即可,
二,随机数的生成在排队模型的数值模拟中,必须按照各种给定的概率分布,得到所需要的随机数,在上面的实例中,按照给定的分布,决定下一零件的到来时间与类型,但是在计算机上,所有的这些数据都是按照一定的算法生成的,因此不是真正随机的,但是,只要其统计特性与真正的随机数很相似,就可以达到我们的目的,我们称它们为 伪随机数,下面我们来介绍 (0,1)区间上的均匀分布,
一般的伪随机数如不特别指明分布时,就是专指这种分布,
均匀分布伪随机数的生成
1.混同余法 适当取定非负整数 a,c,m及 x0,称
x0为种子,按照以下公式生成一个整数序列,
./
.,,2,1,0)(m od1
mxu
nimcaxx
ii
ii
令
则序列 ui即给出了 (0,1)区间上的伪随机数,容易看出,这样产生的伪随机数序列有一个周期,设其长度伪 p,最大可能的周期是 p=m.周期为 m的伪随机数序列称为具有 完全周期,利用数论中的知识来研究它,
现状,当 m=235时,取 a=27+1=129,c=1.
2,乘同余法 混同余法中包含乘法和加法两种运算,下面介绍的算法仅包含乘法,故名,
取定非负整数 a,m及 x0,令
,)( m o d1 maxx ii
仍以序列 ui为 (0,1)区间上的伪随机数,一般来说,乘同余法给出的序列不能达到完全周期 ;
但是,当 m与 x0互素,a满足一定的条件时,可达到一个最大的周期,非常接近 m,
m=235时,取 a=217+3,x0为奇数,
3,伪随机数在微机上的生成算法上面介绍的两种算法中的运算都是整数运算,
不能进行舍入,而现在的微机的整数不能超过 15个二进位 (-32767~32768),因此上面的算法不能在微机上实现,许多人想方设法来克服它,下面介绍 Wichmann和 Hill的公式,
),30323(m o d170
),30307(m o d172),30269(m o d171
1
11
ii
iiii
zz
yyxx
,30323/30307/30269/ iiii zyxT取则伪随机数为 ).(
iii TTR T R U N
从均匀分布随机数生成其他分布的随机数有各种方法,如变换法,舍选法等待,我们不再介绍,
三,过程模拟的更新方法我们在实例中介绍的方法可以视为一次随机实验,目的是得到一些统计推断,但是,由于相继时刻的系统状态是高度相关的,这样就不利于统计分析,还有,传统的方法不能回答应该模拟多长时间以及得到的结果精度如何等问题,
克服这一困难的方法是 更新模拟,
.,
.
,,
),,2,1(
,0
},0);({,
1
210
服从相同的分布度个循环的长为第称即构成一次次的循环一的概率规律所支配而且过去与未来都被同的历史立于过去过程未来的统计性质独对任何随机时间系列由过程本身决定的如果存在一考虑一个随机过程粗略地说
i
iii
i
l
iTTl
iT
TTT
ttx
本节前面介绍的例子中,我们假设输入过程是时齐的,即转移矩阵 Q=(pij)与时间无关,设初始时刻为三台机床都闲置,没有零件等待,当恰有一个零件到来的瞬间我们记为 T0,那么今后所以的这种时刻都是更新时刻,因此模拟更新法的思想是,我们将实验看出是重复实验,
)),(()(,,
)),(()(,
1
1
xfExfN
xfExfxjx
N
j
j
N
j
jj
时一般我们有当然取的取值为在时刻设设 x(t)为一随时间变化的描述系统状态的随机过程,如队列长度,我们关心的是函数 f(x)
的均值,则一般的处理方法是,
对一个具有更新性质的排队系统,我们将时间取离散的整数值,Ti为更新时刻,
.
)(
)(
))((
.}1),,{(,
.,)(
1
1
1
1
1
E
yE
xfEr
iy
TTxfy
ii
iii
T
Tj
ji
i
i
且量是独立同分布的随机向可知由更新过程的性质令此式给出了 r的点估计的令一种方式,但是现在我们可以估计概率,
.),cov (),va r(
),va r(),cov (2)va r(
,0)(,
22
分别表示方差和协方差这里的方差是且则令
iiii
iiiii
ryry
zzEryz
.),1,0(/)(
,,
1
,
1
11
nNryn
n
y
n
y
n
i
i
n
i
i
则由中心极限定理知令由于 z的方差实际上不知道,因此我们用样本估计值来代替,这样
.))((
1
1
,)(
1
1
,)(
1
1
,2,/?
.),1,0(
/
)?(
1
12
1
2
22
1
2
11
22
2
1211
2
n
i
ii
n
i
i
n
i
i
yy
n
s
n
syy
n
s
srsrssyr
nN
s
rrn
其中
).?,?(
1
n
sZ
r
n
sZ
r
r
的置信区间为时如度为从而我们可以求得置信
