Concurrent scheduling and layout of virtual manufacturing cells considering assembly aspects

Abstract

Cell formation, scheduling, and facility layout are three main decisions in designing manufacturing cells. In this paper, we address the integration of these decisions in virtual manufacturing cells considering assembly aspects and process routing. We develop a mathematical model to determine the machine cells, the layout of machines and workstations on the shop floor, the processing route of parts, and the production sequence of operations on the machines. In this mathematical model, material handling costs and cycle time are minimized. To the best of our knowledge, this is the first paper that concurrently addresses the scheduling and layout of virtual manufacturing cells with assembly aspects and so-called criteria. To effectively solve the problem, a Population-based Simulated Annealing (PSA) combined with linear programming is proposed. The practical usability of the developed model is demonstrated in a case study. Finally, instances from the literature are solved to evaluate the performance of the PSA. The comparison results showed the superior performance of the PSA in comparison with CPLEX solver and standard simulated annealing.

Introduction

Cellular Manufacturing (CM) is proposed to maintain the flexibility of job shops and the efficiency of flow shops. ¹ CM seeks to implement group technology principles by centralizing similar production processes in manufacturing cells. In CM, this centralization is achieved through the Cell Formation (CF) process in which machines are grouped into machine cells so as to process a set of similar tasks. CM can contribute to the reduction of in-process inventories, cycle time, and number setup times. ² Also, it can improve product quality, increase flexibility, and simplify production control.^3,4

Despite having several advantages, one of the main disadvantages of CM is that it becomes less responsive as more changes occur in production conditions. ⁵ To deal with this problem, the concept of Virtual Manufacturing Cell (VMC) was introduced by McLean et al. ⁶ in the early 1980s. A traditional cell consists of a set of machines that are physically placed in a specific region of the shop floor. While, a virtual cell is not identifiable as a fixed physical group of machines but rather as data files and computer processes. ⁶ This structure makes VMCs more responsive against the changes in production conditions and prevents the need for costly and time-consuming reconfiguration of shop floor. ⁷

The design of a CM usually starts with the CF and then follows by the Facility Layout (FL) and scheduling. ⁸ The majority of the studies conducted on the CM design have only focused on one of these decisions; examples of studies in this relation can be found.^9–17 However, it should be noted that these decisions are highly inter-dependent. For example, the parameters of the FL and scheduling problems such as the material flow between machines and processing times are highly affected by the CF decision, or the transfer time between machines is directly proportional to the distance between machines. Noticing the importance of these inter-related decisions, some researchers applied an integrative approach to design CM more effectively. The studies conducted in this area can be classified into three main categories: studies addressing the CF jointly either with FL, see for example, in Refs.^8,18–23 or scheduling, see for example, in Refs.^24–30 and studies addressing all of these decisions in an integrated manner, see Table 1. Although the latter approach can produce better design, as Table 1 shows, few studies have been dedicated to the integration of the CF, FL, and scheduling problems. This is due to the excess computational complexity of such an integrated approach since multiple decisions need to be taken at a time. In what follows, we discuss the shortcoming in the studies given in Table 1.

• From Table 1, it is observed that all the related studies have only focused on the design of CM with classical cells, and no research has been done on the integration of the CF, FL, and scheduling problems in VMCs.

• In order to prevent the problem from being more complicated, the parameters such as the dimensions of facilities and the clearance among them have been ignored in the studies shown in Table 1. If the problem can be formulated and solved by considering these parameters, the resulting layout will be more accurate.

• The functionality of a CM is not necessarily limited to the production of parts. In the real world, some factories produce products by assembling parts. According to the survey conducted by Wemmerlov and Johnson ³⁷ at 46 factories in the US, 51% of the cells had joining/assembly operations. After production, the parts are transported to workstations to be used in the assembly process of the products. Despite the assembly operations are common in CM systems, no study has yet investigated this aspect with the scheduling and FL.

• In case assembly operations are taken into consideration, an objective function such as the cycle time can be a better representative of the assembly line’s performance, rather than the objective functions such as makespan or flow time. However, as can be seen in Table 1, this objective function has not been considered in the joint problem.

Table 1.

Design attributes of the studies on the joint CF, FL, and scheduling problem.

Authors	CT	PR	Facilities included		Material weight	Facilities dimension	Aisle width	Distance-based handling		AO	Objectives
			Machines	Workstations				Costs	Times
Wu et al. 31	CC		✓								Flow time
Arkat et al. 32	CC		✓					✓	✓		Makespan and handling cost
Arkat et al. 33	CC		✓					✓	✓		Flow time or handling cost
Fahmy 34	CC		✓					✓	✓		Flow time or makespan
Ebrahimi et al. 35	CC	✓	✓			✓		✓	✓		Makespan, tardiness, and handling cost
Rahimi et al. 36	CC	✓	✓						✓		Flow time
This paper	VC	✓	✓	✓	✓	✓	✓	✓	✓	✓	Handling cost and cycle time

CT: cell type; PR: process routing; AO: assembly operations; CC: classical cells; VC: virtual cells.

To fill the gaps mentioned above, this paper develops a problem to integrate the CF, FL, and scheduling decisions in VMCs, considering assembly aspects and process routing. In what follows, the developed problem is explained.

Problem statement and mathematical model

In this study, we address the scheduling and layout of a VMC in which manufacturing activities consist of two main activities: producing parts and assembling finished parts into final products. Each assembly operation is presented as a task. A directed network is used to define the preceding relation between the tasks. In this network, the nodes are associated with the tasks and the edges show the preceding relations among the tasks. Considering these explanations, we address the following decisions in an integrated problem:

1. CF: this decision is made in order to assign the machines to the cells and determine the process routes of parts.

2. FL: this decision is made to obtain the centroid coordinates of the facilities on a shop floor. In the developed model, the first objective function is the minimization of the material handling costs which is calculated for three activities: transferring parts among machines, transferring parts between machines and workstations, and transferring assembled units among workstations.

3. Scheduling: an operational level scheduling approach is utilized to determine the processing sequence of operations on the machines. The transportation time between facilities is taken into consideration in order to make the model more realistic. As the second objective function of the developed mathematical model, we minimize the cycle time.

It should be mentioned that he following assumptions are considered in the developed model:

• The sequence of operations required to produce each part is known in advance and the operations should be completed in the given order.

• The set of machines that are capable of processing each operation is known and only one machine is allowed to be selected for each operation.

• For each part type, the demand, weight, handling cost per unit distance, and transfer time per unit distance is known.

• The assembly network (explained above) is predefined.

• Facilities (including machines and workstations) are rectangular in shape and their dimensions are pre-specified.

• The number of cells and the maximum number of machines that can be assigned to a cell is known.

• The workstations are assigned to a virtual cell.

• The facilities are laid out on a rectangular-shaped planar shop floor in only one given direction.

• The clearance among facilities is pre-specified.

• The distance between facilities is calculated using the centroid-to-centroid rectilinear function.

• The transfer time between facilities is proportional to the distance between the corresponding facilities.

• The process times are known and setup times are included in the process times.

• Each machine can only process one operation at a time.

• Job splitting is now allowed and once a process is started, it must be completed without interruption.

Mathematical model

According to the explanations given, the problem is formulated as the following mixed-integer program (the notations are defined in the Appendix).

$$\min \rho \frac{H}{H_L}+(1-\rho )\frac{\overline{c}}{{\overline{c}}_L},$$

(1)

subject to:

$$\begin{array}{c}\multicolumn{1}{c}{H=\sum _{p\in P}\sum _{\genfrac{}{}{0ex}{}{o\in O(p)}{o\ne {\overline{O}}_p}}\sum _{m\in M(p,o)}\sum _{\genfrac{}{}{0ex}{}{m\prime \in M(p,o+1)}{m\prime \ne m}}{q}_p{w}_p^P{r}_{\mathrm{pom}}{r}_{p(o+1)m\prime }}\\ \multicolumn{1}{c}{\times (|x_m-x_{m\prime }|+\left|(y_m-y_{m\prime })\right|)}\\ \multicolumn{1}{c}{(c^E+(c^A-c^E)\sum _{c\in C\backslash \{c_0\}}{z}_{\mathrm{mc}}{z}_{m\prime c})}\\ \multicolumn{1}{c}{+\sum _{m\in S}\sum _{t\in T(m)}\sum _{p\in P(t)}\sum _{m\prime \in M(p,{\overline{O}}_p)}}\\ \multicolumn{1}{c}{q_p{w}_p^P{r}_{\mathrm{pom}}{r}_{p,{\overline{O}}_p,m}(|x_m-x_{m\prime }|+\left|(y_m-y_{m\prime })\right|)}\\ \multicolumn{1}{c}{+\sum _{m\in S}\sum _{m\prime \in S}\sum _{t\in T(m)}\sum _{t\prime \in T(m\prime )}}\\ \multicolumn{1}{c}{c^A(f_{\mathrm{tt}\prime }+f_{t\prime t})(|x_m-x_{m\prime }|+\left|(y_m-y_{m\prime })\right|),}\end{array}$$

(2)

$$\sum _{m\in M(p,o)}{r}_{\mathrm{pom}}=1,p\in P,o\in O(p),$$

(3)

$$\sum _{c\in C\backslash \{c_0\}}{z}_{\mathrm{mc}}=1,m\in M,$$

(4)

$$z_{\mathrm{mc}}=1,m\in S,c=c_0,$$

(5)

$$\sum _{m\in M}{z}_{\mathrm{mc}}\le N^M,c\in C\backslash \{c_0\},$$

(6)

$${\alpha }_{\mathrm{mm}\prime }+{\alpha }_{m\prime m}+{\beta }_{\mathrm{mm}\prime }+{\beta }_{m\prime m}=1,m,m\prime \in F,m\prime >m,$$

(7)

$$\begin{array}{c}{x}_{m\prime }-x_m\ge \frac{w_m+w_{m\prime }}{2}+L^E+(L^A-L^E)\sum _{c\in C}{z}_{\mathrm{mc}}{z}_{m\prime c}\hfill \\ -W(1-{\alpha }_{\mathrm{mm}\prime }),m,m\prime \in F,m\prime \ne m,\hfill \end{array}$$

(8)

$$\begin{array}{c}{y}_{m\prime }-y_m\ge \frac{l_m+l_{m\prime }}{2}+L^E+(L^A-L^E)\sum _{c\in C}{z}_{\mathrm{mc}}{z}_{m\prime c}\hfill \\ -L(1-{\alpha }_{\mathrm{mm}\prime }),m,m\prime \in F,m\prime \ne m,\hfill \end{array}$$

(9)

$$\frac{w_m}{2}\le x_m\le W-\frac{w_m}{2},m\in F,$$

(10)

$$\frac{l_m}{2}\le y_m\le L-\frac{l_m}{2},m\in F,$$

(11)

$$\begin{array}{c}{v}_{\mathrm{pp}\prime \mathrm{oo}\prime m}+v_{p\prime \mathrm{po}\prime \mathrm{om}}=1,p,p\prime \in P,p\prime >p,o\in O(p),\hfill \\ o\prime \in O(p\prime ),m\in M,\hfill \end{array}$$

(12)

$$\begin{array}{c}{c}_{p\prime o\prime \mathrm{mn}}^O\ge c_{\mathrm{pomn}}^O+q_{p\prime }{t}_{p\prime o\prime m}-B(3-r_{\mathrm{pom}}-r_{p\prime o\prime m}-v_{\mathrm{pp}\prime \mathrm{oo}\prime m}),\hfill \\ p,p\prime \in P,p\prime \ne p,o\in O(p),o\prime \in O(p\prime ),\hfill \\ m\in M,n\in N,\hfill \end{array}$$

(13)

$$\begin{array}{c}\multicolumn{1}{c}{c_{\mathrm{pomn}}^O\ge c_{p(o-1)m\prime n}^O+q_p{t}_{\mathrm{pom}}}\\ +(T^E+(T^A-T^E)\sum _{c\in C}{z}_{\mathrm{mc}}{z}_{m\prime c})\hfill \\ \times (|x_m-x_{m\prime }|+|y_m-y_{m\prime }|)\hfill \\ -B(2-r_{\mathrm{pom}}-r_{p(o-1)m\prime }),\hfill \\ p\in P,o\in O(p)\backslash \{1\},m,m\prime \in M,n\in N,\hfill \end{array}$$

(14)

$$\begin{array}{c}\multicolumn{1}{c}{c_{\mathrm{pomn}}^O\ge q_p{t}_{\mathrm{pom}}+E_{\mathrm{mn}}-B(1-r_{\mathrm{pom}}),p\in P,o\in O(p),}\\ m\in M,n\in N,\hfill \end{array}$$

(15)

$$\begin{array}{c}\multicolumn{1}{c}{E_{\mathrm{mn}}\ge c_{\mathrm{pom}(n-1)}^O-B(1-r_{\mathrm{pom}}),p\in P,o\in O(p),m\in M,}\\ n\in N\backslash \{1\},\hfill \end{array}$$

(16)

$$\begin{array}{c}{c}_{\mathrm{tn}}^T\ge c_{\mathrm{pomn}}^O+e_t+T^E(|x_m-x_{m\prime }|+|y_m-y_{m\prime }|)-B(1-r_{\mathrm{pom}}),\hfill \\ p\in P(t),o={\overline{O}}_p,m\in M,m\prime \in S,t\in T(m\prime ),n\in N,\hfill \end{array}$$

(17)

$$\begin{array}{c}{c}_{\mathrm{tn}}^T\ge c_{t\prime n}^T+e_t+T^A(|x_m-x_{m\prime }|+|y_m-y_{m\prime }|),m,m\prime \in S,\hfill \\ t\in T(M),t\prime \in (T(m\prime )\cap S(t)),n\in N,\hfill \end{array}$$

(18)

$$c_{\mathrm{tn}}^T\ge c_{t(n-1)}^T+e_t,t\in T,n\in N\backslash \{1\},$$

(19)

$$\overline{c}\ge \frac{c_{\mathrm{tn}}^T}{\overline{N}},t\in T,n=\overline{N},$$

(20)

$$\begin{array}{c}{c}_{\mathrm{pomn}}^P,c_{\mathrm{tn}}^T,E_{\mathrm{mn}}\ge 0,p\in P,o\in O(p),m\in M,\hfill \\ t\in T,n\in N,\hfill \end{array}$$

(21)

$$\begin{array}{c}{r}_{\mathrm{pom}},v_{\mathrm{pp}\prime \mathrm{oo}\prime m}\in \{0,1\},p,p\prime \in P,o\in O(p),o\prime \in O(p),\hfill \\ m\in M,\hfill \end{array}$$

(22)

$$z_{\mathrm{mc}},{\alpha }_{\mathrm{mm}\prime },{\beta }_{\mathrm{mm}\prime }\in \{0,1\},m,m\prime \in F,c\in C.$$

(23)

Objective function (1) minimizes a normalized function of the handling cost and cycle time. Constraint (2) determines the handling cost. Constraint sets (3) represents that each operation is processed by one machine. Constraint sets (4) ensures that each machine is assigned to a cell and constraint sets (5) guarantees that the workstations are assigned to virtual cell $c_0$ . Constraint sets (6) prevents the assignment of more than $N^M$ machines to each cell. Constraints sets (7)–(9) jointly ensure that facilities do not overlap on the shop floor and also there is enough clearance among them. Constraints sets (10) and (11) guarantee that the facilities are laid out within the boundaries of the shop floor. Constraint sets (12) states that each operation is scheduled either before or after another operation. Constraint sets (13) represents that only one operation can be processed by a machine at each moment. Constraint sets (14) states that each operation of a part can be processed after finishing and transferring the previous operation of the same part. Constraint sets (15) makes sure that the production in each cycle is started after finishing the production in the previous cycle. Constraint sets (16) states that the earliest available time of a machine in each production cycle is equal to the completion time of the last operation being processed by the same machine in the previous production cycle. Constraints sets (17) and (18) ensure that a task can be started after finishing and transferring the required parts and assembled units to the corresponding workstation. Constraint (19) represents that a task is started after finishing the same task in the previous production cycle. Constraint sets (20) ensures that the average cycle time in $\overline{N}$ production cycles is correctly calculated. Finally, constraints sets (21)–(23) indicate the type of decision variables.

Proposed solution algorithm

It is well known that the CF, FL, and scheduling are difficult optimization problems. In recent years, Simulated Annealing (SA) algorithm has been effectively applied to solve different types of combinatorial optimization problems, including those problems considered in this study, see for example in Refs.^9,19–23,38 In this section, we propose a Population-based Simulated Annealing (PSA) algorithm to solve the developed problem.

Encoding and decoding approaches

For solution encoding, a two-section structure is proposed. The first section consists of two layers: $\mathbf{R}$ and $\mathbf{S}$ . The information in layer $\mathbf{R}$ shows that which machine is assigned to each operations and the information in layer $\mathbf{S}$ indicates the processing sequence of operations on the machines. The second section consists of three layers: $\mathbf{Z}$ , $\mathbf{A}$ , and $\mathbf{B}$ . The first layer in this section, that is, $\mathbf{Z}$ , indicates which cell is assigned to each facility, and the last two layers, that is, $\mathbf{A}$ and $\mathbf{B}$ , consists of two permutations that are used to specify the relative position of facilities with respect to each other. ³⁹

For further clarification of the encoding approach, an example consisting of 18 operations, eight facilities (six machines and two workstations), and two cells is given in Figure 1. As shown in sections I and II of this figure, the operations are labeled from 1 to 18 and the facilities are labeled from 1 to 8; the first six labels in sections II are used for the machines and last two ones are used for the workstations. The information in layer $\mathbf{R}$ shows that operation #1 is processed on machine $m_5$ , operation #2 is processed on machine $m_3$ , and so on. Also from the information in layer $\mathbf{S}$ it is concluded that operation #15 is scheduled first, then operation #9 is scheduled, and so on. From layer $\mathbf{Z}$ in section sections II, it is seen that the group of machines $m_2$ , $m_3$ , $m_4$ , and $m_5$ are assigned to cell $c_1$ , the group of machines $m_1$ and $m_6$ are assigned to cell $c_2$ , and workstations $s_1$ and $s_2$ are assigned to imaginary cell $c_0$ . Finally, the two permutations shown in layers $\mathbf{A}$ and $\mathbf{B}$ are used to encode the relative position of facilities with respect to each other.

Figure 1.

Example of an encoded solution consisting of eighteen operations, eight facilities (including six machines and two workstations).

To decode layers $\mathbf{A}$ and $\mathbf{B}$ into a feasible layout of non-overlapping facilities, we proposed two algorithms: one based on the placement algorithm presented in Tang et al. ³⁹ and the other based on solving a Linear Program (LP). The following sub-sections explain these decoding algorithms.

First decoding algorithm

To decode layers $\mathbf{A}$ and $\mathbf{B}$ into a feasible layout for our problem, we modify Tang et al.’s ³⁹ algorithm. The modified algorithm, referred to as Algorithm 1, is given below.

Algorithm 1: Decoding algorithm based on Tang et al.’s 39 algorithm.

Input: $(\mathbf{Z},\mathbf{A},\mathbf{B})$ ;Output: $(\mathit{x},\mathit{y})$ ;for $m=1$ to $\overline{M}+\overline{S}$ do  $x_m\leftarrow l_m/2$ ;  $y_m\leftarrow w_m/2$ ;for $m=1$ to $\overline{M}+\overline{S}$ do  $i_1\leftarrow \mathrm{A}(m)$ ; for $m\prime =m+1$ to $\overline{M}+\overline{S}$ do   $i_2\leftarrow \mathrm{A}(m\prime )$ ;  if ${\mathrm{B}}^{-1}(i_2)>{\mathrm{B}}^{-1}(i_1)$ then   if $\mathrm{Z}(i_1)=\mathrm{Z}(i_2)$ then     $x_{i_2}\leftarrow \max \left\{x_{i_2},x_{i_1}+\left(l_{i_1}+l_{i_2}\right)/2+L^A\right\};$    else     $x_{i_2}\leftarrow \max \left\{x_{i_2},x_{i_1}+\left(l_{i_1}+l_{i_2}\right)/2+L^E\right\};$ for $m=\overline{M}+\overline{S}$ down to $1$ do  $i_1\leftarrow \mathrm{A}(m)$ ; for $m\prime =m-1$ down to $1$ do   $i_2\leftarrow \mathrm{A}(m\prime )$ ;  if ${\mathrm{B}}^{-1}(i_2)>{\mathrm{B}}^{-1}(i_1)$ then   if $\mathrm{Z}(i_1)=\mathrm{Z}(i_2)$ then     $\begin{array}{c}{y}_{i_2}\leftarrow max\{y_{i_2},y_{i_1}+(w_{i_1}+w_{i_2})/2+L^A\};\hfill \end{array}$    else     $\begin{array}{c}{y}_{i_2}\leftarrow max\{y_{i_2},y_{i_1}+(w_{i_1}+w_{i_2})/2+L^E\};\hfill \end{array}$

Second decoding algorithm

Although Algorithm 1 can quickly produce a feasible layout, it misses some solutions that might be optimal in terms of the material handling cost. To deal with this issue, we propose Algorithm 2 in which the material handling cost is minimized by solving the following LP.

$$\begin{array}{c}\mathrm{P}1:min\sum _{m\in F}\sum _{m\prime \in G(m)}(g_{\mathrm{mm}\prime }+g_{m\prime m})(v_{m\prime }-v_m)\hfill \\ +\sum _{m\in F}\sum _{m\prime \in F\backslash G(m)}{g}_{\mathrm{mm}\prime }(v_{\mathrm{mm}\prime }^{+}+v_{\mathrm{mm}\prime }^{-})+\gamma \sum _{m\in F}\mathrm{\Delta }{v}_m^{-},\hfill \end{array}$$

subject to:

$$\begin{array}{c}\multicolumn{1}{c}{v_{m\prime }-v_m\ge \frac{h_m+h_{m\prime }}{2}+L^A,m\in F,m\prime \in G(m),}\\ \mathrm{Z}(m)=\mathrm{Z}(m\prime ),\hfill \end{array}$$

$$\begin{array}{c}\multicolumn{1}{c}{v_{m\prime }-v_m\ge \frac{h_m+h_{m\prime }}{2}+L^E,m\in F,m\prime \in G(m),}\\ \multicolumn{1}{c}{\mathrm{Z}(m)\ne \mathrm{Z}(m\prime )}\end{array}$$

$$\begin{array}{c}\multicolumn{1}{c}{v_m-v_{m\prime }=v_{\mathrm{mm}\prime }^{+}-v_{\mathrm{mm}\prime }^{-},m\in F,m\prime \in F\backslash G(m),}\\ \multicolumn{1}{c}{v_m+\mathrm{\Delta }{v}_m^{+}-\mathrm{\Delta }{v}_m^{-}=U-\frac{h_m}{2},m\in F,}\\ \multicolumn{1}{c}{v_m,\mathrm{\Delta }{v}_m^{+},\mathrm{\Delta }{v}_m^{-}\ge 0,m\in F,}\\ \multicolumn{1}{c}{v_{\mathrm{mm}\prime }^{+},v_{\mathrm{mm}\prime }^{-}\ge 0,m\in F,m\prime \in F\backslash G(m).}\end{array}$$

Algorithm 2: LP-based decoding algorithm.

Input: $(\mathbf{R},\mathbf{Z},\mathbf{A},\mathbf{B})$ ; Output: $(\mathit{x},\mathit{y})$ ;Set $g_{\mathrm{mm}\prime },\forall m,m\prime =1,2,\dots ,\overline{M}+\overline{S}$ based on $\mathbf{R}$ ; $G(m)\leftarrow Ø,\forall m=1,2,\dots ,\overline{M}+\overline{S}$ ; for $m=1$ to $\overline{M}+\overline{S}$ do  $i_1\leftarrow \mathrm{A}(m)$ ; for $m\prime =m+1$ to $\overline{M}+\overline{S}$ do   $i_2\leftarrow \mathrm{A}(m\prime )$ ;  if ${\mathrm{B}}^{-1}(i_2)>{\mathrm{B}}^{-1}(i_1)$ then    $G(i)\leftarrow G(i)\cup \{i_2\}$ ; $U\leftarrow W$ ; $\mathit{h}\leftarrow \mathit{w}$ ;Obtain $\mathit{v}$ by solving Problem P1; $\mathit{x}\leftarrow \mathit{v}$ ; $G(m)\leftarrow Ø,\forall m=1,2,\dots ,\overline{M}+\overline{S}$ ; for $m=\overline{M}+\overline{S}$ down to $1$ do   $i_1\leftarrow \mathrm{A}(m)$ ; for $m\prime =m-1$ down to $1$ do  $i_2\leftarrow \mathrm{A}(m\prime )$ ; if ${\mathrm{B}}^{-1}(i_2)>{\mathrm{B}}^{-1}(i_1)$ then   $G(i)\leftarrow G(i)\cup \{i_2\}$ ; $U\leftarrow L$ ; $\mathit{h}\leftarrow \mathit{l}$ ;Obtain $\mathbf{v}$ by solving Problem P1; $\mathit{y}\leftarrow \mathit{v}$ ;

Population-based simulated annealing

The PSA algorithm starts with a randomly generated population of size $N^{\mathrm{Pop}}$ , at an initial temperature that is calculated through equation (24).

$${\theta }_0=\frac{{\sum }_{i=1}^{100}|\Psi ({\omega }_i^1)-\Psi ({\omega }_i^2)|}{100\times \ln (0.5)}+1,$$

(24)

where, ${\omega }_i^1$ and ${\omega }_i^2$ are two solutions randomly generated at $i$ th iteration, and $\mathrm{\Psi }(\omega )$ is the fitness value of solution $\omega$ , which is obtained by equation (25).

$$\begin{array}{l}\mathrm{\Psi }(\omega )=\rho \frac{H}{H_L}+\left(1-\rho \right)\frac{\overline{c}}{{\overline{c}}_L}\\ +{\rho }_1\sum _{c\in C\backslash \left\{c_0\right\}}\max \left\{0,\sum _{m\in M}{z}_{mc}-N^M\right\}\\ +{\rho }_2\sum _{m\in F}(\max \left\{0,x_m+w_m/2-W\right\}\\ +\max \left\{0,y_m+l_m/2-L\right\}),\end{array}$$

(25)

where, ${\rho }_1$ is the penalty cost of violating constraint sets (6) and ${\rho }_2$ is the penalty cost of violating constraints sets (10) and (11).

To create a new population, chromosomes are selected using the roulette wheel procedure, with the fitness scaling function given in equation (26).

$$\exp (-\frac{5\times \mathrm{\Psi }({\omega }_i)}{{\mathrm{\Psi }}^{max}}),$$

(26)

where, ${\mathrm{\Psi }}^{max}$ is the highest fitness value in the population.

To speed up the PSA, Algorithm 1 is used during the solve process. Once a chromosome is selected, one of the mutation operators, illustrated in Figure 2, is applied randomly to create a new chromosome. The new chromosome is accepted if it leads to a better solution. Otherwise, it is accepted with the probability of $e^{-\mathrm{\Delta }/{\theta }_k}$ , where $\mathrm{\Delta }$ is the difference between the fitness value of the new and old chromosomes, and ${\theta }_k$ is the temperature at the current generation. For the selected chromosome, this process is repeated for $N^T$ transitions, and in the end, the resulting chromosome is copied to the new population. We assume that $N^T$ is proportional to the problem size as follows:

$$N^T=N\prime (3\overline{M}+2\overline{S}+\sum _{p\in P}{\overline{O}}_p+\sum _{p\in P}\mathrm{sgn}\{\left|M(p,o)\right|-1\}),$$

(27)

where $N\prime$ is a coefficient to be set.

Figure 2.

Example of the mutation operators and the layers names that each operator is applied to: (a) change ( $\mathbf{R}$ and $\mathbf{Z}$ ), (b) swap ( $\mathbf{S}$ , $\mathbf{Z}$ , $\mathbf{A}$ , and $\mathbf{B}$ ), (c) shift ( $\mathbf{S}$ , $\mathbf{A}$ , and $\mathbf{B}$ ), (d) double swap ( $\mathbf{A}$ and $\mathbf{B}$ ), and (e) double shift ( $\mathbf{A}$ and $\mathbf{B}$ ).

After, generating a new population, the temperature is reduced gradually with a cooling rate that is calculated by equation (28).

$$\mu ={\left(\frac{0.001}{{\theta }_0}\right)}^{\frac{1}{I}}.$$

(28)

This process is repeated for $I$ generations. Finally, the population is re-evaluated by changing the decoding algorithm to Algorithm 2. The flowchart describing the major steps in the PSA algorithm is shown in Figure 3.

Figure 3.

Flowchart of the proposed solution algorithm.

Parameters setting

The PSA was coded using Embarcadero® Delphi 10.3 programming software and CPLEX 12.9 callable library (cplex1290.dll) as a computer program. Based on an initial test, we proposed four different levels for each parameter of the PSA. A Taguchi test with L16 orthogonal array is performed on three instances of different sizes to choose the best parameter combination, see Table 2. Figure 4 shows the effect of the PSA parameters on the solution quality. In this figure, the level selected for each parameter is highlighted in red. Based on these results, the following parameter combination is used in the PSA: $N^{\mathrm{Pop}}=20$ , $I=30$ , and $N\prime =30$ .

Table 2.

Combination of the PSA parameters based on L16 orthogonal array.

Run	$N^{\mathrm{Pop}}$	$I$	$N\prime$	Run	$N^{\mathrm{Pop}}$	$I$	$N\prime$
1	15	20	15	9	25	20	25
2	15	25	20	10	25	25	30
3	15	30	25	11	25	30	15
4	15	35	30	12	25	35	20
5	20	20	20	13	30	20	30
6	20	25	15	14	30	25	25
7	20	30	30	15	30	30	20
8	20	35	25	16	30	35	15

Figure 4.

Effect of the PSA parameters on the solution quality.

Case study

The practical usability of the developed model is demonstrated in a case study selected from Bansee and Chowdary. ⁴⁰ In the selected case, a factory manufactures 20 parts using 15 machines arranged according to a job shop configuration in six departments: Lathe, Drills, Milling, Boring, Grinding, and Shaper. We assume that the final product is made according to the assembly diagram shown in Figure 5. As can be seen in this figure, the assembly process is done in 3 workstations through fulfilling 11 tasks. It is assumed that every task takes 200 s to be completed. The transfer time of parts within the cells and among the cells are respectively assumed to be 10 and 15 s per meter, and the handling cost of the part within the cell and among the cells are respectively set to 1 and 1.5 units of cost per meter. Figure 6(a) shows a schematic view of the current layout in the factory. To be able to compare our solution against Bansee and Chowdary’s ⁴⁰ solutions, we construct a more detailed layout that is shown in Figure 6(b).

Figure 5.

Diagram of assembling parts $\{p_1,p_2,\dots ,p_{20}\}$ into a product by performing tasks $\{t_1,t_2,\dots ,t_{11}\}$ in workstations $\{s_1,s_2,s_3\}$ .

Figure 6.

Layout constructed for the case study based on the schematic layout presented in Bansee and Chowdary ⁴⁰ : (a) schematic layout of the machines and (b) extended layout considering the workstations, dimensions of facilities, and clearance among facilities.

With this information, we solve the problem using the PSA with 100 replications, and during the solving process, we set $n=10$ . The layout corresponding to the best solution is depicted in Figure 7. Also, for the first production cycle, the scheduling Gantt chart corresponding to our solution and Bansee and Chowdary’s ⁴⁰ solution is given in Figure 8.

Figure 7.

Solution obtained for the case study (machines $m_4$ , $m_9$ , and $m_{13}$ are free).

Figure 8.

Gantt charts corresponding to the case study drawn at $n=1$ : (a) Gantt chart of the solution presented in Bansee and Chowdary ⁴⁰ and (b) Gantt chart of the solution obtained in this study.

Based on Bansee and Chowdary’s ⁴⁰ solution, we calculated the handling cost as 833.25 units, and the cycle time equal to 3101 min. While, for our solution, the handling cost is 349.125 units and the cycle time is 3024 min. This shows 58.1% and 2.48% reduction in the handling cost and cycle time, respectively. From Figure 8(a), it is seen that there are fifteen inter-cell moves in Bansee and Chowdary’s ⁴⁰ solution, while there is only one inter-cell move in Figure 8(b). This can justify why our solution is much better than theirs in terms of the handling cost. Furthermore, in Figure 8 it is seen that the transportation times in our solution are much shorter than those in Bansee and Chowdary’s ⁴⁰ solution. This is due to the more efficient process plans and better arrangement of facilities in our solution. Another important finding is that in our solution, machines $m_4$ , $m_9$ , and $m_{13}$ are free, while in Bansee and Chowdary’s ⁴⁰ solution all of the machines are utilized. This demonstrates a more effective machine utilization in our solution, as with less resources, a shorter cycle time is achieved.

Comparison study

Instances selected from the literature are solved to access the performance of the PSA against CPLEX solver and standard SA. Table 3 shows the specification of the instances.

Table 3.

List of the instances selected from the literature.

Instance		Size						Objective functions lower bound
#	Source	$\overline{P}$	$\overline{M}$	$\overline{T}$	$\overline{S}$	$\overline{C}$	$N^M$	$H_L$	${\overline{c}}_L$
1	Feng et al. 30	10	10	7	3	2	5	30,004.75	323.75
2	Feng et al. 30	10	10	7	3	3	4	30,775.25	323.55
3	Jayakumar and Raju 41	12	1	8	4	3	5	139,499.625	1551.55
4	Jayakumar and Raju 41	12	10	8	4	4	4	142,678.125	1565.35
5	Jayakumar and Raju 41	12	15	8	4	3	5	137,837.625	1054.2
6	Jayakumar and Raju 41	12	15	8	4	4	4	145,374.125	1069.2
7	Ramabhatta and Nagi 42	15	15	9	4	3	6	9415	2.03975
8	Ramabhatta and Nagi 42	15	15	9	4	4	4	10,039	2.11
9	Feng et al. 21	15	15	7	3	3	6	124,825.75	46,479.5
10	Feng et al. 21	15	15	7	3	4	4	127,539.75	46,616
11	Chang et al. 43	20	14	10	4	3	6	59,083	822.5
12	Chang et al. 43	20	14	10	4	4	5	61,865	823.3
13	Bansee and Chowdary 40	20	15	11	3	3	5	347.75	3008
14	Bansee and Chowdary 40	20	15	11	3	4	4	353	3008
15	Nagi et al. 44	20	20	10	5	3	7	5935.5	1.615
16	Nagi et al. 44	20	20	10	5	5	5	5964	1.623
17	Chang et al. 43	30	18	13	5	3	7	89,296.5	916.4
18	Chang et al. 43	30	18	13	5	4	5	94,105	916.15

We coded the developed mathematical model in GAMS Studio 0.11.2 modeling software. The time limit for solving each instance is set to 7200 s and CPLEX is chosen as the solver. To make a more precise comparison, each instance is solved 30 times by the SA and PSA. It should be mentioned that all computations are performed on a PC, with an Intel Core i7 4790 k 4.00 GHz CPU and 16 GB of RAM. Table 4 shows a summary of the comparison results; in this table, $F$ is the best objective value of each solution approach, $\overline{t}$ is the average computing time, opt. gap is the relative optimality gap obtained by CPLEX, $F^B$ is the overall best objective function value obtained by any solution approach, and $d$ is the relative difference between $F$ and $F^B$ . Figure 9 also shows the performance profiles ⁴⁵ plotted based on the best objective function values.

Table 4.

Summary of the computational results.

Instance	SA			PSA			CPLEX			$F^B$
#	F	$\overline{t}$	d	F	$\overline{t}$	d	$F$	Optimality gap	d
1	1.091	4.656	1.260	1.077	13.200	0.000	1.298	0.578	20.498	1.077
2	1.054	4.684	0.097	1.053	13.152	0.000	1.598	0.721	51.748	1.053
3	1.138	4.782	1.912	1.117	16.028	0.000	1.408	0.692	26.042	1.117
4 1	1.106	4.749	1.233	1.093	15.974	0.000	1.342	0.676	22.815	1.093
5	1.186	5.205	5.586	1.124	21.290	0.000	Infeasible	–	–	1.124
6	1.164	5.177	0.000	1.172	20.463	0.658	Infeasible	–	–	1.164
7	1.034	10.795	0.000	1.037	60.241	0.314	Infeasible	–	–	1.034
8	1.050	10.903	0.446	1.045	59.921	0.000	Infeasible	–	–	1.045
9	1.097	7.271	0.000	1.105	33.672	0.729	1.569	0.689	42.930	1.097
10	1.118	7.400	1.161	1.105	33.621	0.000	2.020	0.757	82.749	1.105
11	1.012	9.275	0.201	1.010	46.402	0.000	1.726	0.746	70.973	1.010
12	1.012	9.287	0.000	1.019	46.106	0.697	1.364	0.672	34.861	1.012
13	1.015	6.125	0.176	1.013	24.569	0.000	Infeasible	–	–	1.013
14	1.022	6.065	0.763	1.014	24.697	0.000	Infeasible	–	–	1.014
15	1.069	8.828	2.343	1.045	43.067	0.000	Infeasible	–	–	1.045
16	1.039	8.827	1.751	1.021	42.759	0.000	Infeasible	–	–	1.021
17	1.013	12.578	0.000	1.025	87.754	1.160	Infeasible	–	–	1.013
18	1.023	12.414	0.394	1.019	87.046	0.000	Infeasible	–	–	1.019

$d=100(F-F^B)/F^B$ .

1

In this case, CPLEX stopped solving after 7186 s due to insufficient physical memory.

Figure 9.

Performance profiles plotted based on the best objective function values of the algorithms.

CPLEX could not optimally solve any of the instances within the time limit of 7200 s. According to Table 4, out of 18 instances, in 10 instances, CPLEX could not even find a feasible solution in the given time limit. It should be mentioned that CPLEX stopped solving instance #4 after 7186 s due to the lack of physical memory. For CPLEX solutions, it is also seen that the relative deviation from the overall best objective function value, that is, $d$ , is too high. On the other hand, the results of the algorithms show that both the SA and PSA can solve the instances with much better objective function values, in reasonably less computational time. Out of 18 instances, in 13 instances, the overall best solution has been obtained by the PSA, while for the SA, this value is 5. However, it should be noted that the PSA takes more computing time than the SA.

We carried out a Wilcoxon Signed-Rank test on the best objective function values using the Minitab® 19 software to determine if the PSA is meaningfully better than the SA or not. Based on this test, the Wilcoxon statistic was calculated as 131 with a $p$ -value equal to 0.025. So, it is concluded that in terms of the solution quality, the PSA is better than the SA at the 95% confidence level. This conclusion is also supported by the performance profiles plotted in Figure 9. As this figure shows, compared to the SA, the PSA gives better performance at any performance ratio.

Conclusion

In this study, we addressed the joint problem of cell formation, facility layout, and scheduling in virtual manufacturing cells considering assembly aspects and process routing. The problem was formulated as a mixed-integer program with the aim of minimizing the handling costs and cycle time. A Population-based Simulated Annealing (PSA) was applied to solve the problem. The practical usability of the developed model was demonstrated in a case study. The result of case study showed that the by applying the developed joint problem, the factory can reduce the handling costs and cycle time by 58.1% and 2.48%, respectively. Finally, instances from the literature were solved to evaluate the performance of the PSA. The comparison results showed the superior performance of the PSA in comparison with CPLEX solver and standard simulated annealing.