Solving constrained optimization problems¶

In this example, we consider the following constrained optimization problem:

\[ \begin{array}{r c} \text{maximize } & a + b + c \\ \text{subject to } & 2a + 3b \leq 45 \\ & 5a + 2c \leq 75 \\ & 3b + c \leq 50 \\ & -100 \leq a \leq 100 \\ & -100 \leq b \leq 100 \\ & -100 \leq c \leq 100 \\ \end{array} \]

We will now solve this optimization problem using:

Evolutionary computation (using PGPE with constraint penalties)
Penalty method (using gradient-based search with ClipUp, with penalty multipliers iteratively incremented)
Interior points method (with the help of the functional API of PyTorch and the log-barrier function of EvoTorch)

This notebook demonstrates:

How to use functional algorithms of EvoTorch (evolutionary and gradient-based)
How to use constraint penalization utilities of EvoTorch and use them as building blocks for constrained optimization, both for evolutionary and for gradient-based optimization

from evotorch.algorithms.functional import pgpe, pgpe_ask, pgpe_tell, clipup, clipup_ask, clipup_tell
from evotorch.decorators import expects_ndim
from evotorch.tools import penalty, log_barrier

import torch
import torch.func as tfunc
from typing import Union
from functools import partial

from matplotlib import pyplot as plt
import numpy as np

from datetime import datetime

Fitness function implementation¶

We begin by implementing our fitness function. This fitness function has two arguments:

penalty_multiplier: Expected as a scalar, this value represents the multiplier for the negative penalty quantities that will be added onto the fitness. Higher values for penalty_multiplier will result in stronger penalizations.
x: A 1-dimensional tensor, that represents the solution that will be evaluated.

Notice how the fitness function below is decorated by @expects_ndim(0, 1). This decorators declares that the first positional argument (penalty_multipliers) expects a 0-dimensional tensor, and the second positional argument (x) expects a 1-dimensional tensor. If any of these arguments are given with more dimensions, those extra dimensions will be considered as batch dimensions by the decorated function. This auto-batching feature is very helpful, because it allows the decorated function f to work with not just a single solution, but with a population of solutions (or a batch of populations of solutions) when x is given as an n-dimensional tensor with n>1. We will use this auto-batching feature with PGPE.

@expects_ndim(0, 1)
def f(penalty_multiplier: torch.Tensor, x: torch.Tensor) -> torch.Tensor:
    a = x[0]
    b = x[1]
    c = x[2]
    objective = a + b + c

    constraints = [
        [(2 * a) + (3 * b), "<=", 45],
        [(5 * a) + (2 * c), "<=", 75],
        [(3 * b) + c, "<=", 50],
        [a, ">=", -100],
        [a, "<=", 100],
        [b, ">=", -100],
        [b, "<=", 100],
        [c, ">=", -100],
        [c, "<=", 100],
    ]

    penalty_amount = 0.0
    for lhs, op, rhs in constraints:
        # For each constraint, we add a penalty (if there is violation)
        penalty_amount = penalty_amount + penalty(
            # Left-hand-side, comparison operator, and the right-hand-side:
            lhs,
            op,
            rhs,
            #
            # Because this is a function we wish to maximize, the penalties should be in the opposite direction.
            # Therefore, we declare the sign of the penalties as "-", making them negative quantities:
            penalty_sign="-",
            #
            # There will be a penalty in the form: (linear * amount_of_violation)
            linear=1.0,
            #
            # There will also be a penalty in the form: (amount_of_violation ** exp)
            exp=2.0,
            #
            # The exponential penalties are not allowed to exceed this amount:
            exp_inf=5000.0,
            #
            # There will also be a penalty in the form: (step if amount_of_violation > 0 else 0)
            step=10000.0,
        )

    # Scale the accumulated penalty by `penalty_multiplier`, add it onto the objective, then return it
    return objective + (penalty_amount * penalty_multiplier)

Evolutionary computation¶

As the evolutionary algorithm, we use the functional pgpe.

pgpe_state = pgpe(
    # We want to maximize the evaluation results:
    objective_sense="max",

    # The center point of the initial search distribution is given as 0 vector:
    center_init=torch.zeros(3),

    # Learning rates for the center and the standard deviation of the
    # search distribution:
    center_learning_rate=0.1,
    stdev_learning_rate=0.1,

    # The ranking method is "centered", which ranks the worst solution as
    # -0.5, and the best solution as +0.5:
    ranking_method="centered",

    # We use the ClipUp optimizer:
    optimizer="clipup",

    # In the case of this example problem, we use the following max_speed:
    optimizer_config={"max_speed": 1.0},

    # The standard deviation for the initial search distribution:
    stdev_init=100.0,

    # Relative maximum allowed change for the standard deviation.
    # 0.2 means that a standard deviation value cannot change more than 20%
    # of its original value.
    stdev_max_change=0.2,

    # Minimum standard deviation value. The standard deviation is cannot shrink
    # to values lower than this:
    stdev_min=0.01,
)

pgpe_state

Main loop of the evolutionary computation:

num_generations = 1000  # For how many iterations will PGPE work

time_interval = 1  # With this period, we will print the status
last_report_time = datetime.now()

# Initialize the variables that will store the best solution and the best evaluation result
pgpe_best_eval = float("-inf")
pgpe_best = None

for generation in range(1, 1 + num_generations):
    # Ask for a population from PGPE
    population = pgpe_ask(pgpe_state, popsize=1000)

    # Compute the fitnesses
    fitnesses = f(1, population)

    # Inform PGPE of the latest fitnesses, and get its next state
    pgpe_state = pgpe_tell(pgpe_state, population, fitnesses)

    # From the most recent population and fitnesses, update the best solution and the best evaluation result
    pop_best_index = torch.argmax(fitnesses)
    pop_best = population[pop_best_index]
    pop_best_eval = fitnesses[pop_best_index]
    if pop_best_eval > pgpe_best_eval:
        pgpe_best_eval = pop_best_eval
        pgpe_best = population[pop_best_index, :]

    # If it is time to report, print the status
    tnow = datetime.now()
    if ((tnow - last_report_time).total_seconds() > time_interval) or (generation == num_generations):
        print("best solution:", pgpe_best, "best evaluation:", pgpe_best_eval)
        last_report_time = tnow

Best solution:

pgpe_best, pgpe_best_eval

Penalty method¶

As an alternative to evolutionary computation, we now use the gradient-based penalty method below. The main points of the penalty method are as follows:

Use a gradient-based search on a fitness function augmented by penalties (i.e. by penalty terms that are computed according to how much the constraints are penalized)
Periodically increase the multipliers for the penalty terms (when the penalty multiplier is increased, we use the previous optimization's solution as the new starting point)

# Try with these penalty multipliers, ordered from small to large:
penalty_multipliers = [0.1, 1, 4]

# For each penalty multiplier, we do the search for this number of iterations:
num_iterations = 1000

# Initialize the variables that will store the best solution and the best evaluation result
clipup_best_eval = float("-inf")
clipup_best = None

# Initialize the search point as a 0 vector:
x = torch.zeros(3)

for penalty_multiplier in penalty_multipliers:
    # Initialize the ClipUp algorithm for the currently considered penalty multiplier
    clipup_state = clipup(center_init=x, center_learning_rate=0.1, max_speed=1.0)

    # Optimization loop for the current penalty multiplier
    for iteration in range(1, 1 + num_iterations):
        # Ask ClipUp for the current search point
        x = clipup_ask(clipup_state)

        # Compute the gradient and the fitness
        gradient, fitness = tfunc.grad_and_value(partial(f, penalty_multiplier))(x)

        # Update the best-known solution so far
        if fitness > clipup_best_eval:
            clipup_best_eval = fitness
            clipup_best = x

        # Inform ClipUp of the latest gradient, and get its next state
        clipup_state = clipup_tell(clipup_state, follow_grad=gradient)

    # After each optimization loop, print the best known solution
    print("best solution:", clipup_best, " best evaluation:", clipup_best_eval)

Interior points method¶

Although, as its name implies, the main focus of EvoTorch is evolutionary computation, we now demonstrate that it is possible to implement an interior points method by combining:

torch.func.grad
torch.func.hessian
evotorch.tools.log_barrier

The grad(...) and hessian(...) functions are basic building blocks for implementing a Newton-Raphson search. We penalize proximities to the infeasible regions with the help of log_barrier(...).

We begin with a modified implementation of our fitness function, where evotorch.tools.penalty(...) are replaced by evotorch.tools.log_barrier(...).

@expects_ndim(0, 1)
def f_with_log_barriers(sharpness: torch.Tensor, x: torch.Tensor) -> torch.Tensor:
    a = x[0]
    b = x[1]
    c = x[2]
    objective = a + b + c

    constraints = [
        [(2 * a) + (3 * b), "<=", 45],
        [(5 * a) + (2 * c), "<=", 75],
        [(3 * b) + c, "<=", 50],
        [a, ">=", -100],
        [a, "<=", 100],
        [b, ">=", -100],
        [b, "<=", 100],
        [c, ">=", -100],
        [c, "<=", 100],
    ]

    penalty_amount = 0.0
    for lhs, op, rhs in constraints:
        penalty_amount = penalty_amount + log_barrier(lhs, op, rhs, penalty_sign="-", sharpness=sharpness)

    # Return the objective with the log-barrier penalties added onto it.
    # Notice that in the end, we are inverting the sign of the returned quantity.
    # This will allow us to implement the Newton-Raphson's search method from a minimization perspective.
    return -(objective + penalty_amount)

# Start the search from the 0 vector:
x = torch.zeros(3)

# Try with these sharpness values for the log barrier, ordered from small to large:
sharpness_values = [1, 10, 100, 1000, 10000]

# Learning rate for when taking a step
learning_rate = 0.1

# An identity matrix multiplied by the constant below will be added to the Hessian matrix.
# When a diagonal element of the Hessian matrix is 0 because of numerical issues, this trick will allow the
# algorithm to still take a step.
I_multiplier = 0.01

# With this interval (in seconds), we wish to report the status
reporting_interval = 1
last_report_time = datetime.now()

# The interior points method will run for this amount of iterations:
num_iterations = 200

for sharpness in sharpness_values:
    print("sharpness:", sharpness)
    print()

    for iteration in range(1, 1 + num_iterations):
        # Compute the gradient and the solution cost
        g, cost = tfunc.grad_and_value(partial(f_with_log_barriers, sharpness))(x)

        # Compute the Hessian matrix
        H = tfunc.hessian(partial(f_with_log_barriers, sharpness))(x)

        # Add the identity matrix multiplied by a constant to the Hessian
        H = H + (I_multiplier * torch.eye(H.shape[-1]))

        # Take the inverse of the Hessian matrix
        invH = torch.linalg.inv(H)

        # Move the center point
        x = x - learning_rate * (invH @ g)

        # If it is time to report, print the status
        tnow = datetime.now()
        if (tnow - last_report_time).total_seconds() > reporting_interval:
            print("Iteration:", iteration, "  Gradient norm:", torch.linalg.norm(g), "  Solution cost:", cost)
            last_report_time = tnow

    # Print the current search point
    print()
    print("x:", x)
    print()

See this notebook on GitHub