Solving Constrained Problems with Evolutionary Algorithms¶

Every a posteriori evolutionary algorithm in DESDEO accepts constrained problems. This guide shows how to set one up, how to read the results, and — since the four algorithms handle constraints in three different ways — what to watch out for when comparing them.

If you have not used DESDEO's EAs before, read how to use evolutionary algorithms and how to configure them first. This guide assumes the Pydantic interface from the latter.

How DESDEO states a constraint¶

A Constraint is written in standard form, with the expression on the left-hand side and an implied zero on the right. So for a cons_type of <=, a constraint value that is zero or negative is satisfied, and a positive value is the amount by which it is violated. That sign convention is the one thing you need to remember: every selection operator below reads a constraint column exactly that way, and the total violation of a solution is the sum of the positive parts.

Constraints are evaluated alongside the objectives and land in the same output DataFrame, so nothing special is needed to see them.

In [1]:

import numpy as np
import polars as pl

from desdeo.emo import algorithms
from desdeo.emo.operators.selection import ParameterAdaptationStrategy
from desdeo.emo.options.termination import MaxEvaluationsTerminatorOptions
from desdeo.problem.testproblems import car_side_impact

# A real engineering problem: minimize weight, door velocity and passenger load, subject to ten
# safety and manufacturing limits. All seven variables are continuous.
problem = car_side_impact(three_obj=True)

objectives = [objective.symbol for objective in problem.objectives]
constraints = [constraint.symbol for constraint in problem.constraints]

print(f"{len(problem.variables)} variables, {len(objectives)} objectives, {len(constraints)} constraints")
print("constraint types:", {constraint.symbol: constraint.cons_type for constraint in problem.constraints})

import numpy as np import polars as pl from desdeo.emo import algorithms from desdeo.emo.operators.selection import ParameterAdaptationStrategy from desdeo.emo.options.termination import MaxEvaluationsTerminatorOptions from desdeo.problem.testproblems import car_side_impact # A real engineering problem: minimize weight, door velocity and passenger load, subject to ten # safety and manufacturing limits. All seven variables are continuous. problem = car_side_impact(three_obj=True) objectives = [objective.symbol for objective in problem.objectives] constraints = [constraint.symbol for constraint in problem.constraints] print(f"{len(problem.variables)} variables, {len(objectives)} objectives, {len(constraints)} constraints") print("constraint types:", {constraint.symbol: constraint.cons_type for constraint in problem.constraints})

7 variables, 3 objectives, 10 constraints
constraint types: {'g_1': <ConstraintTypeEnum.LTE: '<='>, 'g_2': <ConstraintTypeEnum.LTE: '<='>, 'g_3': <ConstraintTypeEnum.LTE: '<='>, 'g_4': <ConstraintTypeEnum.LTE: '<='>, 'g_5': <ConstraintTypeEnum.LTE: '<='>, 'g_6': <ConstraintTypeEnum.LTE: '<='>, 'g_7': <ConstraintTypeEnum.LTE: '<='>, 'g_8': <ConstraintTypeEnum.LTE: '<='>, 'g_9': <ConstraintTypeEnum.LTE: '<='>, 'g_10': <ConstraintTypeEnum.LTE: '<='>}

In [2]:

def total_violation(outputs: pl.DataFrame) -> np.ndarray:
    """Total constraint violation per solution. Zero means feasible."""
    return np.maximum(outputs[constraints].to_numpy(), 0.0).sum(axis=1)

def total_violation(outputs: pl.DataFrame) -> np.ndarray: """Total constraint violation per solution. Zero means feasible.""" return np.maximum(outputs[constraints].to_numpy(), 0.0).sum(axis=1)

Running the Evolutionary Algorithms¶

Once the problem is set up, the evolutionary algorithms are run in the same way as unconstrained problems.

In [3]:

builders = {
    "RVEA": algorithms.rvea_options,
    "NSGA-III": algorithms.nsga3_options,
    "NSGA-II": algorithms.nsga2_options,
    "IBEA": algorithms.ibea_options,
}

fronts = {}
summary = []

for name, build in builders.items():
    options = build()
    options.template.seed = 0
    options.template.verbosity = 2
    options.template.termination = MaxEvaluationsTerminatorOptions(max_evaluations=20000)
    if name == "RVEA":
        options.template.selection.parameter_adaptation_strategy = ParameterAdaptationStrategy.FUNCTION_EVALUATION_BASED

    solver, _extras = algorithms.emo_constructor(emo_options=options, problem=problem)
    outputs = solver().optimal_outputs

    violation = total_violation(outputs)
    fronts[name] = outputs.filter(violation <= 0)
    summary.append(
        {
            "algorithm": name,
            "population": len(outputs),
            "feasible": int((violation <= 0).sum()),
            "worst violation": float(violation.max()),
        }
    )

print(pl.DataFrame(summary))

builders = { "RVEA": algorithms.rvea_options, "NSGA-III": algorithms.nsga3_options, "NSGA-II": algorithms.nsga2_options, "IBEA": algorithms.ibea_options, } fronts = {} summary = [] for name, build in builders.items(): options = build() options.template.seed = 0 options.template.verbosity = 2 options.template.termination = MaxEvaluationsTerminatorOptions(max_evaluations=20000) if name == "RVEA": options.template.selection.parameter_adaptation_strategy = ParameterAdaptationStrategy.FUNCTION_EVALUATION_BASED solver, _extras = algorithms.emo_constructor(emo_options=options, problem=problem) outputs = solver().optimal_outputs violation = total_violation(outputs) fronts[name] = outputs.filter(violation <= 0) summary.append( { "algorithm": name, "population": len(outputs), "feasible": int((violation <= 0).sum()), "worst violation": float(violation.max()), } ) print(pl.DataFrame(summary))

shape: (4, 4)
┌───────────┬────────────┬──────────┬─────────────────┐
│ algorithm ┆ population ┆ feasible ┆ worst violation │
│ ---       ┆ ---        ┆ ---      ┆ ---             │
│ str       ┆ i64        ┆ i64      ┆ f64             │
╞═══════════╪════════════╪══════════╪═════════════════╡
│ RVEA      ┆ 78         ┆ 78       ┆ 0.0             │
│ NSGA-III  ┆ 91         ┆ 91       ┆ 0.0             │
│ NSGA-II   ┆ 100        ┆ 100      ┆ 0.0             │
│ IBEA      ┆ 100        ┆ 100      ┆ 0.0             │
└───────────┴────────────┴──────────┴─────────────────┘

Three of the four end fully feasible. RVEA usually keeps a handful of infeasible solutions due to how it handles constraints, which the next section explains.

Note also that RVEA's population size fluctuates and sits below NSGA-III's, even though both are built from the same number of reference vectors. RVEA keeps at most one survivor per reference vector and some vectors end up with nothing associated to them, so its realised population size is a result of the search rather than a setting. This is expected behaviour and is documented in the original paper.

Plotting the feasible fronts together¶

Only feasible solutions are plotted. An infeasible solution can sit anywhere, including well past the feasible front.

In [4]:

import plotly.graph_objects as go

fig = go.Figure()
for name, front in fronts.items():
    fig.add_scatter3d(
        x=front["f_1"],
        y=front["f_2"],
        z=front["f_3"],
        mode="markers",
        marker={"size": 3, "opacity": 0.75},
        name=f"{name} ({len(front)} feasible)",
    )

fig.update_layout(
    title="Feasible front approximations for the car-side impact problem",
    scene={
        "xaxis_title": "f_1: weight",
        "yaxis_title": "f_2: door velocity",
        "zaxis_title": "f_3: passenger load",
    },
    legend={"itemsizing": "constant"},
)
fig.show(renderer="notebook", include_plotlyjs="cdn")

import plotly.graph_objects as go fig = go.Figure() for name, front in fronts.items(): fig.add_scatter3d( x=front["f_1"], y=front["f_2"], z=front["f_3"], mode="markers", marker={"size": 3, "opacity": 0.75}, name=f"{name} ({len(front)} feasible)", ) fig.update_layout( title="Feasible front approximations for the car-side impact problem", scene={ "xaxis_title": "f_1: weight", "yaxis_title": "f_2: door velocity", "zaxis_title": "f_3: passenger load", }, legend={"itemsizing": "constant"}, ) fig.show(renderer="notebook", include_plotlyjs="cdn")

The three constraint-handling rules¶

The algorithms agree on the sign convention and on nothing else. Knowing which rule you are using matters, because it decides what happens while the population is still largely infeasible — which is most of the run on a tightly constrained problem.

algorithm	rule	source
NSGA-II	constrained domination — a feasible solution beats an infeasible one, two infeasible solutions are compared by total violation, two feasible ones by ordinary Pareto dominance	Deb, Pratap, Agarwal & Meyarivan (2002)
NSGA-III, IBEA	feasibility-first — while feasible solutions can fill the population the infeasible ones are discarded outright; while they cannot, every feasible solution survives and the remaining places go to the least infeasible	Deb & Jain (2014); Zitzler & Künzli (2004)
RVEA	feasibility-first within each reference vector — each vector keeps its best feasible solution by APD, or, if it has no feasible solution associated with it, its least infeasible one	Cheng, Jin, Olhofer & Sendhoff (2016)

RVEA's rule is the local one, and it is why RVEA finishes with infeasible solutions in the population above. A reference vector pointing into a region where nothing feasible exists still has to keep something, so it keeps the least infeasible candidate it saw..

This behaviour isn't necessarily undesirable. Keeping infeasible solutions may help preserve diversity and can help the search cross a narrow infeasible barrier. However, it also means the final population needs filtering before decision making.

One consequence for benchmarking: comparing algorithms on a constrained problem compares their constraint handling as much as their selection. If that is not the comparison you meant to make, either say so, or hold the constraint handling fixed and vary only the part you are studying.

References¶

1. K. Deb, A. Pratap, S. Agarwal, and T. Meyarivan, "A fast and elitist multiobjective genetic algorithm: NSGA-II," IEEE Transactions on Evolutionary Computation, vol. 6, no. 2, pp. 182-197, 2002.
2. K. Deb and H. Jain, "An evolutionary many-objective optimization algorithm using reference-point-based nondominated sorting approach, part I: solving problems with box constraints," IEEE Transactions on Evolutionary Computation, vol. 18, no. 4, pp. 577-601, 2014.
3. E. Zitzler and S. Künzli, "Indicator-based selection in multiobjective search," in Parallel Problem Solving from Nature - PPSN VIII, LNCS vol. 3242, Springer, 2004, pp. 832-842.
4. R. Cheng, Y. Jin, M. Olhofer, and B. Sendhoff, "A reference vector guided evolutionary algorithm for many-objective optimization," IEEE Transactions on Evolutionary Computation, vol. 20, no. 5, pp. 773-791, 2016.
5. H. Jain and K. Deb, "An evolutionary many-objective optimization algorithm using reference-point-based nondominated sorting approach, part II: handling constraints and extending to an adaptive approach," IEEE Transactions on Evolutionary Computation, vol. 18, no. 4, pp. 602-622, 2014.