Lesson #12: More on Birth Death Simulations and M/G, G/M systems.

From Previous Lesson :

bdSim[birthFunction_,deathFunction_,timeHorizon_]:=Module[{state=0,time=0.0,index=1},
rel=Table[{kuku,bobo},{Round[4*2* (birthFunction[1]+deathFunction[1])*timeHorizon]}];
(*rel={{time,state}};*)
rel[[index++]]={time,state};
(*Now loop untill time > timeHorizon and append to rel as we go*)
While[time<timeHorizon,
birthRate=birthFunction[state];
deathRate=deathFunction[state];
time += RandomReal[ExponentialDistribution[birthRate+deathRate]];
rel[[index++]]={time,state};
(*AppendTo[rel,{time,state}];*)
state+=2RandomInteger[BernoulliDistribution[(birthRate/(birthRate+deathRate))]]-1;
rel[[index++]]={time,state};
(*AppendTo[rel,{time,state}];*)
];
Take[rel,index-1]
]

M/M/1

λ=0.4;
μ=1.0;
mm1Birth[s_]:=λ;
mm1Death[s_]:=If[s>0,μ,0];
rel=bdSim[mm1Birth,mm1Death,10];
ListPlot[rel,Joined->True]

Different Timing Attempts

Table[{t,First[Timing[bdSim[mm1Birth,mm1Death,t];]]},{t,500,5000,500}]

Table[{t,First[Timing[bdSim[mm1Birth,mm1Death,t];]]},{t,500,5000,500}]

Table[{t,First[Timing[bdSim[mm1Birth,mm1Death,t];]]},{t,500,5000,500}]

Table[{t,First[Timing[bdSim[mm1Birth,mm1Death,t];]]},{t,500,5000,500}]

ListPlot[%]

M/M/K/K  (Erlang Loss System)

λ=2.0;
μ=1.0;
KK=8;
mm1Birth[s_]:=If[s<KK,λ,0];
mm1Death[s_]:=μ s;
rel=bdSim[mm1Birth,mm1Death,30];
ListPlot[rel,Joined->True]

Joint Class Exercise

1) Collect Stationary Distribution from Simulation

bdSimWithStat[birthFunction_,deathFunction_,timeHorizon_]:=Module[{state=0,time=0.0},
rel={{time,state}};
histogramTable=Table[0.0,{100}];
(*Now loop untill time > timeHorizon and append to rel as we go*)
While[time<timeHorizon,
birthRate=birthFunction[state];
deathRate=deathFunction[state];
stateDuration=RandomReal[ExponentialDistribution[birthRate+deathRate]];
histogramTable[[state+1]]+=stateDuration;
time += stateDuration;
AppendTo[rel,{time,state}];
state+=2RandomInteger[BernoulliDistribution[(birthRate/(birthRate+deathRate))]]-1;
AppendTo[rel,{time,state}];
];
{rel,(histogramTable/timeHorizon)}
]

M/M/1 Histogram

λ=0.9;
μ=1.0;
mm1Birth[s_]:=λ;
mm1Death[s_]:=If[s>0,μ,0];
out=bdSimWithStat[mm1Birth,mm1Death,4000];
hist=Select[Transpose[{Table[i,{i,0,Length[out[[2]]]-1}],out[[2]]}],#[[2]]>0&];
theoreticalHist=Table[{i,(1-(λ/μ))((((λ/μ)))^i)},{i,0,Round[1.3*Length[hist]]}];

ListPlot[{hist,theoreticalHist},Filling->Axis,Joined->True,PlotRange->All,PlotStyle->PointSize[0.02]]
ListPlot[out[[1]],Joined->True]

M/M/∞ Histogram

λ=3.0;
μ=1.0;
mm1Birth[s_]:=λ;
mm1Death[s_]:=μ s;
out=bdSimWithStat[mm1Birth,mm1Death,500];
hist=Select[Transpose[{Table[i,{i,0,Length[out[[2]]]-1}],out[[2]]}],#[[2]]>0&];
theoreticalHist=Table[{i,(E^(-λ/μ))((((λ/μ))^i)/(i !))},{i,0,Round[1.3*Length[hist]]}];

ListPlot[{hist,theoreticalHist},Filling->Axis,Joined->True,PlotRange->All,PlotStyle->PointSize[0.02]]
ListPlot[out[[1]],Joined->True]

M/M/K/K Histogram

λ=10.0;
μ=1.0;
K=10;
mm1Birth[s_]:=If[s<K,λ,0];
mm1Death[s_]:=μ s;
out=bdSimWithStat[mm1Birth,mm1Death,300];
hist=Select[Transpose[{Table[i,{i,0,Length[out[[2]]]-1}],out[[2]]}],#[[2]]>0&];
theoreticalHist=Table[{i,(1/(Sum[((((λ/μ)))^j)/j!,{j,0,K}]))((((λ/μ))^i)/(i !))},{i,0,K}];

ListPlot[{hist,theoreticalHist},Filling->Axis,Joined->True,PlotRange->All,PlotStyle->PointSize[0.02],AxesOrigin->{0,0}]
ListPlot[out[[1]],Joined->True]

2) Sojourn Times

bdSimWithStatAndSojourns[birthFunction_,deathFunction_,timeHorizon_]:=Module[{state=0,time=0.0},
rel={{time,state}};
soujournList={};
(*queue is a list of arrival times of customers*)
queue={};
histogramTable=Table[0.0,{100}];
(*Now loop untill time > timeHorizon and append to rel as we go*)
While[time<timeHorizon,
birthRate=birthFunction[state];
deathRate=deathFunction[state];
stateDuration=RandomReal[ExponentialDistribution[birthRate+deathRate]];
histogramTable[[state+1]]+=stateDuration;
time += stateDuration;
AppendTo[rel,{time,state}];
change=2RandomInteger[BernoulliDistribution[(birthRate/(birthRate+deathRate))]]-1;
state+=change;
If[change==1,
AppendTo[queue,time]
,
(*Else change==-1*)
AppendTo[soujournList,time-queue[[1]]];
queue=Drop[queue,1]
];
AppendTo[rel,{time,state}];
];
{rel,(histogramTable/timeHorizon),soujournList}
]

Table[
λ=0.3;
μ=1.0;
mm1Birth[s_]:=λ;
mm1Death[s_]:=If[s>0,μ,0];
out=bdSimWithStatAndSojourns[mm1Birth,mm1Death,3000];
out[[3]] //Mean,
{10}]

The M/G/1 Queue - The Pollaczek-Khinchin Formulas

E[W]=(ρ/(1-ρ)) (1/2)(((Subscript[c, s])^2)+1)(1/μ)

Stationary distribution of W and/or  N  is given in form of LSTs,

(SuperStar[W])(s)=(((1-ρ)s)/(λ (SuperStar[G])(s)+s-λ))

Similar formula for stationary distribution of Number in System

The G/M/1 Queue

Here we need to solve the equation

σ=(SuperStar[A])(μ-μ σ)

The Queue length that an arriving customer sees is (Subscript[Geometric, 0])(σ)

Clear[λ,μ,K]

A1[s_]:=(λ/(λ+s))

A2[s_]:=(((((2λ)/(2λ+s))))^2)

-D[A1[s],s]  /.s->0

-D[A2[s],s]  /.s->0

Plot[{A1[s],A2[s]} /.λ->1/2,{s,0,1.1},
AxesOrigin->{0,0},Epilog->{Line[{{0,1},{1,0}}],
Line[{{0,1},{1/2,0}}]}]

A[s_]:=(((((kk λ)/(kk λ+s))))^kk)