**A two-stage stochastic model for designing cellular manufacturing systems with simultaneous multiple processing routes and subcontracting** **Abstract** In recent decades, many researchers have studied the cellular manufacturing system with consideration of various issues such as scheduling, production planning, layout, reliability, etc. However, limited research papers have investigated this problem in an uncertain environment. The present paper addresses a stochastic problem in cellular manufacturing systems considering simultaneous multiple routings and subcontracting. In the developed problem, each part can be simultaneously produced in multiple processing routes. It is also assumed that the unsatisfied part demands as a result of limited machine capacity or high manufacturing cost could be outsourced. A two-stage stochastic programming approach is employed to take the uncertainty into consideration and to formulate the problem. The objective function is to minimize the summation of production, subcontracting, material handling, and machine idleness costs. A sample average approximation method is applied to solve the problem. We also use a numerical example and accomplish some important sensitivity analyses to clarify the problem. Finally, through some numerical examples extracted from related literature, the advantages of constructing a stochastic optimization model for the problem are demonstrated. **Keywords:** Cellular manufacturing system*; *Cell formation; Stochastic programming; Sample average approximation; Multiple processing routes; Subcontracting; **1. Introduction** Cellular Manufacturing system (CMS), as a manufacturing strategy, emanates from group technology concept. In CMSs, each cell is composed of a group of machines that are dedicated to the production of a specific subset of parts called part family. The main advantages expected from production using CMS include a reduction in work-in-process inventories, setup times, lead times, material handling costs, and tool requirements. Production using CMS also cause a noticeable improvement in productivity, product quality, and production control [1,2]. According to what has been presented by Wemmerlöv and Hyer [1], the following decisions could be taken in the design stage of CMSs:

- Cell Formation (CF): grouping parts with similar processing requirements or design features (i.e., making part families) and categorizing machines into machine cells based on the needed operations for part families such that the inter-cell movement of parts is minimized.
- Group layout: designing intra-cell layout (layout of machines inside cells) and inter-cell layout (layout of cells regarding each other).
- Group scheduling: scheduling parts and part families for production.
- Resource allocation: assigning tools, manpower, materials, and other resources.

These issues have been investigated in many research papers. For instance, the CF problem has been addressed in [3–6]. The layout of CMS has been discussed in [7–10]. The cell scheduling problem has been considered in [11–14]. In reality, the aforementioned decisions should be taken in an uncertain environment. However, limited studies have investigated the uncertainty issue in designing CMSs. The uncertainty could appear in demands, operational costs, resource capacities, etc. Dealing with this issue, Wang et al. [15] surveyed a facility layout problem in CMSs in which the demand rate varies over product life cycle. In their problem, the objective was to find the inter- and intra-cell layouts such that the total material handling cost is minimized. A simulated annealing algorithm was applied to solve the presented problem. Jeon and Leep [16] developed a two-phase procedure to configure a CMS. In the first phase, they used a genetic algorithm to obtain part families using a new similarity coefficient. The proposed similarity coefficient considers the possibility of alternative processing routes during machine failure. In the second phase, considering the periodic demand changes, machines are assigned to part families using multiple optimization processes. In these optimization processes, sequential and simultaneous mixed integer programming models were employed to minimize the total costs that are in connection with the operational and scheduling aspects. Tavakkoli-Moghaddam et al. [17] surveyed a facility layout problem in CMSs with normally distributed part demands. The objective was to optimize the total inter- and intra-cell material handling costs. Schaller [18] presented a model for the CF problem considering periodic demand variability. In their problem, the composition of cells is allowed to be changed from period to period. Heuristic algorithms were implemented to minimize the sum of three cost components including the production cost of parts, the amortized cost of machines and the relocating cost of machines during periods. Safaei et al. [19] assumed that part demands and machine capacities are fuzzy numbers. Then, they formulated a dynamic CMS model. The objective of this model was to determine the optimal cell configuration in each period in such a way that the sum of material handling, machine constant/variable, and reconfiguration costs is minimized. They used a fuzzy programming based method to solve the developed problem. Arıkan and Güngör [20] presented a multi-objective CMS model in a fuzzy environment. In this study, part demands, machine capacities, and exceptional elements (EEs) elimination costs were considered as fuzzy numbers. The objective functions that they considered include the minimization of EEs elimination cost, the minimization of inter-cell movements, and the maximization of utilized machine capacity. Ghezavati and Saidi-Mehrabad [21] addressed an integrated mathematical model of CF and group scheduling problems in an uncertain environment. It was assumed that the processing time of parts on machines is a stochastic parameter which is represented by discrete scenarios. The main goal of their model was to minimize the total expected costs of maximum tardiness, EEs subcontracting and resource under-utilization. A hybrid genetic-simulated annealing algorithm was employed as a solution method. Das and Abdul-Kader [22] presented a bi-objective integer-programming model for designing a CMS by considering dynamic changes in machine reliability and parts demands. The first objective function was to minimize the total system costs including the manufacturing, inter-cell material handling, machine under-utilization, and machine duplication costs. The second objective function was to maximize the total system reliability. A *ε*-constraint solution method was used to solve the problem. Ghezavati and Saidi-Mehrabad [23] applied a queuing theory approach to design a CMS with exponentially distributed service and arrival times. It was assumed that each machine works as a server and each part is a customer that should be served by machines. They formulated a mathematical model to maximize the average utilization level of machines. A hybrid method based on genetic and simulated annealing algorithms was exerted to solve the problem. Rabbani et al. [24] proposed a bi-objective CF problem in which part demands are expressed by some probabilistic scenarios. A two-stage stochastic programming model was presented to undertake the uncertain demand of parts. The expected variable cost of all machines and the expected inter-cell material handling cost were considered in the first objective function. The total expected cell load variation was considered as the second objective function. They applied a two-phase fuzzy linear programming approach to solve the presented problem. Forghani et al. [25] applied an interval robust optimization approach to take the uncertainty of part demands into consideration. Then, an integrated CF and layout problem was formulated to minimize the inter- and intra-cell material handling costs. Based on the above survey, the following shortcomings in the developed models can be more investigated:

- Using normal probability distributions: In practice, uncertain parameters are not necessarily normally distributed and may follow other distributions. However, for simplicity, most papers assume that random variables are normally distributed.

- Considering single process routing and other simplifying assumptions: Most related papers usually assume that each part has a unique process routing. While, in practice, each part may have multiple processing routes. Consideration of multiple processing routes in the design of CMSs may enhance planning flexibility and throughput rates [26]. Furthermore, other manufacturing parameters such as operation sequences, processing times, resource capacities and operational costs are not usually addressed in the existing mathematical models.

- Limitation in considering machine duplication: In most current studies, the number of each machine type is known a priori, and there isn’t a capacity limitation on machines. In other words, most developed models do not make a decision about the number of machines in the cell design process. In reality, the capacity of machines is limited. So, each machine can be duplicated so as to cope with the production requirement of parts.

- Neglecting outsourcing costs: Producing all parts in CMSs may not be possible due to either limitation in resource capacities or high manufacturing costs. Accordingly, outsourcing some parts to externals suppliers may be preferred under certain circumstances. However, this issue has only addressed in few research papers.

To deal with the shortcomings mentioned above, in this paper, we propose a stochastic CMS model considering subcontracting and simultaneous multiple processing routes. Furthermore, practical design parameters such as operation sequences, processing times, and machines capacities are taken into account. It is assumed that part demands and outsourcing costs are stochastic parameters with known probability distributions. In the suggested problem, each part can be simultaneously produced in multiple processing routes. Unsatisfied part demands as a result of limited machine capacity or high manufacturing cost is outsourced. A two-stage stochastic programming approach is applied to cope with the uncertainty and formulate the problem. By considering the cell size, machine capacity, and budget constraints, the first stage decisions are the assignment of each machine to a cell and the number of purchased machines of each type. The second stage variables, which are dependent on the various realizations of uncertain parameters and the first stage decisions, are the amount of part demands that should be outsourced and produced in each processing route. The objective function is to minimize the total variable cost, including the production, subcontracting, material handling, and machine idleness costs. A solution method based on Sample Average Approximation (SAA) approach is suggested to find suitable solutions for this stochastic model. To illustrate the problem, a numerical example is solved and sensitivity analyses are conducted. Finally, numerical examples extracted from the related literature are solved to illustrate the efficiency of the model and to demonstrate its advantages compared to other developed models. The remainder of this paper is organized as follows: in Section 2, the proposed problem is explained in detail and a non-linear two-stage stochastic model is presented. Furthermore, a linearization method is applied to linearize the model. In Section 3, the SAA method is presented. In Section 4, in order to clarify the proposed problem, an illustrative numerical example is solved and also, sensitivity analyses are carried out. In Section 5, by means of solving ten numerical examples extracted from the literature, the effectiveness of the solution method, as well as the advantages of the constructed model are examined. Finally, the conclusions and hints for future studies are given in Section 6. **2. Problem Statement**** ** In this section, a two-stage stochastic mathematical programming model is developed for designing a CMS. It is assumed that a set of parts, i=1,…,P , each having an uncertain demand with a known probability distribution, should be produced by a set of machine types, k=1,…,M . Each machine type has a limited capacity which is known a priori. This is assumed that there are Ri processing routes for the production of part i . Each of these routings can be independently implemented in the production of part i . In each routing, the sequence of operations and processing times are known a priori. Unsatisfied part demands which can be as a result of limited machine capacity or high manufacturing costs are outsourced. The suggested subcontracting approach is similar to that proposed by Mohammadi and Forghani [10]. The outsourcing price of each part is also a random variable with a known probability distribution. Machine duplication (i.e., purchasing machines) is allowed and machines of the same type are allocated to the same cell. Given the budget constraint and cell size limit (the maximum number of machines that can be assigned to a cell), machines are grouped into a maximum of Cmax cells. A schematic illustration of the proposed problem is given in Figure 1. **[Please insert Figure 1 about here]** Two-stage stochastic programming is one of the methods which can be used to take the uncertainty into consideration. It has two different types of decision: the first- and second-stage decisions. The first-stage decisions, also called strategic decisions, are made before the realization of the uncertain parameters. In the second stage, where the uncertain parameters are realized, operational and tactical decisions are made [27]. In the suggested problem, we use this method to tackle the uncertainty which arises from part demands and outsourcing costs. **2.1. Notations** The following notations are applied throughout this paper. *Indices*

i | index for parts ( i=1,…,P, where Pis the number of parts) |

j | index for processing routes ( j=1,…,Ri, where Riis the number of processing routes of part i) |

k | index for machine types ( j=1,…,M, where Mis the number of machine types) |

l | index for cells ( l=1,…,Cmax, where Cmaxis the maximum number of cells allowed) |

s | index for scenarios ( s=1,…,S, where Sis the number of scenarios) |

*Parameters*

πs | probability of scenario s |

dis | demand of part iin scenario s |

cisO | unit outsourcing cost of part iin scenario s |

cijP | unit production cost of part iusing processing route j |

cikk’A | unit intra-cell material handling cost for transporting part ibetween machine types kand k’ |

cikk’E | unit inter-cell material handling cost for transporting part ibetween machine types kand k'( cikk’E≥cikk’A) |

fijkk’ | number of times that part iin processing route jis transported between machine types kand k’ |

tijk | processing time of part ion machine type kin processing route j |

ckM | purchase price of one unit of machine type k |

Tk | available time on machine type k |

ckI | idleness cost of machine type kper unit time |

Nkmax | maximum number of machine type kallowed to be purchased |

NM | maximum number of machine types allowed in a cell |

B | budget available for purchasing machines |

wS | estimated total variable cost based on Sscenarios |

*First stage decision variables*

zkl | =1 if machine type kis assigned to cell l; otherwise 0 |

nk | number of machine type kto be purchased |

*Second stage decision variables*

pijs | amount of part ito be produced using processing route jin scenario s |

ois | amount of part ito be outsourced in scenario s |

uks | unused capacity (time) of machine type kin scenario s |

**2.2. Mathematical model** ** ** The two-stage mathematical model of the proposed problem is given below: ** **

wS=min∑s=1Sπs∑i=1P∑j=1RicijPpijs+∑i=1PcisOois+∑k=1MckIuks+∑k=1M-1∑k’=k+1M∑i=1P∑j=1Rifijkk’pijscikk’A∑l=1Cmaxzklzk’l+cikk’E1-∑l=1Cmaxzklzk’l. | (1) |

Subject to: | |

nkNkmax≤∑l=1Cmaxzkl≤1, ∀k, | (2) |

∑l=1Cmaxzkl≤nk, ∀k, | (3) |

∑k=1Mzkl≤NM, ∀l, | (4) |

ois+∑j=1Ripijs=dis, ∀i,s, | (5) |

∑i=1P∑j=1Ritijkpijs+uks=Tknk, ∀k,s, | (6) |

∑k=1MckMnk≤B, | (7) |

nk∈0,…,Nkmax, ∀k, | (8) |

zkl∈0,1, ∀k,l, | (9) |

pijs,ois,uks≥0, ∀i,j,k,s. | (10) |

In the proposed model, objective function (1) minimizes the estimation of the expected total variable cost which composed of production, subcontracting, idleness and material handling costs (i.e., intra- and inter-cell material handling costs) according to different realizations of uncertain parameters. Constraints (2) and (3) jointly represent that if a machine type is purchased, it is only assigned to a single cell. Constraint (4) restricts the number of purchased machines in each cell. Constraint (5) ensures that under various scenarios, the demand of parts is satisfied using production and outsourcing. Constraint (6) guarantees that machine capacity is not exceeded. To be more precise, the amount of used and unused time of machine type k should be equal to its total available time. Constraint (7) presents the budgetary limitations for purchasing machines. Constraint (8) states that the number of each machine type is an integer variable and also limits the number of each machine type. Finally, constraints (9) and (10) indicate the type of other decision variables. **2.3. Mathematical model linearization** The mathematical model presented in subsection 2.2 is a mixed-integer nonlinear program (MINLP) due to the existence of a nonlinear term in the fourth cost component of objective function (1). Therefore, a linearization method is applied to transform the model into a mixed-integer linear program (MIP) which is quite efficient to be solved optimally using high-performance solvers such as GURUBI, XPRESS or CPLEX. In doing so, objective function (1) is rewritten as Eq. (11).

min∑s=1Sπs∑i=1P∑j=1RicijPpijs+∑i=1PcisOois+∑k=1MckIuks+∑k=1M-1∑k’=k+1M∑i=1P∑j=1Rifijkk’cikk’Epijs-∑k=1M-1∑k’=k+1M∑l=1Cmaxzklzk’l∑i=1P∑j=1Ri∑s=1Sπsfijkk’cikk’E-cikk’Apijs. | (11) |

Then, a new auxiliary variable, λkk’l , is defined to be replaced with nonlinear term zklzk’l∑i=1P∑j=1Ri∑s=1Sπsfijkk’cikk’E-cikk’Apijs in Eq. (11). Eq. (11) aims at minimizing a negative coefficient of this nonlinear term, therefore, inequality λkk’l≤minBMzkl,BMzk’l,∑i=1P∑j=1Ri∑s=1Sπsfijkk’cikk’E-cikk’Apijs,∀k’>k is valid, where BM is an upper bound on ∑i=1P∑j=1Ri∑s=1Sπsfijkk’cikk’E-cikk’Apijs (it is reasonable to let BM=∑i=1P∑j=1Ri∑s=1Sπsfijkk’cikk’E-cikk’Adis because, pijs≤dis,∀j , see constraint (5)). From the other side, according to constraints (2) and (9), we know that ∑l=1Cmaxzklzk’l≤1, ∀k’>k . Therefore, inequality ∑l=1Cmaxλkk’l≤∑i=1P∑j=1Ri∑s=1Sπsfijkk’cikk’E-cikk’Apijs,∀k’>k is also valid. Furthermore, we know that cikk’E≥cikk’A , thus, it is correct to derive inequality λkk’l≥0,∀k’>k . Now, based on these three sets of inequalities, constraints (13)–(16) are added to the model, in order to complete the linearization. The linearized mathematical model of the problem is given below.

wS=min∑s=1Sπs∑i=1P∑j=1RicijPpijs+∑i=1PcisOois+∑k=1MckIuks+∑k=1M-1∑k’=k+1M∑i=1P∑j=1Rifijkk’cikk’Epijs-∑k=1M-1∑k’=k+1M∑l=1Cmaxλkk’l. | (12) |

Subject to: (2)–(10), | |

∑l=1Cmaxλkk’l≤∑i=1P∑j=1Ri∑s=1Sπsfijkk’cikk’E-cikk’Apijs, ∀k’>k,l, | (13) |

λkk’l≤∑i=1P∑j=1Ri∑s=1Sπsfijkk’cikk’E-cikk’Adizkl, ∀k’>k,l, | (14) |

λkk’l≤∑i=1P∑j=1Ri∑s=1Sπsfijkk’cikk’E-cikk’Adizk’l, ∀k’>k,l, | (15) |

λkk’l≥0, ∀k’>k,l. | (16) |

It should be mentioned that the model above which is also called as the SAA problem, includes 32MM-1Cmax+3M+Cmax+SP+M+1 constraints, SM+P+∑iRi+M2M-1Cmax positive variables, MCmax binary variables and M integer variables. * * **3. The SAA ****method** SAA is a Monte Carlo simulation-based solution method in which large numbers of generated scenarios are used to find bounds of the objective function of a stochastic problem. In the proposed model, outsourcing costs and part demands are stochastic parameters with known probability distributions. In this section, the stochastic model is formulated in a concise form such that the implementation of the SAA method can be easily explained. To do so, by ignoring the indices of the parameters and the decision variables, they are presented in a bold fashion (e.g., zkl is denoted by z ). The stochastic parameters are presented by ξω=coω,dω , where ω∈Ω ( Ω is the set of all scenarios with a known probability distribution P ) denotes a *scenario *which becomes known while making second-stage recourse decision p,o,u . Also, the objective value under decision z,p,o,u and a particular realization of uncertain parameter co is indicated by wz,p,o,u,coω . Now, assume that X1 denote the set of constraints that are not affected by the uncertainty; this kind of constraints are called as the first-stage constraints. Similarly, assume that X2z,n,ξω denote the set of constraints affected by the uncertainty and the first stage decision; this kind of constraints are also called as the second-stage constraints. According to these definitions, the concise form of the proposed stochastic problem is as follows:

w*=minz,n∈X1EPQz,n,ξω, | (17) |

where

Qz,n,ξω=minwz,p,o,u,cωp,o,u∈X2z,n,dω. | (18) |

In Eq. (17), the vector of binary and integer variables z and n correspond to the first-stage strategic decisions, and Qz,n,ξω represent the optimal total variable cost associated with the given strategic decision z,n and a particular realization of uncertain parameter ξω . Given probability distribution P , S sample scenariosis generated from Ω using Monte Carlo simulation (the samples are denoted by ω1,ω2,…,ωS ) in order to approximate the expected value function EQz,n,ξω . Then, a deterministic optimization problem (also called the SAA problem) specified by the generated samples is solved. The following is the SAA problem in a concise form.

wS=minz,n∈X11S∑s=1SQz,n,ξωs. | (19) |

In fact, the problem above (i.e., Eq. (19)) is a concise form of the model presented in Sub-section 2.3 (where X1 is the set of constraints (2)–(4) and (7)–(9), X2z,n,dω is the set of constraints (5), (6), (10) and (13)–(16) and πs=1/S ). It is worth mentioning that to provide a good statistical estimation of w* , the approximation process should be repeated several times using different independent sample scenarios. At each experiment, the best solution, as well as the estimation of the upper and lower bounds of the objective function, are calculated. This process is repeated until a satisfactory optimality gap is achieved. According to Shapiro and Homem-de-Mello [28], the steps of the SAA method can be summarized as follows: **Step 1.** Randomly generate T independent samples of stochastic parameters each with S scenarios, i.e., ωt1,ωt2,…,ωtS for t=1,…,T . Solve the two-stage problem (i.e., Eq. (19)) T times with S scenarios to obtain objective values wS1,wS2,…,wST and candidate solutions z1,n1,z2,n2,…,zT,nT . **Step 2.** Calculate the average of T optimal objective values by

w̅S=1T∑t=1TwSt, | (20) |

where w̅S is an unbiased estimator of EwS . It is well-known that the expected value of wS is less than or equal to the optimal objective value of the true problem, that is, EwS≤w* ( for example, see [29]). Thus, w̅S provides a statistical estimate for a lower bound on the optimal value of the true problem. The standard deviation of w̅S can be estimated by

σ̂T,S=1T-1T∑t=1TwSt-w̅S2. | (21) |

Calculate an approximate 1001-α% confidence lower bound on EwS or w* by

wL=w̅S-tα,T-1σ̂T,S, | (22) |

where tα,T-1 is the α -critical value of the t -distribution with T–1 degrees of freedom. **Step 3.** Randomly produce S’ ( S’ should be quite larger than S ) independent samples of stochastic parameters, i.e., ω1,ω2,…,ωS’ . Based on these S’ scenarios, compute the estimated objective value of candidate solution zt,nt,∀t=1,…,T , using the following estimator:

ŵS’zt,nt=1S’∑s=1S’ Qzt,nt,ξωs, ∀t=1,…,T. | (23) |

Note that for each candidate solution, this step involves solving S’ independent second-stage sub-problem given in Eq. (18). Choose ẑ*,n̂* as one of the candidate solutions z1,n1,z2,n2,…,zT,nT which has the smallest estimated objective value, that is, ẑ*,n̂*∈zt*,nt*t*∈argmint=1,…,TŵS’zt,nt . Calculate an approximate 1001-α% confidence upper bound on the true objective value of solution ẑ*,n̂* , that is, wẑ*,n̂* , by

wU=ŵS’ẑ*,n̂*+Φ-11-ασ̂S’, | (24) |

where Φ⋅ is the cdf of the standard normal distribution and σ̂S’ is computed by

σ̂S’=1S’-1S’∑s=1S’Qẑ*,n̂*,ξωs-ŵS’ẑ*,n̂*2. | (25) |

**Step 4.** Obtain a statistically valid bound with confidence at least 1001-2α% on the true optimality gap of solution ẑ*,n̂* by

ĝẑ*,n̂*=wU-wL. | (26) |

If the estimated optimality gap is small enough, the obtained solution is almost optimal for the first-stage decision of the true stochastic problem. Otherwise, the values of S or T should be increased and all the above steps should be repeated until a satisfactory optimality gap is achieved. **4. An illustrative example** In this section, a numerical example adopted from [2] is used to illustrate the proposed problem. Also, by using this example, sensitivity analyses are carried out to investigate the behavior of the solution in terms of the budget constraint and the SAA parameters. This numerical example includes 20 parts, 10 machine types, and 36 processing routes. Tables 1 and 2 contain the original data plus the additional ones that we have generated. The maximum number of cells ( Cmax ) is assumed to be 2, and the maximum number of machines allowed in a cell ( NM ) is assumed to be 5. The available budget for purchasing machines ( B ) is limited to $1500. In order to show the capability of the proposed problem in considering various distributions, it is assumed that the demand of part i is normally distributed with mean d̅i , and standard deviation d̅i/3 . Also, it is assumed that the outsourcing cost of part i is uniformly distributed over the interval 1.25maxjcijP,1.75maxjcijP . This example is solved using the SAA method at confidence level 95% ( α=0.025 ) assuming S=30 , S’=2000 . A summary of the results is provided in Table 3. It should be noted that the mathematical model of the SAA problem is coded in the GAMS 24.5, and the GUROBI 6.0 is selected as the default MIP and LP solver. Computations are performed on a PC having Microsoft Windows 10 operating system with Intel(R) Core (TM) i7-4790K 4.00 GHz CPU and 16GB of RAM. In order to use all available CPU cores (4 cores and 8 threads) we set ‘*threads* = 0’ in the GAMS option (this can lead to considerable save in the CPU time). **[Please insert Table 1 about here]** **[Please insert Table 2 about here]** It took almost 351 seconds for the SAA method to solve the illustrative example, with an estimated total variable cost ( ŵS’ẑ*,n̂* ) equal to $65385.243. Refer to Table 3 to see the machine cells and the number of each machine type in each cell. According to Eqs. (22) and (24), the lower bound on the optimum total variable cost ( wL ) and the upper bound on the true total variable cost of the SAA solution ( wU ) equal to $64886.883 and $65617.360, respectively. This also means that the optimal total variable cost of the stochastic model ( w* ) lies in the interval [65100.815, 65226.736] with 95% confidence. Based on these bounds, the estimated optimality gap ( ĝẑ*,n̂* ) is obtained as 730.478. Also, by dividing ĝẑ*,n̂* by wU the relative estimated optimality gap is 1.11%. **[Please insert Table 3 about here]** ** ** In order to verify the solution sensitivity to the available budget, the attempted example is investigated considering different budgets (starting from $0 to $2250 by increments $250) and the results are plotted in Figure 2. As it can be seen in this figure, the total variable cost is very sensitive to the budget within the interval 250,1250 . Such a plot can provide a useful outlook to make decision towards the budget constraint. **[Please insert Figure ****2**** about here]** ** ** The solution quality and the computational time are highly influenced by the number of scenarios in the SAA problem ( S ) and the number of samples ( T ). Therefore, by solving the attempted example, a sensitivity analysis is conducted based on different values of these two parameters. The CPU times and the relative optimality gaps (regarding S= 5, 10, 15, 20, 25 and 30, and T= 5, 10, 15, 20, 25 and 30) are depicted in Figures 3 and 4, respectively. From Figure 4 it is seen that at T=5 , there is a significant difference between the relative optimality gap obtained for S=5 and S=30 . However, as T increases, S loses its significance in the solution quality. On the other hand, it can be seen in Figure 3 that at T=5 , the computation times in terms of various values of S are not much different, while by increasing T to 30 samples, the differences in the computation times increases. As a conclusion, by considering a trade-off between the solution quality and the computation time, it is reasonable to solve the mentioned example with S= 5 or 10 and T=30 . **[Please insert Figure ****3**** about here]** ** ** **[Please insert Figure ****4 ****about here]** ** ** **5. Computational results** In this section, ten numerical examples extracted from the literature are considered to demonstrate the suitability of the proposed stochastic model and the efficiency of the final solution compared to those derived from the literature. The missing data in the source papers are generated and added to the original dataset. In these problems, we assume that the stochastic parameters (i.e., the demands and outsourcing costs) are independent random variables with uniform distribution. Table 4 includes the specifications of the problems, as well as the SAA. In this table, ‘A’ means that the data of the corresponding parameter is available in the source paper. It should be noted that all the computations are performed on the same PC mentioned in Section 4. **[Please insert Table ****4 ****about here]** **5.1. SAA method ****vs.**** the expected value problem** In order to highlight the merit of solving the proposed stochastic problem against the expected value problem (the problem in which the stochastic parameters are substituted with their expected values) a common measure called Value of Stochastic Solution ( VSS ) is taken into consideration. Let z̅,n̅ denote the optimal solution of the expected value problem. So, in our problem, VSS is defined by VSS=wz̅,n̅-w* , where wz̅,n̅ is the true objective value of solution z̅,n̅ and w* is the optimum objective value of the true stochastic problem. Now, we can apply the same S’ scenarios (used in the SAA method) to estimate VSS . Let vss denote the estimated value of VSS ; thus, vss=ŵS’z̅,n̅-ŵS’ẑ*,n̂* , where ẑ*,n̂* is the solution derived from the SAA method. Based on these explanations, we made a comparison between these two approaches by solving the numerical examples given in Table 4. A summary of the comparison results, as well as the number of constraints, positive and discreet variables (including binary and integer variables) in each model is provided in Table 5. In this table, ‘Validation model’ refers to the model with S’ scenarios which is used to estimate the true objective value of a given solution. The estimated optimality gap reported in column ‘ ĝ ’, has been obtained at confidence level 95%, that is, α=0.025 . The CPU time of the SAA method, is the overall time of solving T number of SAA and validation problems. For the sake of time, in some problems we use smaller number of scenarios (to see the SAA parameter values refer to Table 4). Also, the CPU time of the expected value approach is the overall time of solving the expected value and validation problems. **[Please insert Table ****5 ****about here]** According to the results given in Table 5, we can see that the SAA method has solved the problems with satisfactory estimated gap in relatively acceptable time (even when a smaller number of scenarios had been used in the SAA problem). The relative estimated optimality gap of the SAA solution for each problem is also plotted in Figure 5 (relative optimality gap =ĝ/wU , where wU is obtained by Eq. (24)). In this plot, the first and the second numbers in the horizontal axis denotes the problem number and the number of scenarios used in the SAA problem, respectively. According to this plot, we can see that the relative estimated optimality gap of all the problems is below 2.5%. From the other side, according to Table 5, vss in all the problems is a positive number, except for Problem 1. In order to ensure that the true value of the optimality gap (i.e., VSS ) is also a positive number, we examine null hypothesis H0:VSS>0 against H1:VSS≤0 , by obtaining the *p*-value of test statistic z0=vss/σ̂S’SAA2+σ̂S’EV2/S’-1 . According to the test, except for Problem 1, the* p*-value of the remaining problems was almost 0. This implies that solving the stochastic problem using the SAA method gives a better solution in terms of the expected total variable cost, in comparison to the expected value approach. **[Please insert Figure ****5 ****about here]** **5.2. Comparison against the literature solutions** One of the differences between the proposed approach and the ones in the literature is that in this study there isn’t a *routing selection* concept; in fact, parts can be produced using multiple routings. On the other hand, in similar problems extracted from the literature, even though there exist multiple routes in producing a particular part, but a single route is chosen through a routing selection procedure, and finally parts are produced by means of a single routing. Therefore, in order to conduct a comparison, the following two approaches are considered: A) *single routing approach based on the literature solution*: In this approach, the cell formation result ( z̅ ), as well as the routing selection result are extracted from the literature, and then according to this information and the mean value of the stochastic parameters ( c and d ), a mathematical model is solved to determine the number of each machine type in each cell ( n̅ ). Afterwards, based on S’ scenarios used in the SAA method, the validation model is solved to estimate the objective value of z̅,n̅ which is denoted by wS’LS . B) *multiple routings approach based on the literature solution*: In this approach, according to the cell formation result ( z̅ ) in the literature solution and the mean value of the stochastic parameters ( c and d ), a mathematical model (in which the production through multiple routings is allowed) is solved to determine the number of each machine type in each cell ( n̅ ). Afterwards, the estimated objective value of solution z̅,n̅ (for this approach let denoted by wS’LM ) is obtained by solving the validation model. Now a comparison is carried out between the approach proposed in this research, and approaches A and B. The solutions of the ten numerical examples, as well as ensuing improvements in consequence of using any of the intended approaches are reported in Table 6. As it can be inferred from the table, the results of solving the problems using the proposed approach show considerable improvements in comparison to approaches A and B, especially, to A (see columns ‘Imp^{1}’ and ‘Imp^{2}’ in Table 6). Also, we can see that allowing production through multiple routings (approach B) resulted in improvement in the estimated total variable cost in comparison to single routing approach (see ‘Imp^{3}’ in Table 6). **[Please insert Table ****6 ****about here]** ** ** **5. Conclusions** In this paper, a new problem was attempted to design a CMS considering stochastic part demands and stochastic outsourcing costs. The idea of simultaneous multiple processing routes along with subcontracting was addressed in the proposed problem. According to this idea, each part can be produced simultaneously in multiple processing routes and unsatisfied part demands (as a result of limited machine capacity or high manufacturing cost) are outsourced. The problem was formulated as a two-stage stochastic program and a solution procedure based on the SAA method was suggested. The objective function aims at minimizing the summation of the production, outsourcing, material handling and machine idleness costs. To clarify the problem, an illustrative example was investigated in the case that uniform and normal distributions are used for the outsourcing costs and the part demands, respectively. Based on this example, a sensitivity analysis was carried out to study the behavior of the resulted solution in terms of various available budget. Furthermore, some experiments were performed to evaluate the solution quality and the computation time in terms of the number of samples and scenarios in the SAA method. Then, we used ten numerical examples from the literature to demonstrate the performance of the solution method. The computation results indicated that the SAA method could produce efficient solutions with a satisfactory gap in a relatively reasonable computational time. Also, a common measure called VSS was used to highlight the merit of the proposed stochastic approach against the expected value approach. The results demonstrate that solving the stochastic problem using the SAA method is reasonably justified. Moreover, through these numerical examples, a comparison was executed between the obtained solutions and the ones reported in the literature. The comparison results showed that allowing production through multiple routings could lead to an improvement in the estimated total variable cost as compared to the single routing approach. Finally, to provide some directions for future research, the following issues are recommended:

- In this research, we used the SAA method to solve a stochastic CF problem. Even though this method is efficient in solving small and medium-sized instances, for larger problems the computational time is a concern. To overcome such a difficulty, addressing the Bender’s decomposing algorithm in the SAA method could be a possible remedy.
- Subcontracting, simultaneous multiple process routings, and machine duplication were the main issues that we addressed in the CMS design problem. However, these issues could be considered in an integrated problem in which the layout and scheduling problems are also incorporated.
- In order to deal with the uncertainty, we use a two-stage stochastic programming model. However, multi-stage stochastic programming methods could also be used to formulate the problem in an uncertain and dynamic environment.

**References** ** **

- Wemmerlöv, U. and Hyer, N. “Procedures for the part family/machine group identification problem in cellular manufacture”,
*Journal of Operations Management*,**6**, pp. 125–147 (1986). - Solimanpur, M., Vrat, P. and Shankar, R. “A multi-objective genetic algorithm approach to the design of cellular manufacturing systems”,
*International Journal of Production Research*,**42**, pp. 1419–1441 (2004). - Wu, T., Chung, S. and Chang, C. “Hybrid simulated annealing algorithm with mutation operator to the cell formation problem with alternative process routings”,
*Expert Systems with Applications*,**36**, pp. 3652–3661 (2009). - Kao, Y. and Lin, C. “A PSO-based approach to cell formation problems with alternative process routings”,
*International Journal of Production Research*, 50, pp. 4075–4089 (2012). - Arkat, J., Abdollahzadeh, H. and Ghahve, H. “A new branch and bound algorithm for cell formation problem”,
*Applied Mathematical Modelling*,**36**, pp. 5091–5100 (2012). - Boutsinas, B. “Machine-part cell formation using biclustering”,
*European Journal of Operational Research*,**230**, pp. 563–572 (2013). - Wu, X., Chu, Ch-H., Wang, Y and Yue, D. “Genetic algorithms for integrating cell formation with machine layout and scheduling”,
*Computers and Industrial Engineering*,**53**, pp. 277–289 (2007). - Krishnan, K., Mirzaei, S., Venkatasamy, V. and Pillai, V. “A comprehensive approach to facility layout design and cell formation”,
*International Journal of Advanced Manufacturing Technology*,**59**, pp. 737–753 (2012). - Chang, C., Wu, T. and Wu, C. “An efficient approach to determine cell formation, cell layout and intracellular machine sequence in cellular manufacturing systems”,
*Computers and Industrial Engineering*,**66**, pp. 438–450 (2013). - Mohammadi, M. and Forghani, K. “A novel approach for considering layout problem in cellular manufacturing systems with alternative processing routings and subcontracting approach”,
*Applied Mathematical Modelling*,**38**, pp. 3624–3640 (2014). - Solimanpur, M., Vrat, P. and Shankar, R. “A heuristic to minimize makespan of cell scheduling problem”,
*International Journal of Production Economics*,**88**, pp. 231–241 (2004). - Elmi, A., Solimanpur, M., Topaloglu, S. and Elmi, A. “A simulated annealing algorithm for the job shop cell scheduling problem with intercellular moves and reentrant parts”,
*Computers and Industrial Engineering*,**61**, pp. 171–178 (2011). - Arkat, J., Farahani, M. and Ahmadizar, F. “Multi-objective genetic algorithm for cell formation problem considering cellular layout and operations scheduling”,
*International Journal of Computer Integrated Manufacturing*,**25**, pp. 625–635 (2012). - Li, D., Meng, X., Li, M. and Tian, Y., “An ACO-based intercell scheduling approach for job shop cells with multiple single processing machines and one batch processing machine”,
*Journal of Intelligent Manufacturing*,**27**, pp. 283–296 (2016). - Wang, T., Wu, K. and Liu, Y. “A simulated annealing algorithm for facility layout problems under variable demand in Cellular Manufacturing Systems”,
*Computers in Industry*,**46**, pp. 181–188 (2001). - Jeon, G. and Leep, H. “Forming part families by using genetic algorithm and designing machine cells under demand changes”,
*Computers and Operations Research*,**33**, pp. 263–283 (2006). - Tavakkoli-Moghaddam, R., Javadian, N., Javadi, B. and Safaei, N. “Design of a facility layout problem in cellular manufacturing systems with stochastic demands”,
*Applied Mathematics and Computation*,**184**, pp. 721–728 (2007). - Schaller, J. “Designing and redesigning cellular manufacturing systems to handle demand changes”,
*Computers and Industrial Engineering*,**53**, pp. 478–490 (2007). - Safaei, N., Saidi-Mehrabad, M. and Babakhani M. “Designing cellular manufacturing systems under dynamic and uncertain conditions”,
*Journal of Intelligent Manufacturing*.**18**, pp. 383–399 (2007). - Arıkan, F. and Güngör, Z. “Modeling of a manufacturing cell design problem with fuzzy multi-objective parametric programming”,
*Mathematical and Computer Modelling*,**50**, pp. 407–420 (2009). - Ghezavati, V. and Saidi-Mehrabad, M. “Designing integrated cellular manufacturing systems with scheduling considering stochastic processing time”,
*International Journal of Advanced Manufacturing Technology*,**48**, pp. 701–717 (2010). - Das, K. and Abdul-Kader, W. “Consideration of dynamic changes in machine reliability and part demand: a cellular manufacturing systems design model”,
*International Journal of Production Research*,**49**, pp. 2123–2142 (2011). - Ghezavati, V. and Saidi-Mehrabad, M. “An efficient hybrid self-learning method for stochastic cellular manufacturing problem: A queuing-based analysis”,
*Expert Systems with Applications*,**38**, pp. 1326–1335 (2011). - Rabbani, M., Jolai, F., Manavizadeh, N., Radmehr, F. and Javadi, B. “Solving a bi-objective cell formation problem with stochastic production quantities by a two-phase fuzzy linear programming approach”,
*International Journal of Advanced Manufacturing Technology*,**58**, pp. 709–722 (2012). - Forghani, K., Mohammadi, M. and Ghezavati, V. “Designing robust layout in cellular manufacturing systems with uncertain demands”,
*International Journal of Industrial Engineering Computations*,**4**, pp. 215–226 (2013). - Heragu, S. and Chen, J. “Optimal solution of cellular manufacturing system design: Benders’ decomposition approach”,
*European Journal of Operational Research*,**107**, pp. 175–192 (1998). - Ruszczynski, A. and Shapiro, A.
*Handbooks in Operation Research and Management Science: Stochastic Programming*. Elsevier, Amesterdam, (2003). - Shapiro, A. and Homem-de-Mello, T. “On rate of convergence of Monte Carlo approximations of stochastic programs”,
*SIAM Journal**of**Optimization*,**11**, pp. 70–86 (2000). - Mak, W.K., Morton, D.P. and Wood, P.K. “Monte Carlo bounding techniques for determining solution quality in stochastic programs”,
*Operations Research Letters*,**24**, pp. 47–56 (1999). - Kazerooni, M., Luong, H. and Abhary, K. “A genetic algorithm based cell design considering alternative routing”,
*International Journal of Computer Integrated Manufacturing**Systems*,**10**, pp. 93–107 (1997). - Ramabhatta, F.V. and Nagi, R. “An integrated formulation of manufacturing cell formation with capacity planning and multiple routings”,
*Annals of Operations Research*,**77**, pp. 79–95 (1998). - Yin, Y. and Yasuda, K. “Manufacturing cells’ design in consideration of various production factors”,
*International Journal of Production Research*,**40**, 885–906 (2002). - Chan, F.T.S., Lau, K.W., Chan, P.L.Y. and Choy, K.L. “Two-stage approach for machine-part grouping and cell layout problems”,
*Robotics and Computer-Integrated Manufacturing*,**22**, pp. 217–238 (2006).