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.0029671808924744227
```

Code 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 score
```

Code language: JavaScript (javascript)

```
Best Position: [0.999845 1.99783829 2.99998323]
Best Score: -2.4160211942847454e-05
```

Code 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 profit
```

Code 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.00019014543705898434
```

Code 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.