Throughput Maximization of Queueing Networks with Simultaneous Minimization of Service Rates and Buffers

Abstract

The throughput of an acyclic, general-service time queueing network was optimized, and the total number of buffers and the overall service rate was reduced. To satisfy these conflicting objectives, a multiobjective genetic algorithm was developed and employed. Thus, our method produced a set of efficient solutions for more than one objective in the objective function. A comprehensive set of computational experiments was conducted to determine the efficacy and efficiency of the proposed approach. Interesting insights obtained from the analysis of a complex network may assist practitioners in planning general-service queueing networks.

1. Introduction

In this study, the maximization of throughput (Θ) (the number of jobs, parts, clients, etc., served per unit of time) in an acyclic, general-service time queueing network (for an example, see Figure 1) was evaluated. To obtain the maximum Θ, the minimum number of buffers (𝐊={𝐾₁,𝐾₂,…,𝐾_𝑛}) and service rates (𝝁={𝜇₁,𝜇₂,…,𝜇_𝑛}) that must be allocated to a queueing network in a given topology and external arrival rate (Λ={Λ₁,Λ₂,…,Λ_𝑛}) was determined. Potential users of general-service time, finite-queueing network-based optimization models include computer scientists and industrial engineers. Indeed, these models may help to understand and improve various real-life systems, including manufacturing [1–5], production [6–8], and health [9–11] systems, urban or pedestrian traffic [12, 13], computer and communication systems [14–17], web-based applications with tiered configurations [18], and quality-of-service (QoS) requirements measured in terms of response times, throughput, server availability, and security [19].

Figure 1 

A complex network (adapted from Smith and Cruz [20]).

This study focused on single-server queueing networks with exponentially distributed interarrival times and generally distributed service times, configured in an arbitrary acyclic, series-parallel topology. An example of the type of network under consideration is shown in Figure 1. In particular, buffer allocation, server allocation, and throughput tradeoff were optimized in networks of 𝑀/𝐺/1/𝐾 queues. Thus, in Kendall [21] notation, we focused on Markovian arrivals, generally distributed service times, a single server, and the total capacity of 𝐾 items, including the item of service.

Indeed, there is a critical tradeoff between the overall number of buffers and service rates and the resulting throughput. Buffers and service capacities can be very expensive; thus, the total number of buffers and the overall service capacity should be reduced as much as possible. On the other hand, the highest possible network throughput is also desired. Unfortunately, throughput is directly affected by the number of buffers allocated, where an increase in buffers generally leads to a higher throughput. Likewise, the service capacity also directly affects the throughput.

Figure 2 shows the throughput, Θ, for a single 𝑀/𝐺/1/𝐾 queue with 𝑠²=1.5 (squared coefficient of variation of the service time) and Λ=5 users per time unit (external arrival rate), as a function of several values for buffer size, 𝐾, and service rate, 𝜇 (see (2.1) and (2.2)), as well as the respective contour plot. Similar throughput behavior is also observed in a network of queues. The surface of the plot shown in Figure 2 is smooth, and convexity is apparent, which is similar to the results of simple queueing networks [22, 23]. However, the top of the surface plot near the maximum throughput is flat, which creates difficulties for traditional optimization methods. For instance, Smith and Cruz [20] used the Powell method and multiple starts to avoid premature convergence to a local optimum and to derive a successful optimization algorithm.

(a) Throughput (Θ) versus service rate (𝜇) and buffer size (𝐾)

(b) Contour plot

Figure 2 

Results of a single 𝑀/𝐺/1/𝐾 queue for Λ=5.0.

From a modeling point of view, throughput maximization can be defined by a mixed-integer mathematical programming formula, where the total buffer and server costs are minimized, and the throughput, subject to integer buffer allocation and nonnegative service rates, is maximized. By defining a queueing network as a digraph of 𝐺(𝑁,𝐴), where 𝑁 is a finite set of nodes, and 𝐴 is a finite set of arcs, the mixed-integer mathematical programming formula was obtained [24]minimize𝐹(𝐊,𝝁) (1.1) subject to𝐾_𝑖∈{}𝜇1,2,…,∀𝑖∈𝑁,_𝑖≥0,∀𝑖∈𝑁, (1.2) where the decision variables 𝐾_𝑖 and 𝜇_𝑖 indicate the total capacity of the service and the service rate for the 𝑖th 𝑀/𝐺/1/𝐾 queue, respectively. The objective functions, 𝐹(𝐊,𝝁)≡(𝑓₁(𝐊),𝑓₂(𝝁),−𝑓₃(𝐊,𝝁), are the total buffer allocation, 𝑓₁(𝐊)=∑_{∀𝑖∈𝑁}𝐾_𝑖, the overall service allocation, 𝑓₂(𝝁)=∑_{∀𝑖∈𝑁}𝜇_𝑖, and the overall throughput, 𝑓₃(𝐊,𝝁)=Θ(𝐊,𝝁).

Throughput is often modeled as a constraint that must be greater than a target minimum, rather than as an objective that must be maximized [20, 25]. However, to successfully solve the problem, throughput constraints must be relaxed. Thus, parameters such as the threshold throughput (Θ_𝜏) must be determined beforehand. However, establishing an appropriate threshold is not a trivial task. Moreover, it is possible that a small decrease in throughput can result in a significant reduction in buffer allocation (and costs). The tradeoff between throughput and the number of buffers is not apparent in a single-objective formula. Indeed, the weights (𝜔) of a single-objective function have a significant effect on both the objectives and parameters, including errors (𝜀) on performance measure estimates and threshold throughput (Θ_𝜏). Thus, weight determination is difficult, and the results of single objective optimization techniques can be arbitrary.

In this study, an optimization approach was developed that determines the entire set of Pareto-optimal solutions. Thus, our method produces a set of efficient solutions for more than one objective in the objective function [26]. With the proposed approach, the decision maker is able to evaluate the effect of solution replacement. Moreover, the multiobjective approach also allows the user to increase one objective (e.g., throughput) while simultaneously reducing another objective (e.g., buffer and service rate allocation). A multiobjective evolutionary algorithm (MOEA) was used in combination with a generalized expansion method (GEM), which is a well-known method for obtaining accurate approximations of queueing network performance [27–29]. MOEAs are particularly suitable for multiobjective problems and have been shown to perform well in similar multiobjective problems of networks (e.g., see Carrano et al. [30] and references therein).

In this paper, a MOEA, specifically developed to multiobjective optimization, is presented (see Section 2). Additionally, the GEM, a performance evaluation tool used to approximate throughput, is also presented. In Section 3, the results of a comprehensive set of computational experiments are presented to show the efficiency of the approach. Finally, the article is concluded in Section 4, where final remarks and suggestions for future research are discussed.

2. Proposed Algorithms

The exposition of proposed algorithms was conducted in two parts. First, the performance evaluation algorithm was presented, which allowed the overall performance of the system to be estimated in terms of overall throughput, Θ(𝐊,𝝁), for a given configuration of the buffer and service allocation. Then the proposed optimization algorithm was developed, which was applied to obtain the optimal buffer and service allocation.

2.1. Performance Evaluation Algorithm

2.1.1. Single Queues

To maximize the throughput, Θ(𝐊,𝝁) must be estimated. In a single 𝑀/𝐺/1/𝐾 queue, the estimation of Θ(𝐊,𝝁) can be achieved with a computationally efficient and accurate closed-form approximate expression of the blocking probability, 𝑝_𝐾. The method proposed by Smith [31], which is based on a two-moment approximation from Kimura [32], was employed𝑝_𝐾=𝜌^{√((2+𝜌𝑠2−√√𝜌+2(𝐾−1))/(2+𝜌𝑠2−√𝜌))}𝜌−1𝜌^{√(2((2+𝜌𝑠2−√√𝜌+(𝐾−1))/(2+𝜌𝑠2−√𝜌)))}−1, (2.1) where 𝜌<1 is the system utilization, which is the ratio between the total arrival rate and the service rate, 𝜌=𝜆/𝜇.𝑠² is the squared coefficient of variation of the service time, 𝑇_𝑠; thus, 𝑠²=Var(𝑇_𝑠)/E(𝑇_𝑠)². The results indicated that the approximation of 𝑝_𝐾 was accurate for a wide range of values [20, 25, 33].

In order to obtain the throughput in a finite 𝑀/𝐺/1/𝐾 single queue, we need to adjust the arrival rate. In fact a fraction 𝑝_𝐾 of the arrivals cannot join the system because they have come when there is no waiting space left. Thus the actual rate of arrivals to join the system must be adjusted accordingly. Since Poisson arrivals see time averages (the PASTA property), it follows that the effective arrival rate seen by the servers is 𝜆_eff=𝜆(1−𝑝_𝐾) [34]. Thus, the throughput may be given by𝜃=𝜆_eff=𝜆1−𝑝_𝐾. (2.2)

2.1.2. Network of Queues

For a network of queues, the estimation of throughput is considerably more complicated. The generalized expansion method (GEM) is an algorithm that has been successfully used to estimate the performance of arbitrarily configured, finite queueing, and acyclic networks [29]. The GEM is a combination of node-by-node decomposition and repeated trials, where each queue is analyzed separately, and corrections are made to account for interrelated effects between network queues. The GEM uses type I blocking, where the upstream node becomes blocked if the service for an individual customer is complete and the queue at the downstream node is full. This is often referred to as “blocking after service,” which is prevalent in most production, manufacturing, and transportation systems.

With the help of Figure 3, we now describe the GEM. Firstly, we remark that the exponential distribution is a good approximation for the interdeparture times of entities leaving an 𝑀/𝐺/1/𝐾 node. In fact, it is possible to show the quasireversibility of a broader class of finite queues, which are the state-dependent 𝑀/𝐺/𝐶/𝐶 queues [35]. When those entities that are lost are included, the output stream is Poisson. This assumption is supported by several empirical results [7, 8, 13, 20, 25, 36]. The following three stages are involved in the GEM: network reconfiguration, parameter estimation, and feedback elimination.

Figure 3 

Generalized expansion method.

Network Reconfiguration
This stage involves reconfiguring the network. An auxiliary vertex ℎ_𝑗 is created, which is modeled as an 𝑀/𝐺/∞ queue with service rate 𝜇_ℎ and precedes each finite queue 𝑗. When an entity leaves 𝑖, vertex 𝑗 may be blocked with probability 𝑝_𝐾𝑗 or unblocked with probability (1−𝑝_𝐾𝑗). Under blocking, the entities are rerouted to vertex ℎ_𝑗 for a delay while node 𝑗 is busy. After this delay, the entity may be blocked again with probability 𝑝_𝐾𝑗′, for a second delay period. Vertex ℎ_𝑗 accumulates the time an entity has to wait before entering vertex 𝑗 and the effective arrival rate to vertex 𝑗.

Parameter Estimation
In this stage, the parameters 𝑝_𝐾, 𝑝^′_𝐾, and 𝜇_ℎ are estimated (for clarity, we will omit the subscript for node 𝑗).(1)𝑝_𝐾 is obtained by means of a two-moment approximation recently developed by Smith [31]𝑝_𝐾=equation(2.1). (2.3)(2)𝑝^′_𝐾 is obtained with the following approximation from diffusion techniques given by Labetoulle and Pujolle [37]: 𝑝^′_𝐾=𝜇_𝑗+𝜇_ℎ𝜇_ℎ−𝜆𝑟^𝐾₂−𝑟^𝐾₁−𝑟₂^𝐾−1−𝑟₁^𝐾−1𝜇_ℎ𝑟₂^𝐾+1−𝑟₁^𝐾+1−𝑟^𝐾₂−𝑟^𝐾₁⁻¹, (2.4) where 𝑟₁ and 𝑟₂ are the roots to the polynomial 𝜆−𝜆+𝜇_ℎ+𝜇_𝑗𝑥+𝜇_ℎ𝑥²=0, (2.5) with 𝜆=𝜆_𝑗−𝜆_ℎ(1−𝑝^′_𝐾), 𝜆_ℎ is the actual arrival rate to the artificial holding node, and 𝜆_𝑗 is the actual arrival rate to the finite node 𝑗, given by 𝜆_𝑗=𝜆_𝑖1−𝑝_𝐾=𝜆_𝑖−𝜆_ℎ. (2.6)(3)𝜇_ℎ is calculated as follows using renewal theory: 𝜇_ℎ=2𝜇_𝑗1+𝜎²_𝑗𝜇²_𝑗, (2.7) where 𝜎²_𝑗 is the service time variance.

Feedback Elimination
The repeated visits to the holding nodes (due to the feedback) create strong dependence in the arrival process. Therefore, the repeated immediate feedback is eliminated. This is accomplished by giving the customer enough service time during the first passage through the holding node. The adapted service rate for the holding node 𝜇^′_ℎ then becomes 𝜇^′_ℎ=1−𝑝^′_𝐾𝜇_ℎ. (2.8)

Summary 1. The goal of GEM is to provide an approximation scheme to update the service rates of upstream nodes that take into account all blocking after service caused by downstream nodes 𝜇_𝑖⁻¹=𝜇_𝑖⁻¹+𝑝_𝐾𝜇^′_ℎ⁻¹. (2.9)
For each finite node 𝑗 in the network succeeding node 𝑖, we have simultaneous nonlinear equations in variables 𝑝_𝐾, 𝑝^′_𝐾, and 𝜇_ℎ, along with auxiliary variables such as 𝜆 and 𝜆_𝑖. Solving these equations simultaneously, we can compute all the performance measures of the network 𝜆=𝜆_𝑗−𝜆_ℎ1−𝑝^′_𝐾𝜆_𝑗,=𝜆_𝑖1−𝑝_𝐾𝜆_𝑗,=𝜆_𝑖−𝜆_ℎ𝑝_′𝐾,=𝜇_𝑗+𝜇_ℎ𝜇_ℎ−𝜆𝑟^𝐾₂−𝑟^𝐾₁−𝑟₂^𝐾−1−𝑟₁^𝐾−1𝜇_ℎ𝑟₂^𝐾+1−𝑟₁^𝐾+1−𝑟^𝐾₂−𝑟^𝐾₁⁻¹,𝑧=𝜆+2𝜇_ℎ²−4𝜆𝜇_ℎ𝑟₁,=𝜆+2𝜇_ℎ−𝑧^1/22𝜇_ℎ𝑟₂,=𝜆+2𝜇_ℎ+𝑧^1/22𝜇_ℎ𝑝_𝐾,=equation(2.1). (2.10)(2.11)(2.12)(2.13)(2.14)(2.15)(2.16)(2.17)
Equation (2.10) through (2.13) is related to the arrivals and feedback in the holding node. Equation (2.14) through (2.16) is used to solve (2.13) with 𝑧 used as a dummy parameter for simplicity. Lastly, (2.1) gives the blocking probability for the 𝑀/𝐺/1/𝐾 queue. Thus, we essentially have five equations to solve (2.10)–(2.13) and (2.1).

2.2. Optimization Algorithm

For the network under consideration, MOEAs appeared to be a suitable choice for the multiobjective maximization of throughput. MOEAs are optimization algorithms that perform an approximate global search based on information obtained from the evaluation of several points in the search space [38, 39]. The population of points that converge to an optimal value are obtained through the application of genetic operators such as mutation, crossover, selection, and elitism. Each one of these operators characterizes an instance of a MOEA and can be implemented in several different ways. Additionally, MOEA convergence is guaranteed by assigning a value of fitness to each population member and preserving diversity. In fact, recent successful applications of GAs for single-objective applications were reported by Lin [40] and Calvete et al. [41], whereas Carrano et al. [30] employed GAs for multiple-objective applications. Additionally, the efficiency of GAs is well established for multiobjective problems [42, 43]. Many references are provided by the aforementioned authors.

The instance of MOEA used in this study was based upon the elitist nondominated sorting genetic (NSGA-II) algorithm of Deb et al. [44], which is shown in Algorithm 1. In the application of GAs for multiobjective optimization, the selection operator and elitism operator must be specifically structured to correctly identify optimal conditions as shown shortly. Elitism is based on the concept of dominance. Point 𝐱_𝑖=(𝑥_𝑖1,𝑥_𝑖2,…,𝑥_𝑖𝑛) dominates point 𝐱_𝑗=(𝑥_𝑗1,𝑥_𝑗2,…,𝑥_𝑗𝑛) if 𝐱_𝑖 is superior to 𝐱_𝑗 in one objective (𝑓_𝑘(𝐱_𝑖)<𝑓_𝑘(𝐱_𝑗), for minimization) and is not inferior in any other objective (𝑓_ℓ(𝐱_𝑖)≯𝑓_ℓ(𝐱_𝑗), for minimization).

algorithm  read graph, arrival, service rates, 𝐺(𝑁,𝐴),Λ𝑖∀𝑖∈𝑁 𝑃1←𝐆𝐞𝐧𝐞𝐫𝐚𝐭𝐞𝐈𝐧𝐢𝐭𝐢𝐚𝐥𝐏𝐨𝐩𝐮𝐥𝐚𝐭𝐢𝐨𝐧(popSize)  for 𝑖=1 until numGendo   /* generate offspring by crossover and mutation */   𝑄𝑖←𝐌𝐚𝐤𝐞𝐍𝐞𝐰𝐏𝐨𝐩(𝑃𝑖)   /* combine parent and offspring */   𝑅𝑖←𝑃𝑖∪𝑄𝑖   /* find nondominated fronts ℱ=(ℱ1,ℱ2,…) */   ℱ←𝐅𝐚𝐬𝐭𝐍𝐨𝐧𝐃𝐨𝐦𝐢𝐧𝐚𝐭𝐞𝐝𝐒𝐨𝐫𝐭(𝑅𝑖)   /* find new population by */   /* the crowding-distance-assignment */   𝑃𝑖+1←𝐆𝐞𝐧𝐞𝐫𝐚𝐭𝐞𝐍𝐞𝐰𝐏𝐨𝐩𝐮𝐥𝐚𝐭𝐢𝐨𝐧(𝑅𝑖)  end for  𝑃numGen+1←𝐄𝐱𝐭𝐫𝐚𝐜𝐭𝐏𝐚𝐫𝐞𝐭𝐨𝐒𝐞𝐭(𝑃numGen)  write 𝑃numGen+1 end algorithm

Algorithm 1 

Elitist multiobjective genetic algorithm (NSGA-II).

For instance, Figure 4 displays the points for a two-dimension minimization problem. In the figure, point V is dominated by point I, but not by points II, III, and IV. Point VI is dominated by points I, II, and III, but not by point IV. The best front includes points I through IV and is an approximation for the Pareto set, which is the set of points that are superior to other points. To perform elitism, an algorithm commonly referred to as the fast nondominated sorting algorithm was employed (details may be found in Deb et al. [44]). This algorithm separates the individuals in the population into several layers or fronts ℱ_𝑖, such that the solutions in ℱ₁ are nondominated, and every solution in a given front ℱ_𝑖,  𝑖>1, is dominated by at least one solution in ℱ_𝑖−1, and not by any solution in ℱ_𝑗,  𝑗≥𝑖. This can be achieved in 𝒪(𝑛log𝑛) time [44].

Figure 4 

Dominated (▪) and nondominated (•) points.

Selection is performed by sequentially selecting points from each nondominated front (ℱ₁,ℱ₂,…) until the number of required individuals for the next iteration is obtained. Criteria must be applied if, after the addition of a group of individuals from ℱ_𝑖, the maximum number of individuals is exceeded. The algorithm computes a measure of diversity (the crowding distance, as defined by Deb et al. [44]) to ensure the highest possible diversity. Thus, only the points with the largest crowding distance are kept for future iterations, as shown in Figure 5.

Figure 5 

Illustration of the crowding distance.

Crossover and mutation are somewhat independent of the multiobjective nature of the problem but are highly dependent on the application. For the problem at hand, a crossover mechanism known as “uniform” was selected [45]. Uniform crossover is popular in multivariable encodings due to its efficiency in identifying, inheriting, and protecting common genes, as well as recombining noncommon genes [46, 47]. In this mechanism, crossovers were performed for each variable with a probability (rateCro) that is in accordance with the crossover operator. The crossover operator used in the algorithm was the simulated binary crossover operator (SBX) [48, 49], as illustrated in Figure 6. SBX is quite convenient for real-coded GAs because it is able to simulate binary crossover operators but avoids reencoding the variables. The children (𝑥_{𝑖,(⋅,𝑡+1)}) were calculated from the parents (𝑥_{𝑖,(⋅,𝑡)}) according to the following equation:𝑥_{𝑖,(1,𝑡+1)}=0.51+𝛽𝑥_{𝑖,(1,𝑡)}+1−𝛽𝑥_{𝑖,(2,𝑡)},𝑥_{𝑖,(2,𝑡+1)}=0.51−𝛽𝑥_{𝑖,(1,𝑡)}+1+𝛽𝑥_{𝑖,(2,𝑡)}, (2.18) where 𝛽 is a random variable obtained from the following probability distribution function:𝑓𝛽=0.5𝜂+1𝛽^𝜂1,if𝛽≤1,0.5𝜂+1𝛽^𝜂+2,otherwise. (2.19)

Figure 6 

Chromosome representation and simulated binary crossover (SBX).

The function was designed to create a child solution that possesses a similar search power to a single-point crossover of binary-coded GAs [48]. By adjusting 𝜂, several different weights (𝛽) can be generated to produce children that are similar to their parents (i.e., large 𝜂) or not (small 𝜂). Several distributions are shown in Figure 7.

Figure 7 

Probability density function of 𝛽.

For each individual gene (the decision variables 𝐾_𝑖 and 𝜇_𝑖), the mutation scheme occurs with a specific probability (rateMut). As suggested by Deb and Agrawal [48], Gaussian perturbations were added to the decision variables, 𝐾_𝑖+𝜀_𝑖 and 𝜇_𝑖+𝜀_𝑁+𝑖, for all 𝑖∈𝑁, with 𝜀_𝑖∼𝒩(0,1), 𝑖∈{1,2,…,2𝑁}.

After crossover and mutation, constraints (1.2) may no longer apply. To guarantee feasibility, the values of integer variables were rounded accordingly and were readjusted by applying reflection operators𝐾_rfl𝑖=1+𝐾_𝑖,𝜇−1_rfl𝑖=𝜇_lowlim𝑖+𝜇_𝑖−𝜇_lowlim𝑖, (2.20) where 1 is the lower limit of buffer allocation, 𝜇_lowlim𝑖 is the lower limit of service allocation (to ensure that 𝜌<1 is applicable), 𝐾_𝑖 and 𝜇_𝑖 are the resulting values after crossover and mutation, and 𝐾_rfl𝑖 and 𝜇_rfl𝑖 are the results after reflection. The proposed scheme generates feasible solutions without avoiding or favoring any particular solution.

Recently, the stopping criterion of multiobjective optimization evolutionary algorithms has been analyzed in detail (see, e.g., Rudenko and Schoenauer [50] and Martí et al. [51]). Evidently, the maximum number of generations (numGen) plays an important role in the quality of the solutions. However, increasing the number of generation may not be ideal because computational time is wasted when many iterations do not lead to a significant improvement. Thus, Rudenko and Schoenauer [50] suggested that a superior stopping criterion is obtained when a fixed number of iterations are performed without improvement. To demonstrate the complexity of the issue, Rudenko and Schoenauer [50] conducted a comprehensive set of computational experiments. Their results revealed that an obvious stopping criterion, such as the entire population possessing a rank of 1, did not indicate that evolution should be terminated. The authors proposed a local stopping criterion that computes a measure of the stability of nondominated solutions after each iteration. Another global stopping criterion was recently proposed by Martí et al. [51]. This sophisticated method views population evolution as a dynamic system, where the state of the system can be estimated by a Kalman filter. For the sake of simplicity, the criterion of Rudenko and Schoenauer [50] was employed in this study. This criterion is based on the stabilization of the maximal crowding distance, 𝑑_𝑙, measured over 𝐿 generations, and is calculated by the following standard deviation:𝜎_𝐿=1𝐿^𝐿_𝑙=1𝑑_𝑙−𝑑_𝐿². (2.21)

As shown in (2.21), 𝑑_𝐿 is the average of 𝑑_𝑙 over 𝐿 generations. Moreover, (2.21) indicates that the MOEA stops when 𝜎_𝐿<𝛿_lim. Rudenko and Schoenauer [50] suggested that 𝜎_𝐿 does not depend on the actual values of the objective function because crowding distances are normalized. Furthermore, they also suggested that 𝐿 and 𝛿_lim should be set to 40 and 0.02, respectively, which leads to a stopping criteria that is 𝜎₄₀≤0.02. According to empirical evidence, these values are compatible with the network under consideration.

3. Computational Results and Discussion

To use previous implementations of GEM based on the International Mathematics and Statistics Library (IMSL), the optimization algorithm was implemented in Fortran [31, 33]. The code is available from the corresponding author upon request and must be used for educational and research purposes only. The computational experiments were conducted to discover the suboptimal set of parameters that guarantee rapid convergence. Moreover, the analysis of a large and complex network of finite queues was also achieved with the proposed algorithm.

3.1. Setting the Parameters

Unfortunately, to ensure rapid convergence with a minimal amount of computational effort, the suboptimal set of parameters must be determined by trial and error, as indicated by previous studies on GAs. Thus, networks containing 3, 5, and 10 queues were used to set the parameters of MOEA (see Figure 8). For the sake of conciseness, only the results obtained from 3 and 10 nodes are presented (Figures 9–12). Different topologies of acyclic similarly sized networks were also tested, and the results (not presented) were similar to those obtained from series topologies.

(a) A three-node series topology

(b) A five-node series topology

(c) A ten-node series topology

Figure 8 

Tested topologies.

(a) 3-node series network

(b) 10-node series network

Figure 9 

Effect of crossover and mutation.

(a) 3-node series network

(b) 10-node series network

Figure 10 

Effect of population size.

(a) 3-node series network

(b) 10-node series network

Figure 11 

Effect of the mutation rate.

(a) 3-node series network

(b) 10-node series network

Figure 12 

Effect of 𝜂.

In this study, each factor was analyzed independently; specifically, each factor was varied one at a time while the other parameters were held constant. Montgomery [52] reminds us that the major disadvantage of an independent analysis is that it fails to account for interactions between variables. However, recent experiments reported by Cruz et al. [53] indicated that potential interactions were insignificant; thus, interactions between factors were neglected in this study.

Figure 9 presents the convergence speed (in terms of 𝜎_𝐿) as a function of the number of generations. It is possible to conclude that pure mutation could be used to determine the optimal solution (sometimes pure mutation solves the problem, see Mathieu et al. [54]). However, the SBX operator was also utilized because it removed irregularities from the convergence profile. The combination of pure mutation and SBX provided satisfactory results, regardless of the number of queues in the network.

The results in Figure 10 revealed that the population size (popSize) had a significant effect on algorithm convergence. However, the population size cannot be arbitrarily increased because the required computational effort would become prohibitive. Moreover, the performance of the algorithm was not affected by an increase in the number of network nodes.

Figure 11 displays the convergence rate as a function of the parameter rateMut. The results revealed that an increase in the mutation rate accelerated convergence; however, once a specific rate was attained, a further increase did not lead to improved convergence. Under the experimental conditions, mutation rates between 1 and 2% provided superior results, regardless of the number of nodes.

Figure 12 displays the convergence as a function of 𝜂, which controls the dispersion of 𝛽_𝑞 in the SBX operator, (2.19). A further improvement in the convergence speed was not observed for values of 𝜂 greater than 8.

Finally, Figures 13 and 14 illustrate the population evolution, from the starting point to the final generation. They show the population spreading over time to cover an increasing proportion of objective space.

(a) Initial population

(b) After 10 generations

(c) After 100 generations

(d) Final result after 4,000 generations

Figure 13 

Population evolution of a 3-node network.

(a) Initial population

(b) After 10 generations

(c) After 100 generations

(d) Final result after 4,000 generations

Figure 14 

Population evolution of a 10-node network.

In summary all of the attempted problems could be successfully solved by employing the following combination: a combined use of mutation and SBX, a population size of 400 individuals, a mutation rate of 2%, and a dispersion parameter (𝜂) of 8. Moreover, the convergence seemed to be fairly independent of the topology (results not shown, for split topologies, merge topologies, and so on), the external arrival (Λ), the squared coefficient of variation of the service time (𝑠²), and the number of nodes of the network. Additionally, to ensure that the computation is complete within a finite amount of time, the maximum number of generations (numGen) must be predefined. In this study, numGen was set to 4,000 generations.

3.2. Analysis of a Large Complex Network

The complex network presented in Figure 1 was extracted from the literature [20] and analyzed with the proposed method. Three different squared coefficients of variation were analyzed (𝑠²={0.5,1.0,1.5}) at an arrival rate (Λ₁) of 5.0. First, the convergence speed of the genetic algorithm was confirmed to be robust for this type of problem. The experimental setup was identical to the previous analysis. However, the results indicated that convergence was stable at 2,000 iterations. Moreover, as shown in Figure 15, the convergence was largely independent of the squared coefficient of variation.

Figure 15 

Convergence of the 16-node network from Figure 1.

The corresponding profiles are shown in Figure 16, including the contour plot and final surface after convergence. For comparison, an exact contour plot of a single-node queue is presented in Figure 2(b), and the resemblance between the two graphs was encouraging. However, the behavior of a given network cannot be predicted without the use of an algorithmic approach such as the one proposed here. A detailed analysis of the results in Figure 16 revealed that many different pairs of buffers and service rates can be selected for a given throughput. Additionally, it is possible to evaluate the results when the buffer size or service rate is so high that it does not have an effect on the throughput (i.e., when the respective contour lines are parallel to the axes). Moreover, with the proposed algorithm, important insights related to the target throughput are obtained. For example, the results in Figure 16(d) suggest that it is easier to increase the throughput from 2.6 to 3.1 (20% improvement) than 4.1 to 4.5 (10% improvement). Contour lines that are far apart indicate that further improvements can be achieved only by dramatically increasing the buffer size and service rate.

(a) Final surface for 𝑠2=0.5

(b) Contour plot for 𝑠2=0.5

(c) Final surface for 𝑠2=1.5

(d) Contour plot for 𝑠2=1.5

Figure 16 

Final results for the 16-node network from Figure 1.

4. Conclusions and Final Remarks

In this study, a multiobjective approach was developed to maximize the throughput of single server, general queueing networks. By combining the generalized expansion method (GEM) with a multiobjective evolutionary algorithm (MOEA), insightful Pareto curves were obtained. These curves display the tradeoff between throughput, total buffer allocation, and overall service allocation.

Future investigations should be conducted to determine the applicability of this methodology for the determination of other optimal conditions in finite queueing networks. For instance, this method could be applied to optimize throughput in finite general, multiserver queueing networks or queueing networks with loops. Thus, the method could be used to model systems that lead to a reverse stream of products. Moreover, future research should be conducted to evaluate the algorithms in real-life situations.

Acknowledgments

The research of Frederico Cruz has been partially funded by CNPq (Conselho Nacional de Desenvolvimento Científico e Tecnológico; Grants, 201046/1994-6, 301809/1996-8, 307702/2004-9, 472066/2004-8, 304944/2007-6, 561259/2008-9, 553019/2009-0, 550207/2010-4, 501532/2010-2, 303388/2010-2), by CAPES (Coordenação de Aperfeiçoamento de Pessoal de Nível Superior; grant BEX-0522/07-4), and by FAPEMIG (Fundação de Amparo à Pesquisa do Estado de Minas Gerais; Grants, CEX-289/98, CEX-855/98, TEC-875/07, CEX-PPM-00401/08, CEX-PPM-00390-10). The Brazilian government funding agencies mentioned above had no role in the study.