@@ 1,11 1,12 @@
-use std::convert::Infallible;
+use std::{convert::Infallible, env, fs, io::Write, rc::Rc};
use eoa_lib::{
- bounded::BoundedOVector, comparison::MinimizingOperator, constraints::{stochastic_ranking_evolution_algorithm, ConstraintFunction, LowerThanConstraintFunction}, crossover::{ArithmeticCrossover, BoundedCrossover, BoundedCrossoverStrategy}, evolution::EvolutionResult, fitness::FitnessFunction, initializer::{Initializer, RandomInitializer}, pairing::AdjacentPairing, perturbation::{BoundedPerturbation, BoundedPerturbationStrategy, MutationPerturbation, RandomDistributionPerturbation}, population::Population, replacement::{BestReplacement, GenerationalReplacement}, selection::{BestSelection, TournamentSelection}
+ bounded::BoundedOVector, comparison::MinimizingOperator, constraints::{stochastic_ranking_evolution_algorithm, ConstrainedEvalFitness, ConstrainedFitnessFunction, ConstraintFunction, LowerThanConstraintFunction}, crossover::{ArithmeticCrossover, BoundedCrossover, BoundedCrossoverStrategy}, evolution::{EvolutionCandidate, EvolutionResult, EvolutionStats}, fitness::FitnessFunction, initializer::{Initializer, RandomInitializer}, multi_objective_evolution::nsga_2, pairing::AdjacentPairing, perturbation::{BoundedPerturbation, BoundedPerturbationStrategy, MutationPerturbation, RandomDistributionPerturbation}, population::{EvaluatedChromosome, EvaluatedPopulation, Population}, replacement::{BestReplacement, GenerationalReplacement}, selection::{BestSelection, TournamentSelection}
};
-use nalgebra::SVector;
+use nalgebra::{ArrayStorage, Const, Matrix, SVector};
use rand::RngCore;
use rand_distr::Normal;
+use chrono::prelude::*;
pub struct ArbitraryFitness<const SIZE: usize> {
fun: Box<dyn Fn(SVector<f64, SIZE>) -> f64>
@@ 17,6 18,12 @@ impl<const SIZE: usize> ArbitraryFitness<SIZE> {
fun
}
}
+
+ pub fn zero() -> Self {
+ Self {
+ fun: Box::new(|_| 0.0)
+ }
+ }
}
impl<const SIZE: usize> FitnessFunction for ArbitraryFitness<SIZE> {
@@ 30,36 37,72 @@ impl<const SIZE: usize> FitnessFunction for ArbitraryFitness<SIZE> {
}
/// A constrained optimization problem with clear field definitions
-pub struct ConstrainedProblem<const CONSTRAINTS: usize> {
- pub objective: ArbitraryFitness<2>,
- pub constraints: [LowerThanConstraintFunction<SVector<f64, 2>, f64>; CONSTRAINTS],
- pub bounds: (SVector<f64, 2>, SVector<f64, 2>), // (min, max)
+pub struct ConstrainedProblem<const DIM: usize, const CONSTRAINTS: usize> {
+ pub name: String,
+ pub objective: ArbitraryFitness<DIM>,
+ pub constraints: [LowerThanConstraintFunction<SVector<f64, DIM>, f64>; CONSTRAINTS],
+ pub bounds: (SVector<f64, DIM>, SVector<f64, DIM>), // (min, max)
pub optimal_value: f64,
+ pub instantiate_fn: Option<Rc<dyn Fn() -> ConstrainedProblem<DIM, CONSTRAINTS>>>
}
-fn problem_g06() -> ConstrainedProblem<2> {
- ConstrainedProblem {
- objective: ArbitraryFitness::new(
- Box::new(|vec| (vec[0] - 10.0).powi(3) + (vec[1] - 20.0).powi(3))
- ),
- constraints: [
- LowerThanConstraintFunction::new(
- Box::new(|vec| -(vec[0] - 5.0).powi(2) - (vec[1] - 5.0).powi(2) + 100.0)
+impl<const DIM: usize, const CONSTRAINTS: usize> ConstrainedProblem<DIM, CONSTRAINTS> {
+ pub fn new(instantiate: Rc<dyn Fn() -> ConstrainedProblem<DIM, CONSTRAINTS>>) -> Self {
+ let mut problem = instantiate();
+ problem.instantiate_fn = Some(instantiate);
+ problem
+ }
+
+ pub fn clone(&self) -> Self {
+ Self::new(self.instantiate_fn.clone().unwrap())
+ }
+}
+
+/// Configuration for stochastic ranking method
+pub struct StochasticRankingConfig {
+ pub population_size: usize,
+ pub parents_count: usize,
+ pub iterations: usize,
+ pub n_param: usize, // N parameter for stochastic ranking
+ pub p_param: f64, // p parameter for stochastic ranking
+ pub mutation_std_dev: f64,
+}
+
+/// Configuration for NSGA-II method
+pub struct NsgaConfig {
+ pub population_size: usize,
+ pub parents_count: usize,
+ pub iterations: usize,
+ pub mutation_std_dev: f64,
+}
+
+fn problem_g06() -> ConstrainedProblem<2, 2> {
+ ConstrainedProblem::new(Rc::new(||
+ ConstrainedProblem {
+ name: "g06".to_string(),
+ objective: ArbitraryFitness::new(
+ Box::new(|vec| (vec[0] - 10.0).powi(3) + (vec[1] - 20.0).powi(3))
),
- LowerThanConstraintFunction::new(
- Box::new(|vec| (vec[0] - 6.0).powi(2) + (vec[1] - 5.0).powi(2) - 82.81)
+ constraints: [
+ LowerThanConstraintFunction::new(
+ Box::new(|vec| -(vec[0] - 5.0).powi(2) - (vec[1] - 5.0).powi(2) + 100.0)
+ ),
+ LowerThanConstraintFunction::new(
+ Box::new(|vec| (vec[0] - 6.0).powi(2) + (vec[1] - 5.0).powi(2) - 82.81)
+ ),
+ ],
+ bounds: (
+ SVector::<f64, 2>::new(0.0, 0.0), // min bounds
+ SVector::<f64, 2>::new(50.0, 50.0) // max bounds
),
- ],
- bounds: (
- SVector::<f64, 2>::new(0.0, 0.0), // min bounds
- SVector::<f64, 2>::new(100.0, 100.0) // max bounds
- ),
- optimal_value: -6961.8137558015,
- }
+ optimal_value: -6961.8137558015,
+ instantiate_fn: None
+ }))
}
-fn problem_g08() -> ConstrainedProblem<2> {
- ConstrainedProblem {
+fn problem_g08() -> ConstrainedProblem<2, 2> {
+ ConstrainedProblem::new(Rc::new(|| ConstrainedProblem {
+ name: "g08".to_string(),
objective: ArbitraryFitness::new(
Box::new(|vec| {
let num = (2.0 * std::f64::consts::PI * vec[0]).sin().powi(3)
@@ 89,11 132,13 @@ fn problem_g08() -> ConstrainedProblem<2> {
SVector::<f64, 2>::new(10.0, 10.0) // max bounds
),
optimal_value: -0.0958250414180359,
- }
+ instantiate_fn: None
+ }))
}
-pub fn problem_g11(eps: f64) -> ConstrainedProblem<2> {
- ConstrainedProblem {
+pub fn problem_g11(eps: f64) -> ConstrainedProblem<2, 1> {
+ let problem = ConstrainedProblem::new(Rc::new(move || ConstrainedProblem {
+ name: "g11".to_string(),
objective: ArbitraryFitness::new(
Box::new(|vec| {
// Minimize f(x) = x1^2 + (x2 - 1)^2
@@ 102,31 147,28 @@ pub fn problem_g11(eps: f64) -> ConstrainedProblem<2> {
),
constraints: [
// Equality h(x) = x2 - x1^2 = 0
- // Transformed 1: h - eps >= 0 => -(h - eps) <= 0
+ // |h| - eps >= 0.0
LowerThanConstraintFunction::new(
Box::new(move |vec| {
let h = vec[1] - vec[0].powi(2);
- h - eps
- })
- ),
- // Transformed 2: eps - h >= 0 => -(eps - h) <= 0
- LowerThanConstraintFunction::new(
- Box::new(move |vec| {
- let h = vec[1] - vec[0].powi(2);
- eps - h
+ h.abs() - eps
})
),
],
bounds: (
- SVector::<f64, 2>::new(-1.0, -1.0), // min bounds
- SVector::<f64, 2>::new(1.0, 1.0) // max bounds
+ SVector::<f64, 2>::new(-50.0, -50.0), // min bounds
+ SVector::<f64, 2>::new(50.0, 50.0) // max bounds
),
optimal_value: 0.7499, // Best known optimum
- }
+ instantiate_fn: None
+ }));
+
+ problem
}
-pub fn problem_g24() -> ConstrainedProblem<2> {
- ConstrainedProblem {
+pub fn problem_g24() -> ConstrainedProblem<2, 2> {
+ ConstrainedProblem::new(Rc::new(|| ConstrainedProblem {
+ name: "g24".to_string(),
objective: ArbitraryFitness::new(
Box::new(|vec| {
// Minimize f(x) = -x1 - x2
@@ 152,27 194,28 @@ pub fn problem_g24() -> ConstrainedProblem<2> {
SVector::<f64, 2>::new(3.0, 4.0) // max bounds
),
optimal_value: -5.50801327159536, // Best known optimum
- }
+ instantiate_fn: None
+ }))
}
/// Solve a constrained optimization problem using stochastic ranking
///
/// Returns the evolution result with feasible fractions for each iteration
-pub fn solve_with_stochastic_ranking<const CONSTRAINTS: usize>(
- mut problem: ConstrainedProblem<CONSTRAINTS>,
+pub fn solve_with_stochastic_ranking<const DIM: usize, const CONSTRAINTS: usize>(
+ mut problem: ConstrainedProblem<DIM, CONSTRAINTS>,
population_size: usize,
parents_count: usize,
iterations: usize,
stochastic_params: (usize, f64), // (N, p) for stochastic ranking
mutation_std_dev: f64,
rng: &mut dyn RngCore,
-) -> Result<(ConstrainedProblem<CONSTRAINTS>, EvolutionResult<SVector<f64, 2>, f64>, Vec<f64>), Box<dyn std::error::Error>> {
+) -> Result<(EvolutionResult<SVector<f64, DIM>, f64>, Vec<f64>), Box<dyn std::error::Error>> {
// Create initial population
let initializer = RandomInitializer::new(
Box::new(BoundedOVector::new(problem.bounds.0, problem.bounds.1))
);
let initial_population = Population::from_vec(
- initializer.initialize(nalgebra::Const::<2>, population_size, rng)
+ initializer.initialize(nalgebra::Const::<DIM>, population_size, rng)
);
// Setup components as specified
@@ 180,24 223,27 @@ pub fn solve_with_stochastic_ranking<const CONSTRAINTS: usize>(
let mut selection = TournamentSelection::new(5, 0.95);
let mut replacement = GenerationalReplacement;
let mut pairing = AdjacentPairing::new();
- let mut crossover = BoundedCrossover::<nalgebra::Const<2>, 2, _>::new(
- ArithmeticCrossover::new(),
- problem.bounds.0,
- problem.bounds.1,
- BoundedCrossoverStrategy::Trim
- );
+ let mut crossover = ArithmeticCrossover::new();
+ // let mut crossover = BoundedCrossover::<nalgebra::Const<2>, 2, _>::new(
+ // ArithmeticCrossover::new(),
+ // problem.bounds.0,
+ // problem.bounds.1,
+ // BoundedCrossoverStrategy::Retry(5)
+ // );
// Setup bounded random distribution perturbation with Normal distribution
- let normal_perturbation = RandomDistributionPerturbation::<2, Normal<f64>>::normal(mutation_std_dev)?;
- let perturbation = BoundedPerturbation::new(
- normal_perturbation,
- problem.bounds.0,
- problem.bounds.1,
- BoundedPerturbationStrategy::Trim
- );
+ let normal_perturbation = RandomDistributionPerturbation::<DIM, Normal<f64>>::normal(mutation_std_dev)?;
+ let mut perturbation = normal_perturbation;
+ // let perturbation = BoundedPerturbation::new(
+ // normal_perturbation,
+ // problem.bounds.0,
+ // problem.bounds.1,
+ // BoundedPerturbationStrategy::Retry(5)
+ // );
let mut mutation = MutationPerturbation::new(Box::new(perturbation), 0.1);
- // Constraint weights (initially set to 1.0)
+ // The weight is so large mainly because of the g11 that has very small values.
+ // Somehow the higher weights do seem to help, even though I am unsure why exactly.
let constraint_weights = [1.0; CONSTRAINTS];
let (N, p) = stochastic_params;
@@ 222,59 268,676 @@ pub fn solve_with_stochastic_ranking<const CONSTRAINTS: usize>(
&mut replacement,
&better_than,
iterations,
+ rng)
+}
+
+/// Helper function to check if a chromosome is feasible
+fn check_feasibility<const DIM: usize, const CONSTRAINTS: usize>(
+ chromosome: &SVector<f64, DIM>,
+ constraints: &[LowerThanConstraintFunction<SVector<f64, DIM>, f64>; CONSTRAINTS],
+) -> bool {
+ constraints.iter().all(|constraint| {
+ constraint.is_feasible(chromosome).unwrap_or(false)
+ })
+}
+
+/// Individual constraint fitness function - wraps a single constraint to act as a fitness function
+pub struct SingleConstraintFitness<const DIM: usize, TConstraint: ConstraintFunction<Chromosome = SVector<f64, DIM>, Out = f64>> {
+ constraint: TConstraint,
+ capped: bool
+}
+
+impl<const DIM: usize, TConstraint: ConstraintFunction<Chromosome = SVector<f64, DIM>, Out = f64>> SingleConstraintFitness<DIM, TConstraint> {
+ pub fn new(constraint: TConstraint, capped: bool) -> Self {
+ Self { constraint, capped }
+ }
+}
+
+impl<const DIM: usize, TConstraint: ConstraintFunction<Chromosome = SVector<f64, DIM>, Out = f64>> FitnessFunction for SingleConstraintFitness<DIM, TConstraint> {
+ type In = SVector<f64, DIM>;
+ type Out = f64;
+ type Err = Infallible;
+
+ fn fit(&self, inp: &Self::In) -> Result<Self::Out, Self::Err> {
+ Ok(if self.constraint.is_feasible(inp).unwrap() && self.capped {
+ 0.0
+ } else {
+ self.constraint.evaluate(inp).unwrap()
+ })
+ }
+}
+
+/// Solve a constrained optimization problem using NSGA-II
+pub fn solve_with_nsga_ii<const DIM: usize, const CONSTRAINTS: usize>(
+ problem: ConstrainedProblem<DIM, CONSTRAINTS>,
+ population_size: usize,
+ parents_count: usize,
+ iterations: usize,
+ mutation_std_dev: f64,
+ rng: &mut dyn RngCore,
+) -> Result<(EvolutionResult<SVector<f64, DIM>, [f64; 2]>, Vec<f64>), Box<dyn std::error::Error>> {
+ // Create initial population
+ let initializer = RandomInitializer::new(
+ Box::new(BoundedOVector::new(problem.bounds.0, problem.bounds.1))
+ );
+ let initial_population = Population::from_vec(
+ initializer.initialize(nalgebra::Const::<DIM>, population_size, rng)
+ );
+
+ // Setup components
+ let mut pairing = AdjacentPairing::new();
+ let mut crossover = ArithmeticCrossover::new();
+
+ // Setup bounded random distribution perturbation with Normal distribution
+ let normal_perturbation = RandomDistributionPerturbation::<DIM, Normal<f64>>::normal(mutation_std_dev)?;
+ let mut mutation = MutationPerturbation::new(Box::new(normal_perturbation), 0.1);
+
+ let better_than = MinimizingOperator::new();
+
+ // Create objectives: both using ConstrainedFitnessFunction with different fitness components
+ let constraint_weights = [1.0e9; CONSTRAINTS];
+ let constraint_refs = problem.constraints.iter().collect::<Vec<_>>().try_into()
+ .map_err(|_| "Failed to convert constraint references")?;
+
+ let zero_fitness = ArbitraryFitness::zero();
+
+ // Second objective: constraint violation only (zero fitness + constraints)
+ let constrained_fitness_obj2 = ConstrainedFitnessFunction {
+ fitness: &zero_fitness,
+ constraints: constraint_refs,
+ constraint_weights,
+ capped: true
+ };
+
+ let objectives: [Box<dyn FitnessFunction<In = SVector<f64, DIM>, Out = f64, Err = Infallible>>; 2] = [
+ Box::new(problem.objective),
+ Box::new(constrained_fitness_obj2),
+ ];
+
+ let mut feasible_fractions = Vec::with_capacity(iterations);
+
+ let result = nsga_2(
+ initial_population,
+ parents_count,
+ objectives,
+ &mut pairing,
+ &mut crossover,
+ &mut mutation,
+ &better_than,
+ iterations,
rng,
- ).map(|(x, y)| (problem, x, y))
+ |_iteration: usize, _stats: &EvolutionStats<SVector<f64, DIM>, _>, population: &EvaluatedPopulation<SVector<f64, DIM>, _>| {
+ // Calculate feasible fraction based on second objective being close to zero
+ let feasible_count = population.population
+ .iter()
+ .filter(|individual| {
+ individual.evaluation.evaluations[1] <= 0.0
+ })
+ .count();
+
+ // let min_constraint_violation = population.population
+ // .iter()
+ // .map(|individual| {
+ // individual.evaluation.evaluations[1]
+ // })
+ // .min_by(|a, b| a.total_cmp(b)).unwrap();
+ // println!("{}", min_constraint_violation);
+
+ let feasible_fraction = feasible_count as f64 / population.population.len() as f64;
+ feasible_fractions.push(feasible_fraction);
+ },
+ |_, evaluation, best_candidate| {
+ // Do not save infeasible solutions!
+ if evaluation.evaluations[1] > 0.0 {
+ return false;
+ }
+
+ if best_candidate.is_none() {
+ return true;
+ }
+
+ evaluation.evaluations[0] < best_candidate.as_ref().unwrap().evaluated_chromosome.evaluation.evaluations[0]
+ }
+ )?;
+
+ Ok((result, feasible_fractions))
}
-fn main() {
- let mut rng = rand::rng();
+/// Solve a constrained optimization problem using NSGA-II with individual constraint objectives
+/// For simplicity, this function only works with 2-constraint problems
+pub fn solve_with_nsga_multi<const DIM: usize, const CONSTRAINTS: usize, const CONSTRS_PLUS_ONE: usize>(
+ problem: ConstrainedProblem<DIM, CONSTRAINTS>,
+ population_size: usize,
+ parents_count: usize,
+ iterations: usize,
+ mutation_std_dev: f64,
+ rng: &mut dyn RngCore,
+ capped: bool
+) -> Result<(EvolutionResult<SVector<f64, DIM>, [f64; CONSTRS_PLUS_ONE]>, Vec<f64>), Box<dyn std::error::Error>> {
+ // Unfortunately Rustc doesn't support addition in generics...
+ assert_eq!(CONSTRAINTS + 1, CONSTRS_PLUS_ONE);
- // Example usage with problem G06
- let problem = problem_g11(0.01);
+ // Clone the problem to get bounds info first
+ let bounds = (problem.bounds.0, problem.bounds.1);
- let population_size = 500;
- let result = solve_with_stochastic_ranking(
- problem,
- population_size, // population_size
- population_size / 2, // parents_count
- 1000, // iterations
- (population_size * 2, 0.2), // (N, p) for stochastic ranking
- 0.001, // mutation_std_dev
- &mut rng,
+ // Create initial population
+ let initializer = RandomInitializer::new(
+ Box::new(BoundedOVector::new(bounds.0, bounds.1))
);
+ let initial_population = Population::from_vec(
+ initializer.initialize(nalgebra::Const::<DIM>, population_size, rng)
+ );
+
+ // Setup components
+ let mut pairing = AdjacentPairing::new();
+ let mut crossover = ArithmeticCrossover::new();
+ let normal_perturbation = RandomDistributionPerturbation::<DIM, Normal<f64>>::normal(mutation_std_dev)?;
+ let mut mutation = MutationPerturbation::new(Box::new(normal_perturbation), 0.1);
+ let better_than = MinimizingOperator::new();
+
+ // Create objectives: fitness + individual constraints using cloned problem
+ let mut objective = Some(problem.objective);
+ let mut constraints = problem.constraints.into_iter();
+
+ let objectives: [Box<dyn FitnessFunction<In = SVector<f64, DIM>, Out = f64, Err = Infallible>>; CONSTRS_PLUS_ONE] = std::array::from_fn(move |i| {
+ let val: Box<dyn FitnessFunction<In = SVector<f64, DIM>, Out = f64, Err = Infallible>> = if i == 0 {
+ let obj = objective.take().expect("Taken already!");
+ Box::new(obj)
+ } else {
+ Box::new(SingleConstraintFitness::new(constraints.next().unwrap(), capped))
+ };
+
+ val
+ });
- match result {
- Ok((problem, evolution_result, feasible_fractions)) => {
- if let Some(best_candidate) = evolution_result.best_candidate {
- println!("{:?}", feasible_fractions);
- println!("Best solution found:");
- println!(" Chromosome: {:?}", best_candidate.chromosome);
- println!(" Fitness: {}", best_candidate.evaluation);
- println!(" Iterations: {}", evolution_result.iterations);
- println!(" Evaluations: {}", evolution_result.evaluations);
- println!(" Final feasible fraction: {:.2}%", feasible_fractions.last().unwrap_or(&0.0) * 100.0);
-
- // Sanity check: verify the best candidate is feasible
- let best_chromosome = &best_candidate.chromosome;
-
- println!("\nFeasibility check for best solution:");
- for (i, constraint) in problem.constraints.iter().enumerate() {
- match constraint.evaluate(best_chromosome) {
- Ok(value) => {
- let is_feasible = value <= 0.0;
- println!(" Constraint {}: {} ({})", i+1, value, if is_feasible { "FEASIBLE" } else { "INFEASIBLE" });
- }
- Err(e) => {
- println!(" Constraint {}: Error evaluating - {}", i+1, e);
- }
+ let mut feasible_fractions = Vec::with_capacity(iterations);
+
+ let result = nsga_2(
+ initial_population,
+ parents_count,
+ objectives,
+ &mut pairing,
+ &mut crossover,
+ &mut mutation,
+ &better_than,
+ iterations,
+ rng,
+ |_iteration: usize, _stats: &EvolutionStats<SVector<f64, DIM>, _>, population: &EvaluatedPopulation<SVector<f64, DIM>, _>| {
+ // Calculate feasible fraction - all constraints (objectives 1 and 2) must be <= 0
+ let feasible_count: f64 =
+ population.population.iter().filter(
+ |individual| {
+ individual.evaluation.evaluations
+ .iter()
+ .skip(1)
+ .all(|&eval| eval <= 0.0)
}
+ ).count() as f64;
+
+ let feasible_fraction = feasible_count as f64 / population.population.len() as f64;
+ feasible_fractions.push(feasible_fraction);
+ },
+ |_, evaluation, best_candidate| {
+ // Only save feasible solutions (all constraints satisfied)
+ if (1..3).any(|i| evaluation.evaluations[i] > 0.0) {
+ return false;
+ }
+
+ if best_candidate.is_none() {
+ return true;
+ }
+
+ // Compare based on fitness (first objective)
+ evaluation.evaluations[0] < best_candidate.as_ref().unwrap().evaluated_chromosome.evaluation.evaluations[0]
+ }
+ )?;
+
+ Ok((result, feasible_fractions))
+}
+
+const ITERATIONS: usize = 1000;
+const POPULATION: usize = 500;
+const PARENTS_COUNT: usize = 500;
+const G11_EPS: f64 = 0.00015;
+
+fn handle_g06_srank() -> Result<(), Box<dyn std::error::Error>> {
+ let problem = problem_g06();
+ let config = StochasticRankingConfig {
+ population_size: POPULATION,
+ parents_count: PARENTS_COUNT,
+ iterations: ITERATIONS,
+ n_param: 2 * POPULATION,
+ p_param: 0.45,
+ mutation_std_dev: 1.0,
+ };
+ run_stochastic_ranking(problem, config)
+}
+
+fn handle_g08_srank() -> Result<(), Box<dyn std::error::Error>> {
+ let problem = problem_g08();
+ let config = StochasticRankingConfig {
+ population_size: POPULATION,
+ parents_count: PARENTS_COUNT,
+ iterations: ITERATIONS,
+ n_param: 2 * POPULATION,
+ p_param: 0.45,
+ mutation_std_dev: 0.5,
+ };
+ run_stochastic_ranking(problem, config)
+}
+
+fn handle_g11_srank() -> Result<(), Box<dyn std::error::Error>> {
+ let problem = problem_g11(G11_EPS);
+ let config = StochasticRankingConfig {
+ population_size: POPULATION,
+ parents_count: PARENTS_COUNT,
+ iterations: ITERATIONS,
+ n_param: POPULATION * 2,
+
+ // population_size: POPULATION,
+ // parents_count: PARENTS_COUNT,
+ // iterations: ITERATIONS,
+ // n_param: 2 * POPULATION,
+
+ p_param: 0.05,
+ mutation_std_dev: 0.01 / 50.0,
+ };
+ run_stochastic_ranking(problem, config)
+}
+
+fn handle_g24_srank() -> Result<(), Box<dyn std::error::Error>> {
+ let problem = problem_g24();
+ let config = StochasticRankingConfig {
+ population_size: POPULATION,
+ parents_count: PARENTS_COUNT,
+ iterations: ITERATIONS,
+ n_param: 2 * POPULATION,
+ p_param: 0.45,
+ mutation_std_dev: 0.1,
+ };
+ run_stochastic_ranking(problem, config)
+}
+
+/// Generic function to save evolution results to CSV files with date-based naming
+fn save_evolution_results<const DIM: usize, const CONSTRAINTS: usize, TEval: std::fmt::Debug>(
+ method: &str,
+ problem: &ConstrainedProblem<DIM, CONSTRAINTS>,
+ evolution_result: &EvolutionResult<SVector<f64, DIM>, TEval>,
+ feasible_fractions: &[f64],
+) -> Result<(), Box<dyn std::error::Error>> {
+ // Get current date and time for unique file naming
+ let now = Local::now();
+ let timestamp = now.format("%Y%m%d_%H%M%S").to_string();
+
+ // Create output directories
+ let output_dir = format!("solutions/{}/{}", method, problem.name);
+ let feasible_dir = format!("{}/feasible_fraction", output_dir);
+ fs::create_dir_all(&output_dir)?;
+ fs::create_dir_all(&feasible_dir)?;
+
+ // Write best candidates CSV with timestamp
+ let best_candidates_path = format!("{}/best_candidates_{}.csv", output_dir, timestamp);
+ let mut best_file = fs::File::create(&best_candidates_path)?;
+ writeln!(best_file, "iteration,evaluation,fitness")?;
+
+ // Write evolution stats (best candidates through iterations)
+ for candidate in &evolution_result.stats.best_candidates {
+ writeln!(best_file, "{},{:?}", candidate.evaluation, candidate.evaluated_chromosome.evaluation)?;
+ }
+
+ // Write final best candidate and total evaluations
+ if let Some(ref best_candidate) = evolution_result.best_candidate {
+ writeln!(best_file, "{},{:?}", evolution_result.evaluations, best_candidate.evaluation)?;
+ }
+
+ // Write feasible fractions CSV with timestamp
+ let feasible_path = format!("{}/feasible_fractions_{}.csv", feasible_dir, timestamp);
+ let mut feasible_file = fs::File::create(&feasible_path)?;
+ writeln!(feasible_file, "iteration,feasible_fraction")?;
+
+ for (i, fraction) in feasible_fractions.iter().enumerate() {
+ writeln!(feasible_file, "{},{}", i, fraction)?;
+ }
+
+ println!("Results saved to:");
+ println!(" Best candidates: {}", best_candidates_path);
+ println!(" Feasible fractions: {}", feasible_path);
+
+ Ok(())
+}
+
+fn run_stochastic_ranking<const DIM: usize, const CONSTRAINTS: usize>(
+ problem: ConstrainedProblem<DIM, CONSTRAINTS>,
+ config: StochasticRankingConfig,
+) -> Result<(), Box<dyn std::error::Error>> {
+ let mut rng = rand::rng();
+
+ let result = solve_with_stochastic_ranking(
+ problem.clone(),
+ config.population_size,
+ config.parents_count,
+ config.iterations,
+ (config.n_param, config.p_param),
+ config.mutation_std_dev,
+ &mut rng,
+ )?;
+
+ let (evolution_result, feasible_fractions) = result;
+
+ // Save results to CSV files
+ save_evolution_results("srank", &problem, &evolution_result, &feasible_fractions)?;
+
+ if let Some(best_candidate) = evolution_result.best_candidate {
+ println!("Best solution found:");
+ println!(" Chromosome: {:?}", best_candidate.chromosome);
+ println!(
+ " Fitness: {} ({} %)",
+ best_candidate.evaluation,
+ ((best_candidate.evaluation - problem.optimal_value) / problem.optimal_value).abs() * 100.0);
+ println!(" Iterations: {}", evolution_result.iterations);
+ println!(" Evaluations: {}", evolution_result.evaluations);
+ println!(" Final feasible fraction: {:.2}%", feasible_fractions.last().unwrap_or(&0.0) * 100.0);
+
+ // Sanity check: verify the best candidate is feasible
+ let best_chromosome = &best_candidate.chromosome;
+
+ println!("\nFeasibility check for best solution:");
+ for (i, constraint) in problem.constraints.iter().enumerate() {
+ match constraint.evaluate(best_chromosome) {
+ Ok(value) => {
+ let is_feasible = value <= 0.0;
+ println!(" Constraint {}: {} ({})", i+1, value, if is_feasible { "FEASIBLE" } else { "INFEASIBLE" });
+ }
+ Err(e) => {
+ println!(" Constraint {}: Error evaluating - {}", i+1, e);
+ }
+ }
+ }
+ } else {
+ println!("Could not find any feasible solution!")
+ }
+
+ Ok(())
+}
+
+fn run_nsga_ii<const DIM: usize, const CONSTRAINTS: usize>(
+ problem: ConstrainedProblem<DIM, CONSTRAINTS>,
+ config: NsgaConfig,
+) -> Result<(), Box<dyn std::error::Error>> {
+ let mut rng = rand::rng();
+
+ let result = solve_with_nsga_ii(
+ problem.clone(),
+ config.population_size,
+ config.parents_count,
+ config.iterations,
+ config.mutation_std_dev,
+ &mut rng,
+ )?;
+
+ let (evolution_result, feasible_fractions) = result;
+
+ // Save results to CSV files
+ save_evolution_results("nsga", &problem, &evolution_result.clone().map(|x| x[0]), &feasible_fractions)?;
+
+ if let Some(best_candidate) = evolution_result.best_candidate {
+ println!("Best solution found:");
+ println!(" Chromosome: {:?}", best_candidate.chromosome);
+ println!(" Objectives: {:?} ({} %)", best_candidate.evaluation, ((best_candidate.evaluation[0] - problem.optimal_value) / problem.optimal_value).abs() * 100.0);
+ println!(" Iterations: {}", evolution_result.iterations);
+ println!(" Evaluations: {}", evolution_result.evaluations);
+ println!(" Final feasible fraction: {:.2}%", feasible_fractions.last().unwrap_or(&0.0) * 100.0);
+
+ // Sanity check: verify the best candidate is feasible
+ let best_chromosome = &best_candidate.chromosome;
+
+ println!("\nFeasibility check for best solution:");
+ for (i, constraint) in problem.constraints.iter().enumerate() {
+ match constraint.evaluate(best_chromosome) {
+ Ok(value) => {
+ let is_feasible = value <= 0.0;
+ println!(" Constraint {}: {} ({})", i+1, value, if is_feasible { "FEASIBLE" } else { "INFEASIBLE" });
+ }
+ Err(e) => {
+ println!(" Constraint {}: Error evaluating - {}", i+1, e);
}
- } else {
- println!("Could not find any feasible solution!")
}
}
- Err(e) => {
- println!("Error: {}", e);
+ } else {
+ println!("Could not find any feasible solution!")
+ }
+
+ Ok(())
+}
+
+fn run_nsga_multi<const DIM: usize, const CONSTRAINTS: usize, const CONSTRS_PLUS_ONE: usize>(
+ problem: ConstrainedProblem<DIM, CONSTRAINTS>,
+ evolution_result: EvolutionResult<SVector<f64, DIM>, [f64; CONSTRS_PLUS_ONE]>,
+ feasible_fractions: Vec<f64>,
+ capped: bool
+) -> Result<(), Box<dyn std::error::Error>> {
+ // Unfortunately Rustc doesn't support addition in generics...
+ assert_eq!(CONSTRAINTS + 1, CONSTRS_PLUS_ONE);
+
+ // Save results to CSV files
+ if capped {
+ save_evolution_results("nsga_multi", &problem, &evolution_result.clone().map(|x| x[0]), &feasible_fractions)?;
+ } else {
+ save_evolution_results("nsga_multi_noncapped", &problem, &evolution_result.clone().map(|x| x[0]), &feasible_fractions)?;
+ }
+
+ if let Some(best_candidate) = evolution_result.best_candidate {
+ println!("Best solution found:");
+ println!(" Chromosome: {:?}", best_candidate.chromosome);
+ println!(" Objectives: {:?} ({} %)", best_candidate.evaluation, ((best_candidate.evaluation[0] - problem.optimal_value) / problem.optimal_value).abs() * 100.0);
+ println!(" Iterations: {}", evolution_result.iterations);
+ println!(" Evaluations: {}", evolution_result.evaluations);
+ println!(" Final feasible fraction: {:.2}%", feasible_fractions.last().unwrap_or(&0.0) * 100.0);
+
+ // Sanity check: verify the best candidate is feasible
+ let best_chromosome = &best_candidate.chromosome;
+
+ println!("\nFeasibility check for best solution:");
+ for (i, constraint) in problem.constraints.iter().enumerate() {
+ match constraint.evaluate(best_chromosome) {
+ Ok(value) => {
+ let is_feasible = value <= 0.0;
+ println!(" Constraint {}: {} ({})", i+1, value, if is_feasible { "FEASIBLE" } else { "INFEASIBLE" });
+ }
+ Err(e) => {
+ println!(" Constraint {}: Error evaluating - {}", i+1, e);
+ }
+ }
+ }
+ } else {
+ println!("Could not find any feasible solution!")
+ }
+
+ Ok(())
+}
+
+// NSGA-II handler functions
+fn handle_g06_nsga() -> Result<(), Box<dyn std::error::Error>> {
+ let problem = problem_g06();
+ let config = NsgaConfig {
+ population_size: POPULATION,
+ parents_count: PARENTS_COUNT,
+ iterations: ITERATIONS,
+ mutation_std_dev: 1.0,
+ };
+ run_nsga_ii(problem, config)
+}
+
+fn handle_g08_nsga() -> Result<(), Box<dyn std::error::Error>> {
+ let problem = problem_g08();
+ let config = NsgaConfig {
+ population_size: POPULATION,
+ parents_count: PARENTS_COUNT,
+ iterations: ITERATIONS,
+ mutation_std_dev: 0.5,
+ };
+ run_nsga_ii(problem, config)
+}
+
+fn handle_g11_nsga() -> Result<(), Box<dyn std::error::Error>> {
+ let problem = problem_g11(G11_EPS);
+ let config = NsgaConfig {
+ population_size: POPULATION,
+ parents_count: PARENTS_COUNT,
+ iterations: ITERATIONS,
+ mutation_std_dev: 0.01 / 50.0,
+ };
+ run_nsga_ii(problem, config)
+}
+
+fn handle_g24_nsga() -> Result<(), Box<dyn std::error::Error>> {
+ let problem = problem_g24();
+ let config = NsgaConfig {
+ population_size: POPULATION,
+ parents_count: PARENTS_COUNT,
+ iterations: ITERATIONS,
+ mutation_std_dev: 0.1,
+ };
+ run_nsga_ii(problem, config)
+}
+
+// NSGA-Multi handler functions for individual problems
+fn handle_nsga_multi_g06(capped: bool) -> Result<(), Box<dyn std::error::Error>> {
+ let mut rng = rand::rng();
+ let problem = problem_g06();
+ let config = NsgaConfig {
+ population_size: POPULATION,
+ parents_count: PARENTS_COUNT,
+ iterations: ITERATIONS,
+ mutation_std_dev: 0.5,
+ };
+
+ let result = solve_with_nsga_multi::<2, 2, 3>(
+ problem.clone(),
+ config.population_size,
+ config.parents_count,
+ config.iterations,
+ config.mutation_std_dev,
+ &mut rng,
+ capped
+ )?;
+
+ let (evolution_result, feasible_fractions) = result;
+ run_nsga_multi(problem, evolution_result, feasible_fractions, capped)
+}
+
+fn handle_nsga_multi_g08(capped: bool) -> Result<(), Box<dyn std::error::Error>> {
+ let mut rng = rand::rng();
+ let problem = problem_g08();
+ let config = NsgaConfig {
+ population_size: POPULATION,
+ parents_count: PARENTS_COUNT,
+ iterations: ITERATIONS,
+ mutation_std_dev: 0.5,
+ };
+
+ let result = solve_with_nsga_multi::<2, 2, 3>(
+ problem.clone(),
+ config.population_size,
+ config.parents_count,
+ config.iterations,
+ config.mutation_std_dev,
+ &mut rng,
+ capped
+ )?;
+
+ let (evolution_result, feasible_fractions) = result;
+ run_nsga_multi(problem, evolution_result, feasible_fractions, capped)
+}
+
+fn handle_nsga_multi_g11(capped: bool) -> Result<(), Box<dyn std::error::Error>> {
+ let mut rng = rand::rng();
+ let problem = problem_g11(G11_EPS);
+ let config = NsgaConfig {
+ population_size: POPULATION,
+ parents_count: PARENTS_COUNT,
+ iterations: ITERATIONS,
+ mutation_std_dev: 0.01 / 50.0,
+ };
+
+ let result = solve_with_nsga_multi::<2, 1, 2>(
+ problem.clone(),
+ config.population_size,
+ config.parents_count,
+ config.iterations,
+ config.mutation_std_dev,
+ &mut rng,
+ capped
+ )?;
+
+ let (evolution_result, feasible_fractions) = result;
+ run_nsga_multi(problem, evolution_result, feasible_fractions, capped)
+}
+
+fn handle_nsga_multi_g24(capped: bool) -> Result<(), Box<dyn std::error::Error>> {
+ let mut rng = rand::rng();
+ let problem = problem_g24();
+ let config = NsgaConfig {
+ population_size: POPULATION,
+ parents_count: PARENTS_COUNT,
+ iterations: ITERATIONS,
+ mutation_std_dev: 0.1,
+ };
+
+ let result = solve_with_nsga_multi::<2, 2, 3>(
+ problem.clone(),
+ config.population_size,
+ config.parents_count,
+ config.iterations,
+ config.mutation_std_dev,
+ &mut rng,
+ capped
+ )?;
+
+ let (evolution_result, feasible_fractions) = result;
+ run_nsga_multi(problem, evolution_result, feasible_fractions, capped)
+}
+
+fn main() {
+ let args: Vec<String> = env::args().collect();
+
+ if args.len() != 3 {
+ eprintln!("Usage: {} <method> <problem>", args[0]);
+ eprintln!("Methods: srank, nsga, nsga_multi, nsga_multi_noncapped");
+ eprintln!("Problems: g06, g08, g11, g24");
+ std::process::exit(1);
+ }
+
+ let method = &args[1];
+ let problem = &args[2];
+
+ let result = match (method.as_str(), problem.as_str()) {
+ ("srank", "g06") => handle_g06_srank(),
+ ("srank", "g08") => handle_g08_srank(),
+ ("srank", "g11") => handle_g11_srank(),
+ ("srank", "g24") => handle_g24_srank(),
+ ("nsga", "g06") => handle_g06_nsga(),
+ ("nsga", "g08") => handle_g08_nsga(),
+ ("nsga", "g11") => handle_g11_nsga(),
+ ("nsga", "g24") => handle_g24_nsga(),
+ ("nsga_multi", "g06") => handle_nsga_multi_g06(true),
+ ("nsga_multi", "g08") => handle_nsga_multi_g08(true),
+ ("nsga_multi", "g11") => handle_nsga_multi_g11(true),
+ ("nsga_multi", "g24") => handle_nsga_multi_g24(true),
+ ("nsga_multi_noncapped", "g06") => handle_nsga_multi_g06(false),
+ ("nsga_multi_noncapped", "g08") => handle_nsga_multi_g08(false),
+ ("nsga_multi_noncapped", "g11") => handle_nsga_multi_g11(false),
+ ("nsga_multi_noncapped", "g24") => handle_nsga_multi_g24(false),
+ (_, _) => {
+ eprintln!("Invalid method '{}' or problem '{}'", method, problem);
+ eprintln!("Methods: srank, nsga, nsga_multi, nsga_multi_noncapped");
+ eprintln!("Problems: g06, g08, g11, g24");
+ std::process::exit(1);
}
+ };
+
+ if let Err(e) = result {
+ eprintln!("Error: {}", e);
+ std::process::exit(1);
}
}
@@ 4,392 4,175 @@
# Intro
-This is a report for hw01 for EOA course on CTU FEE. The goal
-has been to try to solve traveling salesperson problem by
-means of evolutionary algorithms.
+This is a report for hw02 for EOA course on CTU FEE. The goal
+has been to compare constrained evolutionary algorithms with two approaches:
-The traveling salesperson problem tries to find the shortest possible
-route when traveling through multiple cities, visiting every city
-exactly once, returning to the origin city.
+1. Using stochastic ranking,
+2. Using multi objective evaluation with the original fitness as first objective
+and weighted sum of constraints as second argument.
The report covers the implemented
-algorithms and the results on 10 TSP instances from TSPLIB. All of
-those chosen instances are using euclidean metric and are 2D.
+algorithms and the results on problems taken from "Problem Definitions and Evaluation Criteria for the CEC 2006"
-{latex-placement="ht"}
+# Used algorithms
-# Prepared algorithms
-Three generic algorithms have been used for solving the TSP.
-Then various representations and operators have been tried on them.
+## Stochastic ranking
-## Evolutionary algorithm
-The algorithm keeps a population of N each epoch and generates offsprings
-out of the current population. Specifically it uses the following steps:
+One of the methods to solve constrained problems.
-1. Obtain initial population (i.e., random or specialized heuristic).
-2. Select parents to make offsprings from (selection).
-3. Pair the selected parents (pairing).
-4. Somehow join the parent pairs to make an offspring (crossover).
-5. Replace the current population with a new one, combining the current population and offsprings. (replacement)
-6. Repeat steps 2 to 5 until terminating condition is reached (number of iterations)
+Stochastic ranking uses similar structure as bubble sort. It iterates through the population
+and swaps to elements when one is better than the other (lower fitness).
+But it does that only with probability `p`. With probability `1 - p` it instead compares
+the constraint violations of both elements.
-During its run the algorithm saves the best candidates encountered. It saves them even during
-one epoch for better granularity. This also ensures that if best candidate does not make it
-into the next population, it is still captured by the algorithm and can be returned as
-best candidate at the end. The implementation is in `eoa_lib/src/evolution.rs`.
+Then, the selection takes into account ranking based on the list sorted through stochastic ranking.
+Generational replacement has been used to easily ensure that the stochastic ranking order is still respected.
-For selections those algorithms have been implemented and tried:
+Stochastic ranking is implemented in `eoa_lib/src/constraints.rs`
-- Roulette wheel selection
-- Tournament selection
- - Here two parameters have been used, number of parents to choose from and probability to choose the best parent
- - When probability is set at 1, always the parent with best fitness wins, but if it's lower, it gives chance to worse parents as well
-- Best by fitness selection
+## Unconstrained multi objective evolution (NSGA-II)
-For replacmement, tournament and best by fitness has been implemented
-and tried.
+As second method, NSGA-II has been used to solve the constrained problems.
+It works based on non-dominated fronts.
+All elements in one front are not dominated between each other. A dominates B if
+and only if A is not worse than B in all objectives and A is better than B in at least
+one objective.
-For pairing, only adjacent pairing has been implemented (taking parents
-from selection and pairing first with second, third with fourth and so on)
+NSGA-II first sorts the elements based on the non dominated fronts. Then, within each front,
+crowding distance is calculated. Crowding distance is used to discourage solutions that are
+close to each other. For both selection and replacement the comparison is then done first according
+to the front number and within same front, according to the crowding distance. Higher crowding
+distance is preferred.
-## Local search
-The local search implementation is implementation of the first improving strategy,
-where the current best solution is being perturbated, until a better solution is found.
-And then, the process starts over from this better solution.
-
-The local search implementation supports evolutionary strategies,
-i.e., changing the perturbation parameters during the run, but none
-were used for TSP.
-
-The implementation is available in `eoa_lib/src/local_search/mod.rs`.
+As part of NSGA-II, tournament selection is used and for replacement, best replacement is used,
+keeping only the best individuals. Specific parameters of using 5 candidates for tournament selection
+and probability of 0.99 for choosing the best candidate out of 5.
-## Random search
-Random search initializes new random elements in the search space each
-iteration and saves new elements if they are better than best found.
-
-The implementation is available in `eoa_lib/src/random_search.rs`.
-
-# Representations
-
-Two representations have been chosen for the TSP problem,
-
-1. Node permutation
-2. Binary string half-matrix denoting if city i comes before city j
+Multi objective evolution is implemented in `eoa_lib/src/multi_objective_evolution.rs`
-## Node permutation
-Implemented as a vector of indices of the cities.
-A couple of perturbations and crossovers have been implemented:
+# Used operators
-- Perturbations
- - Move single random city to random position
- - Swap two randomly selected cities
- - Reverse subsequence of cities (random starting point and random length)
-- Crossovers
- - Cycle crossover
- - Partially mapped crossover
- - Edge recombination crossover
-
-Also apart from a random initializer, two additional initializers are implemented,
-one based on minimal spanning tree and second based on nearest neighbors. Detailed
-descriptions of the algorithms follow in next sections.
-
-### Crossovers
-All used crossovers take two parents and it's possible to create two
-offsprings out of the two parents by switching the parent positions.
-(switching first parent and second parent)
-
-#### Cycle crossover
-The cycle crossover creates a cycle of indices, takes parts of the first parent
-from the cycle and moves them to offspring and fills the rest from the second parent.
-
-The cycle is created as follows:
-
-1. Start with index 0, call it current index
-2. Save current index
-3. Look at current index in first parent, let's call it current element
-4. Find the same element in second parent and update current index to this new index
-5. Repeat 2 - 4 until reaching index 0 again
-
-Then the offspring is created as follows:
-
-1. Clone the second parent
-2. Iterate the saved cycle indices, take element at index from first parent and copy it to offspring at the same index
-
-#### Partially mapped crossover
-The partially mapped/matched crossover randomly selects two points. At the end it should
-ensure that the offspring has first parent's elements in between the two cross points.
-
-The way to ensure that, while still ensuring the result is a valid permutation, is to
-always swap the elements.
-
-Offspring is created as follows:
-
-1. Clone second parent
-2. Then, for every index i between the cross points:
-3. Take the element on index i from first parent
-4. Find the city in the offspring, let's call its index j
-5. Swap i with j
-
-#### Edge recombination crossover
-Edge recombination is the most complicated from the three crossover operators.
-
-First, an adjacency list is created for both parents. Let's take an example
-of a permutation: 012345.
-The element 0 has 5 and 1 as adjacent nodes. 1 has 0 and 2 and so on.
-The two adjacency lists are then joined together, while omitting repetitions.
-That means that a single element will have from 2 to 4 elements in its adjacency list.
-
-To make an offspring:
-
-1. First element is taken out of the first parent, this is the current element (city).
-2. The current element is appended to offspring
-3. The current element is removed from adjacency lists of all other elements.
-4. Then, all the cities in adjacency list for current element are taken and one with minimum neighbors in its adjacency list is taken. If there are multiple such cities, random one is chosen.
-5. The found city becomes current element.
-6. Steps 2 to 5 are repeated until all elements are taken.
-
-In the code, the algorithm has been implemented slightly differently.
-Instead of removing any elements, the elements are just marked as removed.
-Also, instead of constructing a list of adjacency nodes progressively,
-a vector of 4 elements is made beforehand. Each element may be unset (an Option).
-Then, the adjacencies from first parent are put to positions 0, 1 and adjacencies from
-second parent to positions 2 and 3. If there would be a repetition, it's left unset.
-This then allows for static allocation only. This is thanks to a library called nalgebra
-that is used and it's possible to tell it the dimensions of adjacency matrix
-beforehand. However for the runs of the algorithms, `Dyn` dimension has been used
-out of convenience, so dynamic allocation is performed.
-
-## Binary string
-Classical perturbations and crossovers have been
-implemented for the binary string representation, specifically:
-
-- Perturbations
- - Flip each bit with probability p
- - Flip single random bit
- - Flip of N bits (N is chosen beforehand, not randomly)
- - Flip of whole binary string
-
-- Crossover
- - N point crossover
-
-As for initialization, random initializer has been implemented.
-
-The fitness function is implemented by form of a wrapper that converts
-the BinaryString into NodePermutation and then the same fitness
-function is used as for node permutation representation.
-
-### N-point crossover
-The N point crossover works on two parents and is
-capable of producing two offsprings.
-The crossover first chooses N cross points randomly.
-
-Then, the cross points are ordered in an ascending order
-and first all bits from first parent are chosen until
-cross point is encountered, then all bits are taken from
-second parent until another cross point is reached. And
-this repeats until the end of the string is reached.
-
-Also, one-point and two-point crossovers
-have been implemented separately for more effective implementation
-than the generic N-point.
-
-# Heuristics
-Instead of starting with random solutions, two heuristics have been tried to
-make the initial populations for the evolutionary algorithms.
-
-## Nearest neighbors
-The heuristic starts at a given node and finds its nearest neighbor. Then
-adds that neighbor to the permutation and moves to it. Then it repeats search
-for the nearest neighbor, making sure to not select nodes twice. The whole
-chromosome is built like this.
-
-Moreover the possibility to select second neighbor instead of first has been
-incorporated in the heuristic as well. Specifically, it's possible to choose
-the probability to choose second neighbor instead of first one.
-
-This makes it possible to generate a lot of initial solutions, specifically it's
-possible to generate nearest neighbors starting from each city and then it's possible
-to tweak the probability of choosing the second neighbor.
-
-To initialize the whole population:
-1. Generate chromosome starting from each city, choosing the first neighbor.
-2. Generate chromosome starting from each city, but always choosing the second neighbor.
-3. The rest of the population is initialized from randomly selected cities with randomly selected probability of choosing second neighbor.
-
-## Minimal spanning tree
-For using a minimal spanning tree the algorithm of Christofides and Serdyukov has been
-chosen as an inspiration. But instead of trying to find best Eulerian tour on the minimal
-spanning tree, a slightly different randomized approach has been taken. Specifically,
-random node is chosen as the starting point and then edges are chosen out of the minimal
-spanning tree to go through. If there are more edges to travel through, random one is chosen.
-If there are no edges left, nearest neighbor is chosen as next node.
-
-To initialize the whole population, the same algorithm is ran multiple times. Since
-it is random, its results should be different, at least slightly.
+Since the chosen problems are real problems, real representation has been used.
-# Results
-To compare all the algorithms on various instances, 10 runs of each algorithm
-have been performed on each instance. All the graphs of probabilities of reaching
-target for given function evaluation were then constructed from
-averaging between the runs and between the instances. The fitness charts sometimes show fewer runs to not be too
-overwhelming. All evolutionary algorithms ran for 10000 iterations (epochs) and the local search
-and random search were adjusted to run on the same number of fitness function evaluations.
-The population size for the evolutionary algorithm is 500.
-The number of parents (and thus offspring) is 250.
-
-On the x-axis of all graphs is the number of evaluations of the fitness function.
-
-The graphs with probability of reaching targets have been computed as follows:
-
-1. Individually computed graphs for targets 1 %, 5 %, 10 % and 20 %
-2. Averaged the graphs
-
-To compute the individual graphs, at each evaluation point, we count how many runs have already achieved better results than the target and divide this by the total number of runs.
-
-Then, standard deviation has been added, clamped between 0 and 1.
-
-The instances chosen from TSPLIB, all 2D Euclidean:
-
-- eil51
-- eil76
-- kroA100
-- eil101
-- pr124
-- ch150
-- kroA150
-- kroA200
-- ts225
-- a280
-
-For all the evolutionary algorithm runs, roulette wheel has been used as selection and truncation of best has
-been selected for replacement. This is because they seemed to provide the best results.
-
-During runs of the EA, the probabilities of mutations were changed throughout the run, by this formula:
-`P = 0.5 * (1.0 + b / e)`, where `b` is iterations since better solution has been found and `e` is number
-of iterations till the end (`max_iterations - current_iteration`), this makes the mutations more aggressive
-later in the algorithm, as long as no better solutions are found, reaching even 100 % for some cases.
+For mutation, the solution is changed with random sample from normal distribution
+with zero mean.
+The standard mean has been chosen differently based on the problems.
-\newpage
+As for crossover, arithmetical crossover is used.
+Arithmetical crossover chooses random `\alpha` between
+0 and 1. Then, the offspring's i-th component is
+`\alpha \cdot p_{1i} + (1 - \alpha) \cdot p_{2i}`,
+where `p_{1i}` (`p_{2i}`) is first (second) parent's i-th
+component respectively.
-## Comparing algorithms
-To compare the algorithms,
-first it has been ensured the algorithms were tweaked to produce the best results (best that the author
-has been capable of). Then, they were run on 10 instances of TSP and averaged in the following chart:
+# Making NSGA-II behave better
-{latex-placement="H"}
-\newpage
+In the default implementation, NSGA-II is not suitable for solving constrained
+problems. This is because all the feasible individuals will be at 0 in the
+objective with constraints sum.
-Here is a an averaged fitness sample on eil76 instance:
+Because of this, a very simple idea can arise to make it better - if multiple feasible solutions
+are sufficiently distant, their constraint function evaluations could also differ
+if they are allowed to go below zero. That feasible solutions can have higher crowding distances.
+But this seems to make sense only for cases where there is objective per constraint. Otherwise
+a constraint that is far below zero would win over and the other constraints might not have chance
+to get satisfied.
-{latex-placement="H"}
-
-From the results, the LS is capable of finding a better solution early on, but then, the EA
-catches up and produces a slightly better solution than LS. LS converges earlier than EA.
-
-\newpage
-
-## Comparing perturbations on LS
-Four perturbations have been evaluated:
-
-- Moving single city
-- Swapping of two cities
-- Reversing subsequence
-- Combination of the above (each cycle, single random perturbation has been chosen to run from the above three)
-
-{latex-placement="H"}
-
-\newpage
-
-{latex-placement="H"}
-
-from the experiment it seems the reversal of a subsequence is the best singular perturbation,
-but the combination is even better. That could be explained by the fact that in cases the subsequence
-reversal gets stuck somewhere, the other perturbations are able to recover.
+# Results
-\newpage
+Every problem with given configuration has been ran 10 times. The initialization
+has been done randomly, but within given bounds different for each problem:
-## Comparing crossovers
-During evaluation of the various crossovers, it has become apparent that with the currently
-chosen parameters of the algorithms, the dominant parts that are responsible for the good
-convergence are: selection, perturbation and replacement, but not the crossover itself.
-It was even tried to use no crossover at all, passing through the parents and the convergence
-has been similar.
-For this reason, the mutations probabilities have been lowered significantly for this comparison.
-This then allows for comparison between the crossovers themselves. Specifically the probability
-of mutation is 0.001.
+- g06 - `0, 0` `50, 50`
+- g08 - `0, 0` `10, 10`
+- g11 - `-1, -1` `1, 1`
+- g24 - `0, 0` `3, 4`
-{latex-placement="H"}
+First also bounding the crossover result and mutation has been tried,
+but it turned out it doesn't even help much, so that has been removed.
+(the method for bounding has been to first try the mutation/crossover
+maximally 5 times if it's out of bounds and after last try, it has
+been bounded forcibly)
-From the runs, the edge recombination has produced the best results, with partially
-mapped crossover being the second and cycle crossover being the worst. This could be explained
-by the fact that edge recombination keeps all of the edges of parents, while mapped crossover makes sure
-to keep edges only inside the two cross point interval. And cycle crossover might miss even more edges.
+For stochastic ranking approach, following parameters have been used:
+arithmetic crossover, normal distribution mutation with probability 0.1,
+tournament selection with 5 individuals and 0.95 probability
+of choosing the best one. Weights of 1e9 have been chosen as with 1.0 there
+have been problems with g11, probably due to very small numbers that were being
+summed together and rounded.
-While at the end CX and PMX are cathing up, I think this is mostly thanks to the mutations rather than
-the crossovers themselves. When mutations are turned off completely, PMX and CX never reach targets as high as ERX.
+As for runtime parameters, population had 500 individuals, there are 1000 iterations,
+500 parents in each iteration. The number of iterations in stochastic ranking is 2 * population.
+Probability of comparing fitness is 0.45. This is all with exception of g11, where different parameters
+have been used to so that feasible solution is found. Specifically epsilon has been chosen to be `0.01`,
+population size 1250, 1000 parents, 500 iterations, 2000 for number of iterations of stochastic
+ranking and the probability 0.05. Still, feasible solution has not been found on every run. The results
+include only runs where solution has been found.
+The following standard deviations have been chosen:
-\newpage
+- g06 - `1.0`
+- g08 - `0.5`
+- g11 - `0.01 / 50` (epsilon divided by 50)
+- g24 - `0.1`
-## Comparing representations
+For NSGA-II
+TODO
-{latex-placement="H"}
+## Comparing the approaches
-{latex-placement="H"}
+The following plot compares:
-On one hand the binary string representation is simpler, because it doesn't require special treatment
-for crossovers. Any crossover leads to valid result. On the other hand since it's not so specific to the
-given problem, it seems to lead to worse results than the node permutation. It seems that mainly the reverse
-subsequence perturbation is making node permutation representation that much better,
-as can be seen from comparison of various LS perturbations.
+- Stochastic ranking
+- NSGA-II
+- NSGA-II with multiple objectives
+- NSGA-II with multiple objectives, allowing constraints under 0
-On top of its worse performance, the whole run is also slower for the same number of iterations, but this
-might be just due to the poor implementation technique chosen for the binary string conversion to node
-permutation.
+As can be seen, NSGA with multiple objectives behaves worse. This is
+possibly because the algorithm can lower second and third objective
+more easily than lowering the first one or lowering all three at the same
+time. This is then even more visible for allowing the constraints to
+go under 0. So the simple idea has failed.
-\newpage
+TODO add the plot
-## Comparing heuristics
-Both of the heuristics are capable of making initial solutions much better
-than random ones. Here is a fitness graph with only the two heuristics.
+## Comparing g11 with g06
-{latex-placement="H"}
+`g11` is special in the regard that it has an
+equality condition. Because of that, it can make
+sense to compare it with other ones. Here's the result:
-From the results it can be seen the nearest neighbor heuristic performs slightly better
-for the initialization, but at the end both heuristics lead to very similar results.
-If the heuristic with minimal spanning tree has been enhanced by searching for better
-solutions instead of randomly, it might outperform nearest neighbor.
+As can be observed, `g11` cannot be optimized easily
+as most of the solutions are going to be infeasible.
+Many of the runs, no feasible solutions are found. And
+since there is only a small range of values to go to when
+in feasible space, it can easily happen mutated solution
+again goes outside of feasible space. Because of this, the
+`p` parameter had to be significantly lowered and this leads
+to optimizing for feasible solutions more than for better solutions.
-The evolutionary algorithm still provides a visible improvement even for the heuristics
-that are already very good starting points.
-Moreover it has to be noted that the results show that with those heuristics, it's possible
-to obtain a better result with the evolutionary algorithm overall that hasn't been reached
-with algorithm with random initialization. This shows possible problem with the configuration,
-the population stagnates. It could be possible to come up with ways to make the population not
-stagnate so much so that the EA algorithm doesn't converge so early.
+Here is a plot with fraction of feasible solutions in iteration for both
+problems:
-# Things done above minimal requirements
+# Tasks done above minimal requirements
-- 1 point - Compared 2 representations (permutation and binary string with precedences)
-- 1 point - Compared 4 LS perturbations (swap, reverse subsequence, move city, combination)
-- 1 point - Compared 3 crossover operators (edge recombination, partially mapped, cycle)
-- 1 point - Compared all algorithms on 10 instances
-- 1 point - Initialized with two constructive heuristics (nearest neighbors, solutions from minimal spanning tree)
+- 1 point - Multi-objective approach with more objectives
+- Tried to improve NSGA with very simple approach, but failed
# Code structure
-Rust has been chosen as the language. There are three subdirectories, `eoa_lib`, `tsp_hw01` and `tsp_plotter`.
+Rust has been chosen as the language. There are four subdirectories, `eoa_lib`, `tsp_hw01`, `constr_hw02` and `tsp_plotter`.
`eoa_lib` is the library with the generic operators defined, with random search, local search and evolution
algorithm functions. It also contains the most common representations, perturbations and crossovers for them.
-`tsp_hw01` contains the TSP implementation and runs all of the algorithms. It then produces CSV results
-with best candidate evaluations for given fitness function evaluation count. This is then utilized by
-`tsp_plotter`. The node permutation TSP representation itself is implemented in `tsp.rs`. The configurations
-of all algorithms used are located in `main.rs`. All the instances used are located in `tsp_hw01/instances` and
-the solutions are put to `tsp_hw01/solutions`.
+`constr_hw02` is the implementation of hw02. It contains the problem definitions and a simple CLI
+program to run any of the problems with either stochastic ranking or with MOE.
+
+`tsp_plotter` plots the graphs. It can now be used forjboth homework.
-`tsp_plotter` contains hard-coded presets for charts to create for the report.
+`tsp_hw01` is folder with source of the first folder, kept for my convenience as I keep them in the same repo
+and would have to be changing Cargo.toml properly to make sure the build doesn't fail if it's not present.
-As part of the work, apart from the standard library, third party libraries have been used, namely:
+As part of the work, apart from the standard library, third party libraries have been used, namely (the list is the same as for hw01):
- nalgebra for working with vertices and matrices,
- plotters for plotting
- rand
@@ 404,39 187,14 @@ For run instructions, including generation of the plots in this report, see READ
# Usage of LLM
-While I was working on this homework, I have used LLM for writing certain parts of the code,
-specifically I have used it very extensively for routine tasks:
-- Loading data,
-- Plotting graphs,
-- Refactoring of a triviality at a lot of places at once (i.e., at first I chose to make perturbation copy the chromosome,
- but I realized this is very ineffective for EAs afterwards and changed it to change in-place instead),
-- Writing some of the tests
-
-As I am not proficient with Rust, sometimes I asked LLM to help with the syntax as well,
-to find a library that will help solve a task or what functions are available for a specific
-task.
-
-I have used LLM only minimally for implementing the algorithms or for deciding on how
-to make the implementation, mainly sometimes in the form of checking if the algorithm
-looks correct. It did find a few issues that I then sometimes let it fix and
-sometimes fixed myself. This is because I believe that I can learn the most by
-writing the implementations myself.
-
-I use Claude from within claude-code, a CLI tool that is capable
-of reading files, executing commands and giving the model live feedback,
-like outputs of ran commands. This means the LLM is able to iterate by
-itself, without my intervention, for example when fixing errors.
-On the other hand it can use a lot of tokens quickly unless I make sure
-to run compact often.
+LLM has been used similarly to how in HW01. The only thing that changed
+is that Gemini has also been utilized, not only Claude.
+(Students got it now for free for a year :))
# Resources
- Lectures of EOA
-- Claude LLM
-- TSPLIB (http://comopt.ifi.uni-heidelberg.de/software/TSPLIB95/)
-- Genetic Algorithm Solution of the TSP Avoiding Special Crossover and Mutation from Gokturk Ucoluk (https://cw.fel.cvut.cz/wiki/_media/courses/a0m33eoa/du/ucoluk2013tspnew.pdf)
-- Segregation of Various Crossover Operators in TSP using G.A by Jyoti Naveen and Rama Chawla (https://cw.fel.cvut.cz/b231/_media/courses/a0m33eoa/du/naween2013crossoversfortsp.pdf)
-- https://en.wikipedia.org/wiki/2-opt
-- https://en.wikipedia.org/wiki/Edge_recombination_operator
-- https://ultraevolution.org/blog/pmx_crossover/
-- https://en.wikipedia.org/wiki/Crossover_(evolutionary_algorithm)
+- Claude, Gemini LLMs
+- https://sciences.uodiyala.edu.iq/uploads/00%20Abdullah%20New%20website/Lectures/comp/%D9%85%D8%AD%D8%A7%D8%B6%D8%B1%D8%A7%D8%AA%20%D8%AF%20%D8%B2%D9%8A%D8%A7%D8%AF%20%D8%B7%D8%A7%D8%B1%D9%82/%D8%A7%D8%AD%D8%AA%D8%B3%D8%A7%D8%A8%20%D8%AA%D8%B7%D9%88%D8%B1%D9%8A%203.pdf
+- https://www.researchgate.net/figure/Blend-crossover-operator-BLX_fig1_226044085
+- https://cw.fel.cvut.cz/wiki/_media/courses/a0m33eoa/cviceni/2006_problem_definitions_and_evaluation_criteria_for_the_cec_2006_special_session_on_constraint_real-parameter_optimization.pdf