Into the Swarm: Enhancing Decision-Making with Particle Swarm Optimization

Optimization is a branch of mathematics concerned with finding the best solution to a problem with questions like: How do I maximize profit, how do I minimize cost, etc. They’re important questions, but sometimes it’s hard to determine the exact best solution. We can turn to heuristic algorithms, which are procedures for more quickly finding a sufficiently good solution. Particle Swarm Optimization is one such procedure.

Let’s start with a problem. An alien race wants to find the best place to put a windmill in a valley to maximize energy production. To assist them, they have a large number (a swarm, even) of simple robots that each have a wind-speed gauge, the ability to remember the position of the highest wind speed they’ve found so far, and a radio to talk to the other robots. They could naively try checking every spot, but it would take too long and they’re very busy avoiding hobbyist telescopes. We’ll try using Particle Swarm instead.

On to the procedure… Because we don’t know where the best windmill placement would be, we’ll start by randomly placing the robots around the valley and let them begin exploring in a random direction, measuring the wind speed as they go, and remembering the best position they’ve each found so far. At each step, we’ll nudge them a little bit towards the best solution each has found individually, hypothesizing that areas of higher wind will probably be close to each other. To allow the robots to collaborate, we’ll also nudge them a little bit towards the best solution they’ve found collectively.

(In this example, purple has found the best spot so far, so blue and green are also nudged toward purple’s best in addition to their own)

Let’s see how they do.

Looks like the robots all worked together and agreed on a best spot and our aliens are happy. Hooray!

Let’s try that example again, this time with some code. We’ll start with a particle and give it some extra parameters to control its behavior. Let’s also give it a random initial position and velocity, and a place to track the best position and score it’s found so far.

import numpy as np 

class Particle: 
  def __init__( 
    self, 
    dimensions, 
    independent_acceleration=0.3, 
    collaborative_acceleration=0.5, 
    inertia=0.9, 
  ): 
    self.dimensions = dimensions 
    self.independent_acceleration = independent_acceleration 
    self.collaborative_acceleration = collaborative_acceleration 
    self.inertia = inertia 

    # random initial n-dimensional position, arbitrarily placed between -5 and +5
    self.position = np.random.rand(self.dimensions) * 10 - 5 
    # random initial n-dimensional  velocity, arbitrarily placed between -5 and +5
    self.velocity = np.random.rand(self.dimensions) * 10 - 5 
    self.local_best_position = np.copy(self.position)
    # start off with a low score, hopefully to be improved over time
    self.local_best_score = float("-inf")

It will also need a way to move through space:

  def _update_position(self): 
    # move position in direction of current velocity
    self.position += self.velocity Code language: PHP (php)

… a way to keep updating its personal locally-best found position, based on an objective function to maximize:

  def _update_local_best(self, objective_function): 
    # assess score for location
    current_score = objective_function(self.position) 
    # if score is best found so far,
    if current_score > self.local_best_score: 
      # set best location and score to current value
      self.local_best_score = current_score 
      self.local_best_position = np.copy(self.position) Code language: PHP (php)

… and something to nudge it towards both the global and local best solutions:

    def _update_velocity(self, global_best_position):
        # Nudge particle towards the best position it's found individually.
        individual_component = (
            self.independent_acceleration
            * np.random.rand(self.dimensions)  # some randomness to encourage exploration
            * (self.local_best_position - self.position)
        )

        # Also nudge particle towards the best position ANY particle has found.
        social_component = (
            self.collaborative_acceleration
            * np.random.rand(self.dimensions)  # some more randomness
            * (global_best_position - self.position)
        )

        # Update velocity according to nudges, weighted by inertia factor
        self.velocity = (
             self.velocity + individual_component + social_component
        ) * self.inertia 
        # not exactly how inertia works, but it serves our purpose.Code language: PHP (php)

For convenience, we’ll wrap those together into a single method.

    def update(self, global_best_position, objective_function):
        self._update_velocity(global_best_position)
        self._update_position()
        self._update_local_best(objective_function)Code language: PHP (php)

Next, let’s add something to keep track of a collection of these particles, the global best position any of them have found, and a loop to iteratively allow the particles to move:

class ParticleSwarmOptimizer: 
  def __init__(self, num_particles, dimensions, iterations): 
    self.iterations = iterations 
    # random initial n-dimensional position, arbitrarily placed between -5 and +5
    self.global_best_position = np.random.rand(dimensions) * 10 - 5 
    # start off with a low score, hopefully to be improved over time
    self.global_best_score = float("-inf") 
    self.particles = [Particle(dimensions) for _ in range(num_particles)] 
 
  def optimize(self, objective_function):
    # limit number of iterations 
    for _ in range(self.iterations): 
      # loop through all particles in swarm
      for particle in self.particles: 
        # move particle, check score
        particle.update(self.global_best_position, objective_function) 
        # compare individual score to global score
        if particle.local_best_score > self.global_best_score: 
          self.global_best_score = particle.local_best_score 
          self.global_best_position = np.copy(particle.local_best_position) 

    return self.global_best_position, self.global_best_score

And voila, we have a simple optimizer that allows our particles to share global information and collaboratively move towards good solutions.

To test it, we’ll make a simple objective function. Here, we’ll just try to find the location that minimizes the distance between a given particle and the three-dimensional point (x=1, y=2, z=3). Ideally, it should find something that looks very close to [1 2 3].

def distance_to_best_location(particle): 
  max_wind_location = np.array([1, 2, 3]) 
  # euclidean distance between particle and our target
  distance = np.linalg.norm(particle - max_wind_location) 
  # multiply by -1 so higher distances return lower scores
  return -1 * distance Code language: PHP (php)

Passing this simple function into our optimizer…

dimensions = 3
num_particles = 3
iterations = 100

pso = ParticleSwarmOptimizer(num_particles, dimensions, iterations)
best_position, best_score = pso.optimize(distance_to_best_location)

print("Best Position:", best_position)
print("Best Score:", best_score)Code language: PHP (php)

Pretty close!

Best Position: [0.99803055 1.99951829 3.00216642]
Best Score: -0.0029671808924744227Code language: CSS (css)

Let’s give it a slightly harder problem, which also happens to have a best answer of [1 2 3].

import math
def harder_objective_function(particle):
    x, y, z = particle
    x_component = math.sin(x / 2 * math.pi) * 12 + (-(x - 1) ** 2) – 12
    y_component = math.cos(y - 2) + (- (y -2) ** 2) - 1
    z_component = (-abs(z - 3) + 1) - 1
    score = x_component + y_component + z_component
    return scoreCode language: JavaScript (javascript)

Best Position: [0.999845   1.99783829 2.99998323]
Best Score: -2.4160211942847454e-05Code language: CSS (css)

As a final test, we’ll give it a totally real and relatable business problem to solve.

Problem statement:
Your company manufactures two products: Twirlyfloppers and Bizzlywigs. 
Each Twirlyflopper sells for $80 and requires the following resources for production: 
    1 hour of quantum computing, 
    1 hour of flux dephasing, 
    1 kg of space dust. 
    
Alternatively, each Bizzlywig sells for $110 and requires the following:
    2 hours of flux dephasing, 
    1 hour of quantum computing, 
    900 grams of space dust.

Your company's weekly resource availability and costs are as follows:

Quantum computing: 40 hours available at a cost of $10 per hour.
Flux dephasing: 50 hours available at a cost of $8 per hour.
Space dust: Unlimited availability at a cost of $4 per kg.

How many Twirlyfloppers and Bizzlywigs should you make this week?

Let’s define these parameters and constraints in code:

## Production Requirements:
# Twirlyflopper
FLUX_DEPHASING_HOURS_PER_TWIRLYFLOPPER = 1
QUANTUM_HOURS_PER_TWIRLYFLOPPER = 1
SPACE_DUST_KG_PER_TWIRLYFLOPPER = 1
# Bizzlywig
FLUX_DEPHASING_HOURS_PER_BIZZLYWIG = 2
QUANTUM_HOURS_PER_BIZZLYWIG = 1
SPACE_DUST_KG_PER_BIZZLYWIG = 0.9

## Constraints and costs
QUANTUM_HOURS_MAX = 40
QUANTUM_HOUR_COST = 10
FLUX_DEPHASING_HOURS_MAX = 50
FLUX_DEPHASING_HOUR_COST = 8
SPACE_DUST_COST_PER_KG = 4

## Revenue
SELL_PRICE_PER_TWIRLYFLOPPER = 80
SELL_PRICE_PER_BIZZLYWIG = 110
Code language: PHP (php)

And we’ll formulate a function for what we want to maximize, in this case profit:

def profit(particle):
    twirlyfloppers, bizzlywigs = particle

    quantum_hours_used = (
        twirlyfloppers * QUANTUM_HOURS_PER_TWIRLYFLOPPER
        + bizzlywigs * QUANTUM_HOURS_PER_BIZZLYWIG
    )
    flux_dephasing_hours_used = (
        twirlyfloppers * FLUX_DEPHASING_HOURS_PER_TWIRLYFLOPPER
        + bizzlywigs * FLUX_DEPHASING_HOURS_PER_BIZZLYWIG
    )
    space_dust_used = (
        twirlyfloppers * SPACE_DUST_KG_PER_TWIRLYFLOPPER
        + bizzlywigs * SPACE_DUST_KG_PER_BIZZLYWIG
    )

    cost = (
        quantum_hours_used * QUANTUM_HOUR_COST
        + flux_dephasing_hours_used * FLUX_DEPHASING_HOUR_COST
        + space_dust_used * SPACE_DUST_COST_PER_KG
    )

    # penalty for exceeding capacities
    penalty = 0
    if quantum_hours_used > QUANTUM_HOURS_MAX:
        penalty += float('inf')

    if flux_dephasing_hours_used > FLUX_DEPHASING_HOURS_MAX:
        penalty += float('inf')

    revenue = (
        twirlyfloppers * SELL_PRICE_PER_TWIRLYFLOPPER
        + bizzlywigs * SELL_PRICE_PER_BIZZLYWIG
    )

    profit = revenue - cost - penalty
    return profitCode language: PHP (php)

And now to plug that into our optimizer. This time we only have two dimensions to represent the numbers of twirlyfloppers and bizzlywigs, respectively.

dimensions = 2
num_particles = 100
iterations = 100

pso = ParticleSwarmOptimizer(num_particles, dimensions, iterations)
(twirlyfloppers, bizzlywigs), best_profit = pso.optimize(profit)
print("Best Production:", f"{twirlyfloppers=}", f"{bizzlywigs=}")
print("Best profit:", best_profit)Code language: PHP (php)

And for our answer..

Best Production: twirlyfloppers=30.010844922626674 bizzlywigs=9.987382774722978
Best profit: 2543.614580600075

Looks like around 30 twirlyfloppers and 10 bizzlywigs is our optimum, which a quick graphical analysis supports:

We’ve done it! All our business problems have been solved with just a few short lines of code.As a bonus, here’s what Dall-E imagines our products look like:

In a closing note… as is common with Python, there’s an existing package for that. I’d recommend pyswarms, but note the objective function expects the entire array of particles, rather than one at a time, and attempts to minimize a cost rather than maximize a score.

import pyswarms as ps

options = {"c1": 0.5, "c2": 0.3, "w": 0.9}

def distances_to_best_location(particles):
    max_wind_location = np.array([1, 2, 3])
    distances = np.sqrt(np.sum((particles - max_wind_location) ** 2, axis=1))
    return distances

optimizer = ps.single.GlobalBestPSO(n_particles=100, dimensions=3, options=options)
cost, pos = optimizer.optimize(distances_to_best_location, iters=100)
print("Best Position:", pos)
print("Best Cost:", cost)Code language: PHP (php)

Not bad, but ours was more fun.

Best Position: [1.00017438 1.9999488  2.99994409]
Best Cost: 0.00019014543705898434Code language: CSS (css)

Due to its flexibility, Particle Swarm Optimization is applicable to other domains as well, such as power system regulation, engineering design, machine learning parameter tuning, and a nearly endless list of other disciplines. If you’d like to see how else data science and software engineering can improve your life, continue your journey with our other blogs.