Multilevel hybrid method for optimal buffer sizing and inspection stations positioning

Designing competitive manufacturing systems with high levels of productivity and quality at a reasonable cost is a complex task. Decision makers must face numerous decision variables which involve multiple and iterative analysis of the estimated cost, quality and productivity of each design alternative. This paper adresses this issue by providing a fast algorithm for solving the buffer sizing and inspection positioning problem of large production lines by combining heuristic and exact algorithms. We develop a multilevel hybrid search method combining a genetic algorithm and tabu search to identify promising locations for the inspection stations and an exact method that optimizes rapidly (in polynomial time) the buffers’ sizes for each location. Our method gives valuable insights into the problem, and its solution time is a small fraction of that required by the exact method on production lines with 10–30 machines.

these two problems are interdependent (see Gershwin 2005, 2008;Colledani and Tolio 2005. We illustrate this by considering the Kanban strategies introduced by Toyota in the 1960s that have since become the paradigm of "lean manufacturing". These strategies essentially advocate areas of limited storage between successive machines in a production line. The storage areas allow a degree of decoupling between machines to increase the productivity of the line. The idea is to limit the impact of shortages on downstream machines and to avoid blocking machines that are upstream from one that is broken down. These areas should be limited because they are associated with frozen capital: they introduce storage costs and extended transit times for parts in the workshop. In an ideal "just in time" scenario, there would be no intermediate storage and finished parts would be pulled from the system as they are produced. In reality, machine breakdowns and shortages of raw material or operators make it impossible to ensure continuous workshop output. The sizes of the storage areas should depend on the likelihood of such events and the estimated costs of the associated service interruptions. If we decide the amount of storage without considering the quality dimension, we risk creating storage areas that contain important quantities of defective parts. These defective parts correspond to misused production-line time, and they undermine the effectiveness of the intermediate storage in terms of increasing productivity. Moreover, the costs associated with storage impact the budget available for improvements in quality and vice versa. In a similar manner, Inman et al. (2013) show that production design impacts quality and vice versa. Also, the location of an inspection station affects both the expected production cost per item and the production rate of the line. The quality control is obviously critical to any production system and has significant managerial impact. As a result, it is crucial to develop strategies to allocate inspection or quality control stations into a manufacturing processes that can help preventing all wastes resulting from unidentified defective items being processed (see Shetwan et al. 2011;Raviv 2013).
In a real manufacturing context, simulation can be used to determine the optimal buffers' sizes and the optimal positions of the inspection stations by considering all possible scenarios. However, this is not practical when the systems are complex: the number of scenarios is ( n m ) × 100 n for a system with n machines and m inspection stations where each buffer of equal size is discretized to 100 levels. Note that the system is required to have an inspection station at the end of the line. Since the position of this station is known, it is not included in the optimization. An alternative to simulation is to develop a realistic model and a fast optimization technique to solve it. The real model is however very complex because many factors impact quality and production. So, we need elaborate some simplifying assumptions in such a way that the resulting model is still realistic and can be solved optimally.
Researchers have investigated the integration of quality aspects and production policies for a single unreliable machine producing a single product type. The system starts in an "in-control" state producing conforming items and then switches to an "out-of-control" state and starts producing nonconforming items. Integrated models can be classified as: (1) integrated production and quality management (e.g., Hajji et al. 2010;Kutzner and Kiesmüller 2013;Mhada et al. 2014b;Naebulharam and Zhang 2014;Matta and Simone 2016; 2) integrated production and maintenance management (e.g., Ben-Daya 2002;Dhouib et al. 2012;Rivera-Gómez et al. 2013;Jafari and Makis 2015;Nourelfath et al. 2016);and (3) integrated production, quality, and maintenance management (e.g., Radhoui et al. 2009Radhoui et al. , 2010Njike et al. 2012;Zequeira et al. 2008;Rivera-Gmez et al. 2016;Bouslah et al. 2016). Some researchers (see Mandroli et al. 2006) focus on determining the optimal inspection-station position in n serial production lines with or without (1) scrapping, (2) reworking, and (3) offline repairs, without considering the buffer sizing.
A large part of published articles related to the inspection suppose the inspection reliable and without errors. They assumed that the produced items are subject to 100% inspection and there are no inspection errors. Some researchers have incorporated quality inspection errors caused, among others, by human failure in their models (see Hsu and Hsu 2013;Duffuaa and Khan 2002;Khan et al. 2014).
In the context of serial production lines with production and inspection machines that follow Bernoulli reliability and quality assumptions, Zhang (2010, 2011) provide important insight into the nature of production and quality bottlenecks. Such systems are encountered in automotive assembly and painting operations where the downtime is relatively short and the defects are a result of uncorrelated random events (Ju et al. 2013). Mhada et al. (2014a) consider a situation where every machine has to satisfy a demand specified as good parts per time unit. Moreover, the inventory size takes into account that the stock is a mixture of good and defective parts, and defective parts are generally eliminated at inspection stations. Mhada et al. (2014a) improve the continuous models of Gershwin (2005, 2008). The former model is more fluid in the sense that quality is treated as a continuous flow, whereas in the Kim and Gershwin models quality is considered to be discrete. They considered a line with 4 machines. Mhada et al. (2014a) develop decomposition methods to reduce the analysis of the line to a series on an equivalent machine that can be isolated and sequentially analyzed. They consider just one inspection station. Ouzineb et al. (2013Ouzineb et al. ( , 2014 generalize the model proposed in Mhada et al. (2014a) by optimizing the number of internal inspection stations (m ≥ 1) and their positioning, and develop an exact search method to solve it. The exact search method (ESM) is exhaustive: the search is guaranteed to generate all possible locations of inspection stations. For each location, the problem is reformulated as a network flow optimization problem that can be efficiently solved by a fast polynomial algorithm. This method is efficient for small instances, but it may need days for reasonably large instances. The solution time increases exponentially with the size of the problem, i.e. with n and m. For larger problems, we expect that the solution time may span over several months. We can have an idea on that from Table 4.
The present paper efficiently adapts a space partitioning technique that decomposes the search space into multiple levels. At each level, we use two well-known metaheuristics to explore some potential regions. A genetic algorithm (GA) is used to identify a population of promising locations of inspection stations. Tabu search (TS) then searches intensively around by exploring locations that are neighbors to some newborns in the population. The population evolves dynamically in a sense that worst locations are discarded and good ones are added during the search. We use a quick shortest path procedure (from ESM) to determine the optimal buffers' sizes and evaluate the objective function (the fitness) for each considered location. From time to time, GA replenishes the population. GA facilitates thus the exploration by guiding the search to unvisited regions with good potential, and this gives a certain diversification in terms of the regions to explore. The new approach is hence a multilevel hybrid heuristic (MHH) that provides certain balance between diversification and intensification. This hybrid approach is effective when the number of solutions is huge: it finds optimal or near-optimal solutions in a fraction of the time required by ESM. We also discuss some important properties of the production line and its sensitivity to system parameters.
Similar ideas have been used to solve some high-dimensional combinatorial optimization problems in other contexts. In Ouzineb et al. (2011), the authors apply a GA to solve an optimal design model for assembly/disassembly manufacturing networks. The objective is to maximize the production rate subject to a total cost constraint, the machines are chosen from a list, and the buffers' sizes are within a predetermined range. Nahas et al. (2014) applies a GA-TS algorithm to optimize the nonhomogeneous redundancy of multistate series-parallel systems. In Ouzineb et al. (2010), the authors use space partitioning to solve two design optimization problems: the first is the expansion scheduling of multistate series-parallel systems, and the second is the redundancy allocation of binary-state series-parallel systems. We adapt space partitioning and GA-TS techniques to efficiently solve the design problem addressed in this paper.
The remainder of this paper is organized as follows. "Problem formulation and discussion" section presents the problem statement and discusses some theoretical properties, and "Methods" section describes our approach. "Numerical results" section presents numerical results for ESM (Ouzineb et al. 2013(Ouzineb et al. , 2014 and our method. "Conclusions and future work" section provides concluding remarks.

Problem formulation and discussion
Figure 1 illustrates the production system studied. It consists of n machines separated by n buffers with m + 1 inspection stations. The machines can be either up or down, starved or blocked. Machine M i is starved if one of the upstream machines is down and all the buffers between that machine and M i are empty. M i is blocked if one of the downstream machines is down and all the buffers between that machine and M i are full. When an operating machine is neither starved nor blocked, it continuously transfers parts from the upstream buffer to the downstream machines. We assume that the first machine can never be starved, and an inspection station is located after the last machine to ensure the quality of the parts received by the customer. The machines M i , i = 1 . . . n can be modeled as a continuous-time Markov chain that produces a single part type with two quality levels, conforming and nonconforming, with a predefined ratio β of nonconforming parts to conforming parts. We assume that the machines have the same maximum production rate k, failure rate p, and repair rate r.
Let d be the demand rate for conforming parts, x i the inventory level for buffer i, d i the long-term average number of parts pulled per unit of time from stock x i , c p the storage cost per time unit and per part, c I the inspection cost per pulled part, and a n des the required availability rate of conforming finished parts.
A binary variable i determines whether or not there is a station before machine M (i+1) to identify nonconforming parts. We assume that there are m inspection stations in the line (m is unknown), i.e., n−1 i=1 i = m, m ∈ {0, 1, . . . , (n − 1)} and n = 1. The problem is to minimize the long-term average global cost per unit time of storage, shortages, and inspection: The constraints are n−1 i=1 i = m, i binary ( n = 1), and m ∈ {0, 1, . . . , (n − 1)}. The objective function (1) can be rewritten as (see Mhada et al. 2014a): with a = (a 1 , a 2 , . . . , a n−1 ) and = ( 1 , 2 , . . . , n−1 ). Here, • a i , i = 1 . . . n − 1 is the total work-in-progress availability coefficient at buffer x i ; this coefficient has a lower bound of max r (r+p) • and finally , and The problem can now be formulated as follows: We aim to find the optimal solution, i.e. the minimal average global cost system structure, (a, , m) that minimizes J and satisfies the constraint (8) on the number of inspection stations, constraints (9) on buffers sizes bounds, and integrality constraints (10)-(11) on admissible values of and m. The proposition below discusses the sensitivity of the production line to c I . A numerical example is given in "Sensitivity analysis" section to illustrate this result.
T F (a n des , a n−1 , ) = ρ n c p Proof The minimum total storage and inspection cost is The optimal total storage and inspection cost as a function of c I is where (a * , * ) = argmin(J (a, )). Observe that does not depend on the value c I , so it is constant on c I and referred to as Constant(c I ). Thus, we can write Eq. (12) as follows: Hence, the total cost as a function of c I is a piecewise linear function; its slope depends only on m and the positions of the inspection stations.
The cost function J is likely to be convex as conjectured by Conjecture 1 and confirmed by the numerical results in "Characteristics of the production line" section. The convexity of the cost function is used to reduce the number of iterations. Actually, the local minimum of a convex function is also a global minimum, and there are many specialized methods and mathematical tools for optimizing convex functions. Conjecture 1 • The minimal total cost J is a convex function of .
• The minimal total cost J is a convex function of the number of inspection stations m.
• If we divide the production line into two parts with the same number of machines and we fix the number of stations in the first part, the minimal total cost J is a convex function of the number of internal stations in the second part.

Methods
In Ouzineb et al. (2013Ouzineb et al. ( , 2014, authors present an exact search methd (ESM) for the problem studied in this paper. For a given location of inspection stations, they show that the optimizing of the buffers' sizes with ESM can be quickly done in polynomial time. Actually, the bottleneck of ESM is the number of possible locations. This number could be reduced by using the property of convexity (Conjecture 1) but nevertheless, it would stay very large. The innovation in this paper is to smartly select potential locations using an evolutionary algorithm combined with tabu search. We summarize ESM in "Exact search method (ESM)" section and then present our new approach in "Multilevel hybrid heuristic (MHH)" section.

Exact search method (ESM)
For a fixed , Ouzineb et al. (2013Ouzineb et al. ( , 2014 shows that the term T (i) (a i , a i−1 , ), and thus the objective function itself, is separable by the variables a i . ESM reformulates then the problem defined by (7)-(11), for a given , as a shortest path problem defined on a network (described below). It uses a standard shortest path algorithm to find the optimal buffers' sizes. For this, we discretize the continuous variables a i to take discrete values 1, 2, ..., 100%. For a fixed , the value of m is determined straightforwardly. So, the constraints (8), (10), and (11) are satisfied. They do not depend on a and consequently could be omitted.
Consider the connected network G (N, A) consisting of a set of nodes N and a set of arcs A, as depicted in Fig. 2. With each machine we associate a set of nodes representing all the possible buffers' sizes, which range from 1 to 100%. The set N consists of 100 × (n − 1) nodes (representing the n − 1 machines) plus two additional nodes N 0 and N n that are adjacent to the nodes representing the first and last machines. The set A consists of 100 2 × (n − 2) + 200 arcs weighted as follows: • N 0 is connected to the nodes N j 1 , j = 1, . . . , 100 representing machine 1. The weight of each arc (N 0 , N j 1 ) is null. N n is connected to the nodes N j n−1 , j = 1, . . . , 100 representing machine n − 1. The weight of each arc (N . . , n − 2 and 1 ≤ j, k ≤ 100. Nodes, and all incident arcs, corresponding to values that are not in the intervals imposed by constraints (9) are removed from the network. Algorithm 1 presents ESM. It finds optimal solutions but takes hours for instances with 20 machines. It is unable to solve some larger instances with 30 machines in one month! Algorithm 1 Exact search method (ESM) Evaluate the optimal cost function J, called also Fitness, and a for the current λ using Algorithm 2.

Algorithm 2 Fitness evaluation
Given λ. Construct the network G. Find a shortest path in G (i.e. the buffers sizes vector a) and set J to its cost. return (J, a).

Multilevel hybrid heuristic (MHH)
Our approach combines ideas from ESM and metaheuristic techniques. The choice of a search space and a neighborhood structure is a critical step in metaheuristic design (see Gendreau 2002). Our search space is the space of all possible line structures, i.e. the vectors (a, , m) that are feasible solutions to the system (8)-(11). As discussed above, the most determinant element of these vectors is because if we know , we can compute m in a straightforward manner and a in a polynomial time. Consequently, our search space S in this section is composed of all possible values. We reduce it by: (1) partitioning this space into a set of disjoint subspaces as explained in "Partitioning the search space" section; (2) applying a hybrid heuristic to selected subspaces in order to find potential solutions; we use Algorithm 2 above to evaluate the fitness of each selected (in a subspace) and to find the optimal buffers' sizes a. See "Hybrid heuristic" section for more details. The multilevel hybrid heuristic pseudocode is then given in "MHH pseudocode" section. This approach finds, as reported in the numerical results, near-optimal solutions in a fraction of the time required by ESM.

Partitioning the search space
We refer hereunder to as a solution vector or simply solution. Each solution vector is assigned an address, also called a level, as explained below.

Definition 1
We divide the solution vector into two equal parts and create its address r as follows: r = Address( ) = ⌊n/2⌋ l=1 l . A search subspace of address r, denoted S r , is defined to be the set of solutions that have the same address, i.e. equal to r.
The lower bound of r is 0 and its upper bound is m. It is clear that (S r ) 1≤r≤m is a partition of the space S, i.e., S r 1 S r 2 = ∅ ∀ r 1 � = r 2 and S = m r=0 S r .
Example Consider a system with n = 20 machines and m = 5 inspection stations: 2 stations in the first part and 3 in the second part: From the above definition, address( ) = 2. So, ∈ S 2 .

Hybrid heuristic
The goal here is to locate promising regions in a given subspace S r . Our hybrid optimization is based on a combination of a GA, which provides diversification, and TS, which provides intensification. The hybrid heuristic pseudocode is presented in Algorithm 3 and a flowchart is given in Fig. 3. N rep and N c are GA tuned by experimentation. The main steps are explained in the paragraphs that follow.
Step 1. Generate a random population of Ns (chromosomes) solutions with the same address r.
Step 2. Select randomly from the population two parent solutions λ1 and λ2 that will "reproduce" using a 1-point crossover operator.
Step 3. Produce a new solution λ (child or offspring) from the selected pair of parent solutions P arent 1 and P arent 2.
Step 4. Find the best solution λ in V (λ), the effective neighborhood of the new child λ, using TS. Evaluate its value J and its corresponding optimal a . if J < J * then J * = J , a * = a , λ * = λ . end if Step 5. If λ is better than the worst solution currently in the population, it replaces that solution; otherwise, it is discarded. if k < Nrep then k = k + 1 and go to Step 2. end if Step 6. Decide if the population needs to be replenished and then to start a new generation: if l < Nc then k = 0, l = l + 1, go to Step 1. end if return (J * , a * , λ * ).

In
Step 3, we use a 1-point crossover operator, which, as illustrated in Table 1, creates a child string for a given pair of parent strings 1 and 2 by: 1. Copying the string elements belonging to the first part from 1 ; 2. Copying the rest of the string elements from 2 .
Observe that the child has the same address as its parents. In Step 4, we use TS, in which we iteratively move from the current solution to a new solution in its neighborhood Neighborhood( ) until some stopping criterion has been satisfied. For a given , we define Neighborhood( ) = {solutions (inspection stations configurations) obtained by applying a single move to }. Moving to a new solution involves changing the positions of two inspection stations; the number of stations does not change. Note that we stay in the same level when we move to a neighbor solution. We use Algorithm 2 to evaluate the fitness of each new solution, and we use theses fitness values to compare different solutions.
TS uses a memory structure: a potential solution is declared tabu so that the algorithm does not visit it again; this prevents cycling. At each iteration, we select the best solution ′ in a subspace V( ) ⊂ Neighborhood( ) as a tabu solution. The subspace V( ), called the effective neighborhood, is obtained by eliminating the tabu solutions from Neighborhood( ). The tabu solutions are stored in a short-term memory, called a tabu list, which contains the solutions that have been visited in the recent past. We use a variable-size tabu list (tabu tenure) because this is generally more efficient (Gendreau 2002;Ouzineb et al. 2008). Finally, our stopping criterion is defined by the maximum number of local iterations (mnli), i.e. the number of visited neighbors that do not improve the bestknown solution. The value of mnli is determined by the experimentation.

MHH pseudocode
In MHH, at each iteration we move from the current level to a new level. We divide the search space into a set of disjoint subspaces (levels) using the partitioning technique of "Partitioning the search space" section and then apply the heuristic presented in "Hybrid heuristic" section to each selected subspace. Algorithm 4 presents the proposed pseudo code and Fig. 4 the corresponding flowchart.  Step 1. Choose the initial solution such that there are no inspection stations (λi = 0 ∀i ≤ n − 1).
Step 2. Randomly increment the address r and calculate the best values of J, a, and λ in Sr using Algorithm 3. if J < J * then J * = J, a * = a, λ * = λ. end if Step 3. Check the termination criterion: the maximum time (maxt) without improving the best known solution. If it is satisfied, return the best solution found. Otherwise, go to step 2. return (J, a, λ).

Numerical results
We use four optimization problems (benchmarks) to test our algorithms: I10 has 10 machines, I20 has 20 machines, I30 has 30 machines and I40 has 40 machines. To get an insight into the algorithms behavior, we fixed the value of m in the model (7)-(11). The optimal objective value can be obtained equivalently by varying m between 1, . . . , n − 1 , sloving the resulting models and take the minimum possible. We implemented MHH and ESM in C ++ under the same conditions (computer architecture, operating system, etc). The tests were performed on an Intel Core i7 at 2.8 GHz with 8 GB of RAM, running Linux. Table 2 shows the parameters we used.

Comparing MHH to ESM
ESM can find the optimal solutions for I10 and I20. We compare the MHH solutions to the optimal solutions to evaluate the performance of MHH. ESM is unable to solve I30 and I40, which has been newly proposed in this paper, for large m. The reduction percentage in terms of CPU time (RPCPU) is defined as follows: For I10, MHH and ESM find the same optimal solutions in comparable times. Tables 3  and 4 give the results for I20 and I30, with m = 1, . . . , 10. The first column indicates the number m of inspection stations. The MHH cost is given in the second column; if the solution is optimal the entry is in italics. MHH finds optimal solutions for all the instances solved by ESM. The column MHH/Best gives the time to find the best solution, and the column MHH gives the total computational time. ESM takes hours on I20 for m ≥ 5. MHH outperforms ESM in terms of the computational time; for I20, it finds optimal solutions in a fraction of ESM time. We also report results for a "pure" genetic algorithm (GA) which is obtained by modifying Algorithm 3 as follows: • Remove Step 4. • Replace Step 5 by: If the new solution is better than the worst solution currently in the population, it replaces that solution; otherwise, it is discarded.
GA is coded in the same technology and runs on the same machine as MHH and ESM. MHH results are signifcantly better in terms of quality and CPU time than GA. For larger instances, we omit GA and present results of MHH only. For I30, MHH finds "optimal" solutions in a reasonable time. ESM cannot solve the instances with m ≥ 7. For example, for m = 7, ESM did not converge after a month of running time. We believe that the MHH solution is optimal or near-optimal for m = 7 and beyond, but we cannot confirm this. Large problems require fast heuristics such as MHH that provide good results in an acceptable time. Table 5 presents results for I40 which consists of very large instances. In this table, we report the cost values for both MHH and ESM because they are slightly different for  The MHH solution is optimal or near-optimal for m = 3 and below. We think that the optimal solution is obtained when m = 5 despite the fact that we cannot solve it by the exact method. Table 6 shows the results for I20 when the computational time is limited to one hour. Both methods find optimal solutions for m ≤ 4; MHH finds better solutions for m ≥ 5. Figure 5 shows the optimal cost as a function of (the number of internal inspection stations) m = 1, 2, . . . , 19. The cost is a convex function of m. The optimal number of internal inspection stations is m = 3, and the optimal positions are 2 = 1, 7 = 1, and 18 = 1. The results obtained by MHH and ESM are identical. Figure 6 shows the best cost found by MHH as a function of m = 1, 2, . . . , 19. The cost is a convex function of m. The optimal value for m is m = 5, and the optimal positions are 1 = 1, 3 = 1, 7 = 1, 15 = 1, and 29 = 1. These results support Conjecture  1. Convexity can play an important role in reducing the solution time. Using convexity, we can conjecture that even I30 is solved optimally by ESM and MHH because the cost starts increasing for m ≥ 6. Figure 7 shows the best cost found by MHH as a function of m = 1, 2, . . . , 10 for I40. The cost is a convex function of m. The optimal value for m is m = 5.

Sensitivity analysis
We now study the sensitivity of the optimal cost and number of inspection stations to the parameters of the model. We study particularily β and c I because they directly affect quality. Table 7 gives the parameters used for this study; we investigate I20.

Inspection cost c I
We begin by varying c I , the inspection cost per time unit and per part (for a fixed β = 0.1 ). The upper plot in Fig. 8 shows that for values of c I close to c p (the storage cost per time unit and per part), the optimal number of stations is high, while for high values   of c I , the optimal number decreases until it reaches zero. This is because the system encourages storage when it is profitable to do so, and it reduces nonconforming parts by imposing more inspection stations when that option is better.
When c I increases, by Proposition 1, m i=1 d(1 + β) n−e i must decrease (m or (n − e i ) decreases). The stations will thus be closer to the end of the line. When a station is placed after machine n − 1, any increase in c I will make it unnecessary, so the optimal number of stations will be reduced by 1 (see Table 8). Figure 9 shows the cases m = 1 (upper plot) and m = 0 (lower plot). For m = 1 and for each station position, the total cost as a function of c I is linear. For m = 0, the total cost as a function of c I is linear. Figure 10 and Table 9 show the behavior of the line as a function of β. As β increases, the quantity of nonconforming parts increases. If we do not eliminate the nonconformity from the system, the machines will work harder and the stocks will fill faster, and the line will not be able to meet the final demand. Here the cost varies exponentially (Fig. 11).

Conclusions and future work
We proposed an efficient hybrid approach, based on ideas from an exact method and metaheuristics, for the buffer sizing problem in unreliable homogeneous production lines with several inspection stations. This is a difficult mixed integer nonlinear program. Our approach combines a genetic algorithm and tabu search to identify profitable configurations (locations of inspection stations). For these locations, we use an exact approach to decide the optimal buffers' sizes. Our final goal is to find an optimal   or near-optimal design as rapidly as possible. This hybrid approach provides a balance between diversification and intensification and works well on homogeneous production lines with up to 40 machines. MHH is significantly faster than ESM and produces solutions that are equally good.   Fig. 9 Optimal total cost as a function of c I (m * = 1 and m * = 0) Future research should focus on developing an optimization method for more realistic nonhomogeneous production lines and solving it by combining simulation and optimization. The model developed here could be used as an approximation model for more complex real life design problems. Fast optimization techniques, like the one proposed   in this paper, would help in rapidly selecting potential scenarios/configurations for a more realistic simulation, and hence reducing the solution time without losing solution applicability in real life context.