forge_core/algo/
pso.rs

1//! Particle Swarm Optimization (Kennedy & Eberhart 1995; Clerc & Kennedy 2002).
2//!
3//! A swarm of particles flies through the search space, each pulled toward its
4//! own best-seen position (*cognitive*) and the swarm's best-seen position
5//! (*social*), with momentum (*inertia*). The default coefficients are the
6//! constriction-equivalent inertia weights `w = 0.7298`, `c1 = c2 = 1.49618`
7//! (Clerc & Kennedy 2002), which keep the swarm convergent without explicit
8//! velocity damping.
9//!
10//! The update is *asynchronous*: each particle reads the current global best at
11//! the moment it moves, so improvements propagate within an iteration. Velocity
12//! is clamped to one full variable range per step, and a particle that hits a
13//! bound has its position clamped and that velocity component zeroed, so it does
14//! not stick to the wall at speed.
15
16use super::{clamp, sample, Evaluator, Optimizer};
17use crate::problem::Problem;
18use crate::rng::Rng;
19use crate::solution::{Report, Solution};
20use crate::termination::Termination;
21
22/// PSO configuration.
23#[derive(Debug, Clone, Copy)]
24pub struct Pso {
25    /// Number of particles (≥ 2).
26    pub swarm_size: usize,
27    /// Inertia weight `w`: momentum carried from the previous velocity.
28    pub inertia: f64,
29    /// Cognitive coefficient `c1`: pull toward the particle's personal best.
30    pub cognitive: f64,
31    /// Social coefficient `c2`: pull toward the swarm's global best.
32    pub social: f64,
33    /// RNG seed; same seed + same problem + same budget ⇒ same result.
34    pub seed: u64,
35}
36
37impl Default for Pso {
38    fn default() -> Self {
39        Pso {
40            swarm_size: 30,
41            inertia: 0.7298,
42            cognitive: 1.49618,
43            social: 1.49618,
44            seed: 42,
45        }
46    }
47}
48
49impl Optimizer for Pso {
50    fn with_seed(&self, seed: u64) -> Self {
51        Pso { seed, ..*self }
52    }
53
54    /// Minimizes `problem` within its bounds using PSO.
55    fn optimize(&self, problem: &dyn Problem, term: &Termination) -> Report {
56        crate::problem::validate(problem).unwrap_or_else(|e| panic!("Pso: invalid problem: {e}"));
57        let bounds = problem.bounds();
58        let dim = bounds.len();
59        let n = self.swarm_size.max(2);
60        let mut rng = Rng::new(self.seed);
61
62        // Per-dimension velocity clamp: one full range per step.
63        let vmax: Vec<f64> = bounds.iter().map(|&(lo, hi)| hi - lo).collect();
64
65        // Initialize positions and velocities, seeding the evaluator with the
66        // first particle.
67        let first_x = sample(bounds, &mut rng);
68        let first_v = finite_or_worst(problem.objective(&first_x));
69        let mut ev = Evaluator::new(
70            problem,
71            term,
72            Solution {
73                x: first_x.clone(),
74                value: first_v,
75            },
76        );
77
78        let mut pos: Vec<Vec<f64>> = Vec::with_capacity(n);
79        let mut vel: Vec<Vec<f64>> = Vec::with_capacity(n);
80        let mut pbest_x: Vec<Vec<f64>> = Vec::with_capacity(n);
81        let mut pbest_f: Vec<f64> = Vec::with_capacity(n);
82
83        // First particle (already evaluated).
84        vel.push(init_velocity(&first_x, bounds, &mut rng));
85        pbest_x.push(first_x.clone());
86        pbest_f.push(first_v);
87        pos.push(first_x);
88
89        for _ in 1..n {
90            let x = sample(bounds, &mut rng);
91            let v = init_velocity(&x, bounds, &mut rng);
92            let f = if ev.done() {
93                f64::INFINITY
94            } else {
95                finite_or_worst(ev.eval(&x))
96            };
97            pbest_x.push(x.clone());
98            pbest_f.push(f);
99            pos.push(x);
100            vel.push(v);
101        }
102
103        while !ev.done() {
104            for i in 0..n {
105                if ev.done() {
106                    break;
107                }
108                // Asynchronous: read the current global best for this move.
109                let gbest = ev.best.x.clone();
110                for d in 0..dim {
111                    let (lo, hi) = bounds[d];
112                    let r1 = rng.uniform();
113                    let r2 = rng.uniform();
114                    let mut v = self.inertia * vel[i][d]
115                        + self.cognitive * r1 * (pbest_x[i][d] - pos[i][d])
116                        + self.social * r2 * (gbest[d] - pos[i][d]);
117                    v = clamp(v, -vmax[d], vmax[d]);
118
119                    let mut x = pos[i][d] + v;
120                    if x < lo {
121                        x = lo;
122                        v = 0.0;
123                    } else if x > hi {
124                        x = hi;
125                        v = 0.0;
126                    }
127                    vel[i][d] = v;
128                    pos[i][d] = x;
129                }
130
131                let f = finite_or_worst(ev.eval(&pos[i]));
132                if f < pbest_f[i] {
133                    pbest_f[i] = f;
134                    pbest_x[i].copy_from_slice(&pos[i]);
135                }
136            }
137        }
138
139        ev.finish()
140    }
141}
142
143#[inline]
144fn finite_or_worst(v: f64) -> f64 {
145    if v.is_finite() {
146        v
147    } else {
148        f64::INFINITY
149    }
150}
151
152/// Initial velocity in the SPSO style: half the gap to a fresh random point,
153/// `v[d] = ½(U(lo, hi) − x[d])`. Keeps early motion on the scale of the box.
154fn init_velocity(x: &[f64], bounds: &[(f64, f64)], rng: &mut Rng) -> Vec<f64> {
155    bounds
156        .iter()
157        .zip(x)
158        .map(|(&(lo, hi), &xi)| 0.5 * (rng.uniform_in(lo, hi) - xi))
159        .collect()
160}
forge_core/algo/pso.rs

forge_core/algo/
pso.rs
