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:

#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;

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.
Trying to download the example script gives an error 404.
Does the script also include the ggplot code for making the plot? If not would also include that code, please.
I just fixed the link. The code to generate the plot is included in the file.
Thanks.
Hi Bart, Thanks for the blog post.
Unfortunately I still get an error when downloading.
ie
Simple dowload monitor error
requested file uploads/lpsolve_example.R not found.
BR
Richard
Hey Richard,
Indeed, the link wasn’t working. Had something to do with a new plugin I installed. It should work now.
Grtz, B.
Hi
I’m trying to get my head round your code enough to generalise it… For the line:
lpmodel<-make.lp(2*NROW(train)+NROW(cargo),12)
I'm guessing this is (pseudocode?!;-):
lpmodel<-make.lp((NCOL(train)-1)*NROW(train)+NROW(cargo),NCOL(train)*NCOL(cargo))
To clarify, I guess I’m after something like:
easy_solve(train,cargo,’max’)
Are there any packages that work at this level of abstraction, i.e. that simple require the data to be shaped in an appropriate way, and then managed (if required) through function parameters?
Hi Tony,
I would suggest you take a look at pulp-or (http://code.google.com/p/pulp-or/). This is a python module which basically takes care of representating problem as a matrix for you. The only thing you have to do is declare the objective and the constraints. E.g. your objective could be to maximize the value of the sum of cargo in wagon1 to I, while adhering to the constraint that every individual wagon can only hold volume X. There are a lot of examples on the pulp-or site.
Grtz,
Bart
This is really a terrific example. One thing I can’t understand is how/why the solution result is integer rather than real number?
It would seem (from reviewing the the lpsolveAPI.pdf on CRAN) that the set.type function would need to be called in order to define the decision variables as type “integer”.
In your code, there is no call to set.type. In fact, when I call get.type, it shows all decision variables as “real”.
Yet, the solution that results matches the answer you provide above (in which all are integers).
So, my question is WHY does this code work as an example of an integer program rather than a simple linear one (e.g. if fractional qty’s were allowed in the solution)? Especially given that the default type of the decision vars is “real” and not “integer”.
Many thanks – Josh
Pingback: glass spigot
Nice post. I enjoyed. And nice blog as well. Keep up the good stuff.
Thanks!