ModelingToolkit Extract conserved quantities

General Usage
1

Hello,

I am interested in a chemical model given by:

using ModelingToolkit
using OrdinaryDiffEq

@parameters k1 k2 k3 k4 k5 k6 k7 k8 k9 k10 k11 M
@variables t O1D(t) O(t) O3(t) O2(t) NO(t) NO2(t) σ(t) r1(t) r2(t) r3(t) r4(t) r5(t) r6(t) r7(t) r8(t) r9(t) r10(t) r11(t)
D = Differential(t)

# Introduce forcing term
# Tr = 4.5
# Ts = 19.5

# σ_fun(tprime) = cos(tprime)#ifelse(mod(tprime/3600, 24) >= Tr &&  mod(t/3600, 24) >= Ts, 
# #                   0.5 + 0.5*cos(π*abs((2*mod(tprime/3600, 24) - Tr - Ts)/(Ts - Tr))*(2*mod(tprime/3600, 24) - Tr - Ts)/(Ts - Tr)),
# #                   0.0)

# @register_symbolic σ_fun(t)
# σ_fun(z)

eqs = [
    # Algebraic relations
       σ ~ cos(6*π*t),
       r1 ~ k1*O2*σ^3,
       r2 ~ k2*O*O2,
       r3 ~ k3*O3*σ,
       r4 ~ k4*O3*O,
       r5 ~ k5*O3*σ^2,
       r6 ~ k6*M*O1D,
       r7 ~ k7*O1D*O3,
       r8 ~ k8*O3*NO,
       r9 ~ k9*NO2*O,
       r10 ~ k10*NO2*σ,
       r11 ~ k11*NO*O,
    # Differential relations
       D(O1D) ~ r5 - r6 -r7,
       D(O) ~ 2*r1 - r2 + r3 - r4 + r6 - r9 + r10 -r11,
       D(O3) ~ r2 - r3 -r4 -r5 -r7 -r8,
       D(O2) ~ -r1 -r2 + r3 + 2*r4 + r5 + 2*r7 + r8 + r9,
       D(NO) ~ -r8 + r9 + r10 - r11,
       D(NO2) ~ r8 -r9 -r10 + r11]

@named symbolic_strato = ODESystem(eqs,t,
                      [O1D, O, O3, O2, NO, NO2, σ, r1, r2, r3, r4, r5, r6, r7, r8, r9, r10, r11],
                      [k1, k2, k3, k4, k5, k6, k7, k8, k9, k10, k11, M],tspan=(0, 10.0))


symbolic_strato = structural_simplify(symbolic_strato)
f_strato = ODEFunction(symbolic_strato)

prob = ODEProblem(f_strato, 1.0 .+ 0.1*randn(17),
    (0.0, 20.0), 0.1*rand(12))

Ut = Matrix([1 1 3 2 1 2;
         0 0 0 0 1 1]')

sol = solve(prob, Rosenbrock23())

In this case, there are two conserved quantities given by the conservation of the oxygen atoms and the nitrogen atoms. Thus U^\top \boldsymbol{x}(t) is a constant, with U^\top = \begin{bmatrix}1 & 1 & 3 & 2 & 1 & 2\\ 0 & 0 & 0 & 0 & 1 & 1] \end{bmatrix} and \boldsymbol{x} the state.

Can we use ModelingToolkit.jl to automatically identify U from the symbolic system of equations?

2

If \dot x = f(x) and U^Tx(t) = c, then

\begin{aligned} \dfrac{U^T x(t)}{dt} &= \dfrac{dc}{dt}\\ U^T f(x) &= 0 \end{aligned}

If you plug in a few different values for x in f(x) you’d get a matrix with a nullspace that matches U?

This method would lose the scaling of U and also the constant c. It appears to almost work. The first 4 rows of N below matches U up to a scaling, but the last two do not. The computed N does satisfy the same condition N^Tx = c_2 though.

I edited your code to not include the algebraic equations in the system of equations to not run a DAE solver when an ODE solver works fine. structural_simplify didn’t seem to simplify the algebraic equations away.

using ModelingToolkit
using OrdinaryDiffEq

@parameters k1 k2 k3 k4 k5 k6 k7 k8 k9 k10 k11 M
@variables t O1D(t) O(t) O3(t) O2(t) NO(t) NO2(t) #σ(t) r1(t) r2(t) r3(t) r4(t) r5(t) r6(t) r7(t) r8(t) r9(t) r10(t) r11(t)
D = Differential(t)

σ = cos(6*π*t)
r1 = k1*O2*σ^3
r2 = k2*O*O2
r3 = k3*O3*σ
r4 = k4*O3*O
r5 = k5*O3*σ^2
r6 = k6*M*O1D
r7 = k7*O1D*O3
r8 = k8*O3*NO
r9 = k9*NO2*O
r10 = k10*NO2*σ
r11 = k11*NO*O
# Differential relations
eqs = [
        (D(O1D) ~ r5 - r6 -r7)
       (D(O) ~ 2*r1 - r2 + r3 - r4 + r6 - r9 + r10 -r11)
       (D(O3) ~ r2 - r3 -r4 -r5 -r7 -r8)
       (D(O2) ~ -r1 -r2 + r3 + 2*r4 + r5 + 2*r7 + r8 + r9)
       (D(NO) ~ -r8 + r9 + r10 - r11)
       (D(NO2) ~ r8 -r9 -r10 + r11)
    ]

@named symbolic_strato = ODESystem(eqs,t,tspan=(0, 10.0))


symbolic_strato = structural_simplify(symbolic_strato)
f_strato = ODEFunction(symbolic_strato)

prob = ODEProblem(f_strato, 1.0 .+ 0.1*randn(6),
    (0.0, 20.0), 0.1*rand(12))

Ut = Matrix([1 1 3 2 1 2;
         0 0 0 0 1 1]')

sol = solve(prob, Rosenbrock23())

X = reduce(hcat, sol.u)
Ut'X

Xdots = reduce(hcat, f_strato.(sol.u, Ref(prob.p), sol.t))

using LinearAlgebra
N = nullspace(Xdots')
NtX = N'X

julia> 4N
6×2 Matrix{Float64}:
  1.00824   0.125311
  1.00824   0.125311
  3.02473   0.375932
  2.01649   0.250621
 -0.852972  2.74418
  0.155272  2.86949

julia> all(isapprox.(NtX[1,1]), NtX[1, :])
true

julia> all(isapprox.(NtX[2,1]), NtX[2, :])
true
1 Like
3

Thank you for your help. I was hoping to get something purely symbolic by variable elimination.

From the Catalyst’s documentation (FAQs · Catalyst.jl):

How to index solution objects using symbolic variables and observables?

One can directly use symbolic variables to index into SciML solution objects. Moreover, observables can also be evaluated in this way. For example, consider the system

using Catalyst, DifferentialEquations, Plots
rn = @reaction_network ABtoC begin
  (k₊,k₋), A + B <--> C
end

# initial condition and parameter values
setdefaults!(rn, [:A => 1.0, :B => 2.0, :C => 0.0, :k₊ => 1.0, :k₋ => 1.0])

Let’s convert it to a system of ODEs, using the conservation laws of the system to eliminate two of the species:

osys = convert(ODESystem, rn; remove_conserved = true)

\frac{dA(t)}{dt} =k_{−}(−A(t)+Γ_2)−k_{+}(A(t)+Γ_1)A(t)

Notice the resulting ODE system has just one ODE, while algebraic observables have been added for the two removed species (in terms of the conservation law constants, Γ[1] and Γ[2])

observed(osys)

\begin{aligned} B(t) & =A(t)+Γ_1\\ C(t) & =−A(t)+Γ_2\\ \end{aligned}​​​

It seems that Catalyst.j can extract the preserved quantities with the routine observed?