Modelling and Solving the Precedence-Constrained Minimum-Cost Arborescence Problem with Waiting-Times

A polynomial-size mixed integer linear programming model for the Precedence-Constrained Minimum-Cost Ar-borescence Problem with Waiting-Times was recently proposed in the literature, that uses a smaller number of variables and constraints compared to previously proposed polynomial-size models. In this work, we extend this model with constraint programming constructs to further enhance its performance. An extensive computational study support that modern constraint programming solvers are the best tool available at solving the models proposed. Several improvements to state-of-the-art results are finally renorted.


I. INTRODUCTION
T HE Minimum-Cost Arborescence (MCA) problem in- volves finding a directed minimum-cost spanning tree, rooted at vertex r, in a given input directed graph.Jack Edmonds [1], and Yoeng-Jin Chu and Tseng-Hong Liu [2] independently introduced the first polynomial time algorithm for solving the problem.Gabow and Tarjan [3] improved the running time of the algorithm by using disjoint-sets and a special implementation of Fibonacci heaps.
Several variations of the MCA problem with different objective function and/or constraints were introduced in the literature since its introduction.Given a finite resource associated with each vertex in the input graph, the Resource-Constrained Minimum-Weight Arborescence problem [4] is an N P-hard problem which asks to find an arborescence with minimum total cost where the sum of the costs of outgoing arcs from each vertex is at most equal to the resource of that vertex.Given an integer Q and non-negative integer vertex demand q j associated with each vertex, the Capacitated Minimum Spanning Tree problem [5] is an N P-hard problem which asks to find a directed minimum spanning tree rooted at r, such that the sum of the weights of the vertices in any subtree off the root is at most Q.Given a weighted directed acyclic graph with each vertex having a specified color from a set of colors, the Maximum Colorful Arborescence problem [6] is an N P-hard problem which asks to find an arborescence of maximum weight, in which no color appears more than once.Given an integer rank associated with each vertex, the Restricted Fathers Tree problem [7] asks to find a minimumcost arborescence rooted at r, such that the path between each vertex and the root contains only vertices with same rank or higher.
Constraint programming (CP) is paradigm for solving combinatorial problems by representing them as constraint satisfaction problems (CSP) [8].A CSP is represented as a set of variables each with a defined domain of values, and a set of relations/constraints on the subsets of these variables.A CP solver takes a CSP and finds an assignment to all the variables that satisfies the constraints, and can also extend the problem to finding optimal solutions according to an optimization criteria.A CP solver searches the solution space systematically using a branch-and-bound algorithm with inference techniques which consists of propagating the information contained in one constraint to the neighboring constraints.Such techniques reduce the size of the solution space that needs to be explored [9].CP has been used to solve a wide range of problems in the literature.Hande [10] proposed a CP model for the Open Vehicle Routing problem with Heterogeneous Vehicle Fleet (HFOVRP).In [10] the CP model is compared with a mixed-integer linear programming (MILP) model of the HFOVRP, and they showed that the CP model is effective for providing good-quality solutions for small-sized instances of the HFOVRP in short computational times compared to the MILP model.Kasapidis et al. [11] presented a MILP model and a CP model for the Multi-Resource Flexible Job-Shop Scheduling problem with Arbitrary Precedence Graphs.The computational experiments conducted in [11] has shown that the CP model is more effective and achieves the best results compared to the MILP model, although more timeconsuming on some instances.Kirac et al. [12] proposed a CP approach for solving the Team Orienteering problem with Time Windows and Mandatory Visits, and they showed that the CP-based approach finds 99 of the best-known solutions and explores 64 new best-known solutions for the benchmark instances.Kizilay et al. [13] proposed a novel CP model for the Mixed-Blocking Permutation Flow Shop Scheduling problem with Batch Delivery that minimizes the total tardiness and batch cost.The results of their study has shown that due to the complexity of the problem, the developed CP model can solve only small-sized instances in reasonable computational time.
Montemanni and Dell'Amico [14] proposed a CP model for the Parallel Drone Scheduling Traveling Salesman problem, and showed that by exploiting multi-threading computation, the method was able to optimally solve all the instances considered in the literature.
The Precedence-Constrained Minimum-Cost Arborescence (PCMCA) problem is an N P-hard problem [15] that was first introduced by Dell'Amico et al. [16].The PCMCA problem is an extension to the MCA problem, in which precedence constraints must be satisfied as follows.Given a set R of ordered pairs of vertices, then for each precedence relationship (s, t) ∈ R, a path in the solution which covers both s and t, must visit vertex s before visiting vertex t.The objective is to find an arborescence of minimum total cost satisfying the precedence constraints.The PCMCA problem has applications in the design of commodity distribution networks where certain paths are not allowed in the network due to logistical constraints [16].Several MILP models of the problem were proposed in [15], [16], [17].
The Precedence-Constrained Minimum-Cost Arborescence Problem with Waiting-Times (PCMCA-WT) is an N P-hard problem that was recently introduced by Chou et al. [15].The PCMCA-WT is a variation on the PCMCA problem characterized by the following differences.Given arc costs indicating the time required to traverse an arc, suppose there is a flow which starts at the root vertex r, that must reach every vertex in an arborescence T .For each precedence relationship (s, t) ∈ R, the flow must enter vertex s at the same time step, or before entering vertex t, which means that the flow can stop at any vertex and wait.The waiting time at vertex t is defined as the difference between the time at which the flow enters s and the time at which the flow reaches t.The objective of the problem is to find an arborescence T of minimum total cost, plus total waiting times, where the flow never enters s after entering t for all (s, t) ∈ R. Several MILP models for solving the problem were proposed in [15].
The PCMCA-WT problem can be formally defined as follows.Given a directed graph G = (V, A, R, r), where V = {1, . . ., n} is the set of vertices, A ⊆ V × V is the set of arcs, R ⊂ V × V is the set of precedence relationships, and r ∈ V is the root of the arborescence.Let c ij be a cost associated with each arc (i, j) ∈ A which represents the time required for the flow to travel from vertex i to vertex j.Let d j be the time step at which the flow enters vertex j ∈ V , and let w j be the waiting time at vertex j ∈ V .The objective of the problem is to find an arborescence T rooted at vertex r, that has a minimum total cost plus total waiting time, where the flow never enters t before entering s for all (s, t) ∈ R (i.e.d t ≥ d s for all (s, t) ∈ R). Figure 1 presents an example of an instance solved as a PCMCA-WT.The instance graph (left graph) shows the precedence relationship (1, 3) ∈ R marked as a dashed arrow, while the solution graph (right graph) shows an optimal solution of that instance, with the corresponding d i and w i value written next to each vertex.The graph on the right shows an optimal PCMCA-WT solution of cost 8 (sum of all the arcs cost plus total waiting time at each vertex), with a resulting waiting time of value 1 at vertex 3, since d 1 = 4, while d 3 = 3 and (1, 3) ∈ R.
The rest of this paper is organized as follows.Section II introduces the MILP model used in this study.Section III introduces a CP model that extends the MILP model introduced in Section II by introducing redundant constraints for a subset of the original inequalities and describing them in terms of CP constructs, in order to further exploit the capabilities of the CP solver.Section IV summarizes computational results, while some conclusions are outlined in Section V.

II. A MIXED INTEGER LINEAR PROGRAMMING MODEL
A polynomial-size MILP model for the PCMCA-WT was recently proposed by Dell'Amico et al. [18].The model extends a classical formulation for the MCA problem [19], through the addition of precedence-enforcing constraints.The precedence-enforcing constraints detect a precedence violating path by propagating a value along all the paths of the solution starting from t for all (s, t) ∈ R [16], [17].
A different version of the model that contains a smaller number of variables and constraints was also proposed in [18].The reduction is achieved by exploiting the special property of the PCMCA-WT, that is for any precedence relationship (s, t) ∈ R, the flow must enter vertex t at the same time step or after entering vertex s.This implies that it is possible to remove a precedence relationship (s, t) ∈ R when the input graph does not contain a zero-cost path that starts from t and ends in s.The reduced model for the PCMCA-WT proposed in [18] is summarized as follows.For further details the reader can refer to [18].
Let x ij be a variable associated with every arc (i, j) ∈ A such that x ij = 1 if (i, j) ∈ T , and 0 otherwise.Let y i be a variable associated with every vertex i ∈ V that indicates the order in which vertex i is visited on the path connecting vertex i to the root r.Let u t j be a variable associated with every vertex j ∈ V , and vertex t ∈ V where t is part of a precedence relationship (i.e.∃(s, t) ∈ R).Let d j be the time at which the flow enters vertex j ∈ V , and let w j be the waiting time before the flow enters vertex j.Let P ij ⊂ A be a simple directed path that starts from i and ends at j, and let c(P ij ) = (i,j)∈P c ij be the cost of that path.For each s ∈ V , let V s = {t ∈ V \{r} | ∃ (s, t) ∈ R, c(P ts ) = 0}.The PCMCA-WT can be formulated as the following MILP model.minimize s.t.: (i,j)∈A The set of constraints (2) impose that every vertex excluding the root must have exactly one parent.Constraints (3) are the subtour elimination constraints, which enforce that any feasible solution is acyclic.The set of constraints ( 2) and ( 3) guarantee that any feasible solution is an arborescence rooted at vertex r ∈ V .Constraints ( 4) and ( 5) fix the values of u t s and u t t to 0 and 1 respectively, for all (s, t) ∈ R, where t ∈ V s .Constraints (6) impose that if x ij = 1 then u t j ≥ u t i (see Figure 2 for further explanation).Constraint ( 7) sets the time step at which the flow enters the root to 0. Constraint (8) sets the waiting time at the root r to be equal to 0. Constraints (9) impose that when arc (i, j) ∈ A is selected to be part of the solution, then the flow enters vertex j at a time step that is greater than or equal to the time step at which the flow enters vertex i plus c ij .Constraints (10) enforce that the waiting time at vertex j is greater than or equal to the difference between the time at which the flow enters vertex j and the time at which the flow enters vertex i plus c ij .Constraints (11) enforce that the time at which the flow enters vertex t is greater than or equal to the time at which the flow enters vertex s for all (s, t) ∈ R. Finally, constraints ( 12)-( 15) define the domain of the variables, and M is an upper bound on the value of an optimal solution.

R
Fig. 2: An example on how a precedence relationship (s, t) ∈ R can be enforced by propagating the value of u t t along every path starting from t, and if the solution contains a path from t to s, then we are propagating a value of one to vertex s and imposing that u t s ≥ 1.However, we enforce u t s = 0, and therefore the solution violates the precedence relationship (s, t) ∈ R.

III. A NEW CONSTRAINT PROGRAMMING MODEL
The CP solver used in this study, CP-SAT [20] is a solver that utilizes integer programming techniques (linear relaxation, presolve, cuts, and branching heuristics) to enhance its performance [21] and has recently been shown to successfully deal with different combinatorial optimization problems [22], [23].Furthermore, the computational results in Section IV show that for the model considered in this work, the CP solver outperforms the MILP solver on a subset of the instances considered in terms of achieved average optimality gap, solution time, and the quality of the solutions obtained.Therefore, we introduce a CP model in this section that extend the MILP model introduced in Section II by adding the set of constraints ( 3) and ( 6) formulated as logical constraints, and merging the two sets of constraints ( 9) and ( 10) into one set of logical constraints.By doing so, we further exploit the capabilities of the CP solver.Since the CP solver used utilizes integer programming techniques, it is beneficial to include both the logical and linear form of the constraints in the model, so that when the logical constraint is not enforced by the SAT solver (i.e. the logical constraint is not included in the model by the solver), their equivalent linear constraint is included in the program when computing its linear relaxation.
Using a set of implication constraints which enforce the implied constraint when the value of the variable is true, the MILP model introduced in Section II can be extended using the following set of constraints.
Constraints ( 16) are the subtour elimination constraints modelling the nonlinear relationship y j = (y i + 1)x ij .Constraints ( 17) are the precedence-enforcing constraints modeling the nonlinear relationship u t j ≥ u t i x ij .Constraints (18) combine the two constraints ( 9) and ( 10) into a single equality constraint that model the nonlinear relationship Note that variables d i and w i are defined as integers (compared to the MILP model), since a CP solver only accepts integer variables and coefficients.This means that c ij for all (i, j) ∈ A should be integer or to be discretized before solving the model.The value of c ij can be discretized by multiplying every c ij by a constant k, and then considering only the integer part of the result.In order to compute the correct solution cost, the objective function value should be divided by k.A higher k value leads to higher numerical precision, whereas a low k value leads to a lower numerical precision and thus faster execution.Therefore, a k value which balances the two factors should be considered.In this study we only consider instances with integer coefficients.However, the interested reader can refer to [14] where the authors show how changing the k value can affect the computation time.

IV. EXPERIMENTAL RESULTS
The computational experiments are based on the benchmark instances of TSPLIB [24], SOPLIB [25], [26], and COMPILERS [27], originally proposed for the Sequential Ordering Problem (SOP) [28], [29], [30].The benchmark instances are the same instances previously adopted in [15], [18] for the PCMCA-WT with the following characteristics.The benchmark sets contain a total of 116 instances (81 open instances) ranging in size between 9 and 700 vertices, with an average of 248 vertices.Finally, all instances have integer coefficients (i.e. the weight of the arcs of the cost graph is integer).All the experiments are performed on an Intel Xeon Platinum 8375C processor with 8 cores running at 2.9 GHz with 16 GB of RAM.For all instances an upper bound on the value of the optimal solution (i.e.M ), is set to the value of the solution cost of solving the instance as a SOP, using a nearest neighbor algorithm [27].This is a valid upper bound for the cost of the optimal solution of the PCMCA-WT, being a feasible solution for the SOP a simple directed path that includes all the vertices of the graph, such that t never precede s for all (s, t) ∈ R.This implies that d t ≥ d s for all (s, t) ∈ R, with a waiting time equal to zero at each vertex by definition.
The computational results are generated using two solvers: a MILP Solver and a CP Solver.The MILP Solver is CPLEX v12.8 [31], and is run with 8 thread standard B&C algorithm, with the two parameters NodeSelect and MIP emphasis are set to BestBound and MIPEmphasisOptimality respectively.The CP Solver is Google OR-Tools [20] v9.5 CP-SAT solver, and is run with its default parameters with all 8 threads available are allocated for the solver.A time limit of 1 hour is set on the computation time of both solvers.For the rest of this section we will be referring to the MILP model introduced in Section II as BM (Basic Model), while the CP model introduced in Section III will be referred to as RM (Reinforced Model).
Tables I, II and III show the complete results of each model and solving method, where we report the following.For each instance, columns Name and Size report the name and size of the instance.Column ρ(R) reports the density of arcs in the set of precedence relationships computed as Column Best-Known reports the best-known bounds on the optimal solution for each instance as [LB, U B], where LB is the lower bound on the optimal solution, and UB is the best-known solution.The best-known solutions are obtained from the results appeared in [18], generated using the same computational setup and configuration used in this study.For each model solved with the corresponding solver, we report the following columns.Columns LB and UB report the lower/upper bound on the optimal solution achieved by the corresponding solving method of that model.Column Gap reports the optimality gap computed as U B−LB U B

. Column
Branches reports the number of branches created in the searchdecision tree, and is only reported when the models are solved with the CP Solver.Finally, column Time [s] reports the solution time in seconds and is only reported for the instances that are solved optimally within the time limit.In the tables, bold numbers indicate that a new best-known lower/upper bound is found.

A. Multi-threading Computation
The performance of CP solvers can often be greatly improved by the use of multi-threading computation, usually more than MILP solvers due to the different approaches used to solve the mathematical model.In this section, we assess the effect of multi-threading on the performance of the CP Solver and MILP Solver at solving the model introduced in Section II.The four instances ft53.1, prob.42,ESC78, and jpeg.4753.54were selected as both solvers are able to optimally solve those instances within the time limit using 8 threads.
Figure 3 reports the time required to optimally solve the different instances considered using a number of threads     A time of 60 minutes reported means that the respective solver was not able to optimally solve the instance within the time limit.
between 1 and 8.In the figure, a time of 60 minuets reported means that the respective solver was not able to optimally solve the instance within the time limit of one hour.
The results reported in Figure 3 show that the CP Solver substantially benefits from the use of multi-threading computation.Furthermore, the results show that the CP Solver is not able to optimally solve three out of four instances within the time limit when less than four threads are allocated for the solver.However, when allocating four or more threads, the CP Solver is able to optimally solve those instances.Furthermore, a drastic change in performance can be observed between four and five threads, reaching a speedup up to 93.5%.On the other hand, the MILP Solver does not seem to benefit as much from multi-threading for the instances considered, possibly due to the overhead of task distribution, and the waiting time incurred by the variety of methods run in parallel.Furthermore, we can notice less consistent gain when increasing the number of threads used by the MILP Solver compared to the CP Solver.It should be noted that the differences between the two solvers might be less extreme when more challenging instances are considered, but this is difficult to investigate as most instances are hard to solve optimally, even with longer computational time limit (hours) is allowed, and with eight threads allocated for the solvers.
In conclusion, the CP Solver appears to greatly benefit from multi-threading computation; therefore, all the experiments reported in this section were run on eight cores.

B. Analysis of the Results
In this section, we first compare and discuss the results achieved by the MILP and CP Solvers by solving the model BM.We then compare and discuss the results achieved by the CP Solver by solving the models BM and RM. to all the instances where both solvers find a feasible/optimal solution before reaching the time limit when solving the model BM.The Average solution time is computed on all the instances that are solved optimally by the both solvers.
The New best-known lower bounds and New best-known upper bounds rows report the number of instances where solving the model by each solver resulted in an improved lower or upper bound.Finally, New optimal solution row reports the number of instances where an optimal solution is found for an instance that was previously open, by each solver.
Considering the model BM, the MILP Solver achieves an average optimality gap of 0.340 across all the instances, but fails to solve a single instance (marked bold in the table) as it runs out of memory while solving the linear relaxation of the model.On the other hand, the CP Solver achieves an average optimality gap of 0.301 (a 11.5% improvement) when excluding the instance that is not solved by the MILP Solver, and an average optimality gap of 0.298 (a 12.4% improvement) across all the instances.By further inspecting the results, we notice that the CP Solver achieves a smaller average optimality gap within the time limit for instances with density less than 0.85 and size smaller than 400.
For a total of 27 instances that are optimally solved by both solvers, the MILP Solver has an average solution time of 297.6 seconds, while the CP Solver has an average solution time of 89.7 seconds (a 69.9% improvement).We should note that the CP Solver generally finds the optimal solution in less time compared to the MILP Solver on small to medium sized instances.
Finally, out of a total of 81 open instances, the CP Solver is able to find an improved lower bound for 9 instances (11.1%), an improved upper bound for 19 instances (23.5%), and finds the optimal solution of one instance that was previously open.On the other hand, the MILP Solver is not able to improve the best-known solution of any instance.Based on the experiments performed on the model BM presented in Section II, we can conclude that the CP Solver has an overall better performance at solving the given MILP model.The rest of this section discusses the results achieved by the CP Solver by solving the two models BM and RM.The results are summarized in Table V where we report the following.The average optimality gap reports the Average optimality gap of all the instances where the solver finds a feasible/optimal solution before reaching the time limit by solving both models.The Average solution time reports the average solution time in seconds of all the instances that are solved optimally by the solver when solving both models.The New best-known lower bounds and New best-known upper bounds rows report the number of instances where solving each model resulted in an improved lower/upper bound.Finally, New optimal solution report the number of instances where an optimal solution is found for an instance that was previously open.
In terms of achieved average optimality gap, the CP Solver achieves an average optimality gap of 0.287 (a 3.7% improvement) when solving the model RM, compared to solving the model BM.Furthermore, for a total of 36 instances that are solved optimally when solving both models, the CP Solver generates 57.9% less branches in the search-decision tree when solving the model RM compared to solving the model BM.By further inspecting the results, we notice that the CP Solver achieves a smaller average optimality gap within the time limit when solving the model RM for instances with density less than 0.89 and size less than 500, which means that the CP Solver performs better on a larger subset of the instances compared to solving the model BM.
For a total of 36 instances that are solved optimaly by the CP Solver when solving both models, the CP Solver has an average solution time of 146.9 seconds when solving the model BM, and an average solution time of 99.4 seconds (a 32.4% improvement) when solving the model RM.
Finally, out of a total of 81 open instances, when solving the model BM the CP Solver finds an improved lower bound for 9 instances (11.1%), an improved upper bound for 19 instances (23.5%), and finds the optimal solution for one instances that was previously open.On the other hand, when solving the model RM the CP Solver finds an improved lower bound for 4 instances (4.9%), an improved upper bound for 27 instances (33.3%), and finds the optimal solution for the same instance that was previously open.Based on the computational experiments and the improvements in the results achieved by the CP Solver when solving the model RM, we can conclude that duplicating the constraints can indeed improve the performance of the solver for the given model.

V. CONCLUSIONS
The computational experiments has shown that the CP Solver outperforms the MILP Solver at solving instances with sizes up to 500 with precedence relationships density that is less than 0.89.Furthermore, the CP Solver achieves a smaller average optimality gap and solution time compared to the MILP Solver.By adding constraint programming constructs to the MILP model, we were able to further exploit the capabilities of the CP Solver, and improve its performance at solving the instances.In terms of solution quality, and out of a total of 81 open instances, the CP Solver was able to find the optimal solution to an instance that was previously open, provide new best-known lower bounds for 13 instances, and establish new best-known solution for 46 instances.Based on the computational experiments performed, we have shown that the CP Solver performs better on average for the given models.Furthermore, duplicating constraints by defining them in their linear form and logical form further pushes the performance of the CP Solver.Future work will consider investigating new valid constraints/inequalities for the PCMCA-WT that can be used within a constraint programming paradigm to further utilize its potential.

Fig. 1 :
Fig. 1: Example of an instance solved as a PCMCA-WT.The graph on left shows the instance with its respective arc costs, and the precedence relationship (1, 3) ∈ R marked as a dashed arrow.The graph on the right shows an optimal PCMCA-WT solution of cost 8.

Fig. 3 :
Fig.3: Time required by the MILP Solver and CP Solver to optimally solve different instances with different number of threads.A time of 60 minutes reported means that the respective solver was not able to optimally solve the instance within the time limit.

TABLE I :
Computational results for TSPLIB instances.
MAURO DELL'AMICO ET AL.: MODELLING AND SOLVING THE PRECEDENCE-CONSTRAINED MINIMUM-COST ARBORESCENCE PROBLEM

TABLE II :
Computational results for SOPLIB instances.
Authorized licensed use limited to the terms of the applicable license agreement with IEEE.Restrictions apply.

TABLE III :
Computational results for COMPILERS instances.

TABLE IV :
Summary of the results achieved by solving the model BM with the MILP Solver and CP Solver.Table IV summarizes the results of solving the model BM by the MILP Solver and CP Solver, where we report the following.The Average optimality gap is computed with respect 428PROCEEDINGS OF THE FEDCSIS.WARSAW, POLAND, 2023Authorized licensed use limited to the terms of the applicable license agreement with IEEE.Restrictions apply.

TABLE V :
Summary of the results achieved by solving each model with the CP Solver.