ALGORITHMS FOR BICRITERIA MINIMIZATION IN THE PERMUTATION FLOW SHOP SCHEDULING PROBLEM

. This paper presents two bi-objective simulated annealing procedures to deal with the classical permutation ﬂow shop scheduling problem considering the makespan and the total completion time as criteria. The proposed methods are based on multi-objective simulated annealing techniques combined with constructive and heuristic algorithms. A computational exper-iment has been carried out and diﬀerent metrics have been computed to check various attributes of each method. For all the tested instances a net set of potentially eﬃcient schedules has been obtained and compared with previously published results. Results indicate that the proposed algorithms provide ef-ﬁcient solutions with little computational eﬀort which can serve as input for interactive procedures.

1. Introduction. We consider the minimization of the makespan and the total completion time for the classic permutation flow shop scheduling problem, and present two algorithms based on a Multi-Objective Simulated Annealing (MOSA) scheme (named as MOSAD and MOSAI) to deal with it. This scheduling problem presents a relatively simple formulation. However it belongs to the set of hard Combinatorial Optimization problems (CO), even considering only one criterion. Moreover, quality is, in real-life, a multidimensional notion. Thus, we have to add the complexity that comes from the multivariant condition of corresponding alternative schedules. A solution which is optimal with respect to a given criterion might be a poor candidate for another. Of course, for multicriteria scheduling problems expecting to find the optimum schedule must usually be discarded. We would be satisfied with finding a set of Pareto optimal alternatives, and we have to let some subjective considerations intervene, such as the decision-maker preferences. It is actually a Multicriterion Decision Making problem (MDM), and at the present time, there is no other rational tool to apply to discard alternatives.
Until the 1990s, most research in the area of multiple criteria scheduling consisted of bi-criteria studies of the single machine case [16]. Only with the breakthrough of metaheursitcs in solving CO problems did researchers begin to adapt metaheuristics to solve Multi-Objective Combinatorial Optimization (MOCO) problems.
MOSA methods are metaheuristics based on Simulated Annealing (SA) to tackle MOCO problems. SA, introduced by Kirkpatrick et al. [23], has demonstrated its ability in solving MOCO problems [42]. Based on a MOSA scheme, we have developed our algorithms to provide the decision-maker with approximate efficient solutions for the permutation flow shop problem, considering the makespan and completion time as objective functions.
Starting with a set of initial sequences, the algorithm samples a point in its neighbourhood. If this generated sequence is dominated, we still accept it with a certain probability. Different constructive and heuristic algorithms are used to compute initial good sequences and lower bounds for the different criteria.
We carried out an intensive computational experiment by making use of the 90 benchmark problems given by Taillard [49]. The procedure is good enough to give efficient solutions with little computational effort. The performance analysis includes a set of metrics specific for evaluating Multi-Objective Optimization algorithms (MOO). We also report a net set of potentially efficient schedules that updates previous published net sets, for the same instances ( [35] and [50]). The presented algorithms have been developed after an extensive study involving the influence of the number of initial solutions, the neighborhood search procedure and the SA parameters, together with the CPU time.
This paper is organized as follows: In the next section the problem statement is presented, and its computational complexity and previous results are looked at. Section 3 describes the new approaches based on the MOSA scheme. Section 4 reports on the computational experiment. We conclude, in section 5, with a summary discussion on further research.
2. Problem statement and review.
2.1. Permutation flow shop scheduling problem. The flow shop problem results in the context where a set of multi-purpose machines are used to manufacture different jobs. So, the aim is to find the schedule that optimizes certain performance measures. In a general flow shop scheduling problem, there are n jobs and m machines or stages. Each job needs to complete one operation on each machine during a fixed processing time. The algorithms presented in this paper are devoted to the permutation flow shop context, where all jobs must pass through all machines in the same order. Thus, the scheduling process involves just finding the sequence of jobs that optimizes criteria. But, the computational complexity grows exponentially, making the problem intractable even if only one criterion should be considered.
2.1.1. Notation. We will use the notation that follows: J: set of n jobs J i (i = 1, . . . , n) M : set of m machines M j (j = 1, . . . , m) p ij : processing time of job J i on machine M j r i : time at which the job J i is ready to be processed w i : priority or weight of job J i C i : completion time of job J i C max : the maximum completion time of all jobs J i . This is the schedule length, which is also called the makespan. F i : flow time of job J i , F i = C i − r i The optimal value of any criterion is denoted with an asterisk, e.g. C * max denotes the optimal makespan value.
We will use the three-parameter notation, α/β/γ , introduced by [13] and extended for [45] to multi-criteria scheduling problems. The first field specifies the machine environment (F represents general permutation flow shop); the second, job characteristics; and the third refers to the chosen optimality criterion for single criteria models, and it extends to cover multicriteria as well as methodology.
2.1.2. Definitions and Assumptions. Consider a set of n independent jobs J i (i = 1, . . . , n) to be processed, each of them on a set of m machines M j (j = 1, . . . , m), that represent m stages of the production process. Every job requires a known, deterministic and non-negative p ij for completion at each machine. Because each machine processes the jobs in the same order, any permutation of the jobs is then a feasible solution for the scheduling problem, regardless the optimization criteria.
The most used criterion is the minimization of the total completion time of the schedule, C max . But there are many performance criteria to be considered when solving scheduling problems. Flow time is applied as a criterion when the cost function is related to the job standing time. Completion time reflects a criterion where the cost depends on the finish time. In the event of all ready times being zero, r i = 0, ∀i, completion time and flow time functions are identical. When some jobs are more important than others, weighted sum measures could be considered. Maximum criteria should be used when interest is focused on the whole system.
We consider C max and i C i simoultaneously, thus, the aim is to find out the set of feasible schedules that are Pareto optimal solutions. Specifically, the makespan criterion is equivalent to the mean utilization performance measure. In contrast, the total flow time criterion is a work in process inventory performance measure. While the makespan criterion is a firm-oriented performance measure, the total flow time is a customer-oriented performance measure. Therefore, the set of efficient solutions to a bi-criteria model that seeks to optimize both measures simultaneously would contain valuable trade-off information, crucial to the decision-maker, who has to identify the most preferable solution, according to her preferences.
Unless explicitly indicated, in the following we assume that: 1. Each job is an entity, composed of m operations, which cannot be processed on more than one machine simultaneously. 2. At every machine, there are no precedence constraints among operations of different jobs. 3. No preemption is allowed. That is to say, once an operation has started, it must be completed before another operation may initiate on the same machine. 4. No cancellation. Each job must be finished. 5. Processing times are independent of sequencing. 6. Job accumulation is allowed. Jobs can be waiting for a free machine. 7. Machine idle time is allowed. Machines can be waiting for jobs or for the end of the total process. 8. No machine can process more than one job simultaneously. 9. Machines never break down and are available throughout the scheduling period. 10. Ready times are zero for all jobs. 11. There is no randomness: (a) the number of jobs, n, is known and fixed; (b) the number of machines, m, is known and fixed; (c) the processing times, p ij (i = 1, . . . , n and j = 1, . . . , m), are known and fixed; (d) all other specifications, needed to define a particular problem, are known and fixed.

2.2.
Computational complexity and state of the art review. The permutation flow shop problem, like almost all deterministic scheduling problems, belongs to the wide class of CO problems, many of which are known to be NP-hard [11].
Thus, heuristic and hyperheuristic methods have been developed, some of them showing an acceptable performance.
Since the early Johnson algorithm [21] that solves F 2//C max in polynomial time, only a few restricted cases have been shown to be efficiently solvable. Minimizing the sum of completion times is still NP-complete for two machines [11]. Real-life scheduling problems require more than one criterion. Nevertheless, the complex nature of flow shop scheduling has prevented the development of models with multiple criteria. Permutation flow shop scheduling research has been mostly restricted to the treatment of one objective at a time. Furthermore, attention was focused on the makespan criterion. However, the total flow time performance measure has also received some attention. These two regular criteria constitute a conflicting pair of objectives [39].

2.2.1.
Multi-Objective Flow Shop Scheduling Problem. The multicriterion nature of the problem at hand leads us to the area of MOO (see [8], for a state of the art review). Due to the NP-hardness of the most scheduling problems just considering only one regular criterion (except for restricted special cases), research work concerning them have been intensively developed. The proliferation of metaheuristic techniques has encouraged researchers to apply them to this highly complex context.
The MOCO problem we are dealing with is a discrete optimization problem which can formulate as follows: min where X is the vector that represents a feasible solution (a sequence), and D is the set of feasible solutions (a discrete set). The criteria are of two different kinds: • Bottleneck function: z 1 = C max = max {C 1 , C 2 , . . . , C n } • Sum function: z 2 = n i=1 C i We call a feasible solution, X (e) ∈ D, efficient, non-dominated, or Pareto optimal, if there is no other feasible solution X ∈ D such that, with at least one strict inequality. The corresponding vector of objective values, is called non-dominated vector.
The set of feasible efficient solutions is denoted by E, and the set of nondominated vectors by N D.
But the discrete nature of the feasible set D, makes the programm non convex. The transformation of the objective function into a linear one (aggregating into weighted sums) does not transform the problem into a convex programm. Except in some special cases, e.g. preemption allowance, or where idle time insertion is advantageous, for which linear programming can be applied, the discrete structure of a MOCO problem persists. An important consequence is the fact that there could be some efficient solutions not optimal for any weighted sum of the objectives. The set of these solutions are named non-supported efficient (N E), whereas the set of the remaining ones are called supported efficient (SE). Concerning computational complexity in obtaining the set E, MOCO problems are in general NP-complete. N E are the main reason why the computation of E is hard [8]. The cardinality of the N E set depends on the number of sum objective functions. For a problem with more than one sum objective function, N E has many more solutions than SE [8]. Despite these results which constitute the essence of the difficulty of MOCO problems, many published works ignore the existence of N E.
The cardinality of E for a MOCO problem may be exponential in the problem size [9], therefore algorithms could determine just an approximation of E in many cases. Thus, methods may be exact or approximate, and metaheuristics are nowadays being applied intensively to MOCO problems.
In the following we will focus on multicriterion flow shop problems. We refer the reader to: [8] for further information on MOCO theory; [25] and [22] for overviews on metaheuristic applications; and [20] for metaheuristics on bi-criteria optimization problems.

2.2.2.
MultiCriteria flow shop scheduling problem review. [36], [33], [14] and [46], present heuristic procedures for the two-machine flow shop problem, where total flow time has to be minimized among the schedules that minimize makespan (lexicographical approach). Local search algorithms based on Ant Colony Optimization (ACO) have been proposed by [47]. [17] presents a technique named Local Dynamic Programming.
[46] presents a Branch and Bound (B&B) algorithm, which can solve problem instances of up to 35 jobs to optimality. [31] and [43] present heuristic and B&B algorithms for the F 2//f ( C i , C max ) problem. For the two-machine flow shop scheduling problem of minimizing makespan and sum of completion times simultaneously, [40] presents an a posteriori approach based on B&B. [6] presents a B&B algorithm for the F 2//f (C max , T max ) , and a heuristic to approximate the set of non-dominated solutions for the more general F//f (C max , T max ) problem. [26] presents B&B algorithms for the F 2//f (C max , U i ) and F 2//f (C max , T i ) problems.
[41] and [51] consider a linear combination of makespan and flow time. [7] presents a comparison of four iterative improvement techniques for flow shop scheduling problems that differ in local search methodology. These techniques are iterative deepening, random search, Tabu Search and Genetic Algorithm (GA). The evaluation function is defined according to the gradual satisfaction of explicitly represented domain constraints and optimization functions. The problem is constrained by a greater variety of antagonistic criteria that are partly contradictory. [15] and [38] propose heuristic procedures for the general m machine case, considering C i , C max , I j . They are based on the idea of minimizing the gaps between the completion times of jobs on adjacent machines (one of our proposed improvement procedures was inspired by this paper). [52] applies ACO to the same problem.
[2] presents a Multi-Objective Genetic Algorithm (MOGA) that improves the previous MOGA presented by [44]. [30] presents a MOGA considering the makespan, total flow time and total tardiness, based on a weighted sum of objective functions with variable weights. This algorithm belongs to the class of Evolutionary Multi-Objective (EMO) algorithms and [19] shows that this algorithm can be improved by adding a local search procedure to the offspring. [3] applies subpopulation GA to the same problems. Artificial chromosomes are created and introduced into the evolution process to improve the efficiency and the quality of the solution. [1] proposes a MOGA algorithm with preservation of dispersion in the population, elitism, and use of a parallel bi-objective local search so as to intensify the search in distinct regions. The algorithm is applied to the makespan-maximum tardiness and makespan-total tardiness problems. [10] investigates both a priori and a posteriori heuristics for makespan and flowime. The last one uncovers non-dominated solutions by varying the weight criteria in an effective way. [4] proposed a GA algorithm for the F//f ( C i , C max ) problem based on the concept of gradual priority weighting (the search process starts in the direction of the first selected objective function, and progresses such that the weight for the first objective function decreases gradually, and the weight for the second objective function increases gradually). [50] applies a similar idea. [35] presents a Pareto GA with local search, based on ranks that are computed by means of crowding distances. [35] and [50] apply the same initial population and improvement schemes.
In [28] an exhaustive review of the multi-objective flow shop scheduling problem can be found. [12] presents a study of the problem structure and the effectiveness of local search neighborhoods within EMO flow shop scheduling problems.
3. MOSAD and MOSAI algorithms. We present two variants of an approximation algorithm to get the Pareto solution set of the bi-criteria permutation flow shop problem, considering the minimization of the makespan and the total completion time simultaneously. The most promising practical approach to general MOCO problems consists in generating efficient solutions with meta and hyper heuristic procedures. In this work we make use of SA techniques.
SA is a generic technique (based on an analogy to physical cooling studied by statistical mechanics), and has to be adapted in the context of the specific problem being studied [23]. It is basically an improvement technique, by which an initial solution is improved by means of local perturbations. MOSA methods are characterized by: • An acceptance rule for new solutions, with some probability that depends on the temperature level. • A scheme of cooling.
• A mechanism for browsing the efficient frontier.
• Information is obtained from the set of solutions.
[42] presents a broad study of the application of SA to MOO.
The proposed algorithms mainly consist in a MOSA scheme which begins with an initial iterate solution, X 0 , that belongs to a set of initial points (feasible solutions of F//{C max , C i }, which are good solutions for one of the two simplified single-criterion versions of the problem). X 0 is then sampled with a point Y in its neighbourhood. But instead of accepting Y , if it is better than the current iterate regarding an aggregated objective function, we accept it, if it is not dominated by the current solution. In this case, we make Y the current iterate, add it to the Potentially Efficient solution set (PE ), and throw out any point in PE that is dominated by Y . On the other hand, if Y is dominated by X 0 , we still make it the current iterate with some probability. This randomization is introduced in the procedure to attempt to reduce the probability of getting stuck in a poor locally optimal solution.
The sampling procedure is iterated according to the MOSA scheme described later. The solutions that are generated, during the optimization process, make iterative updates to the PE point set, to get closer to the Pareto optimal set, E. The only complicated aspect of this algorithm is the necessity of generating solutions in several directions of the bi-objective space search, (C max , C i ) . So, to be able to cover the entire efficient frontier, a diversified set of points must be generated. Improvement and perturbation procedures play a crucial role in the performance of the algorithm. Both of them are carried out by simple neighbourhood exploration heuristics. The objective of these procedures is to approximate the trade-off surface in a more efficient way by using those movements that are more promising according to the current solution. These procedures induce a privileged direction of search to the efficient frontier. So, to be able to cover the entire efficient frontier, a diversified combination of initial solutions and neighbouring generation heuristics must be considered.
Depending on the relation between the optimization criterion for which the iterate solution presents the least deviation (which coincides with the criterion for which the corresponding initial solution has been computed) and the criterion taken into account for the improvement technique, we have differentiated two movement strategies which give rise to the two variants of the proposed method: MOSAD, where the last letter denotes direct: The improvement mechanism search in the same direction of the criterion for which the initial seed of X q has been calculated. Considering that the criterion used in the generation of the initial solution is makespan/total completion time, then, the preferable criterion for improvement mechanism is makespan/total completion time.
MOSAI, where the last letter denotes inverse: The improvement mechanism search in the direction of the criterion for which the initial seed of X q has not been calculated. Considering that the criterion used in the generation of the initial solution is makespan/total completion time, then the improvement technique looks for solutions with less total completion time/makespan.
To avoid the drawback mentioned in section 2.2.1, we do not use any scalarized objective. The objective function values of the MOCO problem play here the role of acceptance rule. In order to avoid the NE solutions to enable entrance into the PE, we have developed a bi-objective model where, simoultaneously, both criteria are included in a domination control. For just bi-criteria models, checking whether a solution is dominated by another is not computationally costly, and besides, updating the non-dominated solution set have to be tackled in any case (for more than two criteria models, the use of aggregated functions may be justified).
3.1. Pseudocode of the proposed algorithms. MOSAD and MOSAI are described with the following pseudocode making use of a variable named ALGO. ALGO=Direct corresponds to MOSAD, while ALGO=Inverse corresponds to MOSAI.
Else return to step 3d. In the following, procedures, specifications and parameters are described in detail.

3.2.
Construction of initial solutions. A set of four initial solutions is built starting from two constructive algorithms devoted to the F//C max problem, and two constructive algorithms devoted to the F// C i problem. The robustness of this set is crucial in reducing the search effort. In previous works [29] we have developed a mechanism of generating a variable number of initial solutions. However, after computational experiments, we have verified that the best trade off (between computational effort and efficiency of solutions) corresponds to four initial solutions.
The first two solutions are obtained by means of two simple but effective constructive algorithms: X 1 : : is built by the improvement by Taillard [48] over NEH [32], looking for the minimum makespan, NEH (C max ), in the following. X 2 : : is built by the algorithm for the F// C i problem presented by Rajendran [37], Raj ( C i ), in the following. While X 1 is considered as a good solution for C max , X 2 is considered as a good solution for C i . X 3 and X 4 are obtained by two list scheduling procedures similar to NEH (C max ) and Raj ( C i ), respectively. While in NEH (C max ) jobs are listed according to descending order of . Re-combining: making list (ordering by p i or w i ), and criterion to be minimized (C max or C i ) we have the following two algorithms, that get X 3 and X 4 . NEH( C i ): Step 1: For each job i, calculate the total processing time p i = m j=1 p ij .
Step 2: List the jobs according to descending order of p i .
Step 3: Schedule the first two jobs (from the list) in order to minimize the partial total completion time (as if there were only these two jobs).
Step 4: For k = 3 to n, insert the job k at the position which minimizes the partial total completion time, among the k possible places. So, X 3 is considered a good solution for the total completion time criterion. By analogy, Rajendran algorithm is applied to obtain X 4 , a good solution for the makespan criterion as follows.
Raj (C max ): Step 2: List the jobs according to ascending order of w i .
Step 3: Schedule the first two jobs (from the list) in order to minimize the partial makespan (as if there were only these two jobs).
Step 4: For k = 3 to n, insert the job k at the position which minimizes the partial makespan, among the k possible places. All of these four generation algorithms share the same four-step structure. In the computation of the initial solutions, the procedure keeps the useful data in order to save computational effort (e.g. job lists, best partial schedules, etc.).
3.3. Improvement scheme. Instead of inducing the search direction by tuning weights, to improve the distribution of non-dominated solutions in PE, we apply different neighbouring search heuristics. In the MOSA scheme described previously, a neighbouring solution of the current permutation must be chosen. The most important neighbourhoods based on a single permutation as an input are: • Exchange, swapping the positions of two jobs at i and k, with i = k. The remaining jobs in the sequence conserve their positions. • Insertion, forward or backward shift, removing a job at i and reinserting it at a different position k, with k > i in forward case, and with k < i in backward case. The remaining jobs in the sequence must be re-arranged in order to keep their relative positions. While insertion has been shown to lead to superior results compared with exchange, for flow shop scheduling problems with C max objective [48], it does not seem to be possible to derive a similar general rule when considering C i . So, for improving makespan, we have only implemented insertion. On the other hand, for completion time, we try with insertion and exchange.
In the valuation of a neighbour, it is very important to save computational effort in order to check a larger neighbourhood. When inserting or exchanging jobs in a schedule, it should be possible to discard some dominated candidate permutations with just small computational requirements examining the corresponding partial schedules. With this in mind, we have developed two neighbouring generation heuristics: one devoted to improve the current solution regarding makespan, another one sampling better solutions regarding total completion time. Thus, depending on the generation procedure employed for the seed, the neighbouring heuristic is chosen. The first one is applied to the seeds X 1 and X 4 (which have been built by minimizing makespan), while the second one is applied to the seeds X 2 and X 3 (which have been built by minimizing total completion time).
3.3.1. Improving Makespan. In order to reduce the search space, we have developed a technique based on elimination by domination conditions. Furthermore, we compute the lower bound for the makespan introduced by [49]. Thus, if we find out a permutation having this makespan value, we stop searching on the C max axis, and concentrate effort in exploring the direction of reducing C i , in the neighborhood of C * max . The lower bound is computed as: NEH (C max ) is actually a procedure to compute the value of the partial makespan when a job i is added in a partial schedule at position k. We employ this algorithm embedded in our neighbourhood search heuristic as a shortcut to evaluate a partial permutation. So, we do not need to compute the objective function for the complete schedule. Based on domination criteria for partial schedules (see [18] and [27]), we have developed our elimination neighbouring search.
Any partial schedule of t jobs, Jp(t) = {J 1 , J 2 , . . . , J t }, where t = 1, 2, . . . , n, is a sequence of the indices corresponding to the jobs in Jp(t), and it could be named as σI Jp(t) . The completion time for a partial schedule σI Jp(t) on machine k, where k = 1, 2, . . . , m, is denoted by C(σI Jp(t) , k). It was proved that: For the case where σII Jp(t) and σI Jp(t−1) are partial schedules of Both theorems allow us to discard any completion of a partial schedule σI Jp(t) or σI Jp(t−1) , because a schedule at least as good exists among the completion of another partial schedule σII Jp(t) .
The improvement scheme proposed in this section is based on the sequential insertion of a job in the current sequence at each possible different position. Since jobs with larger total processing time at the beginning of the schedule bring, in general, schedules with less makespan value, the proposed scheme selects, for insertion, a subset of jobs which are located at the first β% of the total positions in the current sequence. Hence, the set of t jobs scheduled at 1, 2, . . . t positions, in the current permutation, σ Xq , where t = β%n, is selected for exploration consisting in checking whenever a better partial permutation, involving these t jobs, could be built.
Theoretically, we have to check and compare, for each job placed at i on σ Xq , with i = 1 to t, the makespan that results when this job is placed at a different position j, with j = 1 to t. Nevertheless, the elimination criteria described above leads to efficiency gains. The Elimination Criterion 2 will filter any potential permutation generated by moving a job for which, to be placed at a different position with respect to its position in the current schedule, will not yield a sequence with less total completion time. Only for a potential permutation that passes this control, specified for a job to be moved, we check then for different positions. By the Elimination Criterion 1, any potential schedule σ Y , for which the current schedule σ Xq is at least as good, will be discarded. If one partial schedule is not eliminated, then the corresponding complete schedule, σ Y , becomes the generated neighbouring solution, σ Xq = σ Y , and the lower bounds used for computations are updated.
3.3.2. Improving Total Completion time. Following the ideas of [15] and [38], we have developed an improvement heuristic looking for permutations with lower total completion time values, but attempting not to lose the level obtained in makespan. The original idea was to minimize gaps between successive operations that would lead to a better quality solution. The pair of jobs with the most positive gaps has to be placed at the beginning, while the pair of jobs with the most negative gaps has to be placed at the end of the schedule. During the total processing of the whole set of jobs, the larger gaps would have more chance of being compensated with the negative gaps corresponding to the pair of jobs scheduled at the end of the sequence.
The improvement scheme proposed in this section is based on the interchange of adjacent pairs of jobs with positive gaps in the current permutation. Since the objective is to minimize gaps, the jobs are listed in descending order of gaps. Exchanging adjacent jobs with larger gaps is more likely to result in a permutation which yields less completion time value. The procedure selects, for exchanging, a subset of jobs, Jp(t), which are located at the later β% of the total positions in the current sequence σ Xq . The exploration consists in checking whenever a new permutation obtained by exchanging an adjacent pair of jobs of Jp(t), yields a schedule with lower completion time value.
Step 1: The subset of jobs placed at the last t positions of σ Xq , where t = β%n, is selected to constitute the set Jp(t).
Step 2: For the jobs of Jp(t), the gaps between every pair of adjacent jobs in σ Xq , is then computed as Step 3: Jobs in Jp(t) are listed in descending order of gaps G i . Jobs with negative gaps are not included in the list, and ties are broken in descending order of this similar gap: (7) that is computed only in the case of a tie.
Step 4: The first job in the list, Jσ Xq (i), scheduled at position i in the current permutation σ Xq , will be set at i + 1, in a new permutation σII, and its counterpart, the job placed at i + 1 in σ Xq , will be set at position i in σII. Step , then σII is accepted as a new permutation, σ Xq = σII, then return to Step 1. Otherwise proceed to Step 6.
Step 6: Remove Jσ Xq (i) from the list of jobs. If the list is not exhausted, then return to Step 4.

3.4.
Updating potential efficient set. When sampling solutions, only those which pass the domination control are taken into account for listing in the PE solutions set. The rest of the generated solutions are only used as input (for improvement or neighbouring search) and discarded later. When a neighbouring solution Y is accepted and made the current solution, X q = Y , the set of PE solutions should be updated. If X q is a new non-dominated solution, it should be added to the archive set and the set should be updated. Any solution dominated will be removed from the set. For the efficiency of this algorithm the updating PE process is crucial. In order to save computational effort in updating PE, non-dominated solutions are always stored in ascending order of one of the criterion values, thus their other criterion values will be stored in descending order. This arrangement constitutes a fast method of finding out dominated instances with respect to a new solution, and of updating the PE.

3.5.
Perturbation of X k . Starting with the improved solutions, this procedure aims to help in approximating the trade-off surface in a more efficient way. The common procedure explores the neighbourhood of the current solution, X q , and yields a neighbour Y , following the direct search strategy (MOSAD, ALGO=Direct) or the inverse search strategy (MOSAI, ALGO=Inverse), as follows: Let z k be the preferable criterion according to the chosen strategy. For the permutation σ Xq = {I 1 , I 2 , . . . , I n }, where I i , with i = 1, 2, . . . , n, is the index of the job scheduled at the position i in σ Xq , we will check its neighbourhood for finding out a better neighbour with respect to z k . As it is known that insertion brings better improvement than exchange, this procedure generates potential permutations by inserting, forward and backward, removing each job J Ii , and reinserting it in a different position at random. For each job, J Ii , where I i = I 1 , and I i = I n , the procedure will choose randomly two positions for insertion. One position to its right, choosing randomly a number between i + 1 and n, for forward insertion, and another position to its left, choosing randomly a number between 1 and i − 1, for backward insertion. For jobs in extreme positions, I 1 and I n , only one direction of insertion can be chosen. For I 1 only forward insertion is possible to apply, hence, to select a new position, a random number between 2 and n must be generated. In a similar way, I n can only be inserted at positions to its left, so a number between 1 and n−1 has to be chosen. Thus, the z k value of the 2(n−1) potential permutations has to be evaluated and the permutation with minimum z k , is then selected as the neighbouring solution Y .
This mechanism has been inspired on the rule known as Random Insertion Perturbation Scheme, introduced by [34]. The development of MOSAI was motivated for the observation (in previous studies, see [29]) of a tendency of solutions to be concentrated near the extreme points of the efficient frontier. MOSAI intensifies the random generation of seeds, filling in the search space with more central solutions (see Figure 1).
3.6. Simulated annealing parameters. Any SA method requires some specification. Besides the determination of the state space (z 1 , z 2 ), the following parameters have to be specified: the fitness function φ, the acceptance probability function p, the initial and final temperature, T 0 and T f , the temperature reduction factor α, the length of the temperature step, N step , a stopping criterion, etc. To set these appropriate parameters is an important issue of the SA design because of their significant impact on the method's effectiveness. Unfortunately, there is not a set of parameters good for all instances of a problem, and there is no general way to find the best setting for a given problem. After an exhaustive study carried out by varying SA parameters, we have determined the appropriate setting values. Here we point out some aspects of them.
p is the probability, for a dominated alternative solution, of being accepted, and it is computed on the basis of its fitness value and the temperature of the cooling process. In the SA technique, temperature is reduced at every step of iterations. This cooling process decreases the possibility of admitting a dominated solution during the search process. A higher probability, at the beginning of the process, helps in attempting to avoid being trapped in a local optimum. The fitness of a solution is evaluated by means of the normalized deviation function This normalization diminishes the influence of the different dimensions of each criterion, hence we have a dimensionless quantity which indicates the relative deviation of the quality of the generated solution, Y , from that of the current one, X q . Since ∆φ is not dependent upon the instance size, the initial and final temperature values can be fixed more reasonably and accurately to minimize the computational effort without sacrificing the quality of the final solution. Initial temperature should permit the acceptance of inferior quality solutions. After the study including different bound values, we have chosen that the algorithm starts with a temperature value of 475, and finishes when the temperature is below 5 (shorter intervals give worse PE, while longer ones increase the computational time without noticeable improvement in results). These values are set to limit the inferior quality of acceptance of a generated permutation by 50%. This means that the probability of accepting a solution with deviation of performance criteria of 50% is 0.9 at the beginning of the iterations and less than 0.0001 at the later iterations. The temperature will be reduced by the factor α = 0.968. This reduction takes place at every length of step iterations (with or without improvements). Thus, the temperature will be at 140 different level values (T f = α 140 T 0 ). In order to control the computational effort, a stopping criterion must be fixed, thus the number of iterations without improvement, N stop has been fixed. Furthermore, the length of the temperature step, N step , is essential in driving the cooling process. After the aforementioned analysis, we have fixed the N step in 300 (for lower values the algorithms give less efficient solutions, and for higher values the distribution of the solutions is less uniform), and for N stop we have considered 2500 iterations.

4.1.
Metrics. The goal of the proposed MOSA algorithms is to find a good approximation of the set E. For MOO algorithms it is unlikely that the whole set of efficient solutions is fully known (except for small size instances, with non-practical application). While the outcomes from compared algorithms are different, they can still be all equally Pareto efficient. Thus, their performance analysis is quite different from the Optimization techniques for single criterion, and the three following conditions are usually considered as desirable for a good MOO algorithm: 1. The distance between the obtained PE solutions and E should be minimized.
2. The distribution of the solutions in PE should be uniform.
3. The larger the number of obtained solutions, i.e. the cardinality of PE, the better the algorithm. The last two conditions present more weaknesses than strengths. If E does not present a uniform distribution, or [E] = 1, the algorithm that obtains the proper E will not fulfill conditions 2 and 3. Furthermore, an algorithm that just reports a huge number of solutions does not ensure their quality (in terms of efficiency). Thus, at least two strategies can be followed: a) looking for a set of solutions uniformly distributed; b) looking for efficient solutions.
To have an idea of the quality of an algorithm, a Reference set of E, R in the following, could be useful. R is the set of known PE solutions for a given instance problem. The ideal R is the set E. A useful practice is, having a set R as close to E as possible, then getting, with the outcome of the algorithm, a net set. The resultant Net set of non-dominated solutions, N = {X is Pareto efficient in (PE ∪R)}, will be at least as good as R (see Figure 2). Combining the P E yielded by different algorithms, an N set of non-dominated solutions, for an instance problem is obtained. The N set is very useful as a reference for many evaluations of new developments. An important contribution is updating the published N set obtained for benchmark instances.
The following cardinal metrics of the quality of the output can be useful: • The percentage of solutions in PE that survive the filtering process with the R set: • The percentage of reference solutions found by the algorithm: • The percentage of solutions in N found by the algorithm: All of the above metrics are cardinal. However, in the case of real-life MOCO problems it may be impossible to obtain, in a reasonable time, a significant percentage of efficient solutions. Obtaining near-efficient solutions would also be highly appreciated. Following [24], a more general and economic criterion may be to concentrate on evaluating the distance of solutions to the efficient frontier. The Dist1 R and Dist2 R metrics by [5], can serve this purpose. They are non-cardinal measures. Their computations, although more complicated than the above mentioned, do not imply a high complexity. They are based on an achievement scalarizing function: where X ∈ R, Y ∈ PE , and Dist1 R is defined as: While Dist1 R measures the mean distance, over the points in R, of the nearest point in PE, Dist2 R defined as: gives the worst case distance. The lower the values of Dist1 R and Dist2 R the better PE approximates R. Moreover, the lower the ratio Dist2 R /Dist1 R the more uniform the distribution of solutions from PE over R. Dist1 R induces a complete ordering and leads to weak outperformance relations.
Another cardinal measure is the C metric introduced by [53], it compares two sets of PE, A and B. A reference set, R, is not required and it is really easy to compute as: The following statements can aid the understanding of C(A, B): When two algorithms are compared, C(A, B) and C(B, A) must be computed, because they are not necessarily complementary. Unless C(A, B) = 1 and C(B, A) < 1, it is not possible to establish that A weakly outperfoms B. C is not difficult to compute, and it seems to be complementary to Dist1 R , Dist2 R , and Dist2 R /Dist1 R with respect to the properties analyzed by [24].
In order to evaluate the computational complexity of the MOO algorithms, beside the CPU time, the number of candidate solutions enumerated by an algorithm during the entire search process can be useful (obviously, an algorithm is computationally less demanding when it enumerates fewer solutions). An interesting piece of information is the proportion of the enumerated solutions that result in potentially efficient solutions.

Computational experiments.
To know the performance of the proposed algorithms which have been described we have tested on benchmark problems from [49]. The benchmark set contains 90 test instances of 9 different sizes, with the number of jobs equal to 20, 50 and 100, and the number of machines equal to 5, 10 and 20. For each size, a sample of 10 instances is presented. Thus, we have solved these 90 instances with MOSAD and MOSAI. The quality of the obtained approximations has been analyzed regarding the Q 1 (PE ), Q 2 (PE ), C(A, B), Dist1 R and Dist2 R measures described in the previous section.
With the PE for the 90 instances obtained by means of the proposed algorithms, added to the algorithms presented in a previous work [29], we have built a net set to be used as reference (R) for this computational experiment.
We have also updated the net sets for the case published in both papers [50] and [35], size instance 50 × 20; and the net sets published in [50], size instance 20 × 20 and 100 × 20. Therefore, we have updated the net sets based on results from: • Net of MOSA1, MOSA2, GPWGA [4], a posteriori [10], MOGLS [19], ENGA [2], published in [50], size instance 20 × 20, 50 × 20 and 100 × 20. • Net of PGA-ALS, MOGLS, ENGA, GPWGA, published in [35], size instance 50 × 20. • PE of the proposed algorithms, MOSAD and MOSAI, and PE of the previous algorithms [29], size instance 20 × 5, 20 × 10, 20 × 20, 50 × 5, 50 × 10, 50 × 20, 100 × 5, 100 × 10 and 100 × 20. Table 1, 2 and 3 show the net set of non-dominated solutions corresponding to the ten problem instances of size 50x20 (the only size for which both papers [35] and [50], reports on net set). Tables 4 to 8 show the net set corresponding to the ten instances of sizes 20 × 20 and 100 × 20, respectively, combining our results with the net sets published in [50]. The superiority of MOSA2 in yielding the greatest number of nondominated solutions is evident. However, MOSA2 becomes the worst algorithm when looking at computational effort information. A comparative study of the computational effort, in terms of the number of sequences generated during the entire search process, shows that MOSA2 is computationally the most demanding algorithm. Table 9 shows the number of sequences enumerated, for the proposed algorithms: for all the problem instances considered, and the data reported in [50] relative to sequences enumerated. In comparison with the rest of the algorithms MOSA2 demands twice as much effort as MOSA1, 3x that of MOGLS, and thirty times when comparing with the proposed algorithms, on average. Table 1. Net set of non-dominated solutions obtained from various multi-objective flow shop scheduling algorithms for the benchmark problems given by [48], size (50×20) (a).  Table 2. Net set of non-dominated solutions obtained from various multi-objective flow shop scheduling algorithms for the benchmark problems given by [48], size (50×20) (b1). For each problem: In the 1st column makespan is indicated. In the 2nd column total flow time is indicated. In the 3rd column the algorithm that yielded the corresponding solution is indicated as follows: 1 is MOSAI; 2 is MOSAD ; 3 is MOSA1; 4 is MOSA2; 5 is PGA-ALS; 6 is IB8S; 7 is DB4S; 8 is IB4S; 9 is JB4S; 10 is NB4S;11 is DB2A; and 12 is JB2A. Numbers in parenthesis, in the title row, indicate the percentages of MOSAI and MOSAD solutions, respectively, in the Net set. Table 3. Net set of non-dominated solutions obtained from various multi-objective flow shop scheduling algorithms for the benchmark problems given by [48], size (50×20) (b2). For each problem: In the 1st column makespan is indicated. In the 2nd column total flow time is indicated. In the 3rd column the algorithm that yielded the corresponding solution is indicated as follows: 1 is MOSAI; 2 is MOSAD ; 3 is MOSA1; 4 is MOSA2; 5 is PGA-ALS; 6 is IB8S; 7 is DB4S; 8 is IB4S; 9 is JB4S; 10 is NB4S;11 is DB2A; and 12 is JB2A. Numbers in parenthesis, in the title row, indicate the percentages of MOSAI and MOSAD solutions, respectively, in the Net set. Table 4. Net set of non-dominated solutions obtained from various multi-objective flow shop scheduling algorithms for the benchmark problems given by [48], size (20×20) (a). For each problem: In the 1st column makespan is indicated. In the 2nd column total flow time is indicated. In the 3rd column the algorithm that yielded the corresponding solution is indicated as follows: 1 is MOSAI; 2 is MOSA1; 3 is MOSA2; 4 is MOGLS; 5 is IB8S; 6 is CB4S; 7 is IB4S; 8 is NB4S; 9 is IB2S; and 10 is DB4S. Table 5. Net set of non-dominated solutions obtained from various multi-objective flow shop scheduling algorithms for the benchmark problems given by [48], size (20×20) (b). For each problem: In the 1st column makespan is indicated. In the 2nd column total flow time is indicated. In the 3rd column the algorithm that yielded the corresponding solution is indicated as follows: 1 is MOSAI; 2 is MOSA1; 3 is MOSA2; 4 is MOGLS; 5 is IB8S; 6 is CB4S; 7 is IB4S; 8 is NB4S; 9 is IB2S; and 10 is DB4S. Table 6. Net set of non-dominated solutions obtained from various multi-objective flow shop scheduling algorithms for the benchmark problems given by [48], size (100×20) (a1). For each problem: In the 1st column makespan is indicated. In the 2nd column total flow time is indicated. In the 3rd column the algorithm that yielded the corresponding solution is indicated as follows: 1 is MOSAI; 2 is MOSAD, 3 is MOSA1; 4 is MOSA2; 5 is IB8S; 6 is CB4S; 7 is DB4S; 8 is JB4S; 9 is NB4S; 10 is DB2A; 11 is JB2A; and 12 is IB2S. Table 7. Net set of non-dominated solutions obtained from various multi-objective flow shop scheduling algorithms for the benchmark problems given by [48], size (100×20) (a2). For each problem: In the 1st column makespan is indicated. In the 2nd column total flow time is indicated. In the 3rd column the algorithm that yielded the corresponding solution is indicated as follows: 1 is MOSAI; 2 is MOSAD, 3 is MOSA1; 4 is MOSA2; 5 is IB8S; 6 is CB4S; 7 is DB4S; 8 is JB4S; 9 is NB4S; 10 is DB2A; 11 is JB2A; and 12 is IB2S. Table 8. Net set of non-dominated solutions obtained from various multi-objective flow shop scheduling algorithms for the benchmark problems given by [48], size (100×20) (b). For each problem: In the 1st column makespan is indicated. In the 2nd column total flow time is indicated. In the 3rd column the algorithm that yielded the corresponding solution is indicated as follows: 1 is MOSAI; 2 is MOSAD, 3 is MOSA1; 4 is MOSA2; 5 is IB8S; 6 is CB4S; 7 is DB4S; 8 is JB4S; 9 is NB4S; 10 is DB2A; 11 is JB2A; and 12 is IB2S.  We have also computed Q 1 and Q 2 metrics and for pair-wise comparison between different algorithms, C(A, B) have been calculated. Because of limited space only average figures of the obtained results are presented. However, we comment on the most important results (Table 10 presents the average values of Q 1 and Q 2 , respectively). It could be observed that the MOSAI yields more solutions which are kept in the final net set (heuristically efficient) than MOSAD, as much considering Q 1 (relative to [PE ]) as Q 1 (relative to [R]).

ETHEL MOKOTOFF
Going into more depth on the quality relations between these two algorithms, we have the C(MOSAD,MOSAI) and C(MOSAI,MOSAD) measures that make clear how many solutions provided by MOSAD are dominated by solutions from MOSAI, and vice versa. Table 11 shows the comparison between the two algorithms. On average, we can affirm that, considering C, MOSAI is better than MOSAD (only for the instance problems of 100 × 10 MOSAD shows a higher value of C than MOSAI). From the 90 problem instances, while MOSAI weakly outperfoms MOSAD in 27 cases, MOSAD weakly outperfoms MOSAI in only 2 cases. However the outperformance of MOSAI over MOSAD does not apply to all the non-cardinal measures like Dist1 R , Dist2 R , and Dist2 R /Dist1 R (see Table 12). Though figures of Dist1 R , and Dist2 R are still disadvantaged for MOSAD, the values of Dist2 R /Dist1 R make it clear that MOSAD presents a more uniform distribution of solutions from PE over R. Figures 1 and 3 show that, in spite of not yielding as many efficient solutions in R, the set of solutions found by MOSAD achieves two desirable aspects: • It presents better (uniform in most cases) distribution in the bi-objective space. • A wide range of values for each criterion is presented. In spite of having a high percentage of non-efficient solutions, MOSAD gives a wide set of near-efficient solutions. Moreover, another advantage of MOSAD is the  Table 11. Averages values of the cardinal measures C and for comparison between solutions yielded by MOSAD and MOSAI.  fact that it demands less effort than MOSAI. Figures in Table 9 show that, on average values, the number of sequences generated during the entire search process represents 80% of those generated by MOSAI, and it is remarkably less demanding than all the algorithms with solutions in the obtained net sets.
To compute the Q 1 , Q 2 , C(A, B), Dist1 R , Dist2 R and , Dist2 R /Dist1 R metrics, the complete set of PE is required. That is why Tables 10, 11 and 12 show the average results for each size of the 90 benchmark instances obtained only for the proposed algorithms.
From the analysis above, it is possible to conclude that the influence of the improvement technique is crucial for the efficiency of the output. While MOSAD (with a direct search improvement) gives a larger number of solutions which help to improve the diversification of the output with less computational effort, MOSAI (with an inverse search improvement) yields a larger percentage of efficient solutions.

5.
Conclusions. In this paper we have presented two new algorithms based on MOSA techniques for a hard multicriteria scheduling problem. Starting with initial permutations obtained by single criteria constructive algorithms, improvements are made by computing lower bounds on the partial scheduling of neighbors, reducing the objective search space. The selection is made according to the criterion that is preferred at each iteration.
Due to the complexity of evaluating the performance of multi-objective algorithms, a set of different metrics have been computed, considering different attributes of each method. Furthermore, net sets of non-dominated solutions for the benchmarks problems of Taillard [49] have been obtained. After an extensive computational analysis, we can conclude that this kind of approach is appropriate to warrant a quick approximation output, which can serve as input for an interactive procedure. The search process should continue in the direction of the decisionmaker's preferences. The proposed algorithms yield a set of potentially efficient solutions requiring less computational effort than other metaheuristic algorithms that have been published in the last few years.
It is not realistic to hope for general meta-optimization methods for MOCO problems. Results of the analysis above give support to the hypothesis which states that specially-developed algorithms, combining general metaheuristic techniques, for specified combinatorial problems, perform properly.
The development of methaheuristics and hyperheuristics for multi-objective methods is still a growing field. The presented algorithms have been proved to be useful, each having their own advantages and disadvantages. We are working now on the application of these MOO approaches to environment problems.