This gist of this question is to determine what the canonical workflow is for defining systems with changing numbers of components e.g. growth/death of biological organisms. This issue is discussed is the callback section of the DiffEq docs where the intended solution is to resize! the integrator according to your callback logic. This makes sense in the DE universe where you directly define julia functions which return the right hand side of your equation.
With MTK, it apparently makes its own rhs, so it wasn’t immediately clear to me how to define a model with a variable number of equations. More to the point, even if you can do such a thing, the resizing is happening at the level of the ODEProblem, so it’s not obvious what order your parameters are in anymore. One solution I cooked up is to decide a priori what the maximum number of things your system can have, and then multiply the rhs of all of your equations by an “activation” number which is in principle either 0 or 1 and is a parameter in the language of MTK. Then, when your callback logic triggers, you modify the appropriate activation variables to switch from 0 to 1 when something grows and vice versa when it dies. It still runs into the issue where you sort of need to know what order the parameters come in, but it seems like they’re in the order they’re defined? In general, maybe this isn’t a good idea since you’re pushing around possibly many extra equations which do nothing until they activate.
Here is a mwe for the activation method. I use Leaf as the thing that is growing/dying and a Tree as the glue between Leafs. To be clear, this doesn’t have anything to do with real trees and leaves. Of course, the ideal situation would be to be able to just add a new Leaf object directly in the callback logic and have everything Just Work™ but that might be tricky.
using ModelingToolkit, DifferentialEquations
struct LeafParameters{S<:Integer, T<:Real}
x_speed::T
y_speed::T
activated::S
end
struct TreeParameters{T<:Real}
stem_strength::T
end
@variables t
ddt = Differential(t)
function Leaf(xy0::Vector{<:Real}, lp::LeafParameters; name::Symbol)
ps = @parameters act = lp.activated xsp = lp.x_speed ysp = lp.y_speed
@variables x(t) = xy0[1] y(t) = xy0[2]
@variables F(t)[1:2]
eqs = [
ddt(x) ~ act*(xsp + F[1]),
ddt(y) ~ act*(ysp + F[2])
]
return ODESystem(eqs, t, [x, y, F[1:2]...], ps; name)
end
function Tree(xy0::Matrix{<:Real}, lp::Vector{<:LeafParameters}, tp::TreeParameters; name::Symbol)
# a tree is a composition of leaves
N_leaves = size(xy0, 1)
@named leaf 1:N_leaves i -> Leaf(xy0[i,:], lp[i])
ps = @parameters tps = tp.stem_strength
# in practice, these would be complicated interaction forces
forces_x = [
leaf[i].F[1] ~ tps for i = 1:N_leaves
]
forces_y = [
leaf[i].F[2] ~ tps for i = 1:N_leaves
]
eqs = [forces_x ; forces_y]
return compose(ODESystem(eqs, t, [], ps; name = name), leaf[1:N_leaves]...)
end
n_leaves = 3
xy0 = [rand() for i = 1:n_leaves, j = 1:2]
lp = [
LeafParameters(1.1, 2.1, 1),
LeafParameters(1.2, 2.2, 0),
LeafParameters(1.3, 2.3, 0)
]
tp = TreeParameters(3.0)
@named tree = Tree(xy0, lp, tp)
tree = structural_simplify(tree)
t_range = (0.0, 10.0)
prob = ODEProblem(
tree,
[],
t_range
)
# prob.p = [3.0, 1.0, 1.1, 2.1, 0.0, 1.2, 2.2, 0.0, 1.3, 2.3]
function grow_leaf(t1::Real)
condition(u, t, integrator) = t1 - t
function affect!(integrator)
u = integrator.u
u[3] = 0.0 # x
u[4] = 0.0 # y
integrator.p[5] = 1.0 # 5 is the index of the second Leaf's activation
return nothing
end
return ContinuousCallback(condition, affect!)
end
@time sol = solve(prob, saveat = 1.0, callback = grow_leaf(3.0))