In this worksheet we consider the problem of servicing requests, following the example in section 2.2 of the textbook.
> restart;
> with(stats);
1. Distributions for arrival and service time
The arrival time is modeled by the continuous "Poisson flow" with probability density function exp(-t/a)/a. In Maple, this distribution is available as the "exponential distribution" with probability density function alpha*exp(-alpha*(x-a)) for x >= a and zero if x < a, where alpha is a non-negative real number. To relate the Maple exponentional function with the Poisson flow, we take a = 0 and alpha = 1/a.
> a := 10: # mean arrival time
> arrival_time := n -> stats[random,exponential[1/a,0]](n):
> arrival_time(3);
Let us look at the probability density function:
> plot(statevalf[pdf,exponential[1/a,0]],-10..50);
![[Maple Plot]](images/servicing3.gif)
The cumulative densitify function looks like
> plot(statevalf[cdf,exponential[1/a,0]],-10..50);
![[Maple Plot]](images/servicing4.gif)
From the plot above, we can see that the mean occurs where the value 0.5 is reached. Can you now see the meaning of the parameter a? Let us do an experiment.
> samples := arrival_time(1000):
> describe[mean]([samples]);
> describe[standarddeviation]([samples]);
The time to service a customer is normally distributed, with an average of 20 time units, and standard deviation equal to 1. Next we define a Maple procedure to generate n samples from this service time distribution.
> mu := 20: sigma := 1: # defines a normal distribution
> service_time := n -> stats[random,normald[mu,sigma]](n):
> service_time(3);
Now we are ready to perform Monte Carlo experiments.
2. Busy processors reject arriving jobs
Let us first consider the case of one processor. Above we generate three arrival times and three service times. If busy processors reject requests as they arrive, then we have to reject the second request, if it comes within the time it takes to serve the first request.
For several processors, we store for every processor the service time of the job being processed. An idle processor has service time zero.
If we conduct N experiments, we generate a random arrival and service time, according to the Poisson flow and normal distribution. At each arrival, the arrival time is the time elapsed, so we decrease the service times for the busy processors, eventually setting some time to zero. If all processors have positive times, we reject the job and increment the counter for the rejected ones by one.
> m := 3: # number of processors
> N := 10000: # number of simulations
> B := [seq(0,i=1..m)]: # none of the processors is busy
> rejected_ones := 0: # count rejected requests
> for i from 1 to N do # conduct N experiments
> A := arrival_time(1): # new arrival time
> S := service_time(1): # new service time
> for j from 1 to m do # update status of processors
> if B[j] < A
> then B[j] := 0: # release processor
> else B[j] := B[j]-A: # decrease service time
> end if:
> end do:
> reject := 1: # assume we reject the job
> for j from 1 to m do # try to assign job to processor
> if reject = 1 and B[j] = 0
> then B[j] := S: # assign job to j-th processor
> reject := 0: # if assigned, then no reject
> end if;
> end do: # now we can update the counter
> rejected_ones := rejected_ones + reject;
> end do:
Now we can see the ratio of rejected jobs versus the total number of experiments:
> evalf(rejected_ones/N);
Just as in the textbook we find that about 21% of the requests will be rejected.
This Monte Carlo experiment could model the situation in an emergency room treating seriously ill patients. Unless immedate care is given by a limited number of doctors, the patient dies. In a less gruelsome situation, like in a store, our customers can wait.
3. Waiting for busy processors
When all processors are busy, we assign the job to the processor with the smallest service time. The waiting time for this job is then the service time the processor had at the time of assignment. In the Monte Carlo experiment below, we compute the average waiting time.
> m := 3: # number of processors
> N := 10000: # number of simulations
> B := [seq(0,i=1..m)]: # none of the processors is busy
> waiting_time := 0: # sums total waiting time
> for i from 1 to N do # conduct N experiments
> A := arrival_time(1): # new arrival time
> S := service_time(1): # new service time
> for j from 1 to m do # update status of processors
> if B[j] < A
> then B[j] := 0: # release processor
> else B[j] := B[j]-A: # decrease service time
> end if:
> end do:
> min_time := min(op(B)): # look for minimum
> member(min_time,B,'p'): # processor with minimal time
> waiting_time := waiting_time + B[p]:
> B[p] := B[p] + S: # assign job to this processor
> end do:
Now we can ask for the average waiting time:
> evalf(waiting_time/N);
>