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.644342 seconds (4.26 M allocations: 230.701 MiB, 6.67% gc time, 99.90% compilation time)

retcode: Success
Interpolation: specialized 4th order "free" interpolation
t: 2281-element Vector{Float64}:
    0.0
    0.1379730351325356
    0.5564761658770004
    1.1196169933783673
    1.7672551586109673
    2.552567516636946
    3.4680808539852217
    4.564946454459349
    5.876373627177295
    7.483768553665588
    ⋮
 7281.577634096674
 7282.401079280315
 7284.161123198842
 7285.691532471932
 7287.821073788626
 7290.147654177271
 7293.180461584833
 7296.866379309829
 7300.0
u: 2281-element Vector{Vector{Float64}}:
 [400.0]
 [410.58300851936434]
 [437.32809938708135]
 [462.27758078114454]
 [479.52406721317647]
 [490.4502350704483]
 [496.1319303647314]
 [498.7002281857691]
 [499.64738194876736]
 [499.9274082148435]
 ⋮
 [435.6230523006657]
 [422.649927421847]
 [410.67040617502374]
 [407.12322053366086]
 [405.5084068863958]
 [405.0599618315321]
 [404.9429081287553]
 [404.9230484163387]
 [404.919936730924]

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.