Linear Programming with OR-Tools

In this notebook, we do some basic LP solving with Google's OR-Tools. Problems used will be examples in Hamdy Taha's Operations Research: An Introduction, 9th Edition, which I have in paperback.



In [13]:

    
from ortools.linear_solver import pywraplp

Reddy Mikks model

Given the following variables:

$\begin{aligned} x_1 = \textrm{Tons of exterior paint produced daily} \newline x_2 = \textrm{Tons of interior paint produced daily} \end{aligned}$

and knowing that we want to maximize the profit, where \$5000 is the profit from exterior paint and \$4000 is the profit from a ton of interior paint, the Reddy Mikks model is:

$$\textrm{Maximize } z = 5x_1 + 4x_2$$

subject to $$6x_1 + 4x_2 \le 24$$ $$x_1 + 2x_2 \le 6$$ $$-x_1 + x_2 \le 1$$ $$x_2 \le 2$$ $$x_1, x_2 \ge 0$$



In [14]:

    
reddymikks = pywraplp.Solver('Reddy_Mikks', pywraplp.Solver.GLOP_LINEAR_PROGRAMMING)

x1 = reddymikks.NumVar(0, reddymikks.infinity(), 'x1')
x2 = reddymikks.NumVar(0, reddymikks.infinity(), 'x2')

reddymikks.Add(6*x1 + 4*x2 <= 24)
reddymikks.Add(x1 + 2*x2 <= 6)
reddymikks.Add(-x1 + x2 <= 1)
reddymikks.Add(x2 <= 2)

profit = reddymikks.Objective()
profit.SetCoefficient(x1, 5)
profit.SetCoefficient(x2, 4)
profit.SetMaximization()

status = reddymikks.Solve()

if status not in [reddymikks.OPTIMAL, reddymikks.FEASIBLE]:
    raise Exception('No feasible solution found')
    
print("The company should produce",round(x1.solution_value(),2),"tons of exterior paint")
print("The company should produce",round(x2.solution_value(),2),"tons of interior paint")
print("The optimal profit is", profit.Value(), 'thousand USD')









    



The company should produce 3.0 tons of exterior paint
The company should produce 1.5 tons of interior paint
The optimal profit is 21.0 thousand USD

More simple problems

A company that operates 10 hours a day manufactures two products on three sequential processes. The following data characterizes the problem:



In [11]:

    
import pandas as pd

problemdata = pd.DataFrame({'Process 1': [10, 5], 'Process 2':[6, 20], 'Process 3':[8, 10], 'Unit profit':[20, 30]})
problemdata.index = ['Product 1', 'Product 2']

problemdata









    Out[11]:







  
    
      
      Process 1
      Process 2
      Process 3
      Unit profit
    
  
  
    
      Product 1
      10
      6
      8
      20
    
    
      Product 2
      5
      20
      10
      30

Where there are 10 hours a day dedicated to production. Process times are given in minutes per unit while profit is given in USD.

The optimal mix of the two products would be characterized by the following model:

$\begin{aligned} x_1 = \textrm{Units of product 1} \newline x_2 = \textrm{Units of product 2} \end{aligned}$

$$\textrm{Maximize } z = 20x_1 + 30x_2$$

subject to

$$\begin{array}{rcl} 10x_1 + 5x_2 \le 600 \newline 6x_1 + 20x_2 \le 600 \newline 8x_1 + 10x_2 \le 600 \newline x_1, x_2 \ge 0 \end{array}$$

(we will assume that continuous solution values are acceptable for this problem)



In [30]:

    
simpleprod = pywraplp.Solver('Simple_Production', pywraplp.Solver.GLOP_LINEAR_PROGRAMMING)

x1 = simpleprod.NumVar(0, simpleprod.infinity(), 'x1')
x2 = simpleprod.NumVar(0, simpleprod.infinity(), 'x2')

for i in problemdata.columns[:-1]:
    simpleprod.Add(problemdata.loc[problemdata.index[0], i]*x1 + problemdata.loc[problemdata.index[1], i]*x2 <= 600)

profit = simpleprod.Objective()
profit.SetCoefficient(x1, 20)
profit.SetCoefficient(x2, 30)
profit.SetMaximization()

status = simpleprod.Solve()

if status not in [simpleprod.OPTIMAL, simpleprod.FEASIBLE]:
    raise Exception('No feasible solution found')
    
print("The company should produce",round(x1.solution_value(),2),"units of product 1")
print("The company should produce",round(x2.solution_value(),2),"units of product 2")
print("The optimal profit is", round(profit.Value(),2), 'USD')









    



The company should produce 52.94 units of product 1
The company should produce 14.12 units of product 2
The optimal profit is 1482.35 USD