Mathematical models and heuristic methods for the assembly line balancing problem with hierarchical worker assignment

This paper proposes new algorithms for the assembly line balancing problem with hierarchical worker assignment (ALBHW). The ALBHW appears in real industrial contexts, where companies deal with a multi-skilled workforce. It considers task execution times that vary depending on the worker type to whom the task is assigned. Qualification levels among workers are ranked hierarchically, where a lower qualified worker costs less but requires larger execution times then a higher qualified one. The aim is to assign workers and tasks to the stations of an assembly line, in such a way that cycle time and precedence constraints are satisfied, and the total cost is minimised. In this paper, we first present a mathematical model and improve it with preprocessing techniques. Then, we propose a constructive heuristic and a variable neighbourhood descent that are useful to solve large instances. Extensive computational experiments on benchmark instances prove the effectiveness of the algorithms.


Introduction
Assembly lines are essential parts of manufacturing systems for large-scale production (see, e.g. Becker and Scholl 2006, Cortés, Onieva, and Guadix 2010, Otto and Otto 2014Scheduling in Production, Supply Chain and Industry 4.0 Systems by Optimal Control: Fundamentals, Stateof-the-Art and Applications," 2019). An assembling process consists of a fixed set of tasks (indivisible elements), which need to be executed on stations while taking into account technological precedence relations. The simple assembly line balancing problem (SALBP) is the basic problem encountered when optimising assembly systems. It considers straight assembly lines and aims at assigning tasks to stations in such a way that task precedence constraints are satisfied, and the total amount of time for each station does not exceed a maximum productive rate called cycle time. Two main SALBP variants have been considered in the literature: the SALBP-1 aims at minimising the number of stations for a given cycle time, whereas the SALBP-2 aims at minimising the cycle time for a fixed number of stations. We refer the reader interested in the SALBP literature to Baybars (1986), Scholl and Becker (2006), Becker and Scholl (2006), Boysen, Fliedner, and Scholl (2007), and Battaïa and Dolgui (2013), among others. In this paper, we study the assembly line balancing problem with hierarchical worker assignment (ALBHW). The problem has been introduced by Sungur and Yavuz (2015) and generalises the SALBP-1. It considers task processing times that depend on the workers, and workers that are divided into hierarchical skill levels. Workers need to be assigned to a station that contains tasks that are compatible with their skills. A worker with higher qualifications incurs a higher cost but executes tasks in shorter times. The goal is to find a minimumcost assignment of tasks and workers to the stations of the line while satisfying maximum cycle time as well as precedence constraints. The problem is NP-hard in the strong sense. As stressed in Sungur and Yavuz (2015), applications of the ALBHW are numerous and can be found in contexts such as health-care, military, manufacturing, and maintenance systems. In addition, studying hierarchical worker assignment is of interest not only for the ALBHW but also because this characteristic might be found in many other types of assembling and disassembling problems (see, e.g. Ozceylan et al. 2019), just to cite some.
Our contribution to the solution of the ALBHW is threefold. First, we propose a new formulation that strengthens the mathematical model by Sungur and Yavuz (2015). Second, we implement a station-based constructive heuristic that combines a number of task priority and worker priority selection rules. Third, we present a variable neighbourhood descent (VND) algorithm (Mladenović and Hansen 1997) composed of two neighbourhood procedures, both based on the use of mixed integer linear programs (MILP). The computational tests that we performed on the benchmark problem instances show that the new formulation improves the previous one for what concerns both the optimality gap and the execution time. The improvement is sometimes small but consistent on all attempted instances. Instances with up to 50 tasks can be solved to proven optimality in a matter of seconds, whereas for large-sized instances, involving 100 tasks, we need to content us with heuristic solutions. The VND obtained the best overall performance, producing an average improvement of 1% from the best-known solution values and an average optimality gap of 6.35%, requiring short computing times.
The remainder of the paper is organised as follows. In Section 2, we present the formal definition of the ALBHW. Section 3 surveys the related literature, focusing on multi-skilled and heterogeneous worker assignments, and on labour-cost variants. Section 4 presents the mathematical model by Sungur and Yavuz (2015), as well as our improvements. The heuristic procedures that we implemented are described in Section 5, and the computational experiments are discussed in Section 6. Final remarks and conclusions are given in Section 7.

Formal problem definition
We consider a single product straight assembly line. Let T = {1, 2, . . . , n} be the set of indivisible tasks to be allocated in the line, and S = {1, 2, . . . , m} be the set of stations, where m is an upper bound on the number of stations required to process all tasks (e.g.m = n). We use the notation i j to state that task i must be executed before task j in the line, i.e. the index of the station processing i should be lower than or equal to the index of the station processing j. The task precedence network is represented by a topological ordered direct graph G = (T, E), where E = {(i, j) : i ≺ j}. Figure 1 gives an example (inspired from a benchmark instance by Sungur and Yavuz 2015) of a precedence graph involving 13 tasks and 10 precedences.
Let K = H = {1, 2, . . . , l} be two coincident sets denoting, respectively, the set of task types and that of workers types, where l stands for the number of qualification levels. We use index k to denote a task type, and index h to denote a worker type. In this notation, the workers of type 1 are the most qualified, while those of type l are the least qualified. Each task i ∈ T has a type k i ∈ K, and can only be executed by workers of type h satisfying h ≤ k i .
Each task i is also characterised by an execution time t ih that depends from the worker type h that executes i. We assume that the execution time of each task increases when lower-qualified workers are employed, that is, t i1 ≤ t i2 ≤ · · · ≤ t il . The right top part of Figure 1 gives an example of execution times. If a task requires a minimum qualification level higher than the one possessed by  a worker, then the worker cannot execute such task. In our notation, this is simply represented by an extremely large value in the time matrix.
The cost of assigning a worker of type h to a station is defined by c h , for h ∈ H and is proportional to its qualification level, so that c 1 ≥ c 2 ≥ . . . ≥ c l . Given a cycle time C, that is, a maximum total execution time allowed at each station, the ALBHW is to assign workers and tasks to stations in order to minimise the total cost, while satisfying constraints on qualifications levels, cycle time, and task precedences. In Figure 2, we present an optimal solution for the instance given in Figure 1. We use h(s) to denote the worker type assigned to station s, andC(s) to denote the total execution time computed according to the worker/task assignment in the station. In this example, we have C = 600, c 1 = 100, c 2 = 70, and c 3 = 49, and so the solution cost is 319.

Brief literature review
The interest of the scientific community in incorporating heterogeneity in assembly-line problems has increased over the years. Surveying the entire literature produced in the field is out of the scope of our work. Nevertheless, in this section, we provide the reader with a brief literature review that focuses on those aspects that are most relevant to the problem we tackle. In particular, we focus on articles that take into account (i) manufacturing systems with skilled level workers, (ii) heterogeneity over task execution times for each worker, and (iii) hierarchical labour-cost assembly lines. We refer the reader interested in the vast area of assembly line systems to the surveys by Becker and Scholl (2006), Cortés, Onieva, and Guadix (2010), Battaïa and Dolgui (2013), and Otto and Otto (2014). Generalisations of the SALBP involving different task precedence constraints can be found in Dell ' Amico, Díaz, and Iori (2012) and Kramer, Dell' Amico, and Iori (2017), among others. These last works also gave inspiration to the improvement techniques that we present below in Sections 4.1 and 4.2. Mansoor (1968) introduced worker performance skills to define task execution times in assembly line systems. He also proposed a constructive heuristic based on positional weighting criteria. Bartholdi III and Eisenstein (1996) considered a case study at Toyota Sewn Products Management System, where workers were classified according to different speed levels. Their results show that sequencing workers from the slowest to the fastest can produce a reasonable line balance. Gel, Hopp, and Van Oyen (2002) considered the possibility of sharing tasks among workers whose task execution times differ by a multiplicative factor. Their study uses a Markov decision process to find good-quality solutions. This approach was later extended by Hopp, Tekin, and Van Oyen (2004), who included cross-training schemes during task assignments. Techawiboonwong, Yenradee, and Das (2006) presented a mathematical model that considers the presence of both permanent (skilled) and temporary (unskilled) workers, as well as specific constraints on the worker assignment on some stations. In the same context, Corominas, Pastor, and Plans (2008) applied a model to solve a rebalancing problem arising in motorcycle assembly lines. They considered given incompatibilities among groups of tasks and aimed at minimising the number of temporary workers employed. In their experimental tests, optimal solutions were found for some realistic instances within a limited computational effort. Koltai, Tatay, and Kalló (2014) presented a similar approach, but focused, instead, on workforce skills in the context of bicycle assembly lines.
Studies considering paced lines and workforce assignments have also attracted the attention of the scientific community. Battaïa et al. (2015) studied a paced line with products of different types, where a sequence of products is given, and task assignments and sequencings are known for all stations. Workers are identically skilled, and their assignment to stations can be changed dynamically depending on the current state of the line. The objective is to minimise the number of workers. The authors proposed a model and deterministic and randomised heuristics. Computer experiments proved the efficiency of the proposed approaches. Dolgui et al. (2018) studied the possibility of having parallel stations at the line design. The problem they addressed considers an upper bound on the number of workers assigned to an operation, precedence relationships between tasks, and cycle time constraints. The processing time of an operation is dependent on the number of workers assigned to it. The objective is to minimise the number of workers allocated. The authors provided proof of NP-hardness, a mathematical model, a set of heuristics, and a series of computational experiments over real and simulated data. They also described problem similarities with the multi-mode project scheduling and multiprocessor scheduling. Delorme et al. (2019) studied a problem that appears as part of a complex industrial environment. They considered paced mixed-lines, with fully skilled workers, with the objective of minimising the number of workers allocated. Besides a new mathematical formulation, they introduced an NP-hardness proof, a dynamic programming algorithm, and the description of polynomial and pseudo-polynomial algorithms for special cases. Lian et al. (2018) considers a bi-objective approach for a worker assignment problem in seru production systems. These assembly lines have a new type of work-cell-based manufacturing system, where most of the tasks are manual, assisted by hand tools, or lightweight equipment. Workers differ from their skills and proficiency levels. They study proposed a mathematical model, minimising the deviation workload from workers and the serus installed at the line. They also propose an adaptation of the NSGA-II algorithm to solve 100 medium-sized instances. According to the tests, heterogeneous workers with diversified competency perform well in balancing inter-seru workload.  also included workforce satisfaction in their solution methods. They proposed a mathematical model and an -constraint method. The authors reported a Pareto front obtained with a small-sized example with 10 tasks, 6 parallel machines, and 3 workers.
We also mention studies concerning multi-skilled workers at U-shaped lines (Nakade and Nishiwaki 2008), two-sided assembly lines (Tapkan, Özbakır, and Baykasoğlu 2016), project scheduling problems (Karam, Attia, and Duquenne 2017), hybrid flow shop environments (Liu and Yang 2019), automotive contexts (Martignago, Battaïa, and Battini 2017;Yang, Lee, and Kim 2020), and even problems solved as a generalised assignment (Moussavi, Mahdjoub, and Grunder 2018). We refer to De Bruecker et al. (2015) and Dolgui et al., "Workforce Planning and Assignment in Mixed-Model Assembly Lines as a Factor of Line Reconfigurability: State of the Art," (2019) for a survey on workforce planning under different skills. Miralles et al. (2007) introduced the assembly line worker assignment and balancing problem (ALWABP), as a problem arising in the context of sheltered worker centres for disabled (SWD) in Valencia (Spain). In the ALWABP, task execution times vary according to the workers, not all workers can perform all tasks, and the aim is to minimise the required cycle time. The authors solved the problem by utilising a mathematical model. Since the ALWABP is NP-hard, heuristic methods were also developed. Among these, we mention clustering search by Chaves, Miralles, and Lorena (2007), beam search by Blum and Miralles (2011), tabu search by Moreira and Costa (2009), genetic algorithms by Moreira et al. (2012), Mutlu, Polat, and Supciller (2013) and Moreira and Ritt (2019), and variable neighbourhood search by Polat et al. (2016). Concerning exact methods for the ALWABP, we highlight the branch-andbound approaches by Miralles et al. (2008), Borba and Ritt (2014), and Vila and Pereira (2014).
Concerning recent extensions of the ALWABP, Oksuz, Buyukozkan, and Satoglu (2017) proposed a mathematical model, a genetic algorithm, and an artificial bee colony algorithm for the U-shaped ALWABP type-E. Computational tests were performed with three data sets. The authors indicated as further research avenues the study of new algorithms and of a mixed-model variant of the problem. Janardhanan, Li, and Nielsen (2019) addressed the two-sided ALWABP. They proposed a mathematical model that effectively solves (through CPLEX) small-sized instances with less than 16 tasks, and a migrating birds optimisation (MBO) algorithm. The MBO computationally outperformed other metaheuristics as partial swarm optimisation, genetic algorithm, artificial bee colony algorithms, and simulated annealing. The particular ALWABP case in which the heterogeneity derives from the presence of robots is known as the robotic assembly line balancing problem. Cost-oriented approaches, which minimise the robot selection and task allocation costs, have been proposed by Rubinovitz, Bukchin, and Lenz (1993), Bukchin and Tzur (2000), and Pereira, Ritt, and Vásquez (2018).
Hierarchical workforce scheduling (HWS) is a class of problems that has, as its basis, the objective to find the most economical mix of workers that satisfies given labour requirements and worker characteristics. Since the seminal work by Emmons and Burns (1991), several authors have attempted to solve HWS problems. Hung and Emmons (1993) proposed a mathematical model for an HWS variant with compressed workweek arrangements and multiple shifts (rotation and downward worker substitutability). Hung (1994) minimised the labour costs in an HSW problem through a heuristic method, considering the possibility of a single-shift off-day. Still taking a single-shift off-day in seven-daya-week industries, Narasimhan (1997) implemented an algorithm that minimises the labour costs subject to a pre-defined demand of workers and desired work characteristics. Seçkiner, Gökçen, and Kurt (2007) and Pastor and Corominas (2010) extended the approach by Billionnet (1999), who studied an HWS case with variable demand, by including compressed workdays and suitability of task assignments. More recently, Ozguven and Sungur (2013) proposed five mathematical models to solve an HWS variant with divisible and indivisible tasks, by minimising total costs. As far as we know, Sungur and Yavuz (2015) is the only work directly focused on the ALBHW. Apart from introducing the problem, the authors also presented a mathematical model and performed extensive tests on a set of instances adapted from Otto, Otto, and Scholl (2013). Using CPLEX 12.3, their formulation found optimal solutions for small-sized instances with 20 tasks but could not prove optimality for several medium-sized instances involving 50 tasks. These results motivate our work on the development of improved models and new effective heuristics for the problem.

Mathematical models
We first introduce the MILP model by Sungur and Yavuz (2015) and then discuss our improvements. We follow the notation given in Section 2. In addition, we use P i ⊆ T to denote the set of direct predecessors of task i, i ∈ T, and F i ⊆ T to denoted the set of direct successors of task i, i ∈ T. We also compute the transitive closures of these sets, and use P * i and F * i to denote, respectively, all (direct and indirect) predecessors and successors of a task i, i ∈ T. Consider for example Figure 1 7, 8, 9, 11, 12, 13}. The MILP model that we developed is based on the four families of decision variables given in Table 1, three of which are binary and one integer.
The model by Sungur and Yavuz (2015), called for short MSY (Model Sungur Yavuz), is as follows: subject to: h∈H,h≤k i s∈S The objective function (1) minimises the total cost of the assembly line. Constraints (2) guarantee that each task is assigned to a station. Constraints (3) impose the assignment of exactly one worker to a station if such station is opened. Coupling constraints between x ihs and y hs variables are ensured by (4). Constraints (5) determine the total amount of workers of type h that are assigned to the stations. Constraints (6) establish the precedence relations between the tasks. Constraints (7) prevent the total workload of each opened station to exceed the cycle time. Inequalities (8) allow a station s + 1 to be opened only if the previous station s is also opened. Constraints (9)-(12) define the domains of the variables. This model has O(|T||H||S|) variables and O(|T| 2 + |H||S|) constraints.

Preprocessing
We used two preprocessing techniques, both aimed at strengthening the mathematical formulation. In the first technique, we calculate the lower bound on the number of stations that separate every two tasks in E. We then create a new graph and insert arcs having precedence weights equal to the computed lower bound values. The second technique uses the lower bound values computed by the first technique to calculate the earliest and latest possible station for each task in the graph. Recall, indeed, that G does not have precedence weights, so, if (i, j) ∈ E, then i and j can be assigned to the same station. The inclusion of new arcs and precedence weights can be beneficial for creating stronger cuts and valid inequalities (to be described next in Section 4.2). Both our two preprocessing techniques share this scope, but the first technique builds upon the cycle time restriction, whereas the second uses the precedence. The first technique creates an auxiliary graph G , starting from the original graph G, as follows. Let G ij = (N ij , A ij ) be a subgraph of G , where N ij denotes the subset of tasks in all paths connecting two tasks i, j ∈ N having In addition, let A ij denote the subset of arcs that belong to all paths from N ij . Let LB ij be a lower bound on the number of stations required to process all tasks in G ij . We compute LB ij = k∈N ij t k1 /C − 1 as the continuous lower bound on the processing times of the tasks, by considering the case in which they are executed by the fastest worker (of type 1). We then associate each LB ij value with its corresponding arc A ij . If F * i ∩ P * j = ∅, then we set LB ij = 0. These values are used to produce improved inequalities and are the starting point of the next preprocessing technique. Let Z vnd be the cost of the solution found by the VND procedure. We define m ub = Z vnd /c l , which is an upper bound on the number of stations.
In the second technique, we solve the longest path problem using the classical critical path method (CPM) to obtain the earliest and latest station index to which task i ∈ T can be assigned, denoted by e i and l i , respectively. We create a new auxiliary graphḠ = (N,Ā), based on the LB ij values. The arc values of G are first copied intoḠ. LetN = N ∪ {0, n + 1}, where 0 and n+1 are two dummy tasks that represent the beginning and the end of all activities, respectively. To create the new arcs in the graph, we first updateĀ = A ∪ {(0, i)}, ∀i ∈ T having P i = ∅, thus imposing 0 as source vertex. We then updateĀ =Ā ∪ {(j, n + 1)}, ∀j ∈ T having F j = ∅, thus imposing n + 1 as sink. After this update, we have a graph in which all vertices are connected directly or indirectly with the source and the sink. With the first step of the CPM, we are able to define the earliest station e i in which any task i can be processed. With the second step, we compute the latest station l i in which i can be processed. In detail, the value of l i is obtained by first setting l n+1 = m ub . Then, we take the transpose graph of G, named G t , and compute the earliest station of its tasks,ẽ i , by the CPM. The output of the preprocessing procedure is an interval [e i , l i ] of stations, ∀i ∈ T. Figure 3 shows the new graphḠ created by the preprocessing procedures. The dashed arcs were added by the CPM procedure. Each pair of values near a node gives the indices of the earliest and latest stations, respectively, for the node. The red arc is the lower bound number of stations between tasks 2 and 11, given by LB 2,11 = 1 and LB 2,12 = LB 2,13 = 0. Note that we omitted all the arcs that produce LB ij = 0 to simplify the example.

Improvements to the mathematical model
We developed some techniques for improving the linear relaxation of the MSY model and removing useless variables. First, we set the number of stations in the model to m = m ub . Then, we apply two reformulations of constraints (4) and (7). We start by strengthening n in constraints (4) by a sufficient large constant M h ≤ n, thus obtaining: The value of M h is computed by considering the maximum number of tasks that can be assigned to a station, while satisfying the cycle time C without precedence constraints and considering an assignment of only h-type workers to the stations. Such value is obtained, for each h ∈ H, by solving the MILP problem: subject to: where α i is a binary variable that takes the value 1 if task i is selected in the subset of maximum cardinality, and 0 otherwise. Model (14)- (17) is a maximum cardinality knapsack problem with additional incompatibilities imposed by (16). The maximum cardinality knapsack problem is easy, but the addition of the incompatibilities makes it an NP-hard problem. Nevertheless, the above MILP model can be solved very quickly for instances of medium size, as the ones that we face in our tests. We use a similar technique for constraints (7), replacing the original t ih values with new ones, t ih ≥ t ih , resulting in: Each new value is obtained, for i ∈ T and h ∈ H, by solving a subset sum problem (SSP): given a station of capacity C − t ih and a set of tasks T\{i}, determine the subset of tasks not exceeding the capacity but having largest total time. The SSP solution value, denoted by ih , is obtained by solving the following model: subject to: where again α j takes the value 1 if task j is selected, and 0 otherwise. The SSP is known to be NP-hard, but, model (19)-(22) too can be quickly solved for instances of medium size. Once the problem is solved, we set t ih = C − ih . We also modify precedence constraints (6) Furthermore, for each station s that has the index smaller than the earliest possible or greater than the latest possible, we include the following reduction: x ihs = 0 ∀i ∈ T, h ∈ H, s ≤ e i − 1 or s ≥ l i + 1 (24) Let us consider the number z 1 of workers of type 1 (i.e. the highest qualified ones). We limit by below this value by considering those tasks that can only be executed by workers of type 1 (i.e. those for which k i = 1), and computing a continuos lower bound on their processing times. We consequently add to the model the constraint: The resulting MILP model, called MCIM (from Model Campana Iori Moreira) in the following, is to minimise (1) subject to (2), (3), (5), (8)- (13), (18) and (23)-(25). Our idea, confirmed by the computational results in Section 6.2, is to show how we can improve the existing MSY model and obtain an updated model that has better bounds and optimality gaps. Since our aim is to improve the current model and not create a new one using a different concept and structure, it is expected that the two models adopt different parameter configurations and additional inequalities in the preprocessing step.

Heuristic algorithms
In this section, we present a constructive heuristic based on the use of different priority rules, and then improve it by employing a VND algorithm. Our idea is to design fast yet accurate algorithms to solve the ALBHW. Even though we mainly focus on efficiency, computational results (Section 6) show that the procedures obtained good-quality solutions given the time limit imposed.

Constructive heuristic
The constructive heuristic that we built is inspired by Scholl and Voß (1997) and Moreira et al. (2012) and aims at filling one station at a time. For each station, and for each candidate worker type h ∈ H, the heuristic creates a set of tasks by using an input task priority rule while respecting precedence and cycle time constraints. Then, a worker type and the corresponding set of tasks just created are selected to be assigned to the station according to a given worker priority rule. The process is repeated until all tasks have been assigned. We implemented the following task priority rules: (t.1) max successors: select a task i having largest |F * i | value; (t.2) max direct successors: select a task i having largest |F i | value; (t.3) largest min processing time: select a task i with largest t i1 value; (t.4) largest max processing time: select a task i with largest t ik i value; (t.5) largest min positional weight: select a task i with largest t i1 + j∈F * i t j1 value; (t.6) largest max positional weight: select a task i with largest t ik i + j∈F * i t jk i value; (t.7) largest positional weight: select a task i with largest t ih + j∈F * i t jh value; (t.8) shortest min processing time: select a task i with shortest t i1 value; (t.9) largest direct successors processing time: select a task i with largest |F i |/(t ih + j∈F * i t jh ) value; (t.10) largest successors processing time: select a task i with largest |F * i |/t ih value; (t.11) largest processing time: select a task i with largest t ih value; (t.12) largest worker-type processing time: select a task i with largest t ik i , with k i = h; (t.13) shortest worker-type processing time: select a task i with shortest t ik i , with k i = h; where for positional weight we simply mean the processing time of the task and of all its successors. Other techniques that we attempted did not produce improvements in the computational performance, and were thus disregarded.
In addition, we implemented the following worker priority rules: (w.1) look ahead worker: select a worker of type h, by estimating the total cost in case h would not be chosen anymore. We group the tasks not assigned yet, named T , and allocate them to their worker type k i (the one with highest cost that performs task i). Let λ h be the approximated solution cost in case no further worker of type h would be chosen. Then, taking The pseudocode of the constructive heuristic is given in Algorithm 1. It receives in input the set T of tasks, the worker priority rule, w_rule, and the task priority rule, t_rule. Starting from an empty station, each step of the main loop operates as follows. Line 7 defines the setT, composed by candidate tasks that respect the precedence constraints at the current station. Line 9 creates a subset T h ⊆T, for each worker type h, by respecting the cycle time restriction and using the rule t_rule. If the worker of type h cannot execute a task chosen by t_rule, then a new task is chosen with the largest time for the worker h, always considering the precedence or cycle time constraints. The worker priority rule w_rule is used at Line 11 to select the best worker according to its cost and its task set T h . Lines 12-14 perform the necessary updates on the data structures. The algorithm finally returns the set W of selected workers, and the tasks assigned to the stations S.

Improving heuristic solutions
We developed a VND procedure that makes use of two MILP-based neighbourhood procedures. The studies of James and Almada-Lobo (2011) and Moreira, Miralles, and Costa (2015) inspired the idea to embed MILP solvers' resources in our local search procedures. According to Talbi (2009), the VND heuristic iterates neighbourhood structures successively until it reaches a local optimum. Briefly, consider a neighbour structure N k , an initial solution sol, and a set of neighbours of sol denoted by N k (sol). If through N k (sol) there is no possibility to find a better solution than sol, then we replace N k with N k+1 . Every time we find a better solution than the current one, we return to the first neighbourhood and restart the search. This strategy is more effective when neighbourhoods are complementaries since one neighbour's local optimum is not the local optimal of the other neighbour. In our case, we adopt two neighbourhood structures: the first one performs task relocations (Section 5.2.1), and the second one allows worker and task reassignments (Section 5.2.2).
Throughout the procedure, the upper bound on the number of stations is set to the number of stations produced by the solution sol CH found by the constructive heuristic of Section 5.1. Algorithm 2 shows the pseudocode of the VND. A function f (z) that maps a solution z to an objective function value is used to evaluate solutions. The initial solution sol CH is copied into the return sol best 14: end procedure incumbent sol best , which is then modified iteratively by means of k max = 2 neighbourhood procedures N k . Such procedures are described in detail in the next two sections.

Task relocations
The first local search procedure, called MILP-1, allows tasks to be reassigned to stations that are adjacent to the station in which they were assigned to in the input solution. Letm be the number of stations in the input solution. Suppose task i was assigned to station s i , and let S i be the set of stations to which i can be reassigned.
In practice, this is equivalent to the original model above in Section 4, but its size is consistently reduced by replacing constraints (2) by (27).

Worker and task reassignment
The second local search procedure, MILP-2, allows the tasks to be relocated in any position (while respecting the constraints) but requires that the numbers of workers for each type remain unchanged with respect to the input solution. To this aim, we create auxiliary binary variables β s , for s = {1, 2, . . . ,m}, taking the value 1 if and only if station s has no task assigned to it. We then let α s ∈ R, for s = {1, 2, . . . ,m}, be the cost of the worker assigned to s when β s = 1. Letz h be the number of workers of type h selected in the input solution. We obtain the following model: The objective function (28) maximises the savings obtained by excluding a worker from the assembly line. Inequalities (30) and (31) are the coupling constraints between variables β and x, whereas (32) and (33) relate, respectively, α with y and with β.

A numerical example
In this section, we demonstrate how the heuristic procedures perform. We use the same instance of Section 1, with w 1 = 1.1, w 2 = 0.7 and C = 600. Moreover, we choose the combination of worker and task rules for the configuration (w 1 , w 2 ) denoted by w.3 and t.12, respectively, as an example. First, we describe the behaviour of the constructive heuristic (Algorithm 1). We list each step below.
• We start withT = {1, 2, 3} as the available tasks, based on the precedence graph (Line 7). • In Line 9, the algorithm chooses a set of tasks iteratively, based on the task priority rule. In detail, for each worker, the algorithm selects a subset of tasks respecting the precedence constraints and the maximum cycle time. Since the procedure is iterative, when a task i is selected, we remove the precedence arcs, allowing other tasks to be chosen in the next steps. • Let T 1 = {2, 1, 3, 10}, T 2 = ∅, and T 3 = ∅ be the results of the loop of Lines 8-10 in the heuristic.
• The only possible option is to select h * = 1 (Line 11).
• After selecting task 2, the tasks {7, 8, 5, 9} become available, but they are not considered because all these tasks are of type 2. The same happens when task 3 is selected, unlocking tasks {4, 6, 10}. But note that the last chosen task is task 10, which is of type 2 (Line 7). • If the task rule cannot select any task in the available tasks using the t_rule, the procedure picks the task with largest processing time with the worker h (at line 9) respecting the constraints (Line 8). In the VND procedure (Algorithm 2), we start with neighbourhood N 1 , which is the task reallocation MILP (Line 5). We also list the next steps below, as we did for Algorithm 1. where it is again possible to obtain an improvement (Line 8). • The neighbourhood counters k increases when we have no improvement, proceeding to the next neighbourhood N 2 , which is the worker and task reassignment MILP (Line 10). • We impose that the number of workers of each type is the same from the inputz = {2, 1, 1}, but the model can choose to assign them in any station. Note that in N 2 (Line 5), we solve the same original problem but with a reduced solution space. The model can find the optimal solution if (i) the number ofz h ≥ z * h , ∀h ∈ H, which means that the number of workers type h need to be at least equal to the number of workers type h of the optimal solution, and (ii) if there is enough time to solve it.
• At this point, there are no more possible improvements and the procedure ends with the objective function value of 319 (Line 13). Note that the value 319 is indeed optimal.

Computational experiments
In our computational experiments, we evaluate the results obtained by the methods presented in this paper. The evaluation of the mathematical models is given in Section 6.2 and that of the heuristics in Section 6.3. The algorithms have been coded in C++, using Gurobi v9.0 as MILP solver. The tests have been executed on a PC Ubuntu 16.04.6 LTS 64-bit with Intel(R) Xeon(R) CPU E3-1245 v5 3.50GHz and 32GB RAM. The maximum execution time was set to 2 hours for the two models of Section 4, and to 15 s for each execution of the MILP-based neighbourhoods of Section 5.2.
The tests were done using the small-sized and medium-sized instances by Sungur and Yavuz (2015), as well as new large-sized instances that we created to obtain a better evaluation of our methods, as described in Section 6.1. In total, we solved 900 instances. Because of the large size of the testbed that we adopted, we only present aggregate results. All the instances that we used, together with the detailed computational results that we obtained on every single instance, have been made publicly available at https:// github.com/NicolasCampana/ALBHW-Instances.

Benchmark set
To the best of our knowledge, the only benchmark instances for the ALBHW have been proposed by Sungur and Yavuz (2015). They chose 45 small-sized (n = 20) and 45 medium-sized (n = 50) instances among those created by Otto, Otto, and Scholl (2013) for the SALBP-1. For each such instance, they generated 5 ALBHW instances by attempting all combinations of two parameters, w 1 ∈ {1.10, 1.20} and w 2 ∈ {0.70, 0.85} and including an additional instance in which both w 1 and w 2 were set to 1. In this way, they created 45 × 5 = 225 smallsized and 225 medium-sized instances. Parameter w 1 is the time factor that defines the time scaling from a worker type to the next, used in such a way that t i,h+1 = w 1 t ih . Parameter w 2 is used instead to determine the cost scaling among the worker types, in such a way that c h+1 = w 2 c h .
To obtain a better evaluation of our methods, we also created a new set of 225 large-sized ALBHW instances having n = 100 tasks and 225 very large-sized instances having n = 250. These have been obtained following the same procedure adopted by Sungur and Yavuz (2015). We thus adopted a random process to choose the SALBP-1 instances from Otto, Otto, and Scholl (2013) with mixed precedence graphs, balancing the amount of 'extremely tricky', 'very tricky', 'less tricky' and 'tricky' instances. Additional information about the instances can be found in Sungur and Yavuz (2015) and in our public repository. Table 2 presents the results obtained by the two models, MSY by Sungur and Yavuz (2015) and the newly proposed MCIM. Both models do not have any initial solution as an input. For each model, we give the number of proven optimal solutions (#opt), the number of instances for which at least one feasible solution has been found (#feas), the percentage optimality gaps at the root node (%gap root ) and at the end of the execution (%gap), the number of explored nodes (#nodes) and the total computing time in seconds (time(s)). Each line in the table gives aggregate total (for #opt and #feas) or average (for all other columns) values for 45 instances sharing the same value of n, w 1 and w 2 . The lines called 'average/total' give aggregate values for the 225 instances having the same value of n, and the last 'overall' line presents aggregate results for the entire set of 900 instances. The best %gap values are highlighted in bold.

Evaluation of the mathematical models
Both models were very effective in solving the smallsized instances with just 20 items. Some remarkable differences can be noticed, however, already for the cases with n = 50: MSY can only find 148 optima, with an average effort of 2821 seconds an average gap of 1.24%, whereas MCIM obtains 14 more optima, with just 2330 seconds of average effort and producing an optimality gap that is just 0.54%. The difference becomes more evident for the group having n = 100: MSY finds only 33 proven optima and 147 feasible solutions, with an average gap of about 44%, whereas MCIM obtains nine more optima, 24 more feasible solutions, and a smaller gap of 28.89% on average. When n = 250, MCIM found less feasible solutions than MCY, but still more ones. It also obtained better gaps and larger numbers of explored nodes. In any case, it is evident from the results that heuristic algorithms should be developed to tackle the very large-sized instances. Overall, MCIM is more accurate than MSY on all groups of instances. This is achieved at the expense of a slightly larger number of explored nodes, but with a reduction in the computing time. The good MCIM performance can also be imputed to the fact that it has a better root node gap than MSY, especially on the small-sized and medium-sized instances.
Overall, we can conclude that the two models are very effective in solving the small-sized instances but show important drawbacks for the larger-sized instances, for which it is important to invoke heuristic algorithms. We believe it is also important to notice the impact of the hierarchical worker assignment on the difficulty of the instances. This can be observed by comparing the instances having w 1 = 1.00 and w 2 = 1.00, which correspond to the pure SALBP-1 because all workers are equal one another both in terms of cost and of processing times, to the instances in the other groups, where indeed the hierarchical worker assignment has an impact. Such hierarchy increases the difficulty of the problem by increasing the %gap root values and the execution times consistently. This is particularly evident when n = 100, where %gap root increases from less than 10% to around 35% and 90% in the n = 250.

Evaluation of the heuristic algorithms
In Table 3, we present an aggregate evaluation of our heuristic procedures. To facilitate comparison, for MCIM, we resume from Table 2 the number of feasible solutions found (#feas) and the average gap from the lower bound (%gap). For the constructive heuristic, called CH for short in the following, we provide the average solution value (UB) and again the average percentage gap from the lower bound. For the VND, we provide the same information given for CH, but we also report the average improvement with respect to CH, computed as %impr = (UB VND − UB CH )/UB VND * 100, and the average execution time in seconds (time(s)). We execute the VND with an input coming from the best combination task/worker priority rules of the CH, as a greedy strategy to start the method with a good-quality solution, since it is very fast to run all combinations of priority rules. In any case, if one had to choose a single rule that performs well under various conditions, we would suggest rule w.3 and t.3, due to its performance shown in Table 4. The constructive heuristic is 2.62% away from the lower bound on the instances with 20 tasks, and 3.21% on those with 50 tasks. The VND can improve CH by about 1.41% when n = 20, achieving a gap equal to 1.42%, and by 1.35% when n = 50, achieving a gap equal to 2.48%. In these cases, MCIM appears to be a good solution choice. However, as previously noted, MCIM fails in finding feasible solutions for many instances with n = 100 (namely, 54 instances) and for almost all instances with n = 250 (namely, 202). For all such cases, both CH and VND quickly found feasible solutions. In addition, they could reduce the gap consistently from the lower bound, from the 28.89% obtained on average by MCIM, to just 8.15% on average for CH and 7.27% for VND. For the instances with n = 250, CH reduces the average gap from 90.43% to 14.49%, and VND reduces it even more, arriving at 14.24%. The good algorithmic behaviour is also confirmed by the computational efforts, which are quite limited. Indeed, CH never requires more than 0.1 s (not reported in explicit in the  table), and VND about 35 s on average and around 110 seconds for the most difficult group of instances (having n = 250, w 1 = 1.20 and w 2 = 0.70).
To attest to the quality of the solutions provided by the VND heuristic, we apply the well-known benchmark profiler proposed by Dolan and Moré (2002),   which is a useful tool to visualise different algorithmic performances. To generate the graphics, we use the Julia package called BenchmarkProfiles, available at https://github.com/JuliaSmoothOptimizers/Benchmark Profiles.jl. We compare four methods presented in this study: the two mathematical models (MSY and MCIM) and the two heuristic procedures (CH and VND), in which CH corresponds to the best solution from all combinations of task/worker priority rules. In short, let I be a set of instances, A a set of algorithms, andc pa the total cost obtained by algorithm a ∈ A in instance p ∈ I. The performance ratio of algorithm a ∈ A in instance p ∈ I is defined by: The performance profile ρ a (τ ) of an algorithm a ∈ A is given by: The value of ρ a (τ ) indicates the probability for algorithm a ∈ A that a performance ratio r pa is within a factor τ of the best possible ratio. When a set of instances are compared, the algorithm with a larger ρ a (τ ) is preferred. Figures 4-7 present the results for small, medium, large and very large groups of instances, respectively. We only consider instances for which each of the four methods obtained at least a feasible solution. For small-sized instances (Figure 4), we see a similar performance of MSY and MCIM model. If we consider a loss of 10% efficiency, the VND and the CH heuristics achieve almost 100% of the best-known solutions. An identical performance of the four methods is noted when we allow a loss of 30% efficiency. In Figure 5, we observe the dominant behaviour of the MCIM model, which obtained practically all of the best solutions. The model MSY also has a good performance, getting 100% of the benchmark if we take a loss of efficiency greater than 5%. The VND and the CH needs a loss of efficiency of more than 15% to solve all the problems.
In Figures 6 and 7, we note that the VND heuristic starts to outperform the other methods due to the complexity to obtain high-quality solutions found by the models and the simplicity of CH. Concerning solution quality, these analyses confirm a behaviour we expected, that applying heuristic methods (especially the VND) is more appropriate for large and very large instances, while    the mathematical models provide high-quality solutions for small and medium instances.
In Table 4, we try to obtain some insight into the performance of all task and worker priority rules (see Section 5.1) used by CH. Worker rules are summarised as lines in the table, whereas task rules as columns. Differently from the previous tables, in this case, we opted to aggregate instances in groups (w 1 , w 2 ) having the same values of parameters w 1 and w 2 , independently of the value taken by n. As previously noted, the first group (1.00, 1.00) corresponds to the SALBP-1 because all workers are equal both in terms of cost and of processing times. For this group, all worker rules are equivalent, and we thus concentrate only on the task rules. For the other four groups, where hierarchical worker assignment has instead an impact, we show the results obtained by both task and worker priority rules. Each cell in the table gives the number of times in which a pair of rules led CH to find the best solution cost. Total values are provided for each task rule in the last line, and each worker rules in the right-most column. The best values for each group are highlighted in bold. For group (1.00, 1.00), the best task rule is t.3, which consists of selecting the task with the largest processing time using the worker type 1, which is a relaxation to the SALBP-1. For group (1.10, 0.70), the best result is obtained by the combination of rules w.3 and t.12. This corresponds to selecting the tasks with the shortest processing times and the workers that have the shortest time per task assigned. Task rule t.12 is also very effective for the groups (1.1, 0.85) and (1.20, 0.70), both combined with w.3. Task rule t.3 is the best for the last group (1.20, 0.85) and is the one obtaining the highest total number of best results, 296 out of 900. Overall, all rules have a positive impact on the algorithm. For some of them, this is small (as for t.8 and t.11), but for others, it is quite remarkable (as for t.1, t.3, and t.12). There is no clear winner among the worker priority rules, as all of them have a balanced positive impact on the algorithm.  In Figures 8 and 9, we provide a further analysis aiming at evaluating the impact of each heuristic component, selecting only those instances for which MCIM found a feasible solution. In Figure 8, we show the overall gap for each group (w 1 , w 2 ). As the analyse by size is also interesting to show the effectiveness of the procedures, we used Figure 9, organised by each size n its the overall gaps. We then computed the percentage gap of a certain heuristic configuration, say, H, from the upper bound of the model, as (UB H − UB MCIM )/UB H * 100. We tested CH, CH improved by MILP-1 (run once) and CH improved by the VND. We can notice a decrease in the gap that is consistent on all groups and all sizes.

Conclusions
We studied the assembly line balancing problem with hierarchical worker assignment (ALBHW), a problem that requires to assign tasks and workers to the stations of an assembly line, to minimise costs while satisfying precedence and cycle time constraints. Workers are divided into different qualification levels, and lower qualified workers cost less but also require larger execution times for the tasks. We proposed some improvements to an existing mathematical model, and then developed a constructive heuristic and a variable neighbourhood descent algorithm. The improved mathematical model is quite efficient in solving instances with up to 50 tasks, but for larger instances it is convenient to adopt the heuristic algorithms. Several instances, now made publicly available on the Internet, remain unsolved to proven optimality, and future more powerful exact and heuristic algorithms are envisaged for solving the problem.
The inclusion of the hierarchical worker assignment makes the problem much more difficult to solve. Mathematical models incur, indeed, in much larger root node optimality gaps and, consequently, in much larger execution times. Making use of a multi-skill workforce to its full extent may have a critical impact in several production activities not only assembling ones, but also disassembling and scheduling in general. Interesting future research directions include thus the addition of hierarchical worker assignment characteristics to other classes of problems. The implementation of other neighbourhoodbased metaheuristics such as iterated local search and variable neighbourhood search is also a promising future research topic. Concerning the mathematical models, a reasonable idea to feed a feasible solution into MSY and MCIM formulations to guide the search from the given solution, given that Gurobi supports this idea by MIP start. Finally, based on the successful utilisation of branch, bound, and remember algorithm in other assembly line balancing problems (Vila and Pereira 2014;Li, Kucukkoc, and Tang 2020;Li et al. 2020, e.g.), we believe that an adaptation of such an algorithm could be promising for the ALBHW.