Learning Sensitivity of RCPSP by Analyzing the Search Process

1 Introduction

Scheduling problems are important for companies to execute their tasks on time and maximize their productivity. Once a schedule is computed, it is interesting to ask questions and to test different scenarios around this solution. For instance, how much time a company would save by increasing the capacity of a resource? Is it possible to complete a new order in time? What is the company's bottleneck? Which resource causes a slowdown in the production because? Answering these questions is important for continual improvement process. This corresponds to the sensitivity analysis of the solution. "The term sensitivity analysis refers to an analysis of the effect on the optimal solution to a linear programming problem of changes in the input-output coefficients, cost coefficients, and constant terms." [Dantzig, 1963]

The sensitivity analysis can be done by altering the model to encode a "what-if" scenario and solving the new problem. This process requires to solve a variant of the problem that was already solved which can be time consuming, especially if a distinct scenario is required for hundreds of resources or tasks. If the problem is a linear program, it is possible to use the dual solution to find a range of changes for the coefficients of the objective function or the constants on the right-hand side of the inequalities while keeping the optimality of the solution. In constraint programming, one can use the dual inference of an optimal solution to find the range of possible change of a parameter while keeping the optimality of the solution [Hooker, 1999]. The dual inference can be used to do the sensitivity analysis for any constraint satisfaction problem (CSP), integer programming, and mixed integer programming (MIP) [Hooker, 2006].

This article presents a method to predict the impact on the objective value of making changes to the solution or on constraints parameters. We test our method to predict whether increasing the capacity of a cumulative resource improves the solution and to predict the impact on the objective value of finishing a task earlier for a resource-constrained project scheduling problem (RCPSP). We do not just want to predict whether a change keeps the solution optimal, but how much it improves/deteriorates the solution. Particular care is given for making these predictions instantaneous, which prevent using a solver for conducting new optimization.

The method uses the information about the constraint CUMULATIVE [Aggoun and Beldiceanu, 1993] collected during the search and information about the resource usage in the solution found. It also uses information collected offline while solving similar instances obtained from historical data. Once a machine learning classifier is trained, it can quickly predict, for new instances of the problem, whether it is possible to obtain a better solution by increasing the capacity of the resource. This predictor does not need to run the solver a second time to evaluate "what if" scenarios. In fact, executing a trained classifier is instantaneous.

The rest of the paper is divided as follows. First, we present the background of this research. Second, we explain the methodology. Third, we present empirical results using a scheduling benchmark. Fourth, we conclude and present ideas for future work.

2 Background

We present background information about sensitivity analysis (what we want to do), lazy clause generation (with what we work), and the resource-constrained project scheduling problem (our application).

2.1 Sensitivity Analysis

The sensitivity analysis is the study of the impact on the solution of changing the value of a parameter. This problem is widely studied in the literature [Geoffrion and Nauss, 1977; Greenberg, 1998].

Hall et al. [Hall and Posner, 2004] identify 4 questions that sensitivity analysis attempts to answer. What is the range of possible changes to each parameter to keep the solution optimal, i.e. what is the shadow price? What is the new objective value for a specific change of a parameter? What is the new optimal solution following a change to a parameter? What are the answers to the three previous questions if we simultaneously apply several changes to the parameter values? This paper focuses on the second question. Hall et al. present algorithms for the sensitivity analysis for some type of singlemachine or identical parallel machines scheduling problems. They point out that it is more difficult to answer questions when the problem is NP-hard.

Hooker [Hooker, 1999] proposes an approach to do sensitivity analysis on linear or discrete problems. This method involves using the dual inference to obtain a proof of the optimality of the solution. Any optimization problem has a dual inference. Solving the dual inference of a problem consists in inferring the best possible bound on the optimal objective value from the constraints. With the proof found with the inference dual, it is possible to change the data of the problem as long as the proof remains valid. Hooker [Hooker, 1999] shows an example on a 0-1 linear programming problem. Dawande et al. [Dawande and Hooker, 2000] extend this idea to Mixed Integer linear programming problems. Dual inference has two limitations: it requires finding the optimal solution and exploring the entire search tree that proves that a better solution does not exist and it does not tell by how much the objective value changes if it does.

Hadzic et al. [Hadzic and Hooker, 2006] use binary decision diagrams (BDD) to enumerate solutions of a MIP. This can take time for big problems, but it allows to answer questions very quickly afterwards. It allows the user to answer questions about the right-hand side of an inequation but also about the variable domains. For instance, what is the best solution to the problem if a variable is set to a given value?

2.2 Lazy Clause Generation

Constraint solvers such as Chuffed [Chu, 2011] use the lazy clause generation [Stuckey, 2010] to take advantage of the high level modeling and understanding of the structure of the problem as well as the inference graph of SAT solvers. We recall how no-goods are generated since we collect information about these nogoods in order to perform sensitivity analysis.

In a solver using no-goods, propagators generate an explanation for each filtering. This explanation is represented as a clause which is a disjunction of literals. A literal is a Boolean variable which can be represented in the form [[X ≤ v]] or [[X = v]] or its negation where X is a variable and v is a value. The literal [[X ≤ v]] indicates that X is smaller than or equal to v. Consider the constraint X + Y ≤ 10 with dom(X) = [2, 10] and dom(Y ) = [2, 8]. The filtering algorithm prunes the domain of X to [2, 8] and generates the explanation ¬[[Y ≤ 1]] =⇒ [[X ≤ 8]] (equivalent to a clause [[Y ≤ 1]] ∨ [[X ≤ 8]]). When branching on variables and performing constraint propagation, the solver accumulates explanations in an implication graph. Upon failure, a cut in this graph generates a new clause, called no-good, that is added to the model. This new clause prevents the solver from performing a second time the branchings that led to the failure. In summary, a no-good is a clause computed from many filtering explanations generated by multiple constraints.

2.3 Resource-Constrained Project Scheduling Problem

The RCPSP has a set of non-preemptive tasks I and a set of cumulative resources R. A resource r ∈ R has a capacity cr and a task i ∈ I consumes hr,i units of resource r. A task i executes for pi units of time. One needs to find at what time si a task i can start without overloading the resources, i.e. ∑i|si≤t<si+pi hr,i ≤ cr for all time t and resource r. Tasks might be subject to precedence constraints: si + pi ≤ sj . We minimize the makespan. The model uses one CUMULATIVE constraint per resource [Baptiste and Pape, 2000].

3 Problem

We consider the problem of predicting the impact of changing a parameter of the problem on the objective value. Specifically, we answer these questions for a solution of the RCPSP:

1. By how much does the makespan change if one increases/decreases the capacity cr of a resource r?
2. By how much does the makespan change if one forces a task i to complete before a deadline d?

As such questions can be raised for multiple resources and multiple tasks, we consider it impractical to run a solver of on a "what if" scenario for each of these tasks/resources. We allow running the solver once to obtain the initial solution of an instance. We also allow to do off-line computations on historical data. But once the solution is obtained from the solver, we should be able to ask the questions above and instantaneously obtain answers. Since these questions are NPHard to answer, one cannot expect to obtain exact answers in a short computation time, but rather approximations.

Answering these questions explains to the decider why the solver returned such a solution. E.g.: Why is completing a task so late necessary? It also helps determine what should be changed in the parameters of the problem in order to obtain better solutions. E.g.: Should the company sell unnecessary machines/resources in order to purchase missing ones?

4 Methodology

We train a machine learning classifier/predictor to predict the fluctuation of the makespan. We collect information on a collection of historical instances H similar to the ones on which predictions are made (e.g. same resources but different tasks). For each historical instance, we solve "what if" scenarios by changing the capacity of a resource. For an instance H ∈ H, a resource r ∈ R, and a variation δ ∈ D for D = {0.5, 0.6, 0.7, 0.8, 0.9, 1, 1.1, 1.2, 1.3, 1.4, 1.5, 2}, we create a "what if" scenario SδH,r equivalent to the instance H but for the capacity of the resource r set to δ · cr.

We solve each scenario and record three types of features. The search features are collected during the search for a solution. The solution features are computed from the solution. The instance features are related to the instance such as the number of tasks and the capacity of the resources. All these features are presented in detail in the next sections.

We generate |H||R|(|D| 2 ) training examples for a classifier/predictor that predicts the variation of the makespan given the increase/decrease of a resource capacity. For each pair of scenarios Sδ1 H,r and Sδ2 H,r we use the search, solution, and instance features of scenario Sδ1 H,r . As an extra instance feature, we add the variation δ1 - δ2. The target value that needs to be predicted is given by the difference in the objective values (makespan) obtained for Sδ2 H,r and Sδ1 H,r . If one rather wants to train a classifier that predicts whether the objective value changes or not, we use the label 0 when the objective values for Sδ1 r and Sδ2 r are equal and 1 otherwise.

To train a model that predicts the impact of imposing a deadline to a task, we construct the scenarios as follows. For each historical instance H ∈ H, we create the original scenario Z∞ H,i . We solve this instance. For each task i that has no successor (the last task of a job), we randomly and uniformly select a deadline d between the earliest starting time of task i (taking into account the predecessors graph) and the starting time si in the solution of H. We add the constraint si+pi ≤ d to the instance H and call this scenario Z dH,i .

For each pair of scenarios Z∞ H,i and Z dH,i , we find a solution for both scenarios. We create a training example by using the search, solution, and instance features from scenario Z∞ H,i . We use, as an extra instance feature, the difference between the ending time of task i in the solution of S∞ i and the deadline d that is imposed to the task in scenario Sd i . The target value is the difference in the makespan if we train a predictor and a Boolean value if we only want to classify whether the makespan changes or not.

Creating the training data might take time as it requires to solve multiple instances, but once the classifier is trained, it can quickly predict whether increasing the capacity of a resource or imposing a deadline has an effect on the makespan.

4.1 Search Features

As the solver solves an instance, whether it is a scenario created from historical data or the actual instance on which one wants to perform sensitivity analysis, it needs to collect information about the search process. Since the RCPSP can easily be modeled using one CUMULATIVE constraint per resource, it is natural to collect statistics about the activity of this constraint. The intuition is that a CUMULATIVE constraint associated to a scarce resource is more often violated and causes more filtering during the search than an abundant resource.

We use two filtering algorithms for the CUMULATIVE constraint: the time-tabling [Beldiceanu and Carlsson, 2002] and the time-tabling-edge-finder [Vil´ım, 2011]. Other algorithms could also be used, but these two algorithms were already implemented in the solver Chuffed. For each resource r and for each filtering algorithm associated to its CUMULATIVE constraint, we define three search features:

1. the number of filtered values;
2. the number of times a variable domain gets filtered;
3. the average number of filtered values, i.e. the ratio of the two previous features.

It is important to link a no-good to the constraints that generated it. Without this link, information about the constraints is lost since no-goods can end up performing the majority of the filtering. We link a no-good to all the constraints that provided an explanation that contributed to its generation. To achieve this, we modified the solver to create a map between an explanation and the constraint that generated it. We link the no-good to the constraints by checking the explanations used to generate the no-good and the constraints that generated these explanations. If a no-good b is generated from an explanation of a no-good a, we bind the no-good b to the constraints that generated the no-good a. In the end, if a no-good induces a filtering or a backtrack, it is possible to give credits to the constraints that inferred that no-good. The activity of a no-good is the number of times that its explanations are used to produce another no-good, hence this search feature.

1. The sum of the activity of all no-goods linked to the CUMULATIVE constraint of resource r.

We collect information specific to the filtering of the objective value. For each resource r, we create an initialy empty list Lr. If the lower bound of the objective variable gets filtered to v, we analyze the implication graph that contains the explanations for this filtering. If one of these explanations is mapped to the CUMULATIVE constraint of the resource r or a no-good linked to this CUMULATIVE constraint, we insert the value v to the list Lr. Let M be the makespan of the solution returned by the solver at the end of the search. For each resource r, we have the three following search features.

1. The smallest value min(Lr) -M;
2. The largest value max(Lr) -M;
3. The average value avg(Lr) -M;

When training a model to predict how much the makespan changes when increasing/decreasing the capacity of a resource r, we only use the features related to the resource r.

When training a model to predict how much the makespan changes after imposing a deadline to a task i, we use the search features related to all resources. Moreover, we collect information about the filtering of the starting time variable si. These features are similar to features 1, 2, and 3 but specific to the variable si that is the focus of the prediction.

1. The number of times a value from dom(si) is filtered;
2. The number of times dom(si) is filtered;
3. The average number of filtered values from dom(si).

The next three features are inspired from the last three features. Rather than reporting the number of times a filtering occurs, we report how the variable ranks compared to other starting time variables. For instance, if si is the starting time variable that gets filtered the least often, the rank is 1. If it is the second least filtered variable, the rank is 2 and so on.

1. The rank with respect to how often si is filtered;
2. The rank with respect to how many values are filtered;
3. The rank with respect to the average number of filtered values when the variable gets filtered.

4.2 Solution Features

We design 5 features computed from the solution returned by the solver. The utilization rate U rate r of the resource r is the amount of energy pi · hr,i consumed by the tasks over the availability of the resource cr ·M, where M is the makespan.

The saturation duration U max r for a resource r is the time for which it is used at full capacity. Let Ur,t be the usage of r at time t, we have:

The maximum usage ratio U max ratio r of a resource r is the ratio between U max r and the makespanM.

The number of waiting tasks Wr counts the tasks waiting after resource r to start, i.e. the tasks that could have started earlier if r was available. To get this information, we must make sure that the task is not waiting for a predecessor to complete. Let predi be the set of predecessors of task i and ρi the time at which the last predecessor of i finishes.

The waiting time W duration r of resource r is the amount of time the tasks wait after resource r to start. Equation (5) checks if the task could be scheduled earlier, i.e. after its last predecessor but before the current starting time. If it is the case, we conclude the task does not wait after the resource r (it probably waits after another resource) and the waiting time of task i for resource r is set to W duration r,i = 0. Otherwise, we impute a waiting time of si - ρi to the resource r, i.e. the time between the last predecessor finishes and the starting time of task i. Finally, equation (6) computes the total waiting time and it is that quantity that is used as the solution feature.

Like for the search features, we only use solution features related to the resource r when we train a model to predict how much the makespan changes when increasing/decreasing the capacity of a resource r.

When training a model to predict how much the makespan changes after imposing a deadline to a task i, we use the solution features related to all resources in the problem. Moreover we add a solution feature indicating how many tasks finish before the task i. And a feature that is simply the makespan.

4.3 Instance Features

We use as instance features:

1. The number of tasks n;
2. The original capacity cr of each resource r;

When training the model to predict how much the makespan changes when increasing/decreasing the capacity of a resource, we use these features.

1. The new capacity of the resources;
2. The ratio of increase/decrease of the resource capacity.

When training the model to predict how much the makespan changes when imposing a deadline to a task, we add as an instance feature:

1. The difference between the ending time in the original scenario and the deadline in the second scenario.

5 Experiments

We use the RCPSP benchmarks PSPLib [Kolisch and Sprecher, 1997] and Pack [Carlier and N´eron, 2003]. PSPLib has instances of 30, 60, 90, and 120 tasks and 4 resources. The Pack benchmark has instances between 17 and 32 tasks and between 2 and 5 resources. There are more precedence constraints in PSPLib than Pack. We use the model provided by Minizinc [Stuckey et al., 2014]. We use the time-tabling and time-tabling-edge-finder as filtering rules for the CUMULATIVE constraint. We use a timeout of 10 minutes per instance for PSBLib and 3 hours for Pack since Pack has harder instances. Reaching the timeout can add noise to the dataset since the results become solver dependant, but this situation is likely to occur in a real context.

We randomly separate the instances of a benchmark into a training and a testing set with ratio 80/20. Since we produce many scenarios from one instance, we end up producing many training examples from training instances and many testing examples from testing instances. In no case is the testing data contaminated with the training data or vice versa. Solving theses instances require 80 cpu-days for PSPLib and 48 cpu-days for Pack. As the results show later, not that many instances are necessary to obtain good results. We apply a min-max normalization on all features to scale them between 0 and 1 using the relation x′i = xi-min(⃗x) max(⃗x)-min(⃗x) .

We use the random forest classifier and the random forest regressor from Scikit-Learn [Pedregosa et al., 2011]. We use the default parameters except for the number of trees (nEstimators) for which we set the value to 100.

We evaluate the quality of our predictions given which features are used and how many training examples are used. As it can take a lot of time to solve the "what if" scenarios, we want to know the number of examples necessary in the training set to obtain sufficiently good predictions. We compare the accuracy and the f1-score for the classification problem and the mean squared error for the regression problem.

We start by training on 100 randomly chosen examples in the training set and make the predictions on the test set. Then we add 100 more randomly selected examples from the training set and make the predictions on the test set. We do so until we consider all examples in the training set or up to 10,000 examples. We repeat this process 100 times with different training sets and test sets and take the average of the results. We also calculate a 95% confidence interval to ensure that our results do not depend on the training sets or test sets chosen.

5.1 Predicting the MakespanWhen Changing a Resource Capacity

Figure 1 shows the accuracy for the classification problem (i.e. predicting whether the makespan changes or not) according to the number of examples in the training set and the features used for the PSPLib benchmark. For the PSPLib benchmark, 54% of the examples have no change to the makespan (label 0). We compare four groups of features: instance features with solution features (green), instance features with search features (red), all features (blue), and the top 5 features (yellow). The classifier gets better accuracy when using all features, reaching 91.14% with 10,000 examples obtained from solving 892 scenarios. The accuracy of the classifier is similar if it uses the solution features with instance features (88.18%) as it does when using the search features with instance features (88.36%). When using all the features, we obtain an accuracy of 87.32% with 1,000 examples and 90% with 3000 examples. We obtain an f1-score of 90.18% with all features, 86.89% with solution+ instance features and 87.2% with search+instance features.

These results show that the classifier can properly predict whether changing the capacity of the resource affects the makespan. All types of features contribute to the quality of the prediction. We did a feature ranking with recursive feature elimination (RFE) to find the 5 most important features for the random forest classifier. The five most important features are the increase/decrease in the capacity of the resource, the smallest value min(Lr) -M, the utilization rate U rate r , the sum of the activity of all no-goods linked to the CUMULATIVE constraint, and the new capacity of the resource.

Figure 2 shows the accuracy and f1-score for the classification problem for the Pack benchmark. These instances are harder to solve and have 63% of examples with no change to the makespan (label 0). The classifier gets an accuracy of 90.69% with all features and 10,000 training examples which is slightly better than using only the solution+instance features (89.86%) or the search+instance features (85.41%). We obtain a f1-score of 87.41% when using all features, 86.37% with the solution+instance features, and 80.43% with the search+instance features. For Pack, the solution features seem to provide more information than the search features.

Figure 3 shows the mean squared error for the regression problem. For the PSPLib benchmark, with 10,000 examples in the training set, the regressor obtains a mean square error of 38.24 with all features, 41.9 with the solution+instance features and 87.67 with search+instance features. For the Pack benchmark, with 10,000 examples, the regressor obtains a mean square error of 74.58 with all features, 75.41 with solution+instance features and 126.13 with search+instance features. For both benchmarks, the solution features provide more useful information than the search features. Although search features do help obtaining better results. Interestingly, the regressor made no error in its prediction for 48.98 % of the testing examples for PSPLib and 64.47 % for Pack.

5.2 Predicting the MakespanWhen Imposing a Deadline

For the experiments of predicting the impact on makespan of imposing a new deadline to a task, we have fewer training examples than for predicting the makespan when modifying the capacity of a resource. Also, for the Pack benchmark, for a large majority of tasks, the makespan remains unchanged when we constrain the task to finish earlier. For that reason, the Pack benchmark is less interesting and we omit results for this benchmark. This is not the case for the PSPLib benchmark. 67.17% of examples are labeled with 0. Tasks often cannot finish earlier due to the precedence constraints. There is therefore often no impact on the makespan to finish a task earlier. To rebalance the training set, we duplicate the examples with a label of 1. The average change in makespan after imposing a deadline on a task is 0.69 and the maximum is 17. Even though there are many examples with a label of 0, this benchmark is more interesting for evaluating our method.

Figure 4 shows the mean squared error for the regression problem for the PSPLib benchmark. We replicated these experiments 250 times to better see the difference between the different types of features. With 3,250 examples in the training set, the regressor obtains a mean square error of 3.51 when using all features, 3.57 with solution+instance features, and 3.55 with search+instance features. There is not much difference between using a group of features than another. But using all feature provides the best predictions.

Figure 5 shows the distribution of the error when the regressor trains on 3,250 examples. The regressor makes no error or only small errors for a majority of its predictions.

Figure 6 shows the accuracy and the f1-score for the classification problem for the PSPLib benchmark. Predicting a change on the makespan when imposing a deadline on a task seems more difficult than when changing the capacity of the resource. Indeed, the CUMULATIVE constraint provides much information about the impact of the capacity of the resource. No constraint in the original model provides information about the deadline, since there is initially no deadline. Moreover, the smaller dataset and especially the lack of examples with change to the makespan (label 1) makes it harder to train. We notice the search features provide more information than the feature solutions. However, it is by training on all the features that we obtain the best classifying results.

6 Conclusion

We presented a method to predict how the makespan of a schedule is affected when the capacity of a resource is modified or a deadline is imposed on a task. Such predictions can help a decider to sell unnecessary resources to buy scarcer resources. It also tells whether a task could complete earlier. Our method takes advantage of the structure of a solution, but also some events that occur during the search for this solution.