Energy-Saving Oriented Manufacturing Workshop Facility Layout: A Solution Approach Using Multi-Objective Particle Swarm Optimization

Abstract

Low-carbon sustainable development has become the consensus of manufacturing enterprises to fulfill their social responsibilities. Facility layout is an essential part of manufacturing system planning. Current research has demonstrated the advantages of energy saving on the manufacturing system level where operational methods (e.g., energy-efficient production scheduling and path planning) can be utilized and do not require massive investment in the existing legacy system. However, these efforts are mostly based on the existing fixed facility layout. Meanwhile, although facility layout problems have been extensively studied so far, the related work seldom involves the optimization of energy consumption (EC) or other EC-related environmental impact indicators, and does not clearly reveal if EC can be an independent optimization objective in facility layout. Accordingly, whether the energy-saving potential of a manufacturing system can be further tapped through rational facility layout is the gap of the current study. To address this, an investigation into energy-saving oriented manufacturing workshop facility layout is conducted. Correspondingly, an energy-efficient facility layout (EFL) model for the multi-objective optimization problem that minimizes total load transport distance and EC is formulated, and a multi-objective particle swarm optimization-based method is proposed as the solution. Furthermore, experimental studies verify the effectiveness of the presented model and its solution, indicating that EC can be regarded as an independent optimization objective during facility layout, and EFL is a feasible energy-saving approach for a manufacturing system.

1. Introduction

Climate warming is a global environmental problem, and the main reason is the increasing concentration of greenhouse gases in the atmosphere. As an active response to the global climate change crisis, achieving carbon peak and carbon neutralization is a significant development strategy for China and the world’s major economies. The manufacturing industry is the pillar of social and economic development, which is also the major source of resource consumption and waste generation. In China, carbon dioxide emissions mainly stem from energy activities, and statistics indicate that the energy consumed by the manufacturing industry takes up nearly 61.6% of total national energy consumption (EC) in 2019 [1]. Hence, the manufacturing industry has a responsibility to reduce greenhouse gas emissions, and the concept of low-carbon sustainable development has become its consensus.

Meanwhile, with the growth of consumers’ individualized demands and the intensification of market competition, traditional large-scale and single-variety manufacturing modes are being gradually abandoned and flexible manufacturing modes are being adopted by mechanical manufacturing enterprises. Correspondingly, facility layout is an important aspect relevant to flexible manufacturing. In a manufacturing workshop environment, a facility layout problem (FLP) generally involves the rational arrangement of a given number of manufacturing equipment/cells within the given space constraints so that some performance goals of material handling systems can be achieved [2]. However, the objectives mostly concerned about in the current FLP study are transport distance, material handling cost (MHC), effectively-utilized area, etc. [3,4], and energy consumption (EC) or other EC-related environmental impact indicators are seldom touched on. It has been proved by numerous existing research that energy-efficient scheduling is a crucial approach to minimize the EC of machining systems with various forms (e.g., flow shop, job shop, and flexible job shop) [5,6,7], and energy-efficient path planning is beneficial to reducing transport EC [8]. As a typical manufacturing system, a machining system (MS) usually consists of many different elements, such as machine tools, automated material handling systems, and tool storage [9]. The locations of the machines in an MS are directly related to the workpiece transport distances and modes among different equipment, affecting the assignment of jobs and the arrangement of operation sequences in each available machine in the subsequent production scheduling stage. Therefore, the facility layout is also a critical link for manufacturing enterprises to realize low-carbon manufacturing, and there are many answers yet to be sought concerning how to locate manufacturing facilities in an energy-saving/sustainable way. In view of this, the FLP is investigated from the perspective of energy saving in this study, and the corresponding model with transport distance and EC as the optimization objectives, namely the energy-efficient facility layout (EFL) model, is formulated.

The rest of this paper is organized as follows: a review of FLP and its influence on EC in Section 2 is followed by the establishment of the EFL model in Section 3. Then, a solution method based on a multi-objective particle swarm optimization (MOPSO) algorithm is illustrated in Section 4, and Section 5 is devoted to case studies to verify the effectiveness of the presented model and the algorithm. Finally, conclusions are drawn and future work is prospected in Section 6.

2. Literature Review

FLPs are recognized as a crucial class of operations research problems, and numerous efforts have been devoted to them from various perspectives since the 1960s. Generally, facility layouts are directly affected by the specifications of MSs, such as product variety and production volumes. In practice, four common layout forms, namely process layout, product layout, fixed-position layout, and cellular layout [10], can be classified. From the viewpoint of layout configuration, FLPs can be divided into single-row, double-row, multi-row, loop, and multi-floor FLPs [11,12,13]. Meanwhile, based on the layout flexibility and changes over time, FLPs can be categorized into static and dynamic FLPs. A dynamic FLP can be divided into several static FLPs over time, which involves the decisions on whether the layout needs to be modified to another in each period [14]. So far, extensive achievements have been reported on the above layout types in the literature.

To better describe and investigate the FLPs under different conditions, various models have been proposed. The commonly encountered FLP models in the literature include quadratic assignment problem (QAP) models, quadratic set covering problem (QSCP) models, mixed-integer programming (MIP) models, and graph theory-based models [4,10]. Specifically, an FLP with a discrete representation and equal-area facilities tends to be formulated as a QAP model. In a QAP formulation, a plant site is usually divided into rectangular blocks with the same area and shape, and each facility is assigned to exactly one block by disregarding its actual shape and size [15]. However, when the number of facilities is less than that of blocks, the facilities need to be transferred to a larger layout space to guarantee the applicability of the QAP model. The QSCP model is also suitable to formulate discrete layout optimization, which is similar to the QAP model as the total area occupied by all facilities is divided into plenty of small blocks. However, in a QSCP model, the possible area and shape of each facility can be expressed in terms of the occupied blocks [16]. Regarding the MIP model, it consists of the objective functions expressed by mixed integer and non-integer variables, and constraints including a series of equations and inequalities. The MIP model is suitable to formulate an FLP on a continuous surface, leading to widespread utilization in formulating single/double/multi-row FLP, unequal-area FLP, and multi-floor FLP [4,13]. Correspondingly, in an FLP with a continuous representation, the exact positions, orientations, and/or pick-up/drop-off points of facilities can be expressed and acquired by solving the specific model. Furthermore, an FLP can be modeled as a graph in which the facilities are represented by the nodes, and the adjacency relationship between two facilities is represented by an arc connecting two corresponding nodes according to the process routes of related products [17,18]. By analyzing the material flow and activity relationship between different facilities, the graph-theoretic formulation can help generate a spatial relationship diagram to determine the adjacency relationship between facilities, but the final layout may be irregular as the sizes and shapes of facilities are not considered.

Regarding the facility layout approaches, in the early stage, the facility layout mainly depended on the subjective experience and intuitive feeling of planners, lacking scientific and theoretical guidance. With the rapid development of computer and information technology, facility layout methods have developed from the early manual and qualitative methods to the current quantitative, automatic, and intelligent ones. Currently, the facility layout methods reported in the literature can be broadly categorized into two classes: factor assessment approach and mathematical/algorithmic approach [19]. The former can be treated as a decision-making process and usually divides the facility layout process into several fixed and sequential steps. The most representative one is the systematic layout planning (SLP) method [20]. This method comprehensively considers the qualitative and quantitative factors in the facility layout analysis. When the number of facilities involved is small, it can help accurately obtain the adjacency relationship and material flow between facilities and then generate the conceptual layout scheme. Nevertheless, it cannot provide the specific location of each facility. So, its application has certain limitations when the FLP scale is large and the constraints are complex.

In contrast, mathematical/algorithmic approaches are closely related to FLP models, which can be further categorized into exact and approximate approaches. The former aims to seek optimal solutions, and the representative approaches include dynamic programming, semi-definite programming, branch-and-bound, and cutting-plane algorithm [10,16,21,22]. Despite certain advantages in solving small-scale FLPs, the computational effort required tends to increase sharply with the increase in the FLP scale. For example, Kettani and Oral [23] pointed out the number of possible facilities should not be more than 15 to facilitate solving an MIP model with commercial codes. Approximate approaches tend to provide approximate optimal solutions in an acceptable computation time, which can be broadly classified into traditional heuristic methods and meta-heuristic methods. Traditional heuristic methods belong to the local search algorithm and utilize heuristic information and some experience/rules to search solutions along the most feasible direction. Generally, they can be divided into four categories [24]: construction methods (e.g., CORELAP [25]), improvement methods (e.g., CRAFT [26], and MULTIPLE [27]), hybrid methods (e.g., FLAC [28], and DISCON [29]), and graph-theoretic methods. Note that they are still more suitable for solving small-scale FLPs. When the problem scale increases, the solution search is prone to fall into the local optimum. Besides, different heuristic rules often need to be designed for different FLPs, resulting in poor generality. With the rapid development of computing power in the last decade, meta-heuristic methods have been more popular to solve FLPs, especially large-scale FLPs. The typical meta-heuristic methods include genetic algorithm (GA) [3], particle swarm optimization (PSO) [13], simulated annealing (SA) [30], tabu search [31], and ant colony optimization (ACO) [32]. These methods usually allow the emergence of inferior solutions in the search process, and their global search abilities can be enhanced by introducing the mechanism of jumping out the local optimum. Thus, it is likely to obtain global optimization solutions. Moreover, when applying meta-heuristic methods, the algorithm efficiency relies on the individual encoding scheme reflecting feasible solutions to a great extent.

In addition, from the perspective of objectives, a layout can be evaluated qualitatively and quantitatively. Qualitative objectives usually refer to the safety, noise, color, or cleanliness of a layout [33]. Meanwhile, the most significant quantitative indicator is MHC [10], and transport distance, transport time, and material flow are all common MHC-related indicators. Besides, layout area utilization and total production time are some other performance indicators [3]. With the increasing concern on sustainability in the manufacturing industry in recent years, some scholars have begun to study sustainable FLPs. Tayal et al. [34] formulated a sustainable stochastic dynamic FLP and presented a solution method with the hierarchical framework of a meta-heuristic, multiple attribute decision-making techniques, and consensus ranking method. Macroscopically, the sustainability of manufacturing processes involves economic, environmental, and social issues [35]. Accordingly, apart from the economic sustainability represented by manufacturing cost/MHC, attributes such as waste disposal, recycling, noise, maintenance, the safety of human–machine interaction, and EC should also be addressed in sustainable facility layout. In terms of EFL-related studies, Iqbal and Al-Ghamdi [36] took the machining and transport EC as a whole to optimize and presented an energy-saving method by rationally assigning manufacturing processes to various machines and grouping machines in various cells. Specifically, the transport EC was simplified as a function dependent on the part’s weight. Moreover, Wang et al. [37] proposed a facility layout model to minimize the total cost in a large-scale industrial plant with multi-floor structures, and the energy consumed by pumps for overcoming friction and gravity losses in the piping was treated as part of the total cost. Similarly, Lamba et al. [38] formulated a dynamic cellular FLP as a mixed-integer non-linear model, and the sustainability was incorporated by optimizing EC, material handling, and rearrangement cost. Overall, current FLP studies involving EC or other EC-related indicators are rather limited. The related efforts usually treated the EC of material handling equipment as part of transport cost and further considered that reducing the cost was largely equivalent to saving energy. That is, whether EC can be an independent optimization objective in facility layout was not fully answered. Besides, the evaluation of transport EC was also greatly simplified due to the lack of reliable models and methods. A manufacturing system usually consists of multiple energy-consuming equipment, and material handling is an indispensable link in the whole production process. So, transport EC is an essential part of the total production EC. Correspondingly, the workshop facility layout directly affects the transport mode, equipment selection, and route planning, thereby the transport EC and the energy-saving potential of scheduling schemes. Given this, it is of great significance to further study EFL.

3. Problem Description and Energy-Efficient Facility Layout Modeling

To bridge the gap summarized by the above literature review, we conducted a facility layout study considering EC optimization in a manufacturing workshop environment, and an EFL model is formulated in this section. With the increasingly fierce market competition and more diversified consumer demand, multi-variety, small-batch, and flexible production mode, which can quickly respond to market changes and update products, has become a trend for manufacturing enterprises to improve productivity and competitiveness. Therefore, this production mode was our research background. From the perspective of layout types, the cellular layout, with the advantages of the high efficiency of product layout (i.e., flow line) and the flexibility of process layout, was our research focus. Besides, automated guided vehicles (AGVs) have been widely applied in a flexible manufacturing environment, so the material flow between different facilities in a manufacturing cell was assumed to be executed by AGVs. Moreover, when an AGV is assigned material handling tasks, empty transport is often inevitable. Ideally, the empty transport EC should also be taken into account when evaluating the total transport EC. However, owing to the limited AGV resources in a manufacturing workshop, the empty transport EC of each AGV mainly depends on the execution order of the transport tasks assigned to it, and the specific planned path [8], which are the focus of the transport task scheduling and path planning research, respectively. Therefore, we only considered load transport EC in this study. Then, the overall problem is to determine the optimal facility layout scheme such that the load transport distance and EC can both be minimized.

3.1. Model Hypotheses

Based on the existing FLP modeling research, the following assumptions are made to formulate the EFL model:

• The workshop for facility layout has a rectangular/square shape, and its length and width are known in advance.

• The shape of each facility is abstracted as the smallest rectangle/square enveloping the real physical equipment operation area, and its length and width are known and fixed. Moreover, each facility has a safety clearance space.

• Each facility owns a pick-up point and a drop-off point, which are located at the facility boundary and maybe in the same position or different positions, and the positions of such two points relative to the facility centroid are fixed, respectively.

• All facilities are arranged in the same plane, and the rectilinear distance from the pick-up point of one facility to the drop-off point of another facility is utilized to evaluate the transport distance between them.

• All facilities can be arranged freely but must be located in the given workshop layout area. Besides, facilities are not allowed to overlap with each other.

• Each facility has free orientation but can only be placed horizontally (the longer side at the bottom) or vertically (the shorter side at the bottom).

• The shop floor is flat, and the types of AGVs executing all material handling tasks are the same. All AGVs are available when assigned transport tasks. Once an AGV starts the load transport, the material handling process will not be interrupted, and the AGV moves at a constant speed. Moreover, the possible route conflicts among AGVs are ignored.

• The mass, output volume, transport batch size, and process route of each product are known in advance, and the product mass change due to machining is not considered.

3.2. Mathematical Formulation

Based on the research hypotheses, EFL is a static FLP on a two-dimensional continuous plane. To accurately represent the location and orientation of each facility, a reference coordinate system, which takes the bottom-left corner as the coordinate origin and the bottom and left boundaries of the rectangular layout area as the x-axis and y-axis, respectively, is established firstly. The diagram of EFL in the workshop layout area is shown in Figure 1.

The EFL is for the production of n products and involves m facilities in all. Since each facility is abstracted as a regular rectangle, its center can be viewed as the centroid. Let the longer side of a rectangular layout area be the length direction and the shorter side be the width direction, the length and width of a workshop layout area are denoted as L and W, respectively. Then, the definitions of the length and width of each facility follows the length and width directions of the layout area, respectively. As shown in Figure 1, the center, pick-up point, and drop-off point of facility i ($i=1, 2, \dots , m$) are expressed as $(x_i, y_i)$, $(x_i^{\mathrm{out}}, y_i^{\mathrm{out}})$, and $(x_i^{\mathrm{in}}, y_i^{\mathrm{in}})$, separately. To guarantee the safety of production and transport processes, ${\Delta }_i$ is denoted as the safety clearance distance of facility i, and no other facilities are allowed to be placed in its safety clearance area. Furthermore, the length and width of a facility may be exchanged with the change in its orientation. Correspondingly, ${\theta }_i$ is denoted to express the orientation of facility i, and its initial value is 0. The length and width of facility i in the original orientation (i.e., ${\theta }_i=0$) are denoted as $l_{0i}$ and $w_{0i}$, respectively, and $l_{0i}\ge w_{0i}$. The length and width of facility i in the actual layout are denoted as $l_i$ and $w_i$, respectively. Obviously, $l_i$ is equal to $l_{0i}$ or $w_{0i}$, which depends on the specific orientation represented by ${\theta }_i$. Based on the research assumption, ${\theta }_i$ can be set as one element of {0, π/2, π, 3π/2} in the actual layout. Meanwhile, for each facility, since the relative positional relationship between its pick-up (drop-off) point and its center is determined in advance, the pick-up (drop-off) point coordinates depend on the center coordinates and orientation. Specifically, when the center coordinates are determined, the influence of orientation on the coordinates of the pick-up and drop-off points is shown in Figure 2.

It can be observed that the original length and width of facility i exchange only when ${\theta }_i$ is equal to π/2 or 3π/2. Therefore, $r_i$ is defined to indicate whether the original length and width of facility i exchange, and the relationship between $r_i$ and ${\theta }_i$ can be expressed as:

(1)$r_i=\left|\cos {\theta }_i\right|, \forall i$

Accordingly, the value of $r_i$ belongs to {0, 1}. $l_i$ and $w_i$ can be obtained as:

(2)$\{\begin{array}{l}{l}_i=r_i{l}_{0i}+(1-r_i)w_{0i}\\ w_i=r_i{w}_{0i}+(1-r_i)l_{0i}\end{array}, \forall i$

For each facility i in its original orientation, if a plane rectangular coordinate system, with its center as the coordinate origin, the length direction as the x-axis direction, and the width direction as the y-axis direction, is established, the pick-up and drop-off point coordinates can be denoted as $({\widehat{x}}_i^{\mathrm{out}}, {\widehat{y}}_i^{\mathrm{out}})$ and $({\widehat{x}}_i^{\mathrm{in}}, {\widehat{y}}_i^{\mathrm{in}})$ in this coordinate system, respectively. According to the modeling assumption, $({\widehat{x}}_i^{\mathrm{out}}, {\widehat{y}}_i^{\mathrm{out}})$ and $({\widehat{x}}_i^{\mathrm{in}}, {\widehat{y}}_i^{\mathrm{in}})$ are both determined in advance. Then, through coordinate system transformation, $(x_i^{\mathrm{out}}, y_i^{\mathrm{out}})$ and $(x_i^{\mathrm{in}}, y_i^{\mathrm{in}})$ can be acquired as:

(3)$\{\begin{array}{l}{x}_i^{\mathrm{in}}=x_i+{\widehat{x}}_i^{\mathrm{in}}\cos {\theta }_i-{\widehat{y}}_i^{\mathrm{in}}\sin {\theta }_i\\ y_i^{\mathrm{in}}=y_i+{\widehat{y}}_i^{\mathrm{in}}\cos {\theta }_i+{\widehat{x}}_i^{\mathrm{in}}\sin {\theta }_i\end{array}, \forall i$

(4)$\{\begin{array}{l}{x}_i^{\mathrm{out}}=x_i+{\widehat{x}}_i^{\mathrm{out}}\cos {\theta }_i-{\widehat{y}}_i^{\mathrm{out}}\sin {\theta }_i\\ y_i^{\mathrm{out}}=y_i+{\widehat{y}}_i^{\mathrm{out}}\cos {\theta }_i+{\widehat{x}}_i^{\mathrm{out}}\sin {\theta }_i\end{array}, \forall i$

In addition, the total output of product z ($z=1, 2, \dots , n$) is defined as $O_z$, and the average transport batch size in the production process is $B_z$, pieces/time. So, the number of transport batches needed for product z ($N_z$) can be obtained as:

(5)$N_z=\mathrm{ceil}\left(O_z/B_z\right),\forall z$

where ceil is a single variable function used to obtain the smallest integer that is not smaller than the input variable. Obviously, for product z, the number of products in its last batch is not more than $B_z$. Besides, $X_{ijz}$ is an integer indicating the total transfer times from facility i to facility j ($j=1, 2, \dots , m$) in the process route of product z. Owing to the given process routes of relevant products in advance, the value of $X_{ijz}$ can be determined, and $X_{ijz}\ge 0$. Moreover, the mass of product z and the empty AGV are denoted as $m_z$ and $m_0$, respectively. The AGV transport speed is uniformly defined as $v_a$.

3.2.1. Transport Distance Analysis

The load transport distance from facility i to facility j ($d_{ij}$, m) relies on the specific facility layout and can be presented as:

(6)$d_{ij}=\left|x_i^{\mathrm{out}}-x_j^{\mathrm{in}}\right|+\left|y_i^{\mathrm{out}}-y_j^{\mathrm{in}}\right|, \forall i,j$

The frequency of material flow from facility i to facility j ($f_{ij}$) can be calculated as:

(7)$f_{ij}={\displaystyle \sum _{i=1}^m{\displaystyle \sum _{j=1}^m{\displaystyle \sum _{z=1}^n\left(X_{ijz}{N}_z\right)}}}$

Then, the total load transport distance ($D_{\mathrm{total}}$, m) can be formulated as:

(8)$D_{\mathrm{total}}={\displaystyle \sum _{i=1}^m{\displaystyle \sum _{j=1}^m\left(f_{ij}{d}_{ij}\right)}}$

3.2.2. Transport Energy Consumption Analysis

Generally, four kinds of resistance, namely friction, air, slope, and acceleration resistance, need to be overcome to drive an AGV [39]. According to the assumptions of workshop ground and AGV transport speed, slope and acceleration resistance can be ignored. Besides, considering the indoor transport environment and the relatively slow AGV speed, air resistance can also be neglected. Then, from the perspective of motion [8], the AGV transport EC in the load transport process is composed of the standby and uniform motion EC. As long as an AGV is powered on, the standby motion always exists regardless of whether the AGV executes a transport task. Correspondingly, the standby motion power ($P_{\mathrm{sm}}$, W) is usually treated as a constant, which can be measured when the AGV is in standby mode or acquired by referring to the manual. The common energy sources driving an AGV to move at various speeds are motors, and the uniform motion EC can be viewed as the additional EC for igniting an AGV to move at a constant speed on top of the standby motion EC. Therefore, the total load transport EC ($E_{\mathrm{total}}$, J) can be expressed as:

(9)$E_{\mathrm{total}}=E_{\mathrm{sm}}+E_{\mathrm{um}}$

where $E_{\mathrm{sm}}$ and $E_{\mathrm{um}}$ represent the total EC ignited by standby and uniform motions, respectively, J. Owing to the fixed standby motion power, $E_{\mathrm{sm}}$ relies on the load transport time. Based on Equations (6) and (7), the total load transport time from facility i to facility j ($T_{ij}$, s) can be calculated as:

(10)$T_{ij}=\frac{f_{ij}{d}_{ij}}{v_a}, \forall i, j$

Accordingly, $E_{\mathrm{sm}}$ can be acquired as:

(11)$E_{\mathrm{sm}}={\displaystyle \sum _{i=1}^m{\displaystyle \sum _{j=1}^m\left(P_{\mathrm{sm}}{T}_{ij}\right)}}$

In the uniform motion process, the mechanical power output from the AGV’s traveling driving motors is utilized to work against the rolling resistance, which can be calculated as the product of driving force and displacement. The driving force is equal to the rolling resistance, and $E_{\mathrm{um}}$ can be expressed as:

(12)$E_{\mathrm{um}}={\displaystyle \sum _{i=1}^m{\displaystyle \sum _{j=1}^m{\displaystyle \sum _{z=1}^n\left[\frac{1}{\eta }{C}_{\mathrm{r}}(m_z{O}_z+N_z{m}_0)g{d}_{ij}{X}_{ijz}\right]}}}$

where $\eta$ is the overall power factor of driving motors; $C_{\mathrm{r}}$ is the rolling resistance coefficient; and g is the gravity acceleration, 9.81 m/s².

Then, the optimization objectives of EFL are formulated as:

(13)$\{\begin{array}{l}\mathrm{Minimize}[D_{\mathrm{total}}]\\ \mathrm{Minimize}[E_{\mathrm{total}}]\end{array}$

Subject to:

(14)$\frac{1}{2}{l}_i+{\Delta }_i\le x_i\le L-\frac{1}{2}{l}_i-{\Delta }_i, \forall i$

(15)$\frac{1}{2}{w}_i+{\Delta }_i\le y_i\le W-\frac{1}{2}{w}_i-{\Delta }_i, \forall i$

(16)$\left|x_i-x_j\right|\ge \frac{1}{2}(l_i+l_j)+\max \{{\Delta }_i, {\Delta }_j\}, \forall i, j$

(17)$\left|y_i-y_j\right|\ge \frac{1}{2}(w_i+w_j)+\max \{{\Delta }_i, {\Delta }_j\}, \forall i, j$

(18)${\theta }_i\in \{0, \mathsf{\pi }/2, \mathsf{\pi }, 3\mathsf{\pi }/2\}$

Combined with Equations (1) and (2), constraints (14) and (15) ensure that all facilities must be in the assigned layout area. Constraints (16) and (17) confirm that facilities do not overlap in the layout, and a reasonable safety distance should be maintained between any two facilities. Constraint (18) indicates the orientation constraint of a facility to further determine its pick-up and drop-off point coordinates.

4. Model Solution

As shown in the above mathematical formulation, the EFL is also a multi-objective optimization problem (MOP), and $x_i$, $y_i$, and ${\theta }_i$ are key decision variables affecting the optimization objectives. For an MOP, especially when the objectives conflict, there is often no optimal solution that makes all objective functions reach the maximum/minimum value at the same time, and we usually seek its non-inferior solutions (known as non-dominated or Pareto solutions) within the limited computing time. According to the literature review of the FLP solutions, meta-heuristic methods such as ACO, GA, PSO, and SA are suitable for solving large-scale FLPs and own advantages such as high search efficiency, self-adaptation, and self-learning. Because of this, PSO is adopted to solve the EFL model.

4.1. Particle Design

In PSO, the position of each particle represents a feasible problem solution, and all particles move through the search space to find the optimal solution. Since the algorithm design is driven by problem characteristics and the algorithm efficiency relies on the particle encoding scheme, it is the first step to encode a solution of the EFL model into a particle when applying PSO. The value of ${\theta }_i$ is discrete, while $x_i$ and $y_i$ change in continuous space. Therefore, based on the indirect encoding method to solve unequal-area FLPs [40], the values of $x_i$, $y_i$, and ${\theta }_i$ can be determined indirectly by ${\overline{x}}_i$, ${\overline{y}}_i$, and ${\overline{\theta }}_i$, which all belong to [0, 1] and change continuously in this interval. The specific conversion relationship can be expressed as:

(19)$\{\begin{array}{l}{x}_i=x_i^{\min }+{\overline{x}}_i(x_i^{\max }-x_i^{\min }),\forall i\\ y_i=y_i^{\min }+{\overline{y}}_i(y_i^{\max }-y_i^{\min }),\forall i\\ {\theta }_i=\{\begin{array}{c}0\\ \mathsf{\pi }/2\\ \mathsf{\pi }\\ 3\mathsf{\pi }/2\end{array} \begin{array}{c}0\le {\overline{\theta }}_i<0.25\\ 0.25\le {\overline{\theta }}_i<0.5\\ 0.5\le {\overline{\theta }}_i<0.75\\ 0.75\le {\overline{\theta }}_i\le 1\end{array} ,\forall i\end{array}$

where $x_i^{\min }$, $x_i^{\max }$,$y_i^{\min }$, and $y_i^{\max }$ can be determined by constraints (14) and (15), as shown below:

(20)$\{\begin{array}{l}{x}_i^{\min }=\frac{1}{2}{l}_i+{\Delta }_i, \forall i\\ x_i^{\max }=L-\frac{1}{2}{l}_i-{\Delta }_i, \forall i\\ y_i^{\min }=\frac{1}{2}{w}_i+{\Delta }_i, \forall i\\ y_i^{\max }=W-\frac{1}{2}{w}_i-{\Delta }_i, \forall i\end{array}$

Specifically, $l_i$ and $w_i$ rely on the facility orientation and can be calculated by Equations (1) and (2).

Then, each particle can be coded in three layers, and the encoding length of each layer is 2m. Correspondingly, the search space consists of three m-dimensional spaces, and the particle position in each m-dimensional space determines the x coordinate, y coordinate, and orientation of each facility, respectively. The three-layer encoding of each particle can be respectively expressed as the form $({\overline{x}}_1, {\overline{x}}_2, \dots , {\overline{x}}_m,{\overline{vx}}_1, {\overline{vx}}_2, \dots , {\overline{vx}}_m)$, $({\overline{y}}_1, {\overline{y}}_2, \dots , {\overline{y}}_m,{\overline{vy}}_1, {\overline{vy}}_2, \dots , {\overline{vy}}_m)$, and $({\overline{\theta }}_1, {\overline{\theta }}_2, \dots , {\overline{\theta }}_m,{\overline{v\theta }}_1, {\overline{v\theta }}_2, \dots , {\overline{v\theta }}_m)$, where ${\overline{vx}}_i$, ${\overline{vy}}_i$, and ${\overline{v\theta }}_i$ denote the particle velocity in the ith dimension in the corresponding m-dimensional space, respectively. According to the particle position range of each dimension in each m-dimensional search space, the values of ${\overline{vx}}_i$, ${\overline{vy}}_i$, and ${\overline{v\theta }}_i$ are all between –1 and 1. Moreover, for any given particle, when extracting $x_i$ and $y_i$ according to Equations (1), (2), (19), and (20), their position accuracy is all set to the millimeter. So is the accuracy of the pick-up and drop-off point coordinates.

4.2. Steps of the MOPSO

The classical PSO [41] is mainly for single-objective optimization problems (SOPs), and its various improvements for solving MOPs can all be referred to as MOPSO. In PSO, the new position of each particle is determined by the particle velocity, the best position found by the particle, and the position of the group best. For an SOP, it is easy to judge if a particle is the best one in the swarm, and its local best position needs to be updated by comparing the objective values. However, owing to the conflicting optimization objectives in an MOP, the differences between MOPSO and single-objective PSO mainly lie in how to judge if an individual is better than others. Correspondingly, if an MOP is not transformed into an SOP by the common methods such as weighted summation and goal programming, the concept of Pareto domination is usually applied to compare two solutions, which is also adopted in the MOPSO applied in this study. Specifically, based on the objective values of each particle, a particle is defined as dominating the other if its objective values are not worse than the other’s, and it is better than the other at least in one objective.

Moreover, for an FLP in a continuous space, the univariate search with controlled convergence has been reported as an effective strategy in minimizing the size of the block enveloping all facilities and eliminating overlap between facilities simultaneously [42]. Inspired by this, a local search mechanism is designed to execute the local search based on the Pareto solution in the MOPSO. Based on our previous work [43], the flowchart of the MOPSO adopted in this study is shown in Figure 3, and the following illustrates its key processes.

4.2.1. Swarm Initialization

Swarm initialization is to create particles that conform to the designed encoding scheme. Each particle is encoded in three layers, and the initialization method of each layer is the same. Take the layer to determine the x coordinate of each facility as an example, select m random numbers belonging to [0, 1] and assign them to ${\overline{x}}_1, {\overline{x}}_2, \dots , {\overline{x}}_m$ in turn. Then, select m random numbers belonging to [−1, 1] and assign them to ${\overline{vx}}_1, {\overline{vx}}_2, \dots , {\overline{vx}}_m$ in turn.

4.2.2. Fitness Evaluation

The idea of Pareto domination is employed to evaluate the particles in the swarm. For any two particles, it is necessary to obtain their respective objective values before judging the Pareto dominance between them. However, the facility layout extracted from a particle by Formula (19) may be illegal, i.e., constrains (16) and (17) are not satisfied, and facility overlap occurs in the layout scheme. As an evolutionary algorithm, PSO is proposed for unconstrained continuous optimization problems. Considering the penalty-function method is a common method to transform a constrained optimization problem into an unconstrained optimization problem, it is adopted here to handle constraints (16) and (17). Correspondingly, the two original optimization objectives $D_{\mathrm{total}}$ and $E_{\mathrm{total}}$ are transformed into $D_{\mathrm{total}}^{*}$ and $E_{\mathrm{total}}^{*}$, respectively, which are utilized in the Pareto dominance judgment and crowding distance calculation and expressed as:

(21)$D_{\mathrm{total}}^{*}=D_{\mathrm{total}}+\lambda H_{\mathrm{D}}$

(22)$E_{\mathrm{total}}^{*}=E_{\mathrm{total}}+\lambda H_{\mathrm{E}}$

where $H_{\mathrm{D}}$ and $H_{\mathrm{E}}$ represent two positive large numbers; $\lambda$ is the average violation of a facility layout, which can be presented as:

(23)$\lambda =\frac{1}{m^2}{\displaystyle \sum _{i=1}^m{\displaystyle \sum _{j=1}^m{\lambda }_{ij}}}$

where ${\lambda }_{ij}$ is a coefficient reflecting the violation between facility i and facility j, which comprehensively reflects the violations of such two facilities along the x-axis and y-axis, denoted as ${\lambda }_{ij}^x$ and ${\lambda }_{ij}^y$, respectively. Correspondingly, ${\lambda }_{ij}$, ${\lambda }_{ij}^x$, and ${\lambda }_{ij}^y$ can be formulated as follows:

(24)${\lambda }_{ij}=\min \{{\lambda }_{ij}^x, {\lambda }_{ij}^y\}, \forall i, j$

(25)${\lambda }_{ij}^x=\max \{0, 1-\frac{\left|x_i-x_j\right|}{\frac{1}{2}(l_i+l_j)+\max \{{\Delta }_i, {\Delta }_j\}}\}, \forall i, j$

(26)${\lambda }_{ij}^y=\max \{0, 1-\frac{\left|y_i-y_j\right|}{\frac{1}{2}(w_i+w_j)+\max \{{\Delta }_i, {\Delta }_j\}}\}, \forall i, j$

Note that the $l_i$ and $l_j$ in Formula (25), and the $w_i$ and $w_j$ in Formula (26) are determined by Formulas (1) and (2), which rely on the facility orientation information extracted by the particle encoding scheme. Moreover, according to Formulas (24)–(26), if a facility layout scheme is feasible, $\lambda$ is equal to 0, and $D_{\mathrm{total}}^{*}$ and $E_{\mathrm{total}}^{*}$ have practical meaning, respectively.

4.2.3. Crowding Distance Sorting

The concept of crowding distance is often applied in multi-objective optimization algorithms evaluating solutions by comparing the Pareto dominance relationship, e.g., non-dominated sorting genetic algorithm-II (NSGA-II) [44], which reflects the density of other solutions around one solution in solution search space. Through the Pareto-dominance judgment among particles, the particles in the swarm can be divided into different groups according to the number of particles dominated by each particle. Note that the crowding distance calculation can only be executed for particles belonging to the same group. The non-dominated particles are distributed in the Pareto front, and their crowding distances are the focus of the MOPSO. Specifically, all non-dominated particles are sorted by each objective value calculated by Formula (21) or (22); the crowding distance of two boundary particles is set to infinity, respectively, then the crowding distance of each intermediate particle is the sum of the distance difference between its adjacent particles on each objective. Furthermore, after the crowding distance sorting in descending order, the top 10% of non-dominated particles are picked up to form the elite solution set. Correspondingly, in the subsequent particle velocity and position update processes, a particle randomly selected from the elite solution set is viewed as the group best, which is beneficial to guiding the swarm particles to move towards the Pareto front and maintaining the diversity of solutions.

4.2.4. Velocity and Position Update

According to the particle encoding scheme, the velocity of particle k in a swarm with N particles needs to be updated from three aspects, as shown below.

(27)$\begin{array}{l}\{\begin{array}{l}{\overline{vx}}_i^k\left(p+1\right)=\omega {\overline{vx}}_i^k\left(p\right)+C_{\mathrm{L}}{\epsilon }_1\left[{\overline{x}}_i^{k, \mathrm{best}}\left(p\right)-{\overline{x}}_i^k\left(p\right)\right]+C_{\mathrm{G}}{\zeta }_1\left[{\overline{x}}_i^{\mathrm{gbest}}\left(p\right)-{\overline{x}}_i^k\left(p\right)\right]\\ {\overline{vy}}_i^k\left(p+1\right)=\omega {\overline{vy}}_i^k\left(p\right)+C_{\mathrm{L}}{\epsilon }_2\left[{\overline{y}}_i^{k, \mathrm{best}}\left(p\right)-{\overline{y}}_i^k\left(p\right)\right]+C_{\mathrm{G}}{\zeta }_2\left[{\overline{y}}_i^{\mathrm{gbest}}\left(p\right)-{\overline{y}}_i^k\left(p\right)\right]\\ {\overline{v\theta }}_i^k\left(p+1\right)=\omega {\overline{v\theta }}_i^k\left(p\right)+C_{\mathrm{L}}{\epsilon }_3\left[{\overline{\theta }}_i^{k, \mathrm{best}}\left(p\right)-{\overline{\theta }}_i^k\left(p\right)\right]+C_{\mathrm{G}}{\zeta }_3\left[{\overline{\theta }}_i^{\mathrm{gbest}}\left(p\right)-{\overline{\theta }}_i^k\left(p\right)\right]\end{array},\\ k=1, 2, \dots , N, i=1, 2, \dots , m, p=1, 2, \dots , I_{\max }\end{array}$

where p is the iteration index; $\omega$ is the inertia weight coefficient belonging to [0, 1]; C_L and $C_{\mathrm{G}}$ are cognitive and social acceleration coefficients, respectively; in the three m-dimensional search spaces determining the x coordinate, y coordinate, and orientation of each facility, respectively, ${\overline{vx}}_i^k\left(p\right)$, ${\overline{vy}}_i^k\left(p\right)$, and ${\overline{v\theta }}_i^k\left(p\right)$ are the velocity of particle k, and ${\overline{x}}_i^k\left(p\right)$, ${\overline{y}}_i^k\left(p\right)$, and ${\overline{\theta }}_i^k\left(p\right)$ are the position of particle k in the ith dimension of each search space in the pth iteration; ${\overline{x}}_i^{k, \mathrm{best}}\left(p\right)$, ${\overline{y}}_i^{k, \mathrm{best}}\left(p\right)$, and ${\overline{\theta }}_i^{k, \mathrm{best}}\left(p\right)$ are the best position of particle k in the ith dimension of each search space until the pth iteration, and ${\overline{x}}_i^{\mathrm{gbest}}\left(p\right)$, ${\overline{y}}_i^{\mathrm{gbest}}\left(p\right)$, and ${\overline{\theta }}_i^{\mathrm{gbest}}\left(p\right)$ are the best position of the swarm best particle in the ith dimension of each search space until the pth iteration; ${\epsilon }_1$, ${\epsilon }_2$, ${\epsilon }_3$, ${\zeta }_1$, ${\zeta }_2$, and ${\zeta }_3$ are uniform variables between 0 and 1; $I_{\max }$ is the maximum iterations.

As shown in Formula (27), the particle velocity update is composed of three parts. The first part represents that a particle tends to maintain its previous velocity, which is to ensure the global convergence performance of the algorithm; the middle part reflects a particle’s memory of its historical experience and the tendency to approach its best position in history; the last part reflects the group’s historical experience of collaborative cooperation and knowledge sharing among particles, and the particle tends to approach the historical best position of the group. Correspondingly, the second and third parts make the algorithm own local convergence ability. Therefore, the values of $\omega$, $C_{\mathrm{L}}$, and $C_{\mathrm{G}}$ maintain the balance between the global and local search ability together. Specifically, a larger inertia weight is conducive to jumping out of the local optimum and facilitating the global search, while a smaller one is beneficial to local search and promoting algorithm convergence. To avoid premature convergence of the algorithm and the oscillation near the global optimal solution at the later stage of operation, the linearly varying weight [45] is utilized in this study. Accordingly, the $\omega$ in Formula (27) changes with the algorithm iteration, which can be redefined as $\omega (p)$, i.e., the inertia weight in the pth iteration, and presented as:

(28)$\omega \left(p\right)={\omega }_{\max }-p\times \frac{{\omega }_{\max }-{\omega }_{\min }}{I_{\max }}, p=1, 2, \dots , I_{\max }$

where ${\omega }_{\max }$ and ${\omega }_{\min }$ represent the maximum and minimum values of inertia weight, respectively. Moreover, after the velocity update, the new velocity of each particle in each dimension of each search space beyond the pre-set bound is reassigned to its nearest bound, i.e., −1 or 1.

Then, the position of each particle is updated as follows:

(29)$\{\begin{array}{l}{\overline{x}}_i^k\left(p+1\right)={\overline{x}}_i^k\left(p\right)+{\overline{vx}}_i^k\left(p+1\right)\\ {\overline{y}}_i^k\left(p+1\right)={\overline{y}}_i^k\left(p\right)+{\overline{vy}}_i^k\left(p+1\right)\\ {\overline{\theta }}_i^k\left(p+1\right)={\overline{\theta }}_i^k\left(p\right)+{\overline{v\theta }}_i^k\left(p+1\right)\end{array}, k=1, 2, \dots , N, i=1, 2, \dots , m, p=1, 2, \dots , I_{\max }$

Similarly, after each iteration, the new position of each particle in each dimension of each search space exceeding the pre-set bound is set to its nearest bound, i.e., 0 or 1. Accordingly, a new facility layout may be generated by Formula (19), and the corresponding objective values can be acquired by Formulas (21) and (22). Furthermore, whether the local best position of each particle and the global best position need to be updated can be determined through the Pareto-dominance judgment.

4.2.5. Local Search

In each iteration, local search is only executed for the group best particle, i.e., the non-dominated particle, which aims to eliminate the possible facility overlap in the contemporary optimal solution as much as possible and provide a better position near the contemporary optimal particle to guide the swarm particles’ movement. For a randomly selected non-dominated particle in the contemporary swarm, the local search method is performed according to the following steps:

Step 1: Extract the particle position information and sort all facility numbers randomly.

Step 2: Change the particle position relevant to facility i in the qth position, $q=1, 2, \dots , m$, through nine actions. Initially, q is assigned to 1.

• Action 1: Change the value of ${\overline{\theta }}_i$ to 0.1, 0.35, 0.6, and 0.85, respectively, and calculate the values of $D_{\mathrm{total}}^{*}$ and $E_{\mathrm{total}}^{*}$ corresponding to four particles transformed from the original one. Then, judge the Pareto dominance among them, select a non-dominated particle randomly, and record it.

• Action 2: Change the value of ${\overline{y}}_i$ by summing a random number between 0 and 0.5 and the original ${\overline{y}}_i$. Then, execute action 1, and record the best particle obtained.

• Action 3: Change the value of ${\overline{y}}_i$ by subtracting a random number between 0 and 0.5 from the original ${\overline{y}}_i$. Then, execute action 1, and record the best particle obtained.

• Action 4: Change the value of ${\overline{x}}_i$ by summing a random number between 0 and 0.5 and the original ${\overline{x}}_i$. Then, execute action 1, and record the best particle obtained.

• Action 5: Change the value of ${\overline{x}}_i$ by subtracting a random number between 0 and 0.5 from the original ${\overline{x}}_i$. Then, execute action 1, and record the best particle obtained.

• Action 6: Change the value of ${\overline{x}}_i$ by summing a random number between 0 and 0.5 and the original ${\overline{x}}_i$, and the value of ${\overline{y}}_i$ by summing a random number between 0 and 0.5 and the original ${\overline{y}}_i$. Then, execute action 1, and record the best particle obtained.

• Action 7: Change the value of ${\overline{x}}_i$ by subtracting a random number between 0 and 0.5 from the original ${\overline{x}}_i$ and the value of ${\overline{y}}_i$ by summing a random number between 0 and 0.5 and the original ${\overline{y}}_i$. Then, execute action 1, and record the best particle obtained.

• Action 8: Change the value of ${\overline{x}}_i$ by summing a random number between 0 and 0.5 and the original ${\overline{x}}_i$ and the value of ${\overline{y}}_i$ by subtracting a random number between 0 and 0.5 from the original ${\overline{y}}_i$. Then, execute action 1, and record the best particle obtained.

• Action 9: Change the value of ${\overline{x}}_i$ by subtracting a random number between 0 and 0.5 from the original ${\overline{x}}_i$ and the value of ${\overline{y}}_i$ by subtracting a random number between 0 and 0.5 from the original ${\overline{y}}_i$. Then, execute action 1, and record the best particle obtained.

Note that after the position update relevant to ${\overline{x}}_i$ and/or ${\overline{y}}_i$ in actions 2~9, it is necessary to judge whether the updated ${\overline{x}}_i$ and/or ${\overline{y}}_i$ exceeds its position bound. If so, ${\overline{x}}_i$ and/or ${\overline{y}}_i$ is set to its nearest bound, i.e., 0 or 1.

Step 3: Judge the Pareto dominance among the nine particles obtained in Step 2, and select a non-dominated particle randomly.

Step 4: Judge the Pareto dominance between the original particle and the particle acquired in Step 3. If the former is dominated by the latter, its current position is replaced by the latter’s position. Besides, if the local best of the original particle is also dominated by the particle obtained in Step 3, its position is replaced by the position of the particle obtained in Step 3.

Step 5: $q=q+1$. If q is not greater than m, go to Step 2; otherwise, go to the next step.

Step 6: Update the non-dominated particles of the contemporary population.

As noted above, the generation of new particles around the current global best particle mainly depends on the nine actions. Specifically, the first action only involves the change in facility orientation. The other actions may change not only the center coordinates but also the orientation of each facility. When changing the position encoding information related to facility i, the possible facility position adjustments corresponding to such nine actions are described in Table 1.

5. Experiments

To evaluate the relationship between the two optimization objectives in the EFL model and verify the effectiveness of the proposed model solution method, two experiments were executed, respectively.

5.1. Experiment 1

This experiment stemmed from the actual needs of a military manufacturing enterprise in Xi’an, China. The products produced by this enterprise mainly include missiles, guided projectiles, aviation parts, and auto parts. For the sake of confidentiality, the enterprise’s name, specific equipment, and the products involved are hidden in this paper. Due to the special requirements of high quality and high precision of military products, CNC machine tools have become essential tools for military enterprises to produce key parts and develop new products. Correspondingly, with the implementation of the efficiency improvement plan for CNC machine tools and the construction of a digital factory, this enterprise accumulated rich experience in automation and informatization application and owned strong product research and development and manufacturing capacity. In general, this enterprise’s production mode conformed to multi-variety and small/variable-batch. Meanwhile, to adapt to the manufacturing of products with relatively stable annual demand and certain process similarities, the enterprise tried to set up some flexible manufacturing cells (FMCs). Moreover, the manufacturing plants of the enterprise were mainly located in the city. With the increasingly strict environmental protection requirements in the urban area, the enterprise planned to move the plants to the suburbs to guarantee the production capacity and put forward the requirement of the re-layout of existing workshop facilities.

Under this background, we selected an FMC to verify the energy-saving effect of the proposed EFL model. According to the actual measurement and the data provided by the enterprise, the information of the facilities (i.e., various machine tools) belonging to this FMC in their original orientations (i.e., ${\theta }_i=0$) is shown in Table 2, and the information of typical products manufactured by this FMC is presented in Table 3. Regarding the AGV serving this FMC, it had four wheels directly driven by four servo motors, and its traveling speed was often set to 1 m/s when transporting materials. The AGV’s technical manual indicates that the AGV’s net weight is 60 kg, the maximum handling load is 80 kg, and the driving motor efficiency is 0.9.

The AGV contains multiple energy-consuming components. To help analyze the AGV’s EC characteristics, an experimental setup mainly consisting of a DC power sensor module, two Bluetooth communication modules, and a laptop with the EC data acquisition software based on LabVIEW was built. Its architecture and installation on the AGV are shown in Figure 4, and the sampling frequency of power data was set to 40 Hz. The measured $P_{\mathrm{sm}}$ was 25 W, and the relevant energy sources mainly include an STM32 microcomputer, magnetic navigation sensors, motor drivers, fans, and signal lamps. Regarding the AGV uniform motion power ($P_{\mathrm{um}}$, W), it should be stable under the given transport load and speed according to physical kinematics analysis. That is, when the AGV moves at a uniform speed, the motor power output on top of its standby power is utilized to overcome the rolling resistance, so $P_{\mathrm{um}}$ can be expressed as:

(30)$P_{\mathrm{um}}=\frac{1}{\eta }{C}_{\mathrm{r}}{m}_{\mathrm{total}}g{v}_a$

where $m_{\mathrm{total}}$ represents the total mass of the AGV and goods, kg. Correspondingly, $C_{\mathrm{r}}$ is an important parameter that needs to be determined in advance for evaluating $E_{\mathrm{um}}$. Although its value or range is usually acquired by referring to the mechanical design manual in engineering practice, it cannot be accurately obtained through the manual under this experimental condition. In the manufacturing workshop of this enterprise, the ground was usually sprayed with high wear-resistant polyurethane floor paint. Given the AGV, $C_{\mathrm{r}}$ can be regarded as a constant. Consequently, if the AGV travels at the same speed with different loads, $P_{\mathrm{um}}$ should be proportional to the load. Fortunately, the collected AGV power data verified this rule. For instance, when the AGV’s moving speed was 1.1 m/s and $m_{\mathrm{total}}$ was 72.5 kg, 85 kg, 97.5 kg, and 110 kg, respectively, the measured $P_{\mathrm{um}}$ was 39.5 W, 43.7 W, 50.7 W, and 56.9 W, separately. Through the statistical analysis of the measured experimental data, the value of $C_{\mathrm{r}}$ was 0.03. This also indirectly proved that it is somewhat reasonable to treat the power factor of the driving motor as a constant in the AGV’s technical manual when the AGV travels at the specified allowable speed. Furthermore, Formula (12) can be directly utilized to evaluate the $E_{\mathrm{um}}$ corresponding to each layout scheme in the MOPSO running process.

Then, the frequency of material flow, and the mass flow between different facilities could be obtained from Table 3, as shown in Table 4 and Table 5, respectively. Note that the AGV’s weight was also taken into account in the mass flow analysis. The layout area covered by these facilities in the new workshop should not exceed 25 m × 20 m. Besides, the basic parameter setting of the MOPSO, which was implemented in Matlab language on a PC with Intel(R) Core(^TM) i5-6500 CPU @ 3.20 GHz, 16 GB RAM, and Windows 7 OS, is depicted in Table 6. The dimensions of $H_{\mathrm{D}}$ and $H_{\mathrm{E}}$ varied as they were used to evaluate different objectives, but their values were uniformly set to 10²⁰ when calculating the objective values of each particle. After running the MOPSO ten times, it was observed that a unique optimal solution was output each time, and the average search time of ten replications was 83.265 s. The FMC adopted a double-row layout before the workshop was relocated, and the placement and orientation of each facility mainly depended on the experience of process designers and operators. The specific comparison of the layout extracted from the best solution acquired by running the MOPSO ten times, and the original layout is illustrated in Figure 5.

As shown in Figure 5, the rectangle that envelops all facilities including the safety clearance space can be regarded as a whole in the layout planning area. Accordingly, it can be moved along the x-axis or y-axis and rotated like each facility as needed. Though the rectangular area enveloping all facilities was increased by 6.42% after the layout optimization, the optimization objectives concerned in the EFL model, i.e., $D_{\mathrm{total}}$ and $E_{\mathrm{total}}$, were reduced by 38.32% and 39.20%, respectively. Therefore, the energy-saving effect achieved by optimizing the facility layout was significant, and the total load transport distance was also optimized.

In addition, the facilities in the optimized layout were relatively compact, and facility 2 was approximately in the center of the layout due to the highest frequency of the material flow between facility 2 and other facilities. Considering the requirements of facility management and workshop transport network planning, it is difficult to directly apply the optimized layout in practice. Figure 6 depicts a feasible adjustment scheme based on the optimized layout. The positions of some facilities moved along the x-axis so that they were arranged in columns, which can also be interpreted as rows when the envelop rectangle was rotated. Furthermore, it can be found that the rectangular area enveloping all facilities, and the values of $D_{\mathrm{total}}$ and $E_{\mathrm{total}}$, were all increased after adjusting the optimized layout. However, the altered layout was still better than the original ($D_{\mathrm{total}}$ decreased by 27.59% and $E_{\mathrm{total}}$ fell by 29.09%), which means the proposed EFL model can also assist in providing one or more initial layouts for subsequent adjustments in combination with actual situations.

Further, from the viewpoint of motion, the decomposition of the $E_{\mathrm{total}}$ corresponding to the three layout schemes in Figure 5 and Figure 6 is depicted in Figure 7. Comparing the optimized layout with the original layout, it can be observed that $E_{\mathrm{sm}}$ and $E_{\mathrm{um}}$ are both reduced, and the decline of $E_{\mathrm{um}}$ is greater than that of $E_{\mathrm{sm}}$. After adjusting the optimized layout, the distribution of equipment becomes a little loose, leading to the increase of $E_{\mathrm{sm}}$ and $E_{\mathrm{um}}$ of the optimized layout. However, $E_{\mathrm{sm}}$ has a higher percentage of growth than $E_{\mathrm{um}}$. Hence, when optimizing $E_{\mathrm{total}}$, it is inappropriate to focus only on a certain component, and all EC components shall be considered globally.

5.2. Experiment 2

The MOPSO applied in this study was the improvement of our previous work, and its performance had been compared with the common benchmark algorithm NSGA-II in solving some representative MOPs such as AGV path planning and process parameter optimization [8,43]. Theoretically, it is necessary to compare the proposed MOPSO with some other typical meta-heuristic methods to solve the same FLP. However, there were few suitable cases after the literature review, as EC or EC-related factors were rarely involved. Note that when applying an intelligent algorithm to solve an optimization problem, the design of an individual coding scheme reflecting the feasible problem solution directly affects the algorithm efficiency. Correspondingly, the particle encoding scheme proposed in this study should be suitable for solving FLPs, especially unequal-area static FLPs. Hence, Experiment 2 focused on verifying the effectiveness of the particle encoding scheme and evaluating the influence of the local search method on the solution quality. Generally, the main difference between MOPSO and single-objective PSO (SOPSO) lies in the individual fitness evaluation executed to determine the local best position and the group best position. Fortunately, there are massive cases considering the objective of transport distance in the existing research. So, two FLP cases focusing on the optimization of transport distance in [46], including 11 facilities within 15 m × 15 m layout area and 20 facilities within 14 m × 14 m layout area, respectively, were selected, as the scale of them was larger than our actual case. However, in such two cases, the transport distance between any two facilities was the Manhattan distance between their centers. The pick-up and drop-off points ignored, the transport distance was only related to the center position of each facility, as expressed below:

(31)$d_{ij}=\left|x_i-x_j\right|+\left|y_i-y_j\right|, \forall i, j$

Moreover, the original length and width of a facility can still be exchanged, which occurred when the value of ${\theta }_i$ extracted by Formula (19) was equal to π/2 or 3π/2, i.e., $r_i$ equals 1 at this moment.

Then, the SOPSO adopting the particle encoding scheme and the local search method proposed in this study was compared with other methods, which were the modified PSO [40], FLOAT [42], and TOPOPT [46]. Especially, to evaluate the performance of two PSO-based algorithms relatively objectively, the parameter setting of them were the same, except the inertia weight coefficient, as shown in Table 7. Meanwhile, the positioning accuracy was set to four decimal places. After running the SOPSO ten times for each FLP case, the best solutions are depicted in Figure 8, and the comparison of results obtained by various methods is summarized in Table 8. In terms of the best solution, the output result of the SOPSO is better than that of the TOPOPT method, and close to that of the FLOAT method, but worse than that of the modified PSO. Therefore, the particle encoding scheme applied in the MOPSO is effective.

5.3. Discussion

In Experiment 1, as each execution of the algorithm only output a unique optimal solution, the influence relationship between two optimization objectives cannot be directly judged through the Pareto front. Therefore, more data are needed to support the evaluation of their relationship. Correspondingly, plenty of feasible solutions generated in the process of population evolution can be utilized. Specifically, two feasible solutions belonging to the same Pareto rank were found in the algorithm iteration, and the facility layout schemes extracted by decoding the corresponding particles are illustrated in Figure 9. Obviously, the decrease of $D_{\mathrm{total}}$ was accompanied by the growth of $E_{\mathrm{total}}$. Moreover, from the mathematical expressions of the two optimization objectives, they are both linear functions of the transport distance between facilities. Therefore, if such two objectives are linearly-correlated, it can be deduced that all non-zero elements in the same position of the material flow frequency matrix and mass flow matrix, as shown in Table 4 and Table 5, must follow the same proportional relationship. Only in this case, shortening the load transport distance is equivalent to saving EC. However, this situation is a special case. At least in our experiment, it is not tenable. In actual workshop production, it is still difficult to meet this condition strictly due to the differences in product type, batch, and process planning. Therefore, the transport distance and EC can be regarded as two independent optimization objectives in EFL in most cases.

In addition, the existing studies have proved that the effectively utilized area can be an independent optimization goal in traditional FLPs. According to the change in the rectangular area enveloping all facilities with D_total or E_total in Experiment 1, the effectively utilized area can also be treated as an optimization objective in EFL.

Regarding Experiment 2, the average running time of the SOPSO is shorter than that of the modified PSO [40], which is mainly due to the improvement of computer performance. In terms of the output solution quality, although the overall performance of the SOPSO is inferior to that of the modified PSO, it is better than the TOPOPT method, which means the particle encoding scheme is effective, and the proposed MOPSO for solving the EFL model can be acceptable. Furthermore, by analyzing the algorithm process, it can be found that the modified PSO adopts different local search strategies in various iteration stages. The pick-up and drop-off points of each facility ignored, the local search method applied in the MOPSO is similar to the “Local search method 1” in the modified PSO in terms of facility movement actions. Since the local search ability of PSO is generally weak [47], it could be concluded that the lack of flexibility and diversity of local search methods affects the quality of final output solutions. Moreover, by solving the existing cases, it is verified that the particle encoding scheme proposed in this paper is suitable for solving FLPs, especially unequal-area static FLPs.

6. Conclusions

As the manufacturing industry gradually pays attention to sustainability, extensive efforts have been devoted to the energy-efficient manufacturing system. Facility layout is an essential part of manufacturing system planning. To fully excavate the energy-saving potential of a manufacturing system, it is necessary to own energy consciousness in its planning and design stage. Under this background, we investigated energy-efficient facility layout in a manufacturing workshop environment and an EFL model addressing logistics efficiency (i.e., total load transport distance), and total load transport EC was established. Especially, the operation characteristics and layout requirements of the equipment in the actual manufacturing environment are considered as much as possible in the modeling process, i.e., each facility has a pick-up point, a drop-off point, and a safety clearance space. Furthermore, by referring to the common solutions to FLPs, an MOPSO-based approach was put forward to solve the EFL model. Moreover, two experiments, one from the actual engineering requirements and the other from the existing FLP research, were conducted to verify the effectiveness of the proposed model and its solution method. The main research results of this study can be summarized as follows: (1) energy-efficient facility layout is a feasible approach to improve the energy-saving potential of a manufacturing system; (2) EC can be viewed as an independent optimization objective in facility layout; (3) the proposed particle encoding scheme is effective, and the MOPSO utilized in this study can make a trade-off between two optimization objectives objectively and provide a reference for similar research on facility layout methods.

In addition, some defects need attention to guide the follow-up research. Currently, there are few studies on EC modeling for AGVs in a manufacturing workshop environment, which makes it arduous to predict transport EC reliably. In this study, based on our previous work and some simplifications, the load transport EC was predicted from the perspective of motion but only the standby and uniform motions were considered. The AGV motion state is complex, and there are still acceleration, deceleration, and other motions. So, it is necessary to further study the EC characteristics of AGV motions and consider the empty transport EC to predict the AGV EC in the real transport process more reliably. Meanwhile, the current AGV EC data acquisition setup has some limitations in application. Under the Bluetooth communication mode, the effective communication distance was within 10 m, and data transmission was easily disturbed by machine tools and metal products widely existing in the workshop. Consequently, the AGV EC data acquisition setup needs further improvement, and the EC of different AGV energy-consuming components in the transport process also needs further analysis. Besides, the rectilinear distance from the pick-up point of one facility to the drop-off point of another facility is utilized as the transport distance between them, but material handling is usually forbidden to pass through the facility core operation area in practice. For example, the actual transport distance from facility 8 to facility 4 in the optimized layout shown in Figure 5 is inevitably greater than the transport distance calculated by Equation (6). Furthermore, it might be formidable to apply the output EFL scheme in practice directly. How to adjust it quickly to meet the actual needs and guarantee the energy-saving effect as much as possible is also a challenge. With the development of new-generation information technologies such as the Internet of things, big data, and cloud computing, it is gradually becoming possible to realize the integration and cooperation of the digital world and the physical world. Correspondingly, digital twin technology can provide powerful support for evaluating the practical application effect of a facility layout scheme [48]. Therefore, based on our current work, more efforts will be made in: (1) AGV EC characteristics and modeling; (2) improvement of the EFL model, e.g., considering the AGV empty transport EC, obstacle avoidance demand, and the facility orientation constraints for worker communication and human–machine interaction; (3) digital twin-enabled EFL; (4) improvement of the MOPSO.

Author Contributions

Conceptualization, Z.Z.; methodology, Z.Z. and W.Z.; software, Z.Z., W.Z. and T.P.; validation, Z.W., W.Z. and S.J.; formal analysis, Z.W. and S.J.; resources, T.P.; writing—original draft preparation, Z.Z.; writing—review and editing, L.W. and T.P.; visualization, L.W., Z.W. and S.J.; supervision, L.W.; project administration, T.P.; funding acquisition, Z.Z., L.W., Z.W. and S.J. All authors have read and agreed to the published version of the manuscript.

Funding

This research was funded by the National Natural Science Foundation of China (Grant No. U1704156 and 71971130), the Key Science and Technology Program of Henan Province (Grant No. 212102210357), and the Foundation of Henan University of Technology (Grant No. 2016RCJH02).

Institutional Review Board Statement

Not applicable.

Informed Consent Statement

Not applicable.

Data Availability Statement

Not applicable.

Acknowledgments

We sincerely thank all editors and reviewers for their valuable suggestions on the improvement of this paper.

Conflicts of Interest

The authors declare no conflict of interest.

Nomenclature

$$B_z$$	Average transport batch size of product z
$$C_{\mathrm{G}}$$	Social acceleration coefficient
$$C_{\mathrm{L}}$$	Cognitive acceleration coefficient
$$C_{\mathrm{r}}$$	Rolling resistance coefficient
$$d_{ij}$$	Transport distance from facility i to facility j (m)
$$D_{\mathrm{total}}$$	Total load transport distance (m)
$$D_{\mathrm{total}}^{*}$$	Total load transport distance with penalty (m)
$$E_{\mathrm{sm}}$$	AGV standby motion energy consumption (J)
$$E_{\mathrm{total}}$$	Total load transport energy consumption (J)
$$E_{\mathrm{total}}^{*}$$	Total load transport energy consumption with penalty (m)
$$E_{\mathrm{um}}$$	AGV uniform motion energy consumption (J)
$$f_{ij}$$	Frequency of material flow from facility i to facility j
g	Gravity acceleration, 9.81m/s2
i, j	Facility index
$$I_{\max }$$	Maximum iterations
k	Particle index
L	Length of a workshop layout area (m)
$$l_{0i}$$	The original length of facility i (m)
$$l_i$$	The actual length of facility i (m)
m	Number of facilities
$$m_0$$	Empty AGV mass (kg)
$$m_{\mathrm{total}}$$	Total mass of the AGV and cargo (kg)
$$m_z$$	Mass of product z (kg)
n	Number of product types
N	Swarm size
$$N_z$$	Number of transport batches for product z
$$O_z$$	The total output of product z
p	Algorithm iteration index
$$P_{\mathrm{sm}}$$	AGV standby motion power (W)
$$P_{\mathrm{um}}$$	AGV uniform motion power (W)
q	Particle encoding position index
$$r_i$$	Variable to indicate whether the original length and width of facility i exchange
$$T_{ij}$$	Total load transport time from facility i to facility j (s)
$$v_a$$	AGV transport speed [m/s]
$${\overline{vx}}_i$$	Particle velocity in the spatial dimension to which $${\overline{x}}_i$$ belongs
$${\overline{vy}}_i$$	Particle velocity in the spatial dimension to which $${\overline{y}}_i$$ belongs
$${\overline{v\theta }}_i$$	Particle velocity in the spatial dimension to which $${\overline{\theta }}_i$$ belongs
$${\overline{vx}}_i^k\left(p\right)$$	The velocity of particle k in the spatial dimension to which $${\overline{x}}_i^k\left(p\right)$$ belongs
$${\overline{vy}}_i^k\left(p\right)$$	The velocity of particle k in the spatial dimension to which $${\overline{y}}_i^k\left(p\right)$$ belongs
$${\overline{v\theta }}_i^k\left(p\right)$$	The velocity of particle k in the spatial dimension to which $${\overline{\theta }}_i^k\left(p\right)$$ belongs
W	Width of a workshop layout area (m)
$$w_{0i}$$	The original width of facility i (m(
$$w_i$$	The actual width of facility i (m)
$$X_{ijz}$$	Total transfer times from facility i to facility j in the process route of product z
$${\overline{x}}_i$$	Particle position in the dimension reflecting the x-axis coordinate of facility i center
$${\overline{x}}_i^{\mathrm{gbest}}\left(p\right)$$	Position of the swarm best particle in the spatial dimension reflecting the x-axis coordinate of facility i center until the pth iteration of the MOPSO
$${\overline{x}}_i^k\left(p\right)$$	Position of particle k reflecting the x-axis coordinate of facility i center in the pth iteration of the MOPSO
$${\overline{x}}_i^{k, \mathrm{best}}\left(p\right)$$	The best position of particle k in the spatial dimension reflecting the x-axis coordinate of facility i center until the pth iteration of the MOPSO
$$(x_i, y_i)$$	Center coordinates of facility i
$$x_i^{\min }$$,$$x_i^{\max }$$	The x-axis coordinate bounds of facility i center (m)
$$(x_i^{\mathrm{in}}, y_i^{\mathrm{in}})$$	Drop-off point coordinates of facility i
$$(x_i^{\mathrm{out}}, y_i^{\mathrm{out}})$$	Pick-up point coordinates of facility i
$$({\widehat{x}}_i^{\mathrm{in}}, {\widehat{y}}_i^{\mathrm{in}})$$	Coordinates of the drop-off point of facility i relative to its center
$$({\widehat{x}}_i^{\mathrm{out}}, {\widehat{y}}_i^{\mathrm{out}})$$	Coordinates of the pick-up point of facility i relative to its center
$$y_i^{\min }$$,$$y_i^{\max }$$	The y-axis coordinate bounds of facility i center (m)
$${\overline{y}}_i$$	Particle position in the dimension reflecting the y-axis coordinate of facility i center
$${\overline{y}}_i^{\mathrm{gbest}}\left(p\right)$$	Position of the swarm best particle in the spatial dimension reflecting the y-axis coordinate of facility i center until the pth iteration of the MOPSO
$${\overline{y}}_i^k\left(p\right)$$	Position of particle k reflecting the y-axis coordinate of facility i center in the pth iteration of the MOPSO
$${\overline{y}}_i^{k, \mathrm{best}}\left(p\right)$$	The best position of particle k in the spatial dimension reflecting the y-axis coordinate of facility i center until the pth iteration of the MOPSO
z	Product index
$${\Delta }_i$$	Safety clearance distance of facility i
$${\epsilon }_1$$, $${\epsilon }_2$$, $${\epsilon }_3$$,$${\zeta }_1$$,$${\zeta }_2$$,$${\zeta }_3$$	Random uniformly distributed numbers in [0, 1]
$$\eta$$	Overall power factor of AGV driving motors
$${\theta }_i$$	Orientation of facility i
$${\overline{\theta }}_i$$	Particle position in the dimension reflecting the orientation of facility i
$${\overline{\theta }}_i^{\mathrm{gbest}}\left(p\right)$$	Position of the swarm best particle in the spatial dimension reflecting the orientation of facility i until the pth iteration of the MOPSO
$${\overline{\theta }}_i^k\left(p\right)$$	Position of particle k reflecting the orientation of facility i in the pth iteration
$${\overline{\theta }}_i^{k, \mathrm{best}}\left(p\right)$$	The best position of particle k in the spatial dimension reflecting the orientation of facility i until the pth iteration of the MOPSO
$$\lambda$$	Average violation factor of a facility layout
$${\lambda }_{ij}$$	Violation coefficient between facility i and facility j
$${\lambda }_{ij}^x$$	Violation coefficient between facility i and facility j in the x-axis direction
$${\lambda }_{ij}^y$$	Violation coefficient between facility i and facility j in the y-axis direction
$$\omega$$	Inertia weight
$${\omega }_{\min }$$,$${\omega }_{\max }$$	Inertia weight bounds
$$\omega (p)$$	Inertia weight in the pth iteration of the MOPSO

Publisher’s Note: MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Figures and Tables

Figure 1. Diagram of energy-efficient facility layout.

Figure 2. Influence of the facility orientation on the coordinates of the pick-up and drop-off points.

Figure 3. Flowchart of the MOPSO.

Figure 4. Experimental setup for the AGV EC data acquisition.

Figure 5. Comparison of the original layout and the optimized layout.

Figure 6. A feasible adjustment scheme based on the optimized layout.

Figure 7. Decomposition of the AGV load transport EC corresponding to three layout schemes.

Figure 8. Optimal solutions of two FLP cases.

Figure 9. Comparison of two feasible solutions in Experiment 1.

Table 1

Possible facility position adjustments corresponding to nine actions.

Action No.	Affected Decision Variable	Possible Facility Position Adjustment
1	$${\theta }_i$$	rotate to exchange length and width
2	$$y_i$$ $$\mathrm{and} {\theta }_i$$	move upward, and rotate to exchange length and width
3	$$y_i$$ $$\mathrm{and} {\theta }_i$$	move downward, and rotate to exchange length and width
4	$$x_i$$ $$\mathrm{and} {\theta }_i$$	move right, and rotate to exchange length and width
5	$$x_i$$ $$\mathrm{and} {\theta }_i$$	move left, and rotate to exchange length and width
6	$$x_i$$ $$, y_i$$ $$, \mathrm{and} {\theta }_i$$	move toward top-right, and rotate to exchange length and width
7	$$x_i$$ $$, y_i$$ $$, \mathrm{and} {\theta }_i$$	move toward top-left, and rotate to exchange length and width
8	$$x_i$$ $$, y_i$$ $$, \mathrm{and} {\theta }_i$$	move toward bottom-right, and rotate to exchange length and width
9	$$x_i$$ $$, y_i$$ $$, \mathrm{and} {\theta }_i$$	move toward bottom-left, and rotate to exchange length and width

Table 2

The information of facilities in their original orientations.

No.	Original Length (l0i, m)	Original Width (w0i, m)	Safety Clearance Distance (Δi, m)	Drop-Off Point Coordinates ($${\widehat{\mathit{x}}}_{\mathit{i}}^{\mathbf{in}}, {\widehat{\mathit{y}}}_{\mathit{i}}^{\mathbf{in}}$$)	Pick-Up Point Coordinates ($${\widehat{\mathit{x}}}_{\mathit{i}}^{\mathbf{out}}, {\widehat{\mathit{y}}}_{\mathit{i}}^{\mathbf{out}}$$)	Machining Capacity
1	3.4	2.5	2	(−1.020, −1.250)	(1.700, 0.250)	thread machining
2	3.2	2.5	2	(−1.600, 0.750)	(1.600, 0.250)	hole drilling
3	2.8	2	1.5	(0.840, −1.000)	(1.400, 0.200)	plane milling
4	3.2	2.5	1.5	(−0.320, −1.250)	(0.960, 1.250)	slotting; contour milling
5	4.5	3	2.5	(−2.250, 0.900)	(−2.250, −0.900)	external/internal turning
6	3.5	2.5	2	(1.050, −1.250)	(−1.050, 1.250)	hole boring
7	2.5	2	1.5	(0.750, 1.000)	(−1.250, −0.200)	plane grinding
8	2	1.5	1	(−1.000, −0.150)	(1.000, −0.150)	chamfering; deburring

Table 3

The information of products.

No.	Annual Output (Oz, Piece)	Mass of Single Piece (mz, kg)	Average Transport Batch Size (Bz, Pieces/Time)	Number of Transport Batches (Nz)	Machining Process Route
1	720	4.6	15	48	Facility 5→Facility 4→Facility 5→Facility 2→Facility 1→Facility 8
2	450	2.5	15	30	Facility 3→Facility 4→Facility 2→Facility 7→Facility 8
3	300	3.1	15	20	Facility 5→Facility 1→Facility 3→Facility 2→Facility 4→Facility 7→Facility 8
4	700	3.7	15	47	Facility 3→Facility 2→Facility 6→Facility 7→Facility 4→Facility 8
5	400	2.6	15	27	Facility 5→Facility 3→Facility 4→Facility 2→Facility 3→Facility 7→Facility 8
6	900	2	20	45	Facility 3→Facility 7→Facility 2→Facility 6→Facility 1→Facility 8

Table 4

Frequency of material flow between facilities.

Facility	1	2	3	4	5	6	7	8
1	0	0	20	0	0	0	0	93
2	48	0	27	20	0	92	30	0
3	0	67	0	57	0	0	72	0
4	0	57	0	0	48	0	20	47
5	20	48	27	48	0	0	0	0
6	45	0	0	0	0	0	47	0
7	0	45	0	47	0	0	0	77
8	0	0	0	0	0	0	0	0

Table 5

Mass flow between facilities.

Facility	1	2	3	4	5	6	7	8
1	0	0	2130 kg	0	0	0	0	10,692 kg
2	6192 kg	0	2660 kg	2130 kg	0	9910 kg	2925 kg	0
3	0	7540 kg	0	5585 kg	0	0	7160 kg	0
4	0	5585 kg	0	0	6192 kg	0	2130 kg	5410 kg
5	2130 kg	6192 kg	2660 kg	6192 kg	0	0	0	0
6	4500 kg	0	0	0	0	0	5410 kg	0
7	0	4500 kg	0	5410 kg	0	0	0	7715 kg
8	0	0	0	0	0	0	0	0

Table 6

Parameter setting of the MOPSO.

Parameter		Value
Swarm size (N)		50
Maximum iteration (Imax)		500
Cognitive acceleration coefficient (CL)		2
Social acceleration coefficient (CG)		2
Inertia weight coefficient (ω)	Maximum value (ωmax)	0.9
	Minimum value (ωmin)	0.4

Table 7

Parameter setting of the SOPSO for two FLP cases.

Parameter	Value
	Case 1 with 11 Facilities	Case 2 with 20 Facilities
Swarm size (N)	20	50
Maximum iteration $$(I_{\max }$$)	300	400
Cognitive acceleration coefficient($$C_{\mathrm{L}}$$)	0.5	0.5
Social acceleration coefficient ($$C_{\mathrm{G}}$$)	1.2	1.2

Table 8

Comparison of results obtained by various methods.

	Method	The SOPSO	The Modified PSO [40]	FLOAT [42]	TOPOPT [46]
Case Example
FLP with 11 facilities	Best solution [m]	1331.9386	1286.1069	-	-
	Worst solution [m]	1463.5799	1371.3264	-	-
	Mean value of the optimal solution [m]	1413.7281	1335.63845	-	-
	Run time for the best solution [s]	18.674562	888.315646	-	-
	Run time for the worst solution [s]	18.897003	919.736577	-	-
	Average running time [s]	18.873241	-	-	-
FLP with 20 facilities	Best solution [m]	1288.3512	1206.6489	1264.94	1320.72
	Worst solution [m]	1473.6893	1315.2316	-	-
	Mean value of the optimal solution [m]	1348.6483	1264.21306	1333.81	1395.64
	Run time for the best solution [s]	69.381237	2352.12272	-	-
	Run time for the worst solution [s]	70.316711	2250.8654	-	-
	Average running time [s]	69.898404	-	-	-

References

1. Department of Energy Statistics, National Bureau of Statistics, People’s Republic of China. China Energy Statistical Yearbook 2020; China Statistics Press: Beijing, China, 2021.

2. Heragu, S.S.; Kusiak, A. Efficient models for the facility layout problem. Eur. J. Oper. Res.; 1991; 53, pp. 1-13. [DOI: https://dx.doi.org/10.1016/0377-2217(91)90088-D]

3. Sadrzadeh, A. A genetic algorithm with the heuristic procedure to solve the multi-line layout problem. Comput. Ind. Eng.; 2012; 62, pp. 1055-1064. [DOI: https://dx.doi.org/10.1016/j.cie.2011.12.033]

4. Anjos, M.F.; Vieira, M.V.C. Mathematical optimization approaches for facility layout problems: The state-of-the-art and future research directions. Eur. J. Oper. Res.; 2017; 261, pp. 1-16. [DOI: https://dx.doi.org/10.1016/j.ejor.2017.01.049]

5. Dai, M.; Tang, D.; Giret, A.; Salido, M.A.; Li, W.D. Energy-efficient scheduling for a flexible flow shop using an improved genetic-simulated annealing algorithm. Robot. Comput.-Integr. Manuf.; 2013; 29, pp. 418-429. [DOI: https://dx.doi.org/10.1016/j.rcim.2013.04.001]

6. Liu, Y.; Dong, H.; Lohse, N.; Petrovic, S.; Gindy, N. An investigation into minimising total energy consumption and total weighted tardiness in job shops. J. Clean. Prod.; 2014; 65, pp. 87-96. [DOI: https://dx.doi.org/10.1016/j.jclepro.2013.07.060]

7. Zhang, Z.; Wu, L.; Peng, T.; Jia, S. An improved scheduling approach for minimizing total energy consumption and makespan in a flexible job shop environment. Sustainability; 2019; 11, 179. [DOI: https://dx.doi.org/10.3390/su11010179]

8. Zhang, Z.; Wu, L.; Zhang, W.; Peng, T.; Zheng, J. Energy-efficient path planning for a single-load automated guided vehicle in a manufacturing workshop. Comput. Ind. Eng.; 2021; 158, 107397. [DOI: https://dx.doi.org/10.1016/j.cie.2021.107397]

9. Vichare, P.; Nassehi, A.; Kumar, S.; Newman, S.T. A unified manufacturing resource model for representing CNC machining systems. Robot. Comput.-Integr. Manuf.; 2009; 25, pp. 999-1007. [DOI: https://dx.doi.org/10.1016/j.rcim.2009.04.014]

10. Hosseini-Nasab, H.; Fereidouni, S.; Fatemi Ghomi, S.M.T.; Fakhrzad, M.B. Classification of facility layout problems: A review study. Int. J. Adv. Manuf. Technol.; 2018; 94, pp. 957-977. [DOI: https://dx.doi.org/10.1007/s00170-017-0895-8]

11. Asef-Vaziri, A.; Jahandideh, H.; Modarres, M. Loop-based facility layout design under flexible bay structures. Int. J. Prod. Econ.; 2017; 193, pp. 713-725. [DOI: https://dx.doi.org/10.1016/j.ijpe.2017.08.004]

12. Secchin, L.D.; Amaral, A.R.S. An improved mixed-integer programming model for the double row layout of facilities. Optim. Lett.; 2018; 13, pp. 193-199. [DOI: https://dx.doi.org/10.1007/s11590-018-1263-9]

13. Guan, C.; Zhang, Z.; Liu, S.; Gong, J. Multi-objective particle swarm optimization for multi-workshop facility layout problem. J. Manuf. Syst.; 2019; 53, pp. 32-48. [DOI: https://dx.doi.org/10.1016/j.jmsy.2019.09.004]

14. Abedzadeh, M.; Mazinani, M.; Moradinasab, N.; Roghanian, E. Parallel variable neighborhood search for solving fuzzy multi-objective dynamic facility layout problem. Int. J. Adv. Manuf. Technol.; 2013; 65, pp. 197-211. [DOI: https://dx.doi.org/10.1007/s00170-012-4160-x]

15. Balakrishnan, J.; Cheng, C.-H.; Wong, K.-F. FACOPT: A user friendly FACility layout OPTimization system. Comput. Oper. Res.; 2003; 30, pp. 1625-1641. [DOI: https://dx.doi.org/10.1016/S0305-0548(02)00096-5]

16. Bazaraa, M.S. Computerized layout design: A branch and bound approach. AIIE Trans.; 1975; 7, pp. 432-438. [DOI: https://dx.doi.org/10.1080/05695557508975028]

17. Seppänen, J.; Moore, J.M. Facilities planning with graph theory. Manag. Sci.; 1970; 17, pp. B-242-B-253. [DOI: https://dx.doi.org/10.1287/mnsc.17.4.B242]

18. Caccetta, L.; Kusumah, Y.S. Computational aspects of the facility layout design problem. Nonlinear Anal.; 2001; 47, pp. 5599-5610. [DOI: https://dx.doi.org/10.1016/S0362-546X(01)00662-9]

19. Chen, L.; Olhager, J.; Tang, O. Manufacturing facility location and sustainability: A literature review and research agenda. Int. J. Prod. Econ.; 2014; 149, pp. 154-163. [DOI: https://dx.doi.org/10.1016/j.ijpe.2013.05.013]

20. Muther, R. Systematic Layout Planning; 2nd ed. CBI Publishing Co. Inc.: Boston, MA, USA, 1973.

21. Bazaraa, M.S.; Sherali, H.D. On the use of exact and heuristic cutting plane methods for the quadratic assignment problem. J. Oper. Res. Soc.; 1982; 33, pp. 991-1003. [DOI: https://dx.doi.org/10.1057/jors.1982.210]

22. Hungerländer, P.; Anjos, M.F. A semidefinite optimization-based approach for global optimization of multi-row facility layout. Eur. J. Oper. Res.; 2015; 245, pp. 46-61. [DOI: https://dx.doi.org/10.1016/j.ejor.2015.02.049]

23. Kettani, O.; Oral, M. Reformulating quadratic assignment problems for efficient optimization. IIE Trans.; 1993; 25, pp. 97-107. [DOI: https://dx.doi.org/10.1080/07408179308964332]

24. Kusiak, A.; Heragu, S.S. The facility layout problem. Eur. J. Oper. Res.; 1987; 29, pp. 229-251. [DOI: https://dx.doi.org/10.1016/0377-2217(87)90238-4]

25. Lee, R.C.; Moore, J.M. CORELAP-computerized relationship layout planning. J. Ind. Eng.; 1967; 18, pp. 195-200.

26. Armour, G.C.; Buffa, E.S. A heuristic algorithm and simulation approach to relative allocation of facilities. Manag. Sci.; 1963; 9, pp. 294-309. [DOI: https://dx.doi.org/10.1287/mnsc.9.2.294]

27. Bozer, Y.A.; Meller, R.D.; Erlebacher, S.J. An improvement-type layout algorithm for single and multiple-floor facilities. Manag. Sci.; 1994; 40, pp. 918-932. [DOI: https://dx.doi.org/10.1287/mnsc.40.7.918]

28. Scriabin, M.; Vergin, R.C. A cluster-analytic approach to facility layout. Manag. Sci.; 1985; 31, pp. 33-49. [DOI: https://dx.doi.org/10.1287/mnsc.31.1.33]

29. Drezner, Z. DISCON: A new method for the layout problem. Oper. Res.; 1980; 28, pp. 1375-1384. [DOI: https://dx.doi.org/10.1287/opre.28.6.1375]

30. Matai, R. Solving multi objective facility layout problem by modified simulated annealing. Appl. Math. Comput.; 2015; 261, pp. 302-311. [DOI: https://dx.doi.org/10.1016/j.amc.2015.03.107]

31. Kothari, R.; Ghosh, D. Tabu search for the single row facility layout problem using exhaustive 2-opt and insertion neighborhoods. Eur. J. Oper. Res.; 2013; 224, pp. 93-100. [DOI: https://dx.doi.org/10.1016/j.ejor.2012.07.037]

32. Liu, J.; Liu, J. Applying multi-objective ant colony optimization algorithm for solving the unequal area facility layout problems. Appl. Soft Comput.; 2019; 74, pp. 167-189. [DOI: https://dx.doi.org/10.1016/j.asoc.2018.10.012]

33. Ripon, K.S.N.; Glette, K.; Khan, K.N.; Hovin, M.; Torresen, J. Adaptive variable neighborhood search for solving multi-objective facility layout problems with unequal area facilities. Swarm Evol. Comput.; 2013; 8, pp. 1-12. [DOI: https://dx.doi.org/10.1016/j.swevo.2012.07.003]

34. 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, pp. 621-655. [DOI: https://dx.doi.org/10.1007/s10479-016-2351-9]

35. Kumar, R.; Singh, S.P.; Lamba, K. Sustainable robust layout using Big Data approach: A key towards industry 4.0. J. Clean. Prod.; 2018; 204, pp. 643-659. [DOI: https://dx.doi.org/10.1016/j.jclepro.2018.08.327]

36. Iqbal, A.; Al-Ghamdi, K.A. Energy-efficient cellular manufacturing system: Eco-friendly revamping of machine shop configuration. Energy; 2018; 163, pp. 863-872. [DOI: https://dx.doi.org/10.1016/j.energy.2018.08.168]

37. Wang, R.; Zhao, H.; Wu, Y.; Wang, Y.; Feng, X.; Liu, M. An industrial facility layout design method considering energy saving based on surplus rectangle fill algorithm. Energy; 2018; 158, pp. 1038-1051. [DOI: https://dx.doi.org/10.1016/j.energy.2018.06.105]

38. Lamba, K.; Kumar, R.; Mishra, S.; Rajput, S. Sustainable dynamic cellular facility layout: A solution approach using simulated annealing-based meta-heuristic. Ann. Oper. Res.; 2020; 290, pp. 5-26. [DOI: https://dx.doi.org/10.1007/s10479-019-03340-w]

39. Zhang, S.; Gajpal, Y.; Appadoo, S.S.; Abdulkader, M.M.S. Electric vehicle routing problem with recharging stations for minimizing energy consumption. Int. J. Prod. Econ.; 2018; 203, pp. 404-413. [DOI: https://dx.doi.org/10.1016/j.ijpe.2018.07.016]

40. Derakhshan Asl, A.; Wong, K.Y. Solving unequal-area static and dynamic facility layout problems using modified particle swarm optimization. J. Intell. Manuf.; 2015; 28, pp. 1317-1336. [DOI: https://dx.doi.org/10.1007/s10845-015-1053-5]

41. Kennedy, J.; Eberhart, R. Particle swarm optimization. Proceedings of the 1995 IEEE International Conference on Neural Networks; Perth, Australia, 27 November–1 December 1995; pp. 1942-1948.

42. Imam, M.H.; Mir, M. Automated layout of facilities of unequal areas. Comput. Ind. Eng.; 1993; 24, pp. 355-366. [DOI: https://dx.doi.org/10.1016/0360-8352(93)90032-S]

43. Zhang, Z.; Wu, L.; Jia, S.; Peng, T. Multi-objective parameter optimization to support energy-efficient peck deep-hole drilling processes with twist drills. Int. J. Adv. Manuf. Technol.; 2020; 106, pp. 4913-4932. [DOI: https://dx.doi.org/10.1007/s00170-020-04967-x]

44. Deb, K.; Pratap, A.; Agarwal, S.; Meyarivan, T. A fast and elitist multiobjective genetic algorithm: NSGA-II. IEEE Trans. Evol. Comput.; 2002; 6, pp. 182-197. [DOI: https://dx.doi.org/10.1109/4235.996017]

45. Shi, Y.; Eberhart, R. Modified particle swarm optimizer. Proceedings of the 1998 IEEE International Conference on Evolutionary Computation, ICEC’98; Anchorage, AK, USA, 4–9 May 1998; pp. 69-73.

46. Imam, M.H.; Mir, M. Nonlinear programming approach to automated topology optimization. Comput.-Aided Des.; 1989; 21, pp. 107-115. [DOI: https://dx.doi.org/10.1016/0010-4485(89)90146-2]

47. Moslemipour, G.; Lee, T.S.; Rilling, D. A review of intelligent approaches for designing dynamic and robust layouts in flexible manufacturing systems. Int. J. Adv. Manuf. Technol.; 2012; 60, pp. 11-27. [DOI: https://dx.doi.org/10.1007/s00170-011-3614-x]

48. Leng, J.; Zhang, H.; Yan, D.; Liu, Q.; Chen, X.; Zhang, D. Digital twin-driven manufacturing cyber-physical system for parallel controlling of smart workshop. J. Ambient Intell. Humaniz. Comput.; 2019; 10, pp. 1155-1166. [DOI: https://dx.doi.org/10.1007/s12652-018-0881-5]