A simulated annealing based genetic local search algorithm for multi-objective multicast routing problems

This paper presents a new hybrid evolutionary algorithm to solve multi-objective multicast routing problems in telecommunication networks. The algorithm combines simulated annealing based strategies and a genetic local search, aiming at a more flexible and effective exploration and exploitation in the search space of the complex problem to find more non-dominated solutions in the Pareto Front. Due to the complex structure of the multicast tree, crossover and mutation operators have been specifically devised concerning the features and constraints in the problem. A new adaptive mutation probability based on simulated annealing is proposed in the hybrid algorithm to adaptively adjust the mutation rate according to the fitness of the new solution against the average quality of the current population during the evolution procedure. Two simulated annealing based search direction tuning strategies are applied to improve the efficiency and effectiveness of the hybrid evolutionary algorithm. Simulations have been carried out on some benchmark multi-objective multicast routing instances and a large amount of random networks with five real world objectives including cost, delay, link utilisations, average delay and delay variation in telecommunication networks. Experimental results demonstrate that both the simulated annealing based strategies and the genetic local search within the proposed multi-objective algorithm, compared with other multi-objective evolutionary algorithms, can efficiently identify high quality non-dominated solution set for multi-objective multicast routing problems and outperform other conventional multi-objective evolutionary algorithms in the literature.


The multicast routing problem (MRP)
Multicast is a telecommunication technique that simultaneously transfers information (IP datagrams) from a source to a group of destinations in communication networks.Compared to unicast, which relies on the point-to-point transmission, multicast is a more efficient solution which utilises the parallelism in networks.In this work, we consider the Multicast Routing Problem (MRP), which concerns finding the spanning tree while optimising the resource usage within the network.Due to the increasing development of numerous multicast network applications including distance learning, E-commerce and video/audio conferencing, the MRP has become one of the key problems in multimedia telecommunications and received increasing research attention in operational research.
Real world multicast applications generally have some Quality of Service (QoS) parameters or constraints and objectives.For example, an important and common QoS constraint in multicast routing applications is the bounded end-to-end delay.That is, messages must be transmitted from the source to destinations via the multicast tree within a certain limited time; otherwise most customers would cancel their requests.The efficient allocation of network resources to satisfy different QoS requirements, for example, minimising the cost of transmission via the multicast tree, is the primary goal of QoS-based multicast routing.
It is well known that the Steiner tree problem (Hwang and Richards 1992), the underlying model of MRPs, is a NP-hard combinatorial optimisation problem (Garey and Johnson 1979).It has also been proved that finding a feasible multicast tree with two independent path constraints is NP-hard (Chen and Nahrestedt 1998).The constrained Steiner tree problem under various QoS constraints is thus also NP-hard (Kompella et al. 1993).This makes the complex QoS based MRPs one of the challenging optimisation problems.Over the past decade, the problem has attracted increasing attention from the meta-heuristic research community in both computer communications and operational research (Diot et al. 1997;Yeo et al. 2004;Oliveira and Pardalos 2005).A large amount of investigations on meta-heuristic algorithms exist in the literature (Haghighat et al. 2004;Kun et al. 2005;Skorin-Kapov and Kos 2006;Zahrani et al. 2008;Qu et al. 2009).However, at the early stage, the MRPs have been mainly defined and solved as a single-objective optimisation problem subject to certain QoS constraints, i.e. to minimise the tree cost subject to a maximum end-to-end delay restriction.
With a range of inter-dependent and conflicting multiple QoS objectives and constraints (e.g.cost, delay, bandwidth, link utilisation, delay variation, packet loss ratio and hop count) in real world applications, the QoS-based MRPs can be more appropriately defined as multiobjective optimisation problems.Recent multi-objective optimisation algorithms for MRPs have been investigated concerning more realistic constraints and objectives.

Related work
A recent survey in Fabregat et al. (2005) has reviewed a variety of multi-objective multicast routing algorithms.In Table 1, we categorise meta-heuristic algorithms in the literature according to the objectives and constraints considered in problems, where a single multicast tree is constructed.It can be seen that different meta-heuristics, e.g.genetic algorithm, ant colony algorithm, artificial immune algorithm and particle swarm optimisation, have been investigated for multi-objective MRPs with various objectives.Due to the nature of multiobjective optimisation, where a set of alternative solutions is considered, it is not surprising to see that genetic algorithms, one of the mostly studied population based algorithms, have been adapted in most multi-objective multicast routing algorithms.Roy et al. (2002) adapt the widely studied multi-objective NSGA (Non-dominated Sorting based Genetic Algorithm) (Srinivas and Deb 1994) to simultaneously optimise end-toend delay, bandwidth and residential bandwidth utilisation rather than combining them into a single weighted sum objective function for wireless network routing problems.Due to the user mobility and uncertainties in wireless cellular networks, Roy and Das (2004) employ a fast and efficient QoS-based mobile multicast routing protocol based on multi-objective genetic algorithms for dynamic MRPs.In Crichigno andBaran (2004a, 2004b), two multiobjective evolutionary algorithms (MOEA) with an external population of Pareto optimal solutions have been proposed based on the strength Pareto evolutionary algorithm (Zitzler and Thiele 1999).Experimental analysis shows that MOEA1 with a binary tournament selection outperforms MOEA2 with a roulette wheel selection.Other multi-objective genetic algorithms include Koyama et al. (2004), which optimise the cost and delay of the multicast tree and Cui et al. (2003), which develop the algorithm based on Pareto dominance.
A variety of other population based meta-heuristics also appear in the multi-objective multicast routing literature.Two ant colony optimisation algorithms in Diego and Baran (2005) have shown to find more non-dominated solutions than the MOEA2 algorithm in Crichigno and Baran (2004b) on benchmark problems with different features using same computational expenses.Wang et al. (2006) propose a QoS multicast routing model based on an artificial immune system with a gene library and a clone search operator to search for better solutions.The algorithm can effectively identify a set of Pareto optimisation solutions compromising multiple QoS objectives.Particle swarm optimisation has also been investigated in Li et al. (2007) to enhance selected elite individuals before generating the next generation within a hybrid multi-objective genetic algorithm.
In the recent multi-objective optimisation research, various simulated annealing (SA) approaches (Czyzzak and Jaszkiewicz 1998;Ehrgott and Gandibleux 2000;Landa-Silva et al. 2004;Li and Landa-Silva 2008;Martins and Costa 2010;Xu and Qu 2011) have been successfully applied for different multi-objective optimisation problems.Annealing is known as a thermal process, where a solid is melted by increasing its temperature and then followed by a slow progressive temperature decrease aiming at recovering a solid state of lower energy.The SA algorithm simulates the physical annealing process to solve optimisation problems, where a solution corresponds to a state of the physical system and the fitness value of a solution corresponds to the energy of a state.It has the ability in this process to escape from local optima by visiting worse neighbouring solutions, and shows to be very effective when exploring the search space of complex multi-objective optimisation problems.Meanwhile, genetic local search algorithms (Ishibuchi and Murata 1998;Jaszkiewicz 2002;Mendoza et al. 2010) have been investigated for different multi-objective optimisation problems.Due to the ability of local search to find local optima effectively over a relatively small part of the search space, genetic local search algorithms have shown to be very suitable for solving complex multi-objective optimisation problems.Refer to Beume et al. (2007), Zhang and Li (2007), Li and Zhang (2009), Bader andZitzler (2011), Ishibuchi et al. (2011), etc., for some recent multi-objective optimisation algorithms.
To our knowledge, there is no investigation on hybridising SA with genetic local search algorithms to multi-objective MRPs.The only two recent relevant algorithms that we are aware of are applied to single objective MRPs.In Zahrani et al. (2008), a genetic local search utilises a logarithmic simulated annealing in a pre-processing step to analyse the landscape of a single objective MRP subject to multiple constraints in a group multicast scenario.Another genetic simulated annealing algorithm has been proposed in Zhang et al. (2009)   for delay jitter bounded least-cost MRP with bandwidth and delay constraints.Simulated annealing is used to compute the probability of accepting newly generated solutions.These two methods have shown to be effective in solving single objective MRPs.
In our recent work (Xu and Qu 2011), simulated annealing strategies have shown to be effective in driving a population of solutions towards the Pareto front of MRPs with four objectives.However, together with the variable neighbourhoods specially designed for MRPs, they have shown to have much less impact on the algorithm performance.In solving complex problems such as MRPs with special solution structure, specifically designed neighbourhood operators with regard to problem features have shown to be highly effective on improving algorithm performance.
In this work, motivated by the efficiency of both the simulated annealing strategies and genetic local search, we develop the first multi-objective simulated annealing based genetic local search (MOSAGLS) algorithm to solve the multi-objective MRPs.The hybrid MOSAGLS algorithm aims to combine the strengths of both Simulated Annealing and Genetic Algorithm.On the one hand, genetic algorithms have been widely used for solving multi-objective optimisation problems in the literature due to their population-based nature and the ability to simultaneously search different regions of a solution space.On the other hand, Simulated Annealing has the character of escaping from local optima by intelligently accepting worse solutions thus addressing the issue of premature convergence of GAs.In our proposed MOSAGLS, a new genetic local search with genetic operators which are specially designed for MRPs has been developed to simultaneously minimise five real life objectives, namely (1) the cost, (2) the maximum end-to-end delay, (3) the maximum link utilisation, (4) the average delay and (5) the delay variation of the multicast tree.MOSAGLS evolves by using SA-based strategies within the genetic evolutionary process to generate non-dominated solutions.A new SA-based adaptive mutation probability is also used to improve the performance of the hybrid algorithm.The impact of the SA based strategies and the local search within the genetic evolution has been investigated within this new hybrid algorithm.
The rest of the paper is organised as follows.In Sect.2, the multi-objective MRP is formally defined.Sections 3 and 4 present the proposed hybrid algorithm and evaluate its performance by experimental results.Finally, Sect. 5 concludes the paper.

The multi-objective MRP
The multi-objective optimisation problem with n decision variables, k objective functions and q restrictions can be defined as follows (Deb 2005): e(x) = e 1 (x), e 2 (x), . . ., e q (x) ≥ 0 (1) where X: the decision space of feasible regions in the solution space.
x: a vector of decision variables or a solution, x = (x 1 , x 2 , . . ., x n ) ∈ X. f i (x) (i = 1, . . ., k): objective functions with k objectives to be optimised.F (x): the image of x in the k-objective space given by the vector of k objective functions f i (x).e i (x) (i = 1, . . ., q): the set of restrictions which determines the set of feasible solutions.
Multi-objective optimisation generally concerns a set of trade-off optimal solutions, none of which can be considered superior to the others in the search space when all objectives are taken into consideration.The set of all these Pareto-optimal solutions in X is called the Pareto-optimal Set.
To model the general MRP, we denote a communication network as a directed graph G = (V , E) with |V | = n nodes and |E| = l links.The following notations are used in the rest of the paper: (i, j ) ∈ E: the link from node i to node j, i, j ∈ V .c ij ∈ R + : the cost of link (i, j ).d ij ∈ R + : the delay of link (i, j ).z ij ∈ R + : the capacity of link (i, j ) Based on the above definitions, a multi-objective MRP can then be formulated as a multiobjective optimisation problem.In this paper, we consider the multi-objective MRP with more objectives than those defined in our previous work Xu and Qu (2011) and in Crichigno and Baran (2004a).The problem is to find a multicast tree while minimising the values of the following five objectives: • The cost of the multicast tree: • The maximal end-to-end delay of the multicast tree: • The maximal link utilisation: • The average delay of the multicast tree: • Delay variation of the multicast tree: Objective (2) aims to minimise the cost occurred as the multicast tree T occupies certain required bandwidth on links in the network.Objective (3) minimises the maximal delay time of sending the data via the multicast tree so that they arrive all destinations within a shortest bounded time.Objective (4) tries to minimise the maximal link utilisation, i.e. traffic demand over the available bandwidth on the links.Objective (5) minimises the average delay time of sending the data so they arrive all |R| destinations in the shortest average time.Objective (6) minimises the delay variation of the multicast tree, which is defined as the difference between the maximum and minimum delays among all the path delays from the source to all destinations.Note that objective (3) concerns the maximal delay within the multicast tree, while objectives ( 5) and ( 6) minimise the average delay and the delay variation, respectively, thus concerning the delay to all destinations in the network.These five objectives have some correlations.For example, the delay-related objectives (2), ( 4) and (5) (Eqs.(3), ( 5) and ( 6)) which are dependent on the delays from the source to destinations in the tree are strongly correlated.The cost of the multicast tree, i.e., Objective (1) conflicts with these delay-related objective (2), ( 4) and ( 5), since the decrease of the tree cost normally brings the increment of delays.The link utilisation conflicts with the tree cost and the delay-related objectives, since the decrease of the link utilisation causes the increase of the cost and delays.As indicated by the literature, these objectives represent the most common requirements in communications.
It remains interesting future work to formulate a wider range of various objectives based on the above defined problem for different applications with specific requirements.
In communication networks, the total bandwidth of datagrams on a link must not exceed the limited bandwidth available.Hence, the total traffic on link (i, j ), i.e. the traffic demand φ of a multicast request plus the current traffic t ij is subject to the link capacity z ij : Due to the complex real world constraints in multi-objective MRPs, the search space of such problems becomes highly restricted and unpredictable (Xu and Qu 2012).This demands more efficient and effective optimisation techniques to traverse the search space of such problems with many local optimal solutions and disconnected regions of feasible solutions.

The simulated annealing based multi-objective genetic local search (MOSAGLS) algorithm
The proposed multi-objective simulated annealing based genetic local search (MOSAGLS) evolves by using simulated annealing based strategies within the genetic evolutionary process concerning non-dominated solutions with regard to the five objectives defined.Figure 1 shows the flowchart of our proposed MOSAGLS algorithm, details presented in the following subsections.
In MOSAGLS, the initial population of multicast trees is randomly generated.During the evolution, parent solutions are chosen to produce child trees by using the defined crossover and mutation operators.A local search is then applied to the generated child tree to produce a new improved tree.An external solution set NDS is maintained to record the non-dominated solutions obtained during the evolution.The MOSAGLS stops after a certain computational time, or the temperature in the SA drops to the final temperature.The NDS after the evolution is finished is output as the final results.More details of the genetic local search algorithm are given in Sect.3.2.
The proposed MOSAGLS evolves by using the SA strategies to adaptively set the mutation rate, to guide the search directions and to make decision of solution acceptance.A temperature is defined and decreased through generations.Firstly, after crossover, mutation is carried out based on an adaptive rate according to both the current temperature and the fitness of the offspring and current population.Secondly, each solution in the population is associated with a random weight vector.This vector, together with the temperature, takes Fig. 1 The flowchart of the proposed MOSAGLS algorithm part in the solution acceptance and the tuning of search directions.That is, the newly generated trees replace the selected parent based on a probability calculated using the weight vector and the current temperature.When the temperature is decreased to below a threshold, the weight vector is modified to tune the search directions.More details of the SA strategies are given in Sect.3.3.

The representation of the multicast tree
In the proposed MOSAGLS algorithm, we adopt the encoding method in Fabregat et al. (2005) to represent the solutions (multicast trees) for MRPs in both the genetic local search Fig. 2 The NSF (National Science Foundation) network.On each link, d ij , c ij and t ij denote the delay, the cost and the current traffic.The traffic demand φ = 0.2 Mbps, the capacity z ij = 1.5 Mbps.The source node s = 5, the destinations R = {0, 4, 9, 10, 13} Fig. 3 An example multicast tree and its representation for the NSF network in Fig. 2 and SA process.In this simple yet effective representation, a multicast tree is represented by an ordered set of |R| paths from the source node s to each destination r d ∈ R, |R| is the group size.That is, each solution contains |R| components {g 1 , g 2 , . . ., g |R| }, where g i represents a path between the source node s and the d-th destination node r d , d = 1, . . .|R|.
Given the benchmark NSF (National Science Foundation) network (Cui et al. 2003) in Fig. 2, an example multicast tree and its representation in MOSAGLS are shown in Fig. 3.The NSF network is a major part of the early 1990s Internet backbone for mainly academic uses.It has been tested as a benchmark problem in the existing literature by a number of researchers (Crichigno andBaran 2004a, 2004b;Xu and Qu 2011).

The genetic local search in MOSAGLS
Genetic algorithms represent one of the mostly investigated evolutionary algorithms in the literature.It simulates the evolutionary process of the nature to evolve from a population of individuals by using genetic operators (Goldberg 1989).Better individuals of higher fitness have more chance to evolve individuals which inherit good building blocks.
In the evolutionary process of our proposed MOSAGLS, the initial population consists of a fixed number of random multicast trees.They are generated by starting from the source node and randomly selecting the next connected node until all the destination nodes have been added to the tree.In each generation, crossover and mutation operations are carried out on two randomly selected parents from the current population.A local search is used to further explore better neighbouring solutions of the generated tree.A non-dominated set NDS is maintained during the evolution.It stores the newly generated tree if it is not dominated by any tree in the current NDS, and removes the trees which are dominated by the generated tree.
For the multi-objective MRPs being concerned, the strength Pareto based evaluation in Zitzler and Thiele (1999) is adopted in the genetic algorithm in MOSAGLS to calculate the fitness of individuals.It is used to maintain and update the NDS set as well as used in the selection, crossover and mutation operations.

The strength Pareto based evaluation
The value of the five objective functions (2)-( 6) defined in Sect. 2 is calculated for each individual.Based on these values, the fitness of each individual is evaluated by using the evaluation method of the Strength Pareto EA in Zitzler and Thiele (1999) as follows: (1) For each non-dominated solution T i ∈ NDS, a strength q i ∈ [0, 1] is calculated.It is the proportion of the number of solutions T j which are dominated by T i to the population size, i.e.T j is dominated by T i , denoted by T i T j : (2) For each individual T j in the population, the strength q j ∈ [1, 1 + |NDS|] is calculated by summing the strength of all non-dominated solutions T i ∈ NDS, where T i T j , plus one: (3) Finally, the fitness of each individual T j in the population F (T j ) is calculated as the inverse of its strength q j : The strength Pareto based evaluation in our algorithm is similar to that is used in MOEA algorithms in Crichigno andBaran (2004a, 2004b).In this work, we focus on the investigation of genetic local search with SA strategies, so this simple and effective strength Pareto based evaluation method in the literature has been adopted.It also enables us to carry out fair comparisons to evaluate the performance of our proposed hybrid algorithm against the MOEA algorithms in the experiments in Sect. 4.

The selection method
We employ the binary tournament selection method, also used in the MOEA algorithm in Crichigno and Baran (2004a), to select parents.Each time two individuals from the population are randomly selected.The individual with a higher fitness value defined by the strength Pareto based evaluation in (10) wins the tournament and is selected as a parent.Two parents are chosen by applying the tournament selection twice.

The crossover operation
A two-point crossover operator, with a crossover rate of 1, is applied to each selected pair of parents.Based on the representation of the ordered set of paths in Sect.3.1, the paths between two randomly generated points in one parent are selected and replaced by the corresponding paths in the other parent.Note that some selected paths may share the same links with some remaining paths.Such links will not be removed to ensure the tree is connected.To avoid loops in generating the new multicast tree, the selected path will be replaced by adding the new path from its destination node until it connects to an on-tree node.
An example of the crossover operation is shown in Fig. 4, where path g 3 = (5-4-10-12-8-7-13-9) between the selected two crossover points in parent 1 is replaced by the new corresponding path (5-6-1-0-3-10-11-9) of g 3 in parent 2. To avoid loops in the generated multicast tree, the new path (5-6-1-0-3-10-11-9) is added to the tree by starting from the destination node 9, and adding only the path 9-11-10 until it connects to the on-tree node 10.A new tree is then generated as shown in Fig. 4(c).
The simple representation of multicast trees (see Sect. 3.1) facilitates an easy implementation of crossover operations.By adding the selected path(s) in parent 2 from the destination node, the newly generated offspring is guaranteed to be feasible (if the link capacity constraint ( 7) is satisfied).Note that while g 3 = (5-4-10-12-8-7-13-9) in parent 1 is being replaced, the links along the path (5-4-10-12-8-7-13) still remains in the multicast tree as they also appear in g 5 in the original tree of parent 1.In multicast trees, some paths share common links, especially those near the root of the tree, i.e. links near the source node appear multiple times in the solution.This is due to the nature of the multicast tree that all Fig. 5 An example of the mutation operation paths must pass a subset of the d s links from the source node, where d s is the degree of the source node.These (partial) paths in the tree are adaptively selected through the evolution process of MOSAGLS if they present to be good building blocks in the individuals.

The mutation operation
In our MOSAGLS algorithm, mutation is carried out with an adaptive probability by using both the evolutionary information and SA strategies.The adaptive probability is calculated based on not only the fitness of the individual, but also the current temperature.More details of the SA strategies are given in Sect.3.3.1.
For a selected individual, the mutation operation randomly replaces a path by using an alternative path stored in a routing table, which is the same as that is devised in Crichigno and Baran (2004a).The routing table for the destination r d ∈ R of the selected path g d is consists of m least cost, m least delay and m least used paths (least utilisation path) generated by using the k-shortest path algorithm (Eppstein 1998).As objectives (2), (3) and ( 6) defined in Sect. 2 are all related to the delay objective, they share the same least delay paths in the routing table.A new randomly selected path p ∈ {path 1 , . . ., path 3m } in the routing table then replaces the original path g d from the source to the destination r d .

The local search in MOSAGLS
Instead of reproducing the offspring directly to the next generation, a local search is applied to further enhance the offspring, simulating the maturing phenomenon in the nature.To apply the local search, each solution (a multicast tree) is firstly represented by a binary array of |V | = n bits, each corresponding to a node in the multicast tree.Each bit is assigned a value 1 if the corresponding node is on the tree; 0 otherwise.This representation has been widely used in the literature and shows to be effective for local search and genetic operations for MRPs (Skorin-Kapov andKos 2003, 2006).
In our MOSAGLS, the neighbourhood operator of the local search is based on the wellknown Prim's minimum spanning tree algorithm (Betsekas and Gallager 1992) which finds a tree with the minimal total weights of the links spanning a subset of nodes in the graph.The local search repeatedly flips a bit in the binary array which represents a solution until a new better multicast tree is found or a fixed maximum number of nodes have been flipped.This local search method can greatly improve the solutions with regard to the five objectives by using the strength Pareto based evaluation (10).After the local search is applied to the above solution in Fig. 5(b), a new solution is shown in Fig. 6.
This node-based local search has been applied in our previous work (Qu et al. 2009).The selection of the neighbourhood operator is based on our previous observation that the node-based neighbourhood operator is easy to implement and effective for searching better neighbourhood solutions.In this paper, we just investigate this simple yet effective local search operator.More efficient and effective local search methods and the choice of starting solutions for local search (Ishibuchi et al. 2010) may be investigated to reduce the computational time of the hybrid algorithm in our future work.
Based on the procedure described above, Fig. 7 presents the pseudo-code of the hybrid MOSAGLS.In order to improve the performance of MOSAGLS, as shown in Fig. 7, several SA strategies have been applied.We illustrate these strategies in the following subsections.

Simulated annealing strategies in MOSAGLS
Simulated annealing (Kirkpatrick et al. 1983) is one of the mostly studied probabilistic meta-heuristics for global optimisation.The basic idea is inspired from the physical annealing process where the heated material is gradually cooled to reduce the defects and form large size crystals.Based on this physical process, during the search the SA algorithm accepts a non-improved state/solution with a probability depending on the temperature in the cooling schedule at that time.As a result, it is able to escape from local optima by intelligently accepting worse solutions and effectively explore the search space of complex multi-objective optimisation problems (Czyzzak and Jaszkiewicz 1998;Li and Landa-Silva 2008).
In our MOSAGLS, SA strategies have been adopted to (1) set the adaptive mutation rate, (2) make decision of the acceptance of the new offspring in the genetic local search, and (3) guide the search directions at the later stage of the evolution.Notes: Maxgen is the maximum of generation; t max /t min : the starting/final temperatures; t c : the temperature threshold for tuning the weight vector λ; t step : the temperature decrement; P : the current population; NDS: the non-dominated solution set Fig. 7 The pseudo-code of the MOSAGLS algorithm

Simulated annealing based adaptive mutation probability
In MOSAGLS, mutation is always applied to an offspring if it is better than the average of the current population (see more details of the mutation operation in Sect.3.2.4).For a worse new offspring, mutation is applied with an adaptive probability p m based on not only the current temperature but also the difference between the fitness of the new offspring and the average fitness of the current population.The adaptive mutation probability p m is defined as follows: where the average fitness of the current population F avg and the fitness of the new offspring F new are calculated by using the strength Pareto based evaluation function defined in (10); t : the current temperature.On the one hand, better offspring near the Pareto front will be given higher chance to be evolved by the mutation operation, making the evolution more effective by investing computational time on more promising offspring.On the other hand, the mutation rate is higher if the current temperature t is higher.This means at the early stage of the evolution, worse offspring still has the chance to be mutated, giving the whole population the chance to explore more areas of the search space.At the later stage when the mutation rate decreases to a small value along with the temperature decrement, worse solutions are rarely evolved to encourage the convergence of the algorithm.

Simulated annealing based acceptance probability
To assist the decision making of accepting a new generated offspring, and to tune the search directions in the later stage (see Sect. 3.3.3),a weighted sum of a convex combination of different objectives is used.A Pareto-optimal solution is obtained if it is the unique global minimum of the following scalar optimisation problem: where λ is a weight vector, λ i ∈ [0, 1], i = 1, . . ., k, k is the number of objectives and k i=1 λ i = 1; each λ i is associated with an objective function f i (x).Equation ( 12) is called the weighted scalarising function, one of the mostly used scalarising functions in multi-objective algorithms in the literature (see more details of other functions in Miettinen 1999).In our MOSAGLS, the probability of a new solution x to replace a solution x is adopted based on the weighted scalarising function g (ws) in Eq. ( 12): where t > 0 is the current temperature; g(x, x , λ) = g (ws) (x , λ) − g (ws) (x, λ) is the difference between x and x by using the weighted scalarising function in (12).
In the initialisation of our MOSAGLS algorithm, a normalised weight vector λ i is randomly generated for each solution T i .At each generation, a new solution T i is generated after the local search upon the offspring produced from the chosen pair of parent solutions T i and T j .The parent of lower quality is then replaced by the newly generated tree T i based on the weight vector and the current temperature, i.e. parent with a higher probability by comparing P (T i , T i , λ i , t) and P (T j , T i , λ j , t) in Eq. ( 13) is replaced.

Simulated annealing based competition and search direction tuning
Simulated annealing strategies have shown to be effective to tune search directions in solving multi-objective travelling salesman problems (Li and Landa-Silva 2008).To tune the search directions in our proposed MOSAGLS, two similar adaptation strategies to Li and Landa-Silva ( 2008) have been adopted.They are designed to diversify the search directions and avoid solutions being trapped in local optima during the evolutionary process.
(1) The first strategy is defined concerning the competition between close solutions in the current population.Euclidean distance, which is widely used in the multi-objective optimisation literature, has been used in MOSAGLS to measure the distance between two solutions.This is calculated upon the differences between all the objectives values in the two solutions for the problem.After a new solution T i is generated, the closest solution T j measured by the Euclidean distance in the current population is chosen as a neighbour of T i .MOSAGLS then replaces T j with T i if T j is worse than T i by comparing their weighted scalarising function values using Eq. ( 12), i.e. if g (ws) (T i , λ j ) < g(T j , λ j ).  8 The pseudo-code of search directions tuning by using simulated annealing strategies (2) The second strategy aims to tune the search direction when the search is getting closer to the Pareto front at the later stage of the evolution.This strategy is similar to that of multiobjective SA algorithms in Czyzzak and Jaszkiewicz (1998) and Li and Landa-Silva (2008).It adaptively tunes the weight vector according to the closest non-dominated neighbouring solution in the current population when the current temperature is decreased to below a threshold.The pseudo-code of the search direction tuning is presented in Fig. 8.
The latter two strategies in Sects.3.3.2and 3.3.3on the solution acceptance and the tuning of search directions have also been investigated in the simulated annealing based multiobjective algorithm named EMOSA in Xu and Qu (2011) for solving the multi-objective MRPs, where only four objectives (i.e.Eqs. ( 2), ( 3), ( 4) and ( 5), except Eq. ( 6)) have shown effective to find better non-dominated solutions than other multi-objective evolutionary algorithms (Crichigno andBaran 2004a, 2004b) in the literature.Based on our previous work, these two SA based strategies have been adopted to in our proposed MOSAGLS hybridised with the local search operator as the post improvement optimisation strategy for solving the multi-objective MRPs with more objectives in this paper.

Simulation environment and test instances
We have carried out a large amount simulations to test our algorithms on both the benchmark and random networks with different objectives.Two variants of our multi-objective multicast routing algorithms have been evaluated in the following experiments: (1) MOSAGA, the proposed algorithm without local search, and (2) MOSAGLS, the proposed algorithm with local search.
Both of these algorithms are new in solving multi-objective MRPs.We firstly evaluate different components of these algorithms on a set of random networks with two objectives (cost and delay, see Sect. 2).Based on the observations obtained, we then compare our algorithms with the existing algorithms in the literature on two benchmark networks with different number of objectives.Finally, we demonstrate the effectiveness and efficiency of our algorithms upon a set of random networks with all the five objectives defined in Sect. 2.
All simulations have been run on a Windows XP computer with PVI 3.4 GHz, 1 GB RAM, and within a multicast routing simulator developed based on Salama's generator (Salama et al. 1997).The multicast routing simulator generates random network topologies by using the Waxman's graph generation algorithm (Waxman 1988).The nodes of the network are located within a simulated rectangle of size 4000 × 4000 km 2 .The Euclidean metric is used to determine the distance l(i, j ) between pairs of nodes i and j .Links connect nodes (u, v), with a probability where L is the maximum distance between two nodes; α and β are parameters which can be set to different values to create desired characteristics in the network.Therefore, a large value to β creates a high average degree to nodes, and a small value to α creates long connections between nodes.More details of the simulator can be found in Salama et al. (1997) and Qu et al. (2009).In all our simulations, we set α = 0.25 and β = 0.40 to generate random network topologies.Within the network, the link cost c ij is assigned a random value between (0, 100] to define the current total bandwidth reserved on the link.The link delay d ij is defined as the propagation delay depending on the length of the link.The link capacity z ij is set as 1.5 Mbps.The current traffic t ij is randomly loaded around 50 % of its total link capacity, which is set the same as that in Diego and Baran 2005.To encourage scientific comparisons, the detailed information of all the problem instances and the experimental results has been provided at http://www.cs.nott.ac.uk/~rxq/benchmarks.htm.

Parameter settings
In the evolutionary process of MOSAGA, the population size is set to 50, the number of generations is set to 500, and the crossover rate is set to 1.These are the same as those in the other two evolutionary algorithms MOEA1 (Crichigno and Baran 2004a) and MOEA2 (Crichigno and Baran 2004b) for multi-objective MRPs.Instead of using a fixed mutation rate 0.3 in MOEA1 and MOEA2, we employ an adaptive mutation rate in MOSAGA.Based on initial tests, we set m = 25 to generate the routing table of 3 × m alternative paths in the mutation operation.In the SA process of our MOSAGA, we set the initial temperature t max = 50, the final temperature t min = 5, the temperature decrement t step = 5, and the temperature threshold t c = 25.The parameters of the SA process are kept the same as those in the SA based multi-objective algorithm for solving travelling salesman problems (Li and Landa-Silva 2008).For simplicity and fair comparisons, we keep the same parameter settings for the proposed algorithms on all the tested instances in this paper unless otherwise stated.
All the parameter values have been also determined after a set of tests to find the balance between the quality of solutions and the running time.For example, to obtain a proper range of values for the population size and the starting temperature, a set of initial tests have been carried on the NSF network shown in Fig. 2 to compare the performance of MOSAGA with different parameter settings.The optimal Pareto Front (PF) of the NSF benchmark problem that composes of 16 solutions has been found by an exhaustive search method in Crichigno and Baran (2004b), where four objectives (2), (3), ( 4) and ( 5) have been considered in the algorithms.We thus considered the same four objectives in our MOSAGA algorithms in this group of tests.Table 2 presents the maximum, minimum and average number of nondominated solutions found by each variant of MOSAGA in 50 runs.The setting of population size pop = 50 and the starting temperature t max = 50 provides the best solutions and requires less computing time, is thus selected in our algorithms.

Results on random networks with two objectives
In order to visually analyse the impact of SA strategies, the adaptive mutation and the local search on the performance of the proposed algorithms, we carry out a series of experiments

The impact of simulated annealing strategies
In the first group of experiments, we compare MOSAGA with two evolutionary algorithms, MOEA1 (Crichigno and Baran 2004a) and MOEA2 (Crichigno and Baran 2004b), to identify the impact of SA strategies to our multi-objective genetic local search.For a fair comparison, the same parameter setting in MOEA1 and MOEA2 has been used in the evolution process in our MOSAGA.The non-dominated solutions found by the three algorithms after 50 runs are shown in Fig. 9, with their average computing time in Table 3.
The experimental results clearly show that the MOSAGA algorithm found a set of better non-dominated solutions compared with those from the two MOEA algorithms on 3 out of 4 problems, i.e.Fig. 9(a), (c) and (d).For the problem in Fig. 9(b), MOSAGA is able to find better non-dominated solutions except those three solutions found by MOEA2.However, for all other three problems, MOEA2 has performed the worst.With regard to the computing time, MOSAGA requires the least computing time.This demonstrates that the SA strategies in MOSAGA can guide the search direction of genetic algorithm to find better solutions in less computing time in comparison with the MOEA algorithms.

The impact of the adaptive mutation
In the second group of experiments, we test the effect of the adaptive mutation by comparing it with four fixed mutation rates, i.e.P m = 0.3, 0.6, 0.9, 1 in MOSAGA.The nondominated solutions found by MOSAGA with different mutation rates after 50 runs are shown in Fig. 10.In general, MOSAGA with the adaptive mutation rate has the best overall performance among all the other variants with fixed mutation rates, although the Pareto fronts found by variants of the algorithm interweave at certain part of the front for some networks, i.e.Fig. 10(b) and Fig. 10(d).With the fixed mutation rate, MOSAGA occasionally finds several better non-dominated solutions, but fails to obtain the majority of solutions at the Pareto front found by MOSAGA with adaptive mutation rate for all four networks.The adaptive mutation rate makes the search more effective based on both the temperature and the fitness of the new solution and the current population.

The impact of the local search
In the third group of experiments, we compare the MOSAGA and MOSAGLS algorithms to demonstrate the effect of the local search in our proposed algorithms.The local search in MOSAGLS stops after a number of nodes (set as 10 here) has been flipped.The non-dominated solutions found by both algorithms for the four random networks in 50 runs are shown in Fig. 11, clearly showing the improvement made by the local search in MOSAGLS.For the small network in Fig. 11(a), MOSAGLS found four non-dominated solutions, while MOSAGA found only two non-dominated solutions.For the same networks with larger group size and for networks of large size, it is obvious that MOSAGLA outperforms MOSAGA by finding much better non-dominated solutions.This demonstrates  the efficiency of the local search in our proposed MOSAGLS algorithm.Please note the non-dominated solutions obtained by MOSAGA in Fig. 11 are different from the results of MOSAGA in Fig. 10 due to the independent runs of experiments on the random networks.Table 4 presents the average computing time of MOSAGA and MOSAGLS on the four random networks, showing (not surprisingly) that MOSAGLS finds better solutions at the expenses of longer computing time compared with the MOSAGA algorithm.In real world Fig. 11 The non-dominated solutions found by MOSAGLS and MOSAGA on random with two objectives in 50 runs applications, the balance of solution quality and the computational time should be concerned while designing the local search in MOSAGLS.This demonstrates that our proposed local search is able to improve the search results of MOSAGLS while consuming longer computational time.

Results on the NSF network with four objectives
Based on the above observations, we compare the performance of the MOSAGLS and MOSAGA algorithms on the benchmark NSF network with that of MOEA algorithms.As mentioned in Sect.4.2, the optimal Pareto front of 16 solutions for the NSF problem has been found in Crichigno and Baran (2004b), concerning the four objectives (2), (3), (4) and (5) defined in Sect. 2. We thus consider the same four objectives in this group of experiments.
The adaptive mutation rate has been applied in MOSAGLS and MOSAGA.To ensure a fair comparison between algorithms, we have re-implemented the MOEA algorithms, and compared the results obtained by the four algorithms within 60 seconds instead of defining a fixed number of generations.Table 5 presents the maximum, minimum and average number of non-dominated solutions found by each algorithm in 50 runs., 3, 5, 6, 9, 10, 12, 17, 19, 21, 22, 23, 25, 33, 34, 37, 41, 44, 46, 47, 52, 54} 24 For this small NSF problem, our proposed MOSAGLS and MOSAGA algorithms outperform the two MOEA algorithms, i.e. always find all 16 solutions in the optimal Pareto front in 50 runs.This, together with our above observations, indicates that SA based MOSAGA and MOSAGLS with adaptive mutation rate are more effective than the two conventional MOEA algorithms for multi-objective MRPs.

Results on the NTT network with five objectives
We also test our algorithms on the NTT (Nippon Telephone Telegraph) network in Fig. 12.Two different multicast groups have been tested, as shown in Table 6.We set the current traffic t ij of each link as randomly loaded with around 50 % of its total link capacity, which is the same as that in Diego and Baran (2005).All algorithms have been evaluated with respect to all the five objectives defined in Sect. 2.
As the optimal Pareto front for the NTT network with all the five objectives defined in Sect. 2 is not known, we use a six-step procedure devised in Diego and Baran (2005) to obtain an approximation of the Pareto front for each problem.The six-step procedure in Diego and Baran (2005) is presented as follows: (1) Each of the algorithms tested is run 10 times.Based on this approximation of the optimal Pareto solution set Y PF , in order to have an insight of the quality of solutions obtained by each algorithm with regard to Y PF , the following notations are used: (1) Comparisons of MOSAGLS with MOSAGA on the NTT network Firstly, to investigate the effectiveness and efficiency of the local search operator in our proposed MOSAGLS, we compare MOSAGLS with MOSAGA in which no local search is applied.Table 7 presents the total number non-dominated solutions in Y PF obtained by MOSAGA and MOSAGLS for the NTT network in Fig. 12 with two different multicast groups in Table 6 by using the above six-step procedure.For a fair comparison, the computational time for both algorithms for each run is 60 seconds.
Results in both Table 8 and Table 9 show that MOSAGLS found more non-dominated solutions ∈ Y PF than that of MOSAGA, and thus gave a better approximation to the Pareto Front.This becomes more obvious for the NTT network with the large group size (see group 2 in Table 6) in Table 9, where MOSAGLS found 17.27 % of the non-dominated solutions in Y PF , compared with 7.27 % of solutions in Y PF found by MOSAGA.Better experimental results obtained by MOSAGLS demonstrate that the local search operator can improve the performance of MOSAGLS by finding more non-dominated Pareto optimal solutions for solving the multi-objective MRP with five objectives.
(2) Comparisons of MOSAGLS with MOEA1 and MOEA2 on the NTT network Secondly, we compare MOSAGLS with other two multi-objective evolutionary algorithms MOEA1 and MOEA2 on the NTT benchmark problem with five objectives.Each algorithm has been run 10 times on each instance.Table 10 presents the total number non-dominated solutions in Y PF obtained by the three algorithms for the two NTT instances by using the six-step procedure.For a fair and comprehensive comparison, all algorithms tested have been given the same computational time, i.e. 160 seconds and 320 seconds, respectively, for each run.
Tables 11 and 12 present the results obtained for the NTT network with two different multicast groups.Both tables show that MOSAGLS found more non-dominated solutions ∈ Y PF than that of MOEA algorithms.A measurement used in Zitzler and Thiele (1999) is used here to calculate the coverage of the non-dominated solution set obtained by each algorithm.This is represented by the ratio of the number of non-dominated solutions obtained by each algorithm to that of another algorithm.Tables 11 and 12 show that the non-dominated solutions found by MOSAGLS cover the majority of the non-dominated solutions found by MOEA algorithms.This becomes more obvious for the NTT network with larger group size, where MOSAGLS found almost all non-dominated solutions, while the MOEA algorithms found just one or no solution that belongs to MOSAGLS.We set the same running time, i.e. 320 seconds for all algorithms in each run.Table 13 presents the average number of non-dominated solutions (| ȲPF |) obtained from 10 runs of the three algorithms on each of the four types of random networks by using the six-step procedure in Sect.4.4.2.
Tables 14 and 15 present results of the three algorithms on the random networks with different characteristics.Average results from 10 runs on 10 different instances for each type of networks show that MOSAGLS significantly outperforms the two MOEA algorithms by finding more non-dominated solutions in Y PF .The non-dominated solution set obtained by MOSAGLS covers most of the non-dominated solutions found by MOEA1 and MOEA2.This again demonstrates the effectiveness of our proposed MOSAGLS algorithm on the random networks with all five objectives defined.
To summarise, the large amount of experiments on a range of multi-objective MRPs with different features demonstrate the efficiency and effectiveness of our proposed simulated annealing based multi-objective genetic local search algorithm MOSAGLS.The SA strategies improve the performance of the algorithm by finding better non-dominated solutions in less computing time compared with conventional multi-objective evolutionary algorithms.The adaptive mutation contributes a better performance than that of fixed mutation rates.The local search further improves the performance of MOSAGLS, however at a higher computational time.Comparisons on both the benchmark problems and random networks demonstrate that the proposed MOSAGLS has the best performance among variants of algorithms, and show that MOSAGLS is able to find high quality solutions for multi-objective MRPs with different features.

Conclusions
In this paper, we investigate the first simulated annealing based multi-objective genetic local search (MOSAGLS) for solving multi-objective MRPs.A new simulated annealing based adaptive mutation probability is proposed in MOSAGLS, which can adaptively adjust the mutation rate according to the fitness of the new solution against the average quality of the current population during the genetic evolution procedure.In addition, two adaptation strategies based on simulated annealing are adopted to tune the search directions in MOSAGLS.One is the competition between similar members in the current population, another one is the two-phase strategy for tuning weight vectors.Integrated with these simulated annealing strategies, the hybrid multi-objective genetic local search is able to efficiently search towards the Pareto front and diversify the population with regard to the multiple objectives in the problem.
Due to the complex structure of the multicast tree, an ordered set of paths has been used to represent solutions.Based on this simple solution representation, the crossover and mutation operators have been specifically designed concerning the network structure to facilitate efficient and effective operations and improve the objective values of the tree while satisfying the constraint.Finally, the local search further intensifies the search by exploring more promising neighbouring solutions based on a binary string solution representation.
Through a large amount of extensive simulations, we evaluate our proposed MOSAGLS on both benchmarks and a series of random networks with different objectives and different features.Experimental results show that the simulated annealing strategies significantly improve our proposed algorithm compared against another two conventional MOEA algorithms in the literature with respect to both the solution quality and computing time for the multi-objective MRPs.This demonstrates that better control of the evolution by using the simulated annealing strategies significantly improves the efficiency and effectiveness of MOSAGLS.Compared with the simulated annealing based multi-objective genetic algorithm without local search, MOSAGLS obtained significantly better results but at larger computational expenses.
The proposed MOSAGLS algorithm has shown to be an ideal approach for solving the multi-objective MRPs, and is flexible to be extended to solve other multi-objective optimisation problems.In our future work, we also intend to investigate the influence of different selection strategies and adaptive crossover in the hybrid algorithm.More efficient and effective local search methods and the choice of starting solutions for local search may be investigated to reduce the computational time of the hybrid algorithm.Other Pareto dominance methods may also be applied to further improve the search towards the Pareto front in the MOSAGLS algorithm.In addition, some recent powerful multi-objective optimisation algorithms can be investigated in our future work.

Fig. 4
Fig. 4 An example of the two-point crossover operation

Fig. 9
Fig.9The non-dominated solutions found by MOSAGA, MOEA1 and MOEA2 on random networks with two objectives in 50 runs

Fig. 10
Fig. 10 The non-dominated solutions found by MOSAGA with different mutation rates on random networks with two objectives in 50 runs (2) For each algorithm, 10 sets of non-dominated solutions Y 1 , Y 2 , . . ., Y 10 are obtained, one from each run.(3) For each algorithm, an aggregate collection of solutions Y T is obtained, i.e.Y T = 10 i=1 Y i .(4) The non-dominated solutions in Y T are stored for each algorithm, denoted by Y alg .(5) A set of solutions Y is obtained by combining all Y alg , i.e.Y = Y alg .(6) The non-dominated solutions in Y are used as an approximation of the true optimal Pareto solution set, called Y PF .

Fig. 12
Fig. 12 The Nippon Telephone Telegraph (NTT) network of Japan with 55 nodes and 142 links.Costs are shown on each link

Table 1
Multi-objective multicast routing meta-heuristics, categorised by objectives and constraints considered, ordered by the year of publication (MOEA: multi-objective , measured in Mbps.t ij ∈ R + : the current traffic of link (i, j ), measured in Mbps.s ∈ V : the source node of a multicast group.R ⊆ V − {s}: the set of destinations of a multicast group.r d ∈ R: the destinations in a multicast group.|R|: the cardinality of R, i.e. the number of destinations, also called group size.
the network topology; s: the source node; R: the destination set; k: number of objectives Generate the random initial population P with T 1 , . . ., T |P | multicast trees; For T 1 , . . ., T |P | , produce |P | distinct random weight vectors λ 1 , . . ., λ |P | , each with a uniform spread; Select non-dominated solutions with respect to multiple objectives in the population P to form NDS; Calculate the fitness value of all the solutions in NDS and P by using the strength Pareto based evaluation in Eq. (10); Maxgen ) do { Current temperature t = t max .while t ≥ t min do { Randomly select a pair of parent solutions T i and T j from P ; Crossover operation over parents T i and T j ; Mutation with adaptive probability p m to the new solution from crossover; // simulated annealing strategy Generate T i by applying local search to the new solution after mutation; // Sect.3.2.5 if T i is not dominated by parents T i or T j then Update NDS; Replace a selected parent T c by T i with the probability P (T c , T i , λ i , t); // simulated annealing strategy Re-calculate the fitness for all solutions in NDS and P Find the closest solution T k ∈ P of T i if (T k is worse than T i ) then // simulated annealing strategy T k = T i // replace the closest member T k in P t = t − t step // temperature decrement if (t < t c ) then : current temperature; μ: the constant for tuning the search direction foreach T i ∈ P do { Find the closest non-dominated solution T i of T i in the current population by using the Euclidean distance upon λ i and λ i foreach obj ∈ {1, . . ., k} do // k: number of objectives {

Table 2
Comparison of MOSAGA with different parameter settings for solving the NSF network in Fig. 2. |NDS|: the number of non-dominated solutions; Running time = 60 seconds in each run with two objectives, namely the cost and the delay as defined in Sect. 2. Four random networks have been generated with different network sizes (|V | = 50 and |V | = 100) and group sizes (number of destination nodes |R| = 20 % * |V | and |R| = 30 % * |V |).

Table 3
Average computing time of MOSAGA, MOEA1 and MOEA2 on random networks with two objectives

Table 4
Average computing time of MOSAGA and MOSAGLS on random networks with two objectives

Table 5
The number of non-dominated solutions found by different algorithms on the NSF network in Fig.2with four objectives.Computing time = 60 seconds

Table 6
Two different multicast groups used for the experiments on the NTT network in Fig.12

Table 7
Based on the six-step procedure, the total number of non-dominated solutions Y PF obtained by MOSAGLS and MOSAGA for the NTT network with two different multicast groups with respect to all five objectives

Table 8
Comparisons of MOSAGLS and MOSAGA for the NTT network with the small multicast group 1. Computational time for each run is 60 seconds.The best results are in bold

Table 9
Comparisons of MOSAGLS and MOSAGA for solving the NTT network with the large multicast group 2. Computation time for each run is 60 seconds.The best results are in bold

Table 10
Based on the six-step procedure, the total number of non-dominated solutions in Y PF obtained by MOSAGLS, MOEA1 and MOEA2 for the NTT network with two different multicast groups with respect to all five objectives

Table 11
Results of different algorithms for solving the NTT network with small multicast group 1 in Table6

Table 12
Results of different algorithms for solving the NTT network with large multicast group 2 in Table6

Table 13
The average number of non-dominated solutions on the random networks Finally, we test the effectiveness and efficiency of our algorithm upon a set of random networks of different characteristics with respect to all five objectives defined in Sect. 2. We generate 10 random graphs for each network size of |V | = 50 and |V | = 100, and with different group sizes (|R| = 20 % * |V | and |R| = 30 % * |V |).Therefore, there are in total 40 instances (10 for each of the four types of networks) tested.

Table 14
Comparison of different algorithms on random networks of |V | = 50 with different group sizes.Computational time for each run is 320 seconds

Table 15
Comparison of different algorithms on random networks of |V | = 100 with different group sizes.Computational time for each run is 320 seconds Algorithms∈ Y PF Y PF | Ȳ | % Y PF Coverage Y MOSAGLS Y MOEA1 Y MOEA2