1 Introduction

In this paper, a multi-skill resource-constrained multi-project scheduling problem (MSRCMPSP) is modeled and solved. The problem is motivated by a real industrial issue at SNCF, the French national railway company, which carries out the heavy maintenance of its rolling stock in ten different factories. Several rolling stock units are maintained simultaneously, and each unit is considered as a project. To complete each project, a certain number of activities requiring multiple skilled resources must be performed. Different types of resources (maintenance operators and machines) with different characteristics and constraints are taken into account. Two original objective functions, that have rarely been addressed in the literature, are independently considered: (i) Minimization of the sum of the weighted tardiness of the projects (SWTP) and (ii) Minimization of the sum of the weighted duration of the projects (SWDP).

Project management has attracted the attention of many researchers over the years. Planning and scheduling are crucial for controlling the execution time of projects, managing the required resources, and meeting (customer) deadlines. Without a detailed efficient plan, activities would be poorly executed causing resource conflicts, waste of time and thus cost increase. In a world where market competition and customer expectations are very high, companies want to digitize their processes to better control production, optimize resource allocation and reduce costs (Moeuf et al. 2018). However, project scheduling is a very complex problem. The well-known Resource-Constrained Project Scheduling Problem (RCPSP) belongs to the class of hard optimization problems, and instances of more than 60 activities cannot be solved with exact methods in reasonable computational time (Koné et al. 2011). Furthermore, most companies produce several products which generally share the same resources. The survey of Lova et al. (2000) found that 84% of companies work with several projects. The multi-project version of the RCPSP is closer to real-word applications. Yet, most papers study the single project version of the RCPSP. The objective function is often the minimization of the makespan, i.e., the completion time of the last activity. The multi-project extension results in more complex problems because of the scarcity of shared resources, the interactions and the competition among different projects, project specific characteristics, deadlines and more elaborate objective functions (Browning and Yassine 2010).

We first show how to model a complex industrial problem by proposing a time-indexed MILP (Mixed-Integer Linear Programming) model with several constraints such as precedence constraints with lag times, time-dependent resource capacity constraints and machine assignment. In the model, each rolling stock unit is considered as a project, maintenance operations as activities requiring a certain number of human resources with different skills and one specific machine to be executed. Since resources have multiple skills and several rolling stock units are maintained simultaneously, the problem studied in this paper corresponds to the Multi-Skill Resource-Constrained Multi-Project Scheduling Problem (MSRCMPSP) (Pritsker et al. 1969; Bellenguez and Néron 2004). The problem with the considered constraints and objective functions has never been addressed in the literature.

Secondly, several solution approaches are proposed and tested to find the most effective one. To address large industrial instances involving hundreds of projects and thousands of operations, two greedy algorithms are first introduced: a serial scheduling generation scheme and a parallel scheduling generation scheme with different priority rules. These priority rules are determined using a global precedence graph, which is constructed from the precedence graph of each individual project (which is a graph showing the critical paths for each project). The solutions produced by the greedy algorithms are then further optimized using a memetic algorithm. Finally, by conducting various computational experiments, the potential gains are quantified.

The main contributions of this paper are summarized below.

  • To our knowledge, this is the first attempt to tackle the integration of both multi-project and multi-skill versions of the RCPSP problem with an objective function other than the makespan. A time-indexed MILP formulation of the problem is provided.

  • A new Memetic Algorithm (MA) and a new Simulated Genetic Algorithm (hSGA) are proposed for scheduling multiple projects with different resource constraints. The proposed solution approaches are compared, using large industrial instances, to the time-indexed MILP model, and to popular solution approaches of the literature (two constructive heuristics, Simulated Annealing (SA) and a Genetic Algorithm (GA)).

  • The MA is evaluated on benchmark instances, and the computational results show that it stands as one of the most effective methods compared to the existing literature.

  • The gap between theory and practice is reduced by considering more realistic objective functions and constraints. As highlighted by Hartmann and Briskorn (2022), Sánchez et al. (2022) and Rahman et al. (2020) real-world case studies are necessary to motivate more complex models.

The paper is organized as follows. Section 2 gives an overview of the literature on project scheduling problems. Section 3 describes the industrial problem and introduces the associated mathematical model. The solution approaches are presented in Section 4. Section 5 provides the numerical results used to evaluate the performance of the proposed approaches. Finally, Section 6 concludes the paper and gives some perspectives on future work.

2 Literature review

In this section, we review the literature on the multi-project and multi-skill extensions of the RCPSP. Section 2.1 recalls the classical RCPSP and then focuses on the multi-project extension. Section 2.2 reviews the research on project scheduling problems with multi-skilled resources. In Section 2.3, papers integrating both multiple projects and multiple skills are discussed. The limits of the existing literature are also discussed.

2.1 Scheduling problems with multiple projects

The resource-constrained project scheduling problem (RCPSP) is a complex optimization problem that involves scheduling a set of activities subject to precedence constraints and resource availability (Deblaere et al. 2011). It is an extension of the classical job shop scheduling problem and is NP-hard in nature (Blazewicz et al. 1983). Since its formulation by Pritsker et al. (1969), many extensions and solution methods have been proposed. A comprehensive overview of the different RCPSP problems and variants can be found in the literature such as Özdamar and Ulusoy (1995), Brucker et al. (1999), and more recently Habibi et al. (2018) and Hartmann and Briskorn (2022). Kolisch (1996) and Kolisch and Hartmann (2006) examine the performance of various heuristics and priority rules for the classical RCPSP. Lancaster and Ozbayrak (2007) focus on evolutionary algorithms, while Pellerin et al. (2020) provide a comprehensive review on recent hybrid metaheuristics.

A generalization of the RCPSP is to consider the resource-constrained multi-project scheduling problem (RCMPSP). This extension is interesting since it is a step forward on modeling real world problems (Lova et al. 2000). Although the same methods for modeling and solving the RCPSP can be applied to the RCMPSP (Drexl 1991), developing efficient algorithms is more challenging. Dealing with several projects simultaneously significantly increases the size of the problems to solve. The deadlines are relatively easier or harder to meet depending on the tightness of each project, the delay penalties, etc. Additional objective functions, such as the sum of the weighted tardiness of projects (Krüger and Scholl 2009) can be considered and thus exploring the problem structure knowledge may increase the efficiency of solution techniques. In fact, for the multi-project version of the RCPSP, there are two approaches to present the links between activities (Kurtulus and Davis 1982):

  1. 1.

    The single-project approach, that uses two dummy activities and precedence arcs to combine the projects into a single global project. The problem is then reduced from the RCMPSP to the RCPSP and the critical paths of the projects are lost.

  2. 2.

    The multi-project approach, that uses \((P+1)*2\) dummy activities (where P is the number of projects) and where each project has its own critical path(s). Given the objective functions considered in this paper, the multi-project approach is more appropriate and is used to compute priority rules.

Even if most papers on project scheduling focus on solution methods for the single-project version, there exist some work on the multi-project version. Pritsker et al. (1969) are the first to propose a zero-one programming approach for the RCMPSP. Later, Deckro et al. (1991) explore the model of Pritsker et al. (1969) and use a project decomposition approach to solve larger multi-project problems. Similarly, Vercellis (1994) consider a Lagrangian decomposition technique to solve a multi-project planning problem.

Because of the complexity of the RCMPSP, heuristics based on priority rules are widely studied in the literature. Kurtulus and Davis (1982) experiment six new priority rules (PRs). The authors show that priority rules computed using the multi-path method (critical paths for each project) outperform priority rules computed on a single global graph. Browning and Yassine (2010) analyze the performances of 20 PRs and consider various objective functions. The authors help project managers by characterizing the best priority rule based on four problem structure measures: Objective function, network complexity, resource distribution, and resource contention. Lova et al. (2000) use an iterative Forward-Backward heuristic to solve the RCMPSP and consider two time criteria (mean project delay and sum of duration of projects), and four non-time criteria (project splitting, in-process inventory, resource leveling and idle resources). Gonçalves et al. (2008) propose a random key genetic algorithm to solve the RCMPSP. The genes decode the priority of the activities (computed using slack times), the delay of each iteration g (when a new activity is scheduled) and the release date of each project.

The Multi-Mode Resource-Constrained Multi-Project Scheduling Problem (MRCMPSP) is an extension of the multiple project scheduling problem that has attracted the attention of many researchers (Chen et al. 2022). In the MISTA 2013 challenge, various solution methods were proposed to minimize, in lexicographical order, the makespan and the total project delay. For details on the presented methods, the reader can refer to Wauters et al. (2016). Other papers tackle the distributed (or decentralized) resource-constrained multi-project scheduling problem (DRCMPSP) (Confessore et al. 2007). Multi-agent based approaches are the most popular solution methods (Li et al. 2021).

2.2 Multiple skill resource allocation

The multi-skill resource-constrained project scheduling problem (MSRCPSP) is formalized in Bellenguez and Néron (2004). Each activity requires a certain number of (human) resources mastering a given skill. Since resources have multiple skills, not only which resources allocated to each activity must be decided, but the allocated resources must also have the required skills to execute the activity. Generally, skills can be classified into two types: Categorical and hierarchical skills (Snauwaert and Vanhoucke 2023). The categorical class makes no distinction between resources; they either possess the skill or do not. The hierarchical skill class offers additional information about resource efficiency. Resources with higher hierarchical skills can, for instance, process tasks faster. Moreover, certain activities can only be executed by a resource possessing a skill level that meets a specified minimum requirement (Bellenguez and Néron 2004).

These additional decisions and constraints make the problem even harder to solve than the classical RCPSP (Polo-Mejía et al. 2021; Almeida et al. 2019) and exact methods are only considered for small instances. Among the exact approaches that can be found in the literature, Correia and Saldanha-da Gama (2014) propose a MILP model with different valid inequalities. The considered objective functions are the total costs and the makespan. Li and Womer (2009) propose a hybrid MILP/CP Benders decomposition approach to solve the MSRCPSP where activities require only one multi-skill resource. Bellenguez-Morineau and Néron (2007) propose a branch and bound method with instances having at most 32 activities, 10 resources and 5 skills. Montoya et al. (2014) propose a branch and price algorithm with an activity and time decomposition approach and the makespan as objective function.

As for the RCPSP, heuristics based on different priority rules and metaheuristics are also popular approaches to solve the MSRCPSP. Myszkowski et al. (2015) compare state-of-the art priority rules using data from a real world problem. They conclude that complex priority rules are not necessarily better than simple ones. Almeida et al. (2016) propose a parallel scheduling generation scheme (PSGS) with activity priority rules and resource weights to avoid a random resource selection. At each stage of the PSGS, a flow graph is implemented to assign resources. Javanmard et al. (2017) develop a genetic-based algorithm and a particle-swarm-based algorithm for the MSRCPSP with preemptive activities. Lin et al. (2020) implement a genetic programming algorithm (GP) used as a high-level strategy to select a sequence of 10 low-level heuristics (priority rule based heuristics).

For more details about multi-skill resource allocation problems, the reader is referred to the review of Afshar-Nadjafi (2021). The author reviews 160 articles published from 2000 to 2020 and classify the articles based on the objective functions, the mathematical formulations and the solution approaches.

2.3 Multiple skills and multiple projects scheduling problems

In the literature, few papers integrate multiple skills and multiple projects. Cui et al. (2021) consider the multi-mode and multi-skill resource-constrained multi-project scheduling problem in high-end equipment production. A variable neighborhood search metaheuristic is developed to solve the problem. As the objective function is the makespan, the problem can be reduced to a single project with additional project-based precedence constraints.

Some papers consider multiple project scheduling with multi-skilled resources but the only sequencing decisions that are made are the starting order of projects (Heimerl and Kolisch 2010; Felberbauer et al. 2019; Chen et al. 2022). Similarly, Haroune et al. (2022) study a multi-project scheduling and multi-skilled employee assignment problem with preemptive tasks. Each project is broken down into several tasks, but without explicit precedence constraints. A mixed-integer goal programming (MIGP) formulation to optimally solve small instances is implemented. A local search algorithm and a tabu search algorithm are developed to tackle large instances provided by an IT company.

To the best of our knowledge, problems with multiple projects, multi-skilled resources, explicit precedence constraints between activities and other objective functions than the makespan are not studied in the literature.

3 Problem modeling

In Sect. 3.1, the problem of heavy rolling stock maintenance is described in detail. Information about the activities and resources involved in the problem is also provided. In Sect. 3.2, the mathematical model for the MSRCMPSP is presented.

3.1 Problem description

Heavy maintenance refers to the renovation and/or modernization of trains as well as the repair of different components such as electronic cards, bogies, axles and rotors. The main challenges in this process include achieving economic savings and reducing the environmental impact of the railway industry, while also increasing passenger comfort and service quality. During heavy maintenance, the rolling stock is immobilized for several weeks. Since a rolling stock unit is very expensive, the primary objective of SNCF is to minimize the immobilization time of these units. This is highly desirable for clients as it enables the rolling stock units to be operational again as soon as possible. Additionally, by decreasing lead times, idle times for resources are minimized, maintenance costs are reduced and the maintenance process becomes more cost-effective.

The first objective, minimizing the immobilization time of rolling stock units, can be seen as a tactical objective. In fact, each year a plan is made for the next one or two years. The plan is communicated to the clients, and the due dates are fixed. In this way, the clients can build their timetable based on the availability of rolling stock units. But, at the operational level, there are many uncertainties and sometimes due to consecutive or large disturbances (e.g., new crashed rolling stock arrival) it is necessary to reschedule. In this case, respecting the due dates of customers, fixed at the tactical level, is the main objective for SNCF maintenance workshops since, to maintain the train timetable, the rolling stock units must be available on time.

When trains enter into the workshop, some technical tests and observations are performed to re-estimate the workload. Then, trains are uncoupled in several coaches. Components of coaches are dismantled and repaired in parallel. To reassemble the train, the activities of all coaches and components must be finished. Once coaches are coupled, the final activity consists of testing the train to check that safety and quality standards are respected. Figure 1 illustrates the activity precedence graph of a very simplified maintenance procedure of a rolling stock unit composed of only two coaches.

Fig. 1
figure 1

Simplified activity precedence graph of rolling stock heavy maintenance

Full size image
Table 1 Example of the capacities in working hours of a team with two skills \(k_1\) and \(k_2\)
Full size table

In the maintenance workshops, several rolling stock units, with different activities to process and deadlines to respect, are repaired simultaneously and compete for resources. Two kinds of renewable resources are considered: (i) Several teams having several maintenance operators with multiple categorical skills and daily variable capacities. An example of the characteristics of a team is shown in Table 1. The team has a different total capacity for each period and for each skill \(k_1\) and \(k_2\). The capacities of skills \(k_1\) and \(k_2\) are computed according to the availability and the skills of each operator. (ii) Locations in the maintenance center equipped with one or more installations (e.g. a garage pit and a roof access) but, for safety reasons, only one installation at a time can be occupied. A location can be seen as a renewable resource with multiple categorical skills (the skills being the installations at the location).

Fig. 2
figure 2

Example of activity resource requirements

Full size image

To process an activity, several resources are necessary. Each activity is composed of several sub-tasks with a fixed order. Each sub-task requires a workload of team \(r \in \mathcal {R}\) with a skill \(k \in \mathcal {K}\) (during the processing time of the sub-task and not the duration of the activity). To execute the set of sub-tasks, an installation is always necessary (e.g., a garage pit). An example of activity resource requirements is shown in Fig. 2:

  • At period \(t=1\), activity a requires from team \(r_1\) a total of workload of 60 h: 20 h of skill \(k_1\), and 40 h of skill \(k_2\). Human resources become available again at the completion of each sub-task.

  • At period \(t=2\), activity a requires a workload of 30 h of skill \(k_2\) from team \(r_1\) and a workload of 20 h of skill \(k_3\) from team \(r_2\).

  • An installation of type \(i_1\) which, contrary to the human resources (teams), is necessary during the entire processing time of activity a, i.e., the processing time \(p_a\).

Since a team has multiple operators and thus a given capacity for each period and skill, many sub-tasks can be processed in parallel (as it is the case for Task1 and Task2 in Fig. 2) as long as the capacities of the team are not exceeded. In fact, due to long-horizon scheduling (1 to 2 years) and the uncertainty on the number of operators, operators are not assigned. However, we ensure to not exceed a capacity threshold for each team \(r \in \mathcal {R}\) and skill \(k \in \mathcal {K}\) at each period \(t \in \{1,..., H\}\) (constraints (10) and (11)).

Furthermore, the workshop configuration is very unlikely to be modified on the considered scheduling horizon. Thus, we assign to each activity the location with the installation required to process this activity. Since the locations have several installations, this is a multi-skill resource assignment problem (Bellenguez and Néron (2004)).

Fig. 3
figure 3

Example of a solution of a problem with three activities

Full size image

In Fig. 3, a solution of an example with three activities, two machines \(m_1\) and \(m_2\) and one team \(r_1\), is illustrated. A precedence constraint, with a lag time \(l_{a_1,a_2}\) of one time unit, exists between activity \(a_1\) and activity \(a_2\). Machine \(m_1\) includes two installations \(i_1\) and \(i_2\), while machine \(m_2\) contains only installation \(i_1\) (Fig. 3b). The total capacity of \(r_1\) and the capacities of its skills \(k_1\) and \(k_2\) are given in Fig. 3d. To simplify the example, we assume that the capacities remain constant over time.

In this solution, activity \(a_1\) and activity \(a_2\) are assigned to machine \(m_1\) while activity \(a_2\) is assigned to machine \(m_2\) (Fig. 3a). Figure 3c shows the workload and the total capacity of team \(r_1\). The workload and the capacity of skill \(k_1\) of team \(r_1\) is presented in Fig. 3e. Note that activity \(a_3\) cannot start before \(t=2\), since the total capacity of team \(r_1\) would be violated.

The goal is to find a resource feasible solution that minimizes the sum of the weighted tardiness of the projects or the sum of their weighted duration. The weights are given by the planners of SNCF (e.g. the weight of a high-speed train is usually larger than that of a regional train). The proposed mathematical model is presented in the next section.

3.2 Mathematical model

The following notations are considered to model the MSRCMPSP.

Parameters of the model:

  • H, number of periods in the horizon,

  • \(\mathcal {E}\), set of N projects,

  • \(dr_e\), ready date of the project (engine) e,

  • \(dd_e\), due date of project e,

  • \(w_e\), weight of project e,

  • \(\mathcal {A}\), set of activities to schedule,

  • \(p_a\), processing time of activity a,

  • \(\mathcal {P}_a\), set of pairs of activities: \((a, a') \in \mathcal {P}_a\) means that activity \(a \in \mathcal {A}\) must be processed before \(a' \in \mathcal {A}\),

  • \(l_{a,a'}\), positive or negative lag time between \((a, a') \in \mathcal {P}_a \),

  • e(a), function that returns the project of activity a,

  • \(\mathcal {M}\), set of machines (locations in the workshop with one or several installations),

  • \(\mathcal {I}\), set of installations,

  • \(\mathcal {I}_{m}\), set of installations (skills) associated to machine m,

  • \(b_{a,i} = {\left\{ \begin{array}{ll} 1 &{} \text {if installation } i \in \mathcal {I} \text { is necessary to process activity}\,a, \\ 0 &{} \text {otherwise}, \end{array}\right. } \)

  • \(\mathcal {R}\), set of teams ,

  • \(\mathcal {K}\), set of skills,

  • \(W^{\overline{\mathcal {R}}}_{r,t}\), total capacity of team \(r \in \mathcal {R}\) at period \(t \in \{1,..., H\}\),

  • \(W^{\overline{\mathcal {K}}}_{r,k,t}\), total capacity of team \(r \in \mathcal {R}\) for skill \(k \in \mathcal {K}\) at period \(t \in \{1,..., H\}\),

  • \(\alpha _{r,a} = {\left\{ \begin{array}{ll} 1 &{} \text {if team } r \in \mathcal {R} \text { is necessary to process activity}\,a \in \mathcal {A}, \\ 0 &{} \text {otherwise}, \end{array}\right. } \)

  • \(\phi _{k,a}(l)\), workload of skill \(k \in \mathcal {K}\) necessary to process \( a \in \mathcal {A}\) at period \(l \in \{1,..., p_a\}\). Decision variables:

  • \(X_{m,a,t} = {\left\{ \begin{array}{ll} 1 &{} \text {if } a \in \mathcal {A} \text { start at } t\in \{1,..., H\} \text { and } m \in \mathcal {M} \text { is assigned to } a, \\ 0 &{} \text {otherwise}. \end{array}\right. } \)

  • \(Y_{m,i,a} = {\left\{ \begin{array}{ll} 1 &{} \text {if } m \in \mathcal {M} \text { is assigned to } a \in \mathcal {A}\,\text {with installation}\,i \in \mathcal {I}_m, \\ 0 &{} \text {otherwise}. \end{array}\right. } \) To keep a similar structure as the model proposed in Pritsker et al. (1969), but also to ease the understanding of the mathematical model, the following auxiliary variables are defined. Auxiliary variables:

  • \(S_{a,t} = {\left\{ \begin{array}{ll} 1 &{} \text {if activity } a \in \mathcal {A}\text { starts at}\,t \in \{1,..., H\}, \\ 0 &{} \text {otherwise}. \end{array}\right. } \) Let us note that \(S_{a,t}=\sum _{m \in \mathcal {M}} X_{m,a,t}, \hspace{0,5 cm} \forall a \in \mathcal {A},\forall t \in \{1,..., H\}\). Furthermore, the completion time and the tardiness of project \(e \in \mathcal {E}\) are respectively defined as:

    $$\begin{aligned} C_e= & {} \max _{a \in \mathcal {A}; e(a)=e}( \sum _{t=1}^{H}tS_{a,t} + p_a) \\ T_e= & {} \max (0, C_e-dd_{e}) \end{aligned}$$

Using the notations introduced above and based on Pritsker et al. (1969) and Bellenguez and Néron (2004), the following MILP model is proposed.

$$\begin{aligned} Minimize \sum _{e=1}^{N}w_eT_e \end{aligned}$$
(1)

or

$$\begin{aligned} Minimize\sum _{e=1}^{N}w_eC_e \end{aligned}$$
(2)

subject to,

$$\begin{aligned}{} & {} \sum _{t=1}^{dr_{e(a)}-1}S_{a,t}=0 \hspace{1 cm} \forall a \in \mathcal {A} \end{aligned}$$
(3)
$$\begin{aligned}{} & {} \quad \sum _{t=1}^{H}S_{a,t}=1 \hspace{1 cm} \forall a \in \mathcal {A} \end{aligned}$$
(4)
$$\begin{aligned}{} & {} \quad \sum _{t=1}^{H}tS_{a,t}+p_{a}+l_{a,a'}\le \sum _{t=1}^{H}tS_{a',t} \hspace{1 cm} \forall (a,a') \in \mathcal {P}_a \end{aligned}$$
(5)
$$\begin{aligned}{} & {} \quad \sum _{a \in \mathcal {A}}\sum _{t_1=max(1,t-p_a)}^{t}X_{m,a,t_1}\le 1 \hspace{1 cm} \forall m \in \mathcal {M}, \forall t \in \{1,..., H\} \end{aligned}$$
(6)
$$\begin{aligned}{} & {} \quad \sum _{m \in \mathcal {M}}X_{m,a,t}= S_{a,t} \hspace{1 cm} \forall t \in \{1,..., H\} , \forall a \in \mathcal {A} \end{aligned}$$
(7)
$$\begin{aligned}{} & {} \quad \sum _{t=1}^{H}X_{m,a,t}=\sum _{i \in \mathcal {I}_m}Y_{m,i,a}\hspace{1 cm} \forall m \in \mathcal {M}, \forall a \in \mathcal {A} \end{aligned}$$
(8)
$$\begin{aligned}{} & {} \quad \sum _{m \in \mathcal {M}}Y_{m,i,a}= b_{a,i} \hspace{1 cm} \forall i \in \mathcal {I}, \forall a \in \mathcal {A} \end{aligned}$$
(9)
$$\begin{aligned}{} & {} \quad \sum _{a\in \mathcal {A}}\sum _{t_1=max(1, t-p_{a})}^{t}\alpha _{r,a}\phi _{k,a}(t+1-t_1)S_{a,t_1}\le W^{\overline{\mathcal {K}}}_{r,k,t}\hspace{1 cm}\nonumber \\{} & {} \quad \forall r \in \mathcal {R}, \forall k \in \mathcal {K}, \forall t \in \{1,..., H\} \end{aligned}$$
(10)
$$\begin{aligned}{} & {} \quad \sum _{k\in \mathcal {K}}\sum _{a\in \mathcal {A}}\sum _{t_1=max(1,t-p_{a})}^{t}\alpha _{r,a}\phi _{k,a}(t+1-t_1)S_{a,t_1}\le W^{\overline{\mathcal {R}}}_{r,t}\hspace{1 cm} \nonumber \\{} & {} \quad \forall r \in \mathcal {R}, \forall t \in \{1,..., H\} \end{aligned}$$
(11)
$$\begin{aligned}{} & {} \quad S_{a,t}\in \{0,1\}\hspace{1 cm} \forall a \in \mathcal {A},\forall t \in \{1,..., H\} \end{aligned}$$
(12)
$$\begin{aligned}{} & {} \quad X_{m,a,t}\in \{0,1\}\hspace{1 cm} \forall m \in \mathcal {M},\forall a \in \mathcal {A},\forall t \in \{1,..., H\} \end{aligned}$$
(13)
$$\begin{aligned}{} & {} \quad Y_{m,i,a}\in \{0,1\}\hspace{1 cm} \forall m \in \mathcal {M},\forall i \in \mathcal {I},\forall a \in \mathcal {A} \end{aligned}$$
(14)

The objective functions (1) and (2), respectively minimize the weighted sum of the tardiness of the projects and the weighted sum of the duration of the projects. Constraints (3) ensure that activities cannot start before the project start date (note that \(dr_{e(a)}\) can be replaced with the earliest starting date of each activity \(a \in \mathcal {A}\)). Constraints (4) are the non-preemption constraints, meaning that each activity has one and only one possible start date. Constraints (5) are the precedence constraints with minimum lag times (mainly transportation times of coaches from one location to another). An original disaggregated approach to write precedence constraints was proposed by Christofides et al. (1987) and Almeida et al. (2019):

$$\begin{aligned} \sum _{l=1}^{t+p_a-1}lS_{a,l}+p_{a}+l_{a,a'}\le \sum _{l=1}^{H}lS_{a',l} \hspace{1 cm} \forall t \in \{1,..., H\}, \forall (a,a') \in \mathcal {P}_a \end{aligned}$$
(15)

The main advantage of the disaggregated approach is that it has a better linear relaxation and can be used to reinforce the model. Constraints (6) ensure the non-duplication of machines, i.e., that machine (location) \(m \in \mathcal {M}\) is not assigned to more than one activity in a single period \(t \in \{1,..., H\}\). Constraints (7) ensure the synchronization of variables \(X_{m,a,t}\) and auxiliary variables \(S_{a,t}\). Constraints (8) guarantee that if machine \(m \in \mathcal {M}\) is assigned to activity \(a \in \mathcal {A}\), then exactly one of its installations (skills) is used. Constraints (9) ensure that the right installation is assigned to each activity a. Constraints (10) and (11) are the time-varying cumulative capacity constraints. Constraints (10) ensure that the total capacity of each skill \(k \in \mathcal {K}\) of team \(r \in \mathcal {R}\) is not exceeded, and constraints (11) ensure that the total capacity of each team \(r \in \mathcal {R}\) is not exceeded (Table 1). Finally, constraints (12), (13) and (14) are the binary constraints of the decision variables. A similar model for the RCPSP with multiple skills (without constraints (10) and constraints (11)) is proposed in Almeida et al. (2019).

4 Solution approaches

Section 4.1 introduces two greedy algorithms and 7 priority rules. The general framework and the different components of the proposed memetic algorithm are presented in Section 4.2.

4.1 Greedy algorithms

In the literature, there are two methods to determine feasible schedules (also called constructive heuristics) for the RCPSP problem (Kolisch 1996): Serial Scheduling Generation Scheme (SSGS) and Parallel Scheduling Generation Scheme (PSGS). The SSGS performs activity-incrementation while the PSGS performs time-incrementation. In the SSGS method, at each stage g, an activity is selected among a set of eligible activities \(\mathcal {E}_g\) (activities where each predecessor has already been scheduled), and scheduled at its first precedence and resource feasible time. The selection of the activity at stage g is based on one or several priority rules. The algorithm ends when all activities are scheduled, and the solution is a list with the activities in their scheduled order and with their starting times. The procedure of the algorithm is formalized in Algorithm 1.

Algorithm 1
figure a

Greedy algorithm: SSGS + Priority rule

Full size image

In the PSGS method, at each stage g, the smallest completion time \(t_g\) among the active activities (\(\mathcal {A}_{Active}\)) is calculated. Then, the eligible activities at \(t_g\) in terms of precedence constraints are computed and sorted according to a priority rule. If not all eligible activities can be scheduled at \(t_g\) due to lack of resources, activities with the lowest priority are postponed to the next iteration (Algorithm 2). The algorithm ends when all activities are scheduled.

Algorithm 2
figure b

Greedy algorithm: PSGS + Priority rule

Full size image

Because of its relative ease of implementation, SSGS is more popular than PSGS, but the performances of each method seem to be problem related (Kolisch and Hartmann 2006). For more details about complexity, performances and priority rules used in SSGS and PSGS the reader is referred to Hartmann and Kolisch (2000) and Kolisch and Hartmann (2006). In this work, both SSGS and PSGS are adapted and implemented to solve the MSRCMPSP. The main differences compared to the classical SSGS and classical PSGS stand on the assignment of multi-skilled resources and the computation of priority rules since they depend on the tightness of each project. The procedure that allows the availability of resources to be checked so that an activity a may start at time t, is detailed in Algorithm 3. First, the resource availability for each skill \(k \in \mathcal {K}\) and resource \(r \in \mathcal {R}\) is checked for each period \(t' \in \{t,t+1 ,..., t+ p_a\)}. Then, the first available machine \(m \in \mathcal {M}\) is selected and assigned to activity a. Since machines have multiple skills, the selection of machines is based on the number of skills as in Almeida et al. (2019).

Seven priority rules are implemented:

  • EF: The activity with the minimal earliest finish time is selected,

  • ES: The activity with the minimal earliest starting time is selected,

  • LF: The activity with the minimal latest finish time is selected,

  • LS: The activity with the minimal latest start time is selected,

  • Rand: Activities are selected in a uniformly distributed random manner,

  • SA: The activity with the smallest duration is selected,

  • SST: The activity with the smallest slack time is selected.

The priority rules considered in this paper are among the most common ones in the literature. They are computed using a precedence graph for each project.

Algorithm 3
figure c

Implementation of CheckResources(a, t)

Full size image

4.2 Memetic algorithm

The solutions determined by the greedy algorithms in Section 4.1 are improved using a memetic algorithm (MA). A memetic algorithm combines a Genetic Algorithm (GA) and a local search procedure which is usually applied after mutation or as a mutation operator on the new individuals (Moscato and Cotta 2003; Gonçalves et al. 2005).

In this paper, a hybridization of simulated annealing and genetic algorithm (hSGA) is combined with a conditional simulated annealing (SA) used as a local search procedure. The general framework of the proposed MA for the MSRCMPSP is presented in Algorithm 4. The algorithm starts by generating an initial population of solutions. To determine good initial solutions, and based on experimental results conducted in this research (Table 5), the SSGS with a randomized version of priority rule LS (using a geometric distribution with \(p=0.5\)) as a priority rule is used. The randomization of LS ensures the diversification of the initial population. Then, while any of the stopping criteria is not met, two parents are randomly selected from the population. A one-point crossover operator is used to generate two children. Each child has a probability of \(P_m\) to mutate. The mutation operator implemented in this work randomly changes the position of an activity in the sequence by preserving the precedence feasibility. Contrary to the classical generational GA and to avoid the loss of chromosomes with high qualities, a steady-state replacement strategy (Syswerda 1991) based on a simulated annealing method is employed. The incremental replacement can lead to a low population diversity (Essafi et al. 2008) but simulated annealing helps to overcome this drawback. Each child is accepted or not according to the simulated annealing procedure. If the child is accepted, an individual among \(N_{worst}\) individuals is replaced with the accepted child. A common strategy is to replace the worst individual (Chen and Shahandashti 2009) but, in this case, the Metropolis acceptance criterion loses its sense. In fact, there is a high probability that the accepted child with a worse objective than the current population may be replaced before undergoing the crossover process. Accepting worse solutions into the population not only helps in the diversification of the population (which is crucial for GA) but also it helps to escape from a local optimum. To ensure global convergence, a cooling scheme is applied.

Algorithm 4
figure d

Memetic algorithm (MA)

Full size image

The local search procedure of the proposed memetic algorithm is a SA algorithm. SA is known for its ability to escape local optima and many papers found that SA leads to very good results when solving complex scheduling problems (Bouleimen and Lecocq 2003; Yugma et al. 2012; Tamssaouet et al. 2022). The neighborhood function implemented in the SA procedure is the Swap activity move (Asta et al. 2016). To maintain the precedence feasibility of the sequence list, the move starts by randomly selecting an activity \(a_1\). Then, a second activity \(a_2\) between the positions of the latest predecessor and the earliest successor of \(a_1\) (with \(a_1 \ne a_2\)) is selected. The move is accepted or not according to the Metropolis criterion. There are two main methods to apply the local search procedure in a memetic algorithm (Moscato and Cotta 2003): (i) As a mutation operator of the GA algorithm and (ii) After each iteration to all or a part of the population. Many papers underline that the quality of a population based metaheuristic results from the interplay between intensification and diversification (Sörensen and Sevaux 2006). Intensification may conduct the population to be very close to the few first improved individuals (using local search) since the mutation probability is usually small. Consequently the population loses its diversification. Diversification is very time consuming (Gonçalves et al. 2005; Sevaux and Dauzère-Pérès 2003). To overcome these drawbacks, in this work, the SA local search procedure is applied after every \(GA_{iter}\) iterations of the evolutionary process to \(N_{SA}\) random individuals of the population.

All the metaheuristics implemented in this paper work on a precedence feasible activity list. To evaluate the objective function, Algorithm 5 converts the activity list to a time feasible schedule.

Algorithm 5
figure e

Time schedule construction algorithm

Full size image

5 Computational experiments

Section 5.1 presents the design of the computational experiments and how they were conducted. In Section 5.2, the MILP model is compared to MA on 10 small instances. Section 5.3 investigates the performances of SSGS, PSGS and the 7 priority rules. In Sect. 5.4, five solution approaches (Greedy, SA, GA, hSGA and MA) are compared on 27 large instances. The MA is evaluated on benchmark instances in Sect. 5.5. Finally, sensitivity analyses of the MA are discussed in Sect. 5.6.

Table 2 Characteristics of real industrial instances: \(|\) \(\mathcal {E}\) \(|\) is the number of projects, \(|\) \(\mathcal {A}\) \(|\) is the number of activities, SP is the serial/parallel indicator of the global project, avg\(|Pred|\) is the average number of the predecessors for each activity, \(|\) \(\mathcal {R}\) \(|\) is the number of global human resources (teams), avg\(|\) \(\mathcal {K}_r\) \(|\) is the average number of skills per team, max\(|\) \(\mathcal {K}_r\) \(|\) is the maximum number of skills per team, avg\(SS_k\) is the average skill strength of skills, \(|\) \(\mathcal {M}\) \(|\) is the number of machines, avg\(|\) \(\mathcal {I}_m\) \(|\) is the average number of installations (skills) per machine, max\(|\) \(\mathcal {I}_m\) \(|\) is the maximum number of installations (skills) per machine
Full size table

5.1 Design of experiments

The performances of the proposed solution approaches are analyzed using real industrial instances of SNCF. The instances are extracted from the manufacturing execution system (MES) database of the maintenance centers and are feasible. All the numerical experiments are carried out on the two considered objective functions: Minimization of the weighted sum of the completion times of the projects and minimization of the weighted sum of the tardiness of the projects.

First, the MILP model is compared with the proposed MA on 10 small instances defined based on industrial data. Table 2 provides some properties of the small instances. The number of projects varies from 1 to 4 and the number of activities varies from 48 to 156. The serial/parallel indicator SP (defined in Vanhoucke et al. (2008) as \(I_2\)) is lower than 0.5, which means that many activities can be performed in parallel. Note that, since the instances consist of multiple projects, the SP values are computed using the global graph, which is formed by the precedence graphs of each individual project. The number of teams varies from 4 to 8 with approximately 3 skills per team. The average skill strength avg\(SS_k\) (Snauwaert and Vanhoucke 2023) is relatively small, which indicates that all the required skills are scarce, posing a significant challenge when solving the instances. The number of machines varies from 38 to 49 and with an average number of installations per machine from 1.2 to 1.9. For the small instances, very few machines (one or two) were critical. To diversify the instances, data from two different heavy maintenance centers are used.

The MILP model is solved by the standard solver IBM ILOG CPLEX 20.1 with default parameters. The memetic algorithm is implemented in C++ and the numerical experiments were carried out on a personal computer with a 1.60 GHz processor and 16 Gb RAM. To compare the convergence efficiency over time, the computational time is limited first to 600 s and then to 3600  s. The results are summarized in Sect. 5.2 and show that MA outperforms the MILP model on both the quality of solutions and the convergence efficiency.

Then, computational experiments are conducted on 27 large instances with 14 to 380 projects and 1539 to 7119 activities. The large instances correspond to real cases and are provided from three heavy maintenance centers which repair different types of rolling stock units. Since many projects are considered simultaneously, a large number of activities can be executed in parallel. Consequently, the SP values of the large instances are even smaller in comparison to the SP values of the small instances. The number of teams varies from 9 to 21 with an average number of skills per team between 2.9 and 5.2. The average skill strength is very small (lower than 0.14) because, at the tactical level, the resources are sized to reduce the idle times. The number of machines varies from 52 to 107 with an average number of installations per machine between 1.5 and 2.1. Additional properties of the instances can be found in Table 2.

In Sect. 5.3, SSGS is compared to PSGS with the 7 priority rules described in Sect. 4.1. The results show that SSGS outperforms PSGS and that LS is the best priority rule for both objective functions. Five algorithms are compared in Sect. 5.4: Greedy, SA, GA, hSGA and MA. The algorithms are implemented in C++ and empirically parameterized in the following manner:

  • Greedy: Refers to SSGS+LS,

  • SA: Simulated Annealing algorithm with the following parameters:

    • Initial temperature: \(T_{init}=24\) for SWTP and \(T_{init}=48\) for SWDP,

    • Cooling factor: \(C=0.99999\),

  • GA: Genetic Algorithm with the following parameters:

    • Population size: \(N_{pop}=120\),

    • One-point crossover probability: \(P_c=0.95\),

    • Mutation probability: \(P_m=0.75\),

    • Linear ranking selection (Hartmann 1998) method is used for selecting the next generation,

  • hSGA: hybrid Simulated Genetic Algorithm with the same parameters as above and \(N_{rem}=40\),

  • MA: Memetic Algorithm with the same parameters as hSGA and \(GA_{iter}=2500\) iterations, \(N_{SA}=4\) and \(SA_{iter}=2000\).

The parameter tuning is done step by step. First the SA parameters are tuned as in Knopp et al. (2017). By computing the first 100 moves which deteriorates the solution, the 2% percentile on some training instances was calculated. We found that, in most of the cases, the 2% percentile is approximately 24 for the SWTP objective and approximately 48 for the SWDP objective. Note that 24 correspond to 24 hours, which is the duration of most of the activities. Secondly, by empirically testing different values (as in Hartmann (1998)) of the population size, mutation and crossover probabilities, the parameters of the GA were set. The mutation probability (0.75) is notably larger than what is commonly found in papers of the literature that use GA for solving the RCPSP. However, it seems that, for complex and large problems, it is usually better to have a high mutation probability (Elloumi and Fortemps 2010; Murata et al. 1996). For the hGSA, the same parameters of SA and GA are used but; to keep a diversified population and to exploit the advantages of the simulated annealing approach, several values of \(N_{rem}\) (10, 20, 30, 40, 50) were tested. The best results were obtained for \(N_{rem}=40\). Finally, for the MA, the value of \(GA_{iter}\) is the number of iterations required for the GA algorithm to converge towards a good solution. After reaching approximately 2 500 iterations (equivalent to 5 000 schedules), the convergence of the GA significantly decelerates. Adjusting parameters \(N_{SA}\) and \(SA_{iter}\) is challenging, because large values may lead to premature convergence (Sörensen and Sevaux 2006). Consequently, small values are preferred and have been tested by trying different values (2, 4, 6, 8, 10 for \(N_{SA}\) and 2000, 4000, 6000, 8000, 10000 for \(SA_{iter}\)).

The computational time limits were chosen to evaluate the quality of the approaches in a limited amount of time (600 s) and if longer computational times are allowed (1800 and 3600 s). The requirements in terms of computational times from the planners of SNCF depend on how the optimization is used (for generating an initial complete schedule or to perform simulations by varying input parameters).

5.2 Comparison of MILP and MA

In this section, the performances of the MILP model and the proposed MA are compared. Table 3 shows the numerical results and the percentage improvement PI (in %) after 300 s of computational time, while Table 4 shows the results obtained after 3600 s of computational time.

PI is defined as follows:

$$\begin{aligned} PI(\%)=\dfrac{MILP_{Obj}-MA_{Obj}}{MILP_{Obj}}*100, \end{aligned}$$

where \(MA_{Obj}\) is the objective value of the best solution determined by MA and \(MILP_{Obj}\) is the objective value found by CPLEX. Column LB is the lower bound found by CPLEX after 3600 s, column Gap is the gap between LB and the objective value of the MILP, while Gap\(_{MA}\) is the gap between LB and the solution found by MA. The cells filled with “–” indicate that no feasible solution is found by CPLEX, while the cells filled with “*” indicate that the solution found by CPLEX is optimal.

Table 3 MILP vs. MA—Numerical results and percentage of improvement with CPU time limited to 300s
Full size table
Table 4 MILP vs. MA—Objective functions and percentage of improvement with CPU time limited to 3600 s
Full size table

Regarding the SWDP objective function, Table 3 shows a clear dominance of MA, which obtains the best solution for 8 instances out of 10 and equivalent solutions for the 2 remaining instances. The improvement becomes particularly significant when the instances have more than 90 activities. CPLEX obtained an optimal solution in less than 300 s for only 2 instances. For Instance N_3_144, the solver could not find a feasible solution. For Instances N_3_144 and S_4_156, Gap\(_{MA}\) is relatively large (6% and 7.5%), but as CPLEX struggles to find good solutions, it is likely that the lower bound LB is of poor quality.

Regarding the SWTP objective function, MA found the best solution for 8 instances out of 10 and equivalent solution for Instance S_2_58 with a large PI (>78%) for instances with more than 90 activities. For Instance N_2_57, the solver outperforms the proposed MA. Instance N_2_57 has the smallest avg\(SS_k\) and also the average number of skills per machine is larger than in other instances with less than 60 activities, which makes it harder to solve to optimality.

However, after 3600 s of computational time, MA obtains the optimal solution for Instance N_2_57. All instances with less than 60 activities were optimally solved by both CPLEX and MA. For the other instances, MA outperforms CPLEX with a PI larger than 15% for the SWDP objective function and larger than 65.9% for the SWTP objective function. For Instance N_3_144, the solver still could not find a feasible solution while MA reduces its gap from the LB from 6% to 2.3% after 3600 s, which indicates the capability of MA to address hard instances. Concerning the SWTP objective function, large gaps between the solution of MA and the lower bound for the Instances S_3_127, N_3_105 and N_3_144 (respectively 29.5%, 28.6% and 27.8%) can be observed. This is likely due to the small values of the objective function and the poor quality of the lower bounds.

Tables 3 and 4 illustrate the superiority of the proposed MA over CPLEX for the MSRCMPSP. The percentage improvement increases with the size of the instances. The performances of CPLEX could be improved by better tuning the parameters or applying a warm start for example. However, as the real instances include several thousand activities, CPLEX is not adapted to solve the industrial instances. Memory limitations actually prevented us from obtaining even a lower bound for the large real instances.

5.3 Performance of greedy algorithms

In this section, we investigate the performances of SSGS and PSGS tested with the 7 priority rules detailed in Sect. 4.1. Table 5 summarizes the numerical results obtained using the 27 large instances (the detailed results for each instance can be found in Table 11 of Appendix A). Column "#best" gives the number of instances where the corresponding priority rule found the best solution. Column "avg gap (%)" shows the average gap (in %) from the best solution and "max gap (%)" is the maximal gap among the 27 instances. The computational times are not given since both SSGS and PSGS need less than 1 (approximately 5ms) to compute a solution.

Table 5 SSGS vs. PSGS: Performance of priority rules and associated gap from the best solution found by SSGS or PSGS
Full size table

The first observation is that the serial scheduling scheme is superior to the parallel scheduling scheme. Column "#best" of Table 5 shows that SSGS found a better solution than PSGS for the 27 considered instances. The average gap of PSGS is at least 47% for the SWDP objective and at least 90% for the SWTP objective. Kolisch and Hartmann (2006) and Hartmann and Kolisch (2000) concluded that SSGS is better than PSGS for complex problems. In this paper, this significant gap can be explained by the large size of the problem. Having several projects with different characteristics and deadlines makes the problem very complex to solve. The gap is even larger for the SWTP objective function since scheduling in parallel is very myopic regarding the deadline tightness of projects.

Among the 7 considered priority rules, the best one is LS with 22 best solutions found out of 27 for both objective functions. The second best priority rule is LF which found 4 best solutions out of 27 instances for SWDP and 5 best solutions out of 27 instances for SWTP. LS is also the most robust priority rule since LS has the smallest average gap and the smallest maximum gap for both objective functions. Hence, SSGS+LS is used to generate the initial solution of the different metaheuristics. When a population of initial solutions is needed, a randomized version of LS is used as a priority rule to generate the individuals. The numerical results obtained by the metaheuristics are discussed in the next section.

5.4 Comparison of heuristics

Five algorithms are compared in this section: Greedy, SA, GA, hSGA and MA. The results are presented after 600, 1800 and 3600 s for both considered objective functions. The detailed results for each large instance are presented in the 6 tables of Appendix B which have a similar structure (Tables 12, 13, 14, 15, 16, and 17). The tables show the objective value for each algorithm and the associated gap (column "%Gap") from the best solution highlighted in bold.

A summary of the results can be found in Table 6 for the SWDP objective and in Table 8 for the SWTP objective. Column "#best" provides the number of best solutions found by the corresponding algorithm over the 27 large instances. To analyze the efficiency and the robustness of each algorithm, the average, the standard deviation and the maximal gap from the best solution are also provided.

Table 6 Summary of results on the performances of the metaheuristics for objective function SWDP
Full size table

Let us first analyze the results obtained for the SWDP objective. When the computational time is limited to 600 s, from Tables 12 and 6, note that, out of 27 instances MA finds the best solution for 19 instances and hSGA finds the best solution for 9 instances. No best solution is determined by the other algorithms. Consequently, MA and hSGA outperform the Greedy algorithm, GA and SA. MA is slightly better than hSGA since the average gap for MA is equal to 0.1% while the average gap of hSGA is around 0.4%. When compared to the greedy algorithm which is used to generate the initial solutions, MA improves the solutions by around 8%. The maximal gap of the Greedy algorithm is 15.6%. Even if SSGS+LS (i.e., the Greedy algorithm) seems to have good results in Sect. 5.3 compared to PSGS and the other priority rules, still a significant improvement is obtained with the proposed MA.

Regarding the maximal gap, MA is also robust because its worst gap is around 0.8%. The worst gap is obtained for Instance S_21 which has the larger number of activities. In the evolutionary process, many decisions are made randomly. Thus, for large instances it is more difficult to control the convergence. Since GA has worse results compared to SA (with respectively an average gap of 1.7% and 4.0% in Table 6), the fact that the simulated annealing procedure helps to escape local optimum is very important when dealing with large instances and complex problems. Similar conclusions were drawn in Tamssaouet et al. (2022).

Furthermore, after 1800 s of computational time, SA closes the gap from the best solutions. Table 6 shows that SA has an average gap of 1.7% and a maximal gap of 3.5% after 600 s but that after 1800 s the average gap reduces to 1.4% and the maximum gap to 2.8%. SA even finds the best solution for Instance N_b_21_3 (Table 13). However, MA remains the best metaheuristic with 18 best solutions out of 27, followed by hSGA with 8 best solutions. The gap of the Greedy algorithm increases (from 8.1% to 8.6%), meaning that the metaheuristics can still improve the solution after 600 s. In particular, the gap difference between the greedy Algorithm and MA increases from 8% to 8.5%.

Similar observations can be derived when the computational time is limited to 3600 s. SA still closes the gap (from 1.4% to 1.2%) and the number of best solutions found is the same as after 1800 s. However, the best solutions determined by MA, hSGA and SA are not for the same instances as in Table 13. In Table 13, the best solution found by SA was for Instance N_b_21_3 but, in Table 14 the best solution found by SA is for Instance N_22_3. Another difference can be observed for Instances N_b_22_ot and N_b_21_1 where MA overcomes hSGA. However, for Instances S_21_2 and S_22_3, hSGA overcome MA. The gap of the Greedy algorithm still increases from 8.6% after 1800 s to 8.8% after 3600 s. Details on the percentage of improvement (PI) over time are illustrated in Table 7. These results confirm the complexity of the problem and that finding the optimal solution is difficult.

Table 7 Convergence performances of the metaheuristics: Mean of the PI (in %) of the initial solutions provided by SSGS
Full size table

Considering the results presented so far, MA is the most appropriate approach to solve the MSRCMPSP with SWDP as objective function. SA can close the gap over time, which makes SA a very good candidate for the local search procedure of a memetic algorithm.

The same analysis is conducted for the SWTP objective function. Table 8 presents the summary of the results for the 27 instances.

Table 8 Summary of results on the performances of the metaheuristics for objective function SWTP
Full size table

After 600 s of computational time, MA and hSGA both find the best solution for 14 instances. One best solution was found by GA for Instance S_22_2 but the gap with the solution found by MA is only 0.9%, which represents less than one day of total delay for 69 projects (Table 15). MA remains the best approach with an average gap of 2.3%, while the second best approach is hSGA with an average gap of 3.7%. MA also has the best maximal gap (11.6%) and the best standard deviation gap (3.5%), confirming the robustness of MA compared to the other approaches. In particular, the average gap of the Greedy algorithm is 51%, which implies that the tardiness is approximately reduced by a factor of 2. SA and GA have good results as well with an average gap smaller than 15% but their efficiency is not stable since the maximum gap is respectively 33.3% for SA and 29.3% for GA. Once again, this instability can be justified by the complexity of the problem with many projects, activities, teams, machines, skills and constraints. When neighborhood exploration is too random, it may be difficult for the metaheuristics to select promising moves. Table 5 in Section 5.3 shows significant gaps as well, i.e., SSGS and PSGS determine solutions that are far from an optimal solution.

When the computational time is limited to 1800 s, several observations can be made (Table 16):

  • hSGA finds 13 best solutions whereas MA still finds 14 best solutions out of 27 and GA determines one additional best solution, i.e., 2 best solutions.

  • As for the SWDP objective, hSGA converges better than MA for some instances towards the best solution after 1800 s (for example Instance S_21_3). On the contrary, MA overcomes hSGA for Instance N_21_2.

  • SA closes the gap from 14.6% after 600 s to 13% after 1800 s. The other metaheuristics do not considerably improve the solutions since the gap of the Greedy algorithm deteriorates by only 1%. MA and hSGA struggle to improve the solutions (Table 7) and consequently, their average gap slightly deteriorates (from 2.3% to 2.5% for MA and from 3.7% to 3.8% for hSGA).

However, MA improves the maximal gap, in particular when the computational time is limited to 3600 s. The maximal gap of MA decreases from 11.6% after 600 s to 8.8% after 3600 s, which confirms its robustness. SA, GA and hSGA still have significant maximal gaps (respectively 27%, 31.9% and 19.2%). hSGA overcomes GA for Instances S_22_2 and S_21_2 for which GA finds a better solution after 1800 s (Table 16). As for the SWDP objective, accepting worse solutions using the Metropolis criterion allows local optima for the SWTP objective to be escaped.

In view of the results, the intensification in the memetic algorithm, ensured by the local search procedure, increases the average quality of solutions and the robustness of the proposed solution approach. Hence, MA is the best approach to solve the MSRCMPSP. However, in terms of best solutions, hSGA and MA are quite equivalent when SWTP is considered as objective function. These results illustrate that there is still room for improvement and in particular in the intensification mechanism ensured by SA, as SA consistently improves its average gap (for example from 13.0% after 1800 s to 11.9% after 3600 s). An interesting perspective would be to better choose the visited neighborhood. The evaluation of moves is very time-consuming, by discarding uninteresting moves considerable computational time can be saved (Dauzère-Pérès and Paulli 1997; Mati et al. 2011). Hence, with the same time allocated to the intensification part, a more promising neighborhood could be explored. This would probably improve the convergence and the results of the proposed memetic algorithm and the results of SA as well.

5.5 Comparison on benchmark datasets

To the best of our knowledge a benchmark dataset with both multiple projects and multiple skills does not exist in the literature. Since our objective functions are only meaningful in the context of multiple projects, we conducted experiments on the instances of the Multi-Project Scheduling Problem library (MPSPLib: http://www.mpsplib.com/, July 2023). The library contains 20 datasets with a total of 140 instances built by combining several RCPSP instances of Kolisch and Sprecher (1997). The number of projects varies from 2 to 20 and the number of activities varies from 60 to 2 400. For more details about the characteristics of these datasets, the reader is referred to Wauters et al. (2015). The memetic algorithm is compared with the sequential learning-based metaheuristic of Wauters et al. (2015), which is reported to be one of the best (decentralized) solution approaches of the literature (Bredael and Vanhoucke 2022). Their solution approach consists of a sequence learning game played by several project managers. Each manager learns both a local activity list (using reinforcement learning) and a global activity list containing all the activities of all projects. The global activity list is build by adding all the activities of the first project, then the activities of the second project, and so on. Finally, the global list is transformed into a feasible schedule by using the classical serial generation scheme. To ensure a fair comparison, the same stopping criterion (i.e. 100 000 schedules) is used. After preliminary experiments and since the objective values are smaller than the objective values of the real instances, an initial temperature \(T_{init}=2.5\) was set. The other parameters of the proposed memetic algorithm remain unchanged.

Table 9 presents the results of the first configuration (ADP1, Time1 in seconds, and Gap1 in %). The proposed memetic algorithm shows good performance when the average project delay (ADP) is small (less than 30). For example, better solutions are found by MA for Instances MP120_20, MP30_2, MP30_5, MP90_2 and MP90_5. On the contrary, when the ADP is large, MA does not perform well.

Table 9 Comparison of average project delay (ADP) with Wauters et al. (2015)
Full size table

The problem subsets with large ADP have an significant resource Average Utilization Factor (AUF) (Kurtulus and Davis 1982 and the MPSPLib web site). In this case, and as it is underlined in Asta et al. (2016) and Bredael and Vanhoucke (2022), the prioritization of the projects is very important. At the SNCF maintenance centers, the arrival times of the rolling stock (i.e. the starting time of the projects) are decided at a tactical level to avoid overloading the available resources. Hence, the prioritization of the projects was not really necessary. However, we have implemented several neighborhood operators and, among them, a novel one which handles the project priorities (more precisely the priorities of the different phases of a project). The novel neighborhood operator (SortActivities) consists of selecting a sublist from the global activity list, with a random size uniformly distributed between 3 and \(|\) \(\mathcal {A}\) \(|\)/2. The activities of the sublist are then sorted by project delay while conserving the initial order of the activities of a same project. For example, if the following sublist is selected: [1(1), 2(2), 3(2), 4(1)] (where a(n) means that the activity a of project n is in the sublist), and project 1 has a larger delay compared to project 2, after applying the SortActivities neighborhood function, the sublist becomes: [1(1), 4(1), 2(2), 3(2)]. The SortActivities neighborhood function is added as a mutation operator in the memetic algorithm of Sect. 4.2. The memetic algorithm has now two mutation operators with the same probability of being selected. The original one is kept to help the exploration of additional solutions by the SortActivities function. This is the second configuration of our proposed memetic algorithm and the results (ADP2, Time2 (in seconds), and Gap2 (in %)) are detailed in Table 9. The gap is considerably reduced for the instances with large ADP values. For Problem subset MP120_20AC, Gap1 was −108% and the second configuration reduces the gap to −4.4%. Generally speaking, the average gap over all instances is reduced from −32.2% to 6.4%. The memetic algorithm with the second configuration obtains better results for 15 datasets out of 20 compared to the learning approach of Wauters et al. (2015), leading to an average improvement over all instances of 6.4%. In terms of computational times, the proposed memetic algorithm is on average at least 3 times faster for the first configuration and 4 times faster with the second configuration (note that, when the solution quality is poor, Algorithm 3 requires more time, which explains the difference of computational times between Time1 and Time2).

Let us note that, for the instances with small ADP values, the second configuration still finds better solutions compared to the learning approach of Wauters et al. (2015), but does worse than the first configuration (for Instances MP30_2, MP90_2, MP90_5 and MP30_5). In view of the results, the SortActivities mutation operator should be used when the problem has a large AUF. That is why the SortActivities mutation operator is not useful to solve the real instances of SNCF.

The column GapBest gives the best gap between Gap1 and Gap2. If the appropriate configuration of the memetic algorithm is chosen, an overall average improvement of the ADP of 7.7% can be achieved. Even if our algorithm was not initially built to solve the MSPSLib instances, it shows very promising results and could be generalized to solve other related problems.

5.6 Sensitivity analyses

An efficient hybrid metaheuristic requires the configuration of many parameters (Pellerin et al. 2020). Table 10 reports the average gap from the best solution obtained by the different configurations of MA. The computational experiments are conducted on the real instances of SNCF and MA is stopped after 600 s. In this sensitivity analysis, for the sake of brevity and, since they have been identified as the most influential factors during preliminary experiments, we exclusively focus on the parameters \(N_{pop}\), \(T_{init}\), \(N_{rem}\) and \(P_m\).

Table 10 Sensitivity of MA: Impact of input parameters with CPU time limited to 600 s
Full size table

In each row of Table 10, only a single parameter is changed from the initial configuration of MA. Regarding the SWPD objective function, the largest average gaps are obtained for small values of \(P_m\) and \(N_{rem}\) (resp. 3.1% and 2.7%), which are both useful for keeping a diversified population. Similar conclusions can be drawn when the objective function is SWTP.

Regarding the population size, note that, for the SWTP objective function, a larger value (180) is more suitable while, for the SWDP objective function, a lower population (between 80 and 120) performs better. The primary reason could be that numerous individuals yield the same objective values for the SWTP objective function. In fact, a neighborhood move can lead to a different solution but it may not impact the tardiness of the projects and hence, a larger population is necessary to prevent the population to be homogeneous (Van Peteghem and Vanhoucke 2010). This observation could also indicate why the suggested MA demonstrates improved performance with larger values of \(P_m\) and \(N_{rem}\) for the SWTP objective function (resp. 0.85 and 50) compared to the SWDP objective function (resp. 0.75 and 40).

This analysis shows the importance of a proper parameter tuning for the success of a (hybrid) metaheuristic. Except for the initial temperature, we chose a uniform configuration for both of the considered objective functions to facilitate the use of the algorithm.

6 Conclusions

This paper studies an original industrial RCPSP problem with multiple skills and multiple projects, which corresponds to a real case in the heavy maintenance centers of railway rolling stock. Two objective functions are considered: (i) Minimization of the sum of the weighted tardiness of the projects and (ii) Minimization the sum of the weighted duration of the projects. A time-indexed MILP formulation of the problem is provided, and two constructive heuristics (serial and parallel scheduling generation schemes) are implemented and tested with 7 priority rules on 27 large real instances. The solutions provided by SSGS+LS are used as initial solutions and are improved using four metaheuristics: SA, GA, hSGA and MA. Several computational experiments were conducted in this study. The first experiment compared the MILP model and the proposed memetic algorithm (MA) on 10 small instances defined using industrial data. The results showed that MA outperforms the MILP model after 300 and 3600 s of computational time, in particular for instances with more than 90 activities. Larger gaps were observed when the objective function was to minimize the sum of the weighted tardiness of the projects (SWTP). In the second experiment, the performance of 7 priority rules implemented in the serial scheduling generation scheme (SSGS) and the parallel scheduling generation scheme (PSGS) were evaluated using 27 large instances. SSGS found the best solution for the 27 instances, while none were found by PSGS. Among the priority rules, the Latest Start (LS) rule had the best performance, with 22 best solutions found out of 27 for both objective functions.

Finally, the performance of the memetic algorithm was compared to SSGS+LS, simulated annealing (SA), genetic algorithm (GA), and hybrid simulated annealing and genetic algorithm (hSGA) with computational times limited to 600, 1800, and 3600 s. MA has the best average gap and the smallest maximal gap, indicating a high level of stability. In particular, for the SWDP objective function, an average improvement of 8% of the initial solutions provided by the Greedy algorithm was observed. For the SWTP objective function, the total delay was reduced by a factor of 2. Furthermore, an overall average improvement of the ADP of 7.7% was achieved for the benchmark instances. In view of the results, the memetic algorithm is being integrated into the MES of SNCF maintenance centers. The scheduling horizon is decided by the planners, and all the activities starting beyond the starting date of the scheduling horizon are (re)scheduled.

Despite the promising results obtained by the memetic algorithm, several areas for improvement can be identified. One potential direction for future research is to conduct a more advanced sensitivity analysis of the proposed algorithm to optimize its configuration. The proposed approach requires the initialization of numerous parameters, and characterizing the most robust configuration can be a challenging task. A self-adapting component would be of great value for MA to solve other SNCF instances in the future, which may be different from the considered instances. Additionally, as previously discussed, a better exploration of promising neighborhood moves may save computational time and improve the convergence of the memetic algorithm.

Another area of future research is to consider the stochastic version of the problem. In the context of heavy maintenance, many operations are performed by human operators and processing times are uncertain. Moreover, uncertain additional tasks may also be encountered. Therefore, defining robust schedules would be helpful to meet customer deadlines. A possible robustness measure would be to maximize the chances of meeting project deadlines. Similar robustness measures, such as the service level, are proposed in Dauzère-Pérès et al. (2008).