Linear Programming in Python

Linear Programming in Python#

Part 1: Overview and Mathematical Formulation#

A civil engineering company wants to decide on the projects that they should do. Their objective is to minimize the environmental impact of their projects while making enough profit to keep the company running.

They have a portfolio of 6 possible projects to invest in, where A, B , and C are new infrastructure projects (so-called type 1), and D, E, F are refurbishment projects (so-called type 2).

The environmental impact of each project is given by \(I_i\) where \(i \in [1,(...),6]\) is the index of the project. \(I_i=[90,45,78,123,48,60]\)

The profit of each project is given by \(P_i\) where \(i\in [1,(...),6]\) is the index of the project: \(P_i=[120,65,99,110,33,99]\)

The company is operating with the following constraints, please formulate the mathematical program that allows solving the problem:

  • The company wants to do 3 out of the 6 projects

  • the projects of type 2 must be at least as many as the ones of type 1

  • the profit of all projects together must be greater or equal than \(250\) (\(\beta\))

Task 1: Writing the mathematical formulation

Write down every formulation and constraint that is relevant to solve this optimization problem.

Solution
\[ Min(Z) \sum_{i=1}^{6} x_i I_i \]

subject to

\[ \sum_{i=1}^{6} x_i = 3 \]
\[ \sum_{i=4}^{6} x_i \ge \sum_{i=1}^{3} x_i \]
\[ \sum_{i=1}^{6} x_i P_i \ge β \]
\[x\in(1,0), i\in(1,....,6)\]

The constraint related to: Extra Challenge (Task 6) $\( \sum_{i=1}^{6} x_i I_i \le x_1 γ + (1-x_1) M \)$

where M is a sufficiently big number such as 9999

Task 2: Setting up the software

We will use Pyomo (modeling language) together with the open-source solver HiGHS. No licenses are required.

What you need:

  • Python environment (same as previous weeks)
  • Install the packages: pyomo and highspy

Installation (if needed):

pip install pyomo highspy

from pyomo.environ import (
    ConcreteModel, Set, Var, Param, Objective, Constraint,
    Binary, NonNegativeReals, Integers, minimize, maximize,
    SolverFactory, value
)

Task 3: Setting up the problem

Define any inputs you might need to setup your model.

# Project data
I = [1, 2, 3, 4, 5, 6]   # project index set

impact = {1: 90,
          2: 45,
          3: 78,
          4: 123,
          5: 48,
          6: 60}

profit = {1: 120,
          2: 65,
          3: 99,
          4: 110,
          5: 33,
          6: 99}

# Minimum required profit
beta = 250
M = 100000 #Later

type1 = {1, 2, 3}    # type-1 projects
type2 = {4, 5, 6}    # type-2 projects

# Number of projects and types
num_projects = len(I)
num_type1 = len(type1)
num_type2 = len(type2)

Part 2: Create model with Pyomo + HiGHS#

We will follow these steps (Pyomo separates modeling from solving):

  1. Define the model container (ConcreteModel)

  2. Define sets and parameters (Param or Python Dictionary)

  3. Define decision variables (Var)

  4. Define the objective (Objective)

  5. Add constraints (Constraint)

  6. Call the solver (SolverFactory("highs").solve(m, tee=False))

If you need help on a Pyomo component, use Python’s help:

import pyomo.environ as pyo
help(pyo.Var)

Task 4: Create the Pyomo model

Create the Pyomo model, declare the decision variables, objective function, and constraints as specified in the formulation. Use the code skeleton provided to guide you. You can also refer to the Pyomo documentation for additional details of implementation.

# Build the Pyomo model using the inputs from Task 3

m = ConcreteModel()

# Sets
m.I = Set(initialize=I)

# Decision variables
m.x = Var(m.I, domain=Binary)   # 1 if project i is selected

# Objective: minimize total environmental impact
m.obj = Objective(expr=sum(impact[i] * m.x[i] for i in m.I), sense=minimize)

# Constraints
# (1) Select exactly three projects
m.cardinality = Constraint(expr=sum(m.x[i] for i in m.I) == 3)

# (2) Type-2 projects ≥ Type-1 projects
m.type_balance = Constraint(
    expr=sum(m.x[i] for i in type2) >= sum(m.x[i] for i in type1)
)

# (3) Total profit ≥ beta
m.min_profit = Constraint(
    expr=sum(profit[i] * m.x[i] for i in m.I) >= beta
)

Solution.

m = ConcreteModel()
Creates the Pyomo model container. All sets, variables, objective, and constraints are attached to m.

m.I = Set(initialize=I)
Declares the index set of projects using the data vector I defined in Task 3 (e.g., [1,2,3,4,5,6]).

m.x = Var(m.I, domain=Binary)
Creates a binary decision variable for each project i ∈ I.

m.obj = Objective(…, sense=minimize)
Defines the objective to minimize total environmental impact

m.cardinality = Constraint(expr = sum(x[i] for i in I) == 3)
Enforces that exactly three projects are selected in total

m.type_balance = Constraint(expr = sum(x[i] for i in type2) >= sum(x[i] for i in type1))
Ensures at least as many type-2 projects as type-1 projects

m.min_profit = Constraint(expr = sum(profit[i] * x[i] for i in I) >= beta)
Imposes the minimum profit requirement

Notes: All parameters (I, impact, profit, type1, type2, beta) are defined in the Task 3 cell. Solving is executed in a later cell with SolverFactory(“highs”).solve(m, tee=False).

Task 5: Display your results

Display the model in a good way to interpret and print the solution of the optimization.

result = SolverFactory("highs").solve(m, tee=False)
print("Status:", result.solver.status)
print("Termination:", result.solver.termination_condition)

selected = [i for i in m.I if value(m.x[i]) > 0.5] 
total_profit = sum(profit[i] * value(m.x[i]) for i in m.I)
total_impact = sum(impact[i] * value(m.x[i]) for i in m.I)

print("\nSelected projects:", selected)
print("Total profit:", total_profit)
print("Total environmental impact:", total_impact)
print("Type-1 selected:", sum(value(m.x[i]) for i in type1))
print("Type-2 selected:", sum(value(m.x[i]) for i in type2))

Status: ok
Termination: optimal

Selected projects: [1, 5, 6]
Total profit: 252.0
Total environmental impact: 198.0
Type-1 selected: 1.0
Type-2 selected: 2.0

Task 6: Additional constraint

Solve the model with an additional constraint: if project 1 is done then the impact of all projects together should be lower than \(\gamma\) with \(\gamma=130\).

In the first cell you should add the constraint, then in a second cell optimize the model.

gamma = 130
m.impact_conditional = Constraint(
    expr = sum(impact[i] * m.x[i] for i in m.I) <= m.x[1] * gamma + (1 - m.x[1]) * M
)

result2 = SolverFactory("highs").solve(m, tee=False)
print("Status:", result2.solver.status)
print("Termination:", result2.solver.termination_condition)

selected2 = [i for i in m.I if value(m.x[i]) > 0.5]
total_profit2 = sum(profit[i] * value(m.x[i]) for i in m.I)
total_impact2 = sum(impact[i] * value(m.x[i]) for i in m.I)

print("\nSelected projects (new):", selected2)
print("Total profit (new):", total_profit2)
print("Total environmental impact (new):", total_impact2)

Status: ok
Termination: optimal

Selected projects (new): [2, 4, 6]
Total profit (new): 274.0
Total environmental impact (new): 228.0

By Nadia Pourmohammadzia, Delft University of Technology. CC BY 4.0, more info on the Credits page of Workbook.

Contents