Some Linear Programming with Python and Gurobi

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? 

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.

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")

This is the output:

Restricted license - for non-production use only - expires 2026-11-23

Gurobi Optimizer version 12.0.1 build v12.0.1rc0 (win64 - Windows 10.0 (19045.2))

CPU model: Intel(R) Core(TM) i7-6700HQ CPU @ 2.60GHz, instruction set [SSE2|AVX|AVX2]

Thread count: 4 physical cores, 8 logical processors, using up to 8 threads

Optimize a model with 6 rows, 2 columns and 8 nonzeros

Model fingerprint: 0xa3ac08e8

Variable types: 0 continuous, 2 integer (0 binary)

Coefficient statistics:

  Matrix range     [1e+00, 2e+00]

  Objective range  [2e+02, 3e+02]

  Bounds range     [0e+00, 0e+00]

  RHS range        [5e+01, 1e+02]

Found heuristic solution: objective 20000.000000

Presolve removed 4 rows and 0 columns

Presolve time: 0.01s

Presolved: 2 rows, 2 columns, 4 nonzeros

Variable types: 0 continuous, 2 integer (0 binary)

Root relaxation: objective 2.100000e+04, 2 iterations, 0.00 seconds (0.00 work units)

    Nodes    |    Current Node    |     Objective Bounds      |     Work

 Expl Unexpl |  Obj  Depth IntInf | Incumbent    BestBd   Gap | It/Node Time

...

Optimal solution found (tolerance 1.00e-04)

Best objective 2.100000000000e+04, best bound 2.100000000000e+04, gap 0.0000%

Optimal profit: €21000.0

Produce 60.0 budget smartphones and 30.0 premium smartphones

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()

To check the image please got to the Medium 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()

To check the image please got to the Medium story. 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 post was interesting and thank you for taking the time to read it. On my Medium you can find a more in depth story and 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.