GitLab

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 support for the modelling of 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. Moreover, this allows 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 the available 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 this page's 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. It 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 directing 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 vector of variables. 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 port which carries the same state variables, but 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. These are additional inputs to the system, which are supplied from other packages. This means they are not flows, but usually connect to the Modeling Toolkit Standard Library. Often, they can be considered control inputs. An example would be the control input for the pumped flow rate of a FlowPump.
  • Exogenous Outputs: A system can have a set of exogenous_outputs. These are additional outputs of the system which are for connections to other packages. This means they are not flows, but usually connect to the Modeling Toolkit Standard Library. Often, they can be considered some form of measurement outputs. An example would be the measurement of a Sensor.
Exogenous Inputs and Outputs

Note that here, exogenous is not meant that it is exogenous to the modeled system. Instead, exogenous is meant in the sense that they are connecting not to systems in library, but instead leading to other julia packages, such as the Modeling Toolkit Standard Library. This is for including functionality from other packages.

For example, the Sensor has the sensed value as exogenous output which is a RealOutput. Similarly, the Aeration model has an exogenous input (a RealInput) for the aeration strength. With this, a sensor for oxygen could for example be used as input for the PID of the ModelingToolkitStandardLibrary, the output of which drives the aeration. This is how an aeration controller would be implemented.

Thus exogenous here means "exogenous to this julia library" and not "exogenous to the model".

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 in order to directly work with the useful systems provided by the blocks module of the ModelingToolkitStandardLibrary such as PID-controllers or mathematical transformations.

Example

The following example shows how to connect 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 S_O3and the bacteria concentration X_B).

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:

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

When executing the examples in this project, be careful when comparing the numeric output. They may, and most likely will, differ in the last few digits.

This is tied to the representation of floating-point numbers and their operations. If interested, you can check out this wikipedia article as a starting point. Otherwise, just compare the first few digits of numeric output.

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([:S_O3 => 0.5, :X_B => 1]; flowrate = 500)

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

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

# Steady 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)
Example block output

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

\[ \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} \]

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)

\[ \begin{equation} \left[ \begin{array}{c} \mathrm{connect}\left( process_{+}rates, reactor_{+}rates \right) \\ \mathrm{connect}\left( reactor_{+}states, process_{+}states \right) \\ \end{array} \right] \end{equation} \]

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)
Note

Each reactor can also be given a keyword argument initial_states, which is used to set the initial conditions of the CSTR. If not set, as here, initial values for all states are set to zero.

Now the reactor is ready, but the influent is still missing. For this the Influent system can be used, which takes two types of input:

  1. Values for the states: A dictionary or vector of pairs from the symbols to the values of the influent. For each variable it is either a numeric constant value or a function with one time argument (for non-constant values). In the following example, constant values 0.5 and 1 are assumed. Instead of a dictionary, a vector could be passed as well, but then 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.

This yields:

@named influent = Influent([:S_O3 => 0.5, :X_B => 1]; flowrate = 500)

\[ \begin{equation} \left[ \begin{array}{c} \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) \\ \end{array} \right] \end{equation} \]

Note

Instead of specifying the flow rate separately, it could as well be included in the list of variables (the first argument of the list, or with key :q) and the value in the corresponding position of the second argument. In this case, do not use the flowrate keyword.

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_{+}inflow \right) \\ \end{array} \right] \end{equation} \]

The individual systems and the connecting equations are now defined. Therefore, the overall model can be built. For this, the System 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 = System(connect_equations, t; systems = [influent, reactor])

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 = mtkcompile(model)

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

On the Simplification

The number of equations in the original system can be obtained using

julia> length(equations(model))
19

Wow! This system has 19 equations! Way more than needed compared to building this system by hand, where usually only two states (one for each compound in the reactor) and associated equations are considered. So why does this system have so many more equations?

This comes from the inner workings of this package. Looking directly at the equations (see Getting the Equations below on how to display them), it becomes apparent that there are:

  • 3 equations for specifying the functions of each influent component (flowrate + 2 states) in the internally used ModelingToolkitStandardLibrary.Blocks.Constant system to specify the equations.
  • 3 connecting those to a RealInputToOutflowPort
  • 3 to combine those to a FlowPort within the RealInputToOutflowPort
  • 1 connecting the influent to the reactor
  • 3 to compute the mass balances in the reactor (flowrate + 2 states)
  • 2 specifying that the outflow of the CSTR is equal to the states
  • 2 to compute the reaction rates in the process
  • 1 connecting the states of the reactor to the process
  • 1 connecting the reaction rates from the process to the reactor

Most of these equations have not been specified explicitly, but rather implicitly by this framework in order to allow for more flexibility. Details on the inner workings are provided in the section Low Level Interface.

Investigating all these equations (and also from their purpose shortly described here) it can be seen that most of them are simple assignments like, for instance, all connection equations. The simplification process then simplifies the system by various transformations, amongst others, by plugging in all those assignments directly into other equations. This way, a simpler description of the system is found, which is easier to understand and more efficient to compute. Indeed, looking at the number of equations of the simplified system there are only two left, as are required when building by hand:

julia> length(equations(model_simplified))
2

Now we are ready to start the simulation and plotting step.

Simulating the Model and Plotting

To simulate the model, a problem has to be built in a first step, which is then solved using the solve function on it. Various possibilities exist to define the problem. Two of them are of particular interest for this package:

  • 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

The SteadyStateProblem takes the model and the initial state as input. As the initial state is already set by default, an empty vector can be passed in this example. However, this argument could be used for overwriting the set initial condition without rebuilding 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

We solve the problem to get the steady state:

steadystate_sol = solve(steadystate_prob)
retcode: Success
u: 2-element Vector{Float64}:
 0.013806706114398382
 0.023809523809523808

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
Initialization status: FULLY_DETERMINED
Non-trivial mass matrix: false
timespan: (0, 1)
u0: 2-element Vector{Float64}:
 0.0
 0.0

Solving the problem is similar to 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]
 [4.999868754743475e-5, 2.4986879592544387e-5]
 [0.0005497659853782972, 0.000273417971666966]
 [0.0033436315935539878, 0.0016257321982310199]
 [0.008567256406597243, 0.004071374876262846]
 [0.014882472915433462, 0.007150989121001792]
 [0.020758106178211083, 0.010942178652104926]
 [0.022570828846478957, 0.013492554579288276]
 [0.022264107376017417, 0.01595538018390365]
 [0.020707351289634884, 0.01787757255268276]
 ⋮
 [0.01437418470715325, 0.023153879460625304]
 [0.014102235239505101, 0.023456500242681838]
 [0.013949167411311181, 0.023636166844535723]
 [0.01386733432164287, 0.023735192140278777]
 [0.013829166807753743, 0.02378219788784688]
 [0.013814592970047432, 0.023800875499382475]
 [0.013811696705814735, 0.023806954166967997]
 [0.013814351332833576, 0.02380872033733249]
 [0.013810797855815453, 0.02380886334317414]

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

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

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.

Values for variables removed during simplification

As explained in the note above, the system has been simplified to contain only the two states in the reactor. This is a formidable performance boost (solving 2 instead of 19 equations)! However, one might be interested to get the solution for signals that were discarded in the simplification process. Luckily, there is a way to recover all of those from the solution. For example for the rate of S_O3 this can be done as follows:

steadystate_sol[rates(reactor).S_O3]
-0.23809523809523808

Similarly it works for all other variables and as well for the dynamic simulation.

Checking the model connections (displaying the connected model)

When building a large system, errors in the connections between different tanks can happen easily. 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)
Example block output
Warning

If you see the diagram only barely, this is most likely because the background of the generated image is transparent by default. To change this add e.g. a white background using

set_background!(pd, "white")

or directly add $background="white"$ when generating the process diagram.

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.

Displaying the diagram

Displaying the diagram requires TikzPictures, which has it's own dependencies. Check them out directly there.

Alternatively, you can as well generate a LaTeX document to compile externally using get_latex. If you cannot compile LaTeX on your system, you can as well use an online service for this.

For example, under https://www.quicklatex.com you can find such a service. To use it, paste the document preamble (the output of get_latex_preamble) into the Custom LaTeX Document Preamble section on this webpage (you need to expand "Choose Options" first) and then the tikz picture (output of get_tikzpicture) into the field for the latex code (titled Type LaTeX Code on the webpage).

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.states.X\_B}\left( t \right)}{\mathrm{d}t} &= \frac{ - \mathtt{influent.q\_const.k} \mathtt{reactor.reactor.states.X\_B}\left( t \right)}{\mathtt{reactor.reactor.V}} + \frac{\mathtt{influent.X\_B\_const.k} \mathtt{influent.q\_const.k}}{\mathtt{reactor.reactor.V}} - \mathtt{reactor.process.kd} \mathtt{reactor.reactor.states.X\_B}\left( t \right) \mathtt{reactor.reactor.states.S\_O3}\left( t \right) \\ \frac{\mathrm{d} \mathtt{reactor.reactor.states.S\_O3}\left( t \right)}{\mathrm{d}t} &= \frac{\mathtt{influent.S\_O3\_const.k} \mathtt{influent.q\_const.k}}{\mathtt{reactor.reactor.V}} + \frac{ - \mathtt{influent.q\_const.k} \mathtt{reactor.reactor.states.S\_O3}\left( t \right)}{\mathtt{reactor.reactor.V}} - \mathtt{reactor.process.kO3} \mathtt{reactor.reactor.states.S\_O3}\left( t \right) \end{align} \]

Concluding Remarks

This example showed the basic principles of building a model for a 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 include, but are not limited to: