Modeling time-varying loading docks

IN

DEGREE PROJECT MATHEMATICS, SECOND CYCLE, 30 CREDITS

STOCKHOLM SWEDEN 2018,

Modeling time-varying loading docks

CARL FLOGVALL

Modeling time-varying loading docks

FLOGVALL CARL

Degree Projects in Optimization and Systems Theory (30 ECTS credits) Degree Programme in Applied and Computational Mathematics (120 credits) KTH Royal Institute of Technology year 2018

Supervisor at LUP Technologies: Jacob Armé Supervisor at KTH: Per Enqvist

Examiner at KTH: Per Enqvist

TRITA-SCI-GRU 2018:053 MAT-E 2018:19

Royal Institute of Technology School of Engineering Sciences KTH SCI

SE-100 44 Stockholm, Sweden URL: www.kth.se/sci

Modeling time-varying loading docks

Carl Flogvall

May 8, 2018

Abstract

In this paper several Markov models have been created for simulating queues at loading docks with discrete time Markov chains and a optimization algorithm in regard to server capacity has been presented. The queuing systems have been modeled with exponential and Erlang-2 distributed service rate. They have also been modeled for one class or a simplification of two classes. An attempt was made to find a scaling factor for the service rate to use the one class exponential model to replicate the other models. This was sufficiently successful done with derivation and linear regression.

The optimization was done with a greedy approach. It was computationally efficient, but could not guarantee global optimality.

Modellering av tidsvarierande lastkajer

Carl Flogvall

8 maj 2018

Sammanfattning

I den här studien har flera Markov-modeller skapats för att simulera av köer vid lastkajer med hjälp av diskreta tidsmarkovkedjor och en optimeringsalgo- ritm med avseende på betjäningskapacitet har presenterats. Kövsystemen har modellerats med exponentiell och Erlang-2-distribuerad betjäningsintensiteten.

De har även modellerats för en klass eller en förenkling av två klasser. Ett försök gjordes för att hitta en skalfaktor för betjäningsintensiteten för att använda den en klass exponentiella modellen för att replikera de andra modellerna. Detta var tillräckligt framgångsrikt gjort med härledning och linjär regression.

Optimeringen gjordes med en girig inställning. Det var beräkningsmässigt ef- fektivt, men kunde inte garantera global optimalitet.

Acknowledgments

I would like to thank my coworkers at LUP Technologies for the opportunity to write this thesis. I would also like to thank my supervisor Per Enqvist for the support and help to complete this thesis. Lastly, I am thankful for family and friends who have supported me.

Contents

1 Introduction 4

1.1 Problem Description . . . 4

1.2 Literature Review . . . 5

2 Theoretical Background 6 2.1 Probability Distributions . . . 6

2.1.1 Poisson Distribution . . . 6

2.1.2 Exponential Distribution . . . 7

2.1.3 Erlang-k Distribution . . . 7

2.2 Poisson Processes . . . 8

2.3 Discrete-time Markov chain . . . 9

2.4 Queuing Theory . . . 11

2.5 Marginal Allocation . . . 11

2.6 Loading docks . . . 13

3 Queue models 14 3.1 Exponential distributed service rate . . . 15

3.1.1 One class . . . 15

3.1.2 Two classes . . . 15

3.2 Erlang-2 distributed service rate . . . 17

3.2.1 One Class . . . 17

3.2.2 2 Classes . . . 18

3.3 Size differences between models . . . 19

4 Optimization 20 5 Results and analysis 22 5.1 Difference between Exponential and Erlang-2 distributed service rate for one class . . . 23

5.2 Difference between one and two classes for exponential distributed service rate . . . 26

5.3 Difference between one class exponential distributed service rate and two classes Erlang-2 distributed service rate . . . 28

5.4 Optimization . . . 30

6 Conclusions 32

A Results 33

A.1 Scaling Markov Erlang multiple servers . . . 33

B Derivations 35

B.1 Scaling factor for Markov-Erlang one server . . . 35 B.2 One class-Two classes . . . 36

C Bibliography 38

Chapter 1

Introduction

In today’s global world queues are everywhere and affects us daily. Fast food restaurants, call centers, public transports are some examples where queuing systems are present. The area is well researched, but there are many aspects to consider and model after.

1.1 Problem Description

Manufacturers usually do not handle their own shipment, but instead uses a freight forwarder. This mean there is a disconnect of arrival information be- tween the manufacturer and the truck driver. The arrival of trucks is uncertain and the distribution of them is often uneven through out the day. This can create problems on the loading docks with truck queues through out the day of misused personnel on low-intensity periods. The question is what is the best way to solve this problem? The problem is in this way two-folded in modeling and optimization

First, one has to simulate a queue where the arrival rate, service rate and amount of servers can be time-varying to be able to replicate the reality. Further more, there are complications as the service rate’s different probability distributions and modeling the possible build up of trucks during the days.

Secondly, how can one optimize, based on the simulation, to have enough per- sonnel and allocate the capacity during peak periods? How do you optimize integer capacity for queues and does it have a global optimal solution?

1.2 Literature Review

The area of modeling queues is well researched, but very substantial and com- prehensive. Due to the problem formalization of time-varying arrival rate, the main parable was call centers.

The Analysis of Queues with Time-Varying Rates for Telecommunication Mod- els by Massey [2] studies how queues functions during time-varying arrival rates with both delay and loss models for telecommunications. The paper also han- dles capacity of staffing and measuring the time lag between the peak arrivals and peak load for a system.

Kagwai and Takagi’s work on call centers [3] is about modeling time-varying call centers with fluid approximation technique. The analysis also handles after-call work and different distributions like exponential, hyper-exponential and hypo- exponential distribution for the service rate.

Al Hanbali. Alvarez and van der Heijden [8] writes about waiting time dis- tribution in priority queues which are applied methods to real-life data in a case study. They also write about scaling the service time distribution for faster computations.

Both Rolfe [6] and the duo Shanthikumar, J. George, and David D. Yao [6]

wrote papers on optimizing server allocation for a system of multi-server sta- tions. The optimizations are constrained by fixed buffer capacity or manpower and find numerical evidence to support optimal solutions.

Chapter 2

Theoretical Background

2.1 Probability Distributions

2.1.1 Poisson Distribution

The Poisson distribution is a discrete distribution which is mostly used to de- scribe the number of events in a given time period. The distribution is under the assumption the events are independent of each other and occur with a constant average rate λ.

Figure 2.1: Poisson probability distribution function The mean of distribution is µ = λ and the variance is σ²= λ

2.1.2 Exponential Distribution

The exponential distribution describes the time until an event occurs in a Pois- son process. The notation of it is Exp(λ) where λ is the intensity of an event occurring. The exponential distribution is also the only continuous distribution which has the memoryless property.

P (X ≤ t + T |X > t) = P (X ≤ T ). (2.1)

Figure 2.2: Exponential probability distribution function

The mean of the distribution is µ = ¹_λ and the variance is σ²=_λ¹₂.

2.1.3 Erlang-k Distribution

The Erlang distribution is defined by the two parameters, a positive integer k (shape) and a positive λ (rate), where the shape determines the denomination of the distribution. For example, an Erlang-k distribution with k equal to 3 is called an Erlang-3 distribution.

The Erlang distribution is related with the exponential distribution in the fol- lowing way;

If

Xi∈ Exp(λ), (2.2)

then

k

X

i=1

Xi∈ Erlang(k, λ). (2.3)

Figure 2.3: Erlang probability distribution function

The mean of the distribution is µ =^k_λ and the variance is σ²= _λ^k2

2.2 Poisson Processes

A Poisson process is a stochastic process, denoted as N (t), t ≥ 0, for describing the number of events during the time period [0, t]. These events occur for a fixed arrival rate and N(t) has independent increments. The number of events are distributed for Poisson(λt) where λ is the arrival rate. Examples of areas where Poisson processes has been used is queuing theory, insurance rates, and sports [1].

If X_tis a continuous-time stochastic process at the time t ≥ 0, then let

0 ≤ t1< t2< ... < tn < ∞. (2.4) If

X_t2− Xt1, X_t3− Xt2, ..., X_tn− Xtn−1 (2.5) are independent, Xtis a stochastic process with independent increments.

Poisson processes has the attribute of being memoryless. It is defined as P (T > t + s|T > s) = P (T > t) ∀t, s > 0. (2.6) It means the probabilities does not change regardless on how much time has passed. The process’s future is independent of its history

2.3 Discrete-time Markov chain

A discrete-time Markov chain is a discrete stochastic process with regard to state space transitions. For it to be a Markov chain it has to fulfil the Markovian property as follows. Let the process be denoted by X_n and i_n be the state i at the time n if X_n=i_n. Then the following has to hold for i_n, j and n = 0, 1, 2...

P (X_n+1= j|Xn = i, Xn−1= in−1, ..., X₁= i1, X₀= i0) = P (Xn+1= j|Xn= i) (2.7) This shows the Markov chain is memoryless as the probabilities are only depend on what state the system is in and not the past events.

The probability of a state transition is defined as

P (Xn+1= j|Xn = i) = pij (2.8) pij is the probability of moving from state j to state i, whereof pii is the probability of staying in state i. However the following must hold

p_i,j≥ 0∀i, j (2.9)

X

j

pi,j= 1∀j (2.10)

P is a transition matrix with determines the probability of transition between the N states at a given time.

P =















p11 p12 p13 · · · p1N

p21 p22 p23 · · · p2N

p31 p32 p33 · · · p3N

... ... ... . .. ... pN 1 pN 2 pN 3 · · · pN N















(2.11)

x is a state vector which shows the probability of the system being in each state.

If x has N states, then it is defined as

x =













 p1

p2

p3

... p_N















(2.12)

where pj is the probability of the system being in the state j. Like for the transition matrix the two following must hold.

N

X

j=1

pj = 1 (2.13)

pj≥ 0∀j (2.14)

With the state vector and the transition matrix it is possible to calculate the state vector for the next time step

x⁽ⁿ⁺¹⁾= x⁽ⁿ⁾P (2.15)

The generalized equation for calculating the state vector m time steps away given m=1,2,3.. and P is constant.

x^(n+m)= x⁽ⁿ⁾P^m (2.16)

2.4 Queuing Theory

Queuing theory is the application of mathematics for simulating the behavior of waiting lines. Steady state solutions can be used to predict for example ex- pected queue length and waiting times. Markov processes are often used to model them. There are different types of queuing systems, but in this paper only FIFO (First in, first out) systems will be handled.

0 1 2 3 ...

λ

µ

λ

µ

λ

µ

λ

µ Figure 2.4: Queuing model

For a birth-death process where the state variable only increase and decrease by one, the queue can be described by the A/B/C notation. A stands for the distribution of arrival time, B the distribution of service time and C is the number of parallel servers. Systems with exponential distributed arrival and service time is indicated as M/M/C. Such systems can be modeled with Markov chains. The notations for the queuing systems are described in Table (2.1) below.

λ Arrival rate µ Service rate ρ Utilization rate E(W ) Expected queue length

E(n) Expected number of units in system Table 2.1: Properties of queuing system

2.5 Marginal Allocation

The purpose of marginal allocation is to set side discrete resources to a process to maximize/minimize one or several integer-convex objective functions.

For several objective functions the problem is often in compromising and weight- ing between them. Given two integer-convex functions fj(xj) and gj(xj) where f (xj) is strictly decreasing and g(xj) is strictly increasing.

f (x) =

n

X

j=1

fj(xj), (2.17)

g(x) =

n

X

j=1

gj(xj), (2.18)

X = x ∈ Nⁿ|g(x) ≤ g^max. (2.19)

For small xj f (xj) is big but g(xj) small and vice versa for big xj. Often the problem defined as

M = {(g(x), f (x))|x ∈ X} ⊂ R² (2.20)

where M is the finite set of solutions for the problem. To illustrate M, one can plot the points in a coordinate system with g(x) against f (x).

Figure 2.5: Graph of the convex hull

The convex hull of M is defined as the smallest convex set in R²which contains M. The optimal points of problem is in the convex hull, where the optimal solu- tion is depended on the weighting between the functions. The curve in Figure (2.5) is the efficient curve which is the a piecewise linear curve of the convex hull.

The algorithm to calculate the efficient curve is as follows by Krister Svanberg [4]

Step 1:

Calculate the quotients ^−∆f_∆g^j(rj)

j(rj) for each j and r_j = 0, 1, 2, ..., r^max_j where r_j^max is sufficiently big number to reach g^max. Put these quotients in a table with j in each the columns. Set k = 0, x(0) = (0, ..., 0)^T, g(x(0)) = g(0) and f (x(0)) = f (0).

Step 2:

Select the largest quotient in the table with has not been selected before. The

chosen quotient will now be allotted as l.

Step 3:

Let xl(k + 1) = xl(k) + 1 and x(k + 1)j = xj(k) for all j 6= l. Calculate the new values for f (x(k + 1)) and g(x(k + 1)) as f (x(k + 1)) = f (xl(k)) + ∆fl(xl(k)) and g(x(k + 1)) = g(x(k)) + ∆gl(xl(k)).

Step 4:

Terminate the algorithm if g(x(k)) ≥ g_max, else set k = k + 1 and go to step 2.

2.6 Loading docks

This paper will mainly focus on the loading dock setup to handle semi trailer trucks , as it is the most common type and there are many different types of trucks with different setups. The normal procedure of loading a semi trailer is relative simple. Every location have a number of loading docks and each dock has a number of loading zones for where trucks can load. This mean there is a maximum capacity for how most trucks can simultaneity load at each dock.

The two most common types of loading semi truck trailers are loading from back and loading from the side. Back loading is done by opening the hatch on back side, while side loading is done by removing the curtain sided wall and load through it. Most commonly a loading zone has a setup to only handle one type of loading, but there a some mixed capacity.

The process for an arriving truck is comparable to every other queuing sys- tem. Firstly the trucks arrives at the location where they declare their cargo number at the gate. The gatekeeper then tells them which dock to drive to and let them in through the gate. At the dock, the trucks park at an availability loading zone and open the trailer. If all the loading zones are unavailable, the trucks will usually have to wait in a first-come, first -serve queue until a loading zone is available. When it is the trucks turn to be loaded, the truck driver specifies his cargo number to the loader. The loader picks up the cargo and load it onto the truck. When the loading is completed, the truck driver secures the cargo and drive off.

A loading can take any where from couple of minutes to several hours depending on the cargo. For a specified type cargo, the distribution in loading time can also differ depending of the cargo, but an Erlang-n distribution is a common approximation.For example, de Moine’s report[5] found the loading of ship to have an Erlang-2 distribution

Chapter 3

Queue models

Ordinary one could use stationary queuing theory to simulate the queue by dividing the day into intervals and calculating the congestion for each part.

However, it doesn’t work for this problem since the service time is long relative the variation of arrival rate. This mean the queue during the morning affects the queue during the afternoon.

To handle this, a discrete-time Markov process is used. This allows a previ- ous to state affect the current state. The current state is calculated as

x(t + δt) = x(t)P (t, δt), (3.1) where P is built with time dependent arrival rate and depart rate. The are several limitations with this method such as δt has to be small enough for P to only consist of non-negative real numbers and also be a correct equitable approximate of the real continuous system. Furthermore, since the discrete sys- tems has a finite number of states, a maximum capacity of the system has to be set.

The queue models is meant to represent the loading docks where the service is the loading/unloading of a truck. The constant factors between the models will be first-in, first-out for scheduling and a Poisson process for arriving trucks.

The service rate for the models will either be exponential or Erlang-2 distribu- tion. The reason is to examine the difference between them and to see if it is possible to reproduce the result of the simpler exponential distribution model for the more complex model of a Erlang-2 distribution. For the real-life loading dock, the staff will often handle different kind of loadings as it can be different type of goods or loading. Therefore the models will be either one class or two classes. They will also be compared between each other and tried to be reduced for computational speed.

3.1 Exponential distributed service rate

3.1.1 One class

The most common queue model is one class with exponential distributed service rate. This model is, compared to others, relative straightforward. The arrival of a truck increases the state with 1 and the service of a truck decreases the state with 1.

0 1 2 ... c-1 c c+1 ... n-1 n

λ λ

µ

λ

2µ

λ

3µ

λ

(c − 1)µ

λ

cµ

λ

cµ

λ

cµ

λ

cµ cµ

Figure 3.1: One class and exponential distributed service rate queuing model Due to the simplicity of the model, the transition matrix will follow this general form for a M/M/c/n queue.

P =

















1 − λ λ

µ 1 − λ − µ λ

2µ 1 − λ − 2µ λ

.. .

.. .

.. .

(c − 1)µ 1−λ−(c−1)µ λ

cµ 1 − λ − cµ λ

.. .

.. .

.. .

cµ 1 − λ − cµ λ

cµ 1 − cµ

















(3.2)

3.1.2 Two classes

Sometimes the system will handle more than one type of products. These types can have both different arrival and service rates. Ideally one would model the system for two classes were every truck is accounted for. The problem with this is the size of the system becomes very big very quickly. Each line configuration would be a state, which would make the model double in size for every truck in max capacity added. For big system this is unsustainable, the computational effort would be to great. Instead a reduced model is presented were only the amount of trucks in the line is known and not its configuration. This is viable since only the configuration of the trucks being handled/served impact the sys- tem. The trucks in line will only impact in the system in the future.

0

X

Y

XX

XY

YY

XX+1 XX+2

XY+1 XY+2

YY+1 YY+2

λ^X

λY

λ^X

λ_Y µ^X

λ^X

λ_Y µ_Y

λX+λY

2µ^X

λX+λY

pX2µX

p^Y2µ X

pX2µX

p^Y2µ X

λX+λY

µ^X

µY λX+λY

p_X µ_Y

pXµX+pYµY

p^Yµ^X

p_X µ_Y

pXµX+pYµY

p^Yµ^X

λX+λY 2µY

λX+λY pX

2µY

pY2µY p_X

2µ_Y

pY2µY λ^X

Figure 3.2: Two classes and exponential distributed service rate queuing model The reduced model only contain the combinations of trucks served and how many are queuing When a customer is served the probability the next one is a type 1 is calculated as p₁=_λ^λ1

1+λ2. This can create a problem if the arrival rate is varying between the types. The queue can be build up with one type, but later due to the varying arrival rate, the model serves another type from what the queue created from.

For two types the size of the queuesing system becomes ∼ (c + 1)n states where c is the number of servers and n is the capacity. This holds as long as (c << n).

For k classes, the size would increase to ∼ C(c + k − 1, c)n classes.

3.2 Erlang-2 distributed service rate

3.2.1 One Class

For some queuing system the loading/unloading time is not exponentially dis- tributed but have a form of Erlang distribution. Fortunately for the Erlang distribution it is possibly for Markov chains due to the relation between Erlang and exponential distribution. Erlang-2 is chosen for the model since has a suf- ficiently close probability density function to real data and it is the simplest Erlang distribution. For Erlang-2 the simplest way to look it is the service is two folded, as if it goes through the process two times. For X1 the truck is under the first process and X2 is under the second.

0 X1

X2

X1X1

X2X1

X2X2

X1X1

+1 X2X1

+1 X2X2

+1

X1X1

+2 X2X1

+2 X2X2

+2

λX µX

λX µ^X

λX

2µX

λX µ^X

µX

λX 2µ^X

λX

2µX

λX µ^X

µX

λX 2µ^X

λX

2µX µ^X

µX 2µ^X

Figure 3.3: One class and Erlang-2 distributed service rate queuing model The size for a Erlang-2 model becomes (c + 1)n states. For the Erlang-n model the model would have the same structure, but with n process instead of 2 for service.

3.2.2 2 Classes

For two type the model becomes very complicated. Even if with help of the reduced model for two classes, the model becomes very big.

0

X2

Y2

X1

Y1

X2X2

X2Y2

Y2Y2

X2X1 X1X1

Y2Y1 Y1Y1

X1Y2

X2Y1

X1Y1

X2X2

+1

X2Y2

+1

Y2Y2

+1

X2X1

+1

X1X1

+1

Y2Y1

+1

Y1Y1

+1 X1Y2

+1 X2Y1

+1

X1Y1

+1

µ^X

µ_Y

µX

µY 2µ^X

µY

µ^X

2µY

µX

µ^X

2µX

µY µY

2µY µX

µY µY

µ^X

µ^Y µX

p^Y2µ^X pX2µX

p_X µ_Y

pXµX pYµY

p^Yµ^X pX2µY

pY2µY

µX

p^Yµ^X pXµX

2µY

pYµY pXµY

µY 2µY

pX µ_Y

pYµY µ_X

p^Yµ^X p^Xµ^X

µ^Y

µ^Y µX

Figure 3.4: Two class and Erlang-2 distributed service rate queuing model The model is ∼ C(c+3, c)n for two types and for k types it is ∼ C(c+2k −1, c)n.

3.3 Size differences between models

The size difference between the models are important because of computational speed. The table below shows the comparison between the models for the num- ber of states during full service, i.e. the different loading combinations.

Loading distribution No. Classes No. states

M One class 1

M Two classes c+1

E2 One class c+1

E2 Two classes C(c+3,c)

Table 3.1: Size comparison between queuing models

The vastly expanding number of states for two classes and Erlang-2 distributed service rate are demonstrated here. For bigger loading docks the number of states becomes to substantial for the calculations to be applicable. A system with 10 servicing stations would require 286 states for every single truck. Given the fact matrix multiplication has the computational complexity of O(n³) (for two nxn matrices), the matrix size is very important for computational speed.

Chapter 4

Optimization

During a day, most queuing system have peak periods with respect to arrival intensity. Therefore, a manufacturer would ideally have more resources during these periods than others.

The only model which can handle changing server amount is one class with exponential distributed service since the number of states does not depend on the amount of servers. Otherwise it creates a problem as the states cannot correlate between the models when changing capacity. This is demonstrated in table ?? as the models number of states are depended for server capacity except for one class with exponential distributed service.

If one were to optimize the capacity for queue by controlling the amount of servers for every time step, the problem would become very big. Therefore the capacity will have to be set for several time steps as in blocks. A time block would have the same server capacity for a period of time steps e.g having 3 loading dock open between 9:00-10:00.

The marginal allocation algorithm described in the background by Svanberg[4]

0+0 1+0 2+0 3+0 3+1 3+2

0+0 1+0 2+0 2+1 2+2 2+3

Figure 4.1: Demonstration of switching states from 3 servers to 2 servers.

Left number denote trucks in service, right number denote trucks in queue.

does not work for this problem since the objective functions are non-separable, i.e, are not on the form of (2.17) and (2.18).Therefore, it is impossible to create the table in step 1 of section 2.5.

As Rolfe[6] conjectured the expected queuing time in a M/M/c queue is convex to c, Dyer and Proll[12] provided the evidence for the queuing system being convex regarding the number of servers. Since every objective function for each time period is convex regard to numbers of servers, the sum of them could also be conceived to be convex. Due to the convexity of the queue’s waiting time, a modification of the marginal allocation as greedy algorithm for finding the opti- mal solution is used. In other words, the modification turns the algorithm to a steepest algorithm. Murota’s[10] paper investigates the complexity of steepest descent algorithms for discrete convex functions. The goal of the algorithm is to minimize the number of trucks in the system for time blocks while creating a convex hull.

Proposed modified algorithm structure:

1. Create at starting set of servers for the equally long time block and appoint a max capacity for how many servers can work simultaneously.

2. Calculate how the average number of trucks would decrease if you would add one server for each block lower which is below max capacity.

3. Add a server to the time block where it would reduce the average number of costumer in the system the most.

4. Stop if average number of trucks is satisfactory or every time block is at its max capacity, else go back to step 2

Chapter 5

Results and analysis

The examples in this section have the same arrival rates variations for all sim- ulations, which is illustrated in Figure 5.1

Figure 5.1: Arrival rates for the simulations

Where λ1is the arrival rate of the first class, while λ2is of the second. If there is only one class, the arrival rate is noted just as λ, the sum of the two classes,

λ = λ₁+ λ₂. (5.1)

The same notation system holds for the service rate and if not other specified, it will have the following values:

c = 2 µ = 0.04 µ₁ = 0.08 µ₂ = 0.04

Table 5.1: Values for the simulations

The focus of this study analysis is on the expected number of trucks in the system and if it is possible to mimic the chain.

5.1 Difference between Exponential and Erlang- 2 distributed service rate for one class

This section is focused on the difference between the difference between exponen- tial and Erlang-2 distributed service rates for queuing systems. For the models described in the previous chapter, the simulation for M/M/c and M/E₂/c is shown in the figure below.

Figure 5.2: Difference between exponential and Erlang-2 one class There is a significant difference for the expected number of trucks in the system between the exponential and Erlang-2 model. However they follow the same form, so the difference is mostly linearly scaled.

To find out if it is possible to replicate the behavior of the Erlang model with a exponential model, a derivation of queuing theory formulas for the waiting time is made. The derivation only works for one server models since the complexity is too high for multiple servers.

The goal of the set-up was to find the scalar K for µ_E = Kµ, where µ_E is the service rate used for a exponential model to replicate a Erlang model. The derivations are found in Appendix B.1. The scalar was derived as,

K = 1 6

s 9(λ

µ)²− 48λ

µ+ 48 + 3λ µ

!

. (5.2)

For the scaled exponential model with one server and service rate of 0.08 the result is shown below.

Figure 5.3: Scaling exponential to replicate correct Erlang-2 for one server However this only works for one server. As seen by the graph below for three servers and a service rate of 0.04.

Figure 5.4: Scaling exponential to replicate correct Erlang-2 for three servers Since it is possible for a exponential model to replicate the Erlang-2 model with a scaled factor, an examination for multiple servers is made.

The problem is that the expected waiting time for an M/M/c model is a com- plicated formula to derive and for a M/E2/c model it is even worse. It is not feasible to derive the factor the same way for multiple servers as it is for one server. A way around is to do a regression for c and λ/µ.

The scaling is instead calculated for a steady state queue (constant arrival and service rate). The scaling factor between exponential and Erlang-2 for different c, λ and µ was calculated with the help of Newton’s method (See appendix A.1).

A 2-polynomial regression between c and ^λ_µ gave the following scaling equation,

K = 1.079−0.003905c−0.02218λ

µ−0.01816c²+0.05104cλ

µ−0.03211(λ

µ)². (5.3) Scaling with the regression gave a very good fit as shown below

Figure 5.5: Scaling with linear regression exponential to replicate correct Erlang-2 for three servers

The scaling produces a result very close to the original Erlang-2 model. Even though it is a numerical scaling and not a derived, one may conclude it is possible to replicate a Erlang-2 chain with a exponential model for one class.

5.2 Difference between one and two classes for exponential distributed service rate

When using a one class exponential model to replicate a two classes exponential model, the arrival and service rates for the two classes has to be combined. For the arrival rate it is uncomplicated since it is just the sum of the rates of the two classes.

λ = λ1+ λ2 (5.4)

For service rate the same derivation is not feasible since it would be dependent to the arrival rate of each class. Theoretically the service rate should be scaled depending on the arrival rate of the classes. If there is a time period with only arrival of one class λ1, the service rate should be µ1 and vice versa for λ2 and µ2. If one would assume the scaling of the service rate is done directly between the relation of service and arrival rate for the two classes, the result would be,

µ = µ₁λ₁+ µ₂λ₂ λ1+ λ2

. (5.5)

Figure 5.6: Using one class model with assumed relation scaling to replicate correct two class model

The assumption of direct relation is clearly not a good estimate as the difference between them is considerable large as seen in Figure 5.6. In appendix B.2 a derivation for how the service rate should be is made. This was done by comparing the systems for a steady state and specifying the corresponding states as the same, resulting in the equation below.

µ = µ₁µ₂(λ₁+ λ₂) µ1λ2+ µ2λ1

(5.6)

With the derived service rate, the result is shown in figure 5.7.

Figure 5.7: Using one class model with derived relation scaling to replicate two correct class model

The derived service rate is a much better fit than the previous assumed. However it should be noted the problem the derived service rate can have with highly alternating arrival rate. The service rate only depends on the current arrival rate and not the previous, when in reality the current trucks affect the service rate.

5.3 Difference between one class exponential dis- tributed service rate and two classes Erlang- 2 distributed service rate

Since both scaling between one class to two classes and exponential to Erlang-2 exists separately, an attempt is made to combine them to see if it is possible to use a one class Markov chain for a two classes Erlang-2 model. The combination of them gives,

µ = Kµ1µ2(λ1+ λ2) µ1λ2+ µ2λ1

. (5.7)

The comparisons four using the weighted and derived scalar are illustrated be- low.

Figure 5.8: Using both scaling factors

The fit is close to the more complicated two classes Erlang-2 disturbed service rate model when using both the weighted and derived scalar.

5.4 Optimization

As described earlier, the only working model for the optimization algorithm is the queue model with one class and exponential distributed service rate.

The setup was with 8 equally long time blocks for server capacity (as described in chapter 4) and a limitation of a maximum of 5 simultaneously servers at a given time. The result of this with the regular arrival and service rate.

n c1 c2 c3 c4 c5 c6 c7 c8

1 0 1 0 0 0 0 0 0

2 0 1 1 0 0 0 0 0

3 0 1 1 1 0 0 0 0

4 0 1 2 1 0 0 0 0

5 0 1 2 1 1 0 0 0

6 0 1 2 2 1 0 0 0

7 0 1 2 2 1 1 0 0

8 0 2 2 2 1 1 0 0

9 0 2 2 2 1 2 0 0

10 0 2 2 2 1 2 1 0

11 0 2 2 2 1 2 1 1

12 0 2 2 2 1 2 2 1

13 0 2 2 2 1 2 2 2

14 0 2 2 2 2 2 2 2

15 0 2 2 2 2 2 2 3

16 0 2 2 2 2 2 3 3

17 0 3 2 2 2 2 3 3

18 0 3 2 2 3 2 3 3

19 0 3 2 2 3 2 3 4

20 1 3 2 2 3 2 3 4

21 1 3 2 3 3 2 3 4

22 1 3 2 3 3 3 3 4

23 1 3 3 3 3 3 3 4

24 2 3 3 3 3 3 3 4

25 2 3 3 3 3 3 4 4

Table 5.2: Result of optimization where n is the indexed iteration

The algorithm appears to prioritize put servers at the first blocks. This is logi- cal as the objective function is trucks in the system over day. The blocks only affect the trucks in the system during their period and for the blocks after.

The first block then affects the trucks in the system for all the block while the last block only affects itself. This phenomenon is also brought up by Massey[2]

by measuring the time lag between the peak arrivals and peak load for a system.

Given there is enough trucks in the queue, it is better to place servers dur- ing the earlier blocks. In this case, the first block does not seem to have enough arrival rate as the algorithm prioritizes the later blocks. The maximum of 5 simultaneously servers were not reached for a single time block, so it’s impact for this analysis can be ignored.

The objective function’s result of algorithm is presented in the graph below.

Figure 5.9: Convex hull for optimization method

The objective function is strictly decreasing with regards to iterations. This is as since the number of expected trucks in the system is decreasing with regard to servers. It would be impossible to see an increase of trucks in the system in the case of an increase in capacity.

However the graph shows it is not possible to claim global optimality for this algorithm. The graph shows the function is not convex, but quasi-convex. The algorithm can therefore be classified as heuristic.

Nonetheless the benefits of the method is relative computationally effective for finding a heuristic solution.

Chapter 6

Conclusions

The key for modeling Erlang-2 was to have a two folded service, as if the cus- tomer had to go through two processes to be served. This modeling would also work for Erlang-n, but with n processes for every customer. With two classes a simplification was made to combat the large complexity. The simplification was to only have the states of the combinations for trucks being served and how many are trucks queuing. However the model for one class exponential distributed service rate is the only one with the capability to change capacity as the number of states is constant with regard to amount of servers.

Because of this, the attempt was tried to create a scaling factor for the ser- vice rate to replicate the other models. The scaling between one and two classes for exponential distributed service rate was successful with a derivation. Be- tween exponential and Erlang-2 distributed service rate a derivation was also possible in the case of one server, but for multiply servers the scaling factor had to be calculated with the help of linear regression to be successful. The multiplication of these factors also made it possible to replicate the two class Erlang-2 distribute queue with the simpler one class exponential queue.

The scaling factors made it possible to replicate the other models with the one class exponential distributed service rate model and therefore optimize them.

The goal of the optimization was to minimize the number of trucks in the system by controlling the amount of servers. The method chosen was to iterate by al- ways adding servers were they would minimize the objective function the most, aka a steepest descent algorithm. It forms a computational effective method, but a heuristic as global optimality could not be guaranteed from test-data.

Appendix A

Results

A.1 Scaling Markov Erlang multiple servers

The scaling factor K was calculated with the help of Newton’s method by com- paring steady-state between the models. From this table the scaling factor was calculated with linear regression

c λ µ K 1 0.03 0.04 1.0579 1 0.02 0.04 1.0714 1 0.02 0.03 1.0667 2 0.04 0.03 1.0466 2 0.05 0.03 1.0381 2 0.05 0.04 1.0456 3 0.05 0.02 1.0335 3 0.04 0.02 1.0348 3 0.04 0.03 1.0166 4 0.04 0.015 1.0269 4 0.05 0.015 1.0298 4 0.06 0.02 1.0311 1 0.06 0.08 1.0577 1 0.06 0.09 1.0667 1 0.07 0.09 1.0537 2 0.07 0.04 1.0325 2 0.06 0.04 1.0451 2 0.08 0.05 1.0415 3 0.08 0.03 1.0277 3 0.08 0.04 1.0348 3 0.07 0.03 1.0365 4 0.08 0.03 1.0269 4 0.08 0.04 1.0133 4 0.06 0.02 1.0311 1 0.08 0.12 1.0667 2 0.08 0.08 1.0467 10 0.09 0.01 1.0169 5 0.09 0.03 1.0158 Table A.1: Data for linear regression

Appendix B

Derivations

B.1 Scaling factor for Markov-Erlang one server

The goal of the derivation is to find the scaling factor K, which is defined as

µscaled= Kµm. (B.1)

The expected costumers in line for Markov model is, E(Wm) = ρ/µscaled

1 − ρ , (B.2)

where ρ is

ρ = λ

µscaled

< 1. (B.3)

For Erlang-2 the expected costumers in line is, E(We) = ρ

1 − ρ ·r + 1 2 · 1

µ_e, (B.4)

where r is the shape (in this case 2) and ρ is, ρ = λ · r

µe

< 1. (B.5)

The service rate for Erlang-2 has to be twice as high as for the Markov model for them to have the same expected value.

µe= 2µm (B.6)

Finally the last condition is the waiting line should be equal for both the Markov and Erlang model.

E(W_e) = E(W_m) (B.7)

The scale factor can then be solved as K = 1

6 s

9( λ µm

)²− 48 λ µm

+ 48 + 3 λ µm

!

. (B.8)

B.2 One class-Two classes

The purpose of this derivation is to find a service rate for a one class system which is comparable to a two class system. This is done by comparing the systems for a steady state and specifying the corresponding states as the same.

The models represent a M/M/1/1 queue model, the left for two classes and the right for one class. A and D represent the empty states

A

B C E

D

λBA

λCA µBB

µCC

(λB+ λC)D µE

Figure B.1: M/M/1/1 queuing models. Left represents two classes and right represents one class

The probability levels of the system must be equal to one

A + B + C = 1 (B.9)

D + E = 1 (B.10)

For steady state the flow in must equal the flow out

λBA = µBB (B.11)

λCA = µCC (B.12)

(λ_B+ λ_C)D = µE (B.13)

The system’s corresponding states has to equal in the way of they have same probability

A = D (B.14)

With six equations and six unknown variables (the states and µ), it possible to solve µ.

µ = µ_Bµ_C(λB+ λC) µBλC+ µCλB

(B.15) For the general problem the service rate would be.

µ = µ₁µ₂(λ₁+ λ₂)

µ₁λ₂+ µ2λ₁ (B.16)

Appendix C

Bibliography

Bibliography

[1] Tse, K. (2014) Some Applications of the Poisson Process. Applied Mathe- matics, 5, 3011-3017. doi: 10.4236/am.2014.519288.

[2] Massey, W.A. Telecommunication Systems (2002) 21: 173.

doi:10.1023/A:1020990313587

[3] Yosuke Kawai, Hideaki Takagi, Fluid approximation analysis of a call cen- ter model with time-varying arrivals and after-call work, Operations Re- search Perspectives, Volume 2, December 2015, Pages 81-96, ISSN 2214- 7160, http://dx.doi.org/10.1016/j.orp.2015.03.003.

[4] Krister Svanberg, On marginal allocation, Optimization and Systems The- ory, KTH Stockholm Sweden.

[5] G. De Moine, Measuring and evaluating port performance and productivity, UNCTAD Monographs on Port Management, Antwerp Port Engineering and Consulting.

[6] Rolfe, Alan J. “A Note on Marginal Allocation in Multiple-Server Service Systems.” Management Science, vol. 17, no. 9, 1971, pp. 656–658. JSTOR, JSTOR, www.jstor.org/stable/2629050.

[7] Shanthikumar, J. George, and David D. Yao. “Optimal Server Allocation in a System of Multi-Server Stations.” Management Science, vol. 33, no. 9, 1987, pp. 1173–1180. JSTOR, JSTOR, www.jstor.org/stable/2631882.

[8] Al Hanbali, A., Alvarez, E., and van der Heijden, M. C. (2013). Approxima- tions for the waiting time distribution in an M/G/c priority queue. (Beta working paper series; Vol. 411, No. 411). Eindhoven: BETA Research School for Operations Management and Logistics.

[9] Ümit Yüceer, Discrete convexity: convexity for functions defined on discrete spaces, In Discrete Applied Mathematics, Volume 119, Issue 3, 2002, Pages 297-304, ISSN 0166-218X, https://doi.org/10.1016/S0166-218X(01)00191- 3. (http://www.sciencedirect.com/science/article/pii/S0166218X01001913) Keywords: Discrete convexity; Strong discrete convexity; Matrix of second forward differences

[10] Kazuo Murota,On Steepest Descent Algorithms for Discrete Convex Func- tions, SIAM Journal on Optimization 2004 14:3, 699-707

[11] Grassmann, W. “The Convexity of the Mean Queue Size of the M/M/c Queue with Respect to the Traffic Intensity.” Journal of Ap- plied Probability, vol. 20, no. 4, 1983, pp. 916–919. JSTOR, JSTOR, www.jstor.org/stable/3213605.

[12] Dyer, M. E., and L. G. Proll. “On the Validity of Marginal Analysis for Allocating Servers in M/M/c Queues.” Management Science, vol. 23, no. 9, 1977, pp. 1019–1022. JSTOR, JSTOR, www.jstor.org/stable/2630786.s

TRITA -SCI-GRU 2018:053