Simulation Optimization for the Multihoist Scheduling Problem

Although the Multihoist Scheduling Problem (MHSP) can be detailed as a job-shop configuration, the MHSP has additional constraints. Such constraints increase the difficulty and complexity of the schedule. Operation conditions in chemical processes are certainly different from other types of processes. +erefore, in order to model the real-world environment on a chemical production process, a simulation model is built and it emulates the feasibility requirements of such a production system. +e results of the model, i.e., the makespan and the workload of the most loaded tank, are necessary for providing insights about which schedule on the shop floor should be implemented. A new biobjective optimization method is proposed, and it uses the results mentioned above in order to build new scenarios for the MHSP and to solve the aforementioned conflicting objectives. Various numerical experiments are shown to illustrate the performance of this new experimental technique, i.e., the simulation optimization approach. Based on the results, the proposed scheme tackles the inconvenience of the metaheuristics, i.e., lack of diversity of the solutions and poor ability of exploitation. In addition, the optimization approach is able to identify the best solutions by a distance-based ranking model and the solutions located in the first Pareto-front layer contributes to improve the search process of the aforementioned scheme, against other algorithms used in the comparison.


Introduction
Some chemical-production systems use tanks containing treatment baths in order to generate finished products.
ose treatment baths can contain rinsing, acid, or electroplating solutions.
is kind of production system is necessary when a product needs to be treated with a specific treatment bath to enhance its mechanical, electrical, or esthetic properties. e products are soaked in each tank according to a given sequence [1]. Commonly, these kinds of products are loaded on carriers or rack baskets.
As other manufacturing facilities, the handling material of the chemical-production systems is executed by trackmounted hoists, i.e., the rack baskets are moved from facility to facility by one or more hoists. e hoists can be (or not) on the same track. erefore, the hoists are used as means of transport between tanks [2]. e hoists move rack baskets from facility to facility according to the production schedule. It means that the hoists should be available, as an assumption, in order to not compromise the already defined production schedule. However, the hoists frequently operate on shared tracks.
erefore, the hoists also should be scheduled to avoid colliding with each other. Normally, the production schedule is considered as input data for the hoists' schedule. It means that these schedules are usually planned separately. It can be considered a risk to enhance the performance of the chemical-production system. In this case, any production schedule should consider the hoists' availability. e simultaneous schedule of jobs and hoists assumes a central relevance for surface treatments. e treatment processes can be modeled as a job-shop production.
Consequently, a strong need for enhanced hoist scheduling systems emerged, representing a key enabler to maximize productivity and product quality [3].
A classic job-shop configuration can be detailed based on the aforementioned chemical processes [4]. Let a set of n jobs (products) J � J 1 , J 2 , . . . , J n on a set of m machines (tanks) M � M 1 , M 2 , . . . , M m . Each job J i requires a series of p i operations (treatments) O � O i1 , O i2 , . . . , O ip i with precedence constraints. As any job-shop configuration, each machine can process only one job at one time without preemption, and each job can be processed on one and only one machine at one time. However, compared with the classical job-shop configuration, the multihoist scheduling problem (MHSP) includes additional constraints. Such constraints increase the difficulty and complexity of the schedule. Constraints are outlined below: (i) e jobs require a hoist in order to move them between facilities. erefore, buffers are prohibited in the MHSP. e limitations mentioned above make a difference between the classical job-shop configurations and the MHSP. It is a strongly constrained problem, known as an NP-hard one [5]. e problem is indeed NP-hard, and each extension involves an increase of the complexity of both the model and its resolution. is is the reason why we have decided to work on a new approach.
In addition, some industries such as pharmaceutical, chemical and plastic have similar operational conditions, i.e., no-wait in the processes. erefore, building schedules with the no-wait component should be more suitable for modeling in those industries, and also other problems such as train scheduling, aircraft landing scheduling, surgery scheduling, among others [6].
Based on the previous constraints, the assumption on identical jobs, in the simultaneous schedule of jobs and hoists, should be omitted. In addition, this research does not have the objective to find a repetitive hoist scheduling.
ere exist multiple job types and diverse hoists on the chemical-treatment processes. Moreover, this research does not have as objective to preassign to each hoist a precise set of tanks. Although this simplification tackles the collision constraints, it restricts the sequence of treatment of the jobs to follow the layout of the tanks. If the requirements are low-demand, the utilization rate of the tanks would decrease.
In order to identify the difficulty and complexity of the MSHP, an example is provided below.
Assume a schedule with four jobs J � J 1 , J 2 , J 3 , J 4 . Every job J i is formed by a sequence of three operations O i, 1 ; O i,2 ; O i, 3 performed one after another. Table 1 details the alternative machines for every operation and additional data.
Let a sequence be J 1 ≺J 2 ≺J 3 with O 11 in M 1 , O 21 in M 4 , and O 31 in M 2 . Assume the minimum processing time window for all the operations. Figure 1 details the sequence step by step by the hoist. Let the next set of operations be O 12 in M 2 , O 22 in M 4 , and O 32 in M 1 . According to the previous sequence, the hoist has to break free M 2 with O 31 , and break free M 1 with O 11 . However, unlike a normal job-shop configuration, a blocking occurs in the MHSP. e simultaneous schedule of jobs and hoists is a necessity to avoid blocking through the treatment processes. To find a suitable sequence and appropriate assignment of the operations to the machines in order to avoid potential blockings and nowaits in the processes is the problem statement.
In order to simultaneously schedule jobs and hoists, this study proposes a solution for the MHSP. is paper details how a solution can be built by three decisions, i.e., operation scheduling, tank assignment decision, and hoist assignment. In particular, let us discuss vectors, such as permutations, and these vectors can represent the processing sequence of operations.
In this research, the processing sequence of operations is considered as a permutation. Table 2 shows a vector (ranking) of five operations as an example.
In this sense, the members of the population of the proposed algorithm are permutations of elements. Table 3 depicts some permutations, where each element represents an operation, and the processing sequence of the operations is executed according to each permutation.
As previous research, this research also tries to use the space of permutations as solutions. e proposal of this research is to use a probability distribution based on permutations.
For the aforementioned proposal, an Estimation of Distribution Algorithm (EDA) is considered for the MHSP. According to the literature review, the EDA has been scarcely studied for solving the MHSP. EDA is an experimental technique that belongs to the evolutionary computation field. Instead of using traditional evolutionary operators such as crossing and mutation, the EDA produces solutions through a probability model based on the previous solutions. e probability model mentioned above is built with statistical information, i.e., based on the search experience. e EDA makes use of the aforementioned probability model to describe the distribution of the solution space. In this research, the distribution of the permutations, for the MHSP, is the solution space. However, there is no probability model that can produce a feasible solution on a permutation-based representation [7]. Any EDA needs to be reconfigured to tackle permutation-based problems. en, as hypothesis, using a specific probability model for this issue should be effective against other models.
In order to produce feasible solutions from the proposed EDA, this study uses the Mallows model. Mallows [8] establishes the Mallows exponential model. e Mallows model incorporates to each permutation σ, i.e., each solution, a probability. In order to compute the probability mentioned above, the model uses the distance concept between σ and a central permutation σ 0 . A distance metric is computed between permutations, by algebra of permutations. e aforementioned probability will be smaller if the distance between σ and σ 0 is bigger. Figure 2 details an example of such exponential distribution, with five vectors.
en, the Mallows model is considered as the specific probability model for the distribution of the permutations. e search process is based on the exponential model, and any traditional evolutionary operator is omitted. Based on Figure 2, the exponential shape of the model is built using the distance between all the permutations σ n of the solution space and the central permutation σ 0 . Traditionally, the central permutation is estimated or selected between all the members of the population. e best solution, during the evolutionary progress, can be considered as the central permutation. However, the best solution is not a guarantee to be the best estimation. erefore, finding the best estimation should be reached in order to get the best scheduling solution for the MHSP.
Once the central permutation is estimated or selected, generating new offspring, using the exponential model, is the main objective. en, the main idea is to use the probability between permutations to generate a new offspring. However, the Mallows model is not able to generate offspring by itself. e reason is because there can be many permutations with the same distance to the central permutation. It creates confusion when choosing the offspring.
It is solved through the decomposition of the distance between each permutation and the central permutation.
e Generalized Mallows Distribution (GMD) process, detailed in Fligner and Verducci [9], and Fligner and Verducci [10], explains the procedure to get the decomposition of the distance between permutations and thereby generate a new offspring. With this strategy, the inconvenience of probability models in permutation-based problems is solved. e GMD process permits to produce feasible solutions.
Based on the concept of properly modeling the main variables that intervene in the performance of the process has been a priority in the solution of real-world scheduling problems [11]. Such variables or characteristics could be inside or outside of the shop floor and these should be incorporated to efficiently solve the scheduling problem. erefore, this research appends the real environment on the chemical production process mentioned above. Moreover, if the real environment is considered, theoretical assumptions are not necessary to incorporate in the MHSP.
en, a simulation model is built, and it emulates the aforementioned chemical production system. Different resources, in the chemical production process, are considered in the simulation model mentioned above. ose resources perform all the operations required by any job that is scheduled. e aforementioned simulation model, built on Delmia-Quest® platform, included vast details that the chemical production process presents, i.e., the process itself. e results of the simulation model, i.e., the makespan and the workload of the most loaded machine, are necessary for providing insights about which schedule on the shop floor should be implemented. e proposed optimization method uses the results mentioned above in order to build new scenarios for the MHSP and to solve the aforementioned conflicting objectives. e research motivation is to show how the exponential distribution between permutations helps to reduce the deficiencies of the EDA. In addition, it is preferable that the methods used should contain some well-defined math expressions. It helps to understand how the solutions are generated. erefore, the research motivation is to use new methods to explicitly establish the search process of new solutions. e main reason, for doing this research, was to determine what is most useful, i.e., the recent algorithms or the use of exponential distributions to obtain the same or better solutions.
Although diverse methods and strategies have been used to solve the MHSP, this paper contributes to the state of the art as follows:

Mathematical Programming.
Scheduling the movements of a hoist in electroplating facilities is tackled by Varnier et al. [1]. Based on their research, in most works, for the mono-product case, the layout of the shop is generally considered as a fixed data. However, the results of their research show that the layout should be considered in order to improve the productivity of the shop. e authors propose a combination of scheduling and layout design.
Other works, such as Grunder et al. [12], try to identify the interaction between the layout of the shop and the productivity of a treatment line. e authors show that a particular type of configuration of the shop can be maximizing the performance of the line. However, it is not functional for all the studied cases. In addition, the authors propose a branch and bound algorithm to process the optimal layout of a saturated single-hoist production line.
Other types of interaction can be found in Subai et al. [13], where environmental constraints are incorporated in the schedule of a treatment surface process. e authors focus to maximize the shop throughput by means of mathematical models. Such models include a nonlinear global cost function where the environmental cost plays a key role.
A mixed-integer linear programming model for a multiproduct batch plant with nonidentical parallel processing machines is used in Berber et al.'s [14] research. e authors show how to minimize the total production time without using heuristic rules. Diverse numerical instances and one industrial problem are used to test their proposed model. e aforementioned model is able to produce better solutions for the industrial example considered.
An exact solving method is proposed by El Amraoui et al. [15]. e mentioned method is a linear optimization approach.
is work is implemented to optimize the cycle length and the throughput rate in a cyclic hoist scheduling problem where there are r different part-jobs in an electroplating line.
Other mixed-integer linear programming model for multirecipe and multistage material handling processes is developed in Zhao et al. [16]. In this research, the authors simultaneously consider the production line arrangement, and the customized production ratio. Various case studies are used to demonstrate the efficacy of this methodology.
For the two-hoist cyclic scheduling problem, Chtourou and Manier [17] propose a mixed integer linear programming model. e authors consider avoiding collision between the hoists which share a common track whilst the cycle time is minimizing.
For the dynamic hoist scheduling problem with multicapacity machines, Feng et al. [18] detail a mixed-integer programming model. e mentioned model considers jobs randomly arriving throughout the horizon as a dynamic issue.
In a copper production plant, Suominen et al. [19] present a nonlinear optimization and scheduling approach to maximize a smelting furnace production. e proposed approach is presented by simulating the evolution of the process over the optimization horizon.

Applied Computational Intelligence and Soft Computing
For the ethylene cracking process with feedstocks and energy constraints, Su et al. [20] address a scheduling problem by a hybrid mixed-integer nonlinear programming formulation. e problem mentioned above has a special situation, i.e., the facilities require periodic cleanup to restore the performance. e authors make useful suggestions for real cracking process production by means of numerous examples. erefore, the examples mentioned help illustrate the utility of the model.
In an ice cream shop, Wari and Zhu [21] present a mixed-integer linear programming model. e key point of this contribution is that a multiweek production scheduling is considered. Furthermore, the model mentioned above incorporates operational aspects to enhance the solution quality. e aforementioned model is tested and compared with heuristics methods. e results report that the model of this research is able to efficiently and effectively handle the multiweek aspect.
Qu et al. [22] study the effect of integrating scheduling with the optimal design of a two-dimensional production line to maximize production efficiency in the multirecipe and multistage material handling processes. Various case studies are used to demonstrate the efficacy of the integration mentioned above.
Finally, Jianguang et al. [23] consider the cyclic jobshop hoist scheduling with multicapacity reentrant tanks and time-window constraints. As the authors explain, jobs are processed in a series of tanks with a defined processing sequence of operations for all of the jobs. A mixed-integer linear programming model is developed by addressing the time-window constraints and tank capacity constraints.

Heuristics.
For the cyclic hoist scheduling problem considering material as well as resource handling constraints, El Amraoui et al. [24] present a new heuristic where the time windows are maintained for all soaking operations and also overlapping cycles are allowed. As in almost any research, the authors compare the proposed heuristic mentioned above with other existing algorithms to prove its efficiency.
Another example of heuristics for the hoist scheduling problem can be found in Kujawski andŚwiaojtek [2], where the sequence of products is analyzed by changing the ordered items in electroplating production lines in order to reduce the makespan. e items are located in queues where the sequence is monitored.
For the online hoist scheduling problem, Kujawski and Swiątek [25] prepare a set of production scenarios before the real-time system starts. A real-life shop located at Wrocław, Poland, is used to explain carefully the algorithm. In this research, the utilization ratio is considered in order to choose the best schedule.
For the flexible and nonlinear electrochemical processes, Beerbühl et al. [26] propose a heuristic by combining scheduling and capacity planning. e heuristic mentioned above is able to tackle nonconvex and mixed-integer problems, such as the electrolysis of water to produce hydrogen, reformulating these kind of problems into convex and continuous nonlinear problems.
Recently, Laajili et al. [27] detail an adapted variable neighborhood search-based algorithm for the cyclic multihoist design and scheduling problem. e Laajili et al. study considers identical jobs, and therefore identical processing sequence of operations for all the jobs. When identical jobs occur, it permits to create a cyclic multihoist schedule.

Metaheuristics.
For the printed-circuit-board electroplating line, Lim [28] consider to determine the best throughput rate. e authors propose a genetic algorithmbased approach for the cyclic hoist scheduling problem.
rough an experiment, the proposed algorithm is more efficient than the previous mathematical programmingbased algorithms.
In an aluminum casting center, Gravel et al. [29] present an ant colony optimization metaheuristic. e representation of the solution in this research considers different objectives. In addition, the authors implement the proposed method, by introducing software, in the shop casting center.
For the mixed batch and continuous processes, Wang et al. [30] develop a differential evolution algorithm. e solution representation considers capacity constraints. A key characteristic, in the proposed algorithm, is how to compute the crossover probability by using the logistic chaotic map method.
For the single hoist cyclic scheduling problem, El Amraoui et al. [31] consider hard resources and timewindow constraints in their proposed genetic algorithm.
For the multihoist scheduling problem with transportation constraints, Zhang et al. [4] study how to reduce the makespan. e global limitation is that there are no buffers to storage work in process. en, a right assignment is critical in the proposed solution.
e authors detail a modified genetic algorithm to tackle the problem mentioned above. e aforementioned algorithm also uses a modified shifting bottleneck procedure, in order to build a feasible schedule. e results of this research show that the proposed algorithm is able to handle several transport resources.
For other electroplating shops, El Amraoui et al. [32] examine the processing time which is confined within a time window. A genetic algorithm approach is presented to solve the multijobs cyclic hoist scheduling problem with a single transportation resource.
For a single robot and flexible processing times in a robotic flow shop, Lei et al. [33] aim to increase the throughput rate. A hybrid algorithm based on the quantuminspired evolutionary algorithm and genetic operators is presented for solving the cyclic scheduling problem. e algorithm integrates three different decoding strategies to convert quantum individuals into robot move sequences. Besides, crossover and mutation operators with adaptive probabilities are used to increase the population diversity. A repairing procedure is proposed to deal with infeasible individuals. Comparison results on both benchmark and randomly generated instances demonstrate that the proposed algorithm is more effective in solving the studied problem in terms of solution quality and computational time. 6 Applied Computational Intelligence and Soft Computing

Hybrid Approaches.
In a surface treatment system with time window constraints, Chové et al. [34] investigate how to improve throughput rate without loss of treatment quality. e authors propose a new combined approach based on both predictive scheduling and reactive scheduling in an industrial case where only one hoist is used.
Another example of hybrid approach can be found in El Amraoui and Nait-Sidi-Moh [35]. e authors use P-Temporal Petri Net models to describe the behavior of different scenarios of the shop in a specific cyclic hoist scheduling problem. e authors propose a linear programming model to determine the optimal plan, where the exact beginning and ending instants of each task is the output of the model. In addition, a simulation tool validates the hybrid method mentioned above.
For a multicrane scheduling problem, where a set of coils should be transported from a storage location to another side, Xie et al. [36] formulate it as a mixed-integer linear programming model and propose a heuristic algorithm to solve the problem.
Most publications belong to this category, i.e., hybrid approaches, such as the Yan et al.'s [37] research. e authors focus on a bi-objective, i.e., the cycle time and the material handling cost over a cyclic hoist scheduling problem where the parts have different flow patterns. A bi-objective linear programming model is detailed. In addition, a Pareto-front is formulated in this research with respect to the bi-criteria mentioned above. A hybrid discrete differential evolution algorithm computes the Pareto-front. Moreover, the work-inprocess level is used to adjust the exploration and exploitation of the search of the best solution.
Another bi-objective case is found in El Amraoui and Elhafsi [38]; where improving the productivity and quality are studied. Firstly, the authors formulate the problem as a mixed integer linear programming model. In addition, the authors detail an efficient heuristic procedure to obtain the movements of the hoist. e results of the aforementioned heuristic are compared to a lower bound obtained from the model mentioned above and to the best available heuristic in the literature.
For the aerospace and electroplating industries, Basán and Méndez [39] present a hybrid approach using mixedinteger linear programming, heuristics, and a simulation model. e mentioned simulation model contains multiple particularities in a multiproduct multistage production system where a single hoist is analyzed. Real data, of an aircraft manufacturing industry, are used to minimize the operating cost and maximize the productivity of the system.
Another real case is found in Mori and Mahalec [40]. e authors deal with scheduling of the continuous casting of steelmaking. As the authors explain, a mixed-integer linear programming model computationally is intractable. en, the authors produce a production planning by solving a relaxed mixed-integer linear model at the first stage. After that, the authors build schedules by simulated annealing and a shuffled frog-leaping algorithm. As other previous studies, real data are utilized to test the proposed method.
For the scheduling of multicrane operations in an iron and steel enterprise, Xie et al. [41] study how to reduce the makespan, by a mixed-integer linear programming model and a heuristic. e authors identify properties to avoid crane conflicts.
As another study, in a case example of the chemical industry, Hahn and Brandenburg [42] present a linear programming model and an aggregate stochastic queuing network model in order to obtain the best solution. e proposed hybrid approach considers the production-related carbon emission and overtime working hours in the solution.
For a steelmaking-continuous casting manufacturing system, Jiang et al. [43] develop a multiobjective soft scheduling to address the uncertain scheduling problem. ree objectives are tackled, i.e., waiting time, cast-break, and over-waiting. e authors detail a preference-inspired chemical reaction optimization algorithm, and a simulationbased t-test method is used to provide feedback of the solutions. e convergence is handled by a knowledge-based local search embedded to the aforementioned optimization algorithm. Real-world steelmaking-continuous casting instances served as the input parameter to work with the mentioned approach.
Finally, in flexible flow shops, such as the paper mill industry, Zeng et al. [44] construct a multiobject optimization model, where three objectives are analyzed, i.e., makespan, electricity consumption, and material wastage are included in the mentioned model. At the beginning, only two objectives are considered in the solution, i.e., electricity consumption and material wastage. After that, a hybrid nondominated sorting genetic algorithm II method is employed to solve all the objectives. Real-world case study is used in this research.
Other hybrid approaches can found in An et al. [45]; where the authors used simulation and cloud computing, to prevent collision detection, and hoisting path planning, in a three-dimensional hoisting system.
Li et al. [46] present a simulation-based solution for a multicrane-scheduling problem derived from a steelmaking shop. e problem is modeled considering different objectives for the jobs and workload objective for the cranes. e hybrid approach solves the problem by a heuristic. Tamaki et al. [47] propose a simulation-based solution by adopting the metaheuristic methods to solve the cranescheduling problem in manufacturing systems of the jobshop type where semi-products are picked up and delivered by using cranes between the facilities.
Zhang and Oliver [48] include the crane-scheduling problem into the production scheduling environment and combine them together to obtain an integrated schedule. A simulation-based optimization solves this integrated scheduling problem. A genetic algorithm is introduced to determine the allocation of machines and cranes. A simulation model referring to a queuing network is used to evaluate the crane and machine allocation results and provides the fitness value for the genetic algorithm.
Diverse gaps in the current state of the art can be noticed based on the review mentioned above. e main steps in the proposed methods contain greedy procedures to get promising solutions. erefore, the performance of the Applied Computational Intelligence and Soft Computing aforementioned approaches is related to those mentioned procedures. As an example, the genetic algorithms build new solutions by evolutionary operators; however, those operators do not permit to have control in characterizing the solution space explicitly. en, in this research, an explicit probability distribution over the MHSP is offered to characterize the solution space explicitly. In addition, there are currently no publications that use an EDA for the MHSP. It is clear, from the exposed review, that the EDAs have still a gap to improve their performance, and the exponential distributions approach might be a useful way to improve the EDAs. In addition, from the exposed review, it is of interest to determine how much better EDAs can be of the recent algorithms.
Finally, Table 4 depicts the state of the art on the MHSP with other related problems.

Problem Statement
e multihoist scheduling problem shares features and characteristics of a flexible jobshop configuration. Wang et al. [49] and Yan and Wang [50] explain the problem formulation for this configuration. e main constraints are detailed below: (i) For each job, the corresponding operations have to be processed in the given order, that is, the starting time for an operation must not be earlier than the point at which the preceding operation in the sequence of operations of the respective job is completed (ii) Moreover, each operation has to be assigned to exactly one tank (iii) Preemption is not allowed, i.e., each operation must be completed without interruption once it starts (iv) e operations assigned for each tank have to be subsequently established, that is, an operation is only allowed to be assigned to the sequence of a tank if the preceding position on the sequence is already established (v) If the operations i and j are assigned to the same tank k for consecutive positions p − 1 and p, then the starting time of operation j must not be earlier than the completion time of operation i in order to prevent overlapping Additional constraints involved by the MHSP are the following: (i) e processing time in each tank must respect a minimum and maximum limit, i.e., each processing time is bounded, and those limits must be strictly respected to ensure the quality of the products.
(ii) A hoist can perform only one transport operation at a time, and a hoist must have enough time to move between two transport operations.
e mathematical model is based on Pérez-Rodríguez et al. [11] detailed below.
Let J(i) denote the job to which operation i belongs, and let P(i) be the position of operation i in the sequence of operations belonging to job J(i) starting with one, i.e., P(i) � 1 if the operation i is the first operation of a job. Furthermore, the index set I k defined by I k : � i ∈ O|k ∈ M i denotes the indices of operations i ∈ O that can be processed on tank k. Consequently, there are |I k | positions on tank k.
In order to model the assignment of operations to tanks, assignment binary variables x i,k,p for all p k � 1, . . . , |I k |, k � 1, . . . , m, i ∈ O are introduced if x i,k,p � 1 means that the operation i is scheduled for position p on tank k. e processing time of the operation i on tank k is denoted by t i,k . Furthermore, S i is defined as the starting time for operation i. For each job, the corresponding operations have to be processed in the given order, that is, the starting time for an operation must not be earlier than the point at which the preceding operation in the sequence of operations of the respective job is completed. is constraint is imposed simultaneously on all appropriate pairs of operations, aggregated in the set of conjunctions C given by C∶ � (i, j)|i, j ∈ O: J(i) � J(j)∧P(j) � P(i) + 1}. Consequently, the precedence constraints are given by Moreover, each operation has to be assigned to exactly one position, which is ensured by In addition, only one operation can be assigned to each position, due to constraints i∈O x i,k,p ≤ 1, for all p � 1, . . . , I k , k � 1, . . . , m. (3) e positions on each tank have to be subsequently filled, that is, an operation is only allowed to be assigned to a position on a tank if the preceding position is already filled.
In order to interconnect the tank position variables with the starting time variables and to enforce a feasible schedule, nonoverlapping constraints are defined by for all p � 2, . . . , I k , i ≠ j ∈ I k , k � 1, . . . , m. Each operation is bounded within a time window. is condition is ensured by e total time required to conclude all the operations scheduled, C max , is defined by the constraints To minimize, the C max is given by To ensure compliance with these constraints for each hoist and avoid any collision between hoists that share the same track, the Delmia-Quest® simulation language is preferred in this research. Delmia-Quest® avoids any collision between hoists in each simulation run. In addition, the simulation language is able to order each couple of hoist operations. A controller, in the simulation model, establishes the order of movements, in each simulation run, to satisfy the aforementioned constraints. e approach taken in this study combines the key advantages of both MHEDA and event-discrete simulation. Using the Delmia-Quest® simulation language, the main constraints related to the hoists are satisfied. Furthermore, considering Delmia-Quest®, the makespan and the workload of the most loaded machine are obtained directly from the simulation model, and the MHEDA is in charge of modeling the solution space distribution.

MHEDA for the MHSP
To understand the MHEDA methodology, a multihoist scenario is detailed below. Table 5 shows the jobs, operations, machines, and additional information about the multihoist layout configuration. Figure 3 depicts the layout.
MHEDA contains differences and similarities with respect to a recent published algorithm, i.e., the MEDA (Mallows Estimation of Distribution Algorithm). e MEDA algorithm mentioned above can be consulted in Pérez-Rodríguez and Hernández-Aguirre [51]. e details of the MHEDA and also the differences and similarities between the algorithms, MHEDA and MEDA, are detailed below.

Solution Representation.
In this research, three vectors are used to represent a solution, i.e., a task sequence vector, a machine assignment vector, and a hoist assignment vector. e first vector, the task sequence vector, is a permutation-based representation. e task sequence vector models the processing sequence of operations. An example is depicted below; based on the Table 3, there are three jobs and ten operations, indexed from 1 to 10.
A task sequence vector example is where the operation number four should be executed at the beginning; after that, the operation number one, the operation number seven, and so on. Although each element in the task sequence vector shown above can be located in any position along the solution vector, as any permutation-based representation, the precedence between operations of each job must be kept in the vector. For example, the operations 1 and 2, from job number one, must appear in the same order in the vector, from left to right. If the vector solution of the operation sequence satisfies the precedence mentioned above, it satisfies the corresponding original operation sequence from the job mentioned. is representation is based on Gen et al. [52]. As a fixed parameter, 1000 solution vectors are defined for the population. e second vector, the machine assignment vector: its length equals the total number of operations, where each element represents the corresponding selected machine for each operation. To explain the representation, an example is provided by considering the task sequence vector shown below: en, a feasible machine assignment vector can be Based on Table 3, it means that the machine number four must be used for the operation number one, the machine number three for the operation number two, and so on. As a fixed parameter, 1000 solution vectors are defined for the population.
e third vector, the hoist assignment vector: its length equals the total number of operations, where each element represents the corresponding selected hoist for each operation. Based on Table 3, an example is depicted, by considering the task sequence vector shown below: with ten operations, and two hoists, indexed from 1 to 2, a feasible hoist assignment vector can be where the first operation should be executed with the hoist number two, after that, the second operation, with the hoist number one, and so on. As a fixed parameter, 1000 solution vectors are defined for the population.
By using this representation, i.e., through three vectors, infeasible individuals are not generated.

Fitness
Computing. Two objectives to optimize are considered in this research, i.e., the makespan and the maximum workload. For each solution, the corresponding values are obtained from the simulation model, built on Delmia-Quest®. e main details of the simulation model are outlined below: (i) A conveyor system is used to move rack baskets through different hoists (ii) e conveyor system is only unidirectional (iii) When a rack has finished its sequence production, it goes through the conveyor system to exist the system (iv) Different hoists give service to any rack according to predefined sequence (v) Racks can receive service by different hoists based on the predefined sequence (vi) Putting racks in tanks is only possible by hoists (vii) Hoists are also used to put racks on the conveyor With all these features, the simulation model is able to integrate operation times and workflows. Finally, the fitness is used by the MHEDA to build the Pareto-front and to obtain the best solution at the end of the execution.

Pareto-Front.
All members of the population are used in order to build a Pareto-front based on Kacem et al.'s [53] research. Once a Pareto-front is built, the selection process of the best-candidate solutions is based on where the candidates are located on the Pareto-front. Only the members located in the first Pareto-front layer are preferable. Although the MEDA utilizes the Pareto-front for the selection process as the MHEDA, the MEDA requires a tournament process to select the corresponding candidates. In the MHEDA, the candidates are selected without tournament, i.e., the selected candidates must be located in the first Pareto-front layer. Figure 4 depicts a Pareto-optimality approach example. A nondominated set of solutions in each generation is found and used for building the probability model.

Probability Model for Task Sequence Vectors.
As in the MEDA, the MHEDA builds a probability model by the Mallows model for the selected task sequence vectors. Fligner and Verducci [9] formally establish the Mallows model as where θ is a spread parameter, D(σ, σ 0 ) is the distance from σ to the central permutation σ 0 , and ψ(θ) is a normalization constant. In the present work, the Kendall's τ is the distance metric with which the Mallows model is coupled. e Mallows model was initially proposed by Mallows [8] and later improved by Fligner and Verducci [9] through the generalized Mallows distribution (GMD). e GMD is given as follows: where θ j are dispersion parameters, ψ(θ) is a normalization constant, and V j (σ, σ 0 ) can be defined as an auxiliary vector, i.e., it represents the number of positions on the right side of j with values smaller than the current position in the permutation (σ, σ 0 ). Consider a task sequence vector with four tasks, n equals four, as an example.
Let a central permutation be given by σ 0 � 1, 2, 3, 4 { }, i.e., task 1 is located in the first position, task 2 is located in the second position, and so on.
Let a task sequence vector be given by σ � 4, 2, 3, 1 { }. erefore, V j (σ, σ 0 ) for the jth position is calculated as shown below: Applied Computational Intelligence and Soft Computing e corresponding Kendall's τ distance is five, i.e., D t � n−1 j�1 V j (σ, σ 0 ) � 5, i.e., 3 + 1 + 1 � 5. e procedure is carried out for each task sequence vector, i.e., V j (σ, σ 0 ) should be calculated for all the M vectors in the selected population. e next step consists of computing the dispersion parameters θ j , which is given by where V j � N i�1 V j (σ, σ 0 ). Equation (11) can be solved by the Newton-Raphson method. e probability distribution of the random variables in the auxiliary vector V j (σ, σ 0 ) can be written as is means that the possible values for V j (σ, σ 0 ) in the jth position are located between 0 and n − j, where n is the length of the selected task sequence vectors.
For the operation scheduling decision, the offspring are obtained as follows. As an example, let V s (σ, σ 0 ) a sample vector obtained from equation (12). Specifically.

Probability Model for the Machine Assignment Vectors.
As in the MEDA, the MHEDA builds a probability model to determine an estimate of a distribution model to generate new offspring (machine assignment vectors) using the selected members. Again, as in the MEDA, the MHEDA obtains the estimation by the Univariate Marginal Distribution Algorithm (UMDA). e probability model for the selected machine assignment vectors selected previously can be represented by a probability matrix p, i.e., each p i value represents the amount of times where the machine i is elected for a specific position. Each element of the probability matrix p represents the probability that a position is processed on a machine. e value of each element indicates the rationality of a position processed on a certain machine.
Based on Table 3, and as a short example, the vectors shown below are used to build the corresponding probability matrix p for the first position, i.e., p 1 (Table 6). e offspring are obtained as follows: for each position, generate a U[0, 1] value. en, the value mentioned above is interpolated in the cumulative probability matrix p to identify which machine should be selected.

Probability Model for the Hoist Assignment Vectors.
e MHEDA uses the UMDA algorithm to determine an estimate of a distribution model to generate new offspring (hoist assignment vectors) using the selected members. e probability model for the hoist assignment vectors selected previously can be represented by a probability matrix, q, i.e., each q i value represents the amount of times where the hoist i is observed for a specific task. Each element of the probability matrix q represents the probability that a task is processed on a hoist. e value of each element indicates the rationality of a task processed on a certain hoist.
Again, as a short example, the vectors shown below are used to build the corresponding probability matrix q for the first task, i.e., q 1 (Table 7). erefore, the offspring are obtained with the same procedure as that of machine assignment vectors, i.e., for each task, generate a U[0, 1] value. en, the value mentioned above is interpolated in the cumulative probability matrix, q, to identify which hoist should be selected.
In each generation, new candidate solutions are used for building the Pareto-front, and the probability model is updated with the information from the new members located in the first Pareto-front layer.

Multihoist Simulation Model.
e simulation model is built on Delmia-Quest®. is model includes several types of details that the multihoist process presents: set elements, set hoists and tracks, set jobs, set processes, load and unload processing, and transferring of jobs between tanks. ese situations are present in the given process. e model is able to handle any operation of each job that is scheduled. e model is able to identify the sequences of operations for each job. Figure 5 shows a global 3D layout of the production process. e reason to utilize a discreteevent simulation platform is due to the stochasticity of the underlying process.
e main procedure to build the simulation model, with the features mentioned above, is detailed below, and the model is elaborated using batch commands provided by Delmia-Quest®.

Conveyor Structure.
e aforementioned process contains a conveyor system. It is used to transport any basket (job) through the shop floor. e conveyor system is built using 3D linear and arc segments provided by the platform. Each conveyor segment is positioned and connected according to the flow and the real distribution of the shop floor.

Hoists.
Each hoist is used to execute load and unload operations between tanks (machines). Each hoist attends to a specific group of tanks. If any basket requires a specific operation in a specific machine, then the assigned hoist (to that machine) executes the transport task. e hoist's movements, such as forward, return, park, load, and unload, are executed or controlled by logic commands provided by the platform. e aforementioned logic commands are also able to avoid collision between hoists.

Source and Sink.
A source is built in the model. It enables the baskets to enter the model. Also, a sink is built in the model to destroy all the baskets after the production process.

Buffers.
e buffers are used as load and unload locations. Once a basket enters the shop floor, it goes through the conveyor system and waits on the corresponding buffer until the hoist executes the movement.

Tanks.
e tanks are actually machines. ese are built and positioned on the shop floor according to the layout of the process.

Baskets.
e baskets are actually jobs. ese are created by the source. ese baskets go by different routes in the production process according to the requirements.

Setting Processes.
A process is an operation in the MHSP. A process is executed by a previously established tank. However, the process must be detailed, i.e., indicating which basket should be processed, and the processing time window of the process. Once a process is established, it should be linked to the corresponding tank.

Setting Process Sequence.
A process sequence must be established for each basket. is allows the model to verify if the basket has finished all its operations in the production process. If that is not the case, the basket continues in the shop floor by the conveyor system until all its operations have been concluded.

Verification and Validation.
e simulation model is run under different conditions to determine if its computer programming and implementation are correct as an application in verification technique, which is known as the fixed values test [55]. e throughput, as model result, is verified against data provided by the managers of the process. Figure 6 depicts previous descriptions of three months of real production.
Furthermore, the validation of the simulation model is realized statistically. A comparison of the results derived by the simulation model with real production is done under the same initial conditions, satisfying statistical assumptions in the validation. Figure 7 depicts the information below.    Finally, Table 8 depicts the main differences and similarities between the MHEDA and the MEDA for clarify.

Results and Comparison
In order to validate the relevance of this paper, a comparison of the MHEDA results with others is done. A set of standard benchmarking datasets is used for comparison. e Adams et al. [56] instances; the Fisher and ompson [57] instances; the Lawrence [58] instances; the Applegate and Cook [59] instances; the Storer et al. [60] instances; and the Yamada and Nakano [61] instances.
For each instance, 30 trials are executed to account for the stochastic nature of the MHEDA.
ree metrics are used to compare the performance of the algorithms. First, the mean absolute error (MAE) is where c i is the best hyper volume, from the Pareto-front, obtained after running each trial, and c + is the best hyper e MSE measures the amount of error between two datasets, that is, between the values that the algorithm returns and the values that it should obtain.
ird, the relative percentage increase (RPI) is e RPI is used to compare two quantities while taking into account the "sizes" of the things being compared. e comparison is expressed as a ratio.

Comparison with Other Estimations of Distribution
Algorithms. As other previous works, some EDAs have been included as a benchmark for comparison with the MHEDA scheme; the MIMIC by De Bonet et al. [62]; the COMIT by Baluja and Davies [63], and the BOA by Pelikan et al. [64]. Figure 8 indicates the output by the algorithms using equation (13).
rough box and whisper charts, the dispersion of the values obtained, using the MAE metric, is appreciated. e dispersion of the values, using the MAE, shows the results over the instances and over different runs. As we can see, the MHEDA obtains better results than the other algorithms. Based on the MAE, the MHEDA is more accurate with respect to the other algorithms. Figure 9 shows the output by the algorithms using equation (14). e dispersion of the results is similar to Figure 8. Based on the MSE metric, the MHEDA obtains the interval with the smallest error with respect to the other algorithms for the MHSP. Figure 10 depicts the results obtained by the algorithms using the equation (15). e MHEDA obtains better results FitD t−1 ⟵ Evaluate individuals(fitness)through Delmia Quest® Pareto t−1 ⟵ Select best individuals from FitD t−1 σ 0 ⟵ Central permutation computing from Pareto t−1 V t−1 ⟵ Distance computing from D t−1 and σ 0 ∅ ⟵ Spread parameter computing from V t−1 Ds t ⟵ Sampling from V t−1 p t−1 ⟵ p matrix computing from Pareto t−1 Dp t ⟵ Sampling from p t−1 q t−1 ⟵ q matrix computing from Pareto t−1 Dq t ⟵ Sampling from q t−1 D t ⟵ Replacement all old members with new offspring t: � t + 1 Until stopping criterion is met ALGORITHM 1: Pseudocode MHEDA framework.  than the other algorithms.
e MHEDA scheme outperforms all the previous results.
Based on the results, the members located in the first Pareto-front layer contribute to improve the search process of the MHEDA scheme, against other algorithms.
Although the performance of all the algorithms used in the comparison is outstanding, the MHEDA scheme can find the best value in all the trials. In addition, the MHEDA scheme is able to find the best hyper volume for all the instances used in the comparative. It always helps to find the best value in all the trails. e dispersion of MHEDA is less than other algorithms; it means that the solutions found by the MHEDA are more concentrated around the best value, than other algorithms, i.e., the average of solutions of the MHEDA converges better than other approaches to the best found value.

Comparison with Other Multiobjective Algorithms.
Furthermore, as other previous works, two multiobjective algorithms have been considered to evaluate the MHEDA performance, such as the NSGA by Srinivas and Deb [65] and the NSGA-II by Deb et al. [66]. e experiments are executed in the same computer and language specification. Figure 11 details the output for the algorithms using equation (13). e dispersion of the values, using the MAE metric, shows the results over the instances and over different runs. rough box and whisper charts, it is possible to identify that the behavior of the algorithms is similar than the output detailed above. Based on the MAE metric, the MHEDA is more accurate with respect to the other algorithms. Figure 12 indicates the output by the algorithms using equation (14). e dispersion of the results is similar to Figure 11. Based on the MSE metric, the MHEDA obtains the smallest median error with respect to the other algorithms for the MHSP. Figure 13 shows the results obtained for the algorithms using equation (15). e behavior is practically the same between the algorithms, i.e., the MHEDA scheme again outperforms all the previous results. Figure 13 includes the performance of the three algorithms: the NSGA, the NSGA-II, and the MHEDA after running all the instances. Based on equation (15), the MHEDA scheme outperforms all the algorithms used in the comparative. As we can see, the MHEDA is competitive in order to identify the best-candidate solutions. e performances of the algorithms used in the comparative are very similar than the previous comparison. e medians are 0.25% above the best found value. e MHEDA scheme again found the best value in all the trials. e MHEDA scheme is consistently able to find the best hyper volume for all the instances used in the comparative. Again, it

16
Applied Computational Intelligence and Soft Computing always helps to find the best value in all the trails. e dispersion of MHEDA is much less than other algorithms; it means that the solutions found by the MHEDA are more concentrated around the best value, than other algorithms, i.e., the average of solutions of the MHEDA converges better than other approaches to the best found value. Based on the results, the MHEDA scheme is able to identify the best solutions by a distance-based ranking model, i.e., the Mallows model.

Comparison with Recent Algorithms for the MHSP.
Based on the previous results, recent algorithms are proposed as a benchmark for comparison with the MHEDA scheme.
e recent algorithms mentioned above are the mixed-integer programming model with heuristic presented by El Amraoui and Elhafsi [38]; the algorithm proposed by Xie et al. [41]; and the mathematical model with genetic algorithm detailed by Zeng et al. [44]. All the algorithms mentioned are considered recent algorithms for the MHSP. ese algorithms have been implemented by the authors. e experiments are executed with the same parameters and specifications detailed above. Figure 14 shows the output for the algorithms using equation (13). e dispersion of the values, using the MAE metric, depicts the results over the instances and over different runs. rough box and whisper charts, it is possible to identify that the behavior of all the algorithms are similar. Based on the MAE metric, the MHEDA is more accurate with respect to the other algorithms and outperforms all the algorithms used in the comparison. Figure 15 indicates the output by the algorithms using equation (14). e dispersion of the results is similar to Figure 14. Based on the MSE metric, the MHEDA obtains the smallest median error with respect to the other algorithms for the MHSP. Figure 16 presents the results obtained for the algorithms using equation (15). In this case, the MHEDA outperforms all the other recent algorithms. Practically, there exists significant difference based on Figure 17. e performance of recent algorithms and the MHEDA scheme is different.
Based on Figure 16, the medians are 0.10% above the best found value. e MHEDA scheme again can find the best value in all the trials. e MHEDA scheme is consistently able to find the best hyper volume for all the instances used in the comparative. Consistently, it always helps to find the best value in all the trails. e dispersion of MHEDA is much less than other algorithms; it means that the solutions found by the MHEDA are more concentrated around the best value, than other algorithms, i.e., the average of solutions of the MHEDA converges better than other approaches to the best found value.
Based on the results, the MHEDA scheme tackles the inconvenience of the EDAs, i.e., lack of diversity of the solutions and poor ability of exploitation. e proposed  algorithm does not require evolutionary operators to get offspring, such as cross and mutation. e MHEDA makes use of the GMD process to establish a search direction. e EDA scheme considers population size, replacement (also known as generation gap), and selection strategy as key parameters. It is consistently with Grefenstette [67].
(i) e population size; in the current experiments, the population size ranged from 500 to 1000 solutions in increments of 500. (ii) e replacement; the current experiments allowed to vary the percentage of the population to be replaced during each generation between 50% and 100%, in increments of 50%. (iii) e selection; the experiments compared two ways to bubble, i.e., sorting by the makespan and sorting by the maximum workload.
A design of experiment is built to identify the best parameter of each parameter. Parameter tuning is detailed in Table 9.   Finally, the results of the parameter tuning are shown in Figure 18. ere is no statistically significant difference of any of the three controlled parameters (number of generations, initial population size, and selected population size). erefore, the parameters used are the same for all the algorithms.

Conclusions and Future Research
is paper considers the MHSP. It tries to determine the sequence of operations for each hoist to perform so that some performance metric is optimized. e MHEDA scheme is proposed for tackling the problem and simulating a solution.
e aforementioned instances were used as input and test parameters in order to validate the Mallows model as a probability model for the MHSP. e hybridization between the Mallows model and the EDA proposed helps to identify an explicit distribution over a set of permutations. Traditional operators are not considered for building suitable sequences, i.e., operation sequence vectors. ese are obtained from the Mallows model. e proposed exponential model, i.e., the Mallows model, is able to tackle the inconvenience any EDA has with permutation-based problems. e MHEDA scheme does not require to be reconfigured for solving the MHSP. e Mallows model is considered as a specific probability model for this issue. Based on the results, the exponential model is more effective against other algorithms. In addition, the MHEDA scheme considers three probabilistic models instead of only one as almost all the EDAs reported in the literature. e results of the MHEDA are more concentrated around the best Pareto-front obtained in all the trails. Meanwhile the rest of the algorithms have more disperse results. It is consistent in all the experiments detailed above.
Considering only the best solutions, in the selection process through a Pareto-front approach, is suitable to get a better estimation of the central permutation. A better estimate of the central permutation helps to improve the performance of the MHEDA scheme.
Based on the results detailed above, a simulation model is useful to model the critical variables that influence the performance of the process. e simulation model should be considered in almost any proposed solution. Moreover, simulation optimization is an enabling tool to handle diverse manufacturing limitations such as the MHSP. e simulation approach helps to model the difference between the classical job-shop configurations and the MHSP. e simulation language is able to order each couple of hoist operations, and avoid any collision between hoists, to satisfy the aforementioned constrains related to operations of the hoists.
In addition, the replacement step, in the MHEDA scheme, utilizes only the offspring to continue the evolutionary progress. It helps to reduce the dispersion of the results around the best solution.
Finally, as other previous works, this research contributes using the MHEDA as an optimization method for working with any simulation language.
Future research work could consider dynamic aspects such as failures of the hoists, shutdowns, jobs with priority, type of jobs, and type of hoists. erefore, the aforementioned dynamic issues should be integrated for any proposed algorithm.
Other comparisons should be investigated, as future research, in order to know in which other permutation-based problems, the MHEDA has a competitive performance.

Data Availability
All data are included in the manuscript.