Mixed-model moving assembly line material placement optimization for a shorter time-dependent worker walking time

Abstract

Car mass production commonly involves a moving assembly line that mixes several car models. This requires plenty of material supplies at the line side, but available space is scarce. Thus, material is placed apart from ideal positions. Then, picking it up involves walking along the line. This time is non-productive and can encompass 10–15% of total production time. Thus, it is important to estimate and minimize it during production planning. However, the calculations are difficult because the conveyor continuously moves. Therefore, most literature bounds walking time by a constant, but this discards valuable potential. To better approximate it, we use a time-dependent V-shaped function. A comparison indicates that for a majority of instances, constant walking time estimates of 95% confidence are at least 51% higher. Then, we introduce a model to optimize material positions such that the model-mix walking time is minimized. This poses an NP-hard sequencing problem with a recursive and nonlinear objective function. Our key discovery is a lower bound on the objective of partial solutions, established by a Lagrangian relaxation that can be solved in quadratic time. Resulting branch and bound based algorithms allow to quickly and reliably optimize up to the largest real-world sized instances.

1 Introduction

Mass-production of cars was particularly made possible by the moving assembly line. It continuously transports the work pieces from worker to worker. Productivity is highest if the assembly operations are divided equally between all workers. This is an NP-hard bin packing type optimization problem (Álvarez-Miranda and Pereira, 2019) that is called assembly line balancing (Salveson, 1955); surveys are found in Battaïa and Dolgui (2013, 2022), Becker and Scholl (2006), Boysen et al. (2022). Solutions can improve with shorter time buffers given more accurate time estimates, e.g., for a worker’s workload when assigning operations.

A significant time sink can be a worker’s walking time to fetch parts from the line side, amounting for 10–15% of total production time at a major German car manufacturer (Scholl et al., 2013). This time can hardly be estimated with a constant: it is highly variable because the distance is time-dependent. A method to estimate it is introduced in Sedding (2020a) in production of a single model. Walking time is then minimized by placing material supplies at optimal positions along the line side (the operation sequence is fixed). This optimized and more exact walking time estimate allows to better gage the worker’s workload during assembly line balancing, potentially increasing overall line utilization. To employ the estimate in interactive planning software, fast compute times are essential, which are provided by heuristics in Sedding (2020a) with a median runtime of 0.002s. In this paper, we adapt this approach to a model-mix production, which is common in car assembly (Sternatz, 2014).

An assembly worker fetches needed material for an assembly operation by walking to the respective material box. In the typical moving car assembly line, boxes are located at the line side, along which the workpiece moves. Ideally, each box is located close to the workpiece’s position at access time, as this minimizes walking time. However, if the line side space is scarce (Bautista and Pereira, 2007; Bukchin and Meller, 2005; Boysen et al., 2015), it is usually necessary to place a box apart from its ideal location along the line.

The moving workpiece’s distance to a box and the resulting two-way walking time can be modeled by a convex function of time (Sedding, 2020a). However, most models for assembly line planning in the literature are restricted to constant processing time estimates, just like in the classic assembly line balancing problem (Salveson, 1955). Although sequence-dependent nonproductive processing times are included in Andrés et al. (2008), Scholl et al. (2013), Esmaeilbeigi et al. (2016), they cannot be applied to precisely plan walking times at assembly lines with a moving workpiece: Overly large safety factors are needed to upper bound the walking time with a constant. This gives away potential to increase the utilization of an assembly line. In this study, we test the potential of time-dependent walking time estimates.

When assigning operations to a worker, a realistic walking time estimate needs to take possible local optimizations into account. A minimization of a worker’s walking time can be approached in two major ways for a given set of assembly operations. On the one hand, it is possible to optimize the operation sequence (Jaehn and Sedding, 2016; Sedding, 2020b). On the other, the positioning of the respective material boxes can be adjusted (Klampfl et al., 2006; Finnsgård et al., 2011; Sedding, 2017, 2020a). This requires to fix the operation sequence to determine a walking time optimized box placement. Finnsgård et al. (2011) describe a manual optimization and report a walking distance reduction by 52%. Schmid et al. (2021) optimize the placement with a mixed-integer programming approach that considers variable walking costs during material placement, but ignores the workpiece’s continued movement while the worker walks. Klampfl et al. (2006) take this movement into account, using Euclidean distances to calculate time-dependent walking times. In three approaches, they consider placement of boxes along the line. First with overlaps, then without overlaps, and finally with stacking atop each other. Nonlinear programming is used to heuristically find placements. However, they record rather long compute times already for a small instance of five material boxes. A one-dimensional walking time model yields a significantly quicker optimization in Sedding (2017, 2020a, c) for the case of single model placement.

In this paper, we adapt the material placement optimization in Sedding (2020a) to the mixed-model moving assembly line. A mixed-model assembly shares the same line for several product models (Thomopoulos, 1967). Then, the production can better adapt to varying demands of each model. This production method is standard in car assembly (Sternatz, 2014). Line side space can be even more scarce if further material is required for each model (Boysen et al., 2015). In the mixed-model setting, the objective is to minimize total processing time weighted by model share. This allows for longer processing times on some product models, especially rarer ones, provided this can be compensated for on other product models. Allowance is provided by the worker floating up- or downstream the line. However, an overload over the course of several cycles cannot be compensated anymore as increasing walking times exacerbate the overload. Such situations must be prevented during planning, in particular of the production sequence, cf. Boysen et al. (2009c) for a survey.

Our paper’s contribution lies in placing material boxes at the line side for a worker’s assembly operations over a mix of product models such that the model-mix weighted time-dependent walking times to gather material are minimized:

• First, we display all assumptions and introduce an optimization model in Sect. 2. In comparison to Sedding (2020a, c), it permits multiple models and an offset for the placement area (or indirectly the start time).

• We observe that the problem entails a recursive, nonlinear objective function, which impedes evaluation of partial solutions and renders incremental solution procedures difficult. We provide a proof of strong NP-hardness for any number of product models (Sect. 3).

• A mixed-integer program is derived from the single-model case in Sedding (2020a, c).

• The main technical contribution of our paper is a Lagrangian relaxation of the mathematical program, for which we find a partition into two independent subproblems, each of which is solved in quadratic time. This results in a fast lower bound for partial solutions (Sect. 5). This lower bound is the cornerstone of a branch and bound algorithm that incrementally places the boxes. A truncated branch and bound search provides a fast heuristic (Sect. 6).

• Finally, a numerical experiment shows the effectiveness of the approaches (Sect. 7). The total runtime of the best heuristic is negligibly small: for the largest instances, the median runtime is 0.071 s. Such a runtime is suitable for a use in interactive planning software. Potential savings are high in comparison to constant walking time estimates as well as intuitive placements (Sect. 7.5).

• Our model is the first in the literature to efficiently consider and minimize a time-dependent model of walking times for material placement at the mixed-model moving assembly line. As this production mode is standard for car production, our approach has a far-reaching applicability.

Preliminary versions of this paper are available in this special issue’s workshop proceedings (Sedding, 2021) and in the author’s dissertation (Sedding, 2020c, Chapter 6).

Table 1 The model assumptions equal the study of the single-model case in Sedding (2020a, c) except for A3 switching to a model-mix and adding A16

2 Modeling

This section presents the studied optimization problem $\mathscr {P}m$ for placing boxes for a model-mix of m products. After a formal definition of all parameters in a $\mathscr {P}m$ instance, its prerequisites are introduced. The definition of $\mathscr {P}m$ follows.

Fig. 1

On planning the line side placement, boxes are positioned in a sequence at the line side in $[V,W)$, displayed in the figure’s top row. Shown below are assembly operations (jobs) in fixed sequences for three product models. The horizontal axis represents the time; it shows assembly operations performed by the worker while a workpiece moves along the line over time t. Note that the time scale equals the space scale. At the start of a production cycle, the worker is at time and spatial position 0. Blue arrow lines indicate, exemplary for model 3, the worker’s longitudinal movement along the line during one production cycle, alternating between fetching material (walking the curved line) and assembling, during both of which the workpiece moves. Job (3, 1) starts at time 0; with walking time $f_{3,1}(0)$ (indicated by the striped area, to &from the box at $\pi _{3,1}$) and assembly time $\ell _{3,1}$, the job completes at $C_{3,1}(0)$. Job completion times, calculated recursively, are labeled in blue color. A production cycle of model 3 is completed at time (and position) $C_{3,\max }$. The objective is to minimize the average completion time over all models by placing the boxes such that walking times are minimal.

Note that $\mathscr {P}m$ builds upon the optimization problem in Sedding (2020a, c) for a single product assembly $m=1$. While the single-model case is extended to a mixed-model production, all other assumptions remain unchanged. A list of assumptions made for $\mathscr {P}m$ is gathered in Table 1.

Definition 1

An instance of $\mathscr {P}m$ is given by

Before the optimization problem $\mathscr {P}m$ can be defined, we specify how boxes can be placed, see how walking time for fetching part of a boxes can be modeled, and define how to calculate the makespan (completion time) of a production cycle encompassing walking and assembly time. The result is visualized for an example instance in Fig. 1.

2.1 Box placement

Necessary material of an assembly operation is provided in a dedicated box $(i,j)\in J$. It takes up a certain box width $w_{{i,j}}$ along the line side, expressed as a share of the available width (in the dimension along the line). We may assume $w_{{i,j}}\in \mathbb {N}$ by scaling the coordinate system accordingly.

Let set J denote the set of boxes for the the worker. The boxes are placed side-by-side without overlaps. Given a scarce line side space, we assume there are no gaps between boxes. Then, the boxes are placed in a contiguous section at the line side: the space between $V$ and $W=V+\sum _{(i,j)\in J}w_{{i,j}}$.

Then, a box placement 

$$\begin{aligned} \{\pi _{{i,j}}\}_{(i,j)\in J} \end{aligned}$$

(1)

states, for each box $(i,j)\in J$, a rational-valued position ${V}\le \pi _{{i,j}}\le W-w_{{i,j}}$ to place it on interval $[\pi _{{i,j}},\pi _{{i,j}}+w_{{i,j}})$ at the line side such that $\bigcup _{(i,j)\in J} [\pi _{{i,j}},\pi _{{i,j}}+w_{{i,j}})=[{V},W)$. Note that attaining a box placement is equivalent to finding a sequence for the boxes along the line side.

2.2 Walking time

The walking time calculation in Sedding (2020a) is briefly described in the following. In this model, the worker only walks along the moving assembly line; movement orthogonal to the line can be ignored. Three walking strategies are covered: Always walking on the fixed floor, atop the conveyor’s moving floor, or whichever is better in the current movement direction.

For walking up to a box $(i,j){\in J}$, the worker leaves the workpiece at a certain time t, arrives at the box’s position $\pi _{{i,j}}$, and then returns to the workpiece. All the while, the workpiece continues to move. This gives several movement equations, which are solved with a closed formula. Then, the walking time is represented by the piecewise-linear function

$$\begin{aligned} f_{{i,j}}(t)=\max \{-a\cdot (t-\pi _{{i,j}}),\,b\cdot (t-\pi _{{i,j}})\} \end{aligned}$$

(2)

where $\pi _{{i,j}}$ is the box position encoded in time units, and $0\le a\le 1$ and $b\ge 0$ are two slopes that depend on the worker’s velocity and walking strategy. See also Fig. 2.

If the current time equals the box position (case $t=\pi _{{i,j}}$), then the walking time is minimum, which occurs if the workpiece just passes by the box position. Otherwise, the walking time increases linearly. If the workpiece has not yet passed the box position (case $t<\pi _{{i,j}}$), then the walking time corresponds to $-a\cdot (t-\pi _{{i,j}})$, otherwise it is $b\cdot (t-\pi _{{i,j}})$.

The two slopes relate the workpiece’s velocity $v_{\text {conveyor}}$ to the worker’s walking velocity $v_{\text {worker}}>v_{\text {conveyor}}$. For example, if the worker walks on a non-moving floor beside the workpiece, then $a=2/(v+1)$ and $b=2/(v-1)$ where $v=v_{\text {worker}}/v_{\text {conveyor}}$. The same slope values are attained if the worker walks on floor plates that move together with the workpiece. If a worker can always choose the best of both options, then the slopes are $a=(2v+1)/(1+v)^2$ and $b=(2v+1)/v^2$. Note that this walking strategy yields, if $v=13.6$ as in Klampfl et al. (2006), a walking time reduction by 3.5% if $t<\pi _{{i,j}}$, else 4.1%, see Sedding (2020a).

2.3 Makespan calculation

Before an assembly operation can be processed, a walking time occurs to its distinct box $(i,j)\in J$. Together, these two constitute a job, which is also denoted by (i, j). Then, job (i, j) consists of two parts: the nonproductive walking time function $f_{{i,j}}$ in (2), and after that, a productive assembly time $\ell _{{i,j}}\ge 0$ that is a constant nonnegative rational number. Together, they form the job’s processing time $p_{{i,j}}(t)=f_{{i,j}}(t)+\ell _{{i,j}}$. Starting at time t, the job completes at $C_{{i,j}}(t)=t+p_{{i,j}}(t)$. Substituting all components, the job’s completion time is expressed by

$$\begin{aligned} C_{{i,j}}(t)&=t+\ell _{{i,j}}+\max \{-a\cdot (t-\pi _{{i,j}}),\,b\cdot (t-\pi _{{i,j}})\} . \end{aligned}$$

(3)

Each product model $i\in {{I}}$ requires processing of a certain number $n_i$ of jobs in a fixed sequence denoted by

$$\begin{aligned} (i,1),(i,2),\dots ,(i,n_i) \end{aligned}$$

where $n_i$ is the number of jobs in model i. Then, the last completion time $C_{i}^{{\max }}$ in a model i is the composition of job completion times (3), starting the first job (i, 1) at time 0:

$$\begin{aligned} C_{i}^{{\max }}=C_{i,n_i}(\cdots \,C_{i,2}(C_{i,1}(0)))\cdots ). \end{aligned}$$

(4)

Note that a start of the first job at $t_{\min }\ne 0$ is attained by shifting the box placement interval by $-t_{\min }$.

Remark 1

The job sequence is fixed here. As mentioned in the introduction, it is also possible to optimize walking time by permuting the job sequence (Sedding, 2020b, c). This belongs to the field of time-dependent scheduling (Gawiejnowicz, 2020b, a). Related piecewise-linear convex $C_{{i,j}}$ functions are described in Farahani and Hosseini (2013), Kawase et al. (2018), Kononov (1998), while a recent survey is found in Sedding (2020b).

Fig. 2

Visualization of the walking time function $f_{{i,j}}$ of job ${(i,j)\in J}$ in (2) as a function of start time t

2.4 Optimization problem

Minimizing the overall walking time in the mixed-model setting is equal to minimizing the total weighted makespan (completion time) over all product models (Klampfl et al., 2006). This can be attained by a sum of each model i’s last completion time $C_{i}^{{\max }}$ weighted by production share $r_i$. This objective is minimized in the studied optimization problem.

Definition 2

(Problem $\mathscr {P}m$) Given a $\mathscr {P}m$ instance (see Definition 1), minimize the weighted average makespan

$$\begin{aligned}\phi =\sum _{i\in {{I}}} r_i C_{i}^{{\max }}\end{aligned}$$

by determining a box placement $\{\pi _{i,j}\}_{(i,j)\in J}$, which places the boxes, each represented by box width $w_{i,j}$, in a sequence in $[{V},W)$, see (1). This yields, for model $i\in {{I}}$, makespan

$$\begin{aligned}C_{i}^{{\max }}=(C_{i,n_i}\circ \cdots \circ C_{i,2}\circ C_{i,1})(0),\end{aligned}$$

which composes the completion times of model i’s job sequence according to (4). By (3), the completion time $C_{{i,j}}$ of job $(i,j)\in J$ as a function of job start time t is given by

$$\begin{aligned}{C_{{i,j}}(t)=t+\ell _{{i,j}}+\max \{-a\cdot (t-\pi _{{i,j}}),\,b\cdot (t-\pi _{{i,j}})\}.}\end{aligned}$$

In the single-model case, the makespan cannot be larger than the cycle time. Every cycle, a new workpiece arrives, hence a larger makespan would require a line stoppage (or a work overload situation). Herein lies a key advantage of model-mix production: If the weighted average makespan $\phi$ is not above the cycle time, then there is no work overload (on average). Models with a higher makespan let the worker float downstream, others upstream.

Hence, some models may be allowed a makespan longer than the cycle time. However, note that such an overload situation would affect the next production cycle’s start time. Then, walking times are affected. If it occurs repeatedly, however, walking times can grow quickly as the worker is driven away from the boxes. If the worker has floated downstream behind the corresponding boxes, then any delay increases the completion time by a factor of $(1+b)^{n}$ for n remaining jobs, which follows from Sedding (2020b, Corollary 2). In excess, the worker needs assistance and/or the line needs to stop. Such overload situations must prevented in planning of the production sequence.

Although the makespan should equal the cycle time, it can well be different to the placement interval’s end W. Moreover, the model allows a nonzero placement interval start $V$. In both cases, the placement area is incongruent to the assembly station, covering only a part, or extending out of it. This models that the placement area is offset to the assembly station. This is required, e.g., when subdividing the line-side space into different sized regions, which can benefit the overall assembly line balance. A nonzero start time of the first job can depict floating of the worker up- or downstream the line. This is represented in our model by shifting the placement interval back by the same amount.

3 Computational complexity

The main difficulty of optimizing a box placement is caused by the recursive and nonlinear nature of the objective. For example, if the first job’s box is moved, then its walking time changes. As a consequence, succeeding jobs start at a different point in time. This time might be earlier or later than before. As each walking time function is nonlinear, the objective value changes nonlinearly. Moreover, the respective ideal box location of succeeding jobs changes. Then, the placement of these boxes needs further reoptimization.

To highlight the complexity of $\mathscr {P}m$, we prove that it is NP-hard in the strong sense for an arbitrary number of product models $m\ge 1$. For this, we perform a reduction from 3-Partition as defined in Garey and Johnson (1978), which is NP-complete in the strong sense (Garey and Johnson, 1975).

Definition 3

(3-Partition) Given a bound $B\in \mathbb {N}$ and 3z elements in multiset $X=\{x_1,\dots ,x_{3z}\}\subset \mathbb {N}$ with $B/4< {x} < B/2$ for $x\in X$, and $\sum _{x\in X} x=zB$. The question is: does there exist a partition of X into z disjoint multisets $A_{{i}}$, $i=1,\dots ,z$ with $\sum _{x\in A_{{i}}} x=B$?

Theorem 1

$\mathscr {P}m$ is NP-hard in the strong sense for arbitrary $m\ge 1$.

Proof

For $m=2$, we perform a reduction from 3-Partition as of Definition 3. This requires a decision version of $\mathscr {P}m$ that specifies a threshold value $\varPhi$ and asks if a solution exists with objective value $\phi \le \varPhi$.

A corresponding instance is constructed as follows. First, we freely choose any allowed nonzero slope $0<a\le 1$, $b>0$ value. For the two models ${{I}}=\{1,2\}$, we set production share $r_1=0$ and $r_2=1$. Hence, model 1 incurs no walking time in the objective function. However, it occupies space at the line side for its ${n_1=}{}3z$ jobs for which we let $w_{1,j}={x}_j$, $j=1,\dots ,3z$, and choose an arbitrary $\ell _{1,j}$ value. The second model has ${n_2=}{}z+1$ jobs. For each $j=1,\dots ,z+1$, we set $w_{2,j}=1$ and $\ell _{2,j}=B+1$. Finally, we set threshold value

$$\begin{aligned} \varPhi =\sum _{j=1,\dots ,z+1}l_{2,j}=\left( {B+1}\right) \left( {z+1}\right) . \end{aligned}$$

This instance’s objective value is, given the unilateral production shares, $\phi =C_{2,z+1}$. As $C_{2,z+1}\le \varPhi$, there is $\phi \le \varPhi \iff \phi =\varPhi$. Let us assume that $\phi =\varPhi$. This requires in the second product model for each $j=1,\dots ,z+1$ that $p_{2,j}=\ell _{2,j}$. This is the case if and only if the corresponding box is precisely positioned at $(B+1)\cdot j$. These boxes leave gaps of width B. Each gap is closed by the first product model’s boxes if and only if the 3-Partition instance can be solved. Hence, $\mathscr {P}m$ is NP-hard for $m=2$.

We generalize this reduction to $m>2$ by extending the instance with $m-2$ models ${{I}}'$. For each $i\in {{I}}'$, let $r_i=0$. Then, the last completion times of models ${{I}}'$ yield no impact on the objective function. Moreover, we introduce an arbitrary number of jobs for each added model $i\in {{I}}'$, each with the same box width $B+1$ and an arbitrary assembly time. Because these boxes are too wide to be placed between two adjacent boxes of the second product model for $\phi \le \varPhi$, they must be placed after the last one. Hence, they assert no effect on the box placement of the first and the second product model.

For $m=1$, strong NP-hardness is shown in Sedding (2020a). Concluding, a pseudopolynomial reduction of 3-Partition to $\mathscr {P}m$ exists for $m\ge 1$. $\square$

4 Mathematical programming

We describe $\mathscr {P}m$ solutions using mathematical programming, adapting the special case $m=1$ in Sedding (2020a, c) to $m\ge 1$. This includes a makespan calculation for each product model and changes the objective function to a sum of makespans weighted by production share. Then, we derive a mixed-integer program, reformulating box placement constraints.

4.1 Mathematical program

We introduce continuous variables $\pi _{i,j}$ as box position and $C_{i,j}$ as completion time of job $(i,j)\in J$. Then, a mathematical program for $\mathscr {P}m$ can be stated as:

$$\begin{aligned}&\text {minimize } \sum _{i\in {{I}}} r_i \,C_{i,n_i}\nonumber \\&\text {subject to}\nonumber \\&C_{i,0}={}0,\qquad \quad {i\in {{I}},} \end{aligned}$$

(5a)

$$\begin{aligned}&C_{i,j}\ge {}C_{{i,j-1}}+\ell _{i,j}-a\big (C_{i,j-1}-\pi _{i,k}\big ),(i,j)\in J, \end{aligned}$$

(5b)

$$\begin{aligned}&C_{i,j}\ge {}C_{i,j-1}+\ell _{i,j}+b\big (C_{i,j-1}-\pi _{i,k}\big ),(i,j)\in J, \end{aligned}$$

(5c)

$$\begin{aligned}&{\{\pi _{i,j}\}_{(i,j)\in J}\text { being a box placement.}} \end{aligned}$$

(5d)

Completion times are set recursively in constraints (5b) and (5c) starting with (5a). Constraint (5d) restricts the box position variables to a valid box placement as defined in (1).

We observe that attaining a valid box placement corresponds to a job sequencing problem on a single machine, which determines each job’s processing interval from start time to completion time. This is alike to the interval a box is placed upon; the difference being that the job’s interval is typically denoted by its end (the completion time), while we instead describe a box position by the interval start.

On a side note, it is known from job sequencing that idle times add a further layer of complexity, e.g., when optimizing non-regular objectives like earliness and tardiness penalties (Garey et al., 1988). Similarly, we presume that the assumption of placing boxes without gaps provides, besides the practical reason, a computational benefit.

4.2 Mixed-integer program

Model (5) is restricted to a mixed-integer program by substituting the placement constraints (5d). There exists a variety of possible formulations in the related domain of job sequencing, cf. Queyranne and Schulz (1994), Keha et al. (2009), Baker and Keller (2010) for reviews. We use disjunctive sequencing constraints to ensure consistency and comparability with the mixed-integer program and the numerical study in Sedding (2020c) for the single-model case $m=1$. For job sequencing, this method is treated, e.g., in Manne (1960), Balas (1985), Queyranne (1993).

Disjunctive sequencing constraints yield a total order on the boxes to decide the position of each. This is established by disjunctive constraints between each pair of jobs.

Let us ease the notation by introducing $(i,j)\prec (h,k)$, a relation between jobs $(i,j),(h,k)\in J$ that holds if and only if $i<h$, and in case of $i=h$, if $j<k$. Moreover, we abbreviate a job’s pair notation by a single letter, i.e., by writing $x=(i,j)$ or $y=(h,k)$ in this section.

Then, (5d) is substituted by

$$\begin{aligned}&\pi _x+w_x\le {}\pi _y\,\,\, \vee \,\,\, \pi _y+w_y\le {}\pi _x,x,y\in J,\,x\prec y, \end{aligned}$$

(6)

$$\begin{aligned}&{V}\le \pi _x\le W-w_x,\,\qquad \qquad \qquad \quad \ x\in J. \end{aligned}$$

(7)

This can be reformulated as a mixed-integer program using the ‘big-M’ method. This relaxes either of the inequalities in (6) by adding W, because $\pi _x+w_x\le W$ for all $x\in J$. Let $u_{x,y}$ denote a binary variable for each job pair $x,y\in J$ with $x\prec y$. Then, while (7) remains, (6) is replaced with

$$\begin{aligned}&\pi _{x}+w_{x}\le {}\pi _{y}+W\big (1-u_{x,y}\big ),&x,y\in J,\,x\prec y, \end{aligned}$$

(8a)

$$\begin{aligned}&\pi _{y}+w_{y}\le {}\pi _{x}+W\,u_{x,y},&x,y\in J,\,x\prec y, \end{aligned}$$

(8b)

$$\begin{aligned}&u_{x,y}\in \{0,1\},&x,y\in J,\,x\prec y. \end{aligned}$$

(8c)

The resulting mixed-integer program (MIP) encompasses constraints (5a)–(5c), (7), (8a)–(8c), with $n\,(n-1)/2$ binary variables.

5 Lower bound

A lower bound on the minimum objective value $\phi ^*$ of a $\mathscr {P}m$ instance is introduced in the following. We consider a Lagrangian relaxation of the mathematical program (5) and show that it is possible to solve it with a quadratic time algorithm.

The lower bound also accepts a partially solved instance for use within a branch and bound search. A solution can be constructed by placing the boxes in the placement interval $[{V},W)$. It starts at ${V}$ with the first box, and places the next box besides. This is repeated until all boxes are placed. An intermediate, partial solution can be subsumed as follows.

Definition 4

A partial solution provides the box position $\pi _{{{i,j}}}$, ${(i,j)\in J_{\textrm{F}}}$, for a set of fixed jobs $J_{\textrm{F}}\subseteq J$ such that these boxes are placed in $[{V},F)$ where ${F}={V+}\sum _{(i,j)\in J_{\textrm{F}}}w_{{{i,j}}}$, i.e., there is ${V}\le \pi _{{{i,j}}}\le F-w_{{{i,j}}}$ for $(i,j)\in J_{\textrm{F}}$. Then, we call $\{\pi _{{{i,j}}}\}_{(i,j)\in J_{\textrm{F}}}$ a partial box placement. It is completed by placing the boxes of the remaining open jobs $J_{\textrm{O}}=J\setminus J_{\textrm{F}}$ between F and $W$.

5.1 Recurrence solving

In (5), we solve the recurrence relation of the completion time variables to a closed form. Let us introduce, for each job $(i,j)\in J$, a continuous variable walking time $\omega _{i,j}$ and deviation $\delta _{i,j}$ of the job’s start time from its ideal start time (i.e., $\pi _{i,j}$). Each completion time variable is replaced by a sum of all assembly and walking times until then. This replaces constraints (5a) to (5c) with

$$\begin{aligned}&C_{i,j}=\sum _{k=1,\dots ,j}\ell _{i,k}+ \omega _{i,k},&(i,j)\in J, \end{aligned}$$

(9a)

$$\begin{aligned}&\omega _{i,j}\ge -a\delta _{i,j},&(i,j)\in J, \end{aligned}$$

(9b)

$$\begin{aligned}&\omega _{i,j}\ge b\delta _{i,j},&(i,j)\in J, \end{aligned}$$

(9c)

$$\begin{aligned}&\delta _{i,j}=-\pi _{i,j}+\sum _{k=1,\dots ,j-1}\ell _{i,k}+ \omega _{i,k},&(i,j)\in J. \end{aligned}$$

(9d)

Walking time is piecewise linear and, because it is minimized, limited from below in constraints (9b) and (9c). It depends on the deviation of a job’s start time to its position, set in constraints (9d).

The walking time variables’ domain can be limited to strengthen the model: the domain is not less than zero, and at most, it corresponds to the walking time in forwards direction to at most the box position $W-1$, and in the backwards direction, to the lowest possible position ${V}$ (and back). Hence,

$$\begin{aligned} {0\le \omega _{i,j}\le \max \!\left\{ -a\left( {C_{i,j-1}-\left( {W-1}\right) }\right) \!,\,b\left( {C_{i,j-1}-V}\right) \right\} } \end{aligned}$$

for $(i,j)\in J$. Substituting the completion times variables with the sum in (9a) gives, for $(i,j)\in J$, the closed formula

$$\begin{aligned} 0\le \omega _{i,j}\le \max \!&\left\{ \!-a\!\left( {1-W+\!\sum _{k=1,\dots ,j-1} \ell _{i,k}+\omega _{i,k}\!}\right) \!\!,\, \right. \nonumber \\ {}&\quad \,\, \left. b\!\left( {{-V+}\sum _{k=1,\dots ,j-1} \ell _{i,k}+\omega _{i,k}\!}\right) \!\!\right\} \! {}. \end{aligned}$$

(9e)

5.2 Lagrangian relaxation

With the model modifications, we perform a Lagrangian relaxation of constraints (9b) and (9c). This introduces the corresponding Lagrangian multipliers $\lambda _{i,j}\ge 0$, $\lambda '_{i,j}\ge 0$ for $(i,j)\in J$.

Then, the Lagrangian problem is

$$\begin{aligned}\phi ^*_{\text {Lagr}}(L)=\min \phi _{\text {Lagr}}\end{aligned}$$

with

$$\begin{aligned} \phi _{\text {Lagr}}&= \sum _{i\in {{I}}} r_i\,C_{i,n_i} +\sum _{(i,j)\in J} \lambda _{i,j} \left( {-a\delta _{i,j} - \omega _{i,j}}\right) \nonumber \\ {}&\quad + \,\lambda '_{i,j}\left( {b\delta _{i,j} - \omega _{i,j}}\right) \end{aligned}$$

(10)

subject to

$$\begin{aligned}L=\left( \left( \lambda _{i,j},\lambda '_{i,j}\right) \right) _{(i,j)\in J}\ge 0,\end{aligned}$$

as well as constraints (9a) to (9e), and constraint (5d).

The set of multipliers L can be optimized using a standard subgradient optimization, see Fisher (2004). Note that the lower bound inequality $\phi ^*_{\text {Lagr}}(L) \le \phi ^*$ holds for any L.

Substituting the completion time variables in (10) according to (9a) yields

$$\begin{aligned} \phi _{\text {Lagr}}&= \sum _{(i,j)\in J} r_i (\ell _{i,j}\!+\!\omega _{i,j}) \!+\! (b\lambda '_{i,j} \!-\! a\lambda _{i,j})\delta _{i,j} \!-\! (\lambda _{i,j}\!+\!\lambda '_{i,j}) \omega _{i,j}\\&= \sum _{(i,j)\in J} \ell _{i,j} \zeta _{i,j} +\underbrace{\sum _{(i,j)\in J}\omega _{i,j} \theta _{i,j}}_{{\varOmega }} + \underbrace{\sum _{(i,j)\in J}\left( {a\lambda _{i,j} - b\lambda '_{i,j}}\right) \pi _{i,j}}_{{\varPi }} \end{aligned}$$

with constants

$$\begin{aligned} \zeta _{i,j}={}&r_i+ \sum _{\begin{array}{c} k=j+1,\dots ,n_i \end{array}} \left( {b\lambda '_{i,k}-a\lambda _{i,k}}\right) \!, \,&(i,j)\in J, \\ \text {and}\quad \theta _{i,j}={}&\zeta _{i,j} - \left( {\lambda _{i,j}+\lambda '_{i,j}}\right) \!, \,&(i,j)\in J. \end{aligned}$$

Observe that the walking time and box placement variables occur only in different constraints.

Property 1

In the Lagrangian problem, the walking time variables $\omega _{i,j}$, $(i,j)\in J$, and box position variables $\pi _{i,j}$, $(i,j)\in J_{\textrm{O}}$ (from Definition 4), are independent.

Thus, it is possible to separately optimize walking times and box positions. This gives us the partial objective

$$\begin{aligned}{\varOmega }=\sum _{(i,j)\in J}\omega _{i,j} \theta _{i,j}\end{aligned}$$

for the walking times, and

$$\begin{aligned}{\varPi }={\sum _{(i,j)\in J}\left( {a\lambda _{i,j} - b\lambda '_{i,j}}\right) \pi _{i,j}}\end{aligned}$$

for the box positions.

5.3 Optimizing box position values

The boxes of the open jobs $J_{\textrm{O}}$ (cf. Definition 4) are to be placed within a box sequence between ${F}$ and $W$. In partial objective ${\varPi }$, each box $(i,j)\in J_{\textrm{O}}$ adds term $\left( {a\lambda _{i,j} - b\lambda '_{i,j}}\right) \pi _{i,j}$. Hence to minimize ${\varPi }$, we get a classic total weighted completion time scheduling problem of the boxes (as jobs), which is solved in polynomial time by sorting them (Smith, 1956).

Lemma 1

Partial objective ${\varPi }$ is minimum if $\pi _{i,j}$, $(i,\!j)\in J_{\textrm{O}}$, are obtained by sequencing $J_{\textrm{O}}$’s boxes nonincreasingly by

$$\begin{aligned}\frac{a\lambda _{i,j} - b\lambda '_{i,j}}{w_{i,j}}.\end{aligned}$$

Thus for $n_{O }=|J_{\textrm{O}}|$, optimal box position values are attained in $\mathscr {O}\!\left( n_{O }\log n_{O }\right)$ time.

5.4 Optimizing walking time values

The walking time variables $\omega _{i,j}$, $(i,j)\in J$, are optimized by minimizing ${\varOmega }$.

For each $(i,j)\in J$, the value range of $\omega _{i,j}$ is limited by constraints (9e). It can be transformed with constants $q_{i,j} = {-V+}\sum _{k=1,\dots ,j-1} \ell _{i,k}$ and $q'_{i,j} = 1-W+q_{i,j}{+V}$ to the range

$$\begin{aligned} 0\le \;\omega _{i,j}\le \max&\!\left\{ \!\underbrace{-a\!\left( {q'_{i,j}+\!\sum _{k=1,\dots ,j-1}\omega _{i,k}\!}\right) \!}_{\alpha _{i,j}},\, \right. \nonumber \\ {}&\quad \left. \underbrace{b \!\left( {q_{i,j}+\!\sum _{k=1,\dots ,j-1}\omega _{i,k}\!}\right) \!\!}_{\beta _{i,j}}\right\} \!.\!\!\!\! \end{aligned}$$

(11)

We observe for any $i\in {{I}}$, and with increasing j from 1 to $n_i$ that term $\alpha _{i,j}$ (as defined in (11)) is nonincreasing and term $\beta _{i,j}$ (as in (11)) is nondecreasing. Hence, we can find some $\kappa _i\in \{0,\dots ,n_i\}$ such that $\alpha _{i,j}>\beta _{i,j}$ for all $j=\kappa _i+1,\dots ,n_i$. Given such a $\kappa _i$, we can replace constraints (11) by

$$\begin{aligned} \begin{aligned}&0\le \omega _{i,j}\le \alpha _{i,j}\,&\text { if }j\le \kappa _i, \\ {}&0\le \omega _{i,j}\le \beta _{i,j}\,&\text { if }j>\kappa _i. \end{aligned} \end{aligned}$$

Hence, depending on the value of $\kappa _i$, either of the two range constraints is active for job $(i,j)\in J$.

We replace the upper bound by an equality with slack variable $0\le y_{i,j}\le 1$. Then,

$$\begin{aligned} \begin{aligned}&0\le \omega _{i,j}= y_{i,j}\alpha _{i,j}\,&\text { if }j\le \kappa _i, \\ {}&0\le \omega _{i,j}= y_{i,j}\beta _{i,j}\,&\text { if }j>\kappa _i. \end{aligned} \end{aligned}$$

(12)

Property 2

Given $\kappa _i$, $i\in {{I}}$, of an optimal solution. If $\alpha _{i,j}<0$ for any job $(i,j)\in J$ with $j\le \kappa _i$, then it is possible to decrease $\kappa _i$ without changing the objective ${\varOmega }$ such that $\alpha _{i,j}\ge 0$ for each job $(i,j)\in J$ with $j\le \kappa _i$.

Proof

Given the described case, then $\kappa _i\ge 1$ and $\alpha _{i,\kappa _i}<0$ because $\alpha _{i,j}$ is decreasing with j. Hence, $y_{i,\kappa _i}=\omega _{i,\kappa _i}=0$ to fulfill constraints (12).

Let us decrease $\kappa _i$ by one. Then, it is possible to leave $y_{i,\kappa _i+1}=\omega _{i,\kappa _i+1}=0$ in the solution. Thus, ${\varOmega }$ remains unchanged. We repeat this step until $\alpha _{i,\kappa _i}\ge 0$. $\square$

Corollary 1

An optimum solution exists where $\alpha _{i,j}\ge 0$ holds for each $(i,j)\in J$ with $j\le \kappa _i$.

Lemma 2

For each $i\in {{I}}$, let $\kappa _i^*$ be the minimum value for $\kappa _i$ that permits an optimum solution. Then, there exists such a solution where for each job $(i,j)\in J$ there is

$$\begin{aligned}&y_{i,j}={\left\{ \begin{array}{ll} 0,\!\!\!\!\!&{}\quad \textrm{if }\,\,\, \displaystyle \theta _{i,j}+\!\!\!\sum _{k=j+1,\dots ,n_i}\!\! y_{i,k}c_{i,k}\theta _{i,k}>0, \\ 1,\!\!\!\!\!&{} \quad \textrm{else} \end{array}\right. } \nonumber \\&\textrm{with}\quad \! c_{i,k}={\left\{ \begin{array}{ll} -a,\!\!\!\!\!&{}\quad \textrm{if }\,\,\, k\le \kappa _i^*,\\ b,\!\!\!\!\!&{} \quad \textrm{else}. \end{array}\right. } \end{aligned}$$

(13)

Proof

For each model $i\in {{I}}$, we show this by induction for $j=n_i,\dots ,1$. Then,

$$\begin{aligned}&\omega _{i,j}=y_{i,j}\cdot \underbrace{c_{i,j}\left( {d_{i,j}+\sum _{k=1,\dots ,j-1}\omega _{i,k}}\right) }_{\ge 0} \\&\textrm{with}\quad \! c_{i,j}={\left\{ \begin{array}{ll}-a,&{}\!\!\text {if }j\le \kappa _i^*,\\ b,&{}\!\!\text {if }j>\kappa _i^*,\end{array}\right. }\quad d_{i,j}={\left\{ \begin{array}{ll}q'_{i,j},&{}\!\!\text {if }j\le \kappa _i^*,\\ q_{i,j},&{}\!\!\text {if }j>\kappa _i^*.\end{array}\right. } \end{aligned}$$

By choice of $\kappa _i$ and use of Property 2, the slack variable $y_{i,j}$ multiplies a nonnegative value. Hence, $\omega _{i,j}$ is nonnegative. Moreover, $\omega _{i,j}$ influences $\omega _{i,k}$ for each $k=j+1,\dots ,n_i$ unless $y_{i,k}=0$. Thus, $\omega _{i,j}$ contributes to the objective ${\varOmega }$ not only with factor $\theta _{i,j}$, but moreover via $\omega _{i,k}$ with factor $y_{i,k}c_{i,k}\theta _{i,k}$. The total contribution of $\omega _{i,j}$ to ${\varOmega }$ is thus with factor $\theta _{i,j}+\sum _{k=j+1,\dots ,n_i}y_{i,k}c_{i,k}\theta _{i,k}$. If this factor is positive, then the lowest slack value $y_{i,j}=0$ minimizes ${\varOmega }$. If it is negative, then the highest slack value $y_{i,j}=1$ is optimal. If the factor is zero, any value for $y_{i,j}$ is optimal. $\square$

Let $\overline{\kappa _i}\in \{0,\dots ,n_i\}$ for each $i\in {{I}}$ be the maximum $\kappa _i$ for which $-aq'_{i,\overline{\kappa _i}}\ge 0$ holds.

Property 3

An optimum solution exists where $\kappa _i\le \overline{\kappa _i}$ for each $i\in {{I}}$.

Proof

Assume we are given an optimum solution. Naturally, all walking time variables $\omega _{i,j}$, $(i,j)\in J$, are nonnegative. For each $i\in {{I}}$, both $\sum _{k=1,\dots ,j-1}\omega _{i,k}$ and $q'_{i,j}$ are nondecreasing with respect to j, while $-aq'_{i,j}$ is nonincreasing. Hence if $-aq'_{i,{\kappa _i'}}<0$ for any $(i,{\kappa _i'})\in J$ with ${\kappa _i'}\le \kappa _i$, then $\alpha _{i,j}<0$ for $j={\kappa _i'},\dots ,{\kappa _i}$. However it is, according to Property 2, possible to set $\kappa _i<{\kappa _i'}$ such that the solution is optimum and $-aq'_{i,j}<0$ as well as $\alpha _{i,j}\ge 0$ hold for any $(i,j)\in J$ with $j\le \kappa _i$. $\square$

The presented results allow us to describe an algorithm to minimize the walking time variables.

Lemma 3

In an outer loop, set $\kappa _i=0,\dots ,\overline{\kappa _i}$ for $i\in {{I}}$. Given $\kappa _i$, Lemma 2’s recurrence (13) is used to set $y_{i,j}$ for each $(i,j)\in J$, which also sets $\omega _{i,j}$ and objective value ${\varOmega }$. Then, the smallest obtained objective value is optimal.

This algorithm takes quadratic time in terms of $n_i$. Over all m models, the worst case runtime is $\mathscr {O}\!\left( \sum _{i\in {{I}}}n_i^2\right) \le \mathscr {O}\!\left( n^2\right)$.

Combining the algorithm in Lemma 3 with the box sequencing procedure in Lemma 1, we are able to find a solution for the whole Lagrangian problem.

Theorem 2

An optimum solution to $\phi ^*_{\text {Lagr}}(L)$ is computed in $\mathscr {O}\!\left( n^2\right)$ time.

6 Branch and bound methods

For solving $\mathscr {P}m$ instances, we describe a branch and bound search (B &B) and a heuristic version in the following. The upper bound is initialized with basic heuristics. Bounding is performed with the above Lagrangian based lower bound and an additional combinatorial lower bound. The heuristic version introduces a heuristic dominance rule and limits the number of visited descending nodes.

6.1 Basic heuristics

To construct a good upper bound for the branch and bound search we use a constructive heuristic, a local search, and a simulated annealing metaheuristic.

Weighted nearest identity (WNID) As a construction heuristic, an intuitive way to place the boxes is in the same sequence as the jobs, which is already reported in Ford and Crowther (1922, p. 80). For a single model, this strategy is described in Sedding (2020a); it can be called identity sequence placement heuristic. Extending this to multiple models, our strategy is to intersperse the identity sequences of all models ${{I}}$ such that

$$\begin{aligned} \pi _{i,j}<\pi _{i,j'} \text { for all }(i,j),(i,j')\in J\text { with }j<j'. \end{aligned}$$

(14)

Our approach is greedy, it places all boxes iteratively along the line, starting at position 0. In an iteration step, we select the next box to place. To fulfill (14), there are at most m possible boxes to choose from. We rank the boxes by the resulting walking time weighted by the reciprocal of the corresponding model’s production share. Accordingly, we call this method weighted nearest identity (WNID) sequence placement heuristic, see Algorithm 6.1.

Hill Climbing (HC) A local search typically improves a constructive heuristic’s solution. For this, we use the transpose neighborhood that comprises all possible swaps of neighboring boxes. Hence, the number of neighbors is $|J|-1$ for |J| boxes. The local search procedure is initialized with the WNID solution. Of the current solution’s neighborhood, this procedure searches for the neighbor with highest improvement of the objective value. If such a neighbor exists, it is selected as the current solution. This is known as hill climbing (HC) search.

Simulated Annealing (SA) is a metaheuristic that escapes local minima (Kirkpatrick et al., 1983; Černý, 1985). In comparison, the local search procedure stops at some local minimum. It might occasionally correspond to a global minimum, but it is usually sensible to apply SA: Given certain conditions, SA converges to a global optimum (Hajek, 1988). To cross the solution space, SA permits worse solutions with a decreasing probability. In our case, we use the basic and widespread procedure of Press et al. (1992).

6.2 Branch and bound algorithm

An exact solution for $\mathscr {P}m$ can be computed by a branch and bound search (B &B). We use the above described heuristics’ objective value as an initial upper bound.

Branching Our B &B performs a depth first search, discarding nodes that do not lead to an optimal solution. Any node corresponds to a partial box placement $\{\pi _j\}_{j\in J_{\textrm{F}}}$ of some fixed jobs $J_{\textrm{F}}\subseteq J$, as of Definition 4. The root node initializes fixed job set $J_{\textrm{F}}=\emptyset$ and open job set $J_{\textrm{O}}=J$. A leaf node is reached if $J_{\textrm{F}}=J$ and $J_{\textrm{O}}=\emptyset$. A descending node is created by fixing an open job $j\in J_{\textrm{O}}=J{\setminus J_{\textrm{F}}}$, placing it at $\pi _j={V+}\sum _{k\in J_{\textrm{F}}}w_{k}$. A node is discarded if its lower bound is not less than the upper bound. The latter corresponds to the currently best known box placement’s objective value, which is initialized with the described HC or SA heuristics and updated whenever a better solution is reached at a leaf node.

Combinatorial lower bound When our B &B search visits a new node, it is first evaluated by a quickly obtained lower bound. This allows to potentially discard it before computing the more elaborate Lagrangian lower bound. Let us first consider a trivial lower bound: It places each box exactly at its ideal time: the corresponding job’s start time. In a partial solution, this bound can be increased by using the fixed $\pi _{i,j}$ values of all fixed jobs $(i,j)\in J_{\textrm{F}}$. Although the position of the open boxes is not yet known, we know that they will be placed in the unoccupied interval $[F,W)$. Hence, if the start time $t_{i,j}$ of an open job $(i,j)\in J_{\textrm{O}}$ lies within $[F,W-w_{i,j}]$, then we still place the box exactly at this start time. Otherwise, the box can be placed at F or at $W-w_{i,j}$, whichever is closer to $t_{i,j}$. This increases the lower bound further. Concluding, we calculate the start time for each job in a model iteratively. If we encounter an open job $(i,j)\in J_{\textrm{O}}$, then its box position is temporarily set to

$$\begin{aligned}\pi _{i,j}=\max \{F,\,\min \{t,\,W-w_{i,j}\}\}.\end{aligned}$$

Lagrangian lower bound Secondly, our B &B evaluates a partial solution with the Lagrangian based lower bound of Sect. 5. While the lower bound is calculated in quadratic time for a given set of Lagrangian multipliers L, the multipliers are iteratively improved using a subgradient optimization based on Fisher (2004), Held et al. (1974) as follows.

• In our case, the iteration step’s updated set of multipliers $\hat{L}$ for step size $s>0$ is

  $$\begin{aligned}&\hat{\lambda }_{i,j}=\max \!\left\{ \lambda _{i,j}+s\left( {-a\delta _{i,j}-\omega _{i,j}}\right) \!,\,0\right\} \!, \\&\hat{\lambda }'_{i,j}=\max \!\left\{ \lambda '_{i,j}+s\left( {b\delta _{i,j}-\omega _{i,j}}\right) \!,\,0\right\} \!, \end{aligned}$$

  for each $i,j\in J$. For the step size, we employ the common

  $$\begin{aligned}&s= \frac{\displaystyle v\cdot \left( {\textit{UB}-\phi ^*_{\text {Lagr}}\!\left( {L}\right) }\right) }{\displaystyle \sum _{(i,j)\in J} \left( {-a\delta _{i,j}-\omega _{i,j}}\right) ^2+\left( {b\delta _{i,j}-\omega _{i,j}}\right) ^2} \end{aligned}$$

  where $\textit{UB}\ge \phi ^*$ is an upper bound value for $\phi ^*$, and $v\in (0,2]$ is a step size factor. Our initial value for the latter is $v=1$. Then, we divide v by two after 10 iterations of no improvement on the Lagrangian objective $\phi ^*_{\text {Lagr}}$.

• We reduce the number of iterations by reusing multiplier values during the B &B search like in Fisher (2004). This avoids optimizing L from scratch for each partial solution. Exactly one set of L values is memorized. Hence, the last obtained multiplier values L provide a warm-start in the next bound calculation.

• Our B &B performs a higher number of iterations earlier in the search tree. A better bound has a greater utility there. We iterate at most

  $$\begin{aligned}\max \!\left\{ 1,\left\lfloor 4\cdot \sqrt{|J_{\textrm{O}}|} \right\rfloor \right\} \!\end{aligned}$$

  times; except at the root node (with $|J_{\textrm{O}}|=n$), where we allow for up to $10n$ iterations, to initialize L from scratch.

Traversing order The Lagrangian lower bound calculation is additionally used in our B &B search to set a traversing order. This determines, in each node, the sequence for visiting the descending nodes. More promising descending nodes should be visited first. Each descending node places a different box at F. Lemma 1’s artificial position value provides a hint for good posititions of the open boxes $J_{\textrm{O}}$. This motivates our use of $\pi _{i,j}$, $(i,j)\in J_{\textrm{O}}$, as a traversing order value. We use the nondecreasing order of $\pi _{i,j}$ for ranking the open boxes and visit the descending nodes accordingly.

6.3 Truncated branch and bound algorithm

Our truncated branch and bound (TrB &B) leaves out some nodes in the B &B search, which yields a heuristic. Moreover, a node is evaluated with a heuristic dominance rule that compares a currently visited node to previously visited nodes.

Adapting the approach in Sedding (2020a), we limit the number of descending nodes to the maximum branching factor

$$\begin{aligned}\textrm{BF}_{\max }= \min \!\left\{ |J_{\textrm{O}}|,\max \!\left\{ \lceil \psi \rceil ,\left\lfloor \frac{|J_{\textrm{O}}|}{\sigma }\right\rfloor \right\} \right\} \end{aligned}$$

for constant $\psi >0$ and $\sigma >0$. While the number of descending nodes is restricted by $\sigma$ near the root of the search tree, it is restricted by $\psi$ when being deeper in the search tree.

With our heuristic dominance rule, the current node is discarded if the partial solution is (probably) dominated by any previous partial solution. For this, we adapt the approach in Sedding (2020a). Then, the heuristic dominance rule works as follows. Based on a partial box placement, it creates two complete artificial placements:

• (a) All open jobs $j\in J_{\textrm{O}}$ are placed such that $\pi _j=F$,

• (b) all open jobs $j\in J_{\textrm{O}}$ are placed such that $\pi _j=W-w_j$.

This yields two heuristic objective values $\phi _{\text {(a)}}$ and $\phi _{\text {(b)}}$. Then, a partial solution is heuristically dominated (and can be discarded) unless at least one of the two values is less than a previously found value, respectively, for the same set $J_{\textrm{F}}$ of fixed jobs.

7 Numerical results

A quantitative evaluation of the described solution methods is performed in the following numerical study. We create artificial instances and perform a full factorial evaluation with the exact approaches and, secondly, with the heuristics.

Table 2 Exact algorithms’ runtime in seconds, grouped by n and m

7.1 Instance generation

We generate random instances in a variety of parameter settings. Comparability to Sedding (2020c) is ensured by using the same variants for assembly times, box widths, and slopes.

• (1) The number of product models is set to $m\in \{2,4,8\}$.

• (2) The production shares of all models ${{I}}$ need to be positive and sum up to one. We generate these shares randomly analogous to Boysen et al. (2008, 2009a, b) as follows. Let $PC$ denote the total production cycle count. Initialize each model’s demand to one unit. Then, select a model with uniform probability, increase the model’s demand by one unit. Repeat this step until the total demand equals $PC$. This method is equivalent to a uniform sample (with replacement) of $( PC -m)$ values from ${{I}}$. Then, the frequency of each model plus one corresponds to the model’s demand. Finally, dividing the demand of each model by $PC$ gives its production share. We let $PC =1000m$ to avoid quantization. Note that in the limit for $PC \rightarrow \infty$, the expected value for each share is 1/m.

• (3) The number of jobs is set to $n\in \{8,12,\dots ,28\}$. A job’s model is selected with uniform probability. Hence, each model’s expected number of jobs is n/m.

• (4) The jobs’ assembly times are generated in four variants (Jaehn and Sedding, 2016; Sedding, 2020a):

  • (L1) all equal (unitary value 1),

  • (L2) all distinct (a random permutation of $(1,2,\dots ,n)$),

  • (L3) uniform random variates from $\{1,2,\dots , 10\}$,

  • (L4) geometric random variates with mean $\lambda =2$.

• (5) The box widths of the jobs are drawn in either of four variants, and normalized by scaling and rounding to obtain a uniform total width $10\cdot n$ (Sedding, 2020a):

  • (W1) all equal box widths (unitary value 1),

  • (W2) all distinct (a random permutation of $(1,2,\dots ,n)$),

  • (W3) uniform random variates from $\{2^0,\dots ,2^3\}\cup \{3\cdot 2^0,\dots ,3\cdot 2^2\}$, which represent the ISO1-pallet divisions of Euro stacking boxes,

  • (W4) rounded up gamma variates with shape 1.25 and unit scale.

• (6) The slopes a, b are determined by the walking velocity v, which is expressed in multiples of the conveyor velocity, and the walking strategy. Like in Sedding (2020a), we vary the walking velocity $v\in \{2,4,8,16\}$ while the conveyor velocity is kept unitary. Then, the slope values follow from the chosen strategy:

  • (S1) $a=2/(v+1)$ and $b=2/(v-1)$ (walking either besides or atop the moving conveyor),

  • (S2) $a=(2v+1)/(v+1)^2$ and $b=(2v+1)/v^2$ (walking both besides and atop, always choosing the faster option).

• (7) Finally, processing times are harmonized with the box widths with respect to the walking velocity. In the test, we consider the case of a zero box placement interval start $V=0$ and that its end $W$ is approximately equal or larger than the mixed-model weighted average makespan, i.e., the objective value $\phi$. Otherwise, the problem is likely easier, as many start times will be after the box position (if $V\ll 0$ or $W\ll \phi ^*$, cf. Sedding (2020b)) or before (if $V\gg 0$ or $W\gg \phi ^*$). Note it is not necessary to compute the optimum $\phi ^*$ during harmonization because strict equality of W and $\phi ^*$ is not required. Thus, we can use the rather intuitive solution of the WNID heuristic. Then, the deviation $|\phi -W|$ is minimized by linearly scaling all assembly times uniformly with the same factor, using the univariate optimization method in Brent (1971). This procedure is repeated at most 10 times because the WNID’s solution can change due to the scaling.

The parameter setting variants yield a product of $3\cdot 6\cdot 4\cdot 4\cdot 4\cdot 2= 2304$ settings. For each setting, ten instances are generated, which yields 23040 instances total.

7.2 Test setup

Performance comparablity to the single-model study in Sedding (2020c) is ensured by, besides using the same mixed-integer programming approach, using equal hardware and software.

The algorithms are all implemented in the C++11 programming language. C++ STL containers are used for the data structures. Computations remain single-threaded. The code is compiled by GCC to x86-64 binaries. Mixed-integer programs are given to the Gurobi 7.5 C++ API and solved with it. For a fair comparison, the Gurobi computation is limited to one thread. For Simulated Annealing, the reference implementation of Press et al. (1992, pp. 448–451) with default parameters is used. The TrB &B heuristic is paramerized with $\psi =5$, $\sigma =7$.

All programs are ran on Ubuntu Linux. Each instance gets a dedicated CPU core of an Intel Xeon E5-2680.

A computation is terminated after a 1 h time limit. In this case, although the final runtime is unknown, it is lower bounded to at least 1 h. Thus, runtime percentiles (and median) can still be calculated if the respective share of instances finished within the time limit.

7.3 Exact algorithms

Tested exact approaches are the mixed-integer program (MIP, Sect. 4.2) and the branch and bound (B &B, Sect. 6.2). Resulting runtimes are aggregated in Table 2, in groupings by n and m.

Let us first analyze the MIP and the B &B together:

• For both, we can observe an exponential growth of runtime with increasing n, and some increase along with m. This happens irrespective of the number of models m. With the strong NP-hardness proof for $\mathscr {P}m$ in Theorem 1, this behavior is not unexpected. Both approaches exhibit high quartile deviations. The decrease in the quartile deviation in groups of higher n might be explained by the smaller number of solved instances in such a group. The correlation between B &B and MIP runtimes is 0.25, which is weakly positive. Hence both approaches experience some similar difficulty with the instances.

• The B &B is consistently faster than the MIP. The B &B is faster in 18718 (81.24%) of all 20200 solved instances. While B &B is able to solve a clear majority of instances until $n=24$, the MIP can solve the majority only until $n=16$. On an aggregate level, the median speed advantage of the B &B is 66-fold.

Table 3 Median B &B runtime in seconds for $n=24$ and $m=4$ grouped in walking velocity and walking strategy pairs

Table 4 Median B &B runtime in seconds for $n=24$ and $m=4$ grouped in pairs of assembly time and box width settings

Further analysis of B &B performance shows peak difficulty for walking velocity $v=4$ and for box width setting W4. It has lowest difficulty for W1, but remains inconclusive for the different assembly time settings:

• Table 3 lists median runtimes grouped by each (v, S) pair. Both walking strategies yield similar runtimes. The runtime is smallest for $v=16$, rises highest at $v=4$, and falls again for $v=2$. Note that this picture is different from the single-model case in Sedding (2020a): their B &B for $m=1$ had the highest runtime for $v=16$, and the lowest for $v=2$.

• Table 4 groups each (L, W) pair. Compared to the other box width settings W4 is highest with L1, L2, and L4; W3 is highest with L3 (with W4 close up). Hence, W4 is likely the most difficult case.

• Of the assembly time settings in Table 4, L2 runtime is highest for W1 and W2; while L3 is highest with W3 and W4. L4 is lowest for W1, L2 for W2, L3 for W3 and W4. With this diverse ranking it remains unclear which case has the highest or lowest performance.

Concluding, the B &B is clearly the best performing exact approach for solving the $\mathscr {P}m$ with a median 66-fold speed advantage over the MIP. The B &B exhibits relatively uniform performance in all box width variants, is slower with the more realistic assembly times W3 and W4, and is faster for instances of the more realistic velocities $v=8$ and $v=16$.

Table 5 Median (50% percentile), quartile (25%, 75% percentile), and 95% percentile minimum walking time percentage of total weighted work time for $n=24$ and $m=4$ instances with a known optimum, in walking velocity and walking strategy pairs

7.4 Proportion of walking time to total work time

The exact solutions allow us to analyze the share of walking time compared to the total work time in our instances. Scholl et al. (2013) reports that material fetching amounts to about 10–15% of total work time at an assembly line of a large German automobile manufacturer. Our most realistic worker velocity cases are $v=8$ and $v=16$ (Klampfl et al. (2006) document a similar assumption of $v=13.6$). Table 5 shows those cases attain 5.7–11% median walking time for walking strategy S1. It is slightly lower 5.2–11% in strategy S2. Note this is less than the 9–15% mean walking time in Sedding (2020a) for a single product model $m=1$.

7.5 Comparison with constant walking time estimates

In most of the literature, it is assumed that the time to gather material at the line side of a moving conveyor can be estimated by constant estimates. This includes methods time measurements (MTM) and most scientific literature. A compromise might be to represent walking time as costs (Limère et al., 2015; Schmid et al., 2021; Müllerklein et al., 2022). However, this can be imprecise if the space at the line side is scarce and high walking times occur.

In Table 5, a high variability of the minimum walking time can be observed. Thus, using the mean walking time as a constant estimate poses the risk of misguided planning. However, a conservative estimate can be far off. For example for $v=8$ in S1, the 95% percentile is 16.8% walking time. Taking this as a conservative estimate yields an overestimation of walking time in a quarter of the instances (Q25) by at least 75%, and in half of the instances (the median is 11.1%) by at least 51%. Similarly for $v=16$ in S1: here, the 95% percentile is 9.6%; hence the overestimation is over 108% compared to the lower quartile and 68% compared to the median. Therefore, we conclude that constant estimates of walking time are either misguiding or overly conservative because of the observed variability in the instances’ minimum walking times.

An accurate depiction of walking times would henceforth improve planning accuracy, for example, in assembly line balancing for an even distribution of work to workers (Boysen et al., 2022), or in product sequencing, which determines the order of product models at the mixed-model line (Becker and Scholl, 2006; Boysen et al., 2009c, 2012). Both problems would ideally be considered simultaneously (Boysen et al., 2022). An accurate, time-dependent walking time estimate aids to reduce overload situations and minimize idle time. As elucidated in Sect. 2.4, the worker may start a cycle at a varying position due to floating, depending on the previous cycle’s completion time. Hence, walking times are highly dynamic. To minimize them, it may be useful to continuously reorder the worker’s operations where possible, e.g., with methods described in Jaehn and Sedding (2016), Sedding (2020b).

7.6 Heuristics

Tested basic heuristic approaches (cf. Sect. 6.1) are: the weighted nearest identity sequence (WNID) placement heuristic that provides the start for a best improvement local search (HC) and a simulated annealing (SA) metaheuristic. Also, we tested the truncated branch and bound (TrB &B, Sect. 6.3), initialized with the WNID upper bound (denoted TrB &B_UB HC) and also the SA upper bound (TrB &B_UB SA). For a comparison, we let the MIP terminate on time limits of 10, 60, and 3600 s.

We measure solution quality by the percentage increase of walking time, which divides additional walking time by minimum walking time. This is captured by the percentage walking time error

$$\begin{aligned} \textrm{PE}(\phi ,\phi ')=\frac{\phi -\phi '}{\phi '-\sum _{(i,j)\in J}r_i\ell _{{i,j}}}\cdot 100\% \end{aligned}$$

(15)

for two objective values.

We perform a subgroup analysis of the heuristics on the instance subset that is exactly solved by B &B within the 1 h time limit. Let MPE denote mean $\textrm{PE}(\phi ,\phi ^*)$ of the attained objective $\phi$ and an instance’s optimum $\phi ^*$. Please refer to Table 6 to observe the share of optimally solved instances and MPE values grouped by n and m values, and (n, m) pairs. The number of finished instances increases over time, which is plotted for several approaches in Fig. 3 for $n=24$, $m=4$. The resulting PE is shown as box plots in Fig. 4.

Table 6 Heuristic solutions on instances solved optimally by B &B, grouped by n, m, or both

Table 7 Heuristic solutions on instances compared to the heuristic MIP (<1 h), grouped by n, m, or both

Fig. 3

Percentage of finished instances with $n=24$ and $m=4$ in a line plot of algorithm runtime in seconds

Fig. 4

Box plots of percentage walking time error $\textrm{PE}(\phi ,\phi ^*)$ of a heuristic’s $\phi$ for $n=24$ and $m=4$ instances solved by B &B with optimum $\phi ^*$

The WNID yields a weak performance, still it is able to achieve an optimum solution for some small instances. The HC much improves on this, attains an optimum very frequently in 39% of the cases. Its computation time is barely measurable. The SA improves on this by optimally solving 68% of the instances. For $n=24$ and $m=4$, its median runtime is measurable with 0.026 s.

The TrB &B greatly improves on both the HC and the SA’s upper bound. In the former case, it finds the optimum solution for 46% of the instances. It reduces the MPE from 121% to 5.05%. For $n=24$ and $m=4$, the median runtime rises to 0.146 s.

With the SA upper bound, an optimum is found for 72% of the instances, further reducing the MPE to 1.77%. For $n=24$ and $m=4$, the median runtime is 0.140 s, which includes the SA runtime.

Except for small instances, the time-limited MIP at 10 or 60 s has worse PE and MPE values even though the runtime is much higher than of the other heuristics: The MIP’s median runtime corresponds to the time limit because available compute time is mostly used up completely. With a long runtime, only a small PE remains.

An assessment of the full set of instances would require knowledge of an optimal solution of all instances. For 2463 instances, an optimum could not be attained within the 1 h time limit with neither exact algorithm (MIP, B &B). To study the set of all instances in a consistent way, we use the objective value $\phi ^{\text {H}}$ of the long running MIP (<1 h) as a reference, because it returns the smallest overall MPE in the preceding subgroup analysis. The heuristics’ objective value $\phi$ can then be assessed with the percentage of instances that yield $\phi \le \phi ^{\text {H}}$, and with the mean of the percentage walking time error $\textrm{PE}(\phi ,\phi ^{\text {H}})$ denoted by HMPE. We observe similar results for this test as in the subgroup analysis. In particular, executing the MIP for 60 s gives yields relatively high HMPE for large instances. Compared to the MIP runtime of 1 h, the HMPE of the TrB &B is noticeably less for high n and m values, while maintaining a much shorter runtime. See also Table 7.

Concluding, the TrB &B_UB SA provides the best heuristic solutions, makes good use of SA’s initial upper bound, and retains a low runtime. This suggests that the research and implementation effort for the TrB &B is worthwhile, especially in combination with SA. Both the TrB &B and the MIP can be parameterized regarding solution quality, although the latter admits greater walking time in the same runtime.

8 Conclusion

In this paper, we consider the mixed-model placement of boxes to minimize worker walking at the moving assembly line. The time-dependent walking time model allows for a much higher precision according to our numerical experiment: in most instances, constant walking time estimates that include 95% of cases are at least 51% higher than with our time-dependent model (see Sect. 7.5).

We prove that this optimization problem is NP-hard in the strong sense for any number of product models. We observe that moderately sized instances with up to 16 jobs can be solved with a mixed-integer program that employs disjunctive sequencing constraints to avoid box overlapping. For larger instances, we construct branch and bound-based algorithms. A Lagrangian relaxation leads to a lower bound that is solved in quadratic time. This bound is employed in a branch and bound search, for which a truncated search tree yields a heuristic. Also, we describe an intuitive heuristic that is complemented with a local search and simulated annealing to find an initial upper bound.

The numerical results indicate that our exact branch and bound based algorithms provide superior runtime and quality compared to solving a mixed-integer program. The runtime of the best heuristic typically remains below 1 s, contributing to an interactive planning experience that is more accurate and permits less safety buffer time.