Callback for a specific variable in an array while calculating differential equations

Hello everyone,

How would one apply a callback to a variable in an array while calculating differential equations?

I have a complex food web: plants are using nutrients, animals are eating plants,and some animals are being eaten by other animals. I am interested in the change in biomass of this food web (so running differential equations). I want to run the differential equations for this food web until it stabilizes (200 time steps), and then change a harvesting parameter from 0 to 50. But I only want to harvest a specific animal in this food web and then continue running the differential equations.

Is it possible to have such a thing as a callback or is it better for me to forget about the callback and instead run the differential equations for 200 time steps, save the biomass amounts at that time step, manually change the harvest parameter, and then run the differential equations again using the 200 time step biomass amounts as starting values?

I’m not sure that I understand the problem, but do you need a callback for this? Or is your “harvesting parameter”, in fact, a “control input”, i.e., a function of time that you can determine? Consider the following trivial example – is this similar to what you do?

\frac{dN}{dt} = \dot{N}_\mathrm{i} - \dot{N}_\mathrm{e} + \dot{N}_\mathrm{g}.

Here N is the number of, say, animals, \dot{N}_\mathrm{i} is the “immigration” rate of the animals, \dot{N}_\mathrm{e} is the “emmigration” rate, while \dot{N}_\mathrm{g} is the net “generation” rate = births - natural deaths - deaths by “harvesting”.

This is, of course, a much simpler problem than your case, but perhaps it illustrates what you want to do? See below…

In such a description, assume that immigration is given by pressure from external regions, and as such is some function of time, \dot{N}_\mathrm{i}(t) . Emmigration could, typically, be a fraction of the population, say \dot{N}_\mathrm{e} = \epsilon\cdot N. Births and natural deaths are typically proportional to N, e.g., \beta N and \delta N, respectively. Let \dot{N}_\mathrm{h} be the harvesting rate. Then the model can be phrased as:

\frac{dN}{dt} = \dot{N}_\mathrm{i} - \epsilon\cdot N + (\beta - \delta)N - \dot{N}_\mathrm{h}

The harvesting rate could either be given as a time function, \dot{N}_\mathrm{h} = f(t) — say, harvest a certain number of animals in a given time period, or the harvesting rate could be state dependent, e.g., \dot{N}_\mathrm{h} = f(N) (state feedback; the harvesting depends on the population).

If your “harvest parameter” is similar to what I call the “harvesting rate” \dot{N}_\mathrm{h}, then you don’t need callback – you just need a function.

To make the example more concrete…

using DifferentialEquations, Plots; pyplot()
#
# Population model
function Nmod(N,p,t)
    dN = Ndi(t) - p.e*N + (p.b-p.d)*N - Ndh(t)
end
# Model parameters -- named tuple
p = (e=0.1,b=0.05,d=0.01)
# Initial population
N0 = 100.
# Immigration function
function Ndi(t)
    return 5
end
# Harvest function
function Ndh(t)
    if t < 100
        return 0
    else
        return 3
    end
end
#
tspan=(0.,200.)
#
prob = ODEProblem(Nmod,N0,tspan,p);

Solving the problem leads to:

sol = solve(prob)
plot(sol,vars=(0,1))

with result:
image

Here, the drop in N at t=100 is caused by the change in Ndh at t=100.

So — if you by changing the “harvest parameter” mean what I have called the “harvest rate”, you don’t need callback — simply make the harvest rate a function of time.

You could, of course, say that you want to make the harvest rate a function of the number of “animals”. This is what you’d do if you want to approach a “reference population” N_\mathrm{ref} — then you could define a “harvest policy”, say:

\dot{N}_\mathrm{h} = K(N-N_\mathrm{ref})

where K is some suitably chosen constant (“proportional gain” in control speak).

1 Like

Thank you for your reply. The harvest parameter is indeed what you call the harvesting rate and I will try to incorporate it into a harvest function based on your example.

Would you be able to advise me on how I may apply the harvest function to a specific species within an array (in a differential equation)? The array is called bioS and is generated by

using Distributions, DifferentialEquations, Random
# number of nutrients
const num_nutrient = 2
# number of plant species
const num_plant = 15
# number of animal species
const num_animal = 25
# array size
const arr_size = num_nutrient+num_plant+num_animal

function create_bioS(num_animal,num_plant,num_nutrient,arr_size)
    bioS = zeros(arr_size)
    bioS[1:num_nutrient] = rand(Uniform(0,2), num_nutrient)
    bioS[num_nutrient+1:(arr_size-num_animal)] = rand(Uniform(0,10), num_plant)
    bioS[(arr_size-num_animal+1):arr_size] = rand(Uniform(0,10), num_animal)
    return bioS
end
#initial biomass of species with random uniformly distributed values
bioS = create_bioS(num_animal,num_plant,num_nutrient,arr_size)

Then I have separate differential equation for nutrients (dNdT), plants (dPdT), and animals (dAdT). Using your example, how could I specify that I want Ndh(t) to only apply to, for example, bioS[40] ? I am worried that I end up subtracting the harvested biomass of species 40 from each of the 25 animal elements in bioS[(arr_size-num_animal+1):arr_size].

I’ll take a look at it. I’m out of town on a job trip until Saturday, and will respond then. However, my guess would be that (using my notation) if \dot{N}_\mathrm{d}(t) (denoted Nd_h(t)returns an array, you can pick out element 3 (as an example) by the expression Nd_h(t)[3]. But I’ll check.

After thinking about it some more, I think I do need a callback because for my first experiment the rate of harvest does not change - harvesting is either allowed or not, and the amount of harvest is constant.

In the second experiment it will be changing dependent on the population, so @BLI your solution is still very useful. I look forward to your reply.

Just a quick response: I would not use callback here. Instead, I’d do something like:

function Nd_h(t,N,exp1)
    if exp1
        return 10
    else
        K = 2
        Nref = 120
        return K*(N - Nref)
    end
end

If you do:

exp1 = true
Nd_h(t,N,exp1)

then the return value is constant, say 10. If you instead do:

exp1 = false
Nd_h(t,N,exp1)

then the return value is given by the size of the population.

OK… see my modified discussion below. NOTE: my suggestion is based on DifferentialEquations functions calling other functions. I’m not a super expert on Julia, so it is possible that my proposal does not lead to efficient Julia code. On the other hand, it leads to structured code that should be easy to debug.

Importing packages, specifying some plotting “constants”

# Intro packages + define quantities
using Plots; pyplot()
using LaTeXStrings
using DifferentialEquations
#
LW1 = 2.5
LW2 = 1.5
LS1 = :solid
LS2 = :dot
LC1 = :red
LC2 = :blue;

Introducing LW1, etc. is to enforce consistent linewidths, colors, line styles, colors, etc., and make it easy to change these things in one location.

Population model

Comment on notation: in the description below, N is population number within some system boundary, \frac{dN}{dt} is the time derivative of N, while \dot{N} is a rate of N – either (i) the rate of N crossing the system boundary, or (ii) a generation rate within the system boundary. This notation is similar to what is common in some engineering disciplines. Specifically, \dot{N} is not the time derivative of N. In other words, I do not use Newton’s notation for (temporal) derivative.

Consider the population balance

\frac{dN_j}{dt} = \dot{N}_{\mathrm{i},j} - \dot{N}_{\mathrm{e},j} + \dot{N}_{\mathrm{g},j} - \dot{N}_{\mathrm{d},j} + \dot{N}_{\mathrm{m},j}

where N_j is the accumulated value (“population”) of species j, \dot{N}_{\mathrm{i},j} is the influent rate/immigration rate of species j, \dot{N}_{\mathrm{e},j} is the effluent rate/emigration rate of species j, \dot{N}_{\mathrm{g},j} is the growth rate of species j, \dot{N}_{\mathrm{d},j} is the death rate, and \dot{N}_{\mathrm{m},j} is the population management rate.

In the standard Lotka-Volterra model, an unmanaged “island” population is assumed, i.e., \dot{N}_{\mathrm{i},j}\equiv 0, \dot{N}_{\mathrm{e},j}\equiv 0, and \dot{N}_{\mathrm{m},j}\equiv 0. For the predator, growth is based on feeding off the prey, \dot{N}_{\mathrm{g},1} = \gamma_1 N_1 N_2, and death results from natural death/starvation \dot{N}_{\mathrm{d},1} = \delta_1 N_1. For the prey, growth without limiting factors (e.g., unlimited availability of grass) is assumed, \dot{N}_{\mathrm{g},2} = \gamma_2 N_2, and death comes from being consumed by the predator, \dot{N}_{\mathrm{d},1} = \delta_2 N_1 N_2. If prey growth is limited due to resource constraints, a modification is \dot{N}_{\mathrm{g},2} = \gamma_2 N_2\cdot \left( 1-N_2/N_{2,\max} \right); the logistic model.

Some standard parameters for the Lotka-Volterra model are, e.g., \gamma_1 = 0.1, \delta_1 = 0.4, \gamma_2 = 1.1, and \delta_2 = 0.4 (Wikipedia; cheetah and baboon). Typical initial conditions can be N_j(0) = 10.

The DifferentialEquations population model

#
# Population model
# -- population balance
function Nmod(dN,N,p,t)
    dN .= Ndi(t,N,p) .- Nde(N,p) .+ Ndg(N,p) .- Ndd(N,p) .+ Ndm(t,N,p)
end
# -- immigration rate
function Ndi(t,N,p)
    return zeros(size(N))
end
# -- emigration rate
function Nde(N,p)
    return zeros(size(N))
end
# -- growth rate
function Ndg(N,p)
    if p.limited_growth == false
        return [p.g1*N[1]*N[2], p.g2*N[2]]
    else
        return [p.g1*N[1]*N[2], p.g2*N[2]*(1-N[2]/p.N2max)]
    end
end
# -- death rate
function Ndd(N,p)
    return [p.d1*N[1], p.d2*N[1]*N[2]]
end
# -- management rate
function Ndm(t,N,p)
    if p.management_policy == 0
        # zero management
        return [0,0]
    elseif p.management_policy == 1
        # temporal management -- inserting prey for t>= 60
        if t < 60
            return zeros(size(N))
        else
            return [0,1] 
        end
    else
        # state based management
        return [-0.5*max((N[1]-5),0),0]
    end
end
;

Case: unlimited growth, no management

# Model parameters -- named tuple
p = (g1=0.1,g2=1.1,d1=0.4,d2=0.4,N2max=400,limited_growth=false,management_policy=0)
# Initial population
N0 = [10.,10.]
#
tspan=(0.,100.)
#
prob = ODEProblem(Nmod,N0,tspan,p);
#
sol = solve(prob)
plot(sol,vars=(0,1),lc=LC1,lw=LW2,label="Predator")
plot!(sol,vars=(0,2),lc=LC2,lw=LW2,label="Prey")
plot!(title="No management, unlimited prey growth")

The result is as follows:

Case: limited growth, no management

In a more realistic model, the growth rate of the prey is limited by the carrying capacity of the area — e.g., how much grass that is available before the rabbit population starts to eat all of it. The logistic model describes this. Let us set the capacity to N_{2,\max} = 400.

# -- modified growth rate
p = (g1=0.1,g2=1.1,d1=0.4,d2=0.4,N2max=400,limited_growth=true,management_policy=0)
prob = ODEProblem(Nmod,N0,tspan,p)
#
sol = solve(prob)
plot(sol,vars=(0,1),lc=LC1,lw=LW2,label="Predator")
plot!(sol,vars=(0,2),lc=LC2,lw=LW2,label="Prey")
plot!(title="No management, limited prey growth")

with result:

Case: unlimited growth, temporal management (sets in at t\ge 60)

# -- temporal management
p = (g1=0.1,g2=1.1,d1=0.4,d2=0.4,N2max=400,limited_growth=false,management_policy=1)
prob = ODEProblem(Nmod,N0,tspan,p); 
#
sol = solve(prob)
plot(sol,vars=(0,1),lc=LC1,lw=LW2,label="Predator")
plot!(sol,vars=(0,2),lc=LC2,lw=LW2,label="Prey")
plot!(title="Temporal management")

Case: unlimited growth, state based management

Here, the simple policy is used of hunting the predator if its population exceeds 5.

# -- state based management
p = (g1=0.1,g2=1.1,d1=0.4,d2=0.4,N2max=400,limited_growth=false,management_policy=2)
prob = ODEProblem(Nmod,N0,tspan,p); 
#
sol = solve(prob)
plot(sol,vars=(0,1),lc=LC1,lw=LW2,label="Predator")
plot!(sol,vars=(0,2),lc=LC2,lw=LW2,label="Prey")
plot!(title="State based management")

with result:

Remark: Here, the population number N is posed as an integer. The time derivative is meaningless for integers. An “engineering” approach to this is to assume that N is “large” and pretend that N “behaves” as a real number. This is hardly true for N in the order of 10 as above… and the plots show solution values that are non-integers. One rule of thumb seen in the literature is that N should be larger than 10^4 for this approximation to be valid. The requirement for the size of N for this approximation to be true (e.g., pretending that the integer N behaves like a real number) is probably problem dependent. The same difficulty shows up when considering molecules as particles and posing the “number” balance. In many problems, N is so large that this is hardly a problem (i.e., in the order of Avogadro’s number).

1 Like

Oh my… :open_mouth:

[starts reading]