Solving the minimum edge-dilation k-center problem by genetic algorithms

Computers & Industrial Engineering 113 (2017) 282–293

Contents lists available at ScienceDirect

Computers & Industrial Engineering journal homepage: www.elsevier.com/locate/caie

Solving the minimum edge-dilation k-center problem by genetic algorithms a,⁎

b

Dragan Matic , Jozef Kratica , Z. Maksimovic a b c

MARK

c

University of Banjaluka, Faculty of Mathematics and Natural Sciences, Mladena Stojanovica 2, 78000 Banjaluka, Bosnia and Herzegovina Mathematical Institute, Serbian Academy of Sciences and Arts, Kneza Mihaila 36/III, 11000 Belgrade, Serbia Military Academy, University of Defence, Generala Pavla Jurišica Šturma 33, 11000 Belgrade, Serbia

A R T I C L E I N F O

A B S T R A C T

Keywords: Edge-dilation k center Minimum stretch Genetic algorithms Combinatorial optimization

In this paper, we present two genetic algorithms (GA) for solving the minimum edge-dilation k-center (MEDKC) problem. Up to now, MEDKC has been considered only theoretically, and we present the first heuristic approaches for solving it. In both approaches, the assignment of non-center vertices is based on the distance to the center vertices, while modified crossover and mutation genetic operators, specific for each GA, keep individuals feasible during the entire optimization process. The proposed algorithms were empirically tested and compared by using various sets of differently sized test data: some theoretical classes of graphs with known optimal solutions and some special classes of graphs chosen from the literature for the k-minimum spanning tree and Roman domination problems. Obtained results indicate that newly proposed GAs can be reliably used in constructing routing schemes in complex computer networks.

1. Introduction Let G = (V ,E ,ω) be undirected, edge-weighted graph where V = {1,2,…,n} represents set of vertices, E is set of edges, |E| = m and ω is the edge-weight function. Let dω (u,v ) be the shortest path between the vertices u and v and let k > 0 be a given parameter. Let Π ⊂ V be a set of k vertices and let every vertex v ∈ V be assigned to exactly one vertex πv ∈ Π. The center distance between two vertices u,v ∈ V with respect to Π is defined as

d π (u,v ) = dω (u,πu ) + dω (πu,πv ) + dω (πv,v )

(1)

For the given set Π and for a pair of vertices u,v ∈ V the stretch is defined as the ratio of the center distance d π (u,v ) and the shortest distance in the graph, i.e.

stretch (u,v ) = d π (u,v )/ dω (u,v )

(2)

The stretch of a solution (Π,{πv}v ∈ V ) is then defined as the maximum stretch over all pairs of vertices (Konemann, Li, Parekh, & Sinha, 2004):

stretch (Π,{πv}v ∈ V ) =

max (u,v ) ∈ V × V

d π (u,v )/ dω (u,v )

(3)

In other words, the stretch of a solution is defined as the worst-case path-length increase factor relative to the shortest paths. More precisely, for every pair of vertices in the graph, the ratio of the length of the path taken by the chosen assignment to the length of the shortest available path between the same pair of vertices is considered. As it is already mentioned above, the stretch of a solution is then defined as the ⁎

maximum of this ratio among all vertex pairs. The minimum edge-dilation k-center problem (MEDKC) is introduced by Konemann et al. (2004). The aim of MEDKC is to identify the set Π containing k center vertices and to assign every other vertex to exactly one vertex from Π, such that the stretch (Π,{πv}v ∈ V ) is minimal. For a given graph G, let us hereinafter denote the optimal solution with optMEDKC (G ) = min{stretch (Π,{πv}v ∈ V )} over all possible subsets Π and assignments {πv}v ∈ V . Any unweighted graph can be considered as weighted one if we define that any edge-weight is equal to 1. Therefore, the MEDKC can also be considered on unweighted graphs. The main contribution of this paper is the development of two genetic algorithms, which solve the mentioned problem: for a given network and given number k, the GAs obtains the solution, i.e. GAs identify k center vertices (routers) and propose the assignment for all non-center vertices, minimizing the maximal stretch. As it will be shown in Section 5, proposed GAs obtained high-quality results on small, medium and large networks. 1.1. Motivation As it can be seen in Konemann et al. (2004), the motivation for this problem arises from a real application in determining routing schemes in complex computer networks. In a routing scheme, one needs to allow any source vertex to route messages to destination vertices, given the destination’s network identifier. For each host vertex, routing paths to

Corresponding author. E-mail addresses: [email protected] (D. Matic), [email protected] (J. Kratica), [email protected] (Z. Maksimovic).

http://dx.doi.org/10.1016/j.cie.2017.09.029 Received 22 May 2017; Received in revised form 14 September 2017; Accepted 18 September 2017 Available online 20 September 2017 0360-8352/ © 2017 Elsevier Ltd. All rights reserved.

Computers & Industrial Engineering 113 (2017) 282–293

D. Matic et al.

Nomenclature MEDKC GA GA1 GA2 G V n E ω

dω (u,v ) shortest path between the vertices u and v Π set of centers k number of centers (i.e. k = |Π|) πv element from Π assigned to a vertex v d π (u,v ) center distance between two vertices u,v stretch (u,v ) ratio of the center distance d π (u,v ) and the shortest distance between vertices stretch (Π,{πv}v ∈ V ) stretch of a solution optMEDKC (G ) optimal solution of MEDKC problem for a given graph G

minimum edge-dilation k-center problem genetic algorithm GA with binary encoded individuals GA with integer encoded individuals graph set of vertices number of vertices in V (n = |V |) set of edges edge weight function

the problem of long-term forecasting of energy consumption using optimized gene expression programming. The applicability and accuracy of the method are tested using actual data sets from 1971 to 2011. The test results confirmed that the method has higher accuracy than other metaheuristics used for solving this problem. In Abbasi-Pooya and Kashan (2017) the helicopter routing and crew exchange problem is proposed. Two different mathematical models are developed, which solved to optimality test instances of up to 11 jackets and 28 crew members, and for the real case of South Pars Gas Field of up to 11 jackets and 46 passengers. For solving instances of larger sizes a hybrid Grouping Evolution Strategy, is proposed which works based on the composition of the flight sorties rather than working with each nodes in isolation. Experimental results and direct comparison with two existing metaheuristic from the literature Particle Swarm Optimization (PSO) and pure Grouping Evolution Strategy algorithm (GES), show usefulness of proposed method. Mansouri, Kaboli, Ahmadian, and Selvaraj (2012) proposed a hybrid neuro-fuzzy-P.I. speed controller for BLDC. The efficiency and performance of the applied controlling system have been evaluated using various parameters and loads showing better results than those known from the literature. Modiri-Delshad et al. proposed iterated-based algorithms for solving Economic Dispatch problem in Modiri-Delshad, Koohi-Kamali, Taslimi, Kaboli, and Rahim (2013), Modiri-Delshad, Kaboli, Taslimi, Selvaraj, and Rahim (2013). The method is simple to implement and has acceptable speed and accuracy, which is confirmed by the comparison with the exact results obtained with commercial solver CPLEX. A new optimization algorithm, inspired by the rain drops behavior, called Rain-Fall Optimization (RFO) is proposed in Kaboli, Selvaraj, and Rahim (2017). The method is validated by applying on computationally expensive benchmark problems, which proved to offer better quality solutions than compared optimizers, and it is applied for solving Economic Dispatch (ED) problem. Vehicle routing problem with cross docks and split deliveries is solved by Wang, Jagannathan, Zuo, and Murray (2017), which proposed a mixed-integer linear programming formulation and a constructive heuristic with two-layer simulated annealing and tabu search. The first layer optimizes the allocation of trucks to cross docks while the second layer optimizes the visitation order to suppliers and retailers for trucks assigned to each cross dock. Computational experiments demonstrated that the proposed approach effectively solves large-size problems within a reasonable computational time. Another variant of vehicle routing problem, named capacitated vehicle routing problem and re-routing strategies under time-dependent traffic congestion, is solved by the multiple colonies artificial bee colony algorithm by Ng, Lee, Zhang, Wu, and Ho (2017). The computational results indicate the potential of using real time information for data-driven vehicle scheduling. Ant colony optimization is used for solving balancing U-type assembly lines problem in Li, Kucukkoc, and Tang (2017). Comparison with other methods from the literature demonstrate that the proposed model iterates fast and achieves competing results. A decomposition heuristic (DH) and a simulated annealing-based heuristic (SA) were

every other host is kept in its routing table. If one wants to ensure the optimal shortest-path routing for the entire network, for each vertex, in the routing table O (n) memory space size is required for storing information about routing paths to every other vertex. Since modern computer networks contain millions of nodes, it is obvious that storing optimal shortest-paths in routing tables is impossible due to the memory limits. This fact opens an area for research, with the aim to construct a routing scheme as efficient as possible, but having in mind memory and time limitations. There are two main objectives which are considered: minimizing the size of routing tables and/or minimizing the stretch of a routing scheme. In other words, the task of constructing compact routing schemes focuses on the trade-off between the stretch and the size of routing tables. A standard way for dealing with this task is to divide the network into areas. For each vertex, only the set of information about the paths from the same area are stored, while the routing between vertices from different areas are performed via special vertices (usually called routers). Thus, the size of the routing table of each vertex can be significantly reduced, with the hope that the size of the stretch of the entire network is still acceptable. The main idea of the approach presented in this paper follows the principle that some k vertices in a network can be declared as routers, while each of all other vertices is assigned to a router. The strategy of constructing such a compact routing scheme corresponds to solving the MEDKC problem, where routers are exactly the center vertices. 1.2. Recent metaheuristic approaches Metaheuristics proved to be useful in dealing with the problems of combinatorial optimizations on graphs. One broad class of metaheuristics are biologically inspired and have been investigated in recent research. Several techniques have been proposed for dealing with the problems of discrete optimization. An evolutionary technique Backtracking Search Algorithm (BSA) is used in Modiri-Delshad, Kaboli, Taslimi-Renani, and Rahim (2016) for solving economic dispatch problems considering the valve-point loading effects, prohibited operating zones, and multiple fuel options. The BSA method is shown to be both efficient and robust. In Kaboli, Selvaraj, and Rahim (2016) the artificial cooperative search algorithm is applied to the problem of long-term electric energy consumption forecasting. Apart from accuracy and reliability of this method, the balance between exploitation and exploration is achieved using single parameter, which greatly simplifies the fine-tuning of the search process. Rafieerad et al. (2017) developed the Multi Objective Particle Swarm Optimization method (MOPSO) to maximize adhesion and hardness of coating to a surface and an extensive enhancement in wear and corrosion resistance of the nanostructured coating was achieved. The Particle Swarm Optimization meta-heuristic algorithm PSO-PI controller is integrated with a rotational d-q current control method for finding the optimum values of the coefficients in Sebtahmadi, Azad, Kaboli, Islam, and Mekhilef (2017). The experimental results indicate that the proposed system produces high-quality sinusoidal currents at the output. Kaboli, Fallahpour, Selvaraj, and Rahim (2017) discussed 283

Computers & Industrial Engineering 113 (2017) 282–293

D. Matic et al.

problem is to identify the set of k central vertices and to organize the paths over them, in order to minimize the maximal stretch. Since it is proven that the MEDKC is NP hard, exact methods cannot be used for solving large-scaled networks. Moreover, we have found no references involving any heuristic methods for solving MEDKC, which motivates the effort on the present work.

used for solving Transportation-location problem with unknown number of facilities by Carlo, David, and Salvat (2017). Gamma heuristic was used for solving the domestic facility and hub location problem in Khosravi and Jokar (2017). Heuristic was tested in a real-case data-set of Turkish network and validated by CPLEX software in small cases. The remainder of the paper is organized as follows. Some of the main results about different compact routing scheme approaches, as well as other contributions related to this problem are reviewed in the next section. Genetic algorithms are described in Section 3, while some theoretical results are presented in Section 4. Experimental results obtained on various sets of instances are shown in Section 5 and the last section is the conclusion.

3. Genetic algorithms for solving MEDKC Genetic algorithms are well known population-based metaheuristics, which simulate some aspects of natural evolution. Basic concept of the GA has been introduced by Holland (1975) and in the last four decades, GA has been proven to be very successful optimization technique for solving huge classes of various optimization problems on graphs. Since the considered MEDKC problem has not yet been heuristically solved, there is no information about the behavior of any heuristic applied on MEDKC. Therefore, the decision to choose the GA heuristic as optimization method is based on an extensive computational experience in optimization on graphs, indicating that GA often produces high-quality solutions in a reasonable time. In addition, GA has been shown very robust with respect to parameter choice in reasonable bounds on quite different problems (Kratica, Stanimirović, Tošić, & Filipović, 2007; Ljubić & Raidl, 2003; Wolf & Merz, 2007). So, the more comprehensive analysis of the influence of GA control parameters to the overall searching process would be out of the scope of this paper. In the case of MEDKC, in Section 5 Experimental results, it is clearly shown that both GAs are robust with respect to the choice of control parameters on all the considered instances. GA population consists of individuals, and each individual represents one solution of the problem. In the optimization process, GA operators continuously improve the quality of individuals, with the hope that the best individual will achieve the optimal value at the end of the optimization process. GA usually implements genetic operators: crossover, mutation and selection. Crossover typically combines two individuals to produce two new offspring that eventually inherit good genetic material from both parents. The main purpose of the mutation is to diversify the genetic material, in order to avoid premature convergence to local suboptimal solutions. Mutation operator changes only small parts of a genetic code, like changing a single gene with some small probability. Selection procedure is generally used to enable better solutions to pass into the next generation. A numerical value called fitness is calculated for each individual, indicating its quality. Fitness can be considered as the ability of the individual to survive to the next generation. In most GA implementations, better fitness indicates higher probability that the individual will survive to the upcoming generations. Fitness function maps the searching space to the fitness values. It is strongly related to the objective function of the entire optimization problem: the individual with better fitness represents a solution with the objective function closer to the optimal one. The basic object in the GA is a population. The population consists of up to several hundreds of individuals. Each individual is related to one solution of the problem. Individuals are represented by a genetic code which is typically written by binary, integral or some other finite alphabet. Genes are usually represented in a way that allows easier implementation of genetic operators. In order to ensure the diversity of the genetic material, the starting population is ordinarily chosen randomly. Furthermore, there are other approaches, like usage of some other heuristics for generating the starting solutions, typically with the intention to determine some better starting solutions. In the optimization process, genetic operators are iteratively applied to the individuals, until a stopping criterion is satisfied. It is usually based on the limitation of the overall execution time, number of generations, or on their combination. In the rest of this section, two genetic algorithms for solving MEDKC are presented. These two GAs differ in almost all main aspects:

2. Literature review First theoretical results of MEDKC were introduced by Konemann et al. (2004), where the proof of the NP hardness of the problem by using the reduction from minimum vertex-dominating set is given. In the same paper, the authors proposed a polynomial time approximation algorithm that finds an approximate solution of MEDKC, which is not worse than 4·opt + 3, where opt is the value of the optimal stretch. Many papers about compact routing schemes cope with the trade-off between the size of the routing tables and the stretch. In an early work (Peleg & Upfal, 1989), the authors constructed universal compact routing schemes, which can be used on all undirected networks. Routing scheme on weighted graphs was considered in two early works (Awerbuch, Bar-Noy, Linial, & Peleg, 1990; Awerbuch & Peleg, 1990). In Thorup and Zwick (2001), a labeled compact routing scheme is presented, which uses routing tables of size O (1/ k ) -bit at each vertex and the (4k−5) -stretch. This result was later improved by Chechik (2013). By using the Pivot Interval Routing(PIR) strategy, Eilam, Gavoille, and Peleg (2003) constructed routing tables of the size O ( n log3/2 n) bits per node for the stretch factor, at most five and the average stretch factor at most three. In Cowen (2001), a 3-stretch labeled routing scheme is obtained, using local routing tables of size O (n/2/3log 4/3 n) in an arbitrary weighted undirected network. Simulations on real-life networks have indicated that stretch less than 3 can indeed be obtained using sub-linear sized routing tables Krioukov, Fall, and Brady (2007). Following this fact, Enachescu, Wang, and Goel (2008) considered “Internet-like” graphs and presented a landmarkbased routing approach, which achieves stretch values less than 3 with high probability while maintaining O (n) memory requirements. For a given weighted undirected network, a name-independent routing scheme with stretch 3 using a O ( n ) space at each node is presented by Abraham, Gavoille, Malkhi, Nisan, and Thorup (2008). In the recent work (Roditty & Tov, 2015), the authors improved previous results and presented a routing scheme with (5 + ∊) -stretch, 1 1 that uses O ( ∊ n1/3logD) memory at each vertex and O ( ∊ logD) -bit headers, with the label of each vertex of O (logn) size. This results almost matches the optimal stretch/space combination of 5-stretch and O (n1/3) – space routing tables. More detailed overview of the results on compact routing scheme would be out of scope of this paper and can be found in the last-mentioned recent study (Roditty & Tov, 2015). Summarizing the reviewed literature, one can conclude that existing results do not provide any method for solving the MEDKC in practical situations. Shortly, the results from Konemann et al. (2004) are of a theoretical importance, since the proposed approximation algorithm gives only an upper bound for any solution. Furthermore, methods developed for finding the compact routing scheme cannot be directly applied since the natures of the two problems are quite different: while in compact routing scheme problems, one deals with the trade-off between the size of the routing tables and the stretch, in MEDKC there is no memory limit for storing distance information and the core of the 284

Computers & Industrial Engineering 113 (2017) 282–293

D. Matic et al.

operator. Let parent1 and parent2 be two parent individuals. To produce two offspring, we perform the following operations. Starting from the right sides of both parents, we are looking for the pair of corresponding genes, for which is satisfied that parent1 has “1-bit” value at the first position, while parent2 has “0-bit” at the first position of the gene. The individuals then swap their genetic code of those two genes. The similar process is simultaneously performed from the left-hand sides of two parents’ genetic codes. From the left side, the algorithm searches for the pair of genes for which is satisfied that parent1 has “0bit” value at starting position, and parent2 has “1-bit” value at that position. After such a pair of genes is found, the individuals again swaps the genetic code of those two genes. This process of swapping genes from the right and from the left continues until two searches meet each other at some position in the genetic code. By this approach, the number of chosen center vertices in both offsprings is not changed, compared to the parent individuals. Therefore, both offsprings are feasible. This crossover operator is performed with the rate 0.85.

representation, crossover and mutation, while some other aspects are similar. 3.1. GA1 3.1.1. Representation and the objective function Let G = (V ,E ,ω) be an undirected connected weighted graph with n vertices and let k be the cardinality of the center set. In the first variant of GA, the genetic code of each individual consists of n genes. Each gene consists of binary digits (0/1) and has two parts: The first bit of each gene holds the information if the vertex that corresponds to the gene is in the center (in that case, the value of the bit is 1) or not (value 0). The rest of the gene holds the information of the center vertex that is assigned to the current vertex. For each non-center vertex, a list of center vertices is sorted in non-decreasing order with respect to the distance from the current vertex and a center vertex. Center vertices are assigned to themselves. The idea of introducing the sorted list of center vertices for each non-center vertex comes from the observation that non-center vertices are assigned to “closer” vertices in a general case, while assignments to more distance center vertices are rare. By handling the information to which center vertex a vertex is assigned with respect to the distances from the vertex and center vertices in the genetic code, we can better control the overall searching process, directing the search in more promising areas. We should also notice that optimal solutions usually do not involve assignments of each vertex to its nearest center-vertex. In general, it is desirable that in the second part of a gene more zeros than ones appear. Therefore, for each non-center vertex, the assignment to center vertices is based on the approach that “closer” center vertices should be given higher probability than the “far away” ones. To follow this principle, in the initial genetic code, we take care about the probabilities for generating zeros and ones: the digit “1” for the first bit of the second part of the gene is generated with the probability K / n and each following bit will get value 1 with the probability which is two times smaller than the previous one. The value K is given as an input parameter, and it is defined as follows:

K=

{

5, n ⩾ 5 n, n ⩽ 5

3.1.3. Mutation operator This GA implements so called “two level mutation with frozen bits”. This mutation operator changes a randomly selected bit in the genetic code. Since each gene consists of two parts, the algorithm differently handles changes of a bit, depending on the bit position. Recall that first bit of the gene holds the information if the corresponding vertex is chosen as center or not. Therefore, mutation of the first bit influences not only that gene (if “1-bit” mutates to “0-bit”, that center vertex becomes non-center and vice versa), but it can influence many other noncenter vertices. For example, if a vertex was assigned to that center vertex, it will be assigned to another one. Moreover, the entire list of center vertices is updated and a new vertex becomes a candidate for assignments, which can influence many other vertices. Changing bits at the second part of a gene influence the assignment of a center vertex to the corresponding vertex. Recall that if the second part of a gene modulo k is equal to i, then i-th center vertex from the sorted list is assigned to that vertex, change of the bit at the position j,j = 0,1,2…increase or decrease the overall value by 2 j (if we look from the right to the left). So, if mutation changes the bit at the position j from 1 to 0 (respectively from 0 to 1), then the vertex is now assigned to (i−2 j)-th center vertex (with respect to distance), respectively to (i + 2 j)-th center vertex. Because of these reasons, careful refinement of the mutation rates, depending on bit position is of a crucial importance for the overall optimization process. In this implementation, for the first bit position, mutation rate is equal to 0.4/ n . For the second bit position, the mutation rate is 0.1/ n and for each succeeding bit, the mutation rate is two times smaller than previous (0.1/2n,0.1/4n,0.1/8n …). During the iteration process, it may happen that all individuals have the same bit value on certain position. Such a bit is called frozen bit. Any frozen bit can be changed neither by crossover nor selection operator and as a consequence, whole search space can be significantly reduced. From the other hand, basic mutation rate is often too low to be capable to restore lost subregions. Therefore, the mutation operator increases the mutation rates for the frozen bits by the so called “frozen factor”. For the first bit, the frozen factor is equal to 2.5, which means that frozen bits at the first position are mutated with 2.5 times higher rate than basic mutation. For the rest of the gene, the frozen factor is equal to 1.5.

(4)

For each non-center vertex, the second part of its corresponding gene takes values from 0 to k−1 containing information about the assigned center vertex. If for a non-center vertex v the second part of the gene has value j, 0 ⩽ j ⩽ k−1, that means that the vertex v is assigned to the j-th center vertex from the corresponding sorted list of center vertices. Sorting the list of center nodes is performed n−k times in each generation. Although this task requires additional CPU time, it should be noticed that the number k is relatively small in tested instances (k ⩽ 50 ), so it does not significantly affect to the total CPU time. In the initialization phase, as well as during the iterative process, it is possible to happen that total number of center vertices, i.e. the total number of genes with the starting bit equal to 1 is greater or less than k. If that number is greater than k, than first k such vertices still remains center vertices, while the rests are declared as non center vertices. If the number of center vertices is less than k, than the list of center vertices is filled by new vertices starting from the beginning, until k center vertices are identified. This step can always be performed, since k is less than n. The fitness function of an individual is equal to the objective function and is calculated as follows. For each pair of vertices (u,v ) , the assigned center vertices πu and πv are determined and the value stretch (u,v ) is calculated according to the formula (2). The objective function is equal to maximal stretch over all pairs of vertices.

3.1.4. Selection operator Both variants of GA implement so-called fine grained tournament selection (FGTS) (Filipović, 2012). FGTS is an improvement of standard tournament selection that allows fractional average tournament size. If the average tournament size is Ftour , then there are two sizes of

3.1.2. Crossover operator The GA1 implements so called “pair-of-blocks change” crossover 285

Computers & Industrial Engineering 113 (2017) 282–293

D. Matic et al.

tournament: ⌊Ftour ⌋ and ⌊Ftour ⌋ + 1, where ⌊·⌋ is the floor function. In our implementation, Ftour is set to 5.4. As an illustration of this approach, let the FGTS is applied to 50 individuals, then the tournaments of the size 6 are performed 20 times, and the tournaments of the size 5 are performed 30 times. The total number of individuals participating in the tournaments is 6·20 + 5·30 = 270 and the total number of tournaments is 30 + 20 = 50 . Then the average size of tournament is 270/50 = 5.4 .

center vertices assigned to each non-center vertex. It contains n−k genes, each gene corresponding to one non-center vertex, taking value from 0 to k−1. For each non-center vertex, an array of center vertices sorted by distance in non-decreased order is created. If the value of the gene related to the vertex v is equal to j, then the j-th center vertex from the sorted array is assigned to v. Like in the case of GA1, each center vertex is assigned to itself. As it is also the case in the first variant, “distance based” assignment of the center vertices directs the search into the more promising areas, enabling that each non-center vertex can be assigned to a “closer” center. It should be noticed again that non-center vertices are not mandatory assigned to its closest center vertex. From the other hand, experiments indicate that the majority of the vertices are assigned to their “closer” center vertices. During the initialisation phase, the algorithm takes care about the probability for generating “1-bits” in each position. Similarly to the first variant, the probability that “1” is generated at the first bit position is K / n , where K is the input parameter, defined in the same way as in (4). For the second position the probability for “1” is K /2n , for the next one K /4n etc. This strategy enables good starting solutions, ensuring that non-center vertices are frequently assigned to their closer center ones. To calculate the objective function, we first need to decode the genetic code. The list of center vertices is identified from the first segment of the code and then, for each non-center one the assigned center vertex is to be chosen. From the proposed representation, one can conclude that it is possible that duplicate indices of center vertices can appear in the first part of the genetic code. We resolve this situation as follows: If there is an index that is repeating in the genetic code, then we replace it by the next index which was not previously used. If there are no such indices, then the first previous index that is not already used is taken. It should be noticed that we always can perform either first or second way of reparation, since k < n . By this approach, we ensure that exactly k vertices will be considered as center ones. As it is already mentioned, the second segment corresponds to non center vertices. For each such vertex v, the array of established center vertices is non-decreasingly sorted with respect to distances to the current vertex. In order to determine the assigned center vertex for the vertex v, we consider the gene that corresponds to v. If that gene has value j, we take the j-th center vertex from the sorted array. After the algorithm identifies the assigned center vertex for each non-center one, the objective function is simply calculated according to the formula (3).

3.1.5. Other GA1 aspects In this implementation, we address some other GA aspects, also implemented in several other papers. More detailed description of the following GA elements can be seen, for example, in Kratica (2000).

• Number of individuals in population. In this GA implementation, the total number of individuals is set to 150. • Elitistic strategy. The elite individuals preserve highly fitted genes •

•

•

of the population. The objective value of each elite individual is calculated only in the first generation. In this implementation, the number of elite individuals is 2/3 of the total number of individuals. Avoiding duplications of individuals. During the optimization process, it is possible that an individual with same genetic code appears again. This situation should be avoided, since it restricts the overall search and can lead to a premature convergence to a local suboptimum. To discard such an individual, its objective value is set to zero. By this, we ensure that the individual will not survive to the next generation, since the tournament selection operator will remove it from the population. Handling the individuals with the same objective value. If many individuals with different genetic code, but the same objective value appears, it can lead to the situation that not promising genetic code dominates in the population. As a consequence, premature convergence to suboptimal solutions can happen. To avoid this, the appearance of such individuals is limited to a constant number. In this implementation, it is set to 40. Caching. During the execution of the GA, it is possible that the same individual appears in many generations. A standard way for avoiding the multiple calculations of objective values of the same individual is to store genetic codes and reuse already calculated objective values of repeating individuals. The calculated values of the objective function are stored in a hash-queue structure, with a capacity given in advance. When the objective function is to be calculated, the algorithm, firstly looks up if the individual is stored in the cache. If this is the case, the objective value is taken from the hash-queue table, without unnecessary re-calculating. This technique is proved to be very useful in decreasing the overall GA running time (Kratica, 1999). In this GA implementation, the number of objective values stored in the cache is limited to 5000.

3.2.2. Crossover Different representations of individuals in the GA1 and GA2 clearly indicate that different types of genetic operators (crossover and mutation) will be used. Moreover, since the genetic code consists of two parts of different nature, standard one-point or two-point crossover operators cannot be directly used. Therefore, we implement so called “double one-point operator”, which simultaneously refers to both segments: for each segment, a crossover point is randomly chosen, and the genetic code is exchanged after that point. Similarly to the case of GA1, the crossover operator is performed with the rate pcross = 0.85.

3.2. GA2 In this subsection, we describe the second GA variant. Main differences between the proposed two GA variants come from different representations of individuals, which further cause different ways of calculating the objective function, as well as constructing different mutation and crossover operators. In the Sections 3.2.1–3.2.3 we describe specific elements of the GA2. Elements of GA2 which are similar to GA1 are not described in details, and they are only mentioned at the end of this subsection.

3.2.3. Mutation Since the genetic code consists of two segments of different nature, mutation rate differs for the first and the second part of the code. Obviously, mutation performed on the first part of the code, influences on changing the set of center nodes. Mutation of the bits of the second part can have a different impact, depending on the bit position in the gene. Similarly to the case of GA1, if a bit at the position j,j = 0,1,2…is changed, then the overall value is decreased by 2 j . This makes similar effect that is described earlier in the case of GA1, i.e. if the bit at the position j is changed from 1 to 0 (respectively from 0 to 1), then the

3.2.1. Representation and the objective function While in the first variant, each individual consists of n genes (n is the number of vertices), each related to one vertex, in the second variant the genetic code of each individual is divided into two segments: the first segment is a string of the length k, where k is the number of center vertices. Each part of the string, taking values from 0 to k−1 is related to one center vertex. In the second segment, we store information of 286

Computers & Industrial Engineering 113 (2017) 282–293

D. Matic et al.

vertex is now assigned to (i−2 j)-th center vertex (with respect to distance), respectively to (i + 2 j)-th center vertex. Therefore, it is reasonable to consider different mutation rates for the different bit position in the gene. For the first segment of the genetic code, mutation rate is set to 0.7/ n , where n is the number of vertices. For the second segment, mutation rate for the first bit of each gene is 0.1/ n and the following bits are mutated in a rate two times smaller than the previous bit (0.05/ n,0.025/ n …). This mutation operator also addresses frozen bits, enabling higher mutation rates for bit positions having the same value over all individuals. For the first segment, frozen bits are mutated with the frozen factor 2.5, i.e. the basic mutation rate is increased by the value 2.5. The frozen factor for the second segment of the genetic code is 1.5.

V (Wn □Pm) = {(i,vj )| i ∈ {1,2,…,n + 1}, j ∈ {1,2,…,m}} namely and 2mn + (m−1)(n + 1) edges. Two vertices (u,v ) and (u′,v′) are adjacent in Wn □Pm if u = u′ and v and v′ are adjacent vertices in Pm , or v = v′ and u and u′ are adjacent in the wheel Wn . It is known that the Cartesian product of any graph G and a path Pm makes m copies of the graph G, where corresponding vertices of each two “adjacent” copies of G are also adjacent. In our case, the Cartesian product makes m copies Wn1,Wn2,…,Wnm of the wheel Wn , where corresponding vertices of each two “adjacent” wheels are adjacent. It can easily be concluded that if two non-center adjacent vertices u and v are assigned to different centers πu and πv , then d (u,π ) + d (π ,π ) + d (π ,v ) stretch (u,v ) = ω u dω (uu,v)v ω v ⩾ 3, since any of the three ω summands in the numerator is at least 1 and dω (u,v ) = 1. Therefore, we can state the following property.

3.2.4. Other GA2 aspects GA2 also implements FGTS selection operator, already described in the SubSection 3.1.4. Further elements are as well similar to the GA1, described in the SubSection 3.1.5: population size, elitist strategy, avoiding duplications of individuals, handling the individuals of the same objective value and caching.

Property 2. Let G be a graph. If in the optimal solution of the MEDKC there exists two adjacent non-center vertices that are not assigned to the same center, then optMEDKC (G ) ⩾ 3. Let us now consider the graph Wn □Pm and k = m with vertices and edges denoted as it is explained above. We define the center set Π to contain all center vertices of each copy of the wheel Wn . Formally, Π = {(u1,vj )|j = 1,…,m} . Each vertex from the wheel Wnj we assign to the center vertex (u1,vj ) . To calculate stretches for each pair (u,v ) and (u′,v′) , we consider the following cases:

4. Notes on optimality of results Since a heuristic method cannot guarantee the optimality of the obtained solution, we decided to investigate some classes of graphs for which optimal solution can be determined. We determine the optimal solutions for two classes of unweighted graphs: wheel graphs and Cartesian product of the wheel graphs and the paths. In the experiments presented in the SubSection 5.1, we tested the proposed GAs against these results that are theoretically proved as optimal.

(a) Both vertices (u,v ) and (u′,v′) belong to the same copy of the wheel, i.e. v = v′. From Property 1 we have that stretch ((u,v ),(u′,v′)) ⩽ 2 . (b) (u,v ) ∈ Wnj and (u′,v′) ∈ Wnk and without loss of generality suppose that j < k . Then we analyze the following subcases: i. Both vertices (u,v ) and (u′,v′) are centers of the wheels Wnj and k−j Wnk . Then stretch ((u,v ),(u′,v′)) = k − j = 1

4.1. Optimal solution for the wheel graphs

ii. (u,v ) in on the cycle of Wnj and (u′,v′) is the center of the Wnk . 1+k−j Then stretch ((u,v ),(u′,v′)) = 1 + k − j = 1. Similar case appears if

A wheel graph Wn with n + 1, n + 1 ⩾ 4 vertices and 2n edges is a graph that consists of a cycle of the length n and one additional vertex that is adjacent to all others. Mathematically, we denote the set of V (Wn ) = {1,2,…,n + 1} vertices and the set of edges E (Wn ) = {(1,i),(i,i + 1)|i = 2,…,n} ∪ {(1,n + 1),(n + 1,2)} . The vertex 1 is called the center vertex of the wheel. It is obvious that in any unweighted graph, having at least two adjacent non-center vertices, it holds that the stretch is greater or equal 2. So, the minimal stretch of the such entire graph also cannot be less than 2. To reach this lower bound for a wheel Wn, n ⩾ 3 and for any k ⩾ 1, it is enough to declare the vertex 1 of the wheel as a center vertex of the solution of MEDKC and assigns all other vertices to that vertex. For each pair of non-adjacent vertices on the cycle, the stretch is then equal to one, since the shortest path contains the center vertex. If two vertices on the cycle are adjacent, then the stretch is equal to 2, since the distance between them is 1, while the distance over the center vertices is 2. The stretch between the center vertex and any vertex from the cycle is 1. Therefore, we have the following property:

the first vertex is the center of the Wnj and the second one is on the cycle of the Wnk . iii. Neither (u,v ) nor (u′,v′) are center vertices in their wheels. If they are corresponding vertices (u = u′), then we have 1+k−j+1 ⩽ 3, since k−j ⩾ 1. stretch ((u,v ),(u′,v′)) = k−j If (u′,v′) is adjacent to the corresponding vertex of (u,v ) in the cycle of Wnk , then 1+k−j+1 stretch ((u,v ),(u′,v′)) = k − j + 1 , which is less than 3. In the case when (u′,v′) is not adjacent to the corresponding vertex of the (u,v ) in the cycle of Wnk , the stretch is even equal to 1, since one shortest path from (u,v ) to (u′,v′) goes through the center vertices of the wheels. We can see that in any case, the stretch is less or equal to 3. Clearly, the condition of Property 2 is also satisfied, so we can state the following property. Property 3. Let G optMEDKC (Wn □Pm) = 3.

Property 1. Let G be a “wheel” graph Wn, n ⩾ 3 and 1 ⩽ k < n . Then optMEDKC (Wn ) = 2 .

be

the

Cartesian

product Wn □Pm .

Then,

5. Experimental results 4.2. Optimal solution for Cartesian product of the wheels and the paths This section contains experimental results, obtained by both GA implementations. All the tests are executed on the Intel i7-4770 CPU @ 3.40 GHz with 8 GB RAM and Windows 7 operating system. For each execution, only one thread/processor is used. The GAs are implemented in C programming language and compiled with Visual Studio 2013 compiler. For each problem instance, 20 independent runs were performed. Each instance is tested for different values of k ∈ {1,2,3,4,5,10,15,20} . Both GA stop after 5000 generations, or after

We further consider optimal solutions on Cartesian product of the wheel Wn and the path Pm,n ⩾ 3,m > 1 and k = m , where k is the cardinality of the center set. Let Wn be a wheel graph with n + 1 vertices {1,2,3,…,n + 1} , where the vertex 1 is adjacent to all other vertices, while vertices {2,3,…,n + 1} form a cycle of the length n. Let Pm be a path with m vertices {v1,v2,…,vm} . The Cartesian product Wn □Pm is the graph with (n + 1) m vertices, 287

Computers & Industrial Engineering 113 (2017) 282–293

D. Matic et al.

2000 generations without changing the best obtained result. To our knowledge, this problem has not yet been solved by heuristic methods, so we could not compare performances of our GAs to any other method. Because of that reason, there are no public instances specially designated for MEDKC. Therefore, we decided to take some well known publicly available graph instances used for other graph problems. Although the GAs were originally developed for weighted graphs, they can also be applied to unweighted graphs, (each edge weight is assigned to 1). For this purpose, we test the GAs on two classes of graphs for which we determined optimal solutions: wheel graphs and Cartesian product of wheels and paths. Additionally, tests are performed on six random unweighted graphs, named CURRO, created by Currò (2014), originally used for solving the Roman domination problem. The results obtained on wheel instances, as well as the results on Cartesian product of wheels and paths are presented in the SubSection 5.1, while the results obtained on CURRO instances are presented in the SubSection 5.2. The second part of the experiments are performed on weighted graph instances. Like in the case of unweighted graphs, we also decided to use some sets of well known instances from the literature. For this purpose, we took two subsets of weighted instances from Blum and Blesa (2005), used for solving the edge-weighted k-cardinality tree problem. Since each instance is tested for eight values of the parameter k, we decided to take one representative per each instance class. The first set of weighted instances, named BX instances, consists of seven 4regular graphs of different sizes. This set is a subset of total of 35 4regular graph instances provided in Blesa, Blesa, Xhafa, and Girona (2000). The second set of instances, named BB instances is adopted in Blum and Blesa (2005). From the total of twenty available instances, for our purpose we used seven of them. More detailed description of the instances and the results obtained on both sets of weighted instances are shown in the SubSection 5.3.

W50 □P4 , for which GA1 found the result 5.

• GA2 is more successful, since it reaches all optimal results. • For each instance, the average execution time for GA1 is less 28s. •

From the presented results, we can conclude that the GA2 is slightly more successful than GA1, since it reaches all optimal solutions (one more than GA1), even in less time. 5.2. Experimental results on random unweighted instances Beside the theoretical instances, we tested our GAs on the random instances from the literature. From a large set of total of 288 CURRO random instances, we took six of them: three instances with 50 vertices and three instances with 100 vertices, with different density. The size and the density level can be concluded from the name of the instance: if the name is Rand-x-y, then x is the number of vertices and y is the edge’s probability (in percentage) of existence. In our case, x ∈ {50,100} and y ∈ {1,10,50} . For example, the instance Rand-50-10 contains 50 vertices and each edge is presented in the graph with the probability 10%. This percentage indicates that the graph will contain about 112 edges, and it is exactly the case with the instance Rand-50-10. Table 1 contains the experimental results obtained on the CURRO instances. In the first three columns of the table, the name of the instance, the number of vertices of the instance and the cardinality of the center set are given. The next four columns are related to the GA1 and contain: the best obtained result (column best), the average result obtained on 20 runs (column avg), average total execution time in seconds (column total_t[s]) and the time needed to obtain the best result. The last four columns contain data obtained by the GA2 and they are organized in a similar way as for the GA1. From Table 1 one can see that both GAs achieve quality results in a reasonable time, up to 20 s. Both GAs achieve similar results, except for the instance Rand-100-1 and k = 15, where GA1 archives the result 7, while the GA2 achieves 9. The total enumeration could find optimal solutions for all instances for k = 1, since the MEDKC is not NP hard problem for that value of k. Both GAs reach those optimal solutions. With respect to the average solutions, one can conclude that GA1 slightly outperforms GA2 on this class of instances. Comparing total execution times of two GAs, from Table 1 one can calculate that GA2 is slightly faster than GA1: in 30 out of 48 cases, GA2 finds a solution in less time than GA1. Results obtained by both GAs also confirm an assumption that, in case of random graphs, the stretch inversely depends on graph density, i.e. the higher density, the less stretch. In addition, from Table 1 one can see that execution times do not depend on the density. Both conclusions are of a great importance for the usability and reliability of the proposed methods: increased number of connections between nodes in a network will decrease the maximal stretch, while the GA execution time remains nearly the same. For a deeper analysis, we first consider the following observation: For a given instance and given two cardinalities of the center set k1 and k2 , where k1 < k2 , the optimal solution obtained for the cardinality k1 of the center set is greater or equal to the optimal value for the case when the center set is of the cardinality k2 . In order words, in the case of k2 center vertices, we can always take the same k1 center vertices and the same assignments and get the same value. By using additional k2−k1 center vertices, we can just get even better solution. We also expect that the results obtained by GAs should follow this consideration, i.e. the best values for the same instance should decrease as k increases. From Table 1 we can see that this is the case for all instances, with the only one exception of the largest instance and k = 20 . It should be noted that

5.1. Experimental results on the wheel graphs and cartesian product of wheels and paths In Property 1 we proved that the optimal stretch for a wheel graph is two. The first set of experiments were performed to investigate if the proposed GAs can reach the optimal results for the wheels. We tested both variants of GA on 5 wheel instances Wi , where i ∈ {10,20,30,40,50} . The value k, i.e. the cardinality of the center set takes each value from the set {1,2,3,4,5,10,15,20} , with the natural condition that k must be less than i. The overall results on wheel instances can be summarized by the following findings:

• Both GAs reach the optimal solution (which is equal to 2), for each wheel graph W , i ∈ {10,20,30,40,50} and for each tested value of k. • For each instance, the average execution time for GA1 is less than i

•

The largest average execution time is obtained for the largest instance W50 □P5 . GA2 is slightly faster than GA1 on this class of instances. The average execution time for GA2 is less than 25 s, even for the largest instance.

4 s. As expected, the largest average execution time is for the largest instance W50 and k = 20 . Execution time for GA2 is also less than 4 s for all instances; Similarly to the case of GA1, the largest average execution time is for the largest instance W50 and k = 20 .

In the second experiment, we tested the performances of the proposed GA on the Cartesian product of wheels and paths. For this purpose, we generated total of 16 graphs denoted as Wi □Pj , where i ∈ {10,20,30,50},j ∈ {2,3,4,5} and tested both GA variants on them. The cardinality k of the center set is equal to the number of vertices in the path, i.e. for each case, k = j . In Property 3 we stated that the optimal solution for such graphs is equal to 3. Findings obtained by testing the GAs on this class of graphs can be summarized as:

• GA1 reaches all optimal results (equal to 3), except for the instance 288

Computers & Industrial Engineering 113 (2017) 282–293

D. Matic et al.

convenient to overcome this unwanted situation.

this observation also holds for weighted graphs. In several cases of unweighted instances, a so called plateaux problem appears (Davidović, Hansen, & Mladenović, 2005). The plateaux problem is the situation when the algorithm needs to manipulate with many individuals with the same objective value. These situation are generally not welcome, since they usually lead to the premature convergence in suboptimal solutions. In the case of MEDKC and unweighted graphs, this situation appears when the obtained stretch is relatively small. To resolve the problem, we slightly modified the fitness function and involved the additional consideration: if two individuals ind1 and ind2 has the same objective value (same stretch), then the fitness of the individual ind1 is considered better than of the individual ind2 if the number of pairs of vertices of the maximal stretch in ind1 is less than the corresponding number of the individual ind2 . By this approach, we enable the selection operator to discard solutions of the same stretch which cannot easily be “corrected” to the ones with smaller stretch. Experimental results confirmed that this approach is

5.3. Experimental results on weighted instances The first set of weighted instances is the set of BX instances, that consists of seven 4-regular graphs of different sizes, where the number of vertices is from the set {25,50,75,100,200,400,1000} . All these instances are chosen from a larger set of 4-regular graphs, provided by Blesa et al. (2000). The results of the proposed GAs on these instances are shown in Table 2. Columns of the table are organized in a similar way as it is the case of random unweighted instances. The name of the instance is the original name extracted from Blesa et al. (2000). From Table 2 one can see that both GAs obtain results in a reasonable time, even for the largest instances with 1000 vertices. From the total of 56 cases, both GAs achieve the same result in 28 cases. For the rest of instances, GA1 is better in 27 cases, while GA2 is better in only one case (instance g100-4-1 and k = 20 ).

Table 1 Experimental results on CURRO instances. instance

|V |

k

GA1 best

avg

GA2

total_t [s]

time [s]

best

avg

total_t [s]

time [s]

Rand-50-1

50

1 2 3 4 5 10 15 20

35 23 17 15 11 7 5 3

35.00 23.20 17.70 15.30 12.80 7.2 5 3.1

0.170 0.514 0.762 0.869 1.134 2.047 2.907 4.654

0.001 0.054 0.159 0.191 0.262 0.473 0.646 1.200

35 23 17 15 11 7 5 3

35.00 24.60 18.70 15.60 13.90 7.30 5.10 3.10

0.111 0.389 0.752 0.913 1.142 2.775 5.615 9.092

0.001 0.010 0.137 0.156 0.110 0.501 1.344 3.355

Rand-50-10

50

1 2 3 4 5 10 15 20

6 6 6 6 6 5 5 4

6.00 6.00 6.00 6.00 6.00 5.00 5.00 4.80

0.204 0.589 0.912 1.031 1.262 3.406 3.675 4.495

0.001 0.038 0.132 0.157 0.194 1.241 0.702 0.726

6 6 6 6 6 5 5 4

6.00 6.05 6.00 6.00 6.00 5.15 5 4.85

0.119 0.553 0.972 1.003 1.376 4.071 6.205 7.786

0.000 0.032 0.146 0.075 0.158 1.108 1.434 1.608

Rand-50-50

50

1 2 3 4 5 10 15 20

4 4 3 3 3 3 3 3

4.00 4.00 3.2 3.00 3.00 3.70 3.45 3.35

0.220 0.770 1.221 1.819 2.314 3.242 5.032 6.285

0.001 0.189 0.401 0.908 1.237 1.433 3.367 4.482

4 4 3 3 3 3 3 3

4.00 4.00 3.40 3.00 3.05 3.45 3.55 3.50

0.125 0.371 0.784 0.964 1.426 3.681 6.433 9.198

0.001 0.033 0.080 0.283 0.402 1.297 2.721 5.051

Rand-100-1

100

1 2 3 4 5 10 15 20

55 37 27 23 19 11 7 7

55.00 37.00 30.40 24.20 20.50 12.20 8.90 7.10

0.296 1.967 3.351 3.724 4.926 7.772 11.088 15.694

0.002 0.271 0.763 1.036 1.728 2.837 3.833 6.400

55 37 27 23 19 11 9 7

55.00 38.20 31.10 25.60 21.00 11.60 9.10 7.10

0.132 1.222 2.130 2.595 3.234 7.526 12.455 18.049

0.002 0.187 0.261 0.463 0.686 2.263 3.834 5.354

Rand-100-10

100

1 2 3 4 5 10 15 20

5 5 5 5 5 5 5 5

5.00 5.85 5.35 5.45 5.30 5.10 5.00 5.00

0.306 2.097 4.354 4.507 6.703 11.175 12.599 13.160

0.002 0.386 2.240 2.338 4.293 8.167 6.527 5.303

5 5 5 5 5 5 5 5

5 5.45 5.80 5.50 5.85 5.30 5.10 5.00

0.155 1.393 2.042 2.791 3.427 7.739 12.614 15.751

0.001 0.158 0.215 0.946 0.566 2.619 4.735 4.734

Rand-100-50

100

1 2 3 4 5 10 15 20

4 4 4 3 3 3 3 4

4.00 4.00 4.00 3.55 3.005 3.40 3.85 4.00

0.310 3.218 3.812 4.934 7.783 10.909 13.963 14.811

0.001 1.467 1.624 2.994 6.627 8.863 9.596 11.173

4 4 4 3 3 3 3 4

4.00 4.00 4.05 3.90 3.55 3.60 3.90 4.00

0.136 0.764 1.958 2.945 3.290 9.017 13.586 19.363

0.004 0.132 0.340 0.748 0.801 4.172 6.513 14.095

289

Computers & Industrial Engineering 113 (2017) 282–293

D. Matic et al.

containing seven graphs: five grid edge-weighted graphs of different sizes (15 × 15, 33 × 33, 45 × 5, 50 × 50 and 100 × 10), as well as two large Steiner tree benchmarks, named steinc5 (containing 500 vertices) and steind5 (containing 1000 vertices). Names for two Steiner tree graphs are kept from Blum and Blesa (2005). The results are shown in Table 3. Data in the table are organized in a similar way as in previous cases. From the name of a grid graph instance, one can conclude the dimensions of the grid. From Table 3, it is evident that both GAs achieve solutions in a reasonable time, even for the largest instance with 2500 vertices. The

Furthermore, on medium and large instances, GA1 outperforms GA2 with respect to average solution quality and appears to be more stable on average. With respect to the necessary runtime, the picture is different, i.e., for most instances GA2 is faster than GA1. Since the problem is not NP-hard for k = 1, we could check the optimality of the obtained results by the total enumeration. For other values of k, the total enumeration could find optimal solution only for k = 2 and the smallest instance (g25-4-1). Both GAs achieves all these optimal solutions. The second set of weighted instances is the set of BB instances, also Table 2 Experimental results obtained on BX weighted instances. instance

|V |

k

GA1 best

avg

GA2 total_t [s]

time [s]

best

avg

total_t [s]

time [s]

g25-4-01

25

1 2 3 4 5 10 15 20

25.857143 10.333333 8 5.4 5 2.6 2.111111 1.823529

25.857143 10.333333 8 5.4 5.114479 2.651429 2.111111 1.823529

0.509 0.607 0.773 1 1.112 2.023 2.341 2.289

0 0.033 0.022 0.216 0.198 0.684 0.42 0.136

25.857143 10.333333 8 5.4 5.142857 2.6 2.111111 1.823529

25.857143 10.333333 8 5.435135 5.196571 2.821022 2.117111 1.828186

0.522 0.676 0.614 0.746 0.747 1.713 2.746 3.377

0 0.006 0.011 0.056 0.04 0.406 0.556 0.192

g50-4-01

50

1 2 3 4 5 10 15 20

24.75 16 11.222222 9.588235 9 5 4.234043 3.195652

24.75 16.7 11.95235 10.623648 9.177745 5.751302 4.55421 3.655856

1.845 3.683 4.485 4.421 4.672 6.623 9.184 11.288

0.001 0.285 1.04 0.731 1.105 2.267 2.623 3.774

24.75 16 11.222222 9.588235 9 5 4.234043 3.328358

24.75 16.245238 11.911111 10.73225 9.618384 6.118029 4.531934 3.781173

0.328 0.735 0.909 1.019 1.272 4.099 7.547 10.986

0.001 0.06 0.098 0.142 0.203 1.073 2.468 3.824

g75-4-01

75

1 2 3 4 5 10 15 20

31.5 22.666667 17.75 16.2 11.666667 7.666667 5.466667 4.116279

31.5 22.666667 17.7625 16.2025 12.420449 8.489212 6.468417 4.786497

1.623 4.149 4.732 4.397 4.257 7.629 11.043 14.976

0.007 0.303 0.753 0.143 0.722 2.507 4.575 6.772

31.5 22.666667 17.75 16.2 11.666667 8.304348 5.466667 4.482143

31.5 23.214583 17.964167 16.4075 13.66903 9.260038 6.425334 4.969506

0.134 1.287 1.837 2.101 2.951 6.075 10.324 19.892

0.001 0.156 0.137 0.122 0.643 1.086 3.212 11.035

g100-4-01

100

1 2 3 4 5 10 15 20

45.666667 30.25 24.333333 13.333333 12.24 8.9375 7.375 6.052632

45.666667 30.25 24.566666 14.229982 13.3365 9.747968 8.114372 7.159488

0.725 3 4.879 4.946 5.44 10.062 15.902 21.323

0.005 0.025 0.792 0.644 0.745 3.383 6.535 9.689

45.666667 30.25 24.555556 14 13.333333 9.307692 7.545455 5.931818

45.666667 30.404167 24.977778 17.560806 14.051923 10.117313 8.262583 7.09543

0.492 1.562 2.348 2.926 3.508 9.162 16.56 25.843

0.009 0.008 0.019 0.41 0.458 1.573 4.019 8.284

g200-4-01

200

1 2 3 4 5 10 15 20

76 40.142857 31.8 27.111111 22.454545 16.6 11 10.393939

76 40.142857 31.8 28.811111 23.662093 17.532976 12.805148 11.074202

0.829 9.493 11.786 13.348 15.649 25.666 40.44 48.763

0.016 0.245 0.49 1.039 2.353 5.991 17.528 19.736

76 40.142857 31.8 29 23 17 11.181818 10.685714

76 41.125714 32.365714 29 26.813889 18.521628 14.953014 11.22531

0.162 6.945 9.298 11.477 15.804 32.072 40.52 56.671

0.014 0.333 0.283 1.393 4.445 12.229 11.752 18.098

g400-4-01

400

1 2 3 4 5 10 15 20

247 93 71.5 55 51.4 29 21 17.769231

247 93 71.5 55.375 54.64 32.121905 24.448889 20.512649

1.456 36.537 47.988 52.675 60.328 95.382 119.684 193.912

0.179 1.648 5.497 8.245 10.132 32.156 41.624 95.592

247 93 72 55.666667 55 31.5 24.6 19.4

263.3 157.1 82.908333 66.233333 65.275 38.248571 29.68627 25.071429

0.46 31.497 41.13 43.916 46.797 70.535 91.147 105.145

0.24 5.895 8.562 7.844 6.457 22.096 31.325 30.828

g1000-4-01

1000

1 2 3 4 5 10 15 20

183 138 114 90 78 55 45.666667 38.5

183 138.2 114.4 94.05 88.05 63.008333 55.4 46.197619

9.014 467.519 483.743 558.249 594.203 891.103 908.486 988.091

4.351 38.195 33.021 121.116 160.485 371.548 330.472 429.589

183 143 122 107.666667 90 70.5 57 53

183 159.9 142.05 121.35 110.85 79.366667 65.8 64.7

3.725 224.296 248.211 284.886 521.564 758.729 814.247 905.573

1.974 29.203 40.486 84.462 94.691 263.645 241.235 276.545

290

Computers & Industrial Engineering 113 (2017) 282–293

D. Matic et al.

For each instance and the case k = 1, the solutions obtained by both GAs are verified as optimal by the total enumeration. In this case (k = 1), due to the nature of the problem, obtained optimal solutions are pretty high. For other values of k, as it is already mentioned, it is expected that the optimal solutions decrease as k increases. From Tables 2 and 3 one can see that this is the case with the solutions obtained by both GAs, for all BX and BB instances. High quality of the obtained solutions clearly indicates that proposed GAs are capable for real applications of constructing routing schemes in complex computer networks. Both GAs are stable and

execution time for the largest instance is less than 7620 s. In 16 of total of 56 cases, both GAs reach the same best solution. In other cases, GA1 is again more successful than GA2: GA1 is better in 39 cases, while the GA2 is better only in one case (the instance 33x33 and k = 5). Like in the case of BX instances, GA1 also outperforms GA2 with respect to average solution quality. From the total of 56 cases, average solution of GA1 is better in 51 cases, while in the rest of 5 cases, average solutions are equal. With respect to the necessary runtime, the picture is more diverse, but still similar to the case of BX instances: for most instances, GA2 is faster than GA1. Table 3 Experimental results obtained on BB weighted instances. instance

|V |

k

GA1 best

avg

GA2 total_t [s]

time [s]

best

avg

total_t [s]

time [s]

15 × 15

225

1 2 3 4 5 10 15 20

627 403 157 134 98 50.555556 25 20.666667

627 403 159.4 142.4 100.508333 53.266667 31.906667 23.19

0.907 12.207 15.148 15.952 18.017 28.85 41.638 53.142

0.049 0.226 1.636 2.071 3.722 7.223 14.364 20.656

627 423 167 148 98 51 26.25 21

627 423 176.05 152.35 116.593333 57.5 36.463373 23.460833

0.38 6.137 9.865 11.676 21.248 29.562 36.433 52.267

0.013 0.019 0.217 0.517 7.432 10.516 11.321 22.303

33 × 33

1089

1 2 3 4 5 10 15 20

1495 1087 707 609 601 319 195 185

1495 1088 707 613.7 607 368.35 258.416667 206.15

15.172 772.058 910.835 1060.107 1012.91 1354.273 1548.094 1519.992

5.36 105.541 159.031 363.148 195.902 650.386 902.792 741.422

1495 1087 707 643 595 343 227 185

1496.2 1134.3 743 688.7 628.4 383.7 284.316667 215.15

5.413 591.818 673.265 808.344 821.428 968.63 1130.592 1149.519

2.369 32.731 83.901 208.041 148.951 222.947 383.558 348.991

45 × 5

225

1 2 3 4 5 10 15 20

1183 483 208.333333 183 146.333333 59 32 19

1183 483 208.333333 183 150.9 64.926667 37.15 25.608333

0.867 12.429 15.07 16.77 18.221 23.06 32.795 44.755

0.031 0.293 0.638 1.549 3.222 4.193 8.848 19.13

1183 483 208.333333 183 146.333333 61 34 395

1183 483 246.6 198.533333 175.033333 73.94 40.733333 473.3

0.175 10.03 14.231 16.429 17.122 25.9 40.808 6677.6

0.031 0.085 2.285 3.478 1.976 9.116 14.371 2260.595

50 × 50

2500

1 2 3 4 5 10 15 20

1975 1587 1093 927 849 583 459 389

1975 1587 1093 943.9 888.3 632.4 501.666667 444.4

195.191 3674.162 4342.583 6312.908 5936.374 6701.57 7611.624 6799.833

116.199 1043.169 1524.728 3084.247 2530.95 3235.635 4616.252 2901.624

1975 1631 1129 1093 881 641 509 395

2009.3 1713.1 1251 1106 1011.2 735.75 553.266667 473.3

47.672 3517.922 4626.977 4551.593 4598.264 5377.644 5906.175 6677.6

32.98 124.241 680.992 340.345 466.45 781.263 1510.191 2260.595

100 × 10

1000

1 2 3 4 5 10 15 20

3607 1689 1041 867 703 395 265 189

3607 1689 1041.1 869.1 718.7 405.1 297.3 228.123333

11.12 484.535 623.227 623.724 626.312 760.973 886.385 1062.399

4.767 35.178 167.014 163.167 127.257 194.575 333.968 532.296

3607 1689 1121 937 743 397 326 223

3607 1733 1146.3 1041 831.4 487.1 359 261.4

3.278 271.1 351.567 366.924 388.885 492.907 604.22 647.442

1.592 22.128 30.088 59.017 41.486 88.935 169.162 173.97

steinc5

500

1 2 3 4 5 10 15 20

393 257 203 179 169 88 52.428571 47.666667

393 257.7 204.6 182.85 170.55 96.75 62.527381 54.978571

1.979 75.483 82.895 83.491 103.021 138.411 199.472 231.482

0.424 10.194 13.131 8.911 24.268 38.433 88.057 88.144

393 261 221 186 176 99 61.5 49.571429

393 278.6 251.1 218.85 200.45 129.275 97.975 68.400794

0.634 54.861 72.982 68.237 93.12 136.838 161.846 200.101

0.148 5.967 14.068 5.656 24.627 48.709 56.565 71.507

steind5

1000

1 2 3 4 5 10 15 20

545 473 457 367 349 197 139 117

545 473 457.1 378.8 351.6 237.4 162.7 137.85

10.366 485.418 551.608 652.206 606.527 901.829 1072.29 1123.258

3.92 27.734 50.096 144.643 75.588 341.623 630.996 551.433

545 509 459 435 367 271 201 145

551.6 538.8 477.5 449.4 413.8 316.9 247.9 204.6

3.779 395.576 477.755 504.82 530.674 688.88 747.034 854.022

2.032 50.692 33.335 72.868 82.437 205.511 217.661 319.596

291

Computers & Industrial Engineering 113 (2017) 282–293

D. Matic et al.

References

accurate, while execution times still remain relatively small, even for the largest instances. From the experimental results presented in this section, it can be seen that GA1 performs better than GA2 with respect to the best and average solutions. This can be explained in the following way: GA1 uses one-segment representation, while GA2 has two-segment representation. In the integer representation of center vertices in GA2, information, whether a vertex is a center vertex or not and vertex assignments to center vertices is separated in different segments of the genetic code. On the contrary, each bit in the binary representation of center vertices used in GA1 allows the GA1 to join that bit with the index of a center vertex assigned to the current vertex. In this situation, all information about each vertex is grouped in the genetic code, implying shorter building-blocks in GA1 comparing to GA2. The Schema Theorem and the Building Block Hypothesis (Bäck, Fogel, & Michalewicz, 2000) favor shorter building-blocks and binary encoding, which could be a rational explanation for the obtained results.

Abbasi-Pooya, A., & Kashan, A. H. (2017). New mathematical models and a hybrid grouping evolution strategy algorithm for optimal helicopter routing and crew pickup and delivery. Computers & Industrial Engineering, 112, 35–56. Abraham, I., Gavoille, C., Malkhi, D., Nisan, N., & Thorup, M. (2008). Compact nameindependent routing with minimum stretch. ACM Transactions on Algorithms, 4(3), 37:1–37:12. Awerbuch, B., & Peleg, D. (1990). Sparse partitions. In Proceedings of 31st IEEE symposium on foundations of computer science (FOCS), St. Louis (pp. 503–513). Awerbuch, B., Bar-Noy, A., Linial, N., & Peleg, D. (1990). Improved routing strategies with succinct tables. Journal of Algorithms, 11(3), 307–341. Bäck, T., Fogel, D. B., & Michalewicz, Z. (2000). Evolutionary computation 1: Basic algorithms and operators, Vol. 1. CRC Press. Blesa, M. J., Blesa, M. J., Xhafa, F., & Girona, J. (2000). A C++ implementation of tabu search for k-cardinality tree problem based on generic programming and component reuse. In GCSE young researchers workshop 2000 (Part of the Second International Symposium on Generative and Component-based Software Engineering) October 9–12, Net.ObjectDaysForum (pp. 648–652). Blum, C., & Blesa, M. J. (2005). New metaheuristic approaches for the edge-weighted kcardinality tree problem. Computers & Operations Research, 32(6), 1355–1377. Carlo, H. J., David, V., & Salvat, G. (2017). Transportation-location problem with unknown number of facilities. Computers & Industrial Engineering, 112, 212–220. Chechik, S. (2013). Compact routing schemes with improved stretch. Proceedings of the 2013 ACM symposium on principles of distributed computing (pp. 33–41). New York, NY, USA: ACM. Cowen, L. J. (2001). Compact routing with minimum stretch. Journal of Algorithms, 38(1), 170–183. Currò, V. (2014). The roman domination problem on grid graphs (Ph.D. thesis). Davidović, T., Hansen, P., & Mladenović, N. (2005). Permutation-based genetic, tabu, and variable neighborhood search heuristics for multiprocessor scheduling with communication delays. Asia-Pacific Journal of Operational Research, 22(03), 297–326. Eilam, T., Gavoille, C., & Peleg, D. (2003). Compact routing schemes with low stretch factor. Journal of Algorithms, 46(2), 97–114. Enachescu, M., Wang, M., & Goel, A. (2008). Reducing maximum stretch in compact routing. INFOCOM 2008. The 27th conference on computer communications (pp. 977– 985). IEEE. Filipović, V. (2012). Fine-grained tournament selection operator in genetic algorithms. Computing and Informatics, 22(2), 143–161. Holland, J. (1975). Adaption in natural and artificial systems. University of Michigan Press. Kaboli, S. H. A., Fallahpour, A., Selvaraj, J., & Rahim, N. (2017). Long-term electrical energy consumption formulating and forecasting via optimized gene expression programming. Energy, 126, 144–164. Kaboli, S. H. A., Selvaraj, J., & Rahim, N. (2016). Long-term electric energy consumption forecasting via artificial cooperative search algorithm. Energy, 115, 857–871. Kaboli, S. H. A., Selvaraj, J., & Rahim, N. (2017). Rain-fall optimization algorithm: A population based algorithm for solving constrained optimization problems. Journal of Computational Science, 19, 31–42. Khosravi, S., & Jokar, M. R. A. (2017). Facility and hub location model based on gravity rule. Computers & Industrial Engineering, 109, 28–38. Konemann, J., Li, Y., Parekh, O., & Sinha, A. (2004). An approximation algorithm for the edge-dilation k-center problem. Operations Research Letters, 32(5), 491–495. Kratica, J. (1999). Improving performances of the genetic algorithm by caching. Computers and Artificial Intelligence, 18(3), 271–283. Kratica, J. (2000). Parallelization of genetic algorithms for solving some NP-complete problems (in Serbian). Ph.D. thesisFaculty of Mathematics, University of Belgrade. Kratica, J., Stanimirović, Z., Tošić, D., & Filipović, V. (2007). Two genetic algorithms for solving the uncapacitated single allocation p-hub median problem. European Journal of Operational Research, 182(1), 15–28. Krioukov, D., Fall, K., & Brady, A. (2007). On compact routing for the internet. ACM SIGCOMM Computer Communication Review, 37(3), 41–52. Li, Z., Kucukkoc, I., & Tang, Q. (2017). New MILP model and station-oriented ant colony optimization algorithm for balancing U-type assembly lines. Computers & Industrial Engineering, 112, 107–121. Ljubić, I., & Raidl, G. R. (2003). A memetic algorithm for minimum-cost vertex-biconnectivity augmentation of graphs. Journal of Heuristics, 9(5), 401–427. Mansouri, M., Kaboli, S. H. A., Ahmadian, J., & Selvaraj, J. (2012). A hybrid neuro-fuzzy PI speed controller for BLDC enriched with an integral steady state error eliminator. Control system, computing and engineering (ICCSCE), 2012 IEEE international conference on (pp. 234–237). IEEE. Modiri-Delshad, M., Kaboli, S. H. A., Taslimi-Renani, E., & Rahim, N. A. (2016). Backtracking search algorithm for solving economic dispatch problems with valvepoint effects and multiple fuel options. Energy, 116, 637–649. Modiri-Delshad, M., Kaboli, S. H. A., Taslimi, E., Selvaraj, J., & Rahim, N. (2013). An iterated-based optimization method for economic dispatch in power system. Clean energy and technology (CEAT), 2013 IEEE conference on (pp. 88–92). IEEE. Modiri-Delshad, M., Koohi-Kamali, S., Taslimi, E., Kaboli, S. H. A., & Rahim, N. (2013). Economic dispatch in a microgrid through an iterated-based algorithm. Clean energy and technology (CEAT), 2013 IEEE conference on (pp. 82–87). IEEE. Ng, K., Lee, C., Zhang, S., Wu, K., & Ho, W. (2017). A multiple colonies artificial bee colony algorithm for a capacitated vehicle routing problem and re-routing strategies under time-dependent traffic congestion. Computers & Industrial Engineering, 109, 151–168. Peleg, D., & Upfal, E. (1989). A trade-off between space and efficiency for routing tables. Journal of the ACM, 36(3), 510–530.

6. Conclusion In this paper, the MEDKC problem has been solved for the first time by metaheuristic approaches. Two newly developed genetic algorithms which use binary coding and modified genetic operators (crossover, mutation and selection) were proposed. Sorting the center vertices in non-decreasing order with respect to their distances from non center ones, in the process of selection directs both GAs to more promising search regions. In the first variant of genetic algorithm, genetic code consists of n genes, where n is the number of vertices. The first bit of each gene holds information if the vertex is in the set of center vertices or not, while the rest of the gene holds information about the assigned center vertex. In the GA1 the crossover operator was modified using “pair-of-blocks change” technique in order to keep the offspring individuals feasible. The new mutation operator was proposed that differently handles changes of the bits, depending on the bit position in the gene, enabling better diversification of the genetic code. In the GA2, the genetic code consists of two parts: in the first part list of center vertices is kept, while the second part contains information about assigned center vertices for non-center ones. Different representation further causes different crossover and mutation operators. The crossover operator in GA2 was modified to follow “Double-onepoint” technique in order to keep the offspring individuals feasible. The mutation operator was adapted to differently influence on the first and the second segment of the representation of individuals. Experimental results clearly indicate that both proposed methods are robust and produce high-quality results. Tests were performed both on unweighted and weighted instances. On the classes of unweighted instances for which optimal solutions are known, both GAs achieve high-quality results, reaching all optimal solutions except in the case of GA1 and one instance. Results obtained on random unweighted graphs show that both algorithms can find high-quality solutions in a reasonable time, regardless of the graph density. On weighted instances, GA1 outperforms GA2 with respect to the best solution quality and also appears to be more stable on average. Both algorithms achieve same results on a half of 56 BX instances, mainly on small ones, while in all other cases except one, GA1 is better than GA2. Results obtained on the second set give a similar picture: same best results are obtained in 16 out of 56 cases, while in all other cases except one, GA1 outperforms GA2. This research can be extended in several ways. Beside the development of the methods for solving the MEDKC, it could be useful to continue to discover theoretical results related to this problem. Further improvement directions of the proposed GA could be in parallelization of the algorithm and running on the powerful multiprocessor machines.

292

Computers & Industrial Engineering 113 (2017) 282–293

D. Matic et al.

http://dx.doi.org/10.1109/TPEL.2017.2679118. Thorup, M., & Zwick, U. (2001). Compact routing schemes. Proceedings of the thirteenth annual ACM symposium on parallel algorithms and architectures (pp. 1–10). New York, NY, USA: ACM. Wang, J., Jagannathan, A. K. R., Zuo, X., & Murray, C. C. (2017). Two-layer simulated annealing and tabu search heuristics for a vehicle routing problem with cross docks and split deliveries. Computers & Industrial Engineering, 112, 84–98. Wolf, S., & Merz, P. (2007). Evolutionary local search for the super-peer selection problem and the p-hub median problem. Hybrid Metaheuristics, 1–15.

Rafieerad, A., Bushroa, A., Nasiri-Tabrizi, B., Kaboli, S., Khanahmadi, S., Amiri, A., ... Wasa, K. (2017). Toward improved mechanical, tribological, corrosion and in-vitro bioactivity properties of mixed oxide nanotubes on Ti–6Al–7Nb implant using multiobjective PSO. Journal of the Mechanical Behavior of Biomedical Materials, 69, 1–18. Roditty, L., & Tov, R. (2015). New routing techniques and their applications. Proceedings of the 2015 ACM symposium on principles of distributed computing (pp. 23–32). New York, NY, USA: ACM. Sebtahmadi, S. S., Azad, H. B., Kaboli, S. H. A., Islam, M. D., & Mekhilef, S. (2017). A PSODQ current control scheme for performance enhancement of z-source matrix converter to drive im fed by abnormal voltage. IEEE Transactions on Power Electronics.

293