Reactor with multiple Processes: ASM1 with Aeration

This example shows the current procedure of how multiple processes are combined in a single reactor. For this, we use an examplary system of a single aerated reactor with ASM1 and constant input:

Now, let's see how this can be accomplished in three different cases using the high level interface. If you are interested in the inner workings, refer to the same example using the low-level interface.

The three presented cases are:

Combining Variables of the Same Name
Adding a Process to the Default One
Combining Variables with Different Names

Combining Variables of the Same Name

The first example of multiple processes, is concerning the case where two variables have the same name. To do so, two processes are given to the reactor and below it is shown how this is done.

Plug and Play Example

using BioChemicalTreatment # This Package :)
using ModelingToolkit # Basic Modeling Framework. Takes e.g. care of simplification of the equations
using ModelingToolkitStandardLibrary # Standard library for MTK. Provides e.g. simple math operations (adder, product...), controllers (PID) and more
using DifferentialEquations # Solve the equations
using Plots # Plotting

# Processes
@named asm = ASM1() # Use the default parameters
@named aer = Aeration() # The aeration process

# Reactor
@named react = CSTR(1000, [asm, aer]; initial_states=[1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1]) # Volume of 1000, initial_states of all 1 (the default of all 0 does yield infinity, check the model for yourself)

# Influent and aeration input
@named influent = Influent([7.0, 30.0, 6.95, 31.56, 0.0, 0.0, 69.5, 0.0, 28.17, 51.2, 10.59, 0.0, 202.32], unknowns(states(asm)), flowrate = 18446) # BSM1 Stabilization
@named aeration_ctrl = ModelingToolkitStandardLibrary.Blocks.Constant(;k = 240)

# Connection
eqs = [
  connect(outflows(influent, 1), inflows(react, 1))
  connect(exogenous_inputs(react, :aer₊k_La), aeration_ctrl.output)
]

# Model generation
@named model = ODESystem(eqs, t, systems=[influent, react, aeration_ctrl])
model_simplified = structural_simplify(model)

# Simulation
ode_prob = ODEProblem(model_simplified, [], (0, 1))
sol = solve(ode_prob)

# Plotting, only oxygen
plot(sol, idxs=[inflows(react, 1).S_O, states(react).S_O], label=["Inflow" "Reactor"], title="Oxygen concentration in CSTR with ASM1 and aeration", legend=:topright)

# Plotting the rates
plot(sol, idxs=[rates(react).S_O, subsystems(react, :asm).rates.S_O, subsystems(react, :aer).rates.S_O], label=["Combined" "ASM" "Aeration"], title="Oxygen rate in CSTR with ASM1 and aeration", legend=:topright)

Step-by-Step Explanations

As always, the first step is the setup in loading the packages:

using BioChemicalTreatment # This Package :)
using ModelingToolkit # Basic Modeling Framework. Takes e.g. care of simplification of the equations
using ModelingToolkitStandardLibrary # Standard library for MTK. Provides e.g. simple math operations (adder, product...), controllers (PID) and more
using DifferentialEquations # Solve the equations
using Plots # Plotting

And then, the processes (in this case two of them) are created:

@named asm = ASM1() # Use the default parameters
@named aer = Aeration() # The aeration process

Looking at the states of the two processes, it is visible that both have a state S_O, which is the one to be combined:

unknowns(states(asm))' # Transpose for better visibility

1×13 adjoint(::Vector{SymbolicUtils.BasicSymbolic{Real}}) with eltype SymbolicUtils.BasicSymbolic{Real}:
 S_ALK(t)  S_I(t)  S_ND(t)  S_NH(t)  …  X_I(t)  X_ND(t)  X_P(t)  X_S(t)

unknowns(states(aer))' # Transpose for better visibility

1×1 adjoint(::Vector{SymbolicUtils.BasicSymbolic{Real}}) with eltype SymbolicUtils.BasicSymbolic{Real}:
 S_O(t)

(For confirmation, that these are both oxygen and not only equally named variables,see the respective documentation). In this case, generating the reactor with both processes is very easy:

@named react = CSTR(1000, [asm, aer]; initial_states=[1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1]) # Volume of 1000, initial_states of all 1 (the default of all 0 does yield infinity, check the model for yourself)

CSTR with ASM1 and Aeration 'react':
States (13): see states(react)
  S_ALK(t) [defaults to 1.0]: S_ALK
  S_I(t) [defaults to 1.0]: S_I
  S_ND(t) [defaults to 1.0]: S_ND
  S_NH(t) [defaults to 1.0]: S_NH
  ⋮
Inflows (1): see inflows(react)
  :react_reactor₊inflow
Outflows (1): see outflows(react)
  :react_reactor₊outflow
Exogenous Inputs (1): see exogenous_inputs(react)
  :aer₊k_La
Parameters (21): see parameters(react)
  react_reactor₊V [defaults to 1000]
  asm₊Y_H [defaults to 0.67]
  asm₊Y_A [defaults to 0.24]
  asm₊i_XB [defaults to 0.086]
  ⋮

Hereby, the equally named states are automatically connected and the rates of both processes added.

Then, the influent and the input for the aeration (the strength of aeration) are created:

@named influent = Influent([7.0, 30.0, 6.95, 31.56, 0.0, 0.0, 69.5, 0.0, 28.17, 51.2, 10.59, 0.0, 202.32], unknowns(states(asm)), flowrate = 18446) # BSM1 Stabilization
@named aeration_ctrl = ModelingToolkitStandardLibrary.Blocks.Constant(;k = 240)

And everything is connected:

eqs = [
  connect(outflows(influent, 1), inflows(react, 1))
  connect(exogenous_inputs(react, :aer₊k_La), aeration_ctrl.output)
]

\[ \begin{equation} \left[ \begin{array}{c} \mathrm{connect}\left( influent_{+}influent_{to\_outflowport_{+}outflow}, react_{+}react_{reactor_{+}inflow} \right) \\ \mathrm{connect}\left( react_{+}aer_{+}k_{La}, aeration_{ctrl_{+}output} \right) \\ \end{array} \right] \end{equation} \]

And finally, the system can be generated, simulated and plotted:

# Model generation
@named model = ODESystem(eqs, t, systems=[influent, react, aeration_ctrl])
model_simplified = structural_simplify(model)

# Simulation
ode_prob = ODEProblem(model_simplified, [], (0, 1))
sol = solve(ode_prob)

# Plotting, only oxygen
plot(sol, idxs=[inflows(react, 1).S_O, states(react).S_O], label=["Inflow" "Reactor"], title="Oxygen concentration in CSTR with ASM1 and aeration", legend=:topright)

And indeed, the oxygen concentration is not only consumed, but also added by the aeration. Thus both processes are active and it is even possible to look at the individual rates:

plot(sol, idxs=[rates(react).S_O, subsystems(react, :asm).rates.S_O, subsystems(react, :aer).rates.S_O], label=["Combined" "ASM" "Aeration"], title="Oxygen rate in CSTR with ASM1 and aeration", legend=:topright)

Adding a Process to the Default One

The second example of multiple processes, is concerning the case, where a process is to be added to the default ones. This is useful e.g. if a system is built where some of the reactors are aerated, but the others not. Then, the base process can be set as default and then the aeration process can be added to the individual reactors.

Plug and Play Example

using BioChemicalTreatment # This Package :)
using ModelingToolkit # Basic Modeling Framework. Takes e.g. care of simplification of the equations
using ModelingToolkitStandardLibrary # Standard library for MTK. Provides e.g. simple math operations (adder, product...), controllers (PID) and more
using DifferentialEquations # Solve the equations
using Plots # Plotting

# Default process: asm
@set_process asm = ASM1() # Use the default parameters

# Additional process
@named aer = Aeration() # The aeration process

# Reactor with the additional process
@named react = CSTR(1000, aer; initial_states=[1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1]) # Volume of 1000, initial_states of all 1 (the default of all 0 does yield infinity, check the model for yourself)

# Influent and aeration input
@named influent = Influent([7.0, 30.0, 6.95, 31.56, 0.0, 0.0, 69.5, 0.0, 28.17, 51.2, 10.59, 0.0, 202.32], flowrate = 18446) # BSM1 Stabilization
@named aeration_ctrl = ModelingToolkitStandardLibrary.Blocks.Constant(;k = 240)

# Connection
eqs = [
  connect(outflows(influent, 1), inflows(react, 1))
  connect(exogenous_inputs(react, :aer₊k_La), aeration_ctrl.output)
]

# Model generation
@named model = ODESystem(eqs, t, systems=[influent, react, aeration_ctrl])
model_simplified = structural_simplify(model)

# Simulation
ode_prob = ODEProblem(model_simplified, [], (0, 1))
sol = solve(ode_prob)

# Plotting, only oxygen
plot(sol, idxs=[inflows(react, 1).S_O, states(react).S_O], label=["Inflow" "Reactor"], title="Oxygen concentration in CSTR with ASM1 and aeration", legend=:topright)

# Plotting the rates
plot(sol, idxs=[rates(react).S_O, subsystems(react, :react_asm).rates.S_O, subsystems(react, :aer).rates.S_O], label=["Combined" "ASM" "Aeration"], title="Oxygen rate in CSTR with ASM1 and aeration", legend=:topright)

Step-by-Step Explanations

As always, the first step is the setup in loading the packages:

using BioChemicalTreatment # This Package :)
using ModelingToolkit # Basic Modeling Framework. Takes e.g. care of simplification of the equations
using ModelingToolkitStandardLibrary # Standard library for MTK. Provides e.g. simple math operations (adder, product...), controllers (PID) and more
using DifferentialEquations # Solve the equations
using Plots # Plotting

And then, the asm is set as default process:

@set_process asm = ASM1() # Use the default parameters

Then, the reactor with the additional aeration process is created. For this, just the aeration process is created and given to the reactor upon creation.

@named aer = Aeration() # The aeration process

# Reactor with additional aeration
@named react = CSTR(1000, aer; initial_states=[1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1]) # Volume of 1000, initial_states of all 1 (the default of all 0 does yield infinity, check the model for yourself)

Hereby, the default process automatically added, the equally named states are connected and the rates of both processes added.

Note

If the addition of the default process is not wanted, supply use_default_processes=false to the reactor call.

If the new process is to be added to all subsequent reactors, it can be added to the defaults using

@add_process aer = Aeration()

Note the difference between @add_process, which adds a process and @set_process which removes all other processes.

Further note that it is not easily possible to remove default processes, thus it is recommended to use the defaults only for the processes that are really the base processes occurring in every reactor. These should as well define all the components in the flow.

Then, it is proceeded exactly as in the example above:

# Influent and aeration input
@named influent = Influent([7.0, 30.0, 6.95, 31.56, 0.0, 0.0, 69.5, 0.0, 28.17, 51.2, 10.59, 0.0, 202.32], flowrate = 18446) # BSM1 Stabilization
@named aeration_ctrl = ModelingToolkitStandardLibrary.Blocks.Constant(;k = 240)

# Connection
eqs = [
  connect(outflows(influent, 1), inflows(react, 1))
  connect(exogenous_inputs(react, :aer₊k_La), aeration_ctrl.output)
]

# Model generation
@named model = ODESystem(eqs, t, systems=[influent, react, aeration_ctrl])
model_simplified = structural_simplify(model)

# Simulation
ode_prob = ODEProblem(model_simplified, [], (0, 1))
sol = solve(ode_prob)

# Plotting, only oxygen
plot(sol, idxs=[inflows(react, 1).S_O, states(react).S_O], label=["Inflow" "Reactor"], title="Oxygen concentration in CSTR with ASM1 and aeration", legend=:topright)

And the plot of the rates, note that now the asm subsystem of the reactor has a different name, i.e. is prepended by the name of the reactor:

# Plotting the rates
plot(sol, idxs=[rates(react).S_O, subsystems(react, :react_asm).rates.S_O, subsystems(react, :aer).rates.S_O], label=["Combined" "ASM" "Aeration"], title="Oxygen rate in CSTR with ASM1 and aeration", legend=:topright)

Combining Variables with Different Names

The third example of multiple processes, is concerning the case where the variables to be combined do not have the same names. This can occur if different model naming conventions are used, e.g. S_O and S_O2 are both commonly used for the dissolved oxygen concentration.

As seen in the two examples above, the aeration model uses S_O and the directly defined ASM1 (called using ASM1()) as well. However, the matrix-defined version of the ASM1 (called using Process("ASM1")) is using S_O2 instead (as the ASM3 does). In the example below it is shown how this case can be handled using the default processes. If it is rather desired without the defaults, the same as in the first example could be done.

Plug and Play Example

using BioChemicalTreatment # This Package :)
using ModelingToolkit # Basic Modeling Framework. Takes e.g. care of simplification of the equations
using ModelingToolkitStandardLibrary # Standard library for MTK. Provides e.g. simple math operations (adder, product...), controllers (PID) and more
using DifferentialEquations # Solve the equations
using Plots # Plotting

# Default process: asm
@set_process asm = Process("ASM1") # Use the default parameters

# Additional process
@named aer = Aeration() # The aeration process

# Set the state mapping to map both to S_O2
@set_state_mapping S_O2 = states(asm).S_O2, states(aer).S_O

# Reactor with the additional process
@named react = CSTR(1000, aer; initial_states=[1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1]) # Volume of 1000, initial_states of all 1 (the default of all 0 does yield infinity, check the model for yourself)

# Influent and aeration input
@named influent = Influent([7.0, 69.5, 6.95, 0.0, 31.56, 0.0, 0.0, 30.0, 202.32, 10.59, 0.0, 28.17, 0.0, 51.2], flowrate = 18446) # BSM1 Stabilization
@named aeration_ctrl = ModelingToolkitStandardLibrary.Blocks.Constant(;k = 240)

# Connection
eqs = [
  connect(outflows(influent, 1), inflows(react, 1))
  connect(exogenous_inputs(react, :aer₊k_La), aeration_ctrl.output)
]

# Model generation
@named model = ODESystem(eqs, t, systems=[influent, react, aeration_ctrl])
model_simplified = structural_simplify(model)

# Simulation
ode_prob = ODEProblem(model_simplified, [], (0, 1))
sol = solve(ode_prob)

# Plotting, only oxygen
plot(sol, idxs=[inflows(react, 1).S_O2, states(react).S_O2], label=["Inflow" "Reactor"], title="Oxygen concentration in CSTR with ASM1 and aeration", legend=:topright)

# Plotting the rates
plot(sol, idxs=[rates(react).S_O2, subsystems(react, :react_asm).rates.S_O2, subsystems(react, :aer).rates.S_O], label=["Combined" "ASM" "Aeration"], title="Oxygen rate in CSTR with ASM1 and aeration", legend=:topright)

Step-by-Step Explanations

As always, the first step is the setup in loading the packages:

using BioChemicalTreatment # This Package :)
using ModelingToolkit # Basic Modeling Framework. Takes e.g. care of simplification of the equations
using ModelingToolkitStandardLibrary # Standard library for MTK. Provides e.g. simple math operations (adder, product...), controllers (PID) and more
using DifferentialEquations # Solve the equations
using Plots # Plotting

And then, the matrix-defined asm is set as default process and the additional aeration process is created:

@set_process asm = Process("ASM1") # Use the default parameters
@named aer = Aeration() # The aeration process

Then, a mapping between the states S_O2 of the asm and S_O of the aeration is added, to let the program know that they are in fact the same states. To do so @set_state_mapping is used (alternatively @add_state_mapping is possible, with the only difference that the set- version clears the existing mappings).

# Set the state mapping to map both to S_O2
# Left of the equal sign is the name to use, on the right hand side the states that are mapped to it. Note that for the name also any other name not yet existing could have been used, and like this variables could be renamed.
@set_state_mapping S_O2 = states(asm).S_O2, states(aer).S_O

Then, the rest is exactly equal to the example above (possible for both, with and without default processes). The only difference is the name of the states and their order in the influent (it is different for the two ASM1 implementations, as the state names are different and they are ordered alphabetically to have it consistent):

# Reactor with the additional process
@named react = CSTR(1000, aer; initial_states=[1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1]) # Volume of 1000, initial_states of all 1 (the default of all 0 does yield infinity, check the model for yourself)

# Influent and aeration input
@named influent = Influent([7.0, 69.5, 6.95, 0.0, 31.56, 0.0, 0.0, 30.0, 202.32, 10.59, 0.0, 28.17, 0.0, 51.2], flowrate = 18446) # BSM1 Stabilization
@named aeration_ctrl = ModelingToolkitStandardLibrary.Blocks.Constant(;k = 240)

# Connection
eqs = [
  connect(outflows(influent, 1), inflows(react, 1))
  connect(exogenous_inputs(react, :aer₊k_La), aeration_ctrl.output)
]

# Model generation
@named model = ODESystem(eqs, t, systems=[influent, react, aeration_ctrl])
model_simplified = structural_simplify(model)

# Simulation
ode_prob = ODEProblem(model_simplified, [], (0, 1))
sol = solve(ode_prob)

# Plotting, only oxygen
plot(sol, idxs=[inflows(react, 1).S_O2, states(react).S_O2], label=["Inflow" "Reactor"], title="Oxygen concentration in CSTR with ASM1 and aeration", legend=:topright)

And the plot of the rates, note that now the asm subsystem of the reactor has a different name, i.e. is prepended by the name of the reactor:

# Plotting the rates
plot(sol, idxs=[rates(react).S_O2, subsystems(react, :react_asm).rates.S_O2, subsystems(react, :aer).rates.S_O], label=["Combined" "ASM" "Aeration"], title="Oxygen rate in CSTR with ASM1 and aeration", legend=:topright)