数学建模社区-数学中国

标题: 排队论模型(八):Matlab 生成随机数、排队模型的计算机模拟 [打印本页]

作者: 浅夏110    时间: 2020-6-13 09:35
标题: 排队论模型(八):Matlab 生成随机数、排队模型的计算机模拟
1 产生给定分布的随机数的方法
3 m: z" R* M/ D; c5 f) G; J. lMatlab 可以产生常用分布的随机数。下面我们介绍按照给定的概率分布产生随机数的一般方法,这些方法都以U(0,1) 分布的随机变量为基础。
2 q( t) ?/ V9 k# I- M; V" \0 V4 P1 ~* |% P7 {: j; O% j
(i)反变换法
( W( u6 b! D7 u定理 设 X 是一个具有连续分布函数 F(x) 的随机变量,则 F(X ) 在 [0,1] 上服 从均匀分布。
3 B  {  [8 x; j9 `. K3 e0 y2 ]! x, y! ?! W. f7 ^

8 X9 S+ ?1 V3 n
8 j$ a4 |1 f) C6 [3 r' x% B# h3 t* K4 c* s2 @6 ?

, c" ~* g  [4 N1 U7 t" i1 R& o. N7 ^(ii)卷积法% I  J7 ~4 N1 c( s) a: @$ G

/ ]: Z) S2 D: Q# p9 {
. V: T8 `* A1 S
/ J  G# i6 m/ _1 t(iii)取舍法9 S+ ]1 k! d8 W+ [6 w7 c$ p, M
若随机变量 X 在有限区间(a,b) 内变化,但概率密度 f (x)具有任意形式(甚至没 有解析表达式),无法用前面的方法产生时,可用取舍法。一种比较简单的取舍法的步 骤是:0 u, F* D; g7 u8 g1 V* b

" X2 ~' ~; b" O# u; z4 V8 a& ]; P. z: O/ _  t7 t6 t
9 h  l/ {  d6 ?  H" e7 \1 Y0 h* N# g) c
2 排队模型的计算机模拟
$ W7 r- P8 U9 {0 g- p$ r& ~5 E8 h2.1 确定随机变量概率分布的常用方法
' L+ ]7 X7 J, b% K: y2 W, c# U在模拟一个带有随机因素的实际系统时,究竟用什么样的概率分布描述问题中的随 机变量,是我们总是要碰到的一个问题,下面简单介绍确定分布的常用方法:
6 ^9 K4 e) Y$ w' p- [
5 Z+ ^+ m: e' X* Y【1 】根据一般知识和经验,可以假定其概率分布的形式,如顾客到达间隔服从指数 分布 Exp(λ) ;产品需求量服从正态分布   ;订票后但未能按时前往机场登机 的人数服从二项分布 B(n, p) 。然后由实际数据估计分布的参数 λ,μ,σ 等,参数估计 可用极大似然估计、矩估计等方法。& x% b" d1 ?8 c/ m) B7 L
7 q9 T6 l+ O  U: V5 c
【2】 直接由大量的实际数据作直方图,得到经验分布,再通过假设检验,拟合分布 函数,可用  检验等方法。 3 o 既缺少先验知识,又缺少数据时,对区间(a,b) 内变化的随机变量,可选用 Beta 分布(包括均匀分布)。先根据经验确定随机变量的均值 μ 和频率最高时的数值(即密度函数的最大值点)m ,则 Beta 分布中的参数   可由以下关系求出:; {. w. \) Z5 ~4 C; C& {

4 \( q1 z7 k. n1 T. v$ c
" O. Y  _' x3 {* j& G! W( a! V
& M4 r& K: G# L  _' x# y 2 .2  计算机模拟
: C0 |4 u8 p- `5 z6 _: g. d4 t当排队系统的到达间隔时间和服务时间的概率分布很复杂时,或不能用公式给出 时,那么就不能用解析法求解。这就需用随机模拟法求解,现举例说明。1 F# k# G0 R* M% e

* h3 D+ W* Y, y: q5 I; C" G例 14 设某仓库前有一卸货场,货车一般是夜间到达,白天卸货,每天只能卸货 2 车,若一天内到达数超过 2 车,那么就推迟到次日卸货。根据表 3 所示的数据,货车到 达数的概率分布(相对频率)平均为 1.5 车/天,求每天推迟卸货的平均车数。9 `5 R+ h2 p! y5 ?+ W: Q6 z
$ x' P" w3 ?7 [6 D9 z

4 ~" ~/ Q& m3 m- _9 u3 u9 _; Q. m5 t" p2 S% s/ B' Z0 s
解 这是单服务台的排队系统,可验证到达车数不服从泊松分布,服务时间也不服 从指数分布(这是定长服务时间)。 随机模拟法首先要求事件能按历史的概率分布规律出现。模拟时产生的随机数与事 件的对应关系如表 4。
0 q/ g9 B* @% _( P9 B
. f# j) `: v! B) w1 T7 M" v0 p& {$ R% g* _$ U% d
! M5 N0 N7 J) N9 Z+ Q! L  ^5 ]
我们用 a1 表示产生的随机数,a2 表示到达的车数,a3 表示需要卸货车数,a4 表 示实际卸货车数,a5 表示推迟卸货车数。编写程序如下:
" l( H/ d9 a# \; {. N3 w# t" P) x
clear
: o( N/ M# U$ O! |6 L* `rand('state',sum(100*clock));# U5 B: ]: l/ \3 A7 O7 [, i
n=50000;0 D$ ~$ i6 w4 {+ s5 I, s8 o
m=2
5 `+ Z  ]8 p! K  {a1=rand(n,1);; h; s7 G: G$ q. z4 [7 M& z
a2=a1; %a2初始化! M/ G0 g; j+ a7 P
a2(find(a1<0.23))=0;
' K$ }8 Y9 \' W$ va2(find(0.23<=a1&a1<0.53))=1;8 O3 u; a' w$ n' }+ v
a2(find(0.53<=a1&a1<0.83))=2;
) q* k& q* T* Q1 @8 X( la2(find(0.83<=a1&a1<0.93),1)=3;0 m6 \/ ^! X; V2 M# Q8 j
a2(find(0.93<=a1&a1<0.98),1)=4;- }+ `3 w* G+ L" E
a2(find(a1>=0.98))=5;
9 o+ Z/ ?; C: u/ Ga3=zeros(n,1);a4=zeros(n,1);a5=zeros(n,1); %a2初始化+ X: ^/ I$ ]& F6 z* l
a3(1)=a2(1);- _0 ~$ @/ j9 k9 F  i. x9 V, w
if a3(1)<=m% _3 T' e( N) ]. W: k( G  i
    a4(1)=a3(1);a5(1)=0;1 `: t3 f, u- d
else! j2 b" O% |6 r
     a4(1)=m;a5(1)=a2(1)-m;) C: Z7 h: j( q
end1 T$ h% `8 \6 |) t: L* F, ?
for i=2:n
0 V# {0 O0 ]0 K" J    a3(i)=a2(i)+a5(i-1);) ~! z6 }1 F5 a8 S
    if a3(i)<=m
( f; o, N- O8 ~6 Z* p        a4(i)=a3(i);a5(i)=0;
! Y% n  V1 G% E$ w6 V" m' O    else
1 |# `- A. E/ ]/ w3 b( S; U        a4(i)=m;a5(i)=a3(i)-m;
% O- N+ A# G! h    end1 U7 T7 b  V: S
end
$ p* L4 w- H( d3 O! l* t$ Va=[a1,a2,a3,a4,a5];
% H6 q" {1 G% V, f- E9 Xsum(a)/n ' ~% j; P( \5 q3 h- I3 X
例 15 银行计划安置自动取款机,已知 A 型机的价格是 B 型机的 2 倍,而 A 型机 的性能—平均服务率也是 B 型机的 2 倍,问应该购置 1 台 A 型机还是 2 台 B 型机。 为了通过模拟回答这类问题,作如下具体假设,顾客平均每分钟到达 1 位, A 型 机的平均服务时间为 0.9 分钟, B 型机为 1.8 分钟,顾客到达间隔和服务时间都服从 指数分布,2 台 B 型机采取 M / M / 2 模型(排一队),用前 100 名顾客(第 1 位顾客到 达时取款机前为空)的平均等待时间为指标,对 A 型机和 B 型机分别作 1000 次模拟, 进行比较。
  R) g2 \$ ]2 t# C4 ?2 y) U+ A% Z' M( N/ Z" \4 {0 f0 r: ^8 }: {
, Z* H- |! I9 J8 q) |) ]
' P, {/ B6 m0 N
在模拟 A 型机时,我们用cspan表示到达间隔时间,sspan表示服务时间,ctime 表示到达时间,gtime表示离开时间,wtime表示等待时间。我们总共模拟了m 次, 每次n 个顾客。程序如下:5 O+ {% W! }2 m4 i& \9 u
7 x) S- A3 y; P1 H' V
tic0 Q5 b. p0 o8 f8 u
rand('state',sum(100*clock));
# v+ Z/ Z  w  o1 Z1 p  U) T1 Pn=100;m=1000;mu1=1;mu2=0.9;
, j0 [% O, H$ U( Wfor j=1:m  k9 l0 U4 V, q' [. i
    cspan=exprnd(mu1,1,n);sspan=exprnd(mu2,1,n);7 U; g" L9 d, |+ W& [! D9 s6 B
    ctime(1)=cspan(1);% k3 a- w. [4 H: a3 T8 q% e) j
    gtime(1)=ctime(1)+sspan(1);
+ z3 ?0 C, Y, Y* z    wtime(1)=0;: `# s7 G1 G2 c+ A$ n$ N1 P1 w7 O
    for i=2:n( ^, N) B, F/ B+ h/ x
        ctime(i)=ctime(i-1)+cspan(i);" b0 H, D, m( ^
        gtime(i)=max(ctime(i),gtime(i-1))+sspan(i);
5 Q$ }8 B' k% R* i' @  w0 W        wtime(i)=max(0,gtime(i-1)-ctime(i));" F7 \8 J0 n' [' w( f
    end3 B0 g! ~" d9 o) A4 I
    result1(j)=sum(wtime)/n;
3 j# n  g% e0 Uend' c0 k5 r# q9 e( \1 R# u1 D
result_1=sum(result1)/m + F, b+ W% ^& n! m/ W- l
toc6 W% L  H, Z/ a9 S. ]. c
类似地,模拟 B 型机的程序如下:
0 O! _# w4 F5 U* f
  z) H+ k1 Z( B1 ]) utic
  v4 @5 o: c6 v8 C/ M! Xrand('state',sum(100*clock));
& p( k9 G! b! D- p2 D3 ]$ o/ hn=100;m=1000;mu1=1;mu2=1.8;
/ s4 s0 m% Z/ G" W1 {for j=1:m
; X' V  Z7 O6 B8 m  K  i$ l    cspan=exprnd(mu1,1,n);sspan=exprnd(mu2,1,n);
! m5 l: M0 C% ?& d2 v5 P    ctime(1)=cspan(1);ctime(2)=ctime(1)+cspan(2);
0 h$ j, u" ~8 ?" z    gtime(1:2)=ctime(1:2)+sspan(1:2);: J5 M- e) u2 H
    wtime(1:2)=0;flag=gtime(1:2);7 v6 M. L8 N; v4 O3 |5 q5 K
    for i=3:n
3 V( ^6 r, n7 S: @; W0 `9 e        ctime(i)=ctime(i-1)+cspan(i);
. Q# r/ c. f1 P; L: @2 \        gtime(i)=max(ctime(i),min(flag))+sspan(i);
6 }) L* G4 [* W# C; l/ ^9 K        wtime(i)=max(0,min(flag)-ctime(i));: t# B! e/ a1 c: q4 @
        flag=[max(flag),gtime(i)];% S  k1 h9 @  R& R$ M* W" e$ d
    end- @+ O- p; @( Y: e- b4 S
    result2(j)=sum(wtime)/n;' ]/ v6 N4 X- S
end
- n3 j+ y( _; V+ ^( I. L7 Lresult_2=sum(result2)/m
2 K; _' r8 G) }. N' f; n6 Q0 ^& ptoc
; \7 t6 U  p* I读者可以用下面的程序与上面的程序比较了解编程的效率问题。8 ^& p: O: ?& f) h
# a: N% x0 V8 _+ d! I1 E' H2 L7 _
tic' m5 E6 _$ ]5 K3 m
clear! O8 ^& u* x4 ?- Y9 P, o
rand('state',sum(100*clock));
/ h* E9 k; T# S3 i' in=100;m=1000;mu1=1;mu2=0.9;, d  i6 L/ u3 m1 @
for j=1:m  S9 G, D# o' Z5 _/ k
    ctime(1)=exprnd(mu1);
9 }* W. R6 A" E# D% @    gtime(1)=ctime(1)+exprnd(mu2);* @3 x3 O' `- o$ c  R! u& R7 b# p
    wtime(1)=0;
, J% s0 N$ ?: R2 d8 w" p    for i=2:n+ r1 l2 d& u  y( k4 ^
        ctime(i)=ctime(i-1)+exprnd(mu1);
1 }2 P7 v! k4 `        gtime(i)=max(ctime(i),gtime(i-1))+exprnd(mu2);
; I: j7 P3 C: K0 j        wtime(i)=max(0,gtime(i-1)-ctime(i));" r$ L' ]) g' h5 Z
    end( d& e( u# R' P( o* j
    result(j)=sum(wtime)/n;% s" Q. S' s3 Q2 V+ y" i8 K1 f, u
end
% K- l1 C' V' kresult=sum(result)/m
/ w1 o9 o& H9 A* x/ Htoc
! f9 D& ^6 M1 e) D+ A  H1. 一个车间内有10台相同的机器,每台机器运行时每小时能创造4元的利润,且平 均每小时损坏一次。而一个修理工修复一台机器平均需4小时。以上时间均服从指数分 布。设一名修理工一小时工资为6元,试求:( V8 |; I4 A% L, o, I4 R! J! y
1 w  g% z% V+ I$ k# }% A
(i)该车间应设多少名修理工,使总费用为最小;8 _% C3 w8 \4 a; m: [( o& m5 s0 U

; X0 U/ u. [+ l8 d' y* z; B(ii)若要求不能运转的机器的期望数小于4台,则应设多少名修理工;1 J3 n/ B) L, Q" B6 ^% y5 S

: }4 o" Q2 m/ @' G(iii)若要求损坏机器等待修理的时间少于4小时,又应设多少名修理工。
5 H) V- t$ l! s7 O4 }5 x- F" Q# I( Q9 K3 v" j, g
2. 到达某铁路售票处顾客分两类:一类买南方线路票,到达率为λ1 /小时,另一 类买北方线路票,到达率为λ2 /小时,以上均服从泊松分布。该售票处设两个窗口,各窗口服务一名顾客时间均服从参数 μ = 10 的指数分布。试比较下列情况时顾客分别等 待时间Wq :' H' k5 E( d- @% s3 n8 x
, e8 ^: V/ t$ `, d& l4 j( S5 H
(i)两个窗口分别售南方票和北方票;: o2 b( }( |; T/ V  i7 n* T

( j' h6 Y* N7 q5 u9 i(ii)每个窗口两种票均出售。(分别比较 λ1 = λ2 = 2,4,6,8 时的情形)- D5 ~0 Y# H, ?, P
: m4 x( M  I$ h# z
3. 一名修理工负责5台机器的维修,每台机器平均每2h损坏一次,又修理工修复一 台机器平均需时18.75min,以上时间均服从负指数分布。试求:. z/ H7 Y2 S- _' `# \; \, Y
1 ?7 T3 g, i6 z. e1 ?# s
(1)所有机器均正常运转的概率;
4 i# l7 y" C: ~* K* K7 K1 r7 I
(2)等待维修的机器的期望数;4 R' R& j- f: ]1 C% Q
8 W) _- T- W* L
(3)假如希望做到有一半时间所有机器都正常运转,则该修理工最多看管多少台 机器。
) L1 Q! T& m; ^2 e) ], s( p1 U. ]0 I( x
(4)假如维修工工资为8元/h,机器不能正常运转时的损失为40元/h,则该修理工 看管多少台机器较为经济合理。
3 |% D3 p8 R: T" x  {, ]————————————————
- `6 E9 _, p9 r: W& T1 h版权声明:本文为CSDN博主「wamg潇潇」的原创文章,遵循CC 4.0 BY-SA版权协议,转载请附上原文出处链接及本声明。& q9 B4 [# O0 M  Z$ h
原文链接:https://blog.csdn.net/qq_29831163/java/article/details/89738145
6 E% t8 f. ^7 N1 l- @8 X2 i0 M7 S
, Z2 h+ \( h7 ^# U! y' ?( c- X1 V) Y' T- q1 v3 H2 h





欢迎光临 数学建模社区-数学中国 (http://www.madio.net/) Powered by Discuz! X2.5