数学建模社区-数学中国

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

作者: 浅夏110    时间: 2020-6-13 09:35
标题: 排队论模型(八):Matlab 生成随机数、排队模型的计算机模拟
1 产生给定分布的随机数的方法, m- O* x% m9 \) Q7 g2 Z8 ?8 _
Matlab 可以产生常用分布的随机数。下面我们介绍按照给定的概率分布产生随机数的一般方法,这些方法都以U(0,1) 分布的随机变量为基础。
% a. q* M% q5 h: a# }+ A  C9 ~# f
4 P4 n8 ^! K# i" N(i)反变换法
- d+ S3 e, Z0 i7 k' U: W. v定理 设 X 是一个具有连续分布函数 F(x) 的随机变量,则 F(X ) 在 [0,1] 上服 从均匀分布。& X$ n9 f( c9 D. h8 ^7 T# b

/ k* o6 R! ^+ ~8 F+ W% C, @
4 e/ D. i0 U* a% Q9 K2 M8 K
$ ]5 C2 Z( A, W- @
0 {) R4 H$ C) x% G- B" M7 @( c" V2 S8 }1 D
(ii)卷积法
% U2 b/ U6 n" m6 d& a9 ^; B# E* Q) n$ V9 R) C$ t
: j0 Y# D3 Q$ y
4 M* E. K: s0 d0 E) I& T6 q
(iii)取舍法
) s4 I' Q! c3 U- N若随机变量 X 在有限区间(a,b) 内变化,但概率密度 f (x)具有任意形式(甚至没 有解析表达式),无法用前面的方法产生时,可用取舍法。一种比较简单的取舍法的步 骤是:( P* r9 E! J) Q; [! D4 z8 c7 D  q

; X! s+ J" y% b8 |1 x4 ?, U9 w2 d' @! k# s) ?& W5 s
* |* T; R# s1 D, P. n5 ~# u
2 排队模型的计算机模拟
! D5 \# r) }0 u  U' T2.1 确定随机变量概率分布的常用方法
0 q4 U/ j- s4 J2 T2 k在模拟一个带有随机因素的实际系统时,究竟用什么样的概率分布描述问题中的随 机变量,是我们总是要碰到的一个问题,下面简单介绍确定分布的常用方法:5 u0 ^: Y: u& z" x1 W" A

' a% T. P7 p" P5 X4 c: R: ^【1 】根据一般知识和经验,可以假定其概率分布的形式,如顾客到达间隔服从指数 分布 Exp(λ) ;产品需求量服从正态分布   ;订票后但未能按时前往机场登机 的人数服从二项分布 B(n, p) 。然后由实际数据估计分布的参数 λ,μ,σ 等,参数估计 可用极大似然估计、矩估计等方法。
% t: p: \. u2 k2 K& m$ g/ q5 o" F( ~- [- @/ G( o" k  d
【2】 直接由大量的实际数据作直方图,得到经验分布,再通过假设检验,拟合分布 函数,可用  检验等方法。 3 o 既缺少先验知识,又缺少数据时,对区间(a,b) 内变化的随机变量,可选用 Beta 分布(包括均匀分布)。先根据经验确定随机变量的均值 μ 和频率最高时的数值(即密度函数的最大值点)m ,则 Beta 分布中的参数   可由以下关系求出:
# W" o; }4 @' y! y9 E# F7 d# {9 D$ ]! G0 B+ S- c; j6 x6 v

6 p) {2 j  Y4 f1 {
, b/ \* {( l5 X2 A 2 .2  计算机模拟
, Y; m' q4 K' U9 R, k6 n8 k0 L当排队系统的到达间隔时间和服务时间的概率分布很复杂时,或不能用公式给出 时,那么就不能用解析法求解。这就需用随机模拟法求解,现举例说明。# z2 p; W; c& `( T2 `$ A* ?

! j/ Z1 d! y- P* \* O( i- z& p例 14 设某仓库前有一卸货场,货车一般是夜间到达,白天卸货,每天只能卸货 2 车,若一天内到达数超过 2 车,那么就推迟到次日卸货。根据表 3 所示的数据,货车到 达数的概率分布(相对频率)平均为 1.5 车/天,求每天推迟卸货的平均车数。$ T$ f5 l, w3 s# u* m# ]  C

& [  w( N9 N+ n/ C; q$ V2 u
+ V% O( i! V3 |2 g
. \# `2 I2 T" C# j解 这是单服务台的排队系统,可验证到达车数不服从泊松分布,服务时间也不服 从指数分布(这是定长服务时间)。 随机模拟法首先要求事件能按历史的概率分布规律出现。模拟时产生的随机数与事 件的对应关系如表 4。
4 F8 G. x" d1 V# x9 l7 d! {
; W9 r' W0 ]$ d0 @
& @) _& ?2 d3 N% q& W* J# ?& ]/ ]  s( G' T! k4 i  l+ r
我们用 a1 表示产生的随机数,a2 表示到达的车数,a3 表示需要卸货车数,a4 表 示实际卸货车数,a5 表示推迟卸货车数。编写程序如下:
3 x% c, }8 d( p+ w; P$ M! H! N" f2 i% ^) T+ X) f
clear
& v! v* S* ]; I0 Grand('state',sum(100*clock));* @% S  u- e& j
n=50000;. d7 o: C* Y) R! S5 t
m=2
" k$ C- K* H* i1 ~a1=rand(n,1);; E& \" Q) H8 a0 w* X5 K
a2=a1; %a2初始化
9 X6 b5 e$ u. _$ {a2(find(a1<0.23))=0;
% l6 t# b; c; J( h* l/ d( ga2(find(0.23<=a1&a1<0.53))=1;
) M0 d8 f+ n$ l* J# X) c' g+ M% b6 o$ pa2(find(0.53<=a1&a1<0.83))=2;' _8 P7 y! @7 @# ?# B( p) h8 @
a2(find(0.83<=a1&a1<0.93),1)=3;
. a% _# g+ S( @/ Ea2(find(0.93<=a1&a1<0.98),1)=4;; D. u' k8 \4 b. o+ t
a2(find(a1>=0.98))=5;
1 ]/ m4 ^7 |4 E4 |3 T# H: t- t; Ea3=zeros(n,1);a4=zeros(n,1);a5=zeros(n,1); %a2初始化, e; V; B9 s" q4 b& ~* G6 Y
a3(1)=a2(1);( F$ I' g- ?8 u* Q3 R: n
if a3(1)<=m  y! p  {- c" |) P
    a4(1)=a3(1);a5(1)=0;6 i/ k  o4 Z; p: t7 W+ [
else8 p# q. g# k& q7 R% n
     a4(1)=m;a5(1)=a2(1)-m;
( g) `$ z4 {+ H' Fend
' d, n& g6 o- S- u: c  @for i=2:n  y. m1 E( X' X, E
    a3(i)=a2(i)+a5(i-1);
! Z; t. p% ^5 o- x% w# X    if a3(i)<=m
9 v7 e; y; A$ q, F& A' E2 a* _        a4(i)=a3(i);a5(i)=0;
- y2 m9 C% _( A    else! W5 ?5 C' m( H; i4 F
        a4(i)=m;a5(i)=a3(i)-m;: y; l( v* P. i% T* s
    end
. S" D# a: q& aend
1 W9 C8 w( t. B* B# E3 y( [. C# va=[a1,a2,a3,a4,a5];/ k% q6 g1 L8 @/ Z! i8 F/ f! {7 L4 ]
sum(a)/n * Z7 U1 t$ j6 M6 b6 s( ?
例 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 次模拟, 进行比较。3 ]9 `1 W: `4 P: N& R8 I
4 R- n$ c2 d5 |1 v! l/ f
( r4 v& ?' \% }6 ~% q# O4 d
4 b1 |4 ~) Z. H+ ^; N6 y* ^# r
在模拟 A 型机时,我们用cspan表示到达间隔时间,sspan表示服务时间,ctime 表示到达时间,gtime表示离开时间,wtime表示等待时间。我们总共模拟了m 次, 每次n 个顾客。程序如下:3 k/ c* @' {1 j7 g8 u
  W; ^- Z7 z' D- D7 \* j  _) e. [
tic5 S+ [5 Q  t, O% D* y, D1 ?. S+ |& D
rand('state',sum(100*clock));
8 E& W. v) C) qn=100;m=1000;mu1=1;mu2=0.9;
9 y1 s! s8 z1 j$ y1 Qfor j=1:m: k6 A* G1 J0 T2 Z, l, O6 A
    cspan=exprnd(mu1,1,n);sspan=exprnd(mu2,1,n);
" g/ C5 D* E+ e  a6 ]$ U$ d/ I    ctime(1)=cspan(1);
) `) G9 q& v  Y/ q+ d$ S) n# t    gtime(1)=ctime(1)+sspan(1);
2 J( _1 E0 M  y" F: H9 B8 }    wtime(1)=0;
3 K$ b+ K2 y( i% r6 \/ ~    for i=2:n
8 ?* e7 \6 g, D3 X' {& }6 ?% I2 c        ctime(i)=ctime(i-1)+cspan(i);: j1 t, n7 K0 ~, X, k) E
        gtime(i)=max(ctime(i),gtime(i-1))+sspan(i);
. Y2 O) g% j2 r" [- L% o        wtime(i)=max(0,gtime(i-1)-ctime(i));* \- i* U# n+ Q' K( ^% q) @
    end) @% d6 f& O+ D! q3 f
    result1(j)=sum(wtime)/n;
7 A  [4 S* Y( _# j$ \end9 b9 R+ v  e/ s. c5 d
result_1=sum(result1)/m
; K2 x3 T6 _  ~2 F; v  z* I; q, ~5 rtoc
$ a6 w9 J3 N1 I, X2 z# f类似地,模拟 B 型机的程序如下:5 S6 |0 S% z+ B& V. Z# U# G

' S" G6 g" ?' R' k1 utic& }" D. Q7 V# M& d; Y$ u; B+ o  S
rand('state',sum(100*clock));" G! N) V) W% b2 h3 i) t
n=100;m=1000;mu1=1;mu2=1.8;
' F" _5 ~4 C1 q* L9 c/ _& B; qfor j=1:m7 W9 n( _" S- P; J& i5 v' L0 E
    cspan=exprnd(mu1,1,n);sspan=exprnd(mu2,1,n);
; a& t9 Q/ g( Q3 }7 M    ctime(1)=cspan(1);ctime(2)=ctime(1)+cspan(2);
- o( H( Z. a0 Y$ B, U    gtime(1:2)=ctime(1:2)+sspan(1:2);
% I# P7 Z" G) e( R9 J/ {1 q1 ~    wtime(1:2)=0;flag=gtime(1:2);4 f. c- @& i  N2 T+ O; B
    for i=3:n2 X8 b' T8 K# I  }, n4 v9 z
        ctime(i)=ctime(i-1)+cspan(i);' p# g9 a! x- B9 o" T. k3 s
        gtime(i)=max(ctime(i),min(flag))+sspan(i);$ h# B! |- o$ }8 R1 c. m
        wtime(i)=max(0,min(flag)-ctime(i));- f. y! v& m! @0 d2 Z3 J+ [
        flag=[max(flag),gtime(i)];
2 Y) ~" b. m* f$ ^    end% R: z& P9 @9 O. V- P+ ]+ W+ T7 h
    result2(j)=sum(wtime)/n;( T5 S0 O+ d! I; ^! t
end. z9 q: b9 r/ S& D
result_2=sum(result2)/m; w$ l/ D; Z; x. H. q1 t! ]
toc
/ k. K: x; b% V: \% T" |6 @读者可以用下面的程序与上面的程序比较了解编程的效率问题。2 y6 p/ U6 @/ P' k. x6 u0 Q

3 n* R+ y7 y, otic
* B7 t2 j: ?( U0 b5 D& G& t( f/ _8 C6 Gclear
% O  J" {4 W) c1 o% }rand('state',sum(100*clock));, M- S8 f( U5 T
n=100;m=1000;mu1=1;mu2=0.9;/ ~3 N# y% y1 D% A/ _) x- q- N2 u7 @
for j=1:m# Y/ R. l, h$ k8 n- G6 y* E
    ctime(1)=exprnd(mu1);
( \8 s& e. F! Y0 R0 g    gtime(1)=ctime(1)+exprnd(mu2);
. U- I! Y7 i2 V1 W8 _    wtime(1)=0;
8 V6 z& y1 I6 `6 m# _6 D    for i=2:n
7 n1 z- s+ R1 `        ctime(i)=ctime(i-1)+exprnd(mu1);% O- J& M* ]& C
        gtime(i)=max(ctime(i),gtime(i-1))+exprnd(mu2);" c/ b3 [) S0 m  W6 L& n
        wtime(i)=max(0,gtime(i-1)-ctime(i));0 X: C8 }9 H0 U" {  Y
    end
; p* d" c7 q! h; ^% g, L3 w    result(j)=sum(wtime)/n;
  P1 _: B5 O  |+ Xend& r$ L; Y5 o. T  ?
result=sum(result)/m
* ~8 f6 s3 Y2 c1 X9 S% Etoc' f) [, W) b% Q( h0 Z# }
1. 一个车间内有10台相同的机器,每台机器运行时每小时能创造4元的利润,且平 均每小时损坏一次。而一个修理工修复一台机器平均需4小时。以上时间均服从指数分 布。设一名修理工一小时工资为6元,试求:
0 M9 k0 c7 J- N6 a* @+ w: D, s1 y! R+ i3 B" g
(i)该车间应设多少名修理工,使总费用为最小;( O! P% O2 i2 n# q) l- c, Y! q! D
; p5 l( _# ]9 v8 n7 H* O
(ii)若要求不能运转的机器的期望数小于4台,则应设多少名修理工;
, F( ^8 B" W/ S7 A! r5 P# Y$ S; a0 F/ l5 v
(iii)若要求损坏机器等待修理的时间少于4小时,又应设多少名修理工。
% @' r! H# Y3 J/ c" j! |1 }7 l$ R( w7 @( B0 T7 v
2. 到达某铁路售票处顾客分两类:一类买南方线路票,到达率为λ1 /小时,另一 类买北方线路票,到达率为λ2 /小时,以上均服从泊松分布。该售票处设两个窗口,各窗口服务一名顾客时间均服从参数 μ = 10 的指数分布。试比较下列情况时顾客分别等 待时间Wq :7 }' A3 p& f- m
6 P( Y, H7 [: E# F- f
(i)两个窗口分别售南方票和北方票;
7 I+ p' K" Z7 ^1 T! A/ D$ \* B: |6 x" A  J- w- o4 t' `& |; o
(ii)每个窗口两种票均出售。(分别比较 λ1 = λ2 = 2,4,6,8 时的情形)
4 F! x% A" V# `! Z9 Q( g6 C
8 L" s9 w. T/ w1 D3. 一名修理工负责5台机器的维修,每台机器平均每2h损坏一次,又修理工修复一 台机器平均需时18.75min,以上时间均服从负指数分布。试求:
5 C1 }+ _6 M" X% b) T; ]$ m) o0 p; v! E8 z; G/ o2 q
(1)所有机器均正常运转的概率;9 K! Y8 d. m& V( A. o8 ^5 ^
  s6 z/ h7 \# b/ A/ w5 ~/ ], B
(2)等待维修的机器的期望数;$ k" M! L, i# n

7 c5 T5 L8 a0 h1 n& q4 i(3)假如希望做到有一半时间所有机器都正常运转,则该修理工最多看管多少台 机器。
, ?% w/ t# K0 k) S' v! R% i2 h. [- Q! ?7 _: o$ u9 a" v
(4)假如维修工工资为8元/h,机器不能正常运转时的损失为40元/h,则该修理工 看管多少台机器较为经济合理。) r2 J1 ?2 }! _2 q; Q" F% R
————————————————) R8 [6 O3 v# {! ?7 ?$ E
版权声明:本文为CSDN博主「wamg潇潇」的原创文章,遵循CC 4.0 BY-SA版权协议,转载请附上原文出处链接及本声明。
2 z5 o& q8 l# k, i" e原文链接:https://blog.csdn.net/qq_29831163/java/article/details/89738145
$ K2 d$ U1 P8 c* ]9 _. y: J6 r% m7 p7 `% N

1 Y- ?, G1 U  A4 t




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