Operator Composition
There are a lot of cases where there are two different "processes" or "operators" that change the same variable. For example, CO2 in the atmosphere can be emitted by human activity, and it can also be absorbed by the ocean. In models, typically the emission and removal are considered separate processes which are represented by separate model components. However, when we want to combine these two components into a single model, we need to be able to compose them together.
We can use the operator_compose
function for this. It composes to systems of equations together by adding the right-hand side terms together of equations that have matching left-hand sides. The left hand sides of two equations will be considered matching if:
- They are both time derivatives of the same variable.
- The first one is a time derivative of a variable and the second one is the variable itself.
- There is an entry in the optional
translate
dictionary that maps the dependent variable in the first system to the dependent variable in the second system, e.g.Dict(sys1.sys.x => sys2.sys.y)
. - There is an entry in the optional
translate
dictionary that maps the dependent variable in the first system to the dependent variable in the second system, with a conversion factor, e.g.Dict(sys1.sys.x => sys2.sys.y => 6)
.
Perhaps we can make this somewhat clearer with some examples.
Examples
Example with matching variable time derivatives
The example below shows that when we operator_compose
two systems together that are both equal to D(x) = p
, the resulting system is equal to D(x) = 2p
.
using EarthSciMLBase
using ModelingToolkit
using ModelingToolkit: t_nounits, D_nounits
t = t_nounits
D = D_nounits
struct ExampleSysCoupler sys end
function ExampleSys()
@variables x(t)
@parameters p
ODESystem([D(x) ~ p], t; name=:ExampleSys,
metadata=Dict(:coupletype=>ExampleSysCoupler))
end
ExampleSys()
\[ \begin{align} \frac{\mathrm{d} x\left( t \right)}{\mathrm{d}t} &= p \end{align} \]
struct ExampleSys2Coupler sys end
function ExampleSys2()
@variables x(t)
@parameters p
ODESystem([D(x) ~ 2p], t; name=:ExampleSys2,
metadata=Dict(:coupletype=>ExampleSys2Coupler))
end
ExampleSys2()
\[ \begin{align} \frac{\mathrm{d} x\left( t \right)}{\mathrm{d}t} &= 2 p \end{align} \]
sys1 = ExampleSys()
sys2 = ExampleSys2()
function EarthSciMLBase.couple2(sys1::ExampleSysCoupler, sys2::ExampleSys2Coupler)
sys1, sys2 = sys1.sys, sys2.sys
operator_compose(sys1, sys2)
end
combined = couple(sys1, sys2)
combined_mtk = convert(ODESystem, combined)
\[ \begin{align} ExampleSys_{+}ExampleSys2_{ddt\_xˍt}\left( t \right) &= ExampleSys2_{+}ExampleSys2_{ddt\_xˍt}\left( t \right) \\ ExampleSys_{+}x\left( t \right) &= ExampleSys2_{+}x\left( t \right) \\ \frac{\mathrm{d} ExampleSys_{+}x\left( t \right)}{\mathrm{d}t} &= ExampleSys_{+}p + ExampleSys_{+}ExampleSys2_{ddt\_xˍt}\left( t \right) \\ ExampleSys2_{+}ExampleSys2_{ddt\_xˍt}\left( t \right) &= 2 ExampleSys2_{+}p \end{align} \]
The simplified equation should be D(x) = p + sys2_xˍt:
combined_simplified = structural_simplify(combined_mtk)
\[ \begin{align} \frac{\mathrm{d} ExampleSys_{+}x\left( t \right)}{\mathrm{d}t} &= ExampleSys_{+}p + ExampleSys_{+}ExampleSys2_{ddt\_xˍt}\left( t \right) \end{align} \]
where sys2_xˍt is also equal to p:
observed(combined_simplified)
\[ \begin{align} ExampleSys2_{+}ExampleSys2_{ddt\_xˍt}\left( t \right) &= 2 ExampleSys2_{+}p \\ ExampleSys2_{+}x\left( t \right) &= ExampleSys_{+}x\left( t \right) \\ ExampleSys_{+}ExampleSys2_{ddt\_xˍt}\left( t \right) &= ExampleSys2_{+}ExampleSys2_{ddt\_xˍt}\left( t \right) \end{align} \]
Example with non-matching variables
This example demonstrates a case where one variable in the first system is equal to another variable in the second system:
struct ExampleSys3Coupler sys end
function ExampleSys3()
@variables y(t)
@parameters p
ODESystem([D(y) ~ p], t; name=:ExampleSys3,
metadata=Dict(:coupletype=>ExampleSys3Coupler))
end
sys1 = ExampleSys()
sys2 = ExampleSys3()
function EarthSciMLBase.couple2(sys1::ExampleSysCoupler, sys2::ExampleSys3Coupler)
sys1, sys2 = sys1.sys, sys2.sys
operator_compose(sys1, sys2, Dict(sys1.x => sys2.y))
end
combined = couple(sys1, sys2)
combined_simplified = structural_simplify(convert(ODESystem, combined))
\[ \begin{align} \frac{\mathrm{d} ExampleSys_{+}x\left( t \right)}{\mathrm{d}t} &= ExampleSys_{+}p + ExampleSys_{+}ExampleSys3_{ddt\_yˍt}\left( t \right) \end{align} \]
observed(combined_simplified)
\[ \begin{align} ExampleSys3_{+}ExampleSys3_{ddt\_yˍt}\left( t \right) &= ExampleSys3_{+}p \\ ExampleSys3_{+}y\left( t \right) &= ExampleSys_{+}x\left( t \right) \\ ExampleSys_{+}ExampleSys3_{ddt\_yˍt}\left( t \right) &= ExampleSys3_{+}ExampleSys3_{ddt\_yˍt}\left( t \right) \end{align} \]
Example with a non-ODE system
In the second case above, we might have a variable in the second system that is equal to a rate, but it is not a time derivative. This could happen if we are extracting emissions from a file, and those emissions are already in units of kg/s, or something similar. The example below demonstrates this case. (Note that this case can also be used with the conversion factors shown in the last example.)
struct ExampleSysNonODECoupler sys end
function ExampleSysNonODE()
@variables y(t)
@parameters p
ODESystem([y ~ p], t; name=:ExampleSysNonODE,
metadata=Dict(:coupletype=>ExampleSysNonODECoupler))
end
sys1 = ExampleSys()
sys2 = ExampleSysNonODE()
function EarthSciMLBase.couple2(sys1::ExampleSysCoupler, sys2::ExampleSysNonODECoupler)
sys1, sys2 = sys1.sys, sys2.sys
operator_compose(sys1, sys2, Dict(sys1.x => sys2.y))
end
combined = couple(sys1, sys2)
sys_combined = structural_simplify(convert(ODESystem, combined))
\[ \begin{align} \frac{\mathrm{d} ExampleSys_{+}x\left( t \right)}{\mathrm{d}t} &= ExampleSys_{+}p + ExampleSys_{+}ExampleSysNonODE_{y}\left( t \right) \end{align} \]
observed(sys_combined)
\[ \begin{align} ExampleSysNonODE_{+}y\left( t \right) &= ExampleSysNonODE_{+}p \\ ExampleSys_{+}ExampleSysNonODE_{y}\left( t \right) &= ExampleSysNonODE_{+}y\left( t \right) \end{align} \]
Example with non-matching variables and a conversion factor
Finally, this last example shows the fourth case, where a conversion factor is included in the translation dictionary.
struct ExampleSysNonODE2Coupler sys end
function ExampleSysNonODE2()
@variables y(t)
@parameters p
ODESystem([y ~ p], t; name=:Docs₊ExampleSysNonODE2,
metadata=Dict(:coupletype=>ExampleSysNonODE2Coupler))
end
sys1 = ExampleSys()
sys2 = ExampleSysNonODE2()
function EarthSciMLBase.couple2(sys1::ExampleSysCoupler, sys2::ExampleSysNonODE2Coupler)
sys1, sys2 = sys1.sys, sys2.sys
operator_compose(sys1, sys2, Dict(sys1.x => sys2.y => 6.0))
end
combined = couple(sys1, sys2)
combined_simplified = structural_simplify(convert(ODESystem, combined))
\[ \begin{align} \frac{\mathrm{d} ExampleSys_{+}x\left( t \right)}{\mathrm{d}t} &= ExampleSys_{+}p + ExampleSys_{+}Docs_{+}ExampleSysNonODE2_{y}\left( t \right) \end{align} \]
observed(combined_simplified)
\[ \begin{align} Docs_{+}ExampleSysNonODE2_{+}y\left( t \right) &= Docs_{+}ExampleSysNonODE2_{+}p \\ ExampleSys_{+}Docs_{+}ExampleSysNonODE2_{y}\left( t \right) &= 6 Docs_{+}ExampleSysNonODE2_{+}y\left( t \right) \end{align} \]