# Coding the OR problem with cplex

I am new to using OPL. I wrote a CPLEX code for the vaccine distribution from a paper, it doesn't get an error, but it hasn't given the answer. I don't know its problem; please help me; I attached the code in the following. I'm so confused and don't know what can be done?

/*********************************************
* OPL 12.8.0.0 Model
* Author: Ramin
* Creation Date: Feb 12, 2021 at 4:32:21 PM
*********************************************/
//sets--------------------------------------------------------------------------------------------------------------------------------------------------------------------
range D=1..3;   {int} DD=asSet(D);    //Distribution centers.
int c=0;        {int} cc={0};         //center.
range I=4..6;   {int} II={4,5,6};     //Local clinics.
range DI=1..6; //D+I
range cD=0..3; //c+D
range J=0..6;   {int} JJ= cc union DD union II;       //Vertices.
{int} L={0,1,2} ;                      //Levels which central 0, DC 1, local clinic 2.
range T=0..2;   {int} TT=asSet(T);     //Transportation vehicles which cold truck ( = 0), 4 × 4 truck ( = 1), motorbike ( = 2).
range R=0..3;   {int} RR= {0,1,2,3};   //Storage devices which cold room ( = 0), regional device ( = 1), district device ( = 2), local clinic device ( = 3).
range F=0..1;   {int} FF={0,1};        //Replenishment frequency: {quarterly ( = 0), monthly ( = 1)}.
//Define Arcs---------------------------------------------------------------------------------------------------------------------------------------------------------------
tuple Arc{
key int origin;
key int destination;
}
{Arc} A1={<0,1>,<0,2>,<0,3>};          //Arcs between center and Distribution centers.
{Arc} A2={<1,4>,<1,5>,<1,6>,<2,4>,<2,5>,<2,6>,<3,4>,<3,5>,<3,6>,<1,2>,<2,3>};   //Arcs between Distribution centers and locla clinics.
{Arc} A3={<1,2>,<2,3>};               //Arcs between Distribution centers.
{Arc} A={<0,1>,<0,2>,<0,3>,<1,4>,<1,5>,<1,6>,<2,4>,<2,5>,<2,6>,<3,4>,<3,5>,<3,6>,<1,2>,<2,3>};  //total Arcs.
//EXTEND//-----------------------------------------------------------------------------------------------------------------------------------------------------------------
///suppose1:motorbike ( = 2) works by electricity///
{int} chargestation={1,2,3};          //chargestation in route.
range chargestationcopy=1..6;   {int} chargestationcopyy={1,2,3,4,5,6}; //we have three shargestation but it may each chargestation visited more than one.

//parameters-----------------------------------------------------------------------------------------------------------------------------------------------------
float C[A][t in T]=rand(10);          //Transportation cost per kilometer of vehicle type t from location i to location j.
float CS[R]=[8116,1582,600,596];      //Annual storage cost per storage device r.
float C0=40000;                       //Annual level 0 facility (central store) operating cost.
float C1[j in D]=4500;                //Annual level 1 facility (DC) operating cost at location j.
float C2[j in I]=800;                 //Annual level 2 facility (clinic) operating cost at location j.
float P[T]= [9293, 172,5];            //Transportation capacity per trip of vehicle t.
float PS[R]=[18000,1843,76,35];       //Storage capacity of device r.
float G[f in F]=(f==0)?4:12;          //Annual number of replenishments f  ( = 4 if f = 0 and = 12 if f = 1).
float S=0.25;                         //Buffer stock factor for vaccines stored at a location (typically = 0.25).
float H[A]= [2,3,2,2,5,6,7,5,8,9,12,13,12,6]; //Distance (kilometer) between location i and location j.
float M=100000000;                    //Big M.
//EXTEND//--------------------------------------------------------------------------------------------------------------------------------------------------------
///suppose1:motorbike ( = 2) works by electricity///-------------------------------------------------
int Beta=1;                           //(1+ beta) shows the number of copy of shargstations.
int Q =200;                           //motorbike battery capacity.
float e[i in cD][j in DI]=rand(10);   //charge consumption of motorbike from i to j.
int Tmax=200;                         //maximum time that motorbike should drive in two nodes.
///suppose2:cold truck ( = 0), 4 × 4 truck ( = 1) work by diesel, for decrease green gas///----------
int GHL[f in F]=100;
///suppose3:define cost of green house gas in objective function/////--------------------------------
float costco2=20;                     //average cost of co2 on the ration of weight.
float wco2=10;                        //weight of generate co2 by burning fuel per liter.
float fuel[T]=[10,5,0];               //volume of burning fuel per weight of vehicle and distance.
float w[T]=[100,200,0];               //weight of vehicle.
float wcurb[T]=[150,240,0];           //weights of vehicle and its load.
int RVC[T]=[3,2,1];                   //ratio of volume of vehicle on weight (without load).

//variables--------------------------------------------------------------------------------------------------------------------------------
dvar float x[A];                      //Annual flow (volume) of vaccines from location i to location j
dvar float Y[J][R][F];                //Number of storage devices of type r at location j with replenishment frequency f
dvar float Z[A][T][F];                //Number of vehicle trips per replenishment from location i to location j using vehicle type t with replenishment
dvar boolean W[J];                    //1 if a location j is open, 0 otherwise
dvar boolean U[A];                    //1 if vaccines flow from location i to location j, 0 otherwise
dvar boolean V[J];                    //1 if a location j has monthly replenishment frequency, 0 otherwise
//EXTEND//-----------------------------------------------------------------------------------------------------------------------------------------
///suppose1:motorbike ( = 2) works by electricity///
dvar float y[J][T];                   //charge level of motorbike in each node.
dvar float q[chargestationcopy][T];   //charge level of motorbike in each chargestation.
dvar float O[chargestationcopy][T];   //charge of motorbike when depart from chargestation.
dvar float delta[chargestationcopy][T];//period of time being in chargestation.
dvar float timeexitchargestationcopy[chargestationcopy][T];//exit time from chargestation.
dvar float timeenterchargestationcopy[chargestationcopy][T];//enter time to chargstation.
dvar float recievetime[J][T];         //recieve time to a nod
dvar float servicetime[J][T];         //service time for costomer or clinic
dvar float time[cD][DI];              //time period of motorbike when drive from obe node to other.
dvar boolean K[cD][DI];               //if motorebike visit i an j be one and else zero.
///suppose2;cold truck ( = 0), 4 × 4 truck ( = 1) work by diesel, for decrease green gas///
dvar float GHG[F];                    //amount of green gas which generate by vehicle.

//formulation------------------------------------------------------------------------------------------------------------------------------------------
dexpr float objective= sum(i in A,t in T,f in F)2* C[i][t]*G[f]*H[i]*Z[i][t][f]+
sum(j in J,f in F,r in R)CS[r]*Y[j][r][f]+ C0*W+ sum(j in D)C1[j]*W[j]+sum(j in I)C2[j]*W[j]+
sum(i in A,t in T:t!=2)H[i]*costco2*wco2*fuel[t]*w[t]*wload[t]+sum(i in A,t in T:t!=2)x[i]*H[i]*wco2*fuel[t]*wcurb[t]/RVC[t]
;

minimize objective;
subject to{
cos1:sum(i in A)x[i]-sum(j in A)x[j]==0;
forall(i,j in A)cos2:sum(t in T,f in F)P[t]*G[f]*Z[i][t][f]>=x[i];
cos3:sum(t in T,r in R, f in F,j in J)P[t]*G[f]*Y[j][r][f]>=(1+S)*sum(i in A)x[i];
forall(j in DD)cos4:M*W[j]>=sum(i in A)x[i]+sum(j in A)x[j];
forall(i in A)cos5:x[i]<=M*U[i];
cos6:sum(i in A)U[i]<=1;
forall(j in D)cos7:sum(r in R,f in F)Y[j][r][f]<=M*W[j];
forall(j in J)cos8:sum(r in R)Y[j][r]<=M*V[j];
forall(j in J)cos9:sum(r in R)Y[j][r]<=M*(1-V[j]);
forall (i in A3,j in J)cos10:x[i]<=M*V[j];
forall (i in A1,k in D)cos11:V[k]<=2-U[i]-(sum(j in A3)U[j]/abs(3));
forall (i in cc union II)cos12:W[i]==1;
forall (i in I)cos13:V[i]==1;
cos14:Y==1;
forall (i in A)cos15:x[i]>=0;
///////////EXTEND////////////////-----------------------------------------------------------------------------------------------------------------------------------
///suppose1:motorbike ( = 2) works by electricity///
forall(i in cD, t in T:t==2)sum(j in DI)K[i][j]<=1;
forall(i in cD,j in DI,t in T:t==2)e[i][j]*K[i][j]-(1-K[i][j])*Q <=y[i][t]-y[j][t];
forall(i in cD,j in DI,t in T:t==2)y[i][t]-y[j][t]<=e[i][j]*K[i][j]+(1-K[i][j])*Q;
forall(i in cD,j in chargestationcopy,t in T:t==2)e[i][j]*K[i][j]-(1-K[i][j])*Q<= y[i][t]-q[j][t];
forall(i in cD,j in chargestationcopy,t in T:t==2)y[i][t]-q[j][t]<=e[i][j]*K[i][j]-(1-K[i][j])*Q;
forall(i in cc,t in T:t==2)y[i][t]==Q;
forall(i in chargestationcopy,t in T:t==2)y[i][t]==O[i][t];
forall(i in chargestationcopy,t in T:t==2)q[i][t]<=O[i][t];
forall(i in chargestationcopy,t in T:t==2)delta[i][t]==timeexitchargestationcopy[i][t]-timeenterchargestationcopy[i][t];
forall(i in cD,j in DI,t in T:t==2)recievetime[i][t]+(time[i][j]+servicetime[j][t])*K[i][j]-Tmax*(1-K[i][j])<=recievetime[j][t];
forall(i in cD:i!=0,j in chargestationcopy,t in T:t==2)recievetime[i][t]+delta[i][t]+time[i][j]*K[i][j]-Tmax*(1-K[i][j])<=recievetime[j][t];
forall(i in cD,j in DI,t in T:t==2)recievetime[j][t]+time[i][j]<=Tmax;
forall(t in T:t==2)recievetime[t]<=Tmax;
///suppose2:cold truck ( = 0), 4 × 4 truck ( = 1) work by diesel, for decreasing green gas///
forall(f in F,t in T:t!=2)sum(i in A)GHG[f]*H[i]*U[i]<=GHL[f];
}


This is the engine log output!!!!! • Perhaps you can show the log output. Feb 15 at 17:42
• I attached the engine log output.@MarkL.Stone Feb 15 at 18:46
• @Alex Fleischer, could you see the code, thanks. Feb 15 at 18:48
• @ramin, I think you would need to check your model constraints. I tried to reproduce the solution and it seems the model does have many conflicts in the constraints. By disabling a bunch of them the solver log is something like this: Found incumbent of value 1.6560489e+09 after 0.00 sec. (0.02 ticks) Dual infeasible due to empty column 'Z({1,4})(1)(0)'. Presolve time = 0.00 sec. (0.05 ticks) Feb 15 at 19:40
• @A.Omidi thanks, you are right. I should check it. I got this formulation from the paper(onlinelibrary.wiley.com/doi/abs/10.1111/itor.12758) and implement some extension on it. thanks Feb 15 at 19:57

if you name other constraints as:

forall(i in cD, t in T:t==2)  cb1:sum(j in DI)K[i][j]<=1;
forall(i in cD,j in DI,t in T:t==2)  cb2:e[i][j]*K[i][j]-(1-K[i][j])*Q <=y[i][t]-y[j][t];
forall(i in cD,j in DI,t in T:t==2) cb3:y[i][t]-y[j][t]<=e[i][j]*K[i][j]+(1-K[i][j])*Q;
forall(i in cD,j in chargestationcopy,t in T:t==2) cb4:e[i][j]*K[i][j]-(1-K[i][j])*Q<= y[i][t]-q[j][t];
forall(i in cD,j in chargestationcopy,t in T:t==2) cb5:y[i][t]-q[j][t]<=e[i][j]*K[i][j]-(1-K[i][j])*Q;
forall(i in cc,t in T:t==2)y[i][t]==Q;
forall(i in chargestationcopy,t in T:t==2)y[i][t]==O[i][t];
forall(i in chargestationcopy,t in T:t==2) cb6:q[i][t]<=O[i][t];


then you ll get some conflicts to look at: • Thanks. that's so help-full. Feb 16 at 11:03