Optimization - Executable Engineering: Interactive Numerical Methods in Python

import numpy as np
import plotly.graph_objects as go
import scipy.optimize as optimize

Optimization asks us to:

Find $\mathbf{x}$ that minimizes $f(\mathbf{x})$ subject to $\mathbf{g}(\mathbf{x}) = \mathbf{0}$ and $\mathbf{h}(\mathbf{x}) \ge \mathbf{0}$

$f$ is the objective function and typically a scalar
$\mathbf{g}(\mathbf{x})$ is a list of equality constraints
$\mathbf{h}(\mathbf{x})$ is a list of inequality constraints
Maximization for $f$ is merely minimization of $-f$ .
Optimization differs from root finding in the constraints and smoothness of $f$ .

We saw with root finding that using gradient information to help search algorithms has important benefits. That remains true in optimization, but since they are often high dimensional ( $\mathbf{x}$ is a large vector), and $f$ is not necessarily well behaved, great value is placed on methods that do not require (direct) gradient / Hessian calculations.

Local vs Global minima¶

Even in 1D, optimization is confounded by local minima. The lowest of the local minima is the global minimum.

There is no good way to ensure a local minimum is a global one without calculating and comparing them all!

There are lots of algorithms that attempt to find the global minimum. They either:

launch multiple minimizers over the range of arguments
typically involve perturbing good solutions in an attempt to escape from local minima basins.

Interestingly, many of these methods are modelled after nature! (Why?)

Genetic algorithms - mutations, cross-overs, selection
Simulated annealing - metallurgical process that accepts worse solutions
Particle swarm - mimics swarming behaviour to track multiple solutions
Amoeba search (which we will talk about)

Introduction of constraints can further complicate things since sometimes the minimum is at a bound!

Example: local and global minima¶

Minimize $f(x) = x^2 - \cos(5 x)$

x = np.linspace(-2, 2, 200)
f = lambda x: x**2 - np.cos(5 * x)
y = f(x)

opt = optimize.minimize(f, 1)
local_min_x, local_min_y = opt.x[0], opt.fun

fig = go.Figure()
fig.add_trace(go.Scatter(x=x, y=y, mode='lines', name='f(x)'))
fig.add_trace(go.Scatter(x=[0], y=[f(0)], mode='markers', marker=dict(color='red', size=10), name='Global Minimum'))
fig.add_trace(go.Scatter(x=[local_min_x], y=[local_min_y], mode='markers', marker=dict(color='blue', size=10), name='Local Minimum'))
fig.update_layout(title='Function with Local and Global Minima', xaxis_title='x', yaxis_title='f(x)')
fig.show()

Example: Eggholder surfaces¶

Minimize $y \sin(\sqrt{|x + y|}) + x \sin(\sqrt{|x - y|})$

def eggholder(x, y):
    return y * np.sin(np.sqrt(np.abs(x + y))) + x * np.sin(np.sqrt(np.abs(x - y)))

x = np.linspace(-512, 512, 100)
y = np.linspace(-512, 512, 100)
X, Y = np.meshgrid(x, y)
Z = eggholder(X, Y) / 10

fig = go.Figure(data=[go.Surface(z=Z, x=X, y=Y)])
fig.update_layout(
    title='Eggholder Function',
    autosize=False,
    width=600,
    height=500,
    margin=dict(l=65, r=50, b=65, t=90)
)
fig.show()

Example: Constrained minimization¶

Minimize $-x^2 + x$ such that $|x| \le 1$

def f(x):
    return -x**2 + x

bounds = (-1, 1)
x = np.linspace(-1.5, 1.5, 200)
y = f(x)

result = optimize.minimize_scalar(f, bounds=bounds, method='bounded')
x_min, y_min = result.x, result.fun

fig = go.Figure()
fig.add_vrect(x0=-1, x1=1, fillcolor='lightgray', opacity=0.4, line_width=0, annotation_text='Feasible Region')
fig.add_trace(go.Scatter(x=x, y=y, mode='lines', name='f(x)'))
fig.add_trace(go.Scatter(x=[x_min], y=[y_min], mode='markers', marker=dict(color='red', size=10), name=f'Minimum ({x_min:.2f}, {y_min:.2f})'))
fig.update_layout(title='Plot of f(x) = -x^2 + x with |x| <= 1', xaxis_title='x', yaxis_title='f(x)')
fig.show()