Capacity Scalability Planning Algorithms for Job-shop-type Reconfigurable Manufacturing Systems with Dynamic Demands

2.1 Problem description

The system considered in this study is a job-shop-type RMS that consists of computer numerical control machines, loading/unloading (L/U) stations, material transporters and a central buffer. Note that the RMS is fundamentally different from the conventional FMSs in that it has the capability to change system components quickly. Figure 1 shows an RMS with three machines and two L/U stations. A system-level reconfiguration from two to three machines is also illustrated in this figure.

Figure 1.

Reconfigurable Manufacturing System with System-level Reconfiguration: Example

In the RMS, each part is processed through one or more machines in non-unidirectional flows after loaded on a pallet at an L/U station and then released into the central buffer. Similarly, a completed part is unloaded from the pallet at an L/U station and then exits from the RMS. The central buffer, which is an automatic storage/retrieval system with a limited number of storage locations, is used to store the pallets. Each machine can process parts using the required cutting tools stored in its tool magazine with a sufficient tool slot capacity. Parts are produced by one or more operations with precedence relations, each of which can be processed on one of the machines at the pre-specified processing workstation. Finally, one or more automated transporters are used to move pallets among the system components, i.e. L/U station, processing workstations and central buffer.

For given non-decreasing demands of multiple product types over a planning horizon, the basic case is to determine the number of additional components, i.e. machines, transporters, L/U stations and pallets, to satisfy the demands in each period of the planning horizon for the objective of minimizing the sum of component acquisition and configuration change costs. The component acquisition costs are those required to acquire components, and the configuration change costs are those required to install newly acquired components. Due to non-decreasing demands, the system components added in a period are maintained during the remaining periods. Without loss of generality, it is assumed that there are no system components at the beginning of the planning horizon. Besides the demand requirements, other constraints are the limited number of pallets and the minimum allowable station utilization. The number of pallets is limited due to the central buffer capacity, which a practical extension of Yu et al. (2014). Also, the minimum allowable utilization constraint implies that utilizations of processing workstations and L/U stations must not be less than a lower limit, which is needed to avoid over-installations of unnecessary components.

This study considers a deterministic version of the problem, i.e. all data such as process plans, demands and cost values are deterministic and given in advance. Among the data, the process plans, which contain the information on operations, processing workstations and processing/transportation times, is needed to estimate throughputs and station utilizations under a system configuration.

The problem can be formulated as a nonlinear integer programming model, which modifies Yu et al. (2014)’s by adding the limited number of pallets. The notations used are summarized below.

Indices
i part type, i=1, 2, ..., I
t periods, t=1, 2, ..., T
m stations, m=1, 2, ..., M (1, 2, ... M-2 : processing workstations; M-1: L/U station; and M: transportation station)

Parameters
a_mt acquisition cost of a component at station m in period t
c_mt configuration change cost of station m in period t
h acquisition cost of a pallet
d_it demand of part type i in period t
u_min minimum allowable station utilization
q_max maximum number of pallets
G large number

Decision variables
u_mt number of newly acquired components at station m in period t
y_mt = 1 if there is a configuration change at station m in period t, and 0 otherwise
z_t: number of newly acquired pallets in period t

Now, the nonlinear integer programming model is given below. In the model, TH_i(X_t, p_t) and UT_m(X_t, p_t) denote the throughput of part type i and the utilization of station m under configuration (X_t, p_t) in period t, where X_t=(x_1t, x_2t, ..., x_Mt) and $$ p_t=\sum _{l=1}^t z_l$$. Note that $$ x_{mt}=\sum _{l=1}^t u_{ml}$$ and p_t denote the numbers of components at station m and pallets in period t, respectively.

[P-NID] Minimize
$$ \sum _{t=1}^T \left[\sum _{m=1}^M \left(a_{mt}\cdot u_{mt}+c_{mt}\cdot y_{mt}\right)+h\cdot z_t\right]$$
subject to

The objective function denotes the sum of component acquisition and configuration change costs over the planning horizon. Constraints (1) and (2) represent the demand and the minimum allowable station utilization requirements, respectively. In this study, part throughputs and station utilizations are estimated using the closed queuing network model of Yu et al. (2014), originally proposed by Solberg (1977) for performance evaluation of FMSs. The detailed estimation methods will be explained in the next section. Constraint (3) ensures that a component acquisition occurs in a period if there is a configuration change in that period. Constraints (4) and (5) represent the limited number of pallets and there are no system components at the beginning of the planning horizon, respectively. Finally, the remaining constraints represent the conditions of decision variables.

It can be easily seen that the problem [P-NID] is difficult to solve optimally for large-sized instances due to the exponential number of possible configurations and the nonlinear closed queuing network based estimation functions TH_i(X_t, p_t) and UT_m(X_t, p_t). Therefore, instead of the optimal approach with very limited practical applications, the heuristic approach is adopted in this study.

2.2 Solution algorithms

This section presents the heuristic algorithms proposed in this study. Before explaining the heuristics, the closed queuing network (CQN) model is briefly explained.

(1) Estimating throughputs and utilizations

The CQN model represents an RMS configuration (X_t, p_t) as M inter-connected stations with one or more identical components. <Figure 2> shows the CQN model with processing stations (1, 2, ..., M-2), an L/U station (M-1) and a transportation station (M), where the transportation station takes over the role of a central server activated after each individual processing operation is finished on a machine at a processing station. The assumptions are: (a) exponential processing/transportation times; (b) new parts available at all times; (c) first come first served queueing discipline; (d) universal pallets that can load all part types being processed in the RMS; and (e) sufficient central buffer capacity. See Solberg (1977) for more details on the CQN topology.

Figure 2.

Closed Queuing Network (CQN) Model for an RMS Configuration

For a given RMS configuration (X_t, p_t), let n_mt be the number of pallets located at station m in period t and hence $$ p_t=\sum _{m=1}^M n_{mt}$$. Also, let w_mt be the average workload at station m in period t, which can be calculated using the process plans. Then, the probability that the RMS is in state n_t=(n_1t, n_2t, ..., n_Mt) can be represented as

$$ P\left({\mathit{n}}_t\right)=P\left(n_{1t},n_{2t},\dots ,n_{Mt}\right)=\frac{1}{g\left(p_t,M\right)}{\mathrm{\Pi }}_{m=1}^M{f}_m\left(n_{mt}\right),$$

where $$ f_m\left(n_{mt}\right)=w_{mt}^{n_{mt}}/n_{mt}!$$ if n_mt≤x_mt and $$ w_{mt}^{n_{mt}}/x_{mt}!\cdot {\left(x_{mt}\right)}^{n_{mt-x_{mt}}}$$ otherwise. Note that g(p_t, M) is the normalization constants over the state space, which can be calculated by the Buzen algorithm. See Tempelmeier and Kuhn (1993) for more details on the Buzen algorithm.

Then, the throughput of part type i can be estimated as

$$ T{H}_i\left({\mathbf{X}}_t,p_t\right)={\alpha }_i\cdot \frac{\left[g\left(p_t-1,M\right)/g\left(p_t,M\right)\right]}{v_M}$$,

where α_i and v_M denote the production ratio of part type i and the average number of operations per part, respectively. Also, the utilization of station m can be estimated as

$$ U{T}_m\left({\mathbf{X}}_t,p_t\right)=b_m\cdot \frac{r_m\cdot \left[g\left(p_t-1,M\right)/g\left(p_t,M\right)\right]}{x_m}$$,

where b_m and r_m denote the average processing time per operation and the relative arrival frequency (visit ratio) of a part at station, respectively. See Tempelmeier and Kuhn (1993) for more details on the derivations of the above two formulas.

2.3 Computational Results

To test the performance of the heuristics proposed in this study, computational experiments were done and the results are reported in this section. Specifically, the two heuristics, MB-UT and MB-TH, are compared with the forward-throughput (F-TH), forward-utilization (F-UT) and backward-utilization (B-UT) heuristics of Yu et al. (2014) after setting the limited number of pallets. The heuristics were coded in Python 3.1 and the experiments were done on a PC with an Intel core i9 processor at 3.42 GHz clock speed with 64GB RAM memory.

The first experiment was done to show the absolute performance for small sized test instances. For the test, 30 instances with three periods and five stations were generated randomly, i.e. 10 instances for each of three levels for the number of part types (10, 20 and 30) using the data of Yu et al. (2014) and Zhou et al. (2014). The data were generated as follows.

• •Demands in each period (d_it) ~ DU(50,70)+DU(0,10), where DU(0,10) represents the amount of an additional demand
• •Component acquisition costs (c_mt) ~ DU(5000,20000)
• •Configuration changing costs (c_mt) ~ DU(500,2000)
• •Pallet costs (h) ~ DU(200,300)
• •Number of operations for a part type ~ DU(6,16)
• •Operation processing times ~ DU(20,100)
• •Loading/unloading times ~ DU(5,20)
• •Transportation times ~ DU(3,7)

In addition, the limited number of pallets (q_max) was set to 60. The performance measures are: (a) percentage gaps from the optimal solution values that obtained using the full enumeration method and (b) CPU seconds.

Test results are summarized in <Table 1> (a), (b) and (c) that show the average percentage gaps from optimal solution values for three levels of the minimum allowable station utilization (0.6, 0.7 and 0.8). It can be seen from the tables that the heuristics proposed in this study outperform the previous ones in overall averages due to the better last period configurations. Of the heuristics, MB-TH was significantly better than MB-UT because it considers both throughputs and costs. In fact, the overall average gaps of MB-TH were 2.47%, 2.9% and 2.11% when the minimum allowable station utilization was 0.6, 0.7 and 0.8, respectively. To test statistical difference among the heuristics, paired t-tests were done using the average percentage gaps and the results showed that MB-TH is statistically different from the others in the significance level of 0.01. Finally, all the heuristics solve the test instances within 8.2 seconds while the full enumeration method required 2108.0 seconds in overall average.

The second experiment was done for medium-to-large sized test instances. For the test, 240 instances, i.e. 10 instances for each of 24 combinations of three levels for the number of periods (3, 5 and 10), three levels for the number of stations (5, 7 and 9) and three levels for the number of part types (10, 20 and 30), were generated using the data generation method explained earlier. The limited numbers of pallets were set to 60, 80 and 100, and the minimum allowable station utilizations were set to 0.6, 0.7 and 0.8 for the test instances with 5, 7 and 9 stations, respectively. The performance measures are: (a) relative performance ratio (RPR) since the optimal solutions could not be obtained; and (b) CPU seconds, where the RPR of a heuristic for an instance is calculated as

$$ 100\cdot \left(C_a-C_{best}\right)/C_{best }$$,

where C_a is the total cost value obtained from heuristic a and C_best is the best one among those obtained from all the heuristics.

Test results are summarized in <Tables 2> (a), (b) and (c) that show the average RPR values of the heuristics for the three levels of the minimum allowable station utilization. As in the test results for the small sized instances, the new heuristics perform better than the previous ones in overall averages and MB-TH performs the best. Paired t-tests were also done and the results showed that MB-TH is statistically better than the others in the significance level of 0.01, which shows the effectiveness of the throughput/cost ratio based local search that determines the better last period configurations. Finally, all the heuristics gave the solutions within 24.5 seconds.

3.1 Problem Description

The general case is the same as the basic one except that the demands are fluctuating, i.e. increasing or decreasing, over a planning horizon. Therefore, the general case can be defined as follows. For given fluctuating demands of multiple product types, the problem is to determine the number of system components to satisfy the demands in each period of a planning horizon for the objective of minimizing the sum of component acquisition/removal and configuration change costs. As in the basic case, the constraints are demand requirements, the limited number of pallets and the minimum allowable station utilization. It is assumed that the pallets are acquired at the beginning of the planning horizon.

Let s_mt and g_mt denote the removal cost of a component and the number of removed components at station m in period t, respectively. Then, the general case can be formulated as the following nonlinear integer programming model.

[P-FD] Minimize
$$ \begin{array}{c}\sum _{m=1}^M \sum _{t=1}^T \left(a_{mt}\cdot u_{mt}+c_{mt}\cdot y_{mt}+s_{mt}\cdot g_{mt}\right)\\ +h\cdot {\mathrm{m}\mathrm{a}\mathrm{x}}_{\left(t=1,\dots ,T\right)} \left\{p_t\right\}\end{array}$$

subject to

The objective function represents the sum of component acquisition/removal and configuration change costs over the planning horizon. Constraint (3') ensures that a component acquisition/removal occurs in a period if there is a configuration change in that period. Constraints (4') and (5') represent the limited number of pallets and there are no system components at the beginning of the planning horizon, respectively. Constraint (6') ensures that acquisition and removal of the same component cannot occur at the same time in a period. Finally, the remaining constraints (7') and (8) represent the conditions of decision variables.

It can be easily seen that the problem [P-FD] is harder than [P-NID] because components can be removed according to fluctuating demands and hence the solution space gets much larger. Therefore, instead of simple local search heuristics, this study adopts the meta-heuristic approach.

3.2 Solution Algorithms

This section explains the two variable neighborhood search algorithms proposed in this study, i.e. ordinary and hybrid ones. Variable neighborhood search (VNS), proposed by Mladenović and Hansen (1997), is a meta-heuristic that explores a set of predefined neighborhood structures successively to escape from local optimums. See Hansen et al. (2017) for more details on the VNS algorithm.

(1) Ordinary VNS algorithm

An overview of the ordinary VNS algorithm is shown in Figure 3. As can be seen in the figure, the algorithm consists of four main steps: (a) obtaining an initial solution; (b) generating shaking solutions using different neighborhood structures; (c) local search that improves the shaking solution and (d) termination. In the algorithm, a solution is represented by an M×T matrix X=[x_mt], where x_mt denotes the number of components at station m in period t. Therefore, the t th column vector X_t=(x_1t, x_2t, ..., x_Mt)^T represents the configuration with p_t pallets in period t. The details of each step are explained below.

Figure 3.

Ordinary VNS Algorithm: Overview

Step 1. Obtaining an initial solution

An initial solution, denoted as $$ {\mathrm{X}}^i=\left[x_{mt}^i\right]$$, is obtained by determining the feasible configuration in each period using the first step of the MB-UT heuristic, i.e. Step 1 of procedure 1.

Step 2. Shaking

Shaking is done to escape from the local optimums by changing neighborhood structures. More specifically, from the current incumbent solution X, a shaking solution is obtained as follows. First, a neighborhood structure N^' (X) is selected randomly. Then, a feasible solution X^f that satisfies the demand and the minimum allowable utilization constraints is generated randomly using the selected neighborhood structure N^' (X), where X^f∈N^' (X). Finally, a feasible shaking solution X^' is obtained by decreasing the number of pallets in X^f to the minimum possible one.

In this study, the following neighborhood structures are proposed.

One component change in multiple periods (OCC-MP). This method generates a feasible neighborhood solution by selecting P_{OCC_MP} different periods randomly and then changing the number of components at a random station in each of the periods in such a way that the changed number does not exceed the maximum number of components in the initial solution, i.e. $$ {\mathrm{m}\mathrm{a}\mathrm{x}}_{m,t}\left\{x_{mt}^i\right\}$$. <Figure 4 (a)> shows an example with P_{OCC_MP}=2, in which periods 2 and 3 are selected and then the number of components at station 3 (2) in period 2 (3) is changed from 2 (4) to 3 (2).

Figure 4.

Neighborhood Structures: Examples

Multiple component changes in one period (MCC-OP). This method generates a feasible neighborhood solution by selecting a period randomly and then adding or removing one component randomly at each of C_{MCC_OP} random stations in that period. <Figure 4(b)> shows an example with C_{MCC_OP}, in which period 2 is selected and then one component is removed at station 1 while added at stations 3 and 4.

Multiple component replacements in one period (MCR-OP). This method generates a feasible neighborhood solution by selecting a period randomly and then replacing its number of components at each of C_{MCR_OP} random stations by that of the random adjacent period. <Figure 4 (c)> shows an example with C_{MCR_OP}=2, in which period 2 is selected and then the number of components at station 2 (4) is replaced by 4 (3) in the adjacent period 3.

Step 3. Local search improvement

This step improves the shaking solutions by the neighborhood structures explained earlier. Specifically, a neighborhood solution is generated from the current shaking solution X^' using each of OCC_MP, MCC_OP and MCR_OP methods, and the best one X^'' is selected. If X^'' improves the current best solution, the best solution is updated and the next shaking is done. Otherwise, the local search is done for another shaking solution generated by an unused neighborhood structure. Note that the iteration count is reinitialized to 0 if an improved solution is obtained by the consecutive shaking and local search improvement steps, while increased by 1, otherwise.

Step 4. Termination

The ordinary VNS algorithm is terminated when there is no improvement for a certain number P_{O_VNS} of consecutive iterations.

3.3 Computational Results

To test the performance of the VNS algorithms, computational experiments were done and the results are reported in this section. The VNS algorithms were coded in Python 3.1 and the tests were done on a PC with an Intel core i9 processor at 3.42 GHz clock speed with 64GB RAM memory.

Before the main tests, a preliminary test was done to set the parameters of the hybrid VNS algorithm. Specifically, the best one was selected using RPR values after comparing six combinations of three levels for cooling ratio (0.6, 0.7 and 0.8) and two levels for epoch length (6 and 8) when the initial temperature and the termination parameter were set 10000 and 30, respectively. For the test, 10 small sized (3 periods, 5 stations and 10 part types), 10 medium sized (5 periods, 7 stations and 20 part types)and 10 large sized instances (10 periods, 9 stations and 30 part types) were generated randomly. The data were generated as follows.

• •Demands in each period (d_it) ~ DU(10,150)
• •Component acquisition costs (a_mt) ~ DU(10000,20000)
• •Component removal costs (s_mt) ~ DU(5000,10000)
• •Pallet costs (h) ~ DU(200,300)
• •Configuration changing costs (c_mt) ~ DU(500,2000)
• •Number of operations for a part type ~ DU(6,16)
• •Operation processing times ~ DU(20,100)
• •Loading/unloading times ~ DU(5,20)
• •Transportation times ~ DU(3,7)

In addition, the limited numbers of pallets (q_max) were set to 60, 80 and 100 for the instances with 5, 7 and 9 stations, respectively, and the minimum allowable station utilization was set to 0.7. Test results are summarized in Table 3 that shows the average RPR values for all parameter combinations. From the test results, the cooling ratio (β) and the epoch length (Γ) were set to 0.7 and 8, respectively.

Table 3.

Test results for Setting Parameters of hybrid VNS Algorithm

The first test was done to show the absolute performance of the VNS algorithms for small sized test instances. The performance measures are: (a) percentage gaps from the optimal solution values; and (b) CPU seconds, where the optimal solutions were obtained using the full enumeration method. For the test, 30 instances with 3 periods and 5 stations were generated randomly, i.e. 10 instances for each of three levels for the number of part types (10, 20 and 30). The data were generated using the method for the preliminary test and the algorithms were terminated when there was no improvement for 30 consecutive iterations, i.e. P_{O_VNS}=P_{H_VNS}=30.

Test results are summarized in <Table 4(a), (b) and (c)> that show the average percentage gaps from the optimal solution values for the three levels of the minimum allowable station utilization (0.6, 0.7 and 0.8). It can be seen from the tables that the VNS algorithms give optimal or near optimal solutions for all test instances. Of the two algorithms, the hybrid VNS algorithm outperformed the other in overall average and gave more optimal solutions as the minimum allowable station utilization increases. A paired t-test was done and the results showed that the hybrid algorithm is statistically better than the other in the significance level of 0.01. The overall average gaps of the hybrid algorithm were 0.09%, 0.06% and 0.05% for utilization levels 0.6, 0.7 and 0.8, respectively. Also, it gave the 77 optimal solutions out of 90 test instances. This implies that the SA based solution acceptance criterion is an effective method to escape from the local optimums in the local search improvement step. Finally, as can be seen in <Table 5>, the hybrid VNS algorithm required more computation times, but gave the solutions within 854.2 seconds. Note that the full enumeration method required 3203.4 seconds in the overall average.

Table 5.

General Case: CPU Seconds of VNS Algorithms for Small Sized Test Instances

The second main test was done to compare the relative performance of the two VNS algorithms for medium-to-large sized test instances. For the test, 240 instances, i.e. 10 instances for each of 24 combinations of three levels for the number of periods (3, 5 and 10), three levels for the number of stations (5, 7 and 9) and three levels for the number of part types (10, 20 and 30), were generated using the method explained earlier. <Table 6(a), (b) and (c)> summarize the average, minimum and maximum RPR values of the two algorithms for the three levels of the minimum allowable station utilization. As in the results for small sized test instances, the hybrid algorithm outperformed the ordinary one in overall averages and gave smaller RPR values as the minimum allowable station utilization increases, which also implies that the SA based solution acceptance is an effective method. In fact, the result of paired t-test showed that the two algorithms are statistically different in the significance level 0.01. Finally, as can be seen in <Table 7>, there was not much difference in computation times. In fact, the hybrid VNS algorithm required 1443.3 seconds in overall average, which is acceptable for practical applications because capacity scalability planning is a system design problem.

Table 7.

General Case: CPU Seconds of VNS Algorithms for Medium-to-Large Sized Test Instances