Model quantum optimization problems with the `OptimizationProblem` class¶

Introduction¶

This is the initial guide to show how to use the Qiskit Optimization Mapper Addon for quantum optimization modeling. The addon introduces the OptimizationProblem class, which provides a flexible way to handle variables, constraints and objective functions.

The OptimizationProblem class supports higher-order polynomial objectives and constraints, as well as spin variables, unlike the QuadraticProgram class from the legacy qiskit-optimizationpackage objectives and constraints. With OptimizationProblem, you can work with a broader class of optimization models.

OptimizationProblem handles constrained polynomial programs, which you can express in the following general form:

\[\begin{split}\begin{align} \text{minimize}\quad& \sum_{d=1}^{D} \;\; \sum_{S \subseteq \{1,\dots,n\}, |S|=d} c_{S} \prod_{i \in S} x_i \\ \text{subject to}\quad& \sum_{d=1}^{D} \;\; \sum_{S \subseteq \{1,\dots,n\}, |S|=d} a^{(j)}_{S} \prod_{i \in S} x_i \;\; \leq b_j, \quad j=1,\dots,m \\ & l_i \leq x_i \leq u_i, \quad i=1,\dots,n, \end{align}\end{split}\]

where the coefficients \(c_S\) and \(a^{(j)}_S\) may involve terms of any degree up to \(D\). The variables \(x_i\) can be defined as binary, integer, continuous, or spin. In addition to “\(\leq\)” constraints, the OptimizationProblem model also supports “\(\geq\)” and “\(=\)”.

Construct an `OptimizationProblem`¶

This section explains how to build an optimization model by using the OptimizationProblem class. Start by creating an empty model. The model does not yet contain variables, an objective, or constraints. Note that OptimizationProblem has a prettyprint method to generate a comprehensive string representation of the model.

[1]:

from qiskit_addon_opt_mapper import OptimizationProblem

# make an empty problem
mod = OptimizationProblem("my problem")
print(mod.prettyprint())

Problem name: my problem

Minimize
  0

Subject to
  No constraints

  No variables

The OptimizationProblem class supports four types of variables:

Binary variables can have the value 0 or 1.
Integer variables can have any integer value within their bounds.
Continuous variables can have any real value within their bounds.
Spin variables can have the value -1 or 1.

When you define a variable, you specify its name and type. You can also set optional lower and upper bounds.

[2]:

# Add variables
mod.binary_var(name="x")
mod.integer_var(name="y", lowerbound=-1, upperbound=5)
mod.continuous_var(name="z", lowerbound=-1, upperbound=5)
mod.spin_var(name="s")
print(mod.prettyprint())

Problem name: my problem

Minimize
  0

Subject to
  No constraints

  Integer variables (1)
    -1 <= y <= 5

  Continuous variables (1)
    -1 <= z <= 5

  Binary variables (1)
    x

  Spin variables (1)
    -1 <= s <= 1

You set the objective function by calling OptimizationProblem.minimize or OptimizationProblem.maximize.

The objective function can include the following components:

A constant offset.
Linear terms ((c^{\top} x)).
Quadratic terms ((x^{\top} Q x)).
Higher-order terms (\sumS h_S :nbsphinx-math:`prod`{i \in `S} x_i), where (\|S\| :nbsphinx-math:ge 3`).

You can specify these terms by using a list, a matrix, or a dictionary.

The following example shows how to define an objective function by using a dictionary:

For linear terms, the keys are variable names and the values are their coefficients.
For quadratic terms, the keys are pairs of variables and the values are their coefficients.
For higher-order terms, use a nested dictionary with the following structure: Dict[int, Dict[Tuple[str | int, ...], float]]. Where the outer key is the degree of the term and the inner dictionary maps a tuple of variable names (or indices) to its coefficient.

[3]:

# Add objective function using dictionaries
mod.minimize(
    constant=3,
    linear={"x": 1},
    quadratic={("x", "y"): 2, ("z", "z"): -1},
    higher_order={3: {("x", "y", "z"): 4}},
)
print(mod.prettyprint())

Problem name: my problem

Minimize
  4*x*y*z + 2*x*y - z^2 + x + 3

Subject to
  No constraints

  Integer variables (1)
    -1 <= y <= 5

  Continuous variables (1)
    -1 <= z <= 5

  Binary variables (1)
    x

  Spin variables (1)
    -1 <= s <= 1

You can also specify the objective function by using arrays:

For linear terms, use an array that corresponds to the vector \(c\) in the mathematical formulation.
For quadratic terms, use an array that corresponds to the matrix \(Q\).
For higher-order terms, use multi-dimensional arrays. An order-\(k\) tensor represents the coefficients of degree-\(k\) monomials. For example, a third-order tensor \(H^{(3)}\) stores the coefficients \(h_{ijk}\) for terms \(x_i x_j x_k\).

[4]:

import numpy as np

# Add objective function using lists/arrays
# Here, h3 is a 3D array for the cubic terms
h3 = np.zeros((4, 4, 4))  # assuming we have 4 variables
h3[0, 1, 2] = 4  # corresponds to x * y * z
mod.minimize(
    constant=3,
    linear=[1, 0, 0, 0],
    quadratic=[[0, 1, 0, 0], [1, 0, 0, 0], [0, 0, -1, 0], [0, 0, 0, 0]],
    higher_order={3: h3},
)
print(mod.prettyprint())

Problem name: my problem

Minimize
  4*x*y*z + 2*x*y - z^2 + x + 3

Subject to
  No constraints

  Integer variables (1)
    -1 <= y <= 5

  Continuous variables (1)
    -1 <= z <= 5

  Binary variables (1)
    x

  Spin variables (1)
    -1 <= s <= 1

You can access the constant, the linear term, the quadratic term, and the higher-order terms by looking at OptimizationProblem.objective.{constant, linear, quadratic, higher_order}, respectively.

For linear, quadratic, and higher-order terms, you can obtain:

a dense array (to_array),
a sparse representation (coefficients),
or a dictionary (to_dict).

When using dictionaries, you can specify whether to use variable indices or names as keys.

Quadratic terms are stored in a compressed way. For example, {('x', 'y'): 1, ('y', 'x'): 2} is stored as {('x', 'y'): 3}. You can get the quadratic term as a symmetric matrix by calling to_array(symmetric=True).
Higher-order terms are represented as nested dictionaries: Dict[int, Dict[Tuple[str | int, ...], float]], where the outer key is the degree, and the inner dictionary maps a tuple of variables to its coefficient. For example, {3: {('x','y','z'): 4}} corresponds to the cubic term \(4xyz\). Similarly to quadratic terms, you can call to_array(symmetric=True) to aggregate coefficients across permutations of the same monomial.

[5]:

print("constant:\t\t\t", mod.objective.constant)
print("linear dict:\t\t\t", mod.objective.linear.to_dict())
print("linear array:\t\t\t", mod.objective.linear.to_array())
print("linear array as sparse matrix:\n", mod.objective.linear.coefficients, "\n")
print("quadratic dict w/ index:\t", mod.objective.quadratic.to_dict())
print("quadratic dict w/ name:\t\t", mod.objective.quadratic.to_dict(use_name=True))
print(
    "symmetric quadratic dict w/ name:\t",
    mod.objective.quadratic.to_dict(use_name=True, symmetric=True),
)
print("quadratic matrix:\n", mod.objective.quadratic.to_array(), "\n")
print("symmetric quadratic matrix:\n", mod.objective.quadratic.to_array(symmetric=True), "\n")
print("quadratic matrix as sparse matrix:\n", mod.objective.quadratic.coefficients)
print("higher-order dict w/ index:\t", mod.objective.higher_order[3].to_dict())
print("higher-order dict w/ name:\t", mod.objective.higher_order[3].to_dict(use_name=True))
print("higher-order matrix:\n", mod.objective.higher_order[3].to_array(), "\n")
print(
    "symmetric higher-order matrix:\n", mod.objective.higher_order[3].to_array(symmetric=True), "\n"
)
print("higher-order matrix as sparse matrix:\n", mod.objective.higher_order[3].coefficients)

constant:                        3.0
linear dict:                     {np.int32(0): np.int64(1)}
linear array:                    [1 0 0 0]
linear array as sparse matrix:
   (np.int32(0), np.int32(0))   1

quadratic dict w/ index:         {(0, 1): np.float64(2.0), (2, 2): np.float64(-1.0)}
quadratic dict w/ name:          {('x', 'y'): np.float64(2.0), ('z', 'z'): np.float64(-1.0)}
symmetric quadratic dict w/ name:        {('x', 'y'): np.float64(1.0), ('y', 'x'): np.float64(1.0), ('z', 'z'): np.float64(-1.0)}
quadratic matrix:
 [[ 0.  2.  0.  0.]
 [ 0.  0.  0.  0.]
 [ 0.  0. -1.  0.]
 [ 0.  0.  0.  0.]]

symmetric quadratic matrix:
 [[ 0.  1.  0.  0.]
 [ 1.  0.  0.  0.]
 [ 0.  0. -1.  0.]
 [ 0.  0.  0.  0.]]

quadratic matrix as sparse matrix:
   (np.int32(0), np.int32(1))   2.0
  (np.int32(2), np.int32(2))    -1.0
higher-order dict w/ index:      {(0, 1, 2): 4.0}
higher-order dict w/ name:       {('x', 'y', 'z'): 4.0}
higher-order matrix:
 [[[0.         0.         0.         0.        ]
  [0.         0.         0.66666667 0.        ]
  [0.         0.66666667 0.         0.        ]
  [0.         0.         0.         0.        ]]

 [[0.         0.         0.66666667 0.        ]
  [0.         0.         0.         0.        ]
  [0.66666667 0.         0.         0.        ]
  [0.         0.         0.         0.        ]]

 [[0.         0.66666667 0.         0.        ]
  [0.66666667 0.         0.         0.        ]
  [0.         0.         0.         0.        ]
  [0.         0.         0.         0.        ]]

 [[0.         0.         0.         0.        ]
  [0.         0.         0.         0.        ]
  [0.         0.         0.         0.        ]
  [0.         0.         0.         0.        ]]]

symmetric higher-order matrix:
 [[[0. 0. 0. 0.]
  [0. 0. 4. 0.]
  [0. 0. 0. 0.]
  [0. 0. 0. 0.]]

 [[0. 0. 0. 0.]
  [0. 0. 0. 0.]
  [0. 0. 0. 0.]
  [0. 0. 0. 0.]]

 [[0. 0. 0. 0.]
  [0. 0. 0. 0.]
  [0. 0. 0. 0.]
  [0. 0. 0. 0.]]

 [[0. 0. 0. 0.]
  [0. 0. 0. 0.]
  [0. 0. 0. 0.]
  [0. 0. 0. 0.]]]

higher-order matrix as sparse matrix:
 {(0, 1, 2): 4.0}

Add or remove linear, quadratic and higher-order terms to and from the constraints¶

You can add linear constraints by setting name, linear expression, sense and right-hand-side value (rhs). You can use senses '==', '<=', and '>='.

[6]:

# Add linear constraints
mod.linear_constraint(linear={"x": 1, "y": 2}, sense="==", rhs=3, name="lin_eq")
mod.linear_constraint(linear={"x": 1, "y": 2}, sense="<=", rhs=3, name="lin_leq")
mod.linear_constraint(linear={"x": 1, "y": 2}, sense=">=", rhs=3, name="lin_geq")
print(mod.prettyprint())

Problem name: my problem

Minimize
  4*x*y*z + 2*x*y - z^2 + x + 3

Subject to
  Linear constraints (3)
    x + 2*y == 3  'lin_eq'
    x + 2*y <= 3  'lin_leq'
    x + 2*y >= 3  'lin_geq'

  Integer variables (1)
    -1 <= y <= 5

  Continuous variables (1)
    -1 <= z <= 5

  Binary variables (1)
    x

  Spin variables (1)
    -1 <= s <= 1

You can add quadratic constraints by setting name, quadratic expression, sense and right-hand-side value (rhs).

[7]:

# Add quadratic constraints
mod.quadratic_constraint(
    linear={"x": 1, "y": 1},
    quadratic={("x", "x"): 1, ("y", "z"): -1},
    sense="==",
    rhs=1,
    name="quad_eq",
)
mod.quadratic_constraint(
    linear={"x": 1, "y": 1},
    quadratic={("x", "x"): 1, ("y", "z"): -1},
    sense="<=",
    rhs=1,
    name="quad_leq",
)
mod.quadratic_constraint(
    linear={"x": 1, "y": 1},
    quadratic={("x", "x"): 1, ("y", "z"): -1},
    sense=">=",
    rhs=1,
    name="quad_geq",
)
print(mod.prettyprint())

Problem name: my problem

Minimize
  4*x*y*z + 2*x*y - z^2 + x + 3

Subject to
  Linear constraints (3)
    x + 2*y == 3  'lin_eq'
    x + 2*y <= 3  'lin_leq'
    x + 2*y >= 3  'lin_geq'

  Quadratic constraints (3)
    x^2 - y*z + x + y == 1  'quad_eq'
    x^2 - y*z + x + y <= 1  'quad_leq'
    x^2 - y*z + x + y >= 1  'quad_geq'

  Integer variables (1)
    -1 <= y <= 5

  Continuous variables (1)
    -1 <= z <= 5

  Binary variables (1)
    x

  Spin variables (1)
    -1 <= s <= 1

You can also add higher-order constraints by setting name, higher-order expression, sense and right-hand-side value (rhs).

[8]:

# Add higher-order constraints
mod.higher_order_constraint(
    linear={"x": 1},
    quadratic={("x", "x"): 1, ("y", "z"): -1},
    higher_order={3: {("x", "y", "z"): 1}, 4: {("x", "y", "z", "s"): 2}},
    sense="==",
    rhs=1,
    name="hoc_eq",
)
mod.higher_order_constraint(
    linear={"x": 1},
    quadratic={("x", "x"): 1, ("y", "z"): -1},
    higher_order={3: {("x", "y", "z"): 1}, 4: {("x", "y", "z", "s"): 2}},
    sense="<=",
    rhs=1,
    name="hoc_leq",
)
mod.higher_order_constraint(
    linear={"x": 1},
    quadratic={("x", "x"): 1, ("y", "z"): -1},
    higher_order={3: {("x", "y", "z"): 1}, 4: {("x", "y", "z", "s"): 2}},
    sense=">=",
    rhs=1,
    name="hoc_geq",
)
print(mod.prettyprint())

Problem name: my problem

Minimize
  4*x*y*z + 2*x*y - z^2 + x + 3

Subject to
  Linear constraints (3)
    x + 2*y == 3  'lin_eq'
    x + 2*y <= 3  'lin_leq'
    x + 2*y >= 3  'lin_geq'

  Quadratic constraints (3)
    x^2 - y*z + x + y == 1  'quad_eq'
    x^2 - y*z + x + y <= 1  'quad_leq'
    x^2 - y*z + x + y >= 1  'quad_geq'

  Higher-order constraints (3)
    2*s*x*y*z + x*y*z + x^2 - y*z + x == 1  'hoc_eq'
    2*s*x*y*z + x*y*z + x^2 - y*z + x <= 1  'hoc_leq'
    2*s*x*y*z + x*y*z + x^2 - y*z + x >= 1  'hoc_geq'

  Integer variables (1)
    -1 <= y <= 5

  Continuous variables (1)
    -1 <= z <= 5

  Binary variables (1)
    x

  Spin variables (1)
    -1 <= s <= 1

You can access linear and quadratic terms of linear and quadratic constraints as in the same way as the objective function.

[9]:

lin_geq = mod.get_linear_constraint("lin_geq")
print("lin_geq:", lin_geq.linear.to_dict(use_name=True), lin_geq.sense, lin_geq.rhs)
quad_geq = mod.get_quadratic_constraint("quad_geq")
print(
    "quad_geq:",
    quad_geq.linear.to_dict(use_name=True),
    quad_geq.quadratic.to_dict(use_name=True),
    quad_geq.sense,
    lin_geq.rhs,
)
hoc_geq = mod.get_higher_order_constraint("hoc_geq")
print(
    "hoc_geq:",
    hoc_geq.linear.to_dict(use_name=True),
    hoc_geq.quadratic.to_dict(use_name=True),
    {k: v.to_dict(use_name=True) for k, v in hoc_geq.higher_order.items()},
    hoc_geq.sense,
    lin_geq.rhs,
)

lin_geq: {'x': np.float64(1.0), 'y': np.float64(2.0)} ConstraintSense.GE 3
quad_geq: {'x': np.float64(1.0), 'y': np.float64(1.0)} {('x', 'x'): np.float64(1.0), ('y', 'z'): np.float64(-1.0)} ConstraintSense.GE 3
hoc_geq: {'x': np.float64(1.0)} {('x', 'x'): np.float64(1.0), ('y', 'z'): np.float64(-1.0)} {3: {('x', 'y', 'z'): 1.0}, 4: {('x', 'y', 'z', 's'): 2.0}} ConstraintSense.GE 3

You can also remove linear/quadratic/higher-order constraints by using the remove_linear_constraint, remove_quadratic_constraint, and remove_higher_order_constraint methods, respectively.

[10]:

# Remove constraints
mod.remove_linear_constraint("lin_eq")
mod.remove_quadratic_constraint("quad_leq")
mod.remove_higher_order_constraint("hoc_eq")
print(mod.prettyprint())

Problem name: my problem

Minimize
  4*x*y*z + 2*x*y - z^2 + x + 3

Subject to
  Linear constraints (2)
    x + 2*y <= 3  'lin_leq'
    x + 2*y >= 3  'lin_geq'

  Quadratic constraints (2)
    x^2 - y*z + x + y == 1  'quad_eq'
    x^2 - y*z + x + y >= 1  'quad_geq'

  Higher-order constraints (2)
    2*s*x*y*z + x*y*z + x^2 - y*z + x <= 1  'hoc_leq'
    2*s*x*y*z + x*y*z + x^2 - y*z + x >= 1  'hoc_geq'

  Integer variables (1)
    -1 <= y <= 5

  Continuous variables (1)
    -1 <= z <= 5

  Binary variables (1)
    x

  Spin variables (1)
    -1 <= s <= 1

Substitute Variables¶

You can substitute some of the variables with constants or other variables. More precisely, OptimizationProblem has a substitute_variables(constants=..., variables=...) method to deal with the following two cases.

\(x \leftarrow c\): when constants have a dictionary {x: c}.
\(x \leftarrow c y\): when variables have a dictionary {x: (y, c)}.

[11]:

sub = mod.substitute_variables(constants={"x": 1}, variables={"y": ("z", -1)})
print(sub.prettyprint())

Problem name: my problem

Minimize
  -5*z^2 - 2*z + 4

Subject to
  Linear constraints (2)
    -2*z <= 2  'lin_leq'
    -2*z >= 2  'lin_geq'

  Quadratic constraints (2)
    z^2 - z == -1  'quad_eq'
    z^2 - z >= -1  'quad_geq'

  Higher-order constraints (2)
    -2*s*z^2 <= -1  'hoc_leq'
    -2*s*z^2 >= -1  'hoc_geq'

  Continuous variables (1)
    -1 <= z <= 1

  Spin variables (1)
    -1 <= s <= 1

If the resulting problem is infeasible due to lower bounds or upper bounds, the methods returns the status Status.INFEASIBLE. For example, if with the model above you try to replace variable x with -1, but -1 is out of range of x (0 <= x <= 1), it will return Status.INFEASIBLE.

[12]:

sub = mod.substitute_variables(constants={"x": -1})
print(sub.status)

Infeasible substitution for variable: x

OptimizationProblemStatus.INFEASIBLE

You cannot substitute variables multiple times. The method raises an error in such cases.

[13]:

from qiskit_addon_opt_mapper import OptimizationError

try:
    sub = mod.substitute_variables(constants={"x": -1}, variables={"y": ("x", 1)})
except OptimizationError as e:
    print(f"Error: {e}")

Error: Cannot substitute by variable that gets substituted itself: y <- x 1

Model quantum optimization problems with the OptimizationProblem class¶

Introduction¶

Construct an OptimizationProblem¶

Add or remove linear, quadratic and higher-order terms to and from the constraints¶

Substitute Variables¶

Model quantum optimization problems with the `OptimizationProblem` class¶

Construct an `OptimizationProblem`¶