scirs2_optimize/multi_objective/algorithms/
mod.rs

1//! Multi-objective optimization algorithms
2//!
3//! This module provides implementations of various multi-objective optimization algorithms:
4//! - NSGA-II (Non-dominated Sorting Genetic Algorithm II)
5//! - NSGA-III (Non-dominated Sorting Genetic Algorithm III)
6//! - MOEA/D (Multi-Objective Evolutionary Algorithm based on Decomposition)
7//! - SPEA2 (Strength Pareto Evolutionary Algorithm 2)
8
9pub mod moead;
10pub mod nsga2;
11pub mod nsga3;
12pub mod spea2;
13
14pub use moead::MOEAD;
15pub use nsga2::NSGAII;
16pub use nsga3::NSGAIII;
17pub use spea2::SPEA2;
18
19use super::solutions::{MultiObjectiveResult, MultiObjectiveSolution, Population};
20use crate::error::OptimizeError;
21use ndarray::{Array1, ArrayView1};
22
23/// Configuration for multi-objective optimization algorithms
24#[derive(Debug, Clone)]
25pub struct MultiObjectiveConfig {
26    /// Population size
27    pub population_size: usize,
28    /// Maximum number of generations
29    pub max_generations: usize,
30    /// Maximum number of function evaluations
31    pub max_evaluations: Option<usize>,
32    /// Crossover probability
33    pub crossover_probability: f64,
34    /// Mutation probability
35    pub mutation_probability: f64,
36    /// Mutation strength (eta for polynomial mutation)
37    pub mutation_eta: f64,
38    /// Crossover distribution index (eta for SBX)
39    pub crossover_eta: f64,
40    /// Convergence tolerance for hypervolume
41    pub tolerance: f64,
42    /// Reference point for hypervolume calculation
43    pub reference_point: Option<Array1<f64>>,
44    /// Variable bounds (lower, upper)
45    pub bounds: Option<(Array1<f64>, Array1<f64>)>,
46    /// Random seed for reproducibility
47    pub random_seed: Option<u64>,
48    /// Archive size (for algorithms that use archives)
49    pub archive_size: Option<usize>,
50    /// Neighborhood size (for MOEA/D)
51    pub neighborhood_size: Option<usize>,
52}
53
54impl Default for MultiObjectiveConfig {
55    fn default() -> Self {
56        Self {
57            population_size: 100,
58            max_generations: 250,
59            max_evaluations: None,
60            crossover_probability: 0.9,
61            mutation_probability: 0.1,
62            mutation_eta: 20.0,
63            crossover_eta: 15.0,
64            tolerance: 1e-6,
65            reference_point: None,
66            bounds: None,
67            random_seed: None,
68            archive_size: None,
69            neighborhood_size: None,
70        }
71    }
72}
73
74/// Trait for multi-objective optimization algorithms
75pub trait MultiObjectiveOptimizer {
76    /// Optimize the problem
77    fn optimize<F>(&mut self, objective_function: F) -> Result<MultiObjectiveResult, OptimizeError>
78    where
79        F: Fn(&ArrayView1<f64>) -> Array1<f64> + Send + Sync;
80
81    /// Initialize the population
82    fn initialize_population(&mut self) -> Result<(), OptimizeError>;
83
84    /// Perform one generation/iteration
85    fn evolve_generation<F>(&mut self, objective_function: &F) -> Result<(), OptimizeError>
86    where
87        F: Fn(&ArrayView1<f64>) -> Array1<f64> + Send + Sync;
88
89    /// Check convergence criteria
90    fn check_convergence(&self) -> bool;
91
92    /// Get current population
93    fn get_population(&self) -> &Population;
94
95    /// Get current generation number
96    fn get_generation(&self) -> usize;
97
98    /// Get number of function evaluations
99    fn get_evaluations(&self) -> usize;
100
101    /// Get algorithm name
102    fn name(&self) -> &str;
103}
104
105/// Wrapper enum for different multi-objective optimizers
106pub enum MultiObjectiveOptimizerWrapper {
107    NSGAII(NSGAII),
108    NSGAIII(NSGAIII),
109    MOEAD(MOEAD),
110    SPEA2(SPEA2),
111}
112
113impl MultiObjectiveOptimizerWrapper {
114    pub fn optimize<F>(
115        &mut self,
116        objective_function: F,
117    ) -> Result<MultiObjectiveResult, OptimizeError>
118    where
119        F: Fn(&ArrayView1<f64>) -> Array1<f64> + Send + Sync,
120    {
121        match self {
122            Self::NSGAII(opt) => opt.optimize(objective_function),
123            Self::NSGAIII(opt) => opt.optimize(objective_function),
124            Self::MOEAD(opt) => opt.optimize(objective_function),
125            Self::SPEA2(opt) => opt.optimize(objective_function),
126        }
127    }
128}
129
130/// Factory for creating multi-objective optimizers
131pub struct OptimizerFactory;
132
133impl OptimizerFactory {
134    /// Create NSGA-II optimizer
135    pub fn create_nsga2(
136        config: MultiObjectiveConfig,
137        n_objectives: usize,
138        n_variables: usize,
139    ) -> Result<NSGAII, OptimizeError> {
140        NSGAII::new(config, n_objectives, n_variables)
141    }
142
143    /// Create NSGA-III optimizer
144    pub fn create_nsga3(
145        config: MultiObjectiveConfig,
146        n_objectives: usize,
147        n_variables: usize,
148        reference_points: Option<Vec<Array1<f64>>>,
149    ) -> Result<NSGAIII, OptimizeError> {
150        // TODO: Use reference_points when NSGA-III is fully implemented
151        Ok(NSGAIII::new(
152            config.population_size,
153            n_objectives,
154            n_variables,
155        ))
156    }
157
158    /// Create optimizer by name
159    pub fn create_by_name(
160        algorithm: &str,
161        config: MultiObjectiveConfig,
162        n_objectives: usize,
163        n_variables: usize,
164    ) -> Result<MultiObjectiveOptimizerWrapper, OptimizeError> {
165        match algorithm.to_lowercase().as_str() {
166            "nsga2" | "nsga-ii" => Ok(MultiObjectiveOptimizerWrapper::NSGAII(Self::create_nsga2(
167                config,
168                n_objectives,
169                n_variables,
170            )?)),
171            "nsga3" | "nsga-iii" => Ok(MultiObjectiveOptimizerWrapper::NSGAIII(
172                Self::create_nsga3(config, n_objectives, n_variables, None)?,
173            )),
174            _ => Err(OptimizeError::InvalidInput(format!(
175                "Unknown algorithm: {}",
176                algorithm
177            ))),
178        }
179    }
180}
181
182/// Utility functions for multi-objective optimization
183pub mod utils {
184    use super::*;
185    use ndarray::Array2;
186
187    /// Generate Das-Dennis reference points for NSGA-III
188    pub fn generate_das_dennis_points(
189        n_objectives: usize,
190        n_partitions: usize,
191    ) -> Vec<Array1<f64>> {
192        if n_objectives == 1 {
193            return vec![Array1::from_vec(vec![1.0])];
194        }
195
196        let mut points = Vec::new();
197        generate_das_dennis_recursive(
198            &mut points,
199            Array1::zeros(n_objectives),
200            0,
201            n_objectives,
202            n_partitions,
203            n_partitions,
204        );
205
206        // Normalize points
207        for point in &mut points {
208            let sum: f64 = point.sum();
209            if sum > 0.0 {
210                *point /= sum;
211            }
212        }
213
214        points
215    }
216
217    fn generate_das_dennis_recursive(
218        points: &mut Vec<Array1<f64>>,
219        mut current_point: Array1<f64>,
220        index: usize,
221        n_objectives: usize,
222        n_partitions: usize,
223        remaining: usize,
224    ) {
225        if index == n_objectives - 1 {
226            current_point[index] = remaining as f64;
227            points.push(current_point);
228        } else {
229            for i in 0..=remaining {
230                current_point[index] = i as f64;
231                generate_das_dennis_recursive(
232                    points,
233                    current_point.clone(),
234                    index + 1,
235                    n_objectives,
236                    n_partitions,
237                    remaining - i,
238                );
239            }
240        }
241    }
242
243    /// Calculate hypervolume indicator
244    pub fn calculate_hypervolume(
245        pareto_front: &[MultiObjectiveSolution],
246        reference_point: &Array1<f64>,
247    ) -> f64 {
248        if pareto_front.is_empty() {
249            return 0.0;
250        }
251
252        // For 2D case, use simple calculation
253        if reference_point.len() == 2 {
254            return calculate_hypervolume_2d(pareto_front, reference_point);
255        }
256
257        // For higher dimensions, use Monte Carlo approximation
258        calculate_hypervolume_monte_carlo(pareto_front, reference_point, 10000)
259    }
260
261    fn calculate_hypervolume_2d(
262        pareto_front: &[MultiObjectiveSolution],
263        reference_point: &Array1<f64>,
264    ) -> f64 {
265        let mut points: Vec<_> = pareto_front
266            .iter()
267            .map(|sol| (sol.objectives[0], sol.objectives[1]))
268            .collect();
269
270        // Sort by first objective in descending order for correct hypervolume calculation
271        points.sort_by(|a, b| b.0.partial_cmp(&a.0).unwrap());
272
273        let mut hypervolume = 0.0;
274        let mut prev_x = reference_point[0];
275
276        for (x, y) in points {
277            if x < reference_point[0] && y < reference_point[1] {
278                hypervolume += (prev_x - x) * (reference_point[1] - y);
279                prev_x = x;
280            }
281        }
282
283        hypervolume
284    }
285
286    fn calculate_hypervolume_monte_carlo(
287        pareto_front: &[MultiObjectiveSolution],
288        reference_point: &Array1<f64>,
289        n_samples: usize,
290    ) -> f64 {
291        use rand::prelude::*;
292        let mut rng = rand::rng();
293        let n_objectives = reference_point.len();
294
295        // Find bounds for sampling
296        let mut min_bounds = Array1::from_elem(n_objectives, f64::INFINITY);
297        for sol in pareto_front {
298            for (i, &obj) in sol.objectives.iter().enumerate() {
299                min_bounds[i] = min_bounds[i].min(obj);
300            }
301        }
302
303        let mut dominated_count = 0;
304        for _ in 0..n_samples {
305            // Generate random point
306            let mut point = Array1::zeros(n_objectives);
307            for i in 0..n_objectives {
308                point[i] = rng.gen_range(min_bounds[i]..reference_point[i]);
309            }
310
311            // Check if point is dominated by any solution in Pareto front
312            for sol in pareto_front {
313                let mut dominates = true;
314                for (i, &obj) in sol.objectives.iter().enumerate() {
315                    if obj >= point[i] {
316                        dominates = false;
317                        break;
318                    }
319                }
320                if dominates {
321                    dominated_count += 1;
322                    break;
323                }
324            }
325        }
326
327        // Calculate volume
328        let total_volume: f64 = (0..n_objectives)
329            .map(|i| reference_point[i] - min_bounds[i])
330            .product();
331
332        total_volume * (dominated_count as f64 / n_samples as f64)
333    }
334
335    /// Calculate diversity metric (spacing)
336    pub fn calculate_spacing(pareto_front: &[MultiObjectiveSolution]) -> f64 {
337        if pareto_front.len() < 2 {
338            return 0.0;
339        }
340
341        let distances: Vec<f64> = pareto_front
342            .iter()
343            .map(|sol1| {
344                pareto_front
345                    .iter()
346                    .filter(|sol2| sol1 as *const _ != *sol2 as *const _)
347                    .map(|sol2| sol1.objective_distance(sol2))
348                    .fold(f64::INFINITY, f64::min)
349            })
350            .collect();
351
352        let mean_distance = distances.iter().sum::<f64>() / distances.len() as f64;
353        let variance = distances
354            .iter()
355            .map(|d| (d - mean_distance).powi(2))
356            .sum::<f64>()
357            / distances.len() as f64;
358
359        variance.sqrt()
360    }
361
362    /// Calculate convergence metric (average distance to true Pareto front)
363    pub fn calculate_convergence(
364        pareto_front: &[MultiObjectiveSolution],
365        true_pareto_front: &[MultiObjectiveSolution],
366    ) -> f64 {
367        if pareto_front.is_empty() || true_pareto_front.is_empty() {
368            return f64::INFINITY;
369        }
370
371        let total_distance: f64 = pareto_front
372            .iter()
373            .map(|sol1| {
374                true_pareto_front
375                    .iter()
376                    .map(|sol2| sol1.objective_distance(sol2))
377                    .fold(f64::INFINITY, f64::min)
378            })
379            .sum();
380
381        total_distance / pareto_front.len() as f64
382    }
383
384    /// Generate random initial population
385    pub fn generate_random_population(
386        size: usize,
387        n_variables: usize,
388        bounds: &Option<(Array1<f64>, Array1<f64>)>,
389    ) -> Vec<Array1<f64>> {
390        use rand::prelude::*;
391        let mut rng = rand::rng();
392        let mut population = Vec::with_capacity(size);
393
394        let (lower, upper) = match bounds {
395            Some((l, u)) => (l.clone(), u.clone()),
396            None => (Array1::zeros(n_variables), Array1::ones(n_variables)),
397        };
398
399        for _ in 0..size {
400            let mut individual = Array1::zeros(n_variables);
401            for j in 0..n_variables {
402                individual[j] = rng.gen_range(lower[j]..upper[j]);
403            }
404            population.push(individual);
405        }
406
407        population
408    }
409
410    /// Apply bounds to a solution
411    pub fn apply_bounds(individual: &mut Array1<f64>, bounds: &Option<(Array1<f64>, Array1<f64>)>) {
412        if let Some((lower, upper)) = bounds {
413            for (i, value) in individual.iter_mut().enumerate() {
414                *value = value.max(lower[i]).min(upper[i]);
415            }
416        }
417    }
418}
419
420#[cfg(test)]
421mod tests {
422    use super::*;
423    use ndarray::array;
424
425    #[test]
426    fn test_config_default() {
427        let config = MultiObjectiveConfig::default();
428        assert_eq!(config.population_size, 100);
429        assert_eq!(config.max_generations, 250);
430        assert_eq!(config.crossover_probability, 0.9);
431    }
432
433    #[test]
434    fn test_das_dennis_points_2d() {
435        let points = utils::generate_das_dennis_points(2, 3);
436        assert!(!points.is_empty());
437
438        // Check that points are normalized
439        for point in &points {
440            let sum: f64 = point.sum();
441            assert!((sum - 1.0).abs() < 1e-10);
442        }
443    }
444
445    #[test]
446    fn test_hypervolume_2d() {
447        let solutions = vec![
448            MultiObjectiveSolution::new(array![1.0], array![1.0, 3.0]),
449            MultiObjectiveSolution::new(array![2.0], array![2.0, 2.0]),
450            MultiObjectiveSolution::new(array![3.0], array![3.0, 1.0]),
451        ];
452
453        let reference_point = array![4.0, 4.0];
454        let hv = utils::calculate_hypervolume(&solutions, &reference_point);
455        assert!(hv > 0.0);
456    }
457
458    #[test]
459    fn test_spacing_calculation() {
460        let solutions = vec![
461            MultiObjectiveSolution::new(array![1.0], array![1.0, 3.0]),
462            MultiObjectiveSolution::new(array![2.0], array![2.0, 2.0]),
463            MultiObjectiveSolution::new(array![3.0], array![3.0, 1.0]),
464        ];
465
466        let spacing = utils::calculate_spacing(&solutions);
467        assert!(spacing >= 0.0);
468    }
469
470    #[test]
471    fn test_random_population_generation() {
472        let bounds = Some((array![0.0, -1.0], array![1.0, 1.0]));
473        let population = utils::generate_random_population(10, 2, &bounds);
474
475        assert_eq!(population.len(), 10);
476        assert_eq!(population[0].len(), 2);
477
478        // Check bounds
479        for individual in &population {
480            assert!(individual[0] >= 0.0 && individual[0] <= 1.0);
481            assert!(individual[1] >= -1.0 && individual[1] <= 1.0);
482        }
483    }
484
485    #[test]
486    fn test_apply_bounds() {
487        let bounds = Some((array![0.0, -1.0], array![1.0, 1.0]));
488        let mut individual = array![-0.5, 2.0];
489
490        utils::apply_bounds(&mut individual, &bounds);
491
492        assert_eq!(individual[0], 0.0); // Clamped to lower bound
493        assert_eq!(individual[1], 1.0); // Clamped to upper bound
494    }
495
496    #[test]
497    fn test_optimizer_factory() {
498        let config = MultiObjectiveConfig::default();
499
500        let nsga2 = OptimizerFactory::create_by_name("nsga2", config.clone(), 2, 3);
501        assert!(nsga2.is_ok());
502
503        let nsga3 = OptimizerFactory::create_by_name("nsga3", config.clone(), 2, 3);
504        assert!(nsga3.is_ok());
505
506        let unknown = OptimizerFactory::create_by_name("unknown", config, 2, 3);
507        assert!(unknown.is_err());
508    }
509}