Usage

This section gives an introduction to the package, by explaining its basic principles and specialities. These principles are then applied directly in a small example.

BioChemicalTreatment is a Julia package for the simulation of biochemical treatment processes (such as activated sludge processes for wastewater treatment), but it is not limited to them and provides also related physical processes. As such, it allows the user to specify a set of reactors with associated processes and other systems, and connecting them to form a model of a complex system.

To build the models, the ProcessSimulator of BioChemicalTreatment builds upon ModelingToolkit, which provides several convenient features such as automatic simplification of equations for faster solving. Further, this allowed BioChemicalTreatment to be built compatible with the ModelingToolkitStandardLibrary, from where a series of controllers and other helper systems can be used.

The first section gives a general overview of the setup of this package and the caveats when using it, afterwards a small Example is provided on which the use of the package is introduced step by step. Note that this section remains on a high level, explaining what would usually be used when building a model. See the section on the low level interface to have a in depth description of all the elements and of a more flexible (but also less convenient) way to use them.

General Introduction

The general concept behind this package is that provided components (called Systems) are connected to form the model. This model is then used to build the desired problem (e.g. an ODEProblem to integrate the differential equations through time or a SteadyStateProblem to solve for the steady state), which can then be solved automatically by calling solve. In the following, an overview over the types of systems is provided, the different types of connectors are explained, and the detailed process of setting up a basic model is explained in the section Example.

Types of Systems

This package contains several types of systems, each with different properties and purposes. Here the most commonly used ones are described:

• Reactor: Defines the differential equations for the mass balance in the reactor. Has (possibly multiple) inflows and outflows and usually encapsulates one or more Processes.
• Process: Computes the reaction rates to provide them to a reactor. It is usually wrapped by the reactor to which the states are provided.
• FlowElement: Additional systems required for connecting the systems and direct the flow. For example, providing an influent, unifying or separating the flow or an ideal clarifier. These elements are without dynamics i.e. they do not contain differential equations.

All of these systems can be used, together with the ones of the ModelingToolkitStandardLibrary (which provides e.g. basic mathematical functions, controllers...) to build models of complex systems.

Connectors

All the systems in this package have connectors. A connector can be thought of as a port with a series of variables to connect within. For example, a continuously stirred tank reactor (CSTR) has an inflow port which carries the flow rate and the concentrations of the state variables to the inflow. Similarly, it has an outflow which carries the same state variables, just for the outflow of the reactor. These connectors can then be used to conveniently connect (hence the name) two systems. For example, to connect two reactors with outflow outflow1, resp. inflow inflow2 one could use the following code:

connect(outflow1, inflow2)

The classical equivalent to this would be to add equations for the flow rate and the concentration of every state variable in the flow. One can imagine that this would become very tedious for models with many reactors and several compounds involved. Thus the syntax using the connectors is very convenient. Note that the order of inflow1 and outflow2 does not matter and the statement

connect(inflow2, outflow1)

is equivalent to the one above (although probably not so intuitive). For an example on how the connectors are actually used when building a model, see the example below where the connections are made in this subsection.

The important connector types in this package are:

• Flows: Carry a flow rate and concentrations of the state variable in it. Generally systems have separate inflows or outflows (or both).
• Exogenous Inputs: A system can have a set of exogenous_inputs, which are additional inputs to the system which are not in any of the above categories. Often, they can be considered control inputs. An example would be the pumped flow rate of a FlowPump.
• Exogenous Outputs: A system can have a set of exogenous_outputs, which are additional outputs of the system which are not in any of the above categories. An example would be the measurement output of a Sensor.

As a system can have multiple inflows, outflows, exogenous_inputs or exogenous_outputs, these connectors have names by which they can be accessed.

If one would like to investigate single signals in such a port (e.g. one compound in a flow as control input), it is possible to use the Sensor to extract a single or multiple signals into one or multiple RealOutputs

The conversion to the RealInputs and RealOutputs has been chosen, as with these ports, one can directly work with the systems provided by the blocks module of the ModelingToolkitStandardLibrary, which provides useful systems such as PID-controllers or mathematical transformations.

Example

The following example shows how to connect the individual systems to a combined model.

As an example we choose a disinfection and use the predefined process OzoneDisinfection with two state variables (the soluble ozone concentration SO3 and the bacteria concentration XB).

For this simple example, a single CSTR fed by a constant influent of both bacteria and ozone is used. In the CSTR, the provided process for disinfection with ozone is employed. This is schematically depicted as:

Plug and Play version of the example

using BioChemicalTreatment # This Package :)
using ModelingToolkit # Basic Modeling Framework. Takes e.g. care of simplification of the equations
using DifferentialEquations # Solve the equations
using Plots # Plotting

# Set the used process
@set_process process = OzoneDisinfection()

# Generate the systems
@named reactor = CSTR(1000)
@named influent = Influent([0.5, 1]; flowrate = 500)

# Connection
connect_equations = [connect(outflows(influent)[1], inflows(reactor)[1])]

# Preparation for Simulation
@named model = ODESystem(connect_equations, t; systems=[influent, reactor])
model_simplified = structural_simplify(model)

# Setady State solution
steadystate_prob = SteadyStateProblem(model_simplified, [])
steadystate_sol = solve(steadystate_prob)

# Dynamic Solution
ode_prob = ODEProblem(model_simplified, [], (0, 1))
ode_sol = solve(ode_prob)

# Plotting
plot(ode_sol, title = "Disinfection CSTR")
hline!(steadystate_sol, label="Steady States", linestyle=:dash)

Step-by-step explanations

Before starting the modeling, several packages are imported:

using BioChemicalTreatment # This Package :)
using ModelingToolkit # Basic Modeling Framework. Takes e.g. care of simplification of the equations
using DifferentialEquations # Solve the equations
using Plots # Plotting

Creating the Required Systems

The first step to build the model is to create the individual systems needed for it. But even before that, the chosen process model is set as default:

@set_process process = OzoneDisinfection()

Process OzoneDisinfection 'process':
States (2): see states(process)
  S_O3(t) [guess is 0.0]: S_O3
  X_B(t) [guess is 0.0]: X_B
Parameters (2): see parameters(process)
  kO3 [defaults to 10]
  kd [defaults to 1500]

Then, the next step is to create the reactor. This example uses the CSTR, which needs the volume (1000m³ for this example) as input. The process does not need to be provided, as the default has been set. So writing it out:

@named reactor = CSTR(1000)

CSTR with OzoneDisinfection 'reactor':
States (2): see states(reactor)
  S_O3(t) [defaults to 0.0]: S_O3
  X_B(t) [defaults to 0.0]: X_B
Inflows (1): see inflows(reactor)
  :reactor_reactor₊inflow
Outflows (1): see outflows(reactor)
  :reactor_reactor₊outflow
Parameters (3): see parameters(reactor)
  reactor_process₊kO3 [defaults to 10]
  reactor_process₊kd [defaults to 1500]
  reactor_reactor₊V [defaults to 1000]

Here, and for all following systems, the @named macro is used. This macro takes the name of the assigned variable (here reactor) and passes it as keyword argument called name to the called function, i.e. this line of code is equivalent to

reactor = CSTR(1000;name = :reactor)

Now the reactor is ready, but the influent is still missing. For this the Influent system can be used, which takes three values as input (one of which is the flowrate keyword argument):

1. Values for the states: Either a numeric constant value or a function with one time argument (for non-constant values). In ththe following example, constant values 0.5 and 1 are assumed, please note that the order matters!
2. The flow rate keyword argument: Equivalent to the ones for the values for the compounds, just for the flowrate. Here 500 cubic meter per day is assumed.

And making it yields:

@named influent = Influent([0.5, 1]; flowrate = 500)

Influent 'influent':
Outflows (1): see outflows(influent)
  :influent_to_outflowport₊outflow
Parameters (3): see parameters(influent)
  q_const₊k [defaults to 500.0]: Constant output value of block
  S_O3_const₊k [defaults to 0.5]: Constant output value of block
  X_B_const₊k [defaults to 1.0]: Constant output value of block

And with this, all systems are created and the connection of the systems can be started.

Connecting Systems and Building an Overall Model

For connecting the two systems to an overall model, it only needs to be specified that the influent feeds into the reactor:

More explicitly, the outflow of the influent (yes, it is flowing out of the influent) to the inflow of the reactor. As generally all systems can have multiple inflows and outflows, the accessor functions inflows and outflows by default return a vector, however, as the systems here have only each one port, always the first (and only) element of this vector is taken.

The so created equation is then directly added to a vector of equations:

connect_equations = [connect(outflows(influent)[1], inflows(reactor)[1])]

\[ \begin{equation} \left[ \begin{array}{c} \mathrm{connect}\left( influent_{+}influent_{to\_outflowport_{+}outflow}, reactor_{+}reactor_{reactor_{+}inflow} \right) \\ \end{array} \right] \end{equation} \]

Now the individual systems are there and the connecting equations. Thus the overall model can be built. For this the ODESystem function is used , which gets the equations, the time variable (exported from BioChemicalTreatment as t) and a list of all involved systems in the systems keyword argument:

@named model = ODESystem(connect_equations, t; systems=[influent, reactor])

\[ \begin{equation} \left[ \begin{array}{c} \mathrm{connect}\left( influent_{+}influent_{to\_outflowport_{+}outflow}, reactor_{+}reactor_{reactor_{+}inflow} \right) \\ \mathrm{connect}\left( q_{const_{+}output}, influent_{to\_outflowport_{+}q} \right) \\ \mathrm{connect}\left( S_{O3\_const_{+}output}, influent_{to\_outflowport_{+}S\_O3} \right) \\ \mathrm{connect}\left( X_{B\_const_{+}output}, influent_{to\_outflowport_{+}X\_B} \right) \\ \mathtt{influent.influent\_to\_outflowport.q.u}\left( t \right) = \mathtt{influent.influent\_to\_outflowport.outflow.q}\left( t \right) \\ \mathtt{influent.influent\_to\_outflowport.S\_O3.u}\left( t \right) = \mathtt{influent.influent\_to\_outflowport.outflow.S\_O3}\left( t \right) \\ \mathtt{influent.influent\_to\_outflowport.X\_B.u}\left( t \right) = \mathtt{influent.influent\_to\_outflowport.outflow.X\_B}\left( t \right) \\ \mathtt{influent.q\_const.output.u}\left( t \right) = \mathtt{influent.q\_const.k} \\ \mathtt{influent.S\_O3\_const.output.u}\left( t \right) = \mathtt{influent.S\_O3\_const.k} \\ \mathtt{influent.X\_B\_const.output.u}\left( t \right) = \mathtt{influent.X\_B\_const.k} \\ \mathrm{connect}\left( reactor_{process_{+}rates}, reactor_{reactor_{+}rates} \right) \\ \mathrm{connect}\left( reactor_{reactor_{+}states}, reactor_{process_{+}states} \right) \\ \mathtt{reactor.reactor\_process.rates.S\_O3}\left( t \right) = - \mathtt{reactor.reactor\_process.kO3} \mathtt{reactor.reactor\_process.states.S\_O3}\left( t \right) \\ \mathtt{reactor.reactor\_process.rates.X\_B}\left( t \right) = - \mathtt{reactor.reactor\_process.kd} \mathtt{reactor.reactor\_process.states.X\_B}\left( t \right) \mathtt{reactor.reactor\_process.states.S\_O3}\left( t \right) \\ \mathtt{reactor.reactor\_reactor.inflow.q}\left( t \right) = \mathtt{reactor.reactor\_reactor.outflow.q}\left( t \right) \\ \frac{\mathrm{d} \mathtt{reactor.reactor\_reactor.states.S\_O3}\left( t \right)}{\mathrm{d}t} = \frac{ - \mathtt{reactor.reactor\_reactor.outflow.q}\left( t \right) \mathtt{reactor.reactor\_reactor.states.S\_O3}\left( t \right)}{\mathtt{reactor.reactor\_reactor.V}} + \frac{\mathtt{reactor.reactor\_reactor.inflow.q}\left( t \right) \mathtt{reactor.reactor\_reactor.inflow.S\_O3}\left( t \right)}{\mathtt{reactor.reactor\_reactor.V}} + \mathtt{reactor.reactor\_reactor.rates.S\_O3}\left( t \right) \\ \mathtt{reactor.reactor\_reactor.outflow.S\_O3}\left( t \right) = \mathtt{reactor.reactor\_reactor.states.S\_O3}\left( t \right) \\ \frac{\mathrm{d} \mathtt{reactor.reactor\_reactor.states.X\_B}\left( t \right)}{\mathrm{d}t} = \frac{ - \mathtt{reactor.reactor\_reactor.states.X\_B}\left( t \right) \mathtt{reactor.reactor\_reactor.outflow.q}\left( t \right)}{\mathtt{reactor.reactor\_reactor.V}} + \frac{\mathtt{reactor.reactor\_reactor.inflow.q}\left( t \right) \mathtt{reactor.reactor\_reactor.inflow.X\_B}\left( t \right)}{\mathtt{reactor.reactor\_reactor.V}} + \mathtt{reactor.reactor\_reactor.rates.X\_B}\left( t \right) \\ \mathtt{reactor.reactor\_reactor.outflow.X\_B}\left( t \right) = \mathtt{reactor.reactor\_reactor.states.X\_B}\left( t \right) \\ \end{array} \right] \end{equation} \]

Now the model is built!

But before simulating with the model, it is recommended to first simplify for lower simulation times and easier computations:

model_simplified = structural_simplify(model)

\[ \begin{align} \frac{\mathrm{d} \mathtt{reactor.reactor\_reactor.states.S\_O3}\left( t \right)}{\mathrm{d}t} &= \frac{ - \mathtt{reactor.reactor\_reactor.outflow.q}\left( t \right) \mathtt{reactor.reactor\_reactor.states.S\_O3}\left( t \right)}{\mathtt{reactor.reactor\_reactor.V}} + \frac{\mathtt{reactor.reactor\_reactor.inflow.q}\left( t \right) \mathtt{reactor.reactor\_reactor.inflow.S\_O3}\left( t \right)}{\mathtt{reactor.reactor\_reactor.V}} + \mathtt{reactor.reactor\_reactor.rates.S\_O3}\left( t \right) \\ \frac{\mathrm{d} \mathtt{reactor.reactor\_reactor.states.X\_B}\left( t \right)}{\mathrm{d}t} &= \frac{ - \mathtt{reactor.reactor\_reactor.states.X\_B}\left( t \right) \mathtt{reactor.reactor\_reactor.outflow.q}\left( t \right)}{\mathtt{reactor.reactor\_reactor.V}} + \frac{\mathtt{reactor.reactor\_reactor.inflow.q}\left( t \right) \mathtt{reactor.reactor\_reactor.inflow.X\_B}\left( t \right)}{\mathtt{reactor.reactor\_reactor.V}} + \mathtt{reactor.reactor\_reactor.rates.X\_B}\left( t \right) \end{align} \]

julia> length(equations(model))
19

julia> length(equations(model_simplified))
2

Then, the simulation and plotting step can be started.

Simulating the Model and Plotting

To simulate the model, first a problem has to be built, which is then solved using the solve function on it. For the problem, various possibilities exist, two of them are of them with a particular interest for this package and thus this explanation is restricted to those:

• SteadyStateProblem: This problem type computes the steady state of the model
• ODEProblem: This problem type solves a dynamic simulation of the model for a given time span

Starting with the SteadyStateProblem, this function takes the model and the initial state as input. As the initial state is already set by default, this is left an empty vector in this example but this could be used for overwriting the set initial condition without needing to rebuild the model. See Batch Reactor and Initial Condition for an example of setting new initial conditions.

steadystate_prob = SteadyStateProblem(model_simplified, [])

SteadyStateProblem with uType Vector{Float64}. In-place: true
u0: 2-element Vector{Float64}:
 0.0
 0.0

This problem is solved to get the steady state:

steadystate_sol = solve(steadystate_prob)

retcode: Success
u: 2-element Vector{Float64}:
 0.023809523809523808
 0.013806706114398382

It works similarly for the ODEProblem, however this function takes additionally an argument for the time span to integrate the model. For this example, one day is chosen (0, 1):

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

ODEProblem with uType Vector{Float64} and tType Int64. In-place: true
timespan: (0, 1)
u0: 2-element Vector{Float64}:
 0.0
 0.0

The solution works as for the SteadyStateProblem:

ode_sol = solve(ode_prob)

retcode: Success
Interpolation: 3rd order Hermite
t: 23-element Vector{Float64}:
 0.0
 9.999999999999999e-5
 0.0010999999999999998
 0.00673559486483973
 0.017860228919961384
 0.034015522289822075
 0.05860885522576959
 0.07964720700035285
 0.10562328761738204
 0.1323554874799879
 ⋮
 0.3421176067344822
 0.4010795230402114
 0.4688154772190153
 0.5494805988533654
 0.6448428526683345
 0.7545827590183007
 0.8704257910989464
 0.9813336757297095
 1.0
u: 23-element Vector{Vector{Float64}}:
 [0.0, 0.0]
 [2.4986879592544387e-5, 4.999868754743475e-5]
 [0.000273417971666966, 0.0005497659853782972]
 [0.0016257321982310199, 0.0033436315935539878]
 [0.004071374876262846, 0.008567256406597243]
 [0.007150989121001792, 0.01488247291543346]
 [0.010942178652104926, 0.02075810617821108]
 [0.013492554579288276, 0.02257082884647896]
 [0.01595538018390365, 0.02226410737601743]
 [0.01787757255268276, 0.0207073512896349]
 ⋮
 [0.023153879460625304, 0.01437418470715325]
 [0.023456500242681838, 0.014102235239505113]
 [0.023636166844535723, 0.013949167411311162]
 [0.023735192140278777, 0.013867334321642864]
 [0.02378219788784688, 0.013829166807753724]
 [0.023800875499382475, 0.013814592970047373]
 [0.023806954166967997, 0.013811696705814596]
 [0.02380872033733249, 0.013814351332833317]
 [0.02380886334317414, 0.013810797855815323]

Now the solutions can be plotted. Here the two states of the reactor and the corresponding steady states are plotted:

plot(ode_sol, title = "Disinfection CSTR")
hline!(steadystate_sol, label="Steady States", linestyle=:dash)

And as expected, the continuous solution tends to the expected steady state, and also numerically they match:

julia> isapprox(steadystate_sol.u, ode_sol.u[end]; atol=1e-5)
true

The comparison was done only approximately because, 1. never compare floating point numbers for equality and 2. the dynamic solution might not have reached full steady state yet.

steadystate_sol[rates(reactor).S_O3]

-0.23809523809523808

Checking the model connections (displaying the connected model)

When building a large system, errors in the connections between different tanks occur very quickly. Especially similar parts in the reactor configuration, which make it tempting to copy-paste the connections are prone to errors.

Such errors often lead to invalid equations or unrealistic results and are discovered e.g. due to errors during simplification or unrealistic (or even exploding) simulation results. But after noticing this, finding the error by checking the specified connections is quite cumbersome. To facilitate this, it is possible to display the connected system in a graphical way.

For this, a process diagram can be constructed and displayed:

using TikzPictures # Needed for rendering the diagram
pd = ProcessDiagram([influent, reactor], connect_equations)

set_background!(pd, "white")

See the function documentation for all possible arguments for styling the output. Further, the ordering of the systems is based on a simple heuristic and not guaranteed to provide a good result. However, if this needs to be adapted this is possible using the corresponding functions.

Check out the section of the BSM1 example for using this on a more complex system.

Getting the Equations

As this package allows for full transparency, all equations can be obtained.

For example for ozonation process, the equations are obtained as:

equations(process)

\[ \begin{align} \mathtt{rates.S\_O3}\left( t \right) &= - \mathtt{kO3} \mathtt{states.S\_O3}\left( t \right) \\ \mathtt{rates.X\_B}\left( t \right) &= - \mathtt{kd} \mathtt{states.X\_B}\left( t \right) \mathtt{states.S\_O3}\left( t \right) \end{align} \]

Similarly, this works for the connected model and for the simplified version of it. Note that for the simplified model full_equations has to be used instead, as otherwise not all symbols might be plugged in:

full_equations(model_simplified)

\[ \begin{align} \frac{\mathrm{d} \mathtt{reactor.reactor\_reactor.states.S\_O3}\left( t \right)}{\mathrm{d}t} &= \frac{ - \mathtt{influent.q\_const.k} \mathtt{reactor.reactor\_reactor.states.S\_O3}\left( t \right)}{\mathtt{reactor.reactor\_reactor.V}} + \frac{\mathtt{influent.S\_O3\_const.k} \mathtt{influent.q\_const.k}}{\mathtt{reactor.reactor\_reactor.V}} - \mathtt{reactor.reactor\_process.kO3} \mathtt{reactor.reactor\_reactor.states.S\_O3}\left( t \right) \\ \frac{\mathrm{d} \mathtt{reactor.reactor\_reactor.states.X\_B}\left( t \right)}{\mathrm{d}t} &= \frac{ - \mathtt{influent.q\_const.k} \mathtt{reactor.reactor\_reactor.states.X\_B}\left( t \right)}{\mathtt{reactor.reactor\_reactor.V}} + \frac{\mathtt{influent.X\_B\_const.k} \mathtt{influent.q\_const.k}}{\mathtt{reactor.reactor\_reactor.V}} - \mathtt{reactor.reactor\_process.kd} \mathtt{reactor.reactor\_reactor.states.X\_B}\left( t \right) \mathtt{reactor.reactor\_reactor.states.S\_O3}\left( t \right) \end{align} \]

Concluding Remarks

This example showed the basic principles of building a model on this very simple system. However the same principles apply when building more complex models. See the examples section for the ProcessSimulator to find more complex examples. Notable examples there include, but are not limited to:

• Disinfection with Ozone: The same example as in this usage extended by a second reactor
• Activated Sludge System with ASM1: A detailed example on how to include a external recirculation
• Reactor with multiple Processes: ASM1 with Aeration: How to build a reactor with two processes in it
• BSM1: A more complex and realistic plant based on the benchmark simulation model (BSM) number 1