Decision Support Systems

University of Central Florida, Dr. Garibay

Lecture - Linear Programming

Notebook Learning Objectives

After studying this notebook students should be able to:

Visualize a linear relationships

Solve a linear programming problem using Pulp

Overview
Linear programming is a method to achieve the best outcome in a mathematical model whose requirements are represented by linear
relationships. More formally, linear programming is a technique for the optimization of a linear objective function, subject to linear
equality and linear inequality constraints. Its feasible region is a convex polytope, which is a set defined as the intersection of finitely
many half spaces, each of which is defined by a linear inequality. Its objective function is a real-valued affine (linear) function defined
on this polyhedron. A linear programming algorithm finds a point in the polyhedron where this function has the smallest (or largest)
value if such a point exists.

The problem of solving a system of linear inequalities dates back to at least Fourier, who in 1827 published a method for solving it.
Later on John von Neumann imporved the Simplex Method of solving them to use it on game theory. More recently, it was proven that
a system of linear inequalities can be solved in polinomial time and with the introduction in 1984 of the interior-point method there is
now efficient ways to solving linear inequalities numerically.

Linear programming has wide applications in the fields of economics, operations research, and many others. In particular many
problems in operations reserach can be expressed as linear programs.

Standard Form
We can describe a linear programming problems in three parts:

A linear function to be optimized, for example:

f (x1, x2) = c1x1 + c2x2

Problem constraints of the following form, for example:

a11x1 + a12x2 ≤ b1

a21x1 + a22x2 ≤ b2

a31x1 + a32x2 ≤ b3

Non-negative decision variables_, for example:

x1 ≥ 0

x2 ≥ 0

In general, we can express this problem in matrix notation:

T
max{c x | Ax ≤ b ∧ x ≥ 0}

Where, x is a vector of variables to be solved, c and b are vectors of known coefficients, A is a known matrix of coefficients, and ()T is
the matrix transpose operator.

Solving Linear Programs (LP) using PuLP

PuLP is a free open source software written in Python. It is used to describe optimisation problems as mathematical models.

Installing PuLP
To install PuLP in your computer follow this instructions:

Open a terminal window

Execute the following command: conda install -c conda-forge pulp

Now lets load the libraries needed for this notebook, including pulp

In [27]:
import numpy as np
import matplotlib.pyplot as plt

import pulp

A Simple Example
Let’s start using PuPL with a simple example. Lets maximize the following objective function Z :

Z = 4x + 3y

with x and y being our decision variables.

In addition, the objective function is subject to the following constraints:

x ≥ 0

y ≥ 2

2y ≤ 25–x

4y ≥ 2x–8

y ≤ 2x − 5

Before we use PuLP to sove this Linear Program, lets first graph these contrains and the objective function to get a good idea of how
the function we are trying to optimize looks like and how the contrains limit this function.

Visualizing Constrains and Objective

Plotting the constrain boundaries
Lets use standard plotting tools to plot the constrains:

In [28]:
# Construct lines

# x > 0

x = np.linspace(0, 20, 2000) # start=0, stop=20, samples=2,000

# y >= 2

y1 = x*0+2 # we need the "x" so that "y1" also become a vector!

# To graph we just calculate the equality (=)

# at this step, later we will "shade" the area

# included in the inequality (>=) later.

# 2y <= 25 - x

y2 = (25-x)/2.0

# 4y >= 2x - 8

y3 = (2*x-8)/4.0

# y <= 2x - 5
y4 = 2 * x -5

# Make plot
plt.plot(x, y1, label=r'$y\geq2$') #for label in math text use

plt.plot(x, y2, label=r'$y\leq (25-x)/2$') # r'$ $'

plt.plot(x, y3, label=r'$y\geq (2x - 8)/4$')

plt.plot(x, y4, label=r'$y\leq 2x-5$')

plt.xlim((0, 18)) #limits for plot x-axis, left=0, right=18

plt.ylim((0, 11)) #limits for plot y-axis, left=0, right=11

plt.xlabel(r'$x$')

plt.ylabel(r'$y$')

plt.legend(bbox_to_anchor=(1.0, 1.0) ,loc='upper right' )

<matplotlib.legend.Legend at 0x12f1a86d0>
Out[28]:

Plotting the feasible space using .fill_between(x,y1,y2,where)

Fill the area between two horizontal curves

Parameters:
x : array (length N). The x coordinates of the nodes defining the curves.

y1 : array (length N) or scalar. The y coordinates of the nodes defining the first curve.

y2 : array (length N) or scalar, optional, default: 0. The y coordinates of the nodes defining the second curve.

where : array of bool (length N), optional, default: None

Define where to exclude some horizontal regions from being filled. The
filled regions are defined by the coordinates x[where]. More precisely, fill between x[i] and x[i+1] if where[i] and where[i+1]. Note
that this definition implies that an isolated True value between two False values in where will not result in filling. Both sides of the
True position remain unfilled due to the adjacent False values.

In [29]:
# Make plot
plt.plot(x, y1, label=r'$y\geq2$') #for label in math text use

plt.plot(x, y2, label=r'$2y\leq25-x$') # r'$ $'

plt.plot(x, y3, label=r'$4y\geq 2x - 8$')

plt.plot(x, y4, label=r'$y\leq 2x-5$')

plt.xlim((0, 18)) #limits for plot x-axis, left=0, right=18

plt.ylim((0, 11)) #limits for plot y-axis, left=0, right=11

plt.xlabel(r'$x$')

plt.ylabel(r'$y$')

# Fill feasible region

# minimum of

# y <= (25 - x)/2

# y <= 2x - 5

y5 = np.minimum(y2, y4) #Compare two arrays and returns a new

#array containing the element-wise minima

# maximum of

# y >= 2

# y >= (2x - 8)/4

y6 = np.maximum(y1, y3) #Same as above, but return maxima

plt.fill_between(x, y5, y6, where=y5>y6, color='grey', alpha=0.5)

plt.legend(bbox_to_anchor=(1.0, 1.0) ,loc='upper right' )

<matplotlib.legend.Legend at 0x12f2b7590>
Out[29]:

Plotting the feasible space plus the heat map for the objective function
In [30]:
# Make plot
plt.plot(x, y1, label=r'$y\geq2$') #for label in math text use

plt.plot(x, y2, label=r'$2y\leq25-x$') # r'$ $'

plt.plot(x, y3, label=r'$4y\geq 2x - 8$')

plt.plot(x, y4, label=r'$y\leq 2x-5$')

plt.xlim((0, 18)) #limits for plot x-axis, left=0, right=18

plt.ylim((0, 11)) #limits for plot y-axis, left=0, right=11

plt.xlabel(r'$x$')

plt.ylabel(r'$y$')

# Z=4x+3y

def f(x, y):

return 4*x + 3*y

x = np.linspace(0, 20, 2000)

y = np.linspace(0, 20, 2000)

X, Y = np.meshgrid(x, y)

Z = f(X, Y)

plt.contourf(X, Y, Z, 30, cmap = 'coolwarm')

plt.colorbar();

Solving Simple Example using PuLP

Step 1: create linear optimization problem
The first step is to crete a variable for our problem. In our case we are using my_lp_problem , but the name could be anything. We
create this variable using pulp.LpProblem() . This function takes two arguments. The first is the name of our problem, "My LP
Problem" in our example below. The second one can take the values pulp.LpMaximize or pulp.LpMinimize depending if
our optimization problem is a maximiation or a minimization problem.

In [31]:
# Create the 'prob' variable to contain the problem data
my_lp_problem = pulp.LpProblem("My LP Problem", pulp.LpMaximize)

Step 2: setting the decision variables

Optimization problems have variables that we would like to optimize, to find their values that either maximize or minimize the objective
function. These are called the decision variables. The optimal values of these decision variables are the output of a Pulp optimization.
To define these variables we use pulp.LpVariable() . This function takes four arguments. The first one is the name of the
decision variabled, i.e., 'x' . The second and third ones determine the lower bound and upper bounds for the optimization (the lowest
and highest posible values for this variable that we will considered while optimizing), i.e., lowBound=0, upBound=None (explore
from 0 to positive infinity). Finally, the third argument determines if the variable is going to take continuous values or discrete values
using Continuous or Integer respectively, i.e., cat='Continuous' .

In [32]:
x = pulp.LpVariable('x', lowBound=0, upBound=None, cat='Continuous')
y = pulp.LpVariable('y', lowBound=2, upBound=None, cat='Continuous')

Step 3: Input objective function and constrains

Now we start inputing information to our problem. The first function to input is the objective function, the function we are maximizing or
minimizing. We use the operator += to do this. The line of code we use to input the objective function should end with , "Name of
objective function" to indicate that is the objective function and provide a name for it.

In [33]:
# Objective function, first one to input, put lable at end "Z"
my_lp_problem += 4 * x + 3 * y, "Z"

Next, we input all the constrains we may have in our optimization problem. We use += as above follow by a relationship. Usually we
use the following relationship operators: ==, , >=, <= as shown below:

In [34]:
# Constraints
#my_lp_problem += y >= 2 # no needed since enter in variable as lowBound

my_lp_problem += 2 * y <= 25 - x

my_lp_problem += 4 * y >= 2 * x - 8

my_lp_problem += y <= 2 * x - 5

#my_lp_problem += y == 0 # this would make in "infeasible"

Step 4: Solve our optimization using .solve

Once we have setup our problem, input our objective function and input our constrains, finding the solution of our problem requires a
single line of code:

In [35]:
my_lp_problem.solve()

1
Out[35]:

Of course, behind the scenes, Pulp and Python are executing numerical optimization algorithms to find a solution. PuLP choses the
solver to use. We can determine what solver PuPL should use by providing an argument: my_lp_problem.solve(CPLEX()) . We
can see what solvers are available in our installation using pulp.pulpTestAll() .

Step 5: Display Solution

Now we can display the results from the solver as our output. First we ask the solver for he status of the solution. Possible status are:

Not Solved
Infeasible
Unbounded
Undefined
Optimal

In [36]:
print("Status:", pulp.LpStatus[my_lp_problem.status])

('Status:', 'Optimal')

Now, we can print the variables and the solved optimum values as well as the optimal value for the objective function. We use
.name() .varValue() and .objective() as follows:

In [37]:
print x.name, x.varValue
print y.name, y.varValue

print pulp.value(my_lp_problem.objective)

x 14.5

y 5.25

73.75

Example : Car Model Profit Maximization

Lets assume that a company manufactures two models of cars. Lets call these model A and model B. When sold, each of these cars
produce a profit of 30,000 and 45,000 for models A and B respectively. We want to maximaze profit described by the following equation
(our objective function):

prof itmodel = 30, 000A + 45, 000B

But of course, we can not produce infinite amount of cars. The car manufacturing is restricted by factory capacity, supplier capacity,
financial resources available, etc. Furthermore, the deamand for cars is not infinite and the particular market segment this company
operates on could have several competitors. For this example, lets simply assume that the following equations represent the
manufacturing and market constraints:

3A + 4B <= 30000

5A + 6B <= 60000

2A + 3B <= 21000

To solve this LP, we follow the steps described above:

In [38]:
# Instantiate our problem class
model = pulp.LpProblem("Profit maximising problem", pulp.LpMaximize)

# create the decision variables

A = pulp.LpVariable('A', lowBound=0, cat='Integer')

B = pulp.LpVariable('B', lowBound=0, cat='Integer')

# Input objective function

model += 30000 * A + 45000 * B, "Profit"

# Input constraints

model += 3 * A + 4 * B <= 30000

model += 5 * A + 6 * B <= 60000

model += A + 3 * B <= 21000

# Solve our problem

model.solve()

# check status of solution

print pulp.LpStatus[model.status]

# Print our decision variable values

print "Production of Car A = ", A.varValue

print "Production of Car B = ", B.varValue

# Print our objective function value

print "Maximum Profit = ", pulp.value(model.objective)

Optimal

Production of Car A = 1200.0

Production of Car B = 6600.0

Maximum Profit = 333000000.0

Exercises
1. An operations manager is trying to determine a production plan for the next week. There are three products (say, P, Q, and Q) to
produce using four machines (say, A and B, C, and D). Each of the four machines performs a unique process. There is one machine of
each type, and each machine is available for 2400 minutes per week. The unit processing times for each machine is given in Table 1.

Table 1. Unit Processing Time (min)

Machine Product P Product Q Product R Availability (min)

A 20 10 10 2400

B 20 28 16 2400

C 20 6 16 2400

D 20 15 0 2400

Total processing time 57 59 42 9600

The unit revenues and maximum sales for the week are indicated in Table 2. Storage from one week to the next is not permitted. The
operating expenses associated with the plant are $6000 per week, regardless of how many components and products are made. The
$6000 includes all expenses except for material costs.

Table 2. Product Data

Item Product P Product Q Product R

Revenue per unit 20 10 10

Material cost per unit 20 28 16

Profit per unit 20 6 16

Maximum sales 20 15 0

Solution: We formulate the problem as follows:

Defining decision variables:

p , number of units of product P to produce,

q , number of units of product Q to produce,

r , number of units of product R to produce,

choosing an objective function: our objective is to maximize profit:

P rof it = (90 − 45)p + (100 − 40)q + (70 − 20 ) r –6 000 = 45p + 60q + 50r–6000

indentifying the constrains: the amount of time a machine is available and the maximum sales posible restrict the quantities to
be manufactured. Using the unit processing times for each machine, the linear constraints can be written as:

20p + 10q + 10r ≤ 2400 (Machine A)

12p + 28q + 16r ≤ 2400

15p + 6q + 16r ≤ 2400

10p + 15q + 0r ≤ 2400

Market limitations can be written as upper bounds:

p ≤ 100

q ≤ 40

r ≤ 60

Number of products to produce can not be negative:

p ≥ 0

q ≥ 0

r ≥ 0

In [39]:
# Instantiate our problem class
model2 = pulp.LpProblem("Manufacturing problem", pulp.LpMaximize)

# create the decision variables

p = pulp.LpVariable('p', lowBound=0, cat='Continuous')

q = pulp.LpVariable('q', lowBound=0, cat='Continuous')

r = pulp.LpVariable('r', lowBound=0, cat='Continuous')

# Input objective function

model2 += 45*p + 60*q + 50*r-6000, "Profit"

# Input constraints

model2 += 20*p + 10*q + 10*r <= 2400

model2 += 12*p + 28*q + 16*r <= 2400

model2 += 15*p + 6*q + 16*r <= 2400

model2 += 10*p + 15*q + 0*r <= 2400

model2 += p <= 100

model2 += q <= 40

model2 += r <= 60

# Solve our problem

model2.solve()

# check status of solution

print pulp.LpStatus[model2.status]

# Print our decision variable values

print "Production of product p = ", p.varValue

print "Production of product q = ", q.varValue
print "Production of product r = ", r.varValue

# Print our objective function value

print "Maximum Profit = ", pulp.value(model2.objective)

Optimal

Production of product p = 81.818182

Production of product q = 16.363636

Production of product r = 60.0

Maximum Profit = 1663.63635

2. A post office requires different numbers of full- time employees on different days of the week. Each full-time employee must work five
consecutive days and then receive two days off.
In the following table, the number of employees required on each day of the week is
specified.

Day Number of full-time Employees Required

1=Monday 17

2=Tuesday 13

3=Wednesday 15

4=Thursday 19

5=Friday 14

6=Saturday 16

7=Sunday 11