CPLEX: CP Optimizer took forever to run with no solution

Question

I run my model in CPLEX using Constraint Programming for quite a long time, but it still has no solution so I have to stop the model.

Is there something wrong with my model? I tried to run the model with only a few sample in the Excel file but it still failed to give me a solution after hours of run time. Thank you so much in advance!

My mod. file:

using CP;
int Z_val = 8;
int P_val = 5;
int T_val = 20;
int I_val = 4;
int J_val = 5;
int K_val = ...;
int N = ...;// Planning horizon (day);
float alpha = ...;// The average loading rate of delivery trucks;
float gamma = ...;// The maximum of containers waiting in the queue at each block.
int BigM = 1000;
int fact_K_minus_1 = ...;

range Z = 1 .. Z_val;// Number of vessels
range P = 1 .. P_val;// Number of appointment periods
range T = 1 .. T_val;// Number of time intervals
range I = 1 .. I_val;// Number of gate lanes
range J = 1 .. J_val;// Number of yard blocks
range K = 1 .. K_val;// Number of RTGCs

float T_z_A[Z] = ...;
float T_z_D[Z] = ...;
float V_z[Z] = ...;
float beta_zj[Z][J] = ...;

float C_s = ...;
range KSet = 0 .. K_val - 1;
float T_l = ...;
float T_k = ...;

float Y_j[J] = ...;
float c_o = ...;
float c_R = ...;
float factorial[KSet] = ...;

// Define sets
{int} Z_j[J];

execute {

  // Iterate over j and z to update Z sets
  for ( var j in J) {
    for ( var z in Z) {
      if (beta_zj[z][j] > 0) {
        Z_j[j].add(z); // Add z to the set Z_j[j]
      }
    }
  }

  // Print the result (optional)
  for ( var j in J) {
    writeln("Z_j[", j, "] = ", Z_j[j]);
  }
}


int scale = 1000;
dvar int pS[Z] ;// Starting appointment period for delivery trucks related to vessel z
dvar int pE[Z] ;


// Decision variables
dexpr float p_z_S[z in Z] = pS[z] / scale;// Starting appointment period for delivery trucks related to vessel z
dexpr float p_z_E[z in Z] = pE[z] / scale;// Ending appointment period for delivery trucks related to vessel z

// Deived variables
float d = ceil ( T_val / P_val );// Number of time intervals in one appointment period
int m = ftoi ( d );
execute {
  writeln(m);
}


dvar int lam_zp[Z][P]  ;// Appointment quota of export containers related to vessel z at appointment period p
dvar int lit_it[I][T]  ;// Average number of trucks waiting in queue at gate lane i at interval t
dvar int dit_it[I][T]  ;// Actual discharge rate of gate lane i at interval t (truck/min)
dvar int pit_it [I][T]  ;// Capacity utilization rate of gate lane i at interval t
dvar int wp_p [P]  ;// Average waiting time of trucks at terminal gate during appointment period p (min)
dvar int wg_g  ;// Average waiting time of trucks at terminal gate during the planning horizon (min)
dvar int djkt_jkt [J][K][T] ;// Discharge rate of RTGC k deployed to block j at interval t (natural container/min)
dvar int pjkt_jkt [J][K][T] ;// Average utilization rate of RTGC k deployed to block j at interval t
dvar int wjp_jp [J][P] ;// Average waiting time of trucks at block j in appointment period p (min)
dvar int wy_y ;// Average waiting time of trucks at yard in the planning horizon (min)


dexpr float lambda_zp[z in Z][p in P] = lam_zp[z][p] / scale;// Appointment quota of export containers related to vessel z at appointment period p
dexpr float l_it_g[i in I][t in T] = lit_it[i][t] / scale;// Average number of trucks waiting in queue at gate lane i at interval t
dexpr float d_it_g[i in I][t in T] = dit_it[i][t] / scale;// Actual discharge rate of gate lane i at interval t (truck/min)
dexpr float p_it_g[i in I][t in T] = pit_it [i][t] / scale;// Capacity utilization rate of gate lane i at interval t
dexpr float w_p_g[p in P] = wp_p[p] / scale;// Average waiting time of trucks at terminal gate during appointment period p (min)
dexpr float w_g = wg_g/ scale;// Average waiting time of trucks at terminal gate during the planning horizon (min)
dexpr float d_jkt_y[j in J][k in K][t in T] = djkt_jkt[j][k][t] / scale;// Discharge rate of RTGC k deployed to block j at interval t (natural container/min)
dexpr float p_jkt_y[j in J][k in K][t in T] = pjkt_jkt [j][k][t] / scale;// Average utilization rate of RTGC k deployed to block j at interval t
dexpr float w_jp_y[j in J][p in P] = wjp_jp [j][p] / scale;// Average waiting time of trucks at block j in appointment period p (min)
dexpr float w_y = wy_y / scale;// Average waiting time of trucks at yard in the planning horizon (min)
dvar boolean d_zp[Z][P];// Binary variable, 1 if vessel z has departed at appointment period p, 0 otherwise

dvar int  lambda_zt_g[Z][T] ;// Number of trucks related to vessel z arriving at terminal gate at interval t
dvar int lambda_it_g[I][T] ;// Number of trucks arriving at gate lane i at interval t

dvar int lambda_t_y[ T] ;// Number of export containers arriving at yard at interval t
dvar int lambda_jt_y[J][T] ;// Number of export containers arriving at block j at interval t
dvar int l_jt_y[J][T] ;// Number of export containers waiting at block j at interval t

// Objective function
minimize
  ( sum ( t in T ) sum ( i in I ) l_it_g[i][t] + sum ( t in T ) sum ( j in J )
     l_jt_y[j][t] ) * ( 24 * N * c_o / T_val ) + sum ( j in J ) sum ( k in K )
     sum ( t in T ) ( 1 - p_jkt_y[j][k][t] ) * 24 * N * c_R / T_val;
subject to {
  // Constraints 1


    forall ( z in Z )//(1)
    ( p_z_E[z] - p_z_S[z] + 1 ) * 24 * N / P_val >= T_l;

    forall ( z in Z )//(2)
    ( p_z_E[z] - p_z_S[z] + 1 ) * 24 * N / P_val <= T_k;

    forall ( z in Z )//(3)
    p_z_E[z] * 24 * N / P_val <= T_z_A[z];

    forall ( z in Z, p in P )//(4)
    lambda_zp[z][p] * ( p_z_E[z] - p_z_S[z] + 1 ) - V_z[z]<=0.1;

  forall (z in Z, p in P) {
  if ((p - 1) * 24 * N / P_val >= T_z_D[z]) {
    d_zp[z][p] == 1;
  } else {
    d_zp[z][p] == 0;
  }
}

 
    forall ( j in J, p in P )//(6)
    sum ( z in Z_j[j] ) ( sum ( a in 1 .. p ) 
      ( lambda_zp[z][a] - V_z[z] * d_zp[z][a] ) ) * beta_zj[z][j] <= Y_j[j];

  // Constraints at gate
    forall ( z in Z, p in P, t in ( ( p - 1 ) * m + 1 ) .. p * m )//(7)
    lambda_zt_g[z][t] -(lambda_zp[z][p] / ( m * alpha ))<=0.1;

    forall ( i in I, t in T )//(8)
    lambda_it_g[i][t] -( sum ( z in Z ) lambda_zt_g[z][t] / I_val)<=0.1;

    forall ( i in I, t in 1 .. T_val - 1 )//(9)
    l_it_g[i][t + 1] - (maxl ( l_it_g[i][t] + lambda_it_g[i][t] - 
      ( 24 * 60 * N * d_it_g[i][t] / T_val ), 0 ))<=0.1;

    forall ( i in I, t in T )//(10)
    d_it_g[i][t] - ( 19 / 60 ) * p_it_g[i][t]<=0.1;

    forall ( i in I, t in T )//(11)
    //l_it_g[i][t] == p_it_g[i][t] / ( 1 - p_it_g[i][t] );
    l_it_g[i][t] * ( 1 - p_it_g[i][t] ) - p_it_g[i][t] <=0.1;

    forall ( p in P )//(12)
    w_p_g[p] * sum ( t in ( ( p - 1 ) * m + 1 ) .. p * m ) sum ( i in I )
       d_it_g[i][t] - (sum ( t in ( ( p - 1 ) * m + 1 ) .. p * m ) sum 
      ( i in I ) l_it_g[i][t])<=0.1;

  w_g * sum ( t in T ) sum ( i in I ) d_it_g[i][t] - (sum ( t in T ) sum 
  ( i in I ) l_it_g[i][t])<=10;//(13)

  // Constraints at yard
    forall ( t in T )//(14)
    lambda_t_y[t] -( ( 24 * 60 * N / T_val ) * sum ( i in I ) d_it_g[i][t])<=0.1;


    forall ( j in J, t in T )//(15)
    lambda_jt_y[j][t] - (lambda_t_y[t] * alpha * sum ( z in Z_j[j] )
       beta_zj[z][j] * ( lambda_zt_g[z][t] / sum ( i in I ) lambda_it_g[i][t] ))<=0.1;
      

    forall ( j in J, t in 1 .. T_val - 1 )//(16)
    l_jt_y[j][t + 1] - ( maxl ( l_jt_y[j][t] + lambda_jt_y[j][t] - 
      ( 24 * 60 * N / T_val ) * sum ( k in K ) d_jkt_y[j][k][t], 0 ))<=0.1;

  forall ( j in J, t in T, k in K //(17)
  )
    d_jkt_y[j][k][t] - ( 9 / 60 ) * p_jkt_y[j][k][t] <= 0.1;

  forall ( j in J, t in T, k in K //(18)
  )
    l_jt_y[j][t] - (( ( p_jkt_y[j][k][t] * ( 1 + C_s * C_s ) ) / ( 2 * 
      ( K_val - p_jkt_y[j][k][t] ) ) ) * ( 1 + ( fact_K_minus_1 * 
      ( K_val - p_jkt_y[j][k][t] ) ) * ( sum ( n in KSet ) ( 1 / 
      ( factorial[n] * ( p_jkt_y[j][k][t] ^ ( K_val - n ) ) ) ) ) ) ^ ( - 1 )
       + p_jkt_y[j][k][t] )<= 0.1;

    forall ( j in J, t in T )//(19)
    l_jt_y[j][t] <= gamma;

    forall ( j in J, p in P )//(20)
    w_jp_y[j][p] * sum ( t in ( ( p - 1 ) * m + 1 ) .. p * m, k in K )
       d_jkt_y[j][k][t] -( sum ( t in ( ( p - 1 ) * m + 1 ) .. p * m )
       l_jt_y[j][t])<=0.1;

  w_y * sum ( t in T, j in J, k in K ) d_jkt_y[j][k][t] -( sum ( t in T ) sum 
  ( j in J ) l_jt_y[j][t])<=0.1;//(21)



}

My dat. file:

SheetConnection my_sheet("eee.xlsx");
N = 7;
alpha = 1.4;
gamma = 8;
C_s = 0.42687;
T_l= 0.6;
T_k= 2.4;
c_o = 5.728;
c_R = 15.48;

K_val from SheetRead (my_sheet,"Sheet1!B4");
fact_K_minus_1 from SheetRead (my_sheet,"Sheet1!D4");

T_z_A from SheetRead (my_sheet,"Sheet1!B7:B14");
T_z_D from SheetRead (my_sheet,"Sheet1!C7:C14");
V_z from SheetRead (my_sheet,"Sheet1!D7:D14");
Y_j from SheetRead (my_sheet,"Sheet1!L7:L11");
factorial from SheetRead(my_sheet,"Sheet1!I7:I10");

beta_zj from SheetRead (my_sheet,"Sheet1!:B22:I26");

and my excel file here https://docs.google.com/spreadsheets/d/1j4R45RKG53Gaz68DnBP36OLT3vitgcAd/edit#gid=1448362513

Answer 1

I suspect your optimization is unbounded from below. You could check this by limiting the maximum iteration number, noting the objective value, then doubling the maximum iteration number and comparing the newly obtained objective value to the previously noted one.

But just looking at your constraints, I guess the optimizer can take l_it_g to negative infinity (and so the objective value as well) by disabling its constraints the following way:

// maxl chooses 0
forall ( i in I, t in 1 .. T_val - 1 )//(9)
    l_it_g[i][t + 1] - (maxl ( l_it_g[i][t] + lambda_it_g[i][t] - 
      ( 24 * 60 * N * d_it_g[i][t] / T_val ), 0 ))<=0.1;

// p_it_g stays between 0 and 1
forall ( i in I, t in T )//(11)
    //l_it_g[i][t] == p_it_g[i][t] / ( 1 - p_it_g[i][t] );
    l_it_g[i][t] * ( 1 - p_it_g[i][t] ) - p_it_g[i][t] <=0.1;

// w_p_g goes to negative infinity as l_it_g and d_it_g stays around 1
forall ( p in P )//(12)
    w_p_g[p] * sum ( t in ( ( p - 1 ) * m + 1 ) .. p * m ) sum ( i in I )
       d_it_g[i][t] - (sum ( t in ( ( p - 1 ) * m + 1 ) .. p * m ) sum 
      ( i in I ) l_it_g[i][t])<=0.1;

// w_g goes to negative infinity as l_it_g and d_it_g stays around 1
w_g * sum ( t in T ) sum ( i in I ) d_it_g[i][t] - (sum ( t in T ) sum 
  ( i in I ) l_it_g[i][t])<=10;//(13)

Simplifying your model was a good initial step for debugging. The next step is to force the optimizer to stop early (for example by limiting the maximum number of iterations) and try to interpret the partial results whether they make sense...

By the nature of your problem, I guess you are missing a few non-negativity constraints on the decision variables by defining them using int+ and float+.