Performing an Optimisation

Consider the following example to demonstrate how to set up and perform an optimisation with some non-linear constraints.

Suppose that we wish to find

\[\min_{x \in \mathbb{R}^2} \sqrt{x_2} \tag{1}\]

subject to constraints

\[x_2 \ge (a_1x_1 + b_1)^3 \tag{2}\]

\[x_2 \ge (a_2 x_1 + b_2)^3 \tag{3}\]

for parameters \(a_1 = 2\), \(b_1 = 0\), \(a_2 = -1\), \(b_2 = 1\).

This problem expects a minimum at \(\boldsymbol{x} = ( \frac{1}{3}, \frac{8}{27} )\).

Note

This example is ripped straight from the NLOpt docs.

Objective Function

First we must define our objective function; the function we wish to minimise. The function must take a single argument: a 1-D numpy array, where each element in the array is an optimisation parameter. It must return a float.

For our example, our objective function looks like this, note that we don’t actually need to use \(x_1\) (x[0]) here:

def f_objective(x: np.ndarray) -> float:
    """Objective function for optimisation to find solution to eqn (1)."""
    return np.sqrt(x[1])

If using a gradient-based optimisation algorithm, it helps to define an analytical gradient. If you do not give an analytical gradient, a numerical one will be estimated.

As with the objective, the gradient function must take a 1-D numpy array, containing the optimisation parameters as its only argument. It must return a numpy array of the same length, where each index \(i\) contains the partial derivative \(\frac{\partial f}{\partial x_i}\) for the corresponding optimisation parameter.

def df_objective(x) -> np.ndarray:
    """Gradient of the objective function."""
    return np.array([0.0, 0.5 / np.sqrt(x[1])])

Constraints

A constraint function must take a single argument, the array of optimisation parameters. Both vector-valued and scalar constraints are supported, so the function can return a float, or, an array of length \(m\), where \(m\) is the dimensionality of the constraint.

Given a constraint function \(f_c\), and optimisation parameters \(\boldsymbol{x}\), an equality constraint has the form

\[f_c(\boldsymbol{x}) = 0\]

and an inequality constraint has the form

\[f_c(\boldsymbol{x}) \le 0\]

In our example, we have two inequality constraints. We’ll need to adjust these slightly to match the required form. Equation 2 becomes

\[(a_1 x_1 + b_1)^3 - x_2 \le 0\]

and equation 3 becomes

\[(a_2 x_1 + b_2)^3 - x_2 \le 0\]

Notice that these two constraints are similar, so we can define a single function, then use lambdas to set a and b to get the required values.

def f_constraint(x: np.ndarray, a: float, b: float) -> np.ndarray:
    """Inequality constraint."""
    return np.array([(a * x[0] + b) ** 3 - x[1]])

f_constraint_1 = lambda x: f_constraint(x, 2, 0)
f_constraint_2 = lambda x: f_constraint(x, -1, 1)

We can also define the gradient of the constraint. Note that this isn’t strictly necessary, as, if a gradient-based optimiser is used, a numerical approximation is made. However, an analytical gradient will be more reliable.

The constraint’s gradient function takes the array of optimisation parameters, and returns an array with shape \(m \times n\).

The partial derivatives of the constraint in our example are

\begin{gather*} \frac{\partial f_c}{\partial x_1} = 3a(a x_1 + b)^2 \\ \frac{\partial f_c}{\partial x_2} = -1 \end{gather*}

So our Python function will be

def df_constraint(x: np.ndarray, a: float, b: float) -> np.ndarray:
    """Inequality constraint gradient."""
    return np.array([3 * a * (a * x[0] + b) ** 2, -1.0])

Note that we are using two separate constraints here, but it can sometimes be more convenient to express multiple constraints in a single vector-valued one. In this case that vector-valued constraint, and its gradient, could look like this

def vector_constraint(x: np.ndarray, a1: float, b1: float, a2: float, b2: float) -> np.ndarray:
    return np.array([f_constraint(x, a1, b1), f_constraint(x, a2, b2)])


def d_vector_constraint(x: np.ndarray, a1: float, b1: float, a2: float, b2: float) -> np.ndarray:
    return np.vstack([df_constraint(x, a1, b1), df_constraint(x, a2, b2)])

Note

Not all optimisation algorithms support non-linear constraints. They can only be used with SLSQP, COBYLA, and ISRES.

The Optimise Function

Now that we have our objective function and constraints, we can perform the optimisation. To do this, we use the optimise() function.

Within this function, we can select the optimisation algorithm to use, parameter bounds, stopping conditions, the initial guess (if not given, the center of the bounds is used), and whether to record the history of the optimisation parameters (this is off by default, as it impacts run time performance).

>>> from bluemira.optimisation import optimise
>>> result = optimise(
...     f_objective,
...     df_objective=df_objective,
...     algorithm="SLSQP",
...     x0=np.array([1, 1]),
...     opt_conditions={"xtol_rel": 1e-10, "max_eval": 1000},
...     keep_history=True,
...     bounds=(np.array([-np.inf, 0]), np.array([np.inf, np.inf])),
...     ineq_constraints=[
...         {
...             "f_constraint": lambda x: f_constraint(x, 2, 0),
...             "df_constraint": lambda x: df_constraint(x, 2, 0),
...             "tolerance": np.array([1e-8]),
...         },
...         {
...             "f_constraint": lambda x: f_constraint(x, -1, 1),
...             "df_constraint": lambda x: df_constraint(x, -1, 1),
...             "tolerance": np.array([1e-8]),
...         },
...     ],
... )
>>> print(result)
OptimiserResult(f_x=np.float64(0.54...), x=array([0.333..., 0.29...]), n_evals=..., constraints_satisfied=True)

The Optimisation Problem Class

As an alternative to the optimise() function, it is possible to take a class-based approach to performing an optimisation. This can have several benefits, including Liskov Substitution of optimisation problems, shared state between objective functions and constraints, and logical grouping of related functionality.

To define an optimisation problem, inherit from OptimisationProblem and implement the interface.

You must implement the objective() method.

You can optionally override:

df_objective()
Must return the gradient of the objective function at the given parameterisation. If this is not overridden, and a gradient-based algorithm is used, a gradient will be numerically estimated. See also, Objective Function.
eq_constraints()
Must return a list of ConstraintT dictionaries, defining equality constraints. See also, Constraints.
ineq_constraints()
Must return a list of ConstraintT dictionaries, defining inequality constraints. See also, Constraints.
bounds()
Must return the lower and upper bounds of the optimisation parameters. The default is to return (-np.inf, np.inf).

See here for an implemented example of an OptimisationProblem.

Available Optimisation Algorithms

There are several optimisation algorithms that can be used within bluemira. Including gradient and non-gradient based.

SLSQP
COBYLA
SBPLX
MMA
BFGS
DIRECT
DIRECT_L
CRS
ISRES

See the Algorithm enum for a reliably up-to-date list.