**A Comparative Study on Meta Heuristic Algorithms for Solving Multilevel Lot-Sizing Problems**

Ikou Kaku, Yiyong Xiao and Yi Han

Additional information is available at the end of the chapter

http://dx.doi.org/10.5772/55279

## **1. Introduction**

Material requirements planning (MRP) is an old field but still now plays an important role in coordinating replenishment decision for materials/components of complex product structure. As the key part of MRP system, the multilevel lot-sizing (MLLS) problem concerns how to determine the lot sizes forproducing/procuringmultiple items atdifferentlevels withquantita‐ tive interdependences, so as to minimize the total costs of production/procurement setup and inventory holding in the planning horizon. The problem is of great practical importance for the efficient operation of modern manufacturing and assembly processes and has been widely studied both in practice and in academic research over past half century. Optimal solution algorithms exist for the problem; however, only small instances can be solved in reasonable computationaltimebecausetheproblemisNP-hard(SteinbergandNapier,1980).Earlydynamic programmingformulationsusedanetworkrepresentationoftheproblemwithaseriesstructure (Zhangwill, 1968,1969) or an assembly structure (Crowston and Wagner,1973). Other optimal approaches involve the branch and bound algorithms (Afentakis et al., 1984, Afentakis and Gavish, 1986)thatuseda converting approach to change the classicalformulation ofthe general structure into a simple but expanded assembly structure. As the MLLS problem is so common in practice and plays a fundamental role in MRP system, many heuristic approaches have also beendeveloped,consistingfirstofthesequentialapplicationofthesingle-levellot-sizingmodels to each component of the product structure (Yelle,1979, Veral and LaForge,1985), and later, of the application of the multilevel lot-sizing models. The multilevel models quantify item interdependencies and thus perform better than the single-level based models (Blackburn and Millen, 1985, Coleman and McKnew, 1991).

Recently, meta-heuristic algorithms have been proposed to solve the MLLS problem with a low computational load. Examples of hybrid genetic algorithms (Dellaert and Jeunet, 2000, Dellaert et al., 2000), simulated annealing (Tang, 2004, Raza and Akgunduz, 2008), particle

© 2013 Kaku et al.; licensee InTech. This is an open access article distributed under the terms of the Creative Commons Attribution License (http://creativecommons.org/licenses/by/3.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited. © 2013 Kaku et al.; licensee InTech. This is a paper distributed under the terms of the Creative Commons Attribution License (http://creativecommons.org/licenses/by/3.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

swarm optimization (Han et al, 2009, 2012a, 2012b), and soft optimization approach based on segmentation (Kaku and Xu, 2006, Kaku et al, 2010), ant colony optimization system (Pitakaso et al., 2007, Almeda, 2010), variable neighborhood search based approaches (Xiao et al., 2011a, 2011b, 2012), have been developed to solve the MLLS problem of large-scale. Those meta-heuristic algorithms outperform relative simplicity in solving the MLLS problems, together with their cost efficiency, make them appealing tool to industrials, however they are unable to guarantee an optimal solution. Hence those meta-heuristic algorithms that offer a reasonable trade-off between optimality and computational feasibility are highly advisable. It is very reasonable to consider the appropriateness of the algorithm, especially is which algorithm most appropriate for solving the MLLS problems?

finished item is known up to the planning horizon, backlog is not allowed for any items, and the lead time for all production items is zero. Suppose that there are *M* items and the planning horizon is divided into *N* periods. Our purpose is to find the lot sizes of all items so as to minimize the sum of setup and inventory-holding costs, while ensuring that external demands

A Comparative Study on Meta Heuristic Algorithms for Solving Multilevel Lot-Sizing Problems

http://dx.doi.org/10.5772/55279

79

To formulate this problem as an integer optimization problem, we use the same notation of

The objective function is the sum of setup and inventory-holding costs for all items over the

, ,

+

= × +× åå (1)

, ,1 , , , *it it it it I I PD* - = +- (2)

= × " G¹ å (3)

, , 0, *it it P Mx* -× £ (4)

j

( ).

*i it i it*

for the end item are met over the *N*-period planning horizon.

Dellaert and Jeunet (2000a), as follows:

: Unit inventory-holding cost for item *i*

*Di*,*t*: Requirements for item *i* in period *t*

*Pi*,*t*: Production quantity for item *i* in period *t*

*Γ*(*i*): set of immediate successors of items i

*Ii*,*t*: Inventory level of item *i* at the end of period *t*

*xi*,*t*: Binary decision index addressed to capture the setup cost for item *i*

1 1

*TC H I S x*

, ,, | , *j*

*it ij jt l i*

*D CP i*

*i t*

The MLLS problem is to minimize TC under the following constraints:

*i*

*j*

ÎG

= =

*M N*

: Quantity of item *i* required to produce one unit of item *j*

entire planning horizon, denoted by TC (total cost). Then

*i* : Index of items, *i* = 1, 2, …, *M*

: Setup cost for item *i*

*Hi*

*Si*

*Ci*, *<sup>j</sup>*

*M*: A large number

*t* : Index of periods, *t* = 1, 2, …, *N*

In this chapter, We first review the meta-heuristic algorithms for solving the MLLS problems, especially focus on the implemental techniques and their effectives in those meta-heuristic algorithms. For simplicity the MLLS problems are limited with time-varying cost structures and no restrictive assumption on the product structure. Even so the solutions of the MLLS problems are not simply convex but becoming very complex with multi minimums when the cost struc‐ ture is time-varying and the product structure is becoming general. Comparing those imple‐ ment methods used in different algorithms we can find some essential properties of searching bettersolutionoftheMLLSproblems.Usingthepropertiestherefore,wecanspecifythecharacter‐ isticsofthealgorithmsandindicateadirectiononwhichmoreefficientalgorithmwillbedeveloped.

Second, by using these properties as an example, we present a succinct approach—iterated variable neighborhood descent (IVND), a variant of variable neighborhood search (VNS), to efficientlysolve theMLLSproblem.Toexamine theperformanceofthenewalgorithm,different kinds of product structures were considered including the component commonality and multipleend-products,and176benchmarkproblemsunderdifferent scales(small,mediumand large)wereusedtotestagainstinour computational experiments.Theperformanceofthe IVND algorithmwere comparedwiththose ofthreewell-knownalgorithms inliteratures—thehybrid genetic algorithm by Dellaert and Jeunet (2000a), the MMAS algorithm by Pitakaso et al. (2007),andtheparallelgeneticalgorithmbybyHomberger(2008), sincetheyalltackledthesame benchmarkproblems.The results show thatthe IVNDalgorithm is very competitive since it can on average find better solutions in less computing time than other three.

The rest of this chapter is organized as follows. Section 2 describes the MLLS problem. Section 3 gives an overview of related meta-heuristic algorithms firstly, and several implemental techniques used in the algorithms are discussed. Then section 4 explains the initial method and six implemental techniques used in IVND algorithm, and the scheme of the proposed IVND algorithm. In section 5, computational experiments are carried out on three 176 bench‐ mark problems to test the new algorithm of efficiency and effectiveness and compared with existing algorithms. Finally, in section 7, we symmary the chapter.

## **2. The MLLS problems**

The MLLS problem is considered to be a discrete-time, multilevel production/inventory system with an assembly structure and one finished item. We assume that external demand for the finished item is known up to the planning horizon, backlog is not allowed for any items, and the lead time for all production items is zero. Suppose that there are *M* items and the planning horizon is divided into *N* periods. Our purpose is to find the lot sizes of all items so as to minimize the sum of setup and inventory-holding costs, while ensuring that external demands for the end item are met over the *N*-period planning horizon.

To formulate this problem as an integer optimization problem, we use the same notation of Dellaert and Jeunet (2000a), as follows:

*i* : Index of items, *i* = 1, 2, …, *M*

*t* : Index of periods, *t* = 1, 2, …, *N*

*Hi* : Unit inventory-holding cost for item *i*

*Si* : Setup cost for item *i*

swarm optimization (Han et al, 2009, 2012a, 2012b), and soft optimization approach based on segmentation (Kaku and Xu, 2006, Kaku et al, 2010), ant colony optimization system (Pitakaso et al., 2007, Almeda, 2010), variable neighborhood search based approaches (Xiao et al., 2011a, 2011b, 2012), have been developed to solve the MLLS problem of large-scale. Those meta-heuristic algorithms outperform relative simplicity in solving the MLLS problems, together with their cost efficiency, make them appealing tool to industrials, however they are unable to guarantee an optimal solution. Hence those meta-heuristic algorithms that offer a reasonable trade-off between optimality and computational feasibility are highly advisable. It is very reasonable to consider the appropriateness of the algorithm, especially is which

In this chapter, We first review the meta-heuristic algorithms for solving the MLLS problems, especially focus on the implemental techniques and their effectives in those meta-heuristic algorithms. For simplicity the MLLS problems are limited with time-varying cost structures and no restrictive assumption on the product structure. Even so the solutions of the MLLS problems are not simply convex but becoming very complex with multi minimums when the cost struc‐ ture is time-varying and the product structure is becoming general. Comparing those imple‐ ment methods used in different algorithms we can find some essential properties of searching bettersolutionoftheMLLSproblems.Usingthepropertiestherefore,wecanspecifythecharacter‐ isticsofthealgorithmsandindicateadirectiononwhichmoreefficientalgorithmwillbedeveloped. Second, by using these properties as an example, we present a succinct approach—iterated variable neighborhood descent (IVND), a variant of variable neighborhood search (VNS), to efficientlysolve theMLLSproblem.Toexamine theperformanceofthenewalgorithm,different kinds of product structures were considered including the component commonality and multipleend-products,and176benchmarkproblemsunderdifferent scales(small,mediumand large)wereusedtotestagainstinour computational experiments.Theperformanceofthe IVND algorithmwere comparedwiththose ofthreewell-knownalgorithms inliteratures—thehybrid genetic algorithm by Dellaert and Jeunet (2000a), the MMAS algorithm by Pitakaso et al. (2007),andtheparallelgeneticalgorithmbybyHomberger(2008), sincetheyalltackledthesame benchmarkproblems.The results show thatthe IVNDalgorithm is very competitive since it can

algorithm most appropriate for solving the MLLS problems?

78 Recent Advances on Meta-Heuristics and Their Application to Real Scenarios

on average find better solutions in less computing time than other three.

existing algorithms. Finally, in section 7, we symmary the chapter.

**2. The MLLS problems**

The rest of this chapter is organized as follows. Section 2 describes the MLLS problem. Section 3 gives an overview of related meta-heuristic algorithms firstly, and several implemental techniques used in the algorithms are discussed. Then section 4 explains the initial method and six implemental techniques used in IVND algorithm, and the scheme of the proposed IVND algorithm. In section 5, computational experiments are carried out on three 176 bench‐ mark problems to test the new algorithm of efficiency and effectiveness and compared with

The MLLS problem is considered to be a discrete-time, multilevel production/inventory system with an assembly structure and one finished item. We assume that external demand for the *Ii*,*t*: Inventory level of item *i* at the end of period *t*

*xi*,*t*: Binary decision index addressed to capture the setup cost for item *i*

*Di*,*t*: Requirements for item *i* in period *t*

*Pi*,*t*: Production quantity for item *i* in period *t*

*Ci*, *<sup>j</sup>* : Quantity of item *i* required to produce one unit of item *j*

*Γ*(*i*): set of immediate successors of items i

*M*: A large number

The objective function is the sum of setup and inventory-holding costs for all items over the entire planning horizon, denoted by TC (total cost). Then

$$\text{TCC} = \sum\_{i=1}^{M} \sum\_{t=1}^{N} (H\_i \cdot I\_{i,t} + \mathcal{S}\_i \cdot \boldsymbol{\chi}\_{i,t}). \tag{1}$$

The MLLS problem is to minimize TC under the following constraints:

$$I\_{i,t} = I\_{i,t-1} + P\_{i,t} - D\_{i,t'} \tag{2}$$

$$D\_{i,t} = \sum\_{j \neq \Gamma\_i} \mathbb{C}\_{i,j} \cdot \mathbb{P}\_{j, t+l\_j} \qquad \forall i \mid \Gamma\_i \neq \emptyset \,\prime\tag{3}$$

$$P\_{i,t} - M \cdot \boldsymbol{\pi}\_{i,t} \le 0,\tag{4}$$

$$I\_{i,t} \ge 0, \quad P\_{i,t} \ge 0, \quad \propto\_{i,t} \in \{0, 1\}, \quad \forall i, t. \tag{5}$$

**3.1. Soft optimization approach**

samples or by narrowing the range of segment.

It is the point that we try to define and solve in this chapter.

**3.2. Genetic algorithm**

**3.3. Simulated annealing**

Soft optimization approach (SOA) for solving the MLLS problems is based on a general sampling approach (Kaku and Xu, 2006; Kaku et al, 2010). The main merit of soft optimization approach is that it does not require any structure information about the objective function, so it can be used to treat optimization problems with complicated structures. However, it was shown that random sampling (for example simple uniform sampling) method cannot produce a good solution. Several experiments had been derived to find the characteristics of an optimal solution, and as a result applying the solution structure information of the MLLS problem to the sampling method may produce a better result than that arising from the simple uniform sampling method. A heuristic algorithm to segment the solution space with percentage of number of 1s has been developed and the performance improvement of solving MLLS problem was confirmed. It should be pointed that the SOA based on segmentation still remains the essential property of random sampling but limited with the searching ranges, however the adopted new solution(s) does not remain any information of the old solution(s). Therefore the improvement of solution performance can only be achieved by increasing the numbers of

A Comparative Study on Meta Heuristic Algorithms for Solving Multilevel Lot-Sizing Problems

http://dx.doi.org/10.5772/55279

81

Genetic algorithm (GA) has been developed firstly for solving the MLLS problems in (Dellaert and Jeunet, 2000a, Dellaert et al. 2000b, Han et al. 2012a, 2012b). In fact, it firstly created the way that solving MLLS problems by using meta-heuristic algorithms. Several important contributions were achieved. Firstly a very general GA approach was developed and improved by using several specific genetic operators and a roulette rule to gather those operators had been implemented to treat the two dimensional chromosomes. Secondly, comparison studies had been provided to show that better solution could be obtained than those existing heuristic algorithms, based on several benchmarks data collected from literature. Later such bench‐ marks provide a platform of developing meta-heuristic algorithms and evaluating their performance. Because the complexity of MLLS problem and the flexibility and implement ability of GA are matching each other so that GA seems powerful and effective for solving MLLS problem. However, even several operators as single bit mutation; cumulative mutation; inversion; period crossover and product crossover were combined in the GA algorithm but what operators were effective in better solution searching process was not presented clearly.

Simulated annealing (SA) has been also developed to solve the MLLS problem (Tang,2004; Raza and Akgunduz,2008). It starts from a random initial solution and changing neighbouring states of the incumbent solution by a cooling process, in which the new solution is accepted or rejected according to a possibility specified by the Metropolis algorithm. Also parameters used in the algorithmhadbeeninvestigatedbyusingtheanalysisofvarianceapproach.Ithadbeenreported thatthe SAis appropriate for solvingtheMLLSproblemhoweververifiedonlyinverysmalltest problems. Based on our understanding for SA, different from other meta-heuristic algorithms

where Equation 2 expresses the flow conservation constraint for item *i*. Note that, if item *i* is an end item (characterized by *Γ*(*i*)=*φ*), its demand is exogenously given, whereas if it is a component (such that *Γ*(*i*)≠*φ*), its demand is defined by the production of its successors (items belonging to *Γ*(*i*) as stated by Equation 3). Equation 3 guarantees that the amount *P <sup>j</sup>*,*t* of item *j* available in period *t* results from the exact combination of its predecessors (items belonging to *Γ<sup>i</sup>* in period *t*). Equation 4 guarantees that a setup cost is incurred when a production is produced. Equation 5 states that backlog is not allowed, production is either positive or zero, and that decision variables are 0, 1 variables.

Because *xi*,*<sup>t</sup>* ∈{0, 1} is a binary decision variable for item *i* in period *t*, *X* ={*xi*,*t*}*<sup>M</sup>* <sup>×</sup>*N* represents the solution space of the MLLS problem. Searching for an optimal solution of the MLLS problem is equivalent to finding a binary matrix that produces a minimum sum of the setup and inventory-holding costs. Basically, there exists an optimal solution if

$$\mathbf{x}\_{i,t} \cdot I\_{i,t-1} = \mathbf{0} \tag{6}$$

Equation 6 indicates that any optimal lot must cover an integer number of periods of future demand. We set the first column of X to be 1 to ensure that the initial production is feasible because backlog is not allowed for any item and the lead times are zero. Since there is an inner corner property for assembly structure (see Tang (2004)), we need to have *xi*,*<sup>t</sup>* ≥ *xk* ,*t* if item creates internal demand for item k. Thus we need a constraint in order to guarantee that the matrix is feasible.

#### **3. Meta-heuristic algorithms used to solve MLLS problems**

The meta-heuristic algorithms are widely used to refer to simple, hybrid and population-based stochastic local searching (Hoos and Thomas 2005). They transfer the principle of evolution through mutation, recombination and selection of the fittest, which leads to the development of solutions that are better adapted for survival in a given environment, to solving computa‐ tionally hard problems. However, those algorithms often seem to lack the capability of sufficient search intensification, that is, the ability to reach high-quality candidate solutions efficiently when a good starting position is given. Hence, in many cases, the performance of the algorithms for combinatorial problems can be significantly improved by adding some implemental techniques that are used to guide an underlying problem-specific heuristic. In this chapter, our interesting is on the mechanisms of those implemental techniques used to solve a special combinatorial optimization problem, i.e. MLLS problem. Hence we first review several existing meta-heuristic algorithms for solving the MLLS problems.

## **3.1. Soft optimization approach**

,, , 0, 0, {0,1}, , . *it it i t IP x* ³³ Î "*i t* (5)

, ,1 0 *it it x I* - × = (6)

where Equation 2 expresses the flow conservation constraint for item *i*. Note that, if item *i* is an end item (characterized by *Γ*(*i*)=*φ*), its demand is exogenously given, whereas if it is a component (such that *Γ*(*i*)≠*φ*), its demand is defined by the production of its successors (items belonging to *Γ*(*i*) as stated by Equation 3). Equation 3 guarantees that the amount *P <sup>j</sup>*,*t* of item *j* available in period *t* results from the exact combination of its predecessors (items belonging

 in period *t*). Equation 4 guarantees that a setup cost is incurred when a production is produced. Equation 5 states that backlog is not allowed, production is either positive or zero,

Because *xi*,*<sup>t</sup>* ∈{0, 1} is a binary decision variable for item *i* in period *t*, *X* ={*xi*,*t*}*<sup>M</sup>* <sup>×</sup>*N* represents the solution space of the MLLS problem. Searching for an optimal solution of the MLLS problem is equivalent to finding a binary matrix that produces a minimum sum of the setup

Equation 6 indicates that any optimal lot must cover an integer number of periods of future demand. We set the first column of X to be 1 to ensure that the initial production is feasible because backlog is not allowed for any item and the lead times are zero. Since there is an inner corner property for assembly structure (see Tang (2004)), we need to have *xi*,*<sup>t</sup>* ≥ *xk* ,*t* if item creates internal demand for item k. Thus we need a constraint in order to guarantee that the

The meta-heuristic algorithms are widely used to refer to simple, hybrid and population-based stochastic local searching (Hoos and Thomas 2005). They transfer the principle of evolution through mutation, recombination and selection of the fittest, which leads to the development of solutions that are better adapted for survival in a given environment, to solving computa‐ tionally hard problems. However, those algorithms often seem to lack the capability of sufficient search intensification, that is, the ability to reach high-quality candidate solutions efficiently when a good starting position is given. Hence, in many cases, the performance of the algorithms for combinatorial problems can be significantly improved by adding some implemental techniques that are used to guide an underlying problem-specific heuristic. In this chapter, our interesting is on the mechanisms of those implemental techniques used to solve a special combinatorial optimization problem, i.e. MLLS problem. Hence we first review

and inventory-holding costs. Basically, there exists an optimal solution if

**3. Meta-heuristic algorithms used to solve MLLS problems**

several existing meta-heuristic algorithms for solving the MLLS problems.

to *Γ<sup>i</sup>*

matrix is feasible.

and that decision variables are 0, 1 variables.

80 Recent Advances on Meta-Heuristics and Their Application to Real Scenarios

Soft optimization approach (SOA) for solving the MLLS problems is based on a general sampling approach (Kaku and Xu, 2006; Kaku et al, 2010). The main merit of soft optimization approach is that it does not require any structure information about the objective function, so it can be used to treat optimization problems with complicated structures. However, it was shown that random sampling (for example simple uniform sampling) method cannot produce a good solution. Several experiments had been derived to find the characteristics of an optimal solution, and as a result applying the solution structure information of the MLLS problem to the sampling method may produce a better result than that arising from the simple uniform sampling method. A heuristic algorithm to segment the solution space with percentage of number of 1s has been developed and the performance improvement of solving MLLS problem was confirmed. It should be pointed that the SOA based on segmentation still remains the essential property of random sampling but limited with the searching ranges, however the adopted new solution(s) does not remain any information of the old solution(s). Therefore the improvement of solution performance can only be achieved by increasing the numbers of samples or by narrowing the range of segment.

## **3.2. Genetic algorithm**

Genetic algorithm (GA) has been developed firstly for solving the MLLS problems in (Dellaert and Jeunet, 2000a, Dellaert et al. 2000b, Han et al. 2012a, 2012b). In fact, it firstly created the way that solving MLLS problems by using meta-heuristic algorithms. Several important contributions were achieved. Firstly a very general GA approach was developed and improved by using several specific genetic operators and a roulette rule to gather those operators had been implemented to treat the two dimensional chromosomes. Secondly, comparison studies had been provided to show that better solution could be obtained than those existing heuristic algorithms, based on several benchmarks data collected from literature. Later such bench‐ marks provide a platform of developing meta-heuristic algorithms and evaluating their performance. Because the complexity of MLLS problem and the flexibility and implement ability of GA are matching each other so that GA seems powerful and effective for solving MLLS problem. However, even several operators as single bit mutation; cumulative mutation; inversion; period crossover and product crossover were combined in the GA algorithm but what operators were effective in better solution searching process was not presented clearly. It is the point that we try to define and solve in this chapter.

#### **3.3. Simulated annealing**

Simulated annealing (SA) has been also developed to solve the MLLS problem (Tang,2004; Raza and Akgunduz,2008). It starts from a random initial solution and changing neighbouring states of the incumbent solution by a cooling process, in which the new solution is accepted or rejected according to a possibility specified by the Metropolis algorithm. Also parameters used in the algorithmhadbeeninvestigatedbyusingtheanalysisofvarianceapproach.Ithadbeenreported thatthe SAis appropriate for solvingtheMLLSproblemhoweververifiedonlyinverysmalltest problems. Based on our understanding for SA, different from other meta-heuristic algorithms like GA, SA is rather like a local successive search approach from an initial solution. Then almost information of the old solution can be remained which may lead a long time to search better solution if it is far from the original solution. Also several implement methods con be used to improve the effective of SA (see Hoos and Thomas 2005). It is a general point to improve the effective of SA through shortening the cooling time with some other local searching methods.

**3.6. Variable neighbourhood search**

Repeat the following step (1), (2) and (3):

used to terminate the program.

**1.** Single- bit mutation

**2.** Cumulative mutation

of its predecessors.

(2)Repeat (a) and (b) steps:

(1)Find an initial solution(s) by using the *init* function.

(b)Decide whether the new solution should be accepted.

**Figure 1.** General algorithm for solving the MLLS problem

Variable neighborhood search (VNS) is also used to solve the MLLS problem (Xiao et al., 2011a, 2011b, 2012). The main reasoning of this searching strategy, in comparison to most local search heuristics of past where only one neighborhood is used, is based on the idea of a systematic change of predefined and ordered multi neighborhoods within a local search. By introducing a set of distance-leveled neighborhoods and correlated exploration order, the variable neighborhood search algorithm can perform a high efficient searching in nearby

A Comparative Study on Meta Heuristic Algorithms for Solving Multilevel Lot-Sizing Problems

http://dx.doi.org/10.5772/55279

83

There may other different meta-heuristic algorithms have been proposed for solving the MLLS problems, but it can be considered that the algorithms updated above can cover almost fields of the meta-heuristic so that the knowledge obtained in this chapter may has high applicable values. All of the meta-heuristic algorithms used to solve the MLLS problems are generally

neighborhoods where better solutions are more likely to be found.

constructed by the algorithm describe in Fig.1 as follows.

(a)Find a new solution randomly in solution space by using the *step* function;

(3)Terminate the search process by using the *terminate* function and output the best solution found.

In Fig. 1, at the beginning of the search process, an initial solution(s) is generated by using the function *init.* A simple *init* method used to generate the initial solution may be the random sampling method, and often several heuristic concepts of MLLS problem are employed to initialize the solution because they can help obtaining better performance. Moreover, the *init* methods are usually classified into two categories in terms of single solution or multi solutions (usually called population). Function *step* shows the own originality of the algorithm by using different implement techniques. By comparing the total costs of old and new solutions, it accepts usually the solution with smaller one as next incumbent solution. Finally, function *terminate*, a user-defined function such as numbers of calculation, coverage rate and so on, is

All of the implement techniques used in *step* function, which are represented in all of the meta-

Just one position (*i,t*) has randomly selected to change its value (from 0 to 1 or the reverse)

When a single-bit mutation is performed on a given item, it may trigger mutating the value(s)

heuristic algorithms for solving the MLLS problems, can be summarized as below.

#### **3.4. Particle swarm optimization**

Particle swarm optimization (PSO) is also a meta-heuristic algorithm formally introduced (Han et al,2009, 2011). It is a suitable and effective tool for continuous optimization problems. Recently the standard particle swarm optimization algorithm is also converted into a discrete version through redefining all the mathematical operators to solve the MLLS problems (Han et al, 2009). It starts its search procedure with a particle swarm. Each particle in the swarm keeps track of its coordinates in the problem space, which are associated with the best solution (fitness) it has achieved so far. This value is called pBest. Another "best" value tracked by the global version of the particle swarm optimization is the overall best value, and its location, obtained so far by any particle in the population, which is called gBest. Gather those so-called optimal factors into current solutions then they will converge to a better solution. It has been reported that comparing experiments with GA proposed in (Dellaert and Jeunet, 2000a, Dellaert et al. 2000b) had been executed by using the same benchmarks and better performance were obtained. Consider the essential mechanism of PSO, it is clear that those optimal factors (pBest and gBest) follow the information of the particle passed and direct where the better solution being. However, it has not been explained clearly that whether those optimal factors remained the required information when the PSO is converted into a discrete form.

## **3.5. Ant colony optimization**

A special ant colony optimization (ACO) combined with linear program has been developed recently for solving the MLLS problem (Pitakaso et al. 2007, Almeda 2010). The basic idea of ant colony optimization is that a population of agents (artificial ants) repeatedly constructs solutions to a given instance of a combinatorial optimization problem. Ant colony optimization had been used to select the principle production decisions, i.e. for which period production for an item should be schedules in the MLLS problems. It starts from the top items down to the raw materials according to the ordering given by the bill of materials. The ant's decision for production in a certain period is based on the pheromone information as well as on the heuristic information if there is an external (primary) or internal (secondary) demand. The pheromone information represents the impact of a certain production decision on the objective values of previously generated solutions, i.e. the pheromone value is high if a certain produc‐ tion decision has led to good solution in previous iterations. After the selection of the produc‐ tion decisions, a standard LP solver has been used to solve the remaining linear problem. After all ants of an iteration have constructed a solution, the pheromone information is updated by the iteration best as well as the global best solutions. This approach has been reported works well for small and medium-size MLLS problems. However for large instances the solution method leads to high-quality results, but cannot beat highly specialized algorithms.

## **3.6. Variable neighbourhood search**

like GA, SA is rather like a local successive search approach from an initial solution. Then almost information of the old solution can be remained which may lead a long time to search better solution if it is far from the original solution. Also several implement methods con be used to improve the effective of SA (see Hoos and Thomas 2005). It is a general point to improve the effective of SA through shortening the cooling time with some other local searching methods.

Particle swarm optimization (PSO) is also a meta-heuristic algorithm formally introduced (Han et al,2009, 2011). It is a suitable and effective tool for continuous optimization problems. Recently the standard particle swarm optimization algorithm is also converted into a discrete version through redefining all the mathematical operators to solve the MLLS problems (Han et al, 2009). It starts its search procedure with a particle swarm. Each particle in the swarm keeps track of its coordinates in the problem space, which are associated with the best solution (fitness) it has achieved so far. This value is called pBest. Another "best" value tracked by the global version of the particle swarm optimization is the overall best value, and its location, obtained so far by any particle in the population, which is called gBest. Gather those so-called optimal factors into current solutions then they will converge to a better solution. It has been reported that comparing experiments with GA proposed in (Dellaert and Jeunet, 2000a, Dellaert et al. 2000b) had been executed by using the same benchmarks and better performance were obtained. Consider the essential mechanism of PSO, it is clear that those optimal factors (pBest and gBest) follow the information of the particle passed and direct where the better solution being. However, it has not been explained clearly that whether those optimal factors

remained the required information when the PSO is converted into a discrete form.

method leads to high-quality results, but cannot beat highly specialized algorithms.

A special ant colony optimization (ACO) combined with linear program has been developed recently for solving the MLLS problem (Pitakaso et al. 2007, Almeda 2010). The basic idea of ant colony optimization is that a population of agents (artificial ants) repeatedly constructs solutions to a given instance of a combinatorial optimization problem. Ant colony optimization had been used to select the principle production decisions, i.e. for which period production for an item should be schedules in the MLLS problems. It starts from the top items down to the raw materials according to the ordering given by the bill of materials. The ant's decision for production in a certain period is based on the pheromone information as well as on the heuristic information if there is an external (primary) or internal (secondary) demand. The pheromone information represents the impact of a certain production decision on the objective values of previously generated solutions, i.e. the pheromone value is high if a certain produc‐ tion decision has led to good solution in previous iterations. After the selection of the produc‐ tion decisions, a standard LP solver has been used to solve the remaining linear problem. After all ants of an iteration have constructed a solution, the pheromone information is updated by the iteration best as well as the global best solutions. This approach has been reported works well for small and medium-size MLLS problems. However for large instances the solution

**3.4. Particle swarm optimization**

82 Recent Advances on Meta-Heuristics and Their Application to Real Scenarios

**3.5. Ant colony optimization**

Variable neighborhood search (VNS) is also used to solve the MLLS problem (Xiao et al., 2011a, 2011b, 2012). The main reasoning of this searching strategy, in comparison to most local search heuristics of past where only one neighborhood is used, is based on the idea of a systematic change of predefined and ordered multi neighborhoods within a local search. By introducing a set of distance-leveled neighborhoods and correlated exploration order, the variable neighborhood search algorithm can perform a high efficient searching in nearby neighborhoods where better solutions are more likely to be found.

There may other different meta-heuristic algorithms have been proposed for solving the MLLS problems, but it can be considered that the algorithms updated above can cover almost fields of the meta-heuristic so that the knowledge obtained in this chapter may has high applicable values. All of the meta-heuristic algorithms used to solve the MLLS problems are generally constructed by the algorithm describe in Fig.1 as follows.

Repeat the following step (1), (2) and (3): (1)Find an initial solution(s) by using the *init* function. (2)Repeat (a) and (b) steps: (a)Find a new solution randomly in solution space by using the *step* function; (b)Decide whether the new solution should be accepted. (3)Terminate the search process by using the *terminate* function and output the best solution found.

**Figure 1.** General algorithm for solving the MLLS problem

In Fig. 1, at the beginning of the search process, an initial solution(s) is generated by using the function *init.* A simple *init* method used to generate the initial solution may be the random sampling method, and often several heuristic concepts of MLLS problem are employed to initialize the solution because they can help obtaining better performance. Moreover, the *init* methods are usually classified into two categories in terms of single solution or multi solutions (usually called population). Function *step* shows the own originality of the algorithm by using different implement techniques. By comparing the total costs of old and new solutions, it accepts usually the solution with smaller one as next incumbent solution. Finally, function *terminate*, a user-defined function such as numbers of calculation, coverage rate and so on, is used to terminate the program.

All of the implement techniques used in *step* function, which are represented in all of the metaheuristic algorithms for solving the MLLS problems, can be summarized as below.

**1.** Single- bit mutation

Just one position (*i,t*) has randomly selected to change its value (from 0 to 1 or the reverse)

**2.** Cumulative mutation

When a single-bit mutation is performed on a given item, it may trigger mutating the value(s) of its predecessors.

## **3.** Inversion

One position (*i,t*) has been randomly selected, then compare its value to that of its neighbor in position (*i,t+1*). If the two values are different from each other(1 and 0, or 1 and 0), then exchanges the two values. Note the last period *t*=*N* should not be selected for inversion operation.

Note the implement techniques can change the states of the incumbent solution(s) to improve the performance with respect to total cost function, so that which implement technique could do how many contributions to the solution improvement is very important for the computing efficiency of MLLS problems. Here our interesting is on those implement techniques used in the above algorithms, which had been reported appropriate in their calculation situations. We define several criteria for evaluating the performance and effective of the implemental techniques. We firstly define distance-based neighborhood structures of the incumbent solution. The neighborhoods of incumbent solution are sets of feasible solutions associated with different distance measurements from the incumbent solutions. The distance means the

A Comparative Study on Meta Heuristic Algorithms for Solving Multilevel Lot-Sizing Problems

http://dx.doi.org/10.5772/55279

85

**Definition 1**. *Distance from incumbent solution*: For a set of feasible solutions of a MLLS problem, i.e., *X* ={*xi*,*<sup>t</sup>*}, a solution *X* ' belongs to the *k*th-distance of incumbent solution *x*, i.e. *Nk*(*x*), if and only if it satisfies, *x* '∈ *Nk* (*x*)⇔*ρ*(*x*, *x* ')=*k*, where *k* is a positive integer. Distance between any

Where | • \ • | denotes the number of different points between two solutions, i.e.,

For example of a MLLS problem with 3 items and 3 periods, and three feasible solutions: *x*, *x'*


According to **Definition 1** above, we can get their distances such that:

Therefore, creating a solution with *k*th distance, i.e. *Nk*(*x*), to the incumbent solution can be realized by changing the values of *k* different points of incumbent solution. It can be considered that less changes, e.g. *N*1(*x*), can avoid too much damage to the maybe good properties of the incumbent solution, in which some points may be already in its optimal positions. While multiple changes, e.g. *Nk*(*x*) (k>1), lead a new solution to be very different from the incumbent solution, so it may also destroy the goodness of the original solution. However on the other hand local optimization may be overcame by multiple changes. Hence, following hypothesis may be considered to evaluate the performance and effective of an implement technique. Secondly we define range measurements that means the changing points are how far from the incumbent solutions. Eventually, we can evaluate those implemental techniques used in the algorithms to solve the MLLS problems by measuring the distance and range when the solution

1 0 1 1 0 1


( , ') \ ' '\ , ' *xx x x x x xx X* = = "Î (7)

exchanged number of different points of two incumbent solutions.

two elements in in *X*, e.g., *x* and *x'*, is measured by

r



*ρ*(*x*, *x* ')=1, *ρ*(*x* ', *x* ' ')=1, *ρ*(*x*, *x* ' ')=2.

1 0 1 1 1 1


*<sup>x</sup>* <sup>=</sup> <sup>|</sup><sup>1</sup> <sup>0</sup> <sup>0</sup> 1 0 0 1 1 1

has been changed.

*i*=1 *M* ∑ *t*=1 *N*

and *x''*, which are as following,

## **4.** Period crossover

A point *t* in period (1, N) is randomly selected for two parents to produce off-springs by horizontally exchanging half part of their solutions behind period *t*.

## **5.** Product crossover

An item *i* is randomly selected for two parents to produce off-springs by vertically exchanging half part of their solutions below item *i*.

Table 1 summaries what implement techniques and *init* function have been used in the existing meta-heuristic algorithms. From Table 1 it can be observed that all of the algorithms use the singlebit mutation method in their *step* function, but their ways of making the mutations are different. Clearly, SOA creates new solution(s) without using any previous information so that it may is recognized a non-implemental approach. Reversely SA uses single-bit and cumulative muta‐ tion based on incumbent solution therefore it reserves almost all of the previous information. Moreover, a random uphill technique is used in SA to escape from local optima so that the global optimal solution may is obtained. However a very long computing time is needed to get it. On the other hand, PSO uses the single-bit mutation method to create the new candidate in *step* function but starting with multiple initial solutions. According to the concept of original PSO, some information of solutions obtained before may be added into the construction of new solutions in order to increase the probability of generating better solutions. However, it needs a proof of such excellent property is guaranteed in the implemental process but is not appeared in the algorithm proposed before. It may is a challenge work for developing a new PSO algorithm for solving the MLLS problem. Finally, because the techniques of inversion and crossover (which were usually used in other combinational optimization problem such as a Travelling Salesman Problem) only have been used in GA, it is not able to compare them with other algorithms.


**Table 1.** The implement techniques used in various algorithms

Note the implement techniques can change the states of the incumbent solution(s) to improve the performance with respect to total cost function, so that which implement technique could do how many contributions to the solution improvement is very important for the computing efficiency of MLLS problems. Here our interesting is on those implement techniques used in the above algorithms, which had been reported appropriate in their calculation situations. We define several criteria for evaluating the performance and effective of the implemental techniques. We firstly define distance-based neighborhood structures of the incumbent solution. The neighborhoods of incumbent solution are sets of feasible solutions associated with different distance measurements from the incumbent solutions. The distance means the exchanged number of different points of two incumbent solutions.

**Definition 1**. *Distance from incumbent solution*: For a set of feasible solutions of a MLLS problem, i.e., *X* ={*xi*,*<sup>t</sup>*}, a solution *X* ' belongs to the *k*th-distance of incumbent solution *x*, i.e. *Nk*(*x*), if and only if it satisfies, *x* '∈ *Nk* (*x*)⇔*ρ*(*x*, *x* ')=*k*, where *k* is a positive integer. Distance between any two elements in in *X*, e.g., *x* and *x'*, is measured by

$$\varphi(\mathbf{x}, \mathbf{x'}) = \begin{vmatrix} \mathbf{x} \ \mathbf{x'} \end{vmatrix} = \begin{vmatrix} \mathbf{x'} \ \mathbf{x} \end{vmatrix} \qquad \forall \mathbf{x}, \mathbf{x'} \in X \tag{7}$$

Where | • \ • | denotes the number of different points between two solutions, i.e., | • \ • | =∑ *i*=1 *M* ∑ *t*=1 *N* | *xi*,*<sup>t</sup>* − *xi*,*<sup>t</sup>* '|

For example of a MLLS problem with 3 items and 3 periods, and three feasible solutions: *x*, *x'* and *x''*, which are as following,


**3.** Inversion

operation.

**4.** Period crossover

**5.** Product crossover

half part of their solutions below item *i*.

**Single- bit mutation**

One position (*i,t*) has been randomly selected, then compare its value to that of its neighbor in position (*i,t+1*). If the two values are different from each other(1 and 0, or 1 and 0), then exchanges the two values. Note the last period *t*=*N* should not be selected for inversion

A point *t* in period (1, N) is randomly selected for two parents to produce off-springs by

An item *i* is randomly selected for two parents to produce off-springs by vertically exchanging

Table 1 summaries what implement techniques and *init* function have been used in the existing meta-heuristic algorithms. From Table 1 it can be observed that all of the algorithms use the singlebit mutation method in their *step* function, but their ways of making the mutations are different. Clearly, SOA creates new solution(s) without using any previous information so that it may is recognized a non-implemental approach. Reversely SA uses single-bit and cumulative muta‐ tion based on incumbent solution therefore it reserves almost all of the previous information. Moreover, a random uphill technique is used in SA to escape from local optima so that the global optimal solution may is obtained. However a very long computing time is needed to get it. On the other hand, PSO uses the single-bit mutation method to create the new candidate in *step* function but starting with multiple initial solutions. According to the concept of original PSO, some information of solutions obtained before may be added into the construction of new solutions in order to increase the probability of generating better solutions. However, it needs a proof of such excellent property is guaranteed in the implemental process but is not appeared in the algorithm proposed before. It may is a challenge work for developing a new PSO algorithm for solving the MLLS problem. Finally, because the techniques of inversion and crossover (which were usually used in other combinational optimization problem such as a Travelling Salesman Problem) only

**Inversion Period crossover**

SOA ― ― ― ― ― many GA ○ ○ ○ ○ ○ multi SA ○ ○ ― ― ― single PSO ○ ○ ― ― ― multi ACO ○ ― ― ― ― single VNS ○ ○ ○ ― ― single

**Product**

**crossover initial solution**

horizontally exchanging half part of their solutions behind period *t*.

84 Recent Advances on Meta-Heuristics and Their Application to Real Scenarios

have been used in GA, it is not able to compare them with other algorithms.

**Cumulative mutation**

**Table 1.** The implement techniques used in various algorithms

According to **Definition 1** above, we can get their distances such that: *ρ*(*x*, *x* ')=1, *ρ*(*x* ', *x* ' ')=1, *ρ*(*x*, *x* ' ')=2.

Therefore, creating a solution with *k*th distance, i.e. *Nk*(*x*), to the incumbent solution can be realized by changing the values of *k* different points of incumbent solution. It can be considered that less changes, e.g. *N*1(*x*), can avoid too much damage to the maybe good properties of the incumbent solution, in which some points may be already in its optimal positions. While multiple changes, e.g. *Nk*(*x*) (k>1), lead a new solution to be very different from the incumbent solution, so it may also destroy the goodness of the original solution. However on the other hand local optimization may be overcame by multiple changes. Hence, following hypothesis may be considered to evaluate the performance and effective of an implement technique. Secondly we define range measurements that means the changing points are how far from the incumbent solutions. Eventually, we can evaluate those implemental techniques used in the algorithms to solve the MLLS problems by measuring the distance and range when the solution has been changed.

**Definition 2**. *Changing Range of incumbent solution*: Any elements belong to the *k*th-distance of incumbent solution *x*, i.e. *Nk*(*x*), have a range in which the changed items and periods of incumbent solution have been limited.

distance range and searching effective. Then we investigate the performance effective of various implemental techniques. For calculating the distance, set *k*←1, and repeated to find randomly a better solution from *N*k(*x*) and accept it until non better solution could be found (indicated by 100 continuous tries without improvement). And then, search better solutions from next distance (*k*←*k*+1) until *k*>*k*max. The calculation is repeated for 10 runs in different rages and the averages are used. In the *step* of GA, each of the five operators (single-bit mutation, cumulative mutation, inversion, period crossover and product crossover) are used with a fixed probability 20%, to produce new generation. Initial populations are generated and maintained with a population of 50 solutions. In each generation, new 50 solutions are randomly generated and added into the populations, then genetic operation with the five operators was performed and related total cost has been calculated. In the selection, top 50 solutions starting from lowest total cost are remained as seeds for the next generation and the rest solutions are removed from the populations. We use a stop rule of maximum 50 generations in *terminate* function. The successful results are counted in terms of distance and range (maximum deviation item/ period position among changed points). If there is only one point changed, then the range is zero. Finaly, simulation is excuted by using the same banch marks of Dellaert and Jeunet

A Comparative Study on Meta Heuristic Algorithms for Solving Multilevel Lot-Sizing Problems

http://dx.doi.org/10.5772/55279

87

Table 2 and 3 show an example of better solutions with different distances in different ranges. Better solution is defined as that the new solutions are better than the incumbent, so we count the number of better solutions and calculate the ratio of different distances. It can be observed from Table 2 that different distances lead to different number of better solutions found. Most of them (94.82%) are found in the nearest neighbourhood with one point distance and the ratio decreases as the distance increasing, which indicates a less probability of finding better solution among those whose distance are long from the incumbent. However, the incumbent solution did be improved by candidate solutions from longer distance. Even so the results with same tendency can be observed from Table 3 in which a limited range has been used. However very different meanings can be observed from Table 3. Firstly, limiting the search range can lead to a more efficiency for searching better solutions, which is represented by the fact that the total number of better solutions found is about 4% more than that of the non-limited search, and it also leads to obtain a solution with a better total cost (comparing with Table 2). Secondly, even the number of better solutions are almost same in distance 1 (changing just one position of the incumbent solution), but the number of better solutions in distance 2 was about three times of that of non-limited search. That is, longer time and less effect were performed in distance 1 if the range is not limited. This observation can lead a considerable result of implementing a faster searching in distance 1 and a faster changing between distances. That is the superiors of

> **Distance Better solutions Ratio** 1 3333 94.82% 2 135 3.84%

(2000a).

limited searching range.

**Remark 1**: Distance, range and searching effective

Therefore, if the changing result of an incumbent solution by an implement technique falls into a small range, it seems be a local search and hence may give little influence on the total cost function. Otherwise a change of incumbent solution within a large range may be a global search increasing the probability of worse solutions found and resulting lower searching efficiency. Comparing the considerations of constructing mutation in the algorithms above, we can find that only short distance (for cumulative mutation, it is a very special case of long distance) has been executed in GA, SA and PSO with different implement techniques. For example, all of the algorithms use mutation to change the solution but only GA uses the probability of obtaining a good solution to increase the share rate of the operators. While SA uses a probability of accepting a worse solution to create an uphill ability of escaping from the traps of local optima, and PSO uses other two factors (*pBest* and *gBest* ) to point a direction in which the optimal solution maybe exist. SOA and ACO are using basic random sampling principle to select the decisions in all positions of item and period. The difference is that SOA does not remain any information of previous solutions so large numbers of solution should be initialed and long distance (about half of items\*periods) always existed, whereas ACO uses a level by level approach to do single-bit mutation in a solution so goodness of previous solution may has been remained with long distance. Finally, VNS directly use distances to search good solution so that its distance is controllable.

Next, regarding the other three implement techniques, i.e., the inversion, the item crossover, and the period crossover, which are not used in other algorithms except GA, we can find that the inversion is just a special mutation that changes the solution with two points distance. However crossover(s) may change the original solutions with a longer distance so only partial goodness of previous solution has been remained. In GA algorithm, the implement techniques with long and short distance are combined together, therefore it should be considered more effective for solving the MLLS problem. However, we have to say that those implement techniques used in above algorithms (includes GA) seem not effective in finding good solution. Because even a good solution exists near by the original solution, here is no way to guarantee that it will be found by above algorithms since the candidate solution is randomly generated. It can be considered that not only distance but also range should be used in the *step* function to find a better solution. However, two problems here need to answer. Can distance and range evaluate the effective of implement techniques using in the MLLS problem? How do the implement techniques should be applied?

For proofing our consideration above, simulation experiments were executed by using the general algorithm for solving the MLLS problems shown in Fig.1. Initial solution(s) is pro‐ duced with a randomized cumulative Wagner and Whitin (RCWW) algorithm. The *step* function uses a GA algorithm, because all kinds of the implemental techniques had been included in the algorithm. We can control what the implemental technique should be used then evaluate the efficiency of the techniques. However the *terminate*functions used in different problem instances are different. In the experiments we first show the relationship among distance range and searching effective. Then we investigate the performance effective of various implemental techniques. For calculating the distance, set *k*←1, and repeated to find randomly a better solution from *N*k(*x*) and accept it until non better solution could be found (indicated by 100 continuous tries without improvement). And then, search better solutions from next distance (*k*←*k*+1) until *k*>*k*max. The calculation is repeated for 10 runs in different rages and the averages are used. In the *step* of GA, each of the five operators (single-bit mutation, cumulative mutation, inversion, period crossover and product crossover) are used with a fixed probability 20%, to produce new generation. Initial populations are generated and maintained with a population of 50 solutions. In each generation, new 50 solutions are randomly generated and added into the populations, then genetic operation with the five operators was performed and related total cost has been calculated. In the selection, top 50 solutions starting from lowest total cost are remained as seeds for the next generation and the rest solutions are removed from the populations. We use a stop rule of maximum 50 generations in *terminate* function. The successful results are counted in terms of distance and range (maximum deviation item/ period position among changed points). If there is only one point changed, then the range is zero. Finaly, simulation is excuted by using the same banch marks of Dellaert and Jeunet (2000a).

#### **Remark 1**: Distance, range and searching effective

**Definition 2**. *Changing Range of incumbent solution*: Any elements belong to the *k*th-distance of incumbent solution *x*, i.e. *Nk*(*x*), have a range in which the changed items and periods of

Therefore, if the changing result of an incumbent solution by an implement technique falls into a small range, it seems be a local search and hence may give little influence on the total cost function. Otherwise a change of incumbent solution within a large range may be a global search increasing the probability of worse solutions found and resulting lower searching efficiency. Comparing the considerations of constructing mutation in the algorithms above, we can find that only short distance (for cumulative mutation, it is a very special case of long distance) has been executed in GA, SA and PSO with different implement techniques. For example, all of the algorithms use mutation to change the solution but only GA uses the probability of obtaining a good solution to increase the share rate of the operators. While SA uses a probability of accepting a worse solution to create an uphill ability of escaping from the traps of local optima, and PSO uses other two factors (*pBest* and *gBest* ) to point a direction in which the optimal solution maybe exist. SOA and ACO are using basic random sampling principle to select the decisions in all positions of item and period. The difference is that SOA does not remain any information of previous solutions so large numbers of solution should be initialed and long distance (about half of items\*periods) always existed, whereas ACO uses a level by level approach to do single-bit mutation in a solution so goodness of previous solution may has been remained with long distance. Finally, VNS directly use distances to search good

Next, regarding the other three implement techniques, i.e., the inversion, the item crossover, and the period crossover, which are not used in other algorithms except GA, we can find that the inversion is just a special mutation that changes the solution with two points distance. However crossover(s) may change the original solutions with a longer distance so only partial goodness of previous solution has been remained. In GA algorithm, the implement techniques with long and short distance are combined together, therefore it should be considered more effective for solving the MLLS problem. However, we have to say that those implement techniques used in above algorithms (includes GA) seem not effective in finding good solution. Because even a good solution exists near by the original solution, here is no way to guarantee that it will be found by above algorithms since the candidate solution is randomly generated. It can be considered that not only distance but also range should be used in the *step* function to find a better solution. However, two problems here need to answer. Can distance and range evaluate the effective of implement techniques using in the MLLS problem? How do the

For proofing our consideration above, simulation experiments were executed by using the general algorithm for solving the MLLS problems shown in Fig.1. Initial solution(s) is pro‐ duced with a randomized cumulative Wagner and Whitin (RCWW) algorithm. The *step* function uses a GA algorithm, because all kinds of the implemental techniques had been included in the algorithm. We can control what the implemental technique should be used then evaluate the efficiency of the techniques. However the *terminate*functions used in different problem instances are different. In the experiments we first show the relationship among

incumbent solution have been limited.

86 Recent Advances on Meta-Heuristics and Their Application to Real Scenarios

solution so that its distance is controllable.

implement techniques should be applied?

Table 2 and 3 show an example of better solutions with different distances in different ranges. Better solution is defined as that the new solutions are better than the incumbent, so we count the number of better solutions and calculate the ratio of different distances. It can be observed from Table 2 that different distances lead to different number of better solutions found. Most of them (94.82%) are found in the nearest neighbourhood with one point distance and the ratio decreases as the distance increasing, which indicates a less probability of finding better solution among those whose distance are long from the incumbent. However, the incumbent solution did be improved by candidate solutions from longer distance. Even so the results with same tendency can be observed from Table 3 in which a limited range has been used. However very different meanings can be observed from Table 3. Firstly, limiting the search range can lead to a more efficiency for searching better solutions, which is represented by the fact that the total number of better solutions found is about 4% more than that of the non-limited search, and it also leads to obtain a solution with a better total cost (comparing with Table 2). Secondly, even the number of better solutions are almost same in distance 1 (changing just one position of the incumbent solution), but the number of better solutions in distance 2 was about three times of that of non-limited search. That is, longer time and less effect were performed in distance 1 if the range is not limited. This observation can lead a considerable result of implementing a faster searching in distance 1 and a faster changing between distances. That is the superiors of limited searching range.



que. We represent the relationships among distance, range and searching effectives in a three dimensional figure, in which the dimensions are the *distance*, the *range* and the *number of better solutions found*. Then the distribution in the figures can show how the criteria of distance and

A Comparative Study on Meta Heuristic Algorithms for Solving Multilevel Lot-Sizing Problems

http://dx.doi.org/10.5772/55279

89

Firstly, it can be considered that the single-bit mutation and the inversion are very similar since the candidate solutions are always in distance 1 and range ±0 by single-bit mutation and always in distance 2 and range ±1 by inversion, so the better solutions found are always in distance 1 and range ±0 for the former and distance 2 and range ±1 for the latter, which have been verified in Fig.2 and Fig.3. Comparing Fig.2 and Fig.3, we also can find that the searching effective of reversion is little better than that of single-bit mutation because more better solutions can be

Distribution of better solutions by single-bit mutation

<sup>0</sup> <sup>1</sup> <sup>2</sup> <sup>3</sup> <sup>4</sup>

Distribution of better solutions by reverse

<sup>0</sup> <sup>1</sup> <sup>2</sup> <sup>3</sup> <sup>4</sup> <sup>5</sup>

Range

Range

Number of better solutions

Number of better solutions

Dis

Dis

range segmenting the searching effectives.

found by reversion.

**Figure 2.** The result of single-bit mutation

**Figure 3.** The result of reverse

**Table 2.** A non limited search in distances in small MLLS problems (Parameters: *k*max=5, Δ*i* = ± 5, Δ*t* = ± 12)


**Table 3.** A limited search in distances in small MLLS problems (Parameters: *k*max=5, Δ*i* = ± 1, Δ*t* = ± 3)

That means even the implemental techniques used in algorithms such as SA, PSO, and VNS are seemly different, but the results from them may be totally similar. In addition, GA uses other implement techniques like crossover(s), so it may lead a longer distance and improve the searching performance basically. Moreover, distance and range have a very flexible controllability to produce candidate solutions. It gives a very important observation that a special range which includes some heuristic information (such as period heuristics is effective but level heuristic is not effective, and so on) can improve the performance of implemental technique, therefore they should be used as some new implemental methods into the metaheuristic algorithms.

#### **Remark 2**: effective of various implement techniques

Here five implemental techniques, i.e., single-bit mutation, cumulative mutation, inversion, product crossover, and period crossover, were tested by using GA algorithm. These implement techniques are used in the *step* function of our GA algorithm. The RCWW algorithm is used as the initial method to produce a population of 50 solutions as the first generation. Then, each of the five implemental techniques is selected with equal probability, namely, 20%, to produce a new generation, and each generation keeps only a population of 50 solutions after the selections by total cost function. The algorithm stops after 50 generations have been evolved. A solution is said to be a better solution if its objective cost is better than those of its parent and is counted in terms of changing distance and changing range from its parent. Therefore, simulation experiment is able to show the relationships among distance, range and searching effectives (represented by the number of better solutions found) of each implemental techni‐ que. We represent the relationships among distance, range and searching effectives in a three dimensional figure, in which the dimensions are the *distance*, the *range* and the *number of better solutions found*. Then the distribution in the figures can show how the criteria of distance and range segmenting the searching effectives.

Firstly, it can be considered that the single-bit mutation and the inversion are very similar since the candidate solutions are always in distance 1 and range ±0 by single-bit mutation and always in distance 2 and range ±1 by inversion, so the better solutions found are always in distance 1 and range ±0 for the former and distance 2 and range ±1 for the latter, which have been verified in Fig.2 and Fig.3. Comparing Fig.2 and Fig.3, we also can find that the searching effective of reversion is little better than that of single-bit mutation because more better solutions can be found by reversion.

**Figure 2.** The result of single-bit mutation

**Distance Better solutions Ratio** 3 28 0.80% 4 14 0.40% 5 5 0.14%

Total cost=826.23 3515 100.00%

**Distance Better solutions Ratio** 1 3232 88.57% 2 332 9.10% 3 55 1.51% 4 25 0.69% 5 5 0.14%

Total cost=824.37 3649 100%

That means even the implemental techniques used in algorithms such as SA, PSO, and VNS are seemly different, but the results from them may be totally similar. In addition, GA uses other implement techniques like crossover(s), so it may lead a longer distance and improve the searching performance basically. Moreover, distance and range have a very flexible controllability to produce candidate solutions. It gives a very important observation that a special range which includes some heuristic information (such as period heuristics is effective but level heuristic is not effective, and so on) can improve the performance of implemental technique, therefore they should be used as some new implemental methods into the meta-

Here five implemental techniques, i.e., single-bit mutation, cumulative mutation, inversion, product crossover, and period crossover, were tested by using GA algorithm. These implement techniques are used in the *step* function of our GA algorithm. The RCWW algorithm is used as the initial method to produce a population of 50 solutions as the first generation. Then, each of the five implemental techniques is selected with equal probability, namely, 20%, to produce a new generation, and each generation keeps only a population of 50 solutions after the selections by total cost function. The algorithm stops after 50 generations have been evolved. A solution is said to be a better solution if its objective cost is better than those of its parent and is counted in terms of changing distance and changing range from its parent. Therefore, simulation experiment is able to show the relationships among distance, range and searching effectives (represented by the number of better solutions found) of each implemental techni‐

**Table 3.** A limited search in distances in small MLLS problems (Parameters: *k*max=5, Δ*i* = ± 1, Δ*t* = ± 3)

heuristic algorithms.

**Remark 2**: effective of various implement techniques

**Table 2.** A non limited search in distances in small MLLS problems (Parameters: *k*max=5, Δ*i* = ± 5, Δ*t* = ± 12)

88 Recent Advances on Meta-Heuristics and Their Application to Real Scenarios

**Figure 3.** The result of reverse

<sup>0</sup> <sup>1</sup> <sup>2</sup> <sup>3</sup> <sup>4</sup> <sup>5</sup> <sup>6</sup> <sup>7</sup> <sup>8</sup> <sup>9</sup> <sup>10</sup>

Secondly, it seems the cumulative mutation may trigger longer distance and larger range (than the single-bit mutation) to the incumbent solutions, the results of the experiment shown in Fig.4 can illustrate the correctness of this observation. Also it can be observed that only those with small distance and range are more effective in improving the performance of cost

Thirdly, the implemental techniques of item crossover and period crossover are more complex than cumulative mutation and reversion, and the offspring consequently may be associated with very large range and very long distance in comparison to their parents. Still, it can be observed from our experimental results in Fig.5 and Fig.6 that only those offspring with small distance and range are more effective on improving the fitness function. Moreover, comparing with simple implemental techniques (single-bit mutation, reversion and cummulative mutation), crossover techniques can achieve significent performance improvements in

Finally in Fig.7, we show the mixed distribution of all the five implemental techniques. The results of simulation experiments show that a total signifecent effective of searching better solution can be obtained by using the mixed implemental techniques. Also repeatly, almost achievements of searching effective are in smaller distance and range. This phenominon may be used to conduct a more effective implemental technique in GA algorithms than pure free

However, for candidate solutions with larger range and longer distance, they still provided improvements on fitness function with probability so should not be ignored by any metaheuristic algorithms for optimality purpose. Nevertherless, we can develop more effective and more efficient implemental techniques by matching to the ratio of better solution can be found in different ranges and distances. For example, those crossover operations resulting in

offspring with too long distance are non-effective and should be avoided as possible.

Mixed distribution of better solutions by five techniques

Range

Number of better solutions

**Figure 7.** The mixed result of five implemental techniques

function.

searching effective.

crossover.

A Comparative Study on Meta Heuristic Algorithms for Solving Multilevel Lot-Sizing Problems

http://dx.doi.org/10.5772/55279

91

Dis

**Figure 4.** The result of cumulative mutation

**Figure 5.** The result of period crossover

**Figure 6.** The result of item crossover

**Figure 7.** The mixed result of five implemental techniques

<sup>0</sup> <sup>1</sup> <sup>2</sup> <sup>3</sup> <sup>4</sup>

Distribution of better solutions by period crossover

<sup>0</sup> <sup>1</sup> <sup>2</sup> <sup>3</sup> <sup>4</sup> <sup>5</sup> <sup>6</sup> <sup>7</sup> <sup>8</sup> <sup>9</sup>

Distribution of better solutions by item crossover

<sup>0</sup> <sup>1</sup> <sup>2</sup> <sup>3</sup> <sup>4</sup> <sup>5</sup> <sup>6</sup> <sup>7</sup> <sup>8</sup> <sup>9</sup> <sup>10</sup>

Range

Range

Distribution of better solutions by cumulative mutation

Range

Number of better solutions

Number of better solutions

Number of better solutions

90 Recent Advances on Meta-Heuristics and Their Application to Real Scenarios

**Figure 4.** The result of cumulative mutation

**Figure 5.** The result of period crossover

**Figure 6.** The result of item crossover

Dis

Dis

Dis

Secondly, it seems the cumulative mutation may trigger longer distance and larger range (than the single-bit mutation) to the incumbent solutions, the results of the experiment shown in Fig.4 can illustrate the correctness of this observation. Also it can be observed that only those with small distance and range are more effective in improving the performance of cost function.

Thirdly, the implemental techniques of item crossover and period crossover are more complex than cumulative mutation and reversion, and the offspring consequently may be associated with very large range and very long distance in comparison to their parents. Still, it can be observed from our experimental results in Fig.5 and Fig.6 that only those offspring with small distance and range are more effective on improving the fitness function. Moreover, comparing with simple implemental techniques (single-bit mutation, reversion and cummulative mutation), crossover techniques can achieve significent performance improvements in searching effective.

Finally in Fig.7, we show the mixed distribution of all the five implemental techniques. The results of simulation experiments show that a total signifecent effective of searching better solution can be obtained by using the mixed implemental techniques. Also repeatly, almost achievements of searching effective are in smaller distance and range. This phenominon may be used to conduct a more effective implemental technique in GA algorithms than pure free crossover.

However, for candidate solutions with larger range and longer distance, they still provided improvements on fitness function with probability so should not be ignored by any metaheuristic algorithms for optimality purpose. Nevertherless, we can develop more effective and more efficient implemental techniques by matching to the ratio of better solution can be found in different ranges and distances. For example, those crossover operations resulting in offspring with too long distance are non-effective and should be avoided as possible.

## **4. An example: A IVNS algorithm for solving MLLS problems**

The VNS/VND algorithm initiated by Mladenovic and Hansen (1997), Hansen and Mladenovic (2001) and Hansen et al. (2001, 2008), is a top-level methodology for solving the combinatorial optimization problems. Because its principle is simple and easily understood and implement‐ ed, it has been successfully applied to several optimization problems in different fields. The success of VNS is largely due to its enjoying many of the 11 desirable properties of metaheuristic generalized by Hansen et al., (2008), such as simplicity, user-friendliness, efficiency, effectiveness, and so on. Since the MLLS problem is observed to share common characteristics, e.g., a binary decision variable, with those problems successfully solved by VNS/VND based algorithm, it is promising to develop an efficient algorithm for this problem (Xiao et al., 2011a, 2011b, 2012). Here an iterated variable neighborhood descent (IVND) algorithm, which is a variant of VNS, is proposed as an example to show the utility and performance improve‐ ment of considerations descripted above.

additional inventory holding cost from its predecessors because not all the demands for predecessors are met by timely production; some of them may also be satisfied by inventory.

'

A Comparative Study on Meta Heuristic Algorithms for Solving Multilevel Lot-Sizing Problems

*j*

*H*

is the original unit inventory holding cost of product *i* and *q* is a random value

1

å (9)

'

http://dx.doi.org/10.5772/55279

can be

93

Therefore and similarly, the modified unit inventory holding cost of product *i*, i.e., *Hi*

1

*j j*

Here six implemental techniques are used in the IVND algorithm which are integrated together

**1.** *Limit the changing range of incumbent solution within one item*. Limit the changing range of incumbent solution within one item, i.e., *N*1(*x*), when exploring a neighborhood farther than *N*1(*x*). That is, when multiple mutations are done on incumbent solution, they must

**2.** *All mutations are limited between the first period and the last period that have positive demand*. This technique makes the search algorithm to avoid invalid mutations. Nevertheless, the

**3.** *No demand, no setup*. Mutation must not try of arranging setups for products in periods without positive demand, which is obviously non-optimal operation and should be

**4.** *Triggerrecursive mutations*. A mutation of canceling a setup for a product in a period will trigger recursive mutations on all of its ancestors. While a mutation of canceling a setup occurs, e.g. changing the value of bit *xit* from 1 to 0, it withdraws demands for the immediate predecessors in previous periods of leading time. As a consequence, some of these predecessors their demands may drop to zero such that their setups (if they have) in these periods should be canceled at the same time; other predecessors who remain nonzero demands due to product commonality should remain unchanged. The recursive mutations are only triggered for cancellation of a setup for production; they will not occur

**5.** *Shift a setup rather than cancel a setup*. When the setup for product *i* at period *t* need to be canceled, try to shift the setup to the first period with positive demand behind *t*, rather than simply cancel it. For example, when *xit* is to be canceled, find the first period *t\** behind *t* of product *i* that satisfies *Dit\** >0 and arrange a setup by setting *xit\**←1 if *xit\** is 0. Notably,

this arrangement of this setup will also trigger recursive mutations.


*i*

*ii j*

*HHq H*

'

to deliver good performance in the computing experiments.

first period with positive demand should be fixed with 1.

occur on same item but different periods.

banned in the whole searching process.

when arranging a setup.

calculated recursively by

uniformly distributed between [0,1].

**4.2. Implemental techniques**

where *Hi*

#### **4.1. Initial method**

We use a randomized cumulative Wagner and Whitin (RCWW) based approach to initialize the solution for our proposed IVNS algorithm. The RCWW method was introduced by Dellaert and Jeunet(2000a), of which the main idea is based on the fact that lot-sizing a product in the one period will trigger demands for its components in previous periods with leading time corrections( or in same period for the case of zero leading time). Therefore, the real setup cost of an item is in fact greater than its own setup cost and should be modified when using the wellknown sequentiail WW algorithm to generate a initial solution. The time-varying modified setup cost is a improved concept introduced by Dellaert et al(2000b, 2003) and used by Pitakaso et al(2007) which dispos‐ es of using a constant modified setup cost for the whole planning horizon; it suggested the costs might vary from one time period to the next, and reported good in getting better initial solutions.

In the IVND algorithm, we use the sequential WW algorithm based on randomly initialized constant modified setup cost to generate initial solutions. For each item *i*, its modified setup cost *Si* ' can be calculated recursively by

$$\boldsymbol{\mathcal{S}}\_{i} = \boldsymbol{\mathcal{S}}\_{i} + r \left[ \sum\_{j \in \Gamma\_{i}^{-1}} \left( \boldsymbol{\mathcal{S}}\_{j} + \frac{\boldsymbol{\mathcal{S}}\_{j}^{\cdot}}{\left| \Gamma\_{j}^{-1} \right|} \right) \right] \tag{8}$$

where *Si* is the original setup cost of item *i*, *r* is a random value uniformly distributed between [0,1], | *Γ<sup>i</sup>* <sup>−</sup><sup>1</sup> | is the set of immediate predecessors(components) of product *i* and | *<sup>Γ</sup><sup>i</sup>* <sup>−</sup><sup>1</sup> |is its cardinality.

In addition to the modified setup cost, we also use the modified unit inventory holding cost to construct initial solutions. It is simply based on the consideration that a positive inventory balance of one product in one period causes not only its own inventory holding cost but also additional inventory holding cost from its predecessors because not all the demands for predecessors are met by timely production; some of them may also be satisfied by inventory. Therefore and similarly, the modified unit inventory holding cost of product *i*, i.e., *Hi* ' can be calculated recursively by

$$H\_i^\cdot = H\_i + q \left[ \sum\_{j \neq \Gamma\_i^{-1}} \left( H\_j + \frac{H\_j^\cdot}{\left| \Gamma\_j^{-1} \right|} \right) \right] \tag{9}$$

where *Hi* is the original unit inventory holding cost of product *i* and *q* is a random value uniformly distributed between [0,1].

## **4.2. Implemental techniques**

**4. An example: A IVNS algorithm for solving MLLS problems**

92 Recent Advances on Meta-Heuristics and Their Application to Real Scenarios

ment of considerations descripted above.

can be calculated recursively by

'

**4.1. Initial method**

cost *Si* '

where *Si*

[0,1], | *Γ<sup>i</sup>*

cardinality.

The VNS/VND algorithm initiated by Mladenovic and Hansen (1997), Hansen and Mladenovic (2001) and Hansen et al. (2001, 2008), is a top-level methodology for solving the combinatorial optimization problems. Because its principle is simple and easily understood and implement‐ ed, it has been successfully applied to several optimization problems in different fields. The success of VNS is largely due to its enjoying many of the 11 desirable properties of metaheuristic generalized by Hansen et al., (2008), such as simplicity, user-friendliness, efficiency, effectiveness, and so on. Since the MLLS problem is observed to share common characteristics, e.g., a binary decision variable, with those problems successfully solved by VNS/VND based algorithm, it is promising to develop an efficient algorithm for this problem (Xiao et al., 2011a, 2011b, 2012). Here an iterated variable neighborhood descent (IVND) algorithm, which is a variant of VNS, is proposed as an example to show the utility and performance improve‐

We use a randomized cumulative Wagner and Whitin (RCWW) based approach to initialize the solution for our proposed IVNS algorithm. The RCWW method was introduced by Dellaert and Jeunet(2000a), of which the main idea is based on the fact that lot-sizing a product in the one period will trigger demands for its components in previous periods with leading time corrections( or in same period for the case of zero leading time). Therefore, the real setup cost of an item is in fact greater than its own setup cost and should be modified when using the wellknown sequentiail WW algorithm to generate a initial solution. The time-varying modified setup cost is a improved concept introduced by Dellaert et al(2000b, 2003) and used by Pitakaso et al(2007) which dispos‐ es of using a constant modified setup cost for the whole planning horizon; it suggested the costs might vary from one time period to the next, and reported good in getting better initial solutions.

In the IVND algorithm, we use the sequential WW algorithm based on randomly initialized constant modified setup cost to generate initial solutions. For each item *i*, its modified setup

'

*j*

*S*

is the original setup cost of item *i*, *r* is a random value uniformly distributed between

1

å (8)

<sup>−</sup><sup>1</sup> |is its

1

*j j*


<sup>−</sup><sup>1</sup> | is the set of immediate predecessors(components) of product *i* and | *<sup>Γ</sup><sup>i</sup>*

In addition to the modified setup cost, we also use the modified unit inventory holding cost to construct initial solutions. It is simply based on the consideration that a positive inventory balance of one product in one period causes not only its own inventory holding cost but also

*i*

*i i j*

*SSr S*

Here six implemental techniques are used in the IVND algorithm which are integrated together to deliver good performance in the computing experiments.


**6.** *Only affected products and periods are recalculated of their inventory holding costs and setup costs.* Different to other evolutionary algorithms like GA and PSO where a group of incumbent solutions have to be maintained, the IVND algorithm has only one incumbent solution. In fact, when a mutation occurs, most area of the solution states including setup matrix *Yit*, lot-sizing matrix *Xit* and demand matrix *Dit* are remain unchanged. Thus, it just needs to recalculate the affected products and the affected periods of the setup cost and inventory holding cost after a mutation operation. By taking this advantage, the comput‐ ing efficiency of IVND algorithm can be significantly improved since the recalculation of the objective function--the most time-consuming part of IVND algorithm, are avoided.

process loops until the stop condition is met. The stop condition can be a user-specified computing time or a maximum span between two restarts without improvement on the best solution found. In our experiments of the next section, we use the later one, i.e., a fixed times of continuous restarts without improvement, as the stop condition. The basic scheme the

A Comparative Study on Meta Heuristic Algorithms for Solving Multilevel Lot-Sizing Problems

http://dx.doi.org/10.5772/55279

95

Define the set of neighborhood structures *Nk*, *k*=1,…,*k*max, that will be used in the search; choose a stop condition.

There are three parameters, i.e., *P*, *N*, and *K*max, in the IVND algorithm for determining the tradeoff between searching efficiency and the quality of final solution. The first parameter *P* is a positive number which serves as a stop condition indicating the maximum span between two restarts without improvement on best solution found. The second parameter *N* is the maximum number of explorations between two improvements within a neighborhood. If a better solution cannot be found after *N* times of explorations in the neighborhood *Nk*(*x*), it is then deemed as explored and the algorithm goes to explore a farther neighborhood by *k*←*k*+1. The third parameter *K*max is the traditional parameter for VND search indicating the farthest

Three sets of MLLS problem instances under different scales(small, medium and large) are used to test the performance of the proposed IVND algorithm. The first set consists of 96 smallsized MLLS problems involving 5-item assembly structure over a 12-period planning horizon, which was developed by Coleman and McKnew (1991) on the basis of work by Veral and LaForge (1985) and Benton and Srivastava (1985), and also used by Dellaert and Jeunet (2000a). In the 96 small-sized problems, four typical product structures with an one-to-one production ratio are considered, and the lead times of all items are zero. For each product structure, four cost combinations are considered, which assign each individual item with different setup costs and different unit holding costs. Six independent demand patterns with variations to reflect low, medium and high demand are considered over a 12-period planning horizon. Therefore, these combinations produce 4×4×6=96 problems for testing. The optimal

(b)Move or not: if *x'* is better than *x*, then *x x* ' , *k*←1 and *n*←0, go to step (a) ; otherwise, *n*← *n*+1;

proposed IVND algorithm is illustrated in Fig. 9.

(1)Find an initial solution *x*0 by using RCWW algorithm

(3)Until *k*=*k*max repeat (a), (b), and (c) steps: (a)Find at random a solution *x'* in *Nk*(*x*);

(4)Output the best solution found.

**Figure 9.** The IVND algorithm for MLLS problem

neighborhood that the algorithm will go.

**5.1. Problem instances**

**5. Computational experiments and the results**

(2)Set *k*←1, *n*← 0;

Repeat the following step (1), (2) and (3) until the stop condition is met:

(c) If *n*=*N*, then shift to search a farther neighborhood by *k*←*k*+1 and reset *n*←0;

The above six implemental techniques are all used in our proposed IVND algorithm to mutate the incumbent solution into its neighborhood. Steps of implementing these techniques on neighborhood search, e.g., neighborhood *Nk*(*x*), can be explained by Fig.8.

To generate a candidate solution from neighborhood *Nk*(*x*) of the incumbent solution *x*:

(1)Select randomly *k* bits of *<sup>x</sup>*, e.g., 1 2 ,, , , ,..., *<sup>k</sup> it it it xx x* , and sequentially mutate them. The first three implement techniques mentioned above must follow in the selection: the first implemental technique that the to-be bits must be within an identical item; the second implemental technique that the to-be mutated periods should between the first period and last period with positive demand; the third implemental technique that those periods without demand should not be selected for mutation.

(2)For each mutation from 1 to 0(noticeably not including those from 0 to 1), the forth implemental technique must be followed to trigger recursive mutations toward its predecessors with zero demand. In the recursive mutation process, the fifth implemental technique must be implemented to try of shifting the setup to the first sequential period with positive demand, rather than simply removed it.

(3)Whenever a bit of the incumbent solution is changed, the sixth implemental technique is implemented to recalculate the objective function just by recalculating the affected items and their periods.

**Figure 8.** The implementation of implemental techniques

Although the new solutions from *Nk*(*x*) may has a greater than *k* unit distance from the incumbent solution *x* after implemented with these six implemental techniques, it is still considered as a member of *Nk*(*x*). These implemental techniques are only deemed as additional actions implemented on the new solution toward better solution. Moreover, benefiting from these actions, higher efficiency of VNS algorithm could be consequently anticipated, which has been confirmed in the experiments of Section 4.

#### **4.3. The IVND algorithm**

The algorithm IVND is a variant of the basic VND algorithm. It starts from initiating a solution as the incumbent solution, and then launches a VND search. The VND search repeatedly tries of finding a better solution in the nearby neighborhood of the incumbent solution and moves to the better solution found; if a better solution cannot be found in current neighborhood, then go to explore a father neighborhood until the farthest neighborhood is reached. Once the VND process is stopped(characterized by the farthest neighborhood been explored), another initial solution is randomly generated and restarts the VND search again. This simply iterated search process loops until the stop condition is met. The stop condition can be a user-specified computing time or a maximum span between two restarts without improvement on the best solution found. In our experiments of the next section, we use the later one, i.e., a fixed times of continuous restarts without improvement, as the stop condition. The basic scheme the proposed IVND algorithm is illustrated in Fig. 9.

Define the set of neighborhood structures *Nk*, *k*=1,…,*k*max, that will be used in the search; choose a stop condition. Repeat the following step (1), (2) and (3) until the stop condition is met: (1)Find an initial solution *x*0 by using RCWW algorithm (2)Set *k*←1, *n*← 0; (3)Until *k*=*k*max repeat (a), (b), and (c) steps: (a)Find at random a solution *x'* in *Nk*(*x*); (b)Move or not: if *x'* is better than *x*, then *x x* ' , *k*←1 and *n*←0, go to step (a) ; otherwise, *n*← *n*+1; (c) If *n*=*N*, then shift to search a farther neighborhood by *k*←*k*+1 and reset *n*←0; (4)Output the best solution found.

**Figure 9.** The IVND algorithm for MLLS problem

**6.** *Only affected products and periods are recalculated of their inventory holding costs and setup costs.* Different to other evolutionary algorithms like GA and PSO where a group of incumbent solutions have to be maintained, the IVND algorithm has only one incumbent solution. In fact, when a mutation occurs, most area of the solution states including setup matrix *Yit*, lot-sizing matrix *Xit* and demand matrix *Dit* are remain unchanged. Thus, it just needs to recalculate the affected products and the affected periods of the setup cost and inventory holding cost after a mutation operation. By taking this advantage, the comput‐ ing efficiency of IVND algorithm can be significantly improved since the recalculation of the objective function--the most time-consuming part of IVND algorithm, are avoided.

The above six implemental techniques are all used in our proposed IVND algorithm to mutate the incumbent solution into its neighborhood. Steps of implementing these techniques on

(1)Select randomly *k* bits of *<sup>x</sup>*, e.g., 1 2 ,, , , ,..., *<sup>k</sup> it it it xx x* , and sequentially mutate them. The first three implement techniques mentioned above must follow in the selection: the first implemental technique that the to-be bits must be within an identical item; the second implemental technique that the to-be mutated periods should between the first period and last period with positive demand; the third implemental technique that those periods without demand

(2)For each mutation from 1 to 0(noticeably not including those from 0 to 1), the forth implemental technique must be followed to trigger recursive mutations toward its predecessors with zero demand. In the recursive mutation process, the fifth implemental technique must be implemented to try of shifting the setup to the first sequential period

(3)Whenever a bit of the incumbent solution is changed, the sixth implemental technique is implemented to

Although the new solutions from *Nk*(*x*) may has a greater than *k* unit distance from the incumbent solution *x* after implemented with these six implemental techniques, it is still considered as a member of *Nk*(*x*). These implemental techniques are only deemed as additional actions implemented on the new solution toward better solution. Moreover, benefiting from these actions, higher efficiency of VNS algorithm could be consequently anticipated, which

The algorithm IVND is a variant of the basic VND algorithm. It starts from initiating a solution as the incumbent solution, and then launches a VND search. The VND search repeatedly tries of finding a better solution in the nearby neighborhood of the incumbent solution and moves to the better solution found; if a better solution cannot be found in current neighborhood, then go to explore a father neighborhood until the farthest neighborhood is reached. Once the VND process is stopped(characterized by the farthest neighborhood been explored), another initial solution is randomly generated and restarts the VND search again. This simply iterated search

neighborhood search, e.g., neighborhood *Nk*(*x*), can be explained by Fig.8.

To generate a candidate solution from neighborhood *Nk*(*x*) of the incumbent solution *x*:

94 Recent Advances on Meta-Heuristics and Their Application to Real Scenarios

recalculate the objective function just by recalculating the affected items and their periods.

should not be selected for mutation.

**4.3. The IVND algorithm**

with positive demand, rather than simply removed it.

**Figure 8.** The implementation of implemental techniques

has been confirmed in the experiments of Section 4.

There are three parameters, i.e., *P*, *N*, and *K*max, in the IVND algorithm for determining the tradeoff between searching efficiency and the quality of final solution. The first parameter *P* is a positive number which serves as a stop condition indicating the maximum span between two restarts without improvement on best solution found. The second parameter *N* is the maximum number of explorations between two improvements within a neighborhood. If a better solution cannot be found after *N* times of explorations in the neighborhood *Nk*(*x*), it is then deemed as explored and the algorithm goes to explore a farther neighborhood by *k*←*k*+1. The third parameter *K*max is the traditional parameter for VND search indicating the farthest neighborhood that the algorithm will go.

## **5. Computational experiments and the results**

## **5.1. Problem instances**

Three sets of MLLS problem instances under different scales(small, medium and large) are used to test the performance of the proposed IVND algorithm. The first set consists of 96 smallsized MLLS problems involving 5-item assembly structure over a 12-period planning horizon, which was developed by Coleman and McKnew (1991) on the basis of work by Veral and LaForge (1985) and Benton and Srivastava (1985), and also used by Dellaert and Jeunet (2000a). In the 96 small-sized problems, four typical product structures with an one-to-one production ratio are considered, and the lead times of all items are zero. For each product structure, four cost combinations are considered, which assign each individual item with different setup costs and different unit holding costs. Six independent demand patterns with variations to reflect low, medium and high demand are considered over a 12-period planning horizon. Therefore, these combinations produce 4×4×6=96 problems for testing. The optimal solutions of 96 benchmark problem are previously known so that can serve as benchmark for testing against the optimality of the new algorithm.

the column *Best solutions found(%)*. The column *average time(s)* indicates the average computing time in second of one run for each problem. The average result of 10 and the minimum result of 10 are both shown in Table 4 and compared to the HGA algorithm of Dellaert and Jeunet (2000a), the MMAS algorithm of Pitakaso et al.(2007) and the PGA algorithm of Homberger

solutions of 96 benchmark problems. Although the PGAC and the GA3000 can also find 100% optimal solutions, they take longer computing time and also take the advantage of hardware (for PGAC 30 processors were used to make a parallel calculation). After that, we adjust the

Surprisingly, we got 960 optimal solutions (100% optimality) with computing time of 0.27

**Average**

HGA50 810.74 96..88 0.26 5.14s -- 1 Dellaert et al.(2000) MMAS 810.79 92.71 0.26 <1s P4 1.5G 1 Pitakaso et al.(2007) GA100 811.98 75.00 0.68 5s P4 2.4G 1 Homberger(2008) **GA3000 810.67 100.00 0.00 5s P4 2.4G 1 Homberger(2008)** PGAI 810.81 94.79 0.28 5s P4 2.4G 30 Homberger(2008) **PGAC 810.67 100.00 0.00 5s P4 2.4G 30 Homberger(2008)**

100 (Avg. of 10) 810.69 99.58 0.02 0.07s NB 1.6G <sup>1</sup> IVND

**<sup>100</sup>** (Min. Of 10) **810.67 100.00 0.00 0.7s NB 1.6G <sup>1</sup> IVND**

**200(Avg. Of 10) 810.67 100.00 0.00 0.27s NB 1.6G <sup>1</sup> IVND**

**OPT. 810.67 100.00 0.00 -- -- - --**

repeated 10 runs. We summarize the 400 results and compare them with the existing algo‐ rithms in Table 5. More detailed results of 40 problems are listed in Table 6 where the previous best known solutions summarized by Homberger(2008) are also listed for comparison. After

for these 40 medium-sized problems which are listed in column *new best solution* in Table 6.

100 algorithm on 40 medium-sized MLLS benchmark problems and

100 algorithm for several times and updated the best solutions

**time(s) CPU type**

A Comparative Study on Meta Heuristic Algorithms for Solving Multilevel Lot-Sizing Problems

**Mean dev. if best solution not found**

50 uses 0.7 second to find 100% optimal

http://dx.doi.org/10.5772/55279

97

200 on the 96 problems for 10 times again.

**Number of processors**

**Sources**

(2008). It can be observed from Table 4 that IVND200

parameter *P* from 50 to 200 and repeatedly ran IVND200

**Best solutions found (%)**

**Table 4.** Comparing IVND with existing algorithms on 96 small-sized problems

**5.4. Medium-sized MLLS problems**

Secondly, we ran IVND600

that, we repeatedly ran IVND600

second.

IVND 50

**IVND 50**

**IVND 200**

**Method Avg.**

**cost**

The second set consists of 40 medium-sized MLLS problems involving 40/50-item product structure over a 12/24-period planning horizon, which are based on the product structures published by Afentakis et al. (1984), Afentakis and Gavish (1986), and Dellaert and Jeunet (2000). In the 40 medium-sized problems, four product structures with an one-to-one produc‐ tion ratio are constructed. Two of them are 50-item assembly structures with 5 and 9 levels, respectively. The other two are 40-item general structure with 8 and 6 levels, respectively. All lead times were set to zero. Two planning horizons were used: 12 and 24 periods. For each product structure and planning horizon, five test problems were generated, such that a total number of 4×2×5=40 problems could be used for testing.

The third set covers the 40 problem instances with a problem size of 500 products and 36/52 periods synthesized by Dellaert and Jeunet (2000). There are 20 different product structures with one-to-one production ratio and different commonality indices1 . The first 5 instances are pure assembly structures with one end-product. The instances from 6 to 20 are all general structure with five end-products and different communality indices ranges from 1.5 to 2.4. The first 20 instances are all over a 36-period planning horizon; the second 20 instances are of the same product structures of the first 20 instances but over a 52-period planning horizon. The demands are different for each instances and only on end-products.

Since the hybrid GA algorithm developed by Dellaert and Jeunet (2000a) is the first metaheuristic algorithm for solving the MLLS problem, it was always selected as a benchmark algorithm for comparison with newly proposed algorithm. Therefore, we compared the performance of our IVND algorithm with the hybrid GA algorithm on the all instances of three different scales. We also compared our IVND algorithm with the MMAS algorithm developed by Pitakaso et al.(2007), and the parallel GA algorithm developed by Homberger (2008) since both of them used the same three set of instances used in this paper.

## **5.2. Test environment**

The IVND algorithm under examination were coded in VC++6.0 and ran on a notebook computer equipped with a 1.6G CPU operating under Windows XP system. We fixed the parameter *K*max to be 5 for all experiments, and let the parameter *P* and *N* changeable to fit for the different size of problem. The symbol 'IVNDN P ' specifies the IVND algorithm with the parameter *P* and *N*, e.g., IVND200 50 indicates *P*=50, *N*=200, and *Kmax*=5 by default. The effect of individual parameter on the quality of solution was tested in section 5.6.

## **5.3. Small-sized MLLS problems**

We repeatedly ran IVND200 50 on the 96 small-sized MLLS problems for 10 times and got 960 results among which 956 were the optimal results so the optimality was 99.58% indicated by

<sup>1</sup> Commonality index is the average number of successors of all items in a product structure

the column *Best solutions found(%)*. The column *average time(s)* indicates the average computing time in second of one run for each problem. The average result of 10 and the minimum result of 10 are both shown in Table 4 and compared to the HGA algorithm of Dellaert and Jeunet (2000a), the MMAS algorithm of Pitakaso et al.(2007) and the PGA algorithm of Homberger (2008). It can be observed from Table 4 that IVND200 50 uses 0.7 second to find 100% optimal solutions of 96 benchmark problems. Although the PGAC and the GA3000 can also find 100% optimal solutions, they take longer computing time and also take the advantage of hardware (for PGAC 30 processors were used to make a parallel calculation). After that, we adjust the parameter *P* from 50 to 200 and repeatedly ran IVND200 200 on the 96 problems for 10 times again. Surprisingly, we got 960 optimal solutions (100% optimality) with computing time of 0.27 second.


**Table 4.** Comparing IVND with existing algorithms on 96 small-sized problems

#### **5.4. Medium-sized MLLS problems**

solutions of 96 benchmark problem are previously known so that can serve as benchmark for

The second set consists of 40 medium-sized MLLS problems involving 40/50-item product structure over a 12/24-period planning horizon, which are based on the product structures published by Afentakis et al. (1984), Afentakis and Gavish (1986), and Dellaert and Jeunet (2000). In the 40 medium-sized problems, four product structures with an one-to-one produc‐ tion ratio are constructed. Two of them are 50-item assembly structures with 5 and 9 levels, respectively. The other two are 40-item general structure with 8 and 6 levels, respectively. All lead times were set to zero. Two planning horizons were used: 12 and 24 periods. For each product structure and planning horizon, five test problems were generated, such that a total

The third set covers the 40 problem instances with a problem size of 500 products and 36/52 periods synthesized by Dellaert and Jeunet (2000). There are 20 different product structures

pure assembly structures with one end-product. The instances from 6 to 20 are all general structure with five end-products and different communality indices ranges from 1.5 to 2.4. The first 20 instances are all over a 36-period planning horizon; the second 20 instances are of the same product structures of the first 20 instances but over a 52-period planning horizon. The

Since the hybrid GA algorithm developed by Dellaert and Jeunet (2000a) is the first metaheuristic algorithm for solving the MLLS problem, it was always selected as a benchmark algorithm for comparison with newly proposed algorithm. Therefore, we compared the performance of our IVND algorithm with the hybrid GA algorithm on the all instances of three different scales. We also compared our IVND algorithm with the MMAS algorithm developed by Pitakaso et al.(2007), and the parallel GA algorithm developed by Homberger (2008) since

The IVND algorithm under examination were coded in VC++6.0 and ran on a notebook computer equipped with a 1.6G CPU operating under Windows XP system. We fixed the parameter *K*max to be 5 for all experiments, and let the parameter *P* and *N* changeable to fit for

results among which 956 were the optimal results so the optimality was 99.58% indicated by

P '

. The first 5 instances are

specifies the IVND algorithm with the

50 indicates *P*=50, *N*=200, and *Kmax*=5 by default. The effect of

50 on the 96 small-sized MLLS problems for 10 times and got 960

testing against the optimality of the new algorithm.

96 Recent Advances on Meta-Heuristics and Their Application to Real Scenarios

number of 4×2×5=40 problems could be used for testing.

with one-to-one production ratio and different commonality indices1

demands are different for each instances and only on end-products.

both of them used the same three set of instances used in this paper.

individual parameter on the quality of solution was tested in section 5.6.

1 Commonality index is the average number of successors of all items in a product structure

the different size of problem. The symbol 'IVNDN

**5.2. Test environment**

parameter *P* and *N*, e.g., IVND200

**5.3. Small-sized MLLS problems**

We repeatedly ran IVND200

Secondly, we ran IVND600 100 algorithm on 40 medium-sized MLLS benchmark problems and repeated 10 runs. We summarize the 400 results and compare them with the existing algo‐ rithms in Table 5. More detailed results of 40 problems are listed in Table 6 where the previous best known solutions summarized by Homberger(2008) are also listed for comparison. After that, we repeatedly ran IVND600 100 algorithm for several times and updated the best solutions for these 40 medium-sized problems which are listed in column *new best solution* in Table 6.


**Instance IVND400 Prev. best**

**Solution**

A Comparative Study on Meta Heuristic Algorithms for Solving Multilevel Lot-Sizing Problems

**S D I P Avg. of 10 Min. of 10 method**

 2 2 50 12 **155,938 155,938 155,938 155,938** B&B 2 3 50 12 **183,219 183,219 183,219 183,219** B&B 2 4 50 12 136,474 **136,462 136,462 136,462** B&B 2 5 50 12 186,645 **186,597 186,597 186,597** B&B 3 1 40 12 **148,004 148,004 148,004 148,004** PGAC 3 2 40 12 197,727 **197,695 197,695 197,695** PGAC 3 3 40 12 **160,693 160,693 160,693 160,693** MMAS 3 4 40 12 **184,358 184,358 184,358 184,358** PGAC 3 5 40 12 **161,457 161,457 161,457 161,457** PGAC 4 1 40 12 185,507 **185,170 185,170 185,161** PGAC→IVND 4 2 40 12 **185,542 185,542 185,542 185,542** PGAC 4 3 40 12 192,794 192,794 **192,157 192,157** MMAS 4 4 40 12 136,884 **136,757 136,764 136,757** PGAC→IVND 4 5 40 12 166,180 166,122 **166,041 166,041** PGAC 1 6 50 24 344,173 343,855 **343,207 343,207** PGAC 1 7 50 24 293,692 293,373 **292,908 292,908** HGA 1 8 50 24 356,224 355,823 **355,111 355,111** HGA 1 9 50 24 325,701 **325,278 325,304 325,278** PGAC 1 10 50 24 386,322 **386,059 386,082 385,954** HGA→IVND 2 6 50 24 341,087 341,033 **340,686 340,686** HGA 2 7 50 24 378,876 **378,845 378,845 378,845** HGA 2 8 50 24 346,615 346,371 **346,563 346,358** HGA→IVND 2 9 50 24 413,120 412,511 **411,997 411,997** HGA 2 10 50 24 390,385 **390,233 390,410 390,233** PGCA→IVND 3 6 40 24 **344,970 344,970 344,970 344,970** HGA 3 7 40 24 352,641 **352,634 352,634 352,634** PGAC 3 8 40 24 356,626 356,456 **356,427 356,323** PGAC→IVND 3 9 40 24 411,565 **411,438 411,509 411,438** MMAS→IVND 3 10 40 24 **401,732 401,732 401,732 401,732** HGA 4 6 40 24 289,935 **289,846 289,883 289,846** PGAC→IVND 4 7 40 24 339,548 339,299 **337,913 337,913** MMAS 4 8 40 24 320,920 320,426 **319,905 319,905** PGCA

**New best solution**

**New**

99

http://dx.doi.org/10.5772/55279

**Table 5.** Comparing IVND with existing algorithms on 40 medium-sized problems

It can be observed from Table 5 that the PGAC and IVND600 <sup>100</sup> (Min. of 10) algorithm are among the best and very close to each other. Although the optimality of PGAC (65%) is better than that of IVND600 <sup>100</sup> (Min. of 10) (60%) in terms of the baseline of previous best known solutions, it may drop at least 17% if in terms of new best solutions since 12 of 40 problems had been updated their best known solutions by IVND algorithm(see the column *new best solution* in Table 6) and 7 of the 12 updated problems' previous best known solution were previously obtained by PGAC. Furthermore, by taking account into consideration of hardware advantage of the PGAC algorithm(multiple processors and higher CPU speed), we can say that the IVND algorithm performances at least as best as the PGAC algorithm on medium-sized problems, if not better than.


A Comparative Study on Meta Heuristic Algorithms for Solving Multilevel Lot-Sizing Problems http://dx.doi.org/10.5772/55279 99


**Method Avg. cost**

98 Recent Advances on Meta-Heuristics and Their Application to Real Scenarios

IVND 100

IVND 100

that of IVND600

not better than.

**Optimality on previous best solutions (%)**

**Comp.**

HGA250\* 263,931.8 17.50 <60s -- 1 Dellaert et all(2000) MMAS 263,796.3 22.50 <20s P4 1.5G 1 Pitakaso et al.(2007) GA100 271,268.2 0.00 60s P4 2.4G 1 Homberger(2008) GA3000 266,019.8 15.00 60s P4 2.4G 1 Homberger(2008) PGAI 267,881.4 0.00 60s P4 2.4G 30 Homberger(2008) PGAC 263,359.6 65.00 60s P4 2.4G 30 Homberger(2008)

600(Avg. of 10) 263,528.5 30.00 2.67 NB 1.6G <sup>1</sup> IVND

600 (Min. of 10) 263,398.8 60.00 26.7 NB 1.6G <sup>1</sup> IVND

the best and very close to each other. Although the optimality of PGAC (65%) is better than

it may drop at least 17% if in terms of new best solutions since 12 of 40 problems had been updated their best known solutions by IVND algorithm(see the column *new best solution* in Table 6) and 7 of the 12 updated problems' previous best known solution were previously obtained by PGAC. Furthermore, by taking account into consideration of hardware advantage of the PGAC algorithm(multiple processors and higher CPU speed), we can say that the IVND algorithm performances at least as best as the PGAC algorithm on medium-sized problems, if

**Instance IVND400 Prev. best**

<sup>100</sup> (Min. of 10) (60%) in terms of the baseline of previous best known solutions,

**Solution**

**S D I P Avg. of 10 Min. of 10 method**

 1 1 50 12 196,084 196,084 **194,571 194,571** B&B 1 2 50 12 165,682 165,632 **165,110 165,110** B&B 1 3 50 12 **201,226 201,226 201,226 201,226** B&B 1 4 50 12 188,010 188,010 **187,790 187,790** B&B 1 5 50 12 161,561 161,561 **161,304 161,304** B&B 2 1 50 12 **179,761 179,761 179,762 179,761** B&B

**New best solution**

**New**

Prev. best solution 263,199.8 -- -- -- -- -- New best solution 260,678.3 -- -- -- - --

**Table 5.** Comparing IVND with existing algorithms on 40 medium-sized problems

It can be observed from Table 5 that the PGAC and IVND600

**Time(s) CPU type**

**Number of processors**

<sup>100</sup> (Min. of 10) algorithm are among

**Sources**


time to finish their calculation because the interdependencies among items are relatively high for these four problems. The column *Inter D.* in Table 8 is the maximum number of affected

> **Time (m)**

 500 36 500 597,940 597,560 6.3 **595,792 595,792** PGAC 500 36 500 817,615 816,507 6.8 **816,058 816,058** HGA 500 36 500 929,097 927,860 6.4 **911,036 911,036** PGAC 500 36 500 945,317 944,626 6.2 **942,650 942,650** MMAS 500 36 500 1,146,946 **1,145,980** 6.5 **1,149,005 1,145,980** MMAS→IVND 500 36 11036 7,725,323 **7,689,434** 71.9 **7,812,794 7,689,434** PGAC→IVND 500 36 3547 3,928,108 **3,923,336** 22.5 **4,063,248 3,923,336** MMAS→IVND 500 36 1034 2,724,472 2,713,496 17.2 **2,704,332 2,703,004** HGA→IVND 500 36 559 1,898,263 1,886,812 10.3 **1,943,809 1,865,141** PGAC→IVND 500 36 341 1,511,600 1,505,392 6.0 **1,560,030 1,502,371** MMAS→IVND 500 36 193607 59,911,520 **59,842,858** 179.9 **59,866,085 59,842,858** PGAC→IVND 500 36 22973 13,498,853 **13,441,692** 58.6 **13,511,901 13,441,692** PGAC→IVND 500 36 3247 4,751,464 **4,731,818** 34.4 **4,828,331 4,731,818** PGAC→IVND 500 36 914 2,951,232 2,937,914 18.6 **2,910,203 2,910,203** HGA 500 36 708 1,759,976 1,750,611 8.9 **1,791,700 1,740,397** MMAS→IVND 500 36 1099608 472,625,159 472,088,128 106.1 **471,325,517 471,325,517** PGAC 500 36 24234 18,719,243 **18,703,573** 80.5 **18,750,600 18,703,573** MMAS→IVND 500 36 7312 7,321,985 **7,292,340** 33.7 **7,602,730 7,292,340** MMAS→IVND 500 36 1158 3,592,086 **3,550,994** 23.4 **3,616,968 3,550,994** PGAC→IVND 500 36 982 2,326,390 2,293,131 13.8 **2,358,460 2,291,093** MMAS→IVND 500 52 500 1,189,599 1,188,210 22.4 **1,187,090 1,187,090** MMAS 500 52 500 1,343,567 **1,341,412** 14.9 **1,341,584 1,341,412** HGA→IVND 500 52 500 1,403,822 1,402,818 8.7 **1,400,480 1,384,263** MMAS→IVND 500 52 500 1,386,667 1,384,263 9.1 **1,382,150 1,382,150** MMAS 500 52 500 1,660,879 **1,658,156** 8.8 **1,660,860 1,658,156** MMAS→IVND 500 52 11036 12,845,438 12,777,577 117.7 **13,234,362 12,776,833** PGAC→IVND 500 52 3547 7,292,728 **7,246,237** 27.7 **7,625,325 7,246,237** PGAC→IVND

**Prev. best known**

A Comparative Study on Meta Heuristic Algorithms for Solving Multilevel Lot-Sizing Problems

**New best known**

**New source**

101

http://dx.doi.org/10.5772/55279

items when the lot-size of end-product is changed.

**Instance IVND1000**

**I P Inter D. Avg. of 10 Min. of 10**

Note. Boldface type denotes previous best solution; underline type denotes better solution; Boldface&underline denotes the new best solution.

**Table 6.** Results of 40 medium-sized problems and the new best solutions

#### **5.5. Large-sized MLLS problems**

Next, we ran IVND1000 <sup>50</sup> algorithm on 40 large-sized MLLS benchmark problems and repeated 10 runs. We summarize the 400 results and compare them with the existing algorithms in Table 7, and show detailed results of 40 problems in Table 8.


**Table 7.** Comparing IVND with existing algorithms on 40 large-sized problems

It can be observed from Table 7 and Table 8 that the IVND algorithm shows its best optimality among all existing algorithms since 70% of these 40 problems were found new best solution by IVND algorithm. The average computing time for each problem used by IVND algorithm was relatively low. However, four problems, i.e., problem 19, 15, 25, and 50, used much long time to finish their calculation because the interdependencies among items are relatively high for these four problems. The column *Inter D.* in Table 8 is the maximum number of affected items when the lot-size of end-product is changed.

**Instance IVND400 Prev. best**

Avg. computing time 2.67s 26.7s

100 Recent Advances on Meta-Heuristics and Their Application to Real Scenarios

**Table 6.** Results of 40 medium-sized problems and the new best solutions

7, and show detailed results of 40 problems in Table 8.

the new best solution.

Next, we ran IVND1000

IVND 50

IVND 50

**5.5. Large-sized MLLS problems**

**Method Avg. Cost**

**Solution**

<sup>50</sup> algorithm on 40 large-sized MLLS benchmark problems and repeated

**(m) CPU type**

**S D I P Avg. of 10 Min. of 10 method**

38 4 9 40 24 367,531 367,326 **366,872 366,848** PGCA→IVND 39 4 10 40 24 305,729 305,363 **305,172 305,053** PGCA→IVND

Note. Boldface type denotes previous best solution; underline type denotes better solution; Boldface&underline denotes

10 runs. We summarize the 400 results and compare them with the existing algorithms in Table

**Time**

HGA1000\* 40,817,600 10.00 -- 1 Dellaert et all(2000) MMAS 40,371,702 47.50 P4 1.5G 1 Pitakaso et al.(2007) GA100 41,483,590 0.00 60 P4 2.4G 1 Homberger(2008) GA3000 -- -- 60 P4 2.4G 1 Homberger(2008) PGAI 41,002,743 0.00 60 P4 2.4G 30 Homberger(2008) PGAC 39,809,739 52.50 60 P4 2.4G 30 Homberger(2008)

<sup>1000</sup>(Avg. of 10) 40,051,638 65.00 4.44 NB 1.6G <sup>1</sup> IVND

<sup>1000</sup>(Min. of 10) 39,869,210 70.00 44.4 NB 1.6G <sup>1</sup> IVND

It can be observed from Table 7 and Table 8 that the IVND algorithm shows its best optimality among all existing algorithms since 70% of these 40 problems were found new best solution by IVND algorithm. The average computing time for each problem used by IVND algorithm was relatively low. However, four problems, i.e., problem 19, 15, 25, and 50, used much long

**Prev. best solution** 39,792,241 -- -- -- -- **-- New best solution** 39,689,769 **-- -- -- - --**

**Table 7.** Comparing IVND with existing algorithms on 40 large-sized problems

**Optimality on prev. best solutions (%)**

Average 263,529 263,399 263,199.8 260,677.1

**New best solution**

**Number of processors**

**New**

**Sources**



**P N Kmax Avg. Cost of 10 runs Comp. time of 10 runs (s)**

A Comparative Study on Meta Heuristic Algorithms for Solving Multilevel Lot-Sizing Problems

http://dx.doi.org/10.5772/55279

103

5 200 5 264,111 30

10 200 5 263,914 57

20 200 5 263,816 104

30 200 5 263,744 144

50 200 5 263,712 233

70 200 5 263,671 308

100 200 5 263,640 433

130 200 5 263,614 553

160 200 5 263,620 674

200 200 5 263,579 832

50 10 5 266,552 40

50 50 5 264,395 87

50 100 5 263,920 135

50 150 5 263,775 186

50 200 5 263,702 220

50 250 5 263,672 274

50 300 5 263,634 309

50 400 5 263,613 376

50 500 5 263,603 463

50 600 5 263,585 559

50 200 1 263,868 74

50 200 2 263,775 100

50 200 3 263,715 136

50 200 4 263,713 178

50 200 5 263,709 225

50 200 6 263,704 276

50 200 7 263,693 341

Note. Boldface type denotes previous best solution; underline type denotes better solution; Boldface&underline denotes the new best solution.

**Table 8.** Results of 40 large-sized problems and the new best solutions

#### **5.6. The effectiveness of individual parameter of VIND**

Finally, we used the 40 medium-sized problems to test parameters, i.e., *P*, *N* and *K*max, of IVND algorithm their relation between effectiveness and computing load (using medium-sized problems is just for saving computing time). We did three experiments by varying one parameter while fixing other two parameters. First, we fixed *N*=200 and *K*max=5, and increased *P* from 10 to 200, and repeated 10 runs for each *P*. Second, we fixed *P*=50 and *K*max=5, and increase N from 50 to 500. thirdly, *P*, *N* were fixed to 50 and 200, and *K*max was increased from 1 to 10. The average costs gotten by the three experiments against varied parameters are shown in Table 9. A general trend can be observed that increases parameter *P*, *N* or *K*max will all lead to better solutions been found but at the price of more computing time. However, all the parameter may contribute less to the quality of solution when they are increased large enough. Obviously, the best effectiveness-cost combination of these parameters exists for the IVND algorithm which is a worthwhile work to do in future works, but we just set these parameters manually in our experiments.

A Comparative Study on Meta Heuristic Algorithms for Solving Multilevel Lot-Sizing Problems http://dx.doi.org/10.5772/55279 103


**Instance IVND1000**

102 Recent Advances on Meta-Heuristics and Their Application to Real Scenarios

**I P Inter D. Avg. of 10 Min. of 10**

Avg. computing time 4.44m 44.4m

**Table 8.** Results of 40 large-sized problems and the new best solutions

**5.6. The effectiveness of individual parameter of VIND**

the new best solution.

manually in our experiments.

**Prev. best known**

**Time (m)**

 500 52 1034 4,253,400 4,231,896 21.9 **4,320,868 4,199,064** PGAC→IVND 500 52 559 2,905,006 2,889,328 10.7 **2,996,500 2,864,526** MMAS→IVND 500 52 341 2,198,534 **2,186,429** 7.8 **2,277,630 2,186,429** MMAS→IVND 500 52 193607 103,535,103 103,237,091 297.4 **102,457,238 102,457,238** PGAC 500 52 22973 18,160,129 18,104,424 49.4 **18,519,760 18,097,215** PGAC→IVND 500 52 3247 6,932,353 **6,905,070** 44.3 **7,361,610 6,905,070** MMAS→IVND 500 52 914 4,109,712 4,095,109 41.1 **4,256,361 4,080,792** PGAC→IVND 500 52 708 2,602,841 2,573,491 20.5 **2,672,210 2,568,339** MMAS→IVND 500 52 1099608 768,067,039 762,331,081 121.4 **756,980,807 756,980,807** PGAC 500 52 24234 33,393,240 33,377,419 137.7 **33,524,300 33,356,750** MMAS→IVND 500 52 7312 10,506,439 **10,491,324** 52.1 **10,745,900 10,491,324** MMAS→IVND 500 52 1158 5,189,651 5,168,547 29.5 **5,198,011 5,120,701** PGAC→IVND 500 52 982 3,406,764 3,394,470 14.3 **3,485,360 3,381,090** MMAS→IVND Average 40,051,638 39,869,210 44.4 39,792,241 39,689,769 --

Note. Boldface type denotes previous best solution; underline type denotes better solution; Boldface&underline denotes

Finally, we used the 40 medium-sized problems to test parameters, i.e., *P*, *N* and *K*max, of IVND algorithm their relation between effectiveness and computing load (using medium-sized problems is just for saving computing time). We did three experiments by varying one parameter while fixing other two parameters. First, we fixed *N*=200 and *K*max=5, and increased *P* from 10 to 200, and repeated 10 runs for each *P*. Second, we fixed *P*=50 and *K*max=5, and increase N from 50 to 500. thirdly, *P*, *N* were fixed to 50 and 200, and *K*max was increased from 1 to 10. The average costs gotten by the three experiments against varied parameters are shown in Table 9. A general trend can be observed that increases parameter *P*, *N* or *K*max will all lead to better solutions been found but at the price of more computing time. However, all the parameter may contribute less to the quality of solution when they are increased large enough. Obviously, the best effectiveness-cost combination of these parameters exists for the IVND algorithm which is a worthwhile work to do in future works, but we just set these parameters

**New best known**

**New source**


**Acknowledgements**

grant No. 24510192.

**Author details**

**References**

*ence*, , 30, 223-239.

857-866.

Ikou Kaku1\*, Yiyong Xiao2

This work is supported by the Japan Society for the Promotion of Science (JSPS) under the

A Comparative Study on Meta Heuristic Algorithms for Solving Multilevel Lot-Sizing Problems

http://dx.doi.org/10.5772/55279

105

1 Department of Environmental and Information studies, Tokyo City University, Japan

3 College of Business Administration, Zhejiang University of Technology, Hangzhou, China

[1] Afentakis, P, Gavish, B, & Kamarkar, U. (1984). Computationally efficient optimal solutions to the lot-sizing problem in multistage assembly systems, *Management Sci‐*

[2] Afentakis, P, & Gavish, B. (1986). Optimal lot-sizing algorithms for complex product

[3] Almeder, C. (2010). A hybrid optimization approach for multi-level capacitated lot-

[4] Benton, W. C, & Srivastava, R. (1985). Product structure complexity and multilevel

[5] Blackburn, J. D, & Millen, R. A. (1985). An evaluation of heuristic performance in multi-stage lot-sizing systems, *International Journal of Production Research*, , 23,

[6] Coleman, B. J, & Mcknew, M. A. (1991). An improved heuristic for multilevel lot siz‐

[7] Crowston, W. B, & Wagner, H. M. (1973). Dynamic lot size models for multi-stage as‐

sizing problems, *European Journal of Operational Research*, , 200, 599-606.

lot sizing using alternative costing policies, *Decision Sciences*, , 16, 357-369.

ing in material requirements planning. *Decision Sciences*, , 22, 136-156.

2 School of Reliability and System Engineering, Beihang University, Beijing, China

and Yi Han3

\*Address all correspondence to: kakuikou@tcu.ac.jp

structures, *Operations Research*, , 34, 237-249.

sembly system, *Management Science*, , 20, 14-21.

**Table 9.** Experimental results of different parameters for medium-sized problem

#### **6. Summarization**

The consideration of meta-heuristic is widly used in a lot of fields. Deffirent meta-heuristic algorithms are developed for solving deffirent problems, especially combinational optimiza‐ tion problems. In this chapter, we discussed a special case of MLLS problem. First, the general definition of MLLS problem was described. We shown its solution structure and explained its NP completeness. Second, we reviewed the meta-heuristic algorithms which have been use to solve the MLLS problem and pointed their merits and demerits. Based on the recognition, third, we investigated those implement techniques used in the meta-heuristic algorithms for solving the MLLS problems. And two criteria of distance and range were firstly defined to evaluate the effective of those techniques. We brifly discussed the mechanisms and character‐ istics of the techniques by using these two criteria, and provided simulation experiments to prove the correctness of the two criteria and to explain the performance and utility of them. This is one of our contributions. Fourth, we presented a succinct and easily implemented IVND algorithm and six implemental techniques for solving the MLLS problem. The IVND algorithm was evaluated by using 176 benchmark problems of different scales (small, medium and large) from literatures. The results on 96 small-sized benchmark problems showed the IVND algorithm of good optimality; it could find 100% optimal solutions in repeated 10 runs using a very low computing time(less than 1s for each problem). Experiments on other two sets of benchmark problems (40 medium-sized problems and 40 large-sized problems) showed it good efficiency and effectiveness on solving MLLS problem with product structure complex. For the medium-sized problems, the IVND can use 10 repeated runs to reach 40% of the 40 problems of their previous known best solutions and find another 20% of new best known solutions. By more repeated runs, our IVND algorithm actually had updated 30% (12 prob‐ lems) of the best known solutions, and computing efficiency was also very acceptable because the longest computing time for each problem was less than one minute. For the 40 large-sized problems, the IVND algorithm delivered even more exciting results on the quality of solution. Comparison of the best solutions achieved with the new method and those established by previous methods including HGA, MMAS, and PGA shows that the IVND algorithm with the six implemental techniques are till now among the best methods for solving MLLS problem with product structure complexity considered, not only because it is easier to be understood and implemented in practice, but more importantly, it also provides quite good solutions in very acceptable time.

## **Acknowledgements**

**P N Kmax Avg. Cost of 10 runs Comp. time of 10 runs (s)**

50 200 8 263,684 420

50 200 9 263,688 527

50 200 10 263,678 567

The consideration of meta-heuristic is widly used in a lot of fields. Deffirent meta-heuristic algorithms are developed for solving deffirent problems, especially combinational optimiza‐ tion problems. In this chapter, we discussed a special case of MLLS problem. First, the general definition of MLLS problem was described. We shown its solution structure and explained its NP completeness. Second, we reviewed the meta-heuristic algorithms which have been use to solve the MLLS problem and pointed their merits and demerits. Based on the recognition, third, we investigated those implement techniques used in the meta-heuristic algorithms for solving the MLLS problems. And two criteria of distance and range were firstly defined to evaluate the effective of those techniques. We brifly discussed the mechanisms and character‐ istics of the techniques by using these two criteria, and provided simulation experiments to prove the correctness of the two criteria and to explain the performance and utility of them. This is one of our contributions. Fourth, we presented a succinct and easily implemented IVND algorithm and six implemental techniques for solving the MLLS problem. The IVND algorithm was evaluated by using 176 benchmark problems of different scales (small, medium and large) from literatures. The results on 96 small-sized benchmark problems showed the IVND algorithm of good optimality; it could find 100% optimal solutions in repeated 10 runs using a very low computing time(less than 1s for each problem). Experiments on other two sets of benchmark problems (40 medium-sized problems and 40 large-sized problems) showed it good efficiency and effectiveness on solving MLLS problem with product structure complex. For the medium-sized problems, the IVND can use 10 repeated runs to reach 40% of the 40 problems of their previous known best solutions and find another 20% of new best known solutions. By more repeated runs, our IVND algorithm actually had updated 30% (12 prob‐ lems) of the best known solutions, and computing efficiency was also very acceptable because the longest computing time for each problem was less than one minute. For the 40 large-sized problems, the IVND algorithm delivered even more exciting results on the quality of solution. Comparison of the best solutions achieved with the new method and those established by previous methods including HGA, MMAS, and PGA shows that the IVND algorithm with the six implemental techniques are till now among the best methods for solving MLLS problem with product structure complexity considered, not only because it is easier to be understood and implemented in practice, but more importantly, it also provides quite good solutions in

**Table 9.** Experimental results of different parameters for medium-sized problem

104 Recent Advances on Meta-Heuristics and Their Application to Real Scenarios

**6. Summarization**

very acceptable time.

This work is supported by the Japan Society for the Promotion of Science (JSPS) under the grant No. 24510192.

## **Author details**

Ikou Kaku1\*, Yiyong Xiao2 and Yi Han3

\*Address all correspondence to: kakuikou@tcu.ac.jp

1 Department of Environmental and Information studies, Tokyo City University, Japan

2 School of Reliability and System Engineering, Beihang University, Beijing, China

3 College of Business Administration, Zhejiang University of Technology, Hangzhou, China

## **References**


[8] Dellaert, N, & Jeunet, J. (2000). Solving large unconstrained multilevel lot-sizing problems using a hybrid genetic algorithm, *International Journal of Production Re‐ search*, 38(5), 1083-1099.

[22] Pitakaso, R, Almeder, C, Doerner, K. F, & Hartlb, R. F. ant system for unconstrained multi-level lot-sizing problems, *Computer & Operations Research*, , 34, 2533-2552. [23] Raza, A. S, & Akgunduz, A. (2008). A comparative study of heuristic algorithms on Economic Lot Scheduling Problem, *Computers and Industrial Engineering.* 55(1),

A Comparative Study on Meta Heuristic Algorithms for Solving Multilevel Lot-Sizing Problems

http://dx.doi.org/10.5772/55279

107

[24] Steinberg, E, & Napier, H. A. (1980). Optimal multilevel lot sizing for requirements

[25] Tang, O. (2004). Simulated annealing in lot sizing problems, *International Journal of*

[26] Veral, E. A. and LaForge R. L., (1985). The performance of a simple incremental lotsizing rule in a multilevel inventory environment, *Decision Sciences*, , 16, 57-72. [27] Xiao, Y. Y, Kaku, I, Zhao, X. H, & Zhang, R. Q. (2011a). A variable neighborhood search based approach for uncapacitated multilevel lot-sizing problems, *Computers &*

[28] Xiao, Y. Y, Zhao, X. H, Kaku, I, & Zhang, R. Q. (2011b). A reduced variable neighbor‐ hood search algorithm for uncapacitated multilevel lot-sizing problems, *European*

[29] Xiao, Y. Y, Zhao, X. H, Kaku, I, & Zhang, R. Q. (2012). Neighborhood search techni‐ ques for solving uncapacitated multilevel lot-sizing problems, *Computers & Opera‐*

[30] Yelle, L. E. (1979). Materials requirements lot sizing: a multilevel approach, *Interna‐*

[31] Zhangwill, W. I. (1968). Minimum concave cost flows in certain network, *Management*

[32] Zhangwill, W. I. (1969). A backlogging model and a multi-echelon model of a dy‐ namic economic lot size production system-a network approach, *Management Sci‐*

planning systems, *Management Science*, 26(12), 1258-1271.

*Production Economics*, , 88, 173-181.

*Industrial Engineering*, , 60, 218-227.

*tions Research*, 57((3)

*Scienc*e, , 14, 429-450.

*ence*, , 15, 506-527.

*Journal of Operational Research, 214*, 223-231.

*tional Journal of Production Research*, , 17, 223-232.

94-109.


[22] Pitakaso, R, Almeder, C, Doerner, K. F, & Hartlb, R. F. ant system for unconstrained multi-level lot-sizing problems, *Computer & Operations Research*, , 34, 2533-2552.

[8] Dellaert, N, & Jeunet, J. (2000). Solving large unconstrained multilevel lot-sizing problems using a hybrid genetic algorithm, *International Journal of Production Re‐*

[9] Dellaert, N, Jeunet, J, & Jonard, N. (2000). A genetic algorithm to solve the general multi-level lot-sizing problem with time-varying costs, *International Journal of Produc‐*

[10] Han, Y, Tang, J. F, Kaku, I, & Mu, L. F. (2009). Solving incapacitated multilevel lotsizing problem using a particle swarm optimization with flexible inertial weight,

[11] Han, Y, Kaku, I, Tang, J. F, Dellaert, N, Cai, J. H, Li, Y. L, & Zhou, G. G. (2011). A scatter search approach for uncapacitated multilevel lot-sizing problems, *Internation‐*

[12] Han, Y, Cai, J. H, Kaku, I, Lin, H. Z, & Guo, H. D. (2012a). A note on "a genetic algo‐ rithm for the preemptive and non-preemptive multi-mode resource-constrained

[13] Han, Y, Cai, J. H, Kaku, I, Li, Y. L, Chen, Y. Z, & Tang, J. F. (2012b). Evolutionary algorithms for solving unconstrained multilevel lot-sizing problem with series struc‐

[14] Hansen, P, & Mladenovic, N. (2001a). Variable neighborhood search: principles and

[15] Hansen, P, Mladenovic, N, & Perez, D. (2001b). Variable neighborhood decomposi‐

[16] Hansen, P, & Mladenovic, N. and Pe´rez J. A. M., (2008). Variable neighborhood

[17] Homberger, J. (2008). A parallel genetic algorithm for the multilevel unconstrained

[18] Hoos, H. H, & Thomas, S. (2005). Stochastic Local Search-Foundations and Applica‐

[19] Kaku, I, & Xu, C. H. (2006). A soft optimization approach for solving a complicated multilevel lot-sizing problem, *in Proc. 8th Conf. Industrial Management, ICIM'200638*

[20] Kaku, I, Li, Z. S, & Xu, C. H. (2010). Solving Large Multilevel Lot-Sizing Problems with a Simple Heuristic Algorithm Based on Segmentation, *International Journal of In‐*

[21] Mladenovic, N, & Hansen, P. (1997). Variable neighborhood search, *Computers & Op‐*

applications, *European Journal of Operational Research*, 130(3), 449-467.

search, *European Journal of Operational Research*, Editorial., 191, 593-595.

lot-sizing problem, *INFORMS Journal on Computing*, 20(1), 124-132

*novative Computing, Information and Control*, 6(3), 817-827.

*al Journal of Innovative Computing, Information and Control*, 7(7), 4833-4848.

project scheduling problem", *Applied Mechanics and Materials*, , 127, 527-530.

*Computers and Mathematics with Applications*, , 57, 1748-1755.

ture, *Journal of Shanghai Jiaotong University*, 17(1), 39-44.

tion search, *Journal of Heuristics*, , 7, 335-350.

tions, Morgan Kaufmann Publishers.

*erations Research*, , 24, 1097-1100.

*search*, 38(5), 1083-1099.

*tion Economics*, , 68, 241-257.

106 Recent Advances on Meta-Heuristics and Their Application to Real Scenarios


**Chapter 5**

**Provisional chapter**

**A Two-Step Optimisation Method for Dynamic**

**A Two-Step Optimisation Method for Dynamic**

The weapon target assignment (WTA) problem has been designed to match the Command & Control (C2) requirement in military context, of which the goal is to find an allocation plan enabling to treat a specific scenario in assigning available weapons to oncoming targets. The WTA always get into situation weapons defending an area or assets from an enemy aiming to destroy it. Because of the uniqueness of each situation, this problem must be solved in real-time and evolve accordingly to the aerial/ground situation. By the past, the WTA was solved by an operator taking all the decisions, but because of the complexity of the modern warfare, the resolution of the WTA in using the power of computation is inevitable to make possible the resolution in real time of very complex scenarii involving different type of targets. Nowadays, in most of the C2 this process is designed in order to be as a support for a human operator and in helping him in the decision making process. The operator will

The WTA arouses a great interest among the researcher community and many methods have been proposed to cope with this problem. Besides, the WTA has been proved to be NP-complete [1]. There are two families of WTA: the Static WTA (SWTA) and the Dynamic WTA (DWTA). In both of these problems, the optimality of one solution is based either on the minimisation of the target survival after the engagement or the maximisation of the survivability of the defended assets. The main feature of the SWTA stands in its single stage approach. It is considered that all the information about the situations are provided and the problem can be considered as a constrained resource assignment problem. In contrast, the DWTA is a multi-stage problem in which the result of each stage is assessed, then use to update the aerial situation for the upcoming stages. The DWTA can also be expressed as

> ©2012 Siarry et al., licensee InTech. This is an open access chapter distributed under the terms of the Creative Commons Attribution License (http://creativecommons.org/licenses/by/3.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited. © 2013 Leboucher et al.; licensee InTech. This is an open access article distributed under the terms of the Creative Commons Attribution License (http://creativecommons.org/licenses/by/3.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

© 2013 Leboucher et al., licensee InTech. This is a paper distributed under the terms of the Creative Commons Attribution License (http://creativecommons.org/licenses/by/3.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

**Weapon Target Assignment Problem**

**Weapon Target Assignment Problem**

Cédric Leboucher, Hyo-Sang Shin, Patrick Siarry, Rachid Chelouah,

Cédric Leboucher, Hyo-Sang Shin, Patrick Siarry,

Rachid Chelouah, Stéphane Le Ménec and

Additional information is available at the end of the chapter

Stéphane Le Ménec and Antonios Tsourdos

Additional information is available at the end of the chapter

give its final green light to proceed the intervention.

Antonios Tsourdos

10.5772/53606

**1. Introduction**

http://dx.doi.org/10.5772/53606

**Provisional chapter**

## **A Two-Step Optimisation Method for Dynamic Weapon Target Assignment Problem A Two-Step Optimisation Method for Dynamic Weapon Target Assignment Problem**

Cédric Leboucher, Hyo-Sang Shin, Patrick Siarry, Rachid Chelouah, Stéphane Le Ménec and Antonios Tsourdos Cédric Leboucher, Hyo-Sang Shin, Patrick Siarry, Rachid Chelouah, Stéphane Le Ménec and Antonios Tsourdos Additional information is available at the end of the chapter

Additional information is available at the end of the chapter 10.5772/53606

http://dx.doi.org/10.5772/53606

## **1. Introduction**

The weapon target assignment (WTA) problem has been designed to match the Command & Control (C2) requirement in military context, of which the goal is to find an allocation plan enabling to treat a specific scenario in assigning available weapons to oncoming targets. The WTA always get into situation weapons defending an area or assets from an enemy aiming to destroy it. Because of the uniqueness of each situation, this problem must be solved in real-time and evolve accordingly to the aerial/ground situation. By the past, the WTA was solved by an operator taking all the decisions, but because of the complexity of the modern warfare, the resolution of the WTA in using the power of computation is inevitable to make possible the resolution in real time of very complex scenarii involving different type of targets. Nowadays, in most of the C2 this process is designed in order to be as a support for a human operator and in helping him in the decision making process. The operator will give its final green light to proceed the intervention.

The WTA arouses a great interest among the researcher community and many methods have been proposed to cope with this problem. Besides, the WTA has been proved to be NP-complete [1]. There are two families of WTA: the Static WTA (SWTA) and the Dynamic WTA (DWTA). In both of these problems, the optimality of one solution is based either on the minimisation of the target survival after the engagement or the maximisation of the survivability of the defended assets. The main feature of the SWTA stands in its single stage approach. It is considered that all the information about the situations are provided and the problem can be considered as a constrained resource assignment problem. In contrast, the DWTA is a multi-stage problem in which the result of each stage is assessed, then use to update the aerial situation for the upcoming stages. The DWTA can also be expressed as

©2012 Siarry et al., licensee InTech. This is an open access chapter distributed under the terms of the Creative Commons Attribution License (http://creativecommons.org/licenses/by/3.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited. © 2013 Leboucher et al.; licensee InTech. This is an open access article distributed under the terms of the Creative Commons Attribution License (http://creativecommons.org/licenses/by/3.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited. © 2013 Leboucher et al., licensee InTech. This is a paper distributed under the terms of the Creative Commons Attribution License (http://creativecommons.org/licenses/by/3.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

a succession of SWTA, but the optimality of the final solution cannot be guaranteed since it comes to the same as in a greedy optimisation process. One other difference stands in the temporal dimension of the DWTA which does not exist in the SWTA. The weapons can intervene within a certain defined time because of physical, technical and operational constraints. In addition, any DWTA problem has to be solved in using real-time oriented method. By real-time it is assumed that the proposed method has to be fast enough to provide an engagement solution before the oncoming targets reached their goals. Most of the previous work on the WTA was focused on the resolution of the SWTA. Hosein and Athans was among the first to defined a cost function based on the assets [2]. This model was reused in [3] and [4]. Later, a second modelling has been proposed by Karasakal in [5], aiming to maximise the probability to suppress all the oncoming targets. One other variant of the WTA is to take into account a *threatening* value to each target according to its features and the importance of the protected assets. The research of Johansson and Falkman in [6] proposed a good overview of all the possible modelling, taking into account both of the developed models and enabling to take into consideration the value of the defended assets and the threatening index of the incoming target. Kwon *et al.* explored further this principle in assigning a value to the weapon in [7]. The main researches on the SWTA started around the 1950's. Most of the proposals to solve this problem was based on the classic optimisation processes: branch and bound algorithm appears in the survey conducted in 2006 by Cai *et al.* [8]. With the evolution of the new technologies, some more complex methods appeared in [9] in using the neural networks. The genetic algorithms are used in [4], [10] and [11] to solve the SWTA. Cullenbine is using the Tabu Search method in [12]. A different approach angle is used in [13]. In this former approach the WTA problem is treated as a resources management problem and the reactivity of the proposed approach, based on the Tabu Search, was able to deal with real-time requirement. Nowadays, this method is used in many military systems like Rapid Anti-Ship Missile-Integrated Defense System (RAIDS) [14] [13]. Whereas the SWTA had aroused the interest of the researchers first, lately the DWTA had attracted much more attention. The first DWTA was proposed by Hosein and Athans around 1990 [15]. In the proposed approach of Hosein and Athans, a sub-optimal solution was studied in order to determine a solution which was considered as "good enough" [15]. Later they developed exact methods to solve some simplified DWTA [16] [17]. The dynamic programming enables to solve the DWTA in [18], but under the assumption that all the engaged targets are destroyed. Despite its study to decrease the computational time, the problem was still treated in exponential complexity [18]. A more complex DWTA model is designed by Wu *et al.* in [19] where the temporal dimension is included under the form of firing time windows.

10.5772/53606

111

http://dx.doi.org/10.5772/53606

weapons to the targets, then the optimal firing sequence is obtained in using the results obtained from the first step. The optimal assignment is determined in using the graph theory, and more especially the Hungarian method in a bipartite graph. The used of this method in the first step is motivated by the optimality and the polynomial complexity of the method. Then, the computation of the firing sequence is optimised in using a particle swarm optimisation (PSO) process combined to the evolutionary game theory (EGT). This former

A Two-Step Optimisation Method for Dynamic Weapon Target Assignment Problem

The performance index for the evaluation of the assignment is determined by three different criteria: the capacity to propose an early fire, the width of the firing time window and the minimisation of the overflying of the defended area by our own assets for security purpose. The quality of the firing sequence is obtained from the reactivity of the algorithm to treat the targets in the earliest possible way, the respect of the system constraints and the avoidance

The goal of this chapter is to develop an efficient method to solve a target based DWTA problem involving technical and operational constraints. A mission is considered as achieve only if no targets reach the defending area. The contribution of this paper includes the

• Design of a DWTA model taking into account target trajectories and operational and

• A two-step approach based on the graph theory, then a combined swarm intelligence and

• The targets are following a Bezier's curve trajectory in order to sow the confusion among

• The success of one fire is determined by the draw of one random number in [0, 1], then

The rest of this chapter is organised as follows: the second section describes the details of the studied DWTA, the third section introduces the background of Hungarian algorithm, particle swarm optimisation and evolutionary game theory. Then the fourth section details the proposed method before testing and analysing the obtained results by using a dedicated simulator designed for this DWTA problem. The chapter ends with the conclusion of this

The assignment problem arouses the interest of the researchers community for a while. The principle consists of finding a maximum weight matching in a weighted bipartite graph. It is more commonly formulated as: there are two distinct sets, one contain agents, the other one contain tasks. Note that each agent has his own ability to realise a job properly and this capability is represented by a quantitative value. The global objective to assign all the agents to the jobs can be achieved in one optimal way. The Hungarian method published by

evolutionary game method to solve the DWTA in an optimised fashion. • The reducing computational load in order to enable real-time applications.

method has been proved as efficient in general allocation resources problem [20].

of idle time when a firing is possible.

technical constraints on the weapons.

compared to a probability threshold of kill (PK).

**2. Background of the proposed approach**

**2.1. The Hungarian algorithm**

following aspects:

study.

the defending system.

The studied DWTA in this chapter slightly differs from the common defined DWTA in the literature. The proposed model has been designed to fit a specific requirement from industrial application. Whereas the classic problem is considering a multi-stage approach, the solved problem considers a continuous time where the targets are evolving in the space according to their own objectives and features. The targets trajectories are designed in using Bezier's curves defined by 4 control points which the last one is set to the centre of area that we are defending. The choice of this trajectory modelling has been done in order to add more diversity in the tested scenarii. The current situation is updated in real time, which means that the proposed algorithm must be as reactive as possible to cope with the oncoming targets. In order to solve the presented problem in the fastest and the most accurate way, a two-step optimisation method is proposed. The first step optimise the assignment of the weapons to the targets, then the optimal firing sequence is obtained in using the results obtained from the first step. The optimal assignment is determined in using the graph theory, and more especially the Hungarian method in a bipartite graph. The used of this method in the first step is motivated by the optimality and the polynomial complexity of the method. Then, the computation of the firing sequence is optimised in using a particle swarm optimisation (PSO) process combined to the evolutionary game theory (EGT). This former method has been proved as efficient in general allocation resources problem [20].

The performance index for the evaluation of the assignment is determined by three different criteria: the capacity to propose an early fire, the width of the firing time window and the minimisation of the overflying of the defended area by our own assets for security purpose. The quality of the firing sequence is obtained from the reactivity of the algorithm to treat the targets in the earliest possible way, the respect of the system constraints and the avoidance of idle time when a firing is possible.

The goal of this chapter is to develop an efficient method to solve a target based DWTA problem involving technical and operational constraints. A mission is considered as achieve only if no targets reach the defending area. The contribution of this paper includes the following aspects:


The rest of this chapter is organised as follows: the second section describes the details of the studied DWTA, the third section introduces the background of Hungarian algorithm, particle swarm optimisation and evolutionary game theory. Then the fourth section details the proposed method before testing and analysing the obtained results by using a dedicated simulator designed for this DWTA problem. The chapter ends with the conclusion of this study.

## **2. Background of the proposed approach**

## **2.1. The Hungarian algorithm**

2 Recent Advances on Meta-Heuristics and Their Application to Real Scenarios

included under the form of firing time windows.

a succession of SWTA, but the optimality of the final solution cannot be guaranteed since it comes to the same as in a greedy optimisation process. One other difference stands in the temporal dimension of the DWTA which does not exist in the SWTA. The weapons can intervene within a certain defined time because of physical, technical and operational constraints. In addition, any DWTA problem has to be solved in using real-time oriented method. By real-time it is assumed that the proposed method has to be fast enough to provide an engagement solution before the oncoming targets reached their goals. Most of the previous work on the WTA was focused on the resolution of the SWTA. Hosein and Athans was among the first to defined a cost function based on the assets [2]. This model was reused in [3] and [4]. Later, a second modelling has been proposed by Karasakal in [5], aiming to maximise the probability to suppress all the oncoming targets. One other variant of the WTA is to take into account a *threatening* value to each target according to its features and the importance of the protected assets. The research of Johansson and Falkman in [6] proposed a good overview of all the possible modelling, taking into account both of the developed models and enabling to take into consideration the value of the defended assets and the threatening index of the incoming target. Kwon *et al.* explored further this principle in assigning a value to the weapon in [7]. The main researches on the SWTA started around the 1950's. Most of the proposals to solve this problem was based on the classic optimisation processes: branch and bound algorithm appears in the survey conducted in 2006 by Cai *et al.* [8]. With the evolution of the new technologies, some more complex methods appeared in [9] in using the neural networks. The genetic algorithms are used in [4], [10] and [11] to solve the SWTA. Cullenbine is using the Tabu Search method in [12]. A different approach angle is used in [13]. In this former approach the WTA problem is treated as a resources management problem and the reactivity of the proposed approach, based on the Tabu Search, was able to deal with real-time requirement. Nowadays, this method is used in many military systems like Rapid Anti-Ship Missile-Integrated Defense System (RAIDS) [14] [13]. Whereas the SWTA had aroused the interest of the researchers first, lately the DWTA had attracted much more attention. The first DWTA was proposed by Hosein and Athans around 1990 [15]. In the proposed approach of Hosein and Athans, a sub-optimal solution was studied in order to determine a solution which was considered as "good enough" [15]. Later they developed exact methods to solve some simplified DWTA [16] [17]. The dynamic programming enables to solve the DWTA in [18], but under the assumption that all the engaged targets are destroyed. Despite its study to decrease the computational time, the problem was still treated in exponential complexity [18]. A more complex DWTA model is designed by Wu *et al.* in [19] where the temporal dimension is

The studied DWTA in this chapter slightly differs from the common defined DWTA in the literature. The proposed model has been designed to fit a specific requirement from industrial application. Whereas the classic problem is considering a multi-stage approach, the solved problem considers a continuous time where the targets are evolving in the space according to their own objectives and features. The targets trajectories are designed in using Bezier's curves defined by 4 control points which the last one is set to the centre of area that we are defending. The choice of this trajectory modelling has been done in order to add more diversity in the tested scenarii. The current situation is updated in real time, which means that the proposed algorithm must be as reactive as possible to cope with the oncoming targets. In order to solve the presented problem in the fastest and the most accurate way, a two-step optimisation method is proposed. The first step optimise the assignment of the

The assignment problem arouses the interest of the researchers community for a while. The principle consists of finding a maximum weight matching in a weighted bipartite graph. It is more commonly formulated as: there are two distinct sets, one contain agents, the other one contain tasks. Note that each agent has his own ability to realise a job properly and this capability is represented by a quantitative value. The global objective to assign all the agents to the jobs can be achieved in one optimal way. The Hungarian method published by Kuhn in 1955 is inspired by the work of two Hungarian researchers: Dénes Kõnig and Jenõ Egervàry [21]. This method has been proved as optimal and polynomial.

10.5772/53606

113

http://dx.doi.org/10.5772/53606

**2.3. The evolutionary game theory**

in the presence of the others.

*2.3.1. Evolutionary Stable Strategies*

distribution in the population is:

all *p* �= *p*∗,

strategy *q*.

*2.3.2. Replicator dynamics*

described using notations as follows. *ei* is the *i*

(probability distributions over the pure strategies *ei*).

*Aij* = *π*(*ei*,*ej*) is the *m* × *m* payoff matrix.

an individual using strategy *q* in this mixed population, be

for all *ε* > 0 "sufficiently small", for agents minimizing their fitness.

The evolutionary game theory appeared initially in a biologic context. The need to model the evolution phenomena led to the use of mathematical theory of the games to explain the strategic aspect of the evolution. Over the last few decades, the EGT has aroused interest of the economists, sociologists, social scientists, as well as the philosophers. Although the evolutionary game theory found its origin in biologic science, such an expansion to different fields can be explained by three facts. First of all, the notion of evolution has to be understood as the change of beliefs and norms over time. Secondly, the modelling of strategies change provides a social aspect which matches exactly the social system interactions. Finally, it was important to model dynamically the interactions within a population, which was one of the missing elements of the classic game theory. As in this former domain, the evolutionary game theory deals with the equilibrium which is a key point in both of the theories. Here the equilibrium point is called the evolutionary stable strategy. The principle of the EGT is not only based on the strategy performance obtained by itself, but also the performance obtained

A Two-Step Optimisation Method for Dynamic Weapon Target Assignment Problem

An Evolutionary Stable Strategy (ESS) is a strategy such that, if all members of a population adopt it, then no mutant strategy could invade the population under the influence of natural selection. Assume we have a mixed population consisting of mostly *p*<sup>∗</sup> individuals (agents playing optimal strategy *p*∗) with a few individuals using strategy *p*. That is, the strategy

(1 − *ε*)*p*<sup>∗</sup> + *εp* where *ε* > 0 is the small frequency of *p* users in the population. Let the fitness, i.e. payoff of

*π*(*q*,(1 − *ε*)*p*<sup>∗</sup> + *εp*).

Then, an interpretation of Maynard Smith's requirement [25] for *p*<sup>∗</sup> to be an ESS is that, for

*π*(*p*,(1 − *ε*)*p*<sup>∗</sup> + *εp*) > *π*(*p*∗,(1 − *ε*)*p*<sup>∗</sup> + *εp*)

A common way to describe strategy interactions is using matrix games. Matrix games are

∆*<sup>m</sup>* ≡ {*p* = (*p*1, ..., *pm*) | *p*<sup>1</sup> + ... + *pm* = 1, 0 ≤ *pi* ≤ 1} is the set of mixed strategies

Then, *π*(*p*, *q*) = *p* · *Aq<sup>T</sup>* is the payoff of agents playing strategy *p* facing agents playing

*th* unit line vector for *i* = 1, ..., *m*.

Let *G* be a complete bipartite graph composed, one hand by a set of |*A*| agents and one other hand by |*T*| tasks. Then *G* = (*A*, *T*, *E*), where *E* denotes the set of the edges linking the set of Agents with the set of Tasks. Note that each edge from *E* is weighted by a positive cost *c*(*i*, *j*), where *i* ∈ {1, . . . , |*A*|} and *j* ∈ {1, . . . , |*T*|}. The function *P* : (*A* ∪ *T*) −→ **R** represents the potential if *p*(*i*) + *p*(*j*) ≤ *c*(*i*, *j*) for each *i* ∈ *A* and *j* ∈ *T*. The potential value is obtained in summing all the potential from the set *A* ∪ *T*: *p* = ∑*v*∈(*A*∪*T*) *p*(*v*). The Hungarian method enables to find the perfect matching and the potential equalising the cost and the value, which means that both of them are optimal.

## **2.2. The particle swarm optimisation**

Kennedy and Eberhart [22], the founders of the PSO method, was inspired by the behaviour of animals acting in society to achieve a goal. For example, the birds, the fishes, etc. can make up a very efficient collective intelligence in exchanging very basic information about the environment in which they are evolving. From this starting point, the authors have designed the PSO method to solve many optimisation problems over the last few decades. A swarm is composed of particles (representing a solution) flying on the solution space and communicating with the neighbourhood the quality of the current position.

The first step in PSO algorithm is to define the moving rules on the solution space for the particles. Let *X<sup>t</sup> <sup>i</sup>* = [*x<sup>t</sup> <sup>i</sup>*1, *<sup>x</sup><sup>t</sup> <sup>i</sup>*2,..., *<sup>x</sup><sup>t</sup> iD*, ], *<sup>x</sup><sup>t</sup> id* ∈ {0, 1} be a particle in a population of *P* particles and composed of *D* dimensions. The velocity of this particle is denoted as *Vt <sup>i</sup>* = [*v<sup>t</sup> <sup>i</sup>*1, *<sup>v</sup><sup>t</sup> <sup>i</sup>*2,..., *<sup>v</sup><sup>t</sup> iD*, ], *<sup>v</sup><sup>t</sup> id* <sup>∈</sup> **<sup>R</sup>**. Then, as in the PSO method described in [23], the next step is to define the best position for the particle *P<sup>t</sup> <sup>i</sup>* = [*p<sup>t</sup> <sup>i</sup>*1, *<sup>p</sup><sup>t</sup> <sup>i</sup>*2,..., *<sup>p</sup><sup>t</sup> iD*, ], *<sup>p</sup><sup>t</sup> id* <sup>∈</sup> **<sup>R</sup>** , and the best position *P<sup>t</sup> <sup>g</sup>* = [*p<sup>t</sup> <sup>g</sup>*1, *<sup>p</sup><sup>t</sup> <sup>g</sup>*2,..., *<sup>p</sup><sup>t</sup> gD*, ], *<sup>p</sup><sup>t</sup> gd* <sup>∈</sup> **<sup>R</sup>** of the entire population at the iteration *<sup>t</sup>*. The velocity of the particle *i* is adjusted in respect to the direction *d* with:

$$
\sigma\_{id}^{t+1} = \omega\_1 \upsilon\_{id}^t + \omega\_2(\mathbf{x}(t) - \mathbf{x}\_{ind}(t)) + \omega\_3(\mathbf{x}(t) - \mathbf{x}\_{global}(t)).
$$

The parameter *ω*<sup>1</sup> denotes the weight of the particle inertia. *ω*<sup>2</sup> is the coefficient associated to the individual coefficient. Then, *ω*<sup>3</sup> denotes the social coefficient. The final step of one PSO iteration is to update the position of the particles in using the following formula:

$$\mathbf{X}\_{i}^{t+1} = \mathbf{X}\_{i}^{t} + V\_{i}^{t}.$$

This process enables to find an optimal solution in repeating this process. In the classical version of the PSO [24], these coefficients are drawn randomly in order to maximise the exploration of the solution space by the particles. It can be a weakness when the computational time has to be the shortest possible. The studied method proposes to decrease this computational time in using the Evolutionary Game Theory (EGT) to determine the three coefficients *ω*1, *ω*<sup>2</sup> and *ω*3. Since the particles are "jumping" on the solution space, the creators wished to limit the jumped distance to a maximum length determined by the value of *Vmax* usually determined with respect to the solution space.

## **2.3. The evolutionary game theory**

4 Recent Advances on Meta-Heuristics and Their Application to Real Scenarios

which means that both of them are optimal.

**2.2. The particle swarm optimisation**

*<sup>i</sup>* = [*x<sup>t</sup>*

step is to define the best position for the particle *P<sup>t</sup>*

*iD*, ], *<sup>v</sup><sup>t</sup>*

*<sup>g</sup>*1, *<sup>p</sup><sup>t</sup>*

*<sup>i</sup>*1, *<sup>x</sup><sup>t</sup>*

*<sup>g</sup>*2,..., *<sup>p</sup><sup>t</sup>*

of *Vmax* usually determined with respect to the solution space.

the particles. Let *X<sup>t</sup>*

*<sup>i</sup>*1, *<sup>v</sup><sup>t</sup>*

best position *P<sup>t</sup>*

*<sup>i</sup>*2,..., *<sup>v</sup><sup>t</sup>*

*<sup>g</sup>* = [*p<sup>t</sup>*

*vt*+<sup>1</sup> *id* <sup>=</sup> *<sup>ω</sup>*1*v<sup>t</sup>*

*Vt <sup>i</sup>* = [*v<sup>t</sup>*

Egervàry [21]. This method has been proved as optimal and polynomial.

Kuhn in 1955 is inspired by the work of two Hungarian researchers: Dénes Kõnig and Jenõ

Let *G* be a complete bipartite graph composed, one hand by a set of |*A*| agents and one other hand by |*T*| tasks. Then *G* = (*A*, *T*, *E*), where *E* denotes the set of the edges linking the set of Agents with the set of Tasks. Note that each edge from *E* is weighted by a positive cost *c*(*i*, *j*), where *i* ∈ {1, . . . , |*A*|} and *j* ∈ {1, . . . , |*T*|}. The function *P* : (*A* ∪ *T*) −→ **R** represents the potential if *p*(*i*) + *p*(*j*) ≤ *c*(*i*, *j*) for each *i* ∈ *A* and *j* ∈ *T*. The potential value is obtained in summing all the potential from the set *A* ∪ *T*: *p* = ∑*v*∈(*A*∪*T*) *p*(*v*). The Hungarian method enables to find the perfect matching and the potential equalising the cost and the value,

Kennedy and Eberhart [22], the founders of the PSO method, was inspired by the behaviour of animals acting in society to achieve a goal. For example, the birds, the fishes, etc. can make up a very efficient collective intelligence in exchanging very basic information about the environment in which they are evolving. From this starting point, the authors have designed the PSO method to solve many optimisation problems over the last few decades. A swarm is composed of particles (representing a solution) flying on the solution space and

The first step in PSO algorithm is to define the moving rules on the solution space for

*P* particles and composed of *D* dimensions. The velocity of this particle is denoted as

The parameter *ω*<sup>1</sup> denotes the weight of the particle inertia. *ω*<sup>2</sup> is the coefficient associated to the individual coefficient. Then, *ω*<sup>3</sup> denotes the social coefficient. The final step of one PSO iteration is to update the position of the particles in using the following formula:

This process enables to find an optimal solution in repeating this process. In the classical version of the PSO [24], these coefficients are drawn randomly in order to maximise the exploration of the solution space by the particles. It can be a weakness when the computational time has to be the shortest possible. The studied method proposes to decrease this computational time in using the Evolutionary Game Theory (EGT) to determine the three coefficients *ω*1, *ω*<sup>2</sup> and *ω*3. Since the particles are "jumping" on the solution space, the creators wished to limit the jumped distance to a maximum length determined by the value

*Xt*+<sup>1</sup> *<sup>i</sup>* <sup>=</sup> *<sup>X</sup><sup>t</sup>*

*id* ∈ {0, 1} be a particle in a population of

*iD*, ], *<sup>p</sup><sup>t</sup>*

*id* <sup>∈</sup> **<sup>R</sup>** , and the

*<sup>i</sup>*2,..., *<sup>p</sup><sup>t</sup>*

*gd* <sup>∈</sup> **<sup>R</sup>** of the entire population at the iteration *<sup>t</sup>*. The

*id* <sup>∈</sup> **<sup>R</sup>**. Then, as in the PSO method described in [23], the next

*<sup>i</sup>*1, *<sup>p</sup><sup>t</sup>*

*<sup>i</sup>* = [*p<sup>t</sup>*

*id* <sup>+</sup> *<sup>ω</sup>*2(*x*(*t*) <sup>−</sup> *xind*(*t*)) + *<sup>ω</sup>*3(*x*(*t*) <sup>−</sup> *xglobal*(*t*)).

*<sup>i</sup>* <sup>+</sup> *<sup>V</sup><sup>t</sup> i* .

*iD*, ], *<sup>x</sup><sup>t</sup>*

communicating with the neighbourhood the quality of the current position.

*<sup>i</sup>*2,..., *<sup>x</sup><sup>t</sup>*

*gD*, ], *<sup>p</sup><sup>t</sup>*

velocity of the particle *i* is adjusted in respect to the direction *d* with:

The evolutionary game theory appeared initially in a biologic context. The need to model the evolution phenomena led to the use of mathematical theory of the games to explain the strategic aspect of the evolution. Over the last few decades, the EGT has aroused interest of the economists, sociologists, social scientists, as well as the philosophers. Although the evolutionary game theory found its origin in biologic science, such an expansion to different fields can be explained by three facts. First of all, the notion of evolution has to be understood as the change of beliefs and norms over time. Secondly, the modelling of strategies change provides a social aspect which matches exactly the social system interactions. Finally, it was important to model dynamically the interactions within a population, which was one of the missing elements of the classic game theory. As in this former domain, the evolutionary game theory deals with the equilibrium which is a key point in both of the theories. Here the equilibrium point is called the evolutionary stable strategy. The principle of the EGT is not only based on the strategy performance obtained by itself, but also the performance obtained in the presence of the others.

#### *2.3.1. Evolutionary Stable Strategies*

An Evolutionary Stable Strategy (ESS) is a strategy such that, if all members of a population adopt it, then no mutant strategy could invade the population under the influence of natural selection. Assume we have a mixed population consisting of mostly *p*<sup>∗</sup> individuals (agents playing optimal strategy *p*∗) with a few individuals using strategy *p*. That is, the strategy distribution in the population is:

$$(1 - \varepsilon)p^\* + \varepsilon p$$

where *ε* > 0 is the small frequency of *p* users in the population. Let the fitness, i.e. payoff of an individual using strategy *q* in this mixed population, be

$$
\pi(q,(1-\varepsilon)p^\* + \varepsilon p).
$$

Then, an interpretation of Maynard Smith's requirement [25] for *p*<sup>∗</sup> to be an ESS is that, for all *p* �= *p*∗,

$$
\pi(p,(1-\varepsilon)p^\* + \varepsilon p) > \pi(p^\*,(1-\varepsilon)p^\* + \varepsilon p),
$$

for all *ε* > 0 "sufficiently small", for agents minimizing their fitness.

#### *2.3.2. Replicator dynamics*

A common way to describe strategy interactions is using matrix games. Matrix games are described using notations as follows. *ei* is the *i th* unit line vector for *i* = 1, ..., *m*.

*Aij* = *π*(*ei*,*ej*) is the *m* × *m* payoff matrix.

∆*<sup>m</sup>* ≡ {*p* = (*p*1, ..., *pm*) | *p*<sup>1</sup> + ... + *pm* = 1, 0 ≤ *pi* ≤ 1} is the set of mixed strategies (probability distributions over the pure strategies *ei*).

Then, *π*(*p*, *q*) = *p* · *Aq<sup>T</sup>* is the payoff of agents playing strategy *p* facing agents playing strategy *q*.

Another interpretation is *π*(*p*, *q*) being the fitness of a large population of agents playing pure strategies (*p* describing the agent proportion in each behaviour inside a population) with respect to a large *q* population.

10.5772/53606

115

http://dx.doi.org/10.5772/53606

Where t denotes the simulation time and the *Wi*, *i* ∈ {1, 2, 3} and *Tj*, *j* ∈ {1, 2} represent the weapon *i* and the target *j*. The variable *FTi* denotes the Firing Time computed for the weapon *i*. The engagement plan evolves accordingly to the situation and depends on the current simulation time and on the aerial situation. In this application, the engagement plan is recomputed every *P* seconds in order to make up a very reactive engagement plan capable of dealing with the trickier cases in which the targets are constantly changing their

A Two-Step Optimisation Method for Dynamic Weapon Target Assignment Problem

Since the complexity of the presented problem grows exponentially with the number of targets and weapons, to design an algorithm capable of handling the real-time computation, but taking into account very diversified performance indexes, the choice of two different steps was natural. Lloyd [1] proved that the DWTA is a NP-complete problem. Therefore, it is hard to find an exact optimisation method capable of solving the DWTA problem in an exact way within a reasonable time. The reasonable time implies a high frequency which can enable the real-time application of the optimisation method. Note that the system must be able to provide results in real-time in the DWTA problem since the engagement changes as the targets keep evolving in the aerial space during the computation. With these reasons, using a heuristic approach providing suboptimal solutions in real-time could be the best way to handle the DWTA problem. One other problem is to be able to quantify the quality of one proposed solution: the performance index of the assignment, and the firing sequence cannot be evaluated in using the same performance criterion. Whereas the assignment is evaluated from the system point of view, the firing sequence is evaluated from the weapons features. Dividing the problem into two parts could lead to the modification of the solution space and the optimum solution could be not the same as the optimal one if the entire solution space was considered. However regarding the real-time computation, and the heterogeneity of the considered criteria, dividing this problem into two steps makes sense in terms of reality and

applicability of the designed model, and in terms of quality of the found solution.

the initial number of weapons is greater than the number of oncoming targets.

with a target in the earliest possible date in order to avoid any risk.

In order to assign the targets to the available weapons, the Hungarian algorithm is used. The weapons and the targets are modelled as an asymmetric bipartite graph. The Figure 1 shows an example of possible assignment graph used. In the studied problem, it is assumed that

The quality of the proposed assignment is evaluated according to three different criteria: the capacity to propose an early fire, the width of the firing time window and the minimisation of the overflying of the defended area by our own assets for security purpose. These criteria

• the capability of the system to propose an early firing date, and then its ability to cope

• the width of the firing time windows represents the time that we have to cope with one target, then the larger is this firing time windows, the more time we have to propose one

**3.2. The choice of a two-step optimisation method**

**3.3. The weapon-target assignment**

respectively represent:

engagement solution,

trajectories.

The replicator equation (RE) is an Ordinary Differential Equation expressing the difference between the fitness of a strategy and the average fitness in the population. Lower payoffs (agents are minimizers) bring faster reproduction in accordance with Darwinian natural selection process.

$$\dot{p}\_i = -p\_i(e\_i \cdot A p^T - p \cdot A p^T)$$

RE for *i* = 1, ..., *m* describes the evolution of strategy frequencies *pi*. Moreover, for every initial strategy distribution *p*(0) ∈ *δm*, there is an unique solution *p*(*t*) ∈ *δ<sup>m</sup>* for all *t* ≥ 0 that satisfies the replicator equation. The replicator equation is the most widely used evolutionary dynamics. It was introduced for matrix games by Taylor and Jonker [26].

Note that this introducing to the EGT and the PSO comes from one of our previous study in [20].

## **3. The formulation of the DWTA: A target-based model**

A common approach to the DWTA problem based on the capabilities of the defence system to minimise the probability that a target can leak the proposed engagement plan. However, the problem dealt with in this study is slightly different from the classic DWTA. Whereas the classic DWTA is considering a multi-stage approach, the solved problem considers a continuous time where the targets are evolving in the space according to their own objectives and features. The proposed model has been designed to fit a specific requirement from industrial application, which explains this different approach.

The weapon system is defending an area from oncoming targets. This area is represented by a circle. All the weapons are disposed randomly within this range. In order to make the problem as general as possible, it is assumed that each weapon has its own velocity and own range. The targets are aiming the centre of the area to defend. The trajectories of the targets are designed by Bezier's curves in using 4 control points, all randomly drawn on the space, but the last point which is set to the centre of the area to defend. Thus, the problem presents a high diversity and can test the proposed method in the most of possible tricky cases. It is also assumed that the velocity of the targets and of the weapons are constant.

The assignment and the firing time sequence are computed in real-time in order to validate the reactivity of the studied algorithm. Which means that a timer is set at the beginning of the simulation, and the position of the targets evolves accordingly to this time.

#### **3.1. The engagement plan**

The engagement plan represents the solution space. An engagement plan is composed of one assignment weapon/target and completed by a date to fire. For example, if the following situation involves 3 weapons and 2 targets, a possible engagement plan *EP* could be:

$$EP(t) = \{ (\mathcal{W}\_1, T\_2, t + FT\_1); (\mathcal{W}\_3, T\_1, t + FT\_3) \}$$

Where t denotes the simulation time and the *Wi*, *i* ∈ {1, 2, 3} and *Tj*, *j* ∈ {1, 2} represent the weapon *i* and the target *j*. The variable *FTi* denotes the Firing Time computed for the weapon *i*. The engagement plan evolves accordingly to the situation and depends on the current simulation time and on the aerial situation. In this application, the engagement plan is recomputed every *P* seconds in order to make up a very reactive engagement plan capable of dealing with the trickier cases in which the targets are constantly changing their trajectories.

## **3.2. The choice of a two-step optimisation method**

6 Recent Advances on Meta-Heuristics and Their Application to Real Scenarios

with respect to a large *q* population.

selection process.

[20].

Another interpretation is *π*(*p*, *q*) being the fitness of a large population of agents playing pure strategies (*p* describing the agent proportion in each behaviour inside a population)

The replicator equation (RE) is an Ordinary Differential Equation expressing the difference between the fitness of a strategy and the average fitness in the population. Lower payoffs (agents are minimizers) bring faster reproduction in accordance with Darwinian natural

*p*˙*<sup>i</sup>* = −*pi*(*ei* · *Ap<sup>T</sup>* − *p* · *ApT*)

RE for *i* = 1, ..., *m* describes the evolution of strategy frequencies *pi*. Moreover, for every initial strategy distribution *p*(0) ∈ *δm*, there is an unique solution *p*(*t*) ∈ *δ<sup>m</sup>* for all *t* ≥ 0 that satisfies the replicator equation. The replicator equation is the most widely used evolutionary

Note that this introducing to the EGT and the PSO comes from one of our previous study in

A common approach to the DWTA problem based on the capabilities of the defence system to minimise the probability that a target can leak the proposed engagement plan. However, the problem dealt with in this study is slightly different from the classic DWTA. Whereas the classic DWTA is considering a multi-stage approach, the solved problem considers a continuous time where the targets are evolving in the space according to their own objectives and features. The proposed model has been designed to fit a specific requirement from

The weapon system is defending an area from oncoming targets. This area is represented by a circle. All the weapons are disposed randomly within this range. In order to make the problem as general as possible, it is assumed that each weapon has its own velocity and own range. The targets are aiming the centre of the area to defend. The trajectories of the targets are designed by Bezier's curves in using 4 control points, all randomly drawn on the space, but the last point which is set to the centre of the area to defend. Thus, the problem presents a high diversity and can test the proposed method in the most of possible tricky cases. It is

The assignment and the firing time sequence are computed in real-time in order to validate the reactivity of the studied algorithm. Which means that a timer is set at the beginning of

The engagement plan represents the solution space. An engagement plan is composed of one assignment weapon/target and completed by a date to fire. For example, if the following situation involves 3 weapons and 2 targets, a possible engagement plan *EP* could be:

*EP*(*t*) = {(*W*1, *T*2, *t* + *FT*1);(*W*3, *T*1, *t* + *FT*3)}

also assumed that the velocity of the targets and of the weapons are constant.

the simulation, and the position of the targets evolves accordingly to this time.

dynamics. It was introduced for matrix games by Taylor and Jonker [26].

**3. The formulation of the DWTA: A target-based model**

industrial application, which explains this different approach.

**3.1. The engagement plan**

Since the complexity of the presented problem grows exponentially with the number of targets and weapons, to design an algorithm capable of handling the real-time computation, but taking into account very diversified performance indexes, the choice of two different steps was natural. Lloyd [1] proved that the DWTA is a NP-complete problem. Therefore, it is hard to find an exact optimisation method capable of solving the DWTA problem in an exact way within a reasonable time. The reasonable time implies a high frequency which can enable the real-time application of the optimisation method. Note that the system must be able to provide results in real-time in the DWTA problem since the engagement changes as the targets keep evolving in the aerial space during the computation. With these reasons, using a heuristic approach providing suboptimal solutions in real-time could be the best way to handle the DWTA problem. One other problem is to be able to quantify the quality of one proposed solution: the performance index of the assignment, and the firing sequence cannot be evaluated in using the same performance criterion. Whereas the assignment is evaluated from the system point of view, the firing sequence is evaluated from the weapons features. Dividing the problem into two parts could lead to the modification of the solution space and the optimum solution could be not the same as the optimal one if the entire solution space was considered. However regarding the real-time computation, and the heterogeneity of the considered criteria, dividing this problem into two steps makes sense in terms of reality and applicability of the designed model, and in terms of quality of the found solution.

## **3.3. The weapon-target assignment**

In order to assign the targets to the available weapons, the Hungarian algorithm is used. The weapons and the targets are modelled as an asymmetric bipartite graph. The Figure 1 shows an example of possible assignment graph used. In the studied problem, it is assumed that the initial number of weapons is greater than the number of oncoming targets.

The quality of the proposed assignment is evaluated according to three different criteria: the capacity to propose an early fire, the width of the firing time window and the minimisation of the overflying of the defended area by our own assets for security purpose. These criteria respectively represent:


10.5772/53606

117

http://dx.doi.org/10.5772/53606

3, with

(*xt*, *yt*) position and the speed (*vtx* , *vty* ) of the target *t* (respectively (*xw*, *yw*) position and the speed (*vwx* , *vwy* ) of the weapon *w*) in the (*x*,*O*, *y*) plan. The entering point of the target *t* in the capture zone of the weapon *w* and the entering point of the defended area is computed in the same time as the *FTWw*/*<sup>t</sup>* and they are denoted by *Ptin* and *Ptout* . The initial position of

A Two-Step Optimisation Method for Dynamic Weapon Target Assignment Problem

Let *W* be the set of the available weapons and *T* the set of the oncoming targets. If *A* represents the assignments linking the vertices *W* to the targets *T*. *G* = (*W*, *T*, *A*) denotes

The weight of each edge is computed from the linear combination of the three criteria: earliest possible fire, width of the firing time windows and minimising the overfly of the defended

*f*1(*Ew*/*t*) = *EFFw*/*t*,(*w* ∈ *W*),(*t* ∈ *T*)

*f*2(*Ew*/*t*) = *LFFw*/*<sup>t</sup>* − *EFFw*/*t*,(*w* ∈ *W*),(*t* ∈ *T*)

*EFFw*/*<sup>t</sup>* denotes the earliest feasible fire for the weapon *w* on the target *t*. The latest feasible

*f*3(*Ew*/*t*) = *d* (*Ptout* , *Pw*<sup>0</sup> )

Here the function *d*(*P*1, *P*2) represents the Euclidean distance function between the point *P*<sup>1</sup>

Then, the global weight of the assignment *Ew*/*<sup>t</sup>* is the linear combination of the three functions described above: *H*(*Ew*/*t*) = *α*<sup>1</sup> *f*1(*Ew*/*t*) + *α*<sup>2</sup> *f*2(*Ew*/*t*) + *α*<sup>3</sup> *f*3(*Ew*/*t*), where *H*(*Ew*/*t*) denotes the weighting function of the assignment *Ew*/*<sup>t</sup>* and (*α*1, *α*2, *α*3) ∈ [0, 1]

> *E*1/1 *E*2/1 *E*3/1 ... *E*|*W*|/1 *E*1/2 *E*2/2 *E*3/2 ... *E*|*W*|/2

*E*1/|*T*<sup>|</sup> *E*2/|*T*<sup>|</sup> *E*3/|*T*<sup>|</sup> ... *E*|*W*|/|*T*<sup>|</sup>

. ... .

. . The cost matrix used for the Hungarian algorithm has the following form:

. . . . . . . .

As mentioned, *EFFw*/*<sup>t</sup>* denotes the earliest feasible fire for the weapon *w* on the target *t*.

the weapon *w* is denoted by *Pw*<sup>0</sup> = (*xw*<sup>0</sup> , *yw*<sup>0</sup> ).

*3.5.1. The assignment part: Hungarian algorithm*

area. These criteria are represented as follows:

fire for the weapon *w* on the target *t* is denoted by *LFFw*/*t*.

and the point *P*2. This criterion is shown in the Figure 2.

*H* =


the complete bipartite graph.

*α*<sup>1</sup> + *α*<sup>2</sup> + *α*<sup>3</sup> = 1.

**Figure 1.** Example of asymmetric bipartite graph with *w* weapons and *t* targets

• limiting the overfly in our own area enables to cope with security problem in case of material failure.

#### **3.4. The sequencing of the firing time**

As soon as the weapons are assigned to the targets, the sequencing of the firing is computed with respect to the weapons properties (range, velocity) and the firing time windows as well.

In order to evaluate the quality of the proposed solution, the performance index is based on the reactivity of the algorithm, the respect of the system constraints and the avoidance of idle time when a firing is possible.

The system is subject to some technical constraints as a required time between two firing times, which depends on the system. In the designed simulator this time is fixed to 3 seconds.

#### **3.5. Mathematical modelling**

This section describes the mathematical modelling of each step followed to achieve the DWTA. The first step is the assignment of the targets to the weapons, and then the sequencing of the firing time to complete in the best possible way the destroying of all the threatening targets. The weapon-target assignment is done by using the graph theory, especially the Hungarian algorithm. The second part is done by integrating two approaches: the PSO and the EGT to make up an efficient real-time oriented algorithm to solve the firing sequence problem.

In the following section, *FTWw*/*<sup>t</sup>* denotes the set of the firing time windows (time windows in which a weapon *w* can be fired with a given probability to reach the target *t*). *EFFw*/*<sup>t</sup>* denotes the earliest feasible fire for the weapon *w* on the target *t*. The latest feasible fire for the weapon *w* on the target *t* is denoted by *LFFw*/*t*. *Ew*/*<sup>t</sup>* represents the edge linking the weapon *w* with the target *t*. The average speed of the weapon *w* is denoted by *Sw*. *Rt* and *Rw* denote the state of the target *t* (respectively the weapon *w*). The states are composed of the (*xt*, *yt*) position and the speed (*vtx* , *vty* ) of the target *t* (respectively (*xw*, *yw*) position and the speed (*vwx* , *vwy* ) of the weapon *w*) in the (*x*,*O*, *y*) plan. The entering point of the target *t* in the capture zone of the weapon *w* and the entering point of the defended area is computed in the same time as the *FTWw*/*<sup>t</sup>* and they are denoted by *Ptin* and *Ptout* . The initial position of the weapon *w* is denoted by *Pw*<sup>0</sup> = (*xw*<sup>0</sup> , *yw*<sup>0</sup> ).

#### *3.5.1. The assignment part: Hungarian algorithm*

8 Recent Advances on Meta-Heuristics and Their Application to Real Scenarios

**Figure 1.** Example of asymmetric bipartite graph with *w* weapons and *t* targets

**3.4. The sequencing of the firing time**

material failure.

time when a firing is possible.

**3.5. Mathematical modelling**

problem.

• limiting the overfly in our own area enables to cope with security problem in case of

As soon as the weapons are assigned to the targets, the sequencing of the firing is computed with respect to the weapons properties (range, velocity) and the firing time windows as well. In order to evaluate the quality of the proposed solution, the performance index is based on the reactivity of the algorithm, the respect of the system constraints and the avoidance of idle

The system is subject to some technical constraints as a required time between two firing times, which depends on the system. In the designed simulator this time is fixed to 3 seconds.

This section describes the mathematical modelling of each step followed to achieve the DWTA. The first step is the assignment of the targets to the weapons, and then the sequencing of the firing time to complete in the best possible way the destroying of all the threatening targets. The weapon-target assignment is done by using the graph theory, especially the Hungarian algorithm. The second part is done by integrating two approaches: the PSO and the EGT to make up an efficient real-time oriented algorithm to solve the firing sequence

In the following section, *FTWw*/*<sup>t</sup>* denotes the set of the firing time windows (time windows in which a weapon *w* can be fired with a given probability to reach the target *t*). *EFFw*/*<sup>t</sup>* denotes the earliest feasible fire for the weapon *w* on the target *t*. The latest feasible fire for the weapon *w* on the target *t* is denoted by *LFFw*/*t*. *Ew*/*<sup>t</sup>* represents the edge linking the weapon *w* with the target *t*. The average speed of the weapon *w* is denoted by *Sw*. *Rt* and *Rw* denote the state of the target *t* (respectively the weapon *w*). The states are composed of the Let *W* be the set of the available weapons and *T* the set of the oncoming targets. If *A* represents the assignments linking the vertices *W* to the targets *T*. *G* = (*W*, *T*, *A*) denotes the complete bipartite graph.

The weight of each edge is computed from the linear combination of the three criteria: earliest possible fire, width of the firing time windows and minimising the overfly of the defended area. These criteria are represented as follows:

$$f\_1(E\_{w/t}) = EFF\_{w/t} \ (w \in W) \ (t \in T)^\*$$

As mentioned, *EFFw*/*<sup>t</sup>* denotes the earliest feasible fire for the weapon *w* on the target *t*.

$$f\_2(E\_{w/t}) = LFF\_{w/t} - EFF\_{w/t'} (w \in W)\_t (t \in T)$$

*EFFw*/*<sup>t</sup>* denotes the earliest feasible fire for the weapon *w* on the target *t*. The latest feasible fire for the weapon *w* on the target *t* is denoted by *LFFw*/*t*.

$$f\_{\mathfrak{Z}}(E\_{w/t}) = d\left(P\_{t\_{out}}, P\_{w\_0}\right)$$

Here the function *d*(*P*1, *P*2) represents the Euclidean distance function between the point *P*<sup>1</sup> and the point *P*2. This criterion is shown in the Figure 2.

Then, the global weight of the assignment *Ew*/*<sup>t</sup>* is the linear combination of the three functions described above: *H*(*Ew*/*t*) = *α*<sup>1</sup> *f*1(*Ew*/*t*) + *α*<sup>2</sup> *f*2(*Ew*/*t*) + *α*<sup>3</sup> *f*3(*Ew*/*t*), where *H*(*Ew*/*t*) denotes the weighting function of the assignment *Ew*/*<sup>t</sup>* and (*α*1, *α*2, *α*3) ∈ [0, 1] 3, with *α*<sup>1</sup> + *α*<sup>2</sup> + *α*<sup>3</sup> = 1.

The cost matrix used for the Hungarian algorithm has the following form:

$$H = \begin{pmatrix} E\_{1/1} & E\_{2/1} & E\_{3/1} & \dots & E\_{|W|/1} \\ E\_{1/2} & E\_{2/2} & E\_{3/2} & \dots & E\_{|W|/2} \\ \vdots & \vdots & \vdots & \dots & \vdots \\ E\_{1/|T|} & E\_{2/|T|} & E\_{3/|T|} & \dots & E\_{|W|/|T|} \end{pmatrix}$$


10.5772/53606

119

http://dx.doi.org/10.5772/53606

*f*4(*FS*) =

Where, *FT* denotes the firing time of the weapon assigned to the target *t*.

*f*5(*FS*) =

*f*6(*FS*) =

value in order to avoid infeasible solution.

weapon *i*, otherwise *ci* = 0.

available weapons.

current particle.

is obtained as:

with the EGT.

**4. The proposed method**

**4.1. A two step-method**

*T* ∑ *t*=1 *FTt*

A Two-Step Optimisation Method for Dynamic Weapon Target Assignment Problem

The second criterion evaluates the feasibility of the proposed solution to respect the short time delay due to the system constraints. This criterion is based on the presence of constraint violations. When any of the constraints is violated, the proposed solution takes the maximum

> *W* ∑ *w*=1

The vector Conflict = [*ci*], *i* = {1, . . . , |*W*|} with *ci* = 1 if there is a constraint violation by the

The third and last criterion is based on the idle time of the system. This criterion enables to avoid the inactivity of the system if there are possible fires by the current time. In the best case, this value should be reduced to the time constraint multiplied by the number of

> *W*−1 ∑ *w*=1

Note that the *FS* vector is sorted before computing this performance index function to the

When all the criteria are computed, the global performance of the proposed firing sequence

*<sup>F</sup>*(*FS*) = (*f*4(*FS*) + <sup>1</sup>). *<sup>f</sup>*6(*FS*) if *<sup>f</sup>*5(*FS*) = <sup>0</sup>

The proposed method is based on the consecutive use of the Hungarian algorithm to solve the assignment problem before determining the fire sequencing using the PSO combined

As described on the Flowchart 3, the two-step process computes first the optimal assignment of the targets to the weapons, then in a second time the optimal firing sequence is determined.

Conflict(*w*)

(*FSw*<sup>+</sup><sup>1</sup> − *FSw*)

+∞ if *f*5(*FS*) �= 0

**Figure 2.** Representation of the overflying criterion. The used value is the Euclidean distance between the entering point of the target in the area to defend and the initial position of the weapon.

#### *3.5.2. The firing time sequencing: EGPSO*

As described in the section 2.2, the EGPSO process is based on the combination of the PSO algorithm combined to the EGT in order to increase the convergence speed [27]. In this section, *FS* = [*FTi*], *i* = {1, . . . , *w*} denotes a firing sequence for the *w* selected weapons from the previous assignment and *FTi* represents the firing time of the weapon *i* (*i* ≤ |*W*|). In the proposed model, *FS* represents one particle composed by the set of the firing times for each weapon. Since the solution space is composed by the firing time windows, it can be very heterogeneous in terms of length along each dimension. In order to avoid an unequal exploration of the solution space, the normalisation over the solution space is operated. Thus, the solution space is reduced to a [0, 1] <sup>|</sup>*W*<sup>|</sup> hypercube and enables a homogeneous exploring by the particles.

In order to evaluate the performance of a proposed solution, the global performance index is based on the reactivity of the algorithm, the respect of the system constraints and the avoidance of idle time when a firing is possible. The global cost function is obtained in multiplying each criterion. The multiplication is selected to consider evenly all the criteria. Thus, if one criterion is not respected by the proposed engagement plan, the cost function will decrease accordingly to the unsatisfied criterion.

The first performance index based on the time delay enables to quantify the reactivity of the system in summing the firing times. The function *f*<sup>1</sup> enables to express this criterion.

$$f\_4(FS) = \sum\_{t=1}^{T} FT\_t$$

Where, *FT* denotes the firing time of the weapon assigned to the target *t*.

The second criterion evaluates the feasibility of the proposed solution to respect the short time delay due to the system constraints. This criterion is based on the presence of constraint violations. When any of the constraints is violated, the proposed solution takes the maximum value in order to avoid infeasible solution.

$$f\_{\mathbb{S}}(FS) = \sum\_{w=1}^{W} \text{Conflict}(w).$$

The vector Conflict = [*ci*], *i* = {1, . . . , |*W*|} with *ci* = 1 if there is a constraint violation by the weapon *i*, otherwise *ci* = 0.

The third and last criterion is based on the idle time of the system. This criterion enables to avoid the inactivity of the system if there are possible fires by the current time. In the best case, this value should be reduced to the time constraint multiplied by the number of available weapons.

$$f\_6(FS) = \sum\_{w=1}^{W-1} (FS\_{w+1} - FS\_w)^2$$

Note that the *FS* vector is sorted before computing this performance index function to the current particle.

When all the criteria are computed, the global performance of the proposed firing sequence is obtained as:

$$F(FS) = \begin{cases} (f\_4(FS) + 1).f\_6(FS) \text{ if } f\_5(FS) = 0\\ +\infty & \text{if } f\_5(FS) \neq 0 \end{cases}$$

## **4. The proposed method**

10 Recent Advances on Meta-Heuristics and Their Application to Real Scenarios

the target in the area to defend and the initial position of the weapon.

*3.5.2. The firing time sequencing: EGPSO*

the solution space is reduced to a [0, 1]

will decrease accordingly to the unsatisfied criterion.

by the particles.

**Figure 2.** Representation of the overflying criterion. The used value is the Euclidean distance between the entering point of

As described in the section 2.2, the EGPSO process is based on the combination of the PSO algorithm combined to the EGT in order to increase the convergence speed [27]. In this section, *FS* = [*FTi*], *i* = {1, . . . , *w*} denotes a firing sequence for the *w* selected weapons from the previous assignment and *FTi* represents the firing time of the weapon *i* (*i* ≤ |*W*|). In the proposed model, *FS* represents one particle composed by the set of the firing times for each weapon. Since the solution space is composed by the firing time windows, it can be very heterogeneous in terms of length along each dimension. In order to avoid an unequal exploration of the solution space, the normalisation over the solution space is operated. Thus,

In order to evaluate the performance of a proposed solution, the global performance index is based on the reactivity of the algorithm, the respect of the system constraints and the avoidance of idle time when a firing is possible. The global cost function is obtained in multiplying each criterion. The multiplication is selected to consider evenly all the criteria. Thus, if one criterion is not respected by the proposed engagement plan, the cost function

The first performance index based on the time delay enables to quantify the reactivity of the system in summing the firing times. The function *f*<sup>1</sup> enables to express this criterion.

<sup>|</sup>*W*<sup>|</sup> hypercube and enables a homogeneous exploring

The proposed method is based on the consecutive use of the Hungarian algorithm to solve the assignment problem before determining the fire sequencing using the PSO combined with the EGT.

#### **4.1. A two step-method**

As described on the Flowchart 3, the two-step process computes first the optimal assignment of the targets to the weapons, then in a second time the optimal firing sequence is determined.

10.5772/53606

121

http://dx.doi.org/10.5772/53606

matrix is composed by the mean fitness of the particles composing each swarm. Let *S* be the

After one iteration using each strategy successively, the payoff matrix consists of the mean

2

*π*(*s*3) + *π*(*s*2)

The coefficients *π*(*si*), *i* ∈ {1, 2, 3} are the mean value of the swarm after using the pure strategy *si*. The evolutionary game process used to converge to the evolutionary stable strategy is the *replicator dynamic* described in [20]. As soon as the population is stabilised, the proposed algorithm stop running the replicator dynamic. This ESS gives the stable strategy rate, generally composed by a mix of the strategies *s*1, *s*2, and *s*3. Then, the final step uses

The principle of the method is described on the Flowchart 4 and by the following process

(a) Random selection of particles following the classical PSO process (exploration) and

(b) Classic iteration of the PSO in using only one strategy for each swarm (inertial,

(c) Computation of the payoff matrix in computing the mean value of the swarm in using

the particles following the EGPSO (increase computational speed).

(d) Find the evolutionary stable strategy depending on the payoff matrix (e) Classic iteration of the PSO using the previously found coefficients

<sup>2</sup> *<sup>π</sup>*(*s*2) *<sup>π</sup>*(*s*2) + *<sup>π</sup>*(*s*3)

<sup>2</sup> *<sup>π</sup>*(*s*3)

*π*(*s*1) + *π*(*s*3) 2

A Two-Step Optimisation Method for Dynamic Weapon Target Assignment Problem

2

*<sup>π</sup>*(*s*1) *<sup>π</sup>*(*s*1) + *<sup>π</sup>*(*s*2)

set of the available strategies *si*, *i* ∈ {1, 2, 3} which are as follows:

*π*(*s*2) + *π*(*s*1)

*π*(*s*3) + *π*(*s*1) 2

*s*1: Use of the pure strategy *inertia s*2: Use of the pure strategy *individual s*3: Use of the pure strategy *social*

value of the swarm. *A* denotes this payoff matrix:

these rates as coefficients in the PSO algorithm.

2. For a maximum number of iterations

(f) Check if the swarm is stabilised

• If YES, restart the swarm like at the step 1 • If NOT, keep running the algorithm

In the presented simulation, the PSO parameters are defined as:

individual, social)

3. Obtain the optimal solution

the strategies

1. Initialisation of the swarm in position and velocity

step by step:

Π =

**Figure 3.** Representation of the two-step method to solve the DWTA.

## **4.2. The Hungarian algorithm**

The assignment of the targets to the weapons is realised in using the Hungarian algorithm [21]. The section 3.5.1 states all the required details enabling to understand the principles of the used method. Since in real scenarii the number of targets is only rarely the same as the number of weapons, the Hungarian algorithm designed for asymmetric bipartite graphs is used. The following parameters are used to determine the best assignment: the cost matrix has a |*T*|×|*W*| form in order to assign all the targets and the coefficients of this cost matrix are determined in using the equations described in 3.5.1.

#### **4.3. The integration of the particle swarm optimisation with the evolutionary game theory**

There are two main steps in this approach, the first one is the movement of the swarm in using only, first the inertia, then only the individual component, then only the social component. From the obtained results of the movement of the three swarms, the payoff matrix is composed by the mean fitness of the particles composing each swarm. Let *S* be the set of the available strategies *si*, *i* ∈ {1, 2, 3} which are as follows:

*s*1: Use of the pure strategy *inertia*

12 Recent Advances on Meta-Heuristics and Their Application to Real Scenarios

**Figure 3.** Representation of the two-step method to solve the DWTA.

are determined in using the equations described in 3.5.1.

The assignment of the targets to the weapons is realised in using the Hungarian algorithm [21]. The section 3.5.1 states all the required details enabling to understand the principles of the used method. Since in real scenarii the number of targets is only rarely the same as the number of weapons, the Hungarian algorithm designed for asymmetric bipartite graphs is used. The following parameters are used to determine the best assignment: the cost matrix has a |*T*|×|*W*| form in order to assign all the targets and the coefficients of this cost matrix

**4.3. The integration of the particle swarm optimisation with the evolutionary**

There are two main steps in this approach, the first one is the movement of the swarm in using only, first the inertia, then only the individual component, then only the social component. From the obtained results of the movement of the three swarms, the payoff

**4.2. The Hungarian algorithm**

**game theory**


After one iteration using each strategy successively, the payoff matrix consists of the mean value of the swarm. *A* denotes this payoff matrix:

$$
\Pi = \begin{pmatrix}
\pi(s\_1) & \frac{\pi(s\_1) + \pi(s\_2)}{2} & \frac{\pi(s\_1) + \pi(s\_3)}{2} \\
\frac{\pi(s\_2) + \pi(s\_1)}{2} & \pi(s\_2) & \frac{\pi(s\_2) + \pi(s\_3)}{2} \\
\frac{\pi(s\_3) + \pi(s\_1)}{2} & \frac{\pi(s\_3) + \pi(s\_2)}{2} & \pi(s\_3)
\end{pmatrix}
$$

The coefficients *π*(*si*), *i* ∈ {1, 2, 3} are the mean value of the swarm after using the pure strategy *si*. The evolutionary game process used to converge to the evolutionary stable strategy is the *replicator dynamic* described in [20]. As soon as the population is stabilised, the proposed algorithm stop running the replicator dynamic. This ESS gives the stable strategy rate, generally composed by a mix of the strategies *s*1, *s*2, and *s*3. Then, the final step uses these rates as coefficients in the PSO algorithm.

The principle of the method is described on the Flowchart 4 and by the following process step by step:

	- (a) Random selection of particles following the classical PSO process (exploration) and the particles following the EGPSO (increase computational speed).
	- (b) Classic iteration of the PSO in using only one strategy for each swarm (inertial, individual, social)
	- (c) Computation of the payoff matrix in computing the mean value of the swarm in using the strategies
	- (d) Find the evolutionary stable strategy depending on the payoff matrix
	- (e) Classic iteration of the PSO using the previously found coefficients
	- (f) Check if the swarm is stabilised
		- If YES, restart the swarm like at the step 1
		- If NOT, keep running the algorithm

In the presented simulation, the PSO parameters are defined as:

10.5772/53606

123

http://dx.doi.org/10.5772/53606

• 50 particles are used to explore the solution space.

is considered as

**5. Results and comments**

missing fire for example).

Square of 50000 m by 50000 m

on the centre of the space.

16 Weapons vs. 12 Targets. **Condition of engagement success:**

operating a successful shot is 75%.

**The aerial space:**

**Weapons**

**Targets**

(0, 0).

**The initial conditions:**

• The maximum distance travelled by one particle in one iteration is limited to 1/10 along each dimension. Notice that since the solution space has been normalised, the maximum

A Two-Step Optimisation Method for Dynamic Weapon Target Assignment Problem

• In order to be able to be competitive in real-time, the exit criterion is a defined time of

In order to enable a quick convergence to the optimal vector rate of the PSO coefficients, the EGT process is launched in using as payoff matrix Π described in the section 4.3. The replicator equation is computed over 500 hundred generations, and then the obtained result

In this section, the efficiency of the proposed approach is analysed. After running 100 times a simulation, the number of experiences that the mission is successfully achieved is compared to the number of times it fails. Then, in a second time the evolution of the assignment is studied in analysing the target motions and the proposed engagement plan. The study ends with the analysis of the human operator point of view in order to determine if the proposed algorithm can be reliable and stable for the operator. By stable, it is assumed that the operator can have a global overview of the next engagement to execute in advance, and that this plan won't change if there are no major changes in the aerial situation (suppressed enemy or

velocity enables an homogeneous exploration of the solution space.

In the presented simulator, the used parameters are set up as follows:

The speed is randomly drawn between 50 m/s and 900 m/s.

The initial position is within a radius of 3000 m around the central objective

The range of each weapon is randomly drawn between 10000 meters and 15000 meters.

The initial position is set up between 30000 m and 50000 m from the main objective located

The trajectories that the targets are following are modelled in using a Bezier curve defined by 4 control points. The last control point is automatically set as the centre of the space

The success of an engagement one weapon on one target is determined in drawing one random number. If this number is greater than a determined value, then the shoot is considered as a success. Otherwise, it is considered that the target avoids the weapon. In this simulator this value is arbitrary fixed to 0.25, which means that the probability of

2500 ms, after that the best found solution is considered as the optimal one.

**Figure 4.** Details on the method designed to mix EGT and PSO.


In order to enable a quick convergence to the optimal vector rate of the PSO coefficients, the EGT process is launched in using as payoff matrix Π described in the section 4.3. The replicator equation is computed over 500 hundred generations, and then the obtained result is considered as

## **5. Results and comments**

In this section, the efficiency of the proposed approach is analysed. After running 100 times a simulation, the number of experiences that the mission is successfully achieved is compared to the number of times it fails. Then, in a second time the evolution of the assignment is studied in analysing the target motions and the proposed engagement plan. The study ends with the analysis of the human operator point of view in order to determine if the proposed algorithm can be reliable and stable for the operator. By stable, it is assumed that the operator can have a global overview of the next engagement to execute in advance, and that this plan won't change if there are no major changes in the aerial situation (suppressed enemy or missing fire for example).

In the presented simulator, the used parameters are set up as follows:

#### **The aerial space:**

Square of 50000 m by 50000 m

#### **Weapons**

14 Recent Advances on Meta-Heuristics and Their Application to Real Scenarios

**Figure 4.** Details on the method designed to mix EGT and PSO.

The initial position is within a radius of 3000 m around the central objective The range of each weapon is randomly drawn between 10000 meters and 15000 meters.

#### **Targets**

The initial position is set up between 30000 m and 50000 m from the main objective located on the centre of the space.

The trajectories that the targets are following are modelled in using a Bezier curve defined by 4 control points. The last control point is automatically set as the centre of the space (0, 0).

The speed is randomly drawn between 50 m/s and 900 m/s.

#### **The initial conditions:**

16 Weapons vs. 12 Targets.

#### **Condition of engagement success:**

The success of an engagement one weapon on one target is determined in drawing one random number. If this number is greater than a determined value, then the shoot is considered as a success. Otherwise, it is considered that the target avoids the weapon. In this simulator this value is arbitrary fixed to 0.25, which means that the probability of operating a successful shot is 75%.

10.5772/53606

125

http://dx.doi.org/10.5772/53606

**Figure 6.** The upper graph illustrates the variation in the assignment process over the time. The regularity of the proposed assignment can be noticed, especially as long as the aerial situation does not change (no target are suppressed). The black vertical lines highlight these phases. The lower graph shows the evolution of the proposed firing time to engage the target over

A Two-Step Optimisation Method for Dynamic Weapon Target Assignment Problem

been killed, then it denotes a change in the aerial situation. During the different highlighted phases, the assignment presents some interesting features as the regularity over the time when the aerial situation keep being similar. The lower graph on the Figure 6 represents the evolution of the firing time for each target over the time. The vertical lines have the same meaning as the upper graph and denotes a change in the aerial situation like, for example, a suppressed enemy or an unsuccessful fire. This second graph highlights the continuity of the proposed firing sequence over the time. It is shown that the operator can not only approve the firing sequence in executing the firing, but the operator can follow the entire scenario and can anticipate the upcoming events. The Figure 7 focuses on the real time aspect in focusing only on the operator point of view. Indeed this Figure represents a zoom on the 25 last seconds before firing the weapons. The horizontal dash line illustrates a time limit of 5

**Figure 7.** This graph represents a zoom on the final instruction of the operator to execute the firing of the weapon as soon as

the proposed firing time is within 5 seconds of the current time. This limit is illustrated by the horizontal dash line.

seconds from which the operator can execute the firing.

the time.

**Figure 5.** Representation of a possible initialisation of trajectory and weapon position. The triangle marker represents the initial position of the target. The dot line is the trajectory that the target will follow to reach its goal. The continuous line represents the area that we are defending and the cross marker surrounded by a dot line denote the defending weapon and its capture zone.

The Figure 5 shows a possible initialisation of a scenario. Note that if the trajectory is a priori known by the target, the defending side has no information at all but the final point of the target and its current position.

The analysis of the evolution of the assignment of the weapons to the oncoming targets clearly shows stability over the simulation time as long as there are no major change in the scenario. A major change in the scenario can be qualified by the suppression of one enemy which leads to the reconsideration of the entire scenario. Otherwise, the proposed method clearly shows a good stability over the simulation time which is required in the presented case. Considering the presence of a human operator having the final decision making and using this method as a help in the decision making process, it is important for the proposed engagement to be continuous over the time when the aerial situation does not vary dramatically. The upper graph of the Figure 6 displays the assignment of the target *t* over the number of iterations. The vertical lines identify the instants when a target has

16 Recent Advances on Meta-Heuristics and Their Application to Real Scenarios

**Figure 5.** Representation of a possible initialisation of trajectory and weapon position. The triangle marker represents the initial position of the target. The dot line is the trajectory that the target will follow to reach its goal. The continuous line represents the area that we are defending and the cross marker surrounded by a dot line denote the defending weapon and its capture

The Figure 5 shows a possible initialisation of a scenario. Note that if the trajectory is a priori known by the target, the defending side has no information at all but the final point of the

The analysis of the evolution of the assignment of the weapons to the oncoming targets clearly shows stability over the simulation time as long as there are no major change in the scenario. A major change in the scenario can be qualified by the suppression of one enemy which leads to the reconsideration of the entire scenario. Otherwise, the proposed method clearly shows a good stability over the simulation time which is required in the presented case. Considering the presence of a human operator having the final decision making and using this method as a help in the decision making process, it is important for the proposed engagement to be continuous over the time when the aerial situation does not vary dramatically. The upper graph of the Figure 6 displays the assignment of the target *t* over the number of iterations. The vertical lines identify the instants when a target has

zone.

target and its current position.

**Figure 6.** The upper graph illustrates the variation in the assignment process over the time. The regularity of the proposed assignment can be noticed, especially as long as the aerial situation does not change (no target are suppressed). The black vertical lines highlight these phases. The lower graph shows the evolution of the proposed firing time to engage the target over the time.

been killed, then it denotes a change in the aerial situation. During the different highlighted phases, the assignment presents some interesting features as the regularity over the time when the aerial situation keep being similar. The lower graph on the Figure 6 represents the evolution of the firing time for each target over the time. The vertical lines have the same meaning as the upper graph and denotes a change in the aerial situation like, for example, a suppressed enemy or an unsuccessful fire. This second graph highlights the continuity of the proposed firing sequence over the time. It is shown that the operator can not only approve the firing sequence in executing the firing, but the operator can follow the entire scenario and can anticipate the upcoming events. The Figure 7 focuses on the real time aspect in focusing only on the operator point of view. Indeed this Figure represents a zoom on the 25 last seconds before firing the weapons. The horizontal dash line illustrates a time limit of 5 seconds from which the operator can execute the firing.

**Figure 7.** This graph represents a zoom on the final instruction of the operator to execute the firing of the weapon as soon as the proposed firing time is within 5 seconds of the current time. This limit is illustrated by the horizontal dash line.

In order to test the efficiency of the proposed method over different scenario, the designed experience has been launched 100 times and the final result archived. The pie diagram 8 shows the number of times that the proposed method achieved its goal versus the number of time it fails. The analysis of this result shows that the proposed algorithm successfully achieved its mission in 96% of the cases. If we look into details the causes of these failures, we can notice that 3 of the 4 failures was due to the lack of available weapons. Which means that the method does not achieved its goal because of the probability. Indeed, with PK threshold fixed to *PK* = 0.90 and 16 available weapons versus 12 targets, we have an estimate failure rate of approximatively 2 %. This last result comes from the binomial distribution, where the probability of getting exactly *T* success in *W* trials is given by:

10.5772/53606

127

http://dx.doi.org/10.5772/53606

results have shown the efficiency of the proposed two-step approach in various cases. The proposed algorithm achieves its objective in 96% for the given scenarii which include random simulation parameters selected for the generality of the senarii. Note that from a probability study on this application, with the chosen simulation parameters, 2% of the scenarii was expected to be failed simply because of the associated probability laws based on a Binomial

A Two-Step Optimisation Method for Dynamic Weapon Target Assignment Problem

Cédric Leboucher1,<sup>⋆</sup>, Hyo-Sang Shin2, Patrick Siarry3, Rachid Chelouah4,

<sup>⋆</sup> Address all correspondence to: cedric.leboucher@mbda-systems.com

3 Université Paris-Est Créteil (UPEC), LISSI (EA 3956), Créteil, France

2 Cranfield University, School of Engineering, College Road, Cranfield, Bedford, UK

[1] S. P. Lloyd and H. S. Witsenhausen. Weapon allocation is np-complete. In *Proceeding IEEE Summer Simulation Conference*, page 1054 âA ¸S 1058, Reno (USA), 1986. ˘

[2] P. A. Hosein and M. Athans. Preferential defense strategies. part i: The static case. Technical report, MIT Laboratory for Information and Decision Systems with partial

[3] S. Bisht. Hybrid genetic-simulated annealing algorithm for optimal weapon allocation

[4] A. Malhotra and R. K. Jain. Genetic algorithm for optimal weapon allocation in

[5] O. Karasakal. Air defense missile-target allocation models for a naval task group.

[6] F. Johansson and G. Falkman. Sward: System for weapon allocation research & development. *13th Conference on Information Fusion (FUSION)*, 1:1–7, July 2010.

[7] O. Kwon, K. Lee, D. Kang, and S. Park. A branch-and-price algorithm for a targeting

[8] H. Cai, J. Liu, Y. Chen, and H.Wang. Survey of the research on dynamic weapon-target

inmultilayer defence scenario. *Defence Sci. J.*, 54(3):395 – 405, 2004.

multilayer defence scenario. *Defence Sci. J.*, 51(3):285 – 293, 2001.

assignment problem. *J. Syst. Eng. Electron.*, 17(3):559 – 565, 2006.

[9] E. Wacholder. A neural network-based optimization algorithm for the static weapon-target assignment problem. *ORSA J. Comput.*, 1(4):232 – 246, 1989.

1 MBDA France, 1 av. Reaumur, Le Plessis Robinson, France

4 L@ris, EISTI, Avenue du Parc, Cergy-Pontoise, France

Stéphane Le Ménec<sup>1</sup> and Antonios Tsourdos<sup>2</sup>

support, Cambridge (USA), 1990.

*Comput. Oper. Res.*, 35:1759 – 1770, 2008.

problem. *Naval Res. Log.*, 54:732 – 741, 2007.

distribution.

**Author details**

**References**

$$P(T; W, PK) = \frac{W!}{T!(W - T)!} PK^T (1 - PK)^{W - T}$$

Thus, to solve this issue, two possible ways could be explored: first, the increasing of available weapons; second, using more accurate weapons. Although both of the proposed solutions can cope with this issue, it leads to increase the cost of the mission. Controlling this probability enables to optimise the used deployment to protect our area.

**Figure 8.** This bar diagram illustrates the number of time that the simulation is a success versus the number of time that it fails.

## **6. Conclusion**

In this chapter, a two-step optimisation method for the DWTA was proposed. Based on the successive use of the Hungarian algorithm, and a PSO combined with the EGT, the proposed algorithm shows reliable results in terms of performance and real-time computation. The proposed method is verified using one simulator designed to create random scenarii and to follow the normal evolution of the battlefield in real-time. The initialised scenario was composed of 16 weapons versus 12 targets. The stability of the assignment and the continuity of the firing sequence was analysed over the launch of 100 simulations. Regarding the probability of successfully achieved the mission, a short study about the binomial distribution has been done and could be helpful in the mission planning process to determine the optimal number of available weapons before the mission. The simulation results have shown the efficiency of the proposed two-step approach in various cases. The proposed algorithm achieves its objective in 96% for the given scenarii which include random simulation parameters selected for the generality of the senarii. Note that from a probability study on this application, with the chosen simulation parameters, 2% of the scenarii was expected to be failed simply because of the associated probability laws based on a Binomial distribution.

## **Author details**

18 Recent Advances on Meta-Heuristics and Their Application to Real Scenarios

probability of getting exactly *T* success in *W* trials is given by:

**6. Conclusion**

*<sup>P</sup>*(*T*; *<sup>W</sup>*, *PK*) = *<sup>W</sup>*!

this probability enables to optimise the used deployment to protect our area.

In order to test the efficiency of the proposed method over different scenario, the designed experience has been launched 100 times and the final result archived. The pie diagram 8 shows the number of times that the proposed method achieved its goal versus the number of time it fails. The analysis of this result shows that the proposed algorithm successfully achieved its mission in 96% of the cases. If we look into details the causes of these failures, we can notice that 3 of the 4 failures was due to the lack of available weapons. Which means that the method does not achieved its goal because of the probability. Indeed, with PK threshold fixed to *PK* = 0.90 and 16 available weapons versus 12 targets, we have an estimate failure rate of approximatively 2 %. This last result comes from the binomial distribution, where the

*T*!(*W* − *T*)!

Thus, to solve this issue, two possible ways could be explored: first, the increasing of available weapons; second, using more accurate weapons. Although both of the proposed solutions can cope with this issue, it leads to increase the cost of the mission. Controlling

**Figure 8.** This bar diagram illustrates the number of time that the simulation is a success versus the number of time that it fails.

In this chapter, a two-step optimisation method for the DWTA was proposed. Based on the successive use of the Hungarian algorithm, and a PSO combined with the EGT, the proposed algorithm shows reliable results in terms of performance and real-time computation. The proposed method is verified using one simulator designed to create random scenarii and to follow the normal evolution of the battlefield in real-time. The initialised scenario was composed of 16 weapons versus 12 targets. The stability of the assignment and the continuity of the firing sequence was analysed over the launch of 100 simulations. Regarding the probability of successfully achieved the mission, a short study about the binomial distribution has been done and could be helpful in the mission planning process to determine the optimal number of available weapons before the mission. The simulation

*PKT*(<sup>1</sup> − *PK*)

*W*−*T*

Cédric Leboucher1,<sup>⋆</sup>, Hyo-Sang Shin2, Patrick Siarry3, Rachid Chelouah4, Stéphane Le Ménec<sup>1</sup> and Antonios Tsourdos<sup>2</sup>

<sup>⋆</sup> Address all correspondence to: cedric.leboucher@mbda-systems.com

1 MBDA France, 1 av. Reaumur, Le Plessis Robinson, France

2 Cranfield University, School of Engineering, College Road, Cranfield, Bedford, UK

3 Université Paris-Est Créteil (UPEC), LISSI (EA 3956), Créteil, France

4 L@ris, EISTI, Avenue du Parc, Cergy-Pontoise, France

## **References**


[10] K. E. Grant. Optimal resource allocation using genetic algorithms. Technical report, Naval Research Laboratory, Washington (USA), 1993.

10.5772/53606

129

http://dx.doi.org/10.5772/53606

[23] M. Clerc. Discrete particle swarm optimization, illustred by the travelling salesman

algorithm. In *The 1997 IEEE International Conference on Systems, Man, and Cybernetics*,

A Two-Step Optimisation Method for Dynamic Weapon Target Assignment Problem

[25] J. Maynard-Smith. *Evolution and the theory of games*. Cambridge University Press, 1982.

[27] C. Leboucher, R. Chelouah, P. Siarry, and S. Le Ménec. A swarm intelligence method combined to evolutionary game theory applied to ressources allocation problem. In

problem. *New Optimization Techniques in Engineering*, 1:219–239, 2004.

[26] P. Taylor and Jonker L. Evolutionary stable strategies and game dynamics.

*International Conference on Swarm Intelligence*, Cergy (France), June 2011.

volume 5, pages 4104–4108, Orlando (USA), October 1997.

*Mathematical Bioscience*, 40:145–156, 1978.

[24] J. Kennedy and R.C. Eberhart. A discrete binary version of the particle swarm


[23] M. Clerc. Discrete particle swarm optimization, illustred by the travelling salesman problem. *New Optimization Techniques in Engineering*, 1:219–239, 2004.

20 Recent Advances on Meta-Heuristics and Their Application to Real Scenarios

Naval Research Laboratory, Washington (USA), 1993.

*Autom.*, pages 3401 – 3405, Dalian (China), 2006.

Way, WPAFB, OH, 2000.

support, Cambridge (USA), 1990.

Systemswith partial support, 1990.

*Logistics Quarterly*, 2:83 – 97, 1955.

662, 2010.

2008.

[10] K. E. Grant. Optimal resource allocation using genetic algorithms. Technical report,

[11] H. Lu, H. Zhang, X. Zhang, and R. Han. An improved genetic algorithm for target assignment optimization of naval fleet air defense. In *6th World Cong. Intell. Contr.*

[13] D. Blodgett, M. Gendreau, F. Guertin, and J. Y. Potvin. A tabu search heuristic for

[14] B. Xin, J. Chen, J. Zhang, L. Dou, and Z. Peng. Efficient decision makings for dynamic weapon-target assignment by virtual permutation and tabu search heuristics. *IEEE transaction on systems, man, and cybernetics - Part C: Application and reviews*, 40(6):649 –

[15] P. A. Hosein and M. Athans. Preferential defense strategies. part ii: The dynamic case. Technical report, MIT Laboratory for Information and Decision Systems with partial

problems with vulnerable c2 nodes. Technical report, MIT Laboratory for Information

[16] P. A. Hosein, J. T. Walton, and M. Athans. Dynamic weapon-target assignment

[17] P. A. Hosein and M. Athans. Some analytical results for the dynamic weapon-target allocation problem. Technical report, MIT Laboratory for Information and Decision

[18] T. Sikanen. Solving weapon target assignment problem with dynamic programming. Technical report, Mat-2.4108 Independent research projects in applied mathematics,

[20] C. Leboucher, R. Chelouah, P. Siarry, and S. Le Ménec. A swarm intelligence method combined to evolutionary game theory applied to the resources allocation problem.

[21] H.W. Kuhn. The hungarian method for the assignment problem. *Naval Research*

[22] J. Kennedy and R. C. Eberhart. Particle swarm optimization. In *Proceedings of IEEE*

*International Conference on Neural Networks*, pages 1942 – 1948, Piscataway (USA), 1995.

weapon-target allocation problem with decreasing weapons and targets. In *IEEE Congr.*

and Decision Systems with partial support, Cambridge (USA), 1988.

[19] L. Wu, C. Xing, F. Lu, and P. Jia. An anytime algorithm applied to dynamic

*International Journal of Swarm Intelligence Research*, 3(2):20 – 38, 2012.

*Evol. Comput.*, pages 3755 – 3759, Hong Kong (China), 2008.

resource management in naval warfare. *J. Heur.*, 9:145 – 169, 2003.

[12] A. C. Cullenbine. A taboo search approach to the weapon assignment model. Master's thesis, Department of Operational Sciences, Air Force Institute of Technology, Hobson


## *Edited by Javier Del Ser*

This book aims at attracting the interest of researchers and practitioners around the applicability of meta-heuristic algorithms to practical scenarios arising from different knowledge disciplines. Emphasis is placed on evolutionary algorithms and swarm intelligence as computational means to efficiently balance the tradeoff between optimality of the produced solutions and the complexity derived from their estimation. In summary, this book serves as a good start point for early-stage investigators in the initial steps of their research on meta-heuristics, grounded on both a thorough literature review and the practical orientation of its contents.

Recent Advances on Meta-Heuristics and Their Application to Real Scenarios

Recent Advances on

Meta-Heuristics and Their

Application to Real Scenarios

*Edited by Javier Del Ser*

Photo by scyther5 / iStock