Composing open dynamical systems II: Undirected composition

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{Sch}(\mathsf{UWD}), which is defined by the schema below.¹

An undirected wiring diagram is a \mathsf{C}-Set over \mathsf{Sch}(\mathsf{UWD}), i.e. a functor d: \mathsf{Sch}(\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(\box)^{-1}(b) \subseteq P 3. a junction j identifies the ports d(\junc)^{-1}(j) \subseteq P and is exposed by the outer ports d(\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. Choose a composition pattern. Since our composition theory is \mathsf{Sch}(\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: - Adjust the parameters \alpha and \beta. - Change the primitive evaporation system so that it is also sinusoidal. Explore what happens when the periods align like in humid subtropics (rainy summer, dry winter) or anti-align like in Mediterranean climates (rainy winter, dry summer). - Suppose there is a person with a hose and a drain who is attempting to keep precisely 2 liters of water in the pool. Add a third resource sharer to the composite, which models this hose/drain system. Suggested dynamics are \dot w_h = \gamma(2 - w_h).

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,³ it lacks the clarity of the hierarchical approach which we can see from the flattened composition pattern: