Linear programming in R: an lpSolveAPI example

July 14, 2012
By

(This article was first published on Fishy Operations, and kindly contributed to R-bloggers)

First of all, a shout out to R-bloggers for adding my feed to their website!

Linear programming is a valuable instrument when it comes to decision making. This post shows how R in conjunction with the lpSolveAPI package, can be used to build a linear programming model and to analyse its results.

The lpSolveAPI package provides a complete implementation of the lp_solve API.

The example case;

A trading company is looking for a way to maximize profit per transportation of their goods. The company has a train available with 3 wagons. When stocking the wagons they can choose between 4 types of cargo, each with its own specifications. How much of each cargo type should be loaded on which wagon in order to maximize profit?

The following constraints have to be taken in consideration;

  • Weight capacity per train wagon
  • Volume capacity per train wagon
  • Limited availability per cargo type

Let's assume we have the following information at hand:

Train wagon Weight capacity (tonne) Space capacity (m2)
w1 10 5000
w2 8 4000
w3 12 8000
Cargo type Available (tonne) Volume (m2) Profit (per tonne)
c1 18 400 2000
c2 10 300 2500
c3 5 200 5000
c4 20 500 3500

Let's load the necessary libraries and define the datasets to be used.

library(lpSolveAPI)

#used for result visualization
library(ggplot2)
library(reshape)
library(gridExtra)

#define the datasets

train<-data.frame(wagon=c('w1','w2','w3'), weightcapacity=c(10,8,12), spacecapacity=c(5000,4000,8000))

cargo<-data.frame(type=c('c1','c2','c3','c4'), available=c(18,10,5,20), volume=c(400,300,200,500),profit=c(2000,2500,5000,3500))

Next, we start with building the actual lp model, our goal is to end up with the following model:
result_matrix

#create an LP model with 10 constraints and 12 decision variables

lpmodel<-make.lp(2*NROW(train)+NROW(cargo),12)

#I used this to keep count within the loops, I admit that this could be done a lot neater
column<-0
row<-0

#build the model column per column
for(wg in train$wagon){
    row<-row+1
    for(type in seq(1,NROW(cargo$type))){
    column<-column+1

    #this takes the arguments 'column','values' & 'indices' (as in where these values should be placed in the column)
    set.column(lpmodel,column,c(1, cargo[type,'volume'],1), indices=c(row,NROW(train)+row, NROW(train)*2+type))
    }}

#set rhs weight constraints
set.constr.value(lpmodel, rhs=train$weightcapacity, constraints=seq(1,NROW(train)))

#set rhs volume constraints
set.constr.value(lpmodel, rhs=train$spacecapacity, constraints=seq(NROW(train)+1,NROW(train)*2))


#set rhs volume constraints
set.constr.value(lpmodel, rhs=cargo$available, constraints=seq(NROW(train)*2+1,NROW(train)*2+NROW(cargo)))

#set objective coefficients
set.objfn(lpmodel, rep(cargo$profit,NROW(train)))

#set objective direction
lp.control(lpmodel,sense='max')

#I in order to be able to visually check the model, I find it useful to write the model to a text file
write.lp(lpmodel,'model.lp',type='lp')

So, let's have a look at the 'model.lp' file.

/* Objective function */
max: +2000 C1 +2500 C2 +5000 C3 +3500 C4 +2000 C5 +2500 C6 +5000 C7 +3500 C8 +2000 C9 +2500 C10 +5000 C11
 +3500 C12;

/* Constraints */
+C1 +C2 +C3 +C4 <= 10;
+C5 +C6 +C7 +C8 <= 8;
+C9 +C10 +C11 +C12 <= 12;
+400 C1 +300 C2 +200 C3 +500 C4 <= 5000;
+400 C5 +300 C6 +200 C7 +500 C8 <= 4000;
+400 C9 +300 C10 +200 C11 +500 C12 <= 8000;
+C1 +C5 +C9 <= 18;
+C2 +C6 +C10 <= 10;
+C3 +C7 +C11 <= 5;
+C4 +C8 +C12 <= 20;

This seems to be the table which we wanted to create! :)

Now to see if the model will solve:

#solve the model, if this return 0 an optimal solution is found
solve(lpmodel)

#this return the proposed solution
get.objective(lpmodel)

solve(lpmodel) returns a 0, this implies that an optimal solutions was found. The value of the solution is returned by get.objective(lpmodel), this shows a maximum total profit of 107500.

Using ggplot2 we generate the following plot;
results_lp

We can conclude from this exercise what the maximum profit is given the current constraints in transportation resources and available cargo types. Furthermore, in the bottom half of the graph we can see that the proposed configuration leads to a maximum utilization of weight on all train wagons. As we have the model available it is easy to run it for alternative configurations, e.g. an extra trainwagon, new cargo types or different profit levels.

The example R script can be downloaded here.

To leave a comment for the author, please follow the link and comment on his blog: Fishy Operations.

R-bloggers.com offers daily e-mail updates about R news and tutorials on topics such as: visualization (ggplot2, Boxplots, maps, animation), programming (RStudio, Sweave, LaTeX, SQL, Eclipse, git, hadoop, Web Scraping) statistics (regression, PCA, time series, trading) and more...



If you got this far, why not subscribe for updates from the site? Choose your flavor: e-mail, twitter, RSS, or facebook...

Comments are closed.