Efficient Method to Approximately Solve Retrial Systems with Impatience

Abstract

We present a novel technique to solve multiserver retrial systems with impatience. Unfortunately these systems do not present an exact analytic solution, so it is mandatory to resort to approximate techniques. This novel technique does not rely on the numerical solution of the steady-state Kolmogorov equations of the Continuous Time Markov Chain as it is common for this kind of systems but it considers the system in its Markov Decision Process setting. This technique, known as value extrapolation, truncates the infinite state space using a polynomial extrapolation method to approach the states outside the truncated state space. A numerical evaluation is carried out to evaluate this technique and to compare its performance with previous techniques. The obtained results show that value extrapolation greatly outperforms the previous approaches appeared in the literature not only in terms of accuracy but also in terms of computational cost.

1. Introduction

A common assumption when evaluating the performance of communication systems is that users that do not obtain an immediate service leave the system without retrying. However, due to the increasing number of customers and network complexity, the customer behavior in general, and the retrial phenomenon in particular, may have a nonnegligible impact on the system performance. For example, in mobile cellular networks the importance of the retrial phenomenon has been stressed in [1–3]. An extensive bibliography on retrial queues can be found in [4]. The modeling of repeated attempts has been a subject of numerous investigations, because these systems have a nonhomogeneous and infinite state space. However, it is known that the classical theory [5] is developed for random walks on the semistrip {0,…,𝐶}×ℤ₊ (being 𝐶 the number of servers) with infinitesimal transitions subject to conditions of space homogeneity.

When the space-homogeneity condition does not hold, for example, in the case of retrial queues, the problem of calculating the equilibrium distribution has not been solved beyond approximate techniques when the number of servers is higher than two [6]. In particular, Marsan et al. [7] propose a well-known approximate technique for its analysis. In [8], a generalization of the approximate technique in [7] was proposed, showing a substantial improvement in the accuracy at the expense of a marginal increase of the computational cost. Those approximations are based on the reduction of an infinite state space to a finite one by aggregating states. Other solutions maintain the infinite state space but homogenize it beyond a given level in order to solve the system. These later models are known as generalized truncated models [6] and usually present the advantage of providing a much better accuracy than the finite methodologies [9]. In this category we find the models proposed by Falin [10], by Neuts and Rao [11], and by Artalejo and Pozo [6]. All these approaches rely on the numerical solution of the steady-state Kolmogorov equations of the Continuous Time Markov Chain (CTMC) that describes the system under consideration.

Very recently, however, an alternative approach for evaluating infinite state space Markov processes has been introduced by Leino et al. [12–14]. The new technique, named value extrapolation, does not rely on solving the global balance equations. This technique considers the system in its MDP (Markov Decision Process) setting and solves the expected value from the Howard equations written for a truncated state space. Instead of a simple truncation, the relative values of states just outside the truncated state space are estimated using a polynomial extrapolation based on the states inside, obtaining a closed system. Therefore, we can compute any performance parameter as far as we are capable to express it as the expected value of a random variable that is function of the system state.

So far the value extrapolation technique has been applied to multiclass single server queues showing very promising results. It must be noted that a key aspect on the application of value extrapolation lies on the election of the extrapolating function for the relative state values. Indeed, in [14] the authors show that by selecting an appropriate polynomial function the technique yields exact results for the moments of the queue length in a multiclass Discriminatory Processor-Sharing (DPS) system. Unfortunately, the appropriateness of the functional form of the extrapolation depends on the system and also on the revenue function, that is, the performance parameter we are interested in. Hence, there is no universal good choice for the extrapolating function. In this paper we address the application of the value extrapolation technique to an important class of queuing systems, for example, retrial queues, which are essentially different of the type of queues to which this technique has been applied. A potential drawback of value extrapolation compared to conventional state space truncation methods is that, since the stationary state probabilities are not obtained, if one want to compute several performance parameters, the technique has to be applied once per each of them. We apply well-known linear algebra algorithms to compute several performance parameters simultaneously, and through some series of numerical examples we show that, at least for the type of system that we are studying, the relative impact in terms of computational cost is marginal.

The application of the value extrapolation technique has only addressed problems in which relative state values are expected to follow a polynomial tendency. In this paper we develop the value extrapolation technique to solve a multiserver retrial system, addressing also the drawback of computing only a single performance parameter every time the technique is used.

In a first part of the paper, we develop the analytical part of the technique, defining the associated Howard equations of the model and the revenue functions. In a second part, we compare our technique with other previously proposed techniques in terms of accuracy and computational cost. Results show that the proposed technique clearly outperforms the rest of the studied techniques in terms of computational cost, and this improvement is even much higher in terms of accuracy.

The rest of the paper is structured as follows. Section 2 describes the system under study, while Section 3 introduces the solving technique used. In Section 4, the numerical analysis is carried out, evaluating the value extrapolation technique and comparing it with other previous solving techniques proposed in the literature. Final remarks and a summary of results are provided in Section 5.

2. System Model

The system under study is a generic retrial system including user impatience, that is, users leave the system with certain probability after a nonsuccessful retrial. A diagram of the system is shown in Figure 1. Users arrive following a Poisson process with rate 𝜆 to a system with 𝐶 servers and request an exponentially distributed service time with rate 𝜇. Without loss of generality, we consider that each user occupies one server. When a new request finds all servers occupied, it joins the retrial orbit with probability 1. After an exponentially distributed time of rate 𝜇_𝑟, this session retries. The reattempt is successful if it finds a free server. Otherwise, the user leaves the system with probability 𝑃_𝑖 or returns to the retrial orbit with probability (1−𝑃_𝑖), starting the retrial procedure again. Note that we consider an infinite capacity for the retrial orbit.

Figure 1 

Retrial model under study.

The model considered can be represented as a bidimensional CTMC, 𝑆(𝑡)={𝐾(𝑡),𝑀(𝑡)}, where 𝐾(𝑡) is the number of sessions being served and 𝑀(𝑡) the number of users in the retrial orbit at time 𝑡. The state space of the process is defined by𝒮∶=𝑠=(𝑘,𝑚)∶𝑘≤𝐶;𝑚∈ℤ₊. (2.1)

Figure 2 shows the transition diagram of such system, showing two important properties in the dimension corresponding to the number of users in the retrial orbit: on the one hand its infinite cardinality and on the other hand its space-heterogeneity produced by the fact that retrial rate depends on the number of customers in the retrial orbit.

Figure 2 

Transition diagram.

3. Solving Technique

In this section we develop the value extrapolation technique for the model presented in Section 2. Additionally, we present some particularities that should be taken into account when using this technique.

3.1. MDP Settings

As it has been aforementioned, the problem under interest has not a closed form solution when 𝐶>2 [6], so approximation techniques are mandatory. To the best of our knowledge, all the approximate techniques that have appeared in literature compute the steady state probabilities (𝜋(𝑠)) using the balance equations in order to compute the desired performance parameters, that is, solving the linear system of equations:𝜋(𝑠)_{𝑠′≠𝑠}𝑞_𝑠𝑠′=_{𝑠′≠𝑠}𝜋𝑠^′𝑞_𝑠𝑠′∀𝑠∈𝒮, (3.1) along with the normalization condition ∑_𝑠𝜋(𝑠)=1, where 𝑞_𝑠𝑠′ represents the transition rate from state 𝑠 to 𝑠′.

Notwithstanding, value extrapolation is not based on the probability of being in a certain state, but on a new metric called relative state values. Relative state values appear when we consider the system in the setting of an MDP. Formally, an MDP can be defined as a tuple {𝒮,𝒜,𝒫,ℛ}, where 𝒮 is a set of states, 𝒜 is a set of actions, 𝒫 is a state transition function, and ℛ is a revenue function. The state of the system can be controlled by choosing actions 𝑎 from 𝒜, influencing in this way the state transitions. The transition function 𝒫∶𝒮×𝒮×𝒜→ℝ₊ specifies the transition rate between a pair of states when a certain action is taken at the original state. The first characteristic of the value extrapolation technique is the necessity of the definition of a revenue function that must be a function of the system state, that is, 𝑟(𝑠). Following the definition of the revenue function for every state, in the steady state, a mean revenue rate of the entire process can be introduced as 𝑟=∑_𝑠∈𝒮𝜋(𝑠)𝑟(𝑠). In the value extrapolation technique, the revenue function ℛ has to be defined so that the resulting average revenue 𝑟 coincides with the desired performance metric.

Once we have defined the MDP framework as well as the revenue function, we are in a position to define the relative state values. It is obvious that after performing an action in state 𝑠 the system will collect a revenue for that action (𝑟(𝑠)), but as the number of transitions increases, the average revenue collected converges to 𝑟. The relative state value (𝑣(s)) is equal to the difference between the total revenue incurred when the system starts at state 𝑠 and the total revenue incurred in a system for which the revenue rate at all states is 𝑟:𝑣(𝑠)=𝐸^∞₀(𝑟(𝑆(𝑡))−𝑟)𝑑𝑡∣𝑆(0)=𝑠. (3.2)

The Howard equations relate revenues, relative state values, and transition rates:𝑟(𝑠)−𝑟+_𝑠′𝑞_𝑠𝑠′𝑣𝑠^′−𝑣(𝑠)=0∀𝑠∈𝒮. (3.3)

The Howard equations represent the policy evaluation phase of the well-known policy iteration algorithm, the most widespread dynamic programming technique, proposed in [15]. The Howard equations that correspond to the system under study are𝑟(𝑘,𝑚)−𝑟+𝜆(𝑣(𝑘+1,𝑚)−𝑣(𝑘,𝑚))+𝑘𝜇(𝑣(𝑘−1,𝑚)−𝑣(𝑘,𝑚))+𝑚𝜇_𝑟(𝑣()(𝑟()(𝑣()((𝑣()(𝑘+1,𝑚−1−𝑣𝑘,𝑚))=0if𝑘<𝐶,𝐶,𝑚−𝑟+𝜆𝐶,𝑚+1−𝑣𝐶,𝑚))+𝐶𝜇𝐶−1,𝑚−𝑣𝐶,𝑚))+𝑚𝜇_𝑟𝑃_𝑖(𝑣(𝐶,𝑚−1)−𝑣(𝐶,𝑚))=0if𝑘=𝐶. (3.4)

As we can observe the number of states is infinite because 𝑚 can take any value in ℤ₊, thus we need to truncate the state space to 𝒮. In our case, the truncated state space is defined by𝒮∶={𝑠=(𝑘,𝑚)∶𝑘≤𝐶;𝑚≤𝑄}. (3.5)

In general, 𝑄 is known as the truncation level. As we choose a higher value of 𝑄, we can expect a higher accuracy as the system is more similar to the original one, but we will have a higher computation cost too. Therefore, the objective will be to achieve a certain accuracy with the minimum value of 𝑄.

There will be as many Howard equations as number of states, |𝒮|. The number of unknowns will be the |𝒮| relative state values plus the expected revenue 𝑟, that is, |𝒮|+1 unknowns. However, as only the differences in the relative values appear in the Howard equations, we can set 𝑣(0)=0, so we will have a solvable linear system of equations with the same number of equations as unknowns.

3.2. Polynomial Fitting

The traditional truncation sets 𝑞_𝑠𝑠′=0forall𝑠^′∉𝒮, but value extrapolation performs a more efficient truncation. Basically, value extrapolation considers the relative state values outside 𝒮 that appear in the Howard equations as an extrapolation of some relative state values inside 𝒮. As we truncate the retrial orbit dimension beyond a value 𝑄, the value extrapolation technique uses the state value of some states in 𝒮 to approximate 𝑣(𝐶,𝑄+1), which is expected to improve the accuracy significantly, as it is better than ignoring these relative state values. Note that if extrapolation yielded the exact value for those states outside 𝒮, the results obtained by solving the truncated model would be exact. Also note that including value extrapolation neither increases the computational cost nor increases the number of Howard equations, which remains equal to |𝒮|=(𝐶+1)×(𝑄+1).

Summarizing, the objective of value extrapolation is to find an extrapolation function that fits with some points in 𝒮 so that it approximates also points outside 𝒮. It is important to choose a fitting function that makes the Howard equations remain a closed system of linear equations. The most common fitting functions that fulfill that condition are the polynomials. We can use all the states in 𝒮 into the fitting procedure (global fitting) or, what is most commonly used, only a subset (𝒮_𝑓) of them (local fitting).

For the sake of simplicity, in the following description we will assume there exists a mapping 𝒲 from the two-dimensional set of states into a single-dimensional set, for example, the real numbers: 𝒲∶𝒮_𝑓→ℝ. Hence, below we deal with states as if they were real values given as 𝑤=𝑊(𝑠). The specific mapping used for the model under study is specified later on.

The choice of 𝒲 will highly depend on the states we want to extrapolate its relative state value. Note also that the function 𝑓(𝑤) and the set 𝒮_𝑓 need to be chosen so that the parameters in 𝑓(𝑤) have unambiguous values, that is, in the case of choosing a polynomial as the fitting function, the number of different points in 𝒮_𝑓 has to be equal or greater than the number of coefficients in the polynomial. In general, the procedure to compute the coefficients of the fitting polynomial 𝑎_𝑖 consists in minimizing the least mean squared error𝐸=_𝑤∈𝒲𝑓(𝑤)−𝑣(𝑤)². (3.6)

Then the optimal values for the 𝑎_𝑖’s can be computed by solving the equations𝜕𝐸𝜕𝑎_𝑖=0∀𝑖. (3.7)

In our case, we are using as many points as the number of parameters of the fitting polynomial, so the fitting procedure is an ordinary polynomial interpolation and 𝐸=0, that is, all the considered points will lie in the curve of the polynomial. In this case, the problem can be formulated as follows. Given a set of 𝑛=|𝒲(𝒮_𝑓)|=|𝒮_𝑓| points (𝑤₀,𝑣(𝑤₀)),…,(𝑤_𝑛−1,𝑣(𝑤_𝑛−1)), where there are not two identical 𝑤_𝑖, we can determine an (𝑛−1)-th degree polynomial so that 𝑓(𝑤_𝑖)=𝑣(𝑤_𝑖) for 𝑖=0,…,𝑛−1, where𝑓(𝑤)=𝑎₀+𝑎₁𝑤+𝑎₂𝑤²+⋅⋅⋅+𝑎_𝑛−1𝑤^𝑛−1. (3.8)

The interpolating polynomial satisfies the following 𝑛 linear equations:𝑓𝑤_𝑖=𝑎₀+𝑎₁𝑤_𝑖+𝑎₂𝑤²_𝑖+⋅⋅⋅+𝑎_𝑛−1𝑤_𝑖^𝑛−1=𝑣𝑤_𝑖𝑖=0,…,𝑛−1, (3.9) which in a matrix form are𝐴𝑎=1𝑤₀…𝑤₀^𝑛−11𝑤₁…𝑤₁^𝑛−1⋮⋮⋱⋮1𝑤_𝑛−1…𝑤^𝑛−1_𝑛−1𝑎₀𝑎₁⋮𝑎_𝑛−1=𝑣𝑤₀𝑣𝑤₁⋮𝑣𝑤_𝑛−1=𝑏. (3.10)

The matrix of coefficients of this system (𝐴) is a Vandermonde matrix, whose determinant is nonvanishing and therefore 𝐴 is invertible. Thus, there always exists a unique solution to the considered linear system of equations or, equivalently, there exists a unique polynomial that goes through all the 𝑛 points. However, Vandermonde matrices are often badly conditioned, specially if some 𝑤_𝑖 are very close, so the procedure to compute the fitting polynomial is also badly conditioned. It is important to note that the unicity of the fitting polynomial does not mean that it cannot be written in a basis different from the standard basis. More concretely in this work we have used the Lagrange basis.

For the considered interpolation problem, the polynomial in its Lagrange setting is a linear combination𝐿(𝑤)=^𝑛−1_𝑗=0𝑣𝑤_𝑗ℓ_𝑗(𝑤) (3.11) of Lagrange basis polynomialsℓ_𝑗(𝑤)=^𝑛−1_{𝑖=0𝑖≠𝑗}𝑤−𝑤_𝑖𝑤_𝑗−𝑤_𝑖=𝑤−𝑤₀𝑤_𝑗−𝑤₀⋅⋅⋅𝑤−𝑤_𝑗−1𝑤_𝑗−𝑤_𝑗−1𝑤−𝑤_𝑗+1𝑤_𝑗−𝑤_𝑗+1⋅⋅⋅𝑤−𝑤_𝑛−1𝑤_𝑗−𝑤_𝑛−1. (3.12)

For the truncated problem of interest and as shown in Figure 3, we will have a Howard equation in which appears 𝑣(𝐶,𝑄+1), that is a state value of a state that does not belong to 𝒮. Therefore, we must approximate the value 𝑣(𝐶,𝑄+1) by using some relative state values of states belonging to 𝒮. It is important to emphasize that for the extrapolation of 𝑣(𝐶,𝑄+1) we only use states from the last row of the model shown in Figure 3, that is, states of the form 𝑠=(𝐶,𝑚), with varying 𝑚. With this choice, we define the mapping 𝒲 as 𝒲((𝐶,𝑚))=𝑚. Moreover, we use an (𝑛−1)-th degree polynomial that interpolates the 𝑛 points in 𝒮_𝑓∶={𝑠_𝑖=(𝐶,𝑄−𝑖)∣𝑖=0,…,𝑛−1} and then 𝒲(𝒮_𝑓)={𝑤_𝑖=𝑄−𝑖∣𝑖=0,…,𝑛−1}:𝑤₀=𝑄⟶𝑣𝑤₀(),𝑤=𝑣𝐶,𝑄₁=𝑄−1⟶𝑣𝑤₁(),⋮𝑤=𝑣𝐶,𝑄−1_𝑗=𝑄−𝑗⟶𝑣𝑤_𝑗,⋮𝑤=𝑣𝐶,𝑄−𝑗_𝑛−1=𝑄−(𝑛−1)⟶𝑣𝑤_𝑛−1=𝑣(𝐶,𝑄−(𝑛−1)). (3.13)

Figure 3 

Truncated model and states that appear in Howard equations outside the truncated model.

This way, the general form of the extrapolation state when using an (𝑛−1)-th degree polynomial is𝑣^(𝑛)(𝐶,𝑄+1)=𝐿^(𝑛)(𝑄+1)=^𝑛−1_𝑗=0𝑣𝐶,𝑄−𝑗ℓ_𝑗(𝑄+1). (3.14)

For example, in the case of linear extrapolation (𝑛=2), we use (𝑄,𝑣(𝐶,𝑄)) and (𝑄−1,𝑣(𝐶,𝑄−1)), having𝑣⁽²⁾(𝐶,𝑄+1)=𝐿⁽²⁾(𝑄+1)=𝑣(𝐶,𝑄)ℓ₀(𝑄+1)+𝑣(𝐶,𝑄−1)ℓ₁()()()−()𝑄+1=𝑣𝐶,𝑄𝑄+1𝑄−1()()()𝑄−𝑄−1+𝑣𝐶,𝑄−1𝑄+1−𝑄()()().𝑄−1−𝑄=2𝑣𝐶,𝑄−𝑣𝐶,𝑄−1 (3.15)

Following a similar procedure, we can obtain the next relationship for 𝑛=3 and 𝑛=4:𝑣⁽³⁾()()()(),𝑣𝐶,𝑄+1=3𝑣𝐶,𝑄−3𝑣𝐶,𝑄−1+𝑣𝐶,𝑄−2₍₄₎(𝐶,𝑄+1)=4𝑣(𝐶,𝑄)−6𝑣(𝐶,𝑄−1)+4𝑣(𝐶,𝑄−2)−𝑣(𝐶,𝑄−3). (3.16)

In general, for (𝑛−1)-th degree polynomials and using the Lagrange basis to reduce the complexity of the procedure, a simple closed-form expression for the extrapolated value can be obtained by𝑣^(𝑛)(𝐶,𝑄+1)=^𝑛−1_𝑘=0(−1)^𝑘𝑛𝑣()𝑘+1𝐶,𝑄−𝑘, (3.17) where 𝑛 is the number of coefficients taken for Lagrange polynomials.

3.3. Revenue Function

As performance parameters are not computed from the steady state probabilities as usual, it is important to explain more carefully how they are computed. By definition, 𝑟(𝑠) is the revenue rate obtained when the system is in state 𝑠. Therefore, we must define the revenue as the performance parameter we want to compute. The effect of that action is that the computed 𝑟 will be the performance parameter we are looking for. Additionally, the inputs 𝑟(𝑠) in the Howard equations must be properly set. Table 1 gives several examples on how 𝑟(𝑠) can be set in order to obtain certain performance parameters such as: blocking probability 𝑃_b=Prob{𝐾=𝐶}, mean number of users in the retrial orbit 𝑁_ret=𝐸[𝑀], nonservice probability 𝑃_ns (probability of a user leaving the system due to impatience without obtaining service), probability of being in a certain state 𝜋(𝐾,𝑀), and probability of having 𝐾 busy servers 𝐵(𝐾).

Table 1 

Revenue function definition.

As an example, we focus on the blocking probability and we define the revenue function to be 1 in those states in which an attempt is blocked, that is, when 𝑟(𝐶,𝑚)=1,for all 𝑚, and 0 in the rest of states, 𝑟(𝑘,𝑚)=0,𝑘≠𝐶, for all 𝑚.

3.4. Effect of the Value Extrapolation into the Howard Equations

In our problem, and as mentioned above, we will only have to replace 𝑣(𝐶,𝑄+1) by its approximate value in the Howard equation that corresponds to the state 𝑣(𝐶,𝑄). As an example, if we use linear extrapolation (𝑛=2), that equation becomes𝑟(𝐶,𝑄)−𝑟+𝑣(𝐶,𝑄)−𝜆−𝐶𝜇−𝑄𝑃_𝑖𝜇_𝑟+𝜆𝑣(𝐶,𝑄+1)+𝐶𝜇𝑣(𝐶−1,𝑄)+𝑄𝑃_𝑖𝜇_𝑟𝑣(𝐶,𝑄−1)=𝑟(𝐶,𝑄)−𝑟+𝑣(𝐶,𝑄)𝜆−𝐶𝜇−𝑄𝑃_𝑖𝜇_𝑟()++𝐶𝜇𝑣𝐶−1,𝑄𝑄𝑃_𝑖𝜇_𝑟−𝜆𝑣(𝐶,𝑄−1)=0. (3.18)

As 𝑣(𝐶,𝑄+1) no longer appears in the Howard equations, the linear system of equations we have consists of (𝐶+1)×(𝑄+1) equations with the same number of unknowns. This system can be expressed in matrix form for simplicity reasons. Therefore, the system can be seen as 𝐱𝐓=𝐛, where 𝐱 is a vector with the (𝐶+1)×(𝑄+1) unknowns (𝑟 and the relative state values 𝑣(𝑠)) and 𝐛 are the negative revenue rates for the different states:()()𝑣()(),()()().𝐱=𝑟𝑣0,1⋯𝑣0,𝑄1,0⋯𝑣𝐶,𝑄𝐛=−𝑟0,0−𝑟0,1⋯−𝑟𝐶,𝑄 (3.19)

Matrix 𝐓 represents the matrix of coefficients and can be constructed making all the elements in the first row of matrix 𝐓_𝟎 equal to −1.

Matrix 𝐓_𝟎 is given by𝐓_𝟎=𝐀⁰₁𝐀⁰₀𝐀⋯00₁₂𝐀¹₁⋯00⋮⋮⋱⋮⋮00⋯𝐀₁^𝐶−1𝐀₀^𝐶−100⋯𝐀^𝐶₂𝐀^𝐶₁, (3.20) where the submatrices are defined as𝐀^𝑘₀=()(),𝐀𝑘+1𝜇𝐈,for0≤𝑘≤𝐶−1_𝑘2=𝜆𝜇_𝑟0⋯000𝜆2𝜇_𝑟⋯00⋮⋮⋮⋱⋮⋮000⋯𝜆𝑄𝜇_𝑟𝐀000⋯0𝜆,for1≤𝑘≤𝐶,_𝑘1=𝛼00⋯00𝛼−𝜇_𝑟0⋯000𝛼−2𝜇_𝑟⋯0⋮⋮⋮⋱⋮000⋯𝛼−𝑄𝜇_𝑟,for0≤𝑘≤(𝐶−1),𝛼=−𝜆−𝑘𝜇. (3.21)

When 𝑘=𝐶, using linear (𝑛=2) and quadratic (𝑛=3) extrapolation, we obtain, respectively,𝐀^𝐶₁=𝛽𝑃_𝑖𝜇_𝑟⋯00𝜆𝛽−𝑃_𝑖𝜇_𝑟()𝑃⋯000𝜆⋯00⋮⋮⋱⋮⋮00⋯𝛽−𝑄−1_𝑖𝜇_𝑟𝑄𝑃_𝑖𝜇_𝑟−𝜆00⋯𝜆𝜆−𝐶𝜇−𝑄𝑃_𝑖𝜇_𝑟,𝐀^𝐶₁=𝛽𝑃_𝑖𝜇_𝑟⋯00𝜆𝛽−𝑃_𝑖𝜇_𝑟()𝑃⋯000𝜆⋯00⋮⋮⋱⋮⋮00⋯𝑄−1_𝑖𝜇_𝑟𝜆00⋯𝛽−(𝑄−1)𝑃_𝑖𝜇_𝑟𝑄𝑃_𝑖𝜇_𝑟−3𝜆00⋯𝜆2𝜆−𝐶𝜇−𝑄𝑃_𝑖𝜇_𝑟, (3.22) where 𝛽=−𝜆−𝐶𝜇.

In general, if the extrapolation is done with 𝑛≤𝑄+1 points, the matrix 𝐀^𝐶₁ is given as𝐀^𝐶₁=𝛽𝑃_𝑖𝜇_𝑟⋯0𝜆𝑐_𝑄^(𝑛)𝜆𝛽−𝑃_𝑖𝜇_𝑟⋯0𝜆𝑐^(𝑛)_𝑄−10𝜆⋯0𝜆𝑐^(𝑛)_𝑄−2()𝑃⋮⋮⋱⋮⋮00⋯𝑄−1_𝑖𝜇_𝑟𝜆𝑐₂^(𝑛)00⋯𝛽−(𝑄−1)𝑃_𝑖𝜇_𝑟𝑄𝑃_𝑖𝜇_𝑟+𝜆𝑐₁^(𝑛)00⋯𝜆−𝜆−𝐶𝜇−𝑄𝑃_𝑖𝜇_𝑟+𝜆𝑐₀^(𝑛), (3.23) where𝑐_𝑙^(𝑛)=(−1)^𝑙𝑛𝑙+1ifl<𝑛,0if𝑙≥𝑛. (3.24)

Note that the size of matrix 𝐓 does not depend on the order of the polynomial used to perform the extrapolation; only the last column in the matrix 𝐴^𝐶₁ depends on the polynomial adjustment. This characteristic has the advantage that there will not be any difference in the computation cost when using higher order of extrapolation.

The main drawback of the value extrapolation technique is that this technique is only able to compute one performance parameter each time we solve the system. Notwithstanding, we can overcome this drawback in the following way. In a general manner, the solution of the system 𝐱𝐓=𝐛 can be obtained using the inverse matrix of 𝐓 by doing 𝐱=𝐛𝐓^−𝟏. Note also that choosing a different performance parameter to solve will only affect to the values in 𝐛. Therefore, computing a second performance parameter will only increase the computation expenses by the cost of the product 𝐛𝐓^−𝟏, as the rest of the process (specifically the computation of the inverse matrix 𝐓^−𝟏) is solved only once. Similarly, we can compute several performance parameters with a marginal increase in the computation cost using LU factorization, as the first part of the procedure (the factorization, which represents the most computationally expensive part) is done only once for the 𝐓 matrix. This characteristic of the value extrapolation technique can be observed in Figure 4, where we show that the computation time (results have been obtained using Matlab running on an Intel Pentium IV 3 GHz) is only marginally increased when we compute additional performance parameters.

Figure 4 

Computation cost when solving 𝑝 performance parameters simultaneously.

4. Results

In order to evaluate and compare the proposed technique, we have studied its performance in several scenarios. Letting 𝜌=𝜆/(𝐶𝜇), we have studied different system loads by modifying 𝜆 and keeping 𝐶=50 resource units and 𝜇⁻¹=180 s. The retrial phenomenon has been configured with 𝜇_𝑟⁻¹=100 s and 𝑃_𝑖=0.2. Although only one configuration of the retrial orbit has been chosen, there will be fairly different working points, as the system load is widely modified.

For obtaining the results, we have used the relative error of different performance parameters, defined for a generic performance parameter Ψ by 𝜖_Ψ=|Ψ^approx−Ψ^exact|/Ψ^exact. In order to obtain an accurate enough estimate of Ψ which can be used as Ψ^exact, we ran all techniques with increasing and sufficiently high values of 𝑄 so that the value of Ψ had stabilized up to the 14th decimal digit. As expected all techniques converged to the same value in the performance parameters under study, Ψ∈{𝑃_b,𝑃_ns,𝑁_ret}.

4.1. Value Extrapolation Evaluation

Table 2 shows the minimum value of 𝑄 needed to obtain a relative error lower than 10⁻⁸ for different performance parameters and loads (columns) and for different orders of the extrapolation polynomials (rows). Note that VE𝑥 denotes the use of an extrapolation polynomial of order 𝑥=(𝑛−1). The number in bold indicates the lowest truncation level of all the polynomials studied. Finally, the last row of Table 2 shows the exact value of the studied performance parameter for that scenario.

Table 2 

Minimum value of 𝑄 to obtain relative errors (𝜖) lower than 10⁻⁸.

From Table 2 we conclude that there is not a clear choice in the order of the best polynomial. In general, neither the lowest nor the highest order polynomials are recommendable, so we recommend to use the intermediate cases. Furthermore, the fact that using VE𝑥 enforces us to use a model with 𝑄≥𝑥 (see Section 3.2) must be considered in the choice of the polynomial. For that reason we can conclude that, for the problem and scenario of interest and for the relative accuracy we want to achieve, VE8 represents a good tradeoff between accuracy and value of 𝑄 needed. Therefore, hereafter we will use the polynomial of order 8 (VE8) and we will simply denote it as VE.

4.2. Comparison with Other Techniques

In this section we compare the performance of value extrapolation with other techniques based on the traditional approach of solving the steady state probabilities using the balance equations for later computing the performance parameters of interest. Although other approaches exist, we have chosen the technique proposed in [8], referred to hereafter as FM, and the one proposed by Neuts and Rao in [11], referred to as NR. Note that we have not compared the results with the technique proposed by Artalejo and Pozo [6] as this last technique does not include the impatience phenomenon, so it is not directly applicable. A similar reasoning can be done for the technique proposed by Falin [10].

In Table 3 we show the minimum values of 𝑄 needed to obtain a relative error lower than 10⁻⁸ for different performance parameters and for the aforementioned techniques. Results show that value extrapolation clearly outperforms classical techniques as it needs a much lower value of 𝑄 to achieve a certain accuracy in all the scenarios under study and for all the parameters studied. Similarly, in Figures 5, 6 and 7, we plot the relative error for 𝑃_b, 𝑃_ns and 𝑁_ret, respectively, when 𝜌=0.7 and for the different techniques deployed. Results show that, for a same value of 𝑄, VE is able to obtain lower relative errors than NR and FM. The difference in the relative errors is around 4 to 5 orders of magnitude, which supposes a very clear improvement.

Table 3 

Minimum 𝑄 value to obtain relative errors (𝜖) lower than 10⁻⁸.

Figure 5 

Relative error in 𝑃b for different techniques.

Figure 6 

Relative error in 𝑃ns for different techniques.

Figure 7 

Relative error in 𝑁ret for different techniques.

4.3. Computation Cost

Although it is shown that VE clearly outperforms NR and FM techniques, it is interesting to study their associated computation cost. In Figures 8, 9, and 10, we plot the time needed to achieve a certain relative error for 𝑃_b, 𝑃_ns, and 𝑁_ret using the different techniques under study. Note also that, although it has been obtained using VE8, choosing a diferrent order for the extrapolation polynomial would not change the computation cost, as the linear system of equations to be solved remains of the same size. However, results should be interpreted carefully, because computation costs highly depend on the algorithm used to solve the resulting system of equations. More concretely, for solving the systems obtained in the FM and NR techniques, we have made use of the efficient algorithm described in [16] that takes advantage of the block-tridiagonal structure that presents the infinitesimal generator. Unfortunately, the linear system of equations obtained in VE has no longer such a block-tridiagonal structure, and therefore we must use a more general and less efficient algorithm. More concretely, we have used LU factorization. Figures 8–10 show that VE achieves a certain accuracy faster than the other techniques under study.

Figure 8 

Relative error in 𝑃b versus computation cost.

Figure 9 

Relative error in 𝑃ns versus computation cost.

Figure 10 

Relative error in 𝑁ret versus computation cost.

5. Conclusions

Multiserver retrial systems have not an exact solution when the number of servers is higher than two, as their state space presents space heterogeneity along an infinite dimension. For that reason, it is mandatory to develop approximate techniques in order to solve these systems. To the best of our knowledge, all the techniques studied in the literature to solve these systems are based on their steady state probabilities. In this paper we propose an alternative technique based on a different metric: the relative state values and the Howard equations that relate them, instead of the balance equations. With this technique, truncation of the state space can be done in a more efficient way, as the state values outside the truncated state space are extrapolated from some known state values. In order to preserve the linearity of the resulting system of equations, we have only used polynomials as extrapolation functions.

In a first part, we have studied the use of different orders for the extrapolation polynomials. Later, we have compared the new technique with two well-known approaches appeared in the literature [8, 11] in terms of accuracy and computational cost. Results show that the proposed technique highly improves the previous approaches in terms of computational cost and, specially, in terms of accuracy, so its use is highly recommended.

Acknowledgments

This work has been supported by the Spanish government under Projects TIN2010-21378-C02-02 and TIN2008-06739-C04-02/TSI and by Comunidad de Madrid through Project S-2009/TIC-1468.