Fast and Efficient MIP Modelling: A Guide to PyOptInterface

Python has become the language of choice for operations research practitioners due to its versatility and rich ecosystem of libraries. Many of you might be familiar with popular modelling libraries like PuLP and Pyomo, which are often used with open-source solvers like CBC and HiGHS. When you need to use commercial solvers, you might have used gurobipy and docplex libraries for Gurobi and CPLEX solvers respectively. However, they come with a significant cost, making them less accessible for many practitioners and researchers.

The total run time of any optimization problem can be divided into model building time and model solving time. Commercial solvers are undeniably fast compared with open-source solvers, and gurobipy often provides faster model building compared with PuLP and Pyomo. However, you might encounter scenarios where your use case requires the use of open-source solvers but also demands faster model building.

Here comes PyOptInterface, a new and promising library designed to offer the best of both worlds. It provides a user-friendly syntax similar to PuLP and Pyomo, making it easy for you to transition to it. But the standout feature of this library is its performance. Initial benchmarks suggest that PyOptInterface can outperform even gurobipy in terms of model building time.

In the sections that follow, we'll walk you through a step-by-step guide on how to use PyOptInterface to build and solve a standard MIP model.

For this, we'll solve a classic MIP problem: the knapsack problem. Given a set of items, each with a weight and value, the objective is to maximize the total value without exceeding a weight limit. We'll use randomly generated data for our example.

Mathematical Model

Data Preparation

Let's generate sample data for the problem.

from random import randint

n_items = 10
weights = [randint(1, 10) for _ in range(n_items)]
values = [randint(1, 20) for _ in range(n_items)]
max_weight = 30

Model Object

A model is an instance tied to a specific solver. To create a HiGHS model, we need to import the corresponding module and call the constructor of the model class.

import pyoptinterface as poi
from pyoptinterface import highs, quicksum

model = highs.Model()

Decision Variables

Next, we need to add decision variables. A decision variable is defined with three main properties: its type (continuous, integer or binary), its upper bound, and its lower bound. For the above example, we can define decision variables as:

x = {i: model.add_variable(name=f"x_{i}", 
						  domain=poi.VariableDomain.Binary) 
	 for i in range(n_items)}

To add continuous or integer variables, modify domain as poi.VariableDomain.Continuous or poi.VariableDomain.Integer respectively. For such variables you can add upper and lower bounds using ub and lb attributes.

Constraint

model.add_linear_constraint(
	quicksum(weights[i] * x[i] for i in range(n_items)), 
	poi.Leq, max_weight, name="Weight_Constraint")

You can use poi.Geq for greater than or equal to sense, and poi.Eq for equal to sense of the constraints. If you're using v0.2.7 or older version of pyoptinterface, you need to keep right hand side of the constraint as a float value.

Objective Function

model.set_objective(
    quicksum(values[i] * x[i] for i in range(n_items)),
    poi.ObjectiveSense.Maximize)

To minimize the objective function, keep sense as poi.ObjectiveSense.Minimize.

Storing the formulated Model

Sometimes, we need to store the MIP model for debugging purpose. The model can be written to a file in LP, MPS or other formats. The write method of the model can be used to write the model to file.

model.write("model.lp")   # Stores LP file
model.write("model.mps")  # Stores MPS file

Solving the Model

We can solve the model by calling the optimize method of the model. It will invoke the HiGHS solver.

model.optimize()
solution_status = model.get_model_attribute(
					poi.ModelAttribute.TerminationStatus)

Getting Values of Variables

For value of variables and expressions, we can use get_value method of the model.

if solution_status == poi.TerminationStatusCode.OPTIMAL:
    solution = {f"x{i}": model.get_value(x[i]) 
						for i in range(n_items) 
						if model.get_value(x[i])==1}
    print(solution)
else:
    print("No optimal solution found.")

In this blog, we've demonstrated how to build and solve a standard MIP problem using the PyOptInterface library. We only covered high-level modelling in Python. For detailed documentation of PyOptInterface, visit their official page.