Composing open dynamical systems II: Undirected composition

By Sophie Libkind and James Fairbanks on January 25, 2021

In this series of posts on open dynamical systems, we explore various theories of composition for open dynamical systems (shown in the diagram below) and use these theories to build complex systems out of primitive ones. In the previous post, we explored directed theories of composition (shown in purple). In this post, we will focus on undirected theories of composition (shown in blue) and give examples of composing dynamical systems using the AlgebraicDynamics.jl package.





Composing Resource Sharers

To compose dynamical systems, we will need to give them some extra structure in addition to their dynamics. For a dynamical system which composes via directed communication, a machine is the data structure that captures this additional structure. We say that machines are dynamical systems that compose via directed communication. Now we define resources sharers which are dynamical systems that compose via undirected composition. Informally, a resource sharer consists of four components

Throughout this post, we will assume that the evolution rule is an ODE, although discrete-time dynamics is also supported by AlgebraicDynamics.jl.

A pool of water

Imagine two resource sharers that model the evolution of a pool of water over time. One system models the emptying of the pool via evaporation, and the other system models the filling of the pool via rainfall. Each resource sharer has a state variable representing the volume of water in the pool. We can compose these systems by hooking them up to the same pool of water, and they communicate through this shared resource. Since each system may both affect and be affected by the shared pool of water, we consider this style of communication between systems to be undirected. The slogan for composing resource sharers is "add along shared coordinates" since the effect on the pool of water by the composite system is the sum of the effects on the pool of water by the primitive systems.

Recall the general approach for composing dynamical systems:

  1. Choose a theory of composition

  2. Choose a composition pattern in that theory

  3. Define primitive systems

  4. Compose!

Let's sketch the general approach for our pool of water example.

1. Choose a theory of composition. In this case (and throughout this post), we choose the theory of undirected wiring diagrams, \(\mathsf{Th}(\mathsf{UWD})\), which is defined by the schema below.[1]


An undirected wiring diagram is a \(\mathsf{C}\)-Set over \(\mathsf{Th}(\mathsf{UWD})\), i.e. a functor \(d: \mathsf{Th}(\mathsf{UWD}) \to \mathsf{Set}\). It consists of the following data:

  1. a set of boxes \(d(\textrm{Box})\) and a set of junctions \(d(J)\)

  2. every box \(b\) has a set of ports \(d(\operatorname{box})^{-1}(b) \subseteq P\)

  3. a junction \(j\) identifies the ports \(d(\operatorname{junc})^{-1}(j) \subseteq P\) and is exposed by the outer ports \(d(\operatorname{junc}')^{-1}(j) \subseteq P'\)

The theory of undirected wiring diagrams is defined in Catlab.WiringDiagrams, and we can implement an undirected wiring diagram by using the @relation macro.

How will we interpret an undirected wiring diagram as a syntax for composing resource sharers? An undirected wiring diagram specifies how the ports of resource sharers connect to junctions. If multiple ports are connected to the same junction then the variables exposed by the ports are identified in the composite system. Therefore, a junction represents a collection of identified variables, i.e. a shared resource. Outer ports map to junctions and represent the ports of the composite system. [2]

[1] Undirected wiring diagrams form an operad where operadic composition is given by a pushout, shown in (Spivak 2013). Here you can learn more about the operadic composition of undirected wiring diagrams as well as the @relation macro for specifying undirected wiring diagrams.
[2] This interpretation of composition of dynamical systems is mathemetized by the cospan algebra \(\mathsf{Dynam}\) defined in (Baez and Pollard 2017).
2. Choose a composition pattern. Since our composition theory is \(\mathsf{Th}(\mathsf{UWD})\), our composition pattern will be an undirected wiring diagram. We will use the @relation macro in order to implement the composition pattern for the pool of water example. The @relation macro provides a convenient syntax for specifying undirected wiring diagrams via a dialect of Julia code designed for graphical regular logic. The variables are captured as junctions, and function calls represent a system acting on those variables. Ports and connections in the undirected wiring diagram are inferred from the relationships between variables and functions.

using Catlab, Catlab.WiringDiagrams, Catlab.Programs, Catlab.Graphics

pool_pattern = @relation (water,) begin
    evaporation(water)
    rainfall(water)
end

The two boxes designate the two systems to be composed, and the single junction represents the shared resouce — the pool of water.

3. Choose primitive systems to compose. The system modeling the evaporation of water has one variable corresponding to the water in the pool. Let's assume we live in a climate with constant evaporation so the dynamics are given by the ODE

\[\dot w_e(t) = -\alpha w_e(t).\]

The system modeling the rainfall also has one variable corresponding to the water in the pool. Let's assume we live in a climate with distinct rainy and dry seasons. These oscillating dynamics are given by the ODE

\[\dot w_r(t) = \beta (1 + \sin(t)).\]

We will implement these primitive systems as resource sharers, and both resource sharers will have a single port which exposes the water variable. We use a labeled array for our parameter vector p to easily keep track of the systems' parameters.

using AlgebraicDynamics, AlgebraicDynamics.UWDDynam
using LabelledArrays

dotwₑ(u,p,t) = -p.α*u[1]
dotwᵣ(u,p,t) = [p.β * (1 + sin(t))]

evaporation = ContinuousResourceSharer{Float64}(1, 1, dotwₑ, [1])
rainfall    = ContinuousResourceSharer{Float64}(1, 1, dotwᵣ, [1])

4. Compose! Following the slogan "add along shared coordinates," we add the effects of the evaporation and rainfall systems on the pool of water. We get the composite ODE

\[\dot w(t) = -\alpha w(t) + \beta (1 + \sin t).\]
pool_system = oapply(pool_pattern, [evaporation, rainfall])

Let's turn this composite system into an ODEProblem, solve, and plot!

using OrdinaryDiffEq, Plots, Plots.Measures

u0 = [0.0] # initial condition = empty pool of water
params = LVector(α = 0.3, β = 0.5)
nyears = 4
tspan = (0.0, nyears * 2 * pi)
prob = ODEProblem(pool_system, u0, tspan, params)
sol = solve(prob, Tsit5())

For the inspired reader, here are some suggested exercises:

\[\dot w_h = \gamma(2 - w_h).\]

Complex compositions

Careful bookkeeping is endemic to implementing a complicated composition pattern. Previously we applied pullback functorial data migration in order to construct complex composition patterns more easily and with less room for error. In this post, we give two more strategies, which both rely on the rich structure of the composition syntax:

  1. Hierarchical composition, in which a primitive system can itself be a composition of systems. Hierarchical composition is extremely useful when we want to divide and conquer or when we want to build upon a composite system without starting from scratch. This strategy is only available in open composition syntaxes such as \(\mathsf{Th}(\mathsf{UWD})\), \(\mathsf{Th}(\mathsf{DWD})\), \(\mathsf{Th}(\mathsf{OpenCPG})\).

  2. Taking (co)limits of composition patterns. This strategy is available in all of the composition syntaxes discussed so far.

We will demonstrate these two strategies as we construct a complex ecosystem with six species: rabbits, foxes, hawks, little fish, big fish, and sharks.

Divide and conquer

First we will apply hierarchical composition so that we don't have to simultaneously study all of the interactions between the six species. We will separately model the three sub-ecosystems — land, air, and water — which may themselves be compositions of even more primitive system. Once we have constructed these sub-systems, we will compose them according to the following pattern:

eco_pattern = @relation (rabbits, foxes, hawks, littlefish, bigfish, sharks) begin
    land_eco(rabbits, foxes, hawks)
    air_eco(hawks, littlefish)
    ocean_eco(littlefish, bigfish, sharks)
end

In other words, there aren't two independent hawk populations for the land and air sub-ecosystems. So in the final model we identify the hawk population from the land ecosystem and the hawk population from the air ecosystem. Likewise, for little fish.

Land

Let's zoom-in on the land ecosystem. We have five primitive systems which share variables:

  1. rabbit growth: \(\dot r(t) = \alpha_1 r(t)\)

  2. rabbit/fox predation: \(\dot r(t) = -\beta_1 r(t) f(t),\; \dot f(t) = \gamma_1 r(t)f(t)\)

  3. fox decline: \(\dot f(t) = -\delta_1 f(t)\)

  4. rabbit/hawk predation: \(\dot r(t) = -\beta_2 r(t)h(t),\; \dot h(t) = \gamma_2 r(t)h(t)\)

  5. hawk decline: \(\dot h(t) = -\delta_2 h(t)\)

Therefore, the desired composition pattern has five boxes and many ports and wires to keep track of. Instead of implementing this composition pattern by hand, we construct it as a pushout of two simpler composition patterns — one which represents the rabbit/fox interactions and one which represents the rabbit/hawk interactions.

For the rabbit/fox composition pattern note that there are not two independent rabbit populations — one that grows and one that gets eaten by foxes. Likewise, there are not two independent fox populations — one that declines and one that feasts on rabbits. To capture these interactions, we identify the two rabbit populations and identify the two fox populations via the following composition pattern.

rabbitfox_pattern = @relation (rabbits, foxes) begin
    rabbit_growth(rabbits)
    rabbit_fox_predation(rabbits,foxes)
    fox_decline(foxes)
end

The rabbit/hawk composition pattern is identical.

rabbithawk_pattern = @relation (rabbits, hawks) begin
    rabbit_growth(rabbits)
    rabbit_hawk_predation(rabbits,hawks)
    hawk_decline(hawks)
end

Now we construct the complete composition pattern for the land ecosystem by gluing these two composition patterns along the box corresponding to rabbit growth.

using Catlab.CategoricalAlgebra

# Define the composition pattern for rabbit growth
rabbit_pattern = @relation (rabbits,) -> rabbit_growth(rabbits)

# Define transformations between the composition patterns
rabbitfox_transform  = ACSetTransformation((Box=[1], Junction=[1], Port=[1], OuterPort=[1]), rabbit_pattern, rabbitfox_pattern)
rabbithawk_transform = ACSetTransformation((Box=[1], Junction=[1], Port=[1], OuterPort=[1]), rabbit_pattern, rabbithawk_pattern)

# Take the pushout
land_pattern = ob(pushout(rabbitfox_transform, rabbithawk_transform))

Phew! After seeing the complexity of this composition pattern, we are grateful that we used a pushout instead of constructing it by hand. Lastly, we define the primitive systems and compose.

dotr(u,p,t)  = p.α₁*u
dotrf(u,p,t) = [-p.β₁*u[1]*u[2], p.γ₁*u[1]*u[2]]
dotf(u,p,t)  = -p.δ₁*u
dotrh(u, p, t) = [-p.β₂*u[1]*u[2], p.γ₂*u[1]*u[2]]
doth(u, p, t)  = -p.δ₂*u

# Define the primitive systems
rabbit_growth       = ContinuousResourceSharer{Float64}(1, 1, dotr,  [1])
rabbitfox_predation = ContinuousResourceSharer{Float64}(2, 2, dotrf, [1,2])
fox_decline         = ContinuousResourceSharer{Float64}(1, 1, dotf,  [1])
rabbithawk_predation= ContinuousResourceSharer{Float64}(2, 2, dotrh, [1,2])
hawk_decline        = ContinuousResourceSharer{Float64}(1, 1, doth,  [1])

# Compose
land_sys= oapply(land_pattern, [rabbit_growth, rabbitfox_predation, fox_decline, rabbithawk_predation, hawk_decline])

The resource sharer land_sys models the land ecosystem. Although it will play the role of a primitive system in our total ecosystem, we see that it is a composite itself. This structure — primitive systems themselves being composites — is the essence of hierarchical composition.

Air

The air ecosystem is straightforwardly defined by the following resource sharer which models hawk/little fish predation. Here we don't model the decay of hawks or the growth of little fish because these processes are already accounted for in the land and water ecosystems. The air ecosystem is a pure coupling term between the two systems.

dothf(u,p,t) = [p.γ₃*u[1]*u[2], -p.β₃*u[1]*u[2]]
air_sys = ContinuousResourceSharer{Float64}(2, 2, dothf, [1,2])

Water

In the previous post, we defined a ocean ecosystem as the directed composition of machines. In this aquatic foodchain, sharks eat big fish and big fish eat little fish. We can turn the machine representing this ocean ecosystem into a resource sharer as follows.

water_sys = ContinuousResourceSharer{Float64}(3, (u,p,t)->eval_dynamics(ocean_sys, u, [], p, t))

Putting it all together

Now that we have constructed resource sharers for the three sub-ecosystems, we are ready to plug them into the established composition pattern, eco_pattern.

# Compose
eco_system = oapply(eco_pattern, [land_sys, air_sys, water_sys])

# Plot and solve
u0 = [100.0, 50.0, 20.0, 100, 10, 2.0]
params = LVector(
    α₁ = 0.3, β₁ = 0.015, γ₁ = 0.015, δ₁ = 0.7, β₂ = .01, γ₂ = .01, δ₂ = .5,  # land params
    γ₃ = 0.001, β₃ = 0.003, # air params
    α₄ = 0.35, β₄ = 0.015, γ₄ = 0.015, δ₄ = 0.7, β₅ = 0.017, γ₅ = 0.017, δ₅ = 0.35 # water params
)
tspan = (0.0, 75.0)
prob = ODEProblem(eco_system, u0, tspan, params)
sol = solve(prob, Tsit5())
plot(sol, lw=2, label = ["rabbits" "foxes" "hawks" "little fish" "big fish" "sharks"])

A successful example of hierarchical composition!

To really see the benefits of hierarchical composition, let's take a look at how we could have achieved this composition without the layered approach.

# Define the (flattened) composition pattern
flattened_pattern = ocompose(eco_pattern, 1, land_pattern)

#Compose
eco_sys2 = oapply(flattened_pattern, [rabbit_growth, rabbitfox_predation, fox_decline, rabbithawk_predation, hawk_decline, air_sys, water_sys])

While this strategy produces the same composite system,[3] it lacks the clarity of the hierarchical approach which we can see from the flattened composition pattern:

[3] Functoriality of the cospan algebra \(\mathsf{Dynam}\) implies that the resource sharers eco_sys and eco_sys2 have identical dynamics even though the dynamics are described by different expressions.

What next?

You can find more composite dynamical systems (including a multi-city SIR model and a cellular automata) in the documentation for AlgebraicDynamics.jl. As this package develops, a large class of open projects is to implement black and grey boxing functors between operad algebras, such as:

References

Baez and Pollard, 2017. A compositional framework for reaction networks. Reviews in Mathematical Physics. arXiv:1704.02051

Spivak, 2013. The operad of wiring diagrams: formalizing a graphical language for databases, recursion, and plug-and-play circuits. arXiv:1305.0297.