CEC 2006 g17

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. g17 is the seventeenth, the report's equations 35 and 36 (page 10). The report takes it from Himmelblau (1972, Applied Nonlinear Programming, McGraw-Hill).

There are 6 variables:

x1 in [0, 400]    x2 in [0, 1000]    x3, x4 in [340, 420]    x5 in [−1000, 1000]
x6 in [0, 0.5236]

The objective is f(x) = f1(x1) + f2(x2), with two piecewise linear parts:

f1(x1) = 30 x1 for x1 < 300,  31 x1 from 300
f2(x2) = 28 x2 for x2 < 100,  29 x2 for 100 ≤ x2 < 200,  30 x2 from 200

The report defines f1 for x1 < 400 and f2 for x2 < 1000, but its bounds include 400 and 1000: genoxide's G17 extends the last pieces to them. There are four equality constraints:

h1 = −x1 + 300 − (x3 x4 / 131.078) cos(1.48477 − x6) + (0.90798 x3² / 131.078) cos(1.47588) = 0
h2 = −x2 − (x3 x4 / 131.078) cos(1.48477 + x6) + (0.90798 x4² / 131.078) cos(1.47588) = 0
h3 = −x5 − (x3 x4 / 131.078) sin(1.48477 + x6) + (0.90798 x4² / 131.078) sin(1.47588) = 0
h4 = 200 − (x3 x4 / 131.078) sin(1.48477 − x6) + (0.90798 x3² / 131.078) sin(1.47588) = 0

h1, h2 and h3 fix x1, x2 and x5 from x3, x4 and x6, and h4 ties those three together. The feasible set is a surface: two degrees of freedom left of six. 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 G17: the surface becomes a thin layer.

The best known value

The report prints f(x) = 8853.53967480648, but its x evaluates to 8853.53401643571 under its own definition. genoxide's docs trace the printed value to the organizers' code, which evaluates f1 and f2 at the right-hand sides of h1 and h2, the values that x1 and x2 would take with the equalities met exactly, instead of at x1 and x2. Within the tolerance, those differ from x1 and x2 by h1 and h2. At x*, h1 is 0.0000953 and h2 0.0001, and 30 · 0.0000953 + 28 · 0.0001 = 0.00566, the gap between the two values.

Later papers give 8853.5338748065, a little lower. The report's x reaches it with x1 lowered to 201.78446249355, where h1 reaches the tolerance. genoxide uses it as f; its source is unverified, and it isn't proven optimal. It lies at x = (201.78446249355, 99.9999999999999005, 383.071034852773266, 420, −10.9076584514292652, 0.0731482312084287128).

What makes it hard

The objective jumps. f1 jumps by 300 where x1 reaches 300, and f2 by 100 where x2 reaches 100 and by 200 where it reaches 200. A step across a jump changes f far more than steps elsewhere.

The best known solution sits on the edge of a jump: x2 is just below 100, on the piece 28 x2. A step to x2 = 100 costs 100. x4 is at its upper bound 420, and the equalities h1, h2 and h4 are on the edge of the tolerance. The optimum is a corner of the surface, pressed against a jump, a bound and the tolerance at once.

The feasible layer is thin: of 10 million random points in the box, none was feasible. A search must first reach it from outside, guided only by the violation, and then move along a curved surface in small steps.

And there's a local optimum. At f = 8927.5917, 74.06 above f*, x1 ≈ 107.8 and x2 ≈ 196.3, on the pieces 30 x1 and 29 x2, in another region of the surface. A search that settles there stays: the runs of SHADE that reached it never left it within the budget. With Deb's feasibility rules alone, most runs of every algorithm tried here end there.

Representation

A Real genome of 6 genes, x1 to x6, within the report's bounds. genoxide's problems::cec2006::G17 is the fitness: the value f(x) and the total constraint violation, the sum of max(0, |h| − 0.0001) over the four equalities, 0 for a feasible solution.

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 300 and follows Takahama and Sakai's schedule, ε(t) = 300 (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 Deb's rules alone, no algorithm of genoxide solves g17 reliably. With seeds 1 to 25 and the same budget and target:

Algorithm	Met the target	Within 1e-4 of f*
SHADE at an ε level, as here	25	25
SHADE	7	7
SHADE with JADE's strategy, current-to-pbest with p = 0.1	5	5
L-SHADE	8	8
CMA-ES with IPOP restarts	12	14
CMA-ES with BIPOP restarts	15	20
GA, simulated binary crossover and polynomial mutation	0	0

Every run of plain SHADE either met the target or ended at the local optimum, 8927.59, and so did most failed runs of the others, except the GA's. With Deb's rules alone, the first feasible solutions decide where a run ends: the population gathers around them, and whether they lie on the side of x2 < 100 is a matter of luck. With a large ε, much of the box counts as feasible at first, and the value leads the population toward the minimum's side before the layer thins; as ε falls, the population follows the layer down to the surface. The start matters: in 100 to 300 runs each, with ε starting at 100 or 300, 97 to 100 % met the target, at 10, 90 %, and at 1,000 almost none.

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, to 4 significant digits: x2 prints as 100.00, but lies just below 100, on the piece 28 x2, as the sixth line says. The last gives the four equalities. An equality met within the tolerance is "active", else the line gives how far |h| exceeds 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 shows each variable on its range, and each constraint's state. 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. The best's curve begins once the best is feasible, and the median's once half the population is: 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. Measured against the report's printed value, 8853.53967, such a run would be 0.0058 below it.

The recorded run's best, by the ε level, is infeasible for most of the schedule: it uses the tolerance ε allows, and while ε is large, its equalities are far from met. Its violations shrink with ε, and the best is first feasible without it after 147,400 evaluations, when ε is 3·10⁻⁶. It is within 1e-4 of f* after 150,900 evaluations, once ε is 0, and meets its target after 169,300.

The solution is x = (201.8, 100.00, 383.1, 420.0, −10.91, 0.07315), the best known point to 4 digits, with x2 below 100 and x4 at its bound, and all four equalities met.

With seeds 1 to 1,000, every run meets the target, after 167,200 to 178,600 evaluations (at most 171,600 for half of them).

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: 8853.5338748065 (best known, with the equalities met within 0.0001; the report prints 8853.53967480648)

Source: examples/cec2006_g17

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

RustPython

cargo run --release --example cec2006_g17

//! CEC 2006 g17: a piecewise linear function of 6 variables under 4 nonlinear equality
//! constraints, from the CEC 2006 special session on constrained optimization (Liang et al.,
//! 2006). The best known value is 8853.5338748065, with the equalities met within the report's
//! tolerance of 0.0001.
//!
//! genoxide's `G17` 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 300 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, the pieces of the objective it's on, 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_g17
//! ```

mod trace;

use genoxide::prelude::*;
use genoxide::problems::Problem;
use genoxide::problems::cec2006::{EQUALITY_TOLERANCE, G17};
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;
// a constraint within this of its boundary is active
const ACTIVE: f64 = 1e-6;
// the ε level at the start: a violation up to it counts as none
const EPSILON: f64 = 300.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<()> {
    let problem = G17::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 g17, 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, 4)))
        .collect();
    println!("{}", genes.join(", "));
    // the pieces of f1(x1) and f2(x2) that the solution is on
    let f1 = if x[0] < 300.0 { "30 x1" } else { "31 x1" };
    let f2 = if x[1] < 100.0 {
        "28 x2"
    } else if x[1] < 200.0 {
        "29 x2"
    } else {
        "30 x2"
    };
    println!("pieces: f1(x1) = {f1}, f2(x2) = {f2}");
    // h1 to h4, each as its excess over the tolerance, 0 when it's met
    let constraints: Vec<String> = (1..)
        .zip(problem.constraints(x).equalities())
        .map(|(i, &h)| format!("h{i} {}", state((h.abs() - EQUALITY_TOLERANCE).max(0.0))))
        .collect();
    println!("{}", constraints.join(", "));
    trace.write();
    Ok(())
}

// a constraint g(x) <= 0: "active" on its boundary, else its value
fn state(g: f64) -> String {
    if g.abs() <= ACTIVE {
        "active".to_string()
    } else {
        significant(g, 4)
    }
}

// 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_g17/main.py

"""CEC 2006 g17: a piecewise linear function of 6 variables under 4 nonlinear equality
constraints, from the CEC 2006 special session on constrained optimization (Liang et al., 2006).
The best known value is 8853.5338748065, with the equalities met within the report's tolerance of
0.0001.

genoxide's ``G17`` 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 300 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, the pieces of the objective it's on, 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_g17/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
# a constraint within this of its boundary is active
ACTIVE = 1e-6
# the ε level at the start: a violation up to it counts as none
EPSILON = 300.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 state(g):
    """A constraint g(x) <= 0: "active" on its boundary, else its value."""
    return "active" if abs(g) <= ACTIVE else significant(g, 4)


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)


problem = gx.problems.cec2006.G17()
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 g17, 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, 4)}" for i, xi in enumerate(x, 1)))
# the pieces of f1(x1) and f2(x2) that the solution is on
f1 = "30 x1" if x[0] < 300 else "31 x1"
f2 = "28 x2" if x[1] < 100 else "29 x2" if x[1] < 200 else "30 x2"
print(f"pieces: f1(x1) = {f1}, f2(x2) = {f2}")
# h1 to h4, each as its excess over the tolerance, 0 when it's met
h = problem.constraints(result.best_genome).tolist()
print(
    ", ".join(
        f"h{i} {state(max(abs(hi) - EQUALITY_TOLERANCE, 0.0))}" for i, hi in enumerate(h, 1)
    )
)
trace.write()

What it prints, from a seeded run:

SHADE with Deb's feasibility rules at an epsilon level on g17, seed 1
stopped by the target after 169300 evaluations: f(x) - f* < 1e-8, feasible
first feasible after 147400 evaluations, f(x) - f* <= 1e-4 after 150900
f(x) 8853.53, f* 8853.53 (best known)
x1 201.8, x2 100.00, x3 383.1, x4 420.0, x5 -10.91, x6 0.07315
pieces: f1(x1) = 30 x1, f2(x2) = 28 x2
h1 active, h2 active, h3 active, h4 active