CPLEX Studio - Robust Optimization

Question

I'm currently working on a robust optimization of the VRP with CPLEX Studio and at the moment I'm not getting anywhere.

Setting up the deterministic problem is no problem so far, but I am failing at implementing the uncertainty. I am working with uncertainty sets and would like to secure my tour plan against all worst case scenarios. Does anyone have an idea how I can make sure that all possible scenarios are considered? So far, my optimization does calculate a random deviation, but this only affects one scenario and thus does not result in a robust tour plan.

In detail, this means, for example: Transportation time t[1,2] = 300 in the normal case. In addition, there are 30 minutes of delay - but in each scenario only at Γ edges.

The original language of the program is OPL, but if you have any hints for Python or Excel, I would be very grateful too.

// parameters & sets
int n = ...; // amount of customers 

{int} N = asSet(1..n);  // set of customers

{int} N0 = {0} union N union {n+1}; // set of customer incl. depots

int m = ...; // amount of vehicles

{int} M = asSet(1..m); // set of vehicles

float Q = ...; // capacity

int r[N]= ...; // demand

int b[N]= ...;// deadline 

tuple Kanten {int i ; int j;};

{Kanten} A  = {<i,j> | i in N0, j in N0}; // Kantenmenge A

execute {writeln (A);}

//{Kanten} B = {<i,j> | i in N0, j in N0}; // Kantenmenge B
//execute {writeln (B);}

int K = 9999;   // large number
float d[N0,N0]= ...;            // max deviation of travel time

float o [N]= ...;           // max deviation of demand

float G = ...;  // uncertainty budget of travel time                

float L = ...;  // uncertainty budget of demand


float t[N0, N0] = ...;              // travel time matrix

{float} td ={t[i,j]  * (1+d[i,j]) | i,j in N0};                     // Depotverbindung ausblenden
    float proba[N0,N0]=...;                                         // wie i,j == j,i?

int c[N0, N0] = ...;                // distance matrix


 // optimization model
// decision variables
dvar boolean x[N0, N0, M];
dvar float+ s[N0, M];
dvar float+ p[N,M];
dvar float+ z[M];
dvar boolean u[N0,N0];

execute{                                        // indicator function
    for (var i in N0){
      for (var j in N0){
      
        if (t[i][j] == td[i][j]){
          u[i][j] == 0;
        }
        else {
          u[i][j] == 1;
        }
      }
    }
};

// objective function
minimize sum(k in M) sum(i in N0) sum(j in N0) c[i,j]*x[i,j,k];


// constraints
subject to {
  // (2)
  NB2:forall (i in N) {
  sum (k in M) sum (j in N: j!=i) x[i,j,k] == 1;
  };
  
  // (3) 
  NB3:forall (k in M) {
   sum (j in N0: j !=0) x[0,j,k] == 1;
  };

  // (4)
  NB4:forall (i in N, k in M) {
   sum (j in N0) x[i,j,k] - sum (j in N0: j!=i) x[j,i,k] == 0;
  };   
  
  // (5)
  NB5:forall (k in M) {
     sum (i in N0: i!=6) x[i,(n+1),k] == 1;
  };
  
   // (11)
   NB11: forall (k in M) {
    sum (i in N) r[i] * sum (j in N0: j!=i) x[i,j,k] + L*z[k] + sum (i in N) p[i,k] <= Q ;
  };    
  // (12)
  NB12: forall (k in M, i,j in N: i!=j) {
    z[k] + p[i,k] >= o[i] * t[i,j]  * sum (j in N0) x[i,j,k];
  };    
  // (13)
  NB13: forall (i,j in N0: i!=j, k in M) {                                          // if else Befehl für td und t?
    s[i,k] + t[i,j] + d[i,j] * u [i,j] - K * (1 - x [i,j,k]) <= s[j,k]  ;       // Wie A + B berücksichtigen? U berücksichtigen notwendig?
  };
  
  /*execute{                                        // indicator function
    for (var i in N0){
      for (var j in N0){
        for (var k in M){
      
        if (t[i][j] == td[i][j]){
          s[i,k] + t[i,j] - K * (1 - x [i,j,k]) <= s[j,k];
        }
        else {
          s[i,k] + t[i,j] + d[i,j] * u [i,j] - K * (1 - x [i,j,k]) <= s[j,k]
        }
      }
    }
  }    
};*/
  
  
  
  // (14)
  NB14: forall (i in N, k in M) {                                   // N oder N0?
    0 <= s[i,k] <= b[i];  
  };    
  //eigene Ergänzungen
  /*  // (15)                                                                           
    NB15: forall (k in M, i in N0, j in N0: i==j) {
   x[i,j,k] == 0;  
  };
    // (16)                                                     
    NB16: forall (i in N0, k in M) {
      s[0,k] == 0;
    };
    // (17)
    NB17: forall (i,j in N0: i== j) {
      u[i,j] == 0;
    }
    // (18)
    NB18: forall (i,j in N0, k in M){
      x[0,6,k] == 0;
    }*/
};

Answer 1

a tiny robust optimization model to start with https://github.com/AlexFleischerParis/zooopl/blob/master/zoorobust.mod in Easy optimization

int nbKids=300;

{int} nbKidsScenarii={nbKids+i*10 | i in -10..2};
float proba[nbKidsScenarii]=[ 1, 1, 2, 2, 2 ,3 ,3 ,4,  5 ,10 ,50 ,10, 7];

assert sum(s in nbKidsScenarii) proba[s]==100; // total probability is 100

float costBus40=500;
float costBus30=400;

int alpha=80; // We want the constraint to be ok with probability 0.95


dvar int+ nbBus40;
dvar int+ nbBus30;


 
minimize
 costBus40*nbBus40  +nbBus30*costBus30;
 
subject to
{
 ctKids:alpha<=sum(nbKids in nbKidsScenarii)
   (40*nbBus40+nbBus30*30>=nbKids)*proba[nbKids];
}

execute
{
writeln("nbBus40 = ",nbBus40);
writeln("nbBus30 = ",nbBus30);
}