Decision Optimization

Originally posted by: AndyHam

Can you please help me out to understand heightAtEnd? I would like to subtract the cum function value from the starting value as tasks are being completed. The following code failed to keep track of the cum function. Namely, the initial value of 100 does not change. I could not identify my error. Thanks!

using CP;
range M=1..2;
range J=1..5;
dvar interval T[j in J] size j*10;
dvar interval T2M[J][M] optional;
dvar sequence seqM[m in M] in all(j in J) T2M[j][m];

cumulFunction cumf[m in M]= step(0,100)  
            - sum(j in J) stepAtEnd (T2M[j][m], 0, 100); 
     
minimize max(j in J) endOf(T[j]);
constraints {
forall(j in J) alternative(T[j], all(m in M) T2M[j][m] );

forall(m in M, j in J) heightAtEnd(T2M[j][m],cumf[m]) == - sizeOf(T2M[j][m]);
  
forall(m in M) noOverlap(seqM[m]);   
}

Originally posted by: Petr Vilím

Hello Andy,

I think the problem is that capacity of cumul function cumulf[m] is not constrained, there are only heightAtEnd expressions.

Cumul functions are implicitly constrained to be >= 0. However it is done only for those cumul functions that have some another capacity constraint (heightAtEnd is not counted). So I think that the following constraint should fix the problem:

forall(m in M)
  cumf[m] >= 0;

Best regards, Petr

Originally posted by: AndyHam

Thanks! Once I have added the constraint, it worked.
However, I encountered another issue when I changed the heightAtEnd constraint as follow:
forall(m in M, j in J) heightAtEnd(T2M[j][m],cumf[m]) == - 1; //sizeOf(T2M[j][m]);
It should simply subtract by 1 from the initial value as the machine is completing tasks.
Instead, the model generates "conflict error".

using CP;
range M=1..2;
range J=1..5;
dvar interval T[j in J] size j*10;
dvar interval T2M[J][M] optional;
dvar sequence seqM[m in M] in all(j in J) T2M[j][m];

cumulFunction cumf[m in M]= step(0,100)  
            - sum(j in J) stepAtEnd (T2M[j][m], 0, 100); 
     
minimize max(j in J) endOf(T[j]);
constraints {
forall(j in J) alternative(T[j], all(m in M) T2M[j][m] );
forall(m in M, j in J) heightAtEnd(T2M[j][m],cumf[m]) == - 1; //sizeOf(T2M[j][m]);
forall(m in M) noOverlap(seqM[m]);   
forall(m in M) cumf[m] >= 0;
}

Originally posted by: Petr Vilím

Hello Andy,

I exported your model into a cpo file model.cpo. Then I run cpoptimizer executable, imported the model using "read model.cpo", then I computed the conflict by "tools conflict" and finally displayed the conflict by "display conflict". Here is the output:

alternative("T(5)", ["T2M(5)(1)", "T2M(5)(2)"], 1);
-1 == heightAtEnd("T2M(5)(1)", "cumf(1)");
-1 == heightAtEnd("T2M(5)(2)", "cumf(2)");

The problem is that constraint -1 == heightAtEnd("T2M(5)(1)", "cumf(1)") is forcing interval variable T2M(5)(1) to be present because heightAtEnd of an absent interval variable is 0. Similarly T2M(5)(2) must be present. However the constraint "alternative" says that only one of them should be present.

It was working before with sizeOf because sizeOf of an absent interval variable is zero (unless specified otherwise by additional parameter to sizeOf).

Similar to sizeOf, heightAtEnd accepts also additional parameter to specify the value when the interval variable is absent. So if you want to constraint the height only for the case the interval is present you can use that additional parameter:

-1 == heightAtEnd("T2M(5)(2)", "cumf(2)", -1);

This way, if the interval variable is absent then heightAtEnd returns -1 and therefore the constraint will be satisfied.

Best regards, Petr

Originally posted by: AndyHam

I have made many different attempts but all failed. The last thing I tried was to make another interval variable (move) that has the size of travel distance.
Then, I feed the heightAtEnd with sizeOf(move[j][t]).

presenceOf(itvJ2T[j][t])==presenceOf(move[j][t]);
startOf(itvJ2T[j][t])==startOf(move[j][t]);
sizeOf(move[j][t])==dNtoN[j.id][typeOfNext(seqTrk[t], itvJ2T[j][t], j.id, j.id)]; 
heightAtEnd(itvJ2T[j][t],cumBattery[t]) ==  -sizeOf(move[j][t]) ;

If you can shed some light, it will be great!.

Hi Petr,
Thanks for looking into the issue. I think I understand the cause.
Now, I am trying to port the lesson to my battery-charging VRP problem.
I still could not make the code work. I am enclosing the entire model and data files.
The problematic codes are showing below.
I sincerely appreciate your help.

cumulFunction cumBattery[t in Trucks]=  step(0, Q)      // initial Battery level
    - sum(j in aJobs) stepAtStart (itvJ2T[j][t], 0, Q)            //Battery consumption
    + sum(j in aJobs: j.t=="f") stepAtEnd (itvJ2T[j,t], 0, Q); //Battery recharging                     

forall(j in aJobs,t in Trucks) //Battery consumption
    // heightAtStart(itvJ2T[j][t],cumBattery[t],0) == - dNtoN[j.id][typeOfNext(seqTrk[t], itvJ2T[j][t], j.id, j.id)];
     heightAtStart(itvJ2T[j][t],cumBattery[t],-dNtoN[j.id][typeOfNext(seqTrk[t], itvJ2T[j][t], j.id, j.id)]) == - dNtoN[j.id][typeOfNext(seqTrk[t], itvJ2T[j][t], j.id, j.id)];
     heightAtStart(itvJ2T[j][t],cumBattery[t],0) == - presenceOf(itvJ2T[j][t])*dNtoN[j.id][typeOfNext(seqTrk[t], itvJ2T[j][t], j.id, j.id)];
 

forall(j in aJobs,t in Trucks: j.t=="f")  //Battery charging
      heightAtEnd(itvJ2T[j][t],cumBattery[t])  == g*sizeOf(itvJ2T[j][t]);

Thanks,
Andy

Originally posted by: Petr Vilím

Hello Andy,

I'm sorry, I don't understand what is the problem now. The attached model finds a solution easily. You don't like the solution for some reason? Is something wrong with it?

Petr

Originally posted by: AndyHam

Sorry to make you confused.
The cumBattery is not being calculated correctly. After the initial value of 78 is set, the function value should have been updated as the vehicle tour progresses.
I am enclosing the screen capture of the function value.

Originally posted by: Petr Vilím

Hello Andy,

I think the problem is that cumBattery is not constrained in other ways than by heightAtStatrt/End (similarly as cumf was not constrained at the beginning of our discussion). I think the following additional constraint should fix it:

forall(t in Trucks) {
  cumBattery[t] >= 0;
}

However then the problem doesn't have any solution. There are other too strong constraints?

Petr

Originally posted by: AndyHam

The code works well without the battery related codes, meaning the battery is causing issue. So I increased the battery capacity to a big number so that a vehicle does not need to visit a charge station. But the cum function still failed to track the battery consumption at any visit.

I tried all possible trials I can imagine, but all failed.

Thanks,

Andy

Originally posted by: Petr Vilím

Hello Andy,

I'm afraid I still don't fully understand. In the files eVRP.mod and eVRP.dat there is no upper limit of the capacity of cumBattery. Maybe the files are too old now? Could you upload new files?

You wrote: "The code works well without the battery related codes, meaning the battery is causing issue". What issue do you mean? The fact that there is no solution?

Note also that if an interval variable has zero length then heightAtStart and heightAtEnd are the same things (since start and end are happening at the same time).

Best regards, Petr

Originally posted by: Petr Vilím

I was playing with your example and I think that you may run into a bug: in some cases OPL may not evaluate a cumul expression (it may say that it is not fixed) or evaluate it wrongly (typically the height is 0). However the problem is only in OPL, solution in the engine is OK.

To work around this bug it is necessary to store the values in integer variables. The attached model shows the issue, the script at the end prints the values of integer variables and values of corresponding expressions.

Note that I increased the value Q in order to find first solution easily.

I hope it finally helps! Petr

Originally posted by: AndyHam

Please help me. I still cannot figure it out.

Originally posted by: Petr Vilím

Hello Andy,

could you be more specific please? What you cannot figure out? The heightAtStart/End expressions are OK now with the workaround? Or the problem is to find a solution with lower value of Q? Is it still about the eVRP.mod and eVRP.dat files you posted before?

Thanks, Petr

Originally posted by: AndyHam

Sorry! I was too much rush ^^

Q is the initial battery level (hard-constraint) so we should not increase. The solution found by the suggested code is not feasible. The minimum tour distance of 165km should force a vehicle to visit charging station two times. So, I put the original value of Q (77.75) and waited 10 mins, but the code could not find any feasible. Thanks!

 ! ----------------------------------------------------------------------------
 ! Search terminated by limit, no solution found.
 ! Best bound             : 65
 ! ----------------------------------------------------------------------------
 ! Number of branches     : 208,867,509
 ! Number of fails        : 104,431,950
 ! Total memory usage     : 19.8 MB (19.1 MB CP Optimizer + 0.7 MB Concert)
 ! Time spent in solve    : 600.02s (600.01s engine + 0.01s extraction)
 ! Search speed (br. / s) : 348,106.7
 ! ---------------------------------------------------------------------------
 

Schneider, M., Stenger, A., & Goeke, D. (2014). The electric vehicle-routing problem with time windows and recharging stations. Transportation Science, 48(4), 500-520.

Originally posted by: Petr Vilím

Hello Andy,

looking at the data, there are two jobs that could charge the battery (with j.t=="f"): S0 and S15. Both of them has size service time zero, i.e. j.s=0 and therefore size of those interval variables are zero:

dvar interval itvJob[j in aJobs] optional(j.t=="f") size j.s;

The constraint "alternative" propagates the zero size from itvJob to itvJ2T. The charge amount depends on the size of the interval variable:

chargeValues[j][t] == g*sizeOf(itvJ2T[j][t])

In our case sizeOf equals to 0 and therefore there is no charging. I think that's the main issue.

Best regards, Petr

Originally posted by: AndyHam

Thanks for locating the error!
I have addressed the error as following:
dvar interval itvJob[j in aJobs] optional(j.t=="f");
forall(j in aJobs: j.t=="c") sizeOf(itvJob[j])==j.s; 

Even with this change, the cum function still could not track the energy level as a truck travels on both my original code and your work-around.
BTW, your work-around code shows the cum function does not track the values correctly, but it does not provide the work-around of the issue.
I may be misunderstanding your code.
Can you please modify your work-around code to demonstrate that the battery level is successfully updated as a truck travels?
Thanks a lot!

Originally posted by: Petr Vilím

Hello Andy,

I applied your fix to the model file with the workaround. I attach the result.

It prints the following after solve:

OBJECTIVE: 165
From variables: 78  -13  From expressions: 0  0
From variables: 78  0  From expressions: 0  0
From variables: 0  -24  From expressions: 0  0
From variables: 0  0  From expressions: 0  0
From variables: 0  -5  From expressions: 0  0
From variables: 0  -38  From expressions: 0  0
From variables: 0  -36  From expressions: 0  0
From variables: 0  -11  From expressions: 0  0
From variables: 0  -38  From expressions: 0  0
From variables: 0  0  From expressions: 0  0
From variables: 0  0  From expressions: 0  0

It shows the bug: chargeValues and charngeExprs have the same value. Similarly dischargeValues and dischargeExprs.

Values of variables (e.g. chargeValues in our case) are given by the engine. However values of expressions (e.g. chargeExprs) are computed after solve by OPL from values of the variables. And there is a bug in this computation in case of cumulative functions. So what OPL displays as a value of cumulative function cannot be trusted (it applies also to cumBattery). It is necessary to store the values you need in integer variables and compute the cumulative function yourself.

Best regards, Petr 

Originally posted by: AndyHam

Petr,
I can store the values I need in integer variables. However, I was not able to compute the cumulative function by myself. That has been my desperate attempts for the past couple of weeks since I noticed the bug of cum function. If you can demonstrate the possibility, it will be very appreciated. 

This cum function for battery is very critical in many VRP problems related to EV and drone. I was planning to work on E-DARP after I complete E-VRP. Then, move to the drone scheduling problems. I will be looking forward to hearing about the work-around-sample or bug-fix release.
I sincerely appreciate all your help.
Thanks again,
Andy

Originally posted by: Petr Vilím

Hello Andy,

here it is. The script at the end computes the values of cumulative function for each truck and prints it.

Regards, Petr

Originally posted by: AndyHam

I think I confused you.
What I am looking for is to track and maintain the feasibility of tour as a truck visits customers based on the battery level. When the level goes down, the model is supposed to force the truck to visit a recharging station. Once the truck is being charged to an optimal level, the tour continues. The code you offered is the post processing. I need the cumul function-like work-around to maintain the feasibility during the solving. Thanks, Andy

Originally posted by: Petr Vilím

That's what the cumul function cumBattery does, no? The level of the battery never goes below zero..

Petr

Originally posted by: AndyHam

Here is the tour your code generates.
Depot=>S=>65=>98=>20=>24=>57=>S=>Depot (Where S refers the charging station).
This tour is infeasible because the cumulative distance up from S to customer 20 exceeds the maximum battery capacity of 78.

On the other hand, the following is a feasible tour generated by MIP.        
Depot=>65=>S=>98=>S=>20=>24=>57=>S=>Depot
The MIP forces the vehicle to visit charging stations through the tour to maintain the battery level.
Can you modify your work-around code to generate a feasible tour?

The enclosed figure shows the detailed tour.

Originally posted by: Petr Vilím

Hello Andy,

that's because there is only constraint cumBattery >= 0. There should be instead cumBattery <= Q (and then cumBattery >= 0 is added implicitly).

However, if I add this constraint then it is hard to find the initial solution. So I made two more changes in the model:

1) I constrained sizeOf of the last task (the one with j.t=="r" && j.id==2000+t) to be zero. I think that this additional constraint does not limit the solutions.

2) I changed the parameters. In particular I set FailureDirectedSearchEmphasis=1 and Workers=2. This setting allocates one worker to Failure-Directed Search right from the beginning and therefore Failure-Directed Search is used already in the first solution phase. This fixes the problem with initial solution. The second worker runs standard search.

Note that 2) will not work without 1) because we need to limit the search space. Without 2) the size of the last task could be 2^53 and the number of decisions for failure-directed search is too big to fit into memory.

I attach the modified model file.

Note that we may need to continue our discussion here: https://community.ibm.com/community/user/datascience/communities/community-home?communitykey=ab7de0fd-6f43-47a9-8261-33578a231bb7 because this forum will be read-only very soon.

Regards, petr