# Flexible workflow for altering the equations of a dynamical system / diffeq

**URL:** <https://discourse.julialang.org/t/flexible-workflow-for-altering-the-equations-of-a-dynamical-system-diffeq/71390>\
**Category:** Modelling & Simulations\
**Tags:** differentialequation\
**Created:** [November 12, 2021, 3:48pm UTC](https://discourse.julialang.org/t/flexible-workflow-for-altering-the-equations-of-a-dynamical-system-diffeq/71390 "2021-11-12T15:48:26Z")\
**Posts on this page:** 13\
**Page:** 1

<div class="post-metadata">

**Author:** ![Datseris](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/datseris/32/13406_2.png) [@Datseris](https://discourse.julialang.org/u/Datseris)\
**Post date:** [November 12, 2021, 3:48pm UTC](https://discourse.julialang.org/t/flexible-workflow-for-altering-the-equations-of-a-dynamical-system-diffeq/71390/1 "2021-11-12T15:48:26Z")

</div>

Hi all,

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:

```julia
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.

---

<div class="post-metadata">

**Author:** ![baggepinnen](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/baggepinnen/32/693_2.png) [@baggepinnen](https://discourse.julialang.org/u/baggepinnen)\
**Post date:** [November 12, 2021, 3:50pm UTC](https://discourse.julialang.org/t/flexible-workflow-for-altering-the-equations-of-a-dynamical-system-diffeq/71390/2 "2021-11-12T15:50:07Z")

</div>

It sounds like you would benefit from ModelingToolkit.jl, it keeps track of parameters for you

---

<div class="post-metadata">

**Author:** ![Datseris](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/datseris/32/13406_2.png) [@Datseris](https://discourse.julialang.org/u/Datseris)\
**Post date:** [November 12, 2021, 3:56pm UTC](https://discourse.julialang.org/t/flexible-workflow-for-altering-the-equations-of-a-dynamical-system-diffeq/71390/3 "2021-11-12T15:56:22Z")

</div>

I’m aware of ModelingToolkit.jl but I don’t see how it helps me. You still declare the parameters at the top, no? [https://github.com/SciML/ModelingToolkit.jl#high-level-examples](https://github.com/SciML/ModelingToolkit.jl#high-level-examples)

---

<div class="post-metadata">

**Author:** ![Datseris](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/datseris/32/13406_2.png) [@Datseris](https://discourse.julialang.org/u/Datseris)\
**Post date:** [November 12, 2021, 3:58pm UTC](https://discourse.julialang.org/t/flexible-workflow-for-altering-the-equations-of-a-dynamical-system-diffeq/71390/4 "2021-11-12T15:58:55Z")

</div>

Oh, an obvious benefit is that I don’t have to change my “unpacking” inside the equations, as for MTK the parameters are globals anyway.

---

<div class="post-metadata">

**Author:** ![mmiller](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/mmiller/32/20805_2.png) [@mmiller](https://discourse.julialang.org/u/mmiller)\
**Post date:** [November 12, 2021, 4:25pm UTC](https://discourse.julialang.org/t/flexible-workflow-for-altering-the-equations-of-a-dynamical-system-diffeq/71390/5 "2021-11-12T16:25:11Z")

</div>

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 `struct`s may actually be more restrictive.

---

<div class="post-metadata">

**Author:** ![Datseris](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/datseris/32/13406_2.png) [@Datseris](https://discourse.julialang.org/u/Datseris)\
**Post date:** [November 12, 2021, 4:29pm UTC](https://discourse.julialang.org/t/flexible-workflow-for-altering-the-equations-of-a-dynamical-system-diffeq/71390/6 "2021-11-12T16:29:14Z")

</div>

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…

---

<div class="post-metadata">

**Author:** ![Datseris](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/datseris/32/13406_2.png) [@Datseris](https://discourse.julialang.org/u/Datseris)\
**Post date:** [November 12, 2021, 4:33pm UTC](https://discourse.julialang.org/t/flexible-workflow-for-altering-the-equations-of-a-dynamical-system-diffeq/71390/7 "2021-11-12T16:33:28Z")

</div>

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)

---

<div class="post-metadata">

**Author:** ![Raf](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/raf/32/3383_2.png) [@Raf](https://discourse.julialang.org/u/Raf)\
**Post date:** [November 12, 2021, 5:12pm UTC](https://discourse.julialang.org/t/flexible-workflow-for-altering-the-equations-of-a-dynamical-system-diffeq/71390/8 "2021-11-12T17:12:42Z")

</div>

@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.

[https://github.com/rafaqz/ModelParameters.jl](https://github.com/rafaqz/ModelParameters.jl)

---

<div class="post-metadata">

**Author:** ![lamorton](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/lamorton/32/21115_2.png) [@lamorton](https://discourse.julialang.org/u/lamorton)\
**Post date:** [November 13, 2021, 7:34pm UTC](https://discourse.julialang.org/t/flexible-workflow-for-altering-the-equations-of-a-dynamical-system-diffeq/71390/9 "2021-11-13T19:34:40Z")

</div>

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.

---

<div class="post-metadata">

**Author:** ![baggepinnen](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/baggepinnen/32/693_2.png) [@baggepinnen](https://discourse.julialang.org/u/baggepinnen)\
**Post date:** [November 14, 2021, 9:10am UTC](https://discourse.julialang.org/t/flexible-workflow-for-altering-the-equations-of-a-dynamical-system-diffeq/71390/10 "2021-11-14T09:10:23Z")

</div>

With MTK, you can modify your equations and the parameters that are used will be kept track of

```julia
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
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](https://github.com/SciML/ModelingToolkitStandardLibrary.jl/blob/986436310582194c78d3200326e34bc770b2c2bb/src/Blocks/continuous.jl#L21) is used in [the definition for a PID controller](https://github.com/SciML/ModelingToolkitStandardLibrary.jl/blob/986436310582194c78d3200326e34bc770b2c2bb/src/Blocks/continuous.jl#L154), 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
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]

```

---

<div class="post-metadata">

**Author:** ![AshtonSBradley](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/ashtonsbradley/32/11140_2.png) [@AshtonSBradley](https://discourse.julialang.org/u/AshtonSBradley)\
**Post date:** [November 17, 2021, 4:50pm UTC](https://discourse.julialang.org/t/flexible-workflow-for-altering-the-equations-of-a-dynamical-system-diffeq/71390/11 "2021-11-17T16:50:09Z")

</div>

To save one edit you can just do @unpack p  
to access everything. Since it’s just providing references, it shouldnt be a performance hit

---

<div class="post-metadata">

**Author:** ![Datseris](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/datseris/32/13406_2.png) [@Datseris](https://discourse.julialang.org/u/Datseris)\
**Post date:** [November 24, 2021, 2:30pm UTC](https://discourse.julialang.org/t/flexible-workflow-for-altering-the-equations-of-a-dynamical-system-diffeq/71390/12 "2021-11-24T14:30:38Z")

</div>

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:

```julia
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

```julia
@parameters t p=0 d=0
@variables x(t) y(t) S(t)
S_version_1 = ODESystem([S ~ -p*x + y], t)
S_version_2 = ODESystem([S ~ x^2 + y/d], t)

```

and then in my “final” (to-be-integrated) ODESystem do something like:

```julia
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.

---

<div class="post-metadata">

**Author:** ![baggepinnen](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/baggepinnen/32/693_2.png) [@baggepinnen](https://discourse.julialang.org/u/baggepinnen)\
**Post date:** [November 24, 2021, 2:54pm UTC](https://discourse.julialang.org/t/flexible-workflow-for-altering-the-equations-of-a-dynamical-system-diffeq/71390/13 "2021-11-24T14:54:41Z")

</div>

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

```julia
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

```

or any other mechanism for choosing the dynamics