diff --git a/chapters/10.Metaheuristics/PSO.jmd b/chapters/10.Metaheuristics/PSO.jmd
new file mode 100644
index 0000000..56474be
--- /dev/null
+++ b/chapters/10.Metaheuristics/PSO.jmd
@@ -0,0 +1,122 @@
+---
+title : Particle Swarm Optimization
+author : Michiel Stock
+date: 2021-2020
+---
+
+```julia; echo=false
+using Plots
+```
+
+# Description
+
+Particle Swarm Optimization (PSO) is a stochastic metaheuristic to solve non-convex continuous optimization problems. It based on the flocking behavior of birds and insects. PSO holds a list of 'particles', containing a position and a velocity in the design space. In every step, the velocity of each particle is updated according to
+- towards the particle's personal best design point;
+- towards the best design point found over all the groups.
+The particles' velocities are subsequently used to update their positions.
+
+Specifically, the position $\mathbf{x}^{(i)}$ and velocity $\mathbf{v}^{(i)}$ of the $i$-th particle are updated according to
+
+$$
+\mathbf{x}^{(i)} := \mathbf{x}^{(i)} + \mathbf{v}^{(i)}\,,
+$$
+
+$$
+\mathbf{v}^{(i)} := w\mathbf{w}^{(i)} +c_1r_1 (\mathbf{x}^{(i)}_\text{best}-\mathbf{x}^{(i)}) +c_2r_2 (\mathbf{x}_\text{best}-\mathbf{x}^{(i)})\,,
+$$
+with $w$, $c_1$, $c_2$ three parameters dermining the behavior, $r_1$ and $r_2$ two random uniform numbers from the [0,1] interval and $\mathbf{x}^{(i)}_\text{best}$ the best design point of particle $i$ and $\mathbf{x}_\text{best}$ the current global best design point.
+
+Because all particles perform both a partly independent search and share information, PSO exhibits an emerging intelligent swarming behavior.
+
+# Implementation
+
+We will construct a new structure for particles containing their position, velocity, and best point.
+
+```julia
+mutable struct Particle{T}
+    x::T
+    v::T
+    x_best::T
+end
+
+Particle(x) = Particle(x, zero(x), x)
+```
+
+Note that we use parametric typing to infer the type of our design points automatically.
+
+A simple function can generate an intial population:
+
+```julia
+"""Generate an intiation population with `n_particles` randomly positioned between the given limits."""
+init_population(n_particles, lims...) = [Particle([(u - l) * rand() + l for (l, u) in lims]) for i in 1:n_particles]
+```
+
+Then we can move to the bulk of the code.
+
+
+```julia
+"""
+    particle_swarm_optimization!(f, population, k_max;
+            w=1, c1=1, c2=1, tracker=nothing)
+
+Performs Particle Swarm Optimization to minimize a function `f`. Give an initial
+vector of particles (type `Particle`) and the `k_max`, the number of iterations.
+
+Optionally set hyperparameters `w`, `c1` and `c2` (default value of 1).
+"""
+function particle_swarm_optimization!(f, population::Vector{Particle}, k_max;
+        w=1, c1=1, c2=1, tracker=nothing)
+    # find best point
+    y_best, x_best = minimum((((f(part.x_best), part.x_best)) for part in population))
+    for k in 1:k_max
+        # update population
+        for particle in population
+            # For you to complete
+        end
+        # this allows us to keep track of things if we want so.
+        tracker isa Nothing || tracker(population)
+    end
+    return y_best, x_best
+end
+```
+
+# Illustration
+
+We will illustrate it on the Ackley function.
+
+```julia
+using STMO.TestFuns
+
+fun = TestFuns.ackley
+
+x1lims = (-10, 10)
+x2lims = (-10, 10)
+
+pobj = heatmap(-10:0.01:10, -10:0.01:10,
+                fun, color=:speed)
+```
+
+We initialize a population.
+
+```julia
+population = init_population(50, x1lims, x2lims)
+```
+
+Adding the points is easy!
+
+```julia
+pobj = heatmap(-10:0.01:10, -10:0.01:10,
+                fun, color=:speed)
+
+for particle in population
+    x = particle.x
+    scatter!([x[1]], [x[2]], color=:orange, label="",
+            markersize=2)
+end
+pobj
+```
+
+**Assignments**:
+1. Complete the `particle_swarm_optimization!` code.
+2. Minimize the `ackley` function (or a different one). What are the effects of the hyperparameters?
+3. (optional) Make an animation of the swarming behavior of the particles. See the [documentation](https://docs.juliaplots.org/latest/animations/) on how to do this. HINT: you might find it useful to run `particle_swarm_optimization!` for a single iteration `k_max` times.
diff --git a/scripts/pso.jl b/scripts/pso.jl
new file mode 100644
index 0000000..e589db3
--- /dev/null
+++ b/scripts/pso.jl
@@ -0,0 +1,116 @@
+#=
+Created on Thursday 23 July 2020
+Last update: -
+
+@author: Michiel Stock
+michielfmstock@gmail.com
+
+Illustration of Particle Swarm Optimization.
+=#
+
+# first we define a particle with a position `x` and a velocity `v`
+
+using STMO.TestFuns
+
+fun = TestFuns.ackley
+x1lims = (-10, 10)
+x2lims = (-10, 10)
+
+mutable struct Particle{T}
+    x::T
+    v::T
+    x_best::T
+end
+
+Particle(x) = Particle(x, zero(x), x)
+
+init_population(n_particles, lims...) = [Particle([(u - l) * rand() + l for (l, u) in lims]) for i in 1:n_particles]
+
+"""
+    particle_swarm_optimization!(f, population, k_max;
+            w=1, c1=1, c2=1, tracker=nothing)
+
+Performs Particle Swarm Optimization to minimize a function `f`. Give an initial
+vector of particles (type `Particle`) and the `k_max`, the number of iterations.
+
+Optionally set hyperparameters `w`, `c1` and `c2` (default value of 1).
+"""
+function particle_swarm_optimization!(f, population::Vector{Particle}, k_max;
+        w=1, c1=1, c2=1, tracker=nothing)
+    # find best point
+    y_best, x_best = minimum((((f(part.x), part.x)) for part in population))
+    for k in 1:k_max
+        # update population
+        for particle in population
+            r1, r2 = rand(2)
+            particle.v .= w * particle.v + c1 * r1 * (particle.x_best .- particle.x) .+
+                c2 * r2 * (x_best .- particle.x)
+            particle.x .+= particle.v
+            fx = f(particle.x)
+            # update current best
+            fx < f(x_best) && (particle.x_best .= particle.x)
+            # update global best
+            if fx < y_best
+                y_best = fx
+                x_best .= particle.x
+            end
+        end
+        tracker isa Nothing || tracker(population)
+    end
+    return y_best, x_best
+end
+
+population = init_population(50, x1lims, x2lims)
+particle_swarm_optimization!(fun, population, 100, w=1)
+
+function pso_animation(f, population, k_max;
+        w=1, c1=1, c2=1)
+    # find best point
+    y_best, x_best = minimum((((f(part.x), part.x)) for part in population))
+    objective = [y_best]
+    # determine xvals and yvals depending on the spread of the particles
+    x1min = minimum((part.x[1] for part in population))
+    x1max = maximum((part.x[1] for part in population))
+    x2min = minimum((part.x[2] for part in population))
+    x2max = maximum((part.x[2] for part in population))
+    x1steps = (x1max - x1min) / 200
+    x2steps = (x2max - x2min) / 200
+    anim = @animate for k in 1:k_max
+
+        swarmplot = heatmap(x1min:x1steps:x1max, x2min:x2steps:x2max,
+                        (x,y)->f((x,y)), color=:speed)
+        xlims!((x1min, x1max))
+        ylims!((x2min, x2max))
+        # update population
+        for particle in population
+            r1, r2 = rand(2)
+            particle.v .= w * particle.v + c1 * r1 * (particle.x_best .- particle.x) .+
+                c2 * r2 * (x_best .- particle.x)
+            particle.x .+= particle.v
+            fx = f(particle.x)
+            # update current best
+            fx < f(x_best) && (particle.x_best .= particle.x)
+            # update global best
+            if fx < y_best
+                y_best = fx
+                x_best .= particle.x
+            end
+            scatter!(swarmplot, [particle.x[1]], [particle.x[2]],
+                        color=:red, label="")
+        end
+        push!(objective, y_best)
+        pobj = plot(0:k, objective, label="objective")
+        xlabel!("iteration")
+        ylabel!("best objective")
+        plot(swarmplot, pobj)
+    end
+    return anim
+end
+
+population = init_population(50, x1lims, x2lims)
+pso_no_momentum = pso_animation(fun, population, 100, w=1, c1=0.9, c2=0.9)
+gif(pso_no_momentum, "figures/pso_no_momentum.gif", fps=4)
+
+population = init_population(50, x1lims, x2lims)
+pso_momentum = pso_animation(fun, population, 100, w=0.8, c1=0.9, c2=0.9)
+gif(pso_momentum, "figures/pso_momentum.gif", fps=4)
diff --git a/src/testfuns.jl b/src/testfuns.jl
index 2fdcc49..1f7e27c 100644
--- a/src/testfuns.jl
+++ b/src/testfuns.jl
@@ -1,6 +1,6 @@
 #=
 Created on Monday 13 January 2020
-Last update: -
+Last update: Saturday 15 August 2020
 
 @author: Michiel Stock
 michielfmstock@gmail.com
@@ -12,23 +12,40 @@ module TestFuns
 
 using LinearAlgebra
 
+function ackley(x; a=20, b=0.2, c=2π)
+    d = length(x)
+    return -a * exp(-b*sqrt(sum(x.^2)/d)) -
+        exp(sum(cos.(c .* x))/d)
+end
+
+ackley(x...; kwargs...) = ackley(x; kwargs...)
+
 # Branin function
-function branin((x1,x2); a=1, b=5.1/(4pi^2), c=5/pi, r=6, s=10, t=1/8pi)
+function branin((x1, x2); a=1, b=5.1/(4pi^2), c=5/pi, r=6, s=10, t=1/8pi)
     return a * (x2 - b * x1^2 + c * x1 - r)^2 + s * (1 - t) * cos(x1) + s
 end
 
+branin(x1, x2; kwargs...) = branin((x1, x2); kwargs...)
+
+# booth
 booth((x1, x2)) = (x1+2x2-7)^2 + (2x1+x2-5)^2
+booth(x1, x2) = booth((x1, x2))
 
 # Rosenbrock
 rosenbrock((x1, x2); a=1, b=5) = (a-x1)^2 + b*(x2-x1^2)^2
+rosenbrock(x1, x2; kwargs...) = rosenbrock((x1, x2); kwargs...)
 
-rastrigine(x; A=10) = length(x) * A + sum(x.^2 .+ A * cos.(2pi*x))
+rastrigine(x; A=10) = length(x) * A + sum(x.^2 .+ A .* cos.(2pi .* x))
+rastrigine(x...; A=10) = rastrigine(x; A=A)
 
 function flower(x; a=1, b=1, c=4)
-    return a * norm(x) + b * sin(c*atan(x[2], x[1]))
+    return a * norm(x) + b * sin(c * atan(x[2], x[1]))
 end
 
-export branin, rosenbrock, rastrigine, flower, booth
+flower(x1, x2; kwargs...) = flower((x1, x2); kwargs...)
+
+
+export ackley, branin, rosenbrock, rastrigine, flower, booth
 
 # functions for convex optimization
 
diff --git a/test/testfuns.jl b/test/testfuns.jl
index 8d1c8ef..8d58070 100644
--- a/test/testfuns.jl
+++ b/test/testfuns.jl
@@ -1,10 +1,12 @@
 @testset "test functions" begin
 
-    import STMO.TestFuns: branin, rosenbrock, rastrigine, flower, booth, fquadr, fnonquadr
+    import STMO.TestFuns: ackley, branin, rosenbrock, rastrigine, flower, booth, fquadr, fnonquadr
 
-    for fun in [branin, rosenbrock, rastrigine, flower, booth]
+    for fun in [ackley, branin, rosenbrock, rastrigine, flower, booth]
 
         @test fun([3, 4]) isa Real
         @test fun([0.0, 0.0]) isa Real
+        @test fun(3, 4) isa Real
+        @test fun(0.0, 0.0) isa Real
     end
 end