Mixed-model Sequencing with Stochastic Failures: A Case Study for Automobile Industry

In an MMAL, a set of workstations are connected by a conveyor belt. Products, launched with a fixed rate, move along the belt at a constant speed of one time unit (TU). The duration between two consecutive launches is called the cycle time $c$, and we define the station length $l_{k}\geq c$ in TU as the total TU that the workpiece requires to cross the station $k\in K$. Operators work on the assigned tasks and must finish their job within the station length, otherwise, the line is stopped or a so-called utility worker takes over the remaining job. The excess work is called work overload. The sequence of products therefore has a great impact on the efficiency of the assembly line. MMS determines the sequence of a given set of products $V$ by assigning each product $v\in V$ to one of the positions $t\in T$.

Formulating the MMS problem based on the vehicle configurations instead of vehicles is usual (Bolat and Yano, 1992; Bard et al., 1992; Scholl et al., 1998), however, automobile manufacturers offer billions of combinations of options (Pil and Holweg, 2004). When this high level of customization is combined with a short lead time promised to the customers, each vehicle produced in a planning horizon becomes unique. In this study, the vehicles are sequenced instead of configurations since the case study focuses on an automobile manufacturing facility with a high level of customization. In order to do so, we define a binary decision variable $x_{vt}$, which takes value of 1 if vehicle $v\in V$ is assigned to position $t\in T$. The processing time of vehicle $v\in V$ at station $k\in K$ is notated by $p_{kv}$. The starting position and work overload of the vehicle at position $t\in T$ for station $k\in K$ are represented by $z_{kt}$ and $w_{kt}$, respectively. Table 1 lists all the parameters and decision variables used in the proposed model. While second-stage decision variables are scenario-dependent, we drop such dependency for notation simplicity throughout the paper unless it is needed explicitly for clarity.

In this paper, we adopt the side-by-side policy as a work overload handling procedure. A utility worker is assumed to work with the regular worker side-by-side, enabling work to be completed within the station borders. The objective of MMS with the side-by-side policy is to minimize the total duration of work overloads, i.e., the total duration of the remaining tasks that cannot be completed within the station borders. The regular operator stops working on the piece at the station border so they can start working on the next workpiece at position $l_{k}-c$ in the same station.

We illustrate an MMAL with the side-by-side policy in Figure 2 which represents a station that processes five vehicles. The left and right vertical bold lines represent the left and right borders of the station. Assume that the cycle time $c$ is 7 and the station length $l_{k}$ is 10 TU, i.e., it takes 10 TU for the conveyor to flow through the station. This specific station processes two different configurations of vehicles: configurations A and B. While configurations A requires option 1, configurations B does not, so the processing times of configurations A and B are 9 and 5, respectively. Figure 2 illustrates the first five vehicles in the sequence which is [A, B, B, A, A]. The diagonal straight lines represent the position of the vehicle in the station. The worker always starts working on the first vehicle at position zero, left border of the station. The second vehicle is already at position $2=9-c$ when the worker completes working on the first vehicle. Note that the next vehicle enters the station borders a cycle time after the current vehicle enters the station borders. The tasks of the second vehicle are completed when the third vehicle has just entered the station. The worker has 2 TU of idle time when the tasks of the third vehicle are completed. The worker starts working on the fifth vehicle on position 2 and the processing time of the fifth vehicle is 9 TU which causes a work overload of 1 TU, $2+9-l_{k}=1$. The job of processing these five vehicles could not be completed within the station borders but with the help of a utility worker, we assume that the job is completed at the station border at the cost of 1 TU work overload. The worker will continue working on the sixth vehicle at position $3=l_{k}-c$, and this process continues.

We note that the vehicles usually go through the body shop and paint shop in the scheduled sequence before the assembly process. Hence, the failed vehicles must be pulled out of the sequence, and its position cannot be compensated, i.e., resequencing is not an option. It is assumed that each vehicle $v\in V$ (with its unique mix of configurations) has a failure probability $f_{v}$, and failures are independent of each other. The failures are related to the specifications of the vehicle, e.g., the increased production rate of EVs may induce higher failure rates, or painting a vehicle to a specific color may be more problematic. In our numerical experiments in Section 5, we estimate the failure probabilities from the historical data by doing feature analysis and using logistic regression.

In this section, first, we provide a two-stage stochastic program for our problem. Next, we discuss improvements of the proposed formulation.

Motivated by the dynamics of an MMAL, we formulate our problem as a two-stage stochastic program. The sequence of vehicles is decided in the first stage (here-and-now), before the car failures are realized. Once the car failures are realized, the work overload is minimized by determining the second-stage decisions (wait-and-see) given the sequence. First-stage decisions are determined by assigning each vehicle to a position such that the expected work overload in the second stage is minimized.

To formulate the problem, suppose that various realizations of the car failures are represented by a collection of finite scenarios $\Omega$. As each vehicle either exists or fails at a scenario, we have a total of $2^{|V|}$ scenarios. We let $\Omega=\{\omega_{1},\ldots,\omega_{2^{|V|}}\}$, with $\omega$ indicating a generic scenario. To denote a scenario $\omega$, let $e_{v\omega}=1$ if vehicle $v$ exists and $e_{v\omega}=0$ if vehicle $v$ fails at scenario $\omega\in\Omega$. We can then calculate the probability of scenario $\omega$ as $\rho_{\omega}=\prod_{v=1}^{|V|}f_{v}^{1-e_{v\omega}}(1-f_{v})^{e_{v\omega}}$ such that $\sum_{\omega\in\Omega}\rho_{\omega}=1$, where $f_{v}$ denotes the failure probability of vehicle $v\in V$.

A two-stage stochastic program for the full-information problem, where all possible realizations are considered, is as follows:

where

In the first-stage problem (1), the objective function represents the expected work overload, i.e., the cost associated with the second-stage problem. Constraint sets (1b) and (1c) ensure that exactly one vehicle is assigned to each position and each position has exactly one vehicle. respectively. Constraint set (1d) presents the domain of the binary first-stage variable. The second-stage problem (2) minimizes the total work overload throughout the planning horizon, given the sequence and scenario $\omega\in\Omega$. Note that $T_{\omega}$ denotes the set of positions of non-failed vehicles at scenario $\omega\in\Omega$, which is obtained by removing failed vehicles. Constraint set (2b) determines the processing time $b_{kt}$ at station $k$ at position $t$. The starting position and workload of the vehicles at each station are determined by constraint sets (2c) and (2d), respectively. The constraint set (2e) ensures that the first position starts at the left border of the station. Constraint set (2f) builds regenerative production planning, in other words, the first position of the next planning horizon can start at the left border of the station. The constraint set (2g) defines the second-stage variables as continuous and nonnegative.

Several remarks are in order regarding the proposed two-stage stochastic program (1) and (2). First, the number of decision variables and the set of constraints in (2) are scenario-dependent, as the valid positions $T_{\omega}$ are obtained based on the failure scenario $\omega\in\Omega$. Second, the proposed two-stage stochastic program (1) and (2) has a simple recourse. That is, once the sequence is determined and the failures are realized, the work overloads are calculated from the sequence of the existing vehicles, without resequencing.

In the remainder of this section, we first provide two modified models for the second-stage problem so that the number of decision variables and the set of constraints are no longer scenario-dependent. Then, we provide two monolithic MILP formulations for the deterministic equivalent formulation (DEF) of the two-stage stochastic program of MMS with stochastic failures.

For each scenario, we modify the second-stage problem by updating the processing times of failed vehicles instead of removing the failed vehicles. In Figure 3, we demonstrate how the original model (2) and modified models represent car failures. To this end, we consider the example given in Figure 2 and assume that the vehicle at the second position fails. In the original model, the failed vehicles are removed from the sequence and the succeeding vehicles moved forward (Figure 3(a)). The proposed modified models, referred to as standard model and improved model, are explained below. In Section 5.2, we discuss the impact of these modified models on the computational time and solution quality.

Standard Model: In a preprocessing step, the processing time of vehicle $v$ is set to zero for all stations if the vehicle fails in scenario $\omega$. Accordingly, since this modification is conducted by adding uncertainty to the processing times, the scenario index $\omega$ is added to the processing time. That is, $e_{v\omega}=0\Rightarrow p_{kv\omega}=0\enspace k\in K$. Based on this modification, the second-stage problem for a given scenario $\omega$ can be presented as

The objective function (3a) and the constraints (3b)–(3f) are the same as in formulation (2), except the set of positions and the length of the sequence are not scenario-dependent anymore. Constraint set (3g) guarantees that the starting position at station $k$ at position $t+1$ equals the starting position of position $t$ when the vehicle assigned to position $t$ fails. Constraint set (3h) assures that the regenerative production planning is kept in case the vehicle at the end of the sequence fails. The parameter $\beta$ is calculated in a way that $\beta_{k}b_{ktn}>z_{ktn}$ which makes constraint sets (3g) and (3h) non-effective for the positions that have existing vehicles. Hence, $\beta_{k}$ equals the maximum possible starting position divided by the minimum processing time for station $k$, $\beta_{k}=\frac{l_{k}-c}{\min_{v\in V}\{p_{kv}\}}>0$. Note that the processing times in this calculation are the actual processing times before the preprocessing step. Also, $\beta_{k}$ is well-defined as the the minimum processing time is strictly greater than zero. Figure 3(b) demonstrates that in the standard model, the processing time of the second vehicle is set to zero, so the operator starts working on the third vehicle at position two where the operator was going to start working on the second vehicle if it had not failed.

Improved Model: In order to reduce the size of the standard model, we modify this model as follows. During the preprocessing step, the processing time of vehicle $v$ is set to the cycle time for all stations if a vehicle fails at scenario $\omega$. Let us refer to the vehicles with processing time equal to cycle time for all stations as “neutral” because these vehicles do not have any impact on the schedule in terms of work overload (see Proposition 1 and its proof). In other words, we transform failed vehicles into neutral vehicles, i.e., $e_{v\omega}=0\Rightarrow p_{kv\omega}=c\enspace k\in K$.

As a result of Proposition 1, constraints (3g) and (3h) can be removed from the standard model. Hence, the problem size is reduced. Figure 3(c) contains an illustration for Proposition 1. The second vehicle becomes neutral when its processing time is set to cycle time so that the third vehicle starts at the same position as the second vehicle.

Using the standard or improved model, the DEF for MMS with stochastic failures can be obtained by adding the first-stage constraints (1b)–(1d) to the corresponding second-stage formulation, and by adding copies of all second-stage variables and constraints. We skip the details for brevity.

For the ease of exposition, we consider an abstract formulation of the two-stage stochastic program presented in Section 3.2 as follows:

where $x$ denotes the first-stage decisions variables and $X:=\{x\in\{0,1\}^{|V|\times|T|}\colon Ax=b\}$ is the feasible region of decision variables $x$, i.e., the set of points satisfying constraints (1b) - (1d). Moreover, we represent the second-stage problem for the standard or improved model, presented in Section 3.2, as

where $y$ represents the second-stage decision variables and $\xi_{\omega}=(h_{\omega},T_{\omega})$. The expectation of the recourse problem becomes $\mathbb{E}[Q(x,\xi_{\omega})]=\sum_{\omega\in\Omega}\rho_{\omega}Q(x,\xi_{\omega})$.

The L-shaped method is a procedure that has been successfully used to solve large-scale two-stage stochastic programs. Note that for any $\omega\in\Omega$, function $Q(x,\xi_{\omega})$, defined in (5), is convex in $x$ because $x$ appears on the right-hand side of constraints. Hence, we propose to iteratively construct its underestimator. To this end, for each $\omega\in\Omega$ and a given first-stage decision $x\in X$, we consider a subproblem that takes the form of (5). Moreover, we create a relaxed master problem, which contains a partial, but increasingly improving, representation of $Q(x,\xi_{\omega})$, for each $\omega\in\Omega$, through the so-called Benders’ cuts. Recall that our proposed two-stage stochastic programs have a relative complete recourse, that is, for any first-stage decision $x$, there is a feasible second-stage variable $y$. Thus, an underestimator of $Q(x,\xi_{\omega})$ can be constructed by only the so-called Benders’ optimality cuts.

We now describe more details on our proposed L-shaped algorithm. We form the relaxed master problem for formulation (4) and (5) as follows:

where the auxiliary variable $\theta_{\omega}$ approximates the optimal value of the second-stage problem under scenario $\omega\in\Omega$, i.e., $Q(x,\xi_{\omega})$, through cuts $\theta_{\omega}\geq G^{\iota}_{\omega}x+g^{\iota}_{\omega}$ formed up to iteration $l$.

Let $(\hat{x}^{\iota},\hat{\theta}^{\iota})$ be an optimal solution to the relaxed master problem (6). For each scenario $\omega\in\Omega$, we form a subproblem (5) at $\hat{x}^{\iota}$. Suppose that given $\hat{x}^{\iota}$, $\hat{\pi}_{\omega}^{\iota}$ denotes an optimal dual vector associated with the constraints in (5). That is, $\hat{\pi}_{\omega}^{\iota}$ is an optimal extreme point of the dual subproblem (DSP)

where $\pi_{\omega}$ is the associated dual vector. Then, using linear programming duality, we generate an optimality cut as

where $G^{\iota}_{\omega}=-(\hat{\pi}_{\omega}^{\iota})^{\top}T_{\omega}$ and $g^{\iota}_{\omega}=(\hat{\pi}_{\omega}^{\iota})^{\top}h_{\omega}$.

Our proposed L-shaped algorithm iterates between solving the relaxed master problem (6) and subproblems (5) (one for each $\omega\in\Omega$) until a convergence criterion on the upper and lower bounds is satisfied. This algorithm results in an L-shaped method with multiple cuts.

In order to exploit the specific structure of the MMS problem and to provide improvements on the dual problem, let us define variables $\pi^{sp}$, $\pi^{wo}$, $\pi^{fs}$, $\pi^{ch}$, $\pi^{sf}$, and $\pi^{cf}$ corresponding to starting position constraints (3c), work overload constraints (3d), first station starting position constraints (3e), regenerative production planning constraints (3f), starting position of the vehicles following a failed vehicle (3f), and regenerative production planning with failed vehicles (3g), respectively. The DSP for scenario $\omega\in\Omega$ at a candidate solution $\hat{x}^{\iota}$, obtained by solving a relaxed master problem, can be formulated as follows:

We provide improvements to the dual problem in several ways. The dual variables $\pi_{fs}$ and $\pi_{cf}$ are removed since the corresponding subproblem constraints (3f) and (3g) are eliminated in the improved model. The dual variables $\pi^{fs}$ and $\pi^{ch}$ are not in the objective function and are unrestricted, which means that we can remove these variables and the constraints with those variables from the formulation without altering the optimal value of the problem. In our preliminary computational studies, we improved the dual subproblem by removing these variables. However, we observed that most of the DSPs have multiple optimal solutions, and as the number of vehicles and stations increase, it is more likely to have multiple optimal solutions. This naturally raises the question of what optimal dual vector provides the strongest cut, if we add only one cut per iteration per scenario. One can potentially add the cuts corresponding to all optimal dual extreme points, however, this results in an explosion in the size of the relaxed master problem after just a couple of iterations. While there is no reliable way to identify the weak cuts (Rahmaniani et al., 2017), we executed experiments in order to find a pattern for strong cuts. Our findings showed that adding the cut corresponding to the optimal dual extreme point with the most non-zero variables results in the fastest convergence. Thus, we added an $\ell_{1}$ regularization term to the objective function of the DSP, hence, the new objective is encouraged to choose an optimal solution with the most non-zero variables. Accordingly, we propose an improved DSP formulation as follows:

MMS is an NP-hard problem, and stochastic failures of products (cars) increases the computational burden of solving the problem drastically. Hence, it is essential to create efficient heuristic procedures in order to solve industry-sized problems. In this section, we provide a fast and easy-to-implement greedy heuristic to find a good initial feasible first-stage decision (i.e., a sequence of vehicles) and an efficient tabu search (TS) algorithm to improve the solution quality. Although all $|V|!$ vehicle permutations are feasible, the proposed greedy heuristic aims to find a good initial feasible solution. To achieve this, a solution is generated for the deterministic counterpart of the proposed MMS problem, which excludes vehicle failures. We refer to this problem as the one-scenario problem, since the corresponding problem has a single scenario with no failed vehicles. Assuming that the failure probability of each vehicle is less than or equal to 0.5, the scenario with no failed vehicles has the highest probability. Once such a feasible sequence of vehicles is generated, the TS algorithm improves this solution in two parts: first, over the one-scenario problem, and then, over the full-information problem.

In (4), it is assumed that the probability of each scenario is known a priori, which may not hold in practice. In addition, the exponential growth of the number of scenarios causes an explosion in the size of the stochastic program. Hence, we utilize the SAA approach to tackle these issues. Consider the abstract formulation (4). The SAA method approximates the expected value function with an identically and independently distributed (i.i.d) random sample of $N$ realizations of the random vector $\Omega_{N}:=\{\omega_{1},\ldots,\omega_{N}\}\subset\Omega$ as follows:

The optimal value of (12), $z_{N}$, provides an estimate of the true optimal value (Kleywegt et al., 2002). Let $\hat{x}_{N}$ and $x^{*}$ denote an optimal solution to the SAA problem (12) and the true stochastic program (4), respectively. Note that $\mathbb{E}[Q(\hat{x}_{N},\xi_{\omega})]-\mathbb{E}[Q(x^{*},\xi_{\omega})]$ is the optimality gap of solution $\hat{x}_{N}$, where $\mathbb{E}[Q(\hat{x}_{N},\xi_{\omega})]$ is the (true) expected cost of solution $\hat{x}_{N}$ and $\mathbb{E}[Q(x^{*},\xi_{\omega})]$ is the optimal value of the true problem (4). A high quality solution is implied by a small optimality gap. However, as $x^{*}$ (and hence, the optimal value of the true problem) may not be known, one may obtain a statistical estimate of the optimality gap to assess the quality of the candidate solution $\hat{x}_{N}$ (Homem-de Mello and Bayraksan, 2014). That is, given that $\mathbb{E}[z_{N}]\leq\mathbb{E}[Q(x^{*},\xi_{\omega})]$, we can obtain an upper bound on the optimality gap as $\mathbb{E}[Q(\hat{x}_{N},\xi_{\omega})]-\mathbb{E}[z_{N}]$. We employ the multiple replication procedure of Mak et al. (1999) in order to assess the quality of a candidate solution by estimating an upper bound on its optimality gap. A pseudo-code for this procedure is given in Algorithm 1. We utilize MRP, in Section 5, to assess the quality of solutions generated by different approaches.

In addition, we propose an MRP integrated SAA approach for candidate solution generation and quality assessment, given in Algorithm 2. On the one hand, a candidate solution is generated by solving an SAA problem with a sample of $N$ realizations. Then, we use the MRP to estimate an upper bound on the optimality gap of the candidate solution. If the solution is $\epsilon$-optimal, i.e., estimated upper bound on its optimality gap is less than or equal to a $\epsilon$ threshold, the algorithm stops. Otherwise, the sample size increases until a good quality solution is found. The algorithm returns a candidate solution and its optimality gap.

We end this section by noting that each of the DEF, presented in Section 3.2, the L-shaped algorithm, presented in Section 4.1, and the heuristic algorithm, presented in Section 4.2, can be used to solve the SAA problem and obtain a candidate solution. However, the probability of scenarios $\omega\in\Omega$, $\rho_{\omega}$, must change in the formulations so that it reflect the scenarios in a sample $\Omega_{N}$. Let $\hat{N}$ and $n_{\omega}$ represent the set of unique scenarios in $\Omega_{N}$ and the number of their occurrences. Thus, in the described DEF, L-shaped algorithm, and TS algorithm, $\sum_{\omega\in\Omega}\rho_{\omega}(\cdot)$ changes to $\frac{1}{N}\sum_{\omega\in\Omega_{N}}(\cdot)$ or equivalently, $\frac{1}{N}\sum_{\omega\in\hat{N}}n_{\omega}(\cdot)$. Accordingly, in the L-shaped method, we generate one optimality cut for each unique scenario $\omega\in\hat{N}$ by solving $|\hat{N}|$ number of subproblems at each iteration.

We generated real-world inspired instances from our automobile manufacturer partner’s assembly line and planning information. As given in Table 3, we generated three types of instances: (1) small-sized instances with 7-10 vehicles to assess the performance of L-shaped algorithm, (2) medium-sized instances with 40 vehicles to assess the performance of the TS algorithm for the one-scenario problem, (3) large-sized instances with 200, 300, and 400 vehicles to evaluate the performance of the TS algorithm. All instances have five stations, of which the first one is selected as the most restrictive station for EVs, the battery loading station. The rest are selected among other critical stations that conflict with the battery loading station.

The cycle time $c$ is 97 TU, and the station length $l$ is 120 TU for all but the battery loading station, which is two station lengths, 240 TU. The information about the distribution of the processing times is given in Table 4. It can be observed that the average and maximum processing times for each station are lower than the cycle time and the station length, respectively. Moreover, the ratio of the EVs is in the range of [0.25, 0.33] across all instances.

We derived the failure rates from six months of historical data by performing predictive feature analysis on vehicles. Based on the analysis, two groups of vehicles are formed according to their failure probabilities, low-risk and high-risk vehicles, whose failure probabilities are in the range of [0.0, 0.01] and [0.2, 0.35], respectively. The failure probability is mostly higher for recently introduced features, e.g., the average failure probability of EVs is 50% higher than that of other vehicles. High-risk vehicles constitute [0.03, 0.05] of all vehicles. However, this percentage increases to [0.15, 0.25] for the small-sized instances in order to have a higher rate of failed vehicles. We note that the failures are not considered for the medium-sized instances since these instances are used for only the one-scenario problem, which does not involve failures by definition.

The number of failure scenarios, $2^{|V|}$, increases exponentially in the number of vehicles. Thus, we generated an i.i.d random sample of $N$ realizations of the failure scenarios; hence, formed a SAA problem. For each failure scenario and vehicle, we first chose whether the vehicle was high risk or low risk (based on their prevalence). Then, depending on being a high-risk or low-risk vehicle, a failure probability was randomly selected from the respective range. Finally, it was determined whether the vehicle failed or not. In order to have a more representative sample of scenarios for large-sized instances, no low-risk vehicle was allowed to fail at any scenario.

For each parameter configuration, we generated 30 instances. The vehicles of each instance were randomly selected from a production day, respecting the ratios mentioned above. The algorithms were implemented in Python 3. For solving optimization problems we used Gurobi 9.0. The time limit is 600 seconds for all experiments unless otherwise stated. We run our experiments on computing nodes of the Clemson University supercomputer. The experiments with the exact solution approach were run on nodes with a single core and 15 GB of memory, and the experiments with the heuristic solution approach were on on nodes with 16 cores and 125 GB of memory.

In this section, we present results for the solution quality and computational performance of the L-shaped algorithm. We used MRP scheme, explained in Section 4.3, to assess the solution quality. We also compared the computational performance of the L-shaped algorithm with that of solving the DEF. We present the results for 120 small-sized instances consisting of 7 to 10 vehicles. We do not present the results for large-sized instances as our preliminary experiments showed that the number of instances that could be solved to optimality decreases drastically.

We also point that instead of solving a relaxed master problem to optimality at each iteration of the L-shaped algorithm, one can aim for just obtaining a feasible solution $\hat{x}\in X$. This may result in saving a significant amount of computational time that would be otherwise spent on exploring solutions that are already eliminated in previous iterations. This kind of implementation, referred to as branch-and-Benders-cut (B&BC), is studied in the literature, see, e.g., (Hooker, 2011; Thorsteinsson, 2001; Codato and Fischetti, 2006). In our implementation of our proposed L-shaped algorithm, we used Gurobi’s Lazy constraint callback to generate cuts at a feasible integer solution found in the course of the branch-and-bound algorithm.

In this section, we present results for the solution quality and computational performance of the TS algorithm. The solution quality is evaluated by employing the MRP scheme, explained in Section 4.3. We also assess the computational performance of the TS algorithm in various aspects.

We set the operator selection probabilities (weights) based on our preliminary experiments. The weights of swap, forward insertion, backward insertion, and inversion are 0.45, 0.10, 0.15, 0.30, respectively. We set the time limit for the one-scenario problem to 10 seconds, i.e., $\tau_{one}=10$ seconds, which leaves $\tau_{full}=590$ seconds. Finally, based on the results of the quality assessment, we set the sample size $N$ to 1000 for the computational performance assessments.