ABM + DiffEq

https://juliadynamics.github.io/Agents.jl/stable/examples/diffeq/

using Agents
using Distributions
using CairoMakie
using OrdinaryDiffEq
using DiffEqCallbacks
CairoMakie.activate!(px_per_unit = 1.0)

Fisher agents

@agent struct Fisher(NoSpaceAgent)
    competence::Int
    yearly_catch::Float64
end

Set fishing quota

function fish!(integrator, model)
    integrator.p[2] = integrator.u[1] > model.min_threshold ?
        sum(a.yearly_catch for a in allagents(model)) : 0.0
    Agents.step!(model, 1)
end

fish! (generic function with 1 method)

Fish population change

function fish_stock!(ds, s, p, t)
    max_population, h = p
    ds[1] = s[1] * (1 - (s[1] / max_population)) - h
end

fish_stock! (generic function with 1 method)

Agents catch the fish

function agent_cb_step!(agent, model)
    agent.yearly_catch = rand(abmrng(model), Poisson(agent.competence))
end

agent_cb_step! (generic function with 1 method)

ABM model embedded in the callback

function initialise_cb(; min_threshold = 60.0, nagents = 50)
    model = StandardABM(Fisher; agent_step! = agent_cb_step!,
                        properties = Dict(:min_threshold => min_threshold))
    for _ in 1:nagents
        competence = floor(rand(abmrng(model), truncated(LogNormal(), 1, 6)))
        add_agent!(model, competence, 0.0)
    end
    return model
end

initialise_cb (generic function with 1 method)

Setup the problem

modelcb = initialise_cb()
tspan = (0.0, 20.0 * 365.0)
initial_stock = 400.0
max_population = 500.0
prob = OrdinaryDiffEq.ODEProblem(fish_stock!, [initial_stock], tspan, [max_population, 0.0])

# Each Dec 31st, we call fish! that adds our catch modifier to the stock, and steps the ABM model
fish = DiffEqCallbacks.PeriodicCallback(i -> fish!(i, modelcb), 364)
# Stocks are replenished again
reset = DiffEqCallbacks.PeriodicCallback(i -> i.p[2] = 0.0, 365)

@time sol = solve(prob, Tsit5();  callback = CallbackSet(fish, reset))

  1.584149 seconds (2.63 M allocations: 179.848 MiB, 3.36% gc time, 99.90% compilation time)

retcode: Success
Interpolation: specialized 4th order "free" interpolation
t: 2275-element Vector{Float64}:
    0.0
    0.1379730351325356
    0.5564761658770004
    1.1196169933783673
    1.7672552058964655
    2.5525676212598403
    3.4680808264583267
    4.564946465790113
    5.876373568455959
    7.483768743428563
    ⋮
 7281.40838653234
 7282.344654585563
 7284.168706032151
 7285.902426606168
 7288.122506661588
 7290.645031476024
 7293.789405797698
 7297.620970271196
 7300.0
u: 2275-element Vector{Vector{Float64}}:
 [400.0]
 [410.58300851936434]
 [437.32809938708135]
 [462.27758078114454]
 [479.5240681416804]
 [490.45023605017735]
 [496.13193025995974]
 [498.7002282003409]
 [499.647381928803]
 [499.92740822632885]
 ⋮
 [427.18787255980135]
 [407.2447311577636]
 [390.15557608392726]
 [384.2556700810025]
 [381.5906355290433]
 [380.7083947139732]
 [380.44885486693664]
 [380.39430289523807]
 [380.3870229506006]

discrete = vcat(sol(0:365:(365 * 20))[:,:]...)
f = Figure(size = (600, 400))
ax = f[1, 1] = Axis(
        f,
        xlabel = "Year",
        ylabel = "Stock",
        title = "Fishery Inventory",
    )
lines!(ax, discrete, linewidth = 2, color = :blue)
f

---

This notebook was generated using Literate.jl.