Partially Preemptive Multi Skill/Mode Resource-Constrained Project Scheduling with Generalized Precedence Relations and Calendars

Multi skill resource-constrained project scheduling Problems (MS-RCPSP) have been object of studies from many years. Also, preemption is an important feature of real-life scheduling models. However, very little research has been investigated concerning MS-RCPSPs including preemption, and even less research moving out from academic benchmarks to real problem solving. In this paper we present a solution to those problems based on a hybrid method derived from large neighborhood search incorporating constraint programming components tailored to deal with complex scheduling constraints. We also present a constraint programming model adapted to preemption. The methods are implemented in a new open source python library allowing to easily reuse existing modeling languages and solvers. We evaluate the methods on an industrial case study from aircraft manufacturing including additional complicating constraints such as generalized precedence relations, resource calendars and partial preemption on which the standard CP Optimizer solver, even with the preemption-specific model, is unable to provide solutions in reasonable times. The large neighborhood search method is also able to find new best solutions on standard multi-skill project scheduling instances, performing better than a reference method from the literature.


Introduction
Multi skill resource-constrained project scheduling problems (from hereon MS-RCPSP) are critical in business applications like automotive, human resources, aerospace or nuclear industry [3]; managing efficiently the workforce assigned to perform required tasks directly impacts the performance and progress of such businesses.This type of problems have been widely studied in academia from long, as well as the more general yet less compact multi-mode model involving non renewable resources [36].Also, preemption in scheduling is found in many works due to the flexibility it brings to improve the performance of a practical application.However, the scheduling literature including both multi skill and task preemption, all the more grounded in a real business application, is rather sparse [1,34].
The presence of multiple skilled resources adds more complexity to the problem, and new types of goals aside the usual makespan and cost optimization, like assignment score and worker performance [35].Task preemption adds another layer of complexity, and in real scenarios usually is tied to resource availability, with special focus in worker shifts and calendars.Also, it is important to note that whilst a task can be preempted, pausing its progress, maybe only part of the resources can be freed to make them available for other tasks, with others being locked in place.This has been already referred to as partial preemption [25].Finally, practical problems also contain generalized precedence constraints and resource calendars, which have been well known from long ago.These constraints extend simple tasks precedence to detail concrete temporal conditions between tasks.Calendars on the other hand, define resource unavailability periods, hence both features make the problem harder to solve [17].
In this paper, we present a method oriented towards a real assembly line use case, requiring to include all these five features: multi skilled resources, multiple modes, (partially) preemptive tasks, calendars and generalized precedence constraints.Hence, we call our model the partially preemptive-multi-skill/mode resource-constrained project scheduling problem with generalized precedence relations and resource calendars (PP-MS-MM-RCPSP/max-cal).Our solution is aimed for industrial domains, but can solve generic problems in the scope of the model as well.We propose an extended version of the large neighborhood search (LNS) method [21,11], softening constraints initially to find a base solution that will be improved on each iteration.We observed this approach to be faster than other more direct constraint programming approaches.We evaluated our method against different benchmarks, including one using real data from aerospace manufacturing plants, real instances from an hazardous material examination facility [25] and standard multi-skill scheduling instances [38].
The rest of the paper is structured as follows: in Section 2, literature about multi skill scheduling and preemption in scheduling is discussed, along with previous scheduling works oriented towards industrial applications.Section 3 describes the domain of our problem and contains the formulation of our model, including a new CP formulation tailored to preemption.In Section 4 we detail our approach explaining the algorithm we developed to solve PP-MS-RCPSP/max-cal.Following this, Section 5 contains the experiments and benchmarks conducted to evaluate the performance of our solution.Lastly, our conclusions can be seen in Section 6.

Multi-skill/mode, preemption, calendars and generalized precedence in the literature
The multi-mode RCPSP (MM-RCPSP) allows reaching a compromise between resource usage and task duration as it frequently occurs in practice.It also includes non renewable resources to model budget constraints that may prevent from using the fastest modes for all tasks.Many solution methods have been proposed over the years from local search, to sat-based methods using hyperheuristics.[6,9,15,37].CP approaches are currently highly popular and efficient tools to solve RCPSP variants [18] and in particular the MM-RCPSP [36].
MS-RCPSPs are a particular MM-RCPSP based on the concept of skills enabling the compact representation of a possibly large set of modes corresponding to combination of resources managing a set of skills required by a task.There is a variety in the methods and models to solve the problem: tree search [29], genetic algorithms [20], mixed-integer linear programming [3,10,19,31].On the standard MS-RCPSP, we will compare our method with two state-of-the-art methods: the CP approach using no-good learning proposed in [38] and the greedy randomized adaptive search procedure (GRASP) presented in [25].
Preemption can be defined as the capability of stopping the execution of a task in the generated schedule, releasing the resources it was using in order to allow for the execution of a different one, and being able to be continued later.In some cases preemption makes the problem easier to solve when an NP-hard non preemptive scheduling problem has a polynomial preemptive counterpart [5] but in the case of job-shop or resource-constrained project scheduling the problem may become harder to solve due to a much larger search space [7,22].Only a handful of studies have combined multi-skilled resources and preemption in the same work.We can cite [12,24,25,26], where mixed-integer linear programming, constraint programming and metaheuristics methods were proposed.Prominent preemption can cause trouble to CP approaches: in [25] an experiment revealed that MILP obtained better results on a set of highly preemptive MS-RCPSP than the default search of IBM CP Optimizer.For preemptive MS-RCPSP, we will compare our method with the hybrid GRASP-LNS heuristics proposed in [26].For their application, the authors also considered partially preemptive tasks, where only part of the resources are released during preemption.Due to industrial relevance, we also include partial preemption in our model.
Resource calendars are often unavoidable characteristics in industrial contexts.They can be linked to a basic form of preemption in the sense that an on-going task has to be preempted when one of its resource becomes unavailable due to an off-time in the calendar.In [13], calendars were used in a MS-RCPSP context without including task preemption.Our model considers both situations and tasks can be preempted or not by calendars as frequently observed [17].
Generalized precedence constraints, also considered in our model specify more complex temporal relations between tasks than standard precedence constraints.They make the problem more complex in all RCPSP variants as even finding a feasible solution becomes NP-hard in the presence of generalized precedence constraints.Many works have considered such constraints in the standard RCPSP [4,32] and also for multi-mode preemptive RCPSP [27].The RCPSP/max-cal problem involving both generalized precedence and calendar constraints was again successfully solved by CP in [17].
A difference between the performance of generic solvers on academic RCPSP models and on their adaptation to domain problems with very specific requirements is observed in other works and we also could experience it in our benchmark experiments: usually the more special characteristics the problem has, the worse is the performance of generic solvers.Our model falls into this category as it includes all the above-mentioned complicating characteristics.The LNS approach appears as a technique of choice to find the right compromise between genericity and performance.

PP-MS-MM-RCPSP/max-cal: definition and formulation
The PP-MS-MM-RCPSP/max-cal generalizes the partially preemptive MS-RCPSP introduced in [26].An additional feature is the possibility to execute the task in different modes, in the same way it is done in multi-mode RCPSP [15].Generalized precedence constraints and k is a natural number representing the resource demand of activity i for resource k, under mode m i ; 13. ∀i ∈ A, l ∈ L, s i,mi,l is a natural number representing the skill requirement of activity i under mode m i ; 14. P is the set of precedence constraints A × A specifying which activity should precede another one; 15.P sync−start is the set of activity pairs that must start at the same time; 16.P startlag is the set of ordered pairs (i, j) specifying ∆s (i,j) , the minimum time lag between start of i and start of j; 17.P sync−end is the set of activity pairs that must end at the same time; 18. P endlag is the set of ordered pairs (i, j) specifying ∆e (i,j) , the minimum time lag between end of i and start of j; 19.The set A is split in three different subsets : A = A P ∪ A N P ∪ A P P ; A P is the set of fully preemptive activities (activities can be preempted and all the resources are released); A N P is the set of non-preemptive activities; A P P is the set of partially preemptive activities (where at least one resource is not releasable); under mode m i .

Known variants in the literature
Our generic formulation encompasses more classical scheduling problems, that can be therefore solved using the proposed solution, such as:

Constraint Programming formulation
We developed a combinatorial optimisation library containing the CP model for PP-MS-MM-RCPSP/max-cal.It can be found as part of an open source framework to solve discrete optimization problems 1 .
A big focus of the library is on the RCPSP's class of problems described in this paper.All solver approaches, from Local search, rule based heuristics, CP, LP and LNS methods are either coded using the discrete-optimisation library or wrapped into it, allowing the user to easily benchmark methods and hybridize different approaches.
The following sections detail the decision variables and constraints of our formulation for PP-MS-MM-RCPSP/max-cal.

Decision variables
We assume that each task can be preempted at regular discretized time points, resulting in a sequence of small non-preemptive chunks for each activity.A subpart of an activity is thus defined as a subset of adjacent chunks of this activity.Let's define max preemption as an arbitrary input of our problem which represents an upper bound on the number of preemption breaking times allowed for all our activities.We will note J = {1, .., max preemption } the set of preemption breaking time indexes and J − = {2, .., max preemption } 1. Starting time decision : starts is a |A| × max preemption matrix, which contains the starting time of subparts of all the tasks in increasing order.2. Duration decision : durations is a |A|×max preemption matrix, which contains the duration of subparts of all the tasks.3. Mode decision : modes is a |A| vector specifying in which mode a task is executed.

Resource allocation : allocation is a |O| × |A| × max preemption 3D binary matrix which
indicates which worker o is allocated to each subpart of an activity.

Constraints
We will use the notation [:] in matrix indexing to ease the readability of constraints.For e.g starts[i, :] is the vector of starting times for the task i and allocations[o, :, :] is the |A| × max preemption matrix allocation of resource o to all activity subparts.1. Resource consumption constraint : where b modes,k is an array of dimension |A| × max preemption defined by: ∀i ] the resource demand of activity i on the j-th subpart of its execution, which doesn't depend on j in our problem.We use the global cumulative constraint implemented in most CP language to model the cumulative resource consumption in RCPSPs [30].
To take into account variable resource availability, we use a discrete stepped function max t (B kt ) − B ktstep at each t step that includes an artificial fixed task consuming it.This way, whenever availability changes, that task removes the resources from the pool during a period defined by the step function.In our domain, variable resource availability translates into worker's shifts and calendars, also indirectly defining points where a task needs to be preempted.For simplicity we didn't include it in the description of the constraint.

Non-renewable resources
This non-renewable resource constraint is only meaningful when there are several possible modes.In multi-mode settings, some modes value wouldn't satisfy the constraint.

Skills requirement:
∀i Contrary to other multi-skill variants, we consider that Operators assigned to a task contribute to the task skill requirements with all of their skills, similar to what it is seen in [26].

Operator availability:
∀o ∈ O, cumulative(starts, durations, allocation[o, :, :], Ao:) Variable availability of an Operator is dealt with in the same way as Constraint 1 above.

Precedence relation:
Conventions constraints for starts and durations: The following constraints are modeling choices aiming to help the solver to explore the search space.
As soon as there is a 0 duration subtask, all the following ones are 0 too.
: Whenever there is a 0 duration, the remaining starts values are uniquely determined by previous values, pruning redundant solutions.

Partially preemptive and non releasable resources :
We introduce variable blockedduration[i, j, k], ∀i ∈ A, j ∈ J , k ∈ R ρ as the usage duration of resource k for task i and activity subpart j.It is regulated by the following constraint: ∀j ∈ J , This means if the resource is releasable then the duration of the usage is the same as the duration of the subtask.
On the other hand, to describe the consumption of a resource over the total span of the activity for not releasable resources, we have: This will set blockedduration[i, 1, k] to the total span of the activity, including preempted time, and ignore the other subparts.Then, when indicating the cumulative constraint, instead of the normal duration, this resource uses blockedduration: We note that Constraint (13) implementing preemption won't always lead to a feasible solution without backtracking even if there are no time windows or maximum time lags when using CP-Optimizer solver.
Let us take a simple example of a single preemptive activity A 1 of duration 5 task decomposed in 3 subtasks A 1,1 , A 1,2 and A 1,3 .Suppose the solver sets the start time of A 1,1 to 0 and the end time of A 1,1 to 3 (and consequently its duration to 3).Suppose now that the solver sets the start time and the end time of A 1,2 to 3 (duration 0).Then according to constraint (13b) the duration of A 1,3 is set to 0 and a failure occurs due to the impossibility to satisfy constraint (10b).
To deal with this issue, we took advantage of the expressiveness of CP-Optimizer to add an alternative variant in our model: using the concept of interval instead of using variables starts and durations.An interval is defined as follows: ∀i ∈ A, interval[i] is a |max preemption | array of optional intervals variable.Each interval has the attributes present (boolean), duration, start, end (all integers).We also create one unique interval for each task, called spantask.New constraints are applied to this interval variable : [i, j].duration ≥ 1 : We don't consider 0 duration for task intervals, which avoids the above-presented issue.This could also help the solver to remove unnecessary solutions in its search space since otherwise there could be multiple equivalent solutions.
.present = F alse: When one sub-interval is not present, the remaining ones are not present either.This is equivalent to 13b constraint.

Objective function
We will focus exclusively on the makespan, which in terms of our formulation is equal to the ending time of the sink task n.Hence the goal is to minimize makespan = starts[n, max preemption ] + durations[n, max preemption ] However as detailed in section 4.3, due to the difficulty to obtain feasible solution, in our solution approach we opt in a first phase for a relaxation of a set of constraints and a penalization of their violation in the objective function.Thus the actual objective in this phase is a lexicographic optimization of the violation penalty and the makespan (see section 4.3).

A small PP-MS-MM-RCPSP/max-cal instance and its optimal solution
We provide an example PP-MS-MM-RCPSP/max-cal instance to illustrate its output from a simple problem definition contained in Table 1.It contains 6 activities to schedule (including source and sink tasks); it also includes multi-mode tasks, variable operator availability, release and deadlines, complete and partial preemption and finally one synchronisation constraint.This instance can be found in the toy model folder of our open source model's repository 2 .The optimal solution is depicted in Figure 1.We can observe that the tasks A 1 and A 3 are partially preemptive activities : the resource R 1 is therefore still used even though the task is paused.The calendar constraint of Operator 1 is visible in the corresponding Gantt chart where we see no operations assigned to it during its break time.The mode allocation is the following : Table 1 An instance of PP-MS-MM-RCPSP/max-cal.

Generic large neighborhood search algorithm
Our approach to solve the PP-MS-MM-RCPSP/max-cal is to use a generalisation of Large Neighborhood Search (LNS) for scheduling problems [21].In LNS, We iteratively improve the solution quality of a Master Problem MP by solving at each step a Reduced Master Problem RMP, based on the MP and previously found solutions.The way RMP are built are the core of LNS methods.The generic LNS algorithm is described in Algorithm 1, where X iter denotes the solution at iteration iter and Y iter its objective value.
Algorithm 1 Generic Large Neighborhood Search Algorithm.

Initial solution provider
We rely on a generalisation of the serial schedule generation scheme (SGS) procedure [16] to produce an initial solution for PP-MS-MM-RCPSP/max-cal (Algorithm 1, line 2).A similar generalisation was implemented in [26] being the closest one we found in the literature since it includes multi-skill, preemption and partially releasable resources.
Our SGS implementation takes as input a priority ordering of activities order_act, given as a permutation of A and, for each activity i ∈ A, another priority ordering order_res i ∈ S |O| of the workers o ∈ O to be assigned to the activity and the mode id modes i chosen to run the activity, as defined for the CP model.

C P 2 0 2 3
Activities are scheduled incrementally based on their priority as in classic serial SGS, at their earliest possible start date considering resource availability.Worker allocation is also done greedily using the order_res array.The preemptivity feature of our problem allows us to schedule subpart of some activities.Our open source library includes an implementation module handling all the features of PP-MS-MM-RCPSP/max-cal except for the P sync−start/end constraints, release and deadline constraints.The greedy procedure doesn't ensure that the constraints are feasible for these hard constraints.
Using the SGS procedure makes possible to run simple local search algorithms to get one initial solution of the problem.Such solver explores the vectored state space order_act, order_res, modes, and evaluate its fitness using SGS.The output evaluation of SGS is the same than the LNS or CP ones : makespan + potential penalty when deadlines and other constraints are not fulfilled (see details in Section 4.3).The state space being a permutation one, we are using moves chosen randomly from a portfolio of potential moves : partial and total shuffling, swap 2 random activities, 2-opt moves.Only one specific mutation has been also implemented to correct potential deadlines constraints (that put first the late task).In this paper we are using simulated annealing (SA) [14] as a local search solver.SA was parameterized with an initial temperature of T = 3 and exponential cooling schedule of factor α = 0.999.A restart strategy also ensures to roll back to the current best solution when we haven't improved the quality after N = 300 iterations.

Subproblem building
To build the reduced master problem, we rely on the constraint programming formulation and include additional constraints, which makes it simpler to solve.We decompose the subproblem building procedure (Algorithm 1, line 5) in two main methods.The first one, activityselector will split the set of activities A into 2 disjunctive subsets A sub and A sub .Different methods have been implemented and tested to select A sub and A sub as described below: A sub = sample(A, ⌊f • |A|⌋), and A sub = A \ A sub 2. Random selection and neighbors : consist in the previous selection, then add the predecessors of each task in A sub .Formally, given f ∈ R, 0 ≤ f ≤ 1: Cut in equal parts : given c ∈ N, we sort the activities by increasing order of starting times in the current best solution X * , then cut this ordered list in c equal pieces.When called, the function returns one of these c parts.A good strategy we observed consists in returning them in increasing order.4. Generalized precedence adapted selection : when considering the constraints of deadlines, release dates, and the generalized precedence constraints (P sync−start , P startlag , P sync−end , P endlag ) , we can reach a solution X * violating some of the constraints.Therefore we add in priority the tasks responsible for constraints violation in A sub and their predecessors in the precedence graph.

Mixing methods : given a pool containing the previous four methods, it will return A sub
and A sub of one of the pool members (either randomly at each iteration, or based greedily on the effectiveness of the selected method based on previous iterations).
Once we have the activityselector function, its output (A sub and A sub ) can be used in the constraintbuilder function to add constraints to the CP model as described in Algorithm 2. Instead of totally locking the variables corresponding to A sub and keeping free for optimisation the variables of A sub , we consider different ϵ values to constrain start and duration variables.We typically assign high values to parameters ϵ sub, * and small values to ϵ sub, * .

Subproblem solving
To solve the subproblem we use the open source lazy clause generation solver Chuffed3 or alternatively IBM's CP-Optimizer4 backend.Both rely on a CP formulation done in Minizinc language, and worked well in our experiments.At each iteration of the solver, we set an upper time computation h.In case we want to optimize The objective function can be the final activity n's end or the end of the subset of activities A sub we consider the most important in the RMP.Both strategies impact slightly the convergence of the overall algorithm but we didn't assessed the details in our current experiments.
After retrieving the solution of the RMP from the solver, it is post-processed by a left shifting procedure that compresses the full schedule.It is particularly useful in the case where the constraints introduced in constraintbuilder create idle times in the resulting schedule.

Relaxing constraints
The current SGS implementation does not return a feasible solution when certain constraints of PP-MS-MM-RCPSP/max-cal are present.Namely, those are the deadline and generalized precedence constraints P sync−start , P sync−end , P startlag , P endlag .By relaxing constraints ( 6), ( 7), ( 8), ( 9), ( 12) from Subsection 3.2.2we include a violation penalty into the objective function.This violation penalty for P sync−end is illustrated in Figure 2. The objective function that LNS and underlying CP solver minimizes is makespan + M • penalty where M is a large number.
This means that in practice, there will be 2 phases of optimisation : Feasibility phase : In this phase the violation penalty decreases, converging to 0. Optimisation phase : the algorithm optimises the original objective function (makespan) A1 A2 ∆t time 5 Experimental results

Manufacturing domain instances
We ran a series of experiments in order to evaluate our LNS method for PP-MS-MM-RCPSP/max-cal.All experiments have been done in a MacBook Pro with 2,7 GHz Quad-Core Intel Core i7.In the first experiment, our goal is to assess the performance of our method in a practical situation, so we extracted data samples from a real manufacturing use case.The experiment relies on two scheduling problem base instances (A and B) obtained from one of the assembly plants at Airbus, and are described in Table 2. Using these 2 base instances, we built 32 different instances, using all the possible combinations of the following features : Calendar: we can use our definition for resource consumption (Constraint 1 in section 3.2) or consider complete resource availability at any moment, hence making vectors Temporal Constraints: we include or not the deadline times, release times (described in Section 3.2 as the constraints number 11 and 12) and generalized precedence constraints (constraints number 6 to 9) Preemptive : tasks can be preempted if needed or preemption is forbidden; in the non preemptive use case, A N P = A Multiskill : skilled workers are present or not; in the non multiskill use case, workers belong to R and O = ∅ We ran our LNS algorithm with a time limit of 1800 seconds for all of the instances.As a comparison, we also ran the simulated annealing algorithm (SA) and a direct CP solving of Model CP-Base using Minizinc and the CP-Optimizer backend, the Chuffed backend and the CP-SmartPreemption model for the same amount of time.
The results for each instance are detailed in Table 3.Four of those instances (number 8 to 11) are infeasible by design, because they include the calendar constraints but not the preemption ones, so we left out those and will focus on the remaining 28 instances for our experiments.For the sake of clarity, we only show CP-SmartPreemption results (CP-S in the table) as it it consistenly better than the CP-Base model whether the latter is solved by Chuffed or CP-Optimizer.A comparison between CP-SmartPreemption and CP-Base can be found in Appendix A.1.
Table 3 Makespan results obtained for the 28 testing instances (splitted in two tables) using the three solvers: LNS, CP-SmartPreemption (abreviated as CP-S) and SA.For SA penalty values are included when the solution is not satisfying all the constraints (solutions for LNS and CP-SmartPreemption got no penalty).The table also includes each instance base problem (A or B) and its features: M for Multi skill, G for Generalized precedence constraints, P for (partial) Preemption and C for resource Calendars.In bold best results among the three algorithms.As we can see, LNS reaches the best results on 15 out of the 28 instances, CP-SmartPreemption on 16 and SA on 13.By comparing the different instances where each method excels, we observe that instances with high complexity in the type of constraints, i.e a combination of having generalized precedence constraints, calendar, multiskill or preemptive tasks are the ones for which we would rely the most on LNS: instance 15 and 31 are indeed the ones including all the available features.For those instances, CP-SmartPreemption is not returning any result in the allowed time whereas SA fails to find feasible solutions.

ID
On the other hand, SA is generally competitive as long as there is no generalized precedence constraint involved.But if those are present it finds difficulty in yielding a solution, which was expected since the SGS does not fulfill all the generalized precedence constraints.However, in some of the instances SA got better results with some of those features present, thus we theorize that the solution space heavily impacts in this method, helping to overcome the presence of complex constraints and obtaining nonetheless an optimal solution: if we do not take into account instances where no feasible solution was found, SA achieves best results in 13 of 20 instances.
Figures 3a and 3b show makespan and violation penalty from the solutions found using our method for the most complex two variants (containing all the possible features).We can confirm that our strategy of relaxing certain constraints at the cost of adding a penalty is useful to get an initial solution that will be improved from that point, as described in Subsection 4.3.

Experiments on instances with multiskill and partial preemption
We conducted a second experiment using instances of partial preemptive multi-skill problems based on the schedule of operations of a nuclear facility [26].The benchmark is divided in 4 different sets, each having different distribution combinations of A P , A N P and A P P .The multi-skill variant of this benchmark can be directly mapped into our formalisation, except for the minimum number of operators required by activity constraints.To adapt our model we introduced in the problem instances a dummy skill s all mastered by all operators.We then set skill requirement of skill s all for each activity, to the minimum number of operators required.
The comparison results, detailed in Table 4, are based on the average gap to lower bound, like it was done in the original paper.We compared three of our solving methods to the one tested previously.The first one is obtained by running simulated annealing for 300 seconds, working on the permutation of tasks and the SGS algorithm.In the second we ran simulated annealing alone for 240 seconds as second baseline, and finally we ran LN S during 60 seconds.The SGS method combined with SA is competitive with the best dedicated algorithm GRASP+LNS [26].LNS only improves slightly SA results as it can be seen when comparing SA (240 s) and SA (240 s) + LNS (60 s), which may be explained by the performance of SA itself, which is close or even better than the baseline.

Experimental results of LNS on standard multi-skill instances from the literature
Additionally, we did another benchmark comparison using the instances contained in [38].
The goal of running this second experiment is to have again comparative results to the work of [26], which our method can be related with, since it also aims at handling multi-skill and partial preemption in scheduling problems.The mentioned instances are divided in five sets, but we decided only to tests in one of them, set 1B, since it is the set with the lowest solution rates (only 12.5% solved) using the reference solvers.
Set 1B was proposed by [2] and it contains instances with 42 activities, 4 skills and between 20 to 60 resources.For these experiments, we relied on the CP model used in [38] 5 , which models classical multi-skill where operator perform at most 1 skill on a task.We used this CP model in our generic Large Neighborhood Search algorithm to evaluate how well the approach generalizes to different scheduling problem variants.Baseline solutions we compare with are the ones found by running Chuffed solver with the best search strategy.They can also be found in the repository.
The current implementation of SGS prevents us to test simulated annealing as a solver or initial solution provider, because of the constraint stating that an allocated worker can at most perform 1 skill on a task.Therefore, the initial solution used in LN S is computed by using the Chuffed solver directly on the minizinc model from [38] for 30 seconds.Then, if the solution is not optimal, the large neighborhood search is used with the mixing neighborhood strategy for a maximum time of 500 s.Worker allocation is fixed only in 10% of the tasks, which empirically showed good performance in terms of improvement rate during the LN S algorithm.We are improving 151 out of 216 of the Set 1B instances (69.9% of the total) whereas Young et al. results [38] are still better on 29 instances (13.4%), and equal on 36 instances (17.7%).On the 151 improved instances, the average improvement of our solution over the baseline is 4.61%.On the 29 that Young et al.CP method [38] was better, the makespan was worsen in average by 3.46%.In total average, our method gave a 2.76% improvement on the makespan.
Finally, we conducted some preliminary experiments on the new MSLIB [34] benchmark.The authors provided us some baseline best known solutions found by a genetic algorithm (GA) approach [33] on the hardest set of instances MSLIB4.Precisely we ran experiments on the first 404 instances of MSLIB4 using CP-Chuffed 30 seconds as initial solution followed by LN S for 100s.Our computation time is an order of magnitude higher than the GA that has a CPU time of around 15s.
The results indicate a strong impact of the SS (skill strength) parameter of the instances.This ratio roughly gives the scarceness of the skills in the sense that SS = 0 means that the number of resources that master a skill is equal to the maximum demand in operators mastering this skills over the activities, i.e. the minimal value that ensure feasibility.LNS systematically outperforms GA when SS > 0 (90 instances) with an a average improvement of 2.9% but when SS = 0 (164 instances), LNS is largely outperformed by GA with an average 19.2% gap.Additional material on these results are presented in Appendix A.3.

Implementation
We provide an open repository containing scripts to rerun the benchmarks presented in Sections 5.2 and 5.3 to ease replication for further research 6 .Every multiskill variant presented in the paper is using the same python object placeholder with different attribute values.About the algorithm most of it is reused by all variants : 1.The SGS procedure is the same whatever for all benchmark and therefore the metaheuristics methods that only rely on the SGS procedure are reused identically 2. LNS algorithm 1 is the same for all benchmarks 3. Activity selectors mixing method described in 4.2 are used on all benchmarks.Experiments that justify the choice of the mixing method are presented in Appendix A.2.
4. Specific : We are using 2 different minizinc modeling in our benchmark, one containing the most property of PP-MS-MM-RCPSP/max-cal that is used in benchmark 5.1 and 5.2 and one another model for classical MSRCPSP for benchmark 5.3.

Conclusions
In this paper we have presented a large neighborhood search (LNS) method to solve partially preemptive multi-skill/mode resource-constrained project scheduling problem with generalized precedence relations and resource calendars (PP-MS-MM-RCPSP/ max-cal).Solvers including all of those features are scarce, but are usually needed to approach real manufacturing or assembly environments.
In order to validate our method, we performed three experiments to benchmark it against constraint programming and simulated annealing.In the first experiment we used data from an Airbus use case, containing a real world scenario with 32 different variations.To make a fair comparison, we designed the CP-SmartPreemption model that handles preemption in a clever way significantly outperforming the our baseline CP model CP-Base.Despite its good performances, CP-SmartPreemption is unable to find a feasible solution in 6 industrial instances, while the LNS method finds a solution on all instances.From the observations of the samples, we conjecture that our method works better when (almost) all complex features are present in the problem definition (i.e.: calendar, generalized precedence constraints, multiskill and preemption).Then, we tested our method using a benchmark from research work on partially preemptive multi-skill scheduling [26] to evaluate its performance in another real scenario.In this experiment, the simulated annealing methods appears competitive, slightly improved by LNS.Our method was the second best among the tested ones given the benchmark conditions, still being competitive enough to be close to the best one of that benchmark.In the last experiment on standard MS-RCPSP instances, we used our method as a mean of improvement for the solutions given by the CP based method.On the instances from [38] we improved the best known solutions in 69,9% of the instances.On the instances from [34] the performance of our approach is highly correlated with skill scarceness, which open a path for future research.
In conclusion, we found our LNS based method, available in a new discrete optimization open-source library, appropriate to solve scheduling problems containing combinations of complex features like the ones found in industrial instances, and is still reliable to be used for more academic problems.
An interesting perspective is to be able to reuse the subproblem solving information from one iteration the other to save computational time.Current Minizinc implementation has a limitation w.r.t.incrementality.Nevertheless, it would be worth to investigate how solution-based phase saving approaches [8] use nogoods and see if it is applicable in our LNS framework for future work.Adapting propagation guided LNS [23] to our framework is also a promising research direction.

18 .
∀i ∈ A, j∈J (interval[i, j].duration * interval[i, j].present) = p i,modes[i]The sum of duration of present intervals should sum to the duration of the task.19.∀i ∈ A, span(spantask[i], interval[i, :]):We use the native constraint span of CPOptimizer so that the spantask[i] spans over all present intervals in interval[i], ignoring the non-present.The precedence constraints are written using the spantask interval.We call this new model specific to CPOptimizer, Model CP-SmartPreemption while Model (1-14) is called CP-Base.

Figure 1
Figure 1 Operator and resource oriented Gantt chart for the example instance.

Figure 2
Figure 2 ∆t is added to the objective function whenever (A1, A2) ∈ P sync−end .

Figure 3
Figure 3 Evolution of solutions found for LNS algorithm in two of the instances.Makespan shows in green at the left Y axis and penalty in red at the right Y axis for both of the instances.

Table 2
Base instance description of the manufacturing use case.

Table 4
Average gap to lower bound for different algorithms and by set of instances.Best results for each set are bold and second best are underlined.