From d2fb123198a882243cfd55408a80c3248ec45522 Mon Sep 17 00:00:00 2001 From: Rutherther Date: Sun, 30 Nov 2025 18:34:44 +0100 Subject: [PATCH] feat: use stochastic ranking for constr_hw02 --- codes/Cargo.lock | 2 + codes/constr_hw02/Cargo.toml | 2 + codes/constr_hw02/src/main.rs | 202 +++++++++++++++++++++++++++++----- 3 files changed, 179 insertions(+), 27 deletions(-) diff --git a/codes/Cargo.lock b/codes/Cargo.lock index 2f137dc1623518db88d3e7f09d0dbfbbcbd75fe1..7728d46e0173857b917580977f60e91ab6a13afb 100644 --- a/codes/Cargo.lock +++ b/codes/Cargo.lock @@ -103,6 +103,8 @@ version = "0.1.0" dependencies = [ "eoa_lib", "nalgebra", + "rand", + "rand_distr", ] [[package]] diff --git a/codes/constr_hw02/Cargo.toml b/codes/constr_hw02/Cargo.toml index f4e96be17e301d63d66473b43ded5a6232d02826..30807d5b707be04c993c0f4180791a15c031dc55 100644 --- a/codes/constr_hw02/Cargo.toml +++ b/codes/constr_hw02/Cargo.toml @@ -6,3 +6,5 @@ edition = "2024" [dependencies] eoa_lib = { version = "0.1.0", path = "../eoa_lib" } nalgebra = "0.33.2" +rand = "0.9.2" +rand_distr = "0.5.1" diff --git a/codes/constr_hw02/src/main.rs b/codes/constr_hw02/src/main.rs index 54558aea2dda2e8e85a4ff9d33e1608b7f765e8a..2e978bb278843a130909806f62d4431ca9a8116b 100644 --- a/codes/constr_hw02/src/main.rs +++ b/codes/constr_hw02/src/main.rs @@ -1,7 +1,11 @@ use std::convert::Infallible; -use eoa_lib::{constraints::LowerThanConstraintFunction, fitness::FitnessFunction}; -use nalgebra::{OVector, SVector}; +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} +}; +use nalgebra::SVector; +use rand::RngCore; +use rand_distr::Normal; pub struct ArbitraryFitness { fun: Box) -> f64> @@ -25,12 +29,20 @@ impl FitnessFunction for ArbitraryFitness { } } -fn problem_g06() -> (ArbitraryFitness<2>, [LowerThanConstraintFunction, f64>; 2], f64) { - ( - ArbitraryFitness::new( +/// A constrained optimization problem with clear field definitions +pub struct ConstrainedProblem { + pub objective: ArbitraryFitness<2>, + pub constraints: [LowerThanConstraintFunction, f64>; CONSTRAINTS], + pub bounds: (SVector, SVector), // (min, max) + pub optimal_value: f64, +} + +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) ), @@ -38,13 +50,17 @@ fn problem_g06() -> (ArbitraryFitness<2>, [LowerThanConstraintFunction::new(0.0, 0.0), // min bounds + SVector::::new(100.0, 100.0) // max bounds + ), + optimal_value: -6961.8137558015, + } } -fn problem_g08(eps: f64) -> (ArbitraryFitness<2>, [LowerThanConstraintFunction, f64>; 2], f64) { - ( - ArbitraryFitness::new( +fn problem_g08() -> ConstrainedProblem<2> { + ConstrainedProblem { + objective: ArbitraryFitness::new( Box::new(|vec| { let num = (2.0 * std::f64::consts::PI * vec[0]).sin().powi(3) * (2.0 * std::f64::consts::PI * vec[1]).sin(); @@ -52,7 +68,7 @@ fn problem_g08(eps: f64) -> (ArbitraryFitness<2>, [LowerThanConstraintFunction (ArbitraryFitness<2>, [LowerThanConstraintFunction::new(0.0, 0.0), // min bounds + SVector::::new(10.0, 10.0) // max bounds + ), + optimal_value: -0.0958250414180359, + } } -pub fn problem_g11(eps: f64) -> (ArbitraryFitness<2>, [LowerThanConstraintFunction, f64>; 2], f64) { - ( - ArbitraryFitness::new( +pub fn problem_g11(eps: f64) -> ConstrainedProblem<2> { + ConstrainedProblem { + objective: ArbitraryFitness::new( Box::new(|vec| { // Minimize f(x) = x1^2 + (x2 - 1)^2 vec[0].powi(2) + (vec[1] - 1.0).powi(2) }) ), - [ + constraints: [ // Equality h(x) = x2 - x1^2 = 0 // Transformed 1: h - eps >= 0 => -(h - eps) <= 0 LowerThanConstraintFunction::new( @@ -97,19 +117,23 @@ pub fn problem_g11(eps: f64) -> (ArbitraryFitness<2>, [LowerThanConstraintFuncti }) ), ], - 0.7499 // Best known optimum - ) + bounds: ( + SVector::::new(-1.0, -1.0), // min bounds + SVector::::new(1.0, 1.0) // max bounds + ), + optimal_value: 0.7499, // Best known optimum + } } -pub fn problem_g24() -> (ArbitraryFitness<2>, [LowerThanConstraintFunction, f64>; 2], f64) { - ( - ArbitraryFitness::new( +pub fn problem_g24() -> ConstrainedProblem<2> { + ConstrainedProblem { + objective: ArbitraryFitness::new( Box::new(|vec| { // Minimize f(x) = -x1 - x2 -vec[0] - vec[1] }) ), - [ + constraints: [ // g1(x) = -2x1^4 + 8x1^3 - 8x1^2 + x2 - 2 <= 0 LowerThanConstraintFunction::new( Box::new(|vec| { @@ -123,10 +147,134 @@ pub fn problem_g24() -> (ArbitraryFitness<2>, [LowerThanConstraintFunction::new(0.0, 0.0), // min bounds + SVector::::new(3.0, 4.0) // max bounds + ), + optimal_value: -5.50801327159536, // Best known optimum + } +} + +/// Solve a constrained optimization problem using stochastic ranking +/// +/// Returns the evolution result with feasible fractions for each iteration +pub fn solve_with_stochastic_ranking( + mut problem: ConstrainedProblem, + 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, EvolutionResult, f64>, 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::<2>, population_size, rng) + ); + + // Setup components as specified + // let mut selection = TournamentSelection::new(5, 0.9); + let mut selection = TournamentSelection::new(5, 0.95); + let mut replacement = GenerationalReplacement; + let mut pairing = AdjacentPairing::new(); + let mut crossover = BoundedCrossover::, 2, _>::new( + ArithmeticCrossover::new(), + problem.bounds.0, + problem.bounds.1, + BoundedCrossoverStrategy::Trim + ); + + // Setup bounded random distribution perturbation with Normal distribution + let normal_perturbation = RandomDistributionPerturbation::<2, Normal>::normal(mutation_std_dev)?; + let perturbation = BoundedPerturbation::new( + normal_perturbation, + problem.bounds.0, + problem.bounds.1, + BoundedPerturbationStrategy::Trim + ); + let mut mutation = MutationPerturbation::new(Box::new(perturbation), 0.1); + + // Constraint weights (initially set to 1.0) + let constraint_weights = [1.0; CONSTRAINTS]; + + let (N, p) = stochastic_params; + let better_than = MinimizingOperator::new(); + + // Convert constraint array references + let constraint_refs = problem.constraints.iter().collect::>().try_into() + .map_err(|_| "Failed to convert constraint references")?; + + stochastic_ranking_evolution_algorithm( + initial_population, + parents_count, + N, + p, + &mut problem.objective, + constraint_refs, + constraint_weights, + &mut pairing, + &mut selection, + &mut crossover, + &mut mutation, + &mut replacement, + &better_than, + iterations, + rng, + ).map(|(x, y)| (problem, x, y)) } fn main() { - println!("Hello, world!"); + let mut rng = rand::rng(); + + // Example usage with problem G06 + let problem = problem_g11(0.01); + + 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, + ); + + 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); + } + } + } + } else { + println!("Could not find any feasible solution!") + } + } + Err(e) => { + println!("Error: {}", e); + } + } }