5.1.1 Simulation is Imitation

Simulation is one of the key techniques used in scientific research as well in everyday decision making practice. It ranks very high among the most widely used approaches for modeling systems and phenomena and making decisions regarding their design, and prediction of their functioning. It involves using a computer to imitate (simulate) the operation of a process or system (S). It can be very helpful for:

1. risk analysis on financial processes by repeatedly imitating the evolution of the transactions involved to generate a profile of the possible outcomes

2. stochastic systems that will continue operating indefinitely

3. systems where it is too expensive or risky to do live tests. Simulation provides an inexpensive, risk -free way to test changes ranging from a “simple” revision to an existing production line to emulation of a new control system or redesign of an entire supply chain.

4. large or complex systems for which change is being considered. A “best guess” is usually a poor substitute for an objective analysis. Simulation can accurately predict their behavior under changed conditions and reduce the risk of making poor decisions

5. systems where predicting process variability is important. A spreadsheet analysis cannot capture the dynamic aspects of a system, aspects which can have a major impact on system performance. Simulation can help you understand how various components interact with each other and how they affect overall system performance.

6. systems where you have incomplete data. Simulation cannot invent data where it does not exist, but simulation does well at determining sensitivity to unknowns. A high - level model can help you explore alternatives. A more detailed model can help you identify the most important missing data.

5.1.2 Simulation outline

When simulation is used to facilitate modeling and decision making regarding the design/operation of a system then the following steps are followed:

1. firstly a preliminary analysis is done (with an appropriate mathematical model) to develop a rough design of the system (S) or the process (including its operating procedures).

2. Then simulation is ran to experiment with different designs and estimate how well each will perform.

3. After a detailed design is developed, the system is tested in actual use and final touches to the design are made.

5.1.3 Simulation model building blocks

To prepare for simulating a complex system, a detailed simulation model needs to be developed to describe the operation of the system and how it is to be simulated. We use a queueing system M/M/1 to illustrate the simulation model building blocks. Firstly we will briefly introduce/recall what M/M/1 is.

M/M/1 queueing system

Random arrivals/jobs join a service unit with single server, that provides them with a certain type of service. The interarrival times between the job arrivals is assumed to be a random variable \(T\) with exponential distribution, \(T \sim exp(\lambda)\). The service time required by a job is also random variable, say \(S\) and it is assumed that \(S\) is exponential (\(S \sim exp(\mu))\). The service unit is called queueing system and its diagram is depicted above. Therefore

1. the system has a single-server,

2. the interarrival times \(T\) are iid exponential with parameter \(\lambda\) and the mean and the variance of the interarrival time are

\[E(T) = \frac{1}{\lambda},\hspace{1cm}Var(T) = \frac{1}{\lambda^2}.\]

3. exponential service with rate \(\mu\), i.e., the service times \(S\) are iid exponential with parameter \(\mu\) and the mean and the variance of the service time are\

\[ E(S) = \frac{1}{\mu}, \hspace{1cm} Var(S) = \frac{1}{\mu^2}.\]

For M/M/1 queues we can compute their performance measures analytically. For example, the average number of jobs in the system in steady-state is equal to
\[\begin{align*} L &= \frac{\rho}{1-\rho}, &\mbox{where} \hspace{3ex}\rho= \frac{\lambda}{\mu} & \;\;\; \mbox{and} & \rho< 1. \\ \end{align*}\]

Now, we will use the above M/M/1 model to illustrate the building blocks for its simulation model.

1. A definition of the state of (S) - e.g., number of customers in the system

2. Identify the set of possible states of (S) - e.g., if the system cannot have more that 10 customers, so the set of possible states for this queueing system is \(\{0, 1, 2, \ldots, 10 \}\)

3. Identify the possible events that would change the state of (S) - e.g, an arrival or service completion (departure) in a queueing system

4. Simulation clock (measures the simulation time)

5. Method of randomly generating events of various kinds

6. Identification of the state transitions related to the events

5.1.4 Discrete Events Simulation

A change in the state of the (S) occurs instantaneously at random point in time as a result of the occurrence of discrete event.

For example, in a queueing system, as defined earlier, the discrete events that change the state of the system are the arrival of a customer and the departure of a customer due to the completion of its service.

5.1.5 Continuous Simulation

Changes in the state of (S) occur continuously over time

For example, if the system of interest is an airplane in flight and its state is defined as the current position of the airplane, then the state is changing continuously over time. This type of simulation requires differential equations to describe the rate of change of the state variables. The analysis is relatively difficult.

Sometimes, it is possible to approximate the continuous changes in the state of the system by discrete changes, so that a discrete events simulation could be used to approximate the behaviour of a continuous system. This greatly simplifies the analysis.

5.2.1 Simple simulation

You could play the game by yourself many time. One hour spent in repeatedly flipping a coin and recording the earnings or losses might be sufficient to make a decision whether to play this game for money. This is a simulation because you will be imitating the actual play of the game without actually winning or losing any money.

5.2.2 A coin-flipping game: Computer simulation

We need to identify the building blocks of our simulation model.

1. The state \(N(t)\) of the system at time \(t\) is the difference between the number of \(H\) and the number of \(T\) after the \(t^{th}\) coin-flip. \[\begin{equation} \begin{split} N(t) &= \#H - \#T \hspace{2ex}\mbox{after} \hspace{0.5ex} t \hspace{0.5ex} \mbox{flips}, \end{split} \end{equation}\]

2. The set of possible states of (S) \[ N(t) = \{-3, -2, -1,\, 0,\, 1,\, 2,\, 3 \}\]

\(N(t) = -3\) means that the play will end, because there are three more tails, and \(N(t) = 3\) means that the play will end due to having three more heads.

3. The possible events that would change the state of (S) are the outcomes of a single coin-flip \[H \hspace{2ex}\mbox{or} \hspace{2ex} T\]

4. Simulation clock records the number of flips, say \(t\)

5. Method of generating events - generate a uniform random number \((0,1)\). If this number is between \((0.0, 0.5)\) we say that the outcome of the coin-flip is a head (H); if it is within \((0.5, 1.0)\) the outcome is a tail (T), i.e.,

6. State transition - we increase the state by 1 if the outcome of the coin-flip is a \(H\) and reduce it by 1 if the outcome is a \(T\), i.e.,

The game ends when \(N(t) = \pm 3\). Your ``profit’’ is \((8 - t)\) ( it could be positive (you win!) or negative (you lose!)).

Here is a possible R code to simulate the game. So, what is your decision - to play or not to play?

oneplay <- function(){
   nt=0 #state of the sytem
   i=0  # number of flips
   while(-3 < nt && nt < 3){
    rd = runif(1)
    i = i+1 # increase the number of flips
    if(rd < 0.5) nt = nt+1 # the outcome is H, increase by 1
    else nt = nt-1         # the outcome is T, decrease by 1
    }
   payoff= 8-i
   return(payoff)
}

r = rep(0,100000)
for(i in 1:100000){ set.seed(111*i); r[i]<-oneplay()}
av = mean(r)
cat("Average Payoff = ","",av,"\n")

## Average Payoff =   -1.03234

For long runs the average payoff is \(\approx\) -1$. So, should you play this game?

5.3 Monty Hall problem

Suppose you’re on a game show, and you’re given the choice of three doors:

1. Behind one door is a car; behind the others - goats.

2. You pick a door, say No. 1,

3. The host, who knows what’s behind the doors, opens another door, say No. 3, which has a goat.

4. Then he says to you, “Do you want to pick door No. 2?”

5. Is it to your advantage to switch your choice or stick with your initial door choice?

The rules of the game are as follows:

1. You make an initial guess

2. The master opens one of the ``goat’’ doors

3. You have a chance to revise your initial guess before making your final choice.

While revising your initial guess you might select one of the following strategies:

1. strategy A: stick with your initial guess
2. strategy B: switch to another door

Simulate the game! Which strategy (A or B) do you prefer? What is the criterion for your decision?

#no switch strategy
N <- 100000 #number of plays
s <- 0 # counter for wins
for (i in 1:N){
# select the door with the car
DC <- sample(c(1,2,3), 1, replace = FALSE, prob=c(1/3,1/3,1/3))
# select the door of the player choice
C1 <- sample(c(1,2,3), 1, replace = FALSE, prob=c(1/3,1/3,1/3))
# the strategy is "no switch", if the two doors are the same (DC==C1), then the win
if(DC == C1) rez <-1
else rez <- 0

s <- s+rez
}
#computing the probability of win
PW = s/N

#switch strategy

s <- 0
for (i in 1:N){
# select the door with the car
DC <- sample(c(1,2,3), 1, replace = FALSE, prob=c(1/3,1/3,1/3))
# select the door of the player choice
C1 <- sample(c(1,2,3), 1, replace = FALSE, prob=c(1/3,1/3,1/3))
# the strategy is "switch", if the two doors are the not the same (DC!=C1), after switch win

if(DC!= C1) rez <- 1
else rez <- 0

s <- s+rez
}
PW1 <- s/N
cat("ProbWin(Switch) = ", PW1, " ","ProbWin(No Switch) = ", PW, "\n")

## ProbWin(Switch) =  0.666   ProbWin(No Switch) =  0.332

So, the strategy of switching (B) is preferable because the probability of winning the car is much higher.

5.4 The newsboy Freddie problem

The newsboy Freddie faces the following problem. One of the daily newspapers that Freddie sells from his newsstand is the Financial Journal. The distributor brings the day’s copies of the _Financial Journal to the newsstand early each morning. Any copies unsold at the end of the day are returned to the distributor next morning. However, to encourage ordering a large number of copies, the distributor does give a small refund for unsold copies. Here are Freddie’s cost figures:

1.Freddie pays 1.50 per copy delivered;

2.Freddie sells it for 2.50 per copy;

3.Freddie’s refund is 0.50 per unsold copy.

Freddie sells anywhere between 40 and 70 copies inclusively on any given day. The frequency of the numbers between 40 and 70 are roughly equal. The decision that Freddie needs to make is the number of copies to order per day from the distributor, so that to maximize his average daily profit.

You might find the command which.max() useful. It returns the location, i.e., index of the (first) minimum or maximum of a numeric (or logical) vector.

You might want to start by writing a fragment of R code that calculates Freddie’s profit given that he orders, say x number of newspapers, where \(40 \le x \le 70\). Then calculate his average profit for this order of x newspapers. Repeat these calculations for all possible values of x, \(40 \le x \le 70,\) compare the results and make your recommendation on the number of newspapers to order, so to maximise Freddie’s average daily profit.

Here is my R code to simulate this situation:

x - amount ordered, a decision variable

y - the demand - uniform r.v. U[40, 70]

1.5 cost per unit ordered

2.5 sale price per unit

0.5 refund per unsold copy

set.seed(123)
k=5000
rez<-rep(NA,31)

for(x in 40:70)
{S=0
  for (i in 1:k)
  {
    y <-sample(40:70, 1)
    Prof = -1.5*x + 2.5*min(x,y) + 0.5*max(x-y,0)
    S = S + Prof
      }
   rez[x-39] <- S/k
}

#for (i in 1:31)
#{cat(i+39, " ", rez[i], "\n")}

cat("Max profit ",max(rez)," for order ", which.max(rez)+39,"\n")

## Max profit  47.3548  for order  56