Linear Programming with R

linear-programming-with-r.rmd

---
title: "Linear Programming with R"
author: "Abhishek Dubey"
date: "15 June 2020"
output: html_document
---

```{r setup, include=FALSE}
knitr::opts_chunk$set(echo = TRUE)
```

## What is Linear Programming?

Linear Programming is the technique of finding an **optimal solution** for a **linear optimization problem** given an **objective function**, **linear constraints** and **non-negative variables**.

In linear programming, a problem is formulated into a **mathematical model** consisting of the objective function to be maximised or minimised and a system of linear inequalities as constraints. The solution of this system is a polyhedron known as the **feasible region**. Then the objective function is evaluated at all the vertices of the polyhedron. The largest and smallest values are the maximum and minimum respectively.

Linear Programming can be used to solve various scheduling and optimization problems. It can be used for the allocation of limited resources to an operation to several competing activities, finding maximum profit or minimum cost, budgeting, etc.

## Consider the problem given below

A car company produces 2 models, model A and model B. Long-term projections indicate an expected demand of at least 100 model A cars and 80 model B cars each day. Because of limitations on production capacity, no more than 200 model A cars and 170 model B cars can be made daily. To satisfy a shipping contract, a total of at least 200 cars much be shipped each day. If each model A car sold results in a \$2000 loss, but each model B car produces a \$5000 profit, how many of each type should be made daily to maximize net profits?

In this problem, we have to maximize profits given the prices of both types of cars. Let *x* be the number of *A* cars and *y* be the number of *B* cars.

The objective function is:  
<p style="text-align: center;"> **Z = max(-2000x + 5000y)** </p>

Given constraints:  
<p style="text-align: center;"> x + y >= 200    
                                100 <= x <= 200   
                                80 <= y <= 170 </p>
                    
                    
For solving this problem with **R**, we have to use the `lpsolve` library.

```{r include = FALSE}
library(lpSolveAPI)
```
<br>

#### 1. Initializing the Linear Program

Now, we initialize and define the architecture of the problem. 
``` {r include = FALSE}
p = make.lp(0, 2)
```

Here, the parameters passed to `make.lp()` are the number of rows and columns of coefficient matrix. Number of rows represent the number of constraints and number of columns represent the number of variables. Since, we have two variables *x* and *y* in this problem, we set the number of columns to 2. We'll add the constraints later.  
<br>

#### 2. Defining the Objective Function

The objective function is declared using the functions `set.objfc()` and `lpcontrol()`. In both functions, the first argument is the name of the linear optimization interface we are working (in this case, `p`.) `lp.control()` uses the argument `sense`, that can be set to "max" or "min" to reflect that we want to maximize or minimize the objective function, respectively.

`set.objfn()` declares the coefficients of each of the variables in the objective function to be optimized, and takes a vector as an argument into `obj`.

```{r include = FALSE}
lp.control(p, sense = "max")
set.objfn(p, obj = c(-2000, 5000))
```
<br>

#### 3. Defining the Constraints

We can now add *constraining equations*. We can enable ad-hoc addition of constraints to the linear optimization interface using the `row.add.mode()` function. Once enabled, we can now use the `add.constraint()` function to add the constraints. The arguments for `add.constraint()` are:

-> `xt`: The **non-zero coefficients** of the constraining equation. **All coefficients that are not explicitly declared are assigned to zero.**  
-> `type`: The **equality or inequality function** to be compared, as a character string (in double quotes.)  
-> `rhs`: The value of **the right-hand side of the inequality**, which *must be a constant*. Any variables that exist in the constraint must be moved over to the left side of this equation, if they exist.  
-> `indices`: The **indices of the values established by `xt`**, which need not be consecutive but must fall within the range of number of variables.

```{r include = FALSE}
row.add.mode(p, "on")
add.constraint(p, xt = c(1, 1), type = ">=", rhs = 200, indices = c(1, 2))
row.add.mode(p, "off")
```

We can set `row.add.mode` to `"off"` after adding all the constraints. The solver will not work if it is still `"on"`.  
<br>

#### 4. Defining Boundaries

We now add the boundaries. The upper and lower bounds are set using the `set.bounds()` function, in the arguments `lower` and `upper` respectively. As in the objective function, the arguments **must be vectors that contain a value for each variable in the optimization problem.**

``` {r include = FALSE}
set.bounds(p, lower = c(100, 80), upper = c(200, 170))
```
<br>

#### 5. Solving and Viewing Results

We use the `solve()` function to solve the system of equations.
To view the results, we can extract the output of the optimized objective function (in this case, our maximized profit) and the values of each of the variables under optmized conditions using the `get.objective()` and `get.variables()` functions, respectively. The `solve` function returns **0** if a solution is found.


``` {r}
p
solve(p)
get.variables(p)
get.objective(p)
```

Here, the maximum profit is **`r format(get.objective(p), scientific = FALSE)`$** which can be achieved by selling **`r get.variables(p)[1]` A cars** and **`r get.variables(p)[2]` B cars**.
<br><br>