Use Python SciPy Linprog For Linear Programming Optimization

Recently, I was working on an optimization problem where I needed to maximize profits while dealing with multiple constraints in a manufacturing setup. This is when I turned to SciPy’s linear programming functionality.

Linear programming is a powerful mathematical technique for finding the best outcome in a model with linear relationships. And Python’s SciPy library makes this complex task surprisingly accessible.

In this article, I’ll walk you through how to use SciPy’s linprog function to solve real-world optimization problems. I’ll cover everything from basic usage to advanced techniques with practical examples.

What is Linear Programming and SciPy’s Linprog?

Linear programming helps you find the optimal value (maximum or minimum) of a linear objective function, subject to linear constraints.

SciPy’s linprog function is designed specifically for minimization problems. But don’t worry – I’ll show you how to adapt it for maximization problems too.

The general form of a linear programming problem that linprog solves is:

Minimize: c @ x
Subject to: A_ub @ x <= b_ub
           A_eq @ x == b_eq
           bounds[i][0] <= x[i] <= bounds[i][1]

Let’s break down these components and then see them in action.

Read SciPy Stats

Methods to Use Python SciPy Linprog

Now, let’s import the necessary modules:

from scipy.optimize import linprog
import numpy as np

1: Basic Minimization with Linprog

Let’s start with a simple example, imagine you’re managing a small factory that produces two types of products: widgets and gadgets.

Each widget requires 2 hours of labor and 1 hour of machine time
Each gadget requires 1 hour of labor and 3 hours of machine time
You have 100 hours of labor and 90 hours of machine time available
Widgets earn $3 profit each, gadgets earn $5 profit each

Here’s how to find the optimal production plan:

# Define the objective function coefficients (negative for maximization)
c = [-3, -5]  # Negative because linprog minimizes by default

# Define the constraint coefficients
A = [[2, 1],  # Labor constraint
     [1, 3]]  # Machine time constraint

# Define the right-hand side of constraints
b = [100, 90]  # Available resources

# Set bounds for variables (can't produce negative amounts)
x_bounds = [(0, None), (0, None)]

# Solve the linear programming problem
result = linprog(c, A_ub=A, b_ub=b, bounds=x_bounds, method="highs")

print("Optimal solution found:")
print(f"Widgets to produce: {result.x[0]}")
print(f"Gadgets to produce: {result.x[1]}")
print(f"Maximum profit: ${-result.fun}")  # Negate back to get actual profit

Output:

Optimal solution found:
Widgets to produce: 42.0
Gadgets to produce: 16.0
Maximum profit: $206.0

You can refer to the screenshot below to see the output:

The method="highs" parameter specifies that we want to use the HiGHS solver, which is more efficient and robust than the default.

Check out SciPy Misc

2: Use Equality Constraints

Sometimes, you need exact equality constraints. For instance, if you need to produce exactly 30 total units due to packaging requirements:

# Same objective and inequality constraints as before
c = [-3, -5]
A = [[2, 1], [1, 3]]
b = [100, 90]

# Add equality constraint: widgets + gadgets = 30
A_eq = [[1, 1]]
b_eq = [30]

# Solve
result = linprog(c, A_ub=A, b_ub=b, A_eq=A_eq, b_eq=b_eq, 
                bounds=[(0, None), (0, None)], method="highs")

print("Optimal solution with equality constraint:")
print(f"Widgets to produce: {result.x[0]}")
print(f"Gadgets to produce: {result.x[1]}")
print(f"Maximum profit: ${-result.fun}")

Output:

Optimal solution with equality constraint:
Widgets to produce: 0.0
Gadgets to produce: 30.0
Maximum profit: $150.0

You can refer to the screenshot below to see the output:

3: Set Specific Bounds

The bounds parameter in linprog allows you to set specific limits for each variable. This is useful when you have production minimums or maximums:

# Bounds: at least 10 widgets, at most 40 gadgets
x_bounds = [(10, None), (0, 40)]

result = linprog(c, A_ub=A, b_ub=b, bounds=x_bounds, method="highs")

print("Optimal solution with specific bounds:")
print(f"Widgets to produce: {result.x[0]}")
print(f"Gadgets to produce: {result.x[1]}")
print(f"Maximum profit: ${-result.fun}")

Output:

Optimal solution with specific bounds:
Widgets to produce: 42.0
Gadgets to produce: 16.0
Maximum profit: $206.0

You can refer to the screenshot below to see the output:

As noted in Python Guides, bounds are specified as (min, max) pairs, and you can use None to indicate no bound in a particular direction.

Check out SciPy Integrate

4: Portfolio Optimization Example

Let’s apply linprog to a real-world financial problem: portfolio optimization.

Say you have $10,000 to invest across three assets:

Asset A: Expected annual return 5%, risk level 2
Asset B: Expected annual return 8%, risk level 5
Asset C: Expected annual return 12%, risk level 9

You want to maximize returns while keeping risk below a certain threshold:

# Returns (to maximize, so negative for linprog)
c = [-0.05, -0.08, -0.12]

# Risk constraint: total risk score must be <= 50
A = [[2, 5, 9]]
b = [50]

# Constraint: total investment = $10,000
A_eq = [[1, 1, 1]]
b_eq = [10000]

# Non-negative investments
bounds = [(0, None), (0, None), (0, None)]

result = linprog(c, A_ub=A, b_ub=b, A_eq=A_eq, b_eq=b_eq, bounds=bounds, method="highs")

print("Optimal portfolio allocation:")
print(f"Asset A: ${result.x[0]:.2f}")
print(f"Asset B: ${result.x[1]:.2f}")
print(f"Asset C: ${result.x[2]:.2f}")
print(f"Expected annual return: ${-result.fun * 10000:.2f}")

5: Use Different Solver Methods

SciPy’s linprog offers multiple solver methods, each with its advantages:

# Try different methods
methods = ["highs", "highs-ds", "highs-ipm", "interior-point", "simplex", "revised simplex"]

for method in methods:
    try:
        result = linprog(c, A_ub=A, b_ub=b, bounds=x_bounds, method=method)
        print(f"Method {method}: Widgets={result.x[0]:.2f}, Gadgets={result.x[1]:.2f}, Profit=${-result.fun:.2f}")
    except Exception as e:
        print(f"Method {method} failed: {e}")

The HiGHS solver is generally recommended for most problems as it’s faster and more reliable.

Read SciPy Signal

6: Handle Infeasible Problems

Sometimes, your constraints might be too strict, making the problem infeasible. Let’s see how to detect and handle this:

# Impossible constraints
A_impossible = [[2, 1], [1, 3]]
b_impossible = [10, 10]  # Not enough resources

result = linprog(c, A_ub=A_impossible, b_ub=b_impossible, 
                bounds=[(0, None), (0, None)], method="highs")

if result.success:
    print("Solution found:", result.x)
else:
    print("Problem is infeasible:", result.message)

    # You might want to relax constraints and try again
    # For example, increase available resources by 50%
    b_relaxed = [b * 1.5 for b in b_impossible]

    result_relaxed = linprog(c, A_ub=A_impossible, b_ub=b_relaxed, 
                           bounds=[(0, None), (0, None)], method="highs")

    if result_relaxed.success:
        print("Solution with relaxed constraints:", result_relaxed.x)

7: Transportation Problem

The transportation problem is a classic linear programming application. Let’s say you have factories in New York and Los Angeles, and distribution centers in Chicago, Dallas, and Miami.

You need to determine the optimal shipping plan to minimize costs:

# Shipping costs from factories to distribution centers
#              Chicago  Dallas  Miami
costs = np.array([
    [10, 30, 50],  # New York
    [70, 30, 10]   # Los Angeles
])

# Supply at each factory
supply = np.array([100, 150])  # NY, LA

# Demand at each distribution center
demand = np.array([80, 70, 100])  # Chicago, Dallas, Miami

# Create the objective function coefficients
c = costs.flatten()

# Build the constraint matrices
A_eq = []

# Supply constraints
for i in range(len(supply)):
    row = [0] * (len(supply) * len(demand))
    for j in range(len(demand)):
        row[i * len(demand) + j] = 1
    A_eq.append(row)

# Demand constraints
for j in range(len(demand)):
    row = [0] * (len(supply) * len(demand))
    for i in range(len(supply)):
        row[i * len(demand) + j] = 1
    A_eq.append(row)

# Right-hand side of constraints
b_eq = np.concatenate([supply, demand])

# Non-negative shipments
bounds = [(0, None)] * (len(supply) * len(demand))

# Solve
result = linprog(c, A_eq=A_eq, b_eq=b_eq, bounds=bounds, method="highs")

# Reshape the result to match our original problem
shipments = result.x.reshape(len(supply), len(demand))

print("Optimal shipping plan:")
for i, factory in enumerate(["New York", "Los Angeles"]):
    for j, center in enumerate(["Chicago", "Dallas", "Miami"]):
        print(f"Ship {shipments[i, j]:.1f} units from {factory} to {center}")
print(f"Total shipping cost: ${result.fun:.2f}")

Read SciPy Convolve

Understand the Linprog Output

When you run linprog, it returns an object with several attributes:

x: The optimal solution array
fun: The optimal value of the objective function
success: Boolean indicating if the optimizer succeeded
status: Integer indicating the exit status
message: String describing the exit status
nit: Number of iterations performed
slack: Slack in the inequality constraints at the solution

Understanding these outputs helps you interpret your results and troubleshoot issues.

I hope you found this article helpful in understanding how to use SciPy’s linprog function for linear programming. It’s a powerful tool that can optimize countless real-world problems, from production planning to resource allocation and financial optimization.

Whether you’re maximizing profits, minimizing costs, or finding the optimal balance of resources, linprog provides a clean and efficient way to solve these problems in Python.