Optimization Problem¶

Introduction¶

In this tutorial, we introduce how to build optimization problems using the Optimization Add-on module. Optimization Add-on introduces the OptimizationProblem class to model optimization problems. Unlike the traditional QuadraticProgram which was limited to quadratic objectives and constraints, OptimizationProblem supports higher-order polynomial objectives and constraints. More precisely, it deals with constrained polynomial programs given as follows:

\[\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, OptimizationProblem also supports “\(\geq\)” and “\(=\)”.

Constructing an `OptimizationProblem`¶

We now explain how to build a model of an optimization problem using the OptimizationProblem class. Let’s start from an empty 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

OptimizationProblem has a method prettyprint to generate a comprehensive string representation.

The OptimizationProblem 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 value within their bounds
Spin variables can have the value -1 or 1.

When defining variables, you can specify their name, type, and optionally set 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 can set the objective function by invoking OptimizationProblem.minimize or OptimizationProblem.maximize. The objective function can include:

a constant (offset),
linear terms (\(c^{\top} x\)),
quadratic terms (\(x^{\top} Q x\)),
and higher-order terms (\(\sum_S h_S \prod_{i \in S} x_i\), where \(|S| \ge 3\)).

These terms can be specified using a list, a matrix, or a dictionary. For higher-order terms, we use a nested dictionary of the form 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.

The example below shows how to define an objective function using a dictionary. For the linear term, the keys are variable names and the values are their coefficients. For the quadratic term, the keys are pairs of variables, and the values are their coefficients. For higher-order terms, the outer key is the degree of the term, and the inner dictionary maps a tuple of variable names 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

Another way to specify the optimization problem is using arrays.

For the linear term, the array corresponds to the vector \(c\) in the mathematical formulation.
For the quadratic term, the array corresponds to the matrix \(Q\).
For higher-order terms, multi-dimensional arrays are used: an order-\(k\) tensor represents the coefficients of degree-\(k\) monomials. For example, a 3rd-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).

For 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}

Adding/removing linear, quadratic and higher-order terms to/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

Substituting Variables¶

You can substitute some of the variables with constants or other variables. More precisely, OptimizationProblem has a method substitute_variables(constants=..., variables=...) 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. We try to replace variable x with -1, but -1 is out of range of x (0 <= x <= 1). So, it returns 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 a case.

[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

Supporting Docplex models¶

The OptimizationProblem can be constructed from a Docplex model. You can use the from_docplex method to convert a Docplex model to an OptimizationProblem.

[14]:

from docplex.mp.model import Model

docplex_mod = Model(name="docplex model")
x = docplex_mod.binary_var(name="x")
y = docplex_mod.integer_var(name="y", lb=-1, ub=5)
z = docplex_mod.continuous_var(name="z", lb=-1, ub=5)

docplex_mod.minimize(3 + x + 2 * x * y - z * z)
docplex_mod.add_constraint(x + 2 * y == 3, "lin_eq")
docplex_mod.add_constraint(x + 2 * y <= 3, "lin_leq")
docplex_mod.add_constraint(x + 2 * y >= 3, "lin_geq")

print(docplex_mod.export_as_lp_string())

\ This file has been generated by DOcplex
\ ENCODING=ISO-8859-1
\Problem name: docplex model

Minimize
 obj: x + [ 4 x*y - 2 z^2 ]/2 + 3
Subject To
 lin_eq: x + 2 y = 3
 lin_leq: x + 2 y <= 3
 lin_geq: x + 2 y >= 3

Bounds
 0 <= x <= 1
 -1 <= y <= 5
 -1 <= z <= 5

Binaries
 x

Generals
 y
End

Note: When you display your problem as LP format using export_as_lp_string, Binaries denotes binary variables and Generals denotes integer variables. If variables are not included in either Binaries or Generals, such variables are continuous ones with default lower bound = 0 and upper bound = infinity. Note that you cannot use ‘e’ or ‘E’ as the first character of names due to the specification of LP format.

[15]:

from qiskit_addon_opt_mapper.translators import from_docplex_mp

mod2 = from_docplex_mp(docplex_mod)
print(mod2.prettyprint())

Problem name: docplex model

Minimize
  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

You can also convert a OptimizationProblem to a Docplex model using the to_docplex method.

[16]:

from qiskit_addon_opt_mapper.translators import to_docplex_mp

docplex_mod2 = to_docplex_mp(mod2)
print(docplex_mod2.export_as_lp_string())

\ This file has been generated by DOcplex
\ ENCODING=ISO-8859-1
\Problem name: docplex model

Minimize
 obj: x + [ 4 x*y - 2 z^2 ]/2 + 3
Subject To
 lin_eq: x + 2 y = 3
 lin_leq: x + 2 y <= 3
 lin_geq: x + 2 y >= 3

Bounds
 0 <= x <= 1
 -1 <= y <= 5
 -1 <= z <= 5

Binaries
 x

Generals
 y
End