Solving the RE problems¶

The RE problems are a suite of real-world engineering design problems collected and put into a common form by Tanabe and Ishibuchi. Each one is bound-constrained: where the original design problem had constraints, the suite folds their total violation into an extra objective, so an unconstrained solver can be pointed straight at it.

This guide solves the ones DESDEO ships, and checks the result against the approximated Pareto fronts published with the suite.

The problems and where they come from¶

The suite is:

Tanabe, R., & Ishibuchi, H. (2020). An easy-to-use real-world multi-objective optimization problem suite. Applied Soft Computing, 89, 106078. https://doi.org/10.1016/j.asoc.2020.106078

Each problem is drawn from an earlier paper, which is the reference to cite for the design problem itself:

Problem	Design problem	Original source
re21	Four bar truss	Cheng & Li (1999), Engineering Optimization 31(5), 641–661. 10.1080/03052159908941390
re22	Reinforced concrete beam	Amir & Hasegawa (1989), Journal of Structural Engineering 115(3), 626–646. 10.1061/(ASCE)0733-9445(1989)115:3(626)
re23	Pressure vessel	Kannan & Kramer (1994), Journal of Mechanical Design 116(2), 405–411. 10.1115/1.2919393
re24	Hatch cover	Amir & Hasegawa (1989), as above.
re31	Two bar truss	Coello Coello & Pulido (2005), Structural and Multidisciplinary Optimization 30(5), 388–403. 10.1007/s00158-005-0527-z
re32	Welded beam	Ray & Liew (2002), Engineering Optimization 34(2), 141–153. 10.1080/03052150210915
re33	Disc brake	Ray & Liew (2002), as above.
re34	Vehicle crashworthiness	Liao, Li, Yang, Zhang & Li (2008), Structural and Multidisciplinary Optimization 35(6), 561–569. 10.1007/s00158-007-0163-x
re37	Rocket injector	Vaidyanathan, Tucker, Papila & Shyy (2003), 41st AIAA Aerospace Sciences Meeting. 10.2514/6.2003-296
re41	Car side impact	Jain & Deb (2014), IEEE Transactions on Evolutionary Computation 18(4), 602–622. 10.1109/TEVC.2013.2281534
re42	Conceptual marine design	Parsons & Scott (2004), Journal of Ship Research 48(1), 61–76. 10.5957/jsr.2004.48.1.61
re61	Water resource planning	Ray, Tai & Seow (2001), Engineering Optimization 33(4), 399–424. 10.1080/03052150108940926

re22, re23 and re25 have discrete or integer variables and need the mixed-integer algorithm templates, so they are left out below, as is re24, which DESDEO keeps in its constrained form.

In [1]:

import io

import matplotlib.pyplot as plt
import numpy as np
import requests

from desdeo.emo import algorithms
from desdeo.emo.options.termination import MaxEvaluationsTerminatorOptions
from desdeo.problem.testproblems import re21, re31, re32, re33, re34, re37, re41, re42, re61

problems = {
    "re21": re21(),
    "re31": re31(),
    "re32": re32(),
    "re33": re33(),
    "re34": re34(),
    "re37": re37(),
    "re41": re41(),
    "re42": re42(),
    "re61": re61(),
}

# The RE suite orders the rocket injector's objectives as the injector face temperature, the distance
# from the inlet and the post tip temperature; DESDEO's re37 lists the last two the other way round.
# The three functions themselves are the same, so the columns only need matching up by meaning before
# anything is compared against a published front position by position.
re_column_order = {"re37": ["TF_max", "Xcc_max", "TT_max"]}


def objective_columns(name: str) -> list[str]:
    """Objective symbols of a problem, in the order the RE suite lists them."""
    return re_column_order.get(name, [objective.symbol for objective in problems[name].objectives])


print(f"{'problem':<8}{'variables':>10}{'objectives':>12}")
for name, problem in problems.items():
    print(f"{name:<8}{len(problem.variables):>10}{len(problem.objectives):>12}")

import io import matplotlib.pyplot as plt import numpy as np import requests from desdeo.emo import algorithms from desdeo.emo.options.termination import MaxEvaluationsTerminatorOptions from desdeo.problem.testproblems import re21, re31, re32, re33, re34, re37, re41, re42, re61 problems = { "re21": re21(), "re31": re31(), "re32": re32(), "re33": re33(), "re34": re34(), "re37": re37(), "re41": re41(), "re42": re42(), "re61": re61(), } # The RE suite orders the rocket injector's objectives as the injector face temperature, the distance # from the inlet and the post tip temperature; DESDEO's re37 lists the last two the other way round. # The three functions themselves are the same, so the columns only need matching up by meaning before # anything is compared against a published front position by position. re_column_order = {"re37": ["TF_max", "Xcc_max", "TT_max"]} def objective_columns(name: str) -> list[str]: """Objective symbols of a problem, in the order the RE suite lists them.""" return re_column_order.get(name, [objective.symbol for objective in problems[name].objectives]) print(f"{'problem':<8}{'variables':>10}{'objectives':>12}") for name, problem in problems.items(): print(f"{name:<8}{len(problem.variables):>10}{len(problem.objectives):>12}")

problem  variables  objectives
re21             4           2
re31             3           3
re32             4           3
re33             4           3
re34             5           3
re37             4           3
re41             7           4
re42             6           4
re61             3           6

Solving them¶

One call each. NSGA-III is used throughout because it copes with all of the objective counts here, from two up to six.

What is plotted below is the archive rather than the final population. DESDEO keeps a non-dominated archive of everything a run has seen, on the object emo_constructor returns alongside the solver, and it holds far more of the front than the last generation does: a few thousand solutions against a population of about a hundred.

In [2]:

archives = {}

for name, problem in problems.items():
    options = algorithms.nsga3_options()
    options.template.seed = 0
    options.template.termination = MaxEvaluationsTerminatorOptions(max_evaluations=20000)

    solver, extras = algorithms.emo_constructor(emo_options=options, problem=problem)
    final_population = solver().optimal_outputs

    archives[name] = extras.archive.solutions[objective_columns(name)]
    print(f"{name}: {len(archives[name]):>6} archived, {len(final_population):>4} in the final population")

archives = {} for name, problem in problems.items(): options = algorithms.nsga3_options() options.template.seed = 0 options.template.termination = MaxEvaluationsTerminatorOptions(max_evaluations=20000) solver, extras = algorithms.emo_constructor(emo_options=options, problem=problem) final_population = solver().optimal_outputs archives[name] = extras.archive.solutions[objective_columns(name)] print(f"{name}: {len(archives[name]):>6} archived, {len(final_population):>4} in the final population")

re21:   2310 archived,  100 in the final population

re31:   9387 archived,   91 in the final population

re32:   9232 archived,   91 in the final population

re33:   6225 archived,   91 in the final population

/home/runner/work/DESDEO/DESDEO/desdeo/emo/operators/selection.py:1184: RuntimeWarning: divide by zero encountered in divide
  intercepts = 1 / plane

re34:   3886 archived,   91 in the final population

re37:   5138 archived,   91 in the final population

re41:   5879 archived,   84 in the final population

re42:   5383 archived,   84 in the final population

re61:  12984 archived,   56 in the final population

The reference fronts¶

The suite publishes an approximated Pareto front for every problem, obtained by pooling long runs of several algorithms. These are what an implementation should be checked against: if DESDEO's version of a problem were wrong, its archive would sit somewhere else entirely.

The files are fetched from the suite's repository rather than shipped with DESDEO.

In [3]:

REFERENCE_URL = (
    "https://raw.githubusercontent.com/ryojitanabe/reproblems/master/approximated_Pareto_fronts/reference_points_{}.dat"
)

references = {}
for name in problems:
    try:
        response = requests.get(REFERENCE_URL.format(name.upper()), timeout=60)
        response.raise_for_status()
        references[name] = np.loadtxt(io.StringIO(response.text))
    except requests.RequestException as error:
        # The comparison is worth having, but not worth failing a docs build over.
        print(f"{name}: reference front unavailable ({error})")

for name, reference in references.items():
    print(f"{name}: reference front has {len(reference):>5} points")

REFERENCE_URL = ( "https://raw.githubusercontent.com/ryojitanabe/reproblems/master/approximated_Pareto_fronts/reference_points_{}.dat" ) references = {} for name in problems: try: response = requests.get(REFERENCE_URL.format(name.upper()), timeout=60) response.raise_for_status() references[name] = np.loadtxt(io.StringIO(response.text)) except requests.RequestException as error: # The comparison is worth having, but not worth failing a docs build over. print(f"{name}: reference front unavailable ({error})") for name, reference in references.items(): print(f"{name}: reference front has {len(reference):>5} points")

re21: reference front has  1000 points
re31: reference front has  1500 points
re32: reference front has  1500 points
re33: reference front has  1500 points
re34: reference front has  1500 points
re37: reference front has  1500 points
re41: reference front has  2000 points
re42: reference front has  1999 points
re61: reference front has  2999 points

Comparing the archives against the reference fronts¶

Two objectives plot directly, three fit in a 3-D scatter, and anything above that goes into parallel coordinates. The reference front is drawn in grey underneath DESDEO's archive in each case.

In [4]:

def opaque_legend(ax, **kwargs) -> None:
    """Add a legend whose keys are solid, however transparent the points they stand for are."""
    legend = ax.legend(**kwargs)
    for handle in legend.legend_handles:
        handle.set_alpha(1.0)


name = "re21"
front, reference = archives[name], references.get(name)
x, y = front.columns

fig, ax = plt.subplots(figsize=(5.5, 4.5), constrained_layout=True)
if reference is not None:
    ax.scatter(
        reference[:, 0], reference[:, 1], s=22, color="0.55", alpha=0.35, edgecolor="none", label="reference front"
    )
ax.scatter(front[x], front[y], s=7, alpha=0.35, edgecolor="none", label="DESDEO archive")
ax.set(xlabel=x, ylabel=y, title="re21: four bar truss")
opaque_legend(ax)
plt.show()

def opaque_legend(ax, **kwargs) -> None: """Add a legend whose keys are solid, however transparent the points they stand for are.""" legend = ax.legend(**kwargs) for handle in legend.legend_handles: handle.set_alpha(1.0) name = "re21" front, reference = archives[name], references.get(name) x, y = front.columns fig, ax = plt.subplots(figsize=(5.5, 4.5), constrained_layout=True) if reference is not None: ax.scatter( reference[:, 0], reference[:, 1], s=22, color="0.55", alpha=0.35, edgecolor="none", label="reference front" ) ax.scatter(front[x], front[y], s=7, alpha=0.35, edgecolor="none", label="DESDEO archive") ax.set(xlabel=x, ylabel=y, title="re21: four bar truss") opaque_legend(ax) plt.show()

In [5]:

three_objective = ["re31", "re32", "re33", "re34", "re37"]

fig, axes = plt.subplots(2, 3, figsize=(13, 8.5), subplot_kw={"projection": "3d"}, constrained_layout=True)
for index, name in enumerate(three_objective):
    ax = axes.flat[index]
    front, reference = archives[name], references.get(name)
    x, y, z = front.columns

    if reference is not None:
        ax.scatter(
            reference[:, 0],
            reference[:, 1],
            reference[:, 2],
            s=8,
            color="0.55",
            alpha=0.25,
            edgecolor="none",
            label="reference",
        )
    ax.scatter(front[x], front[y], front[z], s=6, alpha=0.3, edgecolor="none", label="DESDEO")

    ax.set(xlabel=x, ylabel=y, zlabel=z, title=name)
    # Every objective here is minimised, so look into the box from the side of the maxima: the
    # corner where all three are smallest sits at the back, and the front faces the reader.
    ax.view_init(elev=25, azim=45)
    ax.tick_params(labelsize=7, pad=0)
    # Pull the axes in a little so that the tick labels are not clipped by the panel edge.
    ax.set_box_aspect(None, zoom=0.85)


opaque_legend(axes.flat[0], loc="upper left", fontsize=8)
for ax in axes.flat[len(three_objective) :]:
    ax.set_axis_off()
plt.show()

three_objective = ["re31", "re32", "re33", "re34", "re37"] fig, axes = plt.subplots(2, 3, figsize=(13, 8.5), subplot_kw={"projection": "3d"}, constrained_layout=True) for index, name in enumerate(three_objective): ax = axes.flat[index] front, reference = archives[name], references.get(name) x, y, z = front.columns if reference is not None: ax.scatter( reference[:, 0], reference[:, 1], reference[:, 2], s=8, color="0.55", alpha=0.25, edgecolor="none", label="reference", ) ax.scatter(front[x], front[y], front[z], s=6, alpha=0.3, edgecolor="none", label="DESDEO") ax.set(xlabel=x, ylabel=y, zlabel=z, title=name) # Every objective here is minimised, so look into the box from the side of the maxima: the # corner where all three are smallest sits at the back, and the front faces the reader. ax.view_init(elev=25, azim=45) ax.tick_params(labelsize=7, pad=0) # Pull the axes in a little so that the tick labels are not clipped by the panel edge. ax.set_box_aspect(None, zoom=0.85) opaque_legend(axes.flat[0], loc="upper left", fontsize=8) for ax in axes.flat[len(three_objective) :]: ax.set_axis_off() plt.show()

Matplotlib has no parallel coordinates plot of its own, so here is a short one. Every axis is scaled to the range the two fronts cover together, so the lines are directly comparable, and the tick labels keep the original values. Both sets are thinned to a few hundred lines, which is about as many as the eye can follow.

In [6]:

def parallel_coordinates(ax, front: np.ndarray, reference: np.ndarray | None, labels: list[str], title: str) -> None:
    """Draw one line per solution, with the reference front behind DESDEO's archive."""
    rng = np.random.default_rng(0)
    scale_from = front if reference is None else np.vstack([front, reference])
    lower, upper = scale_from.min(axis=0), scale_from.max(axis=0)
    # A constant objective would divide by zero, and every solution sits at the same height on it.
    spread = np.where(upper > lower, upper - lower, 1.0)
    positions = np.arange(front.shape[1])

    def draw(values: np.ndarray, colour: str, alpha: float, label: str) -> None:
        if len(values) > 300:
            values = values[rng.choice(len(values), 300, replace=False)]
        for index, row in enumerate((values - lower) / spread):
            ax.plot(positions, row, color=colour, alpha=alpha, linewidth=0.7, label=label if index == 0 else None)

    if reference is not None:
        draw(reference, "0.55", 0.3, "reference front")
    draw(front, "tab:blue", 0.25, "DESDEO archive")

    for position in positions:
        ax.axvline(position, color="0.85", linewidth=0.8, zorder=0)
    ax.set(xticks=positions, xticklabels=labels, yticks=[0, 1], ylim=(-0.05, 1.05), title=title)
    ax.set_yticklabels(["min", "max"])
    for position, low, high in zip(positions, lower, upper, strict=True):
        ax.annotate(f"{low:.3g}", (position, 0), textcoords="offset points", xytext=(0, -14), ha="center", fontsize=7)
        ax.annotate(f"{high:.3g}", (position, 1), textcoords="offset points", xytext=(0, 6), ha="center", fontsize=7)


many_objective = ["re41", "re42", "re61"]

fig, axes = plt.subplots(len(many_objective), 1, figsize=(7.5, 3.4 * len(many_objective)), constrained_layout=True)
for ax, name in zip(axes, many_objective, strict=True):
    parallel_coordinates(ax, archives[name].to_numpy(), references.get(name), archives[name].columns, name)
opaque_legend(axes[0], loc="upper right", fontsize=8)
plt.show()

def parallel_coordinates(ax, front: np.ndarray, reference: np.ndarray | None, labels: list[str], title: str) -> None: """Draw one line per solution, with the reference front behind DESDEO's archive.""" rng = np.random.default_rng(0) scale_from = front if reference is None else np.vstack([front, reference]) lower, upper = scale_from.min(axis=0), scale_from.max(axis=0) # A constant objective would divide by zero, and every solution sits at the same height on it. spread = np.where(upper > lower, upper - lower, 1.0) positions = np.arange(front.shape[1]) def draw(values: np.ndarray, colour: str, alpha: float, label: str) -> None: if len(values) > 300: values = values[rng.choice(len(values), 300, replace=False)] for index, row in enumerate((values - lower) / spread): ax.plot(positions, row, color=colour, alpha=alpha, linewidth=0.7, label=label if index == 0 else None) if reference is not None: draw(reference, "0.55", 0.3, "reference front") draw(front, "tab:blue", 0.25, "DESDEO archive") for position in positions: ax.axvline(position, color="0.85", linewidth=0.8, zorder=0) ax.set(xticks=positions, xticklabels=labels, yticks=[0, 1], ylim=(-0.05, 1.05), title=title) ax.set_yticklabels(["min", "max"]) for position, low, high in zip(positions, lower, upper, strict=True): ax.annotate(f"{low:.3g}", (position, 0), textcoords="offset points", xytext=(0, -14), ha="center", fontsize=7) ax.annotate(f"{high:.3g}", (position, 1), textcoords="offset points", xytext=(0, 6), ha="center", fontsize=7) many_objective = ["re41", "re42", "re61"] fig, axes = plt.subplots(len(many_objective), 1, figsize=(7.5, 3.4 * len(many_objective)), constrained_layout=True) for ax, name in zip(axes, many_objective, strict=True): parallel_coordinates(ax, archives[name].to_numpy(), references.get(name), archives[name].columns, name) opaque_legend(axes[0], loc="upper right", fontsize=8) plt.show()

The archives track the reference fronts on every problem. Where DESDEO's lines stop short of the grey ones — the top of re31's and re33's violation axis, for instance — the reference front reaches extremes that a single 20 000-evaluation run does not, which is a matter of budget rather than a disagreement about what the problem is.