# Efficient Method to Approximately Solve Retrial Systems with Impatience

Authors: Jose Manuel Gimenez-Guzman; M. Jose Domenech-Benlloch; Vicent Pla; Jorge Martinez-Bauset; Vicente Casares-Giner
Journal: Journal of Applied Mathematics
 (2012)

---

## 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.


---

## Body

## 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,…,C}×ℤ+ (being C 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. Section2 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 Figure1. Users arrive following a Poisson process with rate λ to a system with C 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 μr, this session retries. The reattempt is successful if it finds a free server. Otherwise, the user leaves the system with probability Pi or returns to the retrial orbit with probability (1-Pi), 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,S(t)={K(t),M(t)}, where K(t) is the number of sessions being served and M(t) the number of users in the retrial orbit at time t. The state space of the process is defined by(2.1)S:={s=(k,m):k≤C;m∈Z+}.Figure2 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 Section2. 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 whenC>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 (π(s)) using the balance equations in order to compute the desired performance parameters, that is, solving the linear system of equations:(3.1)π(s)∑s′≠sqss′=∑s′≠sπ(s′)qss′∀s∈S,
along with the normalization condition ∑sπ(s)=1, where qss′ represents the transition rate from state s to s′.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 a 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, r(s). 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 r=∑s∈𝒮π(s)r(s). In the value extrapolation  technique, the revenue function ℛ has to be defined so that the resulting average revenue r 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 states the system will collect a revenue for that action (r(s)), but as the number of transitions increases, the average revenue collected converges to r. The relative state value (v(s)) is equal to the difference between the total revenue incurred when the system starts at state s and the total revenue incurred in a system for which the revenue rate at all states is r:(3.2)v(s)=E[∫0∞(r(S(t))-r)dt∣S(0)=s].The Howard equations relate revenues, relative state values, and transition rates:(3.3)r(s)-r+∑s′qss′(v(s′)-v(s))=0∀s∈S.The Howard equations represent thepolicy 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(3.4)r(k,m)-r+λ(v(k+1,m)-v(k,m))+kμ(v(k-1,m)-v(k,m))+mμr(v(k+1,m-1)-v(k,m))=0ifk<C,r(C,m)-r+λ(v(C,m+1)-v(C,m))+Cμ(v(C-1,m)-v(C,m))+mμrPi(v(C,m-1)-v(C,m))=0ifk=C.As we can observe the number of states is infinite becausem 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)Ŝ:={s=(k,m):k≤C;m≤Q}.In general,Q is known as the truncation level. As we choose a higher value of Q, 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 Q.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 r, that is, |𝒮̂|+1 unknowns. However, as only the differences in the relative values appear in the Howard equations, we can set v(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 setsqss′=0foralls′∉𝒮,̂ 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 Q, the value extrapolation technique uses the state value of some states in 𝒮̂ to approximate v(C,Q+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 |𝒮̂|=(C+1)×(Q+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 (𝒮f) 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: 𝒲:𝒮̂f→ℝ. Hence, below we deal with states as if they were real values given as w=W(s). 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 f(w) and the set 𝒮̂f need to be chosen so that the parameters in f(w) have unambiguous values, that is, in the case of choosing a polynomial as the fitting function, the number of different points in 𝒮̂f 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 ai consists in minimizing the least mean squared error(3.6)E=∑w∈W(f(w)-v(w))2.Then the optimal values for theai’s can be computed by solving the equations(3.7)∂E∂ai=0∀i.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 andE=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 n=|𝒲(𝒮̂f)|=|𝒮̂f| points (w0,v(w0)),…,(wn-1,v(wn-1)), where there are not two identical wi, we can determine an (n-1)-th degree polynomial so that f(wi)=v(wi) for i=0,…,n-1, where(3.8)f(w)=a0+a1w+a2w2+⋅⋅⋅+an-1wn-1.The interpolating polynomial satisfies the followingn linear equations:(3.9)f(wi)=a0+a1wi+a2wi2+⋅⋅⋅+an-1win-1=v(wi)i=0,…,n-1,
which in a matrix form are(3.10)Aa=[1w0…w0n-11w1…w1n-1⋮⋮⋱⋮1wn-1…wn-1n-1][a0a1⋮an-1]=[v(w0)v(w1)⋮v(wn-1)]=b.The matrix of coefficients of this system (A) is a Vandermonde matrix, whose determinant is nonvanishing and therefore A 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 n points. However, Vandermonde matrices are often badly conditioned, specially if some wi 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(3.11)L(w)=∑j=0n-1v(wj)lj(w)
of Lagrange basis polynomials(3.12)lj(w)=∏i=0i≠jn-1w-wiwj-wi=w-w0wj-w0⋅⋅⋅w-wj-1wj-wj-1w-wj+1wj-wj+1⋅⋅⋅w-wn-1wj-wn-1.For the truncated problem of interest and as shown in Figure3, we will have a Howard equation in which appears v(C,Q+1), that is a state value of a state that does not belong to 𝒮̂. Therefore, we must approximate the value v(C,Q+1) by using some relative state values of states belonging to 𝒮̂. It is important to emphasize that for the extrapolation of v(C,Q+1) we only use states from the last row of the model shown in Figure 3, that is, states of the form s=(C,m), with varying m. With this choice, we define the mapping 𝒲 as 𝒲((C,m))=m. Moreover, we use an (n-1)-th degree polynomial that interpolates the n points in 𝒮f:={si=(C,Q-i)∣i=0,…,n-1} and then 𝒲(𝒮f)={wi=Q-i∣i=0,…,n-1}:(3.13)w0=Q⟶v(w0)=v(C,Q),w1=Q-1⟶v(w1)=v(C,Q-1),⋮wj=Q-j⟶v(wj)=v(C,Q-j),⋮wn-1=Q-(n-1)⟶v(wn-1)=v(C,Q-(n-1)).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(n-1)-th degree polynomial is(3.14)v(n)(C,Q+1)=L(n)(Q+1)=∑j=0n-1v(C,Q-j)lj(Q+1).For example, in the case of linear extrapolation (n=2), we use (Q,v(C,Q)) and (Q-1,v(C,Q-1)), having(3.15)v(2)(C,Q+1)=L(2)(Q+1)=v(C,Q)l0(Q+1)+v(C,Q-1)l1(Q+1)=v(C,Q)(Q+1)-(Q-1)Q-(Q-1)+v(C,Q-1)(Q+1)-Q(Q-1)-Q=2v(C,Q)-v(C,Q-1).Following a similar procedure, we can obtain the next relationship forn=3 and n=4:(3.16)v(3)(C,Q+1)=3v(C,Q)-3v(C,Q-1)+v(C,Q-2),v(4)(C,Q+1)=4v(C,Q)-6v(C,Q-1)+4v(C,Q-2)-v(C,Q-3).In general, for(n-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(3.17)v(n)(C,Q+1)=∑k=0n-1(-1)k(nk+1)v(C,Q-k),
where n 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,r(s) is the revenue rate obtained when the system is in state s. Therefore, we must define the revenue as the performance parameter we want to compute. The effect of that action is that the computed r will be the performance parameter we are looking for. Additionally, the inputs r(s) in the Howard equations must be properly set. Table 1 gives several examples on how r(s) can be set in order to obtain certain performance parameters such as: blocking probability Pb=Prob{K=C}, mean number of users in the retrial orbit Nret=E[M], nonservice probability Pns (probability of a user leaving the system due to impatience without obtaining service), probability of being in a certain state π(K,M), and probability of having K busy servers B(K).Table 1
Revenue function definition.
Blocking probabilityPbr(k,m)=1 for k=C, forallmr(k,m)=0 otherwiseNonservice probabilityPnsr(k,m)=mμrPi/λ for k=C, forallmr(k,m)=0 otherwiseMean number of users retryingNretr(k,m)=m  forallk, forallmProbability of being in state(K,M)π(K,M)r(k,m)=1 for k=K, m=Mr(k,m)=0 otherwiseProbability of havingK busy serversB(K)r(k,m)=1 for k=K, forallmr(k,m)=0 otherwiseAs 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, whenr(C,m)=1,for all m, and 0 in the rest of states, r(k,m)=0,k≠C, for all m.


### 3.4. Effect of the Value Extrapolation into the Howard Equations

In our problem, and as mentioned above, we will only have to replacev(C,Q+1) by its approximate value in the Howard equation that corresponds to the state v(C,Q). As an example, if we use linear extrapolation (n=2), that equation becomes(3.18)r(C,Q)-r+v(C,Q)(-λ-Cμ-QPiμr)+λv(C,Q+1)+Cμv(C-1,Q)+QPiμrv(C,Q-1)=r(C,Q)-r+v(C,Q)(λ-Cμ-QPiμr)+Cμv(C-1,Q)+(QPiμr-λ)v(C,Q-1)=0.Asv(C,Q+1) no longer appears in the Howard equations, the linear system of equations we have consists of (C+1)×(Q+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 xT=b, where x is a vector with the (C+1)×(Q+1) unknowns (r and the relative state values v(s)) and b are the negative revenue rates for the different states:(3.19)x=[rv(0,1)⋯v(0,Q)v(1,0)⋯v(C,Q)],b=[-r(0,0)-r(0,1)⋯-r(C,Q)].MatrixT represents the matrix of coefficients and can be constructed making all the elements in the first row of matrix T0 equal to −1.MatrixT0 is given by(3.20)T0=[A10A00⋯00A21A11⋯00⋮⋮⋱⋮⋮00⋯A1C-1A0C-100⋯A2CA1C],
where the submatrices are defined as(3.21)A0k=(k+1)μI,for0≤k≤(C-1),A2k=[λμr0⋯000λ2μr⋯00⋮⋮⋮⋱⋮⋮000⋯λQμr000⋯0λ],for1≤k≤C,A1k=[α00⋯00α-μr0⋯000α-2μr⋯0⋮⋮⋮⋱⋮000⋯α-Qμr],for0≤k≤(C-1),α=-λ-kμ.Whenk=C, using linear (n=2) and quadratic (n=3) extrapolation, we obtain, respectively,(3.22)A1C=[βPiμr⋯00λβ-Piμr⋯000λ⋯00⋮⋮⋱⋮⋮00⋯β-(Q-1)PiμrQPiμr-λ00⋯λλ-Cμ-QPiμr],A1C=[βPiμr⋯00λβ-Piμr⋯000λ⋯00⋮⋮⋱⋮⋮00⋯(Q-1)Piμrλ00⋯β-(Q-1)PiμrQPiμr-3λ00⋯λ2λ-Cμ-QPiμr],
where β=-λ-Cμ.In general, if the extrapolation is done withn≤Q+1 points, the matrix A1C is given as(3.23)A1C=[βPiμr⋯0λcQ(n)λβ-Piμr⋯0λcQ-1(n)0λ⋯0λcQ-2(n)⋮⋮⋱⋮⋮00⋯(Q-1)Piμrλc2(n)00⋯β-(Q-1)PiμrQPiμr+λc1(n)00⋯λ-λ-Cμ-QPiμr+λc0(n)],
where(3.24)cl(n)={(-1)l(nl+1)ifl<n,0ifl≥n.Note that the size of matrixT does not depend on the order of the polynomial used to perform the extrapolation; only the last column in the matrix A1C 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 systemxT=b can be obtained using the inverse matrix of T by doing x=bT-1. Note also that choosing a different performance parameter to solve will only affect to the values in b. Therefore, computing a second performance parameter will only increase the computation expenses by the cost of the product bT-1, as the rest of the process (specifically the computation of the inverse matrix T-1) 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 T 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 solvingp performance parameters simultaneously.


## 3.1. MDP Settings

As it has been aforementioned, the problem under interest has not a closed form solution whenC>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 (π(s)) using the balance equations in order to compute the desired performance parameters, that is, solving the linear system of equations:(3.1)π(s)∑s′≠sqss′=∑s′≠sπ(s′)qss′∀s∈S,
along with the normalization condition ∑sπ(s)=1, where qss′ represents the transition rate from state s to s′.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 a 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, r(s). 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 r=∑s∈𝒮π(s)r(s). In the value extrapolation  technique, the revenue function ℛ has to be defined so that the resulting average revenue r 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 states the system will collect a revenue for that action (r(s)), but as the number of transitions increases, the average revenue collected converges to r. The relative state value (v(s)) is equal to the difference between the total revenue incurred when the system starts at state s and the total revenue incurred in a system for which the revenue rate at all states is r:(3.2)v(s)=E[∫0∞(r(S(t))-r)dt∣S(0)=s].The Howard equations relate revenues, relative state values, and transition rates:(3.3)r(s)-r+∑s′qss′(v(s′)-v(s))=0∀s∈S.The Howard equations represent thepolicy 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(3.4)r(k,m)-r+λ(v(k+1,m)-v(k,m))+kμ(v(k-1,m)-v(k,m))+mμr(v(k+1,m-1)-v(k,m))=0ifk<C,r(C,m)-r+λ(v(C,m+1)-v(C,m))+Cμ(v(C-1,m)-v(C,m))+mμrPi(v(C,m-1)-v(C,m))=0ifk=C.As we can observe the number of states is infinite becausem 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)Ŝ:={s=(k,m):k≤C;m≤Q}.In general,Q is known as the truncation level. As we choose a higher value of Q, 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 Q.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 r, that is, |𝒮̂|+1 unknowns. However, as only the differences in the relative values appear in the Howard equations, we can set v(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 setsqss′=0foralls′∉𝒮,̂ 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 Q, the value extrapolation technique uses the state value of some states in 𝒮̂ to approximate v(C,Q+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 |𝒮̂|=(C+1)×(Q+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 (𝒮f) 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: 𝒲:𝒮̂f→ℝ. Hence, below we deal with states as if they were real values given as w=W(s). 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 f(w) and the set 𝒮̂f need to be chosen so that the parameters in f(w) have unambiguous values, that is, in the case of choosing a polynomial as the fitting function, the number of different points in 𝒮̂f 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 ai consists in minimizing the least mean squared error(3.6)E=∑w∈W(f(w)-v(w))2.Then the optimal values for theai’s can be computed by solving the equations(3.7)∂E∂ai=0∀i.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 andE=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 n=|𝒲(𝒮̂f)|=|𝒮̂f| points (w0,v(w0)),…,(wn-1,v(wn-1)), where there are not two identical wi, we can determine an (n-1)-th degree polynomial so that f(wi)=v(wi) for i=0,…,n-1, where(3.8)f(w)=a0+a1w+a2w2+⋅⋅⋅+an-1wn-1.The interpolating polynomial satisfies the followingn linear equations:(3.9)f(wi)=a0+a1wi+a2wi2+⋅⋅⋅+an-1win-1=v(wi)i=0,…,n-1,
which in a matrix form are(3.10)Aa=[1w0…w0n-11w1…w1n-1⋮⋮⋱⋮1wn-1…wn-1n-1][a0a1⋮an-1]=[v(w0)v(w1)⋮v(wn-1)]=b.The matrix of coefficients of this system (A) is a Vandermonde matrix, whose determinant is nonvanishing and therefore A 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 n points. However, Vandermonde matrices are often badly conditioned, specially if some wi 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(3.11)L(w)=∑j=0n-1v(wj)lj(w)
of Lagrange basis polynomials(3.12)lj(w)=∏i=0i≠jn-1w-wiwj-wi=w-w0wj-w0⋅⋅⋅w-wj-1wj-wj-1w-wj+1wj-wj+1⋅⋅⋅w-wn-1wj-wn-1.For the truncated problem of interest and as shown in Figure3, we will have a Howard equation in which appears v(C,Q+1), that is a state value of a state that does not belong to 𝒮̂. Therefore, we must approximate the value v(C,Q+1) by using some relative state values of states belonging to 𝒮̂. It is important to emphasize that for the extrapolation of v(C,Q+1) we only use states from the last row of the model shown in Figure 3, that is, states of the form s=(C,m), with varying m. With this choice, we define the mapping 𝒲 as 𝒲((C,m))=m. Moreover, we use an (n-1)-th degree polynomial that interpolates the n points in 𝒮f:={si=(C,Q-i)∣i=0,…,n-1} and then 𝒲(𝒮f)={wi=Q-i∣i=0,…,n-1}:(3.13)w0=Q⟶v(w0)=v(C,Q),w1=Q-1⟶v(w1)=v(C,Q-1),⋮wj=Q-j⟶v(wj)=v(C,Q-j),⋮wn-1=Q-(n-1)⟶v(wn-1)=v(C,Q-(n-1)).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(n-1)-th degree polynomial is(3.14)v(n)(C,Q+1)=L(n)(Q+1)=∑j=0n-1v(C,Q-j)lj(Q+1).For example, in the case of linear extrapolation (n=2), we use (Q,v(C,Q)) and (Q-1,v(C,Q-1)), having(3.15)v(2)(C,Q+1)=L(2)(Q+1)=v(C,Q)l0(Q+1)+v(C,Q-1)l1(Q+1)=v(C,Q)(Q+1)-(Q-1)Q-(Q-1)+v(C,Q-1)(Q+1)-Q(Q-1)-Q=2v(C,Q)-v(C,Q-1).Following a similar procedure, we can obtain the next relationship forn=3 and n=4:(3.16)v(3)(C,Q+1)=3v(C,Q)-3v(C,Q-1)+v(C,Q-2),v(4)(C,Q+1)=4v(C,Q)-6v(C,Q-1)+4v(C,Q-2)-v(C,Q-3).In general, for(n-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(3.17)v(n)(C,Q+1)=∑k=0n-1(-1)k(nk+1)v(C,Q-k),
where n 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,r(s) is the revenue rate obtained when the system is in state s. Therefore, we must define the revenue as the performance parameter we want to compute. The effect of that action is that the computed r will be the performance parameter we are looking for. Additionally, the inputs r(s) in the Howard equations must be properly set. Table 1 gives several examples on how r(s) can be set in order to obtain certain performance parameters such as: blocking probability Pb=Prob{K=C}, mean number of users in the retrial orbit Nret=E[M], nonservice probability Pns (probability of a user leaving the system due to impatience without obtaining service), probability of being in a certain state π(K,M), and probability of having K busy servers B(K).Table 1
Revenue function definition.
Blocking probabilityPbr(k,m)=1 for k=C, forallmr(k,m)=0 otherwiseNonservice probabilityPnsr(k,m)=mμrPi/λ for k=C, forallmr(k,m)=0 otherwiseMean number of users retryingNretr(k,m)=m  forallk, forallmProbability of being in state(K,M)π(K,M)r(k,m)=1 for k=K, m=Mr(k,m)=0 otherwiseProbability of havingK busy serversB(K)r(k,m)=1 for k=K, forallmr(k,m)=0 otherwiseAs 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, whenr(C,m)=1,for all m, and 0 in the rest of states, r(k,m)=0,k≠C, for all m.


## 3.4. Effect of the Value Extrapolation into the Howard Equations

In our problem, and as mentioned above, we will only have to replacev(C,Q+1) by its approximate value in the Howard equation that corresponds to the state v(C,Q). As an example, if we use linear extrapolation (n=2), that equation becomes(3.18)r(C,Q)-r+v(C,Q)(-λ-Cμ-QPiμr)+λv(C,Q+1)+Cμv(C-1,Q)+QPiμrv(C,Q-1)=r(C,Q)-r+v(C,Q)(λ-Cμ-QPiμr)+Cμv(C-1,Q)+(QPiμr-λ)v(C,Q-1)=0.Asv(C,Q+1) no longer appears in the Howard equations, the linear system of equations we have consists of (C+1)×(Q+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 xT=b, where x is a vector with the (C+1)×(Q+1) unknowns (r and the relative state values v(s)) and b are the negative revenue rates for the different states:(3.19)x=[rv(0,1)⋯v(0,Q)v(1,0)⋯v(C,Q)],b=[-r(0,0)-r(0,1)⋯-r(C,Q)].MatrixT represents the matrix of coefficients and can be constructed making all the elements in the first row of matrix T0 equal to −1.MatrixT0 is given by(3.20)T0=[A10A00⋯00A21A11⋯00⋮⋮⋱⋮⋮00⋯A1C-1A0C-100⋯A2CA1C],
where the submatrices are defined as(3.21)A0k=(k+1)μI,for0≤k≤(C-1),A2k=[λμr0⋯000λ2μr⋯00⋮⋮⋮⋱⋮⋮000⋯λQμr000⋯0λ],for1≤k≤C,A1k=[α00⋯00α-μr0⋯000α-2μr⋯0⋮⋮⋮⋱⋮000⋯α-Qμr],for0≤k≤(C-1),α=-λ-kμ.Whenk=C, using linear (n=2) and quadratic (n=3) extrapolation, we obtain, respectively,(3.22)A1C=[βPiμr⋯00λβ-Piμr⋯000λ⋯00⋮⋮⋱⋮⋮00⋯β-(Q-1)PiμrQPiμr-λ00⋯λλ-Cμ-QPiμr],A1C=[βPiμr⋯00λβ-Piμr⋯000λ⋯00⋮⋮⋱⋮⋮00⋯(Q-1)Piμrλ00⋯β-(Q-1)PiμrQPiμr-3λ00⋯λ2λ-Cμ-QPiμr],
where β=-λ-Cμ.In general, if the extrapolation is done withn≤Q+1 points, the matrix A1C is given as(3.23)A1C=[βPiμr⋯0λcQ(n)λβ-Piμr⋯0λcQ-1(n)0λ⋯0λcQ-2(n)⋮⋮⋱⋮⋮00⋯(Q-1)Piμrλc2(n)00⋯β-(Q-1)PiμrQPiμr+λc1(n)00⋯λ-λ-Cμ-QPiμr+λc0(n)],
where(3.24)cl(n)={(-1)l(nl+1)ifl<n,0ifl≥n.Note that the size of matrixT does not depend on the order of the polynomial used to perform the extrapolation; only the last column in the matrix A1C 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 systemxT=b can be obtained using the inverse matrix of T by doing x=bT-1. Note also that choosing a different performance parameter to solve will only affect to the values in b. Therefore, computing a second performance parameter will only increase the computation expenses by the cost of the product bT-1, as the rest of the process (specifically the computation of the inverse matrix T-1) 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 T 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 solvingp performance parameters simultaneously.


## 4. Results

In order to evaluate and compare the proposed technique, we have studied its performance in several scenarios. Lettingρ=λ/(Cμ), we have studied different system loads by modifying λ and keeping C=50 resource units and μ-1=180 s. The retrial phenomenon has been configured with μr-1=100 s and Pi=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 Q 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, Ψ∈{Pb,Pns,Nret}.


### 4.1. Value Extrapolation Evaluation

Table2 shows the minimum value of Q needed to obtain a relative error lower than 10-8 for different performance parameters and loads (columns) and for different orders of the extrapolation polynomials (rows). Note that VEx denotes the use of an extrapolation polynomial of order x=(n-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 ofQ to obtain relative errors (ϵ) lower than 10-8.
ϵPb<10-8ϵPns<10-8ϵNret<10-8ρ0.50.70.90.50.70.90.50.70.9VE1203261254164223757VE2143153213558173254VE3151848193153162650VE4122547173048142647VE512244412244491843VE6102041142644112239VE77213911244282140VE88173911233681939VE99193810223991334VE10101635102139101735VE11111831111637111837VE12121540122042121742VE13131443141943131843VE14142348262548142448VE15152556152956152556VE16162756182957282757Exact Value3.89·10-60.00450.13536.05·10-81.34·10-40.01105.74·10-50.09814.4789From Table2 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 VEx enforces us to use a model with Q≥x (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 Q 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 Table3 we show the minimum values of Q needed to obtain a relative error lower than 10-8 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 Q 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 Pb, Pns and Nret, respectively, when ρ=0.7 and for the different techniques deployed. Results show that, for a same value of Q, 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
MinimumQ value to obtain relative errors (ϵ) lower than 10-8.
ϵPb<10-8ϵPns<10-8ϵNret<10-8ρ=0.5ρ=0.7ρ=0.9ρ=0.5ρ=0.7ρ=0.9ρ=0.5ρ=0.7ρ=0.9FM233968294670254253NR203161254164223865VE88173911233681939Figure 5
Relative error inPb for different techniques.Figure 6
Relative error inPns for different techniques.Figure 7
Relative error inNret 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 Figures8, 9, and 10, we plot the time needed to achieve a certain relative error for Pb, Pns, and Nret 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 inPb versus computation cost.Figure 9
Relative error inPns versus computation cost.Figure 10
Relative error inNret versus computation cost.


## 4.1. Value Extrapolation Evaluation

Table2 shows the minimum value of Q needed to obtain a relative error lower than 10-8 for different performance parameters and loads (columns) and for different orders of the extrapolation polynomials (rows). Note that VEx denotes the use of an extrapolation polynomial of order x=(n-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 ofQ to obtain relative errors (ϵ) lower than 10-8.
ϵPb<10-8ϵPns<10-8ϵNret<10-8ρ0.50.70.90.50.70.90.50.70.9VE1203261254164223757VE2143153213558173254VE3151848193153162650VE4122547173048142647VE512244412244491843VE6102041142644112239VE77213911244282140VE88173911233681939VE99193810223991334VE10101635102139101735VE11111831111637111837VE12121540122042121742VE13131443141943131843VE14142348262548142448VE15152556152956152556VE16162756182957282757Exact Value3.89·10-60.00450.13536.05·10-81.34·10-40.01105.74·10-50.09814.4789From Table2 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 VEx enforces us to use a model with Q≥x (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 Q 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 Table3 we show the minimum values of Q needed to obtain a relative error lower than 10-8 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 Q 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 Pb, Pns and Nret, respectively, when ρ=0.7 and for the different techniques deployed. Results show that, for a same value of Q, 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
MinimumQ value to obtain relative errors (ϵ) lower than 10-8.
ϵPb<10-8ϵPns<10-8ϵNret<10-8ρ=0.5ρ=0.7ρ=0.9ρ=0.5ρ=0.7ρ=0.9ρ=0.5ρ=0.7ρ=0.9FM233968294670254253NR203161254164223865VE88173911233681939Figure 5
Relative error inPb for different techniques.Figure 6
Relative error inPns for different techniques.Figure 7
Relative error inNret 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 Figures8, 9, and 10, we plot the time needed to achieve a certain relative error for Pb, Pns, and Nret 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 inPb versus computation cost.Figure 9
Relative error inPns versus computation cost.Figure 10
Relative error inNret 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.



---

*Source: 186761-2012-01-22.xml*