Introduction to Linear Optimization in R

Introduction

Linear programming is one of the most extensively used techniques in the toolbox of quantitative methods of optimization. Its origins date as early as 1937, when Leonid Kantorovich published his paper A new method of solving some classes of extremal problems. Kantorovich developed linear programming as a technique for planning expenditures and returns in order to optimize costs to the army and increase losses to the enemy.

The method was kept secret until 1947, when George B. Dantzig published the simplex method for solving linear programming

Roughly speaking, the linear programming problem consists in optimizing (that is, either minimize or maximize) the value of a linear objective function of a vector of decision variables, considering that the variables can only take the values defined by a set of linear constraints. Linear programming is a case of mathematical programming, where objective function and constraints are linear.

A formulation of a linear program in its canonical form of maximum is:

\[max Z = \sum_{i = 1} ^ n c_i x_i\] Subject to

\[\sum_{i = i} ^ n \sum_{j = i} ^ m a_{ij}x_i \le b_j\] \[\sum_{i = i} ^ n \sum_{j = i} ^ m a_{ij}x_i \ge b_k\]

\[\sum_{i = i} ^ n \sum_{j = i} ^ m a_{ij}x_i = b_l\]

The model has the following elements:

1. An objective function of the n decision variables \(x_i\) . Decision variables are affected by the cost coefficients \(c_j\)

2. A set of m constraints, in which a linear combination of the variables affected by coefficients \(a_{ij}\) has to be less or equal than its right hand side value \(b_i\) (constraints with signs greater or equal or equalities are also possible)

3. The bounds of the decision variables. In this case, all decision variables have to be nonnegative

The LP formulation shown above can be expressed in matrix form as follows:

\[max Z = \mathbf{c}^T \mathbf{x}\]

Subject to

\[\mathbf{Ax = b}\]

\[\mathbf{x} \ge 0\]

Modeling Linear Programs

Let’s consider the following situation:

A small business sells two products, named Product 1 and Product 2. Each tonne of Product 1 consumes 30 working hours, and each tonne of Product 2 consumes 20 working hours. The business has a maximum of 2700 working hours for the period considered. As for machine hours, each tonne of Products 1 and 2 consumes 5 and 10 machine hours, respectively. There are 850 machine hours available.

Each tonne of Product 1 yields 20 Me of profit, while Product 2 yields 60 Me for each tonne sold. For technical reasons, the firm must produce a minimum of 95 tonnes in total between both products. We need to know how many tonnes of Product 1 and 2 must be produced to maximize total profit.

This situation is apt to be modeled as a PL model. First, we need to define the decision variables. In this case we have:

1. \(x_1\) number of tonnes produced and sold of Product 1
2. \(x_2\) number of tonnes produced and sold of Product 2

The cost coefficients of these variables are 20 and 60, respectively. Therefore, the objective function is defined multiplying each variable by its corresponding cost coefficient

\[maxZ = 20 x_1 + 60 x_2\]

A constraint WH making that the total amount of working hours used in Product 1 and Product 2, which equals \(30 x_1 + 20 x_2\), is less or equal than 2,700 hours.

A similar constraint MH making that the total machine hours \(5 x_1 + 10 x_2\) are less or equal than 850.

A constraint making that the total units produced and sold \(x_1 + x_2\) are greater or equal than 95

Putting all this together, and considering that the decision variables are nonnegative, the LP that maximizes profit is:

\[maxZ = 20 x_1 + 60 x_2\]

Subject to

\[30 x_1 + 20 x_2 \le 2700\]

\[5 x_1 + 10 x_2 \le 850\] \[x_1 + x_2 \le 95 \] \[x_1, x_2 \ge 0\]

The following code solves that LP with two variables:

prod.sol$objval

[1] 4900

A Trasportation Problem

Let’s consider a transportation problem of two origins a and b, and three destinations 1, 2 and 3. In the Table below are presented the cost \(c_{ij}\) of transporting one unit from the origin \(i\) to destination \(j\), and the maximal capacity of the origins and the required demand in the destinations.

We need to know how we must cover the demand of the destinations at a minimal cost.

Table2 = data.frame(Parameter = c('Origin A', 'Origin B', 'Demand'),
                    Dest.1 = c(8, 2, 40),
                    Dest.2 = c(6, 4, 35),
                    Dest.3 = c(3, 9, 25),
                    Capacity = c(70, 40, NA))
Table2

NA

This situation can be modeled with a LP with the following elements:

1. Decision variables of the form $x_{ij} , representing units transportedfrom origin \(i\) to destination \(j\)
2. An objective function with cost coefficients equal to \(c_{ij}\)
3. Two sets of constraints: a less or equal set of constraints for each origin, limiting the units to be transported, and a greater of equal set of constraints representing that the demand of each destination must be covered.

The resulting LP is:

\[maxZ = 8x_{a1} + 6x_{a2} + 3x_{a3} + 2x_{b1} + 4x_{b2} + 9x_{b3}\]

Subject to:

\[x_{a1} + x_{a2} + x_{a3} \le 70\] \[x_{b1} + x_{b2} + x_{b3} \le 40\] \[x_{a1} + x_{b1} \ge 40\] \[x_{a2} + x_{b2} \ge 35\] \[x_{a3} + x_{b3} \ge 25\]

\[x_{ij} \ge 0\]

sol

     [,1] [,2] [,3]
[1,]    0   35   25
[2,]   40    0    0

Transportation by Truks

In Table3 can be found the quarterly demand (in tonnes) and the acquisition costs per tonne (in ke per tonne) for each quarter of raw materials for a chemical plant. All purchases in a given quarter can be used to cover the demand of the present quarter, or the demand of quarters in the future. The costs of stocking are of 8k USD per tonne stored at the end of each quarter. The stocks at the beginning of first quarter are of 100 tonnes, and it is needed the same amount of stock at the end of the fourth quarter.

In addition to the purchase and storage costs, the transportation costs have to be considered. All the purchased quantity of raw materials has to be transported, using any combination of the two available truck models:

1. Small trucks: cost of 700k uSD, and capacity of 500 tonnes.
2. Large trucks: cost of 1,400k USD, and capacity of 1,200 tonnes.

We need to define a linear programming model that allows the minimization of the total costs: acquisition, storage and transport, obtaining the amount raw materials to purchase, and the amount of trucks of both kinds to be contracted each quarter.

Table3 = data.frame(Quarter = c('Demand', 'Unit Cost'),
                    T1 = c(1000, 20) ,
                    T2 = c(1200, 25),
                    T3 = c(1500, 30),
                    T4 = c(1800, 40))

Table3

The variables to define are:

1. \(q_i\) continuous: tonnes of raw material to purchase in quarter i
2. \(s_i\) continuous: tonnes in stock at the end of quarter i, and \(s_0\) as the initial stock
3. \(t_i\) integer: small trucks to contract in quarter i
4. \(u_i\) integer: large trucks to contract in quarter i

Once defined the variables, two sets of constraints have to be defined:

1. Constraints assuring that the purchase plan meets the demand commitments. These are of the form \(s_{i-1} + q_i - s_i = d_i\), being \(d_i\) the demand of the quarter.
2. Constraints assuring that a sufficient number of each kind of trucks is contracted: \(q_i \le 500t_i + 12004 u_i\)

The resulting model is:

\[minZ = \sum_{i = 1} ^ 4 (c_iq_i + 8s_i + 700t_i + 1400u_i) \] Subject to:

\[s_{i-1} + q_i - s_i = d_i \Rightarrow\forall i = 1,...,4\] \[q_i - 500t_i - 1200ui \le 0 \Rightarrow\forall i = 1,...,4\]

\[s_0 = s_4 = 100\]

\[s_i, q_i \ge 0\]

prod.trans$solution

 [1]  900.000000 1200.000000 3400.000000    0.000000  100.000000    0.000000    0.000000 1900.000000  100.000000    0.000000    0.000000    0.000000
[13]    0.000000    0.750000    1.000000    2.833333    0.000000

In this case, the solution has some practical and technical issues. We supposed to have only integer solutions and in the case of \(u_1 = 0.75\) and \(u_3 = 2.83\) this condition is not true.

To fix this, we use Integer Programming models ensuring all variables to be integer:

Hiring & Firing

In Table4 are listed the needs of pilots able to flight an A320 for the following six months. The cost of a pilot’s salary is $8k USD per month. At the beginning of Month 1 the airline has a staff of 20 pilots, but this staff can be adjusted each month.

Pilots can be hired and fired at the beginning of each month. Newly hired pilots can start working at the same month, and fired pilots stop working the same day they are fired. The cost of firing a pilot is $10k USD, and the hiring cost is of 5K USD per pilot. If it is convenient, the airline can have a staff of pilots larger than the actual needs.

1. Define a linear programming model to obtain the pilots to hire and fire each month to minimize the total cost of pilot staff (costs of salary plus hiring and firing costs).

To model this situation, we’ll have to define the following variables:

1. Variables \(h_i\): pilots hired at the beginning of month i
2. Variables \(f_i\): pilots fired at the beginning of month i
3. Variables \(s_i\): staff of pilots during month i

The model should have the following groups of constraints:

1. Constraints assuring that the staff of pilots at the beginning of month i is equal to \(s_i = h_i - f_i + s_{i-1}\). In this case, we have that \(s_0 = 20\).
2. Constraints assuring that variables \(s_i\) are bigger of equal to the values of staff required di listed in Table4.

table4 = data.frame(Months = c('m1', 'm2', 'm3', 'm4', 'm5', 'm6'),
                    PilotsNeed = c(30, 60, 55, 40, 45, 50))
table4

The linear problem to be solced is:

\[minZ = 5\sum_{i = 1} ^ 6 h_i + 10\sum_{i = 1} ^ 6 f_i + 8\sum_{i = 1} ^ 6 s_i\]

Subject to

\[s_i = h_i - f_i + s_{i-1} \Rightarrow\forall i = 1,...,6\] \[s_i \ge d_i \Rightarrow\forall i = 1,...,6\] \[h_i, f_i, s_i \ge 0\]

#Variable names
var.names = c('h1', 'h2', 'h3', 'h4', 'h5', 'h6',
                  'f1', 'f2', 'f3', 'f4', 'f5', 'f6',
                  's0', 's1', 's2', 's3', 's4', 's5', 's6')


#destinations run j in 1:n
obj.fun <- c(rep(5, 6), rep(10, 6), 0, rep(8, 6))
names(obj.fun) = c('h1', 'h2', 'h3', 'h4', 'h5', 'h6',
                  'f1', 'f2', 'f3', 'f4', 'f5', 'f6',
                  's0', 's1', 's2', 's3', 's4', 's5', 's6')

m <- length(obj.fun) #number of vars
n <- 7 #number of months + 1

constr1 <- matrix (0, nrow = n, ncol = m)
colnames(constr1) = var.names

#Create model for Constraint 1: s[i-1] - s[i] - f[i] + h[i] = 0

constr1[1 : 6, 1 : 6] = diag(x = 1, ncol = 6, nrow = 6)
constr1[1 : 6, 7 : 12] = diag(x = -1, ncol = 6, nrow = 6)
constr1[1 : 6, 13 : 18]  = diag(x = 1, ncol = 6, nrow = 6)
i= 1
j = 14
for(i in 1 : (n - 1)){
    constr1[i, j] = -1
    j = j + 1
}
constr1[7, 13]  = 1

#Create second constraint on minimun requirements s[i] > = d[i]
constr2 <- matrix (0, nrow = n - 1, ncol = m)
constr2[1 : 6, 14 : 19] = diag(x = 1, ncol = 6, nrow = 6)

#Bind Constr1 and Consr2 in a single matrix
constr = rbind(constr1, constr2)

#Create RHS consraintin a vector of 1 row x 13 rows
rhs = c(rep(0, 6), 20, 30, 60, 55, 40, 45, 50)

#Create direction
constr.dir <- c(rep("=", 7), rep(">=", 6))

#Solve Problem
staffing = lp('min', obj.fun, constr, constr.dir, rhs)

#Solution
staffing = staffing$solution
staffing = data.frame(Var.Name = var.names, Solution = staffing)
staffing

NA

2. Modify the linear model to include the constraint that the airline cannot fire pilots if it has hired pilots the previous month.

Looking at the solution of the previous problem, it can be seen that this new constraint does not hold for months 2 and 3: in month 2 are hired 30 pilots, and in month 3 are fired 5 pilots. Then a new model has to be defined to account for this new restriction. To do so, we have to add a new binari variable:

• Variabe \(b_i\): equals one if pilots are hired in month i, and zero otherwise

Then, two new sets of constraints must be added: one set assuring that \(b_i = 0 \Rightarrow f_i = 0\), and another set making that \(b_i = 1 \Rightarrow f_{i+1} = 0\)

The linear problem to be solced is:

\[minZ = 5\sum_{i = 1} ^ 6 h_i + 10\sum_{i = 1} ^ 6 f_i + 8\sum_{i = 1} ^ 6 s_i\]

Subject to

\[s_i = h_i - f_i + s_{i-1} \Rightarrow\forall i = 1,...,6\] \[s_i \ge d_i \Rightarrow\forall i = 1,...,6\] \[f_i \le Mb_i \Rightarrow \forall i = 1,...,5\] \[h_{i+1} \le M(1 - b_i) \Rightarrow \forall i = 1,...,5\]

\[b_i \in [0, 1]\] \[h_i, f_i, s_i \ge 0\]

The letter \(M\) represents an artifitial penalty and is usually asumed to be a big number. In the context of the problem, \(M=1000\) is big.

#Variable names

M = 1000 #Penalty Variable

var.names = c('h1', 'h2', 'h3', 'h4', 'h5', 'h6',
                  'f1', 'f2', 'f3', 'f4', 'f5', 'f6',
                  's0', 's1', 's2', 's3', 's4', 's5', 's6',
              'b1', 'b2', 'b3', 'b4', 'b5')


#destinations run j in 1:n
obj.fun <- c(rep(5, 6), rep(10, 6), 0, rep(8, 6), rep(0, 5))
names(obj.fun) = var.names

m <- length(obj.fun) #number of vars
n <- 7 #number of months + 1

#Create model for Constraint 1: s[i-1] - s[i] - f[i] + h[i] = 0
constr1 <- matrix (0, nrow = n, ncol = m)
colnames(constr1) = var.names


constr1[1 : 6, 1 : 6] = diag(x = 1, ncol = 6, nrow = 6)
constr1[1 : 6, 7 : 12] = diag(x = -1, ncol = 6, nrow = 6)
constr1[1 : 6, 13 : 18]  = diag(x = 1, ncol = 6, nrow = 6)
i= 1
j = 14
for(i in 1 : (n - 1)){
    constr1[i, j] = -1
    j = j + 1
}
constr1[7, 13]  = 1

#Create second constraint on minimun requirements s[i] > = d[i]
constr2 <- matrix (0, nrow = n - 1, ncol = m)
constr2[1 : 6, 14 : 19] = diag(x = 1, ncol = 6, nrow = 6)

#Create the 3rd set of constraints f[i] - M b[i] <= 0
constr3 = matrix(0, nrow = 5, ncol = m)
colnames(constr3) = var.names

constr3[, 7 : 11] = diag(1, nrow = 5, ncol = 5)
constr3[, 20 : 24] = diag(-M, nrow = 5, ncol = 5)

#Create the 4th set of constraints h[i+1] + M b[i] <= M
constr4 = matrix(0, nrow = 5, ncol = m)
colnames(constr4) = var.names

constr4[, 2 : 6] = diag(1, nrow = 5, ncol = 5)
constr4[, 20 : 24] = diag(M, nrow = 5, ncol = 5)

#Bind Constr1 and Consr2 in a single matrix
constr = rbind(constr1, constr2, constr3, constr4)

#Create RHS of ecuations
#Create RHS consraintin a vector of 1 row x 13 rows
rhs = c(rep(0, 6), 20, 30, 60, 55, 40, 45, 50, rep(0, 5), rep(M, 5))

#Create direction
constr.dir <- c(rep("=", 7), rep(">=", 6), rep('<=', 5), rep('<=', 5))

#Solve Problem
staffing = lp('min', obj.fun, constr, constr.dir, rhs, binary.vec = c(20 : 24))

#Solution
staffing = staffing$solution
staffing = data.frame(Var.Name = var.names, Solution = staffing)
staffing

NA
NA
NA