CEC 2006 g16
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. g16 is the sixteenth, the report's equations 32 to 34 (pages 7 to 10). The report takes it from Himmelblau (1972, Applied Nonlinear Programming, McGraw-Hill).
There are 5 variables:
x1 in [704.4148, 906.3855] x2 in [68.6, 288.88] x3 in [0, 134.75]
x4 in [193, 287.0966] x5 in [25, 84.1988]
The objective isn't one formula. The report computes 17 quantities y1 to y17, and 17 helpers c1 to c17, one after the other, each from the variables and the ones before it. The first few:
y1 = x2 + x3 + 41.6 c1 = 0.024 x4 − 4.62 y2 = 12.5 / c1 + 12
c2 = 0.0003535 x1² + 0.5311 x1 + 0.08705 y2 x1 c3 = 0.052 x1 + 78 + 0.002377 y2 x1
y3 = c2 / c3 y4 = 19 y3 ...
and so on to y17, with divisions and products all along. The objective is a weighted sum of some of them:
f(x) = 0.000117 y14 + 0.1365 + 0.00002358 y13 + 0.000001502 y16 + 0.0321 y12 + 0.004324 y5
+ 0.0001 c15 / c16 + 37.48 y2 / c12 − 0.0000005843 y17
There are 38 inequality constraints, each g(x) ≤ 0. Four are direct:
g1 = 0.28 / 0.72 · y5 − y4 g2 = x3 − 1.5 x2
g3 = 3496 y2 / c12 − 21 g4 = 110.6 + y1 − 62212 / c17
and the other 34 keep each of y1 to y17 between a lower and an upper limit: g5 and g6 for y1,
g7 and g8 for y2, and so on to g37 and g38 for y17. genoxide's G16 has the full chain, from the
report's equation 34.
The best known value is f* = −1.90515525853479, at x = (705.174537070090537, 68.6, 102.899999999999991, 282.324931593660324, 37.5841164258054832). It isn't proven optimal.
What makes it hard
The feasible region is small: of 10 million random points in the box, 2,022 were feasible, about 0.02 %. The chain makes most constraints nonlinear functions of all the variables, and couples them: moving one variable changes many y's at once.
The minimum is a corner where many boundaries meet. At the report's x, five constraints are active: g2, g3, g4, g5 and g36, the upper limit of y16. x2 is also at its lower bound, 68.6. That is more conditions than variables. g2 holds when x3 = 1.5 x2, and g5 when y1 = x2 + x3 + 41.6 reaches its lower limit 213.1, x2 + x3 = 171.5: with x2 at 68.6, both give x3 = 102.9. The report's table 3 counts four active constraints at x; genoxide's docs found five.
Other corners nearby are almost as good. At one of them, a little above f*, nearly every step is blocked by a constraint or goes uphill, and a search has to find the narrow edge along which f still falls. Some searches converge there instead.
The constraints also have very different scales. g36 is in the units of y16, which is about 140,000 at the optimum, and g38 in those of y17, in the millions, while g2 is in the units of x2 and x3, below 300.
Representation
A Real genome of 5 genes, x1 to x5, within the report's bounds. genoxide's
problems::cec2006::G16 is the fitness: the value f(x) and the total constraint violation, the
sum of max(0, g(x)) over the 38 inequalities, 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, which matters with constraints of such different scales.
Algorithm
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. It uses genoxide's defaults: a population of 4 + ⌊3 ln 5⌋ = 8, a step size of 0.3 of each gene's range and a random start. With IPOP restarts (Auger and Hansen, 2005, IEEE CEC 2005: 1769-1776), a run that has converged starts again from a random point, with twice the population. A sample outside the bounds is drawn again, up to 100 times, and then clipped to them. Deb's rules rank the samples.
The run has the report's budget of 500,000 evaluations, and stops once 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.
Why CMA-ES: its covariance matrix can stretch the samples along a narrow direction, such as an edge between active constraints, and it needs few evaluations in 5 variables. With seeds 1 to 25, IPOP-CMA-ES met the target on all 25 runs, after a median of 12,216 evaluations (5,712 to 18,712). Without restarts, CMA-ES met it on 17: 6 of the other 8 converged at nearby corners, between 2.1e-5 and 3.3e-3 above f, and 3 of them missed the report's 1e-4; the last 2 converged at the best known point, but stopped 3.3e-8 and 7.5e-7 above f. SHADE (Tanabe and Fukunaga, 2013, IEEE CEC 2013: 71-78), genoxide's default differential evolution, also met the target on all 25, after a median of 37,500 evaluations (35,200 to 39,400), three times as many.
Output
The first line names the run. The second gives what stopped it, after how many evaluations and
restarts, 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 was first feasible, and when its error first
met the report's criterion of success. The fourth compares f(x) with f, to 6 significant digits,
and the fifth gives the solution, to 6 digits. The last names the active constraints, those with
|g| ≤ 1e-4, and counts the others. The threshold is larger than on other pages because of the
scales: the run stops with g4 and g36 about 4e-6 and 1e-5 from their boundaries. In Python, run evaluates
the problem in Rust, so both versions print the same.
The page shows each variable on its range, and each of the 38 constraints' state, as a grid. Its curve shows the error f − f* of the best feasible solution, and of the population's median, on a log scale. The best's curve begins at the first feasible solution, and the median's once half the population is feasible.
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 finds its first feasible solution after 208 evaluations. After about 2,000, it reaches a corner with x1 and x2 at their lower bounds, 704.415 and 68.6, x3 at 102.9 and x5 ≈ 25.01, where g2, g3 and g5 are active but g4 has slack: 3.4e-3 above f, a failure by the report's criterion. It converges there, and after 4,696 evaluations it restarts from a random point with a population of 16. The new run comes to the same corner after about 6,500 evaluations, and then moves on, x5 growing to about 37.58, where g4 is active too. It is within 1e-4 of f after 8,744 evaluations, and meets its target after 10,600, with one restart.
The solution is x = (705.175, 68.6000, 102.900, 282.325, 37.5841), the best known point to 6 digits, with the same five active constraints, g2, g3, g4, g5 and g36, and x2 at its lower bound. The other 33 constraints have slack.
Known optimum: −1.90515525853479 (best known)
Source: examples/cec2006_g16
Interactive run: tachsin.gr/projects/genoxide/examples/cec2006-g16
cargo run --release --example cec2006_g16
//! CEC 2006 g16: a nonlinear function of 5 variables, through a chain of 17 intermediate
//! quantities, under 38 inequality constraints, from the CEC 2006 special session on constrained
//! optimization (Liang et al., 2006). The best known value is −1.90515525853479.
//!
//! genoxide's `G16` gives the value of a solution and its constraint violation, which Deb's
//! feasibility rules compare: a feasible solution beats an infeasible one. CMA-ES with IPOP
//! restarts searches the 5 variables within the report's budget of 500,000 evaluations, and stops
//! once 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_g16
//! ```
mod trace;
use genoxide::prelude::*;
use genoxide::problems::Problem;
use genoxide::problems::cec2006::G16;
// 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: the scales differ, g36 bounds y16, which
// is about 140,000 there, and the run stops with the active constraints within about 1e-5
const ACTIVE: f64 = 1e-4;
fn main() -> Result<()> {
let problem = G16;
let optimum = problem.optimum().expect("known");
let f_star = optimum.value();
let cmaes = Cmaes::builder(problem.representation())
.restarts(cmaes::Restarts::Ipop)
.minimize()
.seed(1)
.build()?;
// 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);
// the restarts, each of which doubles the population
let (mut restarts, mut size) = (0, 0);
let outcome = Engine::new(cmaes, problem)
.stop_when(Stop::target(f_star + ERROR).or(Stop::evaluations(BUDGET)))
.on_generation(|snapshot| {
let progress = snapshot.progress();
let best = progress.best().filter(|best| best.is_feasible());
let error = best.and_then(Fitness::score).map(|value| value - f_star);
if error.is_some() {
feasible.get_or_insert(progress.evaluations());
}
if error.is_some_and(|error| error <= SUCCESS) {
success.get_or_insert(progress.evaluations());
}
let population = snapshot.population().len();
if size != 0 && population > size {
restarts += 1;
}
size = population;
trace.record(snapshot);
})
.run()?;
let best = outcome.best_fitness();
let value = best.score().expect("valid");
let x = outcome.best_genome();
println!("IPOP-CMA-ES with Deb's feasibility rules on g16, seed 1");
let (stop, error) = if outcome.stop_reason() == StopReason::Target {
("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 and {restarts} restarts: 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(", "));
// g1 to g38: the active ones, and how many of the others are satisfied or violated
let constraints = problem.constraints(x);
let g = constraints.inequalities();
let active: Vec<String> = (1..)
.zip(g)
.filter(|(_, g)| g.abs() <= ACTIVE)
.map(|(i, _)| format!("g{i}"))
.collect();
let satisfied = g.iter().filter(|&&g| g < -ACTIVE).count();
let violated = g.iter().filter(|&&g| g > ACTIVE).count();
println!(
"active (|g| <= 1e-4): {}; {satisfied} satisfied, {violated} violated",
active.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_g16/main.py
"""CEC 2006 g16: a nonlinear function of 5 variables, through a chain of 17 intermediate
quantities, under 38 inequality constraints, from the CEC 2006 special session on constrained
optimization (Liang et al., 2006). The best known value is −1.90515525853479.
genoxide's ``G16`` gives the value of a solution and its constraint violation, which Deb's
feasibility rules compare: a feasible solution beats an infeasible one. CMA-ES with IPOP restarts
searches the 5 variables within the report's budget of 500,000 evaluations, and stops once the
error f(x) − f* is at most 1e-8. The example prints the best solution and its constraints.
``run`` evaluates the problem in Rust.
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_g16/main.py
"""
import math
import genoxide as gx
import numpy as np
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: the scales differ, g36 bounds y16, which is
# about 140,000 there, and the run stops with the active constraints within about 1e-5
ACTIVE = 1e-4
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)
problem = gx.problems.cec2006.G16()
optimum = problem.optimum
f_star = optimum.value
cmaes = gx.Cmaes(problem.genome, restarts="ipop", objective=problem.objective, seed=1)
# 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; the restarts, each of which doubles the population
first = {"feasible": None, "success": None, "restarts": 0, "size": 0}
def on_generation(progress):
_, violations = problem.evaluate(progress.best_genome[np.newaxis])
if violations[0] == 0.0:
error = progress.best_fitness - f_star
if first["feasible"] is None:
first["feasible"] = progress.evaluations
if first["success"] is None and error <= SUCCESS:
first["success"] = progress.evaluations
size = len(progress.scores)
if first["size"] and size > first["size"]:
first["restarts"] += 1
first["size"] = size
trace.record(progress)
result = cmaes.run(
problem, target=f_star + ERROR, evaluations=BUDGET, on_generation=on_generation
)
value = result.best_fitness
print("IPOP-CMA-ES with Deb's feasibility rules on g16, seed 1")
if result.stop_reason == "target":
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 and {first['restarts']} restarts: "
f"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)))
# g1 to g38: the active ones, and how many of the others are satisfied or violated
g = problem.constraints(result.best_genome).tolist()
active = [f"g{i}" for i, gi in enumerate(g, 1) if abs(gi) <= ACTIVE]
satisfied = sum(gi < -ACTIVE for gi in g)
violated = sum(gi > ACTIVE for gi in g)
print(f"active (|g| <= 1e-4): {', '.join(active)}; {satisfied} satisfied, {violated} violated")
trace.write()
What it prints, from a seeded run:
IPOP-CMA-ES with Deb's feasibility rules on g16, seed 1
stopped by the target after 10600 evaluations and 1 restarts: f(x) - f* < 1e-8, feasible
first feasible after 208 evaluations, f(x) - f* <= 1e-4 after 8744
f(x) -1.90516, f* -1.90516 (best known)
x1 705.175, x2 68.6000, x3 102.900, x4 282.325, x5 37.5841
active (|g| <= 1e-4): g2, g3, g4, g5, g36; 33 satisfied, 0 violated