Solving a multi-objective manufacturing cell scheduling problem with the consideration of warehouses using a simulated annealing based procedure

Article history: Received January 3


Introduction
The competition faced by manufacturing companies drives the development of innovative methods for better decision making.These methods must be equally fast enough to accomplish highly competitive results.One of the frequent problems tackled by these methods is the scheduling of manufacturing operations.The complexity associated with these problems limits considerably the possibilities of yielding optimal results for real-world cases, investing reasonable computational effort (Błażewicz et al., 2007;Pinedo, 2012;Tirkolaee et al., 2017).Even with these practical limitations, companies need to make continually short, medium and long-term planning and scheduling decisions.As operations scheduling has a significant economic impact on the performance of manufacturing companies, efficient approaches that swiftly provide competitive solutions, i.e., either optimal or "good enough" solutions, to these problems are required (Błażewicz et al., 2007;Pinedo, 2012;Framiñán et al., 2014).
The manufacturing processes include the diverse operations that take part in the transformation from raw materials and basic components to the final products.The general formulation of a given manufacturing scheduling problem can be represented through mathematical programming (Ruiz et al., 2008;Sayadi et al., 2010;Defersha & Chen, 2012).In this work, we specifically address the Flow Shop Scheduling Problem (FSSP), which is a combinatorial optimization scheduling problem in the NP-Complete class (Garey et al., 1976).Given a set of jobs that have to be processed on a set of machines, FSSP can be defined as the problem of finding the best sequence (order) in terms of the relevant objectives in which those jobs should be processed on the machines, considering that this order should be respected for the whole set of machines.The flow shop production setting is appropriate to enhance industrial productivity and, therefore, more than a quarter of real-world production facilities have implemented it (Pan et al., 2011).There are even attempts to extend its benefits to other production processes that, by their very nature, do not have a flow shop configuration such as the case of manufacturing cells, which consist of groups of machines (cells) that process a set of jobs that have common processing characteristics.This allows adopting, within each cell, a similar configuration to a flow shop system.This production strategy is widely used in the industry (Joines et al., 1996;França et al., 2005).However, productivity improvements resulting from the use of this strategy are usually compromised when the cell is not well integrated with the surrounding manufacturing environment.Among other situations, we can mention an inefficient supply of materials required to perform each job or a poor planning in the location of the final products manufactured by the cell.To our knowledge, flow shop optimization literature has not addressed this problem.It is for this reason that, in this work, we conceived the manufacturing cell in a production system, which must provide the necessary materials for the manufacturing process, as well as, receive and store the final products.
In the scheduling research literature, the most popular performance measure used to evaluate the quality of the solutions for the FSSP is the makespan, which is defined as the total completion time required to process all the jobs in the system.Over time, other performance measures have been used and various operational conditions have been added to the basic configuration.For example, Wu et al. (2018) assessed the carbon footprint and the makespan in a multi-objective approach.Another possibility is to consider production and transport costs as it is performed in Ruiz-Torres et al. (2018), where the authors also pursue the minimization of the number of tardy jobs.In Wang and Tang (2017), the authors minimize the makespan, the total tardiness and the total flow time.Han et al. (2017) address a multi objective lot streaming flow shop problem, considering the makespan and the total tardiness as minimization objectives.In our case, the manufacturing system functions as a manufacturing cell environment and incorporates the raw material supply and the final product delivery.This problem considers the handling of the material for the production process, taking into account the supply of initial goods, which is associated with transport operations, and the storage of the final products, that also includes transport operations from the manufacturing cell to the final products warehouse.To the best of our knowledge these transport operations, which are important since they have a direct impact on the system performance, have not been studied in the literature, at least for manufacturing cells.The optimization stategy presented in this work addresses the makespan and the total tardiness in a multi-objective manner.
This paper contributes with a realistic problem, which is new to the literature.For this problem, a mathematical programming formulation is proposed and an appropriate solution methodology, which allows handling efficiently its combinatorial nature, is developed.The solution methodology is a multiobjective metaheuristic procedure that was adapted from the simulated annealing algorithm (Tirkolaee et al., 2016) and is able to generate a set of solutions that provides an approximation to the optimal Pareto front.We conducted a series of experiments in order to evaluate the performance of the proposed technique.
The article's organization is as follows: first, we present the FSSP and detail its configuration.We then formulate a mathematical programming model for the FSSP.Next, we present the Pareto archived simulated annealing procedure and introduce the augmented ε-constraint method as a comparative solution method.We then solve a set of small-sized problems using both approaches in order to compare the results.Subsequently, we apply the Pareto archived simulated annealing algorithm to larger instances to test its performance.Finally, we outline the main conclusions and future research directions.

Problem Presentation
The production scheduling problem introduces a series of elements of interest in real production settings known as flexible manufacturing cells (Joines et al., 1996;França et al., 2005;Hendizaeh et al., 2008;Li et al., 2015;Fichera et al., 2017).The target problem is restricted to the analysis of the scheduling activities in a production environment constituted essentially by a set of machines arranged in a cell and two warehouses: one for raw material and another one for final product.Certain transfer operations of either raw materials or final products must be carried out between the warehouses and the manufacturing cell.These operations have an impact on the performance measures considered to optimize the sequencing problem.Naturally, the time required by these operations will depend on the type of product being processed (quantity, weight, volume, and shape of parts, components and final products).Centobelli et al. (2016) provide further insight on the subject material handling in warehouses.
In the next section, we propose a mathematical programming formulation to solve the problem that can be characterized as a skip flow shop, where jobs do not necessarily need to be processed in all the machines of the cell, with some additional features found in several production environments.Then, we present a multi-objective metaheuristic procedure to solve the sequencing problem.Some examples of metaheuristics for scheduling problems can be found in Minella et al. (2011) and in Shoaardebili and Fattahi (2015).For a deeper understanding on the issue, the review of Ghosh et al. (2011) presents a wide perspective of metaheuristics implementations on cellular manufacturing scheduling problems.
The performance measures used to evaluate the quality of the schedules take into account the utilization of production resources and the fulfillment of due dates.In addition, we considered raw material and final product warehouses transportation times in the model.In particular, we take into account two objectives when analyzing the quality of the candidate solutions: the minimization of the completion time of all jobs, or makespan, and the minimization of the total tardiness of all jobs.Fig. 1 shows an example configuration for the case of a manufacturing cell connected to a raw material warehouse and to a final product warehouse.

Mathematical Formulation
The production configuration considered in this work is a manufacturing cell separated from the raw material warehoused and the final product warehouse.The jobs have a release date, and after their release, the raw materials for producing the job must be transported from the raw material warehouse to the manufacturing cell.Once the raw materials for the job have arrived at the manufacturing cell, the job can start to be processed.The production configuration inside the manufacturing cell is a flow shop setting, thus all jobs have the same machine sequence.However, some jobs do not require processing in all the machines, namely, they can skip some operations.After a job has finished its processing, it must be transported from manufacturing cell to the final product warehouse.This last transportation time affects the delivery time of the job, since it has to be added to the completion time of the job in the manufacturing cell.As we stated before, the performance criteria considered for assessing the schedules are the makespan and the total tardiness.Other assumptions made in the formulation of this problem are the following:  Each machine can process only one job at a time. Each job can be processed by only one machine at a time. Preemption is not allowed. Machines can be idle during the planning horizon.Sets  Jobs: {j}  Machines: {m} Parameters  pj,m: processing time of job j on machine m  : transport time of job j from the raw material warehouse to the manufacturing cell  : release date of job j  : due date of job j  : transport time of job j from the manufacturing cell to the final product warehouse  Ω: a large positive number Variables


, : completion time of job j on machine m  : tardiness time of job j  : makespan  , : binary variable, it has a value of 1 if job j' is processed before job j, and 0 otherwise.Objective functions (1), minimize the makespan and the total tardiness.

∑
(1) Precedence Eq. ( 2) ensures that job j does not start its processing on machine m before it finishes its processing on machine m 1.
, , , ; ∀ , 1 Ordering Eq. ( 3) requires that job j does not start its processing on machine m before all the previously sequenced jobs j' have been processed.(3) Logic Eq. ( 4) enforces that if job j' is sequenced before job j, then the reverse is not valid.
Raw material supply Eq. ( 5), each product cannot start its processing, before its release and raw materials arrive at the manufacturing cell.
, ; ∀ , ∀ Tardiness in Eq. ( 7), for obtaining the tardiness of job j, it is necessary to consider the transportation time from the manufacturing cell to the final product warehouse.
Non-negativity and binary variables value constraints in Eq. ( 8).

Proposed Solution Approaches to Address the Problem
In this section, we describe the algorithm based on simulated annealing used to solve the problem presented before and introduce the augmented ε-constraint method, used as a comparative solution method.The first algorithm has the structure of a simulated annealing type procedure with an archive of Pareto solutions (PASA, Pareto Archived Simulated Annealing).
Simulated Annealing (SA) is a local search based method developed from an analogy with the phenomenon of annealing to solve complex optimization problems (Kirkpatrick et al., 1983).Local search methods look for the solution with the best value of the chosen criterion in the neighborhood of the current solution, accept this new solution as the current solution, and repeat this step until it is not possible to further improve the solution in the explored neighborhood.By systematically applying this procedure, they obtain, in general, a local optimum of the problem.To avoid getting trapped in a local optimum, a diversifying mechanism should be incorporated with the aim of exploring the entire solution space.In the simulated annealing metaheuristic, the diversifying strategy allows moves, with a certain probability, toward solutions that worsen the current value of the objective function.SA has shown the capability of handling regular flow shop environments (Osman & Potts, 1989;Low, 2005;Vahedi Nouri et al., 2013).

PASA Procedure
This method, proposed in Engrand and Mouney (1998) uses an aggregation function of the objective functions, along with a non-dominated solutions file system.It is assumed that the objective functions to be minimized, fp, p = 1, 2,..., P, are positive.Thus, the problem can be transformed into a mono-objective minimization problem using the following aggregation function: where S is a given sequence.Hence, Eq. ( 10): represents the average relative variation of the objective functions between the current solution (SA) and the candidate solution (SC).If z > 0, SC deteriorates the relative mean of the set of objective functions.
If z ≤ 0, SC improves or maintains the relative mean of the set of objective functions.In the first case the solution SC is accepted with a probability given by ∆ / , where T is the control parameter that simulates the role of the temperature in the physical process of annealing.The method takes into account a non-dominated solutions file that is managed as follows:  If at least one of the solutions in the file dominates SC, it is not added to the file. If SC dominates one or more solutions in the file, SC is added, replacing the solutions that it dominates. If SC does not dominate, nor is dominated by any file solution, SC is added without replacing solutions.
In order to obtain an approximation to the entire efficient frontier of non-dominated solutions during the search process, it is necessary to restart the search regularly from one of the archived solutions selected at random.
PASA incorporates the classic parameters of the simulated annealing algorithm:  T: Control parameter (temperature), a positive real value that varies from an initial higher value, T0, to a final lower value, Tf, during the execution of the algorithm. NT: Number of iterations performed by the algorithm for a certain value of T.  α: function of T, α = α (T), which determines the variation of T. In general: α (T) = α T, in practice: α  [0.80, 0.99]. Nstop: Maximum number of iterations allowed without improvement.
The pseudo-code shown in Fig. 2 is applied to generate a set of potentially efficient (Pareto optimal) solutions.

i. Start
A constructive procedure is applied to generate an initial solution, S0.

ii. Iteration t
A solution is randomly generated in the neighborhood of SA, SC  V(SA).
Otherwise, SC is accepted with the following probability: ∆ / .
A random number  uniformly distributed in the interval [0, 1] is generated: If applicable, update SES taking into account Sc.
In the starting stage, the procedure constructs an initial solution by selecting in a random fashion the job that will be added to the sequence.We considered other constructive heuristic procedures in preliminary tests.However, from these tests we concluded that the simulated annealing procedure is competent enough, so it does not require a sophisticated procedure to initiate the algorithm (that can add unnecessary complexity to the solution approach).
Finally, we defined the neighborhood function used to generate a candidate solution as the set of solutions that can be obtained by randomly swapping two jobs in the sequence.

Augmented ε-constraint Method
In order to find the real Pareto Front of each instance and validate the proposed metaheuristic, we applied the augmented ε-constraint method presented in Mavrotas and Florios (2013).This approach is a novel variation of the original method developed by Haimes (1971), having several advantages not only over the traditional ε-constraint but also over another popular multi-objective approach as it is the weighted sum (Mavrotas, 2009;Mavrotas & Florios, 2013;Rossit et al. 2017).Namely, it guarantees the Pareto optimality of the solution and reduces the computing time by avoiding repeated solutions.We present the AUGMECON.Given a set of P objective functions, i.e., f1, f2,…, fp ,...,fP, that we aim to minimize, the AUGMECON method solves the optimization problem that is described in Eqs.(11)(12)(13)(14)(15).These equations were adapted from the ones presented in Mavrotas & Florios (2013), formulated considering a set of objective functions that have to be maximized.
x ϵ X (15) where rp is the range of the p-th objective function over the efficient set, i.e., the distance between the ideal and the nadir value of the p-th objective function over the Pareto Front.Sp is the slack variable of the constraint associate with the p-th objective function.Analogously to the traditional method εp is the parameter for the Right Hand Side (RHS) for the constraint associated with the p-th objective function.Finally, eps is a small positive constant (between 10 -6 and 10 -3 ).x ϵ X represents the structural constraints of the problem that are not directly associated with the application of AUGMECON.In the initial version of the AUGMECON the coefficients a p were all equal to one (Mavrotas, 2009).However, this was modified in the second version, AUGMECON2, presented in Mavrotas and Florios (2013) in which this values are calculated as a =10 -(p-2) , what results in a lexicographic optimization over the objective functions in case there are alternative optima.Thus, the objective functions are optimized sequentially, i.e., first f 1 , then f 2 and so on.The range of the p-th objective function is divided in qp equal intervals.Thus, we obtain in total qp + 1 grid points that are used to systematically vary the RHS of the restriction associated with the p-th objective function (εp).If rp is the range of the p-th objective function, then the discretization step is: Therefore, the RHS of the corresponding constraint associated with the n-th objective in the t-th iteration is: where fmaxp is the nadir value for the p-th criteria.At each iteration, the surplus variable corresponding to the innermost objective function is checked.In the bicriteria problem presented in this paper, this is the objective function with 2. Then we calculate the bypass integer coefficient as: When the surplus variable S2 is larger than step2 in a t-th iteration, the same solution will be obtained in the next iteration t+1 with the only difference of a smaller surplus variable (that will have a value of S2-step2).This makes t+1 iteration redundant and, therefore, it can be bypassed.The bypass coefficient b actually indicates how many consecutive iterations the procedure can jump ahead.
In integer programming problems, as long as we know the efficient ranges of the optimization criteria and the objective functions have integer values, we can generate the complete set of efficient solutions if we set q k r k (Mavrotas & Florios, 2013).For this paper, we obtained the ranges of the objective functions with lexicographic optimization to avoid weakly efficient values.

Computational Experimentation
In order to assess the mathematical formulation and the algorithm proposed in this work, we considered instances of two different sizes: the first group of instances involves up to 20 jobs and 20 machines and the second group of instances involves up to 100 jobs and 20 machines.While the first group is appropriate for validating the PASA algorithm through the comparison with the AUGMECON method, the second group of instances shows the capability of the proposed PASA algorithm for solving larger problems.In Table 1, the combinations of number of machines and number of jobs used to generate the instances are detailed.As we explained in the mathematical formulation section, several parameters must be set to define each scenario within each instance.In the case of processing times, pj,m, the data are obtained from a pseudouniform distribution in the interval [0,100].Since the problem addressed in this paper contemplates skipping operations, the zero value is included in the processing time distribution.Furthermore, and as a special feature, the probability of skipping an operation is set at 3%, that is, not all the possible values have the same probability (that is why there is a reference to a pseudo-uniform distribution).This amplifies the impact of skipping operations in the scheduling problem.For the transportation times from the raw material warehouse to the manufacturing cell and from the manufacturing cell to the final product warehouse, we take the data from a uniform distribution in the interval [10,20].Similarly, we generate the release dates of the jobs from a uniform distribution in the range [1,100].With reference to the due date of the jobs, we followed the guidelines given in Ruiz and Stützle (2008).However, in that article the due dates are calculated using the following expression: ∑ , • 1 • 3 , which contemplates all the operations that job j must pass through for its completion plus an additional time considered by the term random, which is a random number uniformly distributed in [0,1].In our case, each job has more intervals to consider, namely, the transportation times associated with arriving and leaving the manufacturing cell.Thus, in this paper the due date of job j is calculated as follows: , ∑ , • 1 • 0.3 .Moreover, in our due date formulation the constant that multiplies the random number has been reduced from 3 to 0.3 for the small instances (with up to 20 jobs).This produces tighter due dates, giving a higher relevance to the total tardiness with respect to the makespan in the multi-objective evaluation.For each instance, we generated five different sets of parameters (scenarios).For the large instances (more than 50 jobs) the constant that multiplies the random number used in the due date generation was set at 3, respecting the original proposal of the cited paper.
The instances solved using the mathematical programming formulation were addressed using Pyomo (Hart et al., 2011;Hart et al., 2012) as modelling language, and CPLEX 12.5.0 as a solver.The tests were run on a CPU with 8GB of RAM and an I-core 5 processor, and a 64-bit operating system, the run time was limited to 1 hour (3600 seconds).CPLEX could solve instances up to 10 jobs and 10 machines, for the following instance (15 jobs and 10 machines) it could not find a solution with a gap lower than 3%, in average, the gap reported by CPLEX after an hour of computing was larger than 30%.We discarded these solutions.We programmed the instances solved using the simulated annealing metaheuristic in VBA (Visual Basic for Applications) as a modelling language in an MS-Excel setting.The tests were run on a CPU with 8GB of RAM and an I-core 7 processor, with a 64-bit operating system.

Experimental Design
The intention of the experimental design is to assess the ability of the proposed PASA algorithm to solve the problem at hand (Tirkolaee et al., 2018a;Tirkolaee et al., 2018b;Tirkolaee et al., 2018c;Zokaee et al., 2017;Amiri et al., 2018).For this purpose, we need to solve problems of different sizes.However, given that the PASA algorithm is a metaheuristic, we cannot assure the quality of the results by providing only a set of solutions that approximates the Pareto front in the best possible way.Furthermore, given the novelty of the problem, there is no known real Pareto front for the instances, except for the small ones, and due to its hardness, one cannot expect to calculate practically the optimal front.Thus, only the solutions obtained by PASA are considered, then, it remains to evaluate the precision of the solutions obtained with PASA with an exact method (a mathematical programming formulation solved with CPLEX in this case).In this way, it is possible to evaluate whether a given solution obtained with PASA is a real Pareto-optimal solution or not.
The methodology used to evaluate the set of the non-dominated solutions obtained by PASA, consists of assessing whether they are real non-dominated solutions of the problem, i.e., if they are Pareto optimal.
For each solution x, two problems will be solved using the ε-constraint methodology.Each nondominated solution obtained by PASA is defined by a pair of values of the total tardiness and the makespan objectives ∑ * , * , in such a way that PASA cannot find other solution that improves one of the objectives without worsening the other, i.e., a new solution that dominates x.Therefore, we will evaluate the solution x using the ε-constraint method following a hypothesis test procedure: the hypothesis is that x is Pareto optimal.We will perform the hypothesis test demonstrating that it is not possible to improve it.What is done is to take the value ∑ * of x and fixing it as a constraint in the MIP model.Then, we try to minimize in a mono-objective manner the makespan, solving the mathematical programming formulation with CPLEX.If the resulting makespan coincides with * of x, and since CPLEX is an exact method, we verify that this makespan value corresponds to the actual Pareto front.We do the same with the other objective, we fix * as a constraint and we solve the MIP model by minimizing the total tardiness.Similarly if the resulting total tardiness coincides with ∑ * , the value of total tardiness of x belongs to the real Pareto front.If these both value matches occur, then the solution x is an optimal Pareto solution of the problem.Conversely, if one (or both) of these matches do not happen, then the exact method found a solution that dominates x and, thus, x is not Pareto optimal for the problem.

Results
As a first result of the problem, we present Figure 3, which shows the conflicting nature between the total tardiness and makespan objectives.This figure depicts the solutions obtained by PASA for scenario 3 of instance 50J_20M.We can see that when we minimize the makespan, the total tardiness objective deteriorates.Likewise, the same occurs in the other direction of the relation between objective functions: as the total tardiness reduces, the makespan increases.The first stage of the evaluation consists of solving the problem with PASA and then evaluating if these solutions belong to the real Pareto front of the problem.To do this, we assess the PASA's solutions using the procedure described in the experimental design.In Table 2, we present the summary of the results of these assessments (we detail the information of each of the solutions in the Appendix Table ).Table 2 shows for each of the scenarios tested for each instance: the average and the maximum number of solutions obtained with PASA; the percentage proportion of the PASA solutions that were verified to be Pareto optimal solutions of the real problem; and the number of solutions obtained with the augmented ε-constraint method, which represents the complete Pareto Front.In all the scenarios the PASA algorithm manages to find optimal Pareto solutions.
As we show in Table 2, the PASA algorithm was able to effectively address the problem for small instances in a very stable manner since the difference between the maximum number of found solutions and the average number in most of the cases is the same.
In the next stage, we evaluate the PASA algorithm performance for solving in the larger instances, i.e., the ones that have from 15 to 100 jobs and from 10 to 20 machines.For this purpose, we studied the hypervolume indicator obtained in the experimentation.We considered the average hypervolume indicator of the 30 runs that were completed for each scenario.Moreover, we presented the standard deviation and the coefficient of variation to provide a level of representation of the variability of the algorithm in the solution process.We depict all these results in Table 3.The average coefficient of variation for all the runs performed (7 instances × 5 scenarios × 30 runs per scenario: 1050) is 6.53%.Thus, the PASA algorithm is stable, obtaining always very similar solutions in general terms.When comparing the coefficients of variation of the scenarios for a given instance, we observed that there were variations between scenarios of the same problem size in terms of machines and jobs.This may be due to the fact that the data set of one scenario is more difficult to solve than the others, as it is explained in Vallada et al. (2015).Related to this, in order to avoid the case where the only tested scenario is an atypical tricky one, we evaluate more than one scenario for each instance.Nevertheless, the interesting issue to analyze is the convergence of the algorithm regarding the size of the problem, in order to do that we present Table 4. Table 4 shows the average coefficient of variation for the 5 scenarios of each instance.As the size of the instance increases in terms of the number of machines and jobs, the coefficient of variation also tends to increase.This becomes more evident when passing from the instances of 20 jobs to the instances of 50 jobs.We can explain this behavior by the fact that increasing the size of the problem considerably increases the size of the space of feasible solutions.Thus, the algorithm does not always find the same solutions.However, this increase in the coefficient of variation appears to be bounded, since, for the largest instance of 100 jobs and 20 machines, it decreases with respect to the instances of 50 jobs and 20 machines.

Conclusions
In this work, we have studied a particular case of the FSSP (Flow Shop Scheduling Problem).This problem takes into account the material handling for the productive process considering both the supply of initial goods, with its associated transport operations, and the store of the final products, that also includes transport operations from the manufacturing cell to the final products warehouse.We have proposed a mathematical formulation for this realistic problem and presented a multi-objective simulated annealing algorithm (PASA).We have assessed the PASA solutions quality not only in terms of convergence, with a hypothesis test, but also in terms of diversification, with a comparison to the complete Pareto front obtained with the augmented ε-constraint method.These assessments have indicated that the proposed PASA algorithm was capable of addressing this problem, having a high performance on both convergence and Pareto front coverage.
We also addressed problems of larger size using the Pareto archived simulated annealing algorithm.For this set of problems, we observed that the results obtained by means of the algorithm were stable in terms of precision and the solutions were obtained in practical computational times.
Future research will consider larger instances of study and different performance measures with the aim of developing an accurate and fast method for solving this kind of scheduling problems.

Fig. 1 .
Fig. 1.Example configuration of a manufacturing cell with warehouses

Table 1
Detail of the test instances

Table 2
Results of Pareto efficiency tests for PASA's solutions

Table 3
Hypervolume analysis for the larger instances

Table 4
Hypervolume average coefficient of variation for each instance . Handbook on scheduling: from theory to applications.Berlin Heidelberg: Springer Science & Business Media.the authors; licensee Growing Science, Canada.This is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC-BY) license (http://creativecommons.org/licenses/by/4.0/).