Exact and heuristic solutions 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. 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 minimized. In this paper, we first present a mathematical model and improve it with preprocessing techniques. Then, we propose a constructive heuristic and a variable neighborhood 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 2014and Dolgui et al. 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 optimizing 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 minimizing the number of stations for given cycle time, whereas the SALBP-2 aims at minimizing 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 generalizes the SALBP-1. It considers task 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 taking shorter times. The goal is to find a minimum-cost 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., Özceylan 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 neighborhood descent (VND) algorithm (Mladenović and Hansen 1997) composed of two neighborhood 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 produces solutions that are just 2% away, on average, from the best-known solution values, requiring short computing times.
The remainder of the paper is organized 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 labor-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 provided 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 graph G = (T, E), where E = {(i, j) : i j}. Figure 1 gives an example (extracted from a benchmark by Sungur and Yavuz 2015) of a precedence graph involving 13 tasks and 10 precedences.
Let K = H = {1, 2, . . . , l} be two coincident set 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 less 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 characterized 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 a task. In our notation, this is simply represented in the time matrix by a value +∞.
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 minimize 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 assemblyline 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 labor-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). Generalizations of the SALBP involving different task precedence constraints can be found in Dell'Amico, Díaz 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 minimizing 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. We refer to De Bruecker et al. (2015) 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 centers 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 minimize the required cycle time. The authors solved the problem by utilizing a mathematical model. Since the ALWABP is NP-hard, heuristic methods have been developed for its solutions. 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;Moreira and Ritt 2019), and variable neighborhood search by (Polat et al. 2016). Concerning exact methods for the ALWABP, we highlight the branch-and-bound approaches proposed by Miralles et al. (2008), Borba and Ritt (2014), and Vila and Pereira (2014). The particular ALWABP case in which the heterogeneity derives from the presence of robots is known as robotic assembly line balancing problem (RALBP). Cost-oriented approaches, which minimize 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 labor 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) minimized the labor 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-day-a-week industries, Narasimhan (1997) implemented an algorithm that minimizes the labor 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,Özgüven and Sungur (2013) proposed five mathematical models to solve an HWS variant with divisible and indivisible tasks, by minimizing 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, they 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-size instances with 20 tasks but could not prove optimality for a number of 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 all successors of a task i, i ∈ T . Consider for example Figure 1, where we have , 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. takes the value 1 iff station s ∈ S is opened; x ihs ∈ {0, 1} takes the value 1 iff a worker of type h performs task i in station s; y hs ∈ {0, 1} takes the value 1 iff a worker of type h is assigned to station s; z h ∈ Z + gives the total number of workers of type h assigned to the stations.
The model by Sungur and Yavuz (2015), called MSY for short, is as follows: subject to: h∈H,h≤ki s∈S The objective function (1) minimizes 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 pairs of 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 use two preprocessing techniques, both aimed at including additional information, in terms of new arcs and precedence weights on arcs, in the input precedence graph G = (T, E). 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 sequence of precedences. The first technique creates an auxiliary graph G ′ , starting from the original graph G, as follows.
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 station 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 with its corresponding arc A ′ ij . 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 the upper bound of workstations.
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 G = (N ,Ā), based on the LB ij values. The arc values of G ′ are first copied toḠ. 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 workstation e i in which any task i can be processed. With the second step, we compute the latest workstation 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 workstations, ∀i ∈ T .

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 workstations in the model to m = m ub . Then, we apply two reformulations to 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: i∈T :h≤ki 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 (18), 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: 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 in the following, is to minimize (1) subject to (2), (3), (5), (8)-(13), (18) and (23)-(25).

Heuristic algorithms
In this section, we present a constructive heuristic based on the use of different priority rules, and then improve it employing a VND algorithm.

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 iki 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 iki + j∈F * i t jki 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 sucessors processing time: select a task i with largest |F i |/(t ih + j∈S * 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 iki , with k i = h; (t.13) shortest worker-type processing time: select a task i with shortest t iki , with k i = h; where for positional weight we simply mean the processing time of the job 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 divide the tasks not assigned yet by their larger worker type k i . Let λ h be the approximated solution cost in case no further worker of type h would be chosen and each remaining task i would be assigned to a larger worker type h ′ = h, h ′ ≤ k i . Let c curr be the cost of the partial solution under construction and recall c h is the cost of worker type h. We select the worker type h that produces the smallest {c curr + c h + λ h } value. (w.2) smallest cost per task assigned: select a worker type h that has the smallest c h /|T h | value; (w.3) smallest time per task assigned: select a worker type h that has the smallest c h /( i∈Th t ih ) value, where T h is the set of tasks assigned to worker type h; (w.4) largest time per workstation: select a worker type h that has the largest i∈Th t ih value.
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 s, each step of the main loop operates as follows. Line 5 defines the setT , composed by candidate tasks that respect the precedence constraints at the current station. Line 7 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 is used at Line 9 to select the best worker according to its cost and its subset T h . Lines 10-13 perform the necessary updates on the data structures. The algorithm finally returns the number s of stations used, the set W of selected workers, and the task assignment S. return sol CH ={s, W, S} 16: end procedure

Improving heuristic solutions
We developed a VND procedure that makes use of two MILP-based neighborhood procedures. 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. The initial solution sol CH is copied into the incumbent sol best , which is then modified iteratively by means of k max = 2 neighborhood procedures N k . Such procedures are described in detail in the next two sections. return sol best 14: end procedure

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 workstations in the input solution. Suppose task i was assigned to station s i , and let N i be the set of stations to which i can be reassigned. Formally, we set: Procedure MILP-1 finds the reallocations of tasks that is optimal with respect to these limited sets by solving the following model: subject to (3)-(12) and: h∈H,h≤ki s∈N (i) 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 number 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: subject to (3)-(11) and:m s=1 i∈T h∈H,h≤ki The objective function (28) maximizes 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 β.

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 seconds for each execution of the MILP-based neighborhoods of Section 5.2.
The tests were done using the small-size and medium-size instances by Sungur and Yavuz (2015), as well as new large-size 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, we adopted, in this section, 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-size (n = 20) and 45 medium-size (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 small-size and 225 medium-size 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-size ALBHW instances having n = 100 tasks and 225 very large-size 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. 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 small-size 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 just 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-size 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-size and medium-size instances.
Overall, we can conclude that the two models are very effective in solving the smallsize instances, but show important drawbacks for the larger-size 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 = (U B V N D − U B CH )/U B V N D * 100, and the average execution time in seconds (time(s)). 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 very high 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 behavior is also confirmed by the computational efforts, which are quite limited. Indeed, CH never requires more than 0.1 second (not reported in explicit in the table), and VND about 35 seconds on average and around 110 seconds for the most difficult group of instances (having n = 250, w1 = 1.20 and w2 = 0.70). 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 summarized 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 to one another 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 for each worker rule 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 3 and 4, 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 3 we show the overall gap for each group (w1, w2). As the analyse by size is also interesting to show the effectiveness of the procedures, we used the Figure 4, organized 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 (U B H − U B M CIM )/U B 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 minimize 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 neighborhood 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.