CTP4
The problem
Deb, Pratap and Meyarivan (2001) built the CTP problems to test how multi-objective algorithms handle constraints. CTP2 to CTP7 share one form, a generator whose six parameters θ, a, b, c, d and e shape a single constraint:
minimize f₁ = x₁
f₂ = g (1 − √(f₁/g)), g = 1 + x₂
subject to cos θ (f₂ − e) − sin θ f₁ ≥ a |sin(bπ (sin θ (f₂ − e) + cos θ f₁)^c)|^d
x₁, x₂ in [0, 1]
CTP4 is CTP3 with a wave 7.5 times higher: θ = −0.2π, a = 0.75, b = 10, c = 1, d = 0.5 and e = 1. The left side, u, measures the distance above the line f₂ = 1 − 0.7265 f₁, and the right side, a wave of height up to 0.75, touches the line where its sine is 0: at v = 0, 0.1, …, 1.2 along it. Where the wave is higher than the whole objective space above the line, only thin feasible tunnels are left, one down to each touch.
The optimal front is the same as CTP3's: 13 points, (cos θ v, 1 + sin θ v) for v = k/10, from (0, 1) to (0.9708, 0.2947), derived from the definition. At the points themselves the constraint is exactly 0; in floating point, sin(kπ) is about 1e-15 and its square root 3e-8, so the best feasible solutions sit next to them.
The definitions are the paper's: eq. 5 (p. 290) and CTP4's parameters (p. 291). Its preprint, the authors' KanGAL report 200005 (October 2000, p. 8), and Deb's 2001 book (Multi-Objective Optimization Using Evolutionary Algorithms, Wiley, eq. 8.46 on p. 354 and the parameters on p. 356) give the same. All three leave g, the number of variables and their bounds open, and print f₂ as g (1 − f₁/g); their figures draw the unconstrained front as the curve 1 − √f₁, and the authors' NSGA-II code (version 1.1.6, KanGAL) computes g (1 − √(f₁/g)) with g = 1 + x₂, the g the book names on p. 360, and two variables in [0, 1], which genoxide follows.
What makes it hard
The paper puts it this way: "an algorithm now has to travel through a long narrow feasible tunnel in search of the lone Pareto-optimal solution at the end of tunnel". At a height u above the line, a tunnel reaches only about (u/a)²/(bπ) to each side: 6e-6 at u = 0.01, 60 times less than CTP3's. Only 3% of random genomes are feasible, and the tunnels near f₁ = 1 don't connect to the open region above: the population has to land in each of them. Every step closer to a tip needs an offspring inside a narrower sliver than the last, and the steps of crossover and mutation, which change x₁ and x₂ each on its own, rarely follow a tunnel's slant.
The paper found that neither NSGA-II nor Ray et al.'s algorithm got near the 13 points.
Representation
A Real genome of 2 genes in [0, 1]: x₁ and x₂. The problem is genoxide's Ctp4, whose fitness
is the two objectives and the constraint violation.
Solutions compare by constrained dominance: a feasible solution beats an infeasible one, of two infeasible ones the smaller violation wins, and of two feasible ones Pareto dominance, or for MOEA/D the subproblem's value, decides.
Algorithm
Two runs with the same operators: simulated binary crossover with η = 20, at a rate of 0.9, and polynomial mutation with η = 20, at a rate of 1/n per gene for n genes, 0.5.
- NSGA-II with the settings of the paper's experiments, a population of 100 for 500 generations.
- MOEA/D (Zhang and Li, 2007) with 100 subproblems, weight vectors evenly spread on the simplex, Tchebycheff decomposition and genoxide's default neighborhoods of 20, for 70,000 generations: 7 million evaluations, 3.5 seconds on a desktop. Each subproblem keeps its best solution for its direction and lets a child replace a neighbor's only when better, so the solutions in a tunnel creep down it, one small improvement at a time. NSGA-II gets there too with enough generations (an IGD+ of 0.0077 to 0.0089 over seeds 1 to 4 after 100,000), but MOEA/D is faster.
Output
For each run, the first line gives the size of the final front and how many of its solutions are feasible.
The second counts the optimal points that the run reaches: that have a solution within 0.02, in scaled objectives.
The third gives the front's IGD+ (Ishibuchi et al., 2015, EMO 2015, LNCS 9019: 110-125) to 2,000
points of the optimal front (the 13 points, each repeated), and its hypervolume, the area it
dominates up to the reference point (1.1, 1.1), as a share of the optimal front's. Both use
objectives scaled to [0, 1] on the front, by its ideal point (0, 0.2947) and nadir point
(0.9708, 1). IGD+ averages, over the points of the optimal front, the distance to the nearest
solution, counting only the objectives in which the solution is worse. The 13 points'
hypervolume is 0.6683. In Python, run evaluates the problem in Rust, so both versions print the
same.
The run's trace.json, of the MOEA/D run, also has the problem's feasible region, which the page
shades; the tunnels near their tips are narrower than its cells.
The project page plays this run back.
Good results
A good front has a solution next to each of the 13 points, and an IGD+ under 0.01. No finite set of feasible solutions reaches the points' hypervolume, since the points themselves are the limit.
NSGA-II with the paper's budget ends in the tunnels, far from their tips: 22 solutions, none within 0.02 of a point, an IGD+ of 0.097 and 79.18% of the hypervolume. Over seeds 1 to 20 its IGD+ is 0.070 to 0.139, and only one run comes within 0.02 of a single point: the paper's finding again.
MOEA/D reaches all 13: 16 solutions, an IGD+ of 0.0066 and 98.57% of the hypervolume. Over seeds 1 to 20, every run reaches all 13 points, with an IGD+ from 0.0065 to 0.0088.
Known optimum: 13 points on the line f₂ = 1 − tan(0.2π) f₁, from (0, 1) to (0.9708, 0.2947); hypervolume 0.6683 in objectives scaled by the ideal and nadir points (reference point (1.1, 1.1))
Source: examples/ctp4
Interactive run: tachsin.gr/projects/genoxide/examples/ctp4
cargo run --release --example ctp4
//! CTP4: minimize two objectives over two variables subject to one constraint, with
//! a front of 13 separate points, each at the end of a long, narrow feasible tunnel.
//!
//! NSGA-II with the settings of the paper's experiments, a population of 100 for 500
//! generations, and MOEA/D with the same operators for 70,000. Prints, for each, how many
//! solutions of the final front are feasible, how many of the optimal points they reach, their
//! IGD+ to the front and their hypervolume.
//!
//! With `GENOXIDE_TRACE=<file>`, it also writes a trace of the MOEA/D run for the plot on
//! the example's page, with `trace.rs`.
//!
//! ```text
//! cargo run --release --example ctp4
//! ```
mod trace;
use genoxide::Objective::Minimize;
use genoxide::multi::indicator::{hypervolume, igd_plus};
use genoxide::multi::problems::{Ctp4, MultiProblem};
use genoxide::prelude::*;
// how close a solution must come to a piece of the optimal front, in scaled objectives, to reach it
const REACH: f64 = 0.02;
fn main() -> Result<()> {
let problem = Ctp4;
// the settings of the paper's experiments: a population of 100 for 500 generations, SBX
// and polynomial mutation with η = 20, crossover at 0.9 and mutation at 1/n per gene
let nsga2 = Nsga2::builder(problem.representation(), [Minimize; 2])
.population_size(100)
.crossover(SimulatedBinaryCrossover::new(20.0)?)
.mutate(PolynomialMutation::per_gene(0.5, 20.0)?)
.seed(1)
.build()?;
let outcome = MultiEngine::new(nsga2, problem)
.stop_when(Stop::generations(500))
.run()?;
let (front, size) = feasible(outcome.front());
report("NSGA-II", 500, &front, size);
// MOEA/D with 100 subproblems, the same operators and 70,000 generations: 7 million
// evaluations
let moead = Moead::builder(
problem.representation(),
[Minimize; 2],
multi::das_dennis::<2>(99),
)
.crossover(SimulatedBinaryCrossover::new(20.0)?)
.mutate(PolynomialMutation::per_gene(0.5, 20.0)?)
.seed(1)
.build()?;
// with GENOXIDE_TRACE=<file>, a trace of the MOEA/D run for the plot on the example's page
let mut trace = trace::Trace::from_env();
let outcome = MultiEngine::new(moead, problem)
.stop_when(Stop::generations(70_000))
.on_generation(|snapshot| trace.record(snapshot))
.run()?;
let (front, size) = feasible(outcome.front());
report("MOEA/D", 70_000, &front, size);
trace.write();
Ok(())
}
// the objectives scaled to [0, 1] on the optimal front, by its ideal and nadir points
fn scaled(points: &[[f64; 2]]) -> Vec<[f64; 2]> {
let (ideal, nadir) = (
Ctp4.ideal_point().expect("known"),
Ctp4.nadir_point().expect("known"),
);
let scale = |p: &[f64; 2]| [0, 1].map(|j| (p[j] - ideal[j]) / (nadir[j] - ideal[j]));
points.iter().map(scale).collect()
}
// the pieces of a front sorted by f₁, split where neighbors are more than 0.01 apart
fn pieces(front: &[[f64; 2]]) -> Vec<Vec<[f64; 2]>> {
let mut pieces: Vec<Vec<[f64; 2]>> = Vec::new();
for (i, point) in front.iter().enumerate() {
let gap = i == 0 || {
let last = front[i - 1];
((point[0] - last[0]).powi(2) + (point[1] - last[1]).powi(2)).sqrt() > 0.01
};
if gap {
pieces.push(Vec::new());
}
pieces.last_mut().expect("a piece").push(*point);
}
pieces
}
// the feasible solutions of a run's front: how many of them, how many pieces of the optimal
// front they reach, their IGD+ to it and their hypervolume, as a share of the whole front's
fn report(name: &str, generations: u64, front: &[[f64; 2]], size: usize) {
let feasible = if front.len() == size {
"all feasible".to_string()
} else {
format!("{} feasible", front.len())
};
println!("{name}, {generations} generations: {size} solutions on the front, {feasible}");
let optimal = Ctp4.optimal_front(2000).expect("known");
let found = scaled(front);
let pieces = pieces(&optimal);
let reached = pieces
.iter()
.filter(|piece| {
scaled(piece).iter().any(|p| {
found
.iter()
.any(|f| ((f[0] - p[0]).powi(2) + (f[1] - p[1]).powi(2)).sqrt() <= REACH)
})
})
.count();
println!(" optimal points reached: {reached} of {}", pieces.len());
let distance = igd_plus(&found, &scaled(&optimal), &[Minimize; 2]);
let volume = hypervolume(&found, &[1.1, 1.1], &[Minimize; 2]);
let whole = whole_front_hypervolume();
println!(
" IGD+ {distance:.5}, hypervolume {volume:.4}, {:.2}% of the whole front's {whole:.4}",
100.0 * volume / whole
);
}
// the hypervolume of the whole optimal front, from 100,000 of its points, in scaled objectives with
// the reference point (1.1, 1.1)
pub fn whole_front_hypervolume() -> f64 {
let front = scaled(&Ctp4.optimal_front(100_000).expect("known"));
hypervolume(&front, &[1.1, 1.1], &[Minimize; 2])
}
// the objective values of the feasible solutions of a front, and the front's size
fn feasible<G: Genome>(front: &[Individual<G, multi::Scores<2>>]) -> (Vec<[f64; 2]>, usize) {
let values = front
.iter()
.filter_map(|x| {
x.fitness()
.filter(|s| s.is_feasible())
.and_then(|s| s.values())
})
.collect();
(values, front.len())
}
python examples/ctp4/main.py
"""CTP4: minimize two objectives over two variables subject to one constraint, with
a front of 13 separate points, each at the end of a long, narrow feasible tunnel.
NSGA-II with the settings of the paper's experiments, a population of 100 for 500
generations, and MOEA/D with the same operators for 70,000. Prints, for each, how many
solutions of the final front are feasible, how many of the optimal points they reach, their
IGD+ to the front and their hypervolume.
With ``GENOXIDE_TRACE=<file>``, it also writes a trace of the MOEA/D run for the plot on
the example's page, with trace.py.
python examples/ctp4/main.py
"""
import numpy as np
import genoxide as gx
from trace import Trace, scaled, whole_front_hypervolume
problem = gx.problems.Ctp4()
# how close a solution must come to a piece of the optimal front, in scaled objectives, to reach it
REACH = 0.02
def pieces(front):
"""The pieces of a front sorted by f₁, split where neighbors are more than 0.01 apart."""
gaps = np.flatnonzero(np.hypot(*np.diff(front, axis=0).T) > 0.01) + 1
return np.split(front, gaps)
def report(name, generations, result):
"""The feasible solutions of a run's front: how many of them, how many pieces of the optimal
front they reach, their IGD+ to it and their hypervolume, as a share of the whole front's."""
front = result.front_objectives[result.front_violations == 0]
size = len(result.front_objectives)
feasible = "all feasible" if len(front) == size else f"{len(front)} feasible"
print(f"{name}, {generations} generations: {size} solutions on the front, {feasible}")
optimal = problem.optimal_front(2000)
found = scaled(front)
parts = pieces(optimal)
reached = sum(
any(np.hypot(*(found - p).T).min() <= REACH for p in scaled(piece)) for piece in parts
)
print(f" optimal points reached: {reached} of {len(parts)}")
distance = gx.indicators.igd_plus(found, scaled(optimal))
volume = gx.indicators.hypervolume(found, [1.1, 1.1])
whole = whole_front_hypervolume()
print(
f" IGD+ {distance:.5f}, hypervolume {volume:.4f}, "
f"{100 * volume / whole:.2f}% of the whole front's {whole:.4f}"
)
def nsga2():
"""NSGA-II with the settings of the paper's experiments: a population of 100, SBX and
polynomial mutation with η = 20, crossover at 0.9 and mutation at 1/n per gene."""
return gx.Nsga2(
problem.genome,
objectives=problem.objectives,
population_size=100,
crossover=gx.SimulatedBinaryCrossover(20),
mutation=gx.PolynomialMutation(20, rate=0.5),
seed=1,
)
report("NSGA-II", 500, nsga2().run(problem, generations=500))
# MOEA/D with 100 subproblems, the same operators and 70,000 generations: 7 million evaluations
moead = gx.Moead(
problem.genome,
objectives=problem.objectives,
weights=gx.das_dennis(2, 99),
crossover=gx.SimulatedBinaryCrossover(20),
mutation=gx.PolynomialMutation(20, rate=0.5),
seed=1,
)
# with GENOXIDE_TRACE=<file>, a trace of the MOEA/D run for the plot on the example's page
trace = Trace(problem)
result = moead.run(problem, generations=70_000, on_generation=trace.on_generation)
report("MOEA/D", 70_000, result)
trace.write()
What it prints, from a seeded run:
NSGA-II, 500 generations: 22 solutions on the front, all feasible
optimal points reached: 0 of 13
IGD+ 0.09733, hypervolume 0.5292, 79.18% of the whole front's 0.6683
MOEA/D, 70000 generations: 16 solutions on the front, all feasible
optimal points reached: 13 of 13
IGD+ 0.00662, hypervolume 0.6587, 98.57% of the whole front's 0.6683