A decomposition-based algorithm for the double row layout problem

Applied Mathematical Modelling 77 (2020) 963–979

Contents lists available at ScienceDirect

Applied Mathematical Modelling journal homepage: www.elsevier.com/locate/apm

A decomposition-based algorithm for the double row layout problem Jian Guan a,c, Geng Lin b,d,∗, Hui-Bin Feng a, Zhi-Qiang Ruan a a

College of Computer and Control Engineering, Minjiang University, Fuzhou 350121, China College of Mathematics and Data Science, Minjiang University, Fuzhou 350121, China c Modern Educational Technology Center, Minjiang University, Fuzhou 350121, China d Collaborative Innovation Center of IoT Industrialization and Intelligent Production, Minjiang University, Fuzhou 350121, China b

a r t i c l e

i n f o

Article history: Received 15 December 2018 Revised 4 August 2019 Accepted 14 August 2019 Available online 20 August 2019 Keywords: Decomposition Mixed integer programming Double row layout problem Local search Particle swarm optimization

a b s t r a c t The double row layout problem (DRLP) is a common pattern of facility layout problem, which has practical applications in flexible manufacturing systems. The double row layout problem is vital to save transportation cost and enhance productivity. Nevertheless, it is very hard to handle the DRLP because of its characteristic of combination of combinatorial and continuous aspects. In this paper, a decomposition-based algorithm is proposed to solve the DRLP. We decompose the DRLP into two subproblems. In the first subproblem, the adjustable clearances between adjacent facilities are temporarily ignored. A first improvement based local search is applied to optimize the sequences of facilities on double rows. During this process, the facilities of double rows are placed starting at different abscissas rather than starting at the same abscissa for each arrangement. A property of the objective function of the DRLP is used to obtain the optimal difference between two starting abscissas. In the second subproblem, a particle swarm optimization is applied to optimize the adjustable clearances between adjacent facilities under the condition that the sequences of facilities are fixed. Our proposed algorithm is evaluated on 59 test instances and compared with the state-of-the-art methods. The experimental results demonstrate the high competitiveness of our proposed algorithm. © 2019 Elsevier Inc. All rights reserved.

1. Introduction In the flexible manufacturing systems (FMSs), it is efficient to employ automated guided vehicles to handle the material flow on the production line. An automated guided vehicle moves in both directions on a straight line path, where the facilities are placed on both sides in parallel. Since there are unequal distances and material flows between pairs of facilities, the arrangement of facilities has a remarkable impact on the total cost of material handling. Decision makers always try their best to find the optimal arrangement of facilities to minimize the total cost of material handling. This is an example of the double row layout problem (DRLP) [1], which is how to assign a given set of facilities to locations on double rows so that the total cost of material handling among facilities is minimized. The double row layout problem arises in many other real-world applications, such as the arrangement of rooms in buildings [2], the layout of stacker cranes in FMSs [3], the layout of machines in semiconductor manufacturing [4] and the placement of toolsets in LCD fabrication line [5]. The ∗

Corresponding author at: College of Mathematics and Data Science, Minjiang University, Fuzhou 350121, China. E-mail address: [email protected] (G. Lin).

https://doi.org/10.1016/j.apm.2019.08.015 0307-904X/© 2019 Elsevier Inc. All rights reserved.

964

J. Guan, G. Lin and H.-B. Feng et al. / Applied Mathematical Modelling 77 (2020) 963–979

double row layout problem is a particular variant of the facility layout, which plays a key role in a manufacturing system and has a significant contribution towards manufacturing productivity in terms of cost and time [6]. Therefore, it has significant importance in engineering fields. Although the DRLP is commonly encountered in practice, it has received only limited attention in the literature. Chung and Tanchoco [5] first developed a mixed integer programming (MIP) model for the DRLP and applied five heuristic algorithms to solve it. Zhang and Murray [7] pointed out that there were errors in the MIP model formulation of [5]. These errors admitted layouts that did not observe the minimum clearance restrictions. Thus, Zhang and Murray corrected the MIP model formulation for the DRLP. Murray et al. [8] extended the corrected formulation of [7] by considering asymmetric material flows between machines. To solve the problem more efficiently, they improved the constructive heuristics of [5] by embedding a simple local search into the heuristics. Zhang and Cheng [9] decomposed the DRLP into a combinatorial optimization problem and a linear problem. Three types of heuristic and CPLEX were used to solve the combinatorial optimization problem and linear problem, respectively. Anjos et al. [10] presented integer linear programming and semidefinite programming approaches to solve a double row equidistant facility layout problem. Amaral [11] presented a new MIP model of the DRLP based on α -incidence vectors, which were helpful to reduce both the number of continuous variables and the number of binary variables. Secchin and Amaral [12] modified the MIP model of Amaral [11] by reformulating some constraints, obtaining a tighter model. Amaral [13] further established an improved MIP model by considering some valid inequalities and a symmetry breaking constraint. The DRLP with multiobjective optimization has also attracted the attention of researchers: Murray et al. [14] proposed a mixed integer linear programming formulation for an extended DRLP with non-zero aisle widths. The extended DRLP sought to minimize material handling costs and total layout area. They eliminated layout “mirroring” by new constraints and coupled a tabu search heuristic with CPLEX to solve the extended DRLP. Later on, another solution method, a multi-objective tabu search with a linear programming, was proposed by Zuo et al. [4] to solve the extended DRLP. Tang et al. [15] studied a robust DRLP where varying material flows in different periods were considered. A multiobjective evolutionary algorithm based on decomposition (MOEA/D) was proposed for solving the robust DRLP. Considering the dynamic environment, Wang et al. [16] proposed a dynamic DRLP, which is similar to the robust DRLP. They combined simulated annealing and mathematical programming to resolve the dynamic DRLP. Due to the efforts of these researchers, good progress has been made in recent years on solving the DRLP. However, it is still a challenging task to obtain optimal solutions within a reasonable time. The DRLP is related to two other facility layout problems: the corridor allocation problem (CAP) and the parallel row ordering problem (PROP). Similar to the DRLP, the CAP arranges facilities in two parallel rows along a central corridor. The difference between the CAP and the DRLP is that a CAP layout should respect two main conditions: the arrangements in both rows should start from a common point and no clearance is allowed between two adjacent facilities [17]. While these two conditions are not enforced in the DRLP. In the PROP, facilities are also arranged along two parallel rows. Different from the DRLP, the PROP has the following assumptions [18]: facilities are restricted to given rows; the arrangements in both rows start from a common point; no space is allowed between two adjacent facilities; the distance between the two parallel rows is not zero. Due to the assumption that no space is allowed between two adjacent facilities, the CAP and the PROP are combinatorial problems, where only the sequences of machines in double rows need to be determined. While the DRLP has to consider the sequence of facilities in each row and the clearance for each pair of adjacent facilities. The former one is a combinatorial optimization problem, while the latter one is a continuous optimization problem. The complexity of this problem increases because of the incorporation of both combinatorial and continuous aspects, which is one of the novel and thriving researches in optimization theory. Therefore, the DRLP is a more complex problem and has significant importance in academic fields. From our literature review, the following two observations can be made. First, it is common to handle the problem with the incorporation of combinatorial and continuous aspects by decomposition strategies. However, in the decomposition strategies of [8] and [9], the combinatorial subproblem is over simplified just considering that the arrangement of facilities on each row starts at the same abscissa. The optimal starting abscissa of the arrangement of facilities on each row is considered in the continuous subproblem. The over simplified combinatorial subproblem may provide inferior solutions to the continuous subproblem, hindering the search to converge towards the global optimum. Second, CPLEX is able to obtain the optimal solutions. However, it consumes too much time when dealing with the problem with large scale. We aim to propose a modified decomposition, by considering the optimal starting abscissa of the arrangement of facilities on each row in the combinatorial subproblem rather than in the continuous subproblem. Instead of CPLEX, a population-based algorithm is going to be proposed to solve the continuous subproblem. When the DRLP is decomposed into combinatorial and continuous subproblems, different algorithms are separately applied to handle them. Local search is an efficient heuristic approach, which iteratively improves a solution by searching its neighborhood for a better one. It has many variants: basic local search with hill climbing, tabu search, variable neighborhood search, simulated annealing, and so on. Due to its powerful capability of exploitation, local search has been successfully applied for solving different combinatorial optimization problems: traveling thief problem [19], multiple knapsack problem [20], course timetabling problem [21], quantity discounts supply chain planning problem [22], aircraft landing problem [23], max-cut problem [24] and so on. Specifically, it has been successfully implemented for solving facility layout problems such as the single row facility layout problem [25] and the corridor allocation problem [17]. Basic local search with hill climbing is the most simple one. Therefore, we focus on applying a basic local search with hill climbing to solve the combinatorial subproblem of the DRLP. Particle swarm optimization (PSO) is a swarm intelligence algorithm originally developed by

J. Guan, G. Lin and H.-B. Feng et al. / Applied Mathematical Modelling 77 (2020) 963–979

965

Kennedy and Eberhart [26]. It simulates the movement of birds to find food. Because of its fast convergence speed and strong global search ability [27], PSO has been widely used in many fields, especially for continuous problems: global optimization of multimodal functions [28], the optimal power flow problem [29], coverage control of sensor networks [30], structure learning of bayesian networks [31] and so on. What is more, it also owns a merit of easy implementation, which attracts us to apply PSO to solve the continuous subproblem of the DRLP. To the best of our knowledge, there is almost no research work about an algorithm using PSO to solve the continuous subproblem of the DRLP, which is commonly solved by the solver CPLEX in the literature. The main contributions of this paper can be summarized as follows: (1) For the double row layout problem, we introduce a more understandable model than the MIP models in the literature. (2) We have proved a property of the objective function of the double row layout problem: the optimal starting abscissa of the arrangement of facilities in each row can be determined by the sequences in these rows. Therefore, we decompose the double row layout problem into combinatorial and continuous subproblems, considering the optimal starting abscissas in the combinatorial subproblem. It is different from the decomposition strategies in the literature [8,9], where the optimal starting abscissa of the arrangement of facilities in each row is considered in the continuous subproblem. Our combinatorial subproblem is more similar to the double row layout problem and provides superior solutions to the continuous subproblem, helping the search to converge towards the global optimum. (3) We have developed a decomposition-based algorithm (DBA) to solve the double row layout problem. The DBA is able to improve best known results on seven instances. Moreover, to the best of our knowledge, DBA is the first method that is tested on the instances with up to 80 facilities. The remainder of this paper is organized as follows: Section 2 gives a detailed description of the double row layout problem (DRLP). By analyzing the characteristic of the DRLP, we devise a modified decomposition. Section 3 proposes a decomposition-based algorithm (DBA) to solve the DRLP. Section 4 shows the simulation results to demonstrate the efficiency of the proposed DBA. Section 5 outlines our conclusions and directions for future research.

2. Problem statement 2.1. Basic model In the literature, it is common to describe the double row layout problem by a MIP model. In this paper, we give a more understandable description. Given a set B = {b1 , b2 , . . . , bn } of n ≥ 2 facilities, they are in rectangular shape. The length of facility p ∈ B is lp . An n × n matrix C = [c pq ] is given, where cpq is the material flow between facilities p ∈ B and q ∈ B. The value of cpq is non-negative. Another n × n matrix A = [a pq ] is given, where apq is the minimum clearance between facilities p and q. When two facilities p and q are arranged adjacently, the clearance between them should be no less than their minimum clearance. It means that there is an adjustable clearance δ pq ≥ 0 between them. The DRLP consists of arranging the facilities in a layout of double parallel rows, so that the total cost to transport materials among these facilities is minimized. The distance between these two parallel rows is negligible. In a double row layout, m facilities are arranged on the upper row in the sequence π 1 starting at η1 abscissa. The other (n − m ) facilities are arranged on the lower row in the sequence π 2 starting at η2 abscissa. Then, a basic model of the DRLP is given by

min

f =

|π |π1 | 1 |−1  i=1

+

|π |π2 | 2 |−1  i=1

+

cπ1 (i )π1 ( j ) |xπ1 (i ) − xπ1 ( j ) |

j=i+1

cπ2 (i )π2 ( j ) |xπ2 (i ) − xπ2 ( j ) |

j=i+1

|π1 |  |π2 | 

cπ1 (i )π2 ( j ) |xπ1 (i ) − xπ2 ( j ) |,

(1)

i=1 j=1

where:

xπr (i ) = ηr +

i−1  k=1

(lπr (k) + aπr (k)πr (k+1) + δπr (k)πr (k+1) ) +

lπr (i ) , 1 ≤ i ≤ |πr |; r ∈ {1, 2}. 2

(2)

The objective function value f evaluates the total cost to transport materials among facilities. Variable xπr (i ) is the abscissa of the center of the facility at position i on row r, where r = 1 is the upper row and r = 2 is the lower row. The cost includes three parts: the transportation cost among the facilities in the upper row, the transportation cost among the facilities in the lower row and the transportation cost among the facilities in the opposite row. For the objective function value, the sequences (π 1 , π 2 ), their starting abscissas (η1 , η2 ) and the adjustable clearances (δ pq ) are the decision variables. The sequences are the combinatorial aspects, while the starting abscissas and adjustable clearances are the continuous aspects.

966

J. Guan, G. Lin and H.-B. Feng et al. / Applied Mathematical Modelling 77 (2020) 963–979

2.2. Decomposition The combination of combinatorial and continuous aspects makes the DRLP much harder. Murray et al. [8] and Zhang and Cheng [9] decomposed the DRLP into two subproblems. In the first subproblem, the adjustable clearances between adjacent facilities are temporarily ignored and the left most facilities are placed along a virtual wall. The arrangement of facilities on each row can be deemed to start from zero abscissa. In other words, the first subproblem is a special case of DRLP under the condition that η1 = 0, η2 = 0 and δ pq = 0, where p, q ∈ B. We name their first subproblem as DRLP-ηδ . The task of the first subproblem is to determine optimal sequences of facilities. In the second subproblem, the sequences of facilities are fixed. The task of the second subproblem is to find optimal adjustable clearances and starting abscissas (locations of the left most facilities). We name their second subproblem as DRLP+ηδ . 2.3. Modified decomposition The decomposition makes it easier to solve the DRLP by dealing with the combinatorial and continuous aspects separately. However, it has an unsatisfactory performance on the solution quality because DRLP-ηδ is over simplified. The starting abscissa of the arrangement of facilities on each row has a significant influence on the total cost and therefore the solution quality. These motivate us to propose a modified decomposition. We try to determine optimal starting abscissas in the first subproblem rather than in the second subproblem. We name our first subproblem as DRLP-δ . Then, the task of the second subproblem is to find optimal adjustable clearances, named as DRLP+δ . In the first subproblem, the sequences with optimal starting abscissas are evaluated. Subsequently, an important issue to be addressed is how to determine the optimal starting abscissas when the sequences of double rows are given. According to the basic model of DRLP, there are two cases to consider for the DRLP-δ : a material flow occurs between two facilities in the same row and a material flow occurs between two facilities in the opposite row. When a material flow occurs between two facilities in the same row, the distance between them is not related to the starting abscissa. The value of the distance depends on the given sequence π r . When a material flow occurs between two facilities in the opposite row, the distance between them is dπ1 (i )π2 ( j ) = |xπ2 ( j ) − xπ1 (i ) |.

|xπ2 ( j ) − xπ1 (i) |       j−1 i−1   lπ2 ( j ) lπ1 (i )   =  η2 + (lπ2 (k) + aπ2 (k)π2 (k+1) ) + − η1 + (lπ1 (k) + aπ1 (k)π1 (k+1) ) +  2 2   k=1 k=1        j−1 i−1   lπ2 ( j ) lπ1 (i )   =  ( η2 − η1 ) − (lπ1 (k) + aπ1 (k)π1 (k+1) ) + − (lπ2 (k) + aπ2 (k)π2 (k+1) ) +  2 2   k=1 k=1 =

|η − (xπ1 (i) − xπ2 ( j ) )|,

(3)

where i−1 

xπ1 (i ) =

(lπ1 (k) + aπ1 (k)π1 (k+1) ) +

k=1

xπ2 ( j ) =

j−1 

(lπ2 (k) + aπ2 (k)π2 (k+1) ) +

k=1

lπ1 (i ) lπ (i−1) lπ ( i ) = xπ1 (i−1) + 1 + aπ1 (i−1)π1 (i ) + 1 , 2 2 2

(4)

lπ2 ( j ) lπ ( j−1) lπ ( j ) = xπ2 ( j−1) + 2 + aπ2 ( j−1)π2 ( j ) + 2 , 2 2 2

(5)

and η = η2 − η1 . Variable η is the difference between two starting abscissas. The values of xπ (i ) and xπ ( j ) depend on the 1 2 sequences π 1 and π 2 . Given sequences π 1 and π 2 , the cost of material handling between two facilities in the same row is determined, while the cost of material handling between two facilities in the opposite row is not determined because of the variable η. It is a task to find the optimal value of η to minimize the cost of material handing between two facilities in the opposite row. Algorithm 1 provides the pseudo-code of the calculation of the optimal difference between two starting abscissas. First, the cost of material handling between two facilities in the opposite row fopp is written as an expansion formula:

fopp =

|π1 |  |π2 |  i=1 j=1

cπ1 (i )π2 ( j ) |η − (xπ1 (i ) − xπ2 ( j ) )|

= cπ1 (1)π2 (1) |η − (xπ1 (1) − xπ2 (1) )| + cπ1 (1)π2 (2) |η − (xπ1 (1) − xπ2 (2) )| + · · · + cπ1 (i )π2 ( j ) |η − (xπ1 (i ) − xπ2 ( j ) )|

+ · · · + cπ1 (|π1 | )π2 (|π2 | ) |η − (xπ1 (|π1 | ) − xπ2 (|π2 | ) )|.

(6)

J. Guan, G. Lin and H.-B. Feng et al. / Applied Mathematical Modelling 77 (2020) 963–979

967

Algorithm 1 The calculation of the optimal difference between two starting abscissas. Input: A summation of fopp . Output: The optimal difference between two starting abscissas η.  1: sort the elements cπ (i )π ( j ) |η − (x π (i ) − xπ ( j ) )| of the summation of f opp in the ascending order of the value of their 1 2 1

(xπ

2

− xπ ( j ) ); 2 obtain the sequence of the elements as {c1 |η − e1 |, c2 |η − e2 |, . . . , ck |η − ek |, . . . , cs |η − es |}, where e1 < e2 < · · · < ek < · · · < es ; c ← 0 and K ← c1 + c2 + · · · + cs ; i ← 0; while (c ≤ K/2) do i + +; c ← c + ci ; end while η ← ei 1 (i )

2: 3: 4: 5: 6: 7: 8: 9:

Next, the elements of this summation cπ1 (i )π2 ( j ) |η − (xπ (i ) − xπ ( j ) )| are sorted in the ascending order of the value of their 1 2 (xπ (i) − xπ ( j ) ) and fopp is rewritten as: 1

2

fopp = c1 |η − e1 | + c2 |η − e2 | + · · · + ck |η − ek | + · · · + cs |η − es |,

(7)

e1 = min{(xπ (1) − xπ (1) ), . . . , (xπ (i ) − xπ ( j ) ), . . . , (xπ (|π | ) − xπ (|π | ) )}, 1 2 1 2 1 1 2 2

(8)

es = max{(xπ (1) − xπ (1) ), . . . , (xπ (i ) − xπ ( j ) ), . . . , (xπ (|π | ) − xπ (|π | ) )}, 1 2 1 2 1 1 2 2

(9)

where

and e1 < e2 <  < es . It is notable that the elements are removed if their ck = 0 and the elements are merged if their ek are the same. Finally, the variable η with the minimum value of fopp is determined by Property 2.1. Therefore, to minimize the objective value of DRLP-δ , only the sequences π 1 and π 2 need to be considered. According to Eqs. (4) and (5), it needs O(n) to calculate all xπ (i ) and xπ ( j ) . There are s = |π1 | × |π2 | elements in the 1

2

summation of fopp . In general, |π 1 | and |π 2 | are about n/2, then s is about n2 /4. The calculation of all differences (xπ (i ) − 1 xπ ( j ) ) spends O(s). Quick Sort algorithm is used to sort the elements of the summation. Its time complexity is O(s log s) in 2

average and O(s2 ) in the worst case. The calculation of K spends O(s). Then, it requires at most O(s) to find the optimal value for η. Therefore, the complexity of the calculation of η is O(n + 3s + s log s) = O(n + 3 · n2 /4 + n2 /4·log(n2 /4)) = O(n2 logn). Property 2.1. Given a function g(z ) = k1 |z − d1 | + k2 |z − d2 | + · · · + ki |z − di | + · · · + ks |z − ds |, where ki ∈ R+ (R+ is the set of nonnegative real numbers), i ∈ {1, 2, . . . , s}, s is the number of elements to be added and d1 < d2 <  < ds . Let K = k1 + k2 + · · · + ks and K (i ) = k1 + k2 + · · · + ki . When i = θ makes K(θ ) ≤ K/2 and K (θ + 1 ) > K/2, g(z) has a minimum value at z = dθ +1 . Proof. Let z ∈ (di , di+1 ). One can have

g( z ) = k 1 ( z − d 1 ) + k 2 ( z − d 2 ) + · · · + k i ( z − d i ) − ki+1 (z − di+1 ) − ki+2 (z − di+2 ) − · · · − ks (z − ds ) = (k1 + k2 + · · · + ki−1 + ki − ki+1 − · · · − ks )z + (−k1 d1 − k2 d2 − · · · − ki−1 di−1 − ki di + ki+1 di+1 + · · · + ks ds ) = (k1 + k2 + · · · + ki−1 + ki − ki+1 − · · · − ks )z + Q = [2(k1 + k2 + · · · + ki ) − (k1 + k2 + · · · + ki + · · · + ks )]z + Q = (2K (i ) − K )z + Q,

(10)

where Q = −k1 d1 − k2 d2 − · · · − ki−1 di−1 − ki di + ki+1 di+1 + · · · + ks ds , K (i ) = k1 + k2 + · · · + ki and K = k1 + k2 + · · · + ks is a constant. Since ki ∈ R+ , K(i) is a monotonic increasing function, with 0 < K(1) ≤ K(i) ≤ K. Therefore, there is exactly one value for i, which makes K(i) ≤ K/2 and K (i + 1 ) > K/2. Assume the value is θ . According the relationship between i and θ , we distinguish three cases to discuss g(z). Case 1: i < θ . In this case, K(i) < K(θ ) ≤ K/2, i.e. 2K (i ) − K < 0, which implies that g(z) is monotonic decreasing for each zi ∈ (di , di+1 ). Hence, g(d1 ) > g(z1 ) > g(d2 ), g(d2 ) > g(z2 ) > g(d3 ), . . . , g(dθ −1 ) > g(zθ −1 ) > g(dθ ). Therefore, g(d1 ) > g(z1 ) > g(d2 ) > g(z2 ) > g(d3 ) > · · · > g(di ) > g(zi ) > g(di+1 ) > · · · > g(dθ −1 ) > g(zθ −1 ) > g(dθ ), i.e. z ∈ [d1 , dθ ), we have g(z) > g(dθ ). Case 2: i > θ . In this case, K (i ) ≥ K (θ + 1 ) > K/2, i.e. 2K (i ) − K > 0, which implies that g(z) is monotonic increasing for each zi ∈ (di , di+1 ). Hence, g(dθ +1 ) < g(zθ +1 ) < g(dθ +2 ), g(dθ +2 ) < g(zθ +2 ) < g(dθ +3 ), . . . , g(ds−1 ) < g(zs−1 ) < g(ds ). Therefore, g(dθ +1 ) < g(zθ +1 ) < g(dθ +2 ) < g(zθ +2 ) < g(dθ +3 ) < · · · < g(di ) < g(zi ) < g(di+1 ) < · · · < g(ds−1 ) < g(zs−1 ) < g(ds ), i.e. z ∈ (dθ +1 , ds ], we have g(z ) > g(dθ +1 ).

968

J. Guan, G. Lin and H.-B. Feng et al. / Applied Mathematical Modelling 77 (2020) 963–979

Case 3: i = θ . In this case, K (i ) = K (θ ) = K/2 and K (i ) = K (θ ) < K/2 should be considered, respectively. If K (i ) = K (θ ) = K/2, i.e. 2Ki − K = 0, we have g(dθ ) = g(z ) = g(dθ +1 ) for z ∈ (dθ , dθ +1 ). If K (i ) = K (θ ) < K/2, i.e. 2K (i ) − K < 0, g(z) is monotonic decreasing for z ∈ (dθ , dθ +1 ). Hence, g(dθ ) > g(z ) > g(dθ +1 ). These three cases discussed above are just considering z ∈ [d1 , ds ]. In fact, there are the other two cases which have to be discussed. Case 4: z ∈ (−∞, d1 ). In this case, g(z ) = −Kz + Q is monotonic decreasing. Hence, g(z) > g(d1 ). Case 5: z ∈ (ds , ∞). In this case, g(z ) = Kz + Q is monotonic increasing. Hence, g(z) > g(ds ). Based on the observations of Cases 1 and 4, one can have g(z) > g(dθ ) when z ∈ (−∞, dθ ). Based on the observations of Cases 2 and 5, one can have g(z ) > g(dθ +1 ) when z ∈ (dθ +1 , ∞ ). Based on the observation of Case 3, one can have g(dθ ) ≥ g(z ) ≥ g(dθ +1 ) when z ∈ [dθ , dθ +1 ]. In conclusion, g(z ) ≥ g(dθ +1 ), z ∈ R. Therefore, when K(θ ) ≤ K/2 and K (θ + 1 ) > K/2, g(z) has a minimum value at z = dθ +1 . 3. A decomposition-based algorithm The characteristic of DRLP is the incorporation of both combinatorial and continuous aspects. Due to this characteristic, traditional methods are hard to obtain optimal solutions within a reasonable time. In this paper, we decompose the DRLP into two subproblems: a combinatorial optimization problem (DRLP-δ ) and a continuous optimization problem (DRLP+δ ). Based on the decomposition, first improvement based local search is firstly proposed to solve the DRLP-δ and then PSO is proposed to solve the DRLP+δ . 3.1. Solution representation To the DRLP-δ , a solution represents a layout of facilities, which is encoded as an integer vector π , an integer parameter m and a real number η. Vector π = (π (1 ), π (2 ), . . . , π (k ), . . . , π (n )) is a permutation of n facilities and parameter m is the cut point of the permutation π . It means that the first m facilities of π are assigned at the upper row π 1 , i.e. π1 = (π (1 ), π (2 ), . . . , π (m )) and the last (n − m ) facilities of π are assigned at the lower row π 2 , i.e. π2 = (π (m + 1 ), π (m + 2 ), . . . , π (n )). Real number η is the difference between the starting abscissas of the upper row and the lower row. In this paper, we always arrange facilities on the upper row starting at zero axis, i.e. η1 = 0. Then, the starting abscissa of the lower row is η2 = η. To the DRLP+δ , a solution represents the adjustable clearances between adjacent facilities, which is encoded as a vector δ = (δ (1 ), δ (2 ), . . . , δ (k ), . . . , δ (n )) ∈ (R+ )n . A nonnegative real number δ (k) represents the adjustable clearance between the facilities at positions (k − 1 ) and k for a given layout. It is noted that δ (1) is the adjustable clearance between the starting abscissa η1 and the first facility at the upper row, and δ (m + 1 ) is the adjustable clearance between the starting abscissa η2 and the first facility at the lower row. Then, a solution of DRLP is encoded as Z = (π , m, η, δ ), and its objective function value can be written as f(Z) or f(π , m, η, δ ). For example, there is an instance with ten facilities, whose lengths are 5, 3, 8, 4, 3, 5, 5, 7, 4, 6, respectively. The material flow matrix is given:

⎛

0 ⎜0 ⎜8 ⎜ ⎜4 ⎜ ⎜0 C=⎜ ⎜3 ⎜3 ⎜ ⎜2 ⎝4 0

0 0 3 1 0 7 1 5 0 4

8 3 0 0 6 0 2 0 9 1

4 1 0 0 4 6 0 1 0 3

0 0 6 4 0 0 0 5 1 11

3 7 0 6 0 0 2 7 1 1

3 1 2 0 0 2 0 0 0 7

2 5 0 1 5 7 0 0 0 3

4 0 9 0 1 1 0 0 0 0

⎞

0 4 ⎟ 1 ⎟ ⎟ 3 ⎟ ⎟ 11⎟ ⎟. 1 ⎟ 7 ⎟ ⎟ 3 ⎟ ⎠ 0 0

Each pair of facilities requires 0.5 minimum clearance. A solution of the instance consists of cut point m = 5, the difference between two starting abscissas η = 1.5, permutation π = (8, 4, 5, 7, 3, |6, 2, 10, 1, 9), and adjustable clearance vector δ = (0.0, 0.0, 1.0, 0.0, 0.0, | 0.0, 0.0, 0.0, 0.5, 0.0). The solution represents a layout of the instance shown in Fig. 1. Five facilities are arranged on the upper row in the sequence π1 = (8, 4, 5, 7, 3 ) starting at zero abscissa. The other five facilities are arranged on the lower row in the sequence π2 = (6, 2, 10, 1, 9 ) starting at 1.5 abscissa. The adjustable clearance between facilities at positions 2 and 3 on the upper row is 1.0 and the adjustable clearance between facilities at positions 3 and 4 on the lower row is 0.5. While, the adjustable clearances between other adjacent facilities are zero. In addition, the solution of DRLP-δ can be encoded as Z = (π , m, η, δ0 ), where δ0 = (0, 0, . . . , 0 ). Once we obtain the solution code of a double row layout, the abscissa of the center of each facility is calculated by Eq. (2). For each pair of facilities, the transportation cost between two facilities is calculated by multiplying the distance between them with the material flow between them. For example, the abscissa of the center of the facility at position 3 on the upper row is xπ1 (3 ) = x5 = η1 + lπ1 (1 ) + aπ1 (1 )π1 (2 ) + δπ1 (1 )π1 (2 ) + lπ1 (2 ) + aπ1 (2 )π1 (3 ) + δπ1 (2 )π1 (3 ) + lπ1 (3 ) /2 = η1 + l8 + a8,4 +

J. Guan, G. Lin and H.-B. Feng et al. / Applied Mathematical Modelling 77 (2020) 963–979

969

Fig. 1. A solution of DRLP with ten facilities.

δ (2 ) + l4 + a4,5 + δ (3 ) + l5 /2 = 0 + 7 + 0.5 + 0 + 4 + 0.5 + 1 + 3/2 = 14.5. The abscissa of the center of the facility at position 3 on the lower row is xπ2 (3 ) = x10 = η2 + lπ2 (1 ) + aπ2 (1 )π2 (2 ) + δπ2 (1 )π2 (2 ) + lπ2 (2 ) + aπ2 (2 )π2 (3 ) + δπ2 (2 )π2 (3 ) + lπ2 (3 ) /2 = η2 + l6 + a6,2 + δ (m + 2 ) + l2 + a2,10 + δ (m + 3 ) + l10 /2 = 1.5 + 5 + 0.5 + 0 + 3 + 0.5 + 0 + 6/2 = 13.5. The transportation cost between them is cπ1 (3 )π2 (4 ) |xπ1 (3 ) − xπ2 (4 ) | = c5,10 |x5 − x10 | = 11 · |14.5–13.5| = 11. By summing the transportation costs of all pairs, we obtain a total value of 719.

3.2. Cut point selection The cut point selection is a process selecting a cut position in the permutation π , constructing an assignment of the double row layout. The facilities of π before (included) the cut position are assigned to the upper row, and the rest facilities are assigned to the lower row. It has a great impact on the quality of a solution. In the corridor allocation problem, Ghosh and Kothari [32] addressed that as a candidate optimal solution, the sum of the lengths of facilities at the upper row should be close to that at the lower row. Thus, a cut point selection (CPS) similar to their OBTAIN-CAP-SOLUTION algorithm is applied here. The clearances between adjacent facilities are considered in the CPS. Another difference is how to assign the next facility when the sum of lengths of facilities at the upper row is equal to that at the lower row. In the OBTAIN-CAPSOLUTION, the next facility is assigned to the lower row. In the CPS, the next facility is randomly assigned to the upper row or the lower row with equal probability. Algorithm 2 Cut point selection (CPS). Input: A permutation π . Output: The cut point m. 1: m ← 1; k ← n; L1 ← lπ (1 ) ;L2 ← lπ (n ) ; 2: for i from 1 to (n − 2 ) do if (L1 < L2 ) then 3: m + +; L1 ← L1 + lπ (m ) + aπ (m−1 )π (m ) ; 4: else if (L1 > L2 ) then 5: k − −; L2 ← L2 + lπ (k ) + aπ (k )π (k+1 ) ; 6: else 7: 8: r ← generate a uniform random number in (0, 1); if (r < 0.5) then 9: m + +; L1 ← L1 + lπ (m ) + aπ (m−1 )π (m ) ; 10: else 11: k − −; L2 ← L2 + lπ (k ) + aπ (k )π (k+1 ) ; 12: end if 13: end if 14: 15: end for

The pseudo-code of CPS is shown in Algorithm 2. At the beginning, the first facility in the permutation is assigned to the upper row and the the last facility is assigned to the lower row. Variable m increasing from 1 is used to record the position in the permutation, where the facility has been assigned to the upper row. While variable k decreasing from n is used to record the position in the permutation, where the facility has been assigned to the lower row. The sums of lengths of the upper row and the lower row are stored in L1 and L2 , respectively. They contain the lengths of facilities and the clearances between adjacent facilities. If L2 is larger than L1 , the following operator is done: m is increased by one and the facility at m position is assigned to the upper row. If L2 is smaller than L1 , the following operator is done: k is decreased by one and the facility at k position is assigned to the lower row. If L1 equals to L2 , one of above two operators is done at random. The process is repeated until all facilities in the permutation have been assigned. Finally, the cut point m is obtained for a given permutation π .

970

J. Guan, G. Lin and H.-B. Feng et al. / Applied Mathematical Modelling 77 (2020) 963–979

3.3. Local search for the DRLP-δ The purpose of the local search procedure is to seek a local optimum in the search space surrounding a given solution. The main idea of our local search is to look for a local optimum by using a first improvement strategy in the exchange neighborhood. It is a simple and powerful algorithm for combinatorial optimizations. The local search for solving DRLP-δ is shown in Algorithm 3. An exchange of the incumbent permutation π ∗ is performed by randomly selecting a position pair from the set of unchosen position pairs U, generating a new permutation π . The selected position pair is subsequently removed from U. Once a new permutation π is generated, the optimal difference between two starting abscissas η is calculated by Property 2.1. Thus, a new solution of DRLP-δ is obtained. The objective function value of the new solution is compared with the incumbent one. As soon as a new solution is better than the incumbent solution, it is accepted as the new incumbent solution and the search restarts from there. The objective function value of the incumbent solution is stored in flbest . The local search terminates when all exchanges of the incumbent permutation have been performed without improvement. Throughout the local search process, the cut point m remains unchanged, and the adjustable clearances between adjacent facilities δ are ignored. Algorithm 3 Local search. Input: An initial solution Z = (π , m, η, δ0 ). Output: The improved solution Z ∗ = (π ∗ , m, η∗ , δ0 ). 1: π ∗ ← π ,η ∗ ← η , f lbest ← f (Z ); 2: label_Restart: 3: U (unchosen position pairs set)← { (1, 2 ), (1, 3 ), . . . , (i, j ), . . . , (n − 1, n )}, i < j; 4: while (U = ∅) do randomly choose a position pair (such as i and j) from U; 5: 6: U ← U − {(i, j )}; π ← exchange these two facilities at the ith position and the jth position in π ∗ ; 7: η ← obtain the optimal difference between two starting abscissas by Property 2.1; 8: Z ← (π , m, η, δ0 ); 9: 10: if ( f (Z ) < flbest ) then π ∗ ← π ; η∗ ← η; flbest ← f (Z ); 11: goto label_Restart; 12: 13: end if 14: end while 15: Z ∗ ← (π ∗ , m, η ∗ , δ0 ); For a single execution of the “while” loop, there are at most n(n − 1 )/2 candidate solutions to be evaluated in the exchange neighborhood. The complexity of an evaluation mainly depends on the calculation of η (line 8) and f(Z) (line 10). In Section 2.3, the complexity of the calculations of η has been analyzed and it is O(n2 logn). According to Eq. (2), it needs O(n) to calculate the abscissas of all facilities. Using Eq. (1), the calculation of the objective function value spends O(|π1 |2 + |π2 |2 + |π1 | × |π2 | ) = O(n2 /4 + n2 /4 + n2 /4 ) = O(n2 ). In general, |π 1 | and |π 2 | are about n/2. Therefore, a single execution of the “while” loop spends O((n2 logn + n + n2 ) · n(n − 1 )/2 ) = O(n4 logn). When the “while” loop is executed β times, the complexity of local search is O(β n4 logn). 3.4. PSO for the DRLP+δ The particle swarm optimization (PSO) is inspired by fish and birds’ behavior. PSO consists of a group of particles, which have positions and velocities. A solution of an n-dimensional optimization problem is considered as the position vector ui = (ui (1 ), ui (2 ), . . . , ui (k ), . . . , ui (n )) of a particle i. Another vector, namely, velocity vector vi = (vi (1 ), vi (2 ), . . . , vi (k ), . . . , vi (n )) aids particle i to move towards the best position. The velocity is affected by two attractors: the historically best position of particle i (called personal best position pBesti ) and the historically best position of the entire swarm (called global best position gBest). At (t + 1 )th iteration, the velocity and position of each particle i are updated by the following equations:

vt+1 (k ) = wt · vti (k ) + ε1 · σ1 · ( pBesti (k ) − uti (k )) + ε2 · σ2 · (gBest (k ) − uti (k )), i

(11)

ut+1 (k ) = uti (k ) + vt+1 ( k ). i i

(12)

The parameters σ 1 and σ 2 are random numbers in the range of (0, 1). The cognitive parameter ε 1 and social parameter ε 2 are to weigh the relative importance of pBesti and gBest, respectively. The inertia parameter w is to weigh the impact of the previous velocity. Tsoulos and Stavrakoudis [33] suggested that the values for ε 1 and ε 2 should be in [1, 2] and the value for wt should be in [0, 1]. According to the work done by Li et al. [34], we set them to 1.5, 1.5 and 0.7, respectively. To avoid the explosion, the velocity on each dimension is limited in [-Vmax , Vmax ]. We set Vmax to 0.5.

J. Guan, G. Lin and H.-B. Feng et al. / Applied Mathematical Modelling 77 (2020) 963–979

971

For a given layout of facilities, the decision of the adjustable clearances between adjacent facilities δ is a continuous optimization problem. When the adjustable clearance vectors are considered as position vectors of particles, DRLP+δ can be dealt with by PSO. The PSO algorithm for solving DRLP+δ is presented in Algorithm 4. The best solution of DRLP-δ is given as input for the PSO algorithm. The permutation π , the cut point m and the difference between two starting abscissas η remain unchanged throughout the search process of PSO. In the initialization, a swarm is constructed with Pmax particles, whose positions (ui ) and velocities (vi ) are initialized with uniformly distributed random numbers. The initial positions ui of particles are stored in their personal best positions pBesti , where the best one is stored in the global best position gBest. The lower objective function value, the better the position. Then, the positions and velocities of particles are updated iteratively to find better positions. In the iteration, the velocities and positions are updated by Eqs. (11) and (12), respectively. The personal best position pBesti is updated if a new position of particle i is better with lower objective function value. If the new position is also better than the global best position gBest, it is set to be gBest. The updating process is repeated until the number of iterations reaches a predefined integer number Imax . When the process stops, the best solution gBest is produced for the adjustable clearances between adjacent facilities δ ∗ . Algorithm 4 PSO Algorithm. Input: A solution Z = (π , m, η, δ0 ). Output: An improved solution Z ∗ = (π , m, η, δ ∗ ). 1: initialize a swarm with Pmax particles 2: ui (k ) ← U (0, 1 ) and vi (k ) ← U (−0.5, 0.5 ), for every dimension k of each particle i. 3: pBesti ← ui , for each particle i. 4: gBest ← arg min{ f (π , m, η , pBesti )| pBesti = pBest1 , pBest2 , . . . , pBestn } 5: iter ← 1 6: while (iter < Imax ) do for i from 1 to Pmax do 7: Update the velocity vi of particle i by using Eq. (11); 8: 9: Update the position ui of particle i by using Eq. (12); if ( f (π , m, η, ui ) < f (π , m, η, pBesti )) then 10: pBesti ← ui ; 11: if ( f (π , m, η, ui ) < f (π , m, η, gBest )) then 12: 13: gBest ← ui ; end if 14: 15: end if end for 16: 17: it er ← it er +1 18: end while 19: δ ∗ ← gBest The number of particles Pmax is an important factor for the PSO algorithm. A larger number of particles increases the diversity of the swarm and its exploration ability but it increases the computational efforts at the same time. Marini and Walczak [35] mentioned that PSO is not sensitive to increase the exploration ability when the number of particles is larger than 50 in most cases. We set the number of particles Pmax to 30, as done in [36] and [37]. 3.5. Mutation The purpose of the mutation operator is to help the search to escape from local optima. It tries to explore new search regions surrounding local optima. The main idea of the mutation operator is to change the positions of a few facilities in the permutation of the local optimum. An inversion operator is used as the mutation operator. It works by selecting a substring from a permutation and flipping it to form a new one. The substring with a few facilities is randomly selected. For example, consider a permutation (1, 2, 3, 4, 5, 6, 7, 8) and suppose that four successive facilities are randomly selected from the third position. The inversion of these selected facilities gives a new permutation (1, 2, 6, 5, 4, 3, 7, 8). The number of facilities to inverse is a uniform random number in [2, λn] for each mutation operator, where λ represents a percentage of the total number of facilities. 3.6. Overall procedure of the proposed algorithm To solve the DRLP, we propose a decomposition-based algorithm (DBA). The main idea of DBA is to decompose the DRLP into two subproblems: a combinatorial problem DRLP-δ and a continuous problem DRLP+δ . A first improvement based local search is proposed to solve the DRLP-δ and then PSO is applied to solve the DRLP+δ . The DBA balances well between the global exploration and local exploitation. The pseudo code of the proposed DBA is illustrated in Algorithm 5. At first, an initial solution Z is generated (lines 3–6). For the solution Z, the permutation π is randomly generated and its cut point m

972

J. Guan, G. Lin and H.-B. Feng et al. / Applied Mathematical Modelling 77 (2020) 963–979

Algorithm 5 Decomposition-based algorithm. Input: The predefined integer numbers Gmax and Tmax . Output: The global best solution Z ∗ and its value fgbest . 1: f gbest ← ∞. 2: for g from 1 to Gmax do π ← generate a random permutation; 3: m ← CPS(π ); 4: η ← obtain the optimal difference between two starting abscissas by Property 2.1; 5: Z ← (π , m, η, δ0 ); 6: for t from 1 to Tmax do 7: Z r ← LocalSearch(Z); 8: Z a ← P SO(Z r ); 9: if ( f (Z a ) < fgbest ) then 10: fgbest ← f (Z a ); Z ∗ ← Z a ; 11: end if 12: π ← mutation(); 13: m ← CPS(π ); 14: η ← obtain the optimal difference between two starting abscissas by Property 2.1; 15: Z ← (π , m, η, δ0 ); 16: end for 17: 18: end for

is determined by the CPS heuristic. The facilities before (including) the cut point are assigned to the upper row starting at zero abscissa. While the other facilities are assigned to the lower row starting at η abscissa, which is calculated by Property 2.1. Without considering the adjustable clearances, Z is a solution of DRLP-δ . For solving the DRLP-δ , a first improvement based local search is applied to look for a local optimum Zr in the search space surrounding the solution Z (line 8). Next, PSO is applied to solve the DRLP+δ . PSO looks for optimal adjustable clearances for the adjacent facilities of solution Zr and obtains a further improved solution Za (line 9). During the search, the global best solution Z∗ will be updated when an improved solution Za is better (lines 10–12). Afterward, a new solution Z near to the improved solution is generated by the mutation operator, which helps the search to escape from local optima, considering the quality criteria (lines 13–16). Then, local search and PSO continue to improve the new solution Z and try to find a new global best solution. The mutation operator is repeated within Tmax iterations. After that, the search restarts from a new completely random solution Z, which is far away from the improved solution. It helps the search to escape from local optima, considering the diversity criteria. The algorithm is stopped when the search has been restarted Gmax times.

4. Computational results and comparisons To evaluate the performance of the proposed algorithm, three parts of computational experiments have been conducted. First, we tune the parameters of the proposed algorithm. Second, we compare the results obtained by the proposed algorithm with those obtained by other methods in the literature. Third, we compare the proposed algorithm with its two variants to evaluate the most effective component of the proposed algorithm. The proposed DBA is coded in C and performed on a personal computer with Intel core i5-4460 3.20 GHz processor under Windows 7 professional system. To assess the reliability of the heuristic performance, all experiments are run 10 times with random seed. Currently, the state-of-the-art method for the DRLP is studied by Amaral [13]. For a fair comparison of DBA with Amaral’s method, two sets of benchmark instances in his work [13] are also used here: (a) four small instances (S9, S9H, S10, S11) introduced by Simmons [38]; (b) twenty random instances with 9 and 10 facilities and eighteen random instances with n ∈ {11, 12, 13} proposed by Amaral [13]. It should be noted that the instances of these two sets have implicit minimum clearances, which means that minimum clearances apq do not appear in these instances. The minimum clearances are deemed to be all equal and included in the lengths of the facilities. Since the first two instance sets are of small dimension and far from being challenging for our proposed method. The following two sets of larger instances are also used to test our proposed algorithm: eight instances ranging from 15 to 50 facilities created by Murray et al. [8] and four instances ranging from 60 to 80 facilities introduced by Anjos et al. [39]. The former instances have explicit minimum clearances and the latter instances have implicit minimum clearances. It has a risk of getting over fitted results for benchmark instances, if the same instances is used to tune the parameters of an algorithm [40]. In order to avoid this risk, we generated a set of five random instances to tune the parameters of the proposed algorithm. The method of generating a random instance was described in [5]. For five tuning instances with n = 10, 15, 20, 25 and 30 facilities, each factor was created by using a uniform distribution. The lengths of facilities lp were created by U(0, 20). The material flows cpq were created by U(0, 50). The minimum clearances apq were also randomly created

J. Guan, G. Lin and H.-B. Feng et al. / Applied Mathematical Modelling 77 (2020) 963–979

973

Table 1 The parameters and their levels. Parameters

Levels

Symbol

Description

Section

1

2

3

4

Imax Gmax Tmax

Maximum Maximum Maximum Maximum

3.4 3.6 3.6 3.5

log (n) 6log (n) 6log (n) 0.1

2log (n) 8log (n) 8log (n) 0.3

3log (n) 10log (n) 10log (n) 0.5

4log (n) 12log (n) 12log (n) 0.7

λ

iteration of PSO iteration of overall procedure number of mutation executions percentage of facilities to mutate

Table 2 Orthogonal array and results of the experiment. Exp no.

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16

Combinations of parameters

RPD

Imax

Gmax

Tmax

λ

Rt10

Rt15

Rt20

Rt25

Rt30

1 1 1 1 2 2 2 2 3 3 3 3 4 4 4 4

1 2 3 4 1 2 3 4 1 2 3 4 1 2 3 4

1 2 3 4 2 1 4 3 3 4 1 2 4 3 2 1

1 2 3 4 3 4 1 2 4 3 2 1 2 1 4 3

0.37 0.00 0.00 0.00 0.00 0.09 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00

0.72 0.06 0.00 0.00 0.11 0.47 0.00 0.03 0.08 0.03 0.35 0.03 0.05 0.00 0.00 0.43

0.68 0.03 0.04 0.00 0.13 0.50 0.00 0.00 0.17 0.00 0.25 0.00 0.06 0.00 0.00 0.15

0.54 0.23 0.03 0.00 0.30 0.40 0.00 0.04 0.17 0.09 0.26 0.07 0.19 0.11 0.11 0.19

1.00 0.40 0.15 0.10 0.63 0.82 0.09 0.12 0.50 0.21 0.51 0.25 0.49 0.26 0.26 0.33

Table 3 Response value. Level

Imax

Gmax

Tmax

λ

1 2 3 4 Delta Rank

0.218 0.187 0.149 0.132 0.086 3

0.310 0.186 0.103 0.087 0.223 2

0.403 0.131 0.085 0.066 0.338 1

0.206 0.154 0.141 0.184 0.065 4

by U(1, 2). These five tuning instances are named as Rt10, Rt15, Rt20, Rt25 and Rt30. These five instances have explicit minimum clearances. 4.1. Parameter tuning for DBA We first carried out Taguchi based experiments to determine the parameters of the proposed algorithm: the maximum iteration of PSO (Imax ), the maximum iteration of overall procedure (Gmax ), the maximum number of mutation executions (Tmax ), and the maximum percentage of the total number of facilities to mutate (λ). For each parameter, four levels are considered and their values are presented in Table 1. The values of Imax , Gmax and Tmax depend on the size of tested instance n. An orthogonal array L16 (44 ) is selected. Based on the orthogonal array, we tested the DBA with each combination of parameters on five tuning instances. The result of each test reports the relative percentage deviation (RPD) of the objective function value as shown in Table 2. The value of RPD for instance i is calculated by the following formula: nr 1  RP Di = nr j=1





fi j − LBi × 100 . LBi

(13)

where fij is the objective function value obtained in the jth replication for instance i, LBi is the minimum objective function value obtained for instance i, nr is the number of replications. For each instance, nr = 10 independent replications are carried out with a combination of parameters. With the results in Table 2, the response value of each parameter is figured out to analyze the significance rank of each parameter as shown in Table 3. From Table 3, one can observe that Tmax is the most significant parameter, while λ is the least significant one. Furthermore, the main effects for means of factor levels are plotted to illustrate the factor level trend of the DBA as shown in Fig. 2. The DBA exhibits better performance as Tmax

974

J. Guan, G. Lin and H.-B. Feng et al. / Applied Mathematical Modelling 77 (2020) 963–979

Fig. 2. Factor level trend of the DBA.

increases. Similar tendencies appear for Imax and Gmax . The DBA exhibits better performance, when the value of λ trends towards the medium level. A larger value of Gmax or Tmax gives the search more chances to explore new promising regions, increasing the exploration ability. A larger value of Imax helps PSO to converge to local optima, increasing the exploitation ability. A smaller value of λ makes new starting solutions so near to the local optimum that the search is easily trapped in local optimum. Instead, a larger value of λ makes new starting solutions so far away from the local optimum that the exploitation ability may be decreased. Fig. 2 suggests that the DBA will have a better performance when the parameters are set as: Imax = 4 log(n ) (Level 4), Gmax = 12 log(n ) (Level 4), Tmax = 12 log(n ) (Level 4) and λ = 0.5 (Level 3). Thus, these values are adopted in the following experiments. 4.2. Comparison with other methods Currently, the state-of-the-art method for the DRLP is proposed by Amaral, called M2 [13]. Thus, we compared DBA with M2 and another method of Amaral (M1 [11]). The author used the CPLEX 12.7.1.0 solver to deal with two different MIP models of the DRLP. M1 and M2 are exact methods. Table 4 reports the comparison of results on Simmons and Amaral instance sets. The columns “Inst”, “n” and “Best” stand for the instance name, the number of facilities and the best solution obtained by the corresponding method. The columns “t” and “STD” stand for the average execution time and standard deviation of cost over ten independent executions. An additional row is included at the bottom summarizing the average values for each method on all test instances. The best solutions of [11] and [13] are proved optimal solutions because the CPLEX solver obtained the solution with an optimality gap of zero. As seen in the table, the three methods are able to achieve the optimal solution for every test instance. The standard deviation of zero for each instance verifies the robustness of our proposed algorithm. Amaral [13] tested four Simmons instances and the first 20 Amaral instances on a PC with an Intel Core i3-50 05U 2.0 0 GHz CPU, which is around 0.4 times faster than ours according to the CPU benchmarks from http://www.cpubenchmark.net/. While Amaral [13] tested the last 18 Amaral instances on a PC with an Intel Core i7-3770 3.40 GHz CPU, which is around 1.4 times faster than ours. The execution time from [13] is multiplied by the corresponding factor. Considering the average execution time, the DBA runs much faster than M1 and M2, which is obviously displayed in Fig. 3. For M1 and M2, the execution time increases more obviously as the number of facilities increases. The average execution time of DBA over the instances with n = 13 is 2.06 seconds, while those of M1 and M2 are 5752.39 s and 11506.26 s, respectively. We also applied DBA to deal with the double row layout problem with asymmetric material flow and compared it with a heuristic method called LS-minFFasym in [8]. The LS-minFFasym approach offered significantly better performance than the other methods. The LS-minFFasym approach was tested with eight instances ranging from 15 to 50 facilities on an HP 8100 Elite desktop PC. The desktop PC has a quad-core Intel i7-860 processor, which is around 0.75 faster times than ours. We got the original data of [8] from the authors. The comparison of DBA versus LS-minFFasym is given in Table 5. The execution time of LS-minFFasym has been multiplied by 0.75. The LS-minFFasym approach was run with a time limit of 5 minutes, so our DBA is run with a time limit of 3.75 (5 × 0.75) minutes to have a fair comparison. The structure of Table 5 is similar to that of Table 4, except that column “AVG” is added. Column “AVG” reports the average cost over ten independent executions. For each instance, the LS-minFFasym approach is implemented for only one

J. Guan, G. Lin and H.-B. Feng et al. / Applied Mathematical Modelling 77 (2020) 963–979

975

Table 4 Comparison of DBA versus exact methods. Inst

n

S9 S9H S10 S11 Small.09–1 Small.09–2 Small.09–3 Small.09–4 Small.09–5 Small.09–6 Small.09–7 Small.09–8 Small.09–9 Small.09–10 Small.10–1 Small.10–2 Small.10–3 Small.10–4 Small.10–5 Small.10–6 Small.10–7 Small.10–8 Small.10–9 Small.10–10 11a 11b 11c 11d 11e 11f 12a 12b 12c 12d 12e 12f 13a 13b 13c 13d 13e 13f Avg.

9 9 10 11 9 9 9 9 9 9 9 9 9 9 10 10 10 10 10 10 10 10 10 10 11 11 11 11 11 11 12 12 12 12 12 12 13 13 13 13 13 13

M1 [11]

M2 [13]

DBA

Best

t

Best

t

Best

t

STD

1179.0 2293.0 1351.0 3424.5 2393.0 2266.0 2241.0 2231.5 3980.0 4108.0 4089.0 4154.0 935.0 1073.5 1385.0 1437.0 1452.5 1313.5 722.5 792.0 607.5 529.0 929.5 828.0 5559.0 3655.5 3832.5 906.5 578.0 825.5 1493.0 1606.5 2012.5 1107.0 1066.0 997.5 2456.5 2864.0 4136.0 6164.5 6502.5 7699.5 2361.37

3.18 41.89 13.08 130.14 34.77 31.31 41.77 50.02 47.85 59.80 65.91 62.49 18.56 25.62 13.96 19.82 16.35 14.06 8.58 15.09 8.05 4.54 13.53 20.37 176.68 203.91 214.86 141.58 24.14 28.03 217.98 196.52 278.43 77.88 129.43 234.77 6956.49 1334.91 796.43 1855.90 2291.49 21279.10 885.70

1179.0 2293.0 1351.0 3424.5 2393.0 2266.0 2241.0 2231.5 3980.0 4108.0 4089.0 4154.0 935.0 1073.5 1385.0 1437.0 1452.5 1313.5 722.5 792.0 607.5 529.0 929.5 828.0 5559.0 3655.5 3832.5 906.5 578.0 825.5 1493.0 1606.5 2012.5 1107.0 1066.0 997.5 2456.5 2864.0 4136.0 6164.5 6502.5 7699.5 2361.37

2.92 32.85 9.57 180.34 30.08 60.64 81.08 41.52 44.65 123.98 104.74 57.85 14.41 18.46 13.36 15.51 20.04 17.78 10.18 22.88 16.38 6.96 12.90 16.95 481.46 388.53 482.89 96.96 88.91 44.30 524.48 565.49 575.79 161.52 168.63 210.49 17473.64 15309.87 2030.69 4732.08 12802.31 16688.98 1756.74

1179.0 2293.0 1351.0 3424.5 2393.0 2266.0 2241.0 2231.5 3980.0 4108.0 4089.0 4154.0 935.0 1073.5 1385.0 1437.0 1452.5 1313.5 722.5 792.0 607.5 529.0 929.5 828.0 5559.0 3655.5 3832.5 906.5 578.0 825.5 1493.0 1606.5 2012.5 1107.0 1066.0 997.5 2456.5 2864.0 4136.0 6164.5 6502.5 7699.5 2361.37

0.54 0.51 0.75 1.10 0.50 0.49 0.51 0.50 0.50 0.50 0.52 0.51 0.53 0.50 0.77 0.78 0.75 0.76 0.74 0.73 0.74 0.70 0.75 0.74 1.11 1.12 1.09 0.97 0.96 0.98 1.55 1.46 1.50 1.37 1.46 1.34 1.97 1.97 1.95 2.12 2.12 2.23 1.02

0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0

run. Thus, the “STD” and ‘AVG’ of LS-minFFasym are not listed in Table 5. The results marked in bold indicate the best cost among the reference algorithms. Compared with LS-minFFasym, DBA is better on seven instances and equal on one instance in terms of the best cost. The best costs of the last seven Murray instances are reduced ranging from 1.19% to 25.11%. With respect to the average cost, DBA also has a better performance than LS-minFFasym. The results clearly show that DBA outperforms LS-minFFasym. 4.3. Effectiveness of the key feature of DBA The most key feature of DBA is to decompose the DRLP into two subproblems. The decomposition is different from those of [8] and [9]. In the first subproblem of the decomposition in [8] or [9], the left most facilities start at the same abscissa. The optimal starting abscissa of the arrangement of facilities on each row is considered in the second subproblem. While in our decomposition, the optimal starting abscissa of the arrangement of facilities on each row is considered in the first subproblem, where the left most facilities might start at the same abscissa or not. The optimal starting abscissas are determined according to Property 2.1. To analyze the efficiency of our decomposition, two variants of DBA were designed to be compared with DBA. The first one is DBA-ηδ , where our decomposition is replaced by the decomposition method of [9] keeping the other ingredients unchanged. The second one is DBA-PSO, where PSO algorithm is removed from DBA. In other words, its second subproblem is ignored based on our decomposition. To have a fair comparison, these three algorithms were run with the same time limits of (n/10)3 s on a certain instance, where n is the size of the corresponding instance. Experiments were carried out on all instances mentioned above.

976

J. Guan, G. Lin and H.-B. Feng et al. / Applied Mathematical Modelling 77 (2020) 963–979

Fig. 3. Time consumption comparison of DBA with other methods.

Table 5 Comparison of DBA versus a heuristic method. Inst

n

Murray155 Murray206 Murray218 Murray231 Murray243 Murray256 Murray309 Murray322 Avg.

15 20 25 30 35 40 45 50

LS-minFFasym [8]

DBA

Best

t

Best

t

STD

AVG

105578.4 320547.1 636848.9 930012.8 1306111.9 2468886.1 3509780.5 5039125.5 1789611.4

6.97 17.94 101.32 209.16 225.36 225.02 226.32 225.32 154.68

105578.4 311751.1 624653.3 918759.6 1257631.8 2118158.1 2628585.5 4018808.4 1497990.8

4.78 16.53 44.31 101.15 211.11 225.06 225.06 225.10 131.64

0.0 0.0 0.0 655.4 974.1 2096.5 3810.8 3018.0 1319.3

105578.4 311751.1 624653.3 919550.5 1259297.9 2120510.5 2634068.4 4026413.7 1500228.0

Table 6 summarizes the results of experiments, except those instances which are so easy that the three compared algorithms are all able to reach their optimal solutions steadily. Considering the best cost, we can observe that DBA-PSO performs significantly better than DBA-ηδ . DBA-PSO achieves the minimal best cost on 22 out of 30 instances. While DBAηδ achieves the minimal best cost only on 8 out of 30 instances. When PSO is combined into DBA-PSO, a hybrid algorithm DBA becomes further improved. DBA is able to achieve the minimal best cost on 26 out of 30 instances. In order to evaluate the performances of these three algorithms in statistical significance, we also applied a Friedman test to the data in Table 6. The results of Friedman test are listed in Table 7. The average rankings of these three algorithms are reported in the middle three columns and the p-value is reported in the last column. According to the average rankings of ‘Best’, these three algorithms are sorted into the following order: DBA, DBA-PSO and DBA-ηδ . So as ‘AVG’. Thus, DBA exhibits the best performance. With the p-value of 4.84 × 10−4 , there is enough statistical evidence to confirm the significant differences in the best cost between these three algorithms. In terms of average cost, there is also statistically significant difference among these three algorithms, because the p-value is 0.02 at the 0.05 level. With respect to the standard deviations, the results obtained by DBA and DBA-PSO are all zero on the instances with the size smaller than 20. These verify the robustness of DBA and DBA-PSO on small instances. It can be extracted that the algorithms based on our decomposition perform well in terms of the solution quality. It highlights the advantage of our decomposition, considering the optimal starting abscissa of the arrangement of facilities on each row in the first subproblem. Another interesting observation can be obtained from the excellent performance of DBAPSO: without considering the adjustable clearances between adjacent facilities, the solutions of our first subproblem match the best solutions of the DRLP on most test instances. This phenomenon can be explained by the following fact. For the single row facility layout problem (SRFLP), it is clear that an arrangement of facilities with the lowest total transportation cost must be in the condition that there is no adjustable clearances between adjacent facilities. Although the DRLP is different from SRFLP, the characteristic of SRFLP still exists in the upper row or the lower row. Only when the material flows between facilities in the opposite row are significantly larger than those in the same row, adjustable clearances between

J. Guan, G. Lin and H.-B. Feng et al. / Applied Mathematical Modelling 77 (2020) 963–979

977

Table 6 Comparison of DBA versus two variant algorithms. Inst

n

S10 S11 Small.10–3 Small.10–4 Small.10–5 Small.10–6 Small.10–7 11a 11e 12a 12b 12c 13a 13c Rt15 Rt20 Rt25 Rt30 Murray155 Murray206 Murray218 Murray231 Murray243 Murray256 Murray309 Murray322 AKV-60-01 AKV-70-01 AKV-75-01 AKV-80-01 Avg.

10 11 10 10 10 10 10 11 11 12 12 12 13 13 15 20 25 30 15 20 25 30 35 40 45 50 60 70 75 80

DBA-ηδ

DBA-PSO

DBA

Best

AVG

STD

t

Best

AVG

STD

t

Best

AVG

STD

t

1364.5 3427.5 1453.0 1316.5 728.5 793.0 609.5 5621.5 578.5 1524.0 1608.5 2012.5 2457.5 4137.0 43966.3 133281.0 249848.8 400926.4 105578.4 311987.3 624653.3 919239.6 1258926.5 2121176.1 2628451.3 4019903.2 739261.0 764820.0 1197267.5 1035020.5 552731.3

1364.5 3427.5 1453.0 1316.5 728.5 793.0 609.5 5621.5 578.5 1524.0 1608.5 2014.6 2457.5 4137.0 43980.6 133315.6 250022.0 401271.2 105591.0 311987.3 624816.3 919893.9 1260353.0 2122352.8 2631777.4 4024732.9 739384.1 765495.6 1198393.9 1035508.0 553217.0

0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 4.7 0.0 0.0 30.6 76.4 161.9 384.8 39.6 0.0 270.6 497.5 1100.1 1120.7 2594.8 3590.9 122.4 452.8 712.1 318.8 382.6

0.50 0.65 0.50 0.50 0.51 0.51 0.50 0.64 0.65 0.88 0.88 0.86 1.04 1.04 1.57 3.88 8.22 15.65 2.44 6.58 14.53 27.01 42.89 64.02 91.15 125.04 216.04 343.04 421.93 512.15 63.5

1351.0 3424.5 1452.5 1313.5 722.5 792.0 609.0 5559.0 578.0 1493.0 1606.5 2012.5 2456.5 4136.0 43966.3 133281.0 249848.8 400597.9 105578.4 311751.1 624653.3 919657.8 1259002.4 2119755.3 2630090.7 4020367.4 739246.5 764511.0 1197419.5 1034749.5 552732.8

1351.0 3424.5 1452.5 1313.5 722.5 792.0 609.0 5559.0 578.0 1493.0 1606.5 2012.5 2456.5 4136.0 43966.3 133281.0 249873.5 401310.4 105578.4 311751.1 625115.0 920208.0 1260590.4 2122849.2 2636889.5 4028416.8 739842.5 765475.0 1199071.6 1035724.4 553581.6

0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 78.0 379.5 0.0 0.0 542.2 472.4 1121.0 1923.3 4018.1 4440.0 503.2 718.7 1442.8 603.8 541.4

0.29 0.51 0.29 0.30 0.29 0.27 0.29 0.57 0.40 0.79 0.71 0.74 1.08 1.06 2.52 8.01 15.63 27.02 2.86 8.01 15.63 27.02 42.91 64.06 91.26 125.15 216.36 343.47 422.73 513.01 64.4

1351.0 3424.5 1452.5 1313.5 722.5 792.0 607.5 5559.0 578.0 1493.0 1606.5 2012.5 2456.5 4136.0 43966.3 133281.0 249848.8 400597.9 105578.4 311751.1 624653.3 918759.6 1259603.8 2120046.1 2628585.5 4019272.3 739233.0 764554.0 1197240.0 1034735.5 552640.4

1351.0 3424.5 1452.5 1313.5 722.5 792.0 607.5 5559.0 578.0 1493.0 1606.5 2012.5 2456.5 4136.0 43966.3 133281.0 249952.7 401264.6 105578.4 311751.1 624995.2 919890.4 1261059.5 2122665.5 2634717.2 4027539.6 739779.9 765678.2 1199444.0 1035365.6 553481.1

0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 170.2 407.4 0.0 0.0 583.8 666.0 1460.0 1734.6 3436.3 5638.5 579.6 801.0 2153.7 516.3 604.9

0.72 1.06 0.72 0.74 0.73 0.72 0.72 1.12 0.94 1.50 1.41 1.45 1.89 1.89 3.38 8.01 15.63 27.02 3.38 8.01 15.63 27.02 42.92 64.07 91.25 125.14 216.28 343.48 422.55 512.79 64.7

Table 7 Results of Friedmans test applied to the data in Table 6. Average ranking

Best AVG

p-value

DBA-ηδ

DBA-PSO

DBA

2.57 2.40

1.83 1.92

1.60 1.68

4.84 × 10−4 0.02

some adjacent facilities are needed to make several facilities in the opposite row closer. In this condition, the adjustable clearances between these adjacent facilities will be further dealt with by PSO. 4.4. Managerial implications One of the managerial implications of the DRLP model is that a manufacturing company may save the material handling costs when machines are arranged in proper positions. However, it is a challenging task to arrange machines efficiently. Our solution method can be applied to support managerial decisions to determine high quality machine layouts in acceptable time. This study also provides an insight into which factors related to layout affect the the material handling costs significantly. From our experiments, it is observed that the starting abscissas of double rows and the sequences of facilities affect the material handling costs significantly, while the adjustable clearances between adjacent facilities affect the material handling costs slightly. It indicates how to simplify the DRLP model in a context of production process where the speed of layout optimization is also an important objective to be considered. 5. Conclusions and future research We have developed a decomposition-based algorithm (DBA) to solve the double row layout problem (DRLP). The DRLP is decomposed into a combinatorial subproblem and a continuous subproblem, which are dealt with by first improvement based local search and PSO, respectively. Different from the decomposition strategies of [8] and [9], our decomposition considers the sequence of facilities on each row with its optimal starting abscissa in the combinatorial subproblem. Although the starting abscissas of double rows are continuous variables, their optimal difference can be determined by the given

978

J. Guan, G. Lin and H.-B. Feng et al. / Applied Mathematical Modelling 77 (2020) 963–979

sequences of double rows using a property of the objective function of the DRLP. The property has been proved in this paper. The experiments on 59 instances have been carried out to compare the DBA with state-of-the-art methods. The results have shown that the DBA is able to find the optimal solutions in much less time, compared with two exact methods M1 [11] and M2 [13]. In addition, the DBA is able to find better solutions than a heuristic method LS-minFFasym [8], within the same time limit on a certain instance. An extra experiment has been conducted to analyze the effectiveness of the key feature of the proposed DBA. The experiment highlights the advantage of our decomposition in terms of the solution quality. Our combinatorial subproblem considering the sequences and their optimal starting abscissas is more similar to DRLP than one only considering the sequences starting at the same abscissa. Considering the combination of the sequences and their optimal starting abscissas in the first subproblem is more helpful to find the optimal solutions of DRLP. Another interesting observation can be obtained from the experiment: for most DRLP instances, there are no adjustable clearances between adjacent facilities in their optimal solutions. This phenomenon indicates that most DRLP instances belong to the combinatorial subproblem of our decomposition, which can be solved directly by the first improvement based local search. We can conclude that the DBA has excellent performance in terms of solution quality. However, the problem is simplified from a realistic scenario without considering the width of the corridor. In future work, we will further study the problem considering the width of the corridor. For the single row facility layout problem, large-scale instances have recently been solved in [41] with the use of gain techniques for computing objective values. The gain techniques are able to accelerate the process of algorithm. However, there is no effective gain technique for the DRLP till now. Therefore, another interesting topic for future research would be the development of effective gain techniques for the DRLP. Acknowledgments This research was partially supported by the National Natural Science Foundation of China under Grants 61871204, 11301255, the Science and Technology Project of the Education Bureau of Fujian, China, under Grant JAT160396, the Natural Science Foundation of Fujian Province of China under Grant 2017J01767, 2016J01025, the Program for New Century Excellent Talents in Fujian Province University, and the Fund by Collaborative Innovation Center of IoT Industrialization and Intelligent Production, Minjiang University (IIC1703). We wish to thank Chase C. Murray for providing us the test instances with 15–50 facilities and the results of their algorithms. References [1] S.S. Heragu, A. Kusiak, Machine layout problem in flexible manufacturing systems, Oper. Res. 36 (1988) 258–268. [2] H. Ahonen, A.G. de Alvarenga, A.R.S. Amaral, Simulated annealing and tabu search approaches for the corridor allocation problem, Eur. J. Oper. Res. 232 (2014) 221–233. [3] P. Kouvelis, W.C. Chiang, A.S. Kiran, A survey of layout issues in flexible manufacturing systems, Omega 20 (1992) 375–390. [4] X. Zuo, C.C. Murray, A.E. Smith, Solving an extended double row layout problem using multiobjective tabu search and linear programming, IEEE Trans. Autom. Sci. Eng. 11 (2014) 1122–1132. [5] J. Chung, J. Tanchoco, The double row layout problem, Int. J. Prod. Res. 48 (2010) 709–727. [6] M. Mohammadi, K. Forghani, A novel approach for considering layout problem in cellular manufacturing systems with alternative processing routings and subcontracting approach, Appl. Math. Model. 38 (2014) 3624–3640. [7] Z. Zhang, C.C. Murray, A corrected formulation for the double row layout problem, Int. J. Prod. Res. 50 (2012) 4220–4223. [8] C.C. Murray, A.E. Smith, Z. Zhang, An efficient local search heuristic for the double row layout problem with asymmetric material flow, Int. J. Prod. Res. 51 (2013) 6129–6139. [9] Z.-q. Zhang, W.m. Cheng, Decomposition strategies and heuristic for double row layout problem, Comput. Integr. Manuf. Syst. 20 (2014) 559–568. [10] M.F. Anjos, A. Fischer, P. Hungerländer, Solution Approaches for the Double-Row Equidistant Facility Layout Problem, Springer International Publishing, 2016. [11] A.R.S. Amaral, Optimal solutions for the double row layout problem, Optim. Lett. 7 (2013) 407–413. [12] L.D. Secchin, A.R.S. Amaral, An improved mixed-integer programming model for the double row layout of facilities, Optim. Lett. 13 (2019) 193–199. [13] A.R.S. Amaral, A mixed-integer programming formulation for the double row layout of machines in manufacturing systems, Int. J. Prod. Res. 57 (2019) 34–47. [14] C.C. Murray, X. Zuo, A.E. Smith, An extended double row layout problem, in: Proceedings of the 12th International Material Handling Research Colloquium, 2012, pp. 535–550. [15] L. Tang, X. Zuo, C. Wang, X. Zhao, A MOEA/d based approach for solving robust double row layout problem, in: IEEE Congress on Evolutionary Computation, 2015, pp. 1966–1973. [16] S. Wang, X. Zuo, X. Liu, X. Zhao, J. Li, Solving dynamic double row layout problem via combining simulated annealing and mathematical programming, Appl. Soft Comput. 37 (2015) 303–310. [17] A.R.S. Amaral, The corridor allocation problem, Comput. Oper. Res. 39 (2012) 3325–3330. [18] A.R.S. Amaral, A parallel ordering problem in facilities layout, Comput. Oper. Res. 40 (2013) 2930–2939. [19] M.E. Yafrani, B. Ahiod, Efficiently solving the traveling thief problem using hill climbing and simulated annealing, Inf. Sci. 432 (2018) 231–244. [20] M. Avci, S. Topaloglu, A multi-start iterated local search algorithm for the generalized quadratic multiple knapsack problem, Comput. Oper. Res. 83 (2017) 54–65. [21] S.L. Goh, G. Kendall, N.R. Sabar, Improved local search approaches to solve the post enrolment course timetabling problem, Eur. J. Oper. Res. 261 (2017) 17–29. [22] A. Megahed, M. Goetschalckx, A modeling framework and local search solution methodology for a production-distribution problem with supplier selection and time-aggregated quantity discounts, Appl. Math. Model. 68 (2019) 198–218. [23] N.R. Sabar, G. Kendall, An iterated local search with multiple perturbation operators and time varying perturbation strength for the aircraft landing problem, Omega 56 (2015) 88–98. [24] U. Benlic, J.K. Hao, Breakout local search for the max-cut problem, Eng. Appl. Artif. Intell. 26 (2013) 1162–1173. [25] J. Guan, G. Lin, Hybridizing variable neighborhood search with ant colony optimization for solving the single row facility layout problem, Eur. J. Oper. Res. 248 (2016) 899–909. [26] J. Kennedy, R. Eberhart, Particle swarm optimization, in: Proceedings of the IEEE International Conference on Neural Networks, 1995, pp. 1942–1948.

J. Guan, G. Lin and H.-B. Feng et al. / Applied Mathematical Modelling 77 (2020) 963–979

979

[27] Z. Zhan, J. Zhang, Y. Li, H.S. Chung, Adaptive particle swarm optimization, IEEE Trans. Syst. Man Cybern. Part B (Cybern.) 39 (2009) 1362–1381. [28] J.J. Liang, A.K. Qin, P.N. Suganthan, S. Baskar, Comprehensive learning particle swarm optimizer for global optimization of multimodal functions, IEEE Trans. Evol. Comput. 10 (2006) 281–295. [29] R.P. Singh, V. Mukherjee, S.P. Ghoshal, Particle swarm optimization with an aging leader and challengers algorithm for the solution of optimal power flow problem, Appl. Soft Comput. 40 (2016) 161–177. [30] X. Wang, H. Zhang, S. Fan, H. Gu, Coverage control of sensor networks in iot based on RPSO, IEEE Internet Things J. 5 (2018) 3521–3532. [31] S. Gheisari, M.R. Meybodi, BNC-PSO: structure learning of bayesian networks by particle swarm optimization, Inf. Sci. 348 (2016) 272–289. [32] D. Ghosh, R. Kothari, Population Heuristics for the Corridor Allocation Problem, IIMA Working Papers 98(2012) 33–40. [33] I.G. Tsoulos, A. Stavrakoudis, Enhancing PSO methods for global optimization, Appl. Math. Comput. 216 (2010) 2988–3001. [34] Z. Li, W. Wang, Y. Yan, Z. Li, PS-ABC: a hybrid algorithm based on particle swarm and artificial bee colony for high-dimensional optimization problems, Expert Syst. Appl. 42 (2015) 8881–8895. [35] F. Marini, B. Walczak, Particle swarm optimization (PSO). A tutorial, Chemom. Intell. Lab. Syst. 149 (2015) 153–165. [36] Y.-B. Shin, E. Kita, Search performance improvement of particle swarm optimization by second best particle information, Appl. Math. Comput. 246 (2014) 346–354. [37] A. Farahmand, S. Payan, S.M.H. Sarvari, Geometric optimization of radiative enclosures using PSO algorithm, Int. J. Therm. Sci. 60 (2012) 61–69. [38] D.M. Simmons, One-dimensional space allocation: an ordering algorithm, Oper. Res. 17 (1969) 812–826. [39] M.F. Anjos, A. Kennings, A. Vannelli, A semidefinite optimization approach for the single-row layout problem with unequal dimensions, Discret. Optim. 2 (2005) 113–122. [40] Q.-K. Pan, R. Ruiz, An effective iterated greedy algorithm for the mixed no-idle permutation flowshop scheduling problem, Omega 44 (2014) 41–50. [41] G.L. Cravo, A.R.S. Amaral, A GRASP algorithm for solving large-scale single row facility layout problems, Comput. Oper. Res. 106 (2019) 49–61.