数学建模社区-数学中国

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

作者: 浅夏110    时间: 2020-6-13 09:35
标题: 排队论模型(八):Matlab 生成随机数、排队模型的计算机模拟
1 产生给定分布的随机数的方法% g4 `- d- P& j5 t( F4 S" x- W/ X
Matlab 可以产生常用分布的随机数。下面我们介绍按照给定的概率分布产生随机数的一般方法,这些方法都以U(0,1) 分布的随机变量为基础。
3 d: E- w) P2 S0 V0 {. N. X: {7 T* Y& }! A, F2 k: v9 ]1 n
(i)反变换法, p1 H+ O) S% y# Q. r4 o- s) t$ E
定理 设 X 是一个具有连续分布函数 F(x) 的随机变量,则 F(X ) 在 [0,1] 上服 从均匀分布。+ v$ B  r, K" C

+ {- ^* d6 ^0 l# u) j7 W( i, ?) T) t

1 S  r5 H. C5 E8 n* h4 P! G$ n& h# i& y
1 ~1 Q) d; x# [. v$ k9 i
(ii)卷积法# b2 p& K* s! }
/ r  o4 J. p( i6 k$ P; r! \
. o# V. i  n5 N) I0 `: l; f
$ [+ ^: Q& R# |; m" U( t+ k) O
(iii)取舍法
2 K# l1 k$ ?8 o! a4 L若随机变量 X 在有限区间(a,b) 内变化,但概率密度 f (x)具有任意形式(甚至没 有解析表达式),无法用前面的方法产生时,可用取舍法。一种比较简单的取舍法的步 骤是:
( d8 W5 H8 X4 T1 l
: ^  x( F7 A3 M+ n$ l
7 K' U4 y5 \' }/ G$ _! }/ G7 B
8 U9 [* c4 k7 r: Z3 }2 排队模型的计算机模拟
; b6 E. z* T$ }& I2.1 确定随机变量概率分布的常用方法" m, B0 l+ b! a1 f2 U) |
在模拟一个带有随机因素的实际系统时,究竟用什么样的概率分布描述问题中的随 机变量,是我们总是要碰到的一个问题,下面简单介绍确定分布的常用方法:% e; i/ n% z1 c2 l% C1 y+ o9 P
" F* s3 _" x; {# w# Y
【1 】根据一般知识和经验,可以假定其概率分布的形式,如顾客到达间隔服从指数 分布 Exp(λ) ;产品需求量服从正态分布   ;订票后但未能按时前往机场登机 的人数服从二项分布 B(n, p) 。然后由实际数据估计分布的参数 λ,μ,σ 等,参数估计 可用极大似然估计、矩估计等方法。2 q9 a' p$ Y( ]) m' g! x/ s
0 G) \; u7 q2 T9 ~3 E- v
【2】 直接由大量的实际数据作直方图,得到经验分布,再通过假设检验,拟合分布 函数,可用  检验等方法。 3 o 既缺少先验知识,又缺少数据时,对区间(a,b) 内变化的随机变量,可选用 Beta 分布(包括均匀分布)。先根据经验确定随机变量的均值 μ 和频率最高时的数值(即密度函数的最大值点)m ,则 Beta 分布中的参数   可由以下关系求出:
+ _+ t" U. c. T5 K) W1 f7 n1 \8 m* d* `& o

! ~2 ]' Q. @3 k6 N/ {! V, _% `! }/ R7 b: c
2 .2  计算机模拟
$ d1 j. {# k- q% g1 b% s当排队系统的到达间隔时间和服务时间的概率分布很复杂时,或不能用公式给出 时,那么就不能用解析法求解。这就需用随机模拟法求解,现举例说明。7 K$ x+ {* C: E0 B! d

% a6 a' X% Y8 E$ \) c例 14 设某仓库前有一卸货场,货车一般是夜间到达,白天卸货,每天只能卸货 2 车,若一天内到达数超过 2 车,那么就推迟到次日卸货。根据表 3 所示的数据,货车到 达数的概率分布(相对频率)平均为 1.5 车/天,求每天推迟卸货的平均车数。
1 K- g9 V1 r( G: Y$ w2 J) g$ y
1 }( ~/ y9 s) n# N
; u: p8 Q' @+ B3 ~' ~) E  z' n0 N9 f+ w( `8 I* c5 z
解 这是单服务台的排队系统,可验证到达车数不服从泊松分布,服务时间也不服 从指数分布(这是定长服务时间)。 随机模拟法首先要求事件能按历史的概率分布规律出现。模拟时产生的随机数与事 件的对应关系如表 4。
# Q- U" R7 V  @* [: ?% ^
( C8 U" [' h9 W) `4 |
; F& c3 |5 B3 H+ F1 j$ R6 C
7 Z. E- v1 G: c! ]6 I5 u我们用 a1 表示产生的随机数,a2 表示到达的车数,a3 表示需要卸货车数,a4 表 示实际卸货车数,a5 表示推迟卸货车数。编写程序如下:
3 J/ }2 J8 K& A- J* _: y, r( J; |! |8 b8 R  e1 q- ^
clear
- `' x) j0 @* e3 u4 H3 M% Orand('state',sum(100*clock));  D; n  h; R- y: h* N8 L! k( j
n=50000;
  `# X2 {7 }/ Rm=2
% s( X! Y% }5 O) F3 L1 Fa1=rand(n,1);/ I7 v) Q; c6 U% a- K
a2=a1; %a2初始化! l) D( U! w2 \  H! w
a2(find(a1<0.23))=0;
6 [! M( v5 ?" w; X& H- l# _+ ~6 ja2(find(0.23<=a1&a1<0.53))=1;
' R* W$ T: D! N7 {1 ~a2(find(0.53<=a1&a1<0.83))=2;7 t0 W8 L0 e  b
a2(find(0.83<=a1&a1<0.93),1)=3;% W& b! K5 ^( Q' r
a2(find(0.93<=a1&a1<0.98),1)=4;
1 o7 q5 }, L  _( Ba2(find(a1>=0.98))=5;" w5 D" ]. J/ `: V* J  Z
a3=zeros(n,1);a4=zeros(n,1);a5=zeros(n,1); %a2初始化. U/ Q% g3 X! z
a3(1)=a2(1);. S# b7 V. u! B
if a3(1)<=m
8 @# d. E4 M, v6 Q" h! }    a4(1)=a3(1);a5(1)=0;! W3 Y) r; _* A8 E5 e
else
. G. t6 |/ t, m. U/ w; k7 _* J     a4(1)=m;a5(1)=a2(1)-m;
) x8 h$ D( p# [end. b/ }+ }: p+ \! \. O( W6 t* a: g: B, i" ?
for i=2:n
! {( Q2 {. ^! `* n7 G2 [    a3(i)=a2(i)+a5(i-1);8 _3 _. B8 x* X( o' h2 q
    if a3(i)<=m
3 v/ P. S9 Y" z/ B        a4(i)=a3(i);a5(i)=0;
/ Q5 L+ r" @: i) V5 M; v    else3 M3 E1 T# o* q. a1 R0 T
        a4(i)=m;a5(i)=a3(i)-m;% M3 C2 z! X7 X+ I$ f; x
    end9 ]5 }0 l8 D5 D# v
end
: J9 t+ {% c, |+ F# Ea=[a1,a2,a3,a4,a5];
5 P. m( o- @5 y2 C0 Fsum(a)/n
3 u' `  m9 b" W* p# Y例 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 次模拟, 进行比较。
4 T* X6 H2 B- r; N6 q$ O( c/ w+ X& _
) T$ ]  R2 a6 x' t8 `5 z: r/ @+ R/ M$ f# F/ ~# d

; v, p! h! O2 Z在模拟 A 型机时,我们用cspan表示到达间隔时间,sspan表示服务时间,ctime 表示到达时间,gtime表示离开时间,wtime表示等待时间。我们总共模拟了m 次, 每次n 个顾客。程序如下:9 G* J  V9 M8 ?; q, @
& L* o6 K; _2 y
tic
% S. S- P0 V% _7 Q1 R$ h) irand('state',sum(100*clock));* f+ y% B. g5 x
n=100;m=1000;mu1=1;mu2=0.9;
% h6 J) p  G" `. y2 xfor j=1:m
1 y% b1 z- J6 S    cspan=exprnd(mu1,1,n);sspan=exprnd(mu2,1,n);
" X. F( @' ^# t0 B) M    ctime(1)=cspan(1);
6 l4 e  B! x* k3 q6 y! j$ b4 l* z( P    gtime(1)=ctime(1)+sspan(1);
! r# ~3 q3 q& r- H; s5 f0 l    wtime(1)=0;* G  W1 n$ Y$ D3 i# n3 J" b* I
    for i=2:n
7 Q% x7 \3 w: k$ \        ctime(i)=ctime(i-1)+cspan(i);1 k. L/ `5 W; [( `
        gtime(i)=max(ctime(i),gtime(i-1))+sspan(i);
( N; R6 q+ v6 K: D- g) O$ ?; n        wtime(i)=max(0,gtime(i-1)-ctime(i));
! v; B# R0 z  V    end7 J- }7 q0 Z6 {" H' {$ E: {
    result1(j)=sum(wtime)/n;7 O$ o: U% f9 B# q) `5 g
end% a. E* K7 d' w
result_1=sum(result1)/m
, I  J5 U5 T9 q$ U2 Etoc
1 ]. e* m; K! q0 v# z* p类似地,模拟 B 型机的程序如下:/ l5 u; Z5 I6 p# m
2 c8 c* X! T/ t6 g* P: d
tic3 q, p, i: o! G
rand('state',sum(100*clock));
+ [/ d- D" Q( M( y5 I! An=100;m=1000;mu1=1;mu2=1.8;
1 ^3 }; q8 |* c6 hfor j=1:m
, M" s+ {, S, Q+ X) P. ?    cspan=exprnd(mu1,1,n);sspan=exprnd(mu2,1,n);
* g+ I0 m0 X: r/ m- \    ctime(1)=cspan(1);ctime(2)=ctime(1)+cspan(2);' t6 ]# b4 I% ^3 D2 D; k
    gtime(1:2)=ctime(1:2)+sspan(1:2);  W) W% x$ d: l6 i
    wtime(1:2)=0;flag=gtime(1:2);, v$ a6 ~+ D+ ]7 E, R. u/ m
    for i=3:n
: n$ g7 G- l& g: n* Y  o8 D        ctime(i)=ctime(i-1)+cspan(i);
* j7 i9 G$ |7 e8 l! X        gtime(i)=max(ctime(i),min(flag))+sspan(i);
8 {' s8 Y: T0 w! e        wtime(i)=max(0,min(flag)-ctime(i));, a7 R- x2 L7 f' I: h+ C  o8 n
        flag=[max(flag),gtime(i)];2 `9 J) T& n$ `. R, x3 \1 ^
    end
) x9 {: h& \8 |! {5 y    result2(j)=sum(wtime)/n;
- c6 v* Z; m0 o1 Hend3 z# ?8 C$ }% j! o$ z" E, T
result_2=sum(result2)/m
0 y2 I4 ^0 D8 y0 n' p  L0 i/ dtoc
0 @0 Z9 d8 h7 U# x1 \9 L/ D  B读者可以用下面的程序与上面的程序比较了解编程的效率问题。0 W, r" X+ C  N4 R' }

( a5 p3 e5 E+ c0 M' I! X9 ^% ftic, i7 Q: l* a4 h' J* L2 H, Z) J2 |
clear1 i; S  I" w$ f1 Y: }0 }0 X
rand('state',sum(100*clock));# v9 ^; n  k% {2 t" N  H. p: q7 H
n=100;m=1000;mu1=1;mu2=0.9;
3 T) s6 S+ Q" V3 H7 ?for j=1:m/ j% a6 B. e% R8 e
    ctime(1)=exprnd(mu1);& R0 S2 N7 O8 v) U0 S  x( q
    gtime(1)=ctime(1)+exprnd(mu2);
& @( X, s  o. ]    wtime(1)=0;
* _5 D' @4 c# w( O6 M    for i=2:n4 X* }4 X2 e! ]
        ctime(i)=ctime(i-1)+exprnd(mu1);) u& S# p* o; {2 C9 I8 Z  t7 @3 M
        gtime(i)=max(ctime(i),gtime(i-1))+exprnd(mu2);
, V+ T: A5 Q  O% Q- z9 p+ F        wtime(i)=max(0,gtime(i-1)-ctime(i));& L/ O( ^  _% g1 C1 u
    end( o9 s% @4 N0 D( V) I1 \) ?
    result(j)=sum(wtime)/n;/ j- s1 Y0 I1 n% O  H7 j% Q) y
end
/ c" |! |' U+ d4 m: w: Bresult=sum(result)/m2 {% r. B% d6 t1 c( q( ?- R
toc
9 m7 g3 P2 G1 ~1. 一个车间内有10台相同的机器,每台机器运行时每小时能创造4元的利润,且平 均每小时损坏一次。而一个修理工修复一台机器平均需4小时。以上时间均服从指数分 布。设一名修理工一小时工资为6元,试求:
& K4 ~! H! Z; s5 @5 E7 r% A1 u8 |6 c! Y
(i)该车间应设多少名修理工,使总费用为最小;
2 y' C9 V/ D4 u; _% M$ d4 h* X1 ]9 i
(ii)若要求不能运转的机器的期望数小于4台,则应设多少名修理工;
, d1 F8 h7 O; x/ ~( ]+ y% T! G% w# I0 v* F- g% h1 |
(iii)若要求损坏机器等待修理的时间少于4小时,又应设多少名修理工。
' h. _4 T+ M5 N/ E
5 r& c0 U8 A, L: {2. 到达某铁路售票处顾客分两类:一类买南方线路票,到达率为λ1 /小时,另一 类买北方线路票,到达率为λ2 /小时,以上均服从泊松分布。该售票处设两个窗口,各窗口服务一名顾客时间均服从参数 μ = 10 的指数分布。试比较下列情况时顾客分别等 待时间Wq :
( O  |- {- L2 D+ `9 h& z+ F( Q1 ~. G7 t" z; i
(i)两个窗口分别售南方票和北方票;
0 n( l/ L$ T8 P9 I* J9 l4 n9 t  T# J
(ii)每个窗口两种票均出售。(分别比较 λ1 = λ2 = 2,4,6,8 时的情形)/ P; Z( ]$ R$ W3 @( `9 [1 r

( i. {9 g9 N, V& H3. 一名修理工负责5台机器的维修,每台机器平均每2h损坏一次,又修理工修复一 台机器平均需时18.75min,以上时间均服从负指数分布。试求:' j* D  o9 c5 }7 P7 H. I. s6 @' w# D: ?
5 w3 m+ ?# h" I5 d; ]5 |6 {, f- C% V1 a, `
(1)所有机器均正常运转的概率;
  V: P8 v9 z$ m
6 f9 c: q3 p1 A; }' m(2)等待维修的机器的期望数;
+ I) A% M, P! a$ l5 q4 u& a9 K& I& K0 \' f! B1 S8 u
(3)假如希望做到有一半时间所有机器都正常运转,则该修理工最多看管多少台 机器。
( Q  E% o+ d  a& d9 t4 X, r# b. ?: N% D# F7 E
(4)假如维修工工资为8元/h,机器不能正常运转时的损失为40元/h,则该修理工 看管多少台机器较为经济合理。" {1 t7 t) t1 g- W) t' [
————————————————8 \6 d3 N2 x/ [
版权声明:本文为CSDN博主「wamg潇潇」的原创文章,遵循CC 4.0 BY-SA版权协议,转载请附上原文出处链接及本声明。7 i- v5 y" R5 j) |
原文链接:https://blog.csdn.net/qq_29831163/java/article/details/897381458 h1 ]6 [' b8 t+ L& I7 `6 Y. \
* I; c! h6 ]2 _
( c7 q; z" s* _( {





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