Notebooks and Examples

Use Case Tutorials

In the following, you find notebooks with use cases to get you started with the HybridSolver.

• QUBOs Getting Started: Learn how to solve QUBOs to global optimality with Quantagonia’s HybridSolver.

• The Drone Delivery Problem: Discover how to optimize drone delivery models.

• The Facility Layout Planning Problem: Explore how to optimize facility layout planning models.

• The Maintenance Scheduling Problem: Master how to optimize maintenance scheduling models with efficient strategies.

• The Workforce Scheduling Problem: Learn how to optimize workforce scheduling models.

Code Examples

The following examples introduce some of the functionalities of the quantagonia package.

Modeling problems

MIP

This example shows how to model a MIP, solve it with HybridSolver and retrieve the solution status, the optimal solution, and the best objective value computed.

from quantagonia.enums import HybridSolverOptSenses, VarType
from quantagonia.mip.model import Model

# setup model
model = Model()

# setup variables
x0 = model.add_variable(lb=-2, ub=10, name="x_0", var_type=VarType.CONTINUOUS)
x1 = model.add_variable(lb=-7, ub=12, name="x_1", var_type=VarType.INTEGER)
x2 = model.add_variable(lb=-10, ub=17, name="x_2", var_type=VarType.CONTINUOUS)
x3 = model.add_variable(name="x_3", var_type=VarType.BINARY)
x4 = model.add_variable(name="x_4", var_type=VarType.BINARY)

# add constraints
model.add_constraint(x0 + x1 + x2 + x3 + x4 <= 10, "c0")
model.add_constraint(x0 - x1 + x2 - x3 + x4 >= 5, "c1")

# set objective and sense
objective = 2 * x0 - 3 * x1 + 4 * x2 + 5 * x3 - 6 * x4 + 5
model.set_objective(sense=HybridSolverOptSenses.MAXIMIZE, expr=objective)

# set parameters
model.set_absolute_gap(0.0001)
model.set_time_limit(1000)

# solve the model
status = model.optimize()

# print the status
print("Status:", status)

# get the solution
solution = model.solution
print("Solution:")
for var_name, val in solution.items():
    print(f"  {var_name}: {val}")

# get the objective value
objective_value = model.objective_value

# get the runtime
runtime = model.total_time

Quantagonia’s API also provides the possibility to use Gurobipy-like, DOcplex-like, or Python-MIP like syntax to model MIPs. This example shows how to model a MIP problem with the Model class using Gurobipy-like syntax.

from quantagonia.mip.wrappers.gurobipy import GRB, Model

# setup model
m = Model()

# setup variables
x0 = m.addVar(lb=-2, ub=10, name="x_0", vtype=GRB.CONTINUOUS)
x1 = m.addVar(lb=-7, ub=12, name="x_1", vtype=GRB.INTEGER)
x2 = m.addVar(lb=-10, ub=17, name="x_2", vtype=GRB.CONTINUOUS)
x3 = m.addVar(name="x_3", vtype=GRB.BINARY)
x4 = m.addVar(name="x_4", vtype=GRB.BINARY)

# add constraints
m.addConstr(x0 + x1 + x2 + x3 + x4 <= 10, "c0")
m.addConstr(x0 - x1 + x2 - x3 + x4 >= 5, "c1")

objective = 2 * x0 - 3 * x1 + 4 * x2 + 5 * x3 - 6 * x4 + 5
m.setObjective(expr=objective, sense=GRB.MAXIMIZE)

# set parameters
m.Params.MIPGapAbs = 0.0001
m.Params.TimeLimit = 1000

# solve the model
m.optimize()

# print the status
print("Status:", m.Status)

for v in m.getVars():
    print(f"{v.VarName} {v.X:g}")

print(f"Obj: {m.ObjVal:g}")

This example shows how to model a MIP problem with the Model class using CPLEX-like syntax.

from quantagonia.mip.wrappers.cplex import CplexEnum, Model

# setup model
model = Model()

# setup variables
x0 = model.continuous_var(lb=-2, ub=10, name="x_0")
x1 = model.integer_var(lb=-7, ub=12, name="x_1")
x2 = model.continuous_var(lb=-10, ub=17, name="x_2")
x3 = model.binary_var(name="x_3")
x4 = model.binary_var(name="x_4")

# add constraints
model.linear_constraint(x0 + x1 + x2 + x3 + x4, "le", 10, "c0")
model.linear_constraint(x0 - x1 + x2 - x3 + x4, "ge", 5, "c1")

objective = 2 * x0 - 3 * x1 + 4 * x2 + 5 * x3 - 6 * x4 + 5
model.set_objective(expr=objective, sense=CplexEnum.Maximize)

# set parameters
model.parameters.mip.tolerances.absmipgap.set(0.0001)
model.parameters.timelimit.set(300)

# solve the model
status = model.solve()

# print the status
print("Status:", status)

for v in model.get_variables():
    print(f"{v.name} {v.sol_value:g}")

print(f"Obj: {model.objective_value:g}")

Lastly, it is also possible to solve PuLP models with the HybridSolver.

import math
import os

import pulp

from quantagonia import HybridSolverParameters
from quantagonia.mip import pulp_adapter

# retrieve API key from environment
API_KEY = os.getenv("QUANTAGONIA_API_KEY")


###
# We model with vanilla PuLP code
###

# define MIP problem
prob = pulp.LpProblem("test", pulp.LpMaximize)
x1 = pulp.LpVariable("x1", 0, None)
x2 = pulp.LpVariable("x2", 0, None, pulp.LpInteger)
x3 = pulp.LpVariable("x3", 0, None)
prob += 2 * x1 + 4 * x2 + 3 * x3, "obj"
prob += 3 * x1 + 4 * x2 + 2 * x3 <= 60, "c1"
prob += 2 * x1 + 1 * x2 + 2 * x3 <= 40, "c2"
prob += 1 * x1 + 3 * x2 + 2 * x3 <= 80, "c3"

###
# Quantagonia-specific code
###
params = HybridSolverParameters()

# PuLP requires a solver command to be passed to the problem's solve() method
# We get this command as follows
hybrid_solver_cmd = pulp_adapter.HybridSolver_CMD(API_KEY, params, keepFiles=True)

# solve
prob.solve(solver=hybrid_solver_cmd)

# Each of the variables is printed with it's value
for v in prob.variables():
    print(v.name, "=", v.varValue)

# We retrieve the optimal objective value
print("Optimal objective value = ", pulp.value(prob.objective))

QUBO

This example shows how to model a QUBO problem with the QuboModel class. Have a look at our .qubo format for QUBOs to learn about our file format.

import os

from quantagonia import HybridSolver, HybridSolverParameters
from quantagonia.enums import HybridSolverOptSenses
from quantagonia.qubo import QuboModel

API_KEY = os.environ["QUANTAGONIA_API_KEY"]

# setup model
model = QuboModel()

# setup variables
x0 = model.add_variable("x_0", initial=1)
x1 = model.add_variable("x_1", initial=1)
x2 = model.add_variable("x_2", initial=1)
x3 = model.add_variable("x_3", initial=1)
x4 = model.add_variable("x_4", initial=1)

# build objective
model.objective += 2 * x0
model.objective += 2 * x2
model.objective += 2 * x4
model.objective -= x0 * x2
model.objective -= x2 * x0
model.objective -= x0 * x4
model.objective -= x4 * x0
model.objective -= x2 * x4
model.objective -= x4 * x2
model.objective += 3

# set the sense
model.sense = HybridSolverOptSenses.MINIMIZE

print("Problem: ", model)
print("Initial: ", model.eval())

hybrid_solver = HybridSolver(API_KEY)

params = HybridSolverParameters()
params.set_time_limit(10)
res = model.solve(hybrid_solver, params)

print("Runtime:", res["timing"])
print("Status:", res["sol_status"])
print("Objective:", model.eval())


if model.eval() != 3.0:
    error_message = f"Objective value is not correct: {model.eval()} instead of 3.0"
    raise AssertionError(error_message)

Import QUBO problems

The Quantagonia API client offers possibilities to import QUBOs from Qiskit, D-Wave Ocean, and from PyQUBO.

Qiskit

import os

from qiskit_optimization.converters import QuadraticProgramToQubo
from qiskit_optimization.problems.quadratic_program import QuadraticProgram

from quantagonia import HybridSolver, HybridSolverParameters
from quantagonia.qubo.model import QuboModel

API_KEY = os.environ["QUANTAGONIA_API_KEY"]


def build_qiskit_qp() -> QuadraticProgram:
    """Builds a small subset sum quadratic optimization problem."""
    weights = [9, 7, 2, 1]
    capacity = 10

    qp = QuadraticProgram()

    # create binary variables (y_j = 1 if bin j is used, x_i_j = 1 if item i
    # is packed in bin j)
    for j in range(len(weights)):
        qp.binary_var(f"y_{j}")
    for i in range(len(weights)):
        for j in range(len(weights)):
            qp.binary_var(f"x_{i}_{j}")

    # minimize the number of used bins
    qp.minimize(linear={f"y_{j}": 1 for j in range(len(weights))})

    # ensure that each item is packed in exactly one bin
    for i in range(len(weights)):
        qp.linear_constraint(
            name=f"item_placing_{i}", linear={f"x_{i}_{j}": 1 for j in range(len(weights))}, sense="==", rhs=1
        )

    # ensure that the total weight of items in each bin does not exceed the capacity
    for j in range(len(weights)):
        lhs_x = {f"x_{i}_{j}": weights[i] for i in range(len(weights))}
        lhs_y = {f"y_{j}": -capacity}
        qp.linear_constraint(name=f"capacity_bin_{j}", linear={**lhs_x, **lhs_y}, sense="<=", rhs=0)

    return qp


# build qiskit qp
qp = build_qiskit_qp()
# use qiskit to convert qp to qubo
qp2qubo = QuadraticProgramToQubo()
qiskit_qubo = qp2qubo.convert(qp)

# read the qiskit qubo
qubo = QuboModel.from_qiskit_qubo(qiskit_qubo)

# solve with Quantagonia
params = HybridSolverParameters()
hybrid_solver = HybridSolver(API_KEY)
qubo.solve(hybrid_solver, params)

# to be used in tests
obj = qubo.eval()
if obj != -18.0:
    msg = f"Objective value is not correct: {obj} instead of -18.0"
    raise ValueError(msg)

D-Wave Ocean

In this example, we build the map coloring bqm model from here and import it into a QuboModel to solve it with the HybridSolver.

import os

from dimod import BinaryQuadraticModel
from dimod.generators import combinations

from quantagonia import HybridSolver, HybridSolverParameters
from quantagonia.qubo import QuboModel

# build the map coloring bqm model from here:
# https://github.com/dwavesystems/dwave-ocean-sdk/blob/master/docs/examples/map_coloring_full_code.rst

# Represent the map as the nodes and edges of a graph
provinces = ["AB", "BC", "MB", "NB", "NL", "NS", "NT", "NU", "ON", "PE", "QC", "SK", "YT"]
neighbors = [
    ("AB", "BC"),
    ("AB", "NT"),
    ("AB", "SK"),
    ("BC", "NT"),
    ("BC", "YT"),
    ("MB", "NU"),
    ("MB", "ON"),
    ("MB", "SK"),
    ("NB", "NS"),
    ("NB", "QC"),
    ("NL", "QC"),
    ("NT", "NU"),
    ("NT", "SK"),
    ("NT", "YT"),
    ("ON", "QC"),
]

colors = ["y", "g", "r", "b"]

# Add constraint that each node (province) select a single color
bqm_one_color = BinaryQuadraticModel("BINARY")
for province in provinces:
    variables = [province + "_" + c for c in colors]
    bqm_one_color.update(combinations(variables, 1))

# Add constraint that each pair of nodes with a shared edge not both select one color
bqm_neighbors = BinaryQuadraticModel("BINARY")
for neighbor in neighbors:
    v, u = neighbor
    interactions = [(v + "_" + c, u + "_" + c) for c in colors]
    for interaction in interactions:
        bqm_neighbors.add_quadratic(interaction[0], interaction[1], 1)

bqm = bqm_one_color + bqm_neighbors

# read bqm file into QuboModel
qubo = QuboModel.from_dwave_bqm(bqm)

# solve QUBO with Quantagonia's solver
API_KEY = os.environ["QUANTAGONIA_API_KEY"]
params = HybridSolverParameters()
hybrid_solver = HybridSolver(API_KEY)
qubo.solve(hybrid_solver, params)

# print solution vector
print("Optimal solution vector:")
for var_name, var in qubo.vars.items():
    print("\t", var_name, "\t", var.eval())

# in order to use these as test
obj = qubo.eval()
if obj != 0.0:
    msg = f"Objective value is not correct: {obj} instead of 0.0"
    raise ValueError(msg)

PyQUBO

In this example we model the QUBO

\[2x_0 + 2x_2 + 2x_4 - x_0x_2 - x_2x_0 - x_0x_4 - x_4x_0 - x_2x_4 - x_4x_2 + 3\]

using PyQUBO methods and import it into a QuboModel to solve it with the HybridSolver:

import os

import pyqubo as pq

from quantagonia import HybridSolver, HybridSolverParameters
from quantagonia.qubo import QuboModel

API_KEY = os.environ["QUANTAGONIA_API_KEY"]

x0, x1, x2, x3, x4 = pq.Binary("x0"), pq.Binary("x1"), pq.Binary("x2"), pq.Binary("x3"), pq.Binary("x4")
pyqubo_qubo = (2 * x0 + 2 * x2 + 2 * x4 - x0 * x2 - x2 * x0 - x0 * x4 - x4 * x0 - x2 * x4 - x4 * x2 + 3).compile()

# setup model
qubo = QuboModel.from_pyqubo(pyqubo_qubo)

# solve with Quantagonia
params = HybridSolverParameters()
hybrid_solver = HybridSolver(API_KEY)
qubo.solve(hybrid_solver, params)

# to be used in tests
obj = qubo.eval()
if obj != 3.0:
    msg = f"Objective value is not correct: {obj} instead of 3.0"
    raise ValueError(msg)

Solve existing problems

We first import the necessary modules and define the path to the input file, which can be either a QUBO or a MIP in .qubo, .mps, or, .lp file format.

import os

from quantagonia import HybridSolver, HybridSolverParameters

input_file_path = os.path.join(os.path.dirname(os.path.abspath(__file__)), "data", "example.qubo")

Then, we retrieve the API key stored in an environment variable and set up a HybridSolver object.

api_key = os.environ["QUANTAGONIA_API_KEY"]
hybrid_solver = HybridSolver(api_key)

After that, we set solver parameters and solve the problem with the HybridSolver. This is a blocking call, i.e., the solve function blocks until the problem is fully processed.

params = HybridSolverParameters()
res_dict, _ = hybrid_solver.solve(input_file_path, params)

After the problem has been solved, we retrieve and print some information and the solution vector.

# print some results
print("Runtime:", res_dict["timing"])
print("Objective:", res_dict["objective"])
print("Bound:", res_dict["bound"])
print("Solution:")
for idx, val in res_dict["solution"].items():
    print(f"\t{idx}: {val}")

Submit and check

Instead of using the blocking solve method, you can also use the submit method, to submit problem files to the HybridSolver. The method returns a jobid, which you can use to poll and print the progress.

First, do some setup as before:

import os
from time import sleep

from quantagonia import HybridSolver, HybridSolverParameters
from quantagonia.cloud.enums import JobStatus

input_file_path = os.path.join(os.path.dirname(os.path.abspath(__file__)), "data", "air05.mps.gz")

api_key = os.environ["QUANTAGONIA_API_KEY"]

hybrid_solver = HybridSolver(api_key)
params = HybridSolverParameters()

Then, submit the problem and retrieve a job ID:

# submit the job
jobid = hybrid_solver.submit(input_file_path, params)

We use the job ID to query the progress of the job.

# check progress
has_incumbents = False
while not has_incumbents:
    sleep(2)
    # get intermediate progress information
    p = hybrid_solver.progress(jobid)[0]
    has_incumbents = p["num_incumbents"] >= 1
print("Current status:")
print(f" - Found solutions: {p['num_incumbents']}")
print(f" - Objective Value: {p['objective']}")
print(f" - Best Bound:      {p['bound']}")
print(f" - Relative Gap:    {p['rel_gap']}")

We can also poll the status of the job and print the logs, after the job has been finished:

# check status
while hybrid_solver.status(jobid) != JobStatus.finished:
    print("Wait until job is finished")
    sleep(2)

# get logs
print("\nLogs:")
logs = hybrid_solver.logs(jobid)
print(logs[0])

Solve a batch of problems

It is possible to submit several problems in one batched call. They do not have to be the same type, i.e., you can mix QUBOs and MIPs. Using a batched submission may be beneficial for short-running problems, for which submission times comprise a significant portion of the overall processing time. Not that in the example, the batch contains a garbage file, which cannot be read. Still, the job is processed successfully.

import os

from quantagonia import HybridSolver, HybridSolverParameters

mip_path0 = os.path.join(os.path.dirname(os.path.abspath(__file__)), "data", "example.mps")
mip_path1 = os.path.join(os.path.dirname(os.path.abspath(__file__)), "data", "garbage.mps")
qubo_path = os.path.join(os.path.dirname(os.path.abspath(__file__)), "data", "example.qubo")
api_key = os.environ["QUANTAGONIA_API_KEY"]

hybrid_solver = HybridSolver(api_key)

params = HybridSolverParameters()

problems = [mip_path0, mip_path1, qubo_path]
params = [params, params, params]  # use the same parameters for all three problems
results, _ = hybrid_solver.solve(problems, params, suppress_output=True)  # blocks until all problems are solved

for ix, res in enumerate(results):
    print(f"=== PROBLEM {ix}: status {res['status']} ===")
    print(res["solver_log"])

IP as QUBO

Purely integer problems (IPs) can be reformulated to and solved as QUBOs by enabling the as_qubo parameter.

import math
import os

from quantagonia import HybridSolver, HybridSolverParameters

input_file_path = os.path.join(os.path.dirname(os.path.abspath(__file__)), "data", "ip_as_qubo.lp")

api_key = os.environ["QUANTAGONIA_API_KEY"]

# setup the cloud-based solver instance
hybrid_solver = HybridSolver(api_key)

# set solver parameters
params = HybridSolverParameters()
params.set_as_qubo(True)  # this parameterizes the solver to solve the IP as QUBO

res_dict, _ = hybrid_solver.solve(input_file_path, params)

# print some results
print("Runtime:", res_dict["timing"])
print("Objective:", res_dict["objective"])
print("Bound:", res_dict["bound"])
print("Solution:")
for idx, val in res_dict["solution"].items():
    print(f"\t{idx}: {val}")

# in order to use these as test
obj = res_dict["objective"]
if math.fabs(obj - 4.0) > 1e-4:
    msg = f"Objective value is not correct: {obj} instead of 4.0"
    raise ValueError(msg)