A Bi-objective Simulation-optimization Approach for Solving a No-wait two Stages Flexible Flow Shop Scheduling Problem with Rework Ability

The paper suggests a new rule; called no-wait process. The rule has two stages, and is a flexible flow shop scheduling. The process is the subject to maximize tardiness while minimizing the makespan. This hybrid flow shop problem is known to be NP-hard. Therefore, we come to first, Non-dominated Sorting Genetic Algorithm (NSGA-II), then, Multi-Objective Imperialist Competitive Algorithm (MOICA) and finally, Pareto Archive Evolutionary Strategy (PAES) as three multi-objective Pareto based metaheuristic optimization methods. They are developed to solve the problem to approximately figure out optimal Pareto front. The method is investigated in several problems that differed in size and terms of relative percentage deviation of performance metrics. The conclusion, developed by this method is the most efficient and practicable algorithm at the end.


Introduction
Scheduling is arranging and planning sequence of event to complete the work.The goal is to utilize the resources optimally while reaching targets.This field has been attracted by many scholars during recent years (Mosheiov and Sidney 2010, Vallada and Ruiz 2011, Shafaei et al. 2011, Rabiee et al. 2012, Jolai et al. 2013, Ullrich 2013, Jolai et al. 2014, Yang et al. 2014, Tayebi Araghi et al. 2014, Xu et al. 2015, Rabiee et al. 2016, Nesello et al. 2017).The flow shop problem is an important subject in scheduling.One of the reasons is that most of the manufacturing systems follow batch shops, flow shop or semi flow shop routings (Baker 1974, Johnson 1954, Lin et al. 2006).Flexible flow shop or hybrid flow shop, is the one of the most important classes of scheduling which also defined as a flow shop with parallel machines and flexible flow line.For a literature review in this area, the readers are referred to those of (Richard and Zhang 1999, Ruiz and Vazquez-Rodriguez 2010, Ribas et al. 2010) In a no-wait flow shop, the jobs are processed from one machine to the next one without waiting time (Huang et al. 2009).Suppose there are some sequences of jobs, some specified procedures in each that are processed by machine in disciplinary order.No pause or interruption is supposed to occur in the line.In other word, when the process starts, there is no stop for object before or after each machine.No repetition is allowed, one job at a given time by a machine.Therefore, when needed, the start of a job on the first machine must be delayed in order to meet the no-wait requirement (Tasgetiren et al. 2007).Nagano et al. (2013) examines the m machine no-wait flow shop problem with setup times of a job separated from its processing time.The performance measure considered is the makespan.The hybrid metaheuristic Evolutionary Cluster Search (ECS_NSL) is employed to solve this scheduling problem.Nagano et al. (2015) addressed the problem of scheduling jobs in a no -wait flow shop with sequence-dependent setup times with the objective of minimizing the total flow time.As this problem is well-known for being NP-hard, they presented a new constructive heuristic, named QUARTS, in order to obtain good approximate solutions in a short CPU time.Samarghandi and ElMekkawy (2012) studied the problem of no-wait flow shop and proposed two frameworks based on genetic algorithm and particle swarm optimization to deal with the problem.Samarghandi and ElMekkawy (2014) proposed PSO algorithm to solve the problem of scheduling a no-wait flow-shop system with sequence-dependent set-up times.The application of this problem can be found in industries such as chemical industry, steel production, just-in-time manufacturing process, service industries, etc (Rajendran 1994, Grabowski and Pempera 2000, Raaymakers and Hoogeveen 2000, Aldowaisan and Allahverdi 2004).A detailed survey of the research and applications on this topic has been given by Nagano et al. (2016).
Most of literatures about no wait hybrid flow shop scheduling problem have mainly focused on a single objective.For example, Liu et al. (2003) suggests an algorithm known as Least Deviation (LD) in which the focus is on only one machine in each station.The key performance indicator in this point of view is the makespan.The performance is thus high using this algorithm.Also the algorithm is easier in computation and implementation.Having created such values, this is considered favourable.Xie et al. (2004) suggests a new heuristic algorithm named Minimum Deviation Algorithm (MDA) to minimize makespan in a similar method.MDA also performs better than partition method, partition method with LPT, Johnson's and modified Johnson's algorithms.Huang et al. (2009) considered a no-wait two stage flexible flow shop with setup times and with minimum total completion time key performance indicator.The author represents an enhanced programming as well as an Ant Colony approach.The solution was satisfactory with the approach and the results were efficient.
To best use of resources, Jolai et al. (2009) introduced no-wait flexible flow line scheduling problem with time windows and job rejection which is, in turn, an extension of production and delivery scheduling.He presented a similar method known as integer-linear programming model and genetic algorithm process as well.
In comparison with LINGO, studies show that the GA is a better solution in a computational time.Jolai et al. (2012) introduced a new hybrid algorithm with sequence-dependent setup times to minimize the total completion time.They suggest three algorithms.One, Population Based Simulated Annealing (PBSA), second, Adapted Imperialist Competitive Algorithm (AICA) and finally, hybridization of Adapted Imperialist Competitive and Population Based Simulated annealing (AICA+PBSA) for the problem.All the studies support the hybrid algorithm against the others which are applied in literature for related production scheduling problem.Rabiee et al. (2014) proposed the problem with respect to unrelated parallel machines, sequence-dependent setup times, probable reworks and different ready times to actualize the problem.What they proposed is based on imperialist competitive algorithm (ICA), simulated annealing (SA), variable neighborhood search (VNS) and genetic algorithm (GA) to solve the problem.The result revealed the advantages of our algorithm.To reduce makespan, Ramezani et al. (2013) suggested no-wait scheduling problem focused on set up time which is anticipatory and also sequence dependent in a flexible flow shop environment with two sets of same machines in parallel.They introduced it as a novel method since it was NP-hard.Their meta-heuristic method was about invasive weed optimization, adjustable neighborhood exploration and simulated strengthening to attack of the problem.The result showed, hybrid-metaheuristic outperformed in comparison with singular ones.
As said before, most of studies concentrated on single objective problem.But in reality, there is no single objective at all and we need to consider a batch of goals at once.However, there are also some studies working on multi objective in no wait flow shop scheduling problem.Allahverdi and Aldowisan (2004) suggested a method in no wait flow shop scheduling problem in with sum of makespan is important and also the maximum delay as a measure.Their methods for comparison were hybrid simulated annealing and a hybrid genetic heuristics.Also they suggested a dominance relation (DR) and a branch-and-bound algorithm.Herein also, after computation, the heuristic method performs higher and better in comparison with existing heuristics if the makespan and maximum lateness is considered significant.Also the dominance relation and branch and bound algorithm were totally effective.Aiming to minimize average of tardiness and time, Rahimi-Vahed et al. (2008) offered a bi-criteria no-wait flow shop scheduling problem.What they suggest is a new method named multi-objective scatter search as a metaheuristic algorithm for finding near optimal Pareto frontier.They were looking for effectiveness of this approach by solving some experimental problems in comparison with SPEA-II.Here also the better performance was appeared in multi-objective scatter search.
Multi-objective immune algorithm also proposed by Tavakkoli-Moghaddamet al. (2007) to minimize two goals as weighted average completion time and weighted average tardiness for a no-wait flow shop scheduling problem.He compared the algorithm with a conventional multi-objective genetic algorithm, i.e., SPEA-II.The generic algorithm won the game especially for significant and big problems.Pan et al. (2009) offered another algorithm called discrete differential evolution (DDE) for solving the no-wait flow shop scheduling problems with makespan and maximum tardiness measures.The results based on the famous benchmarks and statistical performance comparisons showed that DDE algorithm is much more efficient in comparison to the hybrid differential evolution (HDE) algorithm proposed by Qian et al. (2009).For the same objective, i.e. to minimize makespan and tardiness, Khalili (2012) proposed a multi-objective no-wait hybrid flow shop scheduling problem and suggested a novel Multi-Objective Electromagnetism Algorithm (MOEA) to solve the problem.He formulated the problem with mixed integer-linear programming models and proposed an effective Multi-Objective Electromagnetism Algorithm (MOEA) to reach the goal.
As far as multi-objective approach is considered, there are a few studies about no wait flexible flow shop problem.The author presented three multi-objective based algorithms to discover a no wait two stage flexible flow shop scheduling problem with a number of machines in each stage.The goals were minimizing makespan (i.e.C max ) and maximum tardiness (i.e.T max ).
In the next section; section 2, the multi objective terms are discussed.Then the author wrote about the bi-criteria no wait two stages flexible flow shop.Then in section 3, the multi-objective optimization search techniques are considered.Investigating the proficiency of suggested multi-objective metaheuristic approaches is presented next section‫‬ .And at last, the outcome of the research is induced and directions for further researches are depicted in section 5.

Multi-objective No-wait Two-stage Flexible Flow Shop Optimization
Now we turn to the problem statement.That is the concept of the multi objective optimization and also the structure of the problem.

Multi-objective Optimization
A multi-objective optimization problem formula is as below: 1 2 min ( ) { ( ), ( ), ..., ( )} . .( ) 0 , ( ) 0 s t g x h x  is determined by objective functions.In reality, the aim of a multi-objective optimization is to figure out the entire non-dominated solutions of the problem (any solution that is not able to develop an objective function with no effect on other objective).
If any of below conditions are met in an optimization problem with minimum objectives, solutions 1 x dominates solution 2 x .
(1) For every single objective (2) At least in one objective 1 () fx has a lower value compared with 2 () fx.
Having defined the dominant solutions, the optimal solution of a multi-objective optimization problem is defined as set of non-dominated solutions known as Pareto-optimal set which forms the Pareto front (CoelloCoello et al. 2002).

The Statement of Bi Criteria No-wait two Stage Flexible Flow Shop Problem
The no-wait two stage flexible flow shop scheduling problem is shown below: Given the processing time P ij of job j on stage i(i=1, 2), each of n job will be sequentially processed in stage 1, 2 respectively.At each stage there are m i machines.Also at a given time, each machine can process maximum a single job.Likewise, each action needs to be processed on one machine.Once the order of the action at the first stage is cleared the similar order is done for the second stage.To meet the no-wait boundaries, the end time of a job on a machine must be similar to the start time of the job on the next machine.This way, there is no elapse time in the entire operation.
The aim here is, reduction of makespan (C max ).The matching fitness function is considered at below: Next aim is reduction of maximum tardiness which is calculated below: T is the tardiness of job i, equal to max(0, Cd  , and i d is the due date of job i.

Possible Approaches to Find Solution
There are many classical approaches to solve multi-objective problems among them includes embracing goal programming, integer programming, e-constraint method and weighted sum method.The main features of the classical methods can be described as follows (Deb 2001, CoelloCoello et al. 2002): 1. Changing the problem from multi-objective to single objective 2. Experimentally, the methods might be applied randomly to find out the best solution 3. Each typical method includes some user-defined parameters that are not easy to set in an arbitrary problem.Some meta-heuristics have been developed to eliminate such deficiencies.They are genetic algorithms and evolutionary computation.The ability to figure out a reasonable estimate of Pareto frontier in one operation and good computational time, is one of the benefit of this method.
This paper offers three multi-objective metaheuristic methods to solve the problem; NSGA-II (Non-dominated Sorting Genetic Algorithm), MOICA (Multi-Objective Imperialist Competitive Algorithm) and PAES (Pareto Archive Evolutionary Strategy).They are to examine the output of the algorithms in solving the no -wait two stages flexible flow shop problem.
Metaheuristic algorithms are generally, based on a searching system which is random.Here, the problem altered from a phenotype into a genotype that is informally called chromosome.To discover the best solution, it uses intensification as well as diversification where the first intends to use local search area and the second explores the optimal solution globally.The chromosome, fitness evaluation, related operators and structures of applied system are described elaborated respectively:

Solution Representation and Fitness Evaluation
Some random values that are equals to the number of jobs in length is generated from 0 to 1.This is to show the chromosome.The jobs then are tossed by finding the increasing order of values in vector.See figure 1.

Chromosome
Job sequence Figure 1.An example of chromosome representation and its job sequence After calculation of jobs, procedure of machine assignments is done using a heuristic method.It means in each stage, to assign a job to the machines, the earliest available time machine is chosen and the job with the highest importance is assigned to that machine.If there are two important jobs at the same time, one is chosen randomly.The pseudo code of heuristic procedure is shown in Figure 2.

Crossover
The process in which tow chromosomes are prevented from coupling making progeny is called crossover.The process aims to find a better solution by mixing the chromosomes.The study uses a uniform crossover that first creates a random binary mask with the similar extend as the chromosomes and then substitutes relative gene material of parent chromosome based on created binary mask.This crossover results in a good exploitation of solution space (Syswerda 1989).

Mutation
After crossover, we perform an exchange mutation.For exchange mutation, two different arbitrary genes of the parent chromosome choose and swap the allele values.(Eiben and Smith. 2003).Probability of gaining Pareto-optimal solution using GAs are high, since they work with body of points.This makes it a strong tool for MOOPs too.The Non-dominated Sorting Genetic Algorithm (NSGA-II) is a well-known and extensively used algorithm based on its predecessor NSGA and proposed by Deb et al (2002).
Essentially, NSGA-II differs from non-dominated sorting Genetic Algorithm (NSGA) implementation in a number of ways.First, NSGA-II uses an elite-preserving mechanism, thereby assuring the preservation of previously found good solutions.Second, NSGA-II uses a fast non-dominated sorting procedure.Third, NSGA-II does not require any tunable parameter thereby making the algorithm independent of the user (Sivakumar et al, 2011).NSGA-II is a fast and very efficient Multi-objective evolutionary algorithm (MOEA), which incorporates the features of an elitist archive and a rule for adaptation assignment that takes into account both the rank and the distance of each solution regarding others.Salazar and Kishor have applied and compared the efficiency of NSGA-II with existing methods for reliability optimization problems (Kishor et al, 2008).
NSGA-II is an elitist multi-objective evolutionary algorithm which carries out an approximation of the Pareto front, based on the non-dominance concept.For achieving different Pareto fronts, a ranking procedure is performed at each generation.Also, this algorithm takes advantage of an operator called crowding operator for its diversification.NSGA-II starts from a randomly generated population of chromosomes (solutions), P_0 of size N.The population is sorted based on non-domination.Each solution is assigned a fitness (or rank) equal to its non-domination level (1 is the best level) the minimization of the fitness being assumed.A children population (Q_0) of size N is then created by applying the genetic operators: binary tournament selection, recombination, and mutation (Furtuna, et al, 2011).NSGA-II has been used in several prior studies like: Minella et al. (2008), Behnamian et al (2009), Zandieh and Karimi (20111) and Rabiee et al (2012), Asefi et al. (2014).
The crossover and mutation operators in this algorithm are as mention in NSGAII.The framewo rk of the proposed NSGAII is generally illustrated in the following pseudo code (Figure 3).

PAES Algorithm
The study suggests a simple multi-purpose metaheuristic algorithm called PAES.Offered by Knowles andCorne (1999, 2000), the algorithm applies a local search development strategy to a non-dominated solution in a pool of solutions that already applied.Hereunder, the most alternates of PAES which is (1+1) development strategy is discussed.
A solution is randomly generated, evaluated and saved in archived.Here is when the algorithm starts.Procedure at iteration t is continued with creating a new solution by transform current solution and compare it to current solution for dominance.The one with more dominance is accepted.If both solutions have same dominance priority, the new one is compared to archived solution which is archived.The accepted solution in the archive added there and removes the rejected one.Otherwise the new solution replaced by one of the archive member.At last, if there is no dominate member in the archive and archive has still vacancy, the new solution is added to archive.If there is no space in solution pool, the solution in the busiest area is removed and the new one is added.In when considering intricacy of problem to find non-dominated solutions, we consider a limitless archive scope to maintain non-dominated solutions.This is to obtain more Pareto solutions.The structure of proposed PAES is shown in Figure 4 , where pi is the variable to be optimized.Each variable in a country denotes a socio-political characteristic of a country.From this point of view, the algorithm searches for the best country that is the country with the best combination of socio-political characteristics, such as culture, language and economic policy (Atashpaz-Gargari and Lucas 2007).After country development, a non-dominance technique and a crowding distance are used to shape the fronts and rank member of each front, respectively.At that point, the members of front one are archived.Non-dominance technique and crowding distance described as below: 1: Non-dominance technique: imagine that there are r objective functions.When the following conditions are satisfied, the solution x1 dominates another solution x2.If x1 and x2 do not dominate each other, they are placed in the same front.
(1) For all the objective functions, solution x1 is not poorer than another solution x2.
(2) For at least one of the r objective functions x1 is exactly better than x2.
Solutions that are not dominated by others, constitute in front number 1.Meanwhile, the solutions that are only dominated by solutions in front number 1, organize front number 2. The same order applies to crate the other fronts which is shown in Figure 5.  Next, the prioritized solution is selected in the mass in terms of non-dominance and crowding distance.The selected solution is imperialists and the rests are colonies.
In order to compute the cost of prioritized solution (the imperialist), each target function is calculated.After that each target function is calculated: At last, the total value of each imperialist is calculated through: Where:  r is the quantity of target function After calculation the cost, the strength of each imperialist is obtained as well as the colonies distributed among the imperialist according to power of each imperialist country.
(10) At this point, the primary quantity of colonies is calculated as below: (11) Where:  n NC is the primary quantity of colonies of the n-th imperialist  col N is the number of all colonies n NC colonies are selected randomly and assigned to one imperialist.Apparently, the grater quantity of colonies, the stronger imperialist and the less quantity of colonies, the poorer imperialist.

Total Strength of an Empire
Imperialist country has the major impact on the total strength of an empire.But the strength of its colonies does not have such effect.Therefore, the equation of the total power of an empire is shown below.(Karimi et al. 2010, Shokrollahpour et al. 2011).
 n TP Emp is the total power of the nth empire  zeta (  ) is a positive number which is considered to be less than 1  Total cost of imperialists and colonies are calculated by Eq.8 and Eq.9.First: all of the imperialists and colonies are accumulated and then the non-dominated solutions are chosen.The percentage of the non-dominated solution belonging to each empire is calculated as n QE .
 is small and increasing it will increase the role of the colonies in determining the total power of an empire.

Moving the Colonies of an Empire toward the Imperialist (Assimilating)
Having distributed colonies among imperialists, the imperialist and relevant colonies go together.Figure 7 shows the movement.There, d is the distance between imperialist and colony.X is a random variable with a uniform (or any proper) distribution between 0 and d   and  is a number greater than 1.Direction of the movement is shown by , which is a uniform distribution between   and .A binary competition process is used to find out the best solution for both crossover and mutation operators in this way: Step 1: choose two solutions with the same size Step 2: the lowest front number needs to be chosen if both populations are from different fronts.
If both of them are in same front, choose the solution with highest crowding distance.

Revolution
Some revolutions have been done during last decades on the colonies.To do so, one of the famous policies, swap, reversion and insertion, exchange and local search is randomly chosen.The operators' structures are elaborated as follows: (Pris the revolution rate)  Swap: a colony's initial numbers are chosen randomly (numbers 1.46 & 2.27 in Figure 9) and their substitute each other.
Figure 9. Swap operator Reversion: in the current strategy, in addition to conducting substitution, the number that is placed between the substituted numbers are also changed.

Improve Imperialist
Using the rules describe earlier, the step involves in producing some neighborhoods around each imperialist considering its total power.The less imperialist power, the more neighborhood generation.Similarly, the more imperialist power, the less neighborhood generation.In other word, the quantity of neighborhoods directly depends on the power of the given imperialist.This is a leaner formula ranged between Ne min and Ne max .that is illustrated below: (13) Where:  best TP Emp is the value of total power the most power empire  worst TP Emp is the value of total power of the weakest power empire.
How the number of neighborhood and power are related are shown below: Figure 13.The relationship between power of each imperialist and number of neighborhoods

Colonies Updated
The primary mass of colonies, assimilating, information sharing among colonies, revolution and improve imperialist are combined all together in each decades to shape the empire that is called combined mass.Then, based on combined mass, the archive is updated.Then, for any of the imperialists, the best colony is chosen according to non-domination sorting and crowding distance by size of mass of colonies for a given empire (NC(i)).

Archive Adaption
For the combined mass, the classification is done using non-dominated and crowding distance.To archive, the front one members are chosen.At last, these members are retained and after classifying the solutions in archive, the other members are removed.Meanwhile, the size of archive equals n Archive.

Exchanging Positions of the Imperialist and a Colony
In this step, the total cost of each imperialist is updated.Next, the best imperialist and colony are combined.
Then, this mass is arranged by the non-dominated sorting and crowding distance.At last, the best mass is chosen as imperialist.The step is illustrated below.This term refers to a match among imperialist in which the weaker the territories, the more reduction of power would be and the more powerful the territories, the more power it gains.The competition is to take the hegemony of the weakest colony of the weakest territory.The competition is started by first choosing one or more colonies that are the weakest.Then the hegemony of these weak colonies is taken by a stronger territory through the competition.So far, this does not necessarily mean that the strongest territory is the winner.This means that those that are stronger keep more hegemony.This competition is modeled by just selecting one of the weakest colonies of the weakest territory to formulate the hegemony of each territory first is obtained the normalized total cost as follows.
(14) Where:  n NTP is the normalized total power of nth empire  n TP is the total power of nth empire After calculating normalized total power, the hegemony probability of each territory is obtained by: Next, to allocate the abovementioned colony to a territory, a so called roulette wheel method is used that is shown below: The process that the study is offered is concisely illustrated below: Figure 16.Flowchart of MOICA

Problem Design
The paper investigates the influence of some approaches for 36 test problems.Data are classified in three segments; quantity of operations, quantity of machines in first and second stages and, the dissemination of operation time in those stages.On the other hand, the abovementioned problems has been classified in two major segments; small and large problems.Table 1 depicts the quantity of operations and machine in small and large scale.Small:8,10,14,16,20,24 Large:72,80,88,108,120,132 No Round In addition the due dates are generated using the following formula: )

Parameter Setting
To evaluate the performance the suggested process, we need to set some key success factors.To come to this conclusion, some operations are simulated for both sizes of problems.Table 2 shows the tuned values of the proposed algorithms' values.

Performance Measures
The so-called Pareto based multi-objective optimization algorithm aims to find an estimate of non-dominated front.The key performance indicators of these methods are different from the single-purpose method.Hence, another method is needed to evaluate the key performance indicators of this algorithm.Then to evaluate the output of multi-purpose algorithms quantitatively, below evaluation methods are applied: Number of Pareto solutions (NPS): in this method, the number of non-dominated solutions resulted by an algorithm is computed.
Mean ideal distance (MID): The approximation between Pareto solutions and ideal point is defined.The calculation of MID is: Where n is the quantity of non-dominated solutions and . The lower the value, the better performance the MID has.
 Diversification metric (DM): This performance metric shows the range of solutions resulted by algorithms and is calculated as follows:  The spread of non-dominance solutions (SNS): Regarding to MID, this computes the range of non-dominated solutions.The SNS is defined as follows: (19)  % Domination: this key performance indicator applies a constructed Pareto combination set.Then the percentage of the solution belonging to each algorithm is calculated.

Experimental Results
This section refers to the result of experiments that are done through all the algorithms.The effectivene ss of each algorithm is presented and compared in terms of key performance indicators.All algorithms were coded using MATLAB 2013a and run on personal computer with a 2.66 GHz CPU and 4 GB main memory.
The efficiency of the algorithms was stated by solving 36 variant problems of which 18 are small and 18 are large in scale.The outputs of three algorithms regarding the five key performance indicator for both size are compared and shown in Table 3 to 7 respectively.
Relative Percentage Deviation (RPD) is applied for the best solutions in terms of the key performance indicator.The calculation is shown below: (20) Where Best is the best value between the algorithms Table 4, 5 and 6 show the output with 95% sureness for the percentage of domination, DM and MID key performance indicator for small size respectively.Deep analysis shows that MOICA beats the others in terms of domination percent.Concerning DM, there is no considerable difference between NSGAII and MOICA.Yet the MOICA still beats PAES.The presented facts disclose that for MID key performance indicator, the three algorithms are similar and there is no considerable variance.
Table 7 presents the output of large scale problems.Alike the small scale problems, MOICA beats NSGAII and PEAS.The more size of the problem, the more considerable the advantages are.

Conclusion and Further Researches
The study presents two criteria flow shop that is called two stage flexible flow shop.It aims to shorten the makespan and increase lateness of operations.The paper modeled Pareto optimal solutions and based the process on its similarity.NSGAII, MOICA and PEAS were suggested as the main three metaheuristic Pareto based multi-purpose algorithms.To estimate the efficiency of them, 36 problems, large and small, were answered.In multi-purpose norms, five key performance indicators, NPS, MID, DM, SNS and percentage domination were suggested to disclose the algorithms' efficiency.With similar efficiency, NSGAII beats PAES in both small and large scale of problem, in terms of DM. hence, MOICA is the best algorithm in case of efficiency for all the studied problems.
To Guidance for researchers in similar cases, using other effective metaheuristic algorithm like multi -purpose and colony optimization or multi-purpose invasive weed optimization are suggested to work on.Also there would be valuable result if the stages increase for more than two.
shows the possible solution in n x dimensional search space and () fx is a M dimensional vector of objective values.Map between decision variables of

Figure 2 .
Figure 2. Pseudo code of simulator

Figure 3 .
Figure 3. Proposed NSGA-II algorithm in pseudo code

Figure 5 .
Figure 5. Non-dominance technique 2. Crowding distance: this is a tool to show the quantity of solutions in each step.See below figure.This is an estimate of the solution mass around a given solution.

Figure 6 .
Figure 6.Crowding distanceThe crowding distance measure which is used in MOICA is shown in equation (7).The solutions having a lower value of the crowding distance are prioritized over solutions with upper value of the crowding distance.
, maximum and minimum values of objective function i in each iteration, respectively.
 n QE is the quality of empire nth n QE is determined as below:

Figure 7 .
Figure 7. Moving colonies toward the imperialist with a random angle 3.3.3.4Information Sharing between ColoniesTo improve their position, the colonies share their information.To do so, in this part, one of these operators including one-point, two-point and continuous uniform crossover shown in Figure8are selected randomly.The mass part that are sharing information are shown by Pc.Those colonies with stronger position have more opportunity to share their information since the selection here is by a competition that is describe below:

Figure
Figure 10.Reversion Insertion: for the insertion strategy, like substitution, two numbers of a colony is chosen by chance (Numbers 1.46 & 2.27 in Figure11).After that, the next number is places approximate to the number in the first position.The other number is moved right side consequently.

Figure
Figure 11.Insertion Perturbation: in perturbation strategy, one number is chosen by chance, and another number is generated by chance.Then these two are substituted.

Figure 12 .
Figure 12.Perturbation Local search: the strategy refers to randomly selection of the iteration.After that all of the two -point exchanges are investigated.

Figure 14 .
Figure 14.Exchanging positions of the imperialist and a colony

FigureFigure 15 .
Figure 15.Imperialistic competition 3.3.3.11Eliminating the Powerless Empires Through the competition, weak territories will ruin and their members dispense among other territories.The study refers to this ruined territory as collapses as shown below:
Every single solution in the imperialist competitive algorithm simulates as an array.The arrays include different values that need to be adjusted.What is called chromosome in GA terminology, is named country here.In an N-dimensional optimization problem, a country is a 1×N array.This array is defined by:

Table 1 .
Factors and their levels

Table 2 .
Tuned values of the parameters of the algorithms

Table 8 ,
9 and 10 also shows the outputs of algorithms in terms of domination percent in which DM and MID are shown consequently.By a glance at the table 8, MOICA's superiority in terms of performance in comparison to the next two.Outputs of DM are also illustrated in table 9 and disclose that the performance of NSGAII and MOICA identical and both of them outperform PAES.Moreover, when MID is considered, MOICA beats NSGAII and PEAS as shown in table 10.

Table 3 .
The simulation results for small size problems

Table 4 .
95% confidence interval for % domination in small size problems

Table 4 .
95% confidence interval for % domination in small size problems Individual 95% CIs For Mean Based on