Mathematical Models of Multiserver Queuing System for Dynamic Performance Evaluation in Port

We discuss dynamic system performance evaluation in the river port utilizing queuing models with batch arrivals. The general models of the system are developed. This system is modelled by M/M/n/m queue with finite waiting areas and identical and independent cargo-handling capacities. The models are considered with whole and part batch acceptance or whole and part batch rejections and the interarrival and service times are exponentially distributed. Results related to the batch blocking probability and the blocking probability of an arbitrary vessel in nonstationary and stationary states have been obtained. Numerical results and computational experiments are reported to evaluate the efficiency of the models for the real system.


Introduction
Batch arrival queues with a finite waiting areas or finite-buffer space have wide range of applications in computer networks, telecommunications, transportation, manufacturing, banks, management and logistics systems, and so forth.Many results in queuing theory have been obtained by considering models, where customers arrive one by one and are served individually.However, in numerous real-world situations such as previous; mentioned various practical areas, it is frequently observed that the customers arrive in groups.Consequently, their operation processes can be adequately modelled by batch arrival queue.
Models of this type have also applications in modelling the port systems, showing that initial Markovian models are very accurate in determining a port performance.At seaport

Background in Modelling and Port Terminal Operations
This paper analyses the M X /M/n/m queue, where vessels arrive in batches of random size.Unlike the other batch arrival queue with finite waiting areas anchorage , batches which upon arrival find not enough area at the anchorage are either fully or partially rejected.Some queuing protocols are based on the whole and part batch acceptance models, and it is also known as the total or part batch rejection policy.The stochastic characteristics of the port operation are as follows: time of arrival of single vessel or in batch barge tows in the port cannot be precisely given; the service time is a random variable depending on handling capacities of berths, carriage of barges, the size of an arriving group, and so forth, and the berths are not always occupied in some periods there are no barges-the capacities are underutilized, and there are the time intervals of high utilization when the queues are formed .
The port operation and vessel movement include the following: waiting in the anchorage areas if all anchorages are occupied, the vessels or the barge tows are rejected, or if the vessels or the barge tows are larger in size than the number of available free waiting areas fills the free positions and the remaining vessels of the group are rejected and have to go to another waiting area , vessels move from anchorage to berth, loading or unloading at the berth, and towing of barges after the loading/unloading to the anchorage area or leaving the port.This cycle is called the turnaround time for the vessels in the port.
The ports or, more precisely, the anchorage-ship-berth ASB links are considered as queuing models with batch arrivals, single service, and limited queues at an anchorage.In the port, the ASB link is assumed as follows: the applied queuing system is nonstationary and stationary, with finite waiting area at anchorage; the sources of arriving pattern are not integral parts of ASB link; the service channels are berths with similar or identical and independent cargo-handling capacities; the units arrivals can be single vessels and groups as barge tows; all barge tows and single vessels at anchorage are waiting to be serviced; the service time is a continuous random variable; the size of an arriving group is a random variable; the queue length or the number of waiting anchorage areas are finite and given; interarrival times, the batch sizes, and service times are mutually independent.
The past 30 years show rising interest in the research of port and terminals systems modelling as well as their subsystems involving port equipment.Earlier research related to a bulk port, particularly to the ASB link modelling, using simulation and queuing theory, is summarized in 2-5 .Problems in modelling and simulation in the field of port equipment are shown for instance in 6, 7 .
There are numerous articles dealing with infinite and finite batch arrivals and service queues.For example, in 8 the sojourn time distribution, loss probability of an arbitrary flow, and arbitrary admitted flow in a multiserver loss queue with a flow arrival of customers is analyzed.In the case of the waiting line as a whole, which has a limited capacity, specified by a fixed maximum number of customers or finite-buffer space, queues have been analysed by several researchers in the past.Apart from the classical references 9-12 used for describing the models here, it was necessary to review the following articles.In 13 , the ergodic queue length distribution of a batch service system with finite waiting space by the method of the embedded Markov chain is discussed.The analysis of the M X /G Y /1/N queue is given in 14 .In 15 , the M X /G Y /1/K queue by Cohen's methods is investigated.The M X /M/m/s queue has been analysed in 16 , where several results related to blocking probabilities, the distribution of the number of customers in the system, the cumulative distribution function of the waiting time in the queue, and some numerical results for the single-server system are presented.In 17 , two different methods to study the behavior of a finite queue with batch Poisson inputs, and synchronous server in a computer network are proposed.The analysis M n /G n /1/K queue is given in 18 , while in 19 , the M/G/1/K queue with push-out scheme is analyzed.In 20 , the performance analysis of a discrete-time finite-buffer queue with batch input, general interarrival, and geometric service times GI X /Geom/1/N is presented.In 21 , the GI X /M/c/N queue through a combination of the supplementary variable and the embedded Markov chain techniques is proposed and analyzed.The M X /G Y /1/K B queue with a finite-buffer batch-arrival and batch service queue with variable server capacity has been considered and analysed in 22 .In 23 , a performance analysis of finite-buffer batch-arrival and batch-service Geom X /G Y /1/K B queue is presented.The analysis GI X /MSP/1/N queue is given in 24 .In 25 , the GI/BMSP/1/N queue with a finite-buffer single-server queue with renewal input, where the service is provided in batches of random size according to batch Markovian service process is analyzed.The batch arrival batch service M X /G Y /1/N queue with finite buffer under server's vacation is analysed in 26 .In 27 , the blocking probability in a finitebuffer queue whose arrival process is given by the batch Markovian arrival process was investigated.
One can conclude that the ASB link, as a main port link, has been adequately analyzed and modelled by using different modelling approaches.Various operations research models Mathematical Problems in Engineering and methods in the field of optimizing ASB link modelling are applied more and more in world terminals.In this paper, two models are described by the nonstationary and stationary, multiserver queuing system with finite waiting areas and batch arrival vessels or barge tows into the port system.

Mathematical Models Development
We consider a finite waiting areas batch arrival multiserver queue with identical and independent cargo-handling capacities: M X /M/n/m, where m denotes the number of waiting areas.The vessels or the barge tows arrive at the terminal according to a timehomogeneous Poisson process with mean arrival rate λ.The port terminal has the n berths for the service.The n berths have independent, exponentially distributed service times with common average service time 1/μ the mean cargo-handling rate per berth is μ .The queue discipline is first come first served by tows batch and random within the tow batch.Apart from the possible arrival of n vessels for service, there are m spaces in the waiting queue.The number of vessels X that arrive for service at the same time is a random variable with distribution given by a k P X k , k ≥ 1 whereas k number of vessels in group and mean E X a.The interarrival times, the batch sizes, and service times are mutually independent.The maximum number of vessels allowed in the system at any time is n m.The service times of service batches vessels are independent of the arrival process and the number served.The traffic intensity of the system is θ λa/nμ.
For convenience, three different forms for the distribution of X are considered: uniform distribution with a k 1/ n m , 1 ≤ k ≤ n m, geometric distribution with a k 1 − a a k−1 , k 1, 2, . .., and 0 < a < 1, and shifted Poisson distribution with a k a k−1 / k − 1 !e −a , k 1, 2, . .., and a > 0. Since the waiting area anchorage is finite, the following two batch acceptance models are defined.
I If the group of vessels arriving into the system finds s vessels there, then in the case that k ≤ n m − s, the whole group will be accepted into the system.Otherwise, that is, if k > n m − s, the group is totally rejected.This is called the whole batch acceptance model related to Figure 1.II In the case when group of k vessels arriving into the system finds s vessels there, then if k ≤ n m − s, the whole group will be accepted.If k > n m − s, the system accepts any n m −s vessels from the group, while the rest of k − n m −s vessels will be rejected.This is the part batch acceptance model related to Figure 2.
The state probabilities of the systems for Model I should be determined by the number of vessels within the system, shown in graph in Figure 1, and then, the following system of differential equations for whole batch acceptance model at the moment t are given by

3.1
Analogously, the state probabilities of the systems for Model II should be determined by the number of vessels within the system, shown in graph in Figure 2, and hence, the system of differential equations for part batch acceptance model at t is as follows: . . .

3.2
One of the complexities in the analysis of the models is nonstationary state of work of the systems.The analysis becomes further complicated in the case of multiserver queues as a port systems.The queuing systems for the models are Markovian, and the state probabilities are described by a set of differential equations in the steady state.The first requirement is to determine the state probabilities p k t , k 0, 1, . . ., n m.
Consider the system 3.1 of n m 1 first-order linear differential equations related to the "uniform" random variable X, that is, when a 1 a 2 • • • a m n a the values a i , with i > m n, must be different of a and 0 .By using the normalized condition n m k 0 p k t 1, the last equation of 3.1 becomes The solution of this differential equation is Further by induction based on the recursive method, it can be proved 1 that the solution of system 3.1 is given in the form where C l k are suitable real constants.By substituting the expressions for p i t , 0 ≤ i ≤ m n, given by 3.5 in system 3.1 with a 1 a 2 • • • a m n a, and using the initial conditions p 0 0 1, p i 0 0 for i 1, 2, . . ., m n, we can obtain the recurrence relations between the coefficients C for all s, 0 ≤ s ≤ m are in fact solutions of the system of linear equations which is stationary analogue to the system 3.1 see in detail 1 .Observe that any functions p i t , 0 ≤ i ≤ m n can be written as a linear combination of finite number of exponential functions of the form e −ct , with C > 0.
Unfortunately, the system 3.2 related to the uniform random variable X i.e., a 1 a 2 • • • a m n a cannot be solved using recursive method.Further, if the random variable X has the geometric probability distribution function, with a k 1−a a k−1 , where k 1, 2, . .., and 0 < a < 1, or shifted Poisson distribution with a k a k−1 / k − 1 !e −a , k 1, 2, . .., and a > 0, the systems 3.1 and 3.2 also cannot be solved using recursive method as done in the previous case related to the mentioned random variable X.Hence, MATLAB program is used for solving the corresponding systems 3.1 and 3.2 .

Blocking Probability for Both Models
In order to calculate the blocking probability at the moment t for any considered models in transient nonstationary regime, it is necessary to determine the solutions p k t , k 0, 1, . . ., n m, of the corresponding system of differential equations 3.1 or 3.2 .As previously noticed, it is outlined that closed-form expressions for time-dependent state probabilities for whole batch acceptance model of multiserver queuing system related to the uniform random variable X are obtained in 1 .However, the considered systems 3.1 and 3.2 cannot be solved in the closed-form expressions with respect to other batch size distributions geometric and shifted Poisson distribution for both models and the uniform distribution for Model II .
The batch blocking probability at the moment t is the same for both models, and it is given by the following.

8 Mathematical Problems in Engineering
Proof.Given the batch finds exactly k, 0 ≤ k ≤ n m, vessels in the system the batch will be blocked at arbitrary moment t if this batch is bigger than n m − k in size.Since P X > n m − k ∞ i n m−k 1 p i t , then by the total probability formula, we immediately obtain 3.6 .
Similarly, it can be shown see 16, 20 related to stationary blocking probability the following result.
Proposition 3.2.The blocking probability of an arbitrary vessel at the moment t is given a for Model (I): b for Model (II): where a is the mean batch size.
In the stationary state of work of the system as t → ∞ , using the systems of differential equations 3.1 and 3.2 , the corresponding set of Chapman-Kolmogorov equations can be easily obtained.Note that as t → ∞, the stationary blocking probabilities for Models I and II are p k lim t → ∞ p k t , k 0, 1, . .., where p k , k 0, 1, . . ., n m are solutions of the system 3.1 or 3.2 in the stationary state of work of the system, respectively.Hence, substituting this into 3.6 , 3.7 , and 3.8 , one obtains the corresponding formulae for stationary blocking probabilities as follows:

3.9
For example, in the case of uniform distribution X with a 1 a 2 • • • a m n 1/ n m and hence a i 0 for each i > n m, the first formula of 3.9 becomes P B 1/ n m n m k 1 kp k N ws / n m , where N ws is the average number of vessels in port.
On the other hand, in the case of geometric distribution X with a k 1 − a a k−1 , k 1, 2, . .., and 0 < a < 1, after routine calculation the first formula of 3.9 yields P B n m k 0 a n m−k p k .For example, in the case of uniform distribution X with a 1 a 2 • • • a m n 1/ n m , and hence, a i 0 for each i > n m, the second and third formula of 3.9 become On the other hand, in the case of geometric distribution X with a k 1 − a a k−1 , k 1, 2, . .., and 0 < a < 1, after routine calculations the second and third formula of 3.9 imply It is interesting to note that P B P 2 V for the Model II .

Port Performance Measures
Performance measures are the means to analyse the efficiency of the port queuing system under consideration.When state probabilities are known, performance measures can be easily obtained.In this section, some performance measures are explained, such as the probability that all the berths are occupied-P ob k n 1 p k , probability that berth is busy-P bus 1−p 0 , probability of existing vessels in queue vessels at anchor -P eq m k 1 p n k , average number of vessels at anchor-N w m k 1 k • p n k , average time that vessels spends at anchor-t w N w /λ * , where λ * is the effective arrival rate of vessel 9 , given by λ Recall that the all previous formulae except those on t w and λ * are valid in the same forms for nonstationary state with p k t instead of p k .

Numerical Results and Discussion
The efficiency of operations and processes on anchorage-ship-berth ASB link has been analyzed through the basic operating parameters such as batch blocking probability, blocking probability of an arbitrary vessel, expected number of occupied berths, probability that all the berths are occupied, probability of existence of vessels in queue vessels at anchor , expected number of vessels at anchor, and expected time at anchor.The basic characteristics of the system with the batch arrival of units and the limited waiting queue are shown in previous sections.To demonstrate the applicability of the models presented in this paper, a variety of numerical results have been showed for a combination of various parameters and various performance measures are given for the multiple server port system.All the calculations on probabilities and means were done by MATLAB program, but the results are presented after rounding up after the fourth decimal point.
Recall that both systems 3.1 and 3.2 are first order systems of m n 1 linear homogeneous differential equations with constant coefficients, and so, they may be expressed in matrix notation as dY/dt AY , where Y t p 0 t , p 1 t , . . ., p n m t is a vectorvalued function and A is a square matrix with constant coefficients .Moreover, if λ 1 is an eigenvalue for A i.e., det A − λ 1 I 0 with associated eigenvector a solution of the considered system.We used MATLAB to compute the eigenvalues and eigenvectors of a given matrix A and therefore to calculate the solutions of the corresponding systems.
The first step is to enter the given matrix A: this is done by enclosing in square brackets the rows of A, separated by semicolons.For determining eigenvalues and eigenvectors of A, we enter V, D eig A in order to get two matrices: the matrix V has unit length eigenvectors of A as column vectors, and D is a diagonal matrix with the eigenvalues of A on the diagonal.In our computations it is used MATLAB's numerical procedure "ode45" which is a souped up Runge-Kutta method and it is applied using the syntax t, y ode45 "diffeqn", t0, tf , y0 , 4.1 where t0 is the initial time, tf is the final time and y0 is the initial condition y t0 y0 .

Mathematical Problems in Engineering
In Figures 3 a and 3 b , the parameters on the ASB link of multiple server systems P B , P  In addition, Table 1 gives the numerical results of parameters on the ASB link of multiple-server systems P B , P 1 V , and P 2 V , N w and t w depending on θ 0-1 for the M X /M/2/8 queue with a 4. All of the calculations have been done in MATLAB program.
V , P (2) V a P B , P the M X /M/2/12 queue in Figure 5 a .As can be seen in Figure 5 a , blocking probabilities asymptotically approach to zero when the number of waiting areas increase m 30 .

Experimental Study
This subsection discusses how the developed models can be applied to solve some related problems in river port.In particular, we look at numerical example from 1 in which Model I has been considered to analyze the system performance of bulk cargo terminal in Smederevo city at Danube river in Serbia in the context of nonstationary work order of the ASB link in the port for a long period of time.Because the given models have been developed for analysis of the real system new bulk cargo terminal in Smederevo has been redesigned in 2005 and two new gantries were installed to unload the ore crude materials from barges into trucks , the nonstationary and stationary models, as well as its influence on the characteristics of the system, will be evaluated for the real states of the unloading terminal.The current capacity of the bulk cargo terminal is 2.4 million tons.The plan is to enlarge capacity of the bulk cargo terminal in order to satisfy the increased needs of production for ore crude materials.The purpose of giving models is to show the power of this methodology especially with respect to the various form of number of vessels in a batch distribution.The intent is to show how these models can be used to port performance evaluation.
Duration times of the nonstationary and stationary regimes were defined for the following parameters of the system: the number of berths n 3, the number of areas in the waiting queue anchorage size m 36, and the number of units in a group is given as As m is increased above 32 each model for the same batch size distribution becomes constant with further increase of the anchorage size.
The results presented here support the argument that port operating parameters under the nonstationary working regime have been obtained in the same range as results from stationary state of the port system.Various curves from all figures as a function of constant values of traffic intensity θ 0.67 and mean batch size a 6 for both working regime present always very close results, as can be seen in Table 2.The attained agreement of the results obtained from Table 2 has been also used for validation and verification of considered models.In accordance with that, this correspondence between results gives, in full, the validity both models to be used for port performance evaluation.

Conclusions
Unlike previous studies, this paper is based on real river port systems were presented, and a significant improvement is demonstrated in the operational performance as a result of the M X /M/n/m queue.A real problem of dynamic system performance evaluation at a port has driven the development of these models.Closed-form expressions for state probabilities and blocking probabilities for the whole batch acceptance model are presented.As the problem gets complicated, a closed-form solution may not be possible to derive for the part batch acceptance model.Contribution of this paper is twofold: analytical models development and analysis of dynamic performance measures and the blocking probability P B t for an arbitrary batch and the blocking probability P V t of an arbitrary vessel at the moment t are expressed.It is evident that the influence of the nonstationary work order on operating parameters, expressed by the function P i t , makes use of the queuing theory models to describe the operation of ASB link possible.
The results have revealed that analytical modelling is a very effective method to examine the impact of introducing priority, for certain class of vessels, on the ASB link performance and show that Model I has always lower batch blocking probability and has always higher blocking probability of an arbitrary vessel and hence lower vessel throughput.It was shown that this conclusion does depend on the different distribution of batch size too.

Figure 1 :
Figure 1: Graph of state of system for whole batch acceptance.

Figure 2 :
Figure 2: Graph of state of system for part batch acceptance.
s k .These recurrence relations give the values of each coefficient C s k .Obviously, p n−l C n−l 0 for all l, 0 ≤ l ≤ n and p n m−s C n m−s 0

2 VFigure 3 :
Figure 3: Impact of the nonstationary work order on operating parameters' changes for M X /M/2/12 queue with different distribution of batch size.

1 V and P 2 V 2 V 1 V and P 2 VFigure 4 :
Figure 4: Impact of the traffic intensity on operating parameters' changes for M X /M/2/12 queue with different distribution of batch size.

Figures 4 c
Figures 4 c and 4 dshow N w and t w as a function of θ for the M X /M/2/12 queue with a 4.In addition, Table1gives the numerical results of parameters on the ASB link of multiple-server systems P B , P

Figures 5 a 1 V , and P 2 V
Figures 5 a and 5 b express P B , P 1 V , and P 2 V as a function of anchorage size number of waiting areas for both models.These figures include a family of curves for different forms of the distribution of X.The constant parameters for these curves are θ 0.61 and a 4 for

Figure 5 :
Figure 5: Impact of waiting area size on P B and P V for the M X /M/2/12 queue with different distribution of batch size.

Figure 5 b 1 V and P 2 V
shows P B , P as a function of anchorage size θ 0.92 and a 6 for the M X /M/2/12 queue.

Figure 7 :
Figure 7: Impact of the traffic intensity on operating parameters' changes for M X /M/3/36 queue with different distribution of batch size.

MathematicalFigure 8 :
Figure 8: Impact of waiting area size on operating parameters' changes for M X /M/3/36 queue with different distribution of batch size.

Table 1 :
Impact of the traffic intensity on operating parameters' changes for M X /M/2/8 queue with different distribution of batch size.

Table 2 :
Port performance parameters for M X /M/3/36 queue nonstationary and stationary regime and for M X /M/2/32 queue stationary regime with different distributions of batch size.