## 6.1 Data

The dataset contains 40 instances, with $|P| \in [15, 25]$, $|C| \in [4, 10]$, $T = 12, B = 4$, $|D| \in [20, 50]$. Maximum running time is set to 30 minutes. The dataset is based on past data for the generation of probabilistic distributions. Figure 3 (left) shows an example of a set of distribution for one container. As we can see, the more the building has to make products, the higher the chance to consume.

## 6.2 Model

For each product $p$, variables $p_{p,t} \forall t \in T$ indicates if the product is made by the team $t$. For each delivery $d$, variable $d_{d,c} \forall c \in C$ indicates if the delivery is stored in container $c$. For each container $c$, $y_c$ is the emptying random variable. For each building $b$, $q_b$ is the quantity of product produced in the building. Let $q = (q_1, \dots, q_{|B|})$ be a vector of all $q_b$ variables. Let $x_c \forall c \in C$ be the additional capacity required by container $c$ such that it is no longer overflowing.

$$ \arg \max_{x \in \mathbb{Z}, Y \in \mathbb{P}} \gamma \qquad (26) $$

$$ \text{Such that} \quad \sum_{t=1}^{|T|} p_{i,t} = 1, \quad \forall i \in P \qquad (27) $$

$$ \sum_{i=1}^{|P|} p_{i,t} w_i \le W_t, \quad \forall t \in T \qquad (28) $$

$$ q_b = \sum_{i=1}^{|P|} \sum_{t \in \{t_i | b_{t_i} = b\}} p_{i,t}, \quad \forall b \in B \qquad (29) $$

$$ \sum_{i=1}^{|D|} d_{i,c} q_i \le C_c + x_c, \quad \forall c \in C \qquad (30) $$

$$ \sum_{c=1}^{|C|} d_{i,c} = 1, \quad \forall i \in D \qquad (31) $$

$$ \text{Table}(q, y_c), \quad \forall c \in C \qquad (32) $$

$$ p_{x_c \le y_c} \equiv x_c \le y_c, \quad \forall c \in C \qquad (33) $$

$$ \Pr[(p_{x_1 \le y_1}, \dots, p_{x_{|C|} \le y_{|C|}}) = 1] \ge \gamma \qquad (34) $$

## 6.3 Results

Figure 3 shows the results of the fixed and variable (this paper) methods. First, it is interesting to see that for some instances a confidence of 1 can be found. For these instances, no consumption was required to find a solution. In contrast, we can also see that for some instances, it is hard to find solutions with more than 2 percent confidence. Those instances are hard instances, where the delivered amount is too large. Then in the middle of the plot, it is shown that the confidence of solutions found by the variable method is higher in general than the fixed method. In addition, the variable method makes one order of magnitude less decisions to find better or equivalent solutions. When profiled, the reason why less decisions are made is the repetitive call to the continuous integrated solver. An actual hybrid CP solver would be drastically more efficient. Nevertheless, even with this drawback, in both experiments, the variable showed significant improvements in terms of quality of solution and time.

# 7 Conclusion and Future works

This paper proposed to extend the CP framework to distribution optimization. First, random variables have been extended to the case of distribution variables, then the CP constraints have been extended to deal with this extension. We proposed a generalization of integer
<figure><img src="image_290.png" /><figcaption>Figure 3 Instances are ordered w.r.t. size. (left) Range of probability distributions for a random variable having 9 possible values. The gradient of color indicate the increase of $q$. (middle) Best confidence comparison. (right) Number of decisions (log base 10).</figcaption></figure>

variables and new consistencies for probability distributions. In addition, generic filtering algorithms have been proposed, pruning invalid distributions. We defined an implementation based on a probability mass function decomposition and the minimal arithmetic to model most problems, together with the associated filtering algorithms. Finally, as shown in the experimental section, we used the proposed framework to solve two optimization problems where the distributions of probability are not fixed at the beginning of the problem.

The main future direction is to no longer restrict the search to a finite set of distributions using a table as done in our experiments, but to be able to search directly into the distribution space, which is closer to known methods in machine learning. Doing this would bring CP as one of the main frameworks to do constrained distribution learning, and will require hybridization of CP solvers. Other future directions include the design of specialized global constraints filtering, the generalization of the reuse of existing filtering algorithms [18], and the implementation of different types of encoding of the probability distribution variables. Another important direction is the encoding of dependent random variables. In this paper, we purposely introduced the event space of probability space to encode the dependency between random variables. In future works, constraints on the event spaces will lead to efficient dependency implementation. Finally, this paper is restricted to 2-stage policy, it will be important in the future to extend it to multi-stages.

References

1. Eric Baum and Frank Wilczek. Supervised learning of probability distributions by neural networks. In *Neural information processing systems*, 1987.

2. N Bharanidharan and Harikumar Rajaguru. Improved chicken swarm optimization to classify dementia mri images using a novel controlled randomness optimization algorithm. *International Journal of Imaging Systems and Technology*, 30(3):605-620, 2020.

3. Giuseppe Calafiore and Marco C Campi. Uncertain convex programs: randomized solutions and confidence levels. *Mathematical Programming*, 102(1):25-46, 2005.

4. Gilles Chabert et al. Ibex, an interval-based explorer, 2007.

5. Gilles Chabert and Luc Jaulin. Contractor programming. *Artificial Intelligence*, 173(11):1079-1100, 2009.

6. Abraham Charnes and William W Cooper. Chance-constrained programming. *Management science*, 6(1):73-79, 1959.

7. Abraham Charnes and William W Cooper. Deterministic equivalents for optimizing and satisficing under chance constraints. *Operations research*, 11(1):18-39, 1963.

8. Abraham Charnes, William Wager Cooper, and MJL Kirby. Chance-constrained programming: an extension of statistical method. In *Optimizing methods in statistics*, pages 391-402. Elsevier, 1971.
9 Constantinos Daskalakis, Ilias Diakonikolas, and Rocco A Servedio. Learning poisson binomial distributions. In Proceedings of the forty-fourth annual ACM symposium on Theory of computing, pages 709–728, 2012.

10 Constantinos Daskalakis and Gautam Kamath. Faster and sample near-optimal algorithms for proper learning mixtures of gaussians. In Conference on Learning Theory, pages 1183–1213. PMLR, 2014.

11 Augustin Delecluse, Pierre Schaus, and Pascal Van Hentenryck. Sequence variables for routing problems. In 28th International Conference on Principles and Practice of Constraint Programming (CP 2022), 2022.

12 Mihai Dobrescu, Katerina Argyraki, and Sylvia Ratnasamy. Toward predictable performance in software {Packet-Processing} platforms. In 9th USENIX Symposium on Networked Systems Design and Implementation (NSDI 12), pages 141–154, 2012.

13 Paul Emmerich, Daniel Raumer, Florian Wohlfart, and Georg Carle. Assessing soft-and hardware bottlenecks in pc-based packet forwarding systems. ICN 2015, 90, 2015.

14 Hélène Fargier and Jérôme Lang. Uncertainty in constraint satisfaction problems: a probabilistic approach. In Symbolic and Quantitative Approaches to Reasoning and Uncertainty: European Conference ECSQARU'93 Granada, Spain, November 8–10, 1993 Proceedings 2, pages 97–104. Springer, 1993.

15 Jakob Feldtkeller, David Knichel, Pascal Sasdrich, Amir Moradi, and Tim Güneysu. Randomness optimization for gadget compositions in higher-order masking. Cryptology ePrint Archive, 2022.

16 Juliver Gil Herrera and Juan Felipe Botero. Resource allocation in nfv: A comprehensive survey. IEEE Transactions on Network and Service Management, 13(3):518–532, 2016.

17 Brahim Hnich, Roberto Rossi, S Armagan Tarim, and Steven Prestwich. Synthesizing filtering algorithms for global chance-constraints. In Principles and Practice of Constraint Programming-CP 2009: 15th International Conference, CP 2009 Lisbon, Portugal, September 20-24, 2009 Proceedings 15, pages 439–453. Springer, 2009.

18 Brahim Hnich, Roberto Rossi, S Armagan Tarim, and Steven Prestwich. Filtering algorithms for global chance constraints. Artificial Intelligence, 189:69–94, 2012.

19 JN Hooker. Stochastic decision diagrams. In International Conference on Integration of Constraint Programming, Artificial Intelligence, and Operations Research, pages 138–154. Springer, 2022.

20 Mark H Houck. A chance constrained optimization model for reservoir design and operation. Water Resources Research, 15(5):1011–1016, 1979.

21 Bahareh Kargar, Mir Saman Pishvaee, Hamed Jahani, and Jiuh-Biing Sheu. Organ transportation and allocation problem under medical uncertainty: A real case study of liver transplantation. Transportation Research Part E: Logistics and Transportation Review, 134:101841, 2020.

22 Michael Kearns, Yishay Mansour, Dana Ron, Ronitt Rubinfeld, Robert E Schapire, and Linda Sellie. On the learnability of discrete distributions. In Proceedings of the twenty-sixth annual ACM symposium on Theory of computing, pages 273–282, 1994.

23 Lin Kemeng, Wang Xiaoyan, Xia Weijie, Zhang Jiaming, et al. Optimization of the randomness in einstein which based on monte carlo algorithms. In 2019 Chinese Control And Decision Conference (CCDC), pages 6305–6309. IEEE, 2019.

24 Stefan Kern, Sibylle D Müller, Nikolaus Hansen, Dirk Büche, Jiri Ocenasek, and Petros Koumoutsakos. Learning probability distributions in continuous evolutionary algorithms-a comparative review. Natural Computing, 3:77–112, 2004.

25 Minh Thanh Khong, Christophe Lecoutre, Pierre Schaus, and Yves Deville. Soft-regular with a prefix-size violation measure. In International Conference on the Integration of Constraint Programming, Artificial Intelligence, and Operations Research, pages 333–343. Springer, 2018.

26 Daphne Koller and Nir Friedman. Probabilistic graphical models: principles and techniques. MIT press, 2009.
27 Anna LD Latour, Behrouz Babaki, Daniël Fokkinga, Marie Anastacio, Holger H Hoos, and Siegfried Nijssen. Exact stochastic constraint optimisation with applications in network analysis. *Artificial Intelligence*, 304:103650, 2022.

28 Pu Li, Harvey Arellano-Garcia, and Günter Wozny. Chance constrained programming approach to process optimization under uncertainty. *Computers & chemical engineering*, 32(1-2):25-45, 2008.

29 Xiaoxia Lin, Stacy L Janak, and Christodoulos A Floudas. A new robust optimization approach for scheduling under uncertainty:: I. bounded uncertainty. *Computers & chemical engineering*, 28(6-7):1069-1085, 2004.

30 Alexandre Mercier-Aubin, Ludwig Dumetz, Jonathan Gaudreault, and Claude-Guy Quimper. The confidence constraint: A step towards stochastic cp solvers. In *International Conference on Principles and Practice of Constraint Programming*, pages 759-773. Springer, 2020.

31 Arkadi Nemirovski and Alexander Shapiro. Convex approximations of chance constrained programs. *SIAM Journal on Optimization*, 17(4):969-996, 2007.

32 François Pachet, Pierre Roy, Alexandre Papadopoulos, and Jason Sakellariou. Generating 1/f noise sequences as constraint satisfaction: The voss constraint. In *Twenty-Fourth International Joint Conference on Artificial Intelligence*, 2015.

33 Bernardo K Pagnoncelli, Shabbir Ahmed, and Alexander Shapiro. Sample average approximation method for chance constrained programming: theory and applications. *Journal of optimization theory and applications*, 142(2):399-416, 2009.

34 M Arenas Parra, A Bilbao Terol, B Pérez Gladish, and MV Rodríguez Uría. Solving a multiobjective possibilistic problem through compromise programming. *European Journal of Operational Research*, 164(3):748-759, 2005.

35 Judea Pearl. *Probabilistic reasoning in intelligent systems: networks of plausible inference*. Morgan kaufmann, 1988.

36 Guillaume Perez, Steve Malalel, Gael Glorian, Victor Jung, Alexandre Papadopoulos, Marie Pelleau, Wijnand Suijlen, Jean-Charles Régin, and Arnaud Lallouet. Generalized confidence constraints. In *Proceedings of the AAAI Conference on Artificial Intelligence*, 2023.

37 Guillaume Perez, Brendan Rappazzo, and Carla Gomes. Extending the capacity of 1/f noise generation. In *International Conference on Principles and Practice of Constraint Programming*, pages 601-610. Springer, 2018.

38 Guillaume Perez and Jean-Charles Régin. MDDs are efficient modeling tools: An application to dispersion constraints. In *Integration of AI and OR Techniques in Constraint Programming*, 2017.

39 Guillaume Perez and Jean-Charles Régin. MddS: Sampling and probability constraints. In *International Conference on Principles and Practice of Constraint Programming*, pages 226-242. Springer, 2017.

40 Gilles Pesant. Achieving domain consistency and counting solutions for dispersion constraints. *INFORMS Journal on Computing*, 27(4):690-703, 2015. doi:10.1287/ijoc.2015.0654.

41 Warren B Powell. A unified framework for stochastic optimization. *European Journal of Operational Research*, 275(3):795-821, 2019.

42 Steven D Prestwich, Roberto Rossi, and S Armagan Tarim. Randomness as a constraint. In *Principles and Practice of Constraint Programming: 21st International Conference, CP 2015, Cork, Ireland, August 31-September 4, 2015, Proceedings 21*, pages 351-366. Springer, 2015.

43 Charles Prud'homme, Jean-Guillaume Fages, and Xavier Lorca. Choco solver documentation. *TASC, INRIA Rennes, LINA CNRS UMR*, 6241:13-42, 2016.

44 Jean-Charles Régin. Arc consistency for global cardinality constraints with costs. In *International Conference on Principles and Practice of Constraint Programming*, pages 390-404. Springer, 1999.

45 Francesca Rossi, Peter Van Beek, and Toby Walsh. *Handbook of constraint programming*. Elsevier, 2006.
46 Roberto Rossi, Brahim Hnich, S Armagan Tarim, and Steven Prestwich. Confidence-based reasoning in stochastic constraint programming. *Artificial Intelligence*, 228:129-152, 2015.

47 Christian Schulte and Guido Tack. View-based propagator derivation. *Constraints*, 18(1):75-107, 2013.

48 Kalika Suksomboon, Nobutaka Matsumoto, Shuichi Okamoto, Michiaki Hayashi, and Yusheng Ji. Configuring a software router by the erlang-k-based packet latency prediction. *IEEE Journal on Selected Areas in Communications*, 36(3):422-437, 2018.

49 S Armagan Tarim, Suresh Manandhar, and Toby Walsh. Stochastic constraint programming: A scenario-based approach. *Constraints*, 11:53-80, 2006.

50 Joseph A Tatman and Ross D Shachter. Dynamic programming and influence diagrams. *IEEE transactions on systems, man, and cybernetics*, 20(2):365-379, 1990.

51 Pascal Van Hentenryck and Laurent Michel. Domain views for constraint programming. In *International Conference on Principles and Practice of Constraint Programming*, pages 705-720. Springer, 2014.

52 Willem-Jan Van Hoeve, Gilles Pesant, and Louis-Martin Rousseau. On global warming: Flow-based soft global constraints. *Journal of Heuristics*, 12(4):347-373, 2006.

53 Toby Walsh. Stochastic constraint programming. In *ECAI*, volume 2, pages 111-115, 2002.

54 Weijun Xie and Shabbir Ahmed. On deterministic reformulations of distributionally robust joint chance constrained optimization problems. *SIAM Journal on Optimization*, 28(2):1151-1182, 2018.

55 Hui Zhang and Pu Li. Chance constrained programming for optimal power flow under uncertainty. *IEEE Transactions on Power Systems*, 26(4):2417-2424, 2011.
<figure><img src="image_291.png" /></figure>
# The p-Dispersion Problem with Distance Constraints

Nikolaos Ploskas ⬒

University of Western Macedonia, Kozani, Greece

Kostas Stergiou ⬒

University of Western Macedonia, Kozani, Greece

Dimosthenis C. Tsouros ⬒

KU Leuven, Belgium

## Abstract

In the (maxmin) p-dispersion problem we seek to locate a set of facilities in an area so that the minimum distance between any pair of facilities is maximized. We study a variant of this problem where there exist constraints specifying the minimum allowed distances between the facilities. This type of problem, which we call PDDP, has not received much attention within the literature on location and dispersion problems, despite its relevance to real scenarios. We propose both ILP and CP methods to solve the PDDP. Regarding ILP, we give two formulations derived from a classic and a state-of-the-art model for p-dispersion, respectively. Regarding CP, we first give a generic model that can be implemented within any standard CP solver, and we then propose a specialized heuristic Branch&Bound method. Experiments demonstrate that the ILP formulations are more efficient than the CP model, as the latter is unable to prove optimality in reasonable time, except for small problems, and is usually slower in finding solutions of the same quality than the ILP models. However, although the ILP approach displays good performance on small to medium size problems, it cannot efficiently handle larger ones. The heuristic CP-based method can be very efficient on larger problems and is able to quickly discover solutions to problems that are very hard for an ILP solver.

**2012 ACM Subject Classification** Theory of computation → Constraint and logic programming

**Keywords and phrases** Facility location, distance constraints, optimization

Digital Object Identifier 10.4230/LIPIcs.CP.2023.30

## 1 Introduction

Maximum diversity problems arise in many practical settings from facility location to telecommunications and social network analysis [28]. Arguably, the most famous such problem is the (maxmin) *p-dispersion problem* (PDP) [29]. In the PDP we are given a set of candidate locations $P = \{1, 2, ..., n\}$ for p facilities and an $n \times n$ matrix $(D[i, j])$, $i, j \in P$ with distances between candidate locations *i* and *j*. The goal is to select p items from P to locate the facilities such that the minimum distance between any pair of facilities is maximized.

In practice, the PDP occurs whenever a close proximity of facilities is dangerous or for other reasons undesirable. A standard application is concerned with the location of power plants, where we wish to minimize the risk of losing multiple plants in the event of an accident or an enemy attack. To achieve this, locations for the plants are desired so that interplant distances are as large as possible. Similar applications arise in the military sector, as it is common to scatter military installations in order to make it difficult for the enemy to disarm all of them. In a more peaceful context, we may wish to disperse branches of the same franchise, so that mutual competition between similar shops is minimized, or public facilities which have overlapping areas of service, e.g., schools, hospitals, electoral districts,

<figure><img src="image_292.png" /></figure>

© Nikolaos Ploskas, Kostas Stergiou, and Dimosthenis C. Tsouros;

licensed under Creative Commons License CC-BY 4.0

29th International Conference on Principles and Practice of Constraint Programming (CP 2023).

Editor: Roland H. C. Yap; Article No. 30; pp. 30:1–30:18

<figure><img src="image_293.png" /></figure>

LIPICS

Leibniz International Proceedings in Informatics

Schloss Dagstuhl – Leibniz-Zentrum für Informatik, Dagstuhl Publishing, Germany
30:2 The p-Dispersion Problem with Distance Constraints

etc. [25, 16, 28]. In telecommunications, we may wish to disperse radio transceivers to service
cellular phones so as to minimize interference. Another, more recent, area of application
for the PDP, is if distances are not interpreted physically but as a measure of the diversity
between members of a group [33].

Another dispersion problem that has been deeply studied is the maxsum p-dispersion problem, called p-defence problem in [29], where we seek to locate p facilities so that the sum of distances between the facilities is maximized [25]. Although the maxsum variant also has many applications, it is recognized that the PDP is better suited to model situations where the close proximity of facilities must be definitely avoided (e.g. for safety reasons). This is because maximizing the sum of distances does not guarantee that no two facilities will be placed close to each other [25].

In this paper we are concerned with a variant of the PDP where hard constraints specifying minimum allowed distances exist between facilities. We call this problem *p-dispersion with distance constraints* (PDDP). Although this problem was put forward by Moon and Chaudhry who fist systematically studied location problems with distance constraints and coined the term p-dispersion [29], it has been rather ignored since, despite its relevance to real applications. Distance constraints in dispersion problems can stem from operational needs and regulations, such as clearance distances for safe chemical storage [1], separation distances between packages containing radioactive materials [40], or between portable fire extinguishers in an area [41]. Motivated by such applications, a recent study by Dai et al. considers p-dispersion with distance constraints in the context of circle placement in non-convex polygons [10].

To further motivate our study of the PDDP, consider a scenario where p identical power plants need to be located in an area. Assume that due to safety reasons, any two plants must be more than x km away from one another. If this problem is modeled and solved as a PDP then either of the following two results will hold: 1) The optimal solution places the two closest facilities $y \le x$ km apart. Then the original problem is infeasible, as the safety requirements cannot be satisfied for all pairs of facilities, 2) The optimal solution places the two closest facilities $y > x$ km apart. Then, as y is the minimum distance between any two plants, all the safety requirements are satisfied and the original problem has been solved. But what if all power plants to be located are not identical? What if there are differences in the plants' sizes, the volume of power generated, the waste produced, etc.?

In such cases, the safety requirements regarding the minimum allowed distances between
plants may not be the same for all pairs of plants. For instance, the allowed distance between
smaller and less dangerous plants will probably be smaller than between larger ones. Hence,
the PDP does no longer suffice to model the problem. This is because an optimal solution
that places the closest plants y km apart does not guarantee that the safety distances will
be satisfied for all pairs of plants. The case of non-identical (heterogeneous) facilities has
not received much attention in the p-dispersion literature, as the predominant explicit or
implicit assumption is that the facilities to be located are indistinguishable (homogeneous).
But in practice, not all power plants will be identical, and the same holds for the branches of
a franchise, for public facilities, and almost any type of facilities that we want to disperse.

We start our study of the PDDP by giving two ILP formulations. The first one is based
on the classic formulation for p-dispersion by Kuby [25], while the second is based on a
state-of-the-art model proposed by Sayah and Irnich [34]. Both these formulations are
extended to deal with heterogeneous facilities and to include distance constraints.

Then, we describe a generic CP model for the PDDP that can be implemented within
any standard CP solver. For the purposes of this study, we have implemented this model
in OR-Tools and Choco. Experimental results demonstrate that the ILP formulations,
implemented in Gurobi, are more efficient than the CP solvers, as the latter find it very hard to prove optimality, even for small problems, and are usually slower in locating solutions of the same quality as the ILP solver.

We further explore the applicability of CP technology by introducing a specialized heuristic method based on CP. Specifically, using a simpler model of the problem, with fewer variables and constraints, we have devised a specialized Branch&Bound mechanism, which has been implemented in a custom CP solver. This method tries to prune the search tree early by estimating the best cost that can be achieved if the sub-tree rooted at the currently visited node is explored. If the estimated cost does not improve the cost of the best solution found so far then the current branch is abandoned. The estimation is carried out through a simple greedy assignment of the remaining variables. This reasoning achieves significant search tree pruning, albeit by sacrificing completeness.

In our experimental analysis we first compare the ILP formulations to the CP approaches on randomly generated problems with 5-30 facilities and at most 80 potential location points. Results show that the CP model implemented in the standard CP solvers OR-Tools and Choco cannot compete with the ILP formulations implemented in Gurobi, as the latter solver is quite efficient in all but the largest class that contains 30 facilities and 80 location points. The heuristic CP approach very quickly finds solutions, often optimal ones, on all instances, including those of the hardest class. We then consider harder problems that are generated using the p-dispersion MDPLIB benchmark library as basis [28]. Results demonstrate that the ILP formulations are efficient on problems with 10 facilities and 100 potential locations, but fail to efficiently handle larger problems. On the other hand, the heuristic CP approach can trivially find solutions of good quality on smaller instances, while it can also handle larger instances that are very hard for the ILP solver, finding solutions of much better quality.

## 2 Related Work

The maxmin p-dispersion problem, which is $NP$-hard on general networks for an arbitrary $p$ [18], was originally mentioned by Shier, as far back as 1977 [36]. However, the term p-dispersion first appeared in the analysis of location problems with distance constraints by Moon and Chaudhry [29]. The first ILP solution was proposed by Kuby [25] while the first specialized algorithm was given by Erkut [15]. Kincaid proposed simulated annealing and tabu search methods [24], Ghosh proposed a multi-start heuristic [18] and Resende et al. applied the GRASP methodology to the maxmin problem [32]. More recently, Sayyady & Fathi [35] and Sayah & Irnich [34] proposed ILP approaches to the maxmin problem, which are able to solve large size problems, and as argued in the comprehensive review on OR methods for dispersion problems given in [28], they tend to make heuristic approaches obsolete, as they can handle problems of similar size.

Regarding distance constraints, Moon and Chaudhry were the first to systematically study location problems with distance constraints [29]. The p-dispersion problem with distance constraints was mentioned by them as a problem that can arise in real-life scenarios, but no approaches towards solving it were proposed. Recently, Dai et al. revisited this problem as part of a study on circle (i.e. facility) dispersion in non-convex polygons [10]. A heuristic method, inspired by the mechanics of the n-body problem in physics, was proposed for the plain p-dispersion problem in non-convex polygons, and this method was also adapted to the case where distance constraints exist between pairs of circles.

Distance constraints have also been considered in the context of other location problems. Some early works considered maximum distance constraints between the demand nodes and the facility locations [7, 8, 23, 38]. Tansel et al. studied the distance constrained *p-center*
30:4 The p-Dispersion Problem with Distance Constraints

problem for the case where the network is a tree [38]. Chaudhry et al. proposed heuristics
for selecting a maximum-weight set of locations such that no two are closer than a given
distance from each other [6]. Moon and Papayanopoulos considered the problem of locating
two facilities so as to minimize the maximum of combined Euclidean distances to unweighted
existing points when the facilities must be separated by at least a specified distance [30].
Comley studied the problem of locating a small number of heterogeneous semi-obnoxious
facilities that interact with each other as well as with other existing facilities [9].

Berman and Huang studied the problem of locating homogeneous obnoxious facilities on a network so as to minimize the total demand covered, subject to the condition that no two facilities are allowed to be closer than a pre-specified distance [2]. Drezner et al. proposed the Weber obnoxious facility location problem where we seek to locate one facility so that the weighted sum of distances between the facility and demand points is minimized, with the additional requirement that the facility location is at least a given distance from demand points because it is obnoxious to them [13]. Drezner et al. considered a continuous multiple obnoxious facility location problem where a given number of facilities must be located in a convex polygon with the objective of maximizing the minimum distance between facilities and a given set of communities subject to distance constraints between facilities [14]. Welch and Salhi studied the location of obnoxious facilities with interactions between them [39]. Location problems with distance constraints that restrict the placement of facilities near certain demand points have also been studied, e.g. [31, 4, 27].

There are very few CP-related methods for facility location problems [17, 5, 37] and none of them concerns p-dispersion problems, with or without distance constraints. Regarding our heuristic CP-based method, there are works that follow a similar approach, i.e. sacrificing the completeness of a CP solver to solve optimization problems faster [22, 26, 19]. However, these works typically do this through a more local reasoning, e.g. by adding extra constraints.

Finally, the p-median problem with distance constraints, originally put forward in [29], is
being studied in a paper that is currently under review (details are suppressed to preserve
anonymity). In such a problem, there exist both facilities and clients that are serviced by
the facilities. The goal is to locate p facilities so that the sum of the distances between the
clients and their closest facility is minimized. As here, ILP and CP models for this problem
are proposed and compared. Results demonstrate that the ILP approach is by far the most
efficient on problems with homogeneous facilities, but it is outperformed by a heuristic CP
approach on some classes of problems with heterogeneous facilities.

3 Background

In a p-dispersion with distance constraints problem (PDDP), *p* facilities in a set of facilities *F* are to be placed in an area. We assume that the set *P* of potential location points for the facilities is known. We also assume that the distance between each pair of potential locations (*i*, *j*) is given in a matrix *D*. Between each pair of facilities *f*<sub>*i*</sub> and *f*<sub>*j*</sub> there is a distance constraint *dis*(*f*<sub>*i*</sub>, *f*<sub>*j*</sub>) > *d*<sub>*ij*</sub> specifying that the distance *dis*(*f*<sub>*i*</sub>, *f*<sub>*j*</sub>) between the points where the facilities *f*<sub>*i*</sub> and *f*<sub>*j*</sub> are located must be greater than *d*<sub>*ij*</sub>, where *d*<sub>*ij*</sub> is a constant. To summarize, in a PDDP we have:

- *P*: the set of candidate facility locations.

- *F*: the set of facilities to be located.

- *p*: the number of facilities to be located.

- *D[i,j]: the distance between any two candidate facility locations.*

- *d*<sub>*ij*</sub>: the lower bound in the allowed distance between each pair of facilities (*i*, *j*), where *i*, *j* ∈ *F*.
The distance between two points *i* and *j* can be given by the straight-line (Euclidean) distance, e.g. for the location of hazardous facilities, or by the shortest path in a street network, e.g. for the location of franchises, or by any other suitable metric. The methods we propose do not depend on any particular distance measure because, as is common in the literature, we assume that the pairwise distances between all possible location points are given in a 2-d distance matrix *D*. However, if necessary, instead of precomputing the distances and storing them in a distance matrix, they could be computed "on the fly", under the condition that this operation takes constant time. This holds for Euclidean or Manhattan distances given the coordinates of the points, but it does not hold for the shortest path in a network.

A common assumption in the literature on location problems with distance constraints is that the lower distance bound $d_{ij}$ is fixed to a specific value for all the constraints between facilities. This is a reasonable assumption in the case where the facilities are homogeneous, and therefore in essence indistinguishable, but it is not always realistic, especially when the facilities have different properties, as for example in [1, 40]. In case the facilities are heterogeneous, the distance bound may vary from constraint to constraint.

The goal in a PDDP is to locate each facility to a node so that the minimum distance between any two facilities is maximized subject to the satisfaction of all the distance constraints.

# 4 ILP models

We first give an ILP model for the PDDP, based on the classic formulation of Kuby for p-dispersion [25] and then we give a model based on the state-of-the-art model of Sayah & Irnich [34]. Both formulations are extended to deal with heterogeneous facilities and to include distance constraints between facilities.

## 4.1 Kuby based model

We make use of the following additional notation:

- $C = \{(i, j, f_1, f_2) \mid i, j \in P, f_1, f_2 \in F, D[i, j] \le d_{f_1 f_2}\}, \forall i \in P, \forall j \in P$, and for each pair of facilities $(f_1, f_2)$: the set of quadruples $(i, j, f_1, f_2)$ s.t. facilities $f_1$ and $f_2$ cannot be placed in facility sites $i$ and $j$, respectively, because $i$ and $j$ are not in a safe distance between each other with respect to the allowed distance between $f_1$ and $f_2$.

- $x_{ij} = 1$ if a facility $j \in F$ is located at a facility site $i \in P$ and 0 otherwise.

- $y_i = 1$ if any facility is located at a facility site $i \in P$ and 0 otherwise.

- $b$ is the minimum distance between the facilities that we aim to maximize.

For any $i \in P$, variable $y_i$ shows whether or not facility site $i$ will host any facility. These are the variables that are present in Kuby's formulation for p-dispersion. Given that facilities are considered identical in the p-dispersion literature, in Kuby's model we only need to know whether a site will host a facility or not. But in the case of the PDDP, where facilities can be different and distance constraints exist between them, we also need to know which particular facility will be hosted by a site. Hence, we introduce $|P| \times |F|$ variables, i.e. one variable $x_{ij}, \forall (i, j), i \in P, j \in F$, in order to know whether or not a specific facility $j \in F$ is located at a facility site $i \in P$.
30:6 The p-Dispersion Problem with Distance Constraints

The extension of Kuby's model to capture the PDDP can be expressed as:

$$
\max \qquad b \tag{1}
$$

$$
\text{s.t.} \quad \sum_{i \in P} y_i = p \tag{2}
$$

$$
b \le M (2 - y_i - y_j) + D[i, j] \quad \forall i, j \in P, j > i \quad (3)
$$

$$
\sum_{j \in F} x_{ij} = y_i \qquad \forall i \in P \tag{4}
$$

$$
\sum_{j \in F} x_{ij} \le 1 \qquad \forall i \in P \tag{5}
$$

$$
\sum_{i \in P} x_{ij} = 1 \qquad \forall j \in F \qquad (6)
$$

$$
x_{if_1} + x_{jf_2} \le 1 \qquad \forall (i, j, f_1, f_2) \in C \tag{7}
$$

$$
x_{ij} \in \{0, 1\} \qquad \forall i \in P, \forall j \in F \tag{8}
$$

$$
y_i \in \{0, 1\} \qquad \forall i \in P \tag{9}
$$

$$
b \ge 0 \tag{10}
$$

The objective function 1 aims at maximizing the shortest distance *b* between located
facilities. Constraint 2 specifies that *p* facilities are to be located. Constraint 3 guarantees
that *b* ≤ *D*[*i*, *j*] whenever both locations *i* and *j* are chosen via *y*<sub>*i*</sub> = *y*<sub>*j*</sub> = 1 for the location
of facilities, where *M* represents a Big *M* constant. Variables *b* and *y*<sub>*i*</sub>, *i* ∈ *P*, the objective
function 1 and Constraints 2, 3, 9 and 10 form the original formulation of Kuby for the
p-dispersion problem.

As we explained, in the case of the PDDP, we add variables $x_{ij}$, $\forall(i, j)$, $i \in P, j \in F$, in order to know whether or not a specific facility $j \in F$ is located at a facility site $i \in P$. These variables are linked to the $y_i$ variables via Constraint 4, which specifies that if any facility $j$ is located at a facility site $i$, then variable $y_i$ equals 1 and 0 otherwise. To ensure that no two variables $x_{ij}$ and $x_{ij'}$ are set to 1 (i.e. each facility site must host at most one facility), we add Constraint 5. To ensure that no two variables $x_{ij}$ and $x_{i'j}$ are set to 1 (i.e. each facility must be hosted at exactly one facility site), we add Constraint 6. Finally, Constraint 7 models the distance constraints between facilities. It ensures that each facility is at a safe distance from all other facilities by not allowing two facilities $f_1$ and $f_2$ to be established at sites that are at a distance closer than the allowed distance between $f_1$ and $f_2$. These pairwise constraints are a special case of clique constraints and are an efficient option to model distance constraints in ILP, as demonstrated in [2].

The original Kuby model for the p-dispersion problem has |*P*|+1 variables and Σ|*P*|−1 i=1 i+1 constraints, while our extended Kuby based model for the PDDP has |*P*| × |*F*| + |*P*| + 1 variables and 2 × |*P*| + |*F*| + Σ|*P*|−1 i=1 i + 1 constraints, without considering Constraint 7 which can give (|*P*| × |*F*|)² constraints.

4.2 Sayah and Irnich based model

We now present a model for the PDDP based on the PDP model proposed by Sayah and Irnich, which utilizes the fact that the optimal distance is equal to at least one of the entries of the distance matrix [34]. Let us introduce some additional notation for this model:

- $E = \{(i, j) \in P \times P : i < j\}$: the set of edges between any two candidate facility locations.

- $E(l) = \{(i, j) \in E : D[i, j] < l\}$: a subset of edges given any distance $l$.

- $L^0 < L^1 < \dots < L^{k_{max}}$: the different nonzero values in $D[i, j]$. The associate index sets are $K = \{1, 2, \dots, k_{max}\}$ and $K_0 = \{0\} \cup K$. By definition, $\emptyset = E(L^0) \subsetneq E(L^1) \subsetneq \dots \subsetneq E(L^{k_{max}}) \subsetneq E$ holds.

- $z_k = 1$ if the location decisions satisfy a minimum distance of at least $L^k$, $k \in K$ and 0 otherwise.
Similar to the Kuby model, variable $y_i$, for any $i \in P$, shows whether or not facility site $i$ will host any facility. In addition, variable $z_k$, for any $k \in K$, indicates whether the location decisions satisfy a minimum distance of at least $L^k$. These are the variables that are present in the formulation of Sayah and Irnich for p-dispersion. As we explained in Kuby's model, we also need to introduce $|P| \times |F|$ variables, i.e. one variable $x_{ij}, \forall (i, j), i \in P, j \in F$, in order to know whether or not a specific facility $j \in F$ is located at a facility site $i \in P$.

The model for the PDDP using Sayah and Irnich's model as basis can be expressed as:

$$ \max \quad D^0 + \sum_{k \in K} (L^k - L^{k-1}) z_k \qquad (11) $$

$$ \text{s.t.} \quad \sum_{i \in P} y_i = p \qquad (12) $$

$$ z_k \le z_{k-1} \quad \forall k \in K, k > 1 \qquad (13) $$

$$ y_i + y_j + z_k \le 2 \quad \forall k \in K, (i, j) \in E(L^k) \setminus E(L^{k-1}) \qquad (14) $$

$$ \sum_{j \in F} x_{ij} = y_i \quad \forall i \in P \qquad (15) $$

$$ \sum_{j \in F} x_{ij} \le 1 \quad \forall i \in P \qquad (16) $$

$$ \sum_{i \in P} x_{ij} = 1 \quad \forall j \in F \qquad (17) $$

$$ x_{if_1} + x_{jf_2} \le 1 \quad \forall (i, j, f_1, f_2) \in C \qquad (18) $$

$$ x_{ij} \in \{0, 1\} \quad \forall i \in P, \forall j \in F \qquad (19) $$

$$ y_i \in \{0, 1\} \quad \forall i \in P \qquad (20) $$

$$ z_k \in \{0, 1\} \quad \forall k \in K \qquad (21) $$

The objective function 11 aims at maximizing the shortest distance between located facilities. Constraints 12, 15, 16, 17, 18, 19 and 20 are also present in the Kuby based formulation, while Constraints 13, 14 and 21 replace Constraints 3 and 10. Constraint 13 models the consistency between the $z_k$ variables in the sense that the $z_k$ variables are non-increasing in $k$, while Constraint 14 ensures that no pair $(i, j)$ of locations with distance $D[i, j] < L^k$ is chosen simultaneously. The consistency Constraint 13 specifies that any feasible solution fulfills that there exists a unique $k \in K_0$ with $z_1 = z_2 = \dots = z_k = 1$ and $z_{k+1} = z_{k+2} = \dots = z_{kmax} = 0$. Variables $y_i, i \in P$, and $z_k, k \in K$, the objective function 11 and Constraints 12, 13, 14, 20 and 21 form the original formulation of Sayah and Irnich for the p-dispersion problem.

The addition of variables $x_{ij}, \forall (i, j), i \in P, j \in F$, in the revised Sayah and Irnich model for the PDPP yields Constraints 15, 16, 17, and 19, as already explained in Kuby's model. In addition, Constraint 18 models the distance constraints between facilities.

The original model of Sayah and Irnich for the p-dispersion problem has $|P| + k_{max} - 1$ variables and $k_{max} + (k_{max} - 1) \times \sum_{i=1}^{|P|-1} i - 1$ constraints, while our extended Sayah and Irnich based model for the PDDP has $|P| \times |F| + |P| + k_{max} - 1$ variables and $2 \times |P| + |F| + k_{max} + (k_{max} - 1) \times \sum_{i=1}^{|P|-1} i - 1$ without considering Constraint 18 which can give $(|P| \times |F|)^2$ constraints.

# 5 CP approaches to the PDDP

We first give a generic CP model of the PDDP and we then we describe the mechanics of a specialized heuristic CP solver. The PDDP is modeled as a Constraint Optimization Problem (COP) $(X, Dom, C, O)$, where $X$ is the set of decision variables, $Dom$ is the set
30:8 The p-Dispersion Problem with Distance Constraints

of finite domains, *C* is the set of hard constraints and *O* is the optimization function. The
model is as follows:

1. For each facility $i \in F$ there is a finite domain variable $x_i$. These $p$ variables are the decision variables in the problem, meaning that $|X| = |F| = p$. The domain of each variable $x_i \in X$, denoted by $Dom(x_i)$, includes as values all the points where a facility can be located, i.e. $\forall x_i \in X: Dom(x_i) = P$.

2. $Y$ is a set of auxiliary variables, s.t. for each pair of variables $(x_i, x_j) \in X \times X \mid i < j$, there is a variable $y_{ij} \in Y$ and a constraint $y_{ij} = D[x_i, x_j]$. Hence, each $y_{ij} \in Y$ models the distance between $x_i$ and $x_j$. In CP solvers, this is implemented using the Element global constraint, i.e. $y_{ij} = Element(D, [x_i, x_j])$.

3. For each variable $y_{ij}$, there is a distance constraint $y_{ij} > d_{ij}$.

4. There is a variable $z$, s.t. $z = \min(Y)$.

5. The objective function is $O = \text{maximize}(z)$.

This model contains $p + p \times (p - 1)/2 + 1$ variables, with $p$ being decision variables and the rest auxiliary, and $p \times (p - 1) + 1$ constraints.

We also considered the use of an AllDifferent constraint, to speed up propagation.
However, the distance constraints already propagate the fact that facilities should be placed
at different locations, as they all have bound greater than 0. Experiments with and without
the AllDifferent constraint showed no noticeable difference in run times.

**5.1 A heuristic CP-based method**

We now propose a heuristic technique that tries to prune early the parts of the search tree that do not seem promising, i.e. it is unlikely that exploring them will improve the value of the optimization function. At each node of the search tree this method tries to estimate the best value of the optimization function that can be achieved if we explore the sub-tree rooted at that node. If this value is not better than the cost of the best solution found so far then the current branch is not further explored.

The proposed method can be embedded in a standard CP solver. However, this cannot be done at the modeling level by just specifying variables and posting constraints, because it requires writing specialized code within the solver and possibly a intervention in how the search process works. Naturally, the estimation of the bound at each node cannot always be precise (otherwise we would be able to trivially solve the PDDP), and therefore, a solver that employs this method will not be exact.

To demonstrate our heuristic method, we describe a simple CP solver, specialized for the PDDP, that applies it. This solver uses a simpler model of the problem, dropping the auxiliary variables and relevant constraints. Hence, we now have a model with only the $p$ decision variables and $p \times (p-1)/2$ distance constraints. The optimization function is handled within the solver in the following way: Whenever a new solution is found, its cost is computed so as to determine if this cost is better than the current best cost. If so, then the best cost found so far is updated. If such an update occurs, it will be propagated to the decision variables, as described below.

The heuristic pruning technique works as follows: The cost of the first feasible solution found is used as the initial lower bound denoting the cost of the best solution found so far. Thereafter, at each node of the search tree, an upper bound for the best possible solution under the current assignment is computed. This upper bound gives an estimation of the best possible cost that can be achieved if the sub-tree rooted at the specific node is explored. If this is not higher than the current best cost then the current branch of the search tree is abandoned and the search moves on. Each time a solution with a higher cost than the current lower bound is found, the lower bound is updated.
The computation of the upper bound estimation at each node is performed in a greedy fashion. Specifically, assuming that $x_i$ is the current variable, $(x_1 \leftarrow v_1), \dots, (x_{i-1} \leftarrow v_{i-1})$ is the assignment to past variables and $v_i$ is the value under consideration for $x_i$, we greedily compute the cost of the “best” assignment for the future variables $x_{i+1}, \dots, x_p$. That is, we visit these variables one by one, starting with $x_{i+1}$, and for each variable $x_j, i+1 \le j \le p$, and each value $v_j \in \text{Dom}(x_j)$, we find the minimum distance between $v_j$ and any assignment (location) among variables (facilities) $x_1, \dots, x_{j-1}$. The value that maximizes this distance is then (temporarily) assigned to $x_j$. This is repeated until all variables have been assigned. We then compute the cost of this complete assignment, which gives the upper bound, but in fact, this may be an underestimation of the real cost. If the computed cost is equal to or lower than the current lower bound then the assignment of $v_i$ to $x_i$ is undone and the current branch of the search tree is abandoned, albeit risking to prune the branch leading to the optimal solution. This process is depicted by Algorithm 2.

Algorithm PDDP\_CP\_Solver (Algorithm 1) gives a high level description of the entire solving process.

■ Algorithm 1 PDDP\_CP\_Solver(X, Dom, C, O).

if Propagate(X, Dom, C) = FALSE
    return NULL;
depth ← 1;
best_found ← 0;
select an unassigned variable $x_i$;
while depth ≥ 1
    if all values in Dom($x_i$) have been tried
        depth ← depth-1;
    else
        select a value $a \in \text{Dom}(x_i)$ that has not been tried;
        if depth = n
            cur_cost ← Compute_Solution_Cost(X, Dom, C)
        if cur_cost > best_found
            best_found ← cur_cost;
        else if Propagate(X, Dom, C, $x_i \leftarrow a$) = TRUE
            if best_found ≠ 0 AND Bound(X, Dom, C, $x_i \leftarrow a, best_found$) = TRUE
                depth ← depth+1;
                select an unassigned variable $x_i$;
        if best_found = 0 return NULL;
    return best_found;

Given a PDDP $(X, DOM, C, O)$, where $O$ is the optimization function of the PDDP, the algorithm starts by propagating the hard constraints in $C$, as a typical CP solver does. Function Propagate enforces arc consistency on the distance constraints. If no failure (empty domain) is detected then the algorithm initializes the depth to 1 and the best cost found (*best_found*) to 0 and commences the search by selecting a variable using a variable ordering heuristic. While the depth of search is greater than 0, denoting that the search space has not been exhaustively searched, a branching decision is made, i.e. a value is selected and assigned to the currently selected variable. If all the variables have been assigned (*depth = n*), meaning that a feasible solution has been found, the cost of this solution is computed and if this cost is higher than the best cost found so far then the latter is updated.

If not all variables have been assigned yet then Function Propagate is called to propagate the value assignment just made. If no failure occurs, the heuristic bounding mechanism described above is triggered by calling Function Bound (Algorithm 2), provided that at least one feasible solution has already been found. If this function succeeds, meaning that the
30:10 The p-Dispersion Problem with Distance Constraints

estimated cost is better than the best bound found so far then the algorithm moves forward by increasing the depth of search and selecting a new unassigned variable. On the other hand, if either propagation fails or the estimated bound is not better than the value of best_found then the current branch is abandoned and a new value for the current variable is selected.

■ Algorithm 2 Bound(X, Dom, C, xᵢ ← vᵢ, best_found).

for each xⱼ, i + 1 ≤ j ≤ p
  val ← dis ← -1;
  for each vⱼ ∈ Dom(xⱼ)
    temp ← shortest distance between vⱼ and any assigned variable x₁...xⱼ₋₁;
    if temp > dis
      dis ← temp;
      val ← vⱼ;
    xⱼ ← val;
  bound ← Compute_Solution_Cost(X, Dom, C);
  if bound > best_found return TRUE;
  return FALSE;

If the value *best_found* remains 0 upon termination then the algorithm has proved that the problem is infeasible and the solver returns NULL. Otherwise, the best cost found is returned. In the former case, the heuristic part of the algorithm (i.e. the bounding mechanism) will never be triggered, as no feasible solution will have been found. Hence, the search space will be systematically explored in a typical CP solver fashion until a backtrack to depth -1 occurs, proving that the problem is infeasible.

Function *Propagate* applies arc consistency on the distance constraints, taking into consideration the value of *best_found*, i.e. the best cost found so far, to perform extra pruning, if possible. This is done in typical CP fashion for binary constraints, i.e. using a queue to insert and then process variables that have their domain filtered. Specifically, if a variable *x*ᵢ is removed from the queue then for each variable *x*ⱼ constrained with *x*ᵢ, and each value *v*ⱼ ∈ *Dom*(*x*ⱼ), we check if there exists a value *v*ᵢ in *Dom*(*x*ᵢ) s.t. the two values satisfy the distance constraint between *x*ⱼ and *x*ᵢ and the distance between the two values is greater than *best_found*. In other words, for each possible location *v*ⱼ of *x*ⱼ we search for a location *v*ᵢ for *x*ᵢ s.t. by placing the two facilities at these locations, not only the relevant distance constraint is satisfied, but we can also improve the cost of the optimization function. If no such *v*ᵢ exists then by placing *x*ⱼ at *v*ⱼ there is no way to locate *x*ᵢ so that we can satisfy the distance constraint and improve the cost. Hence, *v*ⱼ can be deleted from *D*(*x*ⱼ), i.e. it can be removed from consideration as a potential location point for *x*ⱼ. If such a value deletion occurs then variable *x*ⱼ is inserted in the queue to propagate the deletion.

Finally, note that the pruning that can be achieved by taking into consideration the current best cost can also be achieved by a standard CP solver that employs the model described further above. Such a solver will typically add a constraint *z* > *best_found* once a new solution with better cost than the previously found solutions is located. The propagation of this constraint may result in the filtering of the *y*ᵢⱼ variables' domains, which in turn will be carried over to the decision variables through the distance constraints *y*ᵢⱼ = *D*[*x*ᵢ, *x*ⱼ].

6 Experimental Results

We experimented with PDDP instances generated in two different ways. The first is a simple method that randomly generates a PDDP with a desired number of facilities and location points. The second uses the p-dispersion benchmark library MDPLIB 2.0 [28] as a basis to
generate PDDPs. Computations were performed on an Intel i7 CPU 8700 with 16 GB of main memory, a clock of 3.2 GHz, an L1 cache of 348 KB, an L2 cache of 2 MB, and an L3 cache of 12 MB, running under CentOS 8.4.

The ILP models were solved using Gurobi 9.0.3 [21]. The exact CP model was written in the CPMpy modeling tool [20] and compiled into OR-Tools [12]. An implementation in Choco [11] was also tried. The heuristic CP method was implemented in a custom solver programmed in C. This solver essentially implements the *CP_Solve* algorithm (Algorithm 1) described above. Choco and the heuristic solver use the dom/wdeg heuristic for variable ordering [3] and lexicographic value ordering, while OR-Tools uses its default options. A time out of 3,600 seconds was set for each instance. We used only one thread on all solvers, to get a fair comparison.

The ILP models are stored in compressed sparse column format since the constraint matrix can sometimes be too large to be stored as a full array in memory. For example, the largest instance considered in the computational study has a constraint matrix of 5,425,061 rows, 105,038 columns and 11,212,542 nonzeros. Its compressed size is only 12MB, while its size as a full matrix is 530GB.

## 6.1 Random problems

For an initial evaluation of the proposed approaches to the PDDP, we ran experiments on problems where we try to locate $p \in \{5, 10, 20, 30\}$ facilities in a $10 \times 10$ grid, with $|P| \in \{30, 80\}$ potential location points selected randomly among the 100 total points. The distances between the points are computed using the Manhattan distance metric. For each distance constraint $dis(f_i, f_j) > d_{ij}$ between facilities $f_i$ and $f_j$, $d_{ij}$ was set to a random integer in the interval $[1, max]/2$, where max is the maximum Manhattan distance between any two potential location points. Ten instances were generated and solved for each combination of parameter values.

Table 1 compares the following: Our two ILP formulations implemented in Gurobi, with Gurobi_k denoting our first model, based on that of Kuby, and Gurobi_s denoting our second model, based on that of Sayah & Irnich, and the CP solvers OR-tools and Choco. For each class, in column $\Sigma$ cpu we give the total cpu time taken by the corresponding solver over all 10 instances, and in brackets we give the number of instances on which the solver timed out. If the solver timed out on at least one instance then $\Sigma$ cpu gives a lower bound on the actual run time. Column cpu_o gives the mean time taken by the solver to locate the optimal solution. A zero means that less than 0.1 seconds were taken on average. In brackets, we give the number of instances for which the solver found the optimal solution. Note that the optimal solutions are known for all instances because at least one solver terminated within the time limit. When a solver managed to find the optimal solution on at least 8 out of the 10 instances, we compute cpu_o, excluding the instances where it timed out. Otherwise, cpu_o is left blank (-), meaning that the solver failed to find the optimal solution on too many instances for this metric to be meaningful.

From Table 1 it is clear that the ILP approach is more efficient than the CP one in this type of instances. Gurobi_k (resp. Gurobi_s) times out on 1 (resp. 2) out of the 70 instances, whereas OR-Tools (resp. Choco) timed out on 40 (resp. 42) instances. With respect to the cases when the optimal solution was located within the time limit, even without proving optimality, Gurobi_k (resp. Gurobi_s) found the optimal solution on 69 (resp. 68 instances), whereas OR-Tools (resp. Choco) on 59 (resp. 45) instances. This indicates that the CP solvers find the proof of optimality especially difficult. Regarding the time required to find the optimal solution (when found), all solvers are quite fast on smaller problems (with 5-10
30:12 The p-Dispersion Problem with Distance Constraints

<table><caption>**Table 1** Comparing exact solvers on random PDDPs. Cpu times are given in seconds.</caption><thead><tr><th>Class<br>(p, |P|)</th><th>Gurobi<sub>k</sub><br>Σcpu</th><th>cpu<sub>o</sub></th><th>Gurobi<sub>s</sub><br>Σcpu</th><th>cpu<sub>o</sub></th><th>OR-Tools<br>Σcpu</th><th>cpu<sub>o</sub></th><th>Choco<br>Σcpu</th><th>cpu<sub>o</sub></th></tr></thead><tbody><tr><td>(5,30)</td><td>0.7 (0)</td><td>0 (10)</td><td>1.6 (0)</td><td>0 (10)</td><td>3 (0)</td><td>0 (10)</td><td>3.5 (0)</td><td>0 (10)</td></tr><tr><td>(10,30)</td><td>2 (0)</td><td>0 (10)</td><td>6 (0)</td><td>0 (10)</td><td>998 (0)</td><td>1.7 (10)</td><td>>11,517 (2)</td><td>50 (10)</td></tr><tr><td>(20,30)</td><td>87 (0)</td><td>7 (10)</td><td>51 (0)</td><td>4 (10)</td><td>>10h (10)</td><td>69 (9)</td><td>>10h (10)</td><td>- (4)</td></tr><tr><td>(5,80)</td><td>7 (0)</td><td>0 (10)</td><td>20 (0)</td><td>0 (10)</td><td>13 (0)</td><td>0 (10)</td><td>3 (0)</td><td>0 (10)</td></tr><tr><td>(10,80)</td><td>12 (0)</td><td>0 (10)</td><td>78 (0)</td><td>0 (10)</td><td>>10h (10)</td><td>0.3 (10)</td><td>>10h (10)</td><td>0 (10)</td></tr><tr><td>(20,80)</td><td>727 (0)</td><td>69 (10)</td><td>725 (0)</td><td>28 (10)</td><td>>10h (10)</td><td>- (4)</td><td>>10h (10)</td><td>- (1)</td></tr><tr><td>(30,80)</td><td>>22,207 (1)</td><td>2,065 (9)</td><td>>10,519 (2)</td><td>414 (8)</td><td>>10h (10)</td><td>- (6)</td><td>>10h (10)</td><td>- (1)</td></tr></tbody></table>

facilities). Problems with 20 facilities and 30 potential locations are also easily manageable by the ILP models, whereas OR-Tools takes more than one minute on average to discover the optimal solution, and even fails to discover it on one instance. Choco fairs even worse, failing to find the optimal solution on 6 instances, and failing to find any solution on one instance.

As the size of the problem grows, the ILP models, and even more so the CP solvers, find it harder to solve the instances. On problems with 30 facilities and 80 location points, the ILP models start giving time outs and take a long time to find the optimal solution, when they manage to do so, whereas the CP solvers, and especially Choco, fail to find the optimal on many instances from the (20,80) and (30,80) classes.

Table 2 compares the exact solvers to the heuristic CP solver (denoted $CP_h$). For this solver, we report the total cpu time taken over the 10 instances of each class ($\Sigma cpu$) and the mean cpu time taken to find the best solution it found within the time limit ($cpu_b$). We also report (in brackets) the number of instances in which the solver managed to find the optimal solution. Then, in the following columns, we give the mean cpu times taken by the exact solvers to find a solution that at least matches the cost of the best solution found by the heuristic solver. In this way, we can evaluate the worth of the heuristic CP method as a heuristic for the PDDP. If the exact solvers manage to quickly match the best solution found by the heuristic one then there is not much point in using the heuristic solver. Whereas, if the heuristic solver quickly discovers a solution that the exact ones take very long to match then it is worth considering this approach for the PDDP. If an exact solver only managed to find a solution as good as that found by the heuristic solver in some instances, we give in brackets the number of times that this occurred, and we consider only these instances for the computation of the mean cpu time.

<table><caption>**Table 2** Comparing solvers on random PDDPs. Cpu times are given in seconds.</caption><thead><tr><th>Class<br>(p, |P|)</th><th>CP<sub>h</sub><br>Σcpu</th><th>cpu<sub>b</sub></th><th>Gurobi<sub>k</sub></th><th>Gurobi<sub>s</sub></th><th>OR-Tools</th><th>Choco</th></tr></thead><tbody><tr><td>(5,30)</td><td>0</td><td>0 (9)</td><td>0</td><td>0</td><td>0</td><td>0</td></tr><tr><td>(10,30)</td><td>0</td><td>0 (10)</td><td>0</td><td>0</td><td>1.7</td><td>50</td></tr><tr><td>(20,30)</td><td>2</td><td>0.2 (9)</td><td>5</td><td>1.7</td><td>65</td><td>160 (4)</td></tr><tr><td>(5,80)</td><td>0.1</td><td>0 (10)</td><td>0</td><td>0</td><td>0</td><td>0</td></tr><tr><td>(10,80)</td><td>1</td><td>0.1 (10)</td><td>0</td><td>0</td><td>0.3</td><td>0</td></tr><tr><td>(20,80)</td><td>2</td><td>0.2 (0)</td><td>3.5</td><td>3</td><td>6</td><td>0</td></tr><tr><td>(30,80)</td><td>36</td><td>3 (10)</td><td>2,218 (9)</td><td>1,051 (8)</td><td>1,760 (6)</td><td>3,453 (1)</td></tr></tbody></table>

As the results in Table 2 demonstrate, $CP_h$ is very fast as it managed to complete all 10 instances of each class in at most 2 seconds, except for the (30,80) class, on which it only took 36 seconds. Importantly, it also discovered the optimal solution in 58 out of the 70
instances, including all instances of the hard class (30,80) class. On the other hand, it did not manage to discover the optimal in any instance of the (20,80) class. Comparing the run times (cpu<sub>b</sub>) of CP<sub>h</sub> to the exact solvers, results from all but the last class show that Gurobi manages to quickly match the best solution found by CP<sub>h</sub>, even if the latter is faster on average. This does not always hold for OR-Tools which takes more than 1 minute on average in the (20,30) class. Choco takes 50 seconds on average in the (10,30) class and it manages to match the best solution of CP<sub>h</sub> in only 4 instances of the (20,30) class. The benefits of the heuristic method in hard problems are demonstrated by the (30,80) class where it takes 3 secs on average to locate the optimal solution, whereas Gurobi<sub>k</sub> and Gurobi<sub>s</sub> require 2,218 and 1,051 secs respectively. OR-Tools found the optimal solution in only 6 out of the 10 instances of this class, taking 1,760 secs on average, while Choco managed to find the optimal in only one instance in 3,453 secs.

Finally, Table 3 takes a closer look at the performance of CP<sub>h</sub>, with respect to the effect of the heuristic pruning method. To investigate this, we report the results obtained by the solver when the heuristic is deactivated (i.e. Function *Bound* is not called). In this case, the solver, denoted as CP<sub>-h</sub>, operates as a typical CP solver. As in Table 1, we give the total cpu time taken, and in brackets the number of time outs, and the mean time taken to locate the optimal solution (the number of instances where the optimal solution was found is given in brackets). In the following column (cpu<sub>h</sub>) we give the mean cpu time taken by CP<sub>-h</sub> to find a solution that at least matches the cost of the best solution found by CP<sub>h</sub>, and in brackets, the number of times that it managed to do so. The next two columns give the mean numbers of visited search tree nodes for CP<sub>-h</sub> and CP<sub>h</sub> (the entry is left blank if there were time-outs). The last two columns give the average number of calls to Function *Bound* in CP<sub>h</sub> and the percentage of fails caused by this function (i.e. the percentage of branches pruned by the heuristic).

<table><caption>**Table 3** A closer look at the performance of the custom solver, with and without the heuristic.</caption><thead><tr><th>(p, |P|)</th><th>CP<sub>-h</sub> &Sigma; cpu</th><th>CP<sub>-h</sub> cpu<sub>o</sub></th><th>CP<sub>-h</sub> cpu<sub>h</sub></th><th>CP<sub>-h</sub> nodes</th><th>CP<sub>h</sub> nodes</th><th>calls</th><th>%fails</th></tr></thead><tbody><tr><td>(5,30)</td><td>0.2 (0)</td><td>0 (10)</td><td>0 (10)</td><td>1,726</td><td>156</td><td>140</td><td>94</td></tr><tr><td>(10,30)</td><td>466 (0)</td><td>5 (10)</td><td>5 (10)</td><td>12.5M</td><td>229</td><td>218</td><td>95</td></tr><tr><td>(20,30)</td><td>&gt;21,827 (3)</td><td>1 (10)</td><td>1 (10)</td><td>-</td><td>4,917</td><td>781</td><td>80</td></tr><tr><td>(5,80)</td><td>6 (0)</td><td>0.6 (10)</td><td>0.6 (10)</td><td>7,624</td><td>501</td><td>477</td><td>97</td></tr><tr><td>(10,80)</td><td>&gt;10h (10)</td><td>- (5)</td><td>520 (5)</td><td>-</td><td>877</td><td>868</td><td>98</td></tr><tr><td>(20,80)</td><td>&gt;10h (10)</td><td>- (0)</td><td>- (0)</td><td>-</td><td>835</td><td>816</td><td>98</td></tr><tr><td>(30,80)</td><td>&gt;10h (10)</td><td>- (7)</td><td>1,044 (7)</td><td>-</td><td>17,207</td><td>12,995</td><td>81</td></tr></tbody></table>

Table 3 demonstrates the pruning power of our proposed heuristic. Without its use, the solver is able to handle the easier classes of problems, but not the harder ones, in accordance with standard CP solvers (Table 1). The solver is very successful compared to OR-Tools and Choco on the (20,30) class, as it times out in only 3 instances and finds the optimal solution in 1 sec on average. However, the solver fares badly on the (20,80) class where it is unable to find the optimal in any instance and actually finds worse solutions that CP<sub>h</sub> in all instances. As for the (30,80) class, CP<sub>-h</sub> finds the optimal solution in 7 instances, which is better than OR-Tools and Choco, but needs 1,044 secs on average to match the solution found by CP<sub>h</sub>. Regarding the pruning achieved by the heuristic when it is activated, it is impressive that in most of the classes, there is a very high percentage of pruned branches over the total calls to the heuristic (up to 98%), which explains its success in speeding up search.
## 6.2 MDPLIB

The MDPLIB collects a large number of dispersion benchmarks (for both `maxmin` and `maxsum` p-dispersion) divided into various classes [28]. In the GKD, MDG, and SOM classes, the distances between the potential facility locations are given by Euclidean distances, random real numbers, and random integers, respectively. We took instances from these classes, having 100-250 potential facility points and 10-20 facilities, and we generated 10 PDDPs for each instance by randomly adding distance constraints between the facilities. For each distance constraint `dis(fᵢ, fⱼ) > dᵢⱼ` between facilities `fᵢ` and `fⱼ`, `dᵢⱼ` was set to a random number in the interval `[1, max]/2`, where `max` is the maximum distance between any two potential location points, as specified in the base MDPLIB instance.

Table 4 compares the exact solvers on the MDPLIB-based problems. We exclude Choco as it is clearly inferior to OR-Tools. For each class we give the number of the MDPLIB instance on which it is based, and in brackets the numbers of potential locations and facilities to be located, e.g. `a1(100,10)` for MDG. For each solver, we report the total cpu time taken over the 10 instances (`∑ cpu columns`), the number of times when the optimal solution was found (`#opt`) columns, and the mean of the optimization function's value for the best solution found within the time limit. In the `∑ cpu` columns we give in brackets the number of instances in which the solvers timed out. In the `#opt` columns we give in brackets the number of times when optimality was proved. Finally, in some cases, a solver did not manage to find any solution within the time limit. In such cases there is a subscript in the value of cost, giving the number of instances in which at least one solution was found. In such a case, the value of cost is computed over these instances only. In GKD classes with 250 points and 20 facilities, OR-Tools suffered memory exhaustion and crashed. This is denoted with `MEM` in the `∑ cpu` column.

The data in Table 4 demonstrates that the PDDPs generated using MDPLIB instances as basis can be very hard for both the ILP and CP approaches. None of the solvers terminates within the time limit on any instance with 20 facilities, while problems with 15 and 10 facilities are also quite hard. In addition, there are some instances of the larger classes (e.g. `GKDd1(250,20)`) where the solvers are unable to discover any solution within 1 hour of cpu time, let alone the optimal one. Comparing ILP to CP, Gurobi, with any of the two formulations, is in general more efficient than OR-Tools. $Gurobi_k$ (resp. $Gurobi_s$) found the optimal solution in 77 (resp. 74) out of the 220 instances, and proved optimality in 73 (resp. 70) instances, whereas OR-Tools did not prove optimality in any instance (as it timed out on all of them) and found the optimal in 10 instances only. However, OR-Tools often managed to find better solutions than Gurobi in hard classes with 20 facilities, as the cost columns indicate, for instance in classes `MDGa2(100,20)` and `MDGb2(100,20)`. Comparing $Gurobi_k$ to $Gurobi_s$, there is no clear winner in terms of run times, but the latter managed to find solutions of better quality than the former in most of the classes. However, $Gurobi_s$ proved optimality or found the optimal solution in slightly fewer instances than $Gurobi_k$.

Table 5 compares $CP_h$ to the exact solvers. For $CP_h$, we report the total cpu time it takes over the 10 instances of each class (and the number of time-outs in brackets), the mean time it takes to find its best solution (`cpu_b`), and the mean of the optimization function's value for the best solution it finds. For the exact solvers, we report the mean times taken to match the value of the best solution found by $CP_h$ (`cpu_h`) columns, and the number of instances on which the solvers managed to find a solution that matches or improves the best solution found by $CP_h$ (in brackets). If this number is 0 or close to 0 then the entry in the `cpu_h` column is left blank (-), as it is impossible or meaningless to compute the value of `cpu_h`. If a value in the `cpu_b` column is blank (-), e.g. `GKDd1(100,20)`, then most of the 10 instances in this class were infeasible (in brackets we give the number of infeasible instances).
<table><caption>**Table 4** Comparing exact solvers on MDPLIB-generated PDDPs. Cpu times are given in seconds. The best mean value of the optimization function for each class is denoted in bold.</caption><thead><tr><th>Class<br>(p,|P|)</th><th>Gurobi<sub>k</sub><br>Σcpu</th><th>#opt</th><th>cost</th><th>Gurobi<sub>s</sub><br>Σcpu</th><th>#opt</th><th>cost</th><th>OR-Tools<br>Σcpu</th><th>#opt</th><th>cost</th></tr></thead><tbody><tr><td>MDG</td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td></tr><tr><td>a1(100,10)</td><td>&gt;25,745 (3)</td><td>9 (7)</td><td>4.67</td><td>11,208 (0)</td><td>10 (10)</td><td><strong>4.68</strong></td><td>&gt;10h (10)</td><td>7 (0)</td><td>4.59</td></tr><tr><td>a1(100,20)</td><td>&gt;10h (10)</td><td>0 (0)</td><td>0.80<sub>s</sub></td><td>&gt;10h (10)</td><td>0 (0)</td><td>1.16</td><td>&gt;10h (10)</td><td>0 (0)</td><td><strong>1.17</strong></td></tr><tr><td>a2(100,10)</td><td>15,327 (0)</td><td>10 (9)</td><td><strong>4.74</strong></td><td>11,412 (0)</td><td>10 (10)</td><td><strong>4.74</strong></td><td>&gt;10h (10)</td><td>1 (0)</td><td>4.36</td></tr><tr><td>a2(100,20)</td><td>&gt;10h (10)</td><td>0 (0)</td><td>0.78<sub>s</sub></td><td>&gt;10h (10)</td><td>0 (0)</td><td>0.83<sub>9</sub></td><td>&gt;10h (10)</td><td>0 (0)</td><td><strong>1.40</strong></td></tr><tr><td>b1(100,10)</td><td>11,852 (0)</td><td>10 (10)</td><td><strong>460.11</strong></td><td>&gt;29,171 (2)</td><td>10 (8)</td><td><strong>460.11</strong></td><td>&gt;10h (10)</td><td>0 (0)</td><td>430.25</td></tr><tr><td>b1(100,20)</td><td>&gt;10h (10)</td><td>0 (0)</td><td>102.42</td><td>&gt;10h (10)</td><td>0 (0)</td><td><strong>109.35</strong></td><td>&gt;10h (10)</td><td>0 (0)</td><td>77.33</td></tr><tr><td>b2(100,10)</td><td>10,305 (2)</td><td>8 (8)</td><td>459.65</td><td>27,894 (2)</td><td>8 (8)</td><td><strong>459.99</strong></td><td>&gt;10h (10)</td><td>1 (0)</td><td>413.85</td></tr><tr><td>b2(100,20)</td><td>&gt;10h (10)</td><td>0 (0)</td><td>106.77<sub>9</sub></td><td>&gt;10h (10)</td><td>0 (0)</td><td>81.09</td><td>&gt;10h (10)</td><td>0 (0)</td><td><strong>113.33</strong></td></tr><tr><td>GKD</td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td></tr><tr><td>d1(100,10)</td><td>4,743 (0)</td><td>10 (10)</td><td><strong>34.06</strong></td><td>5,415 (0)</td><td>10 (10)</td><td><strong>34.06</strong></td><td>&gt;10h (10)</td><td>0 (0)</td><td>33.04</td></tr><tr><td>d1(100,20)</td><td>&gt;10h (10)</td><td>0 (0)</td><td>-</td><td>&gt;10h (10)</td><td>0 (0)</td><td>-</td><td>&gt;10h (10)</td><td>0 (0)</td><td>-</td></tr><tr><td>d1(250,10)</td><td>&gt;15,448 (4)</td><td>6 (6)</td><td><strong>35.98</strong></td><td>34,753 (9)</td><td>3 (1)</td><td>35.27</td><td>&gt;10h (10)</td><td>0 (0)</td><td>-</td></tr><tr><td>d1(250,20)</td><td>&gt;10h (10)</td><td>0 (0)</td><td>9.95<sub>4</sub></td><td>&gt;10h (10)</td><td>0 (0)</td><td>10.55<sub>2</sub></td><td>MEM</td><td>0 (0)</td><td>-</td></tr><tr><td>d2(100,10)</td><td>&gt;5,998 (1)</td><td>9 (9)</td><td>34.34</td><td>2,602</td><td>10 (10)</td><td><strong>34.82</strong></td><td>&gt;10h (10)</td><td>0 (0)</td><td>31.18</td></tr><tr><td>d2(100,20)</td><td>&gt;10h (10)</td><td>0 (0)</td><td>-</td><td>&gt;10h (10)</td><td>0 (0)</td><td>-</td><td>&gt;10h (10)</td><td>0 (0)</td><td>-</td></tr><tr><td>d2(250,10)</td><td>&gt;22,334 (6)</td><td>4 (4)</td><td>35.94</td><td>&gt;33,595 (8)</td><td>2 (2)</td><td><strong>36.31</strong></td><td>&gt;10h (10)</td><td>0 (0)</td><td>-</td></tr><tr><td>d2(250,20)</td><td>&gt;10h (10)</td><td>0 (0)</td><td>-</td><td>&gt;10h (10)</td><td>0 (0)</td><td>-</td><td>MEM</td><td>0 (0)</td><td>-</td></tr><tr><td>SOM</td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td></tr><tr><td>a21(100,10)</td><td>&gt;20,260 (1)</td><td>10 (9)</td><td><strong>5</strong></td><td>9,585 (0)</td><td>10 (10)</td><td><strong>5</strong></td><td>&gt;10h (10)</td><td>1 (0)</td><td>4.1</td></tr><tr><td>a21(100,15)</td><td>&gt;34,887 (9)</td><td>1 (1)</td><td><strong>2.2</strong></td><td>&gt;35,371</td><td>1 (1)</td><td><strong>2.2</strong></td><td>&gt;10h (10)</td><td>0 (0)</td><td>2</td></tr><tr><td>a21(100,20)</td><td>&gt;10h (10)</td><td>0 (0)</td><td>1<sub>6</sub></td><td>&gt;10h (10)</td><td>0 (0)</td><td>1<sub>6</sub></td><td>&gt;10h (10)</td><td>0 (0)</td><td><strong>1<sub>9</sub></strong></td></tr><tr><td>a41(150,15)</td><td>&gt;10h (10)</td><td>0 (0)</td><td>2.6</td><td>&gt;10h (10)</td><td>0 (0)</td><td>2.6</td><td>&gt;10h (10)</td><td>0 (0)</td><td><strong>3</strong></td></tr><tr><td>a41(150,20)</td><td>&gt;10h (10)</td><td>0 (0)</td><td>1<sub>9</sub></td><td>&gt;10h (10)</td><td>0 (0)</td><td>1.11<sub>9</sub></td><td>&gt;10h (10)</td><td>0 (0)</td><td>1</td></tr><tr><td>b5(200,20)</td><td>&gt;10h (10)</td><td>0 (0)</td><td>1.4</td><td>&gt;10h (10)</td><td>0 (0)</td><td><strong>2</strong></td><td>&gt;10h (10)</td><td>0 (0)</td><td>-</td></tr></tbody></table>

Results from Table 5 demonstrate the efficiency of the heuristic CP approach. Regarding run times, $CP_h$ times out in only 7 instances and terminates quickly in all instances of all classes, except for the hard GKDd1(250,20) and GKDd2(250,20). $CP_h$ proved the infeasibility of the 16 infeasible instances of classes GKDd1(100,20) and GKDd2(100,20) and found solutions in the other 4, whereas none of the other solvers managed to do so in any instance. Of course, the fast proof of infeasibility is not due to the bounding mechanism of $CP_h$, as this is not invoked, but most likely due to the lightweight model and mechanics of the custom-written solver.

Regarding the quality of the solutions, as a downside, $CP_h$ finds the optimal in only 2 out of the 80 instances of smaller size, for which the optimal is known (excluding infeasible ones). However, in these classes, the exact solvers (and especially OR-Tools) can be orders of magnitude slower than $CP_h$ in discovering solutions of the same quality as $CP_h$, which reaches its best solution in less than 0.1 secs in most cases. This is evident in class SOMa41(150,15) where $CP_h$ found its best solution in less than 0.1 secs on average, while the exact solvers took more than 1,000 secs to match the solution quality of $CP_h$, on instances where they managed to do this. However, the best solution discovered by these solvers (including OR-Tools) is typically better than the best solution discovered by $CP_h$.
30:16 The p-Dispersion Problem with Distance Constraints

<table><caption>Table 5 Comparing solvers on MDPLIB-generated PDDPs. We denote with bold the mean cost of CP<sub>h</sub> on classes where it was better than the mean cost of all the other solvers.</caption><thead><tr><th>Class<br>(p,|P|)</th><th>CP<sub>h</sub><br>Σcpu</th><th>cpu<sub>b</sub></th><th>cost</th><th>Gurobi<sub>k</sub><br>cpu<sub>h</sub></th><th>Gurobi<sub>s</sub><br>cpu<sub>h</sub></th><th>OR-Tools<br>cpu<sub>h</sub></th></tr></thead><tbody><tr><td>MDG</td><td></td><td></td><td></td><td></td><td></td><td></td></tr><tr><td>a1 (100,10)</td><td>0.5 (0)</td><td>0</td><td>4.35</td><td>243 (10)</td><td>55 (10)</td><td>592 (10)</td></tr><tr><td>a1 (100,20)</td><td>81 (0)</td><td>8</td><td><strong>1.69</strong></td><td>- (0)</td><td>- (1)</td><td>- (0)</td></tr><tr><td>a2 (100,10)</td><td>0.6 (0)</td><td>0</td><td>4.24</td><td>137 (10)</td><td>22 (10)</td><td>873 (7)</td></tr><tr><td>a2 (100,20)</td><td>137 (0)</td><td>12</td><td><strong>1.64</strong></td><td>- (0)</td><td>- (0)</td><td>- (0)</td></tr><tr><td>b1 (100,10)</td><td>0.4 (0)</td><td>0</td><td>428.18</td><td>10 (10)</td><td>345 (10)</td><td>187 (10)</td></tr><tr><td>b1 (100,20)</td><td>54 (0)</td><td>4.5</td><td><strong>181.20</strong></td><td>- (0)</td><td>- (0)</td><td>- (1)</td></tr><tr><td>b2 (100,10)</td><td>0.7 (0)</td><td>0</td><td>428.17</td><td>27 (10)</td><td>41 (10)</td><td>723 (3)</td></tr><tr><td>b2 (100,20)</td><td>92 (0)</td><td>8</td><td><strong>159.93</strong></td><td>- (0)</td><td>- (0)</td><td>- (1)</td></tr><tr><td>GKD</td><td></td><td></td><td></td><td></td><td></td><td></td></tr><tr><td>d1 (100,10)</td><td>1.8 (0)</td><td>0.1</td><td>33.27</td><td>94 (10)</td><td>224 (10)</td><td>717 (5)</td></tr><tr><td>d1 (100,20)</td><td>405 (0)</td><td>- (9)</td><td>-</td><td>- (0)</td><td>- (0)</td><td>- (0)</td></tr><tr><td>d1 (250,10)</td><td>10 (0)</td><td>0.6</td><td>34.24</td><td>800 (8)</td><td>2,074 (6)</td><td>- (0)</td></tr><tr><td>d1 (250,20)</td><td>>18,689 (3)</td><td>1,140</td><td><strong>16.94</strong></td><td>- (0)</td><td>- (0)</td><td>- (0)</td></tr><tr><td>d2 (100,10)</td><td>1.6 (0)</td><td>0</td><td>31.29</td><td>91 (9)</td><td>179 (10)</td><td>1,277 (6)</td></tr><tr><td>d2 (100,20)</td><td>962 (0)</td><td>- (7)</td><td>-</td><td>- (0)</td><td>- (0)</td><td>- (0)</td></tr><tr><td>d2 (250,10)</td><td>8 (0)</td><td>0.5</td><td>35.06</td><td>121 (6)</td><td>997 (8)</td><td>- (0)</td></tr><tr><td>d2 (250,20)</td><td>>23,767 (4)</td><td>1,597</td><td><strong>13.28</strong></td><td>- (0)</td><td>- (0)</td><td>- (0)</td></tr><tr><td>SOM</td><td></td><td></td><td></td><td></td><td></td><td></td></tr><tr><td>a21 (100,10)</td><td>0 (0)</td><td>0</td><td>4</td><td>7 (10)</td><td>3 (10)</td><td>106 (10)</td></tr><tr><td>a21 (100,15)</td><td>0 (0)</td><td>0</td><td>2</td><td>39 (10)</td><td>45 (10)</td><td>206 (10)</td></tr><tr><td>a21 (100,20)</td><td>11 (0)</td><td>1</td><td>1</td><td>41 (6)</td><td>49 (6)</td><td>522 (9)</td></tr><tr><td>a41 (150,15)</td><td>2 (0)</td><td>0</td><td><strong>3</strong></td><td>1,170 (6)</td><td>1,741 (6)</td><td>1,465 (10)</td></tr><tr><td>a41 (150,20)</td><td>12 (0)</td><td>1</td><td><strong>1.9</strong></td><td>- (0)</td><td>- (1)</td><td>- (1)</td></tr><tr><td>b5 (200,20)</td><td>10 (0)</td><td>0.3</td><td>2</td><td>- (4)</td><td>1,272 (10)</td><td>- (0)</td></tr></tbody></table>

But the power of CP<sub>h</sub> as a heuristic method for PDDP is truly evident on the larger classes with 20 facilities where it discovers solutions of (much) better quality than the exact solvers, and excluding the two hard GKD classes with 250 location points and 20 facilities, it does this very fast. Also, CP<sub>h</sub> finds solutions in all 20 instances of the two hard GKD classes, while Gurobi<sub>k</sub> (resp. Gurobi<sub>s</sub>) in only 4 (resp. 2), and OR-Tools in none.

7 Conclusion

We have studied a variant of the p-dispersion problem where distance constraints exist between the facilities to be dispersed. We proposed ILP formulations and a CP model for this problem. We also devised a heuristic CP-based method built around a bounding technique that prunes the search tree by reasoning about the best possible value of the optimization function at each node. Experimental results demonstrated that although the ILP formulations are more efficient than the CP model, they fail to efficiently handle problems with more than 10 facilities, whereas on such problems the heuristic CP method manages to find solutions of better quality than the ILP and CP models in orders of magnitude shorter run times.
References

1 T. Argo and E. Sandstrom. Separation distances in NFPA codes and standards (tech. rep.). Fire Protection Research Foundation. 2014.

2 O. Berman and R. Huang. The minimum weighted covering location problem with distance constraints. *Computers and Operations Research*, 35(12):356-372, 2008.

3 F. Boussemart, F. Heremy, C. Lecoutre, and L. Sais. Boosting systematic search by weighting constraints. In *Proceedings of ECAI'04*, pages 482-486, 2004.

4 J. Brimberg and H. Juel. A minisum model with forbidden regions for locating a semi-desirable facility in the plane. *Location Science*, 6(1):109-120, 1998.

5 H. Cambazard, D. Mehta, B. O'Sullivan, and L. Quesada. A computational geometry-based local search algorithm for planar location problems. In *Proceedings of the 9th International Conference on the Integration of AI and OR Techniques in Constraint Programming for Combinatorial Optimization Problems (CPAIOR 2012)*, pages 97-112, 2012.

6 S. Chaudhry, T. McCormick, and I. D. Moon. Locating independent facilities with maximum weight: Greedy heuristics. *International Journal of Management Science*, 14(5):383-389, 1986.

7 R. L. Church and M. E. Meadows. Results of a new approach to solving the p-median problem with maximum distance constraints. *Geographical Analysis*, 9(4):364-378, 1977.

8 V. Chvatal. A greedy heuristic for the set-covering problem. *Mathematics of Operations Research*, 4(3):233-235, 1979.

9 W. Comley. The location of ambivalent facilities: Use of a quadratic zero-one programming algorithm. *Applied Mathematical Modeling*, 19(1):26-29, 1995.

10 Z. Dai, K. Xu, and M. Ornik. Repulsion-based p-dispersion with distance constraints in non-convex polygons. *Annals of Operations Research*, 307:75-91, 2021.

11 Choco development team. An Open-Source java library for constraint programming. https://choco-solver.org/.

12 OR-Tools development team. OR-Tools, CP-SAT solver. https://developers.google.com/optimization/cp/cp_solver.

13 T. Drezner, Z. Drezner, and A. Schöbel. The weber obnoxious facility location model: A big arc small arc approach. *Computers and Operations Research*, 98:240-250, 2018.

14 Z. Drezner, P. Kalczynski, and S. Salhi. The planar multiple obnoxious facilities location problem: A Voronoi based heuristic. *Omega*, 87:105-116, 2019.

15 E. Erkut. The discrete p-dispersion problem. *European Journal of Operational Research*, 46(1):48-60, 1990.

16 E. Erkut and S. Neuman. Analytical models for locating undesirable facilities. *European Journal of Operational Research*, 40(3):275-291, 1989.

17 M. M. Fazel-Zarandi and J. C. Beck. Solving a location-allocation problem with logic-based benders' decomposition. In *Proceedings of the 15th International Conference on Principles and Practice of Constraint Programming (CP 2009)*, pages 344-351, 2009.

18 J. B. Ghosh. Computational aspects of the maximum diversity problem. *Operations Research Letters*, 19(4):175-181, 1996.

19 C. Gomes and M. Sellmann. Streamlined constraint reasoning. In *Proceedings of the International Conference on Principles and Practice of Constraint Programming (CP 2004)*, pages 274-289, 2004.

20 T. Guns. Increasing modeling language convenience with a universal n-dimensional array, CPpy as python-embedded example. In *Proceedings of the 18th workshop on Constraint Modelling and Reformulation*, 2019.

21 LLC Gurobi Optimization. Gurobi optimizer reference manual. 2023. URL: https://www.gurobi.com.

22 N. Isoart and J.-C. Régin. A k-Opt Based Constraint for the TSP. In *Proceedings of the 27th International Conference on Principles and Practice of Constraint Programming (CP 2021)*, 2021.
30:18 The p-Dispersion Problem with Distance Constraints

23 B. M. Khumawala. An efficient algorithm for the p-median problem with maximum distance constraints. *Geographical Analysis*, 5(4):309–321, 1973.

24 R. K. Kincaid. Good solutions to discrete noxious location problems via meta-heuristics. *Annals of Operations Research*, 40(1):265–281, 1992.

25 M. J. Kuby. Programming models for facility dispersion: The p-dispersion and maxisum dispersion problems. *Mathematical and Computer Modelling*, 10(10):792, 1988.

26 Mikael Zayenz Lagerkvist and Magnus Rattfeldt. Half-checking propagators. In *Proceedings of the 19th workshop on Constraint Modelling and Reformulation*, 2020.

27 A. Maier and H.W. Hamacher. Complexity results on planar multifacility location problems with forbidden regions. *Mathematical Methods of Operations Research*, 89:433–484, 2019.

28 R. Marti, A. Martinez-Gavara, S. Perez-Pelo, and J. Sanchez-Oro. A review on discrete diversity and dispersion maximization from an OR perspective. *European Journal of Operational Research*, 299(3):795–813, 2022.

29 I. D. Moon and S. Chaudhry. An analysis of network location problems with distance constraints. *Management Science*, 30(3):290–307, 1984.

30 I. D. Moon and L. Papayanopoulos. Minimax location of two facilities with minimum separation: Interactive graphical solutions. *Journal of the Operations Research Society*, 42:685–694, 1991.

31 C.S. Orloff. A theoretical model of net accessibility in public facility location. *Geographical Analysis*, 9:244–256, 1977.

32 M. G. Resende, R. Marti, M. Gallego, and A. Duarte. Grasp and path relinking for the max-min diversity problem. *Computers and Operations Research*, 37(3):498–508, 2010.

33 B. Saboonchi, P. Hansen, and S. Perron. MaxMinMin p-dispersion problem: A variable neighborhood search approach. *Computer & Operations Research*, 52:251–259, 2014.

34 D. Sayah and S. Irnich. A new compact formulation for the discrete p-dispersion problem. *European Journal of Operational Research*, 256(1):62–67, 2017.

35 F. Sayyady and Y Fathi. An integer programming approach for solving the p-dispersion problem. *European Journal of Operational Research*, 253(1):216–225, 2016.

36 D. R. Shier. A min-max theorem for p-center problems on a tree. *Transportation Science*, 11(3):243–252, 1977.

37 S.Y.D. Sorkhabi, D.A. Romero, J.C. Beck, and C. Amon. Constrained multi-objective wind farm layout optimization: Novel constraint handling approach based on constraint programming. *Renewable Energy*, 126(C):341–353, 2018.

38 B.C. Tansel, R.L. Francis, T.J. Lowe, and M.L. Chen. Duality and distance constraints for the nonlinear p-center problem and covering problem on a tree network. *Operations Research*, 30(4):725–744, 1982.

39 S.B. Welch and S. Salhi. The obnoxious p facility network location problem with facility interaction. *European Journal of Operations Research*, 102:302–319, 1997.

40 49 C.F.R. §175.701. Separation distance requirements for packages containing class 7 (radio-active) materials in passenger-carrying aircraft. Title 49 Code of Federal Regulations, Part 175. 2021.

41 29 C.F.R. §1910.157. Portable fire extinguishers. Title 29 Code of Federal Regulations, Part 157. 2021.
# Partially Preemptive Multi Skill/Mode Resource-Constrained Project Scheduling with Generalized Precedence Relations and Calendars

Guillaume Povéda ✉️ ID

Airbus (AI Research), Toulouse, France

Nahum Alvarez ✉️ ID

Airbus (AI Research), Toulouse, France

Christian Artigues ✉️ ID

LAAS-CNRS, Université de Toulouse, CNRS, Toulouse, France

## Abstract

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.

**2012 ACM Subject Classification** Computing methodologies → Planning and scheduling; Applied computing → Industry and manufacturing; Theory of computation → Optimization with randomized search heuristics; Theory of computation → Constraint and logic programming

**Keywords and phrases** Large-scale scheduling problem, partial preemption, multi-skill, multi-mode, resource calendars, constraint programming, large neighborhood search

**Digital Object Identifier** 10.4230/LIPIcs.CP.2023.31

**Funding** This work received the support of the Support from the French ANR-3IA Artificial and Natural Intelligence Toulouse Institute (ANITI).

Guillaume Povéda: ANITI

Nahum Alvarez: ANITI

Christian Artigues: ANITI

**Acknowledgements** The authors are grateful to the anonymous reviewers for their constructive comments. In particular, we thank the anonymous reviewer that pointed out the issue linked to preemption in the *CP-Base* model and made inspiring suggestions that leaded to the *CP-SmartPreemption* model.

# 1 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

<figure><img src="image_294.png" /></figure>

© Guillaume Povéda, Nahum Alvarez, and Christian Artigues; licensed under Creative Commons License CC-BY 4.0

29th International Conference on Principles and Practice of Constraint Programming (CP 2023).

Editor: Roland H. C. Yap; Article No. 31; pp. 31:1–31:21

<figure><img src="image_295.png" /></figure>

LIPICS

Leibniz International Proceedings in Informatics

Schloss Dagstuhl – Leibniz-Zentrum für Informatik, Dagstuhl Publishing, Germany
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.

2 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.

### 3 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
31:4 Partially Preemptive Multi Skill/Mode Resource-Constrained Project Scheduling

calendars are also added. Formally a PP-MS-MM-RCPSP/max-cal is defined by:

1. A a set of activities to schedule, $\mathcal{A} = \{1, ..., n\}$, two special activities are created by convention: the source task 1 (predecessor of all other tasks) and sink task n (successor of all other tasks);

2. Each activity $i \in \mathcal{A}$ can be performed in $m_i$ different modes, $m_i \in M_i = \{1, ..., |M_i|\}$; (as defined in [28])

3. $\mathcal{R}^\rho$ a set of renewable resource types, $\mathcal{R}^\rho = \{R_1^\rho, ..., R_{m^\rho}^\rho\}$. $\forall k \in \mathcal{R}^\rho, t \in \mathbb{N}, B_{kt}$ is the discrete amount of available resource $k$ at discrete time $t$ (available between time $t$ and $t+1$);

4. $\mathcal{R}^\nu$ a set of non-renewable resource types, $\mathcal{R}^\nu = \{R_1^\nu, ..., R_{m^\nu}^\nu\}$; $\forall k \in \mathcal{R}^\nu, B_k$ is the total capacity of the non-renewable resource;

5. $\mathcal{O}$ the set of disjunctive resources representing individual skilled workers (from now on, we will refer to them as Operators);

6. $\mathcal{L} = \{1, ..., L\}$ the set of skills;

7. $\forall o \in \mathcal{O}, t \in \mathbb{N}, A_{ot} \in \{0, 1\}$ indicates if the Operator $o$ is present or not at time $t$ (hence we treat the temporal availability of Operators as fixed from the problem definition);

8. $\forall o \in \mathcal{O}, l \in \mathcal{L}, y_{ol} \in \{0, 1\}$ indicates if Operator $o$ masters skill $l$;

9. $\forall i \in \mathcal{A}, r_i$ is the release time of task $i$;

10. $\forall i \in \mathcal{A}, d_i$ is the deadline time of task $i$;

11. $\forall i \in \mathcal{A}, p_{i,m_i}$ is the processing time under mode $m_i \in M_i$;

12. $\forall i \in \mathcal{A}, k \in \mathcal{R}^\rho \cup \mathcal{R}^\nu, b_{i,m_i,k}$ is a natural number representing the resource demand of activity $i$ for resource $k$, under mode $m_i$;

13. $\forall i \in \mathcal{A}, l \in \mathcal{L}, s_{i,m_i,l}$ is a natural number representing the skill requirement of activity $i$ under mode $m_i$;

14. $P$ is the set of precedence constraints $\mathcal{A} \times \mathcal{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 $\Delta 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 $\Delta e_{(i,j)}$, the minimum time lag between end of $i$ and start of $j$;

19. The set $\mathcal{A}$ is split in three different subsets :

- $\mathcal{A} = \mathcal{A}^P \cup \mathcal{A}^{NP} \cup \mathcal{A}^{PP}$;

- $\mathcal{A}^P$ is the set of fully preemptive activities (activities can be preempted and all the resources are released);

- $\mathcal{A}^{NP}$ is the set of non-preemptive activities;

- $\mathcal{A}^{PP}$ is the set of partially preemptive activities (where at least one resource is not releasable);

20. $\forall i \in \mathcal{A}, m_i \in M_i, k \in \mathcal{R}^\rho, \rho_{i,m_i,k} \in \{0, 1\}$ indicates if resource $k$ is releasable for activity $i$ under mode $m_i$.

**3.1 Known variants in the literature**

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

- the classical RCPSP ($\mathcal{L} = \emptyset, \mathcal{R}^\nu = \emptyset, \forall k \in \mathcal{R}^\rho, \forall t, t' \in \mathbb{N}, B_{kt} = B_{kt'}, \mathcal{O} = \emptyset, \mathcal{A}^P = \emptyset, \mathcal{A}^{PP} = \emptyset, \forall i \in \mathcal{A}, M_i = \{1\}$);
- the classical multi-mode RCPSP ($\mathcal{L} = \emptyset, \forall k \in \mathbb{R}^\rho, \forall t, t' \in \mathbb{N}, B_{kt} = B_{kt'}, \mathcal{O} = \emptyset, \mathcal{A}^P = \emptyset, \mathcal{A}^{PP} = \emptyset, \forall i \in \mathcal{A}, M_i = \{1, ..., |M_i|\}$);

- the classical multi-skill RCPSP ($\mathcal{L} = \{1, ..., L\}, \mathcal{O} \neq \emptyset, \mathcal{A}^P = \emptyset, \mathcal{A}^{PP} = \emptyset, \forall i \in \mathcal{A}, M_i = \{1, ..., |M_i|\}$).

## 3.2 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<fnref n="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.

### 3.2.1 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 $\mathcal{J} = \{1, ..., max_{preemption}\}$ the set of preemption breaking time indexes and $\mathcal{J}^- = \{2, ..., max_{preemption}\}$

1. Starting time decision : starts is a $|\mathcal{A}| \times max_{preemption}$ matrix, which contains the starting time of subparts of all the tasks in increasing order.

2. Duration decision : durations is a $|\mathcal{A}| \times max_{preemption}$ matrix, which contains the duration of subparts of all the tasks.

3. Mode decision : modes is a $|\mathcal{A}|$ vector specifying in which mode a task is executed. $\forall i \in \mathcal{A}, modes[i] \in M_i$.

4. Resource allocation : allocation is a $|\mathcal{O}| \times |\mathcal{A}| \times max_{preemption}$ 3D binary matrix which indicates which worker o is allocated to each subpart of an activity.

### 3.2.2 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 $|\mathcal{A}| \times max_{preemption}$ matrix allocation of resource o to all activity subparts.

1. Resource consumption constraint :
   $\forall k \in \mathbb{R}^\rho, \text{cumulative}(\text{starts}, \text{durations}, b_{\text{modes},k}, B_k)$
   where $b_{\text{modes},k}$ is an array of dimension $|\mathcal{A}| \times max_{preemption}$ defined by:

$\forall i \in \mathcal{A}, j \in \mathcal{J}, b_{\text{modes},k}[i, j] = b_{i,\text{modes}[i],k}$, being $b_{\text{modes},k}[i, j]$ 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].

<footnote n="1">¹ https://github.com/airbus/discrete-optimization</footnote>
31:6 Partially Preemptive Multi Skill/Mode Resource-Constrained Project Scheduling

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.

2. Non-renewable resources:
   $\forall k \in \mathcal{R}^v, \sum_{i \in A} b_{i,modes[i],k} \le B_k$. 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.

3. Skills requirement:
   $\forall i \in A, j \in J, l \in L, \\ \text{durations}[i, j] > 0 \rightarrow \sum_{o \in O} \text{allocation}[o, i, j] \cdot y_{ol} \ge s_{i,modes[i],l}$
   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].

4. Operator availability:
   $\forall o \in O, \text{cumulative}(\text{starts}, \text{durations}, \text{allocation}[o, : , :], A_o)$
   Variable availability of an Operator is dealt with in the same way as Constraint 1 above.

5. Precedence relation: $\forall (i, j) \in P, \\ \text{starts}[j, 1] \ge \text{starts}[i, \max_{\text{preemption}}] + \text{durations}[i, \max_{\text{preemption}}]$

6. $P_{\text{sync-start}}$ relations:
   $\forall (i, j) \in P, \text{starts}[i, 1] = \text{starts}[j, 1]$

7. $P_{\text{startlag}}$ relations:
   $\forall (i, j) \in P_{\text{startlag}}, \text{starts}[j, 1] \ge \text{starts}[i, 1] + \Delta s_{(i,j)}$

8. $P_{\text{sync-end}}$ relations: $\forall (i, j) \in P_{\text{sync-end}}, \\ \text{starts}[j, 1] = \text{starts}[i, \max_{\text{preemption}}] + \text{durations}[i, \max_{\text{preemption}}]$

9. $P_{\text{endlag}}$ relations:
   $\forall (i, j) \in P_{\text{endlag}}, \\ \text{starts}[j, 1] \ge \text{starts}[i, \max_{\text{preemption}}] + \text{durations}[i, \max_{\text{preemption}}] + \Delta e_{(i,j)}$

10. Tasks duration:
   a. $\forall i \in A^{NP}, \text{durations}[i, 1] = p_{i,modes[i]} \wedge (\forall j \in [2, \max_{\text{preemption}}], \text{durations}[i, j] = 0)$

b. $\forall i \in A \setminus A^{NP}, \sum_{j \in J} \text{durations}[i, j] = p_{i,modes[i]}$

11. Release time:
   $\forall i \in A, \text{starts}[i, 1] \ge r_i$

12. Deadline time:
   $\forall i \in A, \text{starts}[i, \max_{\text{preemption}}] + \text{durations}[i, \max_{\text{preemption}}] \le d_i$

13. Conventions constraints for starts and durations:
   The following constraints are modeling choices aiming to help the solver to explore the search space.

a. $\forall i \in A, j \in J^-, \text{starts}[i, j] \ge \text{starts}[i, j - 1] + \text{durations}[i, j]$: precedence constraints between each subparts of the task.

b. $\forall i \in A, j \in [1, \max_{\text{preemption}} - 1], \text{durations}[i, j] = 0 \rightarrow \text{durations}[i, j + 1] = 0$: As soon as there is a 0 duration subtask, all the following ones are 0 too.

c. $\forall i \in A, j \in [1, \max_{\text{preemption}}], \text{durations}[i, j] = 0 \rightarrow \text{starts}[i, j] = \text{starts}[i, j - 1] + \text{durations}[i, j - 1]$; Whenever there is a 0 duration, the remaining starts values are uniquely determined by previous values, pruning redundant solutions.
14. Partially preemptive and non releasable resources :

We introduce variable blockedduration[i, j, k], ∀i ∈ A, j ∈ J, k ∈ R p as the usage duration of resource k for task i and activity subpart j. It is regulated by the following constraint: ∀j ∈ J,

$$
\rho_{i, modes[i], k} = 1 \rightarrow \text{blockedduration}[i, j, k] = \text{durations}[i, 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:

$$
\begin{align*}
\rho_{i, modes[i], k} = 0 &\rightarrow (\text{blockedduration}[i, 1, k] = \text{starts}[i, \max_{preemption}] - \text{starts}[i, 1]) \land \\
& (\text{blockedduration}[i, j, k] = 0, \forall j \in \mathcal{J}^{-})
\end{align*}
$$

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*:

$$
\forall k \in \mathcal{R}^p, \text{cumulative}(\text{starts}, \text{blockedduration}, b_{\text{modes},k}, B_k)
$$

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: $\forall i \in 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 :

15. $\forall i \in A, \text{s.t } p_{i, modes[i]} \ge 1, \forall j \in \mathcal{J}, interval[i, j].duration \ge 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.

16. $\forall i \in A, j \in \mathcal{J}^{-}, interval[i, j].start \ge interval[i, j - 1].end$: precedence constraints between sub-intervals of the same task.

17. $\forall i \in A, j \in [1, max_{preemption} - 1], interval[i, j].present = False \rightarrow interval[i, j + 1].present = False$: When one sub-interval is not present, the remaining ones are not present either. This is equivalent to 13b constraint.

18. $\forall i \in A, \sum_{j \in \mathcal{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. $\forall i \in A, span(spantask[i], interval[i, :]):$ We use the native constraint span of CPOptim- izer 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 CPOptimization, Model *CP-SmartPreemption* while Model (1-14) is called *CP-Base*.
### 3.2.3 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).

## 3.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 <fnref n="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 : $mode_{A_0} = 1, mode_{A_1} = 2, mode_{A_2} = 1, mode_{A_3} = 2, mode_{A_4} = 2, mode_{A_5} = 1$.

<table><caption>Table 1 An instance of PP-MS-MM-RCPSP/max-cal.</caption><thead><tr><th>Activity</th><th>Mode</th><th>Duration</th><th>Skills</th><th>Resource</th><th>Deadline</th><th>Release</th><th>Preemption</th><th>Successors</th></tr></thead><tbody><tr><td>A<sub>0</sub></td><td>1</td><td>0</td><td>-</td><td>-</td><td>-</td><td>-</td><td>-</td><td>A<sub>1</sub>, A<sub>2</sub>, A<sub>3</sub>, A<sub>4</sub>, A<sub>5</sub></td></tr><tr><td rowspan="2">A<sub>1</sub></td><td>1</td><td>5</td><td>l<sub>1</sub></td><td>R<sub>1</sub></td><td>-</td><td>-</td><td>A<sup>PP</sup></td><td>A<sub>5</sub></td></tr><tr><td>2</td><td>3</td><td>l<sub>1</sub></td><td>R<sub>1</sub></td><td>5</td><td>-</td><td>A<sup>PP</sup></td><td>A<sub>5</sub></td></tr><tr><td rowspan="2">A<sub>2</sub></td><td>1</td><td>1</td><td>l<sub>3</sub>, l<sub>4</sub></td><td>(R<sub>1</sub>, 1)</td><td>-</td><td>2</td><td>A<sup>NP</sup></td><td>A<sub>3</sub>, A<sub>5</sub></td></tr><tr><td>2</td><td>2</td><td>l<sub>3</sub></td><td>(R<sub>1</sub>, 1)</td><td>-</td><td>2</td><td>A<sup>NP</sup></td><td>A<sub>5</sub></td></tr><tr><td rowspan="2">A<sub>3</sub></td><td>1</td><td>3</td><td>l<sub>2</sub></td><td>(R<sub>1</sub>, 1)</td><td>-</td><td>-</td><td>A<sup>PP</sup></td><td>A<sub>5</sub></td></tr><tr><td>2</td><td>2</td><td>l<sub>3</sub>, l<sub>2</sub></td><td>(R<sub>1</sub>, 1)</td><td>-</td><td>-</td><td>A<sup>PP</sup></td><td>A<sub>5</sub></td></tr><tr><td rowspan="2">A<sub>4</sub></td><td>1</td><td>2</td><td>l<sub>3</sub></td><td>-</td><td>-</td><td>5</td><td>A<sup>P</sup></td><td>A<sub>5</sub></td></tr><tr><td>2</td><td>3</td><td>-</td><td>(R<sub>1</sub>, 1)</td><td>-</td><td>5</td><td>A<sup>P</sup></td><td>A<sub>5</sub></td></tr><tr><td>A<sub>5</sub></td><td>1</td><td>0</td><td>-</td><td>-</td><td>-</td><td>-</td><td>-</td><td>-</td></tr></tbody></table>

<table><thead><tr><th>Operator</th><th>Skills</th><th>Calendar</th></tr></thead><tbody><tr><td>O<sub>1</sub></td><td>l<sub>1</sub>, l<sub>3</sub></td><td>[1, 0, 0, 1, 1, 1, 0, 0, 0, 0, 1, 0, 0, 0]</td></tr><tr><td>O<sub>2</sub></td><td>l<sub>1</sub>, l<sub>2</sub>, l<sub>4</sub></td><td>[1, 1, 0, 0, 1, 1, 0, 0, 0, 0, 1, 0, 0, 1, 1]</td></tr></tbody></table>

<table><thead><tr><th>Resource</th><th>Capacity</th></tr></thead><tbody><tr><td>R<sub>1</sub></td><td>2</td></tr></tbody></table>

<table><tr><td>P<sub>sync-start</sub></td><td>[(A<sub>3</sub>, A<sub>4</sub>)]</td></tr></table>

<footnote n="2">² https://github.com/g-poveda/do_experiments</footnote>
<figure><img src="image_296.png" /><figcaption>Figure 1 Operator and resource oriented Gantt chart for the example instance.</figcaption></figure>

# 4 Algorithms

## 4.1 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.

Begin

1: $Y^* = \infty, X^* = None$

2: $(X^0, Y^0) = (X^*, Y^*) = initial_{solution}(P)$

3: $iter = 0$

4: repeat

5: $RMP = buildsubproblem(MP, X^{iter})$

6: $X^{iter+1}, Y^{iter+1} = solve(RMP)$

7: if $Y^{iter+1} \le Y^*$ then

8: $X^* \leftarrow X^{iter+1}$

9: $Y^* \leftarrow Y^{iter+1}$

10: $iter \leftarrow iter + 1$

11: until stop criterion is met

12: return $X^*, Y^*$

## 4.2 Application to the PP-MS-MM-RCPSP/max-cal

### 4.2.1 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 \in A$, another priority ordering $order\_res_i \in S|_o$ of the workers $o \in O$ to be assigned to the activity and the mode id $modes_i$ chosen to run the activity, as defined for the CP model.
31:10 Partially Preemptive Multi Skill/Mode Resource-Constrained Project Scheduling

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<sub>sync-start/end</sub> 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.

4.2.2 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>sub</sub> and A<sub>sub</sub><sup>−</sup>. Different methods have been implemented and tested to select A<sub>sub</sub> and A<sub>sub</sub><sup>−</sup> as described below:

1. Random selection : given $f \in \mathbb{R}$, $0 \le f \le 1$,
   $A_{sub} = \text{sample}(A, \lfloor f \cdot |A| \rfloor)$, and $A_{sub}^{\text{−}} = A \setminus A_{sub}$

2. Random selection and neighbors : consist in the previous selection, then add the predecessors of each task in $\mathcal{A}_{sub}$. Formally, given $f \in \mathbb{R}$, $0 \le f \le 1$:

$\mathcal{A}_{sub,1} = \text{sample}(\mathcal{A}, \lfloor f \cdot |\mathcal{A}| \rfloor)$, and
$\mathcal{A}_{sub} = \mathcal{A}_{sub,1} \cup \{j \in \mathcal{A} \text{ s.t. } \exists (x, y) \in P, x = j \land y \in \mathcal{A}_{sub,1}\}$

3. Cut in equal parts : given $c \in \mathbb{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 $\mathcal{A}_{sub}$ and their predecessors in the precedence graph.

5. Mixing methods : given a pool containing the previous four methods, it will return $\mathcal{A}_{sub}$ and $\mathcal{A}_{sub}^{\text{−}}$ 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 ($\mathcal{A}_{sub}$ and $\mathcal{A}_{sub}^{\text{−}}$) 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 $\mathcal{A}_{sub}^{\text{−}}$ and keeping free for optimisation
the variables of $\mathcal{A}_{sub}$, we consider different $\epsilon$ values to constrain start and duration variables.
We typically assign high values to parameters $\epsilon_{sub,*}$ and small values to $\epsilon_{sub,*}$.

■ **Algorithm 2** Constraint builder procedure.

1: constraintbuilder($\mathcal{A}_{sub}, \mathcal{A}_{\overline{sub}}, X, \text{model}_{cp}, \text{params}_{\text{constraints}}$):

2: $\epsilon_{sub,+}, \epsilon_{sub,-}, \epsilon_{\overline{sub},+}, \epsilon_{\overline{sub},-}, \epsilon_{\overline{dur}_{sub,+}}, \epsilon_{\overline{dur}_{sub,-}}, \epsilon_{\overline{dur}_{sub,+}}, \epsilon_{\overline{dur}_{sub,-}} = \text{params}_{\text{constraints}}$

3: **for all** $i \in \mathcal{A}$ **do**

4: **for all** $j \in [1, \max_{\text{preemption}}]$ **do**

5: **if** $X.durations[i, j] > 0$ **then**

6: **if** $i \in \mathcal{A}_{sub}$ **then**

7: $(\epsilon_{-}, \epsilon_{+}) \leftarrow (\epsilon_{sub,-}, \epsilon_{sub,+})$

8: $(\epsilon_{\overline{dur}_{-}}, \epsilon_{\overline{dur}_{+}}) \leftarrow (\epsilon_{\overline{dur}_{sub,-}}, \epsilon_{\overline{dur}_{sub,+}})$

9: **else if** $i \in \mathcal{A}_{\overline{sub}}$ **then**

10: $(\epsilon_{-}, \epsilon_{+}) \leftarrow (\epsilon_{\overline{sub},-}, \epsilon_{\overline{sub},+})$

11: $(\epsilon_{\overline{dur}_{-}}, \epsilon_{\overline{dur}_{+}}) \leftarrow (\epsilon_{\overline{dur}_{\overline{sub},-}}, \epsilon_{\overline{dur}_{\overline{sub},+}})$

12: $\text{model}_{cp} \leftarrow X.starts[i, j] - \epsilon_{-} \le \text{model}_{cp}.starts[i, j] \le X.starts[i, j] + \epsilon_{+}$

13: $\text{model}_{cp} \leftarrow \max(0, X.durations[i, j] - \epsilon_{\overline{dur}_{-}}) \le \text{model}_{cp}.durations[i, j] \le X.durations[i, j] + \epsilon_{\overline{dur}_{+}}$

### 4.2.3 Subproblem solving

To solve the subproblem we use the open source lazy clause generation solver Chuffed<fnref n="3" /> or alternatively IBM's CP-Optimizer<fnref n="4" /> 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 $\mathcal{A}_{sub}$ we consider the most important in the $\mathcal{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 $\mathcal{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.

## 4.3 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_{\text{sync-start}}, P_{\text{sync-end}}, P_{\text{startlag}}, P_{\text{endlag}}$. By relaxing constraints (6), (7), (8), (9), (12) from Subsection 3.2.2 we include a violation penalty into the objective function. This violation penalty for $P_{\text{sync-end}}$ is illustrated in Figure 2. The objective function that LNS and underlying CP solver minimizes is `makespan` + $M \cdot$ 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)

<footnote n="3">³ https://github.com/chuffed/chuffed</footnote>

<footnote n="4">⁴ https://www.ibm.com/analytics/cplex-cp-optimizer</footnote>
<figure><img src="image_297.png" /><figcaption>Figure 2 $\Delta t$ is added to the objective function whenever $(A1, A2) \in P_{sync-end}$.</figcaption></figure>

# 5 Experimental results

## 5.1 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 $B_k : \forall k \in \mathcal{R}^\rho$ always constant and $A_o : \forall o \in \mathcal{O}$ always 1.

- 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, $\mathcal{A}^{NP} = \mathcal{A}$

- Multiskill: skilled workers are present or not; in the non multiskill use case, workers belong to $\mathcal{R}$ and $\mathcal{O} = \emptyset$

<table><caption>Table 2 Base instance description of the manufacturing use case.</caption><thead><tr><th>Instance</th><th>|A|</th><th>|R|</th><th>|L|</th><th>|O|</th><th>|P|</th><th>|P<sub>ss</sub>|</th><th>|P<sub>se</sub>|</th><th>|P<sub>sl</sub>|</th><th>|P<sub>el</sub>|</th></tr></thead><tbody><tr><td>A</td><td>291</td><td>8</td><td>23</td><td>27</td><td>842</td><td>3</td><td>70</td><td>0</td><td>4</td></tr><tr><td>B</td><td>199</td><td>3</td><td>6</td><td>6</td><td>676</td><td>1</td><td>34</td><td>6</td><td>4</td></tr></tbody></table>

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 consistently 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* (abbreviated 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.

<table><thead><tr><th>ID</th><th>Features</th><th>LNS</th><th>CP-S</th><th>SA</th><th>p-SA</th><th>ID</th><th>Features</th><th>LNS</th><th>CP-S</th><th>SA</th><th>p-SA</th></tr></thead><tbody><tr><td>0</td><td>A</td><td><b>1080</b></td><td><b>1080</b></td><td><b>1080</b></td><td>0</td><td>18</td><td>B G</td><td><b>3053</b></td><td><b>3053</b></td><td>3275</td><td>2171</td></tr><tr><td>1</td><td>A M</td><td><b>1080</b></td><td><b>1080</b></td><td><b>1080</b></td><td>0</td><td>19</td><td>B G M</td><td><b>3109</b></td><td><b>3109</b></td><td>3530</td><td>2590</td></tr><tr><td>2</td><td>A G</td><td><b>6534</b></td><td><b>6534</b></td><td>6817</td><td>128</td><td>20</td><td>B P</td><td>2340</td><td>2340</td><td><b>2339</b></td><td>0</td></tr><tr><td>3</td><td>A G M</td><td>6553</td><td><b>6534</b></td><td>6543</td><td>104</td><td>21</td><td>B M P</td><td><b>2340</b></td><td>2390</td><td><b>2340</b></td><td>0</td></tr><tr><td>4</td><td>A P</td><td><b>1080</b></td><td><b>1080</b></td><td><b>1080</b></td><td>0</td><td>22</td><td>B G P</td><td>3053</td><td><b>3044</b></td><td>3148</td><td>750</td></tr><tr><td>5</td><td>A M P</td><td><b>1080</b></td><td>2495</td><td><b>1080</b></td><td>0</td><td>23</td><td>B G M P</td><td><b>3117</b></td><td>3180</td><td>3456</td><td>714</td></tr><tr><td>6</td><td>A G P</td><td><b>6534</b></td><td><b>6534</b></td><td><b>6534</b></td><td>0</td><td>24</td><td>B C</td><td>3202</td><td><b>3196</b></td><td>3250</td><td>0</td></tr><tr><td>7</td><td>A G M P</td><td>6553</td><td>6612</td><td><b>6537</b></td><td>0</td><td>25</td><td>B C M</td><td>3289</td><td><b>3208</b></td><td>3263</td><td>0</td></tr><tr><td>12</td><td>A C P</td><td><b>7154</b></td><td>-</td><td><b>7154</b></td><td>0</td><td>26</td><td>B C G</td><td><b>4949</b></td><td><b>4949</b></td><td>5096</td><td>4727</td></tr><tr><td>13</td><td>A C M</td><td><b>7154</b></td><td>-</td><td><b>7154</b></td><td>0</td><td>27</td><td>B C G M</td><td>5096</td><td><b>4949</b></td><td>5096</td><td>4359</td></tr><tr><td>14</td><td>A C G P</td><td>7866</td><td>-</td><td><b>7777</b></td><td>0</td><td>28</td><td>B C P</td><td>3205</td><td>3195</td><td><b>3192</b></td><td>0</td></tr><tr><td>15</td><td>A C G M P</td><td><b>7777</b></td><td>-</td><td>7805</td><td>0</td><td>29</td><td>B C M P</td><td>3193</td><td>3208</td><td><b>3192</b></td><td>0</td></tr><tr><td>16</td><td>B</td><td>2340</td><td><b>2339</b></td><td>2373</td><td>0</td><td>30</td><td>B C G P</td><td>4580</td><td><b>4383</b></td><td>4436</td><td>1167</td></tr><tr><td>17</td><td>B M</td><td>2350</td><td><b>2346</b></td><td>2375</td><td>0</td><td>31</td><td>B C G M P</td><td><b>5010</b></td><td>-</td><td>4840</td><td>1132</td></tr></tbody></table>

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.

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.

## 5.2 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^{NP}$ and $A^{PP}$.
<figure><img src="image_298.png" /><figcaption>(a) Instance n°15.</figcaption></figure>

<caption>(b) Instance n°31.</caption>

■ **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 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.

<table><thead><tr><th>Algorithm</th><th>Set A</th><th>Set B</th><th>Set C</th><th>Set D</th></tr></thead><tbody><tr><td>GRASP + LNS [26]</td><td><strong>1.56%</strong></td><td><strong>1.79%</strong></td><td>6.68%</td><td><strong>2.06%</strong></td></tr><tr><td>SA (300 seconds)</td><td><u>1.59%</u></td><td><u>1.82%</u></td><td><strong>5.65%</strong></td><td>2.33%</td></tr><tr><td>SA (240 s)</td><td>1.82%</td><td>1.83%</td><td>6.57%</td><td>2.17%</td></tr><tr><td>SA (240 s)+LNS (60 s)</td><td>1.82%</td><td>1.83%</td><td><u>6.49%</u></td><td><u>2.12%</u></td></tr></tbody></table>

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*<sub>all</sub> mastered by all operators. We then set skill requirement of skill *s*<sub>all</sub> 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 *LNS* 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.

## 5.3 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 test 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] <fnref n="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 LNS 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 LNS 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 LNS 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 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.

## 5.4 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<fnref n="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.

<footnote n="5">⁵ https://github.com/youngkd/MSPSP-InstLib/blob/master/models/mspsp.mzn</footnote>

<footnote n="6">⁶ https://github.com/g-poveda/do_experiments</footnote>
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.

# 6 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.

# References

1. Behrouz Afshar-Nadjafi. Multi-skilling in scheduling problems: A review on models, methods and applications. *Computers & Industrial Engineering*, 151:107004, 2021.

2. Bernardo F Almeida, Isabel Correia, and Francisco Saldanha-da Gama. Priority-based heuristics for the multi-skill resource constrained project scheduling problem. *Expert Systems with Applications*, 57:91-103, 2016.
3 Bernardo F Almeida, Isabel Correia, and Francisco Saldanha-da Gama. Modeling frameworks for the multi-skill resource-constrained project scheduling problem: a theoretical and empirical comparison. *International Transactions in Operational Research*, 26(3):946–967, 2019.

4 Lucio Bianco and Massimiliano Caramia. An exact algorithm to minimize the makespan in project scheduling with scarce resources and generalized precedence relations. *European Journal of Operational Research*, 219(1):73–85, 2012.

5 P. Brucker. *Complex Scheduling*. Springer, 1999.

6 José Coelho and Mario Vanhoucke. Multi-mode resource-constrained project scheduling using rcpsp and sat solvers. *European Journal of Operational Research*, 213(1):73–82, 2011.

7 Erik L Demeulemeester and Willy S Herroelen. An efficient optimal solution procedure for the preemptive resource-constrained project scheduling problem. *European Journal of Operational Research*, 90(2):334–348, 1996.

8 Emir Demirović, Geoffrey Chu, and Peter J Stuckey. Solution-based phase saving for cp: A value-selection heuristic to simulate local search behavior in complete solvers. In *Principles and Practice of Constraint Programming: 24th International Conference, CP 2018, Lille, France, August 27-31, 2018, Proceedings 24*, pages 99–108. Springer, 2018.

9 Andreas Drexl and Juergen Gruenewald. Nonpreemptive multi-mode resource-constrained project scheduling. *IIE Transactions*, 25(5):74–81, 1993.

10 Davide Giglio, Massimo Paolucci, Abdolreza Roshani, and Flavio Tonelli. Multi-manned assembly line balancing problem with skilled workers: a new mathematical formulation. *IFAC-PapersOnLine*, 50(1):1211–1216, 2017.

11 Daniel Godard, Philippe Laborie, and Wim Nuijten. Randomized large neighborhood search for cumulative scheduling. In *ICAPS*, volume 5, pages 81–89, 2005.

12 Shima Javanmard, Behrouz Afshar-Nadjafi, and Seyed Taghi Akhavan Niaki. Preemptive multi-skilled resource investment project scheduling problem: Mathematical modelling and solution approaches. *Computers & Chemical Engineering*, 96:55–68, 2017.

13 Ahmed Karam, El-Awady Attia, and Philippe Duquenne. A milp model for an integrated project scheduling and multi-skilled workforce allocation with flexible working hours. *IFAC-PapersOnLine*, 50(1):13964–13969, 2017.

14 S. Kirkpatrick, C. D. Gelatt, and M. P. Vecchi. Optimization by simulated annealing. *Science*, 220(4598):671–680, 1983. doi:10.1126/science.220.4598.671.

15 RAINER KOLISCH and ANDREAS DREXL. Local search for nonpreemptive multi-mode resource-constrained project scheduling. *IIE Transactions*, 29(11):987–999, 1997. doi:10.1080/07408179708966417.

16 Rainer Kolisch and Sönke Hartmann. Heuristic algorithms for the resource-constrained project scheduling problem: Classification and computational analysis. In *Project scheduling*, pages 147–178. Springer, 1999.

17 Stefan Kreter, Andreas Schutt, and Peter J Stuckey. Using constraint programming for solving rcpsp/max-cal. *Constraints*, 22(3):432–462, 2017.

18 Philippe Laborie. An update on the comparison of mip, cp and hybrid approaches for mixed resource allocation and scheduling. In Willem-Jan van Hoeve, editor, *Integration of Constraint Programming, Artificial Intelligence, and Operations Research*, pages 403–411, Cham, 2018. Springer International Publishing.

19 Haitao Li and Keith Womer. Scheduling projects with multi-skilled personnel by a hybrid milp/cp benders decomposition algorithm. *Journal of Scheduling*, 12(3):281–298, 2009.

20 Paweł B. Myszkowski, Maciej Laszczyk, Ivan Nikulin, and Marek Skowronski. imopse: a library for bicriteria optimization in multi-skill resource-constrained project scheduling problem. *Soft Computing*, 23:3397–3410, 2019.

21 Mireille Palpant, Christian Artigues, and Philippe Michelon. Lssper: Solving the resource-constrained project scheduling problem with large neighbourhood search. *Annals of Operations Research*, 131(1):237–257, 2004.
31:18 Partially Preemptive Multi Skill/Mode Resource-Constrained Project Scheduling

22 Claude Le Pape and Philippe Baptiste. Heuristic control of a constraint-based algorithm for the preemptive job-shop scheduling problem. *Journal of Heuristics*, 5(3):305-325, 1999.

23 Laurent Perron, Paul Shaw, and Vincent Furnon. Propagation guided large neighborhood search. In *International Conference on Principles and Practice of Constraint Programming*, pages 468-481. Springer, 2004.

24 Oliver Polo-Mejía, Marie-Christine Anselmet, Christian Artigues, and Pierre Lopez. Mixed-integer and constraint programming formulations for a multi-skill project scheduling problem with partial preemption. In *12th International Conference on Modelling, Optimization and Simulation (MOSIM 2018)*, pages 367-374, 2018.

25 Oliver Polo-Mejía, Christian Artigues, Pierre Lopez, and Virginie Basini. Mixed-integer/linear and constraint programming approaches for activity scheduling in a nuclear research facility. *International Journal of Production Research*, 58(23):7149-7166, 2020.

26 Oliver Polo-Mejía, Christian Artigues, Pierre Lopez, Lars Mönch, and Virginie Basini. Heuristic and metaheuristic methods for the multi-skill project scheduling problem with partial preemption. *International Transactions in Operational Research*, 2021.

27 Sacramento Quintanilla, Ángeles Pérez, Pilar Lino, and Vicente Valls. Time and work generalised precedence relationships in project scheduling with pre-emption: An application to the management of service centres. *European Journal of Operational Research*, 219(1):59-72, 2012.

28 Jerzy Rosłon. The multi-mode, resource-constrained project scheduling problem in construction: state of art review and research challenges problem. *Technical Transactions*, 2017. doi: 10.4467/2353737XCT.17.070.6427.

29 Mónica A Santos and Anabela P Tereso. On the multi-mode, multi-skill resource constrained project scheduling problem–a software application. In *Soft computing in industrial applications*, pages 239–248. Springer, 2011.

30 Andreas Schutt, Thibaut Feydy, Peter J Stuckey, and Mark G Wallace. Explaining the cumulative propagator. *Constraints*, 16(3):250-282, 2011.

31 Andreas Schutt, Thibaut Feydy, Peter J. Stuckey, and Mark G. Wallace. Solving rcpsp/max by lazy clause generation. *J. Sched.*, 16(3):273–289, 2013. doi:10.1007/s10951-012-0285-x.

32 Andreas Schutt, Thibaut Feydy, Peter J Stuckey, and Mark G Wallace. Solving rcpsp/max by lazy clause generation. *Journal of scheduling*, 16(3):273–289, 2013.

33 Jakob Snauwaert and Mario Vanhoucke. A new algorithm for resource-constrained project scheduling with breadth and depth of skills. *European Journal of Operational Research*, 292(1):43–59, 2021. doi:10.1016/j.ejor.2020.10.032.

34 Jakob Snauwaert and Mario Vanhoucke. A classification and new benchmark instances for the multi-skilled resource-constrained project scheduling problem. *European Journal of Operational Research*, 2022. doi:10.1016/j.ejor.2022.05.049.

35 Tianshu Song, Ke Xu, Jiangneng Li, Yiming Li, and Yongxin Tong. Multi-skill aware task assignment in real-time spatial crowdsourcing. *GeoInformatica*, 24(1):153-173, 2020.

36 Ria Szeredi and Andreas Schutt. Modelling and solving multi-mode resource-constrained project scheduling. In *International Conference on Principles and Practice of Constraint Programming*, pages 483-492. Springer, 2016.

37 Tony Wauters, Joris Kinable, Pieter Smet, Wim Vancroonenburg, Greet Vanden Berghe, and Jannes Verstichel. The multi-mode resource-constrained multi-project scheduling problem. *Journal of Scheduling*, 19(3):271-283, 2016.

38 Kenneth Young, Thibaut Feydy, and Andreas Schutt. Constraint programming applied to the multi-skill project scheduling problem. In *Principles and Practice of Constraint Programming - 23rd International Conference, CP 2017, Melbourne, VIC, Australia, August 28 - September 1, 2017, Proceedings*, pages 308-317, August 2017. doi:10.1007/978-3-319-66158-2_20.
# A Appendix

## A.1 CP modeling comparisons

We ran new experiments in comparing the two different CP modeling we propose : the default one presented in section 3.2.2 and its upgrade that we called *CP-SmartPreemption*. Computation time are still 1800 s each like in experiments of section 5.1. Results are shown in 4. We interestingly note that *CP-SmartPreemption* is both faster and more efficient than

<table><caption>Table 5 Best solution and time to best solution for CP-Base and *CP-SmartPreemption*.</caption><thead><tr><th>ID</th><th>Features</th><th>CP-B</th><th>CP-B (s)</th><th>CPSmart</th><th>CPSmart (s)</th><th>ID</th><th>Features</th><th>CP-B</th><th>CP-B (s)</th><th>CPSmart</th><th>CPSmart (s)</th></tr></thead><tbody><tr><td>0</td><td>A</td><td>1080</td><td>1.09</td><td>1080</td><td>0.78</td><td>18</td><td>B G</td><td>3053</td><td>0.59</td><td>3053</td><td>0.20</td></tr><tr><td>1</td><td>A M</td><td>-</td><td>-</td><td>1080</td><td>8.20</td><td>19</td><td>B G M</td><td>3109</td><td>81</td><td>3109</td><td>3.20</td></tr><tr><td>2</td><td>A G</td><td>6534</td><td>0.99</td><td>6534</td><td>0.23</td><td>20</td><td>B P</td><td>2342</td><td>133</td><td>2340</td><td>178</td></tr><tr><td>3</td><td>A G M</td><td>-</td><td>-</td><td>6534</td><td>3.90</td><td>21</td><td>B M P</td><td>-</td><td>-</td><td>2390</td><td>196</td></tr><tr><td>4</td><td>A P</td><td>1104</td><td>2.21</td><td>1080</td><td>64.9</td><td>22</td><td>B G P</td><td>3063</td><td>106</td><td>3044</td><td>30.0</td></tr><tr><td>5</td><td>A M P</td><td>-</td><td>-</td><td>2495</td><td>24.4</td><td>23</td><td>B G M P</td><td>-</td><td>-</td><td>3180</td><td>198</td></tr><tr><td>6</td><td>A G P</td><td>6534</td><td>2.40</td><td>6534</td><td>62.3</td><td>24</td><td>B C</td><td>3196</td><td>154</td><td>3196</td><td>170</td></tr><tr><td>7</td><td>A G M P</td><td>-</td><td>-</td><td>6612</td><td>141</td><td>25</td><td>B C M</td><td>3272</td><td>127</td><td>3208</td><td>88</td></tr><tr><td>12</td><td>A C P</td><td>-</td><td>-</td><td>-</td><td>-</td><td>26</td><td>B C G</td><td>4949</td><td>0.52</td><td>4949</td><td>0.17</td></tr><tr><td>13</td><td>A C M</td><td>-</td><td>-</td><td>-</td><td>-</td><td>27</td><td>B C G M</td><td>4949</td><td>2.2</td><td>4949</td><td>0.80</td></tr><tr><td>14</td><td>A C G P</td><td>-</td><td>-</td><td>-</td><td>-</td><td>28</td><td>B C P</td><td>3252</td><td>1.32</td><td>3195</td><td>217</td></tr><tr><td>15</td><td>A C G M P</td><td>-</td><td>-</td><td>-</td><td>-</td><td>29</td><td>B C M P</td><td>-</td><td>-</td><td>3208</td><td>203</td></tr><tr><td>16</td><td>B</td><td>2339</td><td>106</td><td>2339</td><td>22.9</td><td>30</td><td>B C G P</td><td>4949</td><td>1.35</td><td>4383</td><td>59.0</td></tr><tr><td>17</td><td>B M</td><td>2359</td><td>193</td><td>2346</td><td>158</td><td>31</td><td>B C G M P</td><td>-</td><td>-</td><td>-</td><td>-</td></tr></tbody></table>

the default version. Notably, it manages to find feasible solution to 7 additional instances compared to basic CP formulation. It also succeeds in solving instance ID 30 way more efficiently than CP-Base or even the LNS that was presented in Table 3. As of today, the current modeling of *CP-SmartPreemption* still fails at the most complicated instances on the instance A that have the most operators/workers and number of task.

## A.2 Subproblem methods benchmark

We give more hindsight on the general choice done in the paper to use the mixing methods described in section 4.2.2. Naturally we expect that such a portfolio approach that picks a method from a pool of *N* neighborhood methods can be better in average than using only one method. To confirm this intuition, we run LNS methods using different parameters:

1. "Random selection" methods with parameter *f* taking values in [0.1, 0.2, 0.3, 0.4]

2. "Cut in equal parts" methods with parameter *c* taking values in [2, 3, 4, 5, 6]

3. "Mixing method" : portfolio of the previous 9 methods, chosen randomly at each iteration of LNS.

As a starting experiment, we are testing this methods on 5 classical RCPSP instances with 120 tasks, and run the solver 10 times per neighborhood methods and instances to get a better statistical confidence. We limited to 100 the number of iteration of LNS. Detailed results are in Table 6. It shows empirical confidence in using a mixing method, that achieved best performance. We also ran the same experiment on the instances ID=16,17 in Table 7 and got similar behaviour. Deeper hyper-parameters optimisation could be done for the more complex problems as PP-MS-MM-RCPSP/max-cal in a further research.
<table><caption>Table 6 Neighborhood methods performance on a few instances.</caption><thead><tr><th>Method</th><th>j1201_1</th><th>j1201_2</th><th>j1201_3</th><th>j1201_4</th><th>j1201_5</th></tr></thead><tbody><tr><td>RS(0.1)</td><td>116.6</td><td>131.3</td><td>138.6</td><td>106.5</td><td>131.3</td></tr><tr><td>RS(0.2)</td><td>112.2</td><td>129.8</td><td>136.8</td><td>106.0</td><td>128.3</td></tr><tr><td>RS(0.3)</td><td>110.4</td><td>128.4</td><td>134.3</td><td>105.4</td><td>120.4</td></tr><tr><td>RS(0.4)</td><td>108.9</td><td>118.0</td><td>132.6</td><td>103.3</td><td>116.7</td></tr><tr><td>Cut(2)</td><td><strong>105</strong></td><td><u>118.0</u></td><td><u>130</u></td><td><u>102.0</u></td><td><strong>113.0</strong></td></tr><tr><td>Cut(3)</td><td>107</td><td>118.0</td><td>135</td><td><strong>101.0</strong></td><td>118.0</td></tr><tr><td>Cut(4)</td><td>108</td><td>126.0</td><td>131</td><td>106.0</td><td>132.0</td></tr><tr><td>Cut(5)</td><td>111</td><td>128.0</td><td>136</td><td>106.0</td><td>130.0</td></tr><tr><td>Cut(6)</td><td>111</td><td>129.0</td><td>137</td><td>106.0</td><td>130.0</td></tr><tr><td>Mixing</td><td><u>106.5</u></td><td><strong>115.0</strong></td><td><strong>128.3</strong></td><td><u>102.0</u></td><td><strong>113.0</strong></td></tr></tbody></table>

<table><caption>Table 7 Neighborhood methods performance on an industrial use case (ID=16).</caption><thead><tr><th>Method</th><th>ID 16 (B)</th><th>ID 17 (BM)</th></tr></thead><tbody><tr><td>RS(0.1)</td><td>2698.0</td><td>2745.0</td></tr><tr><td>RS(0.2)</td><td>2659.0</td><td>2675.0</td></tr><tr><td>RS(0.3)</td><td>2534.0</td><td>2635.0</td></tr><tr><td>RS(0.4)</td><td>2489.0</td><td>2662.0</td></tr><tr><td>Cut(2)</td><td>2344.0</td><td>2628.0</td></tr><tr><td>Cut(3)</td><td><u>2343.0</u></td><td><strong>2436.0</strong></td></tr><tr><td>Cut(4)</td><td>2346.0</td><td>2480.0</td></tr><tr><td>Cut(5)</td><td>2346.0</td><td>2486.0</td></tr><tr><td>Cut(6)</td><td>2349.0</td><td>2776.0</td></tr><tr><td>Mixing</td><td><strong>2342</strong></td><td><u>2460.0</u></td></tr></tbody></table>

## A.3 MSLIB Experiments

To better understand the results we plot the overcost of the LNS method as a function of different parameters used to generate MSLIB4 instances in Figure 4. We theorize that increasing #R (number of resources) has a negative impact on LNS performance. It also highlights that LNS is underperforming only on instances where SS = 0 and SSα = 0, defined as skill strength and skill strength variability. As future work we want to continue testing on MSLIB benchmark to possibly improve the LNS method. Detailed results comparing the GA and LNS approach over the tested instances are depicted in Tables 8 and 9.

<table><caption>Table 8 Best results over subset of MSLIB4 instances (404 instance).</caption><thead><tr><th>Algorithms</th><th>LNS (ours)</th><th>GA 2021</th><th>Equal</th></tr></thead><tbody><tr><td>#Best Results</td><td>90</td><td>164</td><td>150</td></tr><tr><td>Mean Improvement when better solution than GA</td><td>2.9%</td><td>-</td><td>-</td></tr><tr><td>Mean Degradation when worse solution than GA</td><td>19.2%</td><td>-</td><td>-</td></tr></tbody></table>
Figure 4 Relative overcost to GA of LNS method, function on different features of MSLIB4 instances.

<figure><img src="image_299.png" /></figure>

Table 9 Detailed results comparing GA (from [33]) and our methods, function on instance parameters. The GA, Equal, LNS columns counts the number of instance of some given parameters where each of the algorithms got the best results. (LNS-GA)/GA column store the average overcost of LNS method over the subset of instances described by SP, SS, RA.

<table><thead><tr><th>SP</th><th>SS</th><th>RA</th><th>NbRun</th><th>GA</th><th>Equal</th><th>LNS</th><th>Mean (LNS-GA)/GA (%)</th></tr></thead><tbody><tr><td>0.1</td><td>0.0</td><td>0.0</td><td>79</td><td>41</td><td>34</td><td>4</td><td>12.4</td></tr><tr><td>-</td><td>-</td><td>0.2</td><td>40</td><td>5</td><td>17</td><td>18</td><td>1.89</td></tr><tr><td>-</td><td>-</td><td>0.3</td><td>33</td><td>33</td><td>0</td><td>0</td><td>18.6</td></tr><tr><td>-</td><td>-</td><td>0.4</td><td>25</td><td>9</td><td>9</td><td>7</td><td>5.59</td></tr><tr><td>-</td><td>-</td><td>0.5</td><td>66</td><td>44</td><td>11</td><td>11</td><td>12.4</td></tr><tr><td>-</td><td>-</td><td>0.6</td><td>15</td><td>5</td><td>7</td><td>3</td><td>2.86</td></tr><tr><td>-</td><td>-</td><td>0.7</td><td>18</td><td>11</td><td>5</td><td>2</td><td>9.21</td></tr><tr><td>-</td><td>-</td><td>0.8</td><td>22</td><td>13</td><td>7</td><td>2</td><td>8.04</td></tr><tr><td>-</td><td>-</td><td>0.9</td><td>3</td><td>3</td><td>0</td><td>0</td><td>6.78</td></tr><tr><td>-</td><td>-</td><td>1.0</td><td>28</td><td>0</td><td>18</td><td>10</td><td>-1.01</td></tr><tr><td>-</td><td>0.1</td><td>0.0</td><td>24</td><td>0</td><td>14</td><td>10</td><td>-1.06</td></tr><tr><td>-</td><td>-</td><td>0.2</td><td>6</td><td>0</td><td>2</td><td>4</td><td>-2.71</td></tr><tr><td>-</td><td>-</td><td>0.3</td><td>3</td><td>0</td><td>1</td><td>2</td><td>-3.60</td></tr><tr><td>-</td><td>-</td><td>0.4</td><td>1</td><td>0</td><td>0</td><td>1</td><td>-4.88</td></tr><tr><td>-</td><td>-</td><td>0.5</td><td>1</td><td>0</td><td>0</td><td>1</td><td>-1.67</td></tr><tr><td>-</td><td>-</td><td>0.6</td><td>1</td><td>0</td><td>1</td><td>0</td><td>0.0</td></tr><tr><td>-</td><td>-</td><td>0.7</td><td>2</td><td>0</td><td>1</td><td>1</td><td>-0.76</td></tr><tr><td>-</td><td>0.2</td><td>0.0</td><td>18</td><td>0</td><td>12</td><td>6</td><td>-1.00</td></tr><tr><td>-</td><td>-</td><td>0.2</td><td>6</td><td>0</td><td>2</td><td>4</td><td>-2.67</td></tr><tr><td>-</td><td>-</td><td>0.3</td><td>1</td><td>0</td><td>1</td><td>0</td><td>0.0</td></tr><tr><td>-</td><td>-</td><td>0.4</td><td>2</td><td>0</td><td>1</td><td>1</td><td>-2.38</td></tr><tr><td>-</td><td>0.3</td><td>0.0</td><td>5</td><td>0</td><td>3</td><td>2</td><td>-1.61</td></tr><tr><td>-</td><td>-</td><td>0.2</td><td>1</td><td>0</td><td>0</td><td>1</td><td>-2.63</td></tr><tr><td>-</td><td>0.4</td><td>0.0</td><td>2</td><td>0</td><td>2</td><td>0</td><td>0.0</td></tr><tr><td>-</td><td>0.5</td><td>-</td><td>2</td><td>0</td><td>2</td><td>0</td><td>0.0</td></tr></tbody></table>
<figure><img src="image_300.png" /></figure>
# Assembly Line Preliminary Design Optimization for an Aircraft

Stéphanie Roussel<fnref n="1" /> ⬒ ID

ONERA, ONERA DTIS, Toulouse, Université de Toulouse, France

Thomas Polacsek ⬒ ID

ONERA, ONERA DTIS, Toulouse, Université de Toulouse, France

Anouck Chan ⬒ ID

ONERA, ONERA DTIS, Toulouse, Université de Toulouse, France

## Abstract

In the aeronautics industry, each aircraft family has a dedicated manufacturing system. This system is classically designed once the aircraft design is completely finished, which might lead to poor performance. To mitigate this issue, a strategy is to take into account the production system as early as possible in the aircraft design process. In this work, we define the Assembly Line Preliminary Design Problem, which consists in defining, for a given aircraft design, the best assembly line layout and the type and number of machines equipping each workstation. We propose a Constraint Programming encoding for that problem, along with an algorithm based on epsilon constraint for exploring the set of Pareto solutions. We present experiments run on a set of real industrial data. The results show that the approach is promising and offers support to experts in order to compare aircraft designs with each other.

**2012 ACM Subject Classification** Applied computing → Computer-aided manufacturing

**Keywords and phrases** Assembly line design, Constraint Programming, Multi-objective, Industry 4.0

**Digital Object Identifier** 10.4230/LIPIcs.CP.2023.32

**Supplementary Material** Dataset (Assembly Line Problems): https://github.com/stephroussel/assemblyLine

## 1 Introduction

The aeronautics industry is highly specialized. It requires significant investments in research to design and manufacture a new aircraft. Each new aircraft family, consisting of several models within a given size range, requires the development of a dedicated industrial system. The latter includes the manufacturing and assembly processes, supply chain, tooling and facilities necessary to produce an aircraft at scale. The design of such a system typically involves several steps of high-level objectives refinements, which can take several years and requires collaboration between various stakeholders (aircraft manufacturers, suppliers, and regulatory agencies, etc.). Once established, the system can be used for the entire aircraft family lifecycle, which can last several decades. While aircraft manufacturers can generally leverage existing resources such as buildings and infrastructure, they have to develop new manufacturing and assembly lines to meet the specific requirements of each aircraft family.

An industrial system has its own set of requirements and is designed so as to optimize several criteria. For instance, assembly lines must comply to organizational constraints and ergonomics recommendations while minimizing the overall investment cost. As the manufacturing system depends highly on the aircraft choices (e.g. the aircraft material such as carbon fibers, metal, etc.), the manufacturing system design is generally thought after the

<figure><img src="image_301.png" /></figure>

<footnote n="1">¹ Corresponding author</footnote>
aircraft design. This might result in a low performance of the latter or even the impossibility to meet some requirements. For instance, choosing carbon fibers on the aircraft side requires buying autoclaves on the industrial system side, with an associated cost that exceeds a given investment budget.

Therefore, a nowadays trend in aeronautics is to take into account the industrial system design as early as possible in the aircraft design process ([18, 22]). A way to do so, and the approach we follow in this paper, is to design the best preliminary industrial system for an aircraft design that is not validated yet and assess its performance. This allows to estimate whether expected industrial objectives can be met for a given aircraft. Moreover, if several aircraft designs are candidates, this also allows comparison between them with respect to their industrial suitability. Finally, such an approach offers support for making trade-offs between the aircraft and the manufacturing. More precisely, a modification of the aircraft that downgrades its performance might allow to decrease significantly the cost associated to the production system. For instance, a change in the aircraft shape might increase a little bit its fuel consumption but reduce significantly the time required to build it because the new shape allows to use robots for making junctions.

This work presents contributions produced in the context of an industrial project with an aircraft manufacturer for supporting trade-offs between an aircraft and its industrial system. The objective of this project was to focus on one element of the industrial system, namely the *assembly line*. An aircraft assembly line is responsible for bringing together all of the components of the aircraft, including the fuselage, wings, engines and interiors, into a complete and fully functional aircraft. The assembly lines we consider are pulsed and composed of several workstations that are equipped with specific machines. Assembly lines must be designed to be flexible enough to accommodate changes in the manufacturing process as the aircraft family evolves over time. Moreover, because of the international context, production rates can vary a lot within the aircraft lifecycle. It is therefore essential to choose an aircraft design and an assembly line design that are compatible with such evolution. Other elements of the industrial system, *e.g.* the supply chain, depend also highly from the aircraft design but they were out of scope of our project. Nevertheless, the assembly line by itself represents a quite expensive part of the industrial system. In fact, when an aircraft design allows to have one workstation (or one machine) less than another one, it can correspond to costs of several million euros.

In the preliminary stages of designing an assembly line, aircraft manufacturers can use the PERT (Program Evaluation and Review Technique). PERT is a network diagram that shows the sequence of tasks and their dependencies. It allows to identify critical paths and potential bottlenecks. However, PERT is a very optimistic view of the building process. In fact, because it does not take resources into account, it assumes that all tasks that are not dependent can be performed in parallel. In order to have a more realistic assessment on the assembly line, industrial architects design workstations and equip them, assign tasks to workstations and then compute a fine-grain line balancing that takes into account workers and their profiles. However, such an iterative process can be quite heavy to put in place for each aircraft design candidate.

In this paper, we present an aid-decision tool that we have developed in the context of the previously mentioned project. For a given aircraft, the tool returns the best preliminary assembly line designs with respect to several criteria. These criteria deal with the assembly line layout (workstations and associated equipment), its cost and also with the production duration. We call this problem the Assembly Line Preliminary Design Problem, which can be seen as a Resource-Constrained Project Scheduling Problem (RCPSP) with additional
constraints ([9]). The tool relies on a Constraint Programming encoding of the problem. In fact, Constraint Programming is particularly suited for handling resource and temporal constraints inherent in the assembly line design problem ([16]). In order to explore the Pareto solutions set, we define an epsilon constraint based algorithm ([4]). It basically consists in computing bounds for each criteria, choosing an exploration order and a number of points in the front, and next running a mono-criterion algorithm to generate Pareto solutions.

Our contribution can be summarized as follows:

- we present a new problem, namely the Assembly Line Preliminary Design Problem, and formalize associated constraints and criteria;

- we define a Constraint Programming encoding for this problem;

- we develop an epsilon constraint based algorithmic approach for exploring the Pareto front;

- we present experiments on a set of real data coming from an aircraft manufacturer;

- we make those data available for the community.

The paper is structured as follows. In Section 2, we describe related works. Section 3 is dedicated to an informal description of the problem and some prior work for eliciting constraints and criteria in the case of a preliminary design. In Section 4, we formally model the Assembly Line Preliminary Design Problem and propose an associated Constraint Programming encoding. We detail the algorithmic approach for exploring the Pareto solutions set in Section 5. In Section 6, we describe the real world benchmarks we have used and discuss the results obtained with the paper's approach. Finally, we conclude on future works in Section 7.

## 2 Related Works

A vast body of literature exists on the production and line balancing problems research topic ([5, 8]). The *Simple Assembly Line Balancing Problem* (SALBP) consists in assigning tasks to workstations on a single straight line with respect to some precedence constraint and to a fixed workstation capacity. Such a problem can be seen as a bin packing problem with additional constraints due to precedence. Depending on the objective function, there are two classical extensions: SALBP-1 in which the number of workstations is minimized for a given production rate and SALBP-2 in which the production rate is maximized for a given number of stations.

**Assembly Line Design and Constraint Programming.** When considering the design of assembly lines, additional features are generally taken into account, in particular resources required for the task execution ([24]). Resources are generally used to model machines or equipment that have to be positioned on the assembly line but also to model workers performing assembly operations. Constraint Programming (CP) models have been developed for solving various assembly line problems, and specifically problems in which resources are involved. In [11], the authors present encodings for solving SALBP-1, SALBP-2 and the Task Assignment and Equipment Selection Problem (TAESP). In the latter, one type of equipment has to be assigned to each station and the objective is to minimize the total associated cost. In [2], the authors use CP for solving the Resource-Constrained Assembly Line Balancing Problem (RCALBP) in which there can be several resources candidates for performing each task. The objective is the cycle time minimization, i.e. the time spent in each workstation of the line. This work has been extended in [1] into the Generalized RCALBP, in which several
modes with different resources consumption are considered for each task. These modes are expressed through conjunction normal form formulas that are handled directly by the CP solver through the use of AND and OR operators.

**Multi-criteria Line Design.** As described in [24], the assembly design problem often comes to the optimization of several criteria. In [29], the authors consider the minimization of the number of workstations, the minimization of cycle time and maximization of workload smoothness and of work relatedness. Each criterion is studied independently, except for the last two. The authors of [13] consider the problem of equipping workstations with respect to some exclusion and compatibility constraints between pieces of equipment. They study several criteria: minimization of the investment cost, maximization of the line throughput rate, minimize the occupied space and minimize the workstations complexity. In [19], the assembly line considered is mixed-model, meaning that similar models of a product are produced on the same line. Their objective is to minimize the idle time of each model on the line and minimize the total equipment cost. More recently, in [20], the authors optimize the idle time and the unit product costs. All those works mainly use evolutionary algorithms for computing a set of non-dominated solutions. An extensive comparison between 34 algorithms is provided in [20]. In [15], the authors address the two-sided assembly line balancing problem with multi-operator workstations with respect to several criteria: minimization of the number of mated workstation, minimization of the number of workstation and minimization of the number of workers. Those criteria are optimized sequentially. A CP and an ILP (Integer Linear Programming) approaches, that both consider the time at which tasks start and end, are proposed.

**Line Balancing in the Aeronautical Industry.** In the context of assembly line balancing for the aeronautical industry, literature is less extensive. The authors of [14] consider the efficiency of labor utilization, specifically for tasks that consist in preparing the aircraft for assembly, for instance when it arrives on a workstation. In a generic context of producing in low-volume, which applies for the aircraft industry, the article [6] assigns workers to station and start times to each task so as to minimize the costs of total labor and inventory costs. The authors of [10] consider the local order scheduling for mixed-model assembly lines. In [7] and more recently in [26], exact (Mixed Integer Linear Programming and CP) and heuristic methods are studied for solving a mixed-model assembly line problem in which tasks have been assigned to workstations and the objective is to minimize the overall number of operators. In [3], the authors compare CP and ILP approaches for assigning workers to tasks so as to minimize the overall makespan and while ensuring that some ergonomics constraints are satisfied. Both approaches consider not only the assignment of worker to station but also the precise scheduling of all activities. In [23], a scheduling tool is developed for bridging the gap between aircraft design and aircraft manufacturing. Such a tool relies on a local search approach and is compared with a CP encoding. The problem consists in minimizing the makespan. As done in the present paper, it considers aircraft zones that limit the number of workers that can work simultaneously. However, it does not take into account machines on the assembly line or zone neutralization and does not allow to optimize the takt-time.

Finally, [12] is a previous work in which we present results we have obtained on the same industrial project. In that work, we focus on the elicitation part and present a Goal-Oriented Requirements Engineering methodology that allowed us to obtain constraints and criteria associated with an assembly line performance. We indicate that an operations research approach was used to explore the set of assembly line solutions but this approach is not formally described and only returns a few assembly line designs solutions.
**RCPSP and Multi-Objective Optimization.** As explained in [26], the assembly line balancing problems in the aeronautical industry can be modeled into RCPSPs with some additional constraints, such as temporal constraints or incompatibility constraints. Constraint Programming provides a expressive language for such constraints, along with efficient solvers [16].

Some works of the literature address the Multi-Objective Constrained Optimization Problems (MO-COP). In [25], the authors propose a multi-objective lower bound set that allows to detect inconsistency of a multi-objective problem. The authors of [21] define an interactive algorithm for MO-COP, that gradually refines a set of Pareto solutions by exploring regions chosen by the user in which a Pareto front might exist. The specific case of Multi-Objective in RCPSP has also received attention, as shown in [4].

# 3 Preliminary Assembly Line Design Description

In this section, we describe the assembly lines we have considered, namely pulse assembly lines. Then, we briefly present the work that we have done prior to optimization for eliciting constraints and criteria relevant for a preliminary assembly line design stage.

## 3.1 Pulse Assembly Line

The assembly line we consider in this work is a pulse line composed of several workstations. In aeronautics, pulse assembly lines are common solutions when it comes to high production rates.

The layout of pulse assembly lines we consider is a basic straight line. Each workstation in the line is responsible for performing a specific set of tasks in the assembly process, and the product being assembled moves from one workstation to the next until it is complete. When an aircraft enters the line, it goes to the first workstation. After a duration called *takt-time*<fnref n="2" />, the aircraft goes to the second workstation and a new aircraft enters the first workstation. The aircraft that was on the last workstation is done and leaves the assembly line. Tasks cannot be performed when the aircraft is moving. This implies that all tasks must be assigned to one unique workstation.

The total time that the aircraft remains on the assembly line is called *leadtime*. The leadtime is equal to the multiplication of the takt-time by the number of workstations. Note that the range of values for the takt-time of a given assembly line can be quite large. In fact, when a new line is opened, the production can start with a few aircraft per month, which corresponds to large takt-times. Depending on business considerations, the production rate can reach several dozen aircraft per month (takt-time equal to a few hours).

The set of operations performed on each workstation is always the same and they are performed in the same order for all aircraft instances. Moreover, each workstation is equipped with a set of jigs and tools, or more generally machines, that are dedicated to it. Such machines can be heavy installations, such as scaffolds for supporting the aircraft, dedicated robots or some smaller tools, such as dedicated drilling devices.

Assembly lines we consider in this paper are single-model lines as there is one product to manufacture. In practice, there can be small variants between all the aircraft present on an assembly line but such differences are not taken into account in the preliminary stage design.

<footnote n="2">² Takt-time is generally called cycle time in the literature. However, cycle time sometimes refer to another duration in the aircraft industry.</footnote>
## 3.2 Constraints and Criteria Elicitation

Prior to the tool implementation, a major part of the project we are involved in has been dedicated to the elicitation of preliminary assembly line performance constraints and criteria. To do so, we have followed a goal-oriented requirements engineering based that is described in details in [12].

As expected, one of the criteria for the assembly line design is to minimize the associated investment cost. The costs of designing an aircraft assembly line can vary significantly depending on a variety of factors, including the size and complexity of the aircraft family, the level of automation and technology required, and the specific requirements of the production process. In the case of the preliminary assembly line design, we focus on two types of cost.

**Minimize the number of workstations.** We first consider the workstations cost. In fact, each workstation occupies a floor area in factories and has therefore an associated cost. In reality, the latter depends on the size of machines equipping the workstation and on the volume of the aircraft parts manipulated. In the assembly lines we consider in this work, it is reasonable to consider that the parts manipulated are comparable and therefore that the workstations cost is uniform. Therefore, we try to minimize the number of workstations in the assembly line.

**Minimize the number of machines.** The second type of cost is the machines and tools cost. In fact, to build an assembly line, specialized tooling and equipment must be designed and installed. This may include robots, fixtures and specialized machinery, which can be very costly. As we are in a preliminary design stage, we only consider the sizing machines and tools of the line. Cost range of such machines can be quite large as it depends on the machine features. In this paper, we focus only on highly expensive machines. As each machine is really specific, there is almost one criterion for each type of machine. In this project, we did not have any cost information nor ranking between the machines. In fact, costs can be quite complex to compute (because it implies not only buying the machine but also maintaining it). We have therefore decided to simplify the machines cost and consider only the total number of machines minimization. However, the number of machines of each type should be carefully examined when looking into details solutions in the front.

**Machines assignment and incompatibility.** Machines and tools we consider can be assigned to exactly one workstation. As they can be voluminous, there exists incompatibility between some machines types. Therefore, we have an exclusion relationship between machine.

**Minimize the leadtime.** Leadtime reduction is a critical goal for any manufacturing process. In fact, manufacturers can reduce the amount of cash that is tied up in work-in-progress inventory, which can improve cash flow and financial performance.

**Minimize the takt-time.** As expected, the higher the production rate, the more efficient the assembly line. This corresponds to the takt-time minimization. On that point, we can note that a target takt-time is generally fixed when designing a new aircraft. Nevertheless, it is essential to assess the best takt-time that is compatible with a given aircraft design, in order to prepare future production rate increase.
**Aircraft zones capacity.** The space available for workers inside the aircraft is a major constraint. Workers may need to operate in tight spaces or confined areas to access parts of the aircraft being assembled, which can increase the errors or accidents risk. This is particularly the case for tasks that require precision and accuracy, such as drilling, fastening or installing electrical components. To model that, we consider that the aircraft is divided into zones, and that each zone has a capacity representing the maximum number of workers that can perform assembly tasks simultaneously.

**Zone neutralization.** Because of safety constraints or some assembly task nature, some assembly task can neutralize aircraft zones, meaning that it prevents any other task to be performed in those neutralized zones simultaneously. For instance, if a task requires the removal of a temporary floor during its execution, it prevents access to the associated zone.

**Task features.**

- There is one alternative for each task. Each alternative is characterized by types of machine used, zones occupied and neutralized.

- Duration of tasks is fixed.

- It is not possible to interrupt tasks during their execution (no preemption). However, as we are in preliminary design, tasks we consider are in fact macro ones. Therefore, it is possible for tasks that do not consume any machine to overlap on several workstations.

# 4 Assembly Line Preliminary Design Problem Definition

In this section, we formally define the Assembly Line Preliminary Design Problem (ALPDP) and present a Constraint Programming encoding for it.

## 4.1 Problem Definition

**Assembly Line Preliminary Design Problem.** An Assembly Line Preliminary Design Problem (ALPDP) is formally defined by a tuple $\langle \mathcal{Z}, \mathcal{S}, \mathcal{E}, \mathcal{T}, \mathcal{P} \rangle$ where:

- $\mathcal{Z}$ is the set of zones of the aircraft in which workers perform assembly activities. Each zone $z \in \mathcal{Z}$ has a capacity, denoted $capa_z$, that represents the number of workers that can work simultaneously in zone $z$;

- $\mathcal{S}$ is the set of machine skills that are required by assembly activities. Each skill $s \in \mathcal{S}$ represents a type of jig and tool (i.e. a type of machine) used for the aircraft assembly;

- $\mathcal{E}$ is a symmetric relation included in $\mathcal{S}^2$ that represents skills of machines that are incompatible and cannot belong to the same workstation. Note that for a given skill $s \in \mathcal{S}$, it is possible to have $(s, s) \in \mathcal{E}$: this indicates that two machines with the same skill $s$ cannot be assigned to the same workstation;

- $\mathcal{T}$ is the set of assembly activities that must be performed. $\mathcal{T}$ can be partitioned into two subsets: a set of atomic activities denoted $\mathcal{A}$ and a set of composite activities denoted $\mathcal{C}$. Intuitively, composite activities are a group of atomic activities that allow to write precedence between activities in a compact way. They do not have a fixed duration and do not consume any resource (zone or machine);

- for each composite activity $c \in \mathcal{C}$, $\mathcal{A}_c$ is a subset of $\mathcal{A}$ that represent all the atomic tasks that are part of composite task $c$.

*Assumption.* An atomic task is the child of at most one composite activity. Formally, for $c_1$ and $c_2$ in $\mathcal{C}^2$ such that $c_1 \neq c_2$, we have $\mathcal{A}_{c_1} \cap \mathcal{A}_{c_2} = \emptyset$;

- for each atomic activity $a \in \mathcal{A}$, we consider the following features:
- dur$_{a} \in \mathbb{N}^{+}$ is the duration of a;

- $S_a \subseteq S$ is the set of skills required by $a$. Note that, for a given skill $s$, an activity requires at most one machine with that skill;

- for each zone $z \in \mathcal{Z}$, $occ_{a,z} \in \mathbb{N}^{+}$ is the number of places in $z$ required by activity $a$;

- $\mathcal{Z}_a^{\text{neutr}} \subseteq \mathcal{Z}$ is the set of zones that are neutralized by $a$. If zone $z'$ belongs to $\mathcal{Z}_a^{\text{neutr}}$, then for all activity $a'$ that occupies $z'$ (i.e. such that $occ_{a',z'} > 0$), $a$ and $a'$ cannot temporally overlap.

*Assumption.* For a given activity, the set of zones it occupies and the set of zones it neutralizes have an empty intersection;

- $\mathcal{P}$ is a relation in $\mathcal{T}^2$ that represents precedence requirements between activities. More precisely, $(t, t') \in \mathcal{P}$ if $t$ must be finished before $t'$ can start.

*Assumption 1.* There does not exist precedence between composite activities and their children, i.e. for all $c \in \mathcal{C}$ and $a \in \mathcal{A}_c$, $(a, c) \notin \mathcal{P}$ and $(c, a) \notin \mathcal{P}$.

*Assumption 2.* The directed graph induced by relation $\mathcal{P}$ is acyclic. To formally build such a graph, the first step is to create a node $n_a$ for each atomic activity $a$ and two nodes $s_c$ and $e_c$ for each composite activity $c \in \mathcal{C}$ that respectively represent the start and end of $c$. Then, arcs are added as follows. For each precedence $(a_1, a_2) \in \mathcal{A}^2$, an arc is added between $n_{a_1}$ and $n_{a_2}$. For each precedence $(c_1, c_2) \in \mathcal{C}^2$, an arc is added between $e_{c_1}$ and $s_{c_2}$. For each precedence $(a, c) \in \mathcal{A} \times \mathcal{C}$ (resp. $(c, a) \in \mathcal{C} \times \mathcal{A}$), an arc is added between $n_a$ and $s_c$ (resp. between $e_c$ and $n_a$). Note that this graph construction is similar to the one proposed in [23].

**Alternative.** An alternative for a ALPDP $\langle \mathcal{Z}, \mathcal{S}, \mathcal{E}, \mathcal{T}, \mathcal{P} \rangle$ is described through the following elements:

- a set of workstations $\mathcal{W}$ and for each workstation $w \in \mathcal{W}$, a start date and an end date;

- a takt-time, denoted *takt*, that represents the duration the aircraft stays on each workstation;

- for each skill $s$ and for each workstation $w$, a number of machines with skill $s$ on workstation $w$;

- for each activity $t \in \mathcal{T}$, a start date and an end date.

The leadtime of a solution, denoted $T_{\max}$, is the total time the aircraft spends on workstations. Formally, $T_{\max} = |\mathcal{W}| \cdot \textit{takt}$.

**Solution.** A solution for an ALPDP $\langle \mathcal{Z}, \mathcal{S}, \mathcal{E}, \mathcal{T}, \mathcal{P} \rangle$ is an alternative in which the following constraints are satisfied. Note that we only give an informal definition of constraints here, as the formal definition is detailed through the Constraint Programming encoding later in the paper.

- Start dates and end dates of workstations must represent a consistent assembly line (no temporal hole and no overlap between successive workstations) and consistent with the *takt* value (duration of a workstation is equal to *takt*);

- start dates and end dates of activities are consistent with duration (for atomic activities), with children start and end dates (for composite activities) and with precedence relationship;

- capacity of zones is never exceeded and zone neutralization is respected;

- the number of machines of each skill in each workstation allows the associated activities to be performed;

- skills exclusion is respected;

- atomic activities that consume at least one machine cannot overlap on successive workstations.
<figure><img src="image_302.png" /><figcaption>(a) Activities description.</figcaption></figure>

<figure><img src="image_303.png" /><figcaption>(b) Graph induced by precedence relation $\mathcal{P}$.</figcaption></figure>

Figure 1 Toy example illustration.

**Criteria.** According to the criteria elicitation phase, we consider the following criteria: minimize the takt-time, minimize the leadtime value, minimize the total number of machines on the assembly line and minimize the number of workstations. As described earlier, we do not have any insight on how to aggregate this criteria together.

**► Example 1.** Figure 1 illustrates an ALPDP toy example. In this example, the set of zones $\mathcal{Z}$ is composed of two zones $z_1$ and $z_2$ that respectively have a capacity equal to 1 and 2. The set of skills $\mathcal{S}$ contains three skills, $s_1$, $s_2$ and $s_3$, and we consider that $s_1$ and $s_3$ exclude each other, and so do skills $s_2$ and $s_3$ (i.e. $\mathcal{E} = \{(s_1, s_3), (s_3, s_1), (s_2, s_3), (s_3, s_2)\}$). The set of activities $\mathcal{T}$ contains one composite activity $c_1$ and five atomic activities $a_i$ with $i \in [1..5]$. Composite activity $c_1$ encompasses activities $a_3$ and $a_4$.

Table 1a describes atomic activities features. The X symbol indicates that an activity neutralizes a zone. The ✓ symbol indicates that an activity requires a skill. For instance, the second line of the table indicates that activity $a_2$: has a duration equal to 1, neutralizes zone $z_1$, occupies one place in zone $z_2$, and requires skill $s_2$.

Precedence relation is $\mathcal{P} = \{(a_1, a_2), (a_2, a_5), (c_1, a_5), (a_3, a_4)\}$, as illustrated on Figure 1b$^3$. Figure 2a illustrates a solution, denoted $sol_{3machines}$ for this ALPDP. In this solution, there are two workstations, $w_1$ and $w_2$. The first workstation contains machine $m_1$ with skill $s_1$ and machine $m_2$ with skill $s_2$. Workstation $w_2$ contains machine $m_3$ with skill $s_3$. Start dates of activities are respectively 0, 2, 2, 3, 4 for $a_1, \dots, a_5$. Composite activity $c_1$ starts with its earliest child (here $a_3$ at time 2) and ends with its latest (here $a_4$ at time 4). In this solution, the takt is equal to 4 and the leadtime equal to 8.

A second solution, $sol_{4machines}$, is illustrated on Figure 2b. In this solution, there are two machines with skill $s_1$ in the first workstation. Note that for this solution, composite activity $c_1$ is not continuous in the sense that there is a temporal hole between its children activities. This solution allows to reach a takt equal to 3 and a leadtime equal to 6.

We can notice that $sol_{3machines}$ is better than $sol_{4machines}$ with respect to the total number of machines criteria. Both solutions are equivalent for the number of stations. However, $sol_{4machines}$ is better than $sol_{3machines}$ with respect to the takt and to the leadtime. Thus, no solution dominates the other on all criteria.

<footnote n="3">³ Note that we do not represent explicitly start and end activities associated with $c_1$.</footnote>