A heuristic algorithm combining Pareto optimization and niche technology for multi-objective unequal area facility layout problem

Engineering Applications of Artificial Intelligence 89 (2020) 103453

Contents lists available at ScienceDirect

Engineering Applications of Artificial Intelligence journal homepage: www.elsevier.com/locate/engappai

A heuristic algorithm combining Pareto optimization and niche technology for multi-objective unequal area facility layout problem✩ Jingfa Liu a,b ,∗, Jun Liu c , Xueming Yan a,b , Bitao Peng a a

School of Information Science and Technology, Guangdong University of Foreign Studies, Guangzhou 510006, China Guangzhou Key Laboratory of Multilingual Intelligent Processing, Guangdong University of Foreign Studies, Guangzhou 510006, China c School of Computer & Software, Nanjing University of Information Science & Technology, Nanjing 210044, China b

ARTICLE

INFO

Keywords: Facility layout problem Heuristic algorithm Multi-objective optimization Pareto optimal Global optimization

ABSTRACT The unequal area facility layout problem (UA-FLP) is the problem of placing departments with different areas in a facility so that departments satisfy some given objectives and constraints. In this paper, two objectives including the material handling cost and the closeness rating are optimized. Based on the quasiphysical strategy, we introduce an extrusive elastic potential energy based on the overlapping distance between departments into the layout system. After a novel handling approach of the non-overlapping constraint formed by executing the gradient method with an adaptive step length and subsequent department deformation strategy is developed to deal with the interference among departments and between any department and the facility, the problem is first converted into an optimization problem without the non-overlapping constraint. Then, we use a new heuristic algorithm that combines the local search based on the Pareto optimization and the global optimum search based on the niche technology to obtain Pareto-optimal solutions of the problem. In the proposed heuristic algorithm, in order to overcome the shortcomings of low efficient search toward the diversity of solutions in classical Pareto optimization method, we propose a heuristic layout updating strategy and a niche technology. To improve the convergence of the algorithm to the Pareto front, a mechanism of evolution of population named a feasible layout bank in the algorithm based on the local search and global optimum search is proposed. Two sets of representative instances from the literature with the size of the problem up to 62 departments are tested. The experimental results show that the proposed heuristic algorithm is an effective method for solving the multi-objective UA-FLP.

1. Introduction The facility layout problem (FLP) concerns the most efficient nonoverlapping arrangement of n indivisible departments (such as machines, employee workstations, utilities, customer service areas, restrooms, material storage areas, lunchrooms, offices and so on) in a facility in order to satisfy one or more objectives (Meller and Gau, 1996). In manufacturing organizations, departments usually represent the largest and most expensive assets of the organization and are of crucial importance to the organization. By optimizing the layout design of a facility, a significant reduction in material handling can be realized. It is estimated that between 20% and 50% production costs can be attributed to department planning and material handling, and an effective layout of a facility will help the company to reduce such costs by 10%–30% (Tompkins, 2003). The FLP is one of the most important classical problems of production management and industrial engineering, attracting the attention of

many researchers during recent decades. Broadly, the FLP can be classified into two categories: equal area facility layout problem (EA-FLP) and unequal area facility layout problem (UA-FLP) (Hosseini-Nasab et al., 2018). The EA-FLP focuses on the relative location of equal-area departments, and can be modeled as the quadratic assignment problem (QAP) of assigning equal-area departments to a number of given locations. In most applications and real-world scenarios, equal-area departments are a very poor assumption (Lacksonen, 1994). At present, the majority of approaches focus on the UA-FLPs. However, most of them study either the quantitative single objective (material handling) UA-FLPs or multi-objective UA-FLPs which are usually converted into a single objective problem by using a weighted-sum method, where it is difficult to determine the weight coefficients and standardize objective functions. In comparison, this paper involves the investigation of multiobjective UA-FLPs based on Pareto optimization method that does not need to convert multiple objectives into a single scalar objective. In

✩ No author associated with this paper has disclosed any potential or pertinent conflicts which may be perceived to have impending conflict with this work. For full disclosure statements refer to https://doi.org/10.1016/j.engappai.2019.103453. ∗ Corresponding author. E-mail address: [email protected] (J.F. Liu).

https://doi.org/10.1016/j.engappai.2019.103453 Received 12 August 2018; Received in revised form 22 November 2019; Accepted 23 December 2019 Available online xxxx 0952-1976/© 2019 Elsevier Ltd. All rights reserved.

J.F. Liu, J. Liu, X. Yan et al.

Engineering Applications of Artificial Intelligence 89 (2020) 103453

that combine the FBS or the STR to construct solutions have good performance and outperform other methods. Although the approaches mentioned above have made great progress in the UA-FLP, most of the relevant researches focus only on material handling (MH) costs as an objective function and neglect other important criteria conflicting with it, such as the closeness rating, hazardous movement, safety, and so on. Consequently, to be more realistic, some qualitative factors that have significant influence on the UA-FLP should be also given consideration. When the UA-FLP involves more than one objective, the FLP is treated as a multi-objective optimization (MOO) problem. General methods solving the multiobjective UA-FLP are to convert multiple objectives to a single scalar objective using weighted coefficients (Drira et al., 2007; Liu et al., 2017b). The poor practicability of weighted sum approach is due to the difficulty of normalizing these functions and of quantifying the weights in advance. To overcome this weakness, Pareto optimization (Deb, 2001; Ripon et al., 2009) has become an alternative to the classical weighted sum method. Aiello et al. (2013) employed a multi-objective genetic algorithm (MOGA) and Electre method for the UA-FLP, and the proposed approach was capable of finding a set of Pareto-optimal layouts that optimized the several objective functions simultaneously. Şahin and Türkbey (2009) presented a new hybrid heuristic based on the simulated annealing (SA) algorithm to find the Pareto-optimal solutions for the multi-objective FLP. Ripon et al. (2013) proposed an adaptive VNS based on evolutionary approach for optimizing the material handling cost and the closeness rating. The VNS is an explorative local search method whose basic idea is systematic change of neighborhood within a local search. Liu et al. (2018) proposed a multiobjective particle swarm optimization algorithm based on objective space division for the UA-ULP with multiple objectives. The objectives of the problem aim to optimize the material handling cost, the total adjacency value and the utilization ratio of the shop floor. Jolai et al. (2012) implemented a multi-objective particle swarm optimization to find near optimal solutions for unequal sized dynamic facility layout problem with pickup/drop-off locations. They used two heuristics of shifting the overlapped departments and condensing departments to prevent overlapping of departments. Saraswat et al. (2015) designed a multi-objective meta-heuristic framework that integrated sequence pair representation, 𝜀-accurate model and multi-objective SA algorithm to solve the FLP with three objectives (flow-distance, average workin-process, and the number of required material handling devices). At present, for multi-objective UA-FLPs, although all kinds of MOO methods based on Pareto optimization have been put forward, they still have some shortcomings, such as slow convergence to the Pareto front and low efficient search toward diversity of solutions. Therefore, further improvement and development are still awaited to enhance their effectiveness. For a review of solution methodologies for the UAFLPs including multi-objective UA-FLPs, one can see Anjos and Vieira (2017). In trying to solve the UA-FLPs considered as constrained optimization problems, the non-overlapping constraint handling is another key issue apart from the optimization technique. Constraint handling approaches of the non-overlapping among departments include mainly the flexible bay structure (FBS) (Ulutas and Kulturel-Konak, 2012; Wong and Komarudin, 2010), the slicing tree representation (STR) (Komarudin and Wong, 2010; Chang and Ku, 2013), and penalty function method (PFM) (Tate and Smith, 1995). In FBS formulation, the facility is divided into parallel horizontal or vertical bays where the width of each bay is flexible and dependent on the total area of the departments that are located in that bay. Each bay does not need to have the same width and the same number of departments. In addition, the number of bays and the number of departments in each bay are also flexible. The algorithm tries to find the optimum value for the number of bays, the number of departments contained in each bay, and the department placement order which could minimize the objective function values. Therefore, the problem becomes simpler and easier

addition, in order to overcome the shortcomings of low efficient search toward the diversity of Pareto-optimal solutions in classical Pareto optimization method, we propose a heuristic layout updating strategy and a niche technology. To improve the convergence of the algorithm to the Pareto front, a mechanism of evolution of population named a feasible layout bank in the algorithm based on local search and global optimum search is proposed. UA-FLPs are complex and NP-hard, thus, various kinds of heuristics and meta-heuristics have been proposed to solve them over the past several decades. Armour and Buffa (1963) were the first to state the UA-FLP and proposed a QAP approach by splitting the departments into small blocks with same areas, where the UA-FLP was converted to the EA-FLP. Hereafter, many researchers have studied the various types of the UA-FLPs and put forward different approaches. Bozer and Meller (1997) researched the distance-based facility layout problem, and proposed an alternative distance measure that is based on the expected distance between two departments. Komarudin and Wong (2010) proposed an ant system (AS) which used the slicing tree representation (STR) to easily represent the problems. Their algorithm incorporated several types of local search to improve its search performance and showed encouraging results in solving UA-FLPs. Chang and Ku (2013) developed a heuristic method by combining the STR and a quadratic constraint program model with the harmony search (HS) to solve the UA-FLP. Guan and Lin (2016) proposed a hybrid algorithm based on variable neighborhood search (VNS) and ant colony optimization to solve the single row UA-FLP. The genetic algorithm (GA) has been widely applied to solve the UA-FLPs due to its capability of strong global search (Aiello et al., 2006; Wu and Appleton, 2002). Tate and Smith (1995) presented a GA for solving the UA-FLP using a dynamic or adaptive penalty function to guide the search into feasible solution regions. Gonçalves and Resende (2015) established a linear programming model and proposed a biased random-key GA for the UA-FLP. Liu and Meller (2007) proposed a new formulation to solve the UA-FLP based on the sequence-pairs representation by using the GA and the mix-integer programming (MIP). Scholz et al. (2009) presented a slicing tree-based Tabu search (TS) for the continual plane FLP. They incorporated the possibility to specify various requirements regarding shape and dimensions of each individual department by using bounding curves. Nordin et al. (2009) proposed a hybridized meta-heuristic method for the solution of the UA-FLP by using the hybridization of the GA and the SA. Paes et al. (2017) introduced two algorithmic approaches to address the UA-FLP: a basic GA and a GA combined with a decomposition strategy via partial solution deconstructions and reconstructions. Their algorithms had been tested on classical benchmark instances from the literature, and returned highquality solutions. García-Hernández et al. (2014) proposed a hybrid system by incorporating human expert knowledge into the FLP. A subset of department designs was generated via the GA and then evaluated by human expert’s knowledge. Kulturel-Konak and Konak (2015) presented a large-scale hybrid SA algorithm to solve four different UA-FLPs on the continuous plane. Asl and Wong (2017) developed a modified PSO algorithm which applied two local search methods and the facility swapping method to improve the quality of solutions and prevent local optima. Palomo-Romero et al. (2017) proposed an island model genetic algorithm (IMGA) to solve the UA-FLP, and the parallel evolution of several populations was used to maintain the population diversity and to obtain a wider sampling of the search space. Vitayasak et al. (2017) focused on solving facility layout problems with stochastic demand using GAs and the backtracking search algorithm (BSA). The BSA was a population-based iterative evolutionary algorithm designed to achieve global optimization. Ulutas and Kulturel-Konak (2012) proposed a clonal selection algorithm (CSA) with a new encoding and a novel procedure to copy with dummy departments to solve the UA-FLP with the flexible bay structure (FBS). Kang and Chae (2017) presented a HS-based heuristic algorithm to solve UA-FLPs, and a STR was used as a form of layout structure. In general, the meta-heuristic algorithms 2

J.F. Liu, J. Liu, X. Yan et al.

Engineering Applications of Artificial Intelligence 89 (2020) 103453

fulfill the maximum aspect ratio constraint and the minimum side length restriction for the dimension (i.e., length and width) of each department. Suppose that the origin of Cartesian coordinate system is located at the bottom-left corner of the rectangular facility. Departments may be placed in facility F horizontally or vertically. First, we present some notations for describing the mathematical model of the UA-FLP.

to be solved. However, this rule makes the FBS a restricted layout representation by limiting the number of possible layouts. The slicing tree is a binary tree which is used to represent a slicing structure. For a FLP with N departments, the slicing tree consists of N leaves and N -1 internal nodes, where each leaf represents a department, and each internal node contains information about the direction of cut (horizontal or vertical). Generally, the STR can gain feasible layouts. However, for complex UA-FLPs with large size, the effectiveness of the STR is reduced. In the PFM for handling equality or inequality constraints in optimization problems, the fitness function is defined as the sum of the objective function and a penalty term which depends on the constraint violation. In most approaches solving the UA-FLPs, the PFM prevents from choosing infeasible solutions by considering the penalty for infeasible solutions. The PFM has been the most popular approach, because of their simplicity and ease of implementation. However, the most difficult aspect of the PFM is to find appropriate penalty parameters needed to guide the search towards the constrained optimum. Thus, its performance is not always satisfactory. In this paper, we use a heuristic algorithm that combines the local search based on the Pareto optimization and the global optimum search based on the niche technology to solve the UA-FLP, objectives of which include optimizing the material handling cost and the closeness rating score simultaneously. The main contributions of this paper are embodied in three aspects. (1) A mechanism of evolution of population named a feasible layout bank in the algorithm based on local search and global optimum search is proposed to improve the convergence of the algorithm to Pareto front. (2) A heuristic layout updating strategy and a niche technology are presented to add the diversity of the obtained solutions. (3) A novel non-overlapping constraint handling approach formed by executing the gradient method with an adaptive step length and subsequent department deformation strategy is developed to deal with the interference among departments and between any department and the facility. Furthermore, in application, two groups of instances taken from the literature with the size of the problem up to 62 departments are performed. The numerical results show that the proposed algorithm can efficiently find a set of Pareto-optimal solutions for the UA-FLP with multiple objectives and has better convergence and diversity than other approaches in the literature. The structure of this paper is organized as follows. Section 2 describes the statement of the UA-FLP and Section 3 gives the quasiphysical model of the problem. The heuristic algorithm for searching a set of feasible Pareto-optimal solutions is presented in Section 4. The experimental results and analyses are shown in Section 5. The conclusion and future research are presented in Section 6.

Indexes i, j

Indexes for departments, i, 𝑗 = 1, 2, . . . , N, where N is the number of departments

Parameters 𝐶𝑖𝑗

Cost of a unit material handling per unit distance between departments 𝑅𝑖 and 𝑅𝑗

𝑓𝑖𝑗

Frequency of material flow between departments 𝑅𝑖 and 𝑅𝑗

L W 𝑙𝑖 𝑤𝑖 𝑑𝑖𝑗

Length of facility F Width of facility F Actual length of department 𝑅𝑖 Actual width of department 𝑅𝑖 Euclidean distance between the centers of departments 𝑅𝑖 and 𝑅𝑗

Variables (𝑥𝑖 , 𝑦𝑖 ) Coordinates of the center of department 𝑅𝑖 𝑆𝑖𝑗 , 𝑇𝑖𝑗 Two binary variables used in non-overlapping constraints According to the above descriptions, a mathematical model (called model I) of the UA-FLP can be defined as follows: minimize

𝑁 𝑁 ∑ ∑

𝑓1 (𝑋) =

𝐶𝑖𝑗 𝑓𝑖𝑗 𝑑𝑖𝑗

(1)

𝑟𝑖𝑗

(2)

𝑖=1 𝑗=1

maximize

𝑓2 (𝑋) =

𝑁 ∑ 𝑁 ∑ 𝑖=1 𝑗=1

Subject to 0.5(𝑤𝑖 + 𝑤𝑗 ) − (𝑥𝑖 − 𝑥𝑗 ) ≤ 𝐿(𝑆𝑖𝑗 + 𝑇𝑖𝑗 ),

∀𝑖, 𝑗 ∈ {1, 2, … , 𝑁}

0.5(𝑤𝑖 + 𝑤𝑗 ) − (𝑥𝑗 − 𝑥𝑖 ) ≤ 𝐿(1 − 𝑆𝑖𝑗 + 𝑇𝑖𝑗 ),

2. Problem statement

∀𝑖, 𝑗 ∈ {1, 2, … , 𝑁}

(4)

0.5(𝑙𝑖 + 𝑙𝑗 ) − (𝑦𝑖 − 𝑦𝑗 ) ≤ 𝑊 (1 + 𝑆𝑖𝑗 − 𝑇𝑖𝑗 ),

∀𝑖, 𝑗 ∈ {1, 2, … , 𝑁}

(5)

0.5(𝑙𝑖 + 𝑙𝑗 ) − (𝑦𝑗 − 𝑦𝑖 ) ≤ 𝑊 (2 − 𝑆𝑖𝑗 − 𝑇𝑖𝑗 ),

∀𝑖, 𝑗 ∈ {1, 2, … , 𝑁}

(6)

𝑥𝑖 − 0.5𝑙𝑖 ≥ 0,

∀𝑖 ∈ {1, 2, … , 𝑁}

(7)

𝑥𝑖 + 0.5𝑙𝑖 ≤ 𝐿,

∀𝑖 ∈ {1, 2, … , 𝑁}

(8)

𝑦𝑖 − 0.5𝑤𝑖 ≥ 0,

∀𝑖 ∈ {1, 2, … , 𝑁}

(9)

𝑦𝑖 + 0.5𝑤𝑖 ≤ 𝑊 , ∀𝑖 ∈ {1, 2, … , 𝑁} max{𝑙𝑖 , 𝑤𝑖 } ≤ 𝑑, ∀𝑖 ∈ {1, 2, … , 𝑁} min{𝑙𝑖 , 𝑤𝑖 }

The UA-FLP considered is from Ripon et al. (2013). Suppose that there are N fixed-area departments 𝑅𝑖 (𝑖 = 1, 2, … , 𝑁) and a facility named F, which are regarded as rigid bodies and simplified into rectangles. In this study, two objective functions are set. The first objective function (𝑓1 ), the total material handling (MH) cost, can be represented by minimizing the sum of the product of material flow, distance, and transportation cost per unit distance for each pair of departments. The second objective function (𝑓2 ) which aims to maximize the closeness rating (CR) score is to place departments so as to utilize common materials, personnel, or utilities adjacent to one another, while separating departments for the reasons of safety, noise, or cleanliness and so on. The aim of the proposed multi-objective model is to find a tradeoff layout between the MH cost and the CR score for the UA-FLP and satisfy the following constraints: (1) all departments must be placed inside the given facility, and any two different departments have no interference. Each department must not overlap with the facility; (2) each department’s area is fixed but its length and width may vary within an allowable range. The area of the facility is larger than or equal to the sum of areas of all departments; (3) the layout must

(3)

min{𝑙𝑖 , 𝑤𝑖 } ≥ 𝑠, ⎧ CR value ⎪ 𝑟𝑖𝑗 = ⎨ ⎪ 0 ⎩

∀𝑖 ∈ {1, 2, … , 𝑁}

(10) (11) (12)

if departments 𝑅𝑖 and 𝑅𝑗 are neighbors with a common boundary, ∀𝑖, 𝑗 ∈ {1, 2, … , 𝑁}, 𝑖 ≠ 𝑗 otherwise (13)

𝑆𝑖𝑗 ∈ {0, 1}

∀𝑖, 𝑗 ∈ {1, 2, … , 𝑁}

(14)

𝑇𝑖𝑗 ∈ {0, 1}

∀𝑖, 𝑗 ∈ {1, 2, … , 𝑁}

(15)

Here X = (𝑥1 , 𝑦1 , 𝑙1 , 𝑤1 , 𝑥2 , 𝑦2 , 𝑙2 , 𝑤2 , … , 𝑥𝑁 , 𝑦𝑁 , 𝑙𝑁 , 𝑤𝑁 ) is a layout (or solution) of N departments. Constraints (3)–(6) ensure that there is no interference between any two different departments. 𝑙𝑖 and 𝑤𝑖 denote the actual length and width of department 𝑅𝑖 , respectively. Constraints (7)–(10) make sure that all departments must be placed inside the facility. Constraints (11) and (12) indicate that each of departments 3

J.F. Liu, J. Liu, X. Yan et al.

Engineering Applications of Artificial Intelligence 89 (2020) 103453

Table 1 Departments’ closeness relationship classification.

between points a and b. Thus, the extrusive elastic potential energy of the whole layout is:

Degree of closeness

Mark

CR values

Absolutely necessary Essentially important Important Ordinary Un-important Undesirable

A E I O U X

6 5 4 3 2 1

𝐸(𝑋) =

𝑁−1 ∑

𝑁 ∑

𝑢𝑖𝑗 +

𝑖=1 𝑗=𝑖+1

𝑁 ∑

𝑢𝑖𝐹 =

𝑖=1

𝑁−1 ∑

𝑁 ∑

𝑖=1 𝑗=𝑖+1

2 𝜇𝑙𝑖𝑗 +

𝑁 ∑

2 𝜇𝑙𝑖𝐹

(16)

𝑖=1

Obviously, E(X ) ≥ 0. For a facility F with a given size in advance, we randomly produce a layout X. Starting from layout X, if we can execute a certain algorithm to obtain layout X* so that E(X*) = 0, layout X* satisfies with the constraints of the non-overlapping among departments and between any department and facility F (i.e., satisfies with constraints (3)–(10)). This can be realized by applying two local search methods, the gradient method (GM) with adaptive step length and subsequent heuristic department deformation strategy (DDS), to minimize E(X ). Once the non-overlapping constraints are removed by the quasi-physical strategy, we can reformulate the original model I by the following model II:

must fulfill the maximum aspect ratio constraint and the minimum side length restriction, respectively. d and s which are positive real numbers denote the upper bound of the aspect ratio and the lower bound of the side length of each department, respectively. Expression (13) means that 𝑟𝑖𝑗 (=CR value) is computed if departments 𝑅𝑖 and 𝑅𝑗 are neighbors with a common boundary, 𝑟𝑖𝑗 = 0 otherwise. The CR value 𝑟𝑖𝑗 represents a degree of closeness relationship between departments 𝑅𝑖 and 𝑅𝑗 , and corresponds to adjacency requirement, which is not always quantifiable and is difficult to be given sometimes. In this paper, the most common quantifiable value for 𝑟𝑖𝑗 proposed by Lee et al. (2003) is used as shown in Table 1. For instance, for department 𝑅𝑖 , if we want to separate it from department 𝑅𝑗 due to the presence of noise or the absence of safety, a smaller CR value between them will be set. Thus, maximizing the CR score during the run of the algorithm will oblige department 𝑅𝑖 far away department 𝑅𝑗 and close to other departments that have closeness requirements with department 𝑅𝑖 . This will obtain a greater CR score. Binary constraints (14) and (15) for the two new variables, 𝑆𝑖𝑗 and 𝑇𝑖𝑗 , ensure that only one of constraints (3)–(6) is binding as follows: (a) Constraint (3) becomes active when 𝑆𝑖𝑗 = 𝑇𝑖𝑗 = 0, which implies that department i is placed on the right side of department j. (b) Constraint (4) becomes active when 𝑆𝑖𝑗 = 1 and 𝑇𝑖𝑗 = 0, which implies that department i is placed on the left side of department j. (c) Constraint (5) becomes active when 𝑆𝑖𝑗 = 0 and 𝑇𝑖𝑗 = 1, which implies that department i is placed above department j. (d) Constraint (6) becomes active when 𝑆𝑖𝑗 = 𝑇𝑖𝑗 = 1, which implies that department i is placed below department j.

minimize

𝑓1 (𝑋) =

𝑁 ∑ 𝑁 ∑

𝐶𝑖𝑗 𝑓𝑖𝑗 𝑑𝑖𝑗

(1)

𝑟𝑖𝑗

(2)

𝑖=1 𝑗=1

maximize

𝑓2 (𝑋) =

𝑁 ∑ 𝑁 ∑ 𝑖=1 𝑗=1

Subject to max{𝑙𝑖 , 𝑤𝑖 } ≤ 𝑑, min{𝑙𝑖 , 𝑤𝑖 } min{𝑙𝑖 , 𝑤𝑖 } ≥ 𝑠, ⎧ CR value ⎪ 𝑟𝑖𝑗 = ⎨ ⎪ 0 ⎩

∀𝑖 ∈ {1, 2, … , 𝑁} ∀𝑖 ∈ {1, 2, … , 𝑁}

(11) (12)

if departments 𝑅𝑖 and 𝑅𝑗 are neighbors with a common boundary, ∀𝑖, 𝑗 ∈ {1, 2, … , 𝑁}, 𝑖 ≠ 𝑗 otherwise (13)

4. A heuristic algorithm for UA-FLP To solve the optimization problem based on model II, we first generate randomly n initial layouts where every department satisfies with the maximum aspect ratio constraint and the minimum side length restriction. Then, starting from every layout where constraints (11) and (12) are naturally satisfied when the DDS (see Section 4.3) is applied to change the dimension of departments in an allowable range, a heuristic layout updating strategy is put forward to update the current layout. Afterwards, we incorporate the GM and the heuristic DDS to legalize the layout. Finally, the combination of the local search based on the Pareto optimization and the global optimum search based on the niche technology is presented to obtain Pareto-optimal solutions for the UA-FLP.

3. Quasi-physical model of the problem Based on the quasi-physical strategy (Liu et al., 2009, 2017a) whose working path is to find natural phenomena equivalent to the original mathematical problem in the physical world and then observe the evolution of the motion of matter in it so as to be inspired to obtain a formalistic algorithm for solving the mathematical problem, we introduce the gradient method with adaptive step length which simulates the move of the departments under the action of the extrusive elastic forces by relaxing the constraints and introducing a potential energy based on the overlapping distance to obtain feasible solutions of the optimization problem (1)–(15). Imagine all the rectangular departments and the exterior of the facility as smooth elastic entities. By the elastic mechanics, if there is interference between two different departments 𝑅𝑖 and 𝑅𝑗 , the extrusive elastic potential energy 𝑢𝑖𝑗 between them is proportional to the square of the overlapping depth 𝑙𝑖𝑗 (Liu et al., 2 (𝑖, 𝑗 ∈ {1, 2, … , 𝑁}, 𝑖 ≠ 𝑗), where 2017a) between them, i.e., 𝑢𝑖𝑗 = 𝜇𝑙𝑖𝑗 𝜇 is a physical coefficient. We set 𝜇 = 1 in this study. If there is interference between department 𝑅𝑖 and facility F, the extrusive elastic potential energy 𝑢𝑖𝐹 between them is 𝑢𝑖𝐹 = 𝜇l2𝑖𝐹 (𝑖 ∈ {1, 2, … , 𝑁}), where 𝑙𝑖𝐹 is the overlapping depth between department 𝑅𝑖 and facility F. 𝑙𝑖𝐹 is computed by the following three cases: (i) department 𝑅𝑖 intersects the edge of facility F only in the x-axis direction (Fig. 1(a)); (ii) department 𝑅𝑖 intersects the edge of facility F only in the y-axis direction (Fig. 1(b)); (iii) department 𝑅𝑖 intersects the edge of facility F in directions of both x-axis and y-axis (Fig. 1(c)). We define 𝑙𝑖𝐹 = |𝑏𝑖 𝑝0 |, |𝑐𝑖 𝑝0 |, and |𝑑𝑖 𝑝0 |, respectively, where |𝑎𝑏| is the Euclidean distance

4.1. Heuristic layout updating strategy In the proposed heuristic algorithm, we use the following heuristic layout updating strategy to update the layout and add the diversity of solutions. First, we generate randomly 2N vacant points that are inside facility F but not inside any departments. Then a department to be moved is chosen in accordance with the following two cases. (i) If there exists interference in the current layout, the department with the maximum comparative extrusive elastic energy potential 𝑀𝐶𝐸𝑃 𝑖 is chosen as the department to be moved. ∑𝑁 𝑗=1,𝑗≠𝑖 𝑢𝑖𝑗 + 𝑢𝑖𝐹 𝑀𝐶𝐸𝑃𝑖 = , 𝑖 = 1, 2, … , 𝑁. (17) 𝐵𝑖 where 𝐵𝑖 is the area of department 𝑅𝑖 . (ii) If there exists no interference in the current layout, we select the department with the largest unit material handling cost 𝑀𝐻𝐶 𝑖 . ∑𝑁 𝑗=1,𝑗≠𝑖 𝐶𝑖𝑗 𝑑𝑖𝑗 𝑓𝑖𝑗 𝑀𝐻𝐶𝑖 = ∑𝑁 , 𝑖 = 1, 2, … , 𝑁. (18) 𝑗=1,𝑗≠𝑖 𝐶𝑖𝑗 𝑓𝑖𝑗 4

J.F. Liu, J. Liu, X. Yan et al.

Engineering Applications of Artificial Intelligence 89 (2020) 103453

Fig. 1. The overlapping depth between department 𝑅𝑖 and facility F. (a) 𝑅𝑖 intersecting F only in the x-axis direction. (b) 𝑅𝑖 intersecting F only in the y-axis direction. (c) 𝑅𝑖 intersecting F in directions of both x-axis and y-axis.

department, we transfer similarly the embedded areas to rectangular regions which include the vacant regions around the department in the light of rules of department deformation. If there is not enough space to accommodate the embedded area, we first fill rectangular regions which include the vacant regions using a portion of the embedded area, and transfer the sum of the remaining area and the other embedded areas to the other rectangular regions which include the vacant regions around the department. If there is still interference between departments or between the departments and facility F, we continue to execute the GM with an adaptive step size. These two methods are generally executed within five times in turn until a legal (feasible) layout is obtained. Occasionally, if there exists the situation in which the GM and subsequent DDS happen to be unable to obtain a feasible layout within five alternative runs, we can choose randomly another department in the current layout to produce 4N new layouts according to the heuristic layout updating strategy and continue to run the GM and subsequent DDS to break the deadlock. We give a typical example to describe the process of DDS. Suppose that there is a layout (Fig. 2(a)), where department 𝑅1 intersects departments 𝑅2 and 𝑅3 , and department 𝑅7 intersects department 𝑅8 and facility F simultaneously. We first choose department 𝑅7 to perform deformation operations. We transfer the embedded area between department 𝑅7 and facility F to the rectangular region which includes the vacant region of the left of department 𝑅7 (Fig. 2(b)). Then, we choose department 𝑅1 as the deformation target according to the principle of ascending order of the overlapping depth. We transfer the embedded area between departments 𝑅1 and 𝑅2 to the rectangular region which includes the vacant region around department 𝑅1 (Fig. 2(c)). We continue similar department deformation operations between departments 𝑅1 and 𝑅3 , and the result of this deformation is shown in Fig. 2(d). Next, department 𝑅7 is chosen to execute similar department deformation operations. The result of the deformation is shown in Fig. 2(e), where there is still interference (the shaded region) between departments 𝑅7 and 𝑅9 . We perform the GM with an adaptive step size again, and finally a feasible layout (Fig. 2(f)) is obtained.

Finally, we place the center of selected department 𝑅𝑖 at every vacant point with vertical and horizontal directions. With the positions of the other departments unchanged, we gain 4N new layouts. If there is no any vacant point in the current layout, we generate randomly 2N points that are inside facility F. Perform the same operation as above and gain 4N new layouts. 4.2. Gradient method based on adaptive step length Once a new layout 𝑋 1 is generated from the current layout X by the heuristic layout updating strategy, we use a local search method, the gradient method (GM), to reduce the interference among departments or between the departments and the facility. The GM is also known as the steepest descent method. The search direction is the negative gradient direction. We use an adaptive step length (Liu and Li, 2010; Liu et al., 2016) as the search step size in this study. If the extrusive elastic energy E(𝑋 2 ) of the newly generated layout 𝑋 2 is higher than the energy E(𝑋 1 ) of the previous layout 𝑋 1 after one iteration, which means that the step size h for this step is too large, we decrease h to h*0.8; otherwise, h is kept. The closer the iteration is to the minimum, the slower the potential energy decreases. In order to shorten the running time of the GM, an ‘‘early-escape’’ strategy is introduced. If the energy E(𝑋 2 ) of the newly generated layout 𝑋 2 is rather close to the energy E(𝑋 1 ) of the previous layout 𝑋 1 , i.e., |𝐸(𝑋 2 ) − 𝐸(𝑋 1 )| < 10−4 , a local minimum may be reached. To overcome the defect of the low convergence rate in the GM, the search will leave the GM in advance. If the energy E(𝑋 2 ) is less than 1, we will continue to execute the GM until a feasible layout or the given minimum step size ℎ𝑚𝑖𝑛 = 10−4 in the GM is reached. 4.3. Heuristic department deformation strategy After executing the GM, there may be ‘‘stuck’’ phenomenon where there is still interference between departments or between the departments and the facility in the current layout, but the embedded departments cannot move due to the balance of elastic forces in both axial directions. At that moment, a department deformation strategy (DDS) (Fig. 2) is applied to eliminate the overlapping and obtain a feasible layout as far as possible. The detailed process is outlined below. We first choose the departments that intersect the exterior of the facility according to the ascending order of the overlapping depth, and subsequently deform (i.e. change their lengths and widths) them in turn. This is achieved by transferring the embedded areas between the chosen department and the facility to rectangular regions which include the vacant regions around the department. The rules of department deformation are that the elastic potential energy of layout after deformation operation is lower than that of the previous layout, the area and shape of department are kept unchanged, and the layout must satisfy the maximum aspect ratio constraint and the minimum side length restriction. Then, we select the departments that intersect other departments in the current layout according to the ascending order of the overlapping depth, and subsequently deform them in turn. For a selected

4.4. Local search based on the Pareto optimization Generally, there is not an optimal solution that can dominate all other solutions in a MOO problem. For the current solution 𝑋𝑖 , 4N layouts (solutions) 𝑋𝑗 , 𝑗 = 1, 2, . . . , 4N, are obtained by the heuristic layout updating strategy and the subsequent GM and DDS after an iterative step. The searching process employs the quality function and the distance between the solutions 𝑋𝑖 and 𝑋𝑗 (𝑗 = 1, 2, . . . , 4N ) for determining the search direction of the current solution 𝑋𝑖 . The value of the quality function is determined by Pareto dominance relationship between the obtained solutions 𝑋𝑗 (𝑗 = 1, 2, … , 4𝑁) and the current solution 𝑋𝑖 . If 𝑋𝑗 (𝑗 = 1, 2, … , 4𝑁) is a non-feasible solution, it will have little or no influence on the searching direction of solution 𝑋𝑖 so that the value of the quality function is small. If 𝑋𝑗 is a feasible solution and dominates 𝑋𝑖 , 𝑋𝑗 will be probably selected as the searching direction, which is beneficial for the algorithm to evolve 5

J.F. Liu, J. Liu, X. Yan et al.

Engineering Applications of Artificial Intelligence 89 (2020) 103453

Fig. 2. An example of department deformation. (a) A layout with 𝑅1 intersecting 𝑅2 and 𝑅3 , and 𝑅7 intersecting 𝑅8 and F. (b) Department deformation between 𝑅7 and F. (c) Department deformation between 𝑅1 and 𝑅2 . (d) Department deformation between 𝑅1 and 𝑅3 . (e) Department deformation between 𝑅7 and 𝑅8 . (f) Feasible layout after executing the GM.

in the direction of the Pareto front and feasible solutions. Consequently, the value of the quality function will be higher. For the associated solution 𝑋𝑗 in period t, the quality function 𝜃𝑗 (t ) is defined as follows: ⎧ ⎪ ⎪ 𝜃𝑗 (𝑡) = ⎨ ⎪ ⎪ ⎩

𝜆1 , 𝜆2 , 𝜆3 , 𝜆4 ,

4.5. Global optimum search based on the niche technology Depending on only the local search based on the Pareto optimization, the heuristic algorithm will lack the foresight of global search and the diversity of solutions is not easy to maintain also. In order to make up for these shortcomings, we propose further a global optimum search based on the niche technology. A Pareto solution set named as PF is used to store all non-dominated feasible solutions during the run of the algorithm in this study. The sparsest non-dominated solution in the PF will favor the current solution’s searching direction. Borrowed from the idea of the niche in the classic NSGA (Srinivas and Deb, 1995), a niche count based on the shared function is used to assign the fitness value to the non-dominated solutions in the PF that is collected in the previous searching process. Assume that there are p non-dominated solutions (𝑋1 , 𝑋2 , … , 𝑋𝑝 ) in the PF. For the current solution 𝑋𝑖 , the non-dominated solution 𝑋𝑙 with the smallest niche count signed by niche(𝑋𝑙 ) determines the searching direction of solution. The niche count of non-dominated solution 𝑋𝑙 is given by

if 𝑋𝑗 is non-feasible solution; if 𝑋𝑗 is feasible solution and 𝑋𝑖 ≻ 𝑋𝑗 ; if 𝑋𝑗 is feasible solution, and 𝑋𝑖 and 𝑋𝑗 are not dominated each other; if 𝑋𝑗 is feasible solution and 𝑋𝑗 ≻ 𝑋𝑖 . (19)

where j∈{1, 2, . . . , 4N }. The parameters 𝜆1 , 𝜆2 , 𝜆3 , and 𝜆4 (𝜆1 < 𝜆2 < 𝜆3 < 𝜆4 ) are the four levels for the quality function due to relationship between 𝑋𝑖 and 𝑋𝑗 . 𝑋𝑖 ≻ 𝑋𝑗 means solution 𝑋𝑖 dominates solution 𝑋𝑗 . The searching direction of solution 𝑋𝑖 is related to the quality function and the distance between the current solution 𝑋𝑖 and the produced solution 𝑋𝑗 . The higher the value of quality function is, and the closer the distance is, the greater the probability 𝑃𝑗 (t ) of selecting 𝑋𝑗 as the searching direction is in period t. 𝜃𝑗 (𝑡)𝛿𝑖𝑗 (𝑡) , 𝑗 = 1, 2, … , 4𝑁 𝑃𝑗 (𝑡) = ∑ 4𝑁 𝑘=1 𝜃𝑘 (𝑡)𝛿𝑖𝑘 (𝑡)

(20)

𝑛𝑖𝑐ℎ𝑒(𝑋𝑙 ) =

𝑝 ∑

𝑆(𝐷𝑙𝑗 ),

𝑙 = 1, 2, … , 𝑝.

(22)

𝑗=1,𝑗≠𝑙

where 𝜃𝑗 (t ) is the value of quality function in period t, 𝛿𝑖𝑗 (t )=1/𝐷𝑖𝑗 , and 𝐷𝑖𝑗 is the distance between the objective functions of 𝑋𝑖 and 𝑋𝑗 , and defined by √ √𝑚 √∑ 𝐷𝑖𝑗 = √ (𝑓𝑡 (𝑋𝑖 ) − 𝑓𝑡 (𝑋𝑗 ))2 (21)

where the shared function value S(𝐷𝑙𝑗 ) can be calculated using the following equation { 1 − 𝐷𝑙𝑗 ∕𝑟𝑠ℎ𝑎𝑟𝑒 , if 𝐷𝑙𝑗 < 𝑟𝑠ℎ𝑎𝑟𝑒 ; 𝑆(𝐷𝑙𝑗 ) = (23) 0, else.

𝑡=1

where 𝑟𝑠ℎ𝑎𝑟𝑒 is the given niche radius, and 𝐷𝑙𝑗 is the distance between the objective function values of 𝑋𝑙 and 𝑋𝑗 (see Eq. (21)).

where 𝑗 ∈ {1, 2, … , 4𝑁}, and m is the number of objective functions. 6

J.F. Liu, J. Liu, X. Yan et al.

Engineering Applications of Artificial Intelligence 89 (2020) 103453

The niche technology is introduced into the heuristic algorithm to measure the degree of similarity between solutions in the PF. Through using a niche count based on the shared function to assign the fitness value of solutions in the PF, the algorithm can select a sparsest nondominated solution from the PF. The niche technology can keep the diversity of solutions and improve the global searching capability of the heuristic algorithm. Experimental results verify the effectiveness of the niche technology.

optimum search in a given iterative step number. Finally, a set of Pareto-optimal solutions is returned. Comparing with traditional Pareto optimization methods, which have generally shortcomings, such as slow convergence to the Pareto front and low efficient search toward diversity of Pareto-optimal solutions, the proposed heuristic algorithm can improve the diversity of solutions by using a heuristic layout updating strategy and a niche technology. By means of the heuristic layout updating strategy, the algorithm can produce 4N feasible layouts and obtain a wider sampling of the search space while the niche technology is good for the algorithm to select the sparsest non-dominated solution from the Pareto front. To improve the convergence of the algorithm to the Pareto front, a mechanism of evolution of population (a feasible layout bank) in the algorithm based on local search and global optimum search is proposed, where the combination of the local search based on the Pareto optimization and the global optimum search based on the niche technology is good for the heuristic algorithm to enhance the global searching capability for non-dominated layouts with high quality. In addition, the traditional handling methods of non-overlapping among facilities usually adopt the FBS, the STR, and the PFM. In contrast, in the proposed heuristic algorithm, a novel non-overlapping constraint handling approach formed by executing the GD and subsequent DDS is developed to deal with the interference among departments and between any department and the facility. Experimental results in the following section justify also our heuristic strategies.

4.6. Description of the proposed heuristic algorithm In order to solve the UA-FLP with two objectives, a novel heuristic algorithm based on the Pareto optimization and the niche technology is proposed. First, we generate randomly n initial layouts (solutions) where every department satisfies with constraints of the maximum aspect ratio and minimum side length. Set the maximum iterative step number and relevant parameters. For all n layouts, we legalize the corresponding layouts by executing alternatively the GM and subsequent DDS, and obtain a feasible layout bank C. Then, identify all non-dominated layouts in C and gain the PF. Next, the search process that depends on the two strategies: the global search based on the niche technology and the local search based on the Pareto optimization will guide the search direction of each layout in the feasible layout bank C. The former is used to select a sparsest non-dominated layout from the PF and the latter is applied to select an optimal layout from a feasible layout bank TC obtained by a heuristic layout updating strategy and subsequent GM and DDS which can produce 4N new feasible layouts, where N is the number of departments. In order to enhance the global search effectiveness of the algorithm and keep the diversity of the PF, we give the former more operation chances than the latter by setting a probability parameter 𝑃 = 0.6. After each selection, the layout bank C and the PF will be updated. The above procedure is repeated until the maximum iterative step number is reached, and finally we can obtain a set of Pareto-optimal solutions. The detailed computational procedure of the proposed heuristic algorithm is outlined as follows. (1) Generate randomly n initial layouts (solutions) where every department satisfies with the maximum aspect ratio constraint and the minimum side length restriction. Set the maximum iterative step number T. Let 𝑃 = 0.6, and 𝑡 = 1. (2) For all n layouts, after legalizing the corresponding layouts by executing alternatively the GM and the DDS, we obtain a feasible layout bank, denoted by C. //suppose that |𝐶| = 𝑛. (3) Identify all non-dominated layouts in C, consisting of the PF. (4) Let 𝑖 = 1. (5) If random (0,1) ≤ P, select a non-dominated layout 𝑋𝑙 from the PF by the global optimum search based on the niche technology and go to step (9); otherwise go to step (6). (6) For layout 𝑋𝑖 in C, after implementing the heuristic layout updating strategy and subsequent GM and DDS, we obtain 4N feasible layouts, called test layout bank TC. (7) Select an optimal layout (denoted by 𝑋𝑙 also) from TC by the local search based on the Pareto optimization. (8) If the new layout 𝑋𝑙 is not dominated by any layout in the PF, add layout 𝑋𝑙 into the PF and delete the layouts that are dominated by 𝑋𝑙 ; otherwise, the PF is kept. (9) Let 𝑋𝑖 = 𝑋𝑙 and update layout bank C. (10) Let 𝑖 = 𝑖 + 1. If i ≤ n, go to step (5); otherwise go to step (11). (11) Let 𝑡 = 𝑡 + 1. If t < T, go to step (3); otherwise output the Pareto-optimal solutions in the PF. The corresponding pseudo-code for the procedure of the proposed heuristic algorithm is given in Algorithm 1. In Algorithm 1, steps 1 and 2 produce n legalized layouts called a feasible layout bank (population). Then, step 3 gains an initial non-dominated layout set PF. Steps 4 and 5 initialize parameters. Next, steps 6–22 are applied to execute the evolutions of population and PF based on local search and global

5. Experimental results and discussion To evaluate the efficiency of the proposed heuristic algorithm for the multi-objective UA-FLP, we implement the proposed algorithm in the Java programming language and run it on a PC with 1.5 GHz CPU and 2.0 GB RAM. We test two sets of benchmarks which include seventeen instances from the literature. These benchmark instances vary from small size to large size and from the layouts without spare space to the layouts with redundant space, which are typical representatives. Each instance is run 30 times independently. The experiments are conducted by using 𝑛 = 200 initial solutions and 𝑇 = 100 iterative steps (generations) for instances with up to 14 departments, and the 4 different parameter levels of quality function are 𝜆1 = 0.01, 𝜆2 = 0.1, 𝜆3 = 2, and 𝜆4 = 5. For each of instances with more than 14 departments, the number of the initial solutions and the number of iterative steps (generations) are set to 1000 and 900, respectively. The four corresponding parameter levels of quality function are 𝜆1 = 0.01, 𝜆2 = 0.2, 𝜆3 = 8, and 𝜆4 = 32. For these seventeen instances, the transportation cost of a unit material handling per unit distance 𝐶𝑖𝑗 between any two departments all is set to 1, and we have generated test data sets for CR score as done in Ripon et al. (2013). 5.1. Experimental results and comparisons The first set includes ten benchmark instances which are O7, O8 and O9 in Meller et al. (1998), VC10 in Van Camp et al. (1992), Ba12 and Ba14 in Van Camp (1989), Ab20 in Armour and Buffa (1963), SC30 and SC35 in Scholz et al. (2009), and Du62 in Dunker et al. (2003), respectively. The digit in each instance name represents the number of departments of this instance. The details of these instances including the material flow, the area of each department, the width and length of facility, the maximum aspect ratio constraint and the minimum side length constraint can be found in Komarudin and Wong (2010). For each instance, the numerical results obtained by the heuristic algorithm are listed in Table 2, in comparison with those obtained by the adaptive variable neighborhood search (VNS) (Ripon et al., 2013), the evolutionary algorithm (EA) with conventional local search (EA with LS) and without local search (EA without LS) (Ripon et al., 2011) in the context of the MH costs and CR scores. The best results are printed in bold font in Table 2. 7

J.F. Liu, J. Liu, X. Yan et al.

Engineering Applications of Artificial Intelligence 89 (2020) 103453

From Table 2, it is evident that the best and average values obtained by the heuristic algorithm for the MH cost are smaller than those obtained by the other three approaches in the literature for all ten instances, except instances O7 and Ab20. For instances O7 and Ab20, the heuristic algorithm performs slightly worse than other three algorithms for the best values of the MH cost, while for the average value of the MH cost, the result by the heuristic algorithm is slightly worse than the best result in the literature only for instance O7. The least number of departments in instance O7 and compact structure could limit the diversity and optimality of solutions. In addition, from Table 2, one can see that the best and average values of CR scores by the heuristic algorithm are also desirable, and it finds better results for all ten instances in comparison to the VNS, EA with LS, and EA without LS. About the computational cost, because of the differences in the performance of the running computers and the programming languages, the CPU time of other algorithms in the literature cannot be used in this comparison. However, Table 2 shows that the speed of our algorithm is fast and at least can obtain results in acceptable time for each instance. As Table 2 shows, we run the proposed heuristic algorithm for 60% of the scheduled generations for all test problems to assess their convergence behaviors. The experimental results suggest that at this stage, the proposed approach gets greater proportion of optimal solutions than other methods in the literature, verifying the better convergence of the proposed heuristic algorithm. To demonstrate the optimization performance of the proposed heuristic algorithm over generations, we present the optimization processes of the best and average values for MH cost (minimize) and CR value (maximize) over generations for instances VC10 and SC30 in Figs. 3 and 4. These figures illustrate that the proposed heuristic algorithm has a fast convergence speed and good stability. The best and average values for MH cost and CR value tend to be stable as the iterative step increases. The obtained Pareto-optimal solutions which are composed of the non-dominated solutions in the PF over 30 independent runs for instances Ba14 and SC30 are illustrated in Fig. 5(a) and (b), respectively,

and indicated by the triangular points. From Fig. 5, one can see that the Pareto-optimal solutions by the heuristic algorithm have both good diversity and good convergence, and form an evident front for each instance. It is worth to mention that very few benchmark problems are available for the UA-FLP with multiple objectives, particularly in the case of MH cost and CR score. In order to further test the performance of the proposed heuristic algorithm, we construct another set of instances, which are Z9 from Zhang et al. (2013), Das08, Das10, Das12 from Das (1993), SFLP8, SFLP11, SFLP20 from Asl and Wong (2017), respectively. The first objective of these seven instances in the literature is to minimize material handling cost, but their other optimization objectives do not include maximizing the closeness rating (CR) score. For convenience of test, in this paper we modify the objectives and constraints of these seven instances as those of problem (1)–(15). For these seven instances with different number of departments, the parameter settings of initial solutions, iterative steps, and the four different parameter levels of quality function are the same as above instances. Material flows of each instance, the area of each department, the width and the length of the facility could be found in the literature. The maximum aspect ratio constraint for instances Z9, Das08, SFLP8 and SFLP20 is set to 4, while the minimum side length for instances Das10, Das12 and SFLP11 is set to 1. The best and the worst values, the average values, and the standard deviations of both MH cost and CR score over 30 independent runs are shown in Table 3, waiting for the comparison of other algorithms in the future. From Table 3, one can see that the standard deviations of all seven instances are satisfactory, demonstrating that the heuristic algorithm has good stability and good convergence. To summarize the computational results, the proposed heuristic algorithm is able to optimize both MH cost and CR score successfully. The experimental results show that the over-all property of the proposed 8

J.F. Liu, J. Liu, X. Yan et al.

Engineering Applications of Artificial Intelligence 89 (2020) 103453

Table 2 Comparison of results by the heuristic algorithm with existing three approaches for the first set of ten instances considering objectives of both MH cost and CR score. Instance

MH cost (min)

Algorithm

CR score (max)

Best

Average

Best

Average

Optimal sol. after 60% generations

Time (s)

O7

Heuristic algorithm Adaptive VNS EA with LS EA without LS

103.12 98.69 98.69 98.69

109.29 102.36 102.36 111.009

54 54 54 54

36.652 33.125 33.125 30.825

85% 83% 43% 16%

2900 – – –

O8

Heuristic algorithm Adaptive VNS EA with LS EA without LS

202.29 202.72 202.72 202.72

251.157 270.125 279.69 292.602

70 70 68 64

51.428 42.328 38.00 32.267

80% 76% 21% 7%

3821 – – –

O9

Heuristic algorithm Adaptive VNS EA with LS EA without LS

198.71 201.75 201.75 201.75

209.82 294.36 304.12 444.672

100 96 94 88

82.18 68.16 56.52 39.897

82% 82% 18% 15%

4212 – – –

Heuristic algorithm Adaptive VNS EA with LS EA without LS

19267.37 19901.32 19907.75 19963.74

24082.77 24783.75 25739.038 26065.55

156 150 132 120

120.77 81.016 76.048 61.846

81% 72% 16% 0%

4193 – – –

Heuristic algorithm Adaptive VNS EA with LS EA without LS

7822.73 7982.24 7982.24 8103.85

8532.37 8591.18 8652.94 8711.847

158 128 128 110

100.03 78.09 77.62 71.25

73% 68% 11% 0%

4654 – – –

Heuristic algorithm Adaptive VNS EA with LS EA without LS

3925.01 4712.33 4712.33 4790.83

4114.22 5397.04 5410.04 5500.27

207 198 184 156

136.82 105.96 101.06 98.167

72% 70% 4% 0%

5212 – – –

Heuristic algorithm Adaptive VNS EA with LS EA without LS

4170.18 3978.24 3982.54 4015.25

5792.34 5856.25 5892.946 6131.369

222 216 210 201

112.74 108.86 104.84 99.329

64% 57% 0% 0%

6179 – – –

Heuristic algorithm Adaptive VNS EA with LS EA without LS

3553.53 3707.00 3716.44 3740.75

4101.84 4146.27 4186.04 4444.037

392 380 370 349

287.48 264.62 249.834 241.872

82% 58% 0% 0%

6702 – – –

Heuristic algorithm Adaptive VNS EA with LS EA without LS

3534.64 3636.75 3638.20 3835.78

4034.45 4190.25 4214.075 4305.899

356 340 340 337

295.38 237.65 226.96 216.043

72% 52% 0% 0%

6532 – – –

Heuristic algorithm Adaptive VNS EA with LS EA without LS

2965307.34 2975604.25 2976120.11 2977512.96

3211403.53 3218854.65 3220214.64 3220771.01

322 306 284 248

241.16 201.125 168.17 155.25

64% 56% 0% 0%

14323 – – –

VC10

Ba12

Ba14

Ab20

SC30

SC35

Du62

Fig. 3. Two objectives over generations of instance VC10. (a) MH cost. (b) CR score.

heuristic algorithm for solving the UA-FLP is superior to other three algorithms in the literature. In many cases, it can optimize multiple objectives efficiently and has better searching capability.

to be validated using statistical analysis tools. In order to measure and compare the performance of the final layouts optimizing multiple objectives, various performance metrics are defined and used. Zhou et al. (2011) suggested three desirable aspects of the obtained non-

5.2. Performance analysis

dominated front that can be measured: the convergence measurement of a Pareto optimal front, the uniformity of the obtained solutions, and

The quality of multi-objective optimization is often difficult to define precisely by any single performance metric and the results have

the range of these solutions in the objective space. The commonly used 9

J.F. Liu, J. Liu, X. Yan et al.

Engineering Applications of Artificial Intelligence 89 (2020) 103453

Fig. 4. Two objectives over generations of instance SC30. (a) MH cost. (b) CR score.

Fig. 5. Final Pareto-optimal solutions by the heuristic algorithm. (a) Instance Ba14. (b) Instance SC30.

varies from 0 to 1 and the ideal value of a population is 1. √ √ 𝑝 √ 1 ∑ 𝑆𝑃 = √ (𝑑 − 𝑑)2 𝑝 − 1 𝑖=1 𝑖

Table 3 Results by the heuristic algorithm for the second set of seven instances considering objectives of both MH cost and CR score. Instance

Object

Best

Worst

Average

Standard deviation

Z9

MH cost CR score

1723.23 85

2490.28 54

1913.13 66.21

102.32 20.72

Das08

MH cost CR score

8902.92 72

12321.23 33

9344.24 58.23

1092.61 17.44

Das10

MH cost CR score

17032.37 154

26349.11 61

20972.05 96.53

2923.12 48.23

Das12

MH cost CR score

21938.13 182

28324.52 98

24124.74 123.93

2104.86 41.26

SFLP8

MH cost CR score

188.36 80

230.84 32

210.33 53.46

15.98 12.90

SFLP11

MH cost CR score

1023.35 170

1943.23 118

1463.98 142.11

153.99 29.21

SFLP20

MH cost CR score

1123.95 280

1542.93 192

1323.47 220.12

134.72 52.32

𝑑𝑖 = min 𝑗

𝑚 ∑

|𝑓𝑡 (𝑋𝑖 ) − 𝑓𝑡 (𝑋𝑗 )|,

𝑖, 𝑗 = 1, 2, … , 𝑝, 𝑖 ≠ 𝑗.

(25)

(26)

𝑡=1

where p is the number of the obtained non-dominated solutions, m is the number of objectives, 𝑓𝑡 (𝑋𝑖 ) and 𝑓𝑡 (𝑋𝑗 ) are the values of the objective function 𝑓𝑡 (X ) of solutions 𝑋𝑖 and 𝑋𝑗 in the obtained nondominated solution set, respectively. The distance measure 𝑑𝑖 is the minimum distance between solution 𝑋𝑖 and the others from the obtained non-dominated set, and 𝑑̄ is the mean value of all 𝑑𝑖 , so SP is the standard deviation of p distances. The SP metric measures the uniformity of the obtained solutions with the standard deviations of all 𝑑𝑖 values. When the obtained solutions are distributed more evenly, the SP metric is smaller. 𝑂𝑃 𝑆 =

𝑚 ∏

𝑂𝑃 𝑆𝑡

(27)

𝑡=1

𝑂𝑃 𝑆𝑡 =

corresponding performance metrics are Pareto Ratio (PR) (Collette and Siarry, 2005), Space (SP) (Schott, 1995) and Overall Pareto Spread (OPS) (Zheng et al., 2017), respectively. These three performance metrics can be computed by the following equations. 𝑃𝑅 =

|𝐴𝑃 𝐹 (𝑆)| |𝑆|

| max𝑋∈𝑃 𝐹 𝑓𝑡 (𝑋) − min𝑋∈𝑃 𝐹 𝑓𝑡 (𝑋)| |𝑓𝑡 (𝑃𝐵 ) − 𝑓𝑡 (𝑃𝐺 )|

(28)

where 𝑃𝐵 and 𝑃𝐺 are the Utopia and Nadir design points (Makowski, 2004), which are composed of the best and worst values of objectives in the set of all the solutions, respectively. The PF consists of all nondominated solutions. When the obtained solutions are distributed more widely, the above OPS metric becomes larger. In this paper, the three performance metrics for the first set of ten benchmark instances are shown in Table 4, in comparison with results of the three algorithms in Ripon et al. (2013). The values of the metrics in Table 4 indicate that Pareto-optimal solutions obtained by the heuristic algorithm have better convergence and diversity characteristics than those by the VNS, EA with LS, and EA without LS, except instance O7. The best results are printed in bold font in Table 4. In addition,

(24)

where APF is the set of the approximated Pareto frontier which is close but different from the true Pareto frontier. |𝐴𝑃 𝐹 (𝑆)| represents the number of the solutions set S in the APF. S is the set of all the solutions obtained by the multi-objective optimization method at a given iterative step number, and |𝑆| represents the total number of solutions of S (including solutions from |𝐴𝑃 𝐹 (𝑆)|). The value of PR 10

J.F. Liu, J. Liu, X. Yan et al.

Engineering Applications of Artificial Intelligence 89 (2020) 103453 Table 6 Computational results obtained by the heuristic algorithm with different probability parameter P.

Table 4 Comparison of performance metrics of the four different approaches for the first set of ten benchmark instances. Instance

Algorithm

PR

SP

OPS

O7

Heuristic algorithm Adaptive VNS EA with LS EA without LS

0.98 1 1 0.98

0.11 0.1 0.1 2.431

0.79 0.81 0.81 0.78

O8

Heuristic algorithm Adaptive VNS EA with LS EA without LS

0.90 0.86 0.78 0.58

1.308 1.416 2.586 5.850

0.65 0.64 0.60 0.64

O9

Heuristic algorithm Adaptive VNS EA with LS EA without LS

0.89 0.86 0.77 0.64

2.632 3.014 3.054 4.721

0.65 0.60 0.61 0.59

VC10

Heuristic algorithm Adaptive VNS EA with LS EA without LS

0.78 0.77 0.62 0.62

0.727 1.049 1.390 1.700

0.71 0.61 0.51 0.50

Ba12

Heuristic algorithm Adaptive VNS EA with LS EA without LS

0.87 0.82 0.805 0.7

1.375 1.668 1.654 2.690

0.61 0.59 0.60 0.59

Ba14

Heuristic algorithm Adaptive VNS EA with LS EA without LS

0.81 0.76 0.645 0.505

1.421 1.974 2.040 2.701

0.64 0.60 0.60 0.56

Ab20

Heuristic algorithm Adaptive VNS EA with LS EA without LS

0.74 0.70 0.61 0.505

5.772 5.892 6.942 7.370

0.70 0.66 0.62 0.60

SC30

Heuristic algorithm Adaptive VNS EA with LS EA without LS

0.73 0.61 0.51 0.49

4.823 5.001 5.361 6.971

0.72 0.67 0.62 0.54

SC35

Heuristic algorithm Adaptive VNS EA with LS EA without LS

0.61 0.58 0.505 0.49

5.831 6.035 6.473 8.374

0.69 0.63 0.59 0.54

Du62

Heuristic algorithm Adaptive VNS EA with LS EA without LS

0.75 0.61 0.57 0.47

4.331 4.952 5.587 7.430

0.68 0.62 0.56 0.53

Instance

PR

SP

OPS

Z9 Das08 Das10 Das12 SFLP8 SFLP11 SFLP20

0.71 0.49 0.58 0.81 0.83 0.65 0.72

1.924 2.931 3.128 2.140 1.662 2.230 6.923

0.53 0.63 0.73 0.61 0.54 0.77 0.69

MH cost (min)

CR score (max)

Best

Avg

Best

Avg

O8

0.2 0.4 0.6 0.7 0.8

231.60 202.29 202.29 202.29 221.32

293.12 270.31 251.137 251.20 272.232

64 68 70 70 70

34.343 43.212 51.428 51.248 50.322

SC30

0.2 0.4 0.6 0.7 0.8

3623.12 3612.48 3553.53 3553.53 3671.53

4221.22 4191.76 4101.84 4102.11 4185.71

374 371 392 392 362

262.32 267.76 287.48 286.39 261.92

DU62

0.2 0.4 0.6 0.7 0.8

2979327.52 2975327.78 2965307.34 2965307.34 2966120.12

3221453.02 3215313.13 3211403.53 3211418.57 3220742.21

296 306 322 322 312

228.34 234.82 240.28 241.16 230.27

for each parameter value of P on three different typical instances O8, SC30, and Du62. The best values (Best) and the average values (Avg) of the algorithm with different parameters of P for the MH cost and CR score are listed in Table 6. Table 6 shows that for instances O8, SC30, and Du62, the best values for both objectives can always be obtained while the average values are slightly different when P is set to 0.6 and 0.7, respectively. The results of the algorithm turn worse when P is set to 0.2, 0.4 and 0.8. Overall, when P is set to 0.6 the algorithm can obtain best results, where the local search and the global optimum search make the algorithm reach most profitable trade-off between two objectives. Therefore, we set P to 0.6 in our heuristic algorithm. 6. Conclusions The FLP is one of the most important classical problems of production management and industrial engineering that has attracted the attention of many researches during recent decades. In spite of its significance, the majority of the existing methods solve either the quantitative single objective (material handling) or multi-objective FLPs using a weighted-sum method. In contrast, this paper focuses on a multi-objective unequal area FLP using a novel heuristic algorithm combining the Pareto optimization and the niche technology. The objectives of the problem include optimizing the MH (quantitative aspect) and the CR (qualitative aspect) simultaneously. We convert first the problem into an optimization problem without the non-overlapping constraint using a quasi-physical strategy and then apply the heuristic algorithm to find a set of Pareto-optimal solutions of the problem. The combination of the local search based on the Pareto optimization and the global optimum search based on the niche technology can help the algorithm to improve the searching ability for Pareto-optimal solutions. The heuristic layout updating strategy provides a wide range of alternative layouts so as to add the diversity and optimality of layouts. Specially, the combination of the gradient method based on the adaptive step length and the department deformation strategy is applied to deal with the non-overlapping constraint between departments, which is a key issue when solving the UA-FLP. In order to test the effectiveness of the proposed heuristic algorithm, two sets of instances are used as benchmarks. The experimental results indicate that the heuristic algorithm is an effective method for solving the UAFLP. In the future work, we intend to apply the proposed method to the unequal area dynamic facility layout problem.

Table 5 Performance metrics of the proposed heuristic algorithm for the second set of seven benchmark instances. Instance

P

the three performance metrics of the second set of seven benchmark instances are shown in Table 5, waiting for the comparison of other algorithms in the future. 5.3. Numerical experiments for different parameters In addition, the contrast experiments are designed for testing and verifying the effects of the probability parameter P in the heuristic algorithm that combines the local search based on the Pareto optimization and the global optimum search based on the niche technology to obtain Pareto-optimal solutions of the problem. The parameter P is set to 0.2, 0.4, 0.6, 0.7 and 0.8, respectively, while other parameters in the algorithm are fixed. The algorithm is independently run 30 times

CRediT authorship contribution statement Jingfa Liu: Conceptualization, Methodology, Investigation, Project administration, Writing - review & editing. Jun Liu: Writing - original draft, Software. Xueming Yan: Formal analysis. Bitao Peng: Visualization. 11

J.F. Liu, J. Liu, X. Yan et al.

Engineering Applications of Artificial Intelligence 89 (2020) 103453

Acknowledgments

Liu, Q., Meller, R.D., 2007. A sequence-pair representation and MIP-model-based heuristic for the facility layout problem with rectangular departments. IIE Trans. 39, 377–394. Liu, J.F., Wang, D.W., He, K., Xue, Y., 2017b. Combining Wang-Landau sampling algorithm and heuristics for solving the unequal-area dynamic facility layout problem. European J. Oper. Res. 262 (2), 1052–1063. Liu, J.F., Xue, S.J., Liu, Z.X., Xu, D.H., 2009. An improved energy landscape paving algorithm for the problem of packing circles into a larger containing circle. Comput. Ind. Eng. 57, 1144–1149. Liu, J.F., Zhang, H.Y., He, K., Jiang, S.Y., 2018. Multi-objective particle swarm optimization algorithm based on objective space division for the unequal-area facility layout problem. Expert Syst. Appl. 102, 179–192. Liu, J.F., Zhang, K.W., Yao, Y.L., Xue, Y., Guan, T.Z., 2016. A heuristic quasi-physical algorithm with coarse and fine adjustment for multi-objective weighted circles packing problem. Comput. Ind. Eng. 101, 416–426. Makowski, M., 2004. Multi-objective decision support including sensitivity analysis. In: Encyclopedia of Life Support Systems. pp. 1–24. Meller, R.D., Gau, K.Y., 1996. The facility layout problem: recent and emerging trends and perspectives. J. Oper. Res. Soc. 15, 351–366. Meller, R.D., Narayanan, V., Vance, P.H., 1998. Optimal facility layout design. Oper. Res. Lett. 23, 117–127. Nordin, N.N., Zainuddin, Z.M., Salim, S., Ponnusamy, R.R., 2009. Mathematical modeling and hybrid heuristic for unequal size facility layout problem. J. Fundam. Sci. 5, 87–89. Paes, F.G., Pessoa, A.A., Vidal, T., 2017. A hybrid genetic algorithm with decomposition phases for the unequal area facility layout problem. European J. Oper. Res. 256 (2), 742–756. Palomo-Romero, J.M., Salas-Morera, L., García-Hemández, L., 2017. An island model genetic algorithm for unequal area facility layout problems. Expert Syst. Appl. 68, 151–162. Ripon, K.S.N., Glette, K., Khan, K.N., Hovin, M., Torresen, J., 2013. Adaptive variable neighborhood search for solving multi-objective facility layout problems with unequal area facilities. Swarm Evol. Comput. 8, 1–12. Ripon, K.S.N., Glette, K., Mirmotahari, O., Hovin, O.M., Torresen, J., Pareto optimal based evolutionary approach for solving multi-objective facility layout problem, in: Proceedings of the 16th International Conference on Neural Information Processing, ICONIP, Berlin, Heidelberg, 2009, pp. 159–168. Ripon, K.S.N., Khan, K.N., Glette, K., Hovin, M., Torresen, J., Using pareto-optimality for solving multi-objective unequal area facility layout problem, in: Proceedings of the 13th Annual Conference on Genetic and Evolutionary Computation, GECCO, New York, USA, 2011, pp. 681–688. Saraswat, A., Venkatadri, U., Castillo, I., 2015. A framework for multi-objective facility layout design. Comput. Ind. Eng. 90, 167–176. Scholz, D., Petrick, A., Domschke, W., 2009. STaTS: a slicing tree and tabu search based heuristic for the unequal area facility layout problem. European J. Oper. Res. 197, 166–178. Schott, J.R., 1995. Fault Tolerant Design using Single and Multicriteria Genetic Algorithm Optimization (thesis (MS)). Massachusetts Institute of Technology, Cambridge, MA. Srinivas, N., Deb, K., 1995. Multiobjective optimization using nondominated sorting in genetic algorithms. Evol. Comput. 2, 221–248. Tate, D.M., Smith, A.E., 1995. Unequal area facility layout using genetic search. IIE Trans. 27, 465–472. Tompkins, A., 2003. Facilities Planning. John Wiley & Sons, New York. Ulutas, B.H., Kulturel-Konak, S., 2012. An artificial immune system based algorithm to solve unequal area facility layout problem. Expert Syst. Appl. 39, 5384–5395. Van Camp, D.J., 1989. A Nonlinear Optimization Approach for Solving Facility Layout Problem (thesis). University of Toronto, Canada. Van Camp, D.J., Carter, M.W., Vannelli, A., 1992. A nonlinear optimization approach for solving facility layout problems. European J. Oper. Res. 57, 174–189. Vitayasak, S., Pongcharoen, P., Hicks, C., 2017. A tool for solving stochastic dynamic facility layout problems with stochastic demand using either a genetic algorithm or modified backtracking search algorithm. Int. J. Prod. Econ. 190, 146–157. Wong, K.Y., Komarudin, K., 2010. Solving facility layout problems using flexible Bay structure representation and ant system algorithm. Expert Syst. Appl. 37, 5523–5527. Wu, Y., Appleton, E., 2002. The optimization of block layout and aisle structure by a genetic algorithm. Comput. Ind. Eng. 41 (2), 371–387. Zhang, Y., Lu, C., Zhang, H., Fang, Z.F., 2013. Workshop layout optimization based on differential cellular multi-objective genetic algorithm. Comput. Integr. Manuf. Syst. 19 (2), 727–734. Zheng, K., Yang, R.J., Xu, H., Hu, J.A., 2017. A new distribution metric for comparing Pareto optimal solutions. Struct. Multidiscip. Optim. 55, 1–10. Zhou, A., Qu, B.Y., Li, H., Zhao, S.Z., Suganthan, P.N., Zhang, Q., 2011. Multiobjective volutionary algorithms: a survey of the state of the art. Swarm Evol. Comput. 1, 32–49.

This work is supported by the Natural Science Foundation of Jiangsu Province, China (Grant No. BK20181409), the Major Program of the National Social Science Foundation of China (Grant No. 16ZDA047), the Special Foundation of Guangzhou Key Laboratory of Multilingual Intelligent Processing, China (Grant No. 201905010008), and the Program of Basic and Applied Basic Research of Guangzhou, China. References Aiello, G., Enea, M., Galante, G., 2006. A multi-objective approach to facility layout problem by genetic search algorithm and electre method. Robot. Comput.-Integr. Manuf. 22, 447–455. Aiello, G., Scalia, G.L., Enea, M., 2013. A non-dominated ranking multi objective genetic algorithm and electre method for unequal area facility layout problems. Expert Syst. Appl. 40, 4812–4819. Anjos, M.F., Vieira, M.V.C., 2017. Mathematical optimization approaches for facility layout problems: The state-of-the-art and future research directions. European J. Oper. Res. 261, 1–16. Armour, G.C., Buffa, E.S., 1963. A heuristic algorithm and simulation approach to relative allocation of facilities. Manage. Sci. 9 (2), 294–309. Asl, A.D., Wong, K.Y., 2017. Solving unequal-area static and dynamic facility layout problems using modified particle swarm optimization. J. Intell. Manuf. 28 (2), 1317–1336. Bozer, Y.A., Meller, R.D., 1997. A reexamination of the distance-based facility layout problem. IIE Trans. 29, 549–560. Chang, M.S., Ku, T.C., 2013. A slicing tree representation and QCP-model-based heuristic algorithm for the unequal-area block facility layout problem. Math. Probl. Eng. 2013, 1–19. Collette, Y., Siarry, P., 2005. Three new metrics to measure the convergence of metaheuristics towards the Pareto frontier and the aesthetic of a set solutions in biobjective optimization. Comput. Oper. Res. 32, 773–792. Şahin, R., Türkbey, O., 2009. A simulated annealing algorithm to find approximate pareto optimal solutions for the multi-objective facility layout problem. Int. J. Adv. Manuf. Technol. 41, 1006–1018. Das, S.K., 1993. A facility layout method for flexible manufacturing systems. Int. J. Prod. Res. 31 (2), 279–297. Deb, K., 2001. Multi-Objective Optimization using Evolutionary Algorithms. John Wiley & Sons. Drira, A., Pierreval, H., Hajri-Gabouj, S., 2007. Facility layout problem: A survey. Annu. Rev. Control 31, 255–267. Dunker, T., Radons, G., Westkämper, E., 2003. A coevolutionary algorithm for a facility layout problem. Int. J. Prod. Res. 41, 3479–3500. García-Hernández, L., Pérez-Ortiz, M., Araúzo-Azofra, A., Salas-Morera, L., HervásMartínez, C., 2014. An evolutionary neural system for incorporating expert knowledge into the UA-FLP. Neurocomputing 135, 69–78. Gonçalves, J.F., Resende, M.G.C., 2015. A biased random-key genetic algorithm for the unequal area facility layout problem. European J. Oper. Res. 246, 86–107. Guan, J., Lin, G., 2016. Hybridizing variable neighborhood search with ant colony optimization for solving the single row facility layout problem. European J. Oper. Res. 248 (2), 899–909. Hosseini-Nasab, H., Fereidouni, S., Ghomi, S.M.T.F., Fakhrzad, M.B., 2018. Classification of facility layout problems: a review study. Int. J. Adv. Manuf. Technol. 94 (1–4), 957–977. Jolai, F., Tavakkoli-Moghaddam, R., Taghipour, M., 2012. A multi-objective particle swarm optimization algorithm for unequal-sized dynamic facility layout problem with pickup/drop-off locations. Int. J. Prod. Res. 50, 4279–4293. Kang, S., Chae, J., 2017. Harmony search for the layout design of an unequal area facility. Expert Syst. Appl. 79, 268–281. Komarudin, Wong, K.Y., 2010. Applying ant system for solving unequal area facility layout problems. European J. Oper. Res. 202, 730–746. Kulturel-Konak, S., Konak, A., 2015. A large-scale hybrid simulated annealing algorithm for cyclic facility layout problems. Eng. Optim. 47 (2), 963–978. Lacksonen, T.A., 1994. Static and dynamic layout problems with varying areas. J. Oper. Res. Soc. 45, 59–69. Lee, K.Y., Han, S.N., Roh, M.I., 2003. An improved genetic algorithm for facility layout problems having inner structure walls and passages. Comput. Oper. Res. 30 (2), 117–138. Liu, J.F., Li, G., 2010. Basin filling algorithm for the circular packing problem with equilibrium behavioral constraints. Sci. China Inf. Sci. 53 (2), 885–895. Liu, J.F., Li, J., Lü, Z.P., Xue, Y., 2017a. A quasi-human strategy-based improved basin filling algorithm for the orthogonal rectangular packing problem with mass balance constraint. Comput. Ind. Eng. 107, 196–210.

12