数学建模社区-数学中国

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

作者: 浅夏110    时间: 2020-6-13 09:35
标题: 排队论模型(八):Matlab 生成随机数、排队模型的计算机模拟
1 产生给定分布的随机数的方法8 q9 Q9 |# ~9 H$ ?
Matlab 可以产生常用分布的随机数。下面我们介绍按照给定的概率分布产生随机数的一般方法,这些方法都以U(0,1) 分布的随机变量为基础。3 B# y0 ~/ Q8 H# h+ W
0 A" c" b; Q. r+ o1 H; G, P
(i)反变换法
7 Y0 k, Q/ D1 }' d定理 设 X 是一个具有连续分布函数 F(x) 的随机变量,则 F(X ) 在 [0,1] 上服 从均匀分布。0 X' O0 S/ |: J: X6 H2 J

$ E9 g0 o4 T3 l+ b% p" Q* J( \9 U
( V4 C5 w$ G  n5 Z$ W9 K

; P3 l* Z) ^0 f! L' L, _7 Z) V7 w3 [& _7 j
(ii)卷积法
: C+ h3 n! Z% A! L" d1 D
" n& c" J0 U9 Y3 e- [" x: {  t* ]  \
! [! B( t% V2 O+ n/ W
' l5 x$ S' J; a7 `) t8 {(iii)取舍法
; s+ M2 t# ~- f0 h' d若随机变量 X 在有限区间(a,b) 内变化,但概率密度 f (x)具有任意形式(甚至没 有解析表达式),无法用前面的方法产生时,可用取舍法。一种比较简单的取舍法的步 骤是:
* S9 V" ]0 ]  h$ E5 ]4 a2 }4 e' [' H. u5 `1 M  ~5 u+ Q. H
) n( o, u2 Y- f; k4 f* j' P

9 j  ~% Y8 S6 w; K" i) V% _2 排队模型的计算机模拟
- x# [5 }  ]1 w7 H, S( P2.1 确定随机变量概率分布的常用方法
9 Z5 P: R! `* j! L" R在模拟一个带有随机因素的实际系统时,究竟用什么样的概率分布描述问题中的随 机变量,是我们总是要碰到的一个问题,下面简单介绍确定分布的常用方法:7 K. Q4 z' s  V3 D% G

, v( z' G/ _3 ?! Q" F% ?【1 】根据一般知识和经验,可以假定其概率分布的形式,如顾客到达间隔服从指数 分布 Exp(λ) ;产品需求量服从正态分布   ;订票后但未能按时前往机场登机 的人数服从二项分布 B(n, p) 。然后由实际数据估计分布的参数 λ,μ,σ 等,参数估计 可用极大似然估计、矩估计等方法。6 o" A: B' Y; U0 v1 ~
% r4 `3 N/ V( ?% G7 u4 P9 s0 ?
【2】 直接由大量的实际数据作直方图,得到经验分布,再通过假设检验,拟合分布 函数,可用  检验等方法。 3 o 既缺少先验知识,又缺少数据时,对区间(a,b) 内变化的随机变量,可选用 Beta 分布(包括均匀分布)。先根据经验确定随机变量的均值 μ 和频率最高时的数值(即密度函数的最大值点)m ,则 Beta 分布中的参数   可由以下关系求出:
/ r, x2 D$ V. i' L3 y# J! @" k5 b8 y' K( H
& j! r0 P  }- F! x
& ?# C! @( O8 C8 m4 f4 I
2 .2  计算机模拟
/ E; q; k& r1 d0 o) H当排队系统的到达间隔时间和服务时间的概率分布很复杂时,或不能用公式给出 时,那么就不能用解析法求解。这就需用随机模拟法求解,现举例说明。
) w2 }5 S. U5 d/ m: e. T2 t9 p7 ]0 h* q- }" L% L4 u* D
例 14 设某仓库前有一卸货场,货车一般是夜间到达,白天卸货,每天只能卸货 2 车,若一天内到达数超过 2 车,那么就推迟到次日卸货。根据表 3 所示的数据,货车到 达数的概率分布(相对频率)平均为 1.5 车/天,求每天推迟卸货的平均车数。9 l) g) F7 n9 u* R, @

# @) R7 ^8 a' t! m9 s7 m
4 X* U& _: W+ w/ @2 B
! S& X" W0 U5 K: }% q, I解 这是单服务台的排队系统,可验证到达车数不服从泊松分布,服务时间也不服 从指数分布(这是定长服务时间)。 随机模拟法首先要求事件能按历史的概率分布规律出现。模拟时产生的随机数与事 件的对应关系如表 4。
; n. `7 Y5 t2 ]# S' T1 \* O& ]0 v# U: w2 Y( h; ^3 }
# b) }! r5 r7 Y! P8 Q1 z3 G

$ l  }* T* B6 C& v我们用 a1 表示产生的随机数,a2 表示到达的车数,a3 表示需要卸货车数,a4 表 示实际卸货车数,a5 表示推迟卸货车数。编写程序如下:: u2 [) c  {: U0 H6 W3 g
7 T. W1 O9 x1 ~6 a: j% K
clear
) o/ |( m- v$ s. ]rand('state',sum(100*clock));+ k( E! O& ^7 H- e- A% H, P/ d' I3 b6 U  y
n=50000;
0 X4 }0 D% @0 h. E* F% Dm=27 a" ^4 z/ z* {. q# d
a1=rand(n,1);0 k& E/ l3 g- P7 x: d
a2=a1; %a2初始化
! }1 R. `# s" D1 u" E! M! e4 f1 Aa2(find(a1<0.23))=0;* t1 w5 z8 M7 q8 i$ R# j
a2(find(0.23<=a1&a1<0.53))=1;6 w- ^' z3 ^! ^) b
a2(find(0.53<=a1&a1<0.83))=2;3 x5 H5 F7 p* M% [/ Z
a2(find(0.83<=a1&a1<0.93),1)=3;$ t+ P$ o8 p! {
a2(find(0.93<=a1&a1<0.98),1)=4;) W4 L- d7 c5 K3 F
a2(find(a1>=0.98))=5;# I0 c! s5 l' s8 A$ O0 Q9 N
a3=zeros(n,1);a4=zeros(n,1);a5=zeros(n,1); %a2初始化- p- h/ W2 @( }5 v6 T. m9 f
a3(1)=a2(1);
4 I$ L& L( U  H9 Z9 Cif a3(1)<=m
8 B* K$ q! s* w  H" J) F    a4(1)=a3(1);a5(1)=0;% ]. b" P% K, K9 Z! K& s0 B2 h
else, g5 D4 H7 x  {' R0 T
     a4(1)=m;a5(1)=a2(1)-m;6 k3 I1 g" i$ F0 H! J* {6 |, [
end1 }$ c7 J/ Q/ q2 h8 b* f
for i=2:n
  P3 i0 o8 n% b; j% [4 \    a3(i)=a2(i)+a5(i-1);  X+ [, F4 \* E
    if a3(i)<=m4 N8 L5 i6 W) }8 v# |
        a4(i)=a3(i);a5(i)=0;
( {: S7 {9 d, a+ v( ~    else/ q6 [9 W+ G' a: |3 t8 m
        a4(i)=m;a5(i)=a3(i)-m;4 r8 F  C* _! B! @  H" h
    end
! S( {+ i8 p% _' Kend9 x7 k* [/ A" s! N) F' C2 u
a=[a1,a2,a3,a4,a5];
1 O3 R( U4 @& tsum(a)/n
& A! X* }$ g+ J+ n" Y. I6 o) d3 C例 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 次模拟, 进行比较。  A, |1 ~( p- p, E+ B
8 }7 I* m) y- j& m6 o

5 o  }! Y1 @+ e  u7 Z  H3 d/ i, x) B: G8 O/ v! H
在模拟 A 型机时,我们用cspan表示到达间隔时间,sspan表示服务时间,ctime 表示到达时间,gtime表示离开时间,wtime表示等待时间。我们总共模拟了m 次, 每次n 个顾客。程序如下:1 g; H0 ]% E2 z( [8 z7 o
% l  r  w7 U& \* D
tic0 h+ N/ y8 K% Z1 a5 N& c* u' o2 z
rand('state',sum(100*clock));; L1 Z8 }( y2 b% X& D
n=100;m=1000;mu1=1;mu2=0.9;5 v* {$ l- C/ w$ c
for j=1:m: F" a. O$ E; B& |( [
    cspan=exprnd(mu1,1,n);sspan=exprnd(mu2,1,n);7 X$ ?5 s7 ]- m6 a
    ctime(1)=cspan(1);2 }1 ~0 u  h4 v, I
    gtime(1)=ctime(1)+sspan(1);; c  o4 w- [& L! V) O$ {
    wtime(1)=0;8 Q- R; r* N" L% d9 U/ R  t$ F3 Y
    for i=2:n
7 D3 ^. s  u  F5 U, H        ctime(i)=ctime(i-1)+cspan(i);! Q  t0 I( ^" L! P
        gtime(i)=max(ctime(i),gtime(i-1))+sspan(i);
1 Q! ]- W) B* m) h$ \8 Y) s        wtime(i)=max(0,gtime(i-1)-ctime(i));
8 d$ I/ o9 N% X$ p- B) I    end; E' M6 b2 J& m$ |
    result1(j)=sum(wtime)/n;
! g4 Q% t# z# M% @, Z+ k8 ~end  J) g! Z: D- V# X1 ?- B
result_1=sum(result1)/m 0 f2 U8 J) `  N6 L5 r8 ~$ B
toc
6 d; _, G! X9 J( S+ n& H3 B类似地,模拟 B 型机的程序如下:
2 \" L; l* E  B: s* L. R- m" g
" _" @" d6 D9 i( [$ htic
% E9 b$ j: l- a+ n1 n, irand('state',sum(100*clock));" {& }" h. H' z
n=100;m=1000;mu1=1;mu2=1.8;
* n- n* N8 [8 s0 V( Rfor j=1:m" m' P. C' ?1 z
    cspan=exprnd(mu1,1,n);sspan=exprnd(mu2,1,n);
: H' W5 e, v" M/ I    ctime(1)=cspan(1);ctime(2)=ctime(1)+cspan(2);
  w: ~) f- H* z! y( k9 l    gtime(1:2)=ctime(1:2)+sspan(1:2);
+ S% E* L1 B1 K( y8 J    wtime(1:2)=0;flag=gtime(1:2);
1 G0 y6 U! X7 @7 ~5 v    for i=3:n" `8 O' {1 _- @; x4 W+ j
        ctime(i)=ctime(i-1)+cspan(i);$ d9 I% e) D" W, y5 ~
        gtime(i)=max(ctime(i),min(flag))+sspan(i);& H. W3 \# m- C: j. B
        wtime(i)=max(0,min(flag)-ctime(i));7 m2 |4 e' e; G) I- T
        flag=[max(flag),gtime(i)];3 _8 X7 L) @" ?
    end
* Y3 m1 R) q- z+ v: w    result2(j)=sum(wtime)/n;9 u0 c& E3 G8 R- N' X1 i9 O9 y
end/ }! x# _; V) ^! e
result_2=sum(result2)/m* I) S, v3 l" [* F% G( V) f& _
toc 2 E; o* I; ~, D" ~# ~! i/ F
读者可以用下面的程序与上面的程序比较了解编程的效率问题。
1 h3 B8 u  g2 p: A
1 a- ^0 V4 O1 t! O9 P2 xtic5 W9 ~' E. f2 T# d4 v, i
clear& _" \* G  A! I' Y
rand('state',sum(100*clock));+ a8 X; v+ P9 I% g) P
n=100;m=1000;mu1=1;mu2=0.9;/ @: s/ s# O6 q" P* v
for j=1:m
0 t; e9 W- k, `4 F    ctime(1)=exprnd(mu1);
0 \  z! Y: j; e7 p& u. P& r# j    gtime(1)=ctime(1)+exprnd(mu2);( m: a6 B! ]  R$ b+ p4 H% A
    wtime(1)=0;
. R$ I  A" F6 O2 W' ^    for i=2:n
& l6 f  @  P( }7 L( j6 @8 B        ctime(i)=ctime(i-1)+exprnd(mu1);
1 _/ g; _3 i9 d+ }        gtime(i)=max(ctime(i),gtime(i-1))+exprnd(mu2);$ x0 E9 f; k( t+ j, i
        wtime(i)=max(0,gtime(i-1)-ctime(i));
) u6 a. w  I$ V( {    end0 L0 a$ Y' h% S7 [
    result(j)=sum(wtime)/n;+ T: j1 y3 h* v0 `  l
end
  a1 m$ q2 X$ m  i  Nresult=sum(result)/m
1 ]; m2 D# W# l+ Utoc
1 r" i' a0 i$ O1 q+ [1. 一个车间内有10台相同的机器,每台机器运行时每小时能创造4元的利润,且平 均每小时损坏一次。而一个修理工修复一台机器平均需4小时。以上时间均服从指数分 布。设一名修理工一小时工资为6元,试求:
6 v. U* y2 |+ m/ C7 j6 E* u) r8 d4 r  f9 ~0 ^
(i)该车间应设多少名修理工,使总费用为最小;
5 Q; e1 h& J$ H1 I# a" D
$ ^7 s: V- ~& K' Z; F7 h(ii)若要求不能运转的机器的期望数小于4台,则应设多少名修理工;
6 c+ l' U: M* g( B. I, o
: \! K- ^; H. y4 X/ L(iii)若要求损坏机器等待修理的时间少于4小时,又应设多少名修理工。9 C+ s1 Z/ N8 G

$ Y8 o2 D- u: `9 }& Y/ A2. 到达某铁路售票处顾客分两类:一类买南方线路票,到达率为λ1 /小时,另一 类买北方线路票,到达率为λ2 /小时,以上均服从泊松分布。该售票处设两个窗口,各窗口服务一名顾客时间均服从参数 μ = 10 的指数分布。试比较下列情况时顾客分别等 待时间Wq :4 m- N7 G/ o9 b+ I6 B5 \

8 d# ?" L9 g  P9 r(i)两个窗口分别售南方票和北方票;: [! J8 B% X# x6 O( j+ Z
; _9 ]$ Z' Q# F* V  K
(ii)每个窗口两种票均出售。(分别比较 λ1 = λ2 = 2,4,6,8 时的情形). o' y/ a! A0 f1 S4 e
/ ~  [; m$ a$ c# {
3. 一名修理工负责5台机器的维修,每台机器平均每2h损坏一次,又修理工修复一 台机器平均需时18.75min,以上时间均服从负指数分布。试求:
2 W+ p' B! h( r9 W! \
9 S8 f' g, @' j7 O(1)所有机器均正常运转的概率;
5 R1 q$ D! h% Q! s1 f- t
6 P8 n5 ]! O% e# A(2)等待维修的机器的期望数;1 C$ Z/ m9 ?# Y. F9 j

. [9 x4 r6 g" t4 g2 L; h0 _8 b(3)假如希望做到有一半时间所有机器都正常运转,则该修理工最多看管多少台 机器。
6 I- x- {% _! R$ u
! \- C. B- ^4 C(4)假如维修工工资为8元/h,机器不能正常运转时的损失为40元/h,则该修理工 看管多少台机器较为经济合理。# y5 i! y' N  B* z0 y
————————————————
8 s, g9 O  |1 J7 Z( x版权声明:本文为CSDN博主「wamg潇潇」的原创文章,遵循CC 4.0 BY-SA版权协议,转载请附上原文出处链接及本声明。
" U3 H" i5 y$ u; N! r6 z6 D6 V原文链接:https://blog.csdn.net/qq_29831163/java/article/details/89738145
2 R4 G  H8 g9 `# Q/ p6 E$ [) c& I1 x3 \* f) n$ e
0 ^" z5 Z9 ~. f( Q) |1 u2 ^6 S) l





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