Skip to main content

James Brind

Tabu search – theory and implementation

In this post, I discuss the multi-objective Tabu search for gradient-free solution of typical turbomachinery optimisation problems. There are already many good references outlining this algorithm, but I could not find a reusable Python implementation readily available. After a motivation of why Tabu is a good choice of algorithm, I describe the how the search works and present highlights from my implementation. The article concludes with results from a test problem. The full code is on my sourcehut.

Motivation

Any engineering design process can be formulated as an optimisation problem. We seek a point in the design space that maximises, say, performance or reliability while minimising cost and complexity; while we use the laws of physics to predict each of these factors, external human judgement determines their relative importance in a given situation. There are a large number of possible degrees of freedom available to the designer of a machine, but in practice they will select a much smaller parametrisation that remains flexible enough to navigate the design space and find the optimum.

Moving to the domain of computational fluid dynamics based turbomachinery design, we have more specific considerations. Each evaluation of the objectives, by solving the fluid flow equations, is rather expensive. Aerodynamic performance and mechanical constraints are in direct conflict, leading to a constrained design space.

The key to a successful outcome in computational engineering is not to reinvent the wheel, but to know and deploy the most suitable state-of-the-art algorithms for the problem at hand, developed by other people. Here, a multi-objective Tabu search fits our needs because it,

  • yields the entire Pareto front, enabling flexible judgement of trade-offs;
  • does not require gradients, which are expensive and noisy;
  • naturally handles constraints;
  • scales well to multi-dimensional design spaces.

As an aside, I’m not convinced about adjoint methods or the present fashion for neural network surrogate models. One needs extremely well-converged solutions to get consistent results from the adjoint equations, and I worry about local minima (although the prevalence of local minima in a typical turbine blade design is an open question). Rather than spend PhD-student hours implementing and debugging the additional complexity of a surrogate model, I would prefer to spend core hours brute-forcing a more exhaustive search.

Algorithm

We parametrise the design space using an \(n\)-dimensional vector \(\newcommand{\bv}[1]{\mathbf{#1}}\bv{x}\), and seek to optimise for the minimum of an \(m\)-dimensional objective \(\bv{y} = f(\bv{x})\). The design vector \(\bv{x}\) may be subject to any relevant constraints. In the scalar case of \(m=0\), ranking of candidate optimal \(y\) values is trivial. In the multiobjective case, we are interested in Pareto fronts. When comparing two candidate objectives, \(\bv{y}_0\) and \(\bv{y}_1\), there are three possible cases:

  • If every element of \(\bv{y}_0\) is less than \(\bv{y}_1\), then \(\bv{y}_0\) dominates \(\bv{y}_1\);
  • If every element of \(\bv{y}_0\) is greater than \(\bv{y}_1\), then \(\bv{y}_0\) is dominated by \(\bv{y}_1\);
  • Otherwise, \(\bv{y}_0\) and \(\bv{y}_1\) are Pareto-equivalent and form two points on a Pareto front.

Before starting the search, we initialise the working variables:

  • Local search counter \(i\) holding the number of iterations since the last improvement in the value of the objective;
  • Current step sizes \(\Delta\bv{x}\), a vector with \(n\) elements, one for each design parameter;
  • Short-term memory \(\newcommand{\Mem}[1]{\mathcal{M}_\mathrm{#1}}\Mem{sht}\), to store a fixed number of the most recently-visited \((\bv{x}, \bv{y})\) points which are Tabu;
  • Medium-term memory, \(\Mem{med}\), to store a set of current near-optimal \((\bv{x}, \bv{y})\) points. For \(m=1\), a scalar objective, \(\Mem{med}\) is is a ranked list containing a fixed number of the best points seen so far. For \(m>1\), multidimensional objectives, \(\Mem{med}\) is a Pareto front;
  • Long-term memory \(\Mem{lng}\), to record all \((\bv{x}, \bv{y})\) points visited by the optimiser.

Starting from an initial point \((\bv{x}_0,\bv{y}_0)\), at each iteration:

  1. Generate a set of \(2n\) candidate moves by perturbing the \(n\)th element of \(\bv{x}_0\) by \(\pm\Delta x_n\);
  2. Remove Tabu points, either which have been recently visited and exist in \(\Mem{sht}\), or which violate a constraint;
  3. Evaluate the objective function over non-Tabu candidate moves, reusing previous results from \(\Mem{lng}\) for old points that have been seen before;
  4. Store all new evaluations in \(\Mem{lng}\). If we have found any new near-optimal points, add them to \(\Mem{med}\) and reset \(i=0\), otherwise increment \(i\);
  5. What happens now depends on the value of \(i\):
    • Until \(i\) reaches a certain threshold, the algorithm performs a local search, choosing the next point \(\bv{x_1}\) as the best of the current candidate moves;
    • If \(i = i_\mathrm{diversify} \approx 10\), that is, the optimiser has not found any improvements for 10 moves, we conclude we are stuck in a local minima. To diversify the search, we generate a random \(\bv{x}_1\) in a sparse region of \(\Mem{lng}\);
    • If \(i = i_\mathrm{intensify} \approx 20\), and the optimiser has still found no improvements, select a near-optimal point from \(\Mem{med}\) as \(\bv{x}_1\). This intensification step biases the search towards known-good neighbourhoods, while giving the optimiser flexibility to explore the design space;
    • If \(i = i_\mathrm{restart} \approx 40\), we conclude that the design space has been sufficiently mapped, and no improvement is possible at our current resolution. So reduce the step \(\Delta\bv{x}\) by a factor of, say, 2. Restart the search taking a near-optimal next point \(\bv{x}_1\) from \(\Mem{med}\);
  6. Add the next point, \(\bv{x}_1\) to the Tabu list \(\Mem{sht}\);
  7. Halt if the step size \(\Delta\bv{x}\) reduces below the desired resolution, otherwise replace \(\bv{x}_0 \leftarrow \bv{x}_1\) and return to step 1.

Implementation

Memories

First of all, we need a memory data structure. This is a prime example of something that should be implemented as a class: a tight collection of related data and a few operations that we want to apply to that data. Classes seem most natural to me when they act as more specialised versions of built-in data types, such as the list in Python.

My instinct is to preallocate matrices (not very Pythonic, I know) for the design vectors and objective functions, \(\bv{X}\) and \(\bv{Y}\), where each row represents a new evaluation. Before the memory fills up, we do not want to expose empty rows outside of the class. So we keep track of the number of occupied rows using npts and apply the @property decorator to only return that many rows when accessing the X and Y attributes.

import numpy as np

class Memory:
    def __init__(self, nx, ny, max_points, tol=None):
        """Store a set of design vectors and their objective functions."""

        # Record inputs
        self.nx = nx
        self.ny = ny
        self.max_points = max_points
        self.tol = np.array(tol)

        # Initialise points counter
        self.npts = 0

        # Preallocate matrices for design vectors and objectives
        # Private because the user should not have to deal with empty slots
        self._X = np.empty((max_points, nx))
        self._Y = np.empty((max_points, ny))

    # Public read-only properties for X and Y
    @property
    def X(self):
        """The current set of design vectors."""
        return self._X[: self.npts, :]

    @property
    def Y(self):
        """The current set of objective functions."""
        return self._Y[: self.npts, :]

Elemental operations like Memory.add(), Memory.get(), Memory.delete() are relatively straightforward to implement.

One of the more frequent operations is to check is a given design point is already in the memory, i.e. we have evaluated it before or it is on the Tabu list in short term memory. This amounts to looking for matches of a test vector against each row of a matrix, something like the Matlab ismember function with the rows flag set. There isn’t a Numpy builtin for this, so we roll our own,

def find_rows(A, B, atol=None):
    """Get matching rows in matrices A and B.

    Return:
        logical same shape as A, True where A is in B
        indices same shape as A, of the first found row in B for each A row."""

    # Arrange the A points along a new dimension
    A1 = np.expand_dims(A, 1)

    # nA by nB logical, True where all columns elements match
    if atol is None:
        # Test for equality
        b = (A1 == B).all(axis=-1)
    else:
        # Test for equality within tolerance
        b = np.isclose(A1, B, atol=atol).all(axis=-1)

    # A index is True where it matches any of the B points
    ind_A = b.any(axis=1)

    # Use argmax to find first True along each row
    loc_B = np.argmax(b, axis=1)

    # Where there are no matches, override to a sentinel value -1
    loc_B[~ind_A] = -1

    return ind_A, loc_B

There are a couple of useful tricks here. expand_dims, being the opposite of squeeze, is a useful function to know. Numpy does not really have an equivalent of the find function in Matlab, so we mildly abuse the fact that, where elements are equal, argmax returns the lowest index. We can now implement Memory.contains(), returning true if a design vector is in memory; and Memory.lookup() to return the objective function for a design vector that has already been evaluated.

    def contains(self, Xtest):
        """Boolean index for each row in Xtest, True if x already in memory."""
        if self.npts:
            return find_rows(Xtest, self.X, self.tol)[0]
        else:
            return np.zeros((Xtest.shape[0],), dtype=bool)


    def lookup(self, Xtest):
        """Return objective function for design vector already in mem."""

        # Check that the requested points really are available
        if np.any(~self.contains(Xtest)):
            raise ValueError("The requested points have not been previously evaluated")

        return self.Y[find_rows(Xtest, self.X)[1]]

During the search, we want to keep the medium-term memory as a set of near-optimal points: a fixed number of the best seen scalar objectives, or a Pareto front for vector objectives. Building off our elemental operations, methods for these are not too difficult. Memory.update_front() is a marvel of logical indexing!

    def update_best(self, X, Y):
        """Add or remove points to keep the best N in memory."""

        X, Y = np.atleast_2d(X), np.atleast_2d(Y)

        # Join memory and test points into one matrix
        Yall = np.concatenate((self.Y,Y),axis=0)
        Xall = np.concatenate((self.X,X),axis=0)

        # Sort by objective, truncate to maximum number of points
        isort = np.argsort(Yall[:,0],axis=0)[:self.max_points]
        Xall, Yall = Xall[isort], Yall[isort]

        # Reassign to the memory
        self.npts = len(isort)
        self._X[: self.npts] = Xall
        self._Y[: self.npts] = Yall

        return np.any(self.contains(X))


    def update_front(self, X, Y):
        """Add or remove points to maintain a Pareto front."""
        Yopt = self._Y[: self.npts]

        # Arrange the test points along a new dimension
        Y1 = np.expand_dims(Y, 1)

        # False where an old point is dominated by a new point
        b_old = ~(Y1 < Yopt).all(axis=-1).any(axis=0)

        # False where a new point is dominated by an old point
        b_new = ~(Y1 >= Yopt).all(axis=-1).any(axis=1)

        # False where a new point is dominated by a new point
        b_self = ~(Y1 > Y).all(axis=-1).any(axis=1)

        # We only want new points that are non-dominated
        b_new_self = np.logical_and(b_new, b_self)

        # Delete old points that are now dominated by new points
        self.delete(~b_old)

        # Add new points
        self.add(X[b_new_self], Y[b_new_self])

        # Return true if we added at least one point
        return (np.sum(b_new_self) > 0)

Finally, we need a few methods to sample points from the memory. We can either pick a random point, Memory.sample_random(), or bin each component of the design vector into a histogram and pick the sparsest bin, Memory.sample_sparse(). To really diversify the search, we can generate an entirely new point in a sparse region of the design space by sampling from a uniform distribution of \(x\) in the sparsest bin, Memory.generate_sparse().

Main loop

Once the memory structures are in place, it seems natural to make a TabuSearch class to hold all three memories, the objective function callable, algorithm parameters, etc.

class TabuSearch:
    def __init__(self, objective, constraint, nx, ny, tol):
        """Maximise an objective function using Tabu search."""

        # Store objective and constraint functions
        self.objective = objective
        self.constraint = constraint

        # Store tolerance on x
        self.tol  = tol

        # Default memory sizes
        self.n_short = 20
        self.n_long = 20000
        self.nx = nx
        self.ny = ny
        self.n_med = 2000 if ny>1 else 10

        # Default iteration counters
        self.i_diversify = 10
        self.i_intensify = 20
        self.i_restart = 40
        self.i_pattern = 2

        # Misc algorithm parameters
        self.x_regions = 3
        self.max_fevals = 20000
        self.fac_restart = 0.5
        self.fac_pattern = 2.0

        # Initialise counters
        self.fevals = 0

        # Initialise memories
        self.mem_short = Memory(nx, ny, self.n_short, self.tol)
        self.mem_med = Memory(nx, ny, self.n_med, self.tol)
        self.mem_long = Memory(nx, ny, self.n_long, self.tol)
        self.mem_all = (self.mem_short, self.mem_med, self.mem_long)

I have skipped over most of the TabuSearch methods, as they are not all that interesting, but the main loop in TabuSearch.seach() is the code embodiment of the algorithm detailed above.

    def search(self, x0, dx):
        """Perform a search with given intial point and step size."""

        # Evaluate the objective at given initial guess point, update memories
        y0 = ts.initial_guess(x0)

        # Main loop, until max evaluations reached or step size below tolerance
        i = 0
        while self.fevals < self.max_fevals and np.any(dx > self.tol):

            # Evaluate objective for all permissible candidate moves
            X, Y = ts.evaluate_moves(x0, dx)

            # Put new results into long-term memory
            self.mem_long.add(X, Y)

            # Put Pareto-equivalent results into medium-term memory
            # Flag true if we successfully added a point
            if self.ny==1:
                flag = self.mem_med.update_best(X, Y)
            else:
                flag = self.mem_med.update_front(X, Y)

            # Reset counter if we added to medium memory, otherwise increment 
            i = 0 if flag else i+1

            # Choose next point based on local search counter
            if i == self.i_restart:
                # RESTART: reduce step sizes and randomly select from
                # medium-term
                dx = dx * self.fac_restart
                if self.ny==1:
                    # Pick the current optimum if scalar objective
                    x1, y1 = self.mem_med.get(0)
                else:
                    # Pick from sparse region of Pareto from if multi-objective
                    x1, y1 = self.mem_med.sample_sparse(self.x_regions)
                i = 0
            elif i == self.i_intensify or X.shape[0] == 0:
                # INTENSIFY: Select a near-optimal point
                x1, y1 = self.mem_med.sample_sparse(self.x_regions)
            elif i == self.i_diversify:
                # DIVERSIFY: Generate a new point in sparse design region
                x1 = self.mem_long.generate_sparse(self.x_regions)
                y1 = self.objective(x1)
            else:
                # Normally, choose the best candidate move
                x1, y1 = ts.select_move(x0, y0, X, Y)
                # Check for a pattern move every i_pattern steps
                if np.mod(i, self.i_pattern):
                    x1 = ts.pattern_move(x0, y0, x1, y1)

            # Add chosen point to short-term list (tabu)
            self.mem_short.add(x1)

            # Update current point before next iteration
            x0, y0 = x1, y1

        # After the loop return current point
        return x0, y0

Test problem

To check everything is working correctly, we can apply the algorithm to a benchmark problem from the literature. This one, taken from Deb et al. (2002), is known as ‘CONSTR’. Minimise, \[ \bv{y} = f(\bv{x}) = \begin{bmatrix} x_1 \\ (1+x_2)/x_1\end{bmatrix}\ , \] subject to the constraints, \[ \begin{aligned} &0.1 \le x_1 \le 1.0 \ , \\ &0 \le x_2 \le 5 \ , \\ &9x_1 + x_2 \ge 6 \ , \\ &9x_1 - x_2 \ge 1 \ . \end{aligned} \]

In code, these look like,

def objective(x):
    return np.column_stack((x[:, 0], (1.0 + x[:, 1]) / x[:, 0]))

def constrain(x):
    return np.all(
        (
            x[:, 1] + 9.0 * x[:, 0] >= 6.0,
            -x[:, 1] + 9.0 * x[:, 0] >= 1.0,
            x[:, 0] >= 0.0,
            x[:, 0] <= 1.0,
            x[:, 1] >= 0.0,
            x[:, 1] <= 5.0,
        ),
        axis=0,
    )

and we can run the search with just a few more lines,

x0 = np.atleast_2d([.8,6.])  # Initial guess
dx = np.array([0.2,2.])  # Initial step sizes
tol = dx/64.  # Tolerance to halt optimisation
nx = 2  # Number of design parameters
ny = 2  # Number of objectives
ts = TabuSearch(objective, constrain, nx, ny, tol)
x_opt, y_opt = ts.search(x0, dx)

By plotting the medium and long term memories at each loop iteration, we can observe the progress of the optimiser exploring the Pareto front.

Progress of the optimiser with the 'CONSTR' test problem.

Epilogue

I discussed why the Tabu search is particularly suitable for aerodynamic optimisation problems, outlined the algorithm and presented the more interesting parts of my implementation. This post might be useful as a brief introduction for those new to the topic. More importantly, optimisation is a good tool to have in my toolbox, and by coding it up and writing this post the algorithm is ready for me to deploy on new problems in the future.