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 builds an intermediate model compatible with the combinatorial optimization solver by converting the variables, the polynomial degree of the constraints, and the objective function. If these conversions are impossible, the solver is invoked by generating penalty functions corresponding to the original constraints and automatically adding the penalty functions to the objective function. These procedures allow models with constraints to be solved even by combinatorial optimization solvers that cannot accept constraints as they are.

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,
    FixstarsClient,
    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, for a solver client that can directly handle equality constraints consisting of first-order binary variables, the following transformation from the input model to the intermediate model is performed.

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

Similar to the variable conversion in the objective function, we can see that the constraints included in the intermediate model consist only of binary variables, and the Amplify SDK performs the variable conversion between integer variables and binary variables in 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 constraint expressions by variable conversion or degree reduction, depending on the type and degree of constraints the solver can handle. However, if conversion is impossible, a polynomial expression called a penalty function is assigned to the objective function instead of the constraint condition expression. Because the Amplify SDK generates the penalty function equivalent to the corresponding constraint, the constraint is realized indirectly. When constructing the intermediate model, variable conversion and degree reduction of the penalty function are performed in addition to the objective function, if 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{align*} \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{align*} \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

The value of c.penalty takes the value 0 if c is satisfied, i.e., q[0] + q[1] + q[2] == 1, and a value greater than 0 otherwise. In other words, it satisfies the requirements of the penalty function defined in the previous section.

Penalty function weight

The penalty function added to the objective function behaves like a penalty: it increases the value of the objective function only if the constraints are not satisfied. Since the solver searches for a solution that minimizes the sum of the objective function and penalty values, you can expect both the objective function and the penalty to be small in the obtained solution.

An essential aspect of the penalty method that makes it work well is that you must set the weights of the constraints (\(k_1\), \(k_2\), \(\ldots\), \(k_m\) above) appropriately. This requirement is because the solver minimizes the sum of the objective function and the penalty value. Suppose the value of the penalty for violating a constraint is small. In that case, a solution that does not satisfy the constraint but has a minimal objective function may be “better” than a solution that satisfies the constraint but has a minimal objective function. In other words, the solver must prioritize minimizing the penalty value over minimizing the objective function.

You can set and obtain the weight of the penalty function using the weight property of the Constraint object; the Amplify SDK sets the initial value of the weight to be 1, and the penalty function added to the objective function is the value automatically calculated by the Amplify SDK multiplied by the 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

As described in the next section, by default, the Amplify SDK automatically normalizes and sets the penalty function to take a value of at least 1 as the penalty. The appropriate value for the weight depends on the problem, but a rule of thumb is to set it larger than the value of the objective function involved in the constraint so that a feasible solution is more likely to be obtained. Therefore, if its value is far from 1, you should set the weight of the penalty function accordingly.

You can set the weight of the penalty function directly by assigning it to the weight property, or you can set it by applying a numeric value to the constraint object Constraint or the constraint list ConstraintList. The latter is useful when you want to set the same weight value for multiple constraint objects at once or when you want to set the weights after the model has been constructed.

>>> 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)]

For more information on setting the weights of the penalty function for specific problems, please refer to the following:

• Traveling Salesperson Problem

• Quadratic Allocation Problem

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 with respect to a given constraint when that constraint is not satisfied.

   • For example, in the traveling salesperson problem, the maximum gain from not visiting a city can be estimated as the maximum length of the edges between 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.

4. Run the solver several times, decreasing the weights one at a time.

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

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

Penalty function auto-generation

When you create a constraint object Constraint using helper functions such as equal_to() or less_equal(), the Amplify SDK automatically generates an optimal penalty function as needed.

When the Amplify SDK generates a penalty function for a constraint, it first estimates the upper and lower bounds of possible values for the left-hand side of the constraint expression. For example, for a binary polynomial, the Amplify SDK obtains this estimate by calculating the sum of the negative and positive coefficients in the polynomial.

Then, depending on the type of equality or inequality constraint and the specified algorithm, a penalty function is generated 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\).

For example, for 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 following penalty function is generated.

>>> 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\).

For example, for the product of binary variables q[0] * q[1] = 1, the upper bound of \(f\) and the right-hand side of \(c\) are equal, so the following penalty function is generated.

>>> 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 using the upper and lower bounds on the possible values for the left-hand side of the constraint expression in the form \(a \leq f \leq b\). The Amplify SDK provides the following algorithms to generate penalty functions for inequality constraints of this form. The algorithm is specified using 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 all variables and coefficients in the constraint expression are integers, the Amplify SDK uses the IntegerVariable algorithm; otherwise, it uses the RealVariable algorithm.

IntegerVariable

Generates a penalty using 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\), and the penalty function for the equality constraint \(f - n = 0\) is the penalty function for the inequality constraint \(a \leq f \leq b\).

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

Attention

Note that specifying an IntegerVariable for a constraint that is not an integer value does not result in 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, the larger the range of possible values of an integer variable, the larger the number of auxiliary variables required. Therefore, you can make the variable conversion more efficient by reducing the range of inequality constraints at the time of constraint construction, for example, by dividing both sides by a factor.

For a binary variable sum q[0] + q[1] + q[2] <= 1, the difference between the lower bound \(0\) of the left-hand side \(f\) and the right-hand side \(c=1\) of the constraint is 1, so no integer variables are issued and the following penalty function is generated.

>>> 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 using an auxiliary variable that takes a real value. For inequality constraint \(a \leq f \leq b\), the Amplify SDK issues a real variable \(x\) that takes values between \(a\) and \(b\), and the penalty function for equality constraint \(f - x = 0\) is the penalty function for inequality constraint \(a \leq f \leq b\).

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 if the left-hand side of the constraint \(f\) is quadratic or greater and the lower or upper bound of \(f\) is consistent with the range of the constraint, otherwise uses QuadraticRelaxation.

Attention

Note that the penalty functions generated by the relaxation methods (Relaxation, LinearRelaxation, QuadraticRelaxation) are not exact formulations because they do not meet the requirements for a penalty function as defined in the penalty method.

The penalty is applied in the direction of moving \(f\) closer to \(a\) or \(b\) when LinearRelaxation is applied, and in the direction of moving \(f\) closer to \((a + b) / 2\) when QuadraticRelaxation is applied. However, since the penalty value applied when the constraints are satisfied is not constant, an optimal solution may not be found even if the penalty weights are sufficiently large. Therefore, to find a better solution, the solver should be run many times while changing the weights of the penalty function, and the best solution among the feasible solutions should be selected.

LinearRelaxation

Uses Lagrangian relaxation instead of the penalty method.
For the inequality constraint \(a \leq f \leq b\), if the lower bound of the range of possible values of \(f\) matches \(a\), a normalized version of \(f - a\) is generated; if the upper bound of the range of possible values of \(f\) matches \(b\), a normalized version of \(b - f\) is generated. If neither is the case, QuadraticRelaxation is used.

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))