Originally posted by: AndyHam

Dear IBM,

I am currently testing CP model for DARP TSP/VRP problem. My preliminary result is not promising when it compared to the known efficient MIP model.

For instance, the following two papers discuss rounded capacity inequalities and infeasible path elimination constraints for MIP. So, the MIP model proves an optimal of a medium size instances (50 jobs) within few seconds, while CP finds 97% to optimal after 300 seconds on my computer. Now, my question is how to improve the run-time of CP to be comparable with MIP in this DARP. Any advice will be appreciated.

Ropke, S., Cordeau, J. F., & Laporte, G. (2007). Models and branch‐and‐cut algorithms for pickup and delivery problems with time windows. Networks, 49(4), 258-272.

Braekers, K., Caris, A., & Janssens, G. K. (2014). Exact and meta-heuristic approach for a general heterogeneous dial-a-ride problem with multiple depots. Transportation Research Part B: Methodological, 67, 166-186.

Thanks,

Andy

Originally posted by: PhilippeLaborie

Hi Andy,

The principal way to improve the run time of CP Optimizer is ... to work on the model itself. Could you describe your CP Optimizer model?

Philippe

Originally posted by: AndyHam

Hi Philippe,

Happy new year!

It is a very similar model with the one shared in this forum ("VRP using CP"). The pickup quantity is the size of boarding passengers and the drop quantity is the size of unboarding passengers. There are few other constraints we can add to DARP, for instance, the total acceptable travel hour of passenger. However, I am testing a very basic DARP CP model (attached) against MIP.

Also if you look in to the code, the objective is hierarchical, minimizing number of vehicles used and then minimizing total distance.That seems to make CP slow down a bit. I tried staticLex, but not much improvement. The best known solution is 41 622.8 for the test instance (4 vehicles; 1622.8 travel distance).

Background:
While I am conducting a literature review, I was surprised by MIP's run time by adding valid inequalities constraints. So, I am looking for a corresponding method in CP. Also, I am seeing there are almost no work with CP in DARP. My topic is DARP real-time scheduling. I would like to first setup a baseline demonstrating CP is strong in DARP. Then, extend the model for real-time scheduling. In a given day, a taxi company would have a set of pre-orders given for the day. Then, during the day, new orders dynamically arrive at the system. How to adjust the original schedule based on the new orders would be the key discussion point. Let me seek your advice when I get there.

Regards,
Andy

Originally posted by: PhilippeLaborie

As the dominant aspect of the problem is routing and total transition costs, I don't think you should expect too much of CP Optimizer on optimality proofs for this type of problem.

Things could change if you add additional constraints, like for instance maximum allowed duration between pickup and delivery for the trips, ...

This being said, I played a little bit with your model and in particular added some intervals representing the trip that contribute with a "pulse" on the vehicle load cumul function. This helps a bit, in general, this helps as it explicitly correlates the increase and decrease steps on the cumul function. I also noticed that in this problem (at least on this instance), it helps to constrain the number of vehicles right from the beginning. For instance CPO easily finds a solution with 3 vehicles on this instance if you explicitly constrain the number of vehicles to be less than 3. Here is the corresponding model, it finds a solution with aggregated cost around 31 711.

Originally posted by: AndyHam

It is awesome. The model is much faster. The solution 31 711 is very close to the known optimal value of 31 600 (I just realized this optimal value from the following article. The one I wrote earlier was from the earlier article). The test instance is DD50 from the following article:

Ropke, S., & Cordeau, J. F. (2009). Branch and cut and price for the pickup and delivery problem with time windows. Transportation Science, 43(3), 267-286.

Now, I am trying to meditate your work as usual.
(1) before(seqVeh[v],itvJ2V[j][v],itvJ2V[jh][v]); 

This is valid inequality, right? 
Let me study DARP MIP model to see if we can add more valid inequalities into CP.
Also, let me do more experimentation to see if how much this constraint helps.

(2) additional intervals representing the trip that contribute with a "pulse" on the vehicle load cumul function

I am a little confused about this change, although it improved the performance a lot ^^
Instead of adding new, we could have modified the original as follows:
     cumulFunction Loading[v in Vehicles] = sum(j in Jobs: j.de>0) pulse(itvJ2V[j][v], j.de);

However, it generated the obj value of 31 338, which is lower than the proven optimal value of 31 600, indicating this change is not valid.
I will validate the model with smaller instances.

Regards,
Andy

Originally posted by: PhilippeLaborie

(1) Yes. As both pickup and delivery use the same vehicle, if a vehicle v is used, the pickup occurs before the delivery on that vehicle. It is a valid inequality if you like, in CP, we rather speak of redundant or surrogate constraint (as we do not only handle inequality constraints). 

(2) No, you cannot have a pulse only on the pickup intervals (or delivery). The pulse start to increase the vehicle load at the start time of the pickup and stops at the end time of the delivery. The vehicle load is affected everywhere between these two time-points. That's why you need an additional interval. And indeed if you only use a pulse on the pickup intervals it results in a less constrained (and wrong) problem, which explains why you get "solutions" better than the optimal one.

Note that my model has several differences compared to yours. I must admit I'm not sure which ones are really important. Some of them are probably completely useless and as I quickly worked only with one instance, the model may very well be over-tuned for this particular instance ...

Originally posted by: AndyHam

I think I got a clear picture ^^
I found more valid inequalities in your code:

   endBeforeStart(itvJob[j],itvJob[jh], Dist[j.id][jh.id]);
   sizeOf(useCap[j]) >= Dist[j.id][jh.id] + j.st + jh.st;
 

Regarding the pulse, the additional interval is basically merging two separate intervals (one for pickup and another for delivery) into one. So, we can use "Pulse" function. Cool ^^
I now understand your advice, "this helps as it explicitly correlates the increase and decrease steps on the cumul function.", compared to  my original stepAtStart 

Originally posted by: AndyHam

I would like to see if we can apply the known TSP valid inequalities. For that, we need a binary variable X(i,j)(k) if vehicle k travels from node i to node j. Here is my question. Is there any way to drive X(i,j)(k) from the CP variables during the run-time?

I tried out the following code, but the model says there is no solution after adding the code. I am trying to retrieve the next job ID for each job from the sequence.

forall(j in Jobs, v in Vehicles:  j.id<= nj*2-1)
   nextJob[j]== typeOfNext(seqVeh[v], itvJ2V[j][v], j.id, j.id); 

Originally posted by: PhilippeLaborie

The model is infeasible because yo say that the next of job j is the same on all vehicles v. As only one vehicle will be selected there will be only one v such that "typeOfNext(seqVeh[v], itvJ2V[j][v], j.id, j.id)" is really an actual job (other than j, so not j.id). In all the other cases, the value of the expression will be the "absence" value j.id. But nextJob[j] cannot be at the same time equal to j.id and different from j.id. The model is infeasible.

Concerning your more general question about the valid inequalities, I don't think they will help. As you say, they require to make all these binary variables X(i,j)(k) explicit whereas the main advantage of a CP model is not to handle them explicitly. And if the valid inequalities you are thinking about are the ones one can directly deduce from the time-windows, then, these ones are automatically discovered by the noOverlap constraint (for that you can slightly increase the inference level of the noOverlap: NoOverlapInferenceLevel="Medium").

Originally posted by: AndyHam

Thanks for the explanation. In particular, "And if the valid inequalities you are thinking about are the ones one can directly deduce from the time-windows, then, these ones are automatically discovered by the noOverlap constraint". Your suggestion of increasing the inference level of the noOverlap works well with this model. This option helped to find a better solution quickly.

Can I ask another question? This is about the drone scheduling. I am trying to calculate the total drone flight time of two consecutive jobs: Drop and Pkup.

Originally posted by: PhilippeLaborie

On your drone question: 

1/ Syntax

In OPL you cannot use an integer/float decision expression or constraints in an if (ct1) { ct2 } statement. Instead, you should use the implication constraint (ct1 => ct2) 

2/ Formulation

This being said, in CP it is always good (though I know it is not always possible): (1) to factorize the different elements of the model as much as possible, (2) to avoid this type of composite constraints (*) and (3) to avoid posting a quadratic number of constraints. Here your formulation cumulates the 3 problems. Depending on the context there is probably better ways to formulate the constraint. I see that you are handling pickup and drop activities. I suppose there is some precedences between these activities? Do you have 1 pickup activity for every drop activity ?

* Except for constraints like presenceOf(x)=>presenceOf(y) that are optimized and strongly exploited by the engine.

Originally posted by: AndyHam

"factorize the different elements of the model as much as possible" is something I need to keep it mind.
 

Based on your advice,  I found a temporary solution by using Big-M technique, instead of "IF" condition.

forall(v in Vehicles, j,jh in Jobs: j.type=="Drop"  && jh.type=="Pkup"  && v<=nDrone)

    Dist[j.id][typeOfNext(seqVeh[v], itvJ2V[j][v], j.id, j.id)] + sizeOf(itvJ2V[j][v]) + sizeOf(itvJ2V[jh][v])   

Originally posted by: PhilippeLaborie

The model with big-M will not propagate more than the implication I suggested I think, it will be more or less the same.

But I would at least simplify the "typeOfNext(seqVeh[v], itvJ2V[j][v], j.id, j.id)" that occurs twice.

You want:

( typeOfNext(seqVeh[v], itvJ2V[j][v], j.id, j.id)==jh.id ) => ( Dist[j.id][typeOfNext(seqVeh[v], itvJ2V[j][v], j.id, j.id)] + sizeOf(itvJ2V[j][v]) + sizeOf(itvJ2V[jh][v])   <= DroneEndurance  );

So why not:

( typeOfNext(seqVeh[v], itvJ2V[j][v], j.id, j.id)==jh.id ) => ( Dist[j.id][jh.id] + sizeOf(itvJ2V[j][v]) + sizeOf(itvJ2V[jh][v])   <= DroneEndurance  );

It is always good to replace a variable by a constant !

Originally posted by: AndyHam

( typeOfNext(seqVeh[v], itvJ2V[j][v], j.id, j.id)==jh.id ) => ( Dist[j.id][jh.id] + sizeOf(itvJ2V[j][v]) + sizeOf(itvJ2V[jh][v])   <= DroneEndurance  );

It works like a magic!
Another great lesson for me :)