CEC 2006 g13

The problem

The CEC 2006 special session on constrained optimization (Liang et al., 2006) collected 24 problems, g01 to g24, from the literature, with their best known solutions and rules for comparing algorithms. g13 is the thirteenth, the report's equations 26 and 27 (pages 6-7). The report takes it from Hock and Schittkowski (1981, Test Examples for Nonlinear Programming Codes, Lecture Notes in Economics and Mathematical Systems 187, Springer).

There are 5 variables, x1 and x2 in [−2.3, 2.3], and x3 to x5 in [−3.2, 3.2]. The problem is

minimize   f(x) = exp(x1 x2 x3 x4 x5)
subject to h1(x) = x1² + x2² + x3² + x4² + x5² − 10 = 0
           h2(x) = x2 x3 − 5 x4 x5 = 0
           h3(x) = x1³ + x2³ + 1 = 0

h1 puts x on a sphere of radius √10, h2 ties x4 x5 to x2 x3, and h3 ties x2 to x1. Three equalities in five variables leave a surface of two dimensions. f is least where the product x1 x2 x3 x4 x5 is most negative.

Real-valued samples almost never meet an equality exactly. The report counts an equality as met when |h(x)| ≤ 0.0001, and so does genoxide's G13: the surface becomes a thin shell. The best known value is f = 0.053941514041898, at x = (−1.71714224003, 1.59572124049468, 1.8272502406271, −0.763659881912867, −0.76365986736498), the report's. It meets the equalities within the tolerance only, and genoxide's docs mark it as a best known value, not a proven optimum. As printed, the report's x exceeds the tolerance by 3e-15 in h2, from rounding its digits.

The minimum isn't unique. Changing the signs of two of x3, x4 and x5 leaves the product, h1 and |h2| as they were, so every solution has three twins: x = (−1.717, 1.596, 1.827, 0.764, 0.764) is as good as the report's.

What makes it hard

The feasible region has no volume without the tolerance, and hardly any with it: three shells, each 0.0002 thick in its h, must meet. A random solution is practically never feasible. A search must first reach the shell, guided only by the violation, and then move along it in small steps, since the equalities are nonlinear and a straight step leaves the shell.

The shell also holds a local minimum, with f ≈ 0.4388 at about x = (−0.699, −0.870, 2.790, 0.697, −0.697), and its twins. There, x1 and x2 are both negative, and the product is −0.82 instead of −2.92. A search that has converged there would have to cross worse solutions along the shell to leave it.

Representation

A Real genome of 5 genes, x1 to x5, within the report's bounds. genoxide's problems::cec2006::G13 is the fitness: the value f(x) and the total constraint violation, max(0, |h1| − 0.0001) + max(0, |h2| − 0.0001) + max(0, |h3| − 0.0001), 0 for a feasible solution. The tolerance is the report's, genoxide::problems::cec2006::EQUALITY_TOLERANCE.

genoxide compares fitnesses with Deb's feasibility rules (Deb, 2000, Computer Methods in Applied Mechanics and Engineering 186: 311-338): a feasible solution beats an infeasible one, two feasible ones compare by value, and two infeasible ones by violation. The rules need no penalty weights.

Algorithm

SHADE (Tanabe and Fukunaga, 2013, IEEE CEC 2013: 71-78), a differential evolution that adapts its scale factor and crossover rate from successful trials, with genoxide's defaults: its published population of 100, and a restart after 200 generations without progress.

SHADE compares solutions with Deb's rules at an ε level, the ε constrained method of Takahama and Sakai (2006, "Constrained optimization by the ε constrained differential evolution with gradient-based mutation and feasible elites", IEEE CEC 2006), whose εDE won the CEC 2006 competition: a violation up to ε counts as none. Two solutions within ε compare by value, and the others as Deb's rules have it. ε starts at 20 and follows Takahama and Sakai's schedule, ε(t) = 20 (1 − t / 150,000)⁵ after t evaluations, down to 0 at 150,000. The example applies it with genoxide's features: the fitness function returns the violation beyond ε, and the engine's control lowers ε every 10 generations and scores the population again (reevaluate), since the old violations no longer compare with the new ones. The re-evaluations count in the budget. From 150,000 evaluations on, the rules are Deb's, and the problem the report's.

The run has the report's budget of 500,000 evaluations, and stops once ε is 0 and its best solution is feasible with an error f(x) − f* of at most 1e-8, an absolute error. The report counts a run as successful with an error of at most 1e-4; the example asks for more.

With seeds 1 to 25, in the same budget:

Algorithm	Runs that met the target	Evaluations (median, range)
SHADE at an ε level, as here	25 of 25	150,900 (150,100 to 152,700)
CMA-ES with IPOP restarts	22 of 25	123,808 (19,984 to 443,152)
CMA-ES	7 of 25	57,720 (19,984 to 84,608)
L-SHADE	0 of 25
SHADE	0 of 25

With Deb's rules alone, SHADE and L-SHADE (Tanabe and Fukunaga, 2014, IEEE CEC 2014: 1658-1665), its variant with a shrinking population, build a trial from the difference between two solutions. On a curved shell, that difference points off it, and the trials are infeasible. SHADE ended 12 runs infeasible and 13 with errors from 0.79 to 0.95; L-SHADE ended all 25 feasible, with errors from 0.22 to 0.63. At an ε level, the shell is thick while the population spreads out and converges, and thins as it gathers: its differences shrink with it.

CMA-ES (Hansen and Ostermeier, 2001, Evolutionary Computation 9(2): 159-195) samples a population from a normal distribution, and adapts its mean, step size and covariance matrix, which can learn the directions of the shell. A single run often ends at the local minimum: without restarts, 18 of the 25 runs did, with an error of 0.385. IPOP restarts (Auger and Hansen, 2005, IEEE CEC 2005: 1769-1776) start a new run from a random point, with twice the population, whenever one has converged. In the table, three IPOP runs, with seeds 8, 20 and 24, were still at the local minimum when the budget ran out; with seeds 1 to 100, 13 runs failed, 11 of them there.

Output

The first line names the run. The second gives what stopped it, after how many evaluations, the error f(x) − f and whether the best solution is feasible: "< 1e-8" means the run met its target. The third gives when the best solution, by the ε level, was first feasible without it, and when its error first met the report's criterion of success. The fourth compares f(x) with f, to 6 significant digits. The fifth gives the solution, and the last |h| of each equality, met when it's at most 0.0001. In Python, the fitness function evaluates the problem in Rust, a generation at a time, so both versions print the same.

The page's plot shows each variable on its range, and each constraint's state from its violation, max(0, |h| − 0.0001): violated while ε lets the best solution off the shell, and met from the end of the schedule on. Its curve shows the error f − f* of the best solution, and of the population's median, on a log scale, measured without the ε level: it begins once they are feasible, near the end of the schedule.

The project page plays this run back.

Good results

A good run is feasible and ends within 1e-4 of f, the report's success. A value below f is possible, since f* is only the best known, but the runs here end just above it.

The recorded run's best, by the ε level, lies off the shell for most of the schedule, as far as ε allows. It is first feasible without ε after 147,500 evaluations, already within 1e-4 of f*, and meets its target after 150,200, once ε is 0. The solution is x = (−1.71702, 1.59558, −1.82747, −0.763721, 0.763625), a twin of the report's, with the signs of x3 and x5 changed, and each |h| at 0.0001, on the edge of the tolerance.

With seeds 1 to 1,000, 998 runs meet the target, after 150,100 to 179,800 evaluations (150,900 for half of them). The other 2 end at the local minimum.

Reference: Liang, J. J., Runarsson, T. P., Mezura-Montes, E., Clerc, M., Suganthan, P. N., Coello Coello, C. A. and Deb, K. (2006). Problem Definitions and Evaluation Criteria for the CEC 2006 Special Session on Constrained Real-Parameter Optimization. Technical report, Nanyang Technological University, Singapore.

Known optimum: 0.053941514041898 (best known, with the equalities met within 0.0001)

Source: examples/cec2006_g13

Interactive run: tachsin.gr/projects/genoxide/examples/cec2006-g13

RustPython

cargo run --release --example cec2006_g13

//! CEC 2006 g13: the exponential of a product of 5 variables under 3 nonlinear equality
//! constraints, from the CEC 2006 special session on constrained optimization (Liang et al.,
//! 2006). The best known value is 0.053941514041898, with the equalities met within the report's
//! tolerance of 0.0001.
//!
//! genoxide's `G13` gives the value of a solution and its constraint violation. SHADE, a
//! differential evolution, compares them with Deb's feasibility rules at an ε level (Takahama and
//! Sakai, 2006): a violation up to ε counts as none. ε starts at 20 and falls to 0 over the first
//! 150,000 evaluations, and the population is scored again each time it falls. The run has the
//! report's budget of 500,000 evaluations, and stops once ε is 0 and the error f(x) − f* is at
//! most 1e-8. The example prints the best solution and its constraints.
//!
//! With `GENOXIDE_TRACE=<file>`, it also writes a trace of its run for the plot on the example's
//! page, with `trace.rs`.
//!
//! ```text
//! cargo run --release --example cec2006_g13
//! ```

mod trace;

use genoxide::prelude::*;
use genoxide::problems::Problem;
use genoxide::problems::cec2006::{EQUALITY_TOLERANCE, G13};
use std::sync::Arc;
use std::sync::atomic::{AtomicU64, Ordering};

// the CEC 2006 report's budget of evaluations per run
const BUDGET: u64 = 500_000;
// the run stops once its best is feasible with an error f(x) - f* at most this
const ERROR: f64 = 1e-8;
// the report counts a run as successful once its error is at most this
const SUCCESS: f64 = 1e-4;
// the ε level at the start: a violation up to it counts as none
const EPSILON: f64 = 20.0;
// the evaluations after which ε is 0
const CONTROL: u64 = 150_000;
// ε falls, and the population is scored again, every this many generations
const EVERY: u64 = 10;

// the ε level after `evaluations`: EPSILON (1 - evaluations / CONTROL)^5, then 0
fn epsilon(evaluations: u64) -> f64 {
    if evaluations >= CONTROL {
        return 0.0;
    }
    let rest = 1.0 - evaluations as f64 / CONTROL as f64;
    EPSILON * rest * rest * rest * rest * rest
}

// the ε level in use, shared by the fitness function, the control and the stop condition
#[derive(Clone)]
struct Level(Arc<AtomicU64>);

impl Level {
    fn get(&self) -> f64 {
        f64::from_bits(self.0.load(Ordering::Relaxed))
    }

    fn set(&self, epsilon: f64) {
        self.0.store(epsilon.to_bits(), Ordering::Relaxed);
    }
}

fn main() -> Result<()> {
    // the report's equality tolerance, EQUALITY_TOLERANCE
    let problem = G13::default();
    let optimum = problem.optimum().expect("known");
    let f_star = optimum.value();
    let shade = De::builder(problem.representation())
        .minimize()
        .seed(1)
        .build()?;
    let level = Level(Arc::new(AtomicU64::new(EPSILON.to_bits())));
    // the value, and the violation beyond ε
    let fitness = |x: &Reals| {
        let (value, violation) = problem.evaluate(x);
        (value, (violation - level.get()).max(0.0))
    };
    // the run's target, once ε is 0: a feasible best within ERROR of f*
    let target = {
        let level = level.clone();
        Stop::custom(move |progress| {
            level.get() == 0.0
                && progress.best().is_some_and(|best| {
                    best.is_feasible() && best.score().is_some_and(|value| value <= f_star + ERROR)
                })
        })
    };
    // with GENOXIDE_TRACE=<file>, a trace of the run for the plot on the example's page
    let mut trace = trace::Trace::from_env();
    // the evaluations when the best is first feasible, and when its error first meets the
    // report's criterion of success
    let (mut feasible, mut success) = (None, None);
    let outcome = Engine::new(shade, fitness)
        .stop_when(target.or(Stop::evaluations(BUDGET)))
        .on_generation(|snapshot| {
            let progress = snapshot.progress();
            // the best by the ε level, measured without it
            let (value, violation) = problem.evaluate(snapshot.best().genome());
            if violation == 0.0 {
                feasible.get_or_insert(progress.evaluations());
                if value - f_star <= SUCCESS {
                    success.get_or_insert(progress.evaluations());
                }
            }
            trace.record(snapshot);
        })
        .control(|shade, progress| {
            // every EVERY generations, and once it reaches 0, ε follows its schedule
            let next = epsilon(progress.evaluations());
            let due = progress.generation() % EVERY == 0 || next == 0.0;
            if progress.generation() > 0 && due && next != level.get() {
                level.set(next);
                shade.reevaluate()?;
            }
            Ok(())
        })
        .run()?;

    let best = outcome.best_fitness();
    let value = best.score().expect("valid");
    let x = outcome.best_genome();
    println!("SHADE with Deb's feasibility rules at an epsilon level on g13, seed 1");
    let (stop, error) = if outcome.stop_reason() == StopReason::Custom {
        ("stopped by the target", format!("< {ERROR:.0e}"))
    } else {
        ("stopped", format!("{:.1e}", value - f_star))
    };
    let evaluations = outcome.evaluations();
    let feasibility = if best.is_feasible() {
        "feasible"
    } else {
        "infeasible"
    };
    println!("{stop} after {evaluations} evaluations: f(x) - f* {error}, {feasibility}");
    println!(
        "first feasible after {} evaluations, f(x) - f* <= 1e-4 after {}",
        count(feasible),
        count(success)
    );
    println!(
        "f(x) {}, f* {} ({})",
        significant(value, 6),
        significant(f_star, 6),
        if optimum.is_proven() {
            "proven"
        } else {
            "best known"
        }
    );
    let genes: Vec<String> = (1..)
        .zip(&x[..])
        .map(|(i, xi)| format!("x{i} {}", significant(*xi, 6)))
        .collect();
    println!("{}", genes.join(", "));
    // each equality h = 0 as |h|, met when it's at most 0.0001
    let constraints: Vec<String> = (1..)
        .zip(problem.constraints(x).equalities())
        .map(|(i, h)| {
            let state = if h.abs() <= EQUALITY_TOLERANCE {
                "met"
            } else {
                "violated"
            };
            format!("|h{i}| {:.6} {state}", h.abs())
        })
        .collect();
    println!("{} (|h| <= 0.0001)", constraints.join(", "));
    trace.write();
    Ok(())
}

// the evaluations, or "never"
fn count(evaluations: Option<u64>) -> String {
    evaluations.map_or("never".to_string(), |evaluations| evaluations.to_string())
}

// `digits` significant digits, e.g. 29.9953 or -30665.5 for 6
fn significant(value: f64, digits: i32) -> String {
    let magnitude = value.abs().log10().floor() as i32;
    let decimals = (digits - 1 - magnitude).max(0) as usize;
    format!("{value:.decimals$}")
}

python examples/cec2006_g13/main.py

"""CEC 2006 g13: the exponential of a product of 5 variables under 3 nonlinear equality
constraints, from the CEC 2006 special session on constrained optimization (Liang et al., 2006).
The best known value is 0.053941514041898, with the equalities met within the report's tolerance
of 0.0001.

genoxide's ``G13`` gives the value of a solution and its constraint violation. SHADE, a
differential evolution, compares them with Deb's feasibility rules at an ε level (Takahama and
Sakai, 2006): a violation up to ε counts as none. ε starts at 20 and falls to 0 over the first
150,000 evaluations, and the population is scored again each time it falls. The run has the
report's budget of 500,000 evaluations, and stops once ε is 0 and the error f(x) − f* is at most
1e-8. The example prints the best solution and its constraints. The fitness function evaluates the
problem in Rust, a generation at a time.

With ``GENOXIDE_TRACE=<file>``, it also writes a trace of its run for the plot on the example's
page, with trace.py.

    python examples/cec2006_g13/main.py
"""

import math

import genoxide as gx
import numpy as np
from genoxide.problems.cec2006 import EQUALITY_TOLERANCE

from trace import Trace

# the CEC 2006 report's budget of evaluations per run
BUDGET = 500_000
# the run stops once its best is feasible with an error f(x) - f* at most this
ERROR = 1e-8
# the report counts a run as successful once its error is at most this
SUCCESS = 1e-4
# the ε level at the start: a violation up to it counts as none
EPSILON = 20.0
# the evaluations after which ε is 0
CONTROL = 150_000
# ε falls, and the population is scored again, every this many generations
EVERY = 10


def epsilon(evaluations):
    """The ε level after ``evaluations``: EPSILON (1 - evaluations / CONTROL)^5, then 0."""
    if evaluations >= CONTROL:
        return 0.0
    rest = 1.0 - evaluations / CONTROL
    return EPSILON * rest * rest * rest * rest * rest


def significant(value, digits):
    """``digits`` significant digits, e.g. 29.9953 or -30665.5 for 6."""
    magnitude = math.floor(math.log10(abs(value)))
    return f"{value:.{max(digits - 1 - magnitude, 0)}f}"


def scientific(value, decimals):
    """Scientific notation as Rust writes it, e.g. 1.2e-5 for 1 decimal."""
    mantissa, exponent = f"{value:.{decimals}e}".split("e")
    return f"{mantissa}e{int(exponent)}"


def count(evaluations):
    """The evaluations, or "never"."""
    return "never" if evaluations is None else str(evaluations)


# the report's equality tolerance, EQUALITY_TOLERANCE
problem = gx.problems.cec2006.G13()
optimum = problem.optimum
f_star = optimum.value
shade = gx.De(problem.genome, objective=problem.objective, seed=1)
# the ε level in use
level = {"epsilon": EPSILON}


def fitness(genomes):
    """The values of a generation, and their violations beyond ε."""
    values, violations = problem.evaluate(genomes)
    return values, np.maximum(violations - level["epsilon"], 0.0)


# with GENOXIDE_TRACE=<file>, a trace of the run for the plot on the example's page
trace = Trace(problem)
# the evaluations when the best is first feasible, and when its error first meets the report's
# criterion of success
first = {"feasible": None, "success": None}


def on_generation(progress):
    # the best by the ε level, measured without it
    values, violations = problem.evaluate(progress.best_genome[np.newaxis])
    feasible = violations[0] == 0.0
    if feasible:
        error = values[0] - f_star
        if first["feasible"] is None:
            first["feasible"] = progress.evaluations
        if first["success"] is None and error <= SUCCESS:
            first["success"] = progress.evaluations
    trace.record(progress)
    # the run's target, once ε is 0: a feasible best within ERROR of f*
    return not (level["epsilon"] == 0.0 and feasible and values[0] <= f_star + ERROR)


def control(shade, progress):
    # every EVERY generations, and once it reaches 0, ε follows its schedule
    following = epsilon(progress.evaluations)
    due = progress.generation % EVERY == 0 or following == 0.0
    if progress.generation > 0 and due and following != level["epsilon"]:
        level["epsilon"] = following
        shade.reevaluate()


result = shade.run(
    fitness,
    batch=True,
    evaluations=BUDGET,
    on_generation=on_generation,
    control=control,
)

value = result.best_fitness
print("SHADE with Deb's feasibility rules at an epsilon level on g13, seed 1")
if result.stop_reason == "aborted":
    stop, error = "stopped by the target", f"< {scientific(ERROR, 0)}"
else:
    stop, error = "stopped", scientific(value - f_star, 1)
feasibility = "feasible" if result.violation == 0.0 else "infeasible"
print(f"{stop} after {result.evaluations} evaluations: f(x) - f* {error}, {feasibility}")
print(
    f"first feasible after {count(first['feasible'])} evaluations, f(x) - f* <= 1e-4 after "
    f"{count(first['success'])}"
)
proven = "proven" if optimum.proven else "best known"
print(f"f(x) {significant(value, 6)}, f* {significant(f_star, 6)} ({proven})")
x = result.best_genome.tolist()
print(", ".join(f"x{i} {significant(xi, 6)}" for i, xi in enumerate(x, 1)))
# each equality h = 0 as |h|, met when it's at most 0.0001
equalities = problem.constraints(result.best_genome).tolist()
constraints = [
    f"|h{i}| {abs(h):.6f} {'met' if abs(h) <= EQUALITY_TOLERANCE else 'violated'}"
    for i, h in enumerate(equalities, 1)
]
print(", ".join(constraints) + " (|h| <= 0.0001)")
trace.write()

What it prints, from a seeded run:

SHADE with Deb's feasibility rules at an epsilon level on g13, seed 1
stopped by the target after 150200 evaluations: f(x) - f* < 1e-8, feasible
first feasible after 147500 evaluations, f(x) - f* <= 1e-4 after 147500
f(x) 0.0539415, f* 0.0539415 (best known)
x1 -1.71702, x2 1.59558, x3 -1.82747, x4 -0.763721, x5 0.763625
|h1| 0.000100 met, |h2| 0.000100 met, |h3| 0.000100 met (|h| <= 0.0001)