Next Article in Journal
Network Inference from Gene Expression Data with Distance Correlation and Network Topology Centrality
Next Article in Special Issue
Optimization of the Weighted Multi-Facility Location Problem Using MS Excel
Previous Article in Journal
An Investigation of Alternatives to Transform Protein Sequence Databases to a Columnar Index Schema
 
 
Font Type:
Arial Georgia Verdana
Font Size:
Aa Aa Aa
Line Spacing:
Column Width:
Background:
Article

A Primal-Dual Interior-Point Method for Facility Layout Problem with Relative-Positioning Constraints

1
Graduate School of Creative Science and Engineering, Waseda University, Tokyo 169-8050, Japan
2
School of Creative Science and Engineering, Waseda University, Tokyo 169-8050, Japan
*
Author to whom correspondence should be addressed.
Submission received: 29 December 2020 / Revised: 2 February 2021 / Accepted: 6 February 2021 / Published: 13 February 2021
(This article belongs to the Special Issue Metaheuristic Algorithms and Applications)

Abstract

:
We consider the facility layout problem (FLP) in which we find the arrangements of departments with the smallest material handling cost that can be expressed as the product of distance times flows between departments. It is known that FLP can be formulated as a linear programming problem if the relative positioning of departments is specified, and, thus, can be solved to optimality. In this paper, we describe a custom interior-point algorithm for solving FLP with relative positioning constraints (FLPRC) that is much faster than the standard methods used in the general-purpose solver. We build a compact formation of FLPRC and its duals, which enables us to establish the optimal condition very quickly. We use this optimality condition to implement the primal-dual interior-point method with an efficient Newton step computation that exploit special structure of a Hessian. We confirm effectiveness of our proposed model through applications to several well-known benchmark data sets. Our algorithm shows much faster speed for finding the optimal solution.

1. Introduction

1.1. Facility Layout Problem

Facility Layout Problem (FLP) is one of the most fundamental issues in the design of production system. The goal of FLP is “locating the different facilities or departments in a plant in order to achieve the greatest efficiency in the production of a product or service” (Tompkins et al. 2010 [1]). In an industrial situation, putting in a poor layout may well lead to serious blunders in the use of a company’s available land, in costly re-arrangements in actually tearing down buildings, walls, or major structures which are still usable but which subsequently turn out to be roadblocks to efficiency and low-cost operation [2]. Therefore, it is important for the facility planners to develop the facility layouts that avoid such mistakes.

1.2. Related Research

FLP is a well studied problem and a number of different models and algorithms have been developed. See the surveys Levary and Kalchik(1985) [3], Kusiak and Heragu (1987) [4], Hassan (1994) [5], Meller and Gau (1996) [6], Drira et al.(2007) [7], Arikaran et al. (2010) [8], Anjos and Vieira (2017) [9], and Hosseini-Nasab et al. (2018) [10].
Early research in this field include Koopmans and Beckmann (1957) [11], where the FLP is formulated as QAP (Quadratic Assignment Problem), by approximating the problem as model of assigning all facilities to different locations with the goal of minimizing the sum of the distances multiplied by the corresponding flows. In the 1990s, several exact formulation was presented by representing the candidate layout continuously (called ‘continuous representation’), as opposed to the QAP formulation represents the layout candidate by a grid structure (called ‘discrete representation’). According to Liu and Meller (2007) [12], “By representing the FLP in a discrete fashion, the FLP is simplified, but at the penalty of eliminating many solutions from consideration. Continuous representation is more accurate and realistic than discrete representation, and thus is capable of finding the ‘real optimal’ final layout solution.” However, the continuous representation also increases the complexity of the FLP. In the optimization point of view, FLP belongs to the non-convex problem, and, unfortunately, there are no effective methods for solving the problem in this class. Methods for the general non-convex problem therefore take several different approaches: NLP (NonLinear Programming)-based approach; MH (Meta-Heuristics)-based approach; and MIP (Mixed-Integer-Programming)-based approach, each of which involves some compromise.
In the NLP-based approach, FLP is represented by purely continuous variables and optimized by NLP techniques. The compromise is to give up seeking the optimal solution, and instead they seek a point that is only a locally optimal solution. As a result, the objective value of the local solution obtained is highly dependent on the choice of an initial point. Hence, much of the research is focused on a method for producing a good enough initial point. Those research include Van Camp, Carter and Vaneli (1992) [13], Anjos and Vaneli (2002) [14], Anjos and Vaneli (2006) [15], Jankovits, Luo, Anjos and Vaneli (2011) [16], and Anjos and Vieira (2016) [17]. The general framework is based on the combination of two models: The first is a relaxation of the layout problem and intended to find a good initial point for the iterative algorithm used to solve the second model; the second model is an exact formulation of the facility layout problem and intended to find locally optimal solution. These methods can be fast, can handle large-scale problems, but still have the disadvantage of initial-dependency even with the relaxation step, since the relaxed model that has been proposed is not convex as well.
In the MH-based approach, FLP is represented by a string of finite length integer-variables and optimized by MH techniques (e.g., Simulated Annealing or Genetic Algorithm). The compromise is the quality of solution; an approximate solution is found by drawing some number of candidates from a probability distribution, and taking the best one found as the approximate solution. This approach includes LOGIC and FLEX-BAY. LOGIC is an improvement type algorithm developed by [18], where a candidate layout is represented as a slicing tree. A slicing tree is a binary tree which shows the process of recursive partitioning rectangular in such a way that each rectangular partition corresponds to the space allocated to a facility. In the LOGIC, the slicing tree is optimized using a Genetic Algorithm, and further extensions are conducted by Gau and Meller (1999) [19], Shayan and Chittilappilly (2004) [20], Scholz, Petrick, and Domschke (2009) [21], and Kang and Chae (2017) [22]. FLEX-BAY is an improvement-type algorithm based on a continuous representation developed by Tate and Smith [23]. A layout is represented by a flexible number of vertical bays of varying width, each divided into one or more rectangular departments. Encoding flexible bay layouts is a two-part representation: permutation of the departments and breakpoints for the bays. FLEX-BAY utilizes a genetic algorithm to search the solution space by varying department-to-bay assignments or by adding or removing a bay breakpoint. A family of this approach are Kulturel-Konak, Konak, Norman and Smith(2006) [24], Wong and Komarudin (2010) [25], Kulturel-Konak and Konak (2010) [26], Kulturel-Konak and Konak (2011) [27], Palomo Romero et al. (2017) [28], Garcia-Hernandez et al. (2019) [29], and Garcia-Hernandez et al. (2020) [30]. Recent papers study the multi-objective model. Liu et al. (2018) [31] developed a multi-objective model to minimize material handling cost, to maximize the total adjacency value, and to maximize space utilization. They applied a multi-objective particle swarm optimization. Guan et al. (2019) [32] developed a multi-objective model to minimize material handling cost, to minimize the number of available workshops required for departments, and to maximize space utilization. They also applied a multi-objective particle swarm optimization. Another research trend is the dynamic and stochastic model. Pourvaziri and Pierreval (2017) [33] developed a dynamic facility layout problem based on an open queuing network theory. Turanoğlu and Akkaya (2018) [34] introduced a hybrid heuristic algorithm based on bacterial foraging optimization. Tayal et al. (2017) [35] developed a stochastic dynamic facility layout problem. Several studies integrated these multi-objective and dynamic models. Azevedo et al. (2017) [36] developed a multi-objective model to minimize material handling cost, to minimize reconfiguration cost, and to minimize unsuitability between departments and location. Pourhassan and Raissi (2017) [37] developed a multi-objective model to minimize the total material handling cost and machine rearrangement costs, and to minimize the number of possible transporters accident that can be evaluated by the simulation. They proposed a simulation-based optimization framework in which the genetic algorithm was applied to find an optimal arrangement. While MH-based methods work very well in practice for larger sized problem instances, these approaches have a structural disadvantage in which all these heuristics cannot consider all-feasible solutions due to the encoding scheme that represents solutions as strings for the use of MH techniques.
In the MIP-based approach, FLP is represented by a combination of the binary variables that specify the relative positioning of each pair of departments and the continuous variables that denote the department positioning and shapes. The outline of a general MIP-based-method is as follows. It alternates between two steps: determining a relative position of departments, and the optimizing department positions and the shapes under the specified relative positioning. If the relative positioning of the department is specified, FLP can be reduced to the convex optimization problems such as LP (Linear Programming), QP (Quadratic Programming), or SDP (SemiDefinite Programming), and, thus, can be solved by several convex optimization techniques developed in the past few decades (e.g., simplex method, active-set method, interior-point method). Unlike other two approaches which may miss the optimal solution, MIP-based approach is the true global optimization approach, and, thus, is absolutely reliable. The compromise is, however, the scalability; the first MIP-based method presented by Montreuil in 1990 [38] can only solve the problem with six or less departments. This model is now referred to as FLP0 and several modification for MIP-FLP formulation has been conducted such as FLP1 presented by Meller et al. (1999) [39]; FLP2 by Sherali et al. (2003) [40]; and ε -accurate scheme for controlling the department area by Castillo and Westerlund (2005) [41]; a sequence-pair representation by Liu and Meller (2007) [12]; a graph–pair representation by Bozer and Wang (2012) [42]. Even with those improvements, however, the problem size that can be solved to optimality is still less than twelve (Chae and Regan 2016 [43]).

1.3. Research Purpose

Corresponding to the above-mentioned two steps of MIP, there are two approaches for speeding up MIP: searching the relative position more efficiently, or optimizing department positions and the shapes under the specified relative positioning faster. In this study, we focus primarily on the latter class approach.
In this paper, we describe a custom interior-point method for solving the FLP with relative positioning constraints (FLPRC) that is substantially faster than the standard methods used in the general-purpose solver. We build a compact formation of FLPRC and its dual, which enables us to establish the optimal condition very quickly. We use this optimality condition to implement the primal-dual interior-point method with an efficient Newton step computation that exploits the special structure of a Hessian. The computation effort of our method scales with linear-order of the number of departments, whereas that of standard methods scales with cubic-order of the number of departments.
The outline of this paper is as follows. In Section 2, we describe the FLPRC, formulated as a linear programming in a matrix form. In Section 3, we describe the overall algorithm. We first derive a dual problem and the Karush–Kuhn–Tucker (KKT) optimal condition of the model. We use this KKT condition to implement our optimization method with a primal-dual interior-point. We describe the barrier subproblem associated with a primal-dual interior-point method for the FLPRC, and we show how the special structure of the FLPRC can be exploited to compute the search direction very efficiently. In Section 4, we give numerical results to illustrate the effectiveness of proposed method. Finally, we provide our conclusions and discuss some of the future research in Section 5.

2. Facility Layout Problem with Relative Positioning Constraints

In this section, we first review the FLPRC based on the FLP2 presented by [40]. We then describe our modified model and how our model can be faster than the FLP2.

2.1. The FLPRC Formulation

We consider N departments i = 1 , , N , which have centroid positions ( x i , y i ) R 2 and shapes ( w i , h i ) R 2 of half of their widths and heights, for i = 1 , , N . They must be placed in a rectangle area with width W and height H, and lower left corner at the position (0, 0), under specified relative positioning.
The objective of FLP is to find a set of positions and shapes of department r i = [ x i , y i , w i , h i ] T that yield as small material handling costs as possible. In the FLP2, the material handling cost is measured as the product of distance times flows between departments as in (1)
minimize i = 1 N j = i + 1 N f i j d i j ,
where f i j denotes the material flows from department i to department j and d i j denotes the distance between centroids of department i and department j. The choices of distance measures are usually either rectilinear or Euclidean, but, in the FLP2, the rectilinear distances are used as in (2).
d i j = d i j x + d i j y = | x i x j | + | y i y j |
The absolute values in the distance function (2) can be linearized as in (3)
d i j x x i x j d i j x d i j y y i y j d i j y
We call these inequalities the distance constraints.
In the FLP, it is required that the departments lie inside the bounding rectangle with width W and height H as in (4) and (5).
w i x i W w i , i = 1 , , N
h i y i H h i , i = 1 , , N
These inequalities are called the within-boundary constraints.
The aspect-ratio constraints require that the aspect ratio of a department is bounded within upper and lower limits, i.e.,
l i h i / w i u i , i = 1 , , N ,
where aspect ratio of a department is defined by h i / w i and the upper and lower bound of aspect ratio are defined by u i , l i , respectively. Multiplying both sides of each inequality by w i reduces to these constraints into linear inequalities as in (7).
l i w i h i u i w i , i = 1 , , N ,
In case the dimension of a department is known a priori (often called Machine Layout Problem, or MLP), we can use the fixed aspect ratio with l i = u i , which results in a linear equality constraint.
Let s i be the quarter of minimal area of department i for i = 1 , , N , thus the area constraints require Equation (8):
w i h i s i , i = 1 , , N .
Instead of these exact constraints, we use an outer polyhedra approximation, proposed by [40], to obtain the speed-up as (9):
s i w i + 4 w ¯ i m 2 h i 2 s i w ¯ i m , i = 1 , , N , m = 1 , , M .
where affine support points are defined by w ¯ i m and the number of affine support points are defined by M, respectively. By using the constraints (9), instead of (8), the area constraints can be reduced to the standard linear inequalities, and thus can be solved much faster.
The relative-positioning constraints require that, for each pair of departments ( i , j ) with i j , the department i must be placed to the left, right, above, or below the department j. These relations are specified with the binary variables H (meaning ‘horizontal’) and V (meaning ‘vertical’) as in (10) and (11):
H i j = 1 ( if   i   is   left   of   j ) 0 ( o t h e r w i s e )
V i j = 1 ( if   i   is   below   j ) 0 ( o t h e r w i s e )
Using these variables, the following linear inequalities are imposed as in (12) and (13):
x i + w i x j w j if H i j = 1
y i + h i y j h j if V i j = 1
Although we will not determine the relative positioning, we mention that the following condition must hold to ensure that the departments do not overlap (14):
H i j + H j i + V i j + V j i = 1
These constraints are included in the MIP-based FLP and make the problem a complicated combinational problem.
We summarize the overall problem formulation as in (15).
minimize i = 1 N j = i + 1 N f i j { d i j x + d i j y }
s u b j e c t t o d i j x x i x j d i j x , i , j = 1 , , N .
d i j y y i y j d i j y , i , j = 1 , , N .
w i x i W w i , i = 1 , , N .
h i y i H h i , i = 1 , , N .
l i w i h i u i w i , i = 1 , , N ,
x i + w i x j w j i f H i j = 1 , i , j = 1 , , N .
y i + h i y j h j i f V i j = 1 , i , j = 1 , , N .
s i w i + 4 w ¯ i m 2 h i 2 s i w ¯ i m , i = 1 , , N , m = 1 , , M .
  • N: number of departments;
  • i , j : department indices ( i , j = 1 N );
  • f i j : material flow between department i and department j;
  • W , H : width and height of facility;
  • H i j , V i j : binary constants to specify the relative positioning of department i and j in the horizontal and vertical direction;
  • l i , u i : lower and upper limits on the aspect ratio of department i;
  • s i : quarter of minimum area of department i;
  • d i j x , d i j y : rectilinear distance between department i and j in x and y direction;
  • x i , y i : x and y coordinates of centroid of department i;
  • w i , h i : half length of width and height of department i;

2.2. Reducing Redundant Inequalities by Exploiting Relative-Positioning Constraints Structure

In this section, we describe the redundant inequalities that can be eliminated by exploiting relative-positioning constraints structure. It yields more sparsity of the matrix in the linear equation to compute the Newton step, which is the main effort of the algorithm, and thus yields the speeding-up the algorithm.

2.2.1. Redundancy in Distance Constraints

For each pair of departments ( i , j ) with i j , two inequalities in (12) and (13) can be reduced to a linear equality, if H i j and V i j are specified as follows:
d i j x x i x j d i j x d i j x = x j x i ( if H i j = 1 ) d i j y y i y j d i j y d i j y = y j y i ( if V i j = 1 )
Substituting these equalities into d i j x , d i j y in the objective function (1), these variables and equalities disappeared.
Assuming the non-overlapping condition (14) holds, at least N C 2 variables and 2 × N C 2 inequalities can be reduced. We can further reduce the distance variables and constraints in the similar way, for each pair of departments ( i , j ) with f i j = 0 .

2.2.2. Redundancy in Within-Boundary Constraints

The within-boundary constraints are also reduced by using the relative-positioning structure. The constraints w i x i in (4) only needs to be imposed on the left-most departments, i.e., for i s.t. j = 1 N H i j = 0 . In the similar way, the minimal set of within-boundary constraints can be reduced as follows:
w i x i for   i   s . t . j = 1 N H j i = 0 x i W w i for   i   s . t . j = 1 N H i j = 0 h i y i for   i   s . t . j = 1 N V j i = 0 y i H h i for   i   s . t . j = 1 N V i j = 0
Suppose N l e f t , N r i g h t , N u p p e r , N l o w e r are the number of the most left, right, upper, and lower departments, respectively, the number of inequalities can be reduced to N l e f t + N r i g h t + N u p p e r + N l o w e r . Assuming that the non-overlapping condition (14) holds, this number becomes, at most, less than 2 ( N + 1 ) , which is less than the original 4 N inequalities.

2.2.3. Redundancy in Relative-Positioning Constraints

Finally, we mention the redundant inequalities in relative-positioning constraints. It is obvious that, if the conditions H i j = 1 and H j k = 1 hold, then the condition H i k = 1 must hold. Using these relations, we can reduce the redundant inequalities.
A minimal set of relative-positioning constraints can easily be obtained through eliminating the indirect path in H i j or V i j by the Warshall–Floyd algorithm.

2.3. The FLPRC Formulation in a Matrix Form

In this section, we introduce simpler formulation in a matrix form. This formulation does not reduce the total number of variables nor constraints, but it removes many row/column permutation to solve the linear equation to compute the Newton step.

2.3.1. Notation

First, we describe the notation in this section. 0 m denotes the ( m × 1 ) vector with all entries zero, 1 m denotes the ( m × 1 ) vector with all entries one. I means the identity matrix, and diag(a) means diagonal matrix with diagonal entries a 1 a m , where a i is i-th entry of a .

2.3.2. Decision Vector

Let x , y , w , h R N denote the vector with i t h entry x i , y i , w i , h i , respectively. We define the layout vectorr R 4 N as
r = x y w h R 4 N
and the distance vector d = [ d x T , d y T ] T R N ( N 1 ) as
d x = d 1 x d 2 x d K x = d 12 x d 13 x d ( N 1 ) N x , d y = d 1 y d 2 x d K y = d 12 y d 13 y d ( N 1 ) N y
where d k x = d i j x with k = ( i 1 ) / 2 × ( 2 N i ) + ( j i ) . We denote the overall set of decision variables v R 4 N + ( N 1 ) N as
v = r d .

2.3.3. Objective Function

We define the flow vector f R N ( N 1 ) / 2 as
f = f 1 f 2 f K = f 12 f 13 f ( N 1 ) N
and the overall cost vector c = [ 0 4 N T , f T , f T ] T . Using these vectors, we can write the objective function (1) as in (16):
i = 1 N j = i + 1 N f i j d i j = f 1 d 1 + + f K d K   = f 1 ( d 1 x + d 1 y ) + + f K ( d K x + d K y )   = f T d x + f T d y   = 0 f f T r d x d y   = c T v
If a pair of department ( i , j ) can be reduced with H i j , V i j from the redundancy in distance constraints, then the cost vector becomes c = [ c r T , f T , f T ] T with
( c r ) k = f i j i f k = ( i 1 ) / 2 × ( 2 N i ) + ( j i ) f i j i f k = ( j 1 ) / 2 × ( 2 N j ) + ( i j ) 0 otherwise
and the kth entry corresponding ( i , j ) in f becomes 0.

2.3.4. Distance Constraints

We define the distance constraints matrix A d R N ( N 1 ) × 4 N as
A d = A d x O O O O A d y O O ,
where A d x = A d y are defined as
A d x = A d y = A d 1 A d 2 A d N ,
and A d i is defined as
A d i = O ( N i + 1 ) × ( i 1 ) 1 ( N i + 1 ) I ( N i + 1 ) × ( N i + 1 ) .
Using these matrices, we can write the distance constraints (3) as in (18).
d i j x x i x j d i j x , i < j d i j y y i y j d i j y , i < j x i x j d i j x 0 y i y j d i j y 0 ( x i x j ) d i j x 0 ( y i y j ) d i j y 0 , i < j . A d x x d x 0 A d y y d y 0 A d x x d x 0 A d x y d y 0 A d r d 0 A d r d 0 A d I A d I r d 0
If there are redundant constraints that can be eliminated, the corresponding rows are eliminated. Instead, the corresponding term is represented by the cost vector (17).

2.3.5. Within-Boundary Constraints

We define the within-boundary constraints matrix A w b R 4 N × 4 N and the within-boundary constraints vector b w b R 4 N as follows:
A w b = I I I I , b w b = 0 2 N W × 1 N H × 1 N
Using these matrices, the within-boundary constraints can be written as in (19).
w i x i W w i , i = 1 , , N . h i y i H h i , i = 1 , , N . x i + w i 0 , i = 1 , , N . y i + h i 0 , i = 1 , , N . x i + w i W , i = 1 , , N . y i + h i H , i = 1 , , N . I O I O O I O I I O I O O I O I x y w h 0 N 0 N W × 1 N H × 1 N A w b r b w b .
If there are redundant constraints that can be eliminated, the entries of corresponding department become 0.

2.3.6. Aspect-Ratio Constraints

The aspect-ratio constraints (7) can be written as
l i w i h i 0 , i = 1 , , N , u i w i + h i 0 , i = 1 , , N .
We can express these bounds in a matrix form as
O O L I O O U I x y w h 0 2 N 0 2 N
or simply we can write as in (20)
A a r r b a r
where A a r R 4 N × 4 N , L R N × N , U R N × N , b a r R 4 N are defined as follows:
A a r = O O L I O O U I , L = diag ( l ) , U = diag ( u ) , b a r = 0 4 N

2.3.7. Area Constraints

The approximated area constraints
s i w i 4 w ¯ i m 2 h i 2 s i w ¯ i m , i = 1 , , N , m = 1 M .
can be written as follows:
O O S W ¯ m x y w h 2 × S W ¯ m × 1 N , m = 1 , , M
This can be compactly written as follows:
A s r b s ,
where A s R M N × 4 N , A s m R N × 4 N , b s M R M N , b s m R N , S R N × N , W ¯ m R N × N , s R N , w m ¯ R N are defined as follows:
A s = A s 1 A s M , A s m = O O S W ¯ m , b s M = b s 1 b s M , b s m = 2 × S W ¯ m × 1 N , S = diag ( s ) , W ¯ m = 4 diag ( w ¯ m ) 2 , s = s 1 , , s N T , w ¯ m = w ¯ 1 m , , w ¯ N m T .

2.3.8. Relative-Positioning Constraints

The relative-positioning constraints in a matrix form can be written as
A r p r b r p
where A r p R N ( N 1 ) × 4 N , b r p R N ( N 1 ) are defined as follows:
A r p = A r p x w A r p w h ( A r p x w ) m n = 1 i f ( m = ( i 1 ) / 2 × ( 2 N i ) + ( j i ) ) ( n = i , i + 2 N , j + 2 N ) ( H i j = 1 ) 1 i f ( m = ( i 1 ) / 2 × ( 2 N i ) + ( j i ) ) ( n = j ) ( H i j = 1 ) 0 , o t h e r w i s e ( A r p y h ) m n = 1 i f ( m = ( i 1 ) / 2 × ( 2 N i ) + ( j i ) ) ( n = i + N , i + 3 N , j + 3 N ) ( V i j = 1 ) 1 i f ( m = ( i 1 ) / 2 × ( 2 N i ) + ( j i ) ) ( n = j + N ) ( V i j = 1 ) 0 otherwise b r p = 0 M
For a department pair ( i , j ) with H i j = 1 , in mth row corresponding ( i , j ) , the coefficient of x i , w i , w j becomes 1 and that of x j becomes −1. Similarly, for a department pair ( i , j ) with V i j = 1 , in mth row corresponding ( i , j ) , the coefficient of y i , h i , h j becomes 1 and that of y j becomes −1. We can eliminate the row corresponding the redundant inequalities that can be reduced.

2.3.9. Overall Problem Formulation in Matrix Form

We define the overall constraint matrix A R ( 5 / 2 N 2 + ( 11 / 2 + M ) N ) × ( 4 N + N ( N 1 ) ) and constraint vector b R ( 5 / 2 N 2 + ( 11 / 2 + M ) N ) as
A = A d I A d I A r O , b = 0 0 b r
and A r R ( 1 / 2 N 2 + ( 15 / 2 + M ) N ) × ( 4 N + N ( N 1 ) ) and b r R ( 1 / 2 N 2 + ( 15 / 2 + M ) N ) are defined as follows:
A r = A w b A a r A s A r p , b r = b wb b ar b s b rp
Using these matrices and vectors, we can write the overall constraints as
A d I A d I r d 0 A w b r b w b A a r r b a r A s r b s A r p r b r p A v b
Thus, we can summarize the overall problem formulation as:
m i n i m i z e c T v s u b j e c t t o A v b
This problem (23) is the standard Linear Programming (LP), and, thus, can be solved to optimality.

2.4. Dual Problem and Optimality Conditions for FLPRC

In this section, we derive a dual problem for the FLPRC (23) and outline the optimality condition. We introduce Lagrangian Multiplier vectors λ = [ λ 1 , λ ( 5 / 2 N 2 + ( 11 / 2 + M ) N ) ] T R ( 5 / 2 N 2 + ( 11 / 2 + M ) N ) for the constraint corresponding to a i T v b i , where a i T is i-th row vector of the constraints matrix A. We let n λ = dim ( λ ) = ( 5 / 2 N 2 + ( 11 / 2 + M ) N ) denote the dimensions of λ . The Lagrangian is then
L ( v , λ ) = c T v + i = 1 n λ λ i ( a i T v b i ) = c + i = 1 n λ λ i a i T v λ T b
To obtain the dual function, we minimize L over the primal variables v . We find that the minimum is −∞, unless c + i = 1 n λ λ i a i 0 . The dual function is therefore given by
i n f L ( v , λ ) = inf { c + i = 1 n λ λ i a i T v λ T b } = λ T b i f c + i = 1 n λ λ i a i 0 otherwise
Thus, we can write the dual of the FLPRC as
maximize λ T b subject to A T λ + c = 0 λ 0
The optimality condition, i.e., the Karush–Kuhn–Tucker (KKT) condition, is given as
x L ( v , λ ) = c + A T λ = 0
diag ( λ ) ( b A v ) = 0
A v b 0
λ 0
The first equation (26a) that follows from the gradient of Lagrangian must be zero, since the optimal point v* minimizes the L ( v , λ ) over v . The second equation (26b) is the complementary slickness. The third and the last inequality (26c) and (26d) are the constraints with respect to the primal and dual feasibility, respectively.
If (and only if) there are ( v * , λ * ) that satisfies the KKT condition (26), then the v * is the optimal point. In other words, solving LP (23) is equivalent to the problem of finding the point ( v * , λ * ) that satisfies the KKT conditions. As is the case in most other convex problems, however, it is difficult to directly find such points due to a large number of linear inequalities to maintain the primal and dual feasibility in the KKT conditions (26).
The interior-point method is a technique used to solve those inequality-constrained problems by reducing it to a sequence of problems without inequalities. The details are described in the next section.

3. Custom Interior-Point Method for Solving FLPRC

In this chapter, we describe an efficient method for solving FLPRC (23). The method is primal-dual interior-point method, which is one of the most powerful techniques for solving Linear Programming (LP) or other classes of convex optimization programming. We propose an efficient Newton step computation that exploit the special structure of the problem. Our algorithm is also applicable to other types of interior-point methods with a minor change.

3.1. Barrier Subproblem

The main idea of the interior-point method is casting the inequality constraints into the objective function as
minimize c T v subject   to A v b minimize c T v + 1 t ϕ ( v )
where t is the scaling parameter, and ϕ ( v ) is defined as
ϕ ( v , λ ) = i = 1 n λ ϕ i ( v )
ϕ i ( v , λ ) = l o g ( b i a i T v ) .
The problem is called a barrier subproblem and the function ϕ is called a barrier function. This function is convex and smooth, and, thus, the barrier subproblem becomes convex if the original problem is convex. The first derivative of this function is
ϕ ( v , λ ) = i = 1 n λ a i b i a i T v
We define z with elements z i = 1 / ( b i a i T v ) , and we have
ϕ ( v , λ ) = A T z
The optimality condition for the barrier subproblem is
( c T v + 1 t ϕ ( v ) ) = c + 1 t A T z = 0 .
Comparing this with the gradient of Lagrangian (26a),
x L ( v , λ ) = c + A T λ = 0
we have
λ = 1 t z = 1 t ( b A v ) .
Substituting this into the complementary slackness (26b) and eliminating linear inequalities (26c) and (26d), we have the modified KKT condition:
c + A T λ = 0 diag ( λ ) ( b A v ) = ( 1 / t ) 1
The modified KKT conditions (33) satisfy the following two properties:
  • These conditions are far easier to solve than the KKT condition (26)
  • As t increases, these condition approaches the KKT condition (26).
Combining these two properties, the FLPRC is solved to optimality. The modified KKT conditions (33) are solved by the Newton method, which is described in detail in the next section.

3.2. Newton Method

We rewrite the modified KKT conditions (33) as in (34).
R t ( v , λ ) = R t v ( v , λ ) R t λ ( v , λ ) = c + A T λ diag ( λ ) ( b A v ) ( 1 / t ) 1 = 0
In a primal-dual interior method, we start with some initial point ( v 0 , λ 0 ) and update them in such a way that ( v 0 , λ 0 ) ( v 1 , λ 1 ) ( v k , λ k ) with R t ( v k , λ k ) 0 for k . These updates are made by a Newton method.
At each iteration, primal and dual search directions ( Δ v k , Δ λ k ) and step size γ k are computed, then the new point is obtained as ( v k + 1 , λ k + 1 ) = ( v k + γ k Δ v k , λ k + γ k Δ λ k ) . In the Newton method, the search direction is characterized by the solution of the linear equation
R t ( v k + Δ v k , λ k + Δ λ k ) R t ( v k , λ k ) + D R t ( v k , λ k ) [ Δ v k , Δ λ k ] T = 0
Δ v k Δ λ k = D R t 1 ( v k , λ k ) R t ( v k , λ k )
where
D R t ( v k , λ k ) = v R t v ( v , λ ) λ R t v ( v , λ ) v R t v ( λ , λ ) λ R t λ ( v , λ ) = O A T diag ( λ ) A diag ( b A v )
The first equation (35) follows from the 1st order Taylor approximation. Thus, the solution of (35) may be a good estimate of the solution of Equation (34). After computing the Newton step, we then carry out the back-tracking line search to find the step size γ k .

3.3. Efficient Newton Step Computation

In this section, we show how to compute the Newton step ( Δ v k , Δ λ k ) , i.e., solve Equation (36), efficiently. These equations can be solved using standard methods for linear equations; however, by exploiting the special structure of the equations, they can be solved even faster. The method we describe in this section is based on a sequence of three elimination steps, in which particular blocks of variables are eliminated by the Schur Complement.

3.3.1. 1st Elimination Step

Our first step is to eliminate the variables Δ λ k from Equation (35). From the second equation, we obtain
diag ( λ ) A Δ v k + diag ( b A v ) Δ λ k = R λ Δ λ k = diag ( b A v ) 1 R λ diag ( b A v ) 1 diag ( λ ) A Δ v k Δ λ k = diag ( z ) R λ diag ( z λ ) A Δ v k
The second to third equation is derived using z i = 1 / ( b i a i T v ) . Substituting this into the first equation in (35), we have Δ v k as follows:
A T Δ λ k = R v A T diag ( z ) R λ A T diag ( z λ ) A Δ v k = R v Δ v k = [ A T diag ( z λ ) A ] 1 [ R v + A T diag ( z ) R λ ] Δ v k = [ A T Z L A ] 1 [ R v + A T Z R λ ]
where Z = diag ( z ) , L = diag ( λ ) . These Equations (37) and (38) are alternative for the original Equation (36). Instead of solving the large set of linear Equations (36), we find Δ v k from the smaller Equation (38), which can be easily solved as described in the next two sections, and then calculate the Δ λ from Equation (37).

3.3.2. 2nd Elimination Step

We break down the equation and eliminate some of the blocks further to reduce the size of linear equation. From our definition of v = [ r T , d T ] T , and of A, we can rewrite the Equation (37) as follows:
Δ v k = Δ r k Δ d k = [ A T Z L A ] 1 [ R v + A T Z R λ ] = A d T A d T A r T I I O D d + O O O D d O O O D r A d I A d I A r O 1 g r g d = A d T ( D d + + D d ) A d + A r T D r A r A d T ( D d + D d ) ( D d + D d ) A d D d + + D d 1 g r g d
where D d + , D d , D r are defined as follows
D d + = diag { λ d + T ( A d r d ) 1 } D d = diag { λ d T ( A d r d ) 1 } D r = diag { λ r T ( A r r b r ) 1 }
and g r , g d are defined as follows:
g r = 2 ( A d T λ d + A d T λ d + A r T λ r ) ( 1 / t ) ( A d T z d + T A d T z d T + A r z r T ) g d = [ f T , f T ] T + 2 ( λ d + λ d ) ( 1 / t ) ( z d + z d )
with ( λ d + , λ d , λ r ) are Lagrangian multipliers and ( z d + , z d , z r ) are inverse vectors for the constraints ( A d T r d 0 , A d T r d 0 , A r r b r 0 ) , respectively.
( g r , g d ) are derived from the right-hand-side as (40)
[ R v + A T Z R λ ] = [ ( c + A T λ ) + A T Z ( L z 1 ( 1 / t ) 1 ) ] = [ c + A T ( λ + Z L z 1 ( 1 / t ) z ) ] = [ c + A T ( λ + L Z z 1 ( 1 / t ) z ) ] = [ c + A T ( λ + L 1 ( 1 / t ) z ) ] = [ c + A T ( λ + λ ( 1 / t ) z ) ] = [ c + A T ( 2 λ ( 1 / t ) z ) ] = 0 f f + A d T A d T A r T I I O 2 λ d + λ d λ r ( 1 / t ) z d + z d z r = g r g d
From the second Equation in (40), we have Δ d k as (41).
g d = ( D d + D d ) A d Δ r k + ( D d + + D d ) Δ d k Δ d k = ( D d + + D d ) 1 g d + ( D d + + D d ) 1 ( D d + D d ) A d Δ r k
Substituting this into the first equation in (40), we have Δ r as follows:
g r = A d T ( D d + + D d ) A d + A r T D r A r Δ r k A d T ( D d + D d ) Δ d k = A d T ( D d + + D d ) A d Δ r k + A r T D r A r Δ r k + A d T ( D d + D d ) ( D d + + D d ) 1 g d A d T ( D d + D d ) ( D d + + D d ) 1 ( D d + D d ) A d Δ r k
The sum of the first and forth terms can be summarized as follows:
A d T ( D d + + D d ) ( D d + D d ) ( D d + + D d ) 1 ( D d + D d ) A d r k = A d T ( D d + + D d ) 1 ( D d + + D d ) 2 ( D d + D d ) 2 A d Δ r k = A d T 4 D d + D d ( D d + + D d ) 1 A d Δ r k
We define D d and g as follows:
D d = 4 D d + D d ( D d + + D d ) 1 g = g r + A d T ( D d + D d ) ( D d + + D d ) 1 g d
Using these expressions, we have Δ r as in (42).
( A d T D d A d + A r T D r A r ) Δ r k = g

3.3.3. Third Elimination Step

Equation (42) can be further broken down into smaller pieces, and solved efficiently. We rewrite Equation (42) as in (43)
E Δ r k = g E = E d + E r E d = A d T D d A d E r = A r T D r A r
The term A r T D r A r can be factorized as in Equation (44)
A r T D r A r = A w b T A a r T A s T A r p T D w b O O O O D a r O O O O D s O O O O D r p A w b A a r A s A r p = A w b T D w b A w b + A a r T D a r A a r + A s T D s A s + A r p T D r p A r p = E w b + E a r + E s + E r p ,
where D w b , D a r , D s , D r p are the diagonal matrices with entries with respect to within-boundary, aspect-ratio, area, and relative-positioning constraints. We decompose Δ r k with respect to ( Δ x , Δ y ) and ( Δ w , Δ h ) as Δ r k = [ Δ r ( x y ) k T , Δ r ( w h ) k T ] T . Factorizing (43) into block elements with respect to ( Δ x , Δ y ) and ( Δ w , Δ h ) vectors, we have the following Equation (45)
E 11 E 12 E 21 E 22 Δ r ( x y ) k Δ r ( w h ) k = g x y g w h ,
We discuss the ( m , n ) block element E m n , m = 1 , 2 , n = 1 , 2 . We define the ( m , n ) block element of E d , E r , E w b , E a r , E s , E r p as E m n d , E m n r , E m n w b , E m n a r , E m n s , E m n r p .
First, E m n w b , E m n a r , E m n s , m = 1 , 2 , n = 1 , 2 are all diagonal matrices because the block elements of A w b , A a r , D w b , D a r are diagonal matrices. As for E m n r p , the block elements of A r p is a sparse matrix, and D r p is a diagonal matrix, so E m n r p , m = 1 , 2 , n = 1 , 2 is a sparse symmetric matrix. Since E m n r , m = 1 , 2 , n = 1 , 2 is the sum of these diagonal matrices and a sparse symmetric matrix, E m n r , m = 1 , 2 , n = 1 , 2 is a sparse symmetric matrix. On the other hand, for E m n d , since it has only the components related to ( x , y ) , E 11 d is a dense matrix, and the other block elements are E 12 d , E 21 d , E 22 d = O . From the above, by taking the sum of each block element, E 11 is a dense matrix, and the other block elements E 12 , E 21 , and E 22 are sparse symmetric matrices.
Speed-up is achieved using this structure. From the first equation in (45), we can obtain Δ r w h as Equation (46):
E 11 Δ r ( x y ) k + E 12 Δ r ( w h ) k = g x y Δ r ( w h ) k = E 12 1 ( g x y E 11 Δ r ( x y ) k ) .
Substituting this into the second equation in (45), we have Equation (47) to get Δ r x y :
E 21 Δ r ( x y ) k + E 22 ( E 12 1 ( g x y E 11 Δ r ( x y ) k ) ) = g w h ( E 21 E 22 E 12 1 E 11 ) Δ r ( x y ) k = E 22 E 12 1 g x y g w h .
Using this elimination, we can compute Δ r ( x y ) k efficiently, we then compute the remaining variables Δ r ( x y ) k , Δ d , Δ λ , which is faster than just solving the larger set of equations (33) using a standard factorization method.
Note that all of the matrices are non-singular, if there exists an optimal solution of problem (15).

3.4. Complexity Analysis

In this section, we analyze the calculation complexity. The complexity is measured by the number of floating operations (flops), and all of the calculation in this section is based on Boyd and Vandenberghe [44]. Let n v = dim ( v ) = 4 N + N ( N 1 ) denote the dimension of v. The size of the coefficient matrix of the original Equation (36) is ( n v + n λ ) . Using the standard LU decomposition, the computational complexity of the Equation (36) is ( 2 / 3 ) ( n v + n λ ) 3 . Expanding this, the amount of calculation is expressed as follows:
( 2 / 3 ) ( n v + n λ ) 3 = 343 8 N 6 + ( 2793 8 + 147 4 M ) N 5 + ( 7581 8 + 399 2 M + 21 2 M 2 ) N 4 + ( 6859 3 + 1083 4 M + 57 2 M 2 + M 3 ) N 3
On the other hand, the calculation effort after 1st elimination is performed is 2 n v 2 n λ + ( 2 / 3 ) n v 3 . In this equation, the complexity is proportional to n λ . Since n λ > n v , the amount of calculation is greatly reduced. The amount of calculation after 1st elimination is performed is as follows:
2 n v 2 n λ + ( 2 / 3 ) n v 3 = 17 3 N 6 + ( 47 + 2 M ) N 5 + ( 129 + 12 M ) N 4 + ( 117 + 18 M ) N 3
The calculation effort in the second term can be further reduced by performing 2nd elimination. Let n r = dim ( r ) = 4 N be the dimensions of r, and let n d = dim ( d ) = N ( N 1 ) denote the dimensions of d. The calculation effort after 1st elimination is performed is ( 2 / 3 ) n v 3 = ( 2 / 3 ) ( n r + n d ) 3 . The calculation effort after 2nd elimination is performed is as shown in 2 n r 2 n d + ( 2 / 3 ) n r 3 . In this equation, the complexity is proportional to n d . Since n d > n r in the range of N > 4 , the amount of calculation is greatly reduced. The amount of calculation after 2nd elimination that is performed is as follows:
2 n v 2 n λ + 2 n r 2 n d + ( 2 / 3 ) n r 3 = 4 5 N 6 + ( 79 5 + 2 M ) N 5 + ( 526 5 + 12 M ) N 4 + ( 329 3 + 18 M ) N 3
The calculation effort in the second term can be further reduced by performing 3rd elimination. Let n x y = dim ( x ) + dim ( y ) = 2 N denote the dimensions of x and y, and n w h = dim ( w ) + dim ( h ) = 2 N denote the dimensions of w and h. The calculation effort after 2nd elimination is performed is as shown in ( 2 / 3 ) n r 3 = ( 2 / 3 ) ( n x y + n w h ) 3 . The calculation effort after 3rd elimination is performed is n w h + 2 n x y 2 n w h + ( 2 / 3 ) n x y 3 . Here, the first term is the amount of calculation when sparse factorization is used by utilizing the fact that the matrix E 22 of Equation (47) is a symmetric and sparse matrix. The third term is a block diagonal matrix. Therefore, it can be decomposed into equations related to x and y and the calculation effort is ( 2 / 3 ) N 3 × 2 = ( 4 / 3 ) N 3 . The amount of calculation after 3rd elimination is performed is as follows:
2 n v 2 n λ + 2 n r 2 n d + n w h + 2 n x y 2 n w h + ( 4 / 3 ) N 3 = 4 5 N 6 + ( 79 5 + 2 M ) N 5 + ( 526 5 + 12 M ) N 4 + ( 52 3 + 18 M ) N 3 + 2 N
These complexity is summarized in Table 1. From the above analysis, the problem can be reduced to solving smaller equations by decomposing the matrix rather than solving the original matrix.

3.5. Starting-Point Using Flexible Bay Structure

In this section, we present a technique to create the starting-point for further speeding-up the computation. While our proposed primal-dual interior point method has enough computation effectiveness in itself, there is a potential advantage in utilizing a ‘good’ starting point from a practical perspective. In the context of using interior point method, a starting point should be strictly feasible. Therefore, we define a ‘good’ starting point as points that are close to optimality, yet is sufficiently far from the boundary of the feasible solution. We describe a technique to create such starting point in this section.
Our method is based on a relaxed Flexible-Bay Structure (FBS), in which the departments are located in parallel bays with the varying width associated the total area of the departments in that bay, using two part representation: permutation of the departments and breakpoints for the bays. Figure 1 presents an example illustrating the bay assignment and the obtained layout from the horizontal and vertical relative-positioning graphs.
In our proposal, the permutation and breakpoints are specified using the relative-positioning H i j , V i j , and the shapes of department ( w i , h i ) and the positioning of departments ( x i , y i ) are determined through the relaxed-FBS decoding scheme. However, those obtained from the FBS procedure is near optimality and feasible, but not strictly feasible in general, since each pair of adjacent departments has no separation between them. Therefore, an additional small clearance ρ (say, 10 5 ) is used to make them strictly feasible. The outline of starting-point generation using the relaxed-FBS is described as Algorithm 1.
Algorithm 1:Starting-point generation using the relaxed-FBS
   given H i j , V i j , ρ
   assign department i to the bay k with the horizontal precedence i = 1 N H i j = k .
   sort department i in the bay k in the order of vertical precedence i = 1 N V i j = k
   get bay area B k as: B k = i k S i
   get bay width W k B as: W k B = max i k ( H / B k , w i m i n )
   get department widths and heights as w i = min ( W k B , w i m a x ) , h i = S i / w i
   get department coordinates as w i = min ( W k B , w i m a x ) , h i = S i / w i
                 x i = w i + ρ , i with i = 1 N H i j = 0 (root node)
                 y i = h i + ρ , i with i = 1 N V i j = 0 (root node)
                 x j = x i + ( w i + w j ) + ρ , ( i , j ) with H i j = 1
                 y j = x i + ( h i + h j ) + ρ , ( i , j ) with V i j = 1 ,
   where ρ R is an additional clearance between departments.

3.6. Overall Algorithm

We summarize the overall algorithm is as Algorithm 2.
Algorithm 2:Overall Algorithm

   given t 0 , A, b, c, ϵ
   find initial points ( v 0 , λ 0 ) by warm-start algorithm.
   repeat 1-4:
   1. Find primal and dual search directions [ Δ v k T , Δ λ k T ] T by Newton equation.
   2. Find step size γ k R by back-tracking
   3. Update primal and dual variables [ v k + 1 T , λ k + 1 T ] : = [ v k T , λ k T ] + γ k [ Δ v k T , Δ λ k T ] T
   4. Update t k t k + 1
   until | | R t ( v k , λ k ) | | ϵ
Theoretically, the convergence proof for the primal-dual interior-point algorithm is given in several books and papers such as Boyd and Vandenberghe (2004) [44]. In practice, it is numerically difficult to converge when ϵ is very small. However, the extreme accuracy is not necessary in the context of the intent of this study. The intent of the research is to provide the decision-making support tool for the facility designers, in the phase of the block layout design which the ”macro” flows in the facility is primarily concerned with, while the detailed layout is often concerned with the ”micro” flows (Tompkins et al. (2010) [1]).

4. Numerical Experiments

4.1. Experimental Overview

In this chapter, we give some numerical examples to illustrate the effectiveness of the algorithm described before. All of our examples use the data following from FLP benchmark problem instances, summarized in Liu and Meller (2007) [12]. We generated the randomized relative positioning constraints for each problem, using the following two steps: giving the randomly distributed location first, then specifying the H i j , V i j for i , j with i j as follows:
  • if | x i x j | | y i y j | , H i j = 1
  • otherwise V i j = 1
We abandoned infeasible combinations.
The FLPRC is a convex problem and the optimal solution can be obtained. Therefore, the material handling cost is the same for all algorithms. For this reason, we only compared the results in terms of the computation efficiency. To illustrate the efficiency of our proposed elimination algorithm, we have compared the average CPU TIME for computing a Newton step for each problem with two mathematical solvers (CPLEX, NEOS). We then have compared the average CPU TIME for finding an optimal arrangement under relative positioning constraints for each problem. The parameter settings and the computer specs used to code and run are summarized in Table 2.

4.2. Experimental Results

The results of the comparison of the average CPU TIME for computing a Newton step for each problem with or without eliminations are summarized in Table 3. The results of the comparison of the average CPU TIME for computing an optimal arrangement under relative-constraints are summarized in Table 4.
These examples show that, by exploiting the special structure of the problem, the proposed method reduces the CPU TIME to compute Newton Steps. In our implementation, it corresponds to around [0.1–0.2] µ sec. In particular, the proposed method is much faster for large sized problems. For a FLPRC with N = 35 , the proposed method solves the Newton equations in 0.2 µ sec, approximately 70 times faster than the generic solvers.
The sparsity pattern of the coefficient matrix of the Newton equation for SC35 is shown in Figure 2. In this example, the dimension of the coefficient matrix is n v + n λ = 1330 + 3955 = 5285 based on the (36) equation. However, due to the removal of some constraints from redundancy, the size of the matrix is reduced to n v + n λ = 1330 + 3360 = 4690 . We see that the lower right block is a diagonal matrix. Therefore, by using this structure, the main computational effort is reduced to the effort of solving Equation (38) whose coefficient matrix size is n v = 1330 . Therefore, it is possible to find the solution much more efficiently than solving the original Equation (36).
The sparsity pattern of the coefficient matrix of Equation (38) is shown in Figure 3. As is the same with Figure 2, the lower right block is diagonal, so using this structure can be more efficient than solving Equation (38).
Furthermore, the sparsity pattern of the coefficient matrix of the Equation (42) is as shown in Figure 4. This matrix is a dense matrix only for m = n = 1 , , N component and m = n = N + 1 , , 2 N component, and is a sparse matrix for the others. Therefore, the main computational load is consistent with solving this dense matrix.
Table 5 shows the calculation effort based on the calculation shown in Table 1. The required flops to compute the Newton steps are reduced a lot by the proposed block eliminations. This reduction leads to the faster CPU Time for computing Newton steps as shown in Table 3, and then leads to the faster CPU Time for finding an optimal solution of the FLPRC as shown in Table 4.
These results indicate that the proposed method reduction solves the problem faster, by exploiting the special structure of the matrix. Since the dimension of the final dense matrix matches the number of departments, no further decomposition is possible. Therefore, there is little room for further improvement in this problem.

4.3. Discussion

The method proposed in this study is the optimization of department position and shape given a relative position H i j , V i j . This technique can be used as a subroutine for a solution method using MIP. This approach also includes Liu and Meller [12] and Bozer and Wang [42], which combines MIP and MH. Among these algorithms, FLPRC is solved many times. Therefore, by using the proposed method, the calculation time can be shortened and more time can be devoted to the search of the relative positions.
The proposed method can also be used as a subroutine of the layout technique using FBS, such as [30] because FBS is used to derive the starting point in the proposed method. FBS has an encoding that divides the building into bays. Although this approach is simple, it does not necessarily guarantee optimality due to the layout that cannot be expressed by encoding. On the other hand, in the process of searching, if the proposed method is applied based on the relative position obtained by FBS, the layout of FBS can be improved, and there is a possibility that a better layout can be searched. In recent years, many MH approaches have used FBS and have achieved very good results. Further improvement can be considered by applying this technique.

5. Conclusions

We have considered the problem of computing an optimal arrangement of department within a facility, as measured by the minimum material handling cost, subject to the constraints for departments with respect to within-boundary, aspect-ratio, non-overlapping, and relative positioning. We have developed a primal-dual interior-point method that exploits the special structure of a problem, and is much faster than a standard method. We have proposed an efficient computation of the Newton step, using three types of block eliminations. To illustrate the effectiveness of our algorithm, we have given numerical results, using problem instances from several benchmark problems for FLP with specified relative-positioning given from randomized based rule. For a typical problem with 35 departments, the proposed algorithm is roughly 10 times faster. Using the proposed algorithm, we can accelerate a MIP (Mixed-Integer Programming) based algorithm for solving FLP, which is one of the most efficient methods for a continuous representation based problem.
Future research includes applying the proposed method to the MIP-based approach and MH-based approach. When applied to the MIP-based approach, the relaxed problem in the branch-and-bound method can be solved faster, which may solve a larger problem. The existing MIP-based approach can only solve problems for up to 12 departments. When applied to the MH-based approach, the proposed technique can be applied using the relative positional relationship of the layout obtained by FLEX-BAY or LOGIC. By doing so, the better layout can be obtained at each iteration, and thus there is possibility to get a better solution.

Author Contributions

Conceptualization, S.O. and K.Y.; methodology, S.O.; software, S.O.; validation, S.O.; writing—original draft preparation, S.O.; writing—review and editing, S.O. and K.Y.; supervision, K.Y. All authors have read and agreed to the published version of the manuscript.

Funding

This research received no external funding.

Conflicts of Interest

The authors declare no conflict of interest.

References

  1. Tompkins, J.A.; White, J.A.; Bozer, Y.A.; Tanchoco, J.M.A. Facilities Planning; John Wiley & Sons: Hoboken, NJ, USA, 2010. [Google Scholar]
  2. Richard, M.; Hales, L. Systematic Layout Planning. Available online: http://hpcinc.com/wp-content/uploads/2016/07/Systematic-Layout-Planning-SLP-4th-edition-soft-copy.pdf (accessed on 9 February 2021).
  3. Levary, R.R.; Kalchik, S. Facilities layout: A survey of solution procedures. Comput. Ind. Eng. 1985, 9, 141–148. [Google Scholar] [CrossRef]
  4. Kusiak, A.; Heragu, S.S. The facility layout problem. Eur. J. Oper. Res. 1987, 29, 229–251. [Google Scholar] [CrossRef]
  5. Hassan, M.M. Machine layout problem in modern manufacturing facilities. Int. J. Prod. Res. 1994, 32, 2559–2584. [Google Scholar] [CrossRef]
  6. Meller, R.D.; Gau, K.Y. The facility layout problem: Recent and emerging trends and perspectives. J. Manuf. Syst. 1996, 15, 351–366. [Google Scholar] [CrossRef]
  7. Drira, A.; Pierreval, H.; Hajri-Gabouj, S. Facility layout problems: A survey. Annu. Rev. Control 2007, 31, 255–267. [Google Scholar] [CrossRef]
  8. Arikaran, P.; Jayabalan, V.; Senthilkumar, R. Analysis of unequal areas facility layout problems. Int. J. Eng. 2010, 4, 44–51. [Google Scholar]
  9. Anjos, M.F.; Vieira, M.V. Mathematical optimization approaches for facility layout problems: The state-of-the-art and future research directions. Eur. J. Oper. Res. 2017, 261, 1–16. [Google Scholar] [CrossRef] [Green Version]
  10. Hosseini-Nasab, H.; Fereidouni, S.; Ghomi, S.M.T.F.; Fakhrzad, M.B. Classification of facility layout problems: A review study. Int. J. Adv. Manuf. Technol. 2018, 94, 957–977. [Google Scholar] [CrossRef]
  11. Koopmans, T.C.; Beckmann, M. Assignment Problems and the Location of Econmic Activities. Econometrica 1957, 25, 53–76. [Google Scholar] [CrossRef]
  12. Liu, Q.; Meller, R.D. A sequence-pair representation and MIP-model-based heuristic for the facility layout problem with rectangular departments. IIE Trans. 2007, 39, 377–394. [Google Scholar] [CrossRef]
  13. Van Camp, D.J.; Carter, M.W.; Vannelli, A. A nonlinear optimization approach for solving facility layout problems. Eur. J. Oper. Res. 1992, 57, 174–189. [Google Scholar] [CrossRef]
  14. Anjos, M.F.; Vannelli, A. An attractor-repeller approach to floorplanning. Math. Methods Oper. Res. 2002, 56, 3–27. [Google Scholar] [CrossRef]
  15. Anjos, M.F.; Vannelli, A. A new mathematical-programming framework for facility-layout design. INFORMS J. Comput. 2006, 18, 111–118. [Google Scholar] [CrossRef] [Green Version]
  16. Jankovits, I.; Luo, C.; Anjos, M.F.; Vannelli, A. A convex optimisation framework for the unequal-areas facility layout problem. Eur. J. Oper. Res. 2011, 214, 199–215. [Google Scholar] [CrossRef]
  17. Anjos, M.F.; Vieira, M.V. An improved two-stage optimization-based framework for unequal-areas facility layout. Optim. Lett. 2016, 10, 1379–1392. [Google Scholar] [CrossRef] [Green Version]
  18. Tam, K.Y. Genetic algorithms, function optimization, and facility layout design. Eur. J. Oper. Res. 1992, 63, 322–346. [Google Scholar]
  19. Gau, K.Y.; Meller, R.D. An iterative facility layout algorithm. International J. Prod. Res. 1999, 37, 3739–3758. [Google Scholar] [CrossRef]
  20. Shayan, E.; Chittilappilly, A. Genetic algorithm for facilities layout problems based on slicing tree structure. Int. J. Prod. Res. 2004, 42, 4055–4067. [Google Scholar] [CrossRef]
  21. Scholz, D.; Petrick, A.; Domschke, W. STaTS: A slicing tree and tabu search based heuristic for the unequal area facility layout problem. Eur. J. Oper. Res. 2009, 197, 166–178. [Google Scholar] [CrossRef]
  22. Kang, S.; Chae, J. Harmony search for the layout design of an unequal area facility. Expert Syst. Appl. 2017, 79, 269–281. [Google Scholar] [CrossRef]
  23. Tate, D.M.; Smith, A.E. Unequal-area facility layout by genetic search. IIE Trans. 1995, 27, 465–472. [Google Scholar] [CrossRef]
  24. Kulturel-Konak, S.; Smith, A.E.; Norman, B.A. Multi-objective tabu search using a multinomial probability mass function. Eur. J. Oper. Res. 2006, 169, 918–931. [Google Scholar] [CrossRef]
  25. Wong, K.Y. Solving facility layout problems using flexible bay structure representation and ant system algorithm. Expert Syst. Appl. 2010, 37, 5523–5527. [Google Scholar] [CrossRef]
  26. Kulturel-Konak, S.; Konak, A. A new relaxed flexible bay structure representation and particle swarm optimization for the unequal area facility layout problem. Eng. Optim. 2011, 43, 1263–1287. [Google Scholar] [CrossRef]
  27. Kulturel-Konak, S.; Konak, A. Unequal area flexible bay facility layout using ant colony optimisation. Int. J. Prod. Res. 2011, 49, 1877–1902. [Google Scholar] [CrossRef]
  28. Palomo-Romero, J.M.; Salas-Morera, L.; Garcia-Hernandez, L. An island model genetic algorithm for unequal area facility layout problems. Expert Syst. Appl. 2017, 68, 151–162. [Google Scholar] [CrossRef]
  29. Garcia-Hernandez, L.; Salas-Morera, L.; Garcia-Hernandez, J.A.; Salcedo-Sanz, S.; de Oliveira, J.V. Applying the coral reefs optimization algorithm for solving unequal area facility layout problems. Expert Syst. Appl. 2019, 138, 112819. [Google Scholar] [CrossRef]
  30. Garcia-Hernandez, L.; Salas-Morera, L.; Carmona-Munoz, C.; Garcia-Hernandez, J.A.; Salcedo-Sanz, S. A novel Island Model based on Coral Reefs Optimization algorithm for solving the unequal area facility layout problem. Eng. Appl. Artif. Intell. 2020, 89, 103445. [Google Scholar] [CrossRef]
  31. Liu, J.; Zhang, H.; He, K.; Jiang, S. Multi-objective particle swarm optimization algorithm based on objective space division for the unequal-area facility layout problem. Expert Syst. Appl. 2018, 102, 179–192. [Google Scholar] [CrossRef]
  32. Guan, C.; Zhang, Z.; Liu, S.; Gong, J. Multi-objective particle swarm optimization for multi-workshop facility layout problem. J. Manuf. Syst. 2019, 53, 32–48. [Google Scholar] [CrossRef]
  33. Pourvaziri, H.; Pierreval, H. Dynamic facility layout problem based on open queuing network theory. Eur. J. Oper. Res. 2017, 259, 538–553. [Google Scholar] [CrossRef]
  34. Turanoglu, B.; Akkaya, G. A new hybrid heuristic algorithm based on bacterial foraging optimization for the dynamic facility layout problem. Expert Syst. Appl. 2018, 98, 93–104. [Google Scholar] [CrossRef]
  35. Tayal, A.; Gunasekaran, A.; Singh, S.P.; Dubey, R.; Papadopoulos, T. Formulating and solving sustainable stochastic dynamic facility layout problem: A key to sustainable operations. Ann. Oper. Res. 2017, 253, 621–655. [Google Scholar] [CrossRef] [Green Version]
  36. Azevedo, M.M.; Crispim, J.A.; de Sousa, J.P. A dynamic multi-objective approach for the reconfigurable multi-facility layout problem. J. Manuf. Syst. 2017, 42, 140–152. [Google Scholar] [CrossRef]
  37. Pourhassan, M.R.; Raissi, S. An integrated simulation-based optimization technique for multi-objective dynamic facility layout problem. J. Ind. Inf. Integr. 2017, 8, 49–58. [Google Scholar] [CrossRef]
  38. Montreuil, B. A modelling framework for integrating layout design and flow network design. In Material Handling’ 90; Springer: Berlin/Heidelberg, Germany, 1991; pp. 95–115. [Google Scholar]
  39. Meller, R.D.; Narayanan, V.; Vance, P.H. Optimal facility layout design. Oper. Res. Lett. 1998, 23, 117–127. [Google Scholar] [CrossRef]
  40. Sherali, H.D.; Fraticelli, B.M.; Meller, R.D. Enhanced model formulations for optimal facility layout. Oper. Res. 2003, 51, 629–644. [Google Scholar] [CrossRef]
  41. Castillo, I.; Westerlund, T. An ε-accurate model for optimal unequal-area block layout design. Comput. Oper. Res. 2005, 32, 429–447. [Google Scholar] [CrossRef]
  42. Bozer, Y.A.; Wang, C.T. A graph-pair representation and MIP-model-based heuristic for the unequal-area facility layout problem. Eur. J. Oper. Res. 2012, 218, 382–391. [Google Scholar] [CrossRef]
  43. Chae, J.; Regan, A.C. Layout design problems with heterogeneous area constraints. Comput. Ind. Eng. 2016, 102, 198–207. [Google Scholar] [CrossRef]
  44. Boyd, S.; Vandenberghe, L. Convex Optimization; Cambridge University Press: Cambridge, UK, 2004. [Google Scholar]
Figure 1. Example illustrating the bay assignment (shown at center) and the obtained layout (shown at right) from the horizontal and vertical relative-positioning graphs (shown at left).
Figure 1. Example illustrating the bay assignment (shown at center) and the obtained layout (shown at right) from the horizontal and vertical relative-positioning graphs (shown at left).
Algorithms 14 00060 g001
Figure 2. Sparsity pattern of the coefficient matrix of the Newton equation for SC35.
Figure 2. Sparsity pattern of the coefficient matrix of the Newton equation for SC35.
Algorithms 14 00060 g002
Figure 3. Sparsity pattern of the coefficient matrix after the 1st elimination for SC35.
Figure 3. Sparsity pattern of the coefficient matrix after the 1st elimination for SC35.
Algorithms 14 00060 g003
Figure 4. Sparsity pattern of the coefficient matrix after the 2nd elimination for SC35.
Figure 4. Sparsity pattern of the coefficient matrix after the 2nd elimination for SC35.
Algorithms 14 00060 g004
Table 1. Summary of complexity to solve Newton equation.
Table 1. Summary of complexity to solve Newton equation.
EquationComplexity (# of Flops)
Original Newton Equation (36) ( 2 / 3 ) ( n v + n λ ) 3 343 8 N 6 + ( 2793 8 + 147 4 M ) N 5
+ ( 7581 8 + 399 2 M + 21 2 M 2 ) N 4
+ ( 6859 3 + 1083 4 M + 57 2 M 2 + M 3 ) N 3
After 1st elimination (38) 2 n v 2 n λ + ( 2 / 3 ) n v 3 17 3 N 6 + ( 47 + 2 M ) N 5
+ ( 129 + 12 M ) N 4
+ ( 117 + 18 M ) N 3
After 2nd elimination (42) 2 n v 2 n λ + 2 n r 2 n d 4 5 N 6 + ( 79 5 + 2 M ) N 5
+ ( 2 / 3 ) n r 3 + ( 526 5 + 12 M ) N 4
+ ( 329 3 + 18 M ) N 3
After 3rd elimination (47) 2 n v 2 n λ + 2 n r 2 n d 4 5 N 6 + ( 79 5 + 2 M ) N 5
+ n w h + 2 n x y 2 n w h + ( 526 5 + 12 M ) N 4
+ ( 4 / 3 ) N 3 + ( 52 3 + 18 M ) N 3 + 2 N
Table 2. The outline of experimental conditions.
Table 2. The outline of experimental conditions.
OSMac OS X (ver. 10.7.3)
Processor3.4GHz Intel Core i7
Memory8GB 1333MHz DDR3
LanguageMATLAB 2011Ra
α 0.5
β 0.5
t 0 10 4
ϵ 10 3
Table 3. Average CPU TIME for computing Newton Steps ( µ sec).
Table 3. Average CPU TIME for computing Newton Steps ( µ sec).
Problemwith Eliminationswithout Eliminations
O70.14.2
O80.14.5
O90.14.7
VC100.24.8
Ba100.25.3
M15a0.27.8
M15s0.28.7
AB200.28.9
Tam300.210.1
SC300.212.6
SC350.213.9
Table 4. Average CPU TIME for computing an optimal arrangement under relative-constraints ( µ  sec).
Table 4. Average CPU TIME for computing an optimal arrangement under relative-constraints ( µ  sec).
ProblemProposalCPLEXNEOS
O7141415
O8141322
O9142222
VC10153459
Ba10153789
M15a154794
M15s1566136
AB201577162
Tam301687207
SC3016102225
SC3516117320
Table 5. Summary of complexity to solve Newton equation for SC35.
Table 5. Summary of complexity to solve Newton equation for SC35.
EquationComplexity (# of Flops)
Original Newton Equation (36)68,774,472,667
After 1st elimination (38)133,541,333
After 2nd elimination (42)48,477,333
After 3rd elimination (47)971,903
Publisher’s Note: MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Share and Cite

MDPI and ACS Style

Ohmori, S.; Yoshimoto, K. A Primal-Dual Interior-Point Method for Facility Layout Problem with Relative-Positioning Constraints. Algorithms 2021, 14, 60. https://0-doi-org.brum.beds.ac.uk/10.3390/a14020060

AMA Style

Ohmori S, Yoshimoto K. A Primal-Dual Interior-Point Method for Facility Layout Problem with Relative-Positioning Constraints. Algorithms. 2021; 14(2):60. https://0-doi-org.brum.beds.ac.uk/10.3390/a14020060

Chicago/Turabian Style

Ohmori, Shunichi, and Kazuho Yoshimoto. 2021. "A Primal-Dual Interior-Point Method for Facility Layout Problem with Relative-Positioning Constraints" Algorithms 14, no. 2: 60. https://0-doi-org.brum.beds.ac.uk/10.3390/a14020060

Note that from the first issue of 2016, this journal uses article numbers instead of page numbers. See further details here.

Article Metrics

Back to TopTop