Flexible Job Shop Scheduling Problem with Sequence Dependent Setup Time and Job Splitting: Hospital Catering Case Study

: This paper aims to study a real case of an optimization problem derived from a hospital supply chain. The present work focuses on developing operational decision support models and algorithms for production process scheduling in hospital catering. The addressed production system is considered as a ﬂexible job shop system. The objective is to minimize the total ﬂow time. A novel mathematical model and two metaheuristics for the production scheduling of multi-product and multi-stage food processes are developed. These methods have proven their effectiveness for the scheduling of operations of the food production processes and allowed signiﬁcant improvements in the performance of the studied production system.


Introduction
Nowadays, hospital logistics has become an essential component of healthcare institutions. It allows the synchronization of all the flows inside a hospital to ensure the efficiency of the healthcare system. For many years, the management was commonly focused on improving the quality of medical care, while less attention was usually devoted to operation management. In recent years, the need for containing the costs while increasing the competitiveness along with the new national health service policies for hospital financing forced hospitals to necessarily improve their operational efficiency. It is in this context that the efficient use of resources and the research on optimal service stimulate logistical thinking in hospitals. The difficulties of optimizing flows have led managers to discover new avenues for rationalizing expenses and seeking refined solutions to these difficulties. In this context, optimized logistics solutions allow hospitals to improve inventory management, limit waste, and provide better inventory tracking and traceability of service products. On the other hand, the supply chain is a major source of costs, and its reorganization would make it possible to achieve crucial savings on all hospital expenses.
A hospital's logistics is part of its global performance, where the activities are organized and structured with the aim of patients' satisfaction in terms of quality, quantity, delay, safety, and low cost. The main purpose of this logistics is to control and optimize physical flows from suppliers to patients at the best cost that respects technical, economic, and regulatory conditions for optimal dispensing to patients. Hospital logistics is a complex process characterized by a diversity of needs, users, products, and distribution channels. The coordination of these activities requires logistical expertise that few institutions will be able to develop on their own. This has led researchers to focus for some years on the management and optimization of the supply chain in hospitals. In this context and in order to improve the working conditions of the employees and their well-being, the hospital center of Troyes implements means to improve its daily efficiency. The hospital is carrying out a revision of its supply chain, which must notably consider the management of food flows within the hospital. In the present work, we focused particularly on the scheduling of the food manufacturing process in hospital catering, which is considered as a flexible job shop scheduling problem with a sequence-dependent setup time and job splitting by integrating specific industrial constraints.
The remainder of this paper is organized as follows: Section 2 presents the state of the art regarding the problem of food production process scheduling. The problem statement is defined in Section 3. In Section 4, the mathematical model developed for the studied problem is presented. Sections 5 and 6 present the genetic and the iterated local search algorithms specifically developed for the problem, the different elements of these metaheuristics, and the computational results. Finally, in the last section, an application to a real industrial case and the results obtained are presented.

Literature Review
The production scheduling problem in food industries belongs to a famous class of problems referred to as scheduling with sequence-dependent setups, which are well known to be NP-hard (Sun et al. [1]). In recent years, there has been great interest in the development of intelligent solutions for this problem in various fields of application. The promising results of scheduling methods, such as reduction of production costs, increased throughput and smoother operation of the production equipment, and improvement of working conditions and the well-being of employees, have stimulated considerable research efforts. The existing works in the literature are classified according to the number of products (single products or multiple products), the type of production system, and the expiration dates of products, which may be known or unknown.
Regarding works dealing with a single product, Entrup et al. [2] proposed three different mixed-integer linear programming for scheduling problems in the packing stage of stirred yogurt production in the fresh food industry. They accounted for shelf life issues and fermentation capacity limitations. Doganis and Sarimveis [3] proposed a model that aims to obtain the optimal production scheduling in a single yogurt production line. The model takes into account all the standard constraints encountered in production scheduling (material balances, inventory limitations, machinery capacity). It also considers special features that characterize yogurt production and that are limitations in production sequencing, mainly due to different fat contents and flavors of various products, as well as sequencedependent setup times and costs. However, the model is limited to a single production line. In another study, Doganis and Sarimveis [4] presented a methodology for optimum scheduling of yogurt packaging lines that consist of multiple parallel machines. The methodology incorporates features that allow one to tackle industry-specific problems, such as multiple intermediate due dates, job mixing and splitting, product-specific machine speed, minimum and maximum lot size, and sequence-dependent changeover times and costs. However, the proposed mathematical model does not incorporate multi-stage production decisions and ignores some industry-specific characteristics, such as shelf life. Stefansdottir et al. [5] developed a generic optimization model for lot sizing and scheduling in the typical processing industry setting of flow shops. Sargut and Isık [6] presented a mathematical model for a dynamic economic lot sizing problem with a single machine for a single perishable item under production capacities. They also gave a dynamic programming-based heuristic for the solution of the overall problem.
Many studies have been carried out in the literature on the production scheduling of multi-product food processing industries. Akkerman and van Donk [7] developed a methodology for the analysis of the scheduling problems in food processing. This helps one to understand, describe, and structure scheduling problems in food processing and to evaluate the organizational structures and information flows related to scheduling. Smith, Daniels, and Larry [8] developed a general lot-sizing model for processing industries and applied their method to a representative situation of a food processing facility. Kopanos et al. [9] offered an efficient mathematical framework for detailed production scheduling in the food processing industry. Wauters et al. [10] introduced an integrated approach for production scheduling and demonstrated its applicability to the food processing industry. In this work, the scheduling had to deal with a combination of discrete, continuous, and batch processes, and it was complicated by particular characteristics. Acevedo-Ojeda et al. [11] presented a lot-sizing problem with a single machine that incorporated raw material perishability and analyzed how these considerations enforced specific constraints on a set of fundamental decisions, particularly for multi-level structures. Three variants of the two-level lot-sizing problem incorporating different types of raw-material perishability-fixed shelf life, functionality deterioration, and functionality-volume deterioration-were studied; the authors proposed mixed-integer programming formulations for each variant and performed computational experiments with sensitivity analyses. Copil et al. [12] considered a capacitated dynamic lot-sizing problem with parallel machines for the food industry, in which a given product produced during a specified time period is used to satisfy the related demand. Niaki et al. [13] addressed the integrated lot-sizing and scheduling problem of food production in batch manufacturing systems with multiple shared common resources and proposed a new mixed-integer linear programming formulation with multiple objective functions. Wei et al. [14] proposed a classical multi-level lot-sizing and flow-shop scheduling problem formulation to incorporate perishability issues.
Regarding the shelf life of products, Ahumada and Villalobos [15] reviewed models for the agri-food business, where products may be perishable or not, but their focus was on procurement and harvesting planning. The only goods they were interested in were crops. Sel et al. [16] introduced the planning and scheduling of decisions considering the shelflife restrictions, product-dependent machine speeds, demand due dates, and regular and overtime working hours in the perishable supply chain. Arbib et al. [17] considered a threedimensional matching model for perishable production scheduling, which was studied under two independent aspects: the relative perishability of products and the feasibility of launching/completion time. Basnet et al. [18] described an exact algorithm to solve a scheduling and sequencing problem in the same industry. Chen et al. [19] provided a review of literature on the integration of scheduling and lot sizing for perishable food products, and they categorized the papers by the characteristics of the lot sizing and scheduling that were included in their models, as well as the strategies used to model perishability.
In the present work, the studied production system is considered as a flexible job shop system. Since 1990, the flexible job shop scheduling problem (FJSP) has been extensively investigated by researchers. Liu and MacCarthy [20] discussed the problem in a flexible manufacturing system with transportation times and limited buffers. They developed a mixed-integer linear programming model and heuristics to minimize the makespan, the mean completion time, and the maximum tardiness. Guimaraes and Fernandes [21] proposed a genetic algorithm for the FJSP, where the objective function is used to minimize the makespan and the mean tardiness. Saidi-Mehrabad and Fattahi [22] took into account a special case of FJSP, where each operation could be assigned to one of two parallel machines. A tabu search algorithm for solving the sequencing and assignment problems was developed. The algorithm was compared to a branch and bound algorithm on several random test instances. Defersha and Chen [23] studied the FJSP with attached and detached setup times, machine release dates, and technological time lag constraints. For this problem, a mixed-integer linear programming model was proposed, and a parallel genetic algorithm using multiple threads was introduced. Mati et al. [24] proposed a genetic algorithm for an FJSP in which blocking constraints were taken into consideration. Bagheri and Zandieh [25] developed a variable neighborhood search algorithm. Mousakhani [26] presented a mixedinteger linear programming model and an iterated local search algorithm that minimize the total tardiness. Chaudhry and Khan [27] published a literature review on the methods used to solve the FJSP and the studied objective functions (Figure 1). In this bibliographic study, it was found that almost 59% of the papers used hybrid methods (Rajabinasab and Mansour [28], Geyik and Dosdogru [29], Zhou and Zeng [30]) or evolutionary algorithms. The minimization of the makespan turned out to be the most widely studied criterion. In 88 research papers (44.67%), the makespan was used as the sole objective function, while in another 78 papers (39.59%), makespan was used in combination with another objective function. However, the minimization of the total flow time, which is the target criterion in this work, has been studied very little in the literature. Table 1 represents some research works on flexible job shop problems with minimization of flow time. These works are classified according to the fields of its application in research or an industrial case, the resolution methods used, some workshop details, and the studied objective functions. Finally, it is worth mentioning that, to the best of our knowledge, there are almost no studies addressing the problem of scheduling food production processes in hospital or collective catering. Most existing works in the literature on scheduling food production processes are from the food and dairy industries, where the production systems are parallel machine, flow shop, and single machine, in most cases. Moreover, in the majority of these works, the expiration date of products is unknown, while in this study, it is known. The works found in the literature cannot be adapted to the problem addressed in this study, since the production systems are different. In addition, several constraints specific to the considered problem have not been taken into account in the existing works, and they do not have the same objective of optimization.

Problem Description
The problem of food production process scheduling considered in this study aims to schedule the operations from the pre-treatment of raw materials to the stock of finished products of a meal manufacturing process in hospital catering or, more generally, in collective catering. Figure 2 represents the typical production areas and the different material resources available for realizing the operations of the meal-making process. This problem is considered as a flexible job shop problem with sequence-dependent setup time and job splitting by integrating specific industrial constraints, such as the days of pre-treatment and production of products, the product delivery times, and the amounts of human and material resources available, in addition to the different constraints described below.
The addressed problem is considered as a flexible job shop scheduling problem with sequence-dependent setup time, since each job has its own order of operations and each operation has to be assigned to one among a set of possible machines. This problem can be described by a set of N jobs, where each job i corresponds to the preparation of a dish characterized by a number of portions Q i (quantity) and a set of operations J i for the preparation of the dish (from raw material to finished product). It is worth highlighting that the dishes to be prepared do not have the same order of operating ranges (set of operations necessary for the preparation of the dish). For each job, there is a due date D i to respect.
For each operation of an operating range, there is a set of material resources that are used to realize it, such as ovens, cooking pots, cookers, induction hobs, packaging machines, cooling cells ( Figure 2). These material resources can be classified into three categories: M 1 , M 2 , and M 3 . M 1 ⊂ M is composed of material resources with a capacity of one portion and that can not process several jobs at the same time (material resources that can perform preprocessing and cold production operations). The set M 2 ⊂ M represents the material resources with a capacity greater than one portion and that cannot process several jobs at the same time (ovens, etc.). M 3 ⊂ M correspond to the set of material resources with a capacity greater than one portion and that can process several jobs at the same time (cooling cells). For each material resource, there is a setup time to take into account, which corresponds to the preparation time of the machine before carrying out an operation and the cleaning time of the machine between two consecutive operations. A time window of availability is known for each material resource.
Note that the corresponding machines may not be identical, involving different processing times according to the chosen machine. The setup times of machines are sequence dependent because they depends on the preceding operation on the same machine. The food production process scheduling involves two steps: (i) assignment of operations to machines and (ii) sequencing of operations on machines.
As mentioned previously, in order to respect the production capacities of material resources, a job can be split into smaller sub-lots in such a way that the operations of sub-lots of a job can be performed simultaneously on different machines. This strategy, which is useful when machine capacity does not allow the treatment of the whole job, also enables a more efficient processing scheme.
The criterion to minimize in the present study is the total flow time of jobs in the production system. The choice of this criterion is based on the fact that the respect of the cold chain at each stage of the product life cycle must be ensured. The aim is to constantly maintain a low temperature (positive or negative, depending on the product) to ensure the maintenance of all the food qualities (hygienic, nutritional, and gustatory).

Mathematical Model
In this section, a mixed-integer linear programming model is presented. This model formalizes and solves small-sized instances of the studied problem by using the Cplex solver. The solutions of the small-sized instances can be used to validate the efficiency of the developed metaheuristic methods.

Assumptions
The mathematical model for the scheduling of the food production process inherits its main assumptions from the standard flexible job shop scheduling problem with sequencedependent setup times and additional features due to the job splitting:

•
Jobs are independent of each other, • A job can be split into sub-lots, • The sub-lots of a job can be grouped on the machines to be treated at the same time, • Each sub-lot of a job consists of a set of operations that must be processed consecutively (precedence constraints between operations of sub-lots of jobs), • Each operation of a sub-lot has a given processing time, • The preemption of operations of sub-lots of jobs is not allowed, i.e., operation processing on a machine cannot be interrupted, • Each job has a given due date (finish date of production at latest), • Sub-lot sizes (number of portions) are discrete, • Sub-lots creation is consistent throughout the processing sequence, meaning that job splitting and sub-lot sizes remain constant for all operations, • Machines are independent, • A machine can process, at most, one operation of a job at the same time, • The setup times of machines are dependent on the sequence of operations of sub-lots of jobs, • Material resources have given availability time windows that must be taken into account.
Taking these assumptions into account, the objective is to find a schedule involving sub-lot assignments to machines and sub-lot sequencing for each machine in such a way that each job's demand is fulfilled, different constraints of the problem are respected, and the total flow time of jobs in the production system is minimized.

Notations
The definition of the developed mathematical model's parameters relies on the following sets and indexes: • M: set of all material resources, where m = |M|. • N: set of jobs (dishes to prepare), where n = |N| and {0, n + 1} are two dummy jobs. • J i : set of operations of job i ∈ N, such that the operation j ∈ J i is done before the operation j + 1 ∈ J i and |J 0 | = |J n+1 | = 1.
set of material resources that can perform the operation j ∈ J i of job i ∈ N. • R k : maximum capacity in number of portions of the material resource k ∈ M.
• P ijk : unit processing time of operation j ∈ J i of job i ∈ N on the material resource k ∈ M 1 .

Decision Variables
• X iljk : binary variable that equal to 1 if operation j ∈ J i of sub-lot l ∈ L i of job i ∈ N is assigned to the material resource k ∈ M ij , and 0 otherwise. • F iljhl gk : binary variable that is equal to 1 if operation j ∈ J i of sub-lot l ∈ L i of job i ∈ N directly precedes operation g ∈ J h of sub-lot l ∈ L h of job h ∈ N on the material resource k ∈ M ij ∩ M hg , and 0 otherwise. • Z ill jk : binary variable that is equal to 1 if operation j ∈ J i of sub-lot l ∈ L i of job i ∈ N starts and finishes at the same time as the operation j ∈ J i of sub-lot l ∈ L i of job i ∈ N on the material resource k ∈ M ij , and 0 otherwise.

Mathematical Model
In the literature, there exist different mathematical model formulations for the standard flexible job shop scheduling problem. The mathematical model proposed by Buddala and Mahapatra [31] takes into account constraints on processing times, precedence between operations, and machine capacity in number of operations processed at the same time, which corresponds with the studied problem. However, this model does not take into consideration the constraints on machine capacity in the number of portions, machine availability time windows, sequence-dependent setup time, due dates of jobs, splitting of jobs, and batching of sub-lots.
The developed mathematical model that integrates all these constraints is formulated as indicated through Equations (1) to (22): In the mathematical model presented previously, Equation (1) represents the objective function, which consists in minimizing the total flow time of jobs in the production system. It is defined as the sum of the completion times of all jobs since the release dates are equal to zero. The completion times of jobs are computed as the completion time of the last sub-lot derived from the considered job, as indicated in Equation (2). Note that, due to Equation (3), for given i ∈ N, l ∈ L i , and j ∈ J i , variables S iljk and C iljk are equal to zero if the machine k ∈ M ij is not chosen to realize the operation of the considered sub-lot. On the other hand, when X iljk is equal to 1, Equations (4) and (5) activate the relationship between the starting time and completion time of an operation of a sub-lot. In this case, the processing time does not depend on the quantity of jobs for the material resources M 2 ∪ M 3 (Equation 4), but it depends on the quantity of jobs for the material resources M 1 Equations (5) and (6) consider sequence-dependent setup times between the completion time and starting time of two operations of sub-lots that are processed on a machine one after another. Equation (7) disable Equation (6) if two different sub-lots of the same job are performed at the same time by the same material resource. Equations (8) and (9) require that if two operations of two different sub-lots of the same job are assigned at the same time to a material resource M 2 or M 3 , they must have the same respective starting time and completion time. Equation (10) ensure that only one operation immediately follows the jth operation of sub-lot l ∈ L i of job i ∈ N on machine k ∈ M ij ∩ M hg , and Equation (11) guarantee that only one operation immediately precedes the gth operation of sub-lot l ∈ L h of job h ∈ N on machine k ∈ M ij ∩ M hg . Equation (12) establishes the precedence constraints between two consecutive operations of the same sub-lot. Equation (13) enforce that each operation of each sub-lot should be assigned to exactly one machine among the possible ones. The respect for the due dates of jobs is modeled by Equation (14). Equation (15) ensure that the capacities of material resources in the number of portions are respected. The respect for time windows of availability of material resources is modeled by Equations (16) and (17). Finally, Equations (19)- (22) define the domain of the decision variables.

Computational Results of the Mathematical Model
The mathematical model presented previously was implemented in the Java programming language using the Cplex library. This mathematical model was tested on 150 instances of different types: real instances of the hospital center of Troyes (HCT), randomly generated instances of the HCT type, randomly generated instances, and instances adapted from the literature (Sriboonchandr et al. [32], Nouri et al. [33], Azzouz et al. [34], Mousakhani [26], Lee et al. [35], Bagheri and Zandieh [25], Pezzella et al. [36]). The details of the results of these different types of instances are presented in Section 6. The mathematical model was able to find solutions in less than three hours of execution for the small instances with less than eight jobs, 10 sub-lots, 39 operations, and 29 machines. The execution times of the mathematical model for these instances varied according to the number of jobs, sub-lots, and operations. The computational results of the proposed mathematical model for different types of instances and, specifically, for real instances of HCT show the limits of an exact resolution for the problem of scheduling of the food production process.

Genetic Algorithm
Solving the flexible job shop problem is known to be strongly NP-hard (Xia and Wu [37], Fattahi et al. [38]). The introduction of sequence-dependent setup time and job splitting complicates the already difficult flexible job shop problem. In order to solve this problem efficiently, we developed a hybrid method that combines a genetic algorithm and three local research methods. The elements of this hybrid method are presented in the following subsections.

Solution Representation
By solving a flexible job shop scheduling problem using a genetic algorithm, Kacem [39] used a solution representation by coding both the assignment and the sequencing of operations on different machines. A similar representation is used in this paper to solve the flexible job shop scheduling problem with job splitting by considering each sub-lot as a job. To illustrate this representation, let us consider a small example with three jobs and four machines. The numbers of operations of sub-lots for each job and all the machines eligible for each operation are given in Figure 3. A representation of the assignment of operations to machines and their sequencing is coded in a chromosome, as shown in Figure 3. In this chromosome, each gene is represented by a quadruple (i, l, j, k), designating the assignment of operation j of sub-lot l of job i to machine k. The sequence of genes in the chromosome represents the sequence of operations on machines. The chromosome decoding procedure is carried out by reading it from left to right. For example, the assignment and sequencing of operations on machine 1 can be decoded as follows: (1, 1, 1, 1) → (1, 2, 1, 1) → (3, 1, 3, 1) → (3, 2, 3, 1) → (3, 3, 3, 1). This information is obtained from genes 3, 8, 10, 16, and 17 in the chromosome, where k = 1. In this chromosome, for a given i and l, the gene (i, l, j, k) is always located on the right of all the other genes (i, l, j, k ) with j < j. This ensures that the precedence requirements of the operations of a particular sub-lot are not violated.
operations on machines. The chromosome decoding procedure is carried out by reading it from left to right. For example, the assignment and sequencing of operations on machine 1 can be decoded as follows : (1, 1, 1, 1) → (1, 2, 1, 1) → (3, 1, 3, 1) → (3, 2, 3, 1) → (3, 3, 3, 1). This information is obtained from genes 3, 8, 10, 16 and 17 in the chromosome, where k = 1. In this chromosome, for a given i and l, the gene (i, l, j, k) is always located on the right of all the other genes (i, l, j, k ) with j < j. This ensures that the precedence requirement of the operations of a particular sub-lot are not violated.

Initial population
The initial population plays an important role in the performance of genetic algorithms. An initial population with great diversity between solutions can avoid falling into a premature convergence or a local minimum. Different assignment and sequencing rules have been used to achieve the diversity of solutions : • Random assignment (RA) : the operations are randomly assigned to machines.
• SPT assignment (SPTA) : for each operation, the machine with a smaller processing time is selected to perform this operation.
• LPT assignment (LPTA) : for each operation, the machine with a longer processing time is selected to perform this operation.
• Minimum machine workload assignment (MMWA) : the operations are iteratively assigned to machines based on their processing times and machine workloads. The workload of a machine depends on its type. For the set of machines M 1 , the workload is the sum of processing times of operations assigned to the machine. For the set of machines M 2 and M 3 , the workload is the total machine occupation time. The procedure consists in finding, for each operation, the machine with the minimum workload.
• Random sequence (RS) : this heuristic randomly orders the operations on each machine.
• SPT sequence (SPTS) : operations with the shortest processing time will be firstly processed.  For the construction of the initial population using the heuristics described above, the procedure consists in randomly choosing an assignment and sequencing heuristics and building solutions such as the constraints of precedence between the operations of sub-lots of jobs, the due dates of jobs, the time windows of availability of material resources and the production capacities of machines are respected.

Initial Population
The initial population plays an important role in the performance of genetic algorithms. An initial population with great diversity between solutions can avoid falling into a premature convergence or a local minimum. Different assignment and sequencing rules were used to achieve the diversity of solutions: For the construction of the initial population using the heuristics described above, the procedure consists of randomly choosing an assignment and a sequencing heuristic and building solutions, such as the constraints of precedence between the operations of sub-lots of jobs, the due dates of jobs, the time windows of availability of material resources, and the respect for the production capacities of machines.

Fitness Evaluation
Algorithm 1 decodes the chromosomes to obtain the objective function value and a workable representation of the solutions.
Algorithm 1: Evaluation steps of the fitness function of a chromosome

•
Step 02: Set i, l, j, k, the index values of the gene located at the P position of the chromosome.

•
Step 03: Calculate the completion time C iljk .
• If operation j of sub-lot l of job i is the first operation assigned to the machine k and j = 1, • If operation j of sub-lot l of job i is the first operation assigned to the machine k, j > 1, and the operation j − 1 is assigned to the machine k , • If operation j of sub-lot l of job i is the operation to be processed immediately before the operation j of sub-lot l of job i on the machine k and j = 1, • If operation j of sub-lot l of job i and the operation j of sub-lot l of job i are assigned to the machine k and j = 1, • If operation j of sub-lot l of job i is the operation to be processed immediately before the operation j of sub-lot l of job i on the machine k, j > 1, and the operation j − 1 is assigned to the machine k , • If operation j of sub-lot l of job i and the operation j of sub-lot l of job i are assigned to the machine k, j > 1, and the operation j − 1 is assigned to the machine k , Step 04: If p is less than the total number of operations of sub-lots of jobs, increment its value by 1 and go to Step 2; otherwise, go to Step 05.

•
Step 05: Calculate the fitness function of the solution by using the equation ∑ i C i .
It is worth noting that for a given i and l, the gene (i, l, j, k) is always located on the right of all the other genes (i, l, j , k ) with j < j. Based on this property, when the completion time of operation (i, l, j, k) on machine k is to be calculated, the completion time of operation (i, l, j − 1, k ) is already calculated and available, regardless of the machine to which this preceding operation is assigned. Moreover, the completion time of the operation (i , l , j , k) to be processed on machine k immediately before operation (i, l, j, k) is also calculated and available.
Genetic operators evolve the population to promising regions of the research space. The convergence behavior of the algorithm depends largely on them. These operators are generally categorized as selection, crossover, and mutation operators.

Selection Operator
The process of selecting two parents from the population for reproduction is called selection. The aim of the selection operator is to highlight individuals with the best fitness in hopes that their resulting offspring are more fit individuals. In the proposed genetic algorithm, we used the two-way tournament selection operator introduced in [40]. The selection operator is involved in holding competitions among randomly selected individuals and choosing the one with the best fitness (smallest total flow time). This individual is added to a mating population, which is used to form the next generation. Then, the individuals in the tournament are placed back in the current population, and the process is repeated until the number of individuals added to the mating population is equal to the population size.

Crossover Operator
After the selection of chromosomes for reproduction, a crossover operator is applied to combine the features of the parent chromosomes in order to produce a new child. The individuals in the mating population are randomly paired to create child individuals. For two parents, the algorithm arbitrarily selects one of the available crossover operators and applies it with a certain probability to create two child individuals by exchanging information contained in the parent chromosomes.
Different crossover operators are proposed and can be categorized as assignment or sequence crossover operators. The assignment crossover operators consist of generating offsprings by exchanging the assignment of operations to machines in the parent chromosomes. The role of sequence crossover operators is to produce two new offspring by exchanging partial information of the sequences of operations on machines in the parent chromosomes. The assignment and sequence crossover operators developed in the proposed genetic algorithm are: OMAC (operation to machine assignment crossover), JLOSC (job-level operation sequence crossover), and SLOSC (sub-lot-level operation sequence crossover).
From a given pair of parent chromosomes, OMAC creates two child chromosomes, preserving the sequence of operations in each parent chromosome and modifying the machines' affectation. The creation of a child chromosome using an OMAC crossover operator is illustrated in Figure 4. We assume that in the first step, the operations located in genes 1, 3, 4, 7, 12, 13, 14, and 15 are selected from parent 1 to change their assignments. Then, the second step is to copy all the genetic information of parent 1 to the offspring, except for the assignment information of the selected operations. The final step is to obtain the assignment information of the selected operations from parent 2 and then copy them to the corresponding genes in the offspring to produce a new child chromosome. The same process is repeated to create child 2 by starting to select the operations from the chromosome parent 2. Three versions of OMAC were developed to correspond to different ways of selecting the operations for which the assignment will be modified: The objective of the crossover operator OMAC 2 is to balance the workload between the machines, whereas the goal of the crossover operator OMAC 3 is to reduce the ending dates of jobs with longer completion times. The crossover sequencing operators JLOSC and SLOSC produce two new offspring by exchanging the sequencing information of parent chromosomes while keeping the assignment of machines in a gene unmodified. They are applied with given probabilities. The creation of the child chromosome by JLOSC and SLOSC with the preservation of the assignment information of parent chromosomes is illustrated in Figures 5 and 6. In the JLOSC crossover operator, the first step is to select a set of operations from parent 1. The next step is to copy all of the genetic information of the genes that are associated with the chosen operations to the offspring. The final step is to fill the empty genes in the offspring chromosome, while the machine assignment information of the empty genes is preserved from parent 1. This is done by copying the remaining operations from parent 2 to the empty genes, while the operations retain their appearance order from the second parent. The same procedure is repeated to create child 2 by starting to select the operations from parent 2. SLOSC is identical to JLOSC, with a difference in the second step of SLOSC: Only the genetic information of the chosen genes with the same sub-lot and job indexes is copied to the offspring. Two operators of each kind were developed (JLOSC1, JLOSC2, SLOSC1, and SLOSC2): The SLOSC1 operator is similar to the JLOSC1 operator, but instead of keeping the unmodified operations of all the sub-lots, only the information from the sub-lot of the selected operation is kept and copied to the child chromosome.
The SLOSC2 crossover operator is similar to the JLOSC2 crossover operator except that in the first step, a set of sub-lots of jobs with smaller completion times is chosen. The goal of the JLOSC2 and SLOSC2 crossover operators is to keep the sequencing of operations of jobs and the sequencing of operations of sub-lots of jobs, respectively, with smaller completion times.  It is important to note that, after the application of any crossover operator, in the child chromosome, for a given i and l, the gene (i, l, j, k) is located after any gene (i, l, j , k ), where j < j. This ensures that precedence constraints between operations of a particular sub-lot are not violated in the new created child chromosome.

Mutation Operator
The role of mutation operators is to prevent the algorithm from being trapped in a local optimum and to maintain genetic diversity in the population. The mutation operators are usually applied on chromosomes with given probabilities. Assignment and sequence operators also exist amongst the mutation operators. Assignment mutation operators change the assignment of operations to machines without changing the sequencing of these operations, and sequence operators only change the sequencing property of the chromosome undergoing the mutation, while the assignment property is preserved. The mutation operators used in the developed genetic algorithm are: • ROAM (random operation assignment mutation): This operator is applied with a given probability on a set of operations of a given individual chromosome and randomly changes the assignment property of these operations to another machine. • OSSM (operation sequence shift mutation): Whenever OSSM is applied on an individual, a set of operations is selected. Then, these operations are moved to other positions on the chromosome in such a way that no precedence constraint is violated.

Local Search Methods Local Search through Gene Movement
This local search method improves the quality of any solution built after the application of the crossover and mutation operators. It is repeated as long as the quality of solutions obtained is not improved. The procedure for this method is illustrated through an example in Figure 7. The first step is to choose a set of operations of jobs with longer completion times; the goal is to reduce these completion times. In the second step, the procedure performs a search in the neighborhood of a solution by changing the positions of the operations chosen on the machines in order to improve the quality of this solution. This process is repeated until the iteration criteria are met. The steps of this local search method are as follows:

•
Step 01: Choose a set of operations of jobs with larger completion times in the current chromosome.

•
Step 02: Each operation chosen in Step 01 is positioned just before the previous operation assigned to the same machine (left movement) or just after (right movement) the following operation assigned to the same machine, while respecting the constraints of precedence between operations of sub-lots of jobs. This process allows one to lead research in the neighborhood of a solution by changing the sequencing of operations on the machines.

•
Step 03: This process is repeated until the maximum number of iterations is not reached.

Local Search by Grouping Sub-Lots
This local search method consists of grouping the operations of sub-lots of jobs on machines if the capacities of the machines allow the processing of several operations at the same time. This method can only be applied for machines such as ovens and cooling cells. The procedure for this local search method is illustrated through an example in Figure 8. For each gene, the sub-lot grouping procedure consists of browsing the chromosome from left to right. If the sub-lots of the same job as the considered gene are assigned to the same machine and the latter can process several sub-lots at the same time, all the genes of these sub-lots are repositioned one next to the other in the chromosome, which means that these sub-lots are treated at the same time by the same machine. The grouping of operations of sub-lots of jobs on machines must take into account the precedence between operations of these sub-lots.

Local Search through Intelligent Assignment of Operations
In this local search approach, a set of operations assigned to the most loaded machines is reassigned to the least loaded compatible ones. Table 2 gives the common and different elements between the two developed algorithms. In this table, the heuristics for the generation of the initial population, crossover, and mutation operators of the two genetic algorithms are given. The two algorithms differ in the generation methods of the initial population and in the crossover operators. In both algorithms, several affectation and sequence heuristics were tested for the generation of the initial population. In Table 2, the heuristics allowing one to obtain the best solutions are given. The crossover operators of the first genetic algorithm are based on random choices of operations, while in the second algorithm, different operators that are specific and more adapted to the studied problem have been developed. Regarding the mutation operators, they are common between the two algorithms. The pseudocode of the generic genetic algorithm is given in Algorithm 2. Table 2. Common and different elements between the two genetic algorithms.

Iterated Local Search Algorithm
This section describes an iterated local search (ILS) algorithm for solving the studied problem. ILS is a simple, robust, and highly effective local search procedure that explores local optima in a given neighborhood (Lourenço et al. [41]).
The iterated local search algorithm starts from an initial solution and obtains a local optimum in its neighborhood through a local search procedure. To improve upon the current local optimum, ILS applies perturbation operators to this solution in order to move away from the current local optimum. An acceptance criterion is employed to determine which local optimum will become the current local optimum in the next iteration. The above process is repeated until a termination criterion is satisfied.
To generate the initial solutions of the iterative local search algorithms, the affectation and sequence heuristics for the generation of the initial population of the genetic algorithm are used, except that in the ILS algorithm, only one solution is generated.
The local search procedure is based on two operators: the affectation operator (AO) and sequence operator (SO). The affectation operator changes the assignment of a set of operations to the machines without changing their sequencing on the machines, while the sequence operator changes the sequencing of operations without changing their affectation to to machines.
Three affectation operators, AO1, AO2, and AO3, and three sequence operators, SO1, SO2, and SO3, were developed: • AO1 and SO1: The operations are chosen randomly. • AO2, AO3, SO2, and SO3: The operations are chosen according to rules developed to be specific to the studied problem.
The goal of the AO2 operator is to balance the workload between the machines, while AO3, SO2, and SO3 operators allow the reduction of the ending dates of jobs with larger completion times. The difference between these operators is in the choice of operations to change their assignments and sequences. The first step of the SO operator consists of randomly choosing one of the sequence operators. Then, a set of operations in the current solution to be improved is chosen. The second step is to shift the selected operations to the left until the permutation of these operations with the previous operations assigned to the same machines is effective. If the left shift is blocked (precedence constraints are violated), then, from the solution obtained in the second step, the selected operations are shifted to the right until the permutation of these operations with the next operations assigned to the same machines is effective. This process is repeated until the maximum number of iterations is not reached.
In the AO operator, the first step is to choose one of the assignment operators. Then, a set of operations in the current solution to be improved is chosen. In the second step, the selected operations are assigned to another machine among a set of alternative machines. If the solution obtained after the reassignment of operations is not better than the current solution, then an SO sequence operator is applied on the selected operations with the new machine assignments obtained in the second step. This process is repeated for a given number of iterations.
The local search method can get stuck in a local optimum if the solution space is reduced. To fix this problem, a number of insertion moves to the current local optimum to obtain a perturbation solution are carried out. To determine the appropriate number of perturbations, different strategies have been tested. The best strategy, and that which is used in the algorithm, consists of performing a perturbation if the degree of similarity between the solutions is lower than a given threshold defined beforehand, which ensures a sufficient distance between solutions. The evaluation of the degree of similarity is based on the gaps between solutions. The similarity measures between solutions are evaluated to ensure the stability in the proposed solutions for the production scheduling.
The different steps of the iterative local search method in one iteration are illustrated through an example in Figure 9. The algorithm begins by generating an initial solution using one of the assignments and sequencing heuristics. Then, one of the SO sequence operators is applied. If no improvement of the solution is obtained after the application of the SO operator, the next step consists of applying an AO assignment operator on the operation chosen beforehand. If the application of the AO operator does not improve the initially generated solution, an SO operator with the new affectation obtained after the application of the AO operator is applied. If the solution is still not improved, a perturbation operator is applied on this solution, and the algorithm moves on to the next loop iteration. This operator consists of applying one of the randomly chosen sequencing or assignment operators. This process is repeated until the maximum number of perturbations is not reached. Figure 9. Steps of the iterative local search algorithm for one iteration. Table 3 gives the common and different elements between the two developed algorithms. In this table, the heuristics for the generation of initial population, affectation, and sequence operators of the two iterated local search algorithms are given. The difference between the two algorithms lies in the heuristics used for the generation of initial solutions. In Table 3, the heuristics allowing one to find the best solutions are given. The affectation and sequence operators of the first iterated local search algorithm are based on random choices of operations, while in the second algorithm, the choice of these operations is based on specific rules more suited to the considered problem. The pseudocode of the generic iterated local search algorithm is given in Algorithm 3.

Algorithm 3:
Iterative local search algorithm 1. Initialization: The initial solution is generated by using the previously described assignment and sequencing heuristics ; 1. current solution ← initial solution ; 2. Evaluation: Evaluation of the initial solution ; 1. SimilarDegree ← 0 ; 3. Application of the SO sequence operator; 4. If the solution is improved in 3, then ; 1. current solution ← solution obtained in 3, else go to 5; 5. Application of the AO affectation operator ; 6. If the solution is improved in 5, then; 1. current solution ← solution obtained in 5, else go to 7; 7. Application of the PO operator on the current solution ; 8. Update SimilarDegree by calculating the distance between the current solution and the improved one; 9. If the degree of similarity between the solutions obtained in 8 is less than a predefined threshold,; then go to 3, else go to 10 ; 10. End of algorithm Table 3. Common and different elements between the two iterated local search algorithms. ILS1  SPT, MNOR  SO1  AO1   ILS2 MMWA, MNOR SO2, SO3 AO2, AO3
The real instances were built after collecting data for some examples of production days. The randomly generated instances of the HCT type were built using data (such as the number of machines, due dates of jobs, and time window of availability of machines) from the real instances of HCT. The other parameters were generated from uniform distributions, such as the parameters of the randomly generated instances. To adapt the instances from the literature to the studied problem, the missing data, such as the quantities of dishes, due dates of jobs, machine capacities, and time windows of availability of machines, were randomly generated according to uniform distributions and by taking into account the data on the processing times of operations and the setup times of machines given in these instances. The results of the developed resolution methods for these different categories of instances are presented below.

Optimization of the Metaheuristics' Parameters
To optimize the parameter values of the developed metaheuristics, the Taguchi method was used. The choice of these parameters has a significant effect on the efficiency of the metaheuristic algorithms. The parameter values that are needed for optimization act like controllable factors in the design of experiments. The aim is to find an optimal combination of the parameters such that the total flow time is minimized. To employ the Taguchi Table 4.

Discussion of the Experimental Results
The optimization criteria of the developed metaheuristics are based on the quality of the solutions and rapidity. The quality and efficiency of the metaheuristics were proven by comparing the solutions obtained with these algorithms and the optimal solutions of the mathematical model. Regarding the rapidity, the computation times of the algorithms were compared with those of the mathematical model. From the results presented above (Tables 5-9), we remark that the two genetic algorithms are efficient in terms of quality and rapidity. The performance of the developed resolution methods depends on the types of instances and their sizes, and it also depends on the choice of heuristics for the generation of initial solutions. The first genetic algorithm is more efficient in terms of quality with the MMWA assignment heuristic and the SPT sequence heuristic, while the second genetic algorithm is more efficient with the MMWA and MNOR heuristics. For the iterated local search algorithms, the first algorithm gives better solutions with the SPT assignment heuristic and the MNOR sequence heuristic, while the second algorithm finds good solutions with the MMWA and MNOR heuristics. Regarding the rapidity, the four metaheuristics are more faster with the RA assignment heuristic and the RS sequence heuristic.
The results of the genetic algorithms for different types of instances (Tables 5-9) show that the developed metaheuristics are efficient relative to the quality of solutions. Table 5 represents the results of the different methods developed on a real instance with 82 jobs, 92 sub-lots, 370 operations, and 29 machines. From these instances, several sub-instances were built by increasing the number of jobs, sub-lots, and operations each time to see from what number of jobs, sub-lots, and operations the model is not able to find solutions. From Table 5, we observe that for the real instances of the HCT, the mathematical model quickly found solutions for the small instances with less than five dishes, five sub-lots, 24 operations, and 29 machines. It took a little longer for the instances with more than six dishes, six sub-lots, 29 operations, and 29 machines. However, for the instances with more than nine jobs, 44 operations, and 29 machines, after more than 3 h of execution, the mathematical model did not find solutions. For these instances, the two genetic algorithms succeeded in finding feasible solutions in very short computation times. For example, for the large instance with 82 jobs, 92 sub-lots, 370 operations, and 29 machines, the first genetic algorithm found a feasible solution that respected all the constraints after 5 min of execution.  Tables 6-8 present the performances of the resolution methods in terms of quality and rapidity for randomly generated instances and a comparison between the solutions obtained by these different methods. From the results presented in these tables, we remark that the genetic algorithms found feasible solutions for the large-size instances in reasonable resolution times.
In Table 9, we observe that for the instances adapted from the literature, the genetic algorithms make it possible to improve the solutions of the instances adapted from the literature. For example, for the instance from Lee et al. [35], the genetic algorithms brought about an improvement of 3.92% over the solutions obtained with the methods proposed by Lee et al. [35].   Comparing the two genetic algorithms in terms of rapidity, the first algorithm is faster than the second algorithm, while comparing them in terms of stability both algorithms are stable for the real instances of HCT and randomly generated instances, while for the randomly generated instances of the HCT type, the second algorithm is more stable than the first algorithm. Comparing the two algorithms in terms of quality of solutions, the performance of the algorithms depends on the types of instances. For some instances, the first algorithm is better than the second algorithm, while for other instances, the second algorithm is more efficient than the first algorithm.
By comparing the results obtained with the genetic algorithms and the iterated local search methods for all the tested instances, we find that the ILSs are worse than the GAs in terms of the quality of solutions obtained. However, in terms of rapidity, the ILS methods are faster compared to the GAs.

Application to an Industrial Case
This work was carried out in collaboration with the hospital center of Troyes (HCT). To effectively meet the needs of patients and to improve working conditions and employees' well-being, the hospital center of Troyes implements important measures to improve its daily efficiency. It is in this context that this study of hospital catering activity optimization takes place.
The HCT is carrying out a project to revise its supply chain, which must, in particular, consider the management of food flows. The contribution of the present work is to determine the best plan for meeting customers' demands in terms of food flows and to propose ways to improve the well-being and working conditions of employees of the catering service of the hospital center of Troyes. The aim is to provide methods and tools for scheduling meal-making processes throughout the day.
The hospital center of Troyes is the main member of the South Champagne Hospitals group. Its logistic network ( Figure 10) is composed of a set of customers whose meal needs are met from a warehouse, which is the central kitchen and a set of suppliers. These customers are hospitals, nursing homes, and psychiatric clinics. The central food production kitchen of the HCT produces, on average, 4800 meals per day. Its maximum daily production capacity is 5000 meals per day. The HCT catering service includes a production unit divided into several sectors, from the reception of raw materials to the dispatch of finished products and a coordination unit for menu management. In this study, we are particularly interested in the sectors from the pre-treatment of raw materials to the stock of finished products ( Figure 11). This figure represents the production areas of the studied system and the different machines available for carrying out the operations of the meal-production process.   Table 10 represents the results of the genetic algorithm for some examples of real production days with a comparison between the real solutions, as these production days were organized and the solutions were proposed by the genetic algorithm. In this table, for each instance, the number of dishes, the number of sub-lots of dishes, the total number of operations, the number of machines available, and the average number of meals produced per day are given. The performance indicators between solutions are based on the total flow time and the gaps between them. From these results, we remark that the gaps between the real solutions and those of the genetic algorithm are very important and significant. The performance of the genetic algorithm for these instances depends on their types and sizes. For example, for the instance with 62 dishes, 68 sub-lots, 218 operations, and 29 machines, we brought about a considerable improvement of 18.72% over the real solution, which shows the quality and performance of the developed algorithms.

Conclusions
The present article deals with the study of a new industrial problem. Different resolution methods for the scheduling of production processes in hospital catering were developed. A mathematical model integrating all the constraints of the studied problem was proposed. This model is an improvement of the standard flexible job shop scheduling problem with sequence-dependent setup times by integrating specific industrial constraints. An extensive study confirming the effectiveness of the developed model is presented. The computational results of the mathematical model for different types of instances show the limits of an exact resolution for the problem of scheduling production processes. To solve the large instances of the addressed problem, different metaheuristics were developed and tested on several types of instances. The computational results of these metaheuristics have proven their effectiveness and reliability for the scheduling of operations of the food production process and allowed significant improvements in the real current organization and system performance. As regards future research, this work is limited by the number of human resources available and does not take into account the possible hazards, such as the absence of staff, machinery breakdowns, and the unavailability of raw materials. Further research can be extended by considering these constraints. The present work also opens the way to other perspectives, such as the study of the production planning problem over several days, and our future work will focus on the development of resolution methods for this problem.

Data Availability Statement:
The data presented in this study are available on request from the corresponding author. The data are not publicly available due to confidentiality of the company.