BSM1
This example presents a more complex use of BioChemicalTreatment.jl for modeling water resource recovery facilities, historically called wastewater treatment plants, which can be considered close to real-world applications (see Usage for a simpler example). For this, it is presented how one could proceed to implement the so-called BSM1 model from IWA-MIA (Alex et al., 2008) or this link. There, the structure of this model is depicted as follows: 
Source: http://iwa-mia.org/wp-content/uploads/2019/04/bsm1_layout-1024x478.jpg
Specifically, the open-loop assessment desribed there is performed. See the document there for details. For all variables, where no value is given in the document, the steady-state values from the open-loop assessment is plugged in.
This example shows the open-loop assessment of the BSM1 described in the corresponding article (Alex et al., 2008), with two major differences:
- The rain weather influent is used instead of the dry weather influent: The influent thus consists of the 100 day stabilisation period followed by the rain weather influent.
- Instead of a clarifier using the Takács-Model, an ideal clarifier is used.
Plug and Play Example
using BioChemicalTreatment # Reactors etc.
using ModelingToolkit # Modeling Framework
using ModelingToolkitStandardLibrary.Blocks # Standard library providing several blocks
using HTTP # Load the influent from the BSM1 webpage
using DelimitedFiles # For reading the influent CSV
using DataInterpolations # For interpolating the influent
using DifferentialEquations # Solve the equations
using NonlinearSolve # Solve for initial conditions
using Plots # Plotting
# Bookkeeping variables
sys = []
eqns = Equation[]
##########
# Parameter definitions
##########
## Load the influent data
# Stabilisation period
stabilisation_inf = Dict(
["S_I", "S_S", "X_I", "X_S", "X_BH", "X_BA", "X_P", "S_O", "S_NO", "S_NH", "S_ND", "X_ND", "S_ALK", "q"] .=>
[30, 69.50, 51.2, 202.32, 28.17, 0, 0, 0, 0, 31.56, 6.95, 10.59, 7.0, 18446.0])
# Rain weather from webpage
rain_inf = readdlm(HTTP.get("http://iwa-mia.org/wp-content/uploads/2019/04/Inf_rain_2006.txt").body, header=true)
# Make variable naming match (separated by _ and all uppercase)
rain_inf[2][2:end-1] = uppercase.((x -> x[1]*"_"*x[2:end]).(rain_inf[2][2:end-1]))
# Make influent flowrate lowercase
rain_inf[2][end] = lowercase(rain_inf[2][end])
# Assemble the inflow
inf_t = [0, 100, (rain_inf[1][:, 1] .+ 100)...] # Time, first 100 days stabilisation, then rain inflow
# Values: Dict component -> vector of values (2x stabilization, then rain influent)
inf_val = Dict([comp => [(stabilisation_inf[comp]*[1, 1])..., rain_inf[1][:, i+1]...] for (i, comp) in enumerate(rain_inf[2][2:end])])
## Parameters
# ASM1 Parameters
params = Dict(
:Y_A => 0.24, # autotrophic yield
:Y_H => 0.67, # heterotrophic yield
:f_P => 0.08, # fraction of biomass yielding particulate procucts
:i_XB => 0.08, # mass N / mass COD in biomass
:i_XP => 0.06, # mass N / mass COD in products from biomass
:mu_H => 4.0, # Heterotrophic growth and decay
:K_S => 10.0,
:K_OH => 0.2,
:K_NO => 0.5,
:b_H => 0.3,
:eta_g => 0.8, # correction factor for anoxic growth of heterotrophs
:eta_h => 0.8, # correction factor for anoxic hydrolysis
:kappa_h => 3.0, # hydrolysis
:K_X => 0.1,
:mu_A => 0.5, # Autotrophic growth and decay
:K_NH => 1.0,
:b_A => 0.05,
:K_OA => 0.4,
:kappa_a => 0.05, # ammonification
)
params_aeration = Dict([
:S_O_max => 8 # Oxygen saturation concentration
])
# Volumes
V_na = 1000 # Volume [m^3]
V_a = 1333 # Volume [m^3]
# Aeration Specification (Reactors 3-4, fixed)
KLa = 240 # Aeration Rate [1/d]
# Control variables (Open-Loop assessment)
Qa = 55338 # Flow [m^3/d]
KLa5 = 84 # Aeration rate [1/d]
# Initial conditions
c_init = Dict(
:S_I => 1.0,
:S_S => 1.0,
:X_I => 1.0,
:X_S => 1.0,
:X_BH => 1.0,
:X_BA => 1.0,
:X_P => 1.0,
:S_O => 1.0,
:S_NO => 1.0,
:S_NH => 1.0,
:S_ND => 1.0,
:X_ND => 1.0,
:S_ALK => 1.0
)
##########
# Model Components definitions
##########
# Start by selecting the process model to be used. Here, ASM1, as default process with the parameters defined above
@set_process asm1 = ASM1(;params...)
# Then, we build up the influent.
# For this we use a linear interpolation between the loaded timeseries
interp_fct(y, x) = t -> LinearInterpolation(y, x, extrapolation=ExtrapolationType.Constant)(t) # Generator to build the linear interpolation function
@named influent = Influent(
[k => interp_fct(inf_val[k], inf_t) for k in keys(inf_val)], # The Interpolation function for the corresponding value
)
push!(sys, influent)
# Then define all other needed flow elements
# Unifier for inflow, internal and external recycle
@named unifier = FlowUnifier(;n_in=3)
push!(sys, unifier)
# Separator and pump for internal recycle
@named internal_rec_sep = FlowSeparator()
@named internal_rec_pump = FlowPump(;flowrate=Qa)
append!(sys, [internal_rec_sep, internal_rec_pump])
# Clarifier
@named clarifier = IdealClarifier()
@named clarifier_pump = FlowPump(flowrate=18446.0 + 385.0) # Open-Loop steady-state flow rate from Alex et al. (2018)
append!(sys, [clarifier, clarifier_pump])
# Wastage
@named wastage_separator = FlowSeparator(n_out=2)
@named wastage_pump = FlowPump(flowrate = 385) # Open-Loop steady-state flow rate from Alex et al. (2018)
append!(sys, [wastage_separator, wastage_pump])
# Define the reactors, first two are non-aerated
# First Reactor, non-aerated CSTR with ASM1
@named reactor1 = CSTR(V_na, initial_states = c_init) # CSTR with volume V_na, initial states from c_init
push!(sys, reactor1)
# Second Reactor, non-aerated CSTR with ASM1
@named reactor2 = CSTR(V_na, initial_states = c_init)
push!(sys, reactor2)
# Third Reactor, aerated CSTR with ASM1
@named reactor3 = CSTR(V_a, Aeration(;params_aeration..., name = :reactor3_aer), initial_states = c_init)
@named reactor3_kla = Constant(;k = KLa) # Constant aeration. Replace with controller if added.
append!(sys, [reactor3, reactor3_kla])
push!(eqns, connect(reactor3_kla.output, exogenous_inputs(reactor3, :reactor3_aer₊k_La))) # Connect constant aeration to reactor
# Fourth Reactor, aerated CSTR with ASM1
@named reactor4 = CSTR(V_a, Aeration(;params_aeration..., name = :reactor4_aer), initial_states = c_init)
@named reactor4_kla = Constant(;k = KLa) # Constant aeration. Replace with controller if added.
append!(sys, [reactor4, reactor4_kla])
push!(eqns, connect(reactor4_kla.output, exogenous_inputs(reactor4, :reactor4_aer₊k_La))) # Connect constant aeration to reactor
# Fifth Reactor, aerated CSTR with ASM1, k_La for fifth reactor and open-loop verification
@named reactor5 = CSTR(V_a, Aeration(;params_aeration..., name = :reactor5_aer), initial_states = c_init)
@named reactor5_kla = Constant(;k = KLa) # Constant aeration. Replace with controller if added.
append!(sys, [reactor5, reactor5_kla])
push!(eqns, connect(reactor5_kla.output, exogenous_inputs(reactor5, :reactor5_aer₊k_La))) # Connect constant aeration to reactor
##########
# Connections
##########
## Connect the systems
append!(eqns, [
# Connect the forward flows and add all reactor connection equations
# Connect the influent to the unifier first port
connect(outflows(influent)[1], inflows(unifier)[1]); # Yes the influent has an outflow. Might be unintuitive, but it is flowing out from it...
# Connect Unifier outflow to first reactor
connect(outflows(unifier)[1], inflows(reactor1)[1]);
# Connect reactor 1 to reactor 2
connect(outflows(reactor1)[1], inflows(reactor2)[1]);
# Connect reactor 2 to reactor 3
connect(outflows(reactor2)[1], inflows(reactor3)[1]);
# Connect reactor 3 to reactor 4
connect(outflows(reactor3)[1], inflows(reactor4)[1]);
# Connect reactor 4 to reactor 5
connect(outflows(reactor4)[1], inflows(reactor5)[1]);
# Connect reactor 5 to separator for internal recirculation
connect(outflows(reactor5)[1], inflows(internal_rec_sep)[1]);
# Connect internal recirculation separator to clarifier
connect(outflows(internal_rec_sep)[1], inflows(clarifier)[1]);
# Connect the recycles
# External recirculation
# Connect clarifier to clarifier pump
connect(outflows(clarifier, :sludge), inflows(clarifier_pump)[1]);
# Connect clarifier pump to wastage separator
connect(outflows(clarifier_pump)[1], inflows(wastage_separator)[1])
# Connect wastage separator outflow 1 to unifier inflow 2
connect(outflows(wastage_separator)[1], inflows(unifier)[2]);
# Connect wastage separator outflow 2 to wastage pump
connect(outflows(wastage_separator)[2], inflows(wastage_pump)[1])
# Internal recycle
# Connect internal rec separator to its pump
connect(outflows(internal_rec_sep)[2], inflows(internal_rec_pump)[1])
# Connect internal rec pump back to first unifier inflow 3
connect(outflows(internal_rec_pump)[1], inflows(unifier)[3])
])
##########
# Assembling overall model
##########
@named bsm1_model = ODESystem(
eqns, t, systems = sys
)
bsm1_simplified = structural_simplify(bsm1_model)
##########
# Simulation and Plotting
##########
## Generate the problem to solve (114 days, stabilization + 2 weeks of influent)
t_sim = (0, 114)
prob = ODEProblem(bsm1_simplified, [], t_sim, [], initializealg=ShampineCollocationInit()) # We need to explicitly select the initialization algorithm here. The default doesn't work
# Solve the problem. Use a DAE solver (Rosenbrock methods, see solver docs for more infos on the method)
sol = solve(prob, Rodas5P())
## Plot the output
plot(sol,
xlim = (95, 114), #5 last days of stabilisation and then the influent file
idxs = unknowns(bsm1_model)[occursin.("clarifier₊effluent", string.(unknowns(bsm1_model)))][2:end], # States which are associated with effluent but not flow
label = permutedims(string.(unknowns(outflows(clarifier, :effluent))[2:end])),
legend = :topleft,
title = "BSM1 Effluent",
)
Step-by-Step Explanations
Setup
First, some packages have to be imported:
using BioChemicalTreatment # Reactors etc.
using ModelingToolkit # Modeling Framework
using ModelingToolkitStandardLibrary.Blocks # Standard library providing several blocks
using HTTP # Load the influent from the BSM1 webpage
using DelimitedFiles # For reading the influent CSV
using DataInterpolations # For interpolating the influent
using DifferentialEquations # Solve the equations
using NonlinearSolve # Solve for initial conditions
using Plots # Plotting
# Bookkeeping variables
sys = []
eqns = Equation[]The next step is to load data and define the needed parameters
## Load the influent data
# Stabilisation period
stabilisation_inf = Dict(
["S_I", "S_S", "X_I", "X_S", "X_BH", "X_BA", "X_P", "S_O", "S_NO", "S_NH", "S_ND", "X_ND", "S_ALK", "q"] .=>
[30, 69.50, 51.2, 202.32, 28.17, 0, 0, 0, 0, 31.56, 6.95, 10.59, 7.0, 18446.0])
# Rain weather from webpage
rain_inf = readdlm(HTTP.get("http://iwa-mia.org/wp-content/uploads/2019/04/Inf_rain_2006.txt").body, header=true)
# Make variable naming match (separated by _ and all uppercase)
rain_inf[2][2:end-1] = uppercase.((x -> x[1]*"_"*x[2:end]).(rain_inf[2][2:end-1]))
# Make influent flowrate lowercase
rain_inf[2][end] = lowercase(rain_inf[2][end])
# Assemble the inflow
inf_t = [0, 100, (rain_inf[1][:, 1] .+ 100)...] # Time, first 100 days stabilisation, then rain inflow
# Values: Dict component -> vector of values (2x stabilization, then rain influent)
inf_val = Dict([comp => [(stabilisation_inf[comp]*[1, 1])..., rain_inf[1][:, i+1]...] for (i, comp) in enumerate(rain_inf[2][2:end])])
## Parameters
# ASM1 Parameters
params = Dict(
:Y_A => 0.24, # autotrophic yield
:Y_H => 0.67, # heterotrophic yield
:f_P => 0.08, # fraction of biomass yielding particulate procucts
:i_XB => 0.08, # mass N / mass COD in biomass
:i_XP => 0.06, # mass N / mass COD in products from biomass
:mu_H => 4.0, # Heterotrophic growth and decay
:K_S => 10.0,
:K_OH => 0.2,
:K_NO => 0.5,
:b_H => 0.3,
:eta_g => 0.8, # correction factor for anoxic growth of heterotrophs
:eta_h => 0.8, # correction factor for anoxic hydrolysis
:kappa_h => 3.0, # hydrolysis
:K_X => 0.1,
:mu_A => 0.5, # Autotrophic growth and decay
:K_NH => 1.0,
:b_A => 0.05,
:K_OA => 0.4,
:kappa_a => 0.05, # ammonification
)
params_aeration = Dict([
:S_O_max => 8 # Oxygen saturation concentration
])
# Volumes
V_na = 1000 # Volume [m^3]
V_a = 1333 # Volume [m^3]
# Aeration Specification (Reactors 3-4, fixed)
KLa = 240 # Aeration Rate [1/d]
# Control variables (Open-Loop assessment)
Qa = 55338 # Flow [m^3/d]
KLa5 = 84 # Aeration rate [1/d]
# Initial conditions
c_init = Dict(
:S_I => 1.0,
:S_S => 1.0,
:X_I => 1.0,
:X_S => 1.0,
:X_BH => 1.0,
:X_BA => 1.0,
:X_P => 1.0,
:S_O => 1.0,
:S_NO => 1.0,
:S_NH => 1.0,
:S_ND => 1.0,
:X_ND => 1.0,
:S_ALK => 1.0
)System Creation
Once the data is loaded and the parameters are defined, we start by creating the individual components, which later will be connected to form the overall BSM1 model. Before we come to the reactors, we set the influent and the static flow elements (such as flow splitters and unifiers) up.
There are currently two ways how to call a model in BioChemicalTreatment.jl: First the equation-defined models loaded using the corresponding function (this example uses this way) or through the models that are loaded over the Reactions folder (currently the matrix-defined models). You will find that ASM1 is present at both locations (called as ASM1() vs Process("ASM1")). There is, however, a small difference. The equation-defined model in the ProcessSimulator (ASM1()) is the original published by (Henze et al., 1987) where the mass balance does not come out even, therefore we implemented in the Reactions folder the updated version by (Hauduc et al., 2010) (Process("ASM1")). The idea of the BSM1 model is to compare modelling results with a benchmark. Therefore, we use here the older version, as that is the version used for BSM1.
## Build up model parts
# Start by selecting the process model to be used. Here, ASM1, as default process with the parameters defined abouve
@set_process asm1 = ASM1(;params...)
# Then, we build up the influent.
# For this we use a linear interpolation between the loaded timeseries
interp_fct(y, x) = t -> LinearInterpolation(y, x, extrapolation=ExtrapolationType.Constant)(t) # Generator to build the linear interpolation function
@named influent = Influent(
[k => interp_fct(inf_val[k], inf_t) for k in keys(inf_val)], # The Interpolation function for the corresponding value
)
push!(sys, influent)
# Then define all other needed flow elements
# Unifier for inflow, internal and external recycle
@named unifier = FlowUnifier(;n_in=3)
push!(sys, unifier)
# Separator and pump for internal recycle
@named internal_rec_sep = FlowSeparator()
@named internal_rec_pump = FlowPump(;flowrate=Qa)
append!(sys, [internal_rec_sep, internal_rec_pump])
# Clarifier
@named clarifier = IdealClarifier()
@named clarifier_pump = FlowPump(flowrate=18446.0 + 385.0) # Open-Loop steady-state flow rate from Alex et al. (2018)
append!(sys, [clarifier, clarifier_pump])
# Wastage
@named wastage_separator = FlowSeparator(n_out=2)
@named wastage_pump = FlowPump(flowrate = 385) # Open-Loop steady-state flow rate from Alex et al. (2018)
append!(sys, [wastage_separator, wastage_pump])
nameof.(sys)Then, the continuously stirred tank reactors (CSTR) are set up:
The first two reactors are non-aerated, thus the only added process is ASM1, which is the set default. Thus it is easy to :
# Define the reactors, first two are non-aerated
# First Reactor, non-aerated CSTR with ASM1
@named reactor1 = CSTR(V_na, initial_states = c_init) # CSTR with volume V_na, initial states from c_init
push!(sys, reactor1)
# Second Reactor, non-aerated CSTR with ASM1
@named reactor2 = CSTR(V_na, initial_states = c_init)
push!(sys, reactor2)
nameof.(sys)The third, fourth and fifth reactors are aerated. Thus the aeration process is added to the ASM1. As the oxygen concentration in both cases is called S_O, this is done very conveniently:
# Third Reactor, aerated CSTR with ASM1
@named reactor3 = CSTR(V_a, Aeration(;params_aeration..., name = :reactor3_aer), initial_states = c_init)
@named reactor3_kla = Constant(;k = KLa) # Constant aeration. Replace with controller if added.
append!(sys, [reactor3, reactor3_kla])
push!(eqns, connect(reactor3_kla.output, exogenous_inputs(reactor3, :reactor3_aer₊k_La))) # Connect constant aeration to reactor
# Fourth Reactor, aerated CSTR with ASM1
@named reactor4 = CSTR(V_a, Aeration(;params_aeration..., name = :reactor4_aer), initial_states = c_init)
@named reactor4_kla = Constant(;k = KLa) # Constant aeration. Replace with controller if added.
append!(sys, [reactor4, reactor4_kla])
push!(eqns, connect(reactor4_kla.output, exogenous_inputs(reactor4, :reactor4_aer₊k_La))) # Connect constant aeration to reactor
# Fifth Reactor, aerated CSTR with ASM1, k_La for fifth reactor and open-loop verification
@named reactor5 = CSTR(V_a, Aeration(;params_aeration..., name = :reactor5_aer), initial_states = c_init)
@named reactor5_kla = Constant(;k = KLa) # Constant aeration. Replace with controller if added.
append!(sys, [reactor5, reactor5_kla])
push!(eqns, connect(reactor5_kla.output, exogenous_inputs(reactor5, :reactor5_aer₊k_La))) # Connect constant aeration to reactor
# output
3-element Vector{Equation}:
connect(reactor3_kla.output, reactor3.reactor3_aer.k_La)
connect(reactor4_kla.output, reactor4.reactor4_aer.k_La)
connect(reactor5_kla.output, reactor5.reactor5_aer.k_La)And then we created all individual components and can start to connect them.
System Connection
To connect the systems the connection equations can be created using the connect, inflows and outflows functions:
## Connect the systems
append!(eqns, [
# Connect the forward flows and add all reactor connection equations
# Connect the influent to the unifier first port
connect(outflows(influent)[1], inflows(unifier)[1]); # Yes the influent has an outflow. Might be unintuitive, but it is flowing out from it...
# Connect Unifier outflow to first reactor
connect(outflows(unifier)[1], inflows(reactor1)[1]);
# Connect reactor 1 to reactor 2
connect(outflows(reactor1)[1], inflows(reactor2)[1]);
# Connect reactor 2 to reactor 3
connect(outflows(reactor2)[1], inflows(reactor3)[1]);
# Connect reactor 3 to reactor 4
connect(outflows(reactor3)[1], inflows(reactor4)[1]);
# Connect reactor 4 to reactor 5
connect(outflows(reactor4)[1], inflows(reactor5)[1]);
# Connect reactor 5 to separator for internal recirculation
connect(outflows(reactor5)[1], inflows(internal_rec_sep)[1]);
# Connect internal recirculation separator to clarifier
connect(outflows(internal_rec_sep)[1], inflows(clarifier)[1]);
# Connect the recycles
# External recirculation
# Connect clarifier to clarifier pump
connect(outflows(clarifier, :sludge), inflows(clarifier_pump)[1]);
# Connect clarifier pump to wastage separator
connect(outflows(clarifier_pump)[1], inflows(wastage_separator)[1])
# Connect wastage separator outflow 1 to unifier inflow 2
connect(outflows(wastage_separator)[1], inflows(unifier)[2]);
# Connect wastage separator outflow 2 to wastage pump
connect(outflows(wastage_separator)[2], inflows(wastage_pump)[1])
# Internal recycle
# Connect internal rec separator to its pump
connect(outflows(internal_rec_sep)[2], inflows(internal_rec_pump)[1])
# Connect internal rec pump back to first unifier inflow 3
connect(outflows(internal_rec_pump)[1], inflows(unifier)[3])
])Then, the overall BSM1 system can be built and the equations simplified using ModelingToolkit.jl:
@named bsm1_model = ODESystem(
eqns, t, systems = sys
)
bsm1_simplified = structural_simplify(bsm1_model)@named bsm1_model = ODESystem(
eqns, t, systems = sys
)
bsm1_simplified = structural_simplify(bsm1_model)\[ \begin{align} \frac{\mathrm{d} \mathtt{reactor1.reactor1\_reactor.states.S\_ALK}\left( t \right)}{\mathrm{d}t} &= \frac{ - \mathtt{reactor1.reactor1\_reactor.outflow.q}\left( t \right) \mathtt{reactor1.reactor1\_reactor.states.S\_ALK}\left( t \right)}{\mathtt{reactor1.reactor1\_reactor.V}} + \frac{\mathtt{reactor1.reactor1\_reactor.inflow.S\_ALK}\left( t \right) \mathtt{reactor1.reactor1\_reactor.inflow.q}\left( t \right)}{\mathtt{reactor1.reactor1\_reactor.V}} + \mathtt{reactor1.reactor1\_reactor.rates.S\_ALK}\left( t \right) \\ \frac{\mathrm{d} \mathtt{reactor1.reactor1\_reactor.states.S\_I}\left( t \right)}{\mathrm{d}t} &= \frac{ - \mathtt{reactor1.reactor1\_reactor.states.S\_I}\left( t \right) \mathtt{reactor1.reactor1\_reactor.outflow.q}\left( t \right)}{\mathtt{reactor1.reactor1\_reactor.V}} + \frac{\mathtt{reactor1.reactor1\_reactor.inflow.S\_I}\left( t \right) \mathtt{reactor1.reactor1\_reactor.inflow.q}\left( t \right)}{\mathtt{reactor1.reactor1\_reactor.V}} + \mathtt{reactor1.reactor1\_reactor.rates.S\_I}\left( t \right) \\ \frac{\mathrm{d} \mathtt{reactor1.reactor1\_reactor.states.S\_ND}\left( t \right)}{\mathrm{d}t} &= \mathtt{reactor1.reactor1\_reactor.rates.S\_ND}\left( t \right) + \frac{ - \mathtt{reactor1.reactor1\_reactor.outflow.q}\left( t \right) \mathtt{reactor1.reactor1\_reactor.states.S\_ND}\left( t \right)}{\mathtt{reactor1.reactor1\_reactor.V}} + \frac{\mathtt{reactor1.reactor1\_reactor.inflow.S\_ND}\left( t \right) \mathtt{reactor1.reactor1\_reactor.inflow.q}\left( t \right)}{\mathtt{reactor1.reactor1\_reactor.V}} \\ \frac{\mathrm{d} \mathtt{reactor1.reactor1\_reactor.states.S\_NH}\left( t \right)}{\mathrm{d}t} &= \frac{\mathtt{reactor1.reactor1\_reactor.inflow.S\_NH}\left( t \right) \mathtt{reactor1.reactor1\_reactor.inflow.q}\left( t \right)}{\mathtt{reactor1.reactor1\_reactor.V}} + \frac{ - \mathtt{reactor1.reactor1\_reactor.outflow.q}\left( t \right) \mathtt{reactor1.reactor1\_reactor.states.S\_NH}\left( t \right)}{\mathtt{reactor1.reactor1\_reactor.V}} + \mathtt{reactor1.reactor1\_reactor.rates.S\_NH}\left( t \right) \\ \frac{\mathrm{d} \mathtt{reactor1.reactor1\_reactor.states.S\_NO}\left( t \right)}{\mathrm{d}t} &= \mathtt{reactor1.reactor1\_reactor.rates.S\_NO}\left( t \right) + \frac{\mathtt{reactor1.reactor1\_reactor.inflow.S\_NO}\left( t \right) \mathtt{reactor1.reactor1\_reactor.inflow.q}\left( t \right)}{\mathtt{reactor1.reactor1\_reactor.V}} + \frac{ - \mathtt{reactor1.reactor1\_reactor.states.S\_NO}\left( t \right) \mathtt{reactor1.reactor1\_reactor.outflow.q}\left( t \right)}{\mathtt{reactor1.reactor1\_reactor.V}} \\ \frac{\mathrm{d} \mathtt{reactor1.reactor1\_reactor.states.S\_O}\left( t \right)}{\mathrm{d}t} &= \mathtt{reactor1.reactor1\_reactor.rates.S\_O}\left( t \right) + \frac{ - \mathtt{reactor1.reactor1\_reactor.outflow.q}\left( t \right) \mathtt{reactor1.reactor1\_reactor.states.S\_O}\left( t \right)}{\mathtt{reactor1.reactor1\_reactor.V}} + \frac{\mathtt{reactor1.reactor1\_reactor.inflow.S\_O}\left( t \right) \mathtt{reactor1.reactor1\_reactor.inflow.q}\left( t \right)}{\mathtt{reactor1.reactor1\_reactor.V}} \\ \frac{\mathrm{d} \mathtt{reactor1.reactor1\_reactor.states.S\_S}\left( t \right)}{\mathrm{d}t} &= \mathtt{reactor1.reactor1\_reactor.rates.S\_S}\left( t \right) + \frac{\mathtt{reactor1.reactor1\_reactor.inflow.S\_S}\left( t \right) \mathtt{reactor1.reactor1\_reactor.inflow.q}\left( t \right)}{\mathtt{reactor1.reactor1\_reactor.V}} + \frac{ - \mathtt{reactor1.reactor1\_reactor.outflow.q}\left( t \right) \mathtt{reactor1.reactor1\_reactor.states.S\_S}\left( t \right)}{\mathtt{reactor1.reactor1\_reactor.V}} \\ \frac{\mathrm{d} \mathtt{reactor1.reactor1\_reactor.states.X\_BA}\left( t \right)}{\mathrm{d}t} &= \frac{ - \mathtt{reactor1.reactor1\_reactor.outflow.q}\left( t \right) \mathtt{reactor1.reactor1\_reactor.states.X\_BA}\left( t \right)}{\mathtt{reactor1.reactor1\_reactor.V}} + \mathtt{reactor1.reactor1\_reactor.rates.X\_BA}\left( t \right) + \frac{\mathtt{reactor1.reactor1\_reactor.inflow.X\_BA}\left( t \right) \mathtt{reactor1.reactor1\_reactor.inflow.q}\left( t \right)}{\mathtt{reactor1.reactor1\_reactor.V}} \\ \frac{\mathrm{d} \mathtt{reactor1.reactor1\_reactor.states.X\_BH}\left( t \right)}{\mathrm{d}t} &= \mathtt{reactor1.reactor1\_reactor.rates.X\_BH}\left( t \right) + \frac{ - \mathtt{reactor1.reactor1\_reactor.states.X\_BH}\left( t \right) \mathtt{reactor1.reactor1\_reactor.outflow.q}\left( t \right)}{\mathtt{reactor1.reactor1\_reactor.V}} + \frac{\mathtt{reactor1.reactor1\_reactor.inflow.q}\left( t \right) \mathtt{reactor1.reactor1\_reactor.inflow.X\_BH}\left( t \right)}{\mathtt{reactor1.reactor1\_reactor.V}} \\ \frac{\mathrm{d} \mathtt{reactor1.reactor1\_reactor.states.X\_I}\left( t \right)}{\mathrm{d}t} &= \frac{ - \mathtt{reactor1.reactor1\_reactor.states.X\_I}\left( t \right) \mathtt{reactor1.reactor1\_reactor.outflow.q}\left( t \right)}{\mathtt{reactor1.reactor1\_reactor.V}} + \mathtt{reactor1.reactor1\_reactor.rates.X\_I}\left( t \right) + \frac{\mathtt{reactor1.reactor1\_reactor.inflow.q}\left( t \right) \mathtt{reactor1.reactor1\_reactor.inflow.X\_I}\left( t \right)}{\mathtt{reactor1.reactor1\_reactor.V}} \\ \frac{\mathrm{d} \mathtt{reactor1.reactor1\_reactor.states.X\_ND}\left( t \right)}{\mathrm{d}t} &= \mathtt{reactor1.reactor1\_reactor.rates.X\_ND}\left( t \right) + \frac{ - \mathtt{reactor1.reactor1\_reactor.outflow.q}\left( t \right) \mathtt{reactor1.reactor1\_reactor.states.X\_ND}\left( t \right)}{\mathtt{reactor1.reactor1\_reactor.V}} + \frac{\mathtt{reactor1.reactor1\_reactor.inflow.q}\left( t \right) \mathtt{reactor1.reactor1\_reactor.inflow.X\_ND}\left( t \right)}{\mathtt{reactor1.reactor1\_reactor.V}} \\ \frac{\mathrm{d} \mathtt{reactor1.reactor1\_reactor.states.X\_P}\left( t \right)}{\mathrm{d}t} &= \frac{\mathtt{reactor1.reactor1\_reactor.inflow.X\_P}\left( t \right) \mathtt{reactor1.reactor1\_reactor.inflow.q}\left( t \right)}{\mathtt{reactor1.reactor1\_reactor.V}} + \frac{ - \mathtt{reactor1.reactor1\_reactor.states.X\_P}\left( t \right) \mathtt{reactor1.reactor1\_reactor.outflow.q}\left( t \right)}{\mathtt{reactor1.reactor1\_reactor.V}} + \mathtt{reactor1.reactor1\_reactor.rates.X\_P}\left( t \right) \\ \frac{\mathrm{d} \mathtt{reactor1.reactor1\_reactor.states.X\_S}\left( t \right)}{\mathrm{d}t} &= \frac{ - \mathtt{reactor1.reactor1\_reactor.outflow.q}\left( t \right) \mathtt{reactor1.reactor1\_reactor.states.X\_S}\left( t \right)}{\mathtt{reactor1.reactor1\_reactor.V}} + \mathtt{reactor1.reactor1\_reactor.rates.X\_S}\left( t \right) + \frac{\mathtt{reactor1.reactor1\_reactor.inflow.X\_S}\left( t \right) \mathtt{reactor1.reactor1\_reactor.inflow.q}\left( t \right)}{\mathtt{reactor1.reactor1\_reactor.V}} \\ \frac{\mathrm{d} \mathtt{reactor2.reactor2\_reactor.states.S\_ALK}\left( t \right)}{\mathrm{d}t} &= \frac{ - \mathtt{reactor2.reactor2\_reactor.states.S\_ALK}\left( t \right) \mathtt{reactor2.reactor2\_reactor.outflow.q}\left( t \right)}{\mathtt{reactor2.reactor2\_reactor.V}} + \mathtt{reactor2.reactor2\_reactor.rates.S\_ALK}\left( t \right) + \frac{\mathtt{reactor2.reactor2\_reactor.inflow.S\_ALK}\left( t \right) \mathtt{reactor2.reactor2\_reactor.inflow.q}\left( t \right)}{\mathtt{reactor2.reactor2\_reactor.V}} \\ \frac{\mathrm{d} \mathtt{reactor2.reactor2\_reactor.states.S\_I}\left( t \right)}{\mathrm{d}t} &= \frac{\mathtt{reactor2.reactor2\_reactor.inflow.q}\left( t \right) \mathtt{reactor2.reactor2\_reactor.inflow.S\_I}\left( t \right)}{\mathtt{reactor2.reactor2\_reactor.V}} + \frac{ - \mathtt{reactor2.reactor2\_reactor.states.S\_I}\left( t \right) \mathtt{reactor2.reactor2\_reactor.outflow.q}\left( t \right)}{\mathtt{reactor2.reactor2\_reactor.V}} + \mathtt{reactor2.reactor2\_reactor.rates.S\_I}\left( t \right) \\ \frac{\mathrm{d} \mathtt{reactor2.reactor2\_reactor.states.S\_ND}\left( t \right)}{\mathrm{d}t} &= \frac{\mathtt{reactor2.reactor2\_reactor.inflow.q}\left( t \right) \mathtt{reactor2.reactor2\_reactor.inflow.S\_ND}\left( t \right)}{\mathtt{reactor2.reactor2\_reactor.V}} + \mathtt{reactor2.reactor2\_reactor.rates.S\_ND}\left( t \right) + \frac{ - \mathtt{reactor2.reactor2\_reactor.states.S\_ND}\left( t \right) \mathtt{reactor2.reactor2\_reactor.outflow.q}\left( t \right)}{\mathtt{reactor2.reactor2\_reactor.V}} \\ \frac{\mathrm{d} \mathtt{reactor2.reactor2\_reactor.states.S\_NH}\left( t \right)}{\mathrm{d}t} &= \mathtt{reactor2.reactor2\_reactor.rates.S\_NH}\left( t \right) + \frac{ - \mathtt{reactor2.reactor2\_reactor.outflow.q}\left( t \right) \mathtt{reactor2.reactor2\_reactor.states.S\_NH}\left( t \right)}{\mathtt{reactor2.reactor2\_reactor.V}} + \frac{\mathtt{reactor2.reactor2\_reactor.inflow.S\_NH}\left( t \right) \mathtt{reactor2.reactor2\_reactor.inflow.q}\left( t \right)}{\mathtt{reactor2.reactor2\_reactor.V}} \\ \frac{\mathrm{d} \mathtt{reactor2.reactor2\_reactor.states.S\_NO}\left( t \right)}{\mathrm{d}t} &= \mathtt{reactor2.reactor2\_reactor.rates.S\_NO}\left( t \right) + \frac{ - \mathtt{reactor2.reactor2\_reactor.states.S\_NO}\left( t \right) \mathtt{reactor2.reactor2\_reactor.outflow.q}\left( t \right)}{\mathtt{reactor2.reactor2\_reactor.V}} + \frac{\mathtt{reactor2.reactor2\_reactor.inflow.q}\left( t \right) \mathtt{reactor2.reactor2\_reactor.inflow.S\_NO}\left( t \right)}{\mathtt{reactor2.reactor2\_reactor.V}} \\ \frac{\mathrm{d} \mathtt{reactor2.reactor2\_reactor.states.S\_O}\left( t \right)}{\mathrm{d}t} &= \frac{ - \mathtt{reactor2.reactor2\_reactor.outflow.q}\left( t \right) \mathtt{reactor2.reactor2\_reactor.states.S\_O}\left( t \right)}{\mathtt{reactor2.reactor2\_reactor.V}} + \frac{\mathtt{reactor2.reactor2\_reactor.inflow.q}\left( t \right) \mathtt{reactor2.reactor2\_reactor.inflow.S\_O}\left( t \right)}{\mathtt{reactor2.reactor2\_reactor.V}} + \mathtt{reactor2.reactor2\_reactor.rates.S\_O}\left( t \right) \\ \frac{\mathrm{d} \mathtt{reactor2.reactor2\_reactor.states.S\_S}\left( t \right)}{\mathrm{d}t} &= \frac{\mathtt{reactor2.reactor2\_reactor.inflow.q}\left( t \right) \mathtt{reactor2.reactor2\_reactor.inflow.S\_S}\left( t \right)}{\mathtt{reactor2.reactor2\_reactor.V}} + \mathtt{reactor2.reactor2\_reactor.rates.S\_S}\left( t \right) + \frac{ - \mathtt{reactor2.reactor2\_reactor.outflow.q}\left( t \right) \mathtt{reactor2.reactor2\_reactor.states.S\_S}\left( t \right)}{\mathtt{reactor2.reactor2\_reactor.V}} \\ \frac{\mathrm{d} \mathtt{reactor2.reactor2\_reactor.states.X\_BA}\left( t \right)}{\mathrm{d}t} &= \frac{\mathtt{reactor2.reactor2\_reactor.inflow.X\_BA}\left( t \right) \mathtt{reactor2.reactor2\_reactor.inflow.q}\left( t \right)}{\mathtt{reactor2.reactor2\_reactor.V}} + \mathtt{reactor2.reactor2\_reactor.rates.X\_BA}\left( t \right) + \frac{ - \mathtt{reactor2.reactor2\_reactor.states.X\_BA}\left( t \right) \mathtt{reactor2.reactor2\_reactor.outflow.q}\left( t \right)}{\mathtt{reactor2.reactor2\_reactor.V}} \\ \frac{\mathrm{d} \mathtt{reactor2.reactor2\_reactor.states.X\_BH}\left( t \right)}{\mathrm{d}t} &= \frac{\mathtt{reactor2.reactor2\_reactor.inflow.X\_BH}\left( t \right) \mathtt{reactor2.reactor2\_reactor.inflow.q}\left( t \right)}{\mathtt{reactor2.reactor2\_reactor.V}} + \frac{ - \mathtt{reactor2.reactor2\_reactor.outflow.q}\left( t \right) \mathtt{reactor2.reactor2\_reactor.states.X\_BH}\left( t \right)}{\mathtt{reactor2.reactor2\_reactor.V}} + \mathtt{reactor2.reactor2\_reactor.rates.X\_BH}\left( t \right) \\ \frac{\mathrm{d} \mathtt{reactor2.reactor2\_reactor.states.X\_I}\left( t \right)}{\mathrm{d}t} &= \frac{ - \mathtt{reactor2.reactor2\_reactor.outflow.q}\left( t \right) \mathtt{reactor2.reactor2\_reactor.states.X\_I}\left( t \right)}{\mathtt{reactor2.reactor2\_reactor.V}} + \mathtt{reactor2.reactor2\_reactor.rates.X\_I}\left( t \right) + \frac{\mathtt{reactor2.reactor2\_reactor.inflow.q}\left( t \right) \mathtt{reactor2.reactor2\_reactor.inflow.X\_I}\left( t \right)}{\mathtt{reactor2.reactor2\_reactor.V}} \\ \frac{\mathrm{d} \mathtt{reactor2.reactor2\_reactor.states.X\_ND}\left( t \right)}{\mathrm{d}t} &= \frac{\mathtt{reactor2.reactor2\_reactor.inflow.q}\left( t \right) \mathtt{reactor2.reactor2\_reactor.inflow.X\_ND}\left( t \right)}{\mathtt{reactor2.reactor2\_reactor.V}} + \frac{ - \mathtt{reactor2.reactor2\_reactor.outflow.q}\left( t \right) \mathtt{reactor2.reactor2\_reactor.states.X\_ND}\left( t \right)}{\mathtt{reactor2.reactor2\_reactor.V}} + \mathtt{reactor2.reactor2\_reactor.rates.X\_ND}\left( t \right) \\ \frac{\mathrm{d} \mathtt{reactor2.reactor2\_reactor.states.X\_P}\left( t \right)}{\mathrm{d}t} &= \frac{ - \mathtt{reactor2.reactor2\_reactor.states.X\_P}\left( t \right) \mathtt{reactor2.reactor2\_reactor.outflow.q}\left( t \right)}{\mathtt{reactor2.reactor2\_reactor.V}} + \frac{\mathtt{reactor2.reactor2\_reactor.inflow.q}\left( t \right) \mathtt{reactor2.reactor2\_reactor.inflow.X\_P}\left( t \right)}{\mathtt{reactor2.reactor2\_reactor.V}} + \mathtt{reactor2.reactor2\_reactor.rates.X\_P}\left( t \right) \\ \frac{\mathrm{d} \mathtt{reactor2.reactor2\_reactor.states.X\_S}\left( t \right)}{\mathrm{d}t} &= \frac{ - \mathtt{reactor2.reactor2\_reactor.states.X\_S}\left( t \right) \mathtt{reactor2.reactor2\_reactor.outflow.q}\left( t \right)}{\mathtt{reactor2.reactor2\_reactor.V}} + \mathtt{reactor2.reactor2\_reactor.rates.X\_S}\left( t \right) + \frac{\mathtt{reactor2.reactor2\_reactor.inflow.q}\left( t \right) \mathtt{reactor2.reactor2\_reactor.inflow.X\_S}\left( t \right)}{\mathtt{reactor2.reactor2\_reactor.V}} \\ \frac{\mathrm{d} \mathtt{reactor3.reactor3\_reactor.states.S\_ALK}\left( t \right)}{\mathrm{d}t} &= \frac{\mathtt{reactor3.reactor3\_reactor.inflow.q}\left( t \right) \mathtt{reactor3.reactor3\_reactor.inflow.S\_ALK}\left( t \right)}{\mathtt{reactor3.reactor3\_reactor.V}} + \mathtt{reactor3.reactor3\_reactor.rates.S\_ALK}\left( t \right) + \frac{ - \mathtt{reactor3.reactor3\_reactor.states.S\_ALK}\left( t \right) \mathtt{reactor3.reactor3\_reactor.outflow.q}\left( t \right)}{\mathtt{reactor3.reactor3\_reactor.V}} \\ \frac{\mathrm{d} \mathtt{reactor3.reactor3\_reactor.states.S\_I}\left( t \right)}{\mathrm{d}t} &= \frac{\mathtt{reactor3.reactor3\_reactor.inflow.q}\left( t \right) \mathtt{reactor3.reactor3\_reactor.inflow.S\_I}\left( t \right)}{\mathtt{reactor3.reactor3\_reactor.V}} + \mathtt{reactor3.reactor3\_reactor.rates.S\_I}\left( t \right) + \frac{ - \mathtt{reactor3.reactor3\_reactor.outflow.q}\left( t \right) \mathtt{reactor3.reactor3\_reactor.states.S\_I}\left( t \right)}{\mathtt{reactor3.reactor3\_reactor.V}} \\ \frac{\mathrm{d} \mathtt{reactor3.reactor3\_reactor.states.S\_ND}\left( t \right)}{\mathrm{d}t} &= \frac{ - \mathtt{reactor3.reactor3\_reactor.outflow.q}\left( t \right) \mathtt{reactor3.reactor3\_reactor.states.S\_ND}\left( t \right)}{\mathtt{reactor3.reactor3\_reactor.V}} + \frac{\mathtt{reactor3.reactor3\_reactor.inflow.q}\left( t \right) \mathtt{reactor3.reactor3\_reactor.inflow.S\_ND}\left( t \right)}{\mathtt{reactor3.reactor3\_reactor.V}} + \mathtt{reactor3.reactor3\_reactor.rates.S\_ND}\left( t \right) \\ \frac{\mathrm{d} \mathtt{reactor3.reactor3\_reactor.states.S\_NH}\left( t \right)}{\mathrm{d}t} &= \frac{\mathtt{reactor3.reactor3\_reactor.inflow.S\_NH}\left( t \right) \mathtt{reactor3.reactor3\_reactor.inflow.q}\left( t \right)}{\mathtt{reactor3.reactor3\_reactor.V}} + \frac{ - \mathtt{reactor3.reactor3\_reactor.states.S\_NH}\left( t \right) \mathtt{reactor3.reactor3\_reactor.outflow.q}\left( t \right)}{\mathtt{reactor3.reactor3\_reactor.V}} + \mathtt{reactor3.reactor3\_reactor.rates.S\_NH}\left( t \right) \\ \frac{\mathrm{d} \mathtt{reactor3.reactor3\_reactor.states.S\_NO}\left( t \right)}{\mathrm{d}t} &= \frac{\mathtt{reactor3.reactor3\_reactor.inflow.q}\left( t \right) \mathtt{reactor3.reactor3\_reactor.inflow.S\_NO}\left( t \right)}{\mathtt{reactor3.reactor3\_reactor.V}} + \frac{ - \mathtt{reactor3.reactor3\_reactor.states.S\_NO}\left( t \right) \mathtt{reactor3.reactor3\_reactor.outflow.q}\left( t \right)}{\mathtt{reactor3.reactor3\_reactor.V}} + \mathtt{reactor3.reactor3\_reactor.rates.S\_NO}\left( t \right) \\ \frac{\mathrm{d} \mathtt{reactor3.reactor3\_reactor.states.S\_O}\left( t \right)}{\mathrm{d}t} &= \frac{ - \mathtt{reactor3.reactor3\_reactor.states.S\_O}\left( t \right) \mathtt{reactor3.reactor3\_reactor.outflow.q}\left( t \right)}{\mathtt{reactor3.reactor3\_reactor.V}} + \mathtt{reactor3.reactor3\_reactor.rates.S\_O}\left( t \right) + \frac{\mathtt{reactor3.reactor3\_reactor.inflow.q}\left( t \right) \mathtt{reactor3.reactor3\_reactor.inflow.S\_O}\left( t \right)}{\mathtt{reactor3.reactor3\_reactor.V}} \\ \frac{\mathrm{d} \mathtt{reactor3.reactor3\_reactor.states.S\_S}\left( t \right)}{\mathrm{d}t} &= \frac{ - \mathtt{reactor3.reactor3\_reactor.states.S\_S}\left( t \right) \mathtt{reactor3.reactor3\_reactor.outflow.q}\left( t \right)}{\mathtt{reactor3.reactor3\_reactor.V}} + \frac{\mathtt{reactor3.reactor3\_reactor.inflow.q}\left( t \right) \mathtt{reactor3.reactor3\_reactor.inflow.S\_S}\left( t \right)}{\mathtt{reactor3.reactor3\_reactor.V}} + \mathtt{reactor3.reactor3\_reactor.rates.S\_S}\left( t \right) \\ \frac{\mathrm{d} \mathtt{reactor3.reactor3\_reactor.states.X\_BA}\left( t \right)}{\mathrm{d}t} &= \mathtt{reactor3.reactor3\_reactor.rates.X\_BA}\left( t \right) + \frac{ - \mathtt{reactor3.reactor3\_reactor.outflow.q}\left( t \right) \mathtt{reactor3.reactor3\_reactor.states.X\_BA}\left( t \right)}{\mathtt{reactor3.reactor3\_reactor.V}} + \frac{\mathtt{reactor3.reactor3\_reactor.inflow.q}\left( t \right) \mathtt{reactor3.reactor3\_reactor.inflow.X\_BA}\left( t \right)}{\mathtt{reactor3.reactor3\_reactor.V}} \\ \frac{\mathrm{d} \mathtt{reactor3.reactor3\_reactor.states.X\_BH}\left( t \right)}{\mathrm{d}t} &= \frac{ - \mathtt{reactor3.reactor3\_reactor.outflow.q}\left( t \right) \mathtt{reactor3.reactor3\_reactor.states.X\_BH}\left( t \right)}{\mathtt{reactor3.reactor3\_reactor.V}} + \mathtt{reactor3.reactor3\_reactor.rates.X\_BH}\left( t \right) + \frac{\mathtt{reactor3.reactor3\_reactor.inflow.q}\left( t \right) \mathtt{reactor3.reactor3\_reactor.inflow.X\_BH}\left( t \right)}{\mathtt{reactor3.reactor3\_reactor.V}} \\ \frac{\mathrm{d} \mathtt{reactor3.reactor3\_reactor.states.X\_I}\left( t \right)}{\mathrm{d}t} &= \frac{ - \mathtt{reactor3.reactor3\_reactor.outflow.q}\left( t \right) \mathtt{reactor3.reactor3\_reactor.states.X\_I}\left( t \right)}{\mathtt{reactor3.reactor3\_reactor.V}} + \frac{\mathtt{reactor3.reactor3\_reactor.inflow.q}\left( t \right) \mathtt{reactor3.reactor3\_reactor.inflow.X\_I}\left( t \right)}{\mathtt{reactor3.reactor3\_reactor.V}} + \mathtt{reactor3.reactor3\_reactor.rates.X\_I}\left( t \right) \\ \frac{\mathrm{d} \mathtt{reactor3.reactor3\_reactor.states.X\_ND}\left( t \right)}{\mathrm{d}t} &= \frac{\mathtt{reactor3.reactor3\_reactor.inflow.q}\left( t \right) \mathtt{reactor3.reactor3\_reactor.inflow.X\_ND}\left( t \right)}{\mathtt{reactor3.reactor3\_reactor.V}} + \mathtt{reactor3.reactor3\_reactor.rates.X\_ND}\left( t \right) + \frac{ - \mathtt{reactor3.reactor3\_reactor.states.X\_ND}\left( t \right) \mathtt{reactor3.reactor3\_reactor.outflow.q}\left( t \right)}{\mathtt{reactor3.reactor3\_reactor.V}} \\ \frac{\mathrm{d} \mathtt{reactor3.reactor3\_reactor.states.X\_P}\left( t \right)}{\mathrm{d}t} &= \frac{\mathtt{reactor3.reactor3\_reactor.inflow.q}\left( t \right) \mathtt{reactor3.reactor3\_reactor.inflow.X\_P}\left( t \right)}{\mathtt{reactor3.reactor3\_reactor.V}} + \frac{ - \mathtt{reactor3.reactor3\_reactor.outflow.q}\left( t \right) \mathtt{reactor3.reactor3\_reactor.states.X\_P}\left( t \right)}{\mathtt{reactor3.reactor3\_reactor.V}} + \mathtt{reactor3.reactor3\_reactor.rates.X\_P}\left( t \right) \\ \frac{\mathrm{d} \mathtt{reactor3.reactor3\_reactor.states.X\_S}\left( t \right)}{\mathrm{d}t} &= \frac{\mathtt{reactor3.reactor3\_reactor.inflow.q}\left( t \right) \mathtt{reactor3.reactor3\_reactor.inflow.X\_S}\left( t \right)}{\mathtt{reactor3.reactor3\_reactor.V}} + \frac{ - \mathtt{reactor3.reactor3\_reactor.states.X\_S}\left( t \right) \mathtt{reactor3.reactor3\_reactor.outflow.q}\left( t \right)}{\mathtt{reactor3.reactor3\_reactor.V}} + \mathtt{reactor3.reactor3\_reactor.rates.X\_S}\left( t \right) \\ \frac{\mathrm{d} \mathtt{reactor4.reactor4\_reactor.states.S\_ALK}\left( t \right)}{\mathrm{d}t} &= \frac{\mathtt{reactor4.reactor4\_reactor.inflow.S\_ALK}\left( t \right) \mathtt{reactor4.reactor4\_reactor.inflow.q}\left( t \right)}{\mathtt{reactor4.reactor4\_reactor.V}} + \mathtt{reactor4.reactor4\_reactor.rates.S\_ALK}\left( t \right) + \frac{ - \mathtt{reactor4.reactor4\_reactor.outflow.q}\left( t \right) \mathtt{reactor4.reactor4\_reactor.states.S\_ALK}\left( t \right)}{\mathtt{reactor4.reactor4\_reactor.V}} \\ \frac{\mathrm{d} \mathtt{reactor4.reactor4\_reactor.states.S\_I}\left( t \right)}{\mathrm{d}t} &= \mathtt{reactor4.reactor4\_reactor.rates.S\_I}\left( t \right) + \frac{ - \mathtt{reactor4.reactor4\_reactor.outflow.q}\left( t \right) \mathtt{reactor4.reactor4\_reactor.states.S\_I}\left( t \right)}{\mathtt{reactor4.reactor4\_reactor.V}} + \frac{\mathtt{reactor4.reactor4\_reactor.inflow.S\_I}\left( t \right) \mathtt{reactor4.reactor4\_reactor.inflow.q}\left( t \right)}{\mathtt{reactor4.reactor4\_reactor.V}} \\ \frac{\mathrm{d} \mathtt{reactor4.reactor4\_reactor.states.S\_ND}\left( t \right)}{\mathrm{d}t} &= \frac{ - \mathtt{reactor4.reactor4\_reactor.states.S\_ND}\left( t \right) \mathtt{reactor4.reactor4\_reactor.outflow.q}\left( t \right)}{\mathtt{reactor4.reactor4\_reactor.V}} + \frac{\mathtt{reactor4.reactor4\_reactor.inflow.q}\left( t \right) \mathtt{reactor4.reactor4\_reactor.inflow.S\_ND}\left( t \right)}{\mathtt{reactor4.reactor4\_reactor.V}} + \mathtt{reactor4.reactor4\_reactor.rates.S\_ND}\left( t \right) \\ \frac{\mathrm{d} \mathtt{reactor4.reactor4\_reactor.states.S\_NH}\left( t \right)}{\mathrm{d}t} &= \frac{\mathtt{reactor4.reactor4\_reactor.inflow.q}\left( t \right) \mathtt{reactor4.reactor4\_reactor.inflow.S\_NH}\left( t \right)}{\mathtt{reactor4.reactor4\_reactor.V}} + \frac{ - \mathtt{reactor4.reactor4\_reactor.states.S\_NH}\left( t \right) \mathtt{reactor4.reactor4\_reactor.outflow.q}\left( t \right)}{\mathtt{reactor4.reactor4\_reactor.V}} + \mathtt{reactor4.reactor4\_reactor.rates.S\_NH}\left( t \right) \\ \frac{\mathrm{d} \mathtt{reactor4.reactor4\_reactor.states.S\_NO}\left( t \right)}{\mathrm{d}t} &= \mathtt{reactor4.reactor4\_reactor.rates.S\_NO}\left( t \right) + \frac{\mathtt{reactor4.reactor4\_reactor.inflow.S\_NO}\left( t \right) \mathtt{reactor4.reactor4\_reactor.inflow.q}\left( t \right)}{\mathtt{reactor4.reactor4\_reactor.V}} + \frac{ - \mathtt{reactor4.reactor4\_reactor.outflow.q}\left( t \right) \mathtt{reactor4.reactor4\_reactor.states.S\_NO}\left( t \right)}{\mathtt{reactor4.reactor4\_reactor.V}} \\ \frac{\mathrm{d} \mathtt{reactor4.reactor4\_reactor.states.S\_O}\left( t \right)}{\mathrm{d}t} &= \frac{ - \mathtt{reactor4.reactor4\_reactor.states.S\_O}\left( t \right) \mathtt{reactor4.reactor4\_reactor.outflow.q}\left( t \right)}{\mathtt{reactor4.reactor4\_reactor.V}} + \frac{\mathtt{reactor4.reactor4\_reactor.inflow.S\_O}\left( t \right) \mathtt{reactor4.reactor4\_reactor.inflow.q}\left( t \right)}{\mathtt{reactor4.reactor4\_reactor.V}} + \mathtt{reactor4.reactor4\_reactor.rates.S\_O}\left( t \right) \\ \frac{\mathrm{d} \mathtt{reactor4.reactor4\_reactor.states.S\_S}\left( t \right)}{\mathrm{d}t} &= \mathtt{reactor4.reactor4\_reactor.rates.S\_S}\left( t \right) + \frac{ - \mathtt{reactor4.reactor4\_reactor.states.S\_S}\left( t \right) \mathtt{reactor4.reactor4\_reactor.outflow.q}\left( t \right)}{\mathtt{reactor4.reactor4\_reactor.V}} + \frac{\mathtt{reactor4.reactor4\_reactor.inflow.S\_S}\left( t \right) \mathtt{reactor4.reactor4\_reactor.inflow.q}\left( t \right)}{\mathtt{reactor4.reactor4\_reactor.V}} \\ \frac{\mathrm{d} \mathtt{reactor4.reactor4\_reactor.states.X\_BA}\left( t \right)}{\mathrm{d}t} &= \frac{ - \mathtt{reactor4.reactor4\_reactor.states.X\_BA}\left( t \right) \mathtt{reactor4.reactor4\_reactor.outflow.q}\left( t \right)}{\mathtt{reactor4.reactor4\_reactor.V}} + \mathtt{reactor4.reactor4\_reactor.rates.X\_BA}\left( t \right) + \frac{\mathtt{reactor4.reactor4\_reactor.inflow.X\_BA}\left( t \right) \mathtt{reactor4.reactor4\_reactor.inflow.q}\left( t \right)}{\mathtt{reactor4.reactor4\_reactor.V}} \\ \frac{\mathrm{d} \mathtt{reactor4.reactor4\_reactor.states.X\_BH}\left( t \right)}{\mathrm{d}t} &= \mathtt{reactor4.reactor4\_reactor.rates.X\_BH}\left( t \right) + \frac{ - \mathtt{reactor4.reactor4\_reactor.outflow.q}\left( t \right) \mathtt{reactor4.reactor4\_reactor.states.X\_BH}\left( t \right)}{\mathtt{reactor4.reactor4\_reactor.V}} + \frac{\mathtt{reactor4.reactor4\_reactor.inflow.X\_BH}\left( t \right) \mathtt{reactor4.reactor4\_reactor.inflow.q}\left( t \right)}{\mathtt{reactor4.reactor4\_reactor.V}} \\ \frac{\mathrm{d} \mathtt{reactor4.reactor4\_reactor.states.X\_I}\left( t \right)}{\mathrm{d}t} &= \mathtt{reactor4.reactor4\_reactor.rates.X\_I}\left( t \right) + \frac{\mathtt{reactor4.reactor4\_reactor.inflow.X\_I}\left( t \right) \mathtt{reactor4.reactor4\_reactor.inflow.q}\left( t \right)}{\mathtt{reactor4.reactor4\_reactor.V}} + \frac{ - \mathtt{reactor4.reactor4\_reactor.states.X\_I}\left( t \right) \mathtt{reactor4.reactor4\_reactor.outflow.q}\left( t \right)}{\mathtt{reactor4.reactor4\_reactor.V}} \\ \frac{\mathrm{d} \mathtt{reactor4.reactor4\_reactor.states.X\_ND}\left( t \right)}{\mathrm{d}t} &= \frac{ - \mathtt{reactor4.reactor4\_reactor.states.X\_ND}\left( t \right) \mathtt{reactor4.reactor4\_reactor.outflow.q}\left( t \right)}{\mathtt{reactor4.reactor4\_reactor.V}} + \mathtt{reactor4.reactor4\_reactor.rates.X\_ND}\left( t \right) + \frac{\mathtt{reactor4.reactor4\_reactor.inflow.X\_ND}\left( t \right) \mathtt{reactor4.reactor4\_reactor.inflow.q}\left( t \right)}{\mathtt{reactor4.reactor4\_reactor.V}} \\ \frac{\mathrm{d} \mathtt{reactor4.reactor4\_reactor.states.X\_P}\left( t \right)}{\mathrm{d}t} &= \mathtt{reactor4.reactor4\_reactor.rates.X\_P}\left( t \right) + \frac{\mathtt{reactor4.reactor4\_reactor.inflow.X\_P}\left( t \right) \mathtt{reactor4.reactor4\_reactor.inflow.q}\left( t \right)}{\mathtt{reactor4.reactor4\_reactor.V}} + \frac{ - \mathtt{reactor4.reactor4\_reactor.states.X\_P}\left( t \right) \mathtt{reactor4.reactor4\_reactor.outflow.q}\left( t \right)}{\mathtt{reactor4.reactor4\_reactor.V}} \\ \frac{\mathrm{d} \mathtt{reactor4.reactor4\_reactor.states.X\_S}\left( t \right)}{\mathrm{d}t} &= \frac{ - \mathtt{reactor4.reactor4\_reactor.outflow.q}\left( t \right) \mathtt{reactor4.reactor4\_reactor.states.X\_S}\left( t \right)}{\mathtt{reactor4.reactor4\_reactor.V}} + \frac{\mathtt{reactor4.reactor4\_reactor.inflow.X\_S}\left( t \right) \mathtt{reactor4.reactor4\_reactor.inflow.q}\left( t \right)}{\mathtt{reactor4.reactor4\_reactor.V}} + \mathtt{reactor4.reactor4\_reactor.rates.X\_S}\left( t \right) \\ \frac{\mathrm{d} \mathtt{reactor5.reactor5\_reactor.states.S\_ALK}\left( t \right)}{\mathrm{d}t} &= \mathtt{reactor5.reactor5\_reactor.rates.S\_ALK}\left( t \right) + \frac{\mathtt{reactor5.reactor5\_reactor.inflow.S\_ALK}\left( t \right) \mathtt{reactor5.reactor5\_reactor.inflow.q}\left( t \right)}{\mathtt{reactor5.reactor5\_reactor.V}} + \frac{ - \mathtt{reactor5.reactor5\_reactor.outflow.q}\left( t \right) \mathtt{reactor5.reactor5\_reactor.states.S\_ALK}\left( t \right)}{\mathtt{reactor5.reactor5\_reactor.V}} \\ \frac{\mathrm{d} \mathtt{reactor5.reactor5\_reactor.states.S\_I}\left( t \right)}{\mathrm{d}t} &= \frac{ - \mathtt{reactor5.reactor5\_reactor.states.S\_I}\left( t \right) \mathtt{reactor5.reactor5\_reactor.outflow.q}\left( t \right)}{\mathtt{reactor5.reactor5\_reactor.V}} + \frac{\mathtt{reactor5.reactor5\_reactor.inflow.S\_I}\left( t \right) \mathtt{reactor5.reactor5\_reactor.inflow.q}\left( t \right)}{\mathtt{reactor5.reactor5\_reactor.V}} + \mathtt{reactor5.reactor5\_reactor.rates.S\_I}\left( t \right) \\ \frac{\mathrm{d} \mathtt{reactor5.reactor5\_reactor.states.S\_ND}\left( t \right)}{\mathrm{d}t} &= \frac{ - \mathtt{reactor5.reactor5\_reactor.outflow.q}\left( t \right) \mathtt{reactor5.reactor5\_reactor.states.S\_ND}\left( t \right)}{\mathtt{reactor5.reactor5\_reactor.V}} + \frac{\mathtt{reactor5.reactor5\_reactor.inflow.S\_ND}\left( t \right) \mathtt{reactor5.reactor5\_reactor.inflow.q}\left( t \right)}{\mathtt{reactor5.reactor5\_reactor.V}} + \mathtt{reactor5.reactor5\_reactor.rates.S\_ND}\left( t \right) \\ \frac{\mathrm{d} \mathtt{reactor5.reactor5\_reactor.states.S\_NH}\left( t \right)}{\mathrm{d}t} &= \frac{\mathtt{reactor5.reactor5\_reactor.inflow.q}\left( t \right) \mathtt{reactor5.reactor5\_reactor.inflow.S\_NH}\left( t \right)}{\mathtt{reactor5.reactor5\_reactor.V}} + \frac{ - \mathtt{reactor5.reactor5\_reactor.outflow.q}\left( t \right) \mathtt{reactor5.reactor5\_reactor.states.S\_NH}\left( t \right)}{\mathtt{reactor5.reactor5\_reactor.V}} + \mathtt{reactor5.reactor5\_reactor.rates.S\_NH}\left( t \right) \\ \frac{\mathrm{d} \mathtt{reactor5.reactor5\_reactor.states.S\_NO}\left( t \right)}{\mathrm{d}t} &= \mathtt{reactor5.reactor5\_reactor.rates.S\_NO}\left( t \right) + \frac{ - \mathtt{reactor5.reactor5\_reactor.outflow.q}\left( t \right) \mathtt{reactor5.reactor5\_reactor.states.S\_NO}\left( t \right)}{\mathtt{reactor5.reactor5\_reactor.V}} + \frac{\mathtt{reactor5.reactor5\_reactor.inflow.q}\left( t \right) \mathtt{reactor5.reactor5\_reactor.inflow.S\_NO}\left( t \right)}{\mathtt{reactor5.reactor5\_reactor.V}} \\ \frac{\mathrm{d} \mathtt{reactor5.reactor5\_reactor.states.S\_O}\left( t \right)}{\mathrm{d}t} &= \frac{\mathtt{reactor5.reactor5\_reactor.inflow.S\_O}\left( t \right) \mathtt{reactor5.reactor5\_reactor.inflow.q}\left( t \right)}{\mathtt{reactor5.reactor5\_reactor.V}} + \mathtt{reactor5.reactor5\_reactor.rates.S\_O}\left( t \right) + \frac{ - \mathtt{reactor5.reactor5\_reactor.outflow.q}\left( t \right) \mathtt{reactor5.reactor5\_reactor.states.S\_O}\left( t \right)}{\mathtt{reactor5.reactor5\_reactor.V}} \\ \frac{\mathrm{d} \mathtt{reactor5.reactor5\_reactor.states.S\_S}\left( t \right)}{\mathrm{d}t} &= \frac{ - \mathtt{reactor5.reactor5\_reactor.outflow.q}\left( t \right) \mathtt{reactor5.reactor5\_reactor.states.S\_S}\left( t \right)}{\mathtt{reactor5.reactor5\_reactor.V}} + \frac{\mathtt{reactor5.reactor5\_reactor.inflow.S\_S}\left( t \right) \mathtt{reactor5.reactor5\_reactor.inflow.q}\left( t \right)}{\mathtt{reactor5.reactor5\_reactor.V}} + \mathtt{reactor5.reactor5\_reactor.rates.S\_S}\left( t \right) \\ \frac{\mathrm{d} \mathtt{reactor5.reactor5\_reactor.states.X\_BA}\left( t \right)}{\mathrm{d}t} &= \frac{\mathtt{reactor5.reactor5\_reactor.inflow.q}\left( t \right) \mathtt{reactor5.reactor5\_reactor.inflow.X\_BA}\left( t \right)}{\mathtt{reactor5.reactor5\_reactor.V}} + \mathtt{reactor5.reactor5\_reactor.rates.X\_BA}\left( t \right) + \frac{ - \mathtt{reactor5.reactor5\_reactor.outflow.q}\left( t \right) \mathtt{reactor5.reactor5\_reactor.states.X\_BA}\left( t \right)}{\mathtt{reactor5.reactor5\_reactor.V}} \\ \frac{\mathrm{d} \mathtt{reactor5.reactor5\_reactor.states.X\_BH}\left( t \right)}{\mathrm{d}t} &= \mathtt{reactor5.reactor5\_reactor.rates.X\_BH}\left( t \right) + \frac{ - \mathtt{reactor5.reactor5\_reactor.outflow.q}\left( t \right) \mathtt{reactor5.reactor5\_reactor.states.X\_BH}\left( t \right)}{\mathtt{reactor5.reactor5\_reactor.V}} + \frac{\mathtt{reactor5.reactor5\_reactor.inflow.X\_BH}\left( t \right) \mathtt{reactor5.reactor5\_reactor.inflow.q}\left( t \right)}{\mathtt{reactor5.reactor5\_reactor.V}} \\ \frac{\mathrm{d} \mathtt{reactor5.reactor5\_reactor.states.X\_I}\left( t \right)}{\mathrm{d}t} &= \mathtt{reactor5.reactor5\_reactor.rates.X\_I}\left( t \right) + \frac{ - \mathtt{reactor5.reactor5\_reactor.outflow.q}\left( t \right) \mathtt{reactor5.reactor5\_reactor.states.X\_I}\left( t \right)}{\mathtt{reactor5.reactor5\_reactor.V}} + \frac{\mathtt{reactor5.reactor5\_reactor.inflow.X\_I}\left( t \right) \mathtt{reactor5.reactor5\_reactor.inflow.q}\left( t \right)}{\mathtt{reactor5.reactor5\_reactor.V}} \\ \frac{\mathrm{d} \mathtt{reactor5.reactor5\_reactor.states.X\_ND}\left( t \right)}{\mathrm{d}t} &= \frac{\mathtt{reactor5.reactor5\_reactor.inflow.X\_ND}\left( t \right) \mathtt{reactor5.reactor5\_reactor.inflow.q}\left( t \right)}{\mathtt{reactor5.reactor5\_reactor.V}} + \frac{ - \mathtt{reactor5.reactor5\_reactor.outflow.q}\left( t \right) \mathtt{reactor5.reactor5\_reactor.states.X\_ND}\left( t \right)}{\mathtt{reactor5.reactor5\_reactor.V}} + \mathtt{reactor5.reactor5\_reactor.rates.X\_ND}\left( t \right) \\ \frac{\mathrm{d} \mathtt{reactor5.reactor5\_reactor.states.X\_P}\left( t \right)}{\mathrm{d}t} &= \frac{\mathtt{reactor5.reactor5\_reactor.inflow.X\_P}\left( t \right) \mathtt{reactor5.reactor5\_reactor.inflow.q}\left( t \right)}{\mathtt{reactor5.reactor5\_reactor.V}} + \frac{ - \mathtt{reactor5.reactor5\_reactor.states.X\_P}\left( t \right) \mathtt{reactor5.reactor5\_reactor.outflow.q}\left( t \right)}{\mathtt{reactor5.reactor5\_reactor.V}} + \mathtt{reactor5.reactor5\_reactor.rates.X\_P}\left( t \right) \\ \frac{\mathrm{d} \mathtt{reactor5.reactor5\_reactor.states.X\_S}\left( t \right)}{\mathrm{d}t} &= \mathtt{reactor5.reactor5\_reactor.rates.X\_S}\left( t \right) + \frac{ - \mathtt{reactor5.reactor5\_reactor.outflow.q}\left( t \right) \mathtt{reactor5.reactor5\_reactor.states.X\_S}\left( t \right)}{\mathtt{reactor5.reactor5\_reactor.V}} + \frac{\mathtt{reactor5.reactor5\_reactor.inflow.X\_S}\left( t \right) \mathtt{reactor5.reactor5\_reactor.inflow.q}\left( t \right)}{\mathtt{reactor5.reactor5\_reactor.V}} \\ 0 &= - \mathtt{unifier.inflow\_3.q}\left( t \right) \mathtt{unifier.inflow\_3.S\_ALK}\left( t \right) - \mathtt{unifier.inflow\_2.S\_ALK}\left( t \right) \mathtt{unifier.inflow\_2.q}\left( t \right) + \mathtt{unifier.outflow.q}\left( t \right) \mathtt{unifier.outflow.S\_ALK}\left( t \right) - \mathtt{unifier.inflow\_1.S\_ALK}\left( t \right) \mathtt{unifier.inflow\_1.q}\left( t \right) \\ 0 &= - \mathtt{unifier.inflow\_3.q}\left( t \right) \mathtt{unifier.inflow\_3.S\_I}\left( t \right) + \mathtt{unifier.outflow.S\_I}\left( t \right) \mathtt{unifier.outflow.q}\left( t \right) - \mathtt{unifier.inflow\_2.q}\left( t \right) \mathtt{unifier.inflow\_2.S\_I}\left( t \right) - \mathtt{unifier.inflow\_1.S\_I}\left( t \right) \mathtt{unifier.inflow\_1.q}\left( t \right) \\ 0 &= - \mathtt{unifier.inflow\_3.q}\left( t \right) \mathtt{unifier.inflow\_3.S\_ND}\left( t \right) - \mathtt{unifier.inflow\_2.q}\left( t \right) \mathtt{unifier.inflow\_2.S\_ND}\left( t \right) + \mathtt{unifier.outflow.S\_ND}\left( t \right) \mathtt{unifier.outflow.q}\left( t \right) - \mathtt{unifier.inflow\_1.S\_ND}\left( t \right) \mathtt{unifier.inflow\_1.q}\left( t \right) \\ 0 &= - \mathtt{unifier.inflow\_3.q}\left( t \right) \mathtt{unifier.inflow\_3.S\_NH}\left( t \right) - \mathtt{unifier.inflow\_2.S\_NH}\left( t \right) \mathtt{unifier.inflow\_2.q}\left( t \right) - \mathtt{unifier.inflow\_1.S\_NH}\left( t \right) \mathtt{unifier.inflow\_1.q}\left( t \right) + \mathtt{unifier.outflow.q}\left( t \right) \mathtt{unifier.outflow.S\_NH}\left( t \right) \\ 0 &= - \mathtt{unifier.inflow\_3.q}\left( t \right) \mathtt{unifier.inflow\_3.S\_NO}\left( t \right) - \mathtt{unifier.inflow\_2.S\_NO}\left( t \right) \mathtt{unifier.inflow\_2.q}\left( t \right) + \mathtt{unifier.outflow.S\_NO}\left( t \right) \mathtt{unifier.outflow.q}\left( t \right) - \mathtt{unifier.inflow\_1.S\_NO}\left( t \right) \mathtt{unifier.inflow\_1.q}\left( t \right) \\ 0 &= - \mathtt{unifier.inflow\_3.q}\left( t \right) \mathtt{unifier.inflow\_3.S\_O}\left( t \right) - \mathtt{unifier.inflow\_2.S\_O}\left( t \right) \mathtt{unifier.inflow\_2.q}\left( t \right) + \mathtt{unifier.outflow.S\_O}\left( t \right) \mathtt{unifier.outflow.q}\left( t \right) - \mathtt{unifier.inflow\_1.S\_O}\left( t \right) \mathtt{unifier.inflow\_1.q}\left( t \right) \\ 0 &= - \mathtt{unifier.inflow\_3.q}\left( t \right) \mathtt{unifier.inflow\_3.S\_S}\left( t \right) + \mathtt{unifier.outflow.S\_S}\left( t \right) \mathtt{unifier.outflow.q}\left( t \right) - \mathtt{unifier.inflow\_2.q}\left( t \right) \mathtt{unifier.inflow\_2.S\_S}\left( t \right) - \mathtt{unifier.inflow\_1.S\_S}\left( t \right) \mathtt{unifier.inflow\_1.q}\left( t \right) \\ 0 &= - \mathtt{unifier.inflow\_3.q}\left( t \right) \mathtt{unifier.inflow\_3.X\_BA}\left( t \right) - \mathtt{unifier.inflow\_1.X\_BA}\left( t \right) \mathtt{unifier.inflow\_1.q}\left( t \right) - \mathtt{unifier.inflow\_2.q}\left( t \right) \mathtt{unifier.inflow\_2.X\_BA}\left( t \right) + \mathtt{unifier.outflow.X\_BA}\left( t \right) \mathtt{unifier.outflow.q}\left( t \right) \\ 0 &= - \mathtt{unifier.inflow\_3.q}\left( t \right) \mathtt{unifier.inflow\_3.X\_BH}\left( t \right) - \mathtt{unifier.inflow\_2.X\_BH}\left( t \right) \mathtt{unifier.inflow\_2.q}\left( t \right) - \mathtt{unifier.inflow\_1.X\_BH}\left( t \right) \mathtt{unifier.inflow\_1.q}\left( t \right) + \mathtt{unifier.outflow.q}\left( t \right) \mathtt{unifier.outflow.X\_BH}\left( t \right) \\ 0 &= - \mathtt{unifier.inflow\_3.q}\left( t \right) \mathtt{unifier.inflow\_3.X\_I}\left( t \right) - \mathtt{unifier.inflow\_2.X\_I}\left( t \right) \mathtt{unifier.inflow\_2.q}\left( t \right) + \mathtt{unifier.outflow.q}\left( t \right) \mathtt{unifier.outflow.X\_I}\left( t \right) - \mathtt{unifier.inflow\_1.X\_I}\left( t \right) \mathtt{unifier.inflow\_1.q}\left( t \right) \\ 0 &= - \mathtt{unifier.inflow\_3.q}\left( t \right) \mathtt{unifier.inflow\_3.X\_ND}\left( t \right) - \mathtt{unifier.inflow\_1.X\_ND}\left( t \right) \mathtt{unifier.inflow\_1.q}\left( t \right) - \mathtt{unifier.inflow\_2.X\_ND}\left( t \right) \mathtt{unifier.inflow\_2.q}\left( t \right) + \mathtt{unifier.outflow.q}\left( t \right) \mathtt{unifier.outflow.X\_ND}\left( t \right) \\ 0 &= - \mathtt{unifier.inflow\_3.q}\left( t \right) \mathtt{unifier.inflow\_3.X\_P}\left( t \right) - \mathtt{unifier.inflow\_2.X\_P}\left( t \right) \mathtt{unifier.inflow\_2.q}\left( t \right) + \mathtt{unifier.outflow.q}\left( t \right) \mathtt{unifier.outflow.X\_P}\left( t \right) - \mathtt{unifier.inflow\_1.X\_P}\left( t \right) \mathtt{unifier.inflow\_1.q}\left( t \right) \\ 0 &= - \mathtt{unifier.inflow\_3.q}\left( t \right) \mathtt{unifier.inflow\_3.X\_S}\left( t \right) - \mathtt{unifier.inflow\_2.X\_S}\left( t \right) \mathtt{unifier.inflow\_2.q}\left( t \right) + \mathtt{unifier.outflow.X\_S}\left( t \right) \mathtt{unifier.outflow.q}\left( t \right) - \mathtt{unifier.inflow\_1.X\_S}\left( t \right) \mathtt{unifier.inflow\_1.q}\left( t \right) \\ 0 &= \mathtt{internal\_rec\_sep.outflow\_1.S\_ALK}\left( t \right) \mathtt{internal\_rec\_sep.outflow\_1.q}\left( t \right) - \mathtt{internal\_rec\_sep.inflow.q}\left( t \right) \mathtt{internal\_rec\_sep.inflow.S\_ALK}\left( t \right) + \mathtt{internal\_rec\_sep.outflow\_2.S\_ALK}\left( t \right) \mathtt{internal\_rec\_sep.outflow\_2.q}\left( t \right) \\ 0 &= \mathtt{internal\_rec\_sep.outflow\_2.S\_I}\left( t \right) \mathtt{internal\_rec\_sep.outflow\_2.q}\left( t \right) - \mathtt{internal\_rec\_sep.inflow.q}\left( t \right) \mathtt{internal\_rec\_sep.inflow.S\_I}\left( t \right) + \mathtt{internal\_rec\_sep.outflow\_1.q}\left( t \right) \mathtt{internal\_rec\_sep.outflow\_1.S\_I}\left( t \right) \\ 0 &= - \mathtt{internal\_rec\_sep.inflow.S\_ND}\left( t \right) \mathtt{internal\_rec\_sep.inflow.q}\left( t \right) + \mathtt{internal\_rec\_sep.outflow\_2.q}\left( t \right) \mathtt{internal\_rec\_sep.outflow\_2.S\_ND}\left( t \right) + \mathtt{internal\_rec\_sep.outflow\_1.S\_ND}\left( t \right) \mathtt{internal\_rec\_sep.outflow\_1.q}\left( t \right) \\ 0 &= - \mathtt{internal\_rec\_sep.inflow.q}\left( t \right) \mathtt{internal\_rec\_sep.inflow.S\_NH}\left( t \right) + \mathtt{internal\_rec\_sep.outflow\_2.S\_NH}\left( t \right) \mathtt{internal\_rec\_sep.outflow\_2.q}\left( t \right) + \mathtt{internal\_rec\_sep.outflow\_1.S\_NH}\left( t \right) \mathtt{internal\_rec\_sep.outflow\_1.q}\left( t \right) \\ 0 &= \mathtt{internal\_rec\_sep.outflow\_2.S\_NO}\left( t \right) \mathtt{internal\_rec\_sep.outflow\_2.q}\left( t \right) - \mathtt{internal\_rec\_sep.inflow.q}\left( t \right) \mathtt{internal\_rec\_sep.inflow.S\_NO}\left( t \right) + \mathtt{internal\_rec\_sep.outflow\_1.q}\left( t \right) \mathtt{internal\_rec\_sep.outflow\_1.S\_NO}\left( t \right) \\ 0 &= - \mathtt{internal\_rec\_sep.inflow.S\_O}\left( t \right) \mathtt{internal\_rec\_sep.inflow.q}\left( t \right) + \mathtt{internal\_rec\_sep.outflow\_2.S\_O}\left( t \right) \mathtt{internal\_rec\_sep.outflow\_2.q}\left( t \right) + \mathtt{internal\_rec\_sep.outflow\_1.S\_O}\left( t \right) \mathtt{internal\_rec\_sep.outflow\_1.q}\left( t \right) \\ 0 &= - \mathtt{internal\_rec\_sep.inflow.q}\left( t \right) \mathtt{internal\_rec\_sep.inflow.S\_S}\left( t \right) + \mathtt{internal\_rec\_sep.outflow\_2.q}\left( t \right) \mathtt{internal\_rec\_sep.outflow\_2.S\_S}\left( t \right) + \mathtt{internal\_rec\_sep.outflow\_1.q}\left( t \right) \mathtt{internal\_rec\_sep.outflow\_1.S\_S}\left( t \right) \\ 0 &= \mathtt{internal\_rec\_sep.outflow\_2.X\_BA}\left( t \right) \mathtt{internal\_rec\_sep.outflow\_2.q}\left( t \right) - \mathtt{internal\_rec\_sep.inflow.q}\left( t \right) \mathtt{internal\_rec\_sep.inflow.X\_BA}\left( t \right) + \mathtt{internal\_rec\_sep.outflow\_1.X\_BA}\left( t \right) \mathtt{internal\_rec\_sep.outflow\_1.q}\left( t \right) \\ 0 &= - \mathtt{internal\_rec\_sep.inflow.q}\left( t \right) \mathtt{internal\_rec\_sep.inflow.X\_BH}\left( t \right) + \mathtt{internal\_rec\_sep.outflow\_1.X\_BH}\left( t \right) \mathtt{internal\_rec\_sep.outflow\_1.q}\left( t \right) + \mathtt{internal\_rec\_sep.outflow\_2.q}\left( t \right) \mathtt{internal\_rec\_sep.outflow\_2.X\_BH}\left( t \right) \\ 0 &= - \mathtt{internal\_rec\_sep.inflow.q}\left( t \right) \mathtt{internal\_rec\_sep.inflow.X\_I}\left( t \right) + \mathtt{internal\_rec\_sep.outflow\_1.X\_I}\left( t \right) \mathtt{internal\_rec\_sep.outflow\_1.q}\left( t \right) + \mathtt{internal\_rec\_sep.outflow\_2.q}\left( t \right) \mathtt{internal\_rec\_sep.outflow\_2.X\_I}\left( t \right) \\ 0 &= - \mathtt{internal\_rec\_sep.inflow.q}\left( t \right) \mathtt{internal\_rec\_sep.inflow.X\_ND}\left( t \right) + \mathtt{internal\_rec\_sep.outflow\_2.X\_ND}\left( t \right) \mathtt{internal\_rec\_sep.outflow\_2.q}\left( t \right) + \mathtt{internal\_rec\_sep.outflow\_1.q}\left( t \right) \mathtt{internal\_rec\_sep.outflow\_1.X\_ND}\left( t \right) \\ 0 &= - \mathtt{internal\_rec\_sep.inflow.q}\left( t \right) \mathtt{internal\_rec\_sep.inflow.X\_P}\left( t \right) + \mathtt{internal\_rec\_sep.outflow\_2.q}\left( t \right) \mathtt{internal\_rec\_sep.outflow\_2.X\_P}\left( t \right) + \mathtt{internal\_rec\_sep.outflow\_1.X\_P}\left( t \right) \mathtt{internal\_rec\_sep.outflow\_1.q}\left( t \right) \\ 0 &= - \mathtt{internal\_rec\_sep.inflow.q}\left( t \right) \mathtt{internal\_rec\_sep.inflow.X\_S}\left( t \right) + \mathtt{internal\_rec\_sep.outflow\_2.q}\left( t \right) \mathtt{internal\_rec\_sep.outflow\_2.X\_S}\left( t \right) + \mathtt{internal\_rec\_sep.outflow\_1.X\_S}\left( t \right) \mathtt{internal\_rec\_sep.outflow\_1.q}\left( t \right) \\ 0 &= \mathtt{clarifier.sludge.S\_ALK}\left( t \right) \mathtt{clarifier.sludge.q}\left( t \right) - \mathtt{clarifier.inflow.S\_ALK}\left( t \right) \mathtt{clarifier.inflow.q}\left( t \right) + \mathtt{clarifier.effluent.q}\left( t \right) \mathtt{clarifier.effluent.S\_ALK}\left( t \right) \\ 0 &= \mathtt{clarifier.sludge.q}\left( t \right) \mathtt{clarifier.sludge.S\_I}\left( t \right) + \mathtt{clarifier.effluent.q}\left( t \right) \mathtt{clarifier.effluent.S\_I}\left( t \right) - \mathtt{clarifier.inflow.q}\left( t \right) \mathtt{clarifier.inflow.S\_I}\left( t \right) \\ 0 &= - \mathtt{clarifier.inflow.S\_ND}\left( t \right) \mathtt{clarifier.inflow.q}\left( t \right) + \mathtt{clarifier.sludge.q}\left( t \right) \mathtt{clarifier.sludge.S\_ND}\left( t \right) + \mathtt{clarifier.effluent.q}\left( t \right) \mathtt{clarifier.effluent.S\_ND}\left( t \right) \\ 0 &= \mathtt{clarifier.sludge.q}\left( t \right) \mathtt{clarifier.sludge.S\_NH}\left( t \right) + \mathtt{clarifier.effluent.q}\left( t \right) \mathtt{clarifier.effluent.S\_NH}\left( t \right) - \mathtt{clarifier.inflow.q}\left( t \right) \mathtt{clarifier.inflow.S\_NH}\left( t \right) \\ 0 &= \mathtt{clarifier.sludge.q}\left( t \right) \mathtt{clarifier.sludge.S\_NO}\left( t \right) + \mathtt{clarifier.effluent.q}\left( t \right) \mathtt{clarifier.effluent.S\_NO}\left( t \right) - \mathtt{clarifier.inflow.q}\left( t \right) \mathtt{clarifier.inflow.S\_NO}\left( t \right) \\ 0 &= \mathtt{clarifier.sludge.q}\left( t \right) \mathtt{clarifier.sludge.S\_O}\left( t \right) - \mathtt{clarifier.inflow.S\_O}\left( t \right) \mathtt{clarifier.inflow.q}\left( t \right) + \mathtt{clarifier.effluent.q}\left( t \right) \mathtt{clarifier.effluent.S\_O}\left( t \right) \\ 0 &= - \mathtt{clarifier.inflow.S\_S}\left( t \right) \mathtt{clarifier.inflow.q}\left( t \right) + \mathtt{clarifier.sludge.q}\left( t \right) \mathtt{clarifier.sludge.S\_S}\left( t \right) + \mathtt{clarifier.effluent.S\_S}\left( t \right) \mathtt{clarifier.effluent.q}\left( t \right) \\ 0 &= - \mathtt{clarifier.inflow.X\_BA}\left( t \right) \mathtt{clarifier.inflow.q}\left( t \right) + \mathtt{clarifier.effluent.X\_BA}\left( t \right) \mathtt{clarifier.effluent.q}\left( t \right) + \mathtt{clarifier.sludge.q}\left( t \right) \mathtt{clarifier.sludge.X\_BA}\left( t \right) \\ 0 &= \mathtt{clarifier.sludge.X\_BH}\left( t \right) \mathtt{clarifier.sludge.q}\left( t \right) - \mathtt{clarifier.inflow.X\_BH}\left( t \right) \mathtt{clarifier.inflow.q}\left( t \right) + \mathtt{clarifier.effluent.q}\left( t \right) \mathtt{clarifier.effluent.X\_BH}\left( t \right) \\ 0 &= \mathtt{clarifier.sludge.q}\left( t \right) \mathtt{clarifier.sludge.X\_I}\left( t \right) + \mathtt{clarifier.effluent.q}\left( t \right) \mathtt{clarifier.effluent.X\_I}\left( t \right) - \mathtt{clarifier.inflow.q}\left( t \right) \mathtt{clarifier.inflow.X\_I}\left( t \right) \\ 0 &= - \mathtt{clarifier.inflow.X\_ND}\left( t \right) \mathtt{clarifier.inflow.q}\left( t \right) + \mathtt{clarifier.sludge.X\_ND}\left( t \right) \mathtt{clarifier.sludge.q}\left( t \right) + \mathtt{clarifier.effluent.q}\left( t \right) \mathtt{clarifier.effluent.X\_ND}\left( t \right) \\ 0 &= \mathtt{clarifier.effluent.X\_P}\left( t \right) \mathtt{clarifier.effluent.q}\left( t \right) + \mathtt{clarifier.sludge.q}\left( t \right) \mathtt{clarifier.sludge.X\_P}\left( t \right) - \mathtt{clarifier.inflow.q}\left( t \right) \mathtt{clarifier.inflow.X\_P}\left( t \right) \\ 0 &= - \mathtt{clarifier.inflow.X\_S}\left( t \right) \mathtt{clarifier.inflow.q}\left( t \right) + \mathtt{clarifier.effluent.X\_S}\left( t \right) \mathtt{clarifier.effluent.q}\left( t \right) + \mathtt{clarifier.sludge.q}\left( t \right) \mathtt{clarifier.sludge.X\_S}\left( t \right) \\ 0 &= \mathtt{wastage\_separator.outflow\_1.q}\left( t \right) \mathtt{wastage\_separator.outflow\_1.S\_ALK}\left( t \right) - \mathtt{wastage\_separator.inflow.q}\left( t \right) \mathtt{wastage\_separator.inflow.S\_ALK}\left( t \right) + \mathtt{wastage\_separator.outflow\_2.q}\left( t \right) \mathtt{wastage\_separator.outflow\_2.S\_ALK}\left( t \right) \\ 0 &= \mathtt{wastage\_separator.outflow\_1.q}\left( t \right) \mathtt{wastage\_separator.outflow\_1.S\_I}\left( t \right) - \mathtt{wastage\_separator.inflow.q}\left( t \right) \mathtt{wastage\_separator.inflow.S\_I}\left( t \right) + \mathtt{wastage\_separator.outflow\_2.q}\left( t \right) \mathtt{wastage\_separator.outflow\_2.S\_I}\left( t \right) \\ 0 &= \mathtt{wastage\_separator.outflow\_1.q}\left( t \right) \mathtt{wastage\_separator.outflow\_1.S\_ND}\left( t \right) - \mathtt{wastage\_separator.inflow.q}\left( t \right) \mathtt{wastage\_separator.inflow.S\_ND}\left( t \right) + \mathtt{wastage\_separator.outflow\_2.q}\left( t \right) \mathtt{wastage\_separator.outflow\_2.S\_ND}\left( t \right) \\ 0 &= \mathtt{wastage\_separator.outflow\_1.q}\left( t \right) \mathtt{wastage\_separator.outflow\_1.S\_NH}\left( t \right) - \mathtt{wastage\_separator.inflow.q}\left( t \right) \mathtt{wastage\_separator.inflow.S\_NH}\left( t \right) + \mathtt{wastage\_separator.outflow\_2.q}\left( t \right) \mathtt{wastage\_separator.outflow\_2.S\_NH}\left( t \right) \\ 0 &= \mathtt{wastage\_separator.outflow\_1.q}\left( t \right) \mathtt{wastage\_separator.outflow\_1.S\_NO}\left( t \right) - \mathtt{wastage\_separator.inflow.q}\left( t \right) \mathtt{wastage\_separator.inflow.S\_NO}\left( t \right) + \mathtt{wastage\_separator.outflow\_2.q}\left( t \right) \mathtt{wastage\_separator.outflow\_2.S\_NO}\left( t \right) \\ 0 &= \mathtt{wastage\_separator.outflow\_1.S\_O}\left( t \right) \mathtt{wastage\_separator.outflow\_1.q}\left( t \right) - \mathtt{wastage\_separator.inflow.q}\left( t \right) \mathtt{wastage\_separator.inflow.S\_O}\left( t \right) + \mathtt{wastage\_separator.outflow\_2.q}\left( t \right) \mathtt{wastage\_separator.outflow\_2.S\_O}\left( t \right) \\ 0 &= \mathtt{wastage\_separator.outflow\_1.q}\left( t \right) \mathtt{wastage\_separator.outflow\_1.S\_S}\left( t \right) - \mathtt{wastage\_separator.inflow.q}\left( t \right) \mathtt{wastage\_separator.inflow.S\_S}\left( t \right) + \mathtt{wastage\_separator.outflow\_2.q}\left( t \right) \mathtt{wastage\_separator.outflow\_2.S\_S}\left( t \right) \\ 0 &= \mathtt{wastage\_separator.outflow\_1.q}\left( t \right) \mathtt{wastage\_separator.outflow\_1.X\_BA}\left( t \right) - \mathtt{wastage\_separator.inflow.q}\left( t \right) \mathtt{wastage\_separator.inflow.X\_BA}\left( t \right) + \mathtt{wastage\_separator.outflow\_2.q}\left( t \right) \mathtt{wastage\_separator.outflow\_2.X\_BA}\left( t \right) \\ 0 &= \mathtt{wastage\_separator.outflow\_1.q}\left( t \right) \mathtt{wastage\_separator.outflow\_1.X\_BH}\left( t \right) - \mathtt{wastage\_separator.inflow.q}\left( t \right) \mathtt{wastage\_separator.inflow.X\_BH}\left( t \right) + \mathtt{wastage\_separator.outflow\_2.q}\left( t \right) \mathtt{wastage\_separator.outflow\_2.X\_BH}\left( t \right) \\ 0 &= \mathtt{wastage\_separator.outflow\_1.q}\left( t \right) \mathtt{wastage\_separator.outflow\_1.X\_I}\left( t \right) - \mathtt{wastage\_separator.inflow.q}\left( t \right) \mathtt{wastage\_separator.inflow.X\_I}\left( t \right) + \mathtt{wastage\_separator.outflow\_2.q}\left( t \right) \mathtt{wastage\_separator.outflow\_2.X\_I}\left( t \right) \\ 0 &= \mathtt{wastage\_separator.outflow\_1.q}\left( t \right) \mathtt{wastage\_separator.outflow\_1.X\_ND}\left( t \right) - \mathtt{wastage\_separator.inflow.q}\left( t \right) \mathtt{wastage\_separator.inflow.X\_ND}\left( t \right) + \mathtt{wastage\_separator.outflow\_2.q}\left( t \right) \mathtt{wastage\_separator.outflow\_2.X\_ND}\left( t \right) \\ 0 &= \mathtt{wastage\_separator.outflow\_1.X\_P}\left( t \right) \mathtt{wastage\_separator.outflow\_1.q}\left( t \right) - \mathtt{wastage\_separator.inflow.q}\left( t \right) \mathtt{wastage\_separator.inflow.X\_P}\left( t \right) + \mathtt{wastage\_separator.outflow\_2.q}\left( t \right) \mathtt{wastage\_separator.outflow\_2.X\_P}\left( t \right) \\ 0 &= - \mathtt{wastage\_separator.inflow.X\_S}\left( t \right) \mathtt{wastage\_separator.inflow.q}\left( t \right) + \mathtt{wastage\_separator.outflow\_1.q}\left( t \right) \mathtt{wastage\_separator.outflow\_1.X\_S}\left( t \right) + \mathtt{wastage\_separator.outflow\_2.q}\left( t \right) \mathtt{wastage\_separator.outflow\_2.X\_S}\left( t \right) \end{align} \]
Checking the connections
Now, the overall system has been build and can be simulated. However, are the individual reactors (and other elements) correctly connected? Obviously, for this example this is the case (since otherwise it would not be a good example), but here how to check it graphically (check the corresponding part in the usage first, for the requirements for printing).
A graph can be generated using the ProcessDiagram as follows:
using TikzPictures # Needed for displaying here
# Get the process diagram for displaying
pd = ProcessDiagram(sys, eqns)
And that easily you can get a nice diagram of the flows in the system. (Yes, everything is very rectangular right now and systems only identified by names at the current stage. Might be changed in the future).
But this diagram shows only the flows and not other internal connections, e.g. between the reactors and "controllers". What about errors in the connections there? Luckily, these connections can be shown similarily, it is just disabled by default, to have a cleaner graph. They can be included if the diagram is generated using:
# Get the process diagram for displaying
pd_full = ProcessDiagram(sys, eqns; only_flow=false)
Further, the names of the connectors can be displayed as well, for better debugging (not nicely displayed), using:
print_connector_names!(pd) # or directly pass print_connector_names=true to the ProcessDiagram
Simulation and Plotting
Now, with the overall system at hand and the connections checked, it can be simulated and the output plotted
## Generate the problem to solve (114 days, stabilization + 2 weeks of influent)
t_sim = (0, 114)
prob = ODEProblem(bsm1_simplified, [], t_sim, [], initializealg=ShampineCollocationInit()) # We need to explicitly select the initialization algorithm here. The default doesn't work
# Solve the problem. Use a DAE solver (Rosenbrock methods, see solver docs for more infos on the method)
sol = solve(prob, Rodas5P())
## Plot the output
plot(sol,
xlim = (95, 114), #5 last days of stabilisation and then the influent file
idxs = unknowns(bsm1_model)[occursin.("clarifier₊effluent", string.(unknowns(bsm1_model)))][2:end], # States which are associated with effluent but not flow
label = permutedims(string.(unknowns(outflows(clarifier, :effluent))[2:end])),
legend = :topleft,
title = "BSM1 Effluent",
)