Callback in DifferentialEquation.jl changes the `p` parameter vector. Is this expected behaviour?

Specific Domains Modelling & Simulations
1

This is a MWE showing that p defined outside of the integrator is also modified by the callback.
Is this expected behaviour?

Thanks,

using OrdinaryDiffEq
function rhs(du, u, p, t)
    V, Q, c_in = p
    c = u[1]
    du[1] = (Q/V) * (c_in - c)
end
tspan = (0.0, 10000.0) # Time span for the simulation
u0 = [0.0] # Initial concentrations
c_in = 2e-3 # Inlet concentration
V = 0.1 # Volume of the reactor
Q = 2e-3 # Flow rate in l/s
p = [V, Q, c_in]
@show p
condition(u, t, integrator) = t == 5000.0
affect!(integrator) = integrator.p[3] = 0.0 # Set c_in to 0 at t = 5000s
cb = DiscreteCallback(condition, affect!)
prob = ODEProblem(rhs, u0, tspan, p)
sol = solve(prob, Tsit5(), abstol=1e-8, reltol=1e-8,
 callback=cb, tstops = [5000.0])
@show p
2

Yes it’s expected. A variable does not store an object (unlike some other languages), it’s just a name that can be associated with 1 object at a time. The change you can directly make to a variable is that association, in other words reassigning it. On the other hand, 1 object can be associated with many names or other kinds of references, and a change to an object is mutation.

You never reassign the global p, so it’s always associated with the same vector. You mutated that vector from another reference integrator.p, so when you go to print the object associated with p, you get that mutated vector.

Note that c_in wasn’t reassigned either, so making another vector [V, Q, c_in] contains the original values. Generally when you pass variables into function calls (or expressions like array literals that lower to function calls), the associated objects are inputted, not the variables themselves. In the case of that vector, its elements are assigned to the same objects as the written variables’. Inputting variables somewhat occurs with variable capture by closures (methods in local scopes, comprehensions).

1 Like
3

In the alias controls you can choose for p to not be aliased

4

Thank you for the replies Benny and Chris!
I should have expected that the default behaviour should be to mutate p, but it can be harder to track when we first pass it to a problem and it goes and gets mutated in this integrator object, inside a callback in the solve step.
I am glad I can control this aliasing. This below works:

using OrdinaryDiffEq
using Plots
function rhs(du, u, p, t)
    V, Q, c_in = p
    c = u[1]
    du[1] = (Q/V) * (c_in - c)
end
tspan = (0.0, 10000.0) # Time span for the simulation
u0 = [0.0] # Initial concentrations
c_in = 2e-3 # Inlet concentration
V = 0.1 # Volume of the reactor
Q = 2e-3 # Flow rate in l/s
p = [V, Q, c_in]
@show p
condition(u, t, integrator) = t == 5000.0
affect!(integrator) = integrator.p[3] = 0.0 # Set c_in to 0 at t = 5000s
cb = DiscreteCallback(condition, affect!)

alias = ODEAliasSpecifier(alias_p = false)
prob = ODEProblem(rhs, u0, tspan, p)
sol = solve(prob, Tsit5(), abstol=1e-8, reltol=1e-8,
           callback=cb, tstops = [5000.0], alias = alias)

@show p  # Original p should be unchanged
fig = plot(sol, xlabel="Time (s)", ylabel="Concentration (mol/l)", title="CSTR Concentration Dynamics")
display(fig)