Constraint and Penalty

Constrained Optimization Problem

In mathematical optimization, a constraint is a condition that the solution must satisfy. For example, the following problem is a constrained optimization problem:

$$\begin{aligned} \text{Minimize} & \quad f(x) \\ \text{subject to} & \quad g(x) = 0 \end{aligned}$$

where $f$ and $g$ are functions of the decision variable $x$. The condition $g(x) = 0$ is called an equality constraint. The set of all $x$ that satisfy $g(x) = 0$ is called the feasible set. The constraint may be an equality constraint like $g(x) = 0$ or an inequality constraint like $g(x) \leq 0$. For example, the following problem is also a constrained optimization problem:

$$\begin{aligned} \text{Minimize} & \quad f(x) \\ \text{subject to} & \quad g(x) \leq 0 \end{aligned}$$

In jijmodeling, you can describe both equality and inequality constraints with the Constraint class. For example, an equality constraint $\sum_i x_i = 1$ can be represented as follows:

N = jm.Placeholder("N")
x = jm.BinaryVar("x", shape=(N,))
n = jm.Element('n', belong_to=(0, N))

jm.Constraint("onehot", jm.sum(n, x[n]) == 1)

Please note that in the aforementioned code, there is a string "onehot" as the first argument of jm.Constraint. A constraint object has a name as well as a constraint expression. These names are used during the sampling phase to check whether the constraints are satisfied. The constraint expression must be a boolean expression using only one of the three comparison operators ==, <=, or >=. You can impose multiple constraints on a single problem as follows:

x = jm.BinaryVar("x", shape=(4,))
problem = jm.Problem("constraint_sample")
problem += jm.Constraint("c01", x[0] + x[1] <= 1)
problem += jm.Constraint("c12", x[1] + x[2] <= 1)
problem += jm.Constraint("c23", x[2] + x[3] >= 0)

tip

Other comparison operators, such as >, and logical operators are not supported.

x = jm.BinaryVar("x", shape=(4,))
jm.Constraint("unsupported", (x[0] + x[1] <= 1) | (x[1] + x[2] <= 1))

forall Constraint

Often constraints are indexed by some variables. For example, the following problem is a constrained optimization problem:

$$\begin{aligned} \text{Minimize} & \quad \sum_{i=0}^{N-1} \sum_{j=0}^{M-1} a_{ij} x_{ij} \\ \text{subject to} & \quad \sum_{j = 0}^{M - 1} x_{ij} = 1 \quad \forall i \in \left\{0, \ldots, N - 1\right\} \end{aligned}$$

To represent such a $\forall$ constraint, there is a forall option in the Constraint object. For example, the above problem can be represented as follows:

N = jm.Placeholder("N")
M = jm.Placeholder("M")
a = jm.Placeholder("a", ndim=2)
x = jm.BinaryVar("x", shape=(N, M))
i = jm.Element('i', belong_to=(0, N))
j = jm.Element('j', belong_to=(0, M))

problem = jm.Problem("forall_sample")
problem += jm.sum([i, j], a[i, j] * x[i, j])
problem += jm.Constraint("onehot", jm.sum(j, x[i, j]) == 1, forall=i)

What is Penalty?

The Penalty method and the Lagrange multiplier method are the most popular methods for converting constrained optimization problems into unconstrained ones. In the penalty method, a constraint optimization problem is represented as follows:

$$\begin{aligned} \text{Minimize} & \quad f(x) \\ \text{subject to} & \quad g(x) = 0 \end{aligned}$$

This problem is transformed into an unconstrained optimization problem of the form:

$$\text{Minimize} \quad f(x) + \alpha p(x),$$

where $\alpha$ (penalty coefficient or Lagrangian multiplier) and $p(x)$ (penalty term) play crucial roles. Typically, $p(x)$ is chosen as $p(x) = g(x)^2$. If the minimum of $f(x) + \alpha p(x)$ is found and satisfies $p(x) = 0$, then the minimum of the original constrained optimization problem is found. If the penalty $p(x)$ is positive, increase the value of the penalty coefficient $\alpha$ and rerun the optimization problem.

It's important to note that some solvers only accept unconstrained problems. The "U" in QUBO stands for "Unconstrained." To solve constrained problems defined in jijmodeling using these solvers, they are translated into unconstrained problems by jijmodeling-transpiler or the transpiler layer of JijZept with the default penalty term $p$. However, the choice of $p$ is critical for real-world optimization problems, and jijmodeling provides a way to customize the penalty term using CustomPenaltyTerm.

Translate Constraint to Penalty

The translation of constraints into penalty terms is handled by transpilers, not jijmodeling itself. Here, we explain how jijmodeling-transpiler translates constraints into penalty terms. Let's consider a small problem:

$$\begin{aligned} \text{Minimize} & \quad \sum_{i=0}^{N-1} a_i x_i \\ \text{subject to} & \quad \sum_{i = 0}^{N - 1} x_i = 1 \end{aligned}$$

to see the behavior of jijmodeling-transpiler:

import jijmodeling as jm

a = jm.Placeholder("a", ndim=1)
N = a.len_at(0, latex="N")

i = jm.Element('i', belong_to=(0, N))
x = jm.BinaryVar('x', shape=(N,))

problem = jm.Problem('translate_constraint_to_penalty')
problem += jm.sum(i, a[i]*x[i])
problem += jm.Constraint(
    'onehot',
    jm.sum(i, x[i]) == 1,
)
problem

This problem is converted into an unconstrained problem:

$$\text{Minimize} \quad \sum_{i=0}^{N-1} a_i x_i + \alpha \left(\sum_{i = 0}^{N - 1} x_i - 1\right)^2$$

with $a = [1, 2]$ and $\alpha = 5$ as follows:

import jijmodeling_transpiler as jmt

instance_data = {
    "a": [1, 2],
}

compiled_model = jmt.core.compile_model(problem, instance_data)
pubo_builder = jmt.core.pubo.transpile_to_pubo(
    compiled_model,
    normalize=False  # Disable normalization for simplicity
)
qubo, constant = pubo_builder.get_qubo_dict(multipliers={ 'onehot': 5 })

This yields:

qubo     = {(0, 0): -4.0, (1, 1): -3.0, (0, 1): 10.0}
constant = 5.0

because:

$$\begin{aligned} \sum_{i=0}^{N-1} a_i x_i + \alpha \left(\sum_{i = 0}^{N - 1} x_i - 1\right)^2 &= x_1 + 2 x_2 + 5 (x_1 + x_2 - 1)^2 \\ &= -4 x_1 - 3 x_2 + 10 x_1 x_2 + 5 \\ &= \begin{bmatrix} x_1 & x_2 \end{bmatrix} \begin{bmatrix} -4 & 10 \\ 0 & -3 \end{bmatrix} \begin{bmatrix} x_1 \\ x_2 \end{bmatrix} + 5 \end{aligned}$$

Be sure that $x_i^2 = x_i$ for binary variables.

This instantiation phase is separated into two phases where:

• Convert the Problem object into a CompiledInstance object with instance_data.
• Convert the CompiledInstance object into a QUBO with multipliers.

Since the CompiledInstance represents "what to be solved," and the actual QUBO coefficients with the multiplier represent "how to be solved."

info

jijmodeling-transpiler provides a way to set multipliers for each penalty term using the detail_parameters option of the get_qubo_dict method. Another relaxation method like the Augmented Lagrangian method is also supported. Please see the reference of jijmodeling-transpiler for details.

CustomPenaltyTerm

Although translating constraints into penalty terms is the responsibility of transpilers, one may want to use an original penalty term without creating another transpiler from scratch. jijmodeling provides a way to customize the penalty term using CustomPenaltyTerm. Here, we explain how to use CustomPenaltyTerm to define the same penalty term as in the previous example:

$$\text{Minimize} \quad \sum_{i=0}^{N-1} a_i x_i + \alpha \left(\sum_{i = 0}^{N - 1} x_i - 1\right)^2$$

This problem can be represented using CustomPenaltyTerm as follows:

import jijmodeling as jm

a = jm.Placeholder("a", ndim=1)
N = a.len_at(0, latex="N")

i = jm.Element('i', belong_to=(0, N))
x = jm.BinaryVar('x', shape=(N,))

problem = jm.Problem('penalty_sample')
problem += jm.sum(i, a[i]*x[i])
problem += jm.CustomPenaltyTerm(
    'onehot',
    (jm.sum(i, x[i]) - 1)**2,
)

Note that the multiplier $\alpha$ does not appear here. This can be instantiated in the exact same way as the previous example:

import jijmodeling_transpiler as jmt

instance_data = {
    "a": [1, 2],
}

compiled_model = jmt.core.compile_model(problem, instance_data)
pubo_builder = jmt.core.pubo.transpile_to_pubo(
    compiled_model,
    normalize=False  # Disable normalization for simplicity
)
qubo, constant = pubo_builder.get_qubo_dict(multipliers={ 'onehot': 5 })

This yields the same result:

qubo     = {(0, 0): -4.0, (1, 1): -3.0, (0, 1): 10.0}
constant = 5.0

There is also a forall option for CustomPenaltyTerm:

import jijmodeling as jm

a = jm.Placeholder("a", ndim=2)
N = a.len_at(0, latex="N")
M = a.len_at(1, latex="M")

i = jm.Element('i', belong_to=(0, N))
j = jm.Element('j', belong_to=(0, M))
x = jm.BinaryVar('x', shape=(N, M))

problem = jm.Problem('forall_penalty_sample')
problem += jm.sum([i, j], a[i, j]*x[i, j])
problem += jm.CustomPenaltyTerm(
    'onehot',
    (jm.sum(j, x[i, j]) - 1)**2,
    forall=i
)
problem

$$\begin{array}{cccc} \text{Problem:} & \text{forall\_penalty\_sample} & & \\ & & \min \quad \displaystyle \sum_{i = 0}^{N - 1} \sum_{j = 0}^{M - 1} a_{i, j} \cdot x_{i, j} & \\ \text{{penalty terms}} & & & \\ & \text{onehot} & \displaystyle \left(\left(\sum_{j = 0}^{M - 1} x_{i, j} - 1\right)^{2}\right) & \forall i \in \left\{0,\ldots,N - 1\right\} \\ \text{{where}} & & & \\ & x & 2\text{-dim binary variable}\\ \end{array}$$