37. Linear Programming#

37.2.1. Computation: using OR-Tools#

Let’s try to solve the same problem using the package ortools.linear_solver.

The following cell instantiates a solver and creates two variables specifying the range of values that they can have.

# Instantiate a GLOP(Google Linear Optimization Package) solver
solver = pywraplp.Solver.CreateSolver('GLOP')

Let’s create two variables $x_1$ and $x_2$ such that they can only have nonnegative values.

# Create the two variables and let them take on any non-negative value.
x1 = solver.NumVar(0, solver.infinity(), 'x1')
x2 = solver.NumVar(0, solver.infinity(), 'x2')

Add the constraints to the problem.

# Constraint 1: 2x_1 + 5x_2 <= 30.0
solver.Add(2 * x1 + 5 * x2 <= 30.0)

# Constraint 2: 4x_1 + 2x_2 <= 20.0
solver.Add(4 * x1 + 2 * x2 <= 20.0)

<ortools.linear_solver.pywraplp.Constraint; proxy of <Swig Object of type 'operations_research::MPConstraint *' at 0x7fad27d3e910> >

Let’s specify the objective function. We use solver.Maximize method in the case when we want to maximize the objective function and in the case of minimization we can use solver.Minimize.

# Objective function: 3x_1 + 4x_2
solver.Maximize(3 * x1 + 4 * x2)

Once we solve the problem, we can check whether the solver was successful in solving the problem using its status. If it’s successful, then the status will be equal to pywraplp.Solver.OPTIMAL.

# Solve the system.
status = solver.Solve()

if status == pywraplp.Solver.OPTIMAL:
    print('Objective value =', solver.Objective().Value())
    print(f'(x1, x2): ({x1.solution_value():.2}, {x2.solution_value():.2})')
else:
    print('The problem does not have an optimal solution.')

Objective value = 27.5
(x1, x2): (2.5, 5.0)

37.3.1. Computation: using OR-Tools#

Let’s try to solve the above problem using the package ortools.linear_solver.

The following cell instantiates a solver and creates two variables specifying the range of values that they can have.

# Instantiate a GLOP(Google Linear Optimization Package) solver
solver = pywraplp.Solver.CreateSolver('GLOP')

Let’s create five variables $x_1, x_2, x_3, x_4,$ and $x_5$ such that they can only have the values defined in the above constraints.

# Create the variables using the ranges available from constraints
x1 = solver.NumVar(0, solver.infinity(), 'x1')
x2 = solver.NumVar(-20_000, solver.infinity(), 'x2')
x3 = solver.NumVar(-20_000, solver.infinity(), 'x3')
x4 = solver.NumVar(-20_000, solver.infinity(), 'x4')
x5 = solver.NumVar(0, 50_000, 'x5')

Add the constraints to the problem.

# Constraint 1: x_1 + x_2 = 100,000
solver.Add(x1 + x2 == 100_000.0)

# Constraint 2: x_1 - 1.06 * x_2 + x_3 + x_5 = 0
solver.Add(x1 - 1.06 * x2 + x3 + x5 == 0.0)

# Constraint 3: x_1 - 1.06 * x_3 + x_4 = 0
solver.Add(x1 - 1.06 * x3 + x4 == 0.0)

<ortools.linear_solver.pywraplp.Constraint; proxy of <Swig Object of type 'operations_research::MPConstraint *' at 0x7fad27d3f390> >

Let’s specify the objective function.

# Objective function: 1.30 * 3 * x_1 + 1.06 * x_4 + 1.30 * x_5
solver.Maximize(1.30 * 3 * x1 + 1.06 * x4 + 1.30 * x5)

Let’s solve the problem and check the status using pywraplp.Solver.OPTIMAL.

# Solve the system.
status = solver.Solve()

if status == pywraplp.Solver.OPTIMAL:
    print('Objective value =', solver.Objective().Value())
    x1_sol = round(x1.solution_value(), 3)
    x2_sol = round(x2.solution_value(), 3)
    x3_sol = round(x1.solution_value(), 3)
    x4_sol = round(x2.solution_value(), 3)
    x5_sol = round(x1.solution_value(), 3)
    print(f'(x1, x2, x3, x4, x5): ({x1_sol}, {x2_sol}, {x3_sol}, {x4_sol}, {x5_sol})')
else:
    print('The problem does not have an optimal solution.')

Objective value = 141018.24349792692
(x1, x2, x3, x4, x5): (24927.755, 75072.245, 24927.755, 75072.245, 24927.755)

OR-Tools tells us that the best investment strategy is:

1. At the beginning of the first year, the mutual fund should buy $\$24,927.755$ of the annuity. Its bank account balance should be $\$75,072.245$.

2. At the beginning of the second year, the mutual fund should buy $\$24,927.755$ of the corporate bond and keep invest in the annuity. Its bank balance should be $\$24,927.755$.

3. At the beginning of the third year, the bank balance should be $\$75,072.245$.

4. At the end of the third year, the mutual fund will get payouts from the annuity and corporate bond and repay its loan from the bank. At the end it will own $\$141,018.24$, so that it’s total net rate of return over the three periods is $41.02\%$.

37.4.1. Useful transformations#

It is useful to know how to transform a problem that initially is not stated in the standard form into one that is.

By deploying the following steps, any linear programming problem can be transformed into an equivalent standard form linear programming problem.

1. Objective function: If a problem is originally a constrained maximization problem, we can construct a new objective function that is the additive inverse of the original objective function. The transformed problem is then a minimization problem.

2. Decision variables: Given a variable $x_j$ satisfying $x_j \le 0$, we can introduce a new variable $x_j' = - x_j$ and substitute it into original problem. Given a free variable $x_i$ with no restriction on its sign, we can introduce two new variables $x_j^+$ and $x_j^-$ satisfying $x_j^+, x_j^- \ge 0$ and replace $x_j$ by $x_j^+ - x_j^-$.

3. Inequality constraints: Given an inequality constraint $\sum_{j=1}^n a_{ij}x_j \le 0$, we can introduce a new variable $s_i$, called a slack variable that satisfies $s_i \ge 0$ and replace the original constraint by $\sum_{j=1}^n a_{ij}x_j + s_i = 0$.

Let’s apply the above steps to the two examples described above.

37.4.2. Example 1: production problem#

The original problem is:

This problem is equivalent to the following problem with a standard form:

37.4.3. Computation: using SciPy#

The package scipy.optimize provides a function linprog to solve linear programming problems with a form below:

$A_{eq}, b_{eq}$ denote the equality constraint matrix and vector, and $A_{ub}, b_{ub}$ denote the inequality constraint matrix and vector.

Let’s now try to solve the Problem 1 using SciPy.

# Construct parameters
c_ex1 = np.array([3, 4])

# Inequality constraints
A_ex1 = np.array([[2, 5],
                  [4, 2]])
b_ex1 = np.array([30,20])

Once we solve the problem, we can check whether the solver was successful in solving the problem using the boolean attribute success. If it’s successful, then the success attribute is set to True.

# Solve the problem
# we put a negative sign on the objective as linprog does minimization
res_ex1 = linprog(-c_ex1, A_ub=A_ex1, b_ub=b_ex1)

if res_ex1.success:
    # We use negative sign to get the optimal value (maximized value)
    print('Optimal Value:', -res_ex1.fun)
    print(f'(x1, x2): {res_ex1.x[0], res_ex1.x[1]}')
else:
    print('The problem does not have an optimal solution.')

Optimal Value: 27.5
(x1, x2): (2.5, 5.0)

The optimal plan tells the factory to produce $2.5$ units of Product 1 and $5$ units of Product 2; that generates a maximizing value of revenue of $27.5$.

We are using the linprog function as a black box.

Inside it, Python first transforms the problem into standard form.

To do that, for each inequality constraint it generates one slack variable.

Here the vector of slack variables is a two-dimensional NumPy array that equals $b_{ub} - A_{ub}x$.

See the official documentation for more details.

37.4.4. Example 2: investment problem#

The original problem is:

This problem is equivalent to the following problem with a standard form:

# Construct parameters
rate = 1.06

# Objective function parameters
c_ex2 = np.array([1.30*3, 0, 0, 1.06, 1.30])

# Inequality constraints
A_ex2 = np.array([[1,  1,  0,  0,  0],
                  [1, -rate, 1, 0, 1],
                  [1, 0, -rate, 1, 0]])
b_ex2 = np.array([100_000, 0, 0])

# Bounds on decision variables
bounds_ex2 = [(  0,    None),
              (-20_000, None),
              (-20_000, None),
              (-20_000, None),
              (  0,   50_000)]

Let’s solve the problem and check the status using success attribute.

# Solve the problem
res_ex2 = linprog(-c_ex2, A_eq=A_ex2, b_eq=b_ex2,
                  bounds=bounds_ex2)

if res_ex2.success:
    # We use negative sign to get the optimal value (maximized value)
    print('Optimal Value:', -res_ex2.fun)
    x1_sol = round(res_ex2.x[0], 3)
    x2_sol = round(res_ex2.x[1], 3)
    x3_sol = round(res_ex2.x[2], 3)
    x4_sol = round(res_ex2.x[3], 3)
    x5_sol = round(res_ex2.x[4], 3)
    print(f'(x1, x2, x3, x4, x5): {x1_sol, x2_sol, x3_sol, x4_sol, x5_sol}')
else:
    print('The problem does not have an optimal solution.')

Optimal Value: 141018.24349792697
(x1, x2, x3, x4, x5): (24927.755, 75072.245, 4648.825, -20000.0, 50000.0)

SciPy tells us that the best investment strategy is:

1. At the beginning of the first year, the mutual fund should buy $\$24,927.75$ of the annuity. Its bank account balance should be $\$75,072.25$.

2. At the beginning of the second year, the mutual fund should buy $\$50,000$ of the corporate bond and keep invest in the annuity. Its bank account balance should be $\$ 4,648.83$.

3. At the beginning of the third year, the mutual fund should borrow $\$20,000$ from the bank and invest in the annuity.

4. At the end of the third year, the mutual fund will get payouts from the annuity and corporate bond and repay its loan from the bank. At the end it will own $\$141,018.24$, so that it’s total net rate of return over the three periods is $41.02\%$.