Reactor with Multiple Processes: ASM1/ASM3 with Aeration

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

TikzPicture(
     """;
    options = ">=stealth, auto, node distance=2cm, scale=2, background rectangle/.style={fill=white}, show background rectangle",

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 covers the case where two processes are provided to a reactor and two variables have the same name.

Plug and Play Example

using BioChemicalTreatment # This package :)
using ModelingToolkit # Basic modeling framework. Takes e.g. care of simplification of the equations
using ModelingToolkitStandardLibrary.Blocks # 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 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 = Constant(; k = 240)

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

# Model generation
@mtkcompile model = System(eqs, t, systems = [influent, react, aeration_ctrl])

# Simulation
ode_prob = ODEProblem(model, [], (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, :processes).aer.rates.S_O,
        subsystems(react, :processes).asm.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 by loading the packages:

using BioChemicalTreatment # This package :)
using ModelingToolkit # Basic modeling framework. Takes e.g. care of simplification of the equations
using ModelingToolkitStandardLibrary.Blocks # 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)

The `unknowns` function

The function unknowns, is used to show all the unknown variables of a system. All the states of the system are part of it, but there might be other variables as well (e.g. helper variables for renaming or connecting). Thus the function is not called states directly. But by using unknowns(states(sys)) we get only the states, as the states function returns a system with only the states as unknowns.

(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 yourself)

@named react = CSTR(

\[ \begin{equation} \left[ \begin{array}{c} \mathrm{connect}\left( processes_{+}rate_{assemble_{+}reactionoutput}, reactor_{+}rates \right) \\ \mathrm{connect}\left( reactor_{+}states, processes_{+}state_{split_{+}stateinput} \right) \\ \end{array} \right] \end{equation} \]

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 = Constant(; k = 240)

And everything is connected:

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

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

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

# Model generation
@mtkcompile model = System(eqs, t, systems = [influent, react, aeration_ctrl])

# Simulation
ode_prob = ODEProblem(model, [], (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, :processes).aer.rates.S_O,
        subsystems(react, :processes).asm.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 added to the default ones. This is useful, e.g., if a system is built where some of the reactors are aerated, but 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.Blocks # 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 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 = Constant(; k = 240)

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

# Model generation
@mtkcompile model = System(eqs, t, systems = [influent, react, aeration_ctrl])

# Simulation
ode_prob = ODEProblem(model, [], (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, :processes).asm.rates.S_O,
        subsystems(react, :processes).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 by loading the packages:

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

Next, the ASM is set as default process:

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

Then, the reactor with the additional aeration process is defined. For this, just the aeration process needs to be created and given to the reactor.

@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 is automatically added, the equally named states are connected and the rates of both processes summed up.

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 individual default processes, i.e. you can only use all or none of them (and add additional ones). However, it is not possible to use defaults, except single processes (e.g. if ASM1 and aeration are set as default processes, it is not possible to remove only the aeration from the created reactor. The only option is to not use the defaults/remove them and add only ASM1 to the reactor). 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 = Constant(; k = 240)

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

# Model generation
@mtkcompile model = System(eqs, t, systems = [influent, react, aeration_ctrl])

# Simulation
ode_prob = ODEProblem(model, [], (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 finally, 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, :processes).aer.rates.S_O,
        subsystems(react, :processes).asm.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 shows the case in which 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, ASM3 (called using Process("ASM3")) is using S_O2 instead of S_O. In the example below it is shown how this case can be handled using the default processes. To do the same without default processes, the first example needs to be adjusted accordingly.

Plug and Play Example

using BioChemicalTreatment # This package :)
using ModelingToolkit # Basic modeling framework. Takes e.g. care of simplification of the equations
using ModelingToolkitStandardLibrary.Blocks # 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("ASM3") # 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 yourself)

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

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

# Model generation
@mtkcompile model = System(eqs, t, systems = [influent, react, aeration_ctrl])

# Simulation
ode_prob = ODEProblem(model, [], (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 ASM3 and aeration", legend = :topright)

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

Step-by-Step Explanations

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

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

Next, the matrix-defined ASM3 is set as default process and the additional aeration process is created:

@set_process asm = Process("ASM3") # 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 introduced to let the model 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

The rest is 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:

# 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, 30.0, 0.0, 31.56, 6.95, 0.0, 69.5, 0.0, 28.17, 51.20, 202.32, 216.393, 0.0],
    flowrate = 18446) # BSM1 Stabilization
@named aeration_ctrl = Constant(; k = 240)

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

# Model generation
@mtkcompile model = System(eqs, t, systems = [influent, react, aeration_ctrl])

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

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

And finally, the plot of the rates (note that now the state name of the reactor has changed, however the ones of the process rates are left untouched, but still combined to the effective one):

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