From db937f4c8c792edfa6e98b20bc44f5036bc64699 Mon Sep 17 00:00:00 2001 From: Rutherther Date: Fri, 5 Dec 2025 22:07:18 +0100 Subject: [PATCH] feat: implement improved nsga for constraints with archive of feasible solutions --- codes/constr_hw02/src/main.rs | 363 ++++++++++++++++-- .../eoa_lib/src/multi_objective_evolution.rs | 19 +- 2 files changed, 344 insertions(+), 38 deletions(-) diff --git a/codes/constr_hw02/src/main.rs b/codes/constr_hw02/src/main.rs index a3ab28d4b3c4118cc0e3f6e1230c673c9ecd9ebe..c025c429b40de8e3ddcd0d7fdb421a7b31ee6732 100644 --- a/codes/constr_hw02/src/main.rs +++ b/codes/constr_hw02/src/main.rs @@ -1,4 +1,4 @@ -use std::{convert::Infallible, env, fs, io::Write, rc::Rc}; +use std::{convert::Infallible, env, fs, io::Write, rc::Rc, cell::RefCell}; use eoa_lib::{ bounded::BoundedOVector, comparison::MinimizingOperator, constraints::{stochastic_ranking_evolution_algorithm, ConstrainedEvalFitness, ConstrainedFitnessFunction, ConstraintFunction, LowerThanConstraintFunction}, crossover::{ArithmeticCrossover, BoundedCrossover, BoundedCrossoverStrategy, Crossover}, evolution::{EvolutionCandidate, EvolutionResult, EvolutionStats}, fitness::FitnessFunction, initializer::{Initializer, RandomInitializer}, multi_objective_evolution::{constrained_nsga_2, nsga_2, NSGAEvaluation}, pairing::{AdjacentPairing, ParentPairing}, perturbation::{BoundedPerturbation, BoundedPerturbationStrategy, MutationPerturbation, RandomDistributionPerturbation}, population::{EvaluatedChromosome, EvaluatedPopulation, Population}, replacement::{BestReplacement, GenerationalReplacement}, selection::{BestSelection, TournamentSelection} @@ -117,34 +117,37 @@ impl< let a_feasible = a.evaluation.evaluations.iter().skip(1).all(|&constr| constr <= 0.0); let b_feasible = b.evaluation.evaluations.iter().skip(1).all(|&constr| constr <= 0.0); - match (a_feasible, b_feasible) { - (false, true) => { - if rng.random_bool(self.p_single_replaced) { - pair[0] = rng.random_range(parents_count..full_population); + // Only proceed with replacement if we have archived feasible solutions + if full_population > parents_count { + match (a_feasible, b_feasible) { + (false, true) => { + if rng.random_bool(self.p_single_replaced) { + pair[0] = rng.random_range(parents_count..full_population); + } + }, + (true, false) => { + if rng.random_bool(self.p_single_replaced) { + pair[1] = rng.random_range(parents_count..full_population); + } + }, + (false, false) => { + if rng.random_bool(self.p_double_first_replaced) { + pair[0] = rng.random_range(parents_count..full_population); + + if rng.random_bool(self.p_double_second_replaced) { + pair[1] = rng.random_range(parents_count..full_population); + } + } } - }, - (true, false) => { - if rng.random_bool(self.p_single_replaced) { - pair[1] = rng.random_range(parents_count..full_population); + (true, true) => { + // Do nothing. } - }, - (false, false) => { - if rng.random_bool(self.p_double_first_replaced) { - pair[0] = rng.random_range(parents_count..full_population); - if rng.random_bool(self.p_double_second_replaced) { - pair[1] = rng.random_range(parents_count..full_population); - } - } - } - (true, true) => { - // Do nothing. } - } } - self.crossover(&joined_population, new_pairs.into_iter(), rng) + self.crossover.crossover(&joined_population, new_pairs.into_iter(), rng) } } @@ -677,7 +680,7 @@ pub fn solve_with_nsga_ii( &better_than, iterations, rng, - |_iteration: usize, _stats: &EvolutionStats, _>, population: &EvaluatedPopulation, _>| { + |_iteration: usize, _stats: &EvolutionStats, _>, population: &EvaluatedPopulation, _>, _crossover: &mut _| { // Calculate feasible fraction based on second objective being close to zero let feasible_count = population.population .iter() @@ -775,7 +778,7 @@ pub fn solve_with_nsga_multi, _>, population: &EvaluatedPopulation, _>| { + |_iteration: usize, _stats: &EvolutionStats, _>, population: &EvaluatedPopulation, _>, _crossover: &mut _| { // Calculate feasible fraction - all constraints (objectives 1 and 2) must be <= 0 let feasible_count: f64 = population.population.iter().filter( @@ -893,7 +896,7 @@ pub fn solve_with_nsga_constr( &better_than, iterations, rng, - |_iteration: usize, _stats: &EvolutionStats, _>, population: &EvaluatedPopulation, _>| { + |_iteration: usize, _stats: &EvolutionStats, _>, population: &EvaluatedPopulation, _>, _crossover: &mut _| { // Calculate feasible fraction and average constraint violation let feasible_count = population.population .iter() @@ -941,8 +944,149 @@ pub fn solve_with_nsga_constr( Ok((result, feasible_fractions, avg_constraint_violations)) } +pub fn solve_with_nsga_improved( + problem: &ConstrainedProblem, + population_size: usize, + parents_count: usize, + iterations: usize, + mutation_std_dev: f64, + rng: &mut dyn RngCore, +) -> Result<(EvolutionResult, [f64; 2]>, Vec, Vec), Box> { + // 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::, population_size, rng) + ); + + // Setup components + let mut pairing = AdjacentPairing::new(); + + // Create the wrapped crossover with arithmetic crossover inside + let mut wrapped_crossover = FeasibleCrossoverWrapper { + p_single_replaced: P_SINGLE_REPLACED, + p_double_first_replaced: P_DOUBLE_FIRST_REPLACED, + p_double_second_replaced: P_DOUBLE_SECOND_REPLACED, + archived_count: ARCHIVE_SIZE, + archived_population: Vec::new(), + crossover: ArithmeticCrossover::new(), + }; + + // Setup bounded random distribution perturbation with Normal distribution + let normal_perturbation = RandomDistributionPerturbation::>::normal(mutation_std_dev)?; + let mut mutation = MutationPerturbation::new(Box::new(normal_perturbation), 0.1); + + let better_than = MinimizingOperator::new(); + + // Clone the problem to access its fields + let cloned_problem = problem.clone(); + + // Create objectives: following the nsga algorithm pattern + let constraint_weights = [1.0e9; CONSTRAINTS]; + let constraint_refs = problem.constraints.iter().collect::>().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, Out = f64, Err = Infallible>>; 2] = [ + Box::new(cloned_problem.objective), + Box::new(constrained_fitness_obj2), + ]; + + let mut feasible_fractions = Vec::with_capacity(iterations); + let mut avg_constraint_violations = Vec::with_capacity(iterations); + + let result = nsga_2::<2, SVector, f64, Infallible, 2, _, _, _>( + initial_population, + parents_count, + objectives, + &mut pairing, + &mut wrapped_crossover, + &mut mutation, + &better_than, + iterations, + rng, + |_iteration: usize, _stats: &EvolutionStats, _>, population: &EvaluatedPopulation, _>, crossover: &mut _| { + // Update archive with feasible solutions + let feasible_individuals: Vec<_> = population.population.iter() + .filter(|individual| { + problem.constraints.iter().all(|constraint| { + constraint.is_feasible(&individual.chromosome).unwrap_or(false) + }) + }) + .cloned() + .collect(); + + // Update the crossover's archive now that we have access to it + crossover.update_archive(population); + + // Calculate feasible fraction and average constraint violation + let feasible_count = population.population + .iter() + .filter(|individual| { + // Check if the individual satisfies all constraints using the actual constraint functions + problem.constraints.iter().all(|constraint| { + constraint.is_feasible(&individual.chromosome).unwrap_or(false) + }) + }) + .count(); + + // Calculate average constraint violation + let total_violation: f64 = population.population + .iter() + .map(|individual| { + problem.constraints + .iter() + .map(|constraint| constraint.evaluate(&individual.chromosome).unwrap_or(0.0).max(0.0)) + .sum::() + }) + .sum(); + let avg_constraint_violation = total_violation / population.population.len() as f64; + + let feasible_fraction = feasible_count as f64 / population.population.len() as f64; + feasible_fractions.push(feasible_fraction); + avg_constraint_violations.push(avg_constraint_violation); + }, + |chromosome, evaluation, best_candidate| { + // Only save feasible solutions using actual constraint functions + if !problem.constraints.iter().all(|constraint| { + constraint.is_feasible(chromosome).unwrap_or(false) + }) { + 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] + } + )?; + + // Return both objectives (original objective and constraint sum) + let objective_result = result; + + Ok((objective_result, feasible_fractions, avg_constraint_violations)) +} + const ITERATIONS: usize = 1000; const POPULATION: usize = 500; + +// FeasibleCrossoverWrapper global probability parameters +const P_SINGLE_REPLACED: f64 = 0.4; +const P_DOUBLE_FIRST_REPLACED: f64 = 0.6; +const P_DOUBLE_SECOND_REPLACED: f64 = 0.3; +const ARCHIVE_SIZE: usize = 100; const PARENTS_COUNT: usize = 500; const G11_EPS: f64 = 0.00015; @@ -1402,6 +1546,56 @@ fn run_nsga_constr( Ok(()) } +fn run_nsga_improved( + problem: ConstrainedProblem, + config: NsgaConfig, +) -> Result<(), Box> { + let mut rng = rand::rng(); + + let result = solve_with_nsga_improved( + &problem, + config.population_size, + config.parents_count, + config.iterations, + config.mutation_std_dev, + &mut rng, + )?; + + let (evolution_result, feasible_fractions, avg_constraint_violations) = result; + + // Save results to CSV files + save_evolution_results("nsga_improved", &problem, &evolution_result.clone().map(|x| x[0]), &feasible_fractions, &avg_constraint_violations)?; + + 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-Constr handler functions fn handle_g04_nsga_constr() -> Result<(), Box> { let problem = problem_g04(); @@ -1491,6 +1685,95 @@ fn handle_g24_nsga_constr() -> Result<(), Box> { run_nsga_constr(problem, config) } +// NSGA-Improved handler functions +fn handle_g04_nsga_improved() -> Result<(), Box> { + let problem = problem_g04(); + let config = NsgaConfig { + population_size: POPULATION, + parents_count: PARENTS_COUNT, + iterations: ITERATIONS, + mutation_std_dev: 0.1, + }; + run_nsga_improved(problem, config) +} + +fn handle_g05_nsga_improved() -> Result<(), Box> { + let problem = problem_g05(); + let config = NsgaConfig { + population_size: POPULATION, + parents_count: PARENTS_COUNT, + iterations: ITERATIONS, + mutation_std_dev: 0.1, + }; + run_nsga_improved(problem, config) +} + +fn handle_g06_nsga_improved() -> Result<(), Box> { + let problem = problem_g06(); + let config = NsgaConfig { + population_size: POPULATION, + parents_count: PARENTS_COUNT, + iterations: ITERATIONS, + mutation_std_dev: 0.1, + }; + run_nsga_improved(problem, config) +} + +fn handle_g08_nsga_improved() -> Result<(), Box> { + let problem = problem_g08(); + let config = NsgaConfig { + population_size: POPULATION, + parents_count: PARENTS_COUNT, + iterations: ITERATIONS, + mutation_std_dev: 0.1, + }; + run_nsga_improved(problem, config) +} + +fn handle_g09_nsga_improved() -> Result<(), Box> { + let problem = problem_g09(); + let config = NsgaConfig { + population_size: POPULATION, + parents_count: PARENTS_COUNT, + iterations: ITERATIONS, + mutation_std_dev: 0.1, + }; + run_nsga_improved(problem, config) +} + +fn handle_g11_nsga_improved() -> Result<(), Box> { + let problem = problem_g11(G11_EPS); + let config = NsgaConfig { + population_size: POPULATION, + parents_count: PARENTS_COUNT, + iterations: ITERATIONS, + mutation_std_dev: 0.1, + }; + run_nsga_improved(problem, config) +} + +fn handle_g21_nsga_improved() -> Result<(), Box> { + let problem = problem_g21(); + let config = NsgaConfig { + population_size: POPULATION, + parents_count: PARENTS_COUNT, + iterations: ITERATIONS, + mutation_std_dev: 0.1, + }; + run_nsga_improved(problem, config) +} + +fn handle_g24_nsga_improved() -> Result<(), Box> { + let problem = problem_g24(); + let config = NsgaConfig { + population_size: POPULATION, + parents_count: PARENTS_COUNT, + iterations: ITERATIONS, + mutation_std_dev: 0.1, + }; + run_nsga_improved(problem, config) +} + // NSGA-Multi handler functions for individual problems fn handle_nsga_multi_g06(capped: bool) -> Result<(), Box> { let mut rng = rand::rng(); @@ -1689,11 +1972,27 @@ fn main() { if args.len() != 3 { eprintln!("Usage: {} ", args[0]); - eprintln!("Methods: srank, nsga, nsga_constr, nsga_multi, nsga_multi_noncapped"); + eprintln!("Methods: srank, nsga, nsga_constr, nsga_improved, nsga_multi, nsga_multi_noncapped"); eprintln!("Problems: g04, g05, g06, g08, g09, g11, g21, g24"); std::process::exit(1); } + // Run the main logic with increased stack size + let result = std::thread::Builder::new() + .stack_size(16 * 1024 * 1024) // 16MB stack + .spawn(move || main_with_args(args)) + .expect("Failed to spawn thread") + .join() + .expect("Thread panicked"); + + if let Err(e) = result { + eprintln!("Error: {}", e); + std::process::exit(1); + } +} + +fn main_with_args(args: Vec) -> Result<(), Box> { + let method = &args[1]; let problem = &args[2]; @@ -1722,6 +2021,14 @@ fn main() { ("nsga_constr", "g11") => handle_g11_nsga_constr(), ("nsga_constr", "g21") => handle_g21_nsga_constr(), ("nsga_constr", "g24") => handle_g24_nsga_constr(), + ("nsga_improved", "g04") => handle_g04_nsga_improved(), + ("nsga_improved", "g05") => handle_g05_nsga_improved(), + ("nsga_improved", "g06") => handle_g06_nsga_improved(), + ("nsga_improved", "g08") => handle_g08_nsga_improved(), + ("nsga_improved", "g09") => handle_g09_nsga_improved(), + ("nsga_improved", "g11") => handle_g11_nsga_improved(), + ("nsga_improved", "g21") => handle_g21_nsga_improved(), + ("nsga_improved", "g24") => handle_g24_nsga_improved(), ("nsga_multi", "g04") => handle_nsga_multi_g04(true), ("nsga_multi", "g05") => handle_nsga_multi_g05(true), ("nsga_multi", "g06") => handle_nsga_multi_g06(true), @@ -1740,7 +2047,7 @@ fn main() { ("nsga_multi_noncapped", "g24") => handle_nsga_multi_g24(false), (_, _) => { eprintln!("Invalid method '{}' or problem '{}'", method, problem); - eprintln!("Methods: srank, nsga, nsga_constr, nsga_multi, nsga_multi_noncapped"); + eprintln!("Methods: srank, nsga, nsga_constr, nsga_improved, nsga_multi, nsga_multi_noncapped"); eprintln!("Problems: g04, g05, g06, g08, g09, g11, g21, g24"); std::process::exit(1); } @@ -1750,4 +2057,6 @@ fn main() { eprintln!("Error: {}", e); std::process::exit(1); } + + Ok(()) } diff --git a/codes/eoa_lib/src/multi_objective_evolution.rs b/codes/eoa_lib/src/multi_objective_evolution.rs index ec8d9a061b6486cc63dce671dabb414dd0a85043..5017331d442a7fcf742fd373a4144a149d4d4e90 100644 --- a/codes/eoa_lib/src/multi_objective_evolution.rs +++ b/codes/eoa_lib/src/multi_objective_evolution.rs @@ -171,7 +171,6 @@ impl<'a, } } - #[derive(Debug, Clone, PartialEq)] pub struct NSGAEvaluation { pub evaluations: [TOut; OBJECTIVES], @@ -558,9 +557,8 @@ pub fn nsga_2>, - &EvaluatedPopulation> - - // TODO + &EvaluatedPopulation>, + &mut TCrossover ), better_than_stats: impl Fn(&TChromosome, &NSGAEvaluation, &Option>>) -> bool, ) -> Result, Box> { @@ -576,8 +574,8 @@ pub fn nsga_2>, - &EvaluatedPopulation> - - // TODO + &EvaluatedPopulation>, + &mut TCrossover ), better_than_stats: impl Fn(&TChromosome, &NSGAEvaluation, &Option>>) -> bool, ) -> Result, Box> { @@ -631,8 +628,8 @@ pub fn constrained_nsga_2< &MinimizingOperator::new(), iterations, rng, - |iteration, stats, population, _, _, _, _, _, _| { - evolutionary_strategy(iteration, stats, population); + |iteration, stats, population, _, _, _, crossover, _, _| { + evolutionary_strategy(iteration, stats, population, crossover); }, better_than_stats )?;