Optimization with minimum order quantities using R

I'm new to optimization and I need to implement it in a simple scenario:

There exists a car manufacturer that can produce 5 models of cars/vans. Associated with each model that can be produced is a number of labor hours required and a number of tons of steel required, as well as a profit that is earned from selling one such car/van. The manufacturer currently has a fixed amount of steel and labor available, which should be used in such a way that it optimizes total profit.

Here's the part I'm hung up on - each car also has a minimum order quantity. The company must manufacture a certain number of each model before it becomes economically viable to produce/sell that model. This would be easily sent to optim() if it were not for that final condition because the `lower = ...' argument can be given a vector with the minimum order quantities, but then it does not consider 0 as an option. Could someone help me solve this, taking into account the minimum order, but still allowing for an order of 0? Here's how I've organized the relevant information/constraints:

Dorian <- data.frame(Model = c('SmCar', 'MdCar', 'LgCar', 'MdVan', 'LgVan'),
                   SteelReq = c(1.5,3,5,6,8), LabReq=c(30,25,40,45,55),
                   MinProd = c(1000,1000,1000,200,200),
                   Profit = c(2000,2500,3000,5500,7000))

Materials <- data.frame(Steel=6500,Labor=65000)

NetProfit<-function(x) {
  x[1]->SmCar
  x[2]->MdCar
  x[3]->LgCar
  x[4]->MdVan
  x[5]->LgVan
  np<-sum(Dorian$Profit*c(SmCar,MdCar,LgCar,MdVan,LgVan))
  np
}
LowerVec <- Dorian$MinProd #Or 0, how would I add this option?
UpperVec <- apply(rbind(Materials$Labor/Dorian$LabReq, 
                   Materials$Steel/Dorian$SteelReq),2,min)
# Attempt at using optim()
optim(c(0,0,0,0,0),NetProfit,lower=LowerVec, upper=UpperVec)

Eventually I would like to substitute random variables with known distributions for parameters such as Profit and LabReq (labor required) and wrap this into a function that will take Steel and Labor available as inputs as well as parameters for the random variables. I will want to simulate many times and then find the average solution given specific parameters for the Profit and Labor Required, so ideally this optimization would also be fast so that I could perform the simulations. Thanks in advance for any help!

Solution

If you are not familiar with Linear Programming, start here: http://en.wikipedia.org/wiki/Linear_programming

Also have a look at the part about Mixed-Integer Programming http://en.wikipedia.org/wiki/Mixed_integer_programming#Integer_unknowns. That's when the variables you are trying to solve are not all continuous, but also include booleans or integers.

To all aspects, your problem is a mixed-integer programming (to be exact, an integer programming) as you are trying to solve for integers: the number of vehicles to produce for each model.

There are known algorithms for solving these and thankfully, they are already wrapped into R packages for you. Rglpk is one of them, and I'll show you how to formulate your problem so you can use its Rglpk_solve_LP function.

Let x1, x2, x3, x4, x5 be the variables you are solving for: the number of vehicles to produce for each model.

Your objective is:

Profit = 2000 x1 + 2500 x2 + 3000 x3 + 5500 x4 + 7000 x5.

Your steel constraint is:

1.5 x1 + 3 x2 + 5, x3 + 6 x4 + 8 x5 <= 6500

Your labor constraint is:

30 x1 + 25 x2 + 40 x3 + 45 x4 + 55 x5 <= 65000

Now comes the hard part: modeling the minimum production requirements. Let's take the first one as an example: the minimum production requirement on x1 requires that at least 1000 vehicles be produced (x1 >= 1000) or that no vehicle be produced at all (x1 = 0). To model that requirement, we are going to introduce a boolean variables z1. By boolean, I mean z1 can only take two values: 0 or 1. The requirement can be modeled as follows:

1000 z1 <= x1 <= 9999999 z1

Why does this work? Consider the two possible values for z1:

if z1 = 0, then x1 is forced to 0
if z1 = 1 then x1 is forced to be greater than 1000 (the minimum production requirement) and smaller than 9999999 which I picked as an arbitrarily big number.

Repeating this for each model, you will have to introduce similar boolean variables (z2, z3, z4, z5). In the end, the solver will not only be solving for x1, x2, x3, x4, x5 but also for z1, z2, z3, z4, z5.

Putting all this into practice, here is the code for solving your problem. We are going to solve for the vector x = (x1, x2, x3, x4, x5, z1, z2, z3, z4, z5)

library(Rglpk)

num.models <- nrow(Dorian)

# only x1, x2, x3, x4, x5 contribute to the total profit
objective  <- c(Dorian$Profit, rep(0, num.models))

constraints.mat <- rbind(
    c(Dorian$SteelReq, rep(0, num.models)),                    # total steel used
    c(Dorian$LabReq,   rep(0, num.models)),                    # total labor used
    cbind(-diag(num.models), +diag(Dorian$MinProd)),           # MinProd_i * z_i
    cbind(+diag(num.models), -diag(rep(9999999, num.models)))) # x_i - 9999999 x_i

constraints.dir <- c("<=",
                     "<=",
                     rep("<=", num.models),
                     rep("<=", num.models))

constraints.rhs <- c(Materials$Steel,
                     Materials$Labor,
                     rep(0, num.models),
                     rep(0, num.models))

var.types <- c(rep("I", num.models),  # x1, x2, x3, x4, x5 are integers
               rep("B", num.models))  # z1, z2, z3, z4, z5 are booleans

Rglpk_solve_LP(obj   = objective,
               mat   = constraints.mat,
               dir   = constraints.dir,
               rhs   = constraints.rhs,
               types = var.types,
               max   = TRUE)

# $optimum
# [1] 6408000
#
# $solution
# [1] 1000    0    0  202  471    1    0    0    1    1
#
# $status
# [1] 0

So the optimal solution is to create (1000, 0, 0, 202, 471) vehicles of each respective model, for a total profit of 6,408,000.