14.5 Using R to solve Linear Optimization
The most difficult part about using R to solve a linear optimization problem is to translate the optimization problem into code. Let’s reproduce the table with all the necessary information for the example of Farmer Jean:
| Decision Variables: \(X_1\) = number of plots of parsnips grown \(X_2\) = number of plots of kale grown |
Maximize Profits = 0.15 \(X_1\) + 0.40 \(X_2\) |
|---|---|
| subject to: | - |
| Budget Constraints | \(0.20X_1 + 0.70 X_2 \leq 100\) |
| Space Constraints | \(X_1 + X_2 \leq 200\) |
| Non-negativity | \(X_1 \geq 0\) \(X_2 \geq 0\) |
Here’s how you translate it into code. First, we define the objective function parameters, which are just the coefficients of \(X_1\) and \(X_2\) in the object function: Profits = 0.15 \(X_1\) + 0.40 \(X_2\)
Next, we define the constraints, which are broken up into 3 variables: the constraints matrix, the constraint directions, and the constraint values (or constraint RHS for right-hand-side).
\[0.20X_1 + 0.70 X_2 \leq 100\] \[X_1 + X_2 \leq 200\]
The constraint matrix is simply concatenating all the coefficients here into a matrix. For simplicity, using the convention below, we just read off the matrix from top to bottom, then within each line, from left to right. So we have: (0.20, 0.70, 1, 1). When I construct the matrix, you have to specify the number of columns (ncol) which is simply the number of decision variables we have; in this case it’s 2.
The constraint directions are just a vector that corresponds to each constraint (\(\leq\), \(\leq\)), and the constraint right-hand-side values are just (100, 200)
const.mat <- matrix(c(0.20, 0.70, 1, 1) , ncol=2 , byrow=TRUE)
const.dir <- c("<=", "<=")
const.rhs <- c(100, 200)The important thing to note is to get the order of all the constraints correct. In particular, if you have certain constraints that do not include ALL of the decision variables, to include 0s whenever appropriate.
TIP! For example, if you have an additional constraint that you can only produce a maximum of 500 \(X_1\), this constraint translates to: \(X_1 \leq 500\)
but to be even more helpful to yourself, write it out as: \[ 1 X_1 + 0 X_2 \leq 500 \] This will help you remember to put … 1,0 … into the relevant row of the constraints matrix when you are reading the matrix off from your table.
Finally, we put all of these into a call to the lp function within the lpSolve package. We specify max for maximizing the objective function, pass in the rest of the parameters we just defined, and finally we also ask it to compute.sens=TRUE: we need this for the sensitivity analysis in the next section.
Putting it all together, and looking at the solution (lp.solution$solution) and objective function value (lp.solution$objval)
# defining parameters
objective.fn <- c(0.15, 0.40)
const.mat <- matrix(c(0.20, 0.70, 1, 1) , ncol=2 , byrow=TRUE)
const.dir <- c("<=", "<=")
const.rhs <- c(100, 200)
# solving model
lp.solution <- lp("max", objective.fn, const.mat,
const.dir, const.rhs, compute.sens=TRUE)
# check if it was successful; it also prints out the objective fn value
lp.solution## Success: the objective function is 60## [1] 80 120# Printing it out:
cat("The optimal solution is:", lp.solution$solution, "\nAnd the optimal objective function value is:", lp.solution$objval)## The optimal solution is: 80 120
## And the optimal objective function value is: 60Thus, the optimal solution is \(X_1\) = 80, \(X_2\) = 120, and the optimal profit is 60, which is what we found manually in the previous section.