kind of organizational equation. I have a dynamical system whose equations I change frequently. I add and remove different dynamical components and connections. I’m in a pretty exploratory stage. It becomes really difficult to keep track of the parameter container. Different connections have different parameters and I need to be changing things all over the place each time I want to change a connection. I’m wondering, has anyone gone through a similar phase and found good workflows?
To give you an idea, let’s say I have something like:
function ebm(u, p, t)
T, TD = u
@unpack S, c, ε = p
# C = cloudiness(T, TD)
# ocean_heat_update = γ*(T-TD)
Sf = solar_forcing(S, Sa, t)
α = ice_albedo_feedback(T) #+ C
dT = Sf*(1-α) - ε*(σ₄*T)^4 - eₒ*η*(T-TD)
dTD = η*(T-TD)
return SVector(dT/c, dTD/cD)
end
ice_albedo_feedback(T) = 0.5 - 0.2*tanh((T-263)/4)
solar_forcing(S, Sa, t) = S + Sa*cos(t)
# Parameters
c = 1.0
cD = 10.0
S = 340.0 # solar constant, W/m²
Sa = 0.0
ε = 0.65
eₒ = 1.0
η = 0.7
p = @dict c cD S Sa ε eₒ η
T₀ = 280.0
TD₀ = 260.0
u₀ = [T₀, TD₀]
ds = ContinuousDynamicalSystem(ebm, u₀, p)
# or prob = ODEProblem(ebm, ...)
every time I want to change a connection I need to define my parameter container differently, I need to comment/uncomment different lines in the ebm function, and I need to unpack different parameters.
I’m not quite sure it’ll suit your workflow, but I’ve previously made p a struct which means I can have multiple definitions of ebm that dispatch on p, and it’s then quite easy to switch which set of equations is being solved by just passing in a different type. E.g. you could have a certain struct for a set of components and connections. That may work if you’re just alternating between different sets of components and connections, but if you’re adding new parameters/changing the implementation then using structs may actually be more restrictive.
Honestly, that’s a pretty good idea. I can define all dynamic connections as functions of an abstract parameter. The default dispatch for all dynamic connections would be to return 0 or a standard evaluation. Specialized sets can do things differently. The only point where I see this suffering is when I need to add new dynamic variables (change the ODE problem from 2D to 3D). then this doesn’t help as much…
I wondering if I can use MTK to dispatch on the variable names… If the system state has the “C” variable, do this, if the system state doesn’;t have the “C” variable, do that… (but without using if-else for performance reasons)
@Datseris nested structs how @mmiller describes (with fomulation variants as multiple structs and dispatch methods) is how all of my models work too. And Tuple/NamedTuple for the parts that are in flux or are not duplicated.
ModelParameters.jl is good to get the parameter values out when you need them as a tuple/vector/csv, or if you need metadata for them.
I’m using ModelingToolkit in a similar way – exploratory analysis of system behavior when tinkering with the governing equations. It’s the same issue (although maybe a little better than what you’re facing): needing to edit the code in multiple places consistently (in the system of equations but also in the parameter/variable definitions). It would be nice to auto-populate the list of variables as you add them to the system of equations… but you also need to specify initial values (and potentially units) as well as whether to consider them variables or parameters. So I’m not sure how much one can really do to simplify this.
With MTK, you can modify your equations and the parameters that are used will be kept track of
using ModelingToolkit
@parameters t p=0 s=0
@variables x(t)=0 y(t)=0
D = Differential(t)
eqs = [D(x) ~ -x + p, D(y) ~ -y + s]
@named sys = ODESystem(eqs, t)
parameters(sys) # returns `p, s`
@named sys2 = ODESystem(eqs[1], t) # only add the first equation
parameters(sys2) # returns `p`
You can also give your parameters default values (like above) and provide initial conditions for state variables.
julia> sys
Model sys with 2 equations
States (2):
x(t) [defaults to 0]
y(t) [defaults to 0]
Parameters (2):
p [defaults to 0]
s [defaults to 0]
A common pattern to defining model components using MTK is to create constructors for ODESystem, e.g., the definition of an Integrating block is used in the definition for a PID controller, if you modify the integrator definition, it will be updated everywhere where the constructor is called. In the PID controller below, the state x of the integrator is renamed to I.x and similarly with it’s parameters I.k.
julia> using ModelingToolkitStandardLibrary.Blocks: PID
julia> @named pid = PID(k=1, Ti=1)
Model pid with 9 equations
States (11):
e(t) [defaults to 0]
u_r(t) [defaults to 0]
u_y(t) [defaults to 0]
ep(t) [defaults to 0]
ea(t) [defaults to 0]
sat₊y(t) [defaults to 0]
sat₊u(t) [defaults to 0]
I₊u(t) [defaults to 0]
I₊y(t) [defaults to 0]
y(t)
I₊x(t) [defaults to 0]
Parameters (6):
wp [defaults to 1]
Ni [defaults to 0.001]
k [defaults to 1]
I₊k [defaults to 1.0]
sat₊y_max [defaults to Inf]
sat₊y_min [defaults to -Inf]
Thanks @baggepinnen , this seems very promising and I’d like to try it out. I am noticing a couple of key differences however. The main being, I don’t want to be modifying existing terms in the equations, like e.g., your integrator, but rather swapping them over for different ones. Let’s assume that my Left Hand Side is fixed, then I have equations looking like:
D(x) ~ S + L - T
D(y) ~ +T
in a over-simplified manner. All letters here depend on (x,y) in some way. I want to be able to swap different equations for the capital letters. Does this still make sense with your ODESystem approach? I.e., would you recommend that I define something like
and then in my “final” (to-be-integrated) ODESystem do something like:
eqs = [D(x) ~ S_version_1.S + L - T,
D(y) ~ +T]
# or
eqs = [D(x) ~ S_version_2.S + L - T,
D(y) ~ +T]
EDIT: Okay so, even though this initially much more complex that using standard function and numbers, I do see some clear benefits. I do not need to keep track of what is the “independent variables x, y”, what is the “independent variable t”, and what are parameters “p, d”.
So the only part that remains as the “difficult” part is that I would still need to use multiple dispatch or some other method to differentiate in my final equations whether I should use version 2 of S or version 1.
You seem to have the right idea of how to use it. It’s often convenient to define your own “constructors” for ODESystems that can take arguments that affect the dynamics. For instance, you could define something like
function my_system(version = 1)
@parameters t p=0 d=0
@variables x(t) y(t) S(t)
rhs = version == 1 ? -p*x + y : x^2 + y/d
S = ODESystem([S ~ rhs], t)
eqs = [D(x) ~ S.S + L - T,
D(y) ~ +T]
ODESystem(eqs, t, systems=[S])
end