Topic
• 3 replies
• Latest Post - ‏2013-05-17T12:44:52Z by Petr Vilím
JorisK
5 Posts

# Pinned topic using constraint programming to determine tight variable domains

‏2013-05-06T12:51:55Z | constraint domain propagation queries

Dear,

I'm currently trying to solve a problem involving a complex mixture of a routing and a scheduling problem.

Some (but not all) features of my problem are:

-heterogeneous fleet of vehicles

-split deliveries (i.e. customer is visited by multiple vehicles, a single vehicle may visit the same customer more than once)

-time windows, time lags and sequence dependent setup times

One of our solution approaches utilizes a Mixed Integer Programming (MIP) model. I would like to strengthen the model by tightening bounds and adding a number of valid inequalities. For example:

-The number of required deliveries to a single customer are unknown and depend on the capacity of the vehicles. A weak upper bound on the number of required deliveries is the demand of the customer divided by the capacity of the smallest vehicle.

-A time window is associated with each customer. Currently, each individual delivery for a customer inherits this time window. It would however be better to tighten the time windows associated with individual deliveries. E.g. setting the time windows of deliveries d1, d2 for customer c1 (time window [0,10]) to d1=[0,5], d2=[5,10] is better than setting d1=d2=[0,10].

-Certain customers are incompatible with each other, i.e. it is not possible to satisfy the demand of all customers due to capacity and time constraints. For two incompatible customers c1, c2, this could lead to a cut like: c1+c2 <= 1

Reasoning about the variable domains in the context of my problem is notoriously difficult. Hence I was wondering whether I could use constraint programming to do the job for me. Would it be possible to do the following:

a. model my problem using Ilog constraint programming (java interface)

b. apply constraint propagation (but NO branching)

c. query the domains of the variables after all constraints have been propagated and adjust the variable domains of my MIP model accordingly

1. Are steps b, c possible using IloCP+Java? Is there a way to query the variable domains after constraint propagation?

2. Does this make sense, or is this a rather ridiculous idea (I'm for example not sure whether or not Cplex already implements domain reductions or whether cplex would benefit from the strengthened domains)? As I'm not an expert in constraint programming, it would take me a considerable amount of time to implement my problem, so basically I'm trying to determine whether it is worth the effort :)

• Petr Vilím
29 Posts

#### Re: using constraint programming to determine tight variable domains

‏2013-05-06T20:46:59Z

Hello,

what you suggest makes perfect sense and the API of CP Optimizer supports it. And possible to do even more, in particular you can inject solution found by CPLEX a starting point for CP Optimizer search and try to improve it.

You can get the maximum from both solvers by modelling the problem in completely different ways. In CP Optimizer you can use interval variables for actions/activities, noOverlap constraint with setup times for activities that cannot overlap. There are also optional interval variables for activities that may but doesn't have to be executed what allows to model alternatives in the schedule (using constraint 'alternative'). A good starting point to read about modeling scheduling problems using CP Optimizer may be the following in the documentation:

CP Optimizer > CP Optimizer C++ API Reference Manual > Concepts (on the top of the right page) > Introduction to Scheduling Concepts in CP Optimizer

Those mathematical modeling concepts are the same for C++ and Java, so don't mind that it is inside C++ documentation.

That would be for the part (a). For (b), to apply propagation but do not search, just call propagate() on IloCP. It returns true if the model was not found infeasible during the propagation. After propagate, to answer part (c), you can call functions on IloCP to check variable domains, for example getStartMin(intervalVariable), getMax(integerVariable), isPresent(intervalVariable) etc.

Finally, give a try to CP Optimizer search. I understand that your goal is to improve your MIP model, but I believe that CP Optimizer may find some pretty good solutions for this kind of problem too. Or it can improve a solution found by CPLEX if the solution is passed to IloCP using setStartingPoint function.

Best, Petr

• JorisK
5 Posts

#### Re: using constraint programming to determine tight variable domains

‏2013-05-16T10:02:21Z

Hello,

what you suggest makes perfect sense and the API of CP Optimizer supports it. And possible to do even more, in particular you can inject solution found by CPLEX a starting point for CP Optimizer search and try to improve it.

You can get the maximum from both solvers by modelling the problem in completely different ways. In CP Optimizer you can use interval variables for actions/activities, noOverlap constraint with setup times for activities that cannot overlap. There are also optional interval variables for activities that may but doesn't have to be executed what allows to model alternatives in the schedule (using constraint 'alternative'). A good starting point to read about modeling scheduling problems using CP Optimizer may be the following in the documentation:

CP Optimizer > CP Optimizer C++ API Reference Manual > Concepts (on the top of the right page) > Introduction to Scheduling Concepts in CP Optimizer

Those mathematical modeling concepts are the same for C++ and Java, so don't mind that it is inside C++ documentation.

That would be for the part (a). For (b), to apply propagation but do not search, just call propagate() on IloCP. It returns true if the model was not found infeasible during the propagation. After propagate, to answer part (c), you can call functions on IloCP to check variable domains, for example getStartMin(intervalVariable), getMax(integerVariable), isPresent(intervalVariable) etc.

Finally, give a try to CP Optimizer search. I understand that your goal is to improve your MIP model, but I believe that CP Optimizer may find some pretty good solutions for this kind of problem too. Or it can improve a solution found by CPLEX if the solution is passed to IloCP using setStartingPoint function.

Best, Petr

Dear,

Thank you for your reply! I have implemented my model into cp optimizer (v12.5.1) via the java interface. Unfortunately, the results are poor (only very small instances can be solved), basically due to weak domain propagation.
Here I'll post some of the propagation issues I encounter using a highly simplified version of my model (java implementation attached). Hopefully you can tell me what I'm doing wrong.

Note: in my experiments, I have set: cp.setParameter(IloCP.IntParam.DefaultInferenceLevel, IloCP.ParameterValues.Extended);

First, a simplified version of my model:

================================
tuple Customer{
key int id; //unique identifier 1..n
int demand; //demand of customer
int release; //Time window during which customer can be serviced
int max_deliveries;
}

tuple Delivery{
Customer customer;
int id;
int duration; //Time required to perform delivery
int capacity; //Amount delivered
}

//Not all customers can be served, hence optional
dvar interval customers[c in Customers] in c.release ... c.deadline optional
forall(d in Deliveries)
dvar interval deliveries[d] in d.customer.release ... d.customer.deadline size d.duration optional

//The deliveries must be made in order, i.e first delivery 1, then delivery 2 etc:
//CONSTRAINTS 1a,b:
forall(d1,d2 in Deliveries: d1.id+1==d2.id){
endBeforeStart(delivery[d1],delivery[d2]);
presenceOf(delivery[d2])=>presenceOf(delivery[d1]);
}

//The customer is 'satisfied' iff sufficient deliveries are made such that its demand is fulfilled:
//Constraints 2:
forall(c in Customers){
precenceOf(customer[c]) ==
(sum(d in Deliveries: d.customer==c) (precenceOf(delivery[d])*d.capacity) >= c.demand);
}

//The objective simply maximizes the number of satisfied customers:
maximize sum(c in Customers)precenceOf(c)

================================

1. Instance 1
One customer with demand 10, time window [100,200]
Two deliveries with capacity 5, duration 10

Expected propagation:
a. Propagation on time windows of deliveries due to constraints 1a,b. Delivery d1 must start within the interval [100,180]. Delivery d2 must start within [110,190]

Result after invoking cp.propagate():
c1[0..1: 100..200 -- (0..100)0..100 --> 100..200]
d1[0..1: 100..190 -- (10)10 --> 110..200]
d2[0..1: 110..190 -- (10)10 --> 120..200]

Conclusion:
The starts of the intervals d1 and d2 have not been propagated well. To have a non-zero objective, interval c1 must be present. To have c1 present, both d1 and d2 must be present (constraint 2). Having d1 and d2 present implies (constraints 1a,1b):
d1[0..1: 100..180 -- (10)10 --> 110..190]
d2[0..1: 110..190 -- (10)10 --> 120..200]

Any suggestions how to tighten the propagation in this instance?
In my real problem, I have a heterogeneous fleet; hence the number of deliveries required to satisfy a customer is unknown and depend on the vehicles that perform the deliveries. However a simple lower bound is: Demand of customer divided by the largest vehicle, rounded up.

2. Instance 2
One customer with demand 10, time window [100,200]
Two deliveries with capacity 5, duration 60

Expected propagation:
a. To satisfy the customer, 2 deliveries must be made. However, since each delivery takes 60 time units, it is impossible to perform 2 deliveries within the time window of the customer. Hence, this customer cannot be satisfied.

Result after invoking cp.propagate():
c1[0]0 -- (0)0 --> 0]
d1[0..1: 100..140 -- (60)60 --> 160..200]
d2[0]0 -- (0)0 --> 0]

Conclusion:
Excellent propagation. CP correctly identified that customer c1 cannot be satisfied.

Question: so in this case, CP is capable of identifying the fact that customer c1 could not be scheduled. Imagine that I have 2 customers which can be scheduled individually, but not together (due to resource restrictions). These 'sets' of infeasible customers are highly valuable to me as I can use these sets to generate cuts in my Integer Programming model. Is there a smart way to produce these sets using CP? One approach would be to do the following:

for every combination of customers C:

forall(c in C)

c.setPresent()

forall(c not in C)

c.setOptional();

Next you can propagate the model again and see whether the result is infeasible. This approach seems however rather expensive, especially since you would have to run this approach for every combination of customers? Is there a smarter approach?

br,

Joris

#### Attachments

Updated on 2013-05-16T13:04:38Z at 2013-05-16T13:04:38Z by JorisK
• Petr Vilím
29 Posts

#### Re: using constraint programming to determine tight variable domains

‏2013-05-17T12:44:52Z
• JorisK
• ‏2013-05-16T10:02:21Z

Dear,

Thank you for your reply! I have implemented my model into cp optimizer (v12.5.1) via the java interface. Unfortunately, the results are poor (only very small instances can be solved), basically due to weak domain propagation.
Here I'll post some of the propagation issues I encounter using a highly simplified version of my model (java implementation attached). Hopefully you can tell me what I'm doing wrong.

Note: in my experiments, I have set: cp.setParameter(IloCP.IntParam.DefaultInferenceLevel, IloCP.ParameterValues.Extended);

First, a simplified version of my model:

================================
tuple Customer{
key int id; //unique identifier 1..n
int demand; //demand of customer
int release; //Time window during which customer can be serviced
int max_deliveries;
}

tuple Delivery{
Customer customer;
int id;
int duration; //Time required to perform delivery
int capacity; //Amount delivered
}

//Not all customers can be served, hence optional
dvar interval customers[c in Customers] in c.release ... c.deadline optional
forall(d in Deliveries)
dvar interval deliveries[d] in d.customer.release ... d.customer.deadline size d.duration optional

//The deliveries must be made in order, i.e first delivery 1, then delivery 2 etc:
//CONSTRAINTS 1a,b:
forall(d1,d2 in Deliveries: d1.id+1==d2.id){
endBeforeStart(delivery[d1],delivery[d2]);
presenceOf(delivery[d2])=>presenceOf(delivery[d1]);
}

//The customer is 'satisfied' iff sufficient deliveries are made such that its demand is fulfilled:
//Constraints 2:
forall(c in Customers){
precenceOf(customer[c]) ==
(sum(d in Deliveries: d.customer==c) (precenceOf(delivery[d])*d.capacity) >= c.demand);
}

//The objective simply maximizes the number of satisfied customers:
maximize sum(c in Customers)precenceOf(c)

================================

1. Instance 1
One customer with demand 10, time window [100,200]
Two deliveries with capacity 5, duration 10

Expected propagation:
a. Propagation on time windows of deliveries due to constraints 1a,b. Delivery d1 must start within the interval [100,180]. Delivery d2 must start within [110,190]

Result after invoking cp.propagate():
c1[0..1: 100..200 -- (0..100)0..100 --> 100..200]
d1[0..1: 100..190 -- (10)10 --> 110..200]
d2[0..1: 110..190 -- (10)10 --> 120..200]

Conclusion:
The starts of the intervals d1 and d2 have not been propagated well. To have a non-zero objective, interval c1 must be present. To have c1 present, both d1 and d2 must be present (constraint 2). Having d1 and d2 present implies (constraints 1a,1b):
d1[0..1: 100..180 -- (10)10 --> 110..190]
d2[0..1: 110..190 -- (10)10 --> 120..200]

Any suggestions how to tighten the propagation in this instance?
In my real problem, I have a heterogeneous fleet; hence the number of deliveries required to satisfy a customer is unknown and depend on the vehicles that perform the deliveries. However a simple lower bound is: Demand of customer divided by the largest vehicle, rounded up.

2. Instance 2
One customer with demand 10, time window [100,200]
Two deliveries with capacity 5, duration 60

Expected propagation:
a. To satisfy the customer, 2 deliveries must be made. However, since each delivery takes 60 time units, it is impossible to perform 2 deliveries within the time window of the customer. Hence, this customer cannot be satisfied.

Result after invoking cp.propagate():
c1[0]0 -- (0)0 --> 0]
d1[0..1: 100..140 -- (60)60 --> 160..200]
d2[0]0 -- (0)0 --> 0]

Conclusion:
Excellent propagation. CP correctly identified that customer c1 cannot be satisfied.

Question: so in this case, CP is capable of identifying the fact that customer c1 could not be scheduled. Imagine that I have 2 customers which can be scheduled individually, but not together (due to resource restrictions). These 'sets' of infeasible customers are highly valuable to me as I can use these sets to generate cuts in my Integer Programming model. Is there a smart way to produce these sets using CP? One approach would be to do the following:

for every combination of customers C:

forall(c in C)

c.setPresent()

forall(c not in C)

c.setOptional();

Next you can propagate the model again and see whether the result is infeasible. This approach seems however rather expensive, especially since you would have to run this approach for every combination of customers? Is there a smarter approach?

br,

Joris

Hello Joris,

I'll start with a detailed explanation of propagation in CP Optimizer (CPO for short). CPO gathers all constraints such as "presenceOf(x) => presenceOf(y)" into "logical network". Logical network contains all kind of implications between presence statuses, in particular negations are also allowed. Logical network is basically a 2-SAT and it is also propagated this way.

Precedence constraints fooBeforeFoo (for example endBeforeStart) are gathered in a similar way into "temporal network". Propagation of temporal network is using logical network to know what intervals variables influence each other. You can find more details about this propagation in the paper by P. Laborie and J.Rogerie: Reasoning with Conditional Time-intervals. This is the propagation that is able to see that customer cannot be satistfied in the instance 2. The paper can help you to better understand what's going on behind the scene.

The important fact is that for a constraint to be included in logical network, it must be written in the form: [!]presenceOf(x)=>[!]presenceOf(y). So the constraint to satisfy customers: precenceOf(customer[c]) ==  (sum(d in Deliveries: d.customer==c) (precenceOf(delivery[d])*d.capacity) >= c.demand) is not included in logical network. And currently there is no presolve step that would derive some redundant implications form this constraint.

So, if you think that there may be some redundant implications derived from this constraint, then it may help to add those implications into the model manually. The advantage is that temporal network will be able to take those implications into account and therefore propagate better.

To the instance 1:

As you say to have a non-zero objective, interval c1 must be present. However there is no constraint that forces objective to be non-zero. The purpose of objective is to propagate that the objective value must be better than objective value of the best solution found so far. As only propagate() and not solve() is called, there is no previous solution and therefore objective is completely ignored.

So, to get better propagation, it would be necessary to add a constraint saying that objective expression is > 0. The more you constraint the objective value the more propagation you may get. So it may help to wait until some initial solution is found. And also it may pay off to redo the propagation when a better solution was found.

To the instance 2:

Your idea to set some intervals as present and propagate is perfect. However I'm not sure whether you have to do it for every combination of customers C. Lets say that we try to set only one optional interval X to be present and call propagate(). If it fails then X must be absent. If it doesn't fail, then some other optional interval Y may become absent and yet another interval Z may become present. Then it is possible to add the following implications both into MIP and CP models:

presenceOf(X) => !presenceOf(Y)

presenceOf(X) => presenceOf(Z)

Of course, it will capture only binary disjunctions. It is possible to go further, but it will be more costly: set two interval variables X and Y to be present at the same time. And depending on the outcome add into the model:

if propagation fails: !presenceOf(X) ||  !presenceOf(Y) (will be included in logical network)

if Z becomes absent: (presenceOf(X) && presenceOf(Y)) => !presenceOf(Z)

if Z becomes present: (presenceOf(X) && presenceOf(Y)) => presenceOf(Z)

Best, Petr