数学建模社区-数学中国
标题:
排队论模型(八):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 E
8 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$ }& I
2.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) W
1 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% O
rand('state',sum(100*clock));
D; n h; R- y: h* N8 L! k( j
n=50000;
`# X2 {7 }/ R
m=2
% s( X! Y% }5 O) F3 L1 F
a1=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 j
a2(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 _( B
a2(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
else
3 M3 E1 T# o* q. a1 R0 T
a4(i)=m;a5(i)=a3(i)-m;
% M3 C2 z! X7 X+ I$ f; x
end
9 ]5 }0 l8 D5 D# v
end
: J9 t+ {% c, |+ F# E
a=[a1,a2,a3,a4,a5];
5 P. m( o- @5 y2 C0 F
sum(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) i
rand('state',sum(100*clock));
* f+ y% B. g5 x
n=100;m=1000;mu1=1;mu2=0.9;
% h6 J) p G" `. y2 x
for 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
end
7 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 E
toc
1 ]. e* m; K! q0 v# z* p
类似地,模拟 B 型机的程序如下:
/ l5 u; Z5 I6 p# m
2 c8 c* X! T/ t6 g* P: d
tic
3 q, p, i: o! G
rand('state',sum(100*clock));
+ [/ d- D" Q( M( y5 I! A
n=100;m=1000;mu1=1;mu2=1.8;
1 ^3 }; q8 |* c6 h
for 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 H
end
3 z# ?8 C$ }% j! o$ z" E, T
result_2=sum(result2)/m
0 y2 I4 ^0 D8 y0 n' p L0 i/ d
toc
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 ^% f
tic
, i7 Q: l* a4 h' J* L2 H, Z) J2 |
clear
1 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:n
4 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: B
result=sum(result)/m
2 {% r. B% d6 t1 c( q( ?- R
toc
9 m7 g3 P2 G1 ~
1. 一个车间内有10台相同的机器,每台机器运行时每小时能创造4元的利润,且平 均每小时损坏一次。而一个修理工修复一台机器平均需4小时。以上时间均服从指数分 布。设一名修理工一小时工资为6元,试求:
& K4 ~! H! Z; s
5 @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 P
9 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& H
3. 一名修理工负责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& a
9 K& I& K0 \' f! B1 S8 u
(3)假如希望做到有一半时间所有机器都正常运转,则该修理工最多看管多少台 机器。
( Q E% o+ d a& d9 t
4 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/89738145
8 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