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

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 difﬁcult 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% conﬁdence 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.


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 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 manufac-  (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, 2017Sedding, , 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 (2017Sedding ( , 2020aSedding ( , 2020c 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 timedependent walking times to gather material are minimized: -First, we display all assumptions and introduce an optimization model in Sect. 2. In comparison to Sedding (2020aSedding ( , 2020c, 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 singlemodel case in Sedding (2020aSedding ( , 2020c. -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.
formal definition of all parameters in a Pm instance, its prerequisites are introduced. The definition of Pm follows. Note that Pm builds upon the optimization problem in Sedding (2020aSedding ( , 2020c 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 Pm is gathered in Table 1. as walking time slope (box behind), − I as set of product models, − m = |I | as number of product models, as production rate of model i ∈ I s.t. i∈I r i = 1, as set of all jobs, j ∈ {1, . . . , n i }} − n i ∈ N as number of jobs of model i ∈ I , which remain in the fixed sequence (i, 1), (i, 2), . . . , (i, n i ), − n = |J | = i∈I n i as total number of jobs, Before the optimization problem Pm 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.

Box placement
Necessary material of an assembly operation is provided in a dedicated box (i, j) ∈ 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 ∈ 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 Then, a box placement A3 Multiple product variants We consider a model-mix of several product variants, each with a different production rate and its own set of assembly operations and containers. While a production sequence of the product variants is not yet determined, it is assumed that the sequence fulfills the production rates on average A4 Fixed operation sequence A worker has an immutable sequence of operations albeit it may be changed in practice A5 Single cycle A production cycle begins when the workpiece enters the worker's zone, and it ends when all operations have been finished on this workpiece. If that takes shorter or longer than the average cycle time, the worker can start the next cycle early or late. This may cause a different starting point, highly depending on the production sequence. However, the sequence is not available during material placement. Thus, we assume an average cycle time and disregard floating A6 Single work point The work point is at one location at the workpiece. We assume that other work points are assigned to different workers on assembly line balancing (Becker and Scholl, 2009) A7 One box per operation One single box aggregates all material containers of an operation, e.g., as a stack or a shelf A8 One-dimensional box placement Boxes are placed along a line parallel to the assembly line. Hence, the box width along this dimension is the decisive size to measure a box's occupied interval at the line side A9 No space between boxes There are no gaps between adjacent boxes in our model. This assumption is realistic as space at the line side is typically scarce (Sternatz, 2015). Note that more elaborate logistics operations like material sequencing or kitting may reduce required space at the line side (Boysen et al., 2015;Schmid and Limère, 2019) A10 Stationary box positions Box positions are fixed during production and optimized beforehand A11 Uniform walking velocity Walking velocity is constant and the same for all operations A12 One-dimensional walking Only walking in parallel to the assembly line is counted: Boxes are located in close tangency the conveyor to reduce the orthogonal margin to a gripping distance A13 One walk per operation The worker fetches necessary parts all and only before the start of an assembly operation from a single box as in A7, An operation can aggregate several sub-operations, for all of which material is gathered in one single run A14 No pick time Pick times are neglected. Finnsgård et al. (2011) reports that pick times encompass just 6% of nonproductive time A15 Picking at upstream side We simplistically assume picking at the upstream (left) edge, albeit several pick points may exist for a box A16 Placement area offset The box placement interval start can be given with a shift off the station start along the line 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 π 3,1 ) and assembly time 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.

states, for each box
Note that attaining a box placement is equivalent to finding a sequence for the boxes along the line side.

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)∈ J , the worker leaves the workpiece at a certain time t, arrives at the box's position π 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 where π i, j is the box position encoded in time units, and 0 ≤ a ≤ 1 and b ≥ 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 = π i, j ), then the walking time is minimum, which occurs if the work-piece just passes by the box position. Otherwise, the walking time increases linearly. If the workpiece has not yet passed the box position (case t < π i, j ), then the walking time cor- The two slopes relate the workpiece's velocity v conveyor to the worker's walking velocity v worker > v 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 worker /v 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 < π i, j , else 4.1%, see Sedding (2020a).

Makespan calculation
Before an assembly operation can be processed, a walking time occurs to its distinct box (i, j) ∈ 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 i, j ≥ 0 that is a constant nonnegative rational number. Together, they form the job's processing time p i, j (t) = f i, j (t) + i, j . Starting at time t, the job completes at C i, j (t) = t + p i, j (t). Substituting all components, Each product model i ∈ I requires processing of a certain number n i of jobs in a fixed sequence denoted by where n i is the number of jobs in model i. Then, the last completion time C max i in a model i is the composition of job completion times (3), starting the first job (i, 1) at time 0: Note that a start of the first job at t min = 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(Sedding, , 2020c. This belongs to the field of time-dependent scheduling (Gawiejnowicz, 2020b, 2020a). 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).

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 max i weighted by production share r i . This objective is minimized in the studied optimization problem.
Definition 2 (Problem Pm) Given a Pm instance (see Definition 1), minimize the weighted average makespan φ = i∈I r i C max i by determining a box placement {π i, j } (i, j)∈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 ∈ I , makespan 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) ∈ J as a function of job start time t is given by 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 φ 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.

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 Pm, we prove that it is NPhard in the strong sense for an arbitrary number of product models m ≥ 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).
Theorem 1 Pm is NP-hard in the strong sense for arbitrary m ≥ 1.
Proof For m = 2, we perform a reduction from 3-Partition as of Definition 3. This requires a decision version of Pm that specifies a threshold value Φ and asks if a solution exists with objective value φ ≤ Φ.
A corresponding instance is constructed as follows. First, we freely choose any allowed nonzero slope 0 < a ≤ 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, . . . , 3z, and choose an arbitrary 1, j value. The second model has n 2 =z +1 jobs. For each j = 1, . . . , z +1, we set w 2, j = 1 and 2, j = B + 1. Finally, we set threshold value This instance's objective value is, given the unilateral production shares, φ = C 2,z+1 . As This requires in the second product model for each j = 1, . . . , z + 1 that p 2, j = 2, j . This is the case if and only if the corresponding box is precisely positioned at (B + 1) · 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, Pm 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 ∈ 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 ∈ 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 φ ≤ Φ, 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 Pm exists for m ≥ 1.

Mathematical programming
We describe Pm solutions using mathematical programming, adapting the special case m = 1 in Sedding (2020aSedding ( , 2020c to m ≥ 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.

Mathematical program
We introduce continuous variables π i, j as box position and C i, j as completion time of job (i, j) ∈ J . Then, a mathematical program for Pm can be stated as: 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.

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 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 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 π x + w x ≤ W for all x ∈ J . Let u x,y denote a binary variable for each job pair x, y ∈ J with x ≺ y. Then, while (7) remains, (6) is replaced with The resulting mixed-integer program (MIP) encompasses constraints (5a)-(5c), (7), (8a)-(8c), with n (n − 1)/2 binary variables.

Lower bound
A lower bound on the minimum objective value φ * of a Pm 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.

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) ∈ J , a continuous variable walking time ω i, j and deviation δ i, j of the job's start time from its ideal start time (i.e., π 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 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, Substituting the completion times variables with the sum in (9a) gives, for (i, j) ∈ J , the closed formula

Lagrangian relaxation
With the model modifications, we perform a Lagrangian relaxation of constraints (9b) and (9c). This introduces the corresponding Lagrangian multipliers λ i, j ≥ 0, λ i, j ≥ 0 for (i, j) ∈ J . Then, the Lagrangian problem is subject to 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 φ * Lagr (L) ≤ φ * holds for any L. Substituting the completion time variables in (10) according to (9a) yields Observe that the walking time and box placement variables occur only in different constraints. Property 1 In the Lagrangian problem, the walking time variables ω i, j , (i, j) ∈ J , and box position variables π i, j , (i, j) ∈ J O (from Definition 4), are independent.
Thus, it is possible to separately optimize walking times and box positions. This gives us the partial objective for the walking times, and for the box positions.

Optimizing box position values
The boxes of the open jobs J O (cf. Definition 4) are to be placed within a box sequence between F and W . In partial objective Π , each box (i, j) ∈ J O adds term aλ i, j − bλ i, j π i, j . Hence to minimize Π , 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).

Optimizing walking time values
The walking time variables ω i, j , (i, j) ∈ J , are optimized by minimizing Ω.
For each (i, j) ∈ J , the value range of ω i, j is limited by constraints (9e). It can be transformed with constants q i, j = −V + k=1,..., j−1 i,k and q i, j = 1 − W + q i, j +V to the range We observe for any i ∈ I , and with increasing j from 1 to n i that term α i, j (as defined in (11)) is nonincreasing and term β i, j (as in (11)) is nondecreasing. Hence, we can find some κ i ∈ {0, . . . , n i } such that α i, j > β i, j for all j = κ i + 1, . . . , n i . Given such a κ i , we can replace constraints (11) by Hence, depending on the value of κ i , either of the two range constraints is active for job (i, j) ∈ J .
We replace the upper bound by an equality with slack variable 0 ≤ y i, j ≤ 1. Then, Property 2 Given κ i , i ∈ I , of an optimal solution. If α i, j < 0 for any job (i, j) ∈ J with j ≤ κ i , then it is possible to decrease κ i without changing the objective Ω such that α i, j ≥ 0 for each job (i, j) ∈ J with j ≤ κ i .
Proof Given the described case, then κ i ≥ 1 and α i,κ i < 0 because α i, j is decreasing with j. Hence, y i,κ i = ω i,κ i = 0 to fulfill constraints (12). Let us decrease κ i by one. Then, it is possible to leave y i,κ i +1 = ω i,κ i +1 = 0 in the solution. Thus, Ω remains unchanged. We repeat this step until α i,κ i ≥ 0.

Corollary 1 An optimum solution exists where
Lemma 2 For each i ∈ I , let κ * i be the minimum value for κ i that permits an optimum solution. Then, there exists such a solution where for each job (i, j) ∈ J there is Proof For each model i ∈ I , we show this by induction for j = n i , . . . , 1. Then, By choice of κ i and use of Property 2, the slack variable y i, j multiplies a nonnegative value. Hence, ω i, j is nonnegative. Moreover, ω i, j influences ω i,k for each k = j + 1, . . . , n i unless y i,k = 0. Thus, ω i, j contributes to the objective Ω not only with factor θ i, j , but moreover via ω i,k with factor y i,k c i,k θ i,k . The total contribution of ω i, j to Ω is thus with factor θ i, j + k= j+1,...,n i y i,k c i,k θ i,k . If this factor is positive, then the lowest slack value y i, j = 0 minimizes Ω. 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.
Let κ i ∈ {0, . . . , n i } for each i ∈ I be the maximum κ i for which −aq i,κ i ≥ 0 holds.

Property 3 An optimum solution exists where
Proof Assume we are given an optimum solution. Naturally, all walking time variables ω i, j , (i, j) ∈ J , are nonnegative. For each i ∈ I , both k=1,..., j−1 ω i,k and q i, j are nondecreasing with respect to j, while −aq i, j is nonincreasing.
. . , κ i . However it is, according to Property 2, possible to set κ i < κ i such that the solution is optimum and −aq i, j < 0 as well as α i, j ≥ 0 hold for any The presented results allow us to describe an algorithm to minimize the walking time variables.

Lemma 3
In an outer loop, set κ i = 0, . . . , κ i for i ∈ I . Given κ i , Lemma 2's recurrence (13) is used to set y i, j for each (i, j) ∈ J , which also sets ω i, j and objective value Ω. 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 O i∈I n 2 i ≤ O n 2 . 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 φ * Lagr (L) is computed in O n 2 time.

Branch and bound methods
For solving Pm 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.

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 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.
Algorithm 1 WNID: Weighted nearest identity sequence placement heuristic 1: ( j i ) i∈I ← 1 set the next box of each model 2: (t i ) i∈I ← 0 set the current time in each model 3: F ← V set the next position to place a box 4: loop |J | times 5: F ← F + w i, j increase F to the next free space 10: 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).

Branch and bound algorithm
An exact solution for Pm 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 {π j } j∈J F of some fixed jobs J 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 π i, j values of all fixed jobs (i, j) ∈ J 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) ∈ J 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) ∈ J O , then its box position is temporarily set to 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 multipli-ersL for step size s > 0 iŝ for each i, j ∈ J . For the step size, we employ the common where UB ≥ φ * is an upper bound value for φ * , and v ∈ (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 φ * 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 max 1, 4 · |J O | times; except at the root node (with |J 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 O . This motivates our use of π i, j , (i, j) ∈ J O , as a traversing order value. We use the nondecreasing order of π i, j for ranking the open boxes and visit the descending nodes accordingly.

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 for constant ψ > 0 and σ > 0. While the number of descending nodes is restricted by σ near the root of the search tree, it is restricted by ψ 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: This yields two heuristic objective values φ (a) and φ (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 F of fixed jobs.

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.

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.
(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. (2008Boysen et al. ( , 2009aBoysen et al. ( , 2009b 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 → ∞, the expected value for each share is 1/m. (3) The number of jobs is set to n ∈ {8, 12, . . . , 28}. A job's model is selected with uniform probability. Hence, each model's expected number of jobs is n/m. Md, median runtime in seconds; Q75, upper quartile in seconds; solved, percentage of instances solved in 1 h * aggregate of all variants of the respective parameter (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, . . . , n)), (L3) uniform random variates from {1, 2, . . . , 10}, (L4) geometric random variates with mean λ = 2.
(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 φ. Otherwise, the problem is likely easier, as many start times will be after the box position (if V 0 or W φ * , cf. Sedding (2020b)) or before (if V 0 or W φ * ). Note it is not necessary to compute the optimum φ * during harmonization because strict equality of W and φ * is not required. Thus, we can use the rather intuitive solution of the WNID heuristic. Then, the deviation |φ−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.

Test setup
Performance comparablity to the single-model study in Sedding (2020c) is ensured by, besides using the same mixedinteger 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 ψ = 5, σ = 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.

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 Pm in Theo- rem 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.
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: - Concluding, the B&B is clearly the best performing exact approach for solving the Pm 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.

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.

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., 2009cBoysen et al., , 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).

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 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 PE(φ, φ * ) of the attained objective φ and an instance's optimum φ * . 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.
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 φ 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 φ can then be assessed with the percentage of instances that yield φ ≤ φ H , and with the mean of the percentage walking time error PE(φ, φ 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.

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 Table 6 Heuristic solutions on instances solved optimally by B&B, grouped by n, m, 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.
Funding Open access funding provided by ZHAW Zurich University of Applied Sciences.

Conflict of interest
The author has no competing interests to declare that are relevant to the content of this article.
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecomm ons.org/licenses/by/4.0/.