Constraints and Penalty Functions

The Amplify SDK allows you to create models with any variable and polynomial degree constraints. However, each combinatorial optimization solver handles different types and orders of variables as constraints, and in particular, some QUBO solvers do not accept constraints themselves.

The Amplify SDK converts the variables, the polynomial degree of the constraints, and the objective function. It builds an intermediate model that the combinatorial optimization solver accepts. If a conversion is impossible, the SDK generates a penalty function for each original constraint and adds that function to the objective function. It then calls the solver. A solver that cannot accept a constraint directly can therefore solve a model that holds one.

Constraints in intermediate models

Each solver client defines the types of variables and degrees of equality and inequality constraints, as well as the objective function the solver can handle.

Here are examples of input model conversions involving equality constraints on integer variables to an intermediate model.

from amplify import (
    VariableGenerator,
    Model,
    AmplifyAEClient,
    AcceptableDegrees,
    equal_to,
)

gen = VariableGenerator()
n = gen.scalar("Integer", bounds=(-10, 10))  # Issue integer variable
c = equal_to(n, 1)  # Create a constraint, n = 1

model = Model(c)

For example, take a solver client that handles an equality constraint of first-order binary variables directly. The SDK then transforms the input model into the intermediate model as follows.

bqbl = AcceptableDegrees(
    objective={"Binary": "Quadratic"},
    equality_constraints={"Binary": "Linear"}
)
im, mapping = model.to_intermediate_model(bqbl)

>>> print(im)
minimize:
  0
subject to:
  q_0 + 2 q_1 + 4 q_2 + 8 q_3 + 5 q_4 - 10 == 1 (weight: 1)
>>> print(mapping[n])
q_0 + 2 q_1 + 4 q_2 + 8 q_3 + 5 q_4 - 10

The constraints of the intermediate model hold binary variables only, as the variable conversion in the objective function does. The Amplify SDK converts between the integer variables and the binary variables inside the intermediate model.

Similarly, for solvers that cannot handle constraints at all, you can still solve problems involving constraints with the Amplify SDK. The following example shows how to transform an input model with integer variables into an intermediate model for the QUBO solver.

bq = AcceptableDegrees(objective={"Binary": "Quadratic"})
im, mapping = model.to_intermediate_model(bq)

>>> print(im)
minimize:
  0
subject to:
  q_0 + 2 q_1 + 4 q_2 + 8 q_3 + 5 q_4 - 10 == 1 (weight: 1)
>>> print(mapping[n])
q_0 + 2 q_1 + 4 q_2 + 8 q_3 + 5 q_4 - 10

At first glance, there appears to be no difference from the previous example. Still, a polynomial called the penalty function, which we will explain in the next section, is generated and variable conversions are performed. Because the penalty function is part of the objective function, it undergoes the same variable conversion and degree reduction as the objective function. You can see this process as follows.

>>> # Penalty function for the constraint of the input model
>>> print(model.constraints[0].penalty)
n_0^2 - 2 n_0 + 1

>>> # Penalty function for the constraint of the intermediate model
>>> print(im.constraints[0].penalty)
4 q_0 q_1 + 8 q_0 q_2 + 16 q_0 q_3 + 10 q_0 q_4 + 16 q_1 q_2 + 32 q_1 q_3
  + 20 q_1 q_4 + 64 q_2 q_3 + 40 q_2 q_4 + 80 q_3 q_4 - 21 q_0 - 40 q_1
  - 72 q_2 - 112 q_3 - 85 q_4 + 121

Penalty method

The Amplify SDK attempts to convert a constraint expression by variable conversion or by degree reduction. The type and the degree of the constraints that the solver handles decide which one it uses. If no conversion is possible, the SDK assigns a polynomial to the objective function in place of the constraint expression. That polynomial is a penalty function. That penalty function is equivalent to the constraint, so the constraint holds indirectly. While the SDK constructs the intermediate model, it applies variable conversion and degree reduction to the penalty function as well. It does so when they are necessary.

Penalty method

For a given constraint \(c\) on variables \(x_1, x_1, \ldots, x_n\), if a real-valued function \(p\) satisfies:

\[\begin{split} p(x_1, x_2, \ldots, x_n) \begin{cases} = 0 \quad & \text{if } c \text{ is satisfied} \\ > 0 \quad & \text{otherwise} \end{cases}, \end{split}\]

then \(p\) is called the penalty function.

The penalty method is a method of transforming a combinatorial optimization problem:

\[\begin{split} \begin{aligned} \text{minimize} & \quad f(x) \\ \text{subject to} & \quad g_1(x) \leq c_1, \, g_2(x) \leq c_2, \, \ldots, \, g_m(x) \leq c_m \end{aligned} \end{split}\]

into a combinatorial optimization problem without constraints by computing \(m\) sets of penalty functions \(p_1, p_2, \ldots, p_m\) for each of the constraint conditions:

\[ \text{minimize} \quad f(x) + k_1 p_1(x) + k_2 p_2(x) + \cdots + k_m p_m(x). \]

Here, \(k_1, k_2, \ldots, k_m\) are sufficiently large hyperparameters.

For instances of the Constraint class, You can obtain its penalty function using the penalty property.

For example, the penalty function of a constraint object created using the equal_to() helper function is as follows.

from amplify import VariableGenerator, equal_to

gen = VariableGenerator()
q = gen.array("Binary", 6)
c = equal_to(q[0] + q[1] + q[2], 1)

>>> print(c.penalty)
2 q_0 q_1 + 2 q_0 q_2 + 2 q_1 q_2 - q_0 - q_1 - q_2 + 1

c.penalty takes the value 0 if c is satisfied, i.e., if q[0] + q[1] + q[2] == 1. It takes a value greater than 0 otherwise. It therefore satisfies the requirements of the penalty function that the previous section defines.

Penalty function weight

Note

With its default setting, AmplifyAEClient needs no weight from you. Amplify AE adjusts the weight of each constraint by itself, and the Amplify SDK does not send weight to it. This section applies to a solver that uses a penalty function for a constraint. It also applies to AmplifyAEClient when you set its solver property to Pubo or to Qubo.

The penalty function that the SDK adds to the objective function behaves like a penalty. Only an unsatisfied constraint makes it increase the value of the objective function. The solver searches for a solution that minimizes the sum of the objective function and the penalty. You can therefore expect both to be small in the solution that you obtain.

The penalty method works only if you set the weights of the constraints (\(k_1\), \(k_2\), \(\ldots\), \(k_m\) above) to values that match the objective function. The reason is that the solver minimizes the sum of the objective function and the penalty value. Suppose that the penalty for a violated constraint is small. A solution that breaks the constraint and has a smaller objective function then looks “better” than a solution that keeps it. The solver must therefore put the penalty value before the objective function.

The weight property of the Constraint object sets and gives the weight of the penalty function. The Amplify SDK sets its initial value to 1. The penalty function that the SDK adds to the objective function is the function that the SDK computes, multiplied by weight.

gen = VariableGenerator()
q = gen.array("Binary", shape=(2, 3))
c_list = equal_to(q, 1, axis=1)

>>> print(c_list)
[q_{0,0} + q_{0,1} + q_{0,2} == 1 (weight: 1),
 q_{1,0} + q_{1,1} + q_{1,2} == 1 (weight: 1)]
>>> c_list[0].weight
1.0

The next section shows that the Amplify SDK normalizes the penalty function by default, so that the penalty takes a value of at least 1. The value that weight needs depends on the problem. A rule of thumb is to set it larger than the part of the objective function that the constraint touches. A feasible solution is then more likely. If that part is far from 1, set the weight of the penalty function accordingly.

You can set the weight of the penalty function in two ways. The first way assigns the value to the weight property. The second way multiplies a constraint object Constraint or a constraint list ConstraintList by a number. The second way sets the same weight on more than one constraint object at once. It also sets a weight after you build the model.

>>> c_list *= 2.0
>>> print(c_list)
[q_{0,0} + q_{0,1} + q_{0,2} == 1 (weight: 2),
 q_{1,0} + q_{1,1} + q_{1,2} == 1 (weight: 2)]

Hint

For example, the following procedure can be used to adjust the weights of the penalty function.

  1. Estimate the maximum gain of the objective function that a violated constraint gives.

    • In the traveling salesperson problem, you can estimate the maximum gain from not visiting a city as the longest edge between two cities.

  2. Set the weights of the penalty function for the constraint of interest to a value greater than the estimated value.

  3. If no feasible solution is obtained, increase the weights and rerun the solver.

    • For example, you might increase the weights by a factor of 2.

If you want to increase the accuracy of the solution, follow the steps below after performing the above steps.

  1. Run the solver again and again, and decrease the weights one step at a time.

    • For example, you can decrease the weight by a factor of 0.9.

  2. Of all the solutions obtained, the one with the smallest objective function value is the final solution.

Penalty function auto-generation

A helper function such as equal_to() or less_equal() makes a constraint object Constraint. The Amplify SDK generates an optimal penalty function for that object as needed.

To generate a penalty function, the Amplify SDK first estimates the bounds of the left-hand side of the constraint expression. For a binary polynomial, it adds the negative coefficients and the positive coefficients of the polynomial separately.

The type of the constraint and the algorithm that you name then decide the penalty function, as follows.

Equality constraint

For equality constraint \(f(x) = c\), a penalty function \(p\) is generated as follows.

(i) If the lower bound of the range of possible values of the left-hand side \(f\) of the constraint expression is \(c\)

The penalty function \(p\) is set to be \(f - c\).

Take the product of binary variables q[0] * q[1] = 0. The lower bound of \(f\) and the right-hand side \(c\) are equal, so the SDK generates this penalty function.

>>> c = equal_to(q[0] * q[1], 0)
>>> print(c.penalty)
q_0 q_1
(ii) When the upper bound of the range of possible values of the left-hand side \(f\) of the constraint equation is \(c\)

The penalty function \(p\) is set to be \(c - f\).

Take the product of binary variables q[0] * q[1] = 1. The upper bound of \(f\) and the right-hand side \(c\) are equal, so the SDK generates this penalty function.

>>> c = equal_to(q[0] * q[1], 1)
>>> print(c.penalty)
- q_0 q_1 + 1
(iii) Other cases

The penalty function \(p\) is set to be \((f - c)^2\).

For example, for a binary variable sum q[0] + q[1] + q[2] = 2, the following penalty function is generated.

>>> c = equal_to(q[0] + q[1] + q[2], 2)
>>> print(c.penalty)
2 q_0 q_1 + 2 q_0 q_2 + 2 q_1 q_2 - 3 q_0 - 3 q_1 - 3 q_2 + 4

Note

The above applies to an inequality constraint where the upper and lower bounds are equal, and you can regard such constraint as an equality constraint.

Inequality constraints

To generate a penalty function for an inequality constraint, the Amplify SDK first rewrites the constraint into the form \(a \leq f \leq b\). It takes \(a\) from the lower bound and \(b\) from the upper bound of the left-hand side of the constraint expression. The Amplify SDK gives the following algorithms for a constraint of that form. You specify the algorithm with the penalty_formulation keyword argument of the less_equal(), greater_equal(), and clamp() inequality constraint-generating helper functions. The default is Default.

Default (Default)

If every variable and every coefficient in the constraint expression is an integer, the Amplify SDK uses the IntegerVariable algorithm. It uses the RealVariable algorithm otherwise.

IntegerVariable

Generates a penalty with auxiliary variables that take integer values. For an inequality constraint \(a \leq f \leq b\), the Amplify SDK issues an integer variable \(n\) that takes a value between \(a\) and \(b\). The penalty function of the equality constraint \(f - n = 0\) is then the penalty function of the inequality constraint.

However, as an exception, if \(b - a = 1\), then \((f - a)(f - b) / 2\) is the penalty function without issuing an integer variable.

Attention

An IntegerVariable on a constraint that is not an integer value does not give an exact formulation.

For example, for a binary variable sum q[0] + q[1] + q[2] <= 2, the following penalty function is generated

>>> c = less_equal(q[0] + q[1] + q[2], 2, penalty_formulation="IntegerVariable")
>>> print(c.penalty)
2 q_0 q_1 + 2 q_0 q_2 - 2 q_0 n_0 + 2 q_1 q_2 - 2 q_1 n_0 - 2 q_2 n_0
  + n_0^2 + q_0 + q_1 + q_2
>>> print(gen.variables[3])
{name: n_0, id: 3, type: Integer, lower_bound: -0, upper_bound: 2}

Here, n_0 is an integer auxiliary variable output when the Amplify SDK generates the penalty function.

Note

The Amplify SDK performs variable conversions during intermediate model construction for solvers that cannot handle integer variables, such as the QUBO solver. See “Variable Conversion and Degree Reduction” for details.

In this variable conversion, a wider range of an integer variable needs more auxiliary variables. To make the conversion more efficient, narrow the range of the inequality constraint while you build it. Dividing both sides by a common factor is one way.

Take the binary variable sum q[0] + q[1] + q[2] <= 1. The lower bound \(0\) of the left-hand side \(f\) and the right-hand side \(c=1\) differ by 1, so the SDK issues no integer variable. It generates the following penalty function.

>>> c = less_equal(q[0] + q[1] + q[2], 1, penalty_formulation="IntegerVariable")
>>> print(c.penalty)
q_0 q_1 + q_0 q_2 + q_1 q_2
RealVariable

Generates a penalty function with an auxiliary variable that takes a real value. For an inequality constraint \(a \leq f \leq b\), the Amplify SDK issues a real variable \(x\) that takes a value between \(a\) and \(b\). The penalty function of the equality constraint \(f - x = 0\) is then the penalty function of the inequality constraint.

For example, for the sum of the binary variables 0.1 * q[0] + 0.2 * q[1] + 0.4 * q[2] <= 0.5, the following penalty function is generated.

>>> c = less_equal(0.1 * q[0] + 0.2 * q[1] + 0.4 * q[2], 0.5, penalty_formulation="RealVariable")
>>> print(c.penalty)
0.04 q_0 q_1 + 0.08 q_0 q_2 - 0.2 q_0 x_0 + 0.16 q_1 q_2 - 0.4 q_1 x_0 - 0.8 q_2 x_0 + x_0^2 + 0.01 q_0 + 0.04 q_1 + 0.16 q_2
>>> print(gen.variables[3])
{name: x_0, id: 3, type: Real, lower_bound: 0, upper_bound: 0.5}

Here, x_0 is the real auxiliary variable that is issued when the penalty function is generated.

Note

The Amplify SDK performs variable conversions from real to binary variables for solvers that cannot handle real variables, such as the QUBO solver. The real_encoding_method keyword argument of the solve() function then specifies the conversion method. However, if the number of binary variables used for the conversion is small, the constraint conditions may not be accurately expressed.

In such cases, the following approach can be considered to avoid using real variables.

  • Multiply both sides of the constraint equation by a constant and use IntegerVariable.

  • Try approximation by Relaxation described next.

Relaxation

Uses LinearRelaxation when both of these hold. The left-hand side \(f\) of the constraint is quadratic or of a higher degree. The lower bound or the upper bound of \(f\) agrees with the range of the constraint. It uses QuadraticRelaxation in every other case.

Attention

The relaxation methods (Relaxation, LinearRelaxation, QuadraticRelaxation) do not give an exact formulation. Their penalty functions do not meet the requirements that the penalty method defines.

LinearRelaxation applies the penalty in the direction that moves \(f\) closer to \(a\) or to \(b\). QuadraticRelaxation applies it in the direction that moves \(f\) closer to \((a + b) / 2\). The penalty value is not constant while the constraints hold, so the solver may not find an optimal solution even with sufficiently large penalty weights. To find a better solution, run the solver many times and change the weights of the penalty function between the runs. Then take the best solution among the feasible ones.

LinearRelaxation

Uses Lagrangian relaxation instead of the penalty method.
For the inequality constraint \(a \leq f \leq b\), if the lower bound of \(f\) matches \(a\), the Amplify SDK generates a normalized \(f - a\). If the upper bound of \(f\) matches \(b\), it generates a normalized \(b - f\). It uses QuadraticRelaxation if neither bound matches.

For example, for a binary variable sum q[0] + q[1] + q[2] <= 2, the following penalty function is generated.

>>> c = less_equal(q[0] + q[1] + q[2], 2, penalty_formulation="LinearRelaxation")
>>> print(c.penalty)
0.5 q_0 + 0.5 q_1 + 0.5 q_2
QuadraticRelaxation

Uses Lagrangian relaxation instead of the penalty method. For inequality constraint \(a \leq f \leq b\), a normalized version of \((f - \left(a + b\right)/2)^2\) is generated.

For example, for a sum of binary variables q[0] + q[1] + q[2] <= 2, the following penalty function is generated.

>>> c = less_equal(q[0] + q[1] + q[2], 2, penalty_formulation="QuadraticRelaxation")
>>> print(c.penalty)
2 q_0 q_1 + 2 q_0 q_2 + 2 q_1 q_2 - q_0 - q_1 - q_2 + 1

Specifying a penalty function

For constraints created with helper functions, the Amplify SDK automatically specifies a penalty function. On the other hand, if you want to set a custom penalty function, you can call the Constraint constructor.

As shown below, you can create a constraint object by specifying a constraint expression and a penalty function.

  • Creating an equality constraint with constraint expression \(q_0 + q_1 = 2\) and penalty \(- q_0 - q_1\):

>>> c = Constraint(q[0] + q[1], eq=2, penalty=-q[0] - q[1])

  • Creating an inequality constraint with a constraint expression of \(q_0 + q_1 \leq 1\) and a penalty of \(q_0 q_1\):

>>> c = Constraint(q[0] + q[1], le=1, penalty=q[0] * q[1])

  • Creating an inequality constraint with constraint expression \(q_0 + q_1 \geq 1\) and penalty \(q_0 q_1 - q_0 - q_1\):

>>> c = Constraint(q[0] + q[1], ge=1, penalty=q[0] * q[1] - q[0] - q[1])

  • Creating an inequality constraint with constraint expression \(1 \leq q_0 + q_1 + q_2 \leq 2\) and penalty \((q_0 + q_1 + q_2 - 1)(q_0 + q_1 + q_2 - 2)\):

>>> f = q[0] + q[1] + q[2]
>>> c = Constraint(f, bounds=(1, 2), penalty=(f - 1) * (f - 2))