Implementation in Python

Implementation in Python#

In this notebook we build and solve the sand & clay LP using Pyomo and HiGHS. You will learn how to declare the model in Pyomo (algebraic style), solve it with HiGHS, and inspect solutions.

Problem description#

A company extracts sand and clay from a site which when sold gives a profit of 57 and 60 monetary units per thousand units of product. For this extraction, sand needs a manpower of 50 men x hours to extract 1000 units of product 1 whilst clay needs 13 men x hours for 1000 units.
4h of backhoe work are needed to extract one thousand units of sand and 5h of backhoe work per thousand units of clay. The number of hours needed of truck transport is 8h and 4h respectively for sand and clay for each 1000 units of product transported.
The company has a work schedule of 40 hours per week for the men but also for the equipment (truck and backhoe). There are 5 men who can be used interchangeably between the transport of the two products. There is only one truck and one backhoe.
What should be the extraction plan of this company in a week? What should their objective be? What decisions need to be made?

Setting up the software

In this notebook we use Pyomo (modeling language) together with the open-source solver HiGHS. The second of these does not work as live code in the book, so in order to execute this notebook, you need to download the file and run it locally. You may first need to install the pyomo and highspy packages, using pip install or conda install:

  conda install pyomo highspy

For more information on these packages see:

Pyomo + HiGHS tutorial#

In this notebook, we are going to build a simple model using Pyomo + HiGHS.

Pyomo is a modeling language: you write math-like expressions in Python.
HiGHS is the solver: it optimizes the model (LP/MIP) you give it.
You can inspect the Pyomo class/method docstrings for full details. You can also Use Python’s help() to explore different objects:

import pyomo.environ as pyo
help(pyo.Var)          # decision variables
help(pyo.Objective)    # objective
help(pyo.Constraint)   # constraints

Import Pyomo + HiGHS package#

After preparing the environment for Pyomo + HiGHS, you should be able to run the following code to import the Pyomo + HiGHS package, where pyomo.environ is the main module of Pyomo.

import pyomo.environ as pyo

# (Optional) If you want a terser import style:

# from pyomo.environ import (
#     ConcreteModel, 
#     Var, 
#     Objective, 
#     Constraint,
#     NonNegativeReals, 
#     maximize, 
#     SolverFactory, 
#     value
# )

Construct the model#

Step 1 — Problem data (parameters)#

Pyomo separates data from the model structure so we clearly specify parameters before building the model. We first name the parameters from the problem. Writing them here (not hard-coding in expressions) makes the model readable & easy to change.

	Profit	Manpower	Backhoe work	Truck transport
Sand	57	50	4	8
Clay	60	13	5	4

# Profits per 1000 units
profit = {
    "sand": 57,
    "clay": 60,
}

# Hours needed per 1000 units
truck_hours   = {"sand": 8, "clay": 4}  # 1 truck available, 40h
backhoe_hours = {"sand": 4, "clay": 5}  # 1 backhoe available, 40h
worker_hours  = {"sand": 50, "clay": 13}# 5 workers * 40h = 200h total

# Weekly capacities
cap_truck   = 40
cap_backhoe = 40
cap_workers = 200

Tip: If you later add more materials, keep the same pattern and loop over sets.

Step 2 — Create a model container#

Think of this as an empty box where we put variables, objective, and constraints.

A ConcreteModel is an object that holds model components. It is “concrete” because all sets/parameters are defined in Python immediately.

m = pyo.ConcreteModel()

Step 3 — Declare decision variables#

To solve an optimization problem computationally, we must explicitly tell Pyomo what the unknowns are — i.e., what quantities the solver is allowed to adjust in order to optimize the objective while respecting the constraints.

In this sand & clay problem, our decisions are:

$x_1$ : thousands of units of sand produced and transported per week
$x_2$ : thousands of units of clay produced and transported per week

These quantities can take any real value ≥ 0 (e.g., 1.5 thousand units). Because fractional production is allowed, this is a continuous linear program (LP) rather than an integer (MILP) problem.

Variable domains in Pyomo#

When declaring a variable, we must specify its domain, i.e. the allowed type and range:

Pyomo domain	Meaning	Example use
`NonNegativeReals`	Any real number ≥ 0	Production, time
`Reals`	Any real number	Balance variables
`NonNegativeIntegers`	Integer ≥ 0	Integer production counts
`Integers`	Any integer	Net flow
`Binary`	0 or 1	Yes/No decisions

m.x1 = pyo.Var(domain=pyo.NonNegativeReals)   # sand (in 1000 units)
m.x2 = pyo.Var(domain=pyo.NonNegativeReals)   # clay (in 1000 units)

Step 4 — Add the objective#

We want to maximize weekly profit from selling sand and clay. $$ \begin{align} Max \quad L = 57 x_1 + 60 x_2 \end{align} $$

In Pyomo, we express this with Objective(expr=...). We also need to specify whether we are maximizing or minimizing. Pyomo defines:

sense=pyo.maximize or
sense=pyo.minimize

m.obj = pyo.Objective(
    expr = profit["sand"]*m.x1 + profit["clay"]*m.x2,
    sense = pyo.maximize
)

Step 5 — Add constraints#

Next, we specify limitations. These describe what combinations of $x_1$ and $x_2$ are allowed.

In Pyomo, we add a constraint using the Constraint()component:
Note that non-negativity constraints are not always necessary in Pyomo + HiGHS since we defined it in the domain of variables.

We have three set of constraints here:

Truck hours

\[\begin{split} \begin{align} & 8 x_1 + 4 x_2 \leq 40 \\ \end{align} \end{split}\]

Backhoe hours

\[\begin{split} \begin{align} & 4 x_1 + 5 x_2 \leq 40 \\ \end{align} \end{split}\]

Worker hours

\[\begin{split} \begin{align} & 50 x_1 + 13 x_2 \leq 200 \\ \end{align} \end{split}\]

# Trucks hours
m.truck = pyo.Constraint(
    expr = truck_hours["sand"]*m.x1 
         + truck_hours["clay"]*m.x2 <= cap_truck
)

# Backhoe hours
m.backhoe = pyo.Constraint(
    expr = backhoe_hours["sand"]*m.x1
         + backhoe_hours["clay"]*m.x2 <= cap_backhoe
)

# Worker hours
m.workers = pyo.Constraint(
    expr = worker_hours["sand"]*m.x1
         + worker_hours["clay"]*m.x2 <= cap_workers
)

Step 6 — Solve#

Use HiGHS through Pyomo’s SolverFactory. The tee flag prints the solver log

result = pyo.SolverFactory("highs").solve(m, tee=False)  # or "appsi_highs"
print("\nSolver status:", result.solver.status)
print("Termination  :", result.solver.termination_condition)

Step 7 — Read solution#

pyo.value() reads the numeric value of variables/expressions after solving.

x1 = pyo.value(m.x1)
x2 = pyo.value(m.x2)
obj = pyo.value(m.obj)

print(f"x1 (sand, 1000s units) = {x1:.3f}")
print(f"x2 (clay, 1000s units) = {x2:.3f}")
print(f"Max weekly profit      = {obj:.2f}")

Optional: Resource usage & binding constraints (slacks)#

Which constraints are tight/binding at the optimum? Here you can compute actual usage and slacks.

truck_use   = truck_hours["sand"]  * x1 + truck_hours["clay"]  * x2
backhoe_use = backhoe_hours["sand"]* x1 + backhoe_hours["clay"]* x2
worker_use  = worker_hours["sand"] * x1 + worker_hours["clay"] * x2

print("\nResource usage:")
print(f"  Truck   : {truck_use:.2f} / {cap_truck}")
print(f"  Backhoe : {backhoe_use:.2f} / {cap_backhoe}")
print(f"  Workers : {worker_use:.2f} / {cap_workers}")

print("\nBinding? (usage close to capacity):")
def binding(used, cap, tol=1e-6):
    return abs(used - cap) <= tol

print("  Truck   :", "YES" if binding(truck_use, cap_truck) else "NO")
print("  Backhoe :", "YES" if binding(backhoe_use, cap_backhoe) else "NO")
print("  Workers :", "YES" if binding(worker_use, cap_workers) else "NO")

Exercise: The sand and clay extraction problem#

Now that you have seen how to model and solve a simple problem in Pyomo + HiGHS, let’s try to make some changes in the problem. Imagine in the meantime the company has purchased another truck and has hired an extra (full time) worker to increase its profit. So you have to model and solve a very similar problem with minor changes to reflect the new employee and the new truck.

You need to make changes when defining attributes.