Some Linear Programming with Python and Gurobi

Let’s use gurobipy to solve a simple product mix problem

Smartphones Production Optimization

Linear programming (LP) is a mathematical method used to optimize a linear objective function, subject to a set of linear constraints. In simpler terms, it’s about maximizing profit or minimizing cost in a problem where variables have linear relationships.

For example, let’s imagine you are running a small factory that manufactures two smartphone models (budget and premium) and produces the components needed to build them.

These fantastic smartphones are composed of only two components:

• 1 screen
• 1 battery

The budget smartphones are composed of budget components and the premium ones are composed of premium components.

The screens department can’t work more than 120 hours per week and 1 hour is needed to produce 1 budget screen while 2 hours are needed to produce 1 premium screen.

The batteries department can’t work more than 90 hours per week and 1 hour is needed to produce 1 battery regardless of whether it is a budget or premium one.

The factory has 2 assembly lines:

• Line A assembles the budget smartphones and can’t assemble more than 70 units per week
• Line B assembles premium smartphones and can’t assemble more than 50 units per week

Each budget smartphone earns 200 Euro profit, and each premium smartphone earns 300 Euro. How many smartphones of each model should you produce weekly to maximize the profit?

Smartphone factory schema

These kinds of problems in the literature are referred to as product mix problems.

So, the basic structure of our linear programming problem is:

• Decision Variables (the things you can control): x is how many units we produce weekly of budget smartphones and y is how many units we produce weekly of premium smartphones
• Objective Function (what you want to optimize): P is the weekly profit, and we have that P = 200 * x + 300 * y
• Constraints (the limitations you have): our constraints are:

a. screens department work: 1*x + 2*y ≤ 120

b. batteries department work: 1*x + 1*y ≤ 90

c. line A production: x ≤ 70

d. line B production: y ≤ 50

e. we also must add the not negative constraint on x and y: x ≥ 0 and y ≥ 0

Instinctively we could say that since premium smartphones give us higher profits, we should maximize their production, so:

y=50

This implies that:

1. The screens department must work 100 hours to produce 50 premium screens per week, meaning it can work only 20 hours to produce the budget screens therefore it can produce 20 budget screens per week.
2. The batteries department must work 50 hours to produce 50 premium batteries per week, meaning it can work only 20 hours to produce the budget batteries therefore it can produce 20 budget batteries per week.
3. Line A assembles 20 budget smartphones per week
4. Line B assembles 50 premium smartphones per week

and we obtain the weekly profit P = 200 * 20 + 300 * 50 = 19000 Euro

Why do we need Linear Programming since we got our answer without using it? Because our instinct led us to a suboptimal solution.

I actually spoilered you the optimal solution in the first image of this story: by producing 60 budget smartphones (x = 60) and 30 premium smartphones (y = 30) we can get a higher profit:

Weekly profit P = 200 * 60 + 300 * 30 = 21000 Euro

You can check that these values of x and y respect the constraints, and they are indeed the optimal solution.

Since the problem is simple, we could do the math by hand but let’s stick to the purpose of this story and use Gurobi.

https://www.gurobi.com/

Gurobi is a commercial optimization solver designed to tackle linear programming (LP), mixed-integer programming (MIP), and quadratic programming (QP) problems. Developed by Gurobi Optimization LLC, it’s one of the fastest and most reliable solvers available and it is used by companies like Uber, Walmart, and the NFL for real-world decision-making.

Gurobi isn’t free, but they offer an academic license for students and researchers, and you can get a trial to test it out. More information on this page Academic Program and Licenses.

The gurobipy package comes with a trial license that allows you to solve problems of limited size.

import gurobipy as gp

# Create a new model
model = gp.Model("SmartphoneProduction")

# Add the decision variables
x = model.addVar(vtype=gp.GRB.INTEGER, name="budget_smartphones")
y = model.addVar(vtype=gp.GRB.INTEGER, name="premium_smartphones")

# Set objective function: maximize profit
model.setObjective(200 * x + 300 * y, gp.GRB.MAXIMIZE)

# Add constraints
model.addConstr(x + 2 * y <= 120, "screen_hours")  # Screens department
model.addConstr(x + y <= 90, "battery_hours")      # Batteries department 
model.addConstr(x <= 70, "linea_limit")            # Line A: budget smartphones limit
model.addConstr(y <= 50, "lineb_limit")            # Line B: premium smartphones limit
model.addConstr(x >= 0, "nonneg_x")                
model.addConstr(y >= 0, "nonneg_y")

# Save the model before optimization
model.write("smartphone_production_model.lp")  # LP file is a human-readable linear programming format
model.write("smartphone_production_model.mps") # MPS file is a standard optimization format

# Optimize the model
model.optimize()

if model.status == gp.GRB.OPTIMAL:
    opt_x, opt_y = x.x, y.x
    print(f"Optimal profit: €{model.objVal}")
    print(f"Produce {opt_x} budget smartphones and {opt_y} premium smartphones")
    model.write("smartphone_production_solution.sol") # SOL file is the solution after optimization
else:
    print("No optimal solution found")

Code execution output

LP format

MPS format

SOL format

So, we found that the optimal solution is producing 60 budget smartphones (x = 60) and 30 premium smartphones (y = 30) and the profit is 21000 Euro.

Let’s now visualize what is happening:

import numpy as np
import matplotlib.pyplot as plt

x_vals = np.linspace(0, 120, 400)

# Screens constraint: x + 2y <= 120 ---> y <= (120 - x)/2
y1 = (120 - x_vals) / 2

# Batteries constraint: x + y <= 90 ---> y <= 90 - x
y2 = 90 - x_vals

# Line A constraint: x <= 70
x_limit = 70

# Line B constraint: y <= 50
y_limit = 50

plt.figure(figsize=(8, 6))
plt.plot(x_vals, y1, label="x + 2y <= 120 (Screens)")
plt.plot(x_vals, y2, label="x + y <= 90 (Batteries)")
plt.axvline(x_limit, color='red', label="x <= 70 (Line A)")
plt.axhline(y_limit, color='purple', label="y <= 50 (Line B)")
plt.axvline(0, color='gray', linestyle='--', label="x >= 0")
plt.axhline(0, color='gray', linestyle='--', label="y >= 0")

# Feasible region
plt.fill_between(x_vals, np.minimum(np.minimum(y1, y2), y_limit), where=(y1 >= 0) & (y2 >= 0) & (x_vals <= x_limit), color='lightgreen', alpha=0.5)

# Optimal point
if model.status == gp.GRB.OPTIMAL:
    plt.plot(opt_x, opt_y, 'ro', label=f"Optimal ({opt_x}, {opt_y})")
else:
    pass

plt.xlabel("Budget smartphones (x)")
plt.ylabel("Premium smartphones (y)")
plt.title("Smartphones Production Optimization")
plt.legend()
plt.grid(True)
plt.savefig("smartphone_production.png")
plt.show()

Smartphones Production Optimization

This is the same image you can see at the beginning of the story. In light green we can see the feasible region also called solution space which is the set of all possible points of an optimization problem that satisfy the problem’s constraints. The colored lines are the constraints, the dashed lines are the not negative constraints, and the red dot is the optimal solution.

Let’s now visualize the relationship between our feasible region and the profit:

An isocost line shows all combinations of inputs which cost the same total amount. By setting different values of the cost we get parallels lines. In our case an isocost line will be all the points that will get us the same profit.

# Plot isocost lines
plt.figure(figsize=(8, 6))

# Feasible region
plt.fill_between(x_vals, np.minimum(np.minimum(y1, y2), y_limit), where=(y1 >= 0) & (y2 >= 0) & (x_vals <= x_limit), color='lightgreen', alpha=0.5)

# Optimal point
if model.status == gp.GRB.OPTIMAL:
    plt.plot(opt_x, opt_y, 'ro', label=f"Optimal ({opt_x}, {opt_y})")
else:
    pass

X, Y = np.meshgrid(np.linspace(0, 125, 100), np.linspace(0, 85, 100))
Z = 200 * X + 300 * Y  # Profit function
profits = [6000, 14000, 19000, 21000]  # Example profit levels, including optimal
isocost_lines = plt.contour(X, Y, Z, levels=profits, colors='blue', linestyles='dashed')
plt.clabel(isocost_lines, inline=True, fontsize=10, fmt="€%d")

plt.xlabel("Budget smartphones (x)")
plt.ylabel("Premium smartphones (y)")
plt.title("Smartphones Production Optimization")
plt.legend()
plt.grid(True)
plt.savefig("smartphone_production_isocost_lines.png")
plt.show()

Isocost lines

We have the same feasible region in light green as before and each dashed blue line is an isocost line, i.e C6000, C14000, C19000 and C21000.

Note that C19000 represents the profit we got using the first suboptimal solution (x=20 and y=50).

C21000 represents the optimal profit, and it intersects the feasible region in the optimal point (x=60 and y=30).

Outro

I hope the story was interesting and thank you for taking the time to read it. On my Blogspot you can find the same post in Italian. Let me know if you have any question and if you like the content that I create feel free to buy me a coffee.