Low-Level Reactor with multiple Processes: ASM1 with Aeration

Note

This example is exactly the same as the corresponding the high level example. Thus the two interfaces can be compared.

For this system it can already be seen how the high-level interface simplifies the definition and one can imagine how it goes on for more complex systems.

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, combining two processes in a single reactor is exactly the point where the split between the Reactor and the Process such that the reaction rate is provided as an external signal to the Reactor comes in handy (see Low Level: Reactors vs. Processes for details). This is, because the reaction rates of the different processes can simply be added, and thus this can be done using external connections. Now this can be represented as:

For setting this up however, an adder block is needed. Luckily, this can be taken from the ModelingToolkitStandardLibrary as ModelingToolkitStandardLibrary.Blocks.Add where the compatibility systems provided in this library come in handy.

There is only one caveat left: The two processes do not need the same states and also do not provide the same rates. Instead, the Aeration needs only the dissolved oxygen state as input and also provides only a rate for the dissolved oxygen. This makes the combination more complex, but also this can be solved using the standard library compatibility systems, as they act as kind of a multiplexer: they split the ports into their individual signals and they can also assemble individual signals into the corresponding port.

Now, let's start the example of how to actually do this.

Setup

As always, the first step is the setup. Note that now also the ModelingToolkitStandardLibrary is loaded for the adder system and the bookkeeping variables are also introduced.

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

sys = ODESystem[]
eqs = Equation[]

Creating the Processes and Reactor

Then, as usual, first the processes (in this case two of them) are created.

@named asm = Process("ASM1") # Use the matrix-defined version with default parameters
@named aer = Aeration() # The aeration process
append!(sys, [asm, aer]) # Bookkeeping

@named react = CSTR(1000, unknowns(states(asm)); 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)
push!(sys, react) # Bookkeeping
nameof.(sys) # Show all added systems

# output

3-element Vector{Symbol}:
 :asm
 :aer
 :react

For the reactor, only the states of the ASM1 are used, as those superset the states of the aeration.

Creating the Systems for the combination

For the combination of the processes, a series of other helper systems is needed. First the state of the reactor needs to be split into single ports

@named reactor_state_to_single_ports = StateInputToRealOutput(unknowns(states(react))) # note the use of the states of the reactor here as the reactor state port is to be split
push!(sys, reactor_state_to_single_ports)
nameof.(sys) # Show all added systems

# output

4-element Vector{Symbol}:
 :asm
 :aer
 :react
 :reactor_state_to_single_ports

And then the aeration state needs to be reassembled from it:

@named single_ports_to_aeration_states = RealInputToStateOutput(unknowns(states(aer))) # Similarly here, aeration state is to be built -> use unknowns(states(aer))
push!(sys, single_ports_to_aeration_states)
nameof.(sys) # Show all added systems

# output

5-element Vector{Symbol}:
 :asm
 :aer
 :react
 :reactor_state_to_single_ports
 :single_ports_to_aeration_states

Now the state inputs for both processes are available (for the ASM1 directly the state from the reactor can be taken). So it's time to look at the other side of the rates. There, first the rate outputs of both processes have to be split apart

@named asm_rates_to_single_ports = ReactionInputToRealOutput(unknowns(rates(asm))) # asm rates is to be split -> unknowns(rates(asm))
@named aeration_rates_to_single_ports = ReactionInputToRealOutput(unknowns(rates(aer))) # aeration rates is to be split -> unknowns(rates(aer))
append!(sys, [asm_rates_to_single_ports, aeration_rates_to_single_ports])
nameof.(sys) # Show all added systems

# output

7-element Vector{Symbol}:
 :asm
 :aer
 :react
 :reactor_state_to_single_ports
 :single_ports_to_aeration_states
 :asm_rates_to_single_ports
 :aeration_rates_to_single_ports

then, an adder is needed to add the rates of the dissolved oxigen

@named add_aeration_rate = ModelingToolkitStandardLibrary.Blocks.Add()
push!(sys, add_aeration_rate)
nameof.(sys) # Show all added systems

# output

8-element Vector{Symbol}:
 :asm
 :aer
 :react
 :reactor_state_to_single_ports
 :single_ports_to_aeration_states
 :asm_rates_to_single_ports
 :aeration_rates_to_single_ports
 :add_aeration_rate

and finally the rates port has to be reassembled to be able to provide it to the reactor

@named single_ports_to_reaction_output = RealInputToReactionOutput(unknowns(rates(react))) # reactor rates is to be built -> unknowns(rates(react))
push!(sys, single_ports_to_reaction_output)
nameof.(sys) # Show all added systems

# output

9-element Vector{Symbol}:
 :asm
 :aer
 :react
 :reactor_state_to_single_ports
 :single_ports_to_aeration_states
 :asm_rates_to_single_ports
 :aeration_rates_to_single_ports
 :add_aeration_rate
 :single_ports_to_reaction_output

Creating the other needed systems

Now the most of the needed systems have been created, but two more are missing. Do you figure out which?

They are the influent and a aeration strength input for the aeration. The aeration process models the aeration, but it needs an input of the oxygen transfer coefficient, which allows it to be controlled (e.g. to hold a constant level of dissolved oxygen). See the docs of the aeration for details.

As influent and this control input, constant values are assumed for simplicity. For the influent the stabilization from the BSM1 is taken:

@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], unknowns(states(asm)), flowrate = 18446) # This is the stabilization influent from BSM1
push!(sys, influent)
nameof.(sys) # Show all added systems

# output

10-element Vector{Symbol}:
 :asm
 :aer
 :react
 :reactor_state_to_single_ports
 :single_ports_to_aeration_states
 :asm_rates_to_single_ports
 :aeration_rates_to_single_ports
 :add_aeration_rate
 :single_ports_to_reaction_output
 :influent

And for the control input, a value of 240 1/d is assumed. This port is a ModelingToolkitStandardLibrary.Blocks.RealInput as this port can directly be connected to the systems from there, e.g. the PID controller. Luckily, the library there also features a constant, which is used to specify this constant input:

@named aeration_ctrl = ModelingToolkitStandardLibrary.Blocks.Constant(;k = 240)
push!(sys, aeration_ctrl)
nameof.(sys) # Show all added systems

# output

11-element Vector{Symbol}:
 :asm
 :aer
 :react
 :reactor_state_to_single_ports
 :single_ports_to_aeration_states
 :asm_rates_to_single_ports
 :aeration_rates_to_single_ports
 :add_aeration_rate
 :single_ports_to_reaction_output
 :influent
 :aeration_ctrl

Connecting

It is started by connecting the current state of the reactor to the processes. As the ASM1 has the same states as the reactor this works as usual:

push!(eqs, connect(states(react), states(asm)))

For the Aeration, which only needs the dissolved oxygen concentration, first the port is split into its individual values to extract the dissolved oxygen, which is then fed to the system assembling the state port for the Aeration process and then to the process itself:

push!(eqs, connect(states(react), states(reactor_state_to_single_ports))) # First the reactor states to the splitting system
push!(eqs, connect(exogenous_outputs(reactor_state_to_single_ports, :S_O2), exogenous_inputs(single_ports_to_aeration_states, :S_O))) # Then the split up dissolved oxygen to the reassembly system
push!(eqs, connect(states(single_ports_to_aeration_states), states(aer))) # Then the reassembly system to the aeration states

Warning

Note the different names of the dissolved oxygen in the two processes, S_O vs S_O2. This is unfortunate, but comes from the fact that there are multiple naming conventions out there, the one from the original ASM (Henze et al., 2006) (used by the Aeration system) and the revised one from (Corominas et al., 2010) (Used by the matrix-defined ASM1). If the original ASM1 (using ASM1()) would have been used, both variables would be called S_O. The variable naming in this package is still work in progress and it might be unified at some point in the future. However, this does not hinder the combination of two processes.

This implies, however, that arbitrary process rates can be added and combined like this. Thus the user must pay attention that the correct rates are connected!

Then, the other side of adding and connecting the rates is tackled. For this, first the rate outputs of both systems have to be split up:

push!(eqs, connect(rates(asm), rates(asm_rates_to_single_ports)))
push!(eqs, connect(rates(aer), rates(aeration_rates_to_single_ports)))

And then, the dissolved oxygen rates can be summed up by connection to the addition system:

push!(eqs, connect(exogenous_outputs(asm_rates_to_single_ports, :S_O2), add_aeration_rate.input1))
push!(eqs, connect(exogenous_outputs(aeration_rates_to_single_ports, :S_O), add_aeration_rate.input2))

It is continued by assembling the combined rates port by connecting all single signals from the split ASM1 rates port and the summed up rate for dissolved oxygen to the reassembly system:

append!(eqs, [
    # First the summed up aeration rate
    connect(add_aeration_rate.output, exogenous_inputs(single_ports_to_reaction_output, :S_O2))
    # Then all others directly from the splitted asm rates 
    connect(exogenous_outputs(asm_rates_to_single_ports, :S_Alk), exogenous_inputs(single_ports_to_reaction_output, :S_Alk))
    connect(exogenous_outputs(asm_rates_to_single_ports, :S_B), exogenous_inputs(single_ports_to_reaction_output, :S_B))
    connect(exogenous_outputs(asm_rates_to_single_ports, :S_BN), exogenous_inputs(single_ports_to_reaction_output, :S_BN))
    connect(exogenous_outputs(asm_rates_to_single_ports, :S_N2), exogenous_inputs(single_ports_to_reaction_output, :S_N2))
    connect(exogenous_outputs(asm_rates_to_single_ports, :S_NHx), exogenous_inputs(single_ports_to_reaction_output, :S_NHx))
    connect(exogenous_outputs(asm_rates_to_single_ports, :S_NOx), exogenous_inputs(single_ports_to_reaction_output, :S_NOx))
    connect(exogenous_outputs(asm_rates_to_single_ports, :S_U), exogenous_inputs(single_ports_to_reaction_output, :S_U))
    connect(exogenous_outputs(asm_rates_to_single_ports, :XC_B), exogenous_inputs(single_ports_to_reaction_output, :XC_B))
    connect(exogenous_outputs(asm_rates_to_single_ports, :XC_BN), exogenous_inputs(single_ports_to_reaction_output, :XC_BN))
    connect(exogenous_outputs(asm_rates_to_single_ports, :X_ANO), exogenous_inputs(single_ports_to_reaction_output, :X_ANO))
    connect(exogenous_outputs(asm_rates_to_single_ports, :X_OHO), exogenous_inputs(single_ports_to_reaction_output, :X_OHO))
    connect(exogenous_outputs(asm_rates_to_single_ports, :X_UE), exogenous_inputs(single_ports_to_reaction_output, :X_UE))
    connect(exogenous_outputs(asm_rates_to_single_ports, :X_UInf), exogenous_inputs(single_ports_to_reaction_output, :X_UInf))
    ])

And finally, the assembled rates to the reactor:

push!(eqs, connect(rates(single_ports_to_reaction_output), rates(react)))

With this, the reactor and process are connected and ready. Now adding the influent and the control input, and then the system is ready to be simulated:

push!(eqs, connect(outflows(influent, 1), inflows(react, 1)))
push!(eqs, connect(exogenous_inputs(aer, :k_La), aeration_ctrl.output))

Simulation and plotting

Before simulating, the model must be built and simplified, as usual:

@named model = ODESystem(eqs, t, systems=sys)
model_simplified = structural_simplify(model)

Note

Incredibly, the simplification was able to simplify all of this to 14 differential equation, one for each state in the ASM1. And the dissolved oxygen has indeed the rates terms of the ASM and the aeration!

Check it for yourself:

Number of equations of the unsimplified model:

julia> length(equations(model))
168

Number of equations of the simplified model:

julia> length(equations(model_simplified))
14

Check the equations of the simplified model for yourself using

full_equations(model_simplified)

\[ \begin{align} \frac{\mathrm{d} \mathtt{react.states.S\_Alk}\left( t \right)}{\mathrm{d}t} &= \frac{ - \mathtt{asm.K\_O2OHO} \left( \frac{\left( 1 - \mathtt{asm.Y\_OHO} \right) \mathtt{asm.i\_ChargeSNOx}}{\mathtt{asm.Y\_OHO} \mathtt{asm.i\_NO3N2}} + \mathtt{asm.i\_ChargeSNHx} \mathtt{asm.i\_NXBio} \right) \mathtt{asm.m\_OHOMax} \mathtt{asm.n\_mOHOAx} \mathtt{react.states.S\_NOx}\left( t \right) \mathtt{react.states.X\_OHO}\left( t \right) \mathtt{react.states.S\_NHx}\left( t \right) \mathtt{react.states.S\_B}\left( t \right)}{\left( \mathtt{asm.K\_NHxOHO} + \mathtt{react.states.S\_NHx}\left( t \right) \right) \left( \mathtt{asm.K\_NOxOHO} + \mathtt{react.states.S\_NOx}\left( t \right) \right) \left( \mathtt{asm.K\_O2OHO} + \mathtt{react.states.S\_O2}\left( t \right) \right) \left( \mathtt{asm.K\_SBOHO} + \mathtt{react.states.S\_B}\left( t \right) \right)} + \frac{ - \left( \frac{ - \mathtt{asm.i\_ChargeSNOx}}{\mathtt{asm.Y\_ANO}} - \mathtt{asm.i\_ChargeSNHx} \left( - \mathtt{asm.i\_NXBio} + \frac{-1}{\mathtt{asm.Y\_ANO}} \right) \right) \mathtt{asm.m\_ANOMax} \mathtt{react.states.S\_O2}\left( t \right) \mathtt{react.states.X\_ANO}\left( t \right) \mathtt{react.states.S\_NHx}\left( t \right)}{\left( \mathtt{asm.K\_NHxANO} + \mathtt{react.states.S\_NHx}\left( t \right) \right) \left( \mathtt{asm.K\_O2ANO} + \mathtt{react.states.S\_O2}\left( t \right) \right)} + \frac{ - \mathtt{asm.i\_ChargeSNHx} \mathtt{asm.i\_NXBio} \mathtt{asm.m\_OHOMax} \mathtt{react.states.S\_O2}\left( t \right) \mathtt{react.states.X\_OHO}\left( t \right) \mathtt{react.states.S\_NHx}\left( t \right) \mathtt{react.states.S\_B}\left( t \right)}{\left( \mathtt{asm.K\_NHxOHO} + \mathtt{react.states.S\_NHx}\left( t \right) \right) \left( \mathtt{asm.K\_O2OHO} + \mathtt{react.states.S\_O2}\left( t \right) \right) \left( \mathtt{asm.K\_SBOHO} + \mathtt{react.states.S\_B}\left( t \right) \right)} + \frac{\mathtt{influent.S\_Alk\_const.k} \mathtt{influent.q\_const.k}}{\mathtt{react.V}} + \frac{ - \mathtt{influent.q\_const.k} \mathtt{react.states.S\_Alk}\left( t \right)}{\mathtt{react.V}} + \mathtt{asm.i\_ChargeSNHx} \mathtt{asm.q\_am} \mathtt{react.states.X\_OHO}\left( t \right) \mathtt{react.states.S\_BN}\left( t \right) \\ \frac{\mathrm{d} \mathtt{react.states.S\_B}\left( t \right)}{\mathrm{d}t} &= \frac{\mathtt{influent.S\_B\_const.k} \mathtt{influent.q\_const.k}}{\mathtt{react.V}} + \frac{ - \mathtt{asm.m\_OHOMax} \mathtt{react.states.S\_O2}\left( t \right) \mathtt{react.states.X\_OHO}\left( t \right) \mathtt{react.states.S\_NHx}\left( t \right) \mathtt{react.states.S\_B}\left( t \right)}{\left( \mathtt{asm.K\_NHxOHO} + \mathtt{react.states.S\_NHx}\left( t \right) \right) \left( \mathtt{asm.K\_O2OHO} + \mathtt{react.states.S\_O2}\left( t \right) \right) \left( \mathtt{asm.K\_SBOHO} + \mathtt{react.states.S\_B}\left( t \right) \right) \mathtt{asm.Y\_OHO}} + \frac{ - \mathtt{influent.q\_const.k} \mathtt{react.states.S\_B}\left( t \right)}{\mathtt{react.V}} + \frac{\mathtt{asm.q\_XCBSBhyd} \mathtt{react.states.XC\_B}\left( t \right) \left( \frac{\mathtt{asm.K\_O2OHO} \mathtt{asm.n\_qhydAx} \mathtt{react.states.S\_NOx}\left( t \right)}{\left( \mathtt{asm.K\_NOxOHO} + \mathtt{react.states.S\_NOx}\left( t \right) \right) \left( \mathtt{asm.K\_O2OHO} + \mathtt{react.states.S\_O2}\left( t \right) \right)} + \frac{\mathtt{react.states.S\_O2}\left( t \right)}{\mathtt{asm.K\_O2OHO} + \mathtt{react.states.S\_O2}\left( t \right)} \right)}{\mathtt{asm.K\_XCBhyd} + \frac{\mathtt{react.states.XC\_B}\left( t \right)}{\mathtt{react.states.X\_OHO}\left( t \right)}} + \frac{ - \mathtt{asm.K\_O2OHO} \mathtt{asm.m\_OHOMax} \mathtt{asm.n\_mOHOAx} \mathtt{react.states.S\_NOx}\left( t \right) \mathtt{react.states.X\_OHO}\left( t \right) \mathtt{react.states.S\_NHx}\left( t \right) \mathtt{react.states.S\_B}\left( t \right)}{\left( \mathtt{asm.K\_NHxOHO} + \mathtt{react.states.S\_NHx}\left( t \right) \right) \left( \mathtt{asm.K\_NOxOHO} + \mathtt{react.states.S\_NOx}\left( t \right) \right) \left( \mathtt{asm.K\_O2OHO} + \mathtt{react.states.S\_O2}\left( t \right) \right) \left( \mathtt{asm.K\_SBOHO} + \mathtt{react.states.S\_B}\left( t \right) \right) \mathtt{asm.Y\_OHO}} \\ \frac{\mathrm{d} \mathtt{react.states.S\_BN}\left( t \right)}{\mathrm{d}t} &= \frac{\mathtt{asm.q\_XCBSBhyd} \left( \frac{\mathtt{asm.K\_O2OHO} \mathtt{asm.n\_qhydAx} \mathtt{react.states.S\_NOx}\left( t \right)}{\left( \mathtt{asm.K\_NOxOHO} + \mathtt{react.states.S\_NOx}\left( t \right) \right) \left( \mathtt{asm.K\_O2OHO} + \mathtt{react.states.S\_O2}\left( t \right) \right)} + \frac{\mathtt{react.states.S\_O2}\left( t \right)}{\mathtt{asm.K\_O2OHO} + \mathtt{react.states.S\_O2}\left( t \right)} \right) \mathtt{react.states.XC\_BN}\left( t \right)}{\mathtt{asm.K\_XCBhyd} + \frac{\mathtt{react.states.XC\_B}\left( t \right)}{\mathtt{react.states.X\_OHO}\left( t \right)}} + \frac{\mathtt{influent.S\_BN\_const.k} \mathtt{influent.q\_const.k}}{\mathtt{react.V}} + \frac{ - \mathtt{influent.q\_const.k} \mathtt{react.states.S\_BN}\left( t \right)}{\mathtt{react.V}} - \mathtt{asm.q\_am} \mathtt{react.states.X\_OHO}\left( t \right) \mathtt{react.states.S\_BN}\left( t \right) \\ \frac{\mathrm{d} \mathtt{react.states.S\_N2}\left( t \right)}{\mathrm{d}t} &= \frac{\mathtt{asm.K\_O2OHO} \left( 1 - \mathtt{asm.Y\_OHO} \right) \mathtt{asm.m\_OHOMax} \mathtt{asm.n\_mOHOAx} \mathtt{react.states.S\_NOx}\left( t \right) \mathtt{react.states.X\_OHO}\left( t \right) \mathtt{react.states.S\_NHx}\left( t \right) \mathtt{react.states.S\_B}\left( t \right)}{\left( \mathtt{asm.K\_NHxOHO} + \mathtt{react.states.S\_NHx}\left( t \right) \right) \left( \mathtt{asm.K\_NOxOHO} + \mathtt{react.states.S\_NOx}\left( t \right) \right) \left( \mathtt{asm.K\_O2OHO} + \mathtt{react.states.S\_O2}\left( t \right) \right) \left( \mathtt{asm.K\_SBOHO} + \mathtt{react.states.S\_B}\left( t \right) \right) \mathtt{asm.Y\_OHO} \mathtt{asm.i\_NO3N2}} + \frac{ - \mathtt{influent.q\_const.k} \mathtt{react.states.S\_N2}\left( t \right)}{\mathtt{react.V}} + \frac{\mathtt{influent.S\_N2\_const.k} \mathtt{influent.q\_const.k}}{\mathtt{react.V}} \\ \frac{\mathrm{d} \mathtt{react.states.S\_NHx}\left( t \right)}{\mathrm{d}t} &= \frac{ - \mathtt{influent.q\_const.k} \mathtt{react.states.S\_NHx}\left( t \right)}{\mathtt{react.V}} + \frac{ - \mathtt{asm.i\_NXBio} \mathtt{asm.m\_OHOMax} \mathtt{react.states.S\_O2}\left( t \right) \mathtt{react.states.X\_OHO}\left( t \right) \mathtt{react.states.S\_NHx}\left( t \right) \mathtt{react.states.S\_B}\left( t \right)}{\left( \mathtt{asm.K\_NHxOHO} + \mathtt{react.states.S\_NHx}\left( t \right) \right) \left( \mathtt{asm.K\_O2OHO} + \mathtt{react.states.S\_O2}\left( t \right) \right) \left( \mathtt{asm.K\_SBOHO} + \mathtt{react.states.S\_B}\left( t \right) \right)} + \frac{ - \mathtt{asm.K\_O2OHO} \mathtt{asm.i\_NXBio} \mathtt{asm.m\_OHOMax} \mathtt{asm.n\_mOHOAx} \mathtt{react.states.S\_NOx}\left( t \right) \mathtt{react.states.X\_OHO}\left( t \right) \mathtt{react.states.S\_NHx}\left( t \right) \mathtt{react.states.S\_B}\left( t \right)}{\left( \mathtt{asm.K\_NHxOHO} + \mathtt{react.states.S\_NHx}\left( t \right) \right) \left( \mathtt{asm.K\_NOxOHO} + \mathtt{react.states.S\_NOx}\left( t \right) \right) \left( \mathtt{asm.K\_O2OHO} + \mathtt{react.states.S\_O2}\left( t \right) \right) \left( \mathtt{asm.K\_SBOHO} + \mathtt{react.states.S\_B}\left( t \right) \right)} + \frac{\mathtt{influent.S\_NHx\_const.k} \mathtt{influent.q\_const.k}}{\mathtt{react.V}} + \frac{\left( - \mathtt{asm.i\_NXBio} + \frac{-1}{\mathtt{asm.Y\_ANO}} \right) \mathtt{asm.m\_ANOMax} \mathtt{react.states.S\_O2}\left( t \right) \mathtt{react.states.X\_ANO}\left( t \right) \mathtt{react.states.S\_NHx}\left( t \right)}{\left( \mathtt{asm.K\_NHxANO} + \mathtt{react.states.S\_NHx}\left( t \right) \right) \left( \mathtt{asm.K\_O2ANO} + \mathtt{react.states.S\_O2}\left( t \right) \right)} + \mathtt{asm.q\_am} \mathtt{react.states.X\_OHO}\left( t \right) \mathtt{react.states.S\_BN}\left( t \right) \\ \frac{\mathrm{d} \mathtt{react.states.S\_NOx}\left( t \right)}{\mathrm{d}t} &= \frac{\mathtt{asm.K\_O2OHO} \left( -1 + \mathtt{asm.Y\_OHO} \right) \mathtt{asm.m\_OHOMax} \mathtt{asm.n\_mOHOAx} \mathtt{react.states.S\_NOx}\left( t \right) \mathtt{react.states.X\_OHO}\left( t \right) \mathtt{react.states.S\_NHx}\left( t \right) \mathtt{react.states.S\_B}\left( t \right)}{\left( \mathtt{asm.K\_NHxOHO} + \mathtt{react.states.S\_NHx}\left( t \right) \right) \left( \mathtt{asm.K\_NOxOHO} + \mathtt{react.states.S\_NOx}\left( t \right) \right) \left( \mathtt{asm.K\_O2OHO} + \mathtt{react.states.S\_O2}\left( t \right) \right) \left( \mathtt{asm.K\_SBOHO} + \mathtt{react.states.S\_B}\left( t \right) \right) \mathtt{asm.Y\_OHO} \mathtt{asm.i\_NO3N2}} + \frac{\mathtt{influent.S\_NOx\_const.k} \mathtt{influent.q\_const.k}}{\mathtt{react.V}} + \frac{ - \mathtt{influent.q\_const.k} \mathtt{react.states.S\_NOx}\left( t \right)}{\mathtt{react.V}} + \frac{\mathtt{asm.m\_ANOMax} \mathtt{react.states.S\_O2}\left( t \right) \mathtt{react.states.X\_ANO}\left( t \right) \mathtt{react.states.S\_NHx}\left( t \right)}{\left( \mathtt{asm.K\_NHxANO} + \mathtt{react.states.S\_NHx}\left( t \right) \right) \left( \mathtt{asm.K\_O2ANO} + \mathtt{react.states.S\_O2}\left( t \right) \right) \mathtt{asm.Y\_ANO}} \\ \frac{\mathrm{d} \mathtt{react.states.S\_O2}\left( t \right)}{\mathrm{d}t} &= \frac{\mathtt{influent.S\_O2\_const.k} \mathtt{influent.q\_const.k}}{\mathtt{react.V}} + \frac{ - \mathtt{influent.q\_const.k} \mathtt{react.states.S\_O2}\left( t \right)}{\mathtt{react.V}} + \mathtt{add\_aeration\_rate.k1} \left( \frac{\left( -1 + \mathtt{asm.Y\_OHO} \right) \mathtt{asm.m\_OHOMax} \mathtt{react.states.S\_O2}\left( t \right) \mathtt{react.states.X\_OHO}\left( t \right) \mathtt{react.states.S\_NHx}\left( t \right) \mathtt{react.states.S\_B}\left( t \right)}{\left( \mathtt{asm.K\_NHxOHO} + \mathtt{react.states.S\_NHx}\left( t \right) \right) \left( \mathtt{asm.K\_O2OHO} + \mathtt{react.states.S\_O2}\left( t \right) \right) \left( \mathtt{asm.K\_SBOHO} + \mathtt{react.states.S\_B}\left( t \right) \right) \mathtt{asm.Y\_OHO}} + \frac{\left( \mathtt{asm.Y\_ANO} + \mathtt{asm.i\_CODNO3} \right) \mathtt{asm.m\_ANOMax} \mathtt{react.states.S\_O2}\left( t \right) \mathtt{react.states.X\_ANO}\left( t \right) \mathtt{react.states.S\_NHx}\left( t \right)}{\left( \mathtt{asm.K\_NHxANO} + \mathtt{react.states.S\_NHx}\left( t \right) \right) \left( \mathtt{asm.K\_O2ANO} + \mathtt{react.states.S\_O2}\left( t \right) \right) \mathtt{asm.Y\_ANO}} \right) + \mathtt{add\_aeration\_rate.k2} \mathtt{aeration\_ctrl.k} \left( \mathtt{aer.S\_O\_max} - \mathtt{react.states.S\_O2}\left( t \right) \right) \\ \frac{\mathrm{d} \mathtt{react.states.S\_U}\left( t \right)}{\mathrm{d}t} &= \frac{ - \mathtt{influent.q\_const.k} \mathtt{react.states.S\_U}\left( t \right)}{\mathtt{react.V}} + \frac{\mathtt{influent.S\_U\_const.k} \mathtt{influent.q\_const.k}}{\mathtt{react.V}} \\ \frac{\mathrm{d} \mathtt{react.states.XC\_B}\left( t \right)}{\mathrm{d}t} &= \frac{ - \mathtt{influent.q\_const.k} \mathtt{react.states.XC\_B}\left( t \right)}{\mathtt{react.V}} + \frac{ - \mathtt{asm.q\_XCBSBhyd} \mathtt{react.states.XC\_B}\left( t \right) \left( \frac{\mathtt{asm.K\_O2OHO} \mathtt{asm.n\_qhydAx} \mathtt{react.states.S\_NOx}\left( t \right)}{\left( \mathtt{asm.K\_NOxOHO} + \mathtt{react.states.S\_NOx}\left( t \right) \right) \left( \mathtt{asm.K\_O2OHO} + \mathtt{react.states.S\_O2}\left( t \right) \right)} + \frac{\mathtt{react.states.S\_O2}\left( t \right)}{\mathtt{asm.K\_O2OHO} + \mathtt{react.states.S\_O2}\left( t \right)} \right)}{\mathtt{asm.K\_XCBhyd} + \frac{\mathtt{react.states.XC\_B}\left( t \right)}{\mathtt{react.states.X\_OHO}\left( t \right)}} + \frac{\mathtt{influent.XC\_B\_const.k} \mathtt{influent.q\_const.k}}{\mathtt{react.V}} + \mathtt{asm.b\_ANO} \left( 1 - \mathtt{asm.f\_XUBiolys} \right) \mathtt{react.states.X\_ANO}\left( t \right) + \mathtt{asm.b\_OHO} \left( 1 - \mathtt{asm.f\_XUBiolys} \right) \mathtt{react.states.X\_OHO}\left( t \right) \\ \frac{\mathrm{d} \mathtt{react.states.XC\_BN}\left( t \right)}{\mathrm{d}t} &= \frac{ - \mathtt{asm.q\_XCBSBhyd} \left( \frac{\mathtt{asm.K\_O2OHO} \mathtt{asm.n\_qhydAx} \mathtt{react.states.S\_NOx}\left( t \right)}{\left( \mathtt{asm.K\_NOxOHO} + \mathtt{react.states.S\_NOx}\left( t \right) \right) \left( \mathtt{asm.K\_O2OHO} + \mathtt{react.states.S\_O2}\left( t \right) \right)} + \frac{\mathtt{react.states.S\_O2}\left( t \right)}{\mathtt{asm.K\_O2OHO} + \mathtt{react.states.S\_O2}\left( t \right)} \right) \mathtt{react.states.XC\_BN}\left( t \right)}{\mathtt{asm.K\_XCBhyd} + \frac{\mathtt{react.states.XC\_B}\left( t \right)}{\mathtt{react.states.X\_OHO}\left( t \right)}} + \frac{ - \mathtt{influent.q\_const.k} \mathtt{react.states.XC\_BN}\left( t \right)}{\mathtt{react.V}} + \frac{\mathtt{influent.XC\_BN\_const.k} \mathtt{influent.q\_const.k}}{\mathtt{react.V}} + \mathtt{asm.b\_ANO} \left( \mathtt{asm.i\_NXBio} - \mathtt{asm.f\_XUBiolys} \mathtt{asm.i\_NXUE} \right) \mathtt{react.states.X\_ANO}\left( t \right) + \mathtt{asm.b\_OHO} \left( \mathtt{asm.i\_NXBio} - \mathtt{asm.f\_XUBiolys} \mathtt{asm.i\_NXUE} \right) \mathtt{react.states.X\_OHO}\left( t \right) \\ \frac{\mathrm{d} \mathtt{react.states.X\_ANO}\left( t \right)}{\mathrm{d}t} &= \frac{\mathtt{asm.m\_ANOMax} \mathtt{react.states.S\_O2}\left( t \right) \mathtt{react.states.X\_ANO}\left( t \right) \mathtt{react.states.S\_NHx}\left( t \right)}{\left( \mathtt{asm.K\_NHxANO} + \mathtt{react.states.S\_NHx}\left( t \right) \right) \left( \mathtt{asm.K\_O2ANO} + \mathtt{react.states.S\_O2}\left( t \right) \right)} + \frac{ - \mathtt{influent.q\_const.k} \mathtt{react.states.X\_ANO}\left( t \right)}{\mathtt{react.V}} + \frac{\mathtt{influent.X\_ANO\_const.k} \mathtt{influent.q\_const.k}}{\mathtt{react.V}} - \mathtt{asm.b\_ANO} \mathtt{react.states.X\_ANO}\left( t \right) \\ \frac{\mathrm{d} \mathtt{react.states.X\_OHO}\left( t \right)}{\mathrm{d}t} &= \frac{\mathtt{influent.X\_OHO\_const.k} \mathtt{influent.q\_const.k}}{\mathtt{react.V}} + \frac{\mathtt{asm.K\_O2OHO} \mathtt{asm.m\_OHOMax} \mathtt{asm.n\_mOHOAx} \mathtt{react.states.S\_NOx}\left( t \right) \mathtt{react.states.X\_OHO}\left( t \right) \mathtt{react.states.S\_NHx}\left( t \right) \mathtt{react.states.S\_B}\left( t \right)}{\left( \mathtt{asm.K\_NHxOHO} + \mathtt{react.states.S\_NHx}\left( t \right) \right) \left( \mathtt{asm.K\_NOxOHO} + \mathtt{react.states.S\_NOx}\left( t \right) \right) \left( \mathtt{asm.K\_O2OHO} + \mathtt{react.states.S\_O2}\left( t \right) \right) \left( \mathtt{asm.K\_SBOHO} + \mathtt{react.states.S\_B}\left( t \right) \right)} + \frac{ - \mathtt{influent.q\_const.k} \mathtt{react.states.X\_OHO}\left( t \right)}{\mathtt{react.V}} + \frac{\mathtt{asm.m\_OHOMax} \mathtt{react.states.S\_O2}\left( t \right) \mathtt{react.states.X\_OHO}\left( t \right) \mathtt{react.states.S\_NHx}\left( t \right) \mathtt{react.states.S\_B}\left( t \right)}{\left( \mathtt{asm.K\_NHxOHO} + \mathtt{react.states.S\_NHx}\left( t \right) \right) \left( \mathtt{asm.K\_O2OHO} + \mathtt{react.states.S\_O2}\left( t \right) \right) \left( \mathtt{asm.K\_SBOHO} + \mathtt{react.states.S\_B}\left( t \right) \right)} - \mathtt{asm.b\_OHO} \mathtt{react.states.X\_OHO}\left( t \right) \\ \frac{\mathrm{d} \mathtt{react.states.X\_UE}\left( t \right)}{\mathrm{d}t} &= \frac{\mathtt{influent.X\_UE\_const.k} \mathtt{influent.q\_const.k}}{\mathtt{react.V}} + \frac{ - \mathtt{influent.q\_const.k} \mathtt{react.states.X\_UE}\left( t \right)}{\mathtt{react.V}} + \mathtt{asm.b\_ANO} \mathtt{asm.f\_XUBiolys} \mathtt{react.states.X\_ANO}\left( t \right) + \mathtt{asm.b\_OHO} \mathtt{asm.f\_XUBiolys} \mathtt{react.states.X\_OHO}\left( t \right) \\ \frac{\mathrm{d} \mathtt{react.states.X\_UInf}\left( t \right)}{\mathrm{d}t} &= \frac{\mathtt{influent.X\_UInf\_const.k} \mathtt{influent.q\_const.k}}{\mathtt{react.V}} + \frac{ - \mathtt{influent.q\_const.k} \mathtt{react.states.X\_UInf}\left( t \right)}{\mathtt{react.V}} \end{align} \]

Then, the problem can be built and solved:

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

And the solution plotted:

plot(sol, title="CSTR with ASM1 and aeration", legend=:topright)

Indeed it is visible that there is some oxygen in the tank and as well some heterotrophs (which need oxygen to grow) seem to be present. But, let's have a look at the oxygen concentrations:

plot(sol, idxs=[inflows(react, 1).S_O2, outflows(react, 1).S_O2], label=["Reactor Inflow" "Reactor Outflow"], title="Dissolved oxygen concentration", legend=:right)

This makes it evident, that the aeration works and the oxygen is supplied. However, it might be interesting to have a look at the rate of the dissolved oxygen, the total one in the reactor and as well the contributions of the two processes. It is also possible to plot those:

plot(sol, idxs=[rates(react).S_O2, rates(asm).S_O2, rates(aer).S_O], label=["Total Rate" "Rate from ASM1" "Rate from Aeration"], title="Dissolved oxygen reaction rates", legend=:topright)

And as expected, the two individual rates sum up to the total rate in the reactor.