# Multiobjective Simulation-Based Optimization Based on Artificial Immune Systems for a Distribution Center

Authors: Chris S. K. Leung; Henry Y. K. Lau
Journal: Journal of Optimization
 (2018)

---

## Abstract

Competitive market factors, such as more stringent government regulations, larger number of competitors, and shorter product life cycle, in recent years have created more significant pressure on the management in all supply chain parties. To this end, the ability of analyzing and evaluating systems and related operations involving the deployment of complex multiobjective material handling systems is vital for distribution practitioners. In this respect, simulation modeling techniques together with optimization have emerged as a very useful tool to facilitate the effective analysis of these complex operations and systems. In this paper, we apply a multiobjective simulation-based optimization framework consisting of a hybrid immune-inspired algorithm named Suppression-controlled Multiobjective Immune Algorithm (SCMIA) and a simulation model for solving a real-life multiobjective optimization problem. The results show that the framework is able to solve large scale problems with a large number of parameters, operators, and equipment involved.


---

## Body

## 1. Introduction

Simulation modeling is indeed a powerful industrial engineering technique for studying the functioning, performance, and operation of complex systems. As such, it becomes an extremely useful tool for stakeholders and decision-makers in various industries and domains including multiobjective optimization. By changing input data and operating parameters of a system being studied with simulation, predictions about the system’s behaviors can be obtained through computer-aided simulation for helping management make the right decisions. Unlike a mathematical model, simulation can handle a variety of complex factors that are commonly found in real world. In addition, simulation is a cost-effective means for existing process redesign and new system design because alternative solutions can be studied and evaluated for correctness and feasibility before actual implementation. More importantly, the accuracy of the performance measures of the complex systems obtained from simulation models is normally higher than that of analytical methods because analytical methods in general involve making unrealistic assumptions for the systems or problems under investigation [1].In real world, many problems no matter whether they are in the domain of engineering, finance, business, or science can be formulated into different forms of optimization problems. These problems are characterized by the requirement of finding the best possible solution that fulfills certain criteria under certain constraints. Most of the real-world optimization problems normally involve multiple objectives rather than one single objective, in which some objectives conflict with others. Solving this kind of problems is never an easy task because objectives of such problems are often found to be noncommensurable and conflicting. Very often, there is no single best solution to the multiobjectives optimization problems, but rather a set of optimal solutions that exists among the objectives. However, using simulation modeling alone cannot provide us with optimal solutions to these optimization problems. Therefore, an optimization algorithm is needed to guide the search process to the optimal solutions. Over the past decades, different algorithms have been developed for solving multiobjective optimization problems. However, some algorithms, such as Simulated Annealing (SA) [2] and Tabu Search (TS) [3–5], do not easily solve multiobjective optimization problems since they are initially developed for solving single objective optimization problems without having the ability to generate a set of optimal trade-off solutions. One of the possible approaches for them to solve the multiobjective optimization problems is to transform the multiple-objective problems into the single objective problems by emphasizing one particular optimal solution at each run [6]. Some other algorithms, such as Genetic Algorithm (GA) [7], Evolutionary Strategy (ES) [8], and Artificial Immune Systems (AIS) [9] which are inspired by natural processes, have been proved to be the effective algorithms in both single objective and multiobjective optimization contexts. Among these algorithms, AIS based on the concepts of biological immune system have received special attention among the research community because the biological immune system provides a rich source of stimulation and inspiration to the research community with their interesting characteristics: distributed nature, self-organization, memory, and learning capabilities. Motivated by its great potential for solving multiple-objective optimization problems, the study reported in this paper is to demonstrate how a real-life multiobjective simulation-based optimization problem in logistics industry where material handling system (MHS) is involved can be solved by a simulation-based optimization framework comprising a simulation tool together with an AIS-based algorithm.This paper proceeds as follows: Section2 starts by giving a brief overview of material handling and multiobjective optimization. Section 3 introduces the multiobjective optimization framework for simulation-based optimization. Section 4 evaluates the performance of the framework through a real-life case study on a MHS by benchmarking with some well-known optimization algorithms. Finally, concluding remarks are presented in Section 5.


## 2. Background

### 2.1. Importance of Material Handling in Distribution Industry

Material handling, which is the total management of material concerns in an operation, is a vital element of industrial processes. Material handling involves a variety of operations including the movement, storage, protection, and control of materials, products, and wastes throughout the processes of manufacturing, distribution, and disposal. Having efficient MHS is of great importance to various industries to maintain and facilitate a continuous flow of material through the workplace and guarantees that required materials are available when needed. It is especially important to logistics and manufacturing industries as it accounts for a large percentage of the operation in these industries. In the manufacturing sector, the time spent on different kinds of product handling and transportation can be as much as the time used on the value-added processes. Banks et al. [10] claimed that the time of material handling accounts for about 85% of the total manufacturing time. In addition to time, the money used on material handling activities is equally high because the material handling equipment and systems require large capital investment in terms of system design, installation, operations, maintenance, and so forth. A number of studies conducted in different industries show that the cost of material handling alone is about 20–25% of the total manufacturing cost [11].There are different kinds of material handling equipment and systems available that range from simple hand truck and pallet rack to complex conveyor system and Automated Storage and Retrieval System (AS/RS). A typical MHS is composed of different smaller components closely working together, thus making business activities more efficient and cost-effective. Over the past decade, MHS and its function have undergone a big change because of new advances in technologies, such as the development and applications of automation techniques and robotics, by which a large number of manual handling jobs are replaced by machines. Since the entire production or distribution process is automated, MHS has to respond just-in-time to the requirements of different processing activities. These new technologies, today, increasingly become prevalent in different industries as they help ease drudgery for manual labors and some of the mechanized or automated handling jobs are physically impossible to be done by workers.Norman [12] claimed that equipment capacity, speed, and arrangement are the most critical considerations when modeling and optimizing MHS. Capacity under this context is the maximum quantity of products handled by the equipment. Speed is the average operating speed of the equipment, which may include acceleration, deceleration, and speeds of various equipment components. Configuration is the layout and structure of the MHS or its moving paths.


### 2.2. Optimization of MHS via Simulation

In the literature, there are a number of studies that dedicatedly contribute to the optimization of MHS via simulation. For example, Ebbesen et al. [13] studied the baggage handling system at airports. They developed an approach to optimize the design of conveyor systems, that is, the design of tracks suited for baggage handling systems with the use of a time domain simulation model of the entire system. Sergueyevich et al. [14] proposed a simulation model of the overhead monorail conveyor system coupled with statistical methods for analyzing and solving the multiobjective optimization problem regarding the manufacturing process. The aim was to determine the optimum speed of conveyor, lengths of queues, time in system, capacity of machines, and so forth under certain limitations. Elahi et al. [15] studied the General Motors paint shop conveyor system by developing a simulation model. The model works firmly with a decision optimizer incorporating integer linear programming model and dynamic programming model at critical points such as the beginning and end of buffer conveyors in the system in order to regroup batches of different color cars. Leung and Lau [16] proposed a simulation-based optimization framework that combines the processes of optimization and simulation for solving typical linear optimization problems related to logistics and production operation. The framework integrates an AIS-based algorithm with a simulation tool for the evaluation of optimal system parameters and to reveal the performance of systems. Subulan and Cakmakci [17] made use of ARENA simulation program and Taguchi experimental design method to build a solution model for effectively designing material handling–transfer systems and optimizing the performance of automation technologies in automobile industry. Chang et al. [18] proposed a framework that integrates simulation optimization and data envelopment analysis techniques to find out the optimal vehicle fleet size for a multiobjective problem in automated materials handling systems. Lin and Huang [19] extended the optimal computing budget allocation by adding Genetic Algorithm together with the help of a simulation model for optimizing the vehicle allocation for the automated material handling system in semiconductor industry.


### 2.3. Multiobjective Simulation-Based Optimization Approaches

In practice, optimization problems involve several objectives that often conflict with each other and must be simultaneously optimized so that a possibly uncountable set of trade-off solutions rather than a single optimal point is found with respect to the contradicting objectives. Therefore, the aim of these problems is to find out the global trade-off solutions that effectively spread over the Pareto front. No solution from the Pareto front is worse than any other solution because it is better in at least one objective. These problems are normally termed as multiobjective problems, which were first studied in an economic context and then extended to the fields of science and engineering [20]. Since multiobjective optimization problems involve several objectives, the view towards optimum has changed, hence changing the aim from finding a single solution to obtaining a set of compromised solutions. Today, the notion of optimum for a multiobjective optimization problem is frequently called Pareto optimum, which was first proposed by Edgeworth [21] and then generalized by Pareto [22, 23].As is known, the notion of optimality for multiobjective optimization problems is different from that of single objective optimization problems because the aim of multiobjective optimization problems is to find a set of optimal trade-off solutions rather than a single optimal solution. Thus, in the absence of preference information on the objectives, the concept of Pareto optimality is adopted in this study for solving multiobjective optimization problems [24] and several definitions about Pareto optimality [20, 25, 26] are considered and stated below.Definition 1 (Pareto optimality).
A pointX∗∈Ω in the search space is said to be Pareto optimal/nondominated with respect to Ω if and only if there are no other solutions X∈Ω for which f→(X) dominates f→(X∗).In other words, X∗ is Pareto optimal with regard to the whole decision variable space if it cannot be improved in any one objective without resulting in a simultaneous degradation in other objectives.Definition 2 (Pareto dominance).
Consider, without loss of generality, two decision vectorsXa and Xb∈X of a minimization problem. Then, a vector of decision variables Xa is said to dominate another vector Xb, which is denoted by Xa≺Xb, if and only if Xa is partially less than Xb. That is, the following two conditions must be satisfied:(1)
Xa is not worse than Xb in all objectives(1)fiXa≤fiXb,∀i=1,2,…,r(2)
Xa is better than Xb in at least one objective(2)fiXa<fiXb,∃i=1,2,…,r,where i is the number of objectives.However, when any of these conditions are violated, the two solutionsXa and Xb are said to be indifferent to each other instead of dominating the other or being dominated by the other. Based on the above relations, Pareto optimal set, Pareto front, and nondominated set are defined below.Definition 3 (Pareto optimal set).
Pareto optimal set of solutions is a collection of all Pareto optimal solutions, which is defined as(3)PS≔X∈Ω∣¬∃X∗∈Ω,X∗≺X.
Pareto optimal solutions are those solutions in the decision variable space whose corresponding objective vector elements cannot be all simultaneously improved [27]. The solutions inside the Pareto optimal set may have no apparent relationship except their membership in the set. Pareto optimal solutions are classified as such based on their values being evaluated through whatever means.Definition 4 (Pareto front).
A surface or line containing all nondominated solutions is called Pareto front, which is represented by(4)PF≔f→X∣X∈PS.
According to the literature, finding an analytical expression of the Pareto front is a very difficult task. Therefore, a common approach for Pareto front generation is to find out the points within Ω and their corresponding valuef(X),  X∈Ω. When this procedure is repeated a sufficient number of times, the nondominated points and hence the Pareto front are most likely to be found in the objective space [20].Definition 5 (nondominated set).
Of a solution setS, a nondominated set of solutions S∗ comprises solutions that are not dominated by other solutions in the set S. It is worth mentioning that while Pareto optimal solutions are always nondominated solutions, nondominated solutions may include both non-Pareto optimal solutions and Pareto optimal solutions, thus revealing that the true Pareto optimal solutions could hardly be represented by the nondominated solutions obtained from running an optimizer. Thus, the idea, stated by van Veldhuizen and Lamont [28], about the true Pareto front PFtrue distinguishing it from the final set of nondominated solutions found by an algorithm is called known Pareto front PFknown.Modern optimization approaches are very often population-based and evolutionary in nature. In such methods, the search for the global optima essentially comprises an iterative process that replaces the candidate solutions in the population by newly generated ones with an aim of achieving continuous improvement in the performance of the best candidate solutions through the help of mechanisms that guide the search to find a set of nondominated solutions. The use of the modern optimization approaches, especially population-based evolutionary algorithms, to solving the complex multiobjective optimization problems has been motivated mainly because of the following critical reasons. First, population-based evolutionary algorithms can recognize the specificity of multiobjective optimization problems by working simultaneously on all objectives and finally generating a group of optimal trade-off solutions, thus forming a uniformly distributed Pareto front. Second, as the name implies, the population-based approaches can deal with a population of candidate solutions simultaneously, allowing the generation of several elements of the Pareto optimal set in a single run of an optimizer instead of performing many separate runs when using classical mathematical programming methods [29]. This allows the decision-makers to simply pick the one that best fits the problem at hand, thus preventing the need to reconfigure and to rerun the optimization tool for finding other alternative Pareto optimal solutions. Third, their capability of maintaining diversity among the candidate solutions in the population is important to prevent the search from premature convergence to a specific region of the solution space, thus allowing a better exploration of the solution space and minimizing the susceptibility of the search to the presence of poor local optima in the optimization problems [30]. The final reason is that the population-based approaches are less susceptible to the continuity or shape of the Pareto front as they can deal with concave and discontinuous Pareto fronts without difficulty [27, 31]. Nevertheless, on the other hand, the weaknesses of these population-based methods are that they do not guarantee the identification of optimal trade-off solutions and the solutions obtained are likely to be stuck at some “good” approximations [32].During the past few decades, a large number of publications have been done in population-based evolutionary algorithms and proved to be effective for solving multiobjective optimization problems since the first multiobjective evolutionary algorithm has been developed by Schaffer [33]. These algorithms include PAES [34], NSGA-II [35], PESA [36], SPEA2 [37], PESA-II [38], micro-GA [39], micro-GA2 [40], omni-aiNet [41], NNIA [42], omni-AIOS [43], MTLBO [44], and SCMIA [45].


## 2.1. Importance of Material Handling in Distribution Industry

Material handling, which is the total management of material concerns in an operation, is a vital element of industrial processes. Material handling involves a variety of operations including the movement, storage, protection, and control of materials, products, and wastes throughout the processes of manufacturing, distribution, and disposal. Having efficient MHS is of great importance to various industries to maintain and facilitate a continuous flow of material through the workplace and guarantees that required materials are available when needed. It is especially important to logistics and manufacturing industries as it accounts for a large percentage of the operation in these industries. In the manufacturing sector, the time spent on different kinds of product handling and transportation can be as much as the time used on the value-added processes. Banks et al. [10] claimed that the time of material handling accounts for about 85% of the total manufacturing time. In addition to time, the money used on material handling activities is equally high because the material handling equipment and systems require large capital investment in terms of system design, installation, operations, maintenance, and so forth. A number of studies conducted in different industries show that the cost of material handling alone is about 20–25% of the total manufacturing cost [11].There are different kinds of material handling equipment and systems available that range from simple hand truck and pallet rack to complex conveyor system and Automated Storage and Retrieval System (AS/RS). A typical MHS is composed of different smaller components closely working together, thus making business activities more efficient and cost-effective. Over the past decade, MHS and its function have undergone a big change because of new advances in technologies, such as the development and applications of automation techniques and robotics, by which a large number of manual handling jobs are replaced by machines. Since the entire production or distribution process is automated, MHS has to respond just-in-time to the requirements of different processing activities. These new technologies, today, increasingly become prevalent in different industries as they help ease drudgery for manual labors and some of the mechanized or automated handling jobs are physically impossible to be done by workers.Norman [12] claimed that equipment capacity, speed, and arrangement are the most critical considerations when modeling and optimizing MHS. Capacity under this context is the maximum quantity of products handled by the equipment. Speed is the average operating speed of the equipment, which may include acceleration, deceleration, and speeds of various equipment components. Configuration is the layout and structure of the MHS or its moving paths.


## 2.2. Optimization of MHS via Simulation

In the literature, there are a number of studies that dedicatedly contribute to the optimization of MHS via simulation. For example, Ebbesen et al. [13] studied the baggage handling system at airports. They developed an approach to optimize the design of conveyor systems, that is, the design of tracks suited for baggage handling systems with the use of a time domain simulation model of the entire system. Sergueyevich et al. [14] proposed a simulation model of the overhead monorail conveyor system coupled with statistical methods for analyzing and solving the multiobjective optimization problem regarding the manufacturing process. The aim was to determine the optimum speed of conveyor, lengths of queues, time in system, capacity of machines, and so forth under certain limitations. Elahi et al. [15] studied the General Motors paint shop conveyor system by developing a simulation model. The model works firmly with a decision optimizer incorporating integer linear programming model and dynamic programming model at critical points such as the beginning and end of buffer conveyors in the system in order to regroup batches of different color cars. Leung and Lau [16] proposed a simulation-based optimization framework that combines the processes of optimization and simulation for solving typical linear optimization problems related to logistics and production operation. The framework integrates an AIS-based algorithm with a simulation tool for the evaluation of optimal system parameters and to reveal the performance of systems. Subulan and Cakmakci [17] made use of ARENA simulation program and Taguchi experimental design method to build a solution model for effectively designing material handling–transfer systems and optimizing the performance of automation technologies in automobile industry. Chang et al. [18] proposed a framework that integrates simulation optimization and data envelopment analysis techniques to find out the optimal vehicle fleet size for a multiobjective problem in automated materials handling systems. Lin and Huang [19] extended the optimal computing budget allocation by adding Genetic Algorithm together with the help of a simulation model for optimizing the vehicle allocation for the automated material handling system in semiconductor industry.


## 2.3. Multiobjective Simulation-Based Optimization Approaches

In practice, optimization problems involve several objectives that often conflict with each other and must be simultaneously optimized so that a possibly uncountable set of trade-off solutions rather than a single optimal point is found with respect to the contradicting objectives. Therefore, the aim of these problems is to find out the global trade-off solutions that effectively spread over the Pareto front. No solution from the Pareto front is worse than any other solution because it is better in at least one objective. These problems are normally termed as multiobjective problems, which were first studied in an economic context and then extended to the fields of science and engineering [20]. Since multiobjective optimization problems involve several objectives, the view towards optimum has changed, hence changing the aim from finding a single solution to obtaining a set of compromised solutions. Today, the notion of optimum for a multiobjective optimization problem is frequently called Pareto optimum, which was first proposed by Edgeworth [21] and then generalized by Pareto [22, 23].As is known, the notion of optimality for multiobjective optimization problems is different from that of single objective optimization problems because the aim of multiobjective optimization problems is to find a set of optimal trade-off solutions rather than a single optimal solution. Thus, in the absence of preference information on the objectives, the concept of Pareto optimality is adopted in this study for solving multiobjective optimization problems [24] and several definitions about Pareto optimality [20, 25, 26] are considered and stated below.Definition 1 (Pareto optimality).
A pointX∗∈Ω in the search space is said to be Pareto optimal/nondominated with respect to Ω if and only if there are no other solutions X∈Ω for which f→(X) dominates f→(X∗).In other words, X∗ is Pareto optimal with regard to the whole decision variable space if it cannot be improved in any one objective without resulting in a simultaneous degradation in other objectives.Definition 2 (Pareto dominance).
Consider, without loss of generality, two decision vectorsXa and Xb∈X of a minimization problem. Then, a vector of decision variables Xa is said to dominate another vector Xb, which is denoted by Xa≺Xb, if and only if Xa is partially less than Xb. That is, the following two conditions must be satisfied:(1)
Xa is not worse than Xb in all objectives(1)fiXa≤fiXb,∀i=1,2,…,r(2)
Xa is better than Xb in at least one objective(2)fiXa<fiXb,∃i=1,2,…,r,where i is the number of objectives.However, when any of these conditions are violated, the two solutionsXa and Xb are said to be indifferent to each other instead of dominating the other or being dominated by the other. Based on the above relations, Pareto optimal set, Pareto front, and nondominated set are defined below.Definition 3 (Pareto optimal set).
Pareto optimal set of solutions is a collection of all Pareto optimal solutions, which is defined as(3)PS≔X∈Ω∣¬∃X∗∈Ω,X∗≺X.
Pareto optimal solutions are those solutions in the decision variable space whose corresponding objective vector elements cannot be all simultaneously improved [27]. The solutions inside the Pareto optimal set may have no apparent relationship except their membership in the set. Pareto optimal solutions are classified as such based on their values being evaluated through whatever means.Definition 4 (Pareto front).
A surface or line containing all nondominated solutions is called Pareto front, which is represented by(4)PF≔f→X∣X∈PS.
According to the literature, finding an analytical expression of the Pareto front is a very difficult task. Therefore, a common approach for Pareto front generation is to find out the points within Ω and their corresponding valuef(X),  X∈Ω. When this procedure is repeated a sufficient number of times, the nondominated points and hence the Pareto front are most likely to be found in the objective space [20].Definition 5 (nondominated set).
Of a solution setS, a nondominated set of solutions S∗ comprises solutions that are not dominated by other solutions in the set S. It is worth mentioning that while Pareto optimal solutions are always nondominated solutions, nondominated solutions may include both non-Pareto optimal solutions and Pareto optimal solutions, thus revealing that the true Pareto optimal solutions could hardly be represented by the nondominated solutions obtained from running an optimizer. Thus, the idea, stated by van Veldhuizen and Lamont [28], about the true Pareto front PFtrue distinguishing it from the final set of nondominated solutions found by an algorithm is called known Pareto front PFknown.Modern optimization approaches are very often population-based and evolutionary in nature. In such methods, the search for the global optima essentially comprises an iterative process that replaces the candidate solutions in the population by newly generated ones with an aim of achieving continuous improvement in the performance of the best candidate solutions through the help of mechanisms that guide the search to find a set of nondominated solutions. The use of the modern optimization approaches, especially population-based evolutionary algorithms, to solving the complex multiobjective optimization problems has been motivated mainly because of the following critical reasons. First, population-based evolutionary algorithms can recognize the specificity of multiobjective optimization problems by working simultaneously on all objectives and finally generating a group of optimal trade-off solutions, thus forming a uniformly distributed Pareto front. Second, as the name implies, the population-based approaches can deal with a population of candidate solutions simultaneously, allowing the generation of several elements of the Pareto optimal set in a single run of an optimizer instead of performing many separate runs when using classical mathematical programming methods [29]. This allows the decision-makers to simply pick the one that best fits the problem at hand, thus preventing the need to reconfigure and to rerun the optimization tool for finding other alternative Pareto optimal solutions. Third, their capability of maintaining diversity among the candidate solutions in the population is important to prevent the search from premature convergence to a specific region of the solution space, thus allowing a better exploration of the solution space and minimizing the susceptibility of the search to the presence of poor local optima in the optimization problems [30]. The final reason is that the population-based approaches are less susceptible to the continuity or shape of the Pareto front as they can deal with concave and discontinuous Pareto fronts without difficulty [27, 31]. Nevertheless, on the other hand, the weaknesses of these population-based methods are that they do not guarantee the identification of optimal trade-off solutions and the solutions obtained are likely to be stuck at some “good” approximations [32].During the past few decades, a large number of publications have been done in population-based evolutionary algorithms and proved to be effective for solving multiobjective optimization problems since the first multiobjective evolutionary algorithm has been developed by Schaffer [33]. These algorithms include PAES [34], NSGA-II [35], PESA [36], SPEA2 [37], PESA-II [38], micro-GA [39], micro-GA2 [40], omni-aiNet [41], NNIA [42], omni-AIOS [43], MTLBO [44], and SCMIA [45].


## 3. A Multiobjective Simulation-Based Optimization Framework

### 3.1. Mechanisms of the Multiobjective Simulation-Based Optimization Framework

The optimization framework adopted in this study is developed by taking advantage of the idea of separation between the optimization method and the simulation model. For this reason, the optimization framework can remain the same or require only minor modifications such as changing range of parameters, data type of decision variables, and number of decision variables, to optimize the simulation model that incorporates new requirements. This framework in fact is a modified version of Leung and Lau’s work [16], which incorporates a multiobjective optimization algorithm instead of a single objective algorithm (Figure 1). The multiobjective algorithm adopted in the framework is Suppression-controlled Multiobjective Immune Algorithm (SCMIA) proposed by Leung and Lau [45]. Concepts of the algorithm are briefly discussed in the next section.Figure 1
The framework for multiobjective simulation-based optimization (source: modified from Leung and Lau’s work [16]).The framework consists of two critical components, namely,(1) the simulation model and (2) the immunity-inspired multiobjective optimization algorithm. A problem or a system to be optimized is represented by a discrete-event simulation model, which is developed separately from the optimization algorithm. The simulation model under study takes inputs (solutions) produced from the optimization algorithm (except the initial inputs which are entered by user). Based on these inputs, the simulation is executed automatically with various parameter settings. In other words, the simulation of the system with the parameters under consideration is iteratively conducted in an automatic manner. And then the simulation model generates several output performance metrics as feedback on the quality of solutions obtained from the optimization algorithm. During the process of simulation-based optimization, it is worth emphasizing that the general guidelines of discrete-event simulation should be followed, such as performing multiple replications of the simulation with the same input parameters but different random number streams to minimize the variability of the simulation outputs.In essence, the AIS-based multiobjective optimization algorithm is a search algorithm. The simulation results (output performance metrics) direct the algorithm to search for a new set of input parameters that takes the system towards its optimal setting with respect to certain criteria. Hence a feedback mechanism is incorporated to direct the search for optimal solution in a controlled manner. In other words, the generation of a new set of input parameters for the simulation model is based on the simulation results of past evaluations. As such, the input parameters and the output performance metrics shown in the figure in fact serve as the feedback to the simulation model and the optimization algorithm, respectively. This iterative process repeats until prespecified termination criteria are met. The criteria could include the following: for example, no improving solutions can be found or the predefined number of iterations is reached.In this framework, the optimization algorithm was implemented with Excel VBA, whereas the simulation model was developed by using the FlexSim simulation tool [46].


### 3.2. Multiobjective Simulation-Based Optimization Algorithm

The fundamental of the hybrid AIS-based optimization algorithm [45] adopted in the framework was inspired from mechanisms of biological immune system and biological evolution. It was developed by hybridizing the clonal selection principle and immune network theory with the idea from Genetic Algorithm (GA). The algorithm makes use of the Pareto dominance for fitness assignment and some common AIS-based algorithm’s features for guiding the search process, such as clonal selection and expansion, affinity maturation, antibody concentration, metadynamics, and immune memory. The interesting feature of this algorithm is the introduction of an innovative suppression operator, which is used to help eliminate similar antibodies, hence significantly minimizing the number of unnecessary searches and increasing the population diversity. The similarity among antibodies is determined in terms of both the objective space and the decision variable space to ensure that only similar antibodies are eliminated in the suppression operation. Moreover, a modified crossover operator originated from the biological evolution was also developed to help further enhance the diversity of the clone population and the convergence of the algorithm because some good genes from the elite parents can be passed to the offspring for facilitating the search of optimal solutions, otherwise it may take a longer time to converge towards the Pareto front [27] especially in simulation-based optimization context. The mapping between the biological immune system and the proposed artificial one is given in Table 1.Table 1
Mapping between the biological immune system and SCMIA.
Biological Immune System SCMIA Antigen (Ag) Objective function (simulation model in simulation study) to be optimized Antibody (Ab) Candidate solution (a set of decision variables) to be optimized Ag-Ab affinity Fitness value of each candidate solution evaluated based on Pareto dominance Ab-Ab affinity Crowding-distance proposed by Deb et al. [6] working as a measure of population diversity Immune suppression Mechanism to control the number of nearby candidate solutions based on similarity among candidate solutions in both the objective space and decision variable space Memory cell Current best nondominated solutionThe algorithm comprises five immune operators: cloning operator, hypermutation operator, suppression operator, selection and receptor editing operator, and memory updating operator, and one genetic operator: crossover operator. Each of them takes responsibility for different tasks for the purpose of finding uniformly distributed Pareto front. The block diagram showing the computational steps for SCMIA is presented in Figure2 and the details of these operators are provided as below.Figure 2
Computational steps for SCMIA.Nondominated Neighbor-Based Selection.Based on the idea proposed by Gong et al. [42], only nondominated antibodies are selected to form an active parent population A(t) (where t is the iteration counter) with the size of NA for undergoing cloning, crossover, and hypermutation. This process mimics the biological immune system in which only those B-cells (antibodies) that are capable of binding with foreign antigens will undergo clonal expansion and then hypermutation. However, if the number of nondominated antibodies is larger than NA, an antibody density measure called crowding-distance [6] analogous with Ab-Ab affinity in biological immune system is employed for selecting the nondominated antibodies. The crowding-distance or Ab-Ab affinity for a nondominated antibody is computed as follows:(5)dabi∗=∑j=1rabi+1∗fj-abi-1∗fjabfjmax-abfjmin,where dabi∗ is the crowding-distance of the ith nondominated antibody ab∗, r is the number of objectives, ab(fjmax) and ab(fjmin) are the maximum and minimum fitness values of the jth objective, and abi+1∗fj and abi-1∗fj are the fitness values of the nearest neighboring antibodies from both sides.Cloning Operator. It enlarges the population by generating a number of copies of each antibody in the active parent populationA(t) and the number of copies is directly proportional to its Ab-Ab affinity, thus forming a clone population, which is denoted as C(t). Hence the size of the population now is NA+Nc and Nc is obtained by(6)ci=roundcmax×dabi∗,Nc=∑ci,where Nc is the total number of copies produced, ci is the clone size for the antibody abi∗∈A(t), cmax is the predefined maximum clone size of each antibody, and round() is an operator for rounding its argument to the closest integer. Thus, the higher the Ab-Ab affinity an antibody has, the more the number of copies it generates.Hypermutation Operator. It brings multipoint mutations to the pool of clones, hoping for producing better offspring. The mutations are dependent on the Ab-Ab affinity of their parents for the purpose of preventing the crowding of antibodies and maintaining the population diversity. The mechanism is that if the Ab-Ab affinity of an active parent is lower, the mutation rate for the clones of the parent should be higher to obtain mutated children in sparse locations in order to diversify the population. The clones are mutated proportionally as follows:(7)α=e-p×dab∗Cmt=Ct+α×R,where α represents the mutation rate inversely proportional to the Ab-Ab affinity d(ab∗), p is an exponential coefficient controlling the decay of α,R∈[-1,1] is a m-dimensional random vector obtained with uniform distribution, and Cm(t) is the mutated clone population.Crossover Operator. It works on the mutated clone population to generate a mature clone population denoted as Cc(t) with the size of Nc. It is a modified single point crossover operator through which offspring is generated by randomly selecting a single crossover point on a clone and then swapping its content beyond that point with that of a parent antibody randomly selected from A(t). This operator is introduced to control the diversity of the clone population and the convergence of the algorithm. That being said, the diversity can be enhanced while the quick convergence can be ensured because some good genes from the active parent can be passed to the offspring.Suppression Operator. It works on the whole population including the parent cells and mutated clones A(t)∪Cc(t) to eliminate similar individuals so as to avoid a particular search space being over exploited and acquire the uniformly distributed Pareto front based on the idea of immune network theory [47] such that B-cells are stimulated and suppressed by not only non-self-antigens but the interacted B-cells. Different from other AIS-based multiobjective optimization algorithms, the similarity among antibodies in this algorithm is determined in terms of both the objective space and the decision variable space so as to determine whether to retain or discard individual antibodies. Therefore, the suppression operation in this algorithm has two phases: in first phase, the suppression will be applied to all antibodies and the similarity between two antibodies is defined as follows:(8)dOaba,abbj=abafj-abbfj≤δj,where dOaba,abbj is the distance between antibodies aba and abb in terms of jth objective and δj refers to the threshold value for jth objective.In this phase, if the distances for all objectives between two antibodies are smaller than the thresholds, the two antibodies are said to be similar and hence the cell with poorer Pareto fitnesspf will be suppressed and eliminated from the population.In second phase, the suppression will only be applied to the similarity between nondominated cells and dominated cells and the similarity between two antibodies is defined as follows:(9)dVaba,abb=∑j=1mabaxj-abbxj2≤ε,where dV(aba,abb) is the Euclidean distance in decision variable space between the two antibodies aba (dominated cell) and abb (nondominated cell) ∈A(t)∪Cc(t), m is the length of each antibody, and ε refers to the threshold value for the decision variable space.In this phase, if the distance between two cells is smaller than the threshold in decision variable space, the two cells are said to be similar and hence the dominated cell will be suppressed and eliminated from the population. Eventually, surviving populationsAs(t)∪Cs(t) are obtained and then enter into the selection and receptor editing process and memory updating process simultaneously.To enhance the population diversity and facilitate the search of uniformly distributed nondominated solutions, the threshold values for the decision variable space and the objective space are adaptive values that are dynamically calculated based on the maximum and minimum values found at each iteration. The adaptive threshold values are computed according to the following equations:(10)ε=abxmax-abxminn,δj=abfjmax-abfjminn,where ab(xmax) and ab(xmin) are the maximum and minimum distance values in the decision variable space, ab(fjmax) and ab(fjmin) are the maximum and minimum fitness values of the jth objective, and n is the number of antibodies in the population under investigation.Selection and Receptor Editing Operator. It works like a director to guide the search towards the promising regions by selecting all nondominated antibodies with respect to the Pareto fitnesspf from the populations As(t)∪Cs(t) to form a new population Ab(t+1) with the size of N for the next generation t+1. If the new population Ab(t+1) is not full, dominated antibodies with a better Pareto fitnesspf are selected and some genes of these antibodies are then randomly selected to be replaced by randomly generated genes. These restructured antibodies are finally added to the new population until the population is full in order to further enhance the population diversity. This process actually mimics the process of receptor editing in the biological immune system. However, if the number of nondominated antibodies found exceeds the limit of the population size N, only N nondominated antibodies with higher Ab-Ab affinity d(ab) are selected.Memory Updating Operator.It works as an elitist mechanism for helping preserve the best solutions that represent the Pareto front found over the search process. The memory set is updated at each iteration through the following procedure:(1)
If the size ofP(t) is equal to Nm, there is no memory archiving procedure to proceed.(2)
If the size ofP(t) is smaller than Nm, an archiving procedure is triggered to copy some of the dominated antibodies that have the best Pareto fitnesspf and Ab-Ab affinity d(ab) to P(t) for filling the vacant space up.(3)
If the size ofP(t) is larger than Nm, an archiving procedure is triggered to compress P(t) based on the antibody similarity (the measure used in the suppression operator) until its size becomes Nm.After performing the final memory updating procedure at the last iteration, the memory setP(t) obtained represents the resulting solution set of the algorithm.The algorithm is conducted by applying these heuristic and stochastic operators on the antibody population for balancing both the local and global search capabilities. For the detailed discussion of the hybrid multiobjective algorithm, one can refer to [45].


## 3.1. Mechanisms of the Multiobjective Simulation-Based Optimization Framework

The optimization framework adopted in this study is developed by taking advantage of the idea of separation between the optimization method and the simulation model. For this reason, the optimization framework can remain the same or require only minor modifications such as changing range of parameters, data type of decision variables, and number of decision variables, to optimize the simulation model that incorporates new requirements. This framework in fact is a modified version of Leung and Lau’s work [16], which incorporates a multiobjective optimization algorithm instead of a single objective algorithm (Figure 1). The multiobjective algorithm adopted in the framework is Suppression-controlled Multiobjective Immune Algorithm (SCMIA) proposed by Leung and Lau [45]. Concepts of the algorithm are briefly discussed in the next section.Figure 1
The framework for multiobjective simulation-based optimization (source: modified from Leung and Lau’s work [16]).The framework consists of two critical components, namely,(1) the simulation model and (2) the immunity-inspired multiobjective optimization algorithm. A problem or a system to be optimized is represented by a discrete-event simulation model, which is developed separately from the optimization algorithm. The simulation model under study takes inputs (solutions) produced from the optimization algorithm (except the initial inputs which are entered by user). Based on these inputs, the simulation is executed automatically with various parameter settings. In other words, the simulation of the system with the parameters under consideration is iteratively conducted in an automatic manner. And then the simulation model generates several output performance metrics as feedback on the quality of solutions obtained from the optimization algorithm. During the process of simulation-based optimization, it is worth emphasizing that the general guidelines of discrete-event simulation should be followed, such as performing multiple replications of the simulation with the same input parameters but different random number streams to minimize the variability of the simulation outputs.In essence, the AIS-based multiobjective optimization algorithm is a search algorithm. The simulation results (output performance metrics) direct the algorithm to search for a new set of input parameters that takes the system towards its optimal setting with respect to certain criteria. Hence a feedback mechanism is incorporated to direct the search for optimal solution in a controlled manner. In other words, the generation of a new set of input parameters for the simulation model is based on the simulation results of past evaluations. As such, the input parameters and the output performance metrics shown in the figure in fact serve as the feedback to the simulation model and the optimization algorithm, respectively. This iterative process repeats until prespecified termination criteria are met. The criteria could include the following: for example, no improving solutions can be found or the predefined number of iterations is reached.In this framework, the optimization algorithm was implemented with Excel VBA, whereas the simulation model was developed by using the FlexSim simulation tool [46].


## 3.2. Multiobjective Simulation-Based Optimization Algorithm

The fundamental of the hybrid AIS-based optimization algorithm [45] adopted in the framework was inspired from mechanisms of biological immune system and biological evolution. It was developed by hybridizing the clonal selection principle and immune network theory with the idea from Genetic Algorithm (GA). The algorithm makes use of the Pareto dominance for fitness assignment and some common AIS-based algorithm’s features for guiding the search process, such as clonal selection and expansion, affinity maturation, antibody concentration, metadynamics, and immune memory. The interesting feature of this algorithm is the introduction of an innovative suppression operator, which is used to help eliminate similar antibodies, hence significantly minimizing the number of unnecessary searches and increasing the population diversity. The similarity among antibodies is determined in terms of both the objective space and the decision variable space to ensure that only similar antibodies are eliminated in the suppression operation. Moreover, a modified crossover operator originated from the biological evolution was also developed to help further enhance the diversity of the clone population and the convergence of the algorithm because some good genes from the elite parents can be passed to the offspring for facilitating the search of optimal solutions, otherwise it may take a longer time to converge towards the Pareto front [27] especially in simulation-based optimization context. The mapping between the biological immune system and the proposed artificial one is given in Table 1.Table 1
Mapping between the biological immune system and SCMIA.
Biological Immune System SCMIA Antigen (Ag) Objective function (simulation model in simulation study) to be optimized Antibody (Ab) Candidate solution (a set of decision variables) to be optimized Ag-Ab affinity Fitness value of each candidate solution evaluated based on Pareto dominance Ab-Ab affinity Crowding-distance proposed by Deb et al. [6] working as a measure of population diversity Immune suppression Mechanism to control the number of nearby candidate solutions based on similarity among candidate solutions in both the objective space and decision variable space Memory cell Current best nondominated solutionThe algorithm comprises five immune operators: cloning operator, hypermutation operator, suppression operator, selection and receptor editing operator, and memory updating operator, and one genetic operator: crossover operator. Each of them takes responsibility for different tasks for the purpose of finding uniformly distributed Pareto front. The block diagram showing the computational steps for SCMIA is presented in Figure2 and the details of these operators are provided as below.Figure 2
Computational steps for SCMIA.Nondominated Neighbor-Based Selection.Based on the idea proposed by Gong et al. [42], only nondominated antibodies are selected to form an active parent population A(t) (where t is the iteration counter) with the size of NA for undergoing cloning, crossover, and hypermutation. This process mimics the biological immune system in which only those B-cells (antibodies) that are capable of binding with foreign antigens will undergo clonal expansion and then hypermutation. However, if the number of nondominated antibodies is larger than NA, an antibody density measure called crowding-distance [6] analogous with Ab-Ab affinity in biological immune system is employed for selecting the nondominated antibodies. The crowding-distance or Ab-Ab affinity for a nondominated antibody is computed as follows:(5)dabi∗=∑j=1rabi+1∗fj-abi-1∗fjabfjmax-abfjmin,where dabi∗ is the crowding-distance of the ith nondominated antibody ab∗, r is the number of objectives, ab(fjmax) and ab(fjmin) are the maximum and minimum fitness values of the jth objective, and abi+1∗fj and abi-1∗fj are the fitness values of the nearest neighboring antibodies from both sides.Cloning Operator. It enlarges the population by generating a number of copies of each antibody in the active parent populationA(t) and the number of copies is directly proportional to its Ab-Ab affinity, thus forming a clone population, which is denoted as C(t). Hence the size of the population now is NA+Nc and Nc is obtained by(6)ci=roundcmax×dabi∗,Nc=∑ci,where Nc is the total number of copies produced, ci is the clone size for the antibody abi∗∈A(t), cmax is the predefined maximum clone size of each antibody, and round() is an operator for rounding its argument to the closest integer. Thus, the higher the Ab-Ab affinity an antibody has, the more the number of copies it generates.Hypermutation Operator. It brings multipoint mutations to the pool of clones, hoping for producing better offspring. The mutations are dependent on the Ab-Ab affinity of their parents for the purpose of preventing the crowding of antibodies and maintaining the population diversity. The mechanism is that if the Ab-Ab affinity of an active parent is lower, the mutation rate for the clones of the parent should be higher to obtain mutated children in sparse locations in order to diversify the population. The clones are mutated proportionally as follows:(7)α=e-p×dab∗Cmt=Ct+α×R,where α represents the mutation rate inversely proportional to the Ab-Ab affinity d(ab∗), p is an exponential coefficient controlling the decay of α,R∈[-1,1] is a m-dimensional random vector obtained with uniform distribution, and Cm(t) is the mutated clone population.Crossover Operator. It works on the mutated clone population to generate a mature clone population denoted as Cc(t) with the size of Nc. It is a modified single point crossover operator through which offspring is generated by randomly selecting a single crossover point on a clone and then swapping its content beyond that point with that of a parent antibody randomly selected from A(t). This operator is introduced to control the diversity of the clone population and the convergence of the algorithm. That being said, the diversity can be enhanced while the quick convergence can be ensured because some good genes from the active parent can be passed to the offspring.Suppression Operator. It works on the whole population including the parent cells and mutated clones A(t)∪Cc(t) to eliminate similar individuals so as to avoid a particular search space being over exploited and acquire the uniformly distributed Pareto front based on the idea of immune network theory [47] such that B-cells are stimulated and suppressed by not only non-self-antigens but the interacted B-cells. Different from other AIS-based multiobjective optimization algorithms, the similarity among antibodies in this algorithm is determined in terms of both the objective space and the decision variable space so as to determine whether to retain or discard individual antibodies. Therefore, the suppression operation in this algorithm has two phases: in first phase, the suppression will be applied to all antibodies and the similarity between two antibodies is defined as follows:(8)dOaba,abbj=abafj-abbfj≤δj,where dOaba,abbj is the distance between antibodies aba and abb in terms of jth objective and δj refers to the threshold value for jth objective.In this phase, if the distances for all objectives between two antibodies are smaller than the thresholds, the two antibodies are said to be similar and hence the cell with poorer Pareto fitnesspf will be suppressed and eliminated from the population.In second phase, the suppression will only be applied to the similarity between nondominated cells and dominated cells and the similarity between two antibodies is defined as follows:(9)dVaba,abb=∑j=1mabaxj-abbxj2≤ε,where dV(aba,abb) is the Euclidean distance in decision variable space between the two antibodies aba (dominated cell) and abb (nondominated cell) ∈A(t)∪Cc(t), m is the length of each antibody, and ε refers to the threshold value for the decision variable space.In this phase, if the distance between two cells is smaller than the threshold in decision variable space, the two cells are said to be similar and hence the dominated cell will be suppressed and eliminated from the population. Eventually, surviving populationsAs(t)∪Cs(t) are obtained and then enter into the selection and receptor editing process and memory updating process simultaneously.To enhance the population diversity and facilitate the search of uniformly distributed nondominated solutions, the threshold values for the decision variable space and the objective space are adaptive values that are dynamically calculated based on the maximum and minimum values found at each iteration. The adaptive threshold values are computed according to the following equations:(10)ε=abxmax-abxminn,δj=abfjmax-abfjminn,where ab(xmax) and ab(xmin) are the maximum and minimum distance values in the decision variable space, ab(fjmax) and ab(fjmin) are the maximum and minimum fitness values of the jth objective, and n is the number of antibodies in the population under investigation.Selection and Receptor Editing Operator. It works like a director to guide the search towards the promising regions by selecting all nondominated antibodies with respect to the Pareto fitnesspf from the populations As(t)∪Cs(t) to form a new population Ab(t+1) with the size of N for the next generation t+1. If the new population Ab(t+1) is not full, dominated antibodies with a better Pareto fitnesspf are selected and some genes of these antibodies are then randomly selected to be replaced by randomly generated genes. These restructured antibodies are finally added to the new population until the population is full in order to further enhance the population diversity. This process actually mimics the process of receptor editing in the biological immune system. However, if the number of nondominated antibodies found exceeds the limit of the population size N, only N nondominated antibodies with higher Ab-Ab affinity d(ab) are selected.Memory Updating Operator.It works as an elitist mechanism for helping preserve the best solutions that represent the Pareto front found over the search process. The memory set is updated at each iteration through the following procedure:(1)
If the size ofP(t) is equal to Nm, there is no memory archiving procedure to proceed.(2)
If the size ofP(t) is smaller than Nm, an archiving procedure is triggered to copy some of the dominated antibodies that have the best Pareto fitnesspf and Ab-Ab affinity d(ab) to P(t) for filling the vacant space up.(3)
If the size ofP(t) is larger than Nm, an archiving procedure is triggered to compress P(t) based on the antibody similarity (the measure used in the suppression operator) until its size becomes Nm.After performing the final memory updating procedure at the last iteration, the memory setP(t) obtained represents the resulting solution set of the algorithm.The algorithm is conducted by applying these heuristic and stochastic operators on the antibody population for balancing both the local and global search capabilities. For the detailed discussion of the hybrid multiobjective algorithm, one can refer to [45].


## 4. Multiobjective Simulation-Based Optimization Study

In this study, a set of experiments based on a real-life multiobjective optimization problem was performed to evaluate the performance and capability of the optimization framework and algorithm. In order to provide a better view on the performance of each method, both the graphical presentations of the simulation results and the statistical results of performance metrics on the problem under investigation were analyzed. All these experiments were conducted using a computer with Xeon E5-2620 2 GHz CPU with 2 GB RAM.


### 4.1. Performance Metrics

In this simulation-based optimization study, two performance metrics, namely, error ratio (ER) [48] and spacing (S) [49], are adopted to examine the quality of solution set in terms of the optimality and diversity.Error Ratio (ER).The PFtrue was originally used to compute this metric. However, since the PFtrue can hardly be found for the real-life problem under investigation in this paper, we use the reference Pareto front “PFref”, that is, the best approximation of the true Pareto front for measuring the optimality of each solution set. The PFref is found by using all of the algorithms used in this study. To achieve this, a set of nondominated solutions, that is, a Pareto front for each algorithm, is firstly generated by running a very large number of iterations, e.g., 100 over 20 trials, and then all the fronts obtained by all the trials of all the compared algorithms were merged together to form a reference Pareto front. Hence, this metric is mathematically defined as(11)ER≔∑i=1nein,ei=0,xi∈PFref,1,xi∉PFref,where n is the number of solutions in the PFknown found by the algorithm being evaluated, ei = 0 if solution i belongs to the PFref, and ei = 1 otherwise. ER=0 indicates that all the generated solutions belong to the PFref.Spacing (S).This metric is used for determining the diversity of the generated solutions. That being said, it shows how well the solutions are distributed in the PFknown. This metric is mathematically defined as(12)S≔∑i=1nd--di2n-1,di=minj⁡∑k=1mf1ix-f1jx,i,j=1,…,n,where d- is the mean of all di, k is the number of objective functions, and n is the number of solutions in the PFknown. S=0 indicates that all the nondominated solutions found are equally spaced and uniformly distributed in objective space.


### 4.2. Experimental Setup

In this section, a real-life multiobjective optimization problem, that is, the distribution operation of the material handling system (MHS) implemented at the distribution center (DC) of S.F. Express (Hong Kong) Limited, was studied through the simulation-based optimization approach. S.F. Express launched its service in Hong Kong in 1993. At present, its service network is covering almost all areas of Hong Kong, which is mainly served by the DC located in Tin Shui Wai (TSW) in northern New Territories (Figure3). Its service stores are located in 30 areas of Hong Kong [50].Figure 3
The distribution flow of S.F. Express’s imported goods in Hong Kong.The DC has two major conveyor systems: one for handling the exported goods and another for imported goods. In this study, we focus on the physical goods flow at the DC, where the items are imported from China and then distributed to all parts of Hong Kong. At the DC, most of the items received at inbound docks (Figure4) are directly transferred to outbound docks (Figure 5) and then shipped to the final destinations with little material handling in between, such as deconsolidation and sortation. As such, the expensive physical distribution functions of putaway, inventory holding, and order picking from the DC are eliminated. This approach is called cross-docking. To implement cross-docking effectively and efficiently, timely distribution of freight and better synchronization of all inbound and outbound shipments are required by making use of the information system infrastructure and advanced automated MHS in the supply chain, such as automated conveyor system, Warehouse Management System (WMS), and real-time material identification and tracking system (e.g., barcode).Figure 4
Inbound docks connecting to the conveyor system.Figure 5
Outbound docks connecting to the conveyor system.


#### 4.2.1. Current Physical Layout and Labor Deployment of the MHS

The MHS is a circular shaped automated conveyor system, which comprises a number of interconnected 5-, 10-, and 20-meter long straight modular conveyor units and 2.5-meter radius curved conveyor units. The layout of the MHS is depicted in Figure7. Each conveyor unit has a programmable logic controller to control the movement (speed, direction, acceleration, etc.) of the items being put on it and to communicate with the central computer. The conveyor system has 4 entrances connecting to 4 inbound docks and 16 exits connecting 30 outbound docks (each outbound dock serves trucks distributing parcels to one destination) (Figure 6).Figure 6
The docks and trucks at the DC.Figure 7
The simulation model of the MHS implemented with FlexSim.At each conveyor entrance, 7 workers are deployed to deconsolidate the incoming bulky consolidated parcels uploaded from the big inbound truck (16 tonnes) for facilitating the subsequent sortation process. To enable the distribution process to go well and items to be accurately sorted according to customer requirements, 4 workers are assigned to each conveyor exit serving 2 destinations (except the 2 exits close to the entrances serving 1 destination that requires 2 workers) and equipped a barcode reader for confirming the destination of each small parcel and helping to load the parcel to the small outbound truck (5.5 tonnes).


#### 4.2.2. Operation of the MHS

Each step of the operation is performed sequentially in the simulation process and the operation is given as follows:(1)
When the operation starts, each of the incoming bulky consolidated parcels is unloaded by the workers from the 4 inbound trucks to the staging area of the 4 inbound docks.(2)
And then each bulky consolidated parcel is deconsolidated into multiple small parcels and transported to the conveyor entrances by the workers.(3)
After the parcels arrive at the conveyor, they are transported to the different exits of the conveyor system according to their destinations.(4)
When the parcels reach the corresponding exits, they are moved to the outbound docks by the workers and ready for delivery.


#### 4.2.3. Assumptions

To focus our study on the behavior of the MHS, certain real-world factors were simplified. Thus, the following assumptions were made:(1)
Every day, the DC handles around 3 to 4 batches of parcels. Since the operation is similar in each batch, only one batch in the simulation model is simulated and studied. There are 400 bulky consolidated parcels in each batch coming from 4 supply sources (each contributes 100 bulky parcels) in which the processing time (cycle time) in one batch and the workers’ utilization are studied.(2)
The time for unloading bulky parcels follows a uniform distribution.(3)
The time for deconsolidation of the bulky parcels follows a normal distribution in which each bulky parcel is separated into 10 small parcels. Therefore, there are 4000 small parcels in total being handled by the MHS.(4)
The system only processes two types of parcels (namely, bulky consolidated parcels and small parcels) because they are packed similarly in terms of volume and weight.(5)
The demand for each destination is similar which follows a uniform distribution, that is, the parcels are uniformly distributed among the 30 destinations.


#### 4.2.4. Simulation Model

The simulation model with specific system configuration and behavior is implemented with FlexSim (Figure7). Details of model components and initial model settings (Table 2) are as follows:(i)
Entities are the 400 bulky consolidated parcels and the 4000 small parcels deconsolidated from the bulky parcels, which are processed through the system causing changes in the system state over time.(ii)
Activities are the deconsolidation of bulky consolidated parcels unloaded from each incoming truck (represented by a source in the simulation model), the picking of small parcels from the circular conveyor to each outgoing truck (represented by a sink in the simulation model) according to their destinations.(iii)
Item attribute is the product type associated with its destination (one product type corresponds to one specific destination).Table 2
Initial model settings.
Item Value Conveyor speed 2.5 m/sec (limit: 1–2.5 m/sec)  Conveyor spacing 1 parcel  Number of workers deployed for each conveyor entrance 7 workers (limit: 1–9 workers)  Number of workers deployed for each conveyor exit (serving 1 destination) 2 workers (limit: 1–4 workers)  Number of workers deployed for each conveyor exit (serving 2 destinations) 4 workers (limit: 1–6 workers)  Handling Capacity of Worker 1 parcel  Arrival pattern for each source Uniform distribution with a min. of 5 sec and a max. of 10 sec  Processing time of deconsolidation process Normal distribution with a mean of 30 sec and a standard deviation of 2 sec  Demand for each destination Uniform distribution with a min. of 220 units and a max. of 280 units The above model parameters are set based on the real system settings and observation.


#### 4.2.5. Parameter Setting

Several parameters, including number of generations, initial population size, active population size, maximum clone size, and mutation factor are studied through sensitivity analysis so as to examine the individual parameters’ impact on the overall performance of the system and to determine the value for each parameter. Given with the results of the sensitivity analysis, the parameters of SCMIA were set as follows.SCMIA.Initial population size is N=30; size of active population, NA = 12; size of the memory population, Nm = 30; maximum number of clones for each cell, max_clone = 6; exponential distribution coefficient, ρ=0.05; number of simulation replications per fitness evaluation, replication = 10.To allow a fair comparison among the algorithms compared, the parameters of the benchmarking algorithms were set with similar values and the values suggested by the authors.MISA.Initial population size is N=30; size of clone population, Nc = 180; size of the memory population, Nm = 30; number of grid subdivisions, subd_size = 25; initial mutation rate,ω= 0.6 (it decreases linearly over time until reaching the rate of 1/m, where m is the number of decision variables); number of simulation replications per fitness evaluation, replication = 10.NNIA.Size of dominant population is Np = 30; size of active population is NA = 8; size of clone population is Nc = 30; crossover probability is pc = 0.9; mutation probability is pm = 1/m; distribution indexes for crossover and mutation operators are nc and nm = 20 (crossover probability, pc = 1; mutation probability, pm = 1/m); number of simulation replications per fitness evaluation is replication = 10.NSGA-II.Initial population size is N=30; crossover probability is pc = 0.9; mutation probability is pm = 1/m; distribution indexes for crossover and mutation operators are nc and nm = 20; number of simulation replications per fitness evaluation is replication = 10.SPEA2.Initial population size is N=30; archive size, Nm = 30; crossover probability, pc = 1; mutation probability, pm = 0.006; number of simulation replications per fitness evaluation, replication = 10.Antibody Definition.An antibody ab (a set of decision variables) that has the direct impact on the system’s performance in terms of cycle time (CT) and workers’ utilization (WU) is defined as follows: x1 is taken to be the conveyor speed, x2 is the number of workers deployed for each conveyor entrance, x3 is the number of workers deployed for each conveyor exit (serving 2 destinations), and x4 is the number of workers deployed for each conveyor exit (serving 1 destination). Since the objective of the study is to optimize the performance of the MHS by minimizing the system cycle time and maximizing the workers’ utilization, the optimization problem can be given by(13)f→abab∈ΩOptimize=ECT,WU,ωsubject to(14)1≤x1≤2.5,1≤x2≤9,1≤x3≤6,1≤x4≤4, where a set of objective functions f→(ab) to be optimized in objective space are the expected values of the random output variables [CT,WU,ω] that are obtained from running the simulation model,ωis a sample path (i.e., the sequence of random numbers used in a simulation run), and (14) define a set of physical constraints.


#### 4.2.6. Experimental Results and Analysis

We conducted two experiments to evaluate the performance of the SCMIA and the simulation-based optimization framework, that is,(1) to compare the results of integrating simulation and optimization with the results without using any optimization algorithm and (2) to benchmark SCMIA against two immune-inspired algorithms, MISA [20] and NNIA [42], and two other evolutionary algorithms, NSGA-II [35] and SPEA2 [37], under the framework. The first experiment was conducted to evaluate the effectiveness of optimization algorithms when being employed in the simulation study while the second one was used to emphasize the comparison among the algorithms under investigation with respect to the optimality and the diversity using the same simulation-based optimization framework. All algorithms were run for 30 generations over 30 trials to obtain the average performance of each algorithm on the same condition.(1) Simulation without Optimization versus Simulation with Optimization.The results shown in Table 3 are the optimized results obtained by making use of all optimization algorithms studied in this research.Table 3
Performance comparison between the results obtained from the simulation alone and that from the simulation-based optimization approach by using different optimization algorithms (the best results are bolded).
Cycle Time (the improvement in % compared with the one without optimization) Workers’ Utilization (the improvement in % compared with the one without optimization) Simulation without optimization 6788.64 sec 45.04% Simulation-based optimization with SCMIA 5763.86 sec (15.10%) 65.46% (45.34%) Simulation-based optimization with MISA 5754.88 sec (15.23%) 63.45% (40.87%) Simulation-based optimization with NNIA 5672.87 sec (16.44 % ) 62.00% (37.66%) Simulation-based optimization with NSGA-II 5680.23 sec (16.33%) 62.90% (39.65%) Simulation-based optimization with SPEA2 5948.53 sec (12.38%) 67.10 % (48.76 % )The table shows that the cycle time of the whole distribution system at the DC is reduced by about 12–16% and the workers’ utilization is increased by 38–49% when optimization algorithms are deployed in the simulation process. This proves that the use of the optimizers can enhance the performance in the system’s cycle time and the workers’ utilization. However, the higher the utilization achieved, the longer the cycle time spent, and vice versa. When comparing SCMIA with other benchmark algorithms, SCMIA is able to produce comparable results in both the cycle time and workers’ utilization.(2) Performance Comparison between SCMIA and Other Benchmark Algorithms.The performance of the algorithms studied in this study was compared by using the graphical representation together with the application of the metrics mentioned in Section 4.1. The comparison between the known Pareto fronts shown in Figure 8 suggests that the overall Pareto front patterns generated by the five algorithms are largely similar in which the cycle time ranges from around 5300 sec to 7000 sec and the workers’ utilization ranges from around 50% to 75%. From Figure 8, it is shown that the higher the utilization achieved, the longer the cycle time spent, and vice versa. This implies that increasing the utilization does not mean that the efficiency of the system can be increased. Therefore, according to the figure, although the optimal solution with 75% utilization while spending almost 7000 sec and another case achieving less than 50% utilization while spending only around 5300 sec both fall within the reference Pareto front PFref; the latter case is preferred in practice because the cycle time is much shorter and hence the efficiency and productivity of the system are much higher. The utilization becoming lower in the latter case is mainly due to the increased values in the decision variables, namely, conveyor speed and number of workers. This implies that in order to further enhance both objectives at the same time for the DC, the company may need to do something other than just changing the system parameters, such as redesigning the layout of the material handling systems particularly the conveyor system.Figure 8
Graphical comparisons of the known Pareto fronts generated by the five algorithms.The performance regarding the optimality and diversity of SCMIA in this multiobjective optimization problem was then examined by applying the metrics, namely, spacing and error ratio as shown in Table4.Table 4
Spacing and error ratio values generated by the five algorithms (the best results are bolded).
SCMIA MISA NNIA NSGA-II SPEA2 Mean (Standard Deviation) Error Ratio (ER) 0.71 (0.10) 0.85 (0.16) 0.78 (0.18) 0.76 (0.14) 0.74 (0.06) Spacing (S) 38.74 (10.97) 58.02 (55.82) 49.23 (26.10) 51.85 (17.58) 39.32 (37.16)In this experiment, we compared the results of the mean and standard deviation of the two metrics over 30 trials obtained by SCMIA with that of the other benchmark algorithms. From the results shown in Table4, we found that SCMIA generally is able to provide the best results in terms of the diversity and optimality because it generates the lowest values in the metrics of ER (0.71) and S (38.74) and the latter metric is significantly lower than most of the other algorithms. This implies that the generated front is very close to the PFref. In terms of the stability, SCMIA is also the best one among these five algorithms in error ratio and spacing because it has much lower standard deviations (0.10) and (10.97), respectively, than other algorithms except SPEA2, implying that SCMIA is able to provide a relatively consistent result for each trial.(3) Comparison of Computational Time.The computational times of the investigated algorithms are provided in Table 5, which are the mean times of 30 trials of each algorithm running for 30 generations. It can be observed that the computational time of MISA is the longest (104.19 hours) because MISA spends very much time in solution evaluation through running the complex simulation model due to the largest number of candidate solutions being generated and evaluated at each iteration. This is very time-consuming. SCMIA comes second (33.75 hours), although the performance in terms of computational time is much better than MISA. The reason is that SCMIA also generates a large number of solutions at each iteration, but the suppression operator adopted in the algorithm helps reduce the number of similar solutions and hence reduce the number of solution evaluations. On the contrary, NNIA, NSGA-II, and SPEA2 consume less and similar computational resources largely due to the much smaller number of solution evaluations being conducted at each iteration.Table 5
Computational times of the five algorithms.
SCMIA MISA NNIA NSGA-II SPEA2 Time (hours) 33.75 104.19 24.46 24.54 23.75


## 4.1. Performance Metrics

In this simulation-based optimization study, two performance metrics, namely, error ratio (ER) [48] and spacing (S) [49], are adopted to examine the quality of solution set in terms of the optimality and diversity.Error Ratio (ER).The PFtrue was originally used to compute this metric. However, since the PFtrue can hardly be found for the real-life problem under investigation in this paper, we use the reference Pareto front “PFref”, that is, the best approximation of the true Pareto front for measuring the optimality of each solution set. The PFref is found by using all of the algorithms used in this study. To achieve this, a set of nondominated solutions, that is, a Pareto front for each algorithm, is firstly generated by running a very large number of iterations, e.g., 100 over 20 trials, and then all the fronts obtained by all the trials of all the compared algorithms were merged together to form a reference Pareto front. Hence, this metric is mathematically defined as(11)ER≔∑i=1nein,ei=0,xi∈PFref,1,xi∉PFref,where n is the number of solutions in the PFknown found by the algorithm being evaluated, ei = 0 if solution i belongs to the PFref, and ei = 1 otherwise. ER=0 indicates that all the generated solutions belong to the PFref.Spacing (S).This metric is used for determining the diversity of the generated solutions. That being said, it shows how well the solutions are distributed in the PFknown. This metric is mathematically defined as(12)S≔∑i=1nd--di2n-1,di=minj⁡∑k=1mf1ix-f1jx,i,j=1,…,n,where d- is the mean of all di, k is the number of objective functions, and n is the number of solutions in the PFknown. S=0 indicates that all the nondominated solutions found are equally spaced and uniformly distributed in objective space.


## 4.2. Experimental Setup

In this section, a real-life multiobjective optimization problem, that is, the distribution operation of the material handling system (MHS) implemented at the distribution center (DC) of S.F. Express (Hong Kong) Limited, was studied through the simulation-based optimization approach. S.F. Express launched its service in Hong Kong in 1993. At present, its service network is covering almost all areas of Hong Kong, which is mainly served by the DC located in Tin Shui Wai (TSW) in northern New Territories (Figure3). Its service stores are located in 30 areas of Hong Kong [50].Figure 3
The distribution flow of S.F. Express’s imported goods in Hong Kong.The DC has two major conveyor systems: one for handling the exported goods and another for imported goods. In this study, we focus on the physical goods flow at the DC, where the items are imported from China and then distributed to all parts of Hong Kong. At the DC, most of the items received at inbound docks (Figure4) are directly transferred to outbound docks (Figure 5) and then shipped to the final destinations with little material handling in between, such as deconsolidation and sortation. As such, the expensive physical distribution functions of putaway, inventory holding, and order picking from the DC are eliminated. This approach is called cross-docking. To implement cross-docking effectively and efficiently, timely distribution of freight and better synchronization of all inbound and outbound shipments are required by making use of the information system infrastructure and advanced automated MHS in the supply chain, such as automated conveyor system, Warehouse Management System (WMS), and real-time material identification and tracking system (e.g., barcode).Figure 4
Inbound docks connecting to the conveyor system.Figure 5
Outbound docks connecting to the conveyor system.


### 4.2.1. Current Physical Layout and Labor Deployment of the MHS

The MHS is a circular shaped automated conveyor system, which comprises a number of interconnected 5-, 10-, and 20-meter long straight modular conveyor units and 2.5-meter radius curved conveyor units. The layout of the MHS is depicted in Figure7. Each conveyor unit has a programmable logic controller to control the movement (speed, direction, acceleration, etc.) of the items being put on it and to communicate with the central computer. The conveyor system has 4 entrances connecting to 4 inbound docks and 16 exits connecting 30 outbound docks (each outbound dock serves trucks distributing parcels to one destination) (Figure 6).Figure 6
The docks and trucks at the DC.Figure 7
The simulation model of the MHS implemented with FlexSim.At each conveyor entrance, 7 workers are deployed to deconsolidate the incoming bulky consolidated parcels uploaded from the big inbound truck (16 tonnes) for facilitating the subsequent sortation process. To enable the distribution process to go well and items to be accurately sorted according to customer requirements, 4 workers are assigned to each conveyor exit serving 2 destinations (except the 2 exits close to the entrances serving 1 destination that requires 2 workers) and equipped a barcode reader for confirming the destination of each small parcel and helping to load the parcel to the small outbound truck (5.5 tonnes).


### 4.2.2. Operation of the MHS

Each step of the operation is performed sequentially in the simulation process and the operation is given as follows:(1)
When the operation starts, each of the incoming bulky consolidated parcels is unloaded by the workers from the 4 inbound trucks to the staging area of the 4 inbound docks.(2)
And then each bulky consolidated parcel is deconsolidated into multiple small parcels and transported to the conveyor entrances by the workers.(3)
After the parcels arrive at the conveyor, they are transported to the different exits of the conveyor system according to their destinations.(4)
When the parcels reach the corresponding exits, they are moved to the outbound docks by the workers and ready for delivery.


### 4.2.3. Assumptions

To focus our study on the behavior of the MHS, certain real-world factors were simplified. Thus, the following assumptions were made:(1)
Every day, the DC handles around 3 to 4 batches of parcels. Since the operation is similar in each batch, only one batch in the simulation model is simulated and studied. There are 400 bulky consolidated parcels in each batch coming from 4 supply sources (each contributes 100 bulky parcels) in which the processing time (cycle time) in one batch and the workers’ utilization are studied.(2)
The time for unloading bulky parcels follows a uniform distribution.(3)
The time for deconsolidation of the bulky parcels follows a normal distribution in which each bulky parcel is separated into 10 small parcels. Therefore, there are 4000 small parcels in total being handled by the MHS.(4)
The system only processes two types of parcels (namely, bulky consolidated parcels and small parcels) because they are packed similarly in terms of volume and weight.(5)
The demand for each destination is similar which follows a uniform distribution, that is, the parcels are uniformly distributed among the 30 destinations.


### 4.2.4. Simulation Model

The simulation model with specific system configuration and behavior is implemented with FlexSim (Figure7). Details of model components and initial model settings (Table 2) are as follows:(i)
Entities are the 400 bulky consolidated parcels and the 4000 small parcels deconsolidated from the bulky parcels, which are processed through the system causing changes in the system state over time.(ii)
Activities are the deconsolidation of bulky consolidated parcels unloaded from each incoming truck (represented by a source in the simulation model), the picking of small parcels from the circular conveyor to each outgoing truck (represented by a sink in the simulation model) according to their destinations.(iii)
Item attribute is the product type associated with its destination (one product type corresponds to one specific destination).Table 2
Initial model settings.
Item Value Conveyor speed 2.5 m/sec (limit: 1–2.5 m/sec)  Conveyor spacing 1 parcel  Number of workers deployed for each conveyor entrance 7 workers (limit: 1–9 workers)  Number of workers deployed for each conveyor exit (serving 1 destination) 2 workers (limit: 1–4 workers)  Number of workers deployed for each conveyor exit (serving 2 destinations) 4 workers (limit: 1–6 workers)  Handling Capacity of Worker 1 parcel  Arrival pattern for each source Uniform distribution with a min. of 5 sec and a max. of 10 sec  Processing time of deconsolidation process Normal distribution with a mean of 30 sec and a standard deviation of 2 sec  Demand for each destination Uniform distribution with a min. of 220 units and a max. of 280 units The above model parameters are set based on the real system settings and observation.


### 4.2.5. Parameter Setting

Several parameters, including number of generations, initial population size, active population size, maximum clone size, and mutation factor are studied through sensitivity analysis so as to examine the individual parameters’ impact on the overall performance of the system and to determine the value for each parameter. Given with the results of the sensitivity analysis, the parameters of SCMIA were set as follows.SCMIA.Initial population size is N=30; size of active population, NA = 12; size of the memory population, Nm = 30; maximum number of clones for each cell, max_clone = 6; exponential distribution coefficient, ρ=0.05; number of simulation replications per fitness evaluation, replication = 10.To allow a fair comparison among the algorithms compared, the parameters of the benchmarking algorithms were set with similar values and the values suggested by the authors.MISA.Initial population size is N=30; size of clone population, Nc = 180; size of the memory population, Nm = 30; number of grid subdivisions, subd_size = 25; initial mutation rate,ω= 0.6 (it decreases linearly over time until reaching the rate of 1/m, where m is the number of decision variables); number of simulation replications per fitness evaluation, replication = 10.NNIA.Size of dominant population is Np = 30; size of active population is NA = 8; size of clone population is Nc = 30; crossover probability is pc = 0.9; mutation probability is pm = 1/m; distribution indexes for crossover and mutation operators are nc and nm = 20 (crossover probability, pc = 1; mutation probability, pm = 1/m); number of simulation replications per fitness evaluation is replication = 10.NSGA-II.Initial population size is N=30; crossover probability is pc = 0.9; mutation probability is pm = 1/m; distribution indexes for crossover and mutation operators are nc and nm = 20; number of simulation replications per fitness evaluation is replication = 10.SPEA2.Initial population size is N=30; archive size, Nm = 30; crossover probability, pc = 1; mutation probability, pm = 0.006; number of simulation replications per fitness evaluation, replication = 10.Antibody Definition.An antibody ab (a set of decision variables) that has the direct impact on the system’s performance in terms of cycle time (CT) and workers’ utilization (WU) is defined as follows: x1 is taken to be the conveyor speed, x2 is the number of workers deployed for each conveyor entrance, x3 is the number of workers deployed for each conveyor exit (serving 2 destinations), and x4 is the number of workers deployed for each conveyor exit (serving 1 destination). Since the objective of the study is to optimize the performance of the MHS by minimizing the system cycle time and maximizing the workers’ utilization, the optimization problem can be given by(13)f→abab∈ΩOptimize=ECT,WU,ωsubject to(14)1≤x1≤2.5,1≤x2≤9,1≤x3≤6,1≤x4≤4, where a set of objective functions f→(ab) to be optimized in objective space are the expected values of the random output variables [CT,WU,ω] that are obtained from running the simulation model,ωis a sample path (i.e., the sequence of random numbers used in a simulation run), and (14) define a set of physical constraints.


### 4.2.6. Experimental Results and Analysis

We conducted two experiments to evaluate the performance of the SCMIA and the simulation-based optimization framework, that is,(1) to compare the results of integrating simulation and optimization with the results without using any optimization algorithm and (2) to benchmark SCMIA against two immune-inspired algorithms, MISA [20] and NNIA [42], and two other evolutionary algorithms, NSGA-II [35] and SPEA2 [37], under the framework. The first experiment was conducted to evaluate the effectiveness of optimization algorithms when being employed in the simulation study while the second one was used to emphasize the comparison among the algorithms under investigation with respect to the optimality and the diversity using the same simulation-based optimization framework. All algorithms were run for 30 generations over 30 trials to obtain the average performance of each algorithm on the same condition.(1) Simulation without Optimization versus Simulation with Optimization.The results shown in Table 3 are the optimized results obtained by making use of all optimization algorithms studied in this research.Table 3
Performance comparison between the results obtained from the simulation alone and that from the simulation-based optimization approach by using different optimization algorithms (the best results are bolded).
Cycle Time (the improvement in % compared with the one without optimization) Workers’ Utilization (the improvement in % compared with the one without optimization) Simulation without optimization 6788.64 sec 45.04% Simulation-based optimization with SCMIA 5763.86 sec (15.10%) 65.46% (45.34%) Simulation-based optimization with MISA 5754.88 sec (15.23%) 63.45% (40.87%) Simulation-based optimization with NNIA 5672.87 sec (16.44 % ) 62.00% (37.66%) Simulation-based optimization with NSGA-II 5680.23 sec (16.33%) 62.90% (39.65%) Simulation-based optimization with SPEA2 5948.53 sec (12.38%) 67.10 % (48.76 % )The table shows that the cycle time of the whole distribution system at the DC is reduced by about 12–16% and the workers’ utilization is increased by 38–49% when optimization algorithms are deployed in the simulation process. This proves that the use of the optimizers can enhance the performance in the system’s cycle time and the workers’ utilization. However, the higher the utilization achieved, the longer the cycle time spent, and vice versa. When comparing SCMIA with other benchmark algorithms, SCMIA is able to produce comparable results in both the cycle time and workers’ utilization.(2) Performance Comparison between SCMIA and Other Benchmark Algorithms.The performance of the algorithms studied in this study was compared by using the graphical representation together with the application of the metrics mentioned in Section 4.1. The comparison between the known Pareto fronts shown in Figure 8 suggests that the overall Pareto front patterns generated by the five algorithms are largely similar in which the cycle time ranges from around 5300 sec to 7000 sec and the workers’ utilization ranges from around 50% to 75%. From Figure 8, it is shown that the higher the utilization achieved, the longer the cycle time spent, and vice versa. This implies that increasing the utilization does not mean that the efficiency of the system can be increased. Therefore, according to the figure, although the optimal solution with 75% utilization while spending almost 7000 sec and another case achieving less than 50% utilization while spending only around 5300 sec both fall within the reference Pareto front PFref; the latter case is preferred in practice because the cycle time is much shorter and hence the efficiency and productivity of the system are much higher. The utilization becoming lower in the latter case is mainly due to the increased values in the decision variables, namely, conveyor speed and number of workers. This implies that in order to further enhance both objectives at the same time for the DC, the company may need to do something other than just changing the system parameters, such as redesigning the layout of the material handling systems particularly the conveyor system.Figure 8
Graphical comparisons of the known Pareto fronts generated by the five algorithms.The performance regarding the optimality and diversity of SCMIA in this multiobjective optimization problem was then examined by applying the metrics, namely, spacing and error ratio as shown in Table4.Table 4
Spacing and error ratio values generated by the five algorithms (the best results are bolded).
SCMIA MISA NNIA NSGA-II SPEA2 Mean (Standard Deviation) Error Ratio (ER) 0.71 (0.10) 0.85 (0.16) 0.78 (0.18) 0.76 (0.14) 0.74 (0.06) Spacing (S) 38.74 (10.97) 58.02 (55.82) 49.23 (26.10) 51.85 (17.58) 39.32 (37.16)In this experiment, we compared the results of the mean and standard deviation of the two metrics over 30 trials obtained by SCMIA with that of the other benchmark algorithms. From the results shown in Table4, we found that SCMIA generally is able to provide the best results in terms of the diversity and optimality because it generates the lowest values in the metrics of ER (0.71) and S (38.74) and the latter metric is significantly lower than most of the other algorithms. This implies that the generated front is very close to the PFref. In terms of the stability, SCMIA is also the best one among these five algorithms in error ratio and spacing because it has much lower standard deviations (0.10) and (10.97), respectively, than other algorithms except SPEA2, implying that SCMIA is able to provide a relatively consistent result for each trial.(3) Comparison of Computational Time.The computational times of the investigated algorithms are provided in Table 5, which are the mean times of 30 trials of each algorithm running for 30 generations. It can be observed that the computational time of MISA is the longest (104.19 hours) because MISA spends very much time in solution evaluation through running the complex simulation model due to the largest number of candidate solutions being generated and evaluated at each iteration. This is very time-consuming. SCMIA comes second (33.75 hours), although the performance in terms of computational time is much better than MISA. The reason is that SCMIA also generates a large number of solutions at each iteration, but the suppression operator adopted in the algorithm helps reduce the number of similar solutions and hence reduce the number of solution evaluations. On the contrary, NNIA, NSGA-II, and SPEA2 consume less and similar computational resources largely due to the much smaller number of solution evaluations being conducted at each iteration.Table 5
Computational times of the five algorithms.
SCMIA MISA NNIA NSGA-II SPEA2 Time (hours) 33.75 104.19 24.46 24.54 23.75


## 4.2.1. Current Physical Layout and Labor Deployment of the MHS

The MHS is a circular shaped automated conveyor system, which comprises a number of interconnected 5-, 10-, and 20-meter long straight modular conveyor units and 2.5-meter radius curved conveyor units. The layout of the MHS is depicted in Figure7. Each conveyor unit has a programmable logic controller to control the movement (speed, direction, acceleration, etc.) of the items being put on it and to communicate with the central computer. The conveyor system has 4 entrances connecting to 4 inbound docks and 16 exits connecting 30 outbound docks (each outbound dock serves trucks distributing parcels to one destination) (Figure 6).Figure 6
The docks and trucks at the DC.Figure 7
The simulation model of the MHS implemented with FlexSim.At each conveyor entrance, 7 workers are deployed to deconsolidate the incoming bulky consolidated parcels uploaded from the big inbound truck (16 tonnes) for facilitating the subsequent sortation process. To enable the distribution process to go well and items to be accurately sorted according to customer requirements, 4 workers are assigned to each conveyor exit serving 2 destinations (except the 2 exits close to the entrances serving 1 destination that requires 2 workers) and equipped a barcode reader for confirming the destination of each small parcel and helping to load the parcel to the small outbound truck (5.5 tonnes).


## 4.2.2. Operation of the MHS

Each step of the operation is performed sequentially in the simulation process and the operation is given as follows:(1)
When the operation starts, each of the incoming bulky consolidated parcels is unloaded by the workers from the 4 inbound trucks to the staging area of the 4 inbound docks.(2)
And then each bulky consolidated parcel is deconsolidated into multiple small parcels and transported to the conveyor entrances by the workers.(3)
After the parcels arrive at the conveyor, they are transported to the different exits of the conveyor system according to their destinations.(4)
When the parcels reach the corresponding exits, they are moved to the outbound docks by the workers and ready for delivery.


## 4.2.3. Assumptions

To focus our study on the behavior of the MHS, certain real-world factors were simplified. Thus, the following assumptions were made:(1)
Every day, the DC handles around 3 to 4 batches of parcels. Since the operation is similar in each batch, only one batch in the simulation model is simulated and studied. There are 400 bulky consolidated parcels in each batch coming from 4 supply sources (each contributes 100 bulky parcels) in which the processing time (cycle time) in one batch and the workers’ utilization are studied.(2)
The time for unloading bulky parcels follows a uniform distribution.(3)
The time for deconsolidation of the bulky parcels follows a normal distribution in which each bulky parcel is separated into 10 small parcels. Therefore, there are 4000 small parcels in total being handled by the MHS.(4)
The system only processes two types of parcels (namely, bulky consolidated parcels and small parcels) because they are packed similarly in terms of volume and weight.(5)
The demand for each destination is similar which follows a uniform distribution, that is, the parcels are uniformly distributed among the 30 destinations.


## 4.2.4. Simulation Model

The simulation model with specific system configuration and behavior is implemented with FlexSim (Figure7). Details of model components and initial model settings (Table 2) are as follows:(i)
Entities are the 400 bulky consolidated parcels and the 4000 small parcels deconsolidated from the bulky parcels, which are processed through the system causing changes in the system state over time.(ii)
Activities are the deconsolidation of bulky consolidated parcels unloaded from each incoming truck (represented by a source in the simulation model), the picking of small parcels from the circular conveyor to each outgoing truck (represented by a sink in the simulation model) according to their destinations.(iii)
Item attribute is the product type associated with its destination (one product type corresponds to one specific destination).Table 2
Initial model settings.
Item Value Conveyor speed 2.5 m/sec (limit: 1–2.5 m/sec)  Conveyor spacing 1 parcel  Number of workers deployed for each conveyor entrance 7 workers (limit: 1–9 workers)  Number of workers deployed for each conveyor exit (serving 1 destination) 2 workers (limit: 1–4 workers)  Number of workers deployed for each conveyor exit (serving 2 destinations) 4 workers (limit: 1–6 workers)  Handling Capacity of Worker 1 parcel  Arrival pattern for each source Uniform distribution with a min. of 5 sec and a max. of 10 sec  Processing time of deconsolidation process Normal distribution with a mean of 30 sec and a standard deviation of 2 sec  Demand for each destination Uniform distribution with a min. of 220 units and a max. of 280 units The above model parameters are set based on the real system settings and observation.


## 4.2.5. Parameter Setting

Several parameters, including number of generations, initial population size, active population size, maximum clone size, and mutation factor are studied through sensitivity analysis so as to examine the individual parameters’ impact on the overall performance of the system and to determine the value for each parameter. Given with the results of the sensitivity analysis, the parameters of SCMIA were set as follows.SCMIA.Initial population size is N=30; size of active population, NA = 12; size of the memory population, Nm = 30; maximum number of clones for each cell, max_clone = 6; exponential distribution coefficient, ρ=0.05; number of simulation replications per fitness evaluation, replication = 10.To allow a fair comparison among the algorithms compared, the parameters of the benchmarking algorithms were set with similar values and the values suggested by the authors.MISA.Initial population size is N=30; size of clone population, Nc = 180; size of the memory population, Nm = 30; number of grid subdivisions, subd_size = 25; initial mutation rate,ω= 0.6 (it decreases linearly over time until reaching the rate of 1/m, where m is the number of decision variables); number of simulation replications per fitness evaluation, replication = 10.NNIA.Size of dominant population is Np = 30; size of active population is NA = 8; size of clone population is Nc = 30; crossover probability is pc = 0.9; mutation probability is pm = 1/m; distribution indexes for crossover and mutation operators are nc and nm = 20 (crossover probability, pc = 1; mutation probability, pm = 1/m); number of simulation replications per fitness evaluation is replication = 10.NSGA-II.Initial population size is N=30; crossover probability is pc = 0.9; mutation probability is pm = 1/m; distribution indexes for crossover and mutation operators are nc and nm = 20; number of simulation replications per fitness evaluation is replication = 10.SPEA2.Initial population size is N=30; archive size, Nm = 30; crossover probability, pc = 1; mutation probability, pm = 0.006; number of simulation replications per fitness evaluation, replication = 10.Antibody Definition.An antibody ab (a set of decision variables) that has the direct impact on the system’s performance in terms of cycle time (CT) and workers’ utilization (WU) is defined as follows: x1 is taken to be the conveyor speed, x2 is the number of workers deployed for each conveyor entrance, x3 is the number of workers deployed for each conveyor exit (serving 2 destinations), and x4 is the number of workers deployed for each conveyor exit (serving 1 destination). Since the objective of the study is to optimize the performance of the MHS by minimizing the system cycle time and maximizing the workers’ utilization, the optimization problem can be given by(13)f→abab∈ΩOptimize=ECT,WU,ωsubject to(14)1≤x1≤2.5,1≤x2≤9,1≤x3≤6,1≤x4≤4, where a set of objective functions f→(ab) to be optimized in objective space are the expected values of the random output variables [CT,WU,ω] that are obtained from running the simulation model,ωis a sample path (i.e., the sequence of random numbers used in a simulation run), and (14) define a set of physical constraints.


## 4.2.6. Experimental Results and Analysis

We conducted two experiments to evaluate the performance of the SCMIA and the simulation-based optimization framework, that is,(1) to compare the results of integrating simulation and optimization with the results without using any optimization algorithm and (2) to benchmark SCMIA against two immune-inspired algorithms, MISA [20] and NNIA [42], and two other evolutionary algorithms, NSGA-II [35] and SPEA2 [37], under the framework. The first experiment was conducted to evaluate the effectiveness of optimization algorithms when being employed in the simulation study while the second one was used to emphasize the comparison among the algorithms under investigation with respect to the optimality and the diversity using the same simulation-based optimization framework. All algorithms were run for 30 generations over 30 trials to obtain the average performance of each algorithm on the same condition.(1) Simulation without Optimization versus Simulation with Optimization.The results shown in Table 3 are the optimized results obtained by making use of all optimization algorithms studied in this research.Table 3
Performance comparison between the results obtained from the simulation alone and that from the simulation-based optimization approach by using different optimization algorithms (the best results are bolded).
Cycle Time (the improvement in % compared with the one without optimization) Workers’ Utilization (the improvement in % compared with the one without optimization) Simulation without optimization 6788.64 sec 45.04% Simulation-based optimization with SCMIA 5763.86 sec (15.10%) 65.46% (45.34%) Simulation-based optimization with MISA 5754.88 sec (15.23%) 63.45% (40.87%) Simulation-based optimization with NNIA 5672.87 sec (16.44 % ) 62.00% (37.66%) Simulation-based optimization with NSGA-II 5680.23 sec (16.33%) 62.90% (39.65%) Simulation-based optimization with SPEA2 5948.53 sec (12.38%) 67.10 % (48.76 % )The table shows that the cycle time of the whole distribution system at the DC is reduced by about 12–16% and the workers’ utilization is increased by 38–49% when optimization algorithms are deployed in the simulation process. This proves that the use of the optimizers can enhance the performance in the system’s cycle time and the workers’ utilization. However, the higher the utilization achieved, the longer the cycle time spent, and vice versa. When comparing SCMIA with other benchmark algorithms, SCMIA is able to produce comparable results in both the cycle time and workers’ utilization.(2) Performance Comparison between SCMIA and Other Benchmark Algorithms.The performance of the algorithms studied in this study was compared by using the graphical representation together with the application of the metrics mentioned in Section 4.1. The comparison between the known Pareto fronts shown in Figure 8 suggests that the overall Pareto front patterns generated by the five algorithms are largely similar in which the cycle time ranges from around 5300 sec to 7000 sec and the workers’ utilization ranges from around 50% to 75%. From Figure 8, it is shown that the higher the utilization achieved, the longer the cycle time spent, and vice versa. This implies that increasing the utilization does not mean that the efficiency of the system can be increased. Therefore, according to the figure, although the optimal solution with 75% utilization while spending almost 7000 sec and another case achieving less than 50% utilization while spending only around 5300 sec both fall within the reference Pareto front PFref; the latter case is preferred in practice because the cycle time is much shorter and hence the efficiency and productivity of the system are much higher. The utilization becoming lower in the latter case is mainly due to the increased values in the decision variables, namely, conveyor speed and number of workers. This implies that in order to further enhance both objectives at the same time for the DC, the company may need to do something other than just changing the system parameters, such as redesigning the layout of the material handling systems particularly the conveyor system.Figure 8
Graphical comparisons of the known Pareto fronts generated by the five algorithms.The performance regarding the optimality and diversity of SCMIA in this multiobjective optimization problem was then examined by applying the metrics, namely, spacing and error ratio as shown in Table4.Table 4
Spacing and error ratio values generated by the five algorithms (the best results are bolded).
SCMIA MISA NNIA NSGA-II SPEA2 Mean (Standard Deviation) Error Ratio (ER) 0.71 (0.10) 0.85 (0.16) 0.78 (0.18) 0.76 (0.14) 0.74 (0.06) Spacing (S) 38.74 (10.97) 58.02 (55.82) 49.23 (26.10) 51.85 (17.58) 39.32 (37.16)In this experiment, we compared the results of the mean and standard deviation of the two metrics over 30 trials obtained by SCMIA with that of the other benchmark algorithms. From the results shown in Table4, we found that SCMIA generally is able to provide the best results in terms of the diversity and optimality because it generates the lowest values in the metrics of ER (0.71) and S (38.74) and the latter metric is significantly lower than most of the other algorithms. This implies that the generated front is very close to the PFref. In terms of the stability, SCMIA is also the best one among these five algorithms in error ratio and spacing because it has much lower standard deviations (0.10) and (10.97), respectively, than other algorithms except SPEA2, implying that SCMIA is able to provide a relatively consistent result for each trial.(3) Comparison of Computational Time.The computational times of the investigated algorithms are provided in Table 5, which are the mean times of 30 trials of each algorithm running for 30 generations. It can be observed that the computational time of MISA is the longest (104.19 hours) because MISA spends very much time in solution evaluation through running the complex simulation model due to the largest number of candidate solutions being generated and evaluated at each iteration. This is very time-consuming. SCMIA comes second (33.75 hours), although the performance in terms of computational time is much better than MISA. The reason is that SCMIA also generates a large number of solutions at each iteration, but the suppression operator adopted in the algorithm helps reduce the number of similar solutions and hence reduce the number of solution evaluations. On the contrary, NNIA, NSGA-II, and SPEA2 consume less and similar computational resources largely due to the much smaller number of solution evaluations being conducted at each iteration.Table 5
Computational times of the five algorithms.
SCMIA MISA NNIA NSGA-II SPEA2 Time (hours) 33.75 104.19 24.46 24.54 23.75


## 5. Conclusion

This study applies a multiobjective simulation-based optimization framework incorporating a hybrid immune optimization algorithm SCMIA for the evaluation of the optimality of the distribution system with respect to two criteria: system cycle time and workers’ utilization through simulation modeling. Based on the results of the simulation-based optimization study, the following conclusions and implications can be drawn regarding the performance of the framework and algorithm.SCMIA generally performs better than other benchmark algorithms especially in the diversity aspect. This is largely attributed to the operators employed in the algorithm. For example, the selection operator incorporates the crowding-distance as a measure to select nondominated antibodies for undergoing the subsequent evolutionary processes so that the antibodies in less crowded regions will have a higher priority to be selected. The cloning operator and hypermutation operator are based on the same measure to generate a number of copies for exploring the solution space and bringing variation to the clone population, respectively, where less crowded individuals are given more chances for cloning and hypermutation in order to hopefully produce better offspring and increase population diversity. The crossover operator helps further enhance the diversity of the clone population and the convergence of the algorithm because some good genes from the active parent can be passed to the offspring. The suppression operator helps reduce antibody redundancy by eliminating similar individuals, hence significantly minimizing the number of unnecessary searches and increasing the population diversity. The memory updating operator takes account of the antibody similarity in terms of both the objective space and the decision variable space to formulate the memory population. As a result, SCMIA is able to generate a well-distributed set of solutions while it is a good approximation to the reference Pareto front. Other than the optimality and diversity, the stability of each algorithm can also be observed based on the standard deviations of the metrics shown in Table4. It can be seen that the standard deviation of spacing obtained by SCMIA is much smaller than the other four algorithms. In a word, SCMIA is the most stable one among the five algorithms in terms of population diversity in simulation-based optimization.The results overall demonstrate the ability of the multiobjective simulation-based optimization framework to serve as a decision support tool for helping management to effectively and efficiently find near optimal system operating parameters such as the speed of machines, the number of workers, or any other decision variables of interest in order to fulfill different objectives including system’s cycle time and workers’ unitization. As a result, significant savings in money, energy, and so forth, are achieved through the effective and efficient deployment of material handling systems and well-coordinated activities based on the optimized results. The research also sheds light into important issues of logistics system operation and management based on the case study, such as manpower allocation and design capacity of facilities.Based on the findings of the current undertaking, it is worthwhile to extend the framework to tackle other complex problems involving many objectives to be solved in an efficient and effective manner in the future. Future research could also extend this approach to solve other real-world complex business problems other than the problems associated with material handling systems so as to establish the practical value of the framework in the simulation-based optimization context.



---

*Source: 5852469-2018-05-21.xml*