Intoduction

In this project we want to optimize a distribution network. We will create a script where the user enters

and the script compiles the CPLEX code by itself and then runs it. Alternatively, the user can also run the CPLEX code produced in IBM’s ILOG CPLEX Optimization Studio.

The problem

The problem is as follows.

A company has 3 factories. The location of each factory and the maximum monthly supply is given below:

  • Tilbury 600 tons/month
  • Florence 270 tons/month
  • Barcelona 350 tons/ month

Furthermore the company has warehouses in multiple locations. The location of each supply node, and the demand is given below:

  • Manchester 75 tons/month
  • Tilbury 200 tons/month
  • Paris 250 tons/month
  • Madrid 150 tons/month
  • Milan 100 tons/month
  • Munich 125 tons/month
  • Reykjavik 25 tons/month
  • Sofia 47 tons/month
  • Frankfurt 111 tons/month
  • Prague 31 tons/month

Determine the optimal quantities (and associated costs) to be shipped between plants and distribution centres.

Load the locations of the Supply and Demand Nodes

initial_data<-read.csv("initial_data.csv")
initial_data
##       SupplyNodes DemandNodes
## 1         Tilbury  Manchester
## 2        Florence     Tilbury
## 3 Barcelona+Spain       Paris
## 4                      Madrid
## 5                      Munich
## 6                   Reykjavik
## 7                       Sofia
## 8                   Frankfurt
## 9                      Prague

Calculate the driving distances between each supply and demand node

We need an API key for Google Maps.

Set up your own, correct, API key on the second line of the code below. (Notice the errors due to wrong API key)

require(gmapsdistance)
set.api.key("MY-GOOGLE-MAPS-API-KEY")

Factories<-as.character(initial_data[!(initial_data$SupplyNodes)=="",1])
DistributionCentres<-as.character(initial_data[!(initial_data$DemandNodes)=="",2])
results <- gmapsdistance(origin = Factories,destination = DistributionCentres, mode="driving")
## Error in gmapsdistance(origin = Factories, destination = DistributionCentres, : Google API returned an error: The provided API key is invalid.
distances<- results$Distance
## Error in eval(expr, envir, enclos): object 'results' not found
write.csv(distances,"distances.csv")
## Error in is.data.frame(x): object 'distances' not found

This is how the distances file should look like

distances
##                or Distance.Manchester Distance.Tilbury Distance.Paris
## 1         Tilbury              364719                0         456583
## 2        Florence             1934033          1593591        1153877
## 3 Barcelona+Spain             1825109          1484667        1037337
##   Distance.Madrid Distance.Milan Distance.Munich Distance.Reykjavik
## 1         1716265        1255619         1134963            3665339
## 2         1683893         305015          647098            4198205
## 3          623860         982069         1369384            4574952
##   Distance.Sofia Distance.Frankfurt Distance.Prague
## 1        2478338             762476         1269039
## 2        1400361             968701         1035411
## 3        2375835            1331419         1720865

Cost matrix & Constraints

After doing our research we have calculated the cost per ton per 100km in Europe. This is roughly equal to 2.29 euros.

In the following script we:

  • Create a distance and a cost matrix, using the calculated distances. Two .csvs will be created (distance_matrix_for_report.csv & cost_matrix_for_report.csv)
  • Create two files for the constraints (SupplierConstraints.csv & DemandConstraints.csv), one for the supply nodes and one for the demand nodes, Afterwards we ask the user to fill in the constraints and press enter, so the CPLEX code can be written.
    • Each of these files has 3 columns: NodeName, Demand or Supply and Direction. The user needs to feel the demand or supply, and the direction of the constraint (<, >, <= or >=).
#Here we declare the filenames so we can easily change them, without having to look all over the code
DISTANCES_FILE<-"distances.csv"
FACTORYCONSTRAINTS_FILE<-"SupplierConstraints.csv"
DEMANDNODECONSTRAINTS_FILE<-"DemandConstraints.csv"
CPLEX_FILENAME<-"filename3.lp"
SOLUTION_FILENAME<-"solution.txt"


# cost of freight per ton per 100km in euros
cost<-2.29

#### Create the cost matrix
distances<-read.csv(DISTANCES_FILE)
costmatrix<-distances[, 2:ncol(distances)]/100000*cost


#Create vectors with the correct variable names,eliminate countries (also gmapdistance inserts "Distance." to the names)
headers<-as.vector(names(costmatrix))
newheaders<-sapply(strsplit(headers, split='.', fixed=TRUE), function(x) (x[2]))
newheaders<-sapply(strsplit(newheaders, split='+', fixed=TRUE), function(x) (x[1]))

Factories<-as.vector(distances[,1])
Factories<-sapply(strsplit(Factories, split='+', fixed=TRUE), function(x) (x[1]))
DemandNodes<-as.vector(newheaders)

# renaming and writing the matrixes for our report 
names(distances)[2:ncol(distances)]<-c(newheaders)
names(costmatrix)<-c(newheaders)
costmatrixtowrite<-cbind(distances[,1],costmatrix)
names(distances)[1]<-c("to")
names(costmatrixtowrite)[1]<-c("to")
write.csv(distances,"distance_matrix_for_report.csv",row.names=FALSE)
write.csv(costmatrixtowrite,"cost_matrix_for_report.csv",row.names=FALSE)


for (i in 1:length(Factories)){
  for (j in 1:length(DemandNodes)){
assign(paste(Factories[i], "_to_", DemandNodes[j],sep=""),costmatrix[i,j])
  }
}

#Create Constraints Input Files
SupplierConstraints<-data.frame(Factories)
Names<-c("Production","Direction")
SupplierConstraints[,Names]<-NA

DemandConstraints<-data.frame(DemandNodes)
Names<-c("Demand","Direction")
DemandConstraints[,Names]<-NA

write.csv(SupplierConstraints,"SupplierConstraints.csv",row.names=FALSE)
write.csv(DemandConstraints,"DemandConstraints.csv",row.names=FALSE)

#Ask User to Fill the Constraints Files and Press Enter

readline(prompt="Please fill the constraint values correctly and press [Enter] to continue")
## Please fill the constraint values correctly and press [Enter] to continue
## [1] ""

We now wait for the user to fill in the files created, and then press enter.

This is how the files currently look like

SupplierConstraints
##   Factories Production Direction
## 1   Tilbury         NA        NA
## 2  Florence         NA        NA
## 3 Barcelona         NA        NA
DemandConstraints
##    DemandNodes Demand Direction
## 1   Manchester     NA        NA
## 2      Tilbury     NA        NA
## 3        Paris     NA        NA
## 4       Madrid     NA        NA
## 5        Milan     NA        NA
## 6       Munich     NA        NA
## 7    Reykjavik     NA        NA
## 8        Sofia     NA        NA
## 9    Frankfurt     NA        NA
## 10      Prague     NA        NA

And this is how they should look like when they are filled.

constraints1
##   Factories Production Direction
## 1   Tilbury        600        <=
## 2  Florence        270        <=
## 3 Barcelona        350        <=
constraints2
##    DemandNodes Demand Direction
## 1   Manchester     75         =
## 2      Tilbury    200         =
## 3        Paris    250         =
## 4       Madrid    150         =
## 5        Milan    100         =
## 6       Munich    125         =
## 7    Reykjavik     25         =
## 8        Sofia     47         =
## 9    Frankfurt    111         =
## 10      Prague     31         =

Write the CPLEX code and create an .LP file

Writing in CPLEX is quite tedious. However we can program R to compile the CPLEX code for us. We will be using the sink and cat commands for this.

Isn’t it wonderful?

##### Create CIPLEX_LP file
##### Form the objective function

#First create the variable names
varnames<-data.frame()
for (i in 1:length(Factories)){
  for (j in 1:length(DemandNodes)){
    varnames[i,j]<-(paste(Factories[i], "_to_", DemandNodes[j],sep=""))
  }
}

#Then create the function
minthis<-NA
for (i in 1:length(Factories)){
  for (j in 1:length(DemandNodes)){
minthis<-paste(minthis,"+",costmatrix[i,j],varnames[i,j])
  }
}
minthis<-sapply(strsplit(minthis, split='NA +', fixed=TRUE), function(x) (x[2]))

##### Create the Constraints
# Factory Constraints
constraints1<-read.csv(FACTORYCONSTRAINTS_FILE)
constraint<-as.vector(NA)
for (i in 1:length(Factories)){
  for (j in 1:length(DemandNodes)){
constraint[i]<-  (paste(constraint[i]," + ",Factories[i], "_to_", DemandNodes[j],sep=""))
 }
  }
for (i in 1:length(Factories)){
constraint[i]<- paste(constraint[i],constraints1[i,3],constraints1[i,2])
}

constraint<-sapply(strsplit(constraint, split="NA + ", fixed=TRUE), function(x) (x[2]))

# Supply Node Constraint
constraints2<-read.csv(DEMANDNODECONSTRAINTS_FILE)
constraint2<-as.vector(NA)
for (i in 1:length(DemandNodes)){
  for (j in 1:length(Factories)){
    constraint2[i]<-  (paste(constraint2[i]," + ",Factories[j], "_to_", DemandNodes[i],sep=""))
  }
}
for (i in 1:length(DemandNodes)){
  constraint2[i]<- paste(constraint2[i],constraints2[i,3],constraints2[i,2])
}

constraint2<-sapply(strsplit(constraint2, split="NA + ", fixed=TRUE), function(x) (x[2]))


constraints<-as.vector(rbind(constraint,constraint2))



#Output to file
sink(CPLEX_FILENAME)
cat("Minimize\n")
cat(paste(" obj:",minthis,"\n"))
cat("Subject To \n")
c<-as.vector(0)
for (i in 1:length(constraints)){
  c[i]<-as.vector(paste("c",i,":",sep=""))
  cat(paste(" ",c[i],constraints[i],"\n"))
}
cat("End")
sink(NULL)

LP File

This is how the .LP file looks

Minimize
 obj:  8.3520651 Tilbury_to_Manchester + 0 Tilbury_to_Tilbury + 10.4557507 Tilbury_to_Paris + 39.3024685 Tilbury_to_Madrid + 28.7536751 Tilbury_to_Milan + 25.9906527 Tilbury_to_Munich + 83.9362631 Tilbury_to_Reykjavik + 56.7539402 Tilbury_to_Sofia + 17.4607004 Tilbury_to_Frankfurt + 29.0609931 Tilbury_to_Prague + 44.2893557 Florence_to_Manchester + 36.4932339 Florence_to_Tilbury + 26.4237833 Florence_to_Paris + 38.5611497 Florence_to_Madrid + 6.9848435 Florence_to_Milan + 14.8185442 Florence_to_Munich + 96.1388945 Florence_to_Reykjavik + 32.0682669 Florence_to_Sofia + 22.1832529 Florence_to_Frankfurt + 23.7109119 Florence_to_Prague + 41.7949961 Barcelona_to_Manchester + 33.9988743 Barcelona_to_Tilbury + 23.7550173 Barcelona_to_Paris + 14.286394 Barcelona_to_Madrid + 22.4893801 Barcelona_to_Milan + 31.3588936 Barcelona_to_Munich + 104.7664008 Barcelona_to_Reykjavik + 54.4066215 Barcelona_to_Sofia + 30.4894951 Barcelona_to_Frankfurt + 39.4078085 Barcelona_to_Prague 
Subject To 
  c1: Tilbury_to_Manchester + Tilbury_to_Tilbury + Tilbury_to_Paris + Tilbury_to_Madrid + Tilbury_to_Milan + Tilbury_to_Munich + Tilbury_to_Reykjavik + Tilbury_to_Sofia + Tilbury_to_Frankfurt + Tilbury_to_Prague <= 600 
  c2: Tilbury_to_Manchester + Florence_to_Manchester + Barcelona_to_Manchester = 75 
  c3: Florence_to_Manchester + Florence_to_Tilbury + Florence_to_Paris + Florence_to_Madrid + Florence_to_Milan + Florence_to_Munich + Florence_to_Reykjavik + Florence_to_Sofia + Florence_to_Frankfurt + Florence_to_Prague <= 270 
  c4: Tilbury_to_Tilbury + Florence_to_Tilbury + Barcelona_to_Tilbury = 200 
  c5: Barcelona_to_Manchester + Barcelona_to_Tilbury + Barcelona_to_Paris + Barcelona_to_Madrid + Barcelona_to_Milan + Barcelona_to_Munich + Barcelona_to_Reykjavik + Barcelona_to_Sofia + Barcelona_to_Frankfurt + Barcelona_to_Prague <= 350 
  c6: Tilbury_to_Paris + Florence_to_Paris + Barcelona_to_Paris = 250 
  c7: Tilbury_to_Manchester + Tilbury_to_Tilbury + Tilbury_to_Paris + Tilbury_to_Madrid + Tilbury_to_Milan + Tilbury_to_Munich + Tilbury_to_Reykjavik + Tilbury_to_Sofia + Tilbury_to_Frankfurt + Tilbury_to_Prague <= 600 
  c8: Tilbury_to_Madrid + Florence_to_Madrid + Barcelona_to_Madrid = 150 
  c9: Florence_to_Manchester + Florence_to_Tilbury + Florence_to_Paris + Florence_to_Madrid + Florence_to_Milan + Florence_to_Munich + Florence_to_Reykjavik + Florence_to_Sofia + Florence_to_Frankfurt + Florence_to_Prague <= 270 
  c10: Tilbury_to_Milan + Florence_to_Milan + Barcelona_to_Milan = 100 
  c11: Barcelona_to_Manchester + Barcelona_to_Tilbury + Barcelona_to_Paris + Barcelona_to_Madrid + Barcelona_to_Milan + Barcelona_to_Munich + Barcelona_to_Reykjavik + Barcelona_to_Sofia + Barcelona_to_Frankfurt + Barcelona_to_Prague <= 350 
  c12: Tilbury_to_Munich + Florence_to_Munich + Barcelona_to_Munich = 125 
  c13: Tilbury_to_Manchester + Tilbury_to_Tilbury + Tilbury_to_Paris + Tilbury_to_Madrid + Tilbury_to_Milan + Tilbury_to_Munich + Tilbury_to_Reykjavik + Tilbury_to_Sofia + Tilbury_to_Frankfurt + Tilbury_to_Prague <= 600 
  c14: Tilbury_to_Reykjavik + Florence_to_Reykjavik + Barcelona_to_Reykjavik = 25 
  c15: Florence_to_Manchester + Florence_to_Tilbury + Florence_to_Paris + Florence_to_Madrid + Florence_to_Milan + Florence_to_Munich + Florence_to_Reykjavik + Florence_to_Sofia + Florence_to_Frankfurt + Florence_to_Prague <= 270 
  c16: Tilbury_to_Sofia + Florence_to_Sofia + Barcelona_to_Sofia = 47 
  c17: Barcelona_to_Manchester + Barcelona_to_Tilbury + Barcelona_to_Paris + Barcelona_to_Madrid + Barcelona_to_Milan + Barcelona_to_Munich + Barcelona_to_Reykjavik + Barcelona_to_Sofia + Barcelona_to_Frankfurt + Barcelona_to_Prague <= 350 
  c18: Tilbury_to_Frankfurt + Florence_to_Frankfurt + Barcelona_to_Frankfurt = 111 
  c19: Tilbury_to_Manchester + Tilbury_to_Tilbury + Tilbury_to_Paris + Tilbury_to_Madrid + Tilbury_to_Milan + Tilbury_to_Munich + Tilbury_to_Reykjavik + Tilbury_to_Sofia + Tilbury_to_Frankfurt + Tilbury_to_Prague <= 600 
  c20: Tilbury_to_Prague + Florence_to_Prague + Barcelona_to_Prague = 31 
End

Solution

We now got our file, named “filename3.lp”. It contains the code needed to solve our problem, written in CPLEX, assembled in R.

We need to load it either in IBM’s optimizer and hit run OR use the package “Rglpk” in R.

We will use R.

# Solve the problem

require(Rglpk)
x<-Rglpk_read_file(file=CPLEX_FILENAME, type = "CPLEX_LP")
solution<-Rglpk_solve_LP(x$objective,mat=x$constraints[[1]],dir=x$constraints[[2]],rhs=x$constraints[[3]],x$bounds, x$types, x$maximum)

# Solution Output 

require("gdata")
sink(SOLUTION_FILENAME)
vars<-as.vector(unmatrix(varnames,byrow=T))
costs<-as.vector(unmatrix(costmatrix,byrow=T))

cat(paste("Total Cost (euros):",solution$optimum,"\n\n"))
for (i in 1:length(vars)){
cat(paste(vars[i],"=",solution$solution[i],"cost:",solution$solution[i]*costs[i],"\n"))
}
sink(NULL)

Clean solution

That’s it our solution is dumbed into a text file.

However it also contains the arcs/routes with 0 freight.

If we want a clean solution, then we can slightly change the above code and only include the arcs with a freight different than 0:

vars<-as.vector(unmatrix(varnames,byrow=T))
costs<-as.vector(unmatrix(costmatrix,byrow=T))

only_include<-which(solution$solution!=0)

vars<-vars[only_include]
costs<-costs[only_include]
clean_solution<-solution$solution[only_include]

cat(paste("Total Cost (euros):",solution$optimum,"\n\n"))
for (i in 1:length(vars)){
cat(paste(vars[i],"=",clean_solution[i],"cost:",clean_solution[i]*costs[i],"\n"))
}

Results

Total Cost (euros): 15519.3013521 

Tilbury_to_Manchester = 75 cost: 626.4048825 
Tilbury_to_Tilbury = 200 cost: 0 
Tilbury_to_Paris = 250 cost: 2613.937675 
Tilbury_to_Reykjavik = 25 cost: 2098.4065775 
Tilbury_to_Frankfurt = 50 cost: 873.03502 
Florence_to_Milan = 67 cost: 467.9845145 
Florence_to_Munich = 125 cost: 1852.318025 
Florence_to_Sofia = 47 cost: 1507.2085443 
Florence_to_Prague = 31 cost: 735.0382689 
Barcelona_to_Madrid = 150 cost: 2142.9591 
Barcelona_to_Milan = 33 cost: 742.1495433 
Barcelona_to_Frankfurt = 61 cost: 1859.8592011