Usage of the Low-Level Interface

This section gives an introduction to the low-level interface of the package, by explaining the basic principles and specialities of the framework. These principles are then applied directly in a small example.

The low-level interface is more verbose than the high-level interface and thus requires more code and more effort to write. In turn, however, it allows for some more flexibility, e.g. transformation of process rates in a reactor.

Thus, using low-level interface is intended for advanced users with special requirements. But this introduction can also be helpful for other users, to better understand the inner workings of the overall framework and might thus help to debug.

The first section gives a general overview of the setup of this package and the caveats when using it, afterwards a small Low-Level: Example is provided on which the use of the package is introduced step by step.

Low-Level: 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 of the types of systems is provided, the very important difference between reactors and processes is highlighted, the different types of connectors are explained, and the process of setting up a model is detailed in the section Low-Level: Example.

Low-Level: Types of Systems

This package contains four different types of systems, each with different properties and purposes:

• Reactor: Defines the differential equations for the mass balance in the reactor. It has (multiple) inflows and (multiple) outflows and additionally gets the reaction rates from a process and provides the states as output.
• Process: Computes the reaction rates to provide them to a reactor. Takes the states and possibly other external signals as input and returns the reaction rates for the corresponding process as output.
• 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. This element is without dynamics, i.e. it contains no differential equations.
• Standard Library Compatibility: Systems required for converting the signals from this library to ports that can be used in the ModelingToolkitStandardLibrary. The models from there can be used, e.g. for including controllers and other signal processing.

Note that three of them are also part of the high-level interfiace, with the Standard Library Compatibility-type systems added here. These system can of course also be used in the high-level interface (as the two interfaces can be mixed), but there their use is generally not needed, as individual states can be extracted using the Sensor and the additional systems mentioned here are intended for advanced use.

All of these systems can be used together with the ones of the ModelingToolkitStandardLibrary to build models of complex systems.

Low-Level: Reactors vs. Processes

It may be confusing that in this library a distinction between Reactors and Processes is made, since usually the processes take place within reactors and thus the two are inherently combined. Here, this design decision is explained.

Looking at the mathematical model of a reactor, a reactor is generally described by the following mass balance equation for the concentration of every compound C:

\[\frac{dC}{dt} = \frac{q_{in}}{V}C_{in} - \frac{q_{out}}{V}C_{out} + \frac{r_C}{V}\]

where $C_{in}$ is the concentration of C of the inflow to the reactor, $C_{out}$ is the concentration in the outflow and $r_C$ is the rate term describing the process rate. Typically, the outflow of the reactor is considered equal to the current state (the so-called "lumped parameter assumption"):

\[C_{out} = C\]

Note that this assumption is not necessary (and this framework is not restricted to it), but it is made in most cases.

If we now consider a reactor of a certain type, the mass flow (relationship of $q_{in}$, $q_{out}$, and $C_{out}$) are always equal. The processes within the reactor influence only the rate term $r_C$. Therefore, the mass balance is a property of the reactor, while the process rate is determined by the processes.

To make the reactor independent of the process and with it usable with any (combination of) process(es), this framework moved the computation of the reaction rates into an individual system, the Process.

This has several advantages:

1. Each Process and Reactor can be defined independent of each other.
2. Any combination of Processand Reactor is possible.
3. The combination of multiple reactions is achieved by adding the process rates before feeding them to the reactor.

The high-level interface presented in the Usage links Processes and the Reactors automatically. Here, in the low-level, however, the process and the reactor have to be defined individually. For example, for a Completely Stirred Tank Reactor (CSTR) with the Activated Sludge Model No. 1 this definition works as follows:

@named process = ASM1() # defines the process
@named reactor = CSTR(1000, unknowns(states(process))) # The CSTR requires an argument that defines the reactor volume and the states to generate equations for, the second can always be obtained using unknowns(states(process)) with the corresponding process.

eqs = [
    connect(states(reactor), states(process)),
    connect(rates(process), rates(reactor))
] # Equations to connect the states and rates of the reactor and corresponding process.

Low-Level: Connectors

All the systems in this package have connectors. A connector can be thought of as a port with a vector of variables to connect with. For example, a reactor 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 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).
• State Ports: Carry the information about the states of the system. Usually systems have only one port for states used either as input or output depending on the system.
• Rate Ports: Carry the process rates. Usually systems have one (either input or output) port for rates.
• Exogenous Inputs: A system can have a set of exogenous_inputs. These are additional inputs to the system which are not in any of the above categories. An example would be 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 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), one has two possibilities:

1. Use the Sensor to extract a single or multiple signals into one or multiple RealOutputs
2. Unpack all signals into RealOutputs using InflowPortToRealOutput, ReactionInputToRealOutput or StateInputToRealOutput

Similarly, a series of RealInputs can be combined to one of the above ports using RealInputToOutflowPort, RealInputToReactionOutput or RealInputToStateOutput.

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.

Low-Level: 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_O3 and 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",

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

Low-Level: Creating the Required Systems

The first step to build the model is to create the individual systems needed for it. The chosen process model is created as follows:

@named process = OzoneDisinfection()
nameof(process)

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

process = OzoneDisinfection(; name = :process)

This macro is just used for convenience and it will be the same in the following.

The next step is to create the reactor. This example uses the CSTR. From the docs it can be seen that this needs two arguments (besides the name):

1. The reactor volume: For this example 1000 cubic meters is assumed.
2. The states which are in the flows: This depends on the process. Therefore, it can be extracted from the process using unknowns(states(process)), which works independent of the process. This is also the reason why the building of the model started with the process. A variety of systems need this argument, however, for all of them it can be obtained in the same way.

So writing it out:

@named reactor = CSTR(1000, unknowns(states(process)))
nameof(reactor)

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

1. The states in the inflow: As for the reactor, it is taken from the process as unknowns(states(process)).
2. Values for these states: 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, please note that the order matters!
3. 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([0.5, 1], unknowns(states(process)); flowrate = 500)
nameof(influent)

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

Low-Level: Connecting Systems and Building an Overall Model

For connecting the single systems to an overall model, there are two types of connection to be considered: Firstly, the connection of the flow that goes from the influent to the reactor. Secondly, the reactor and the process need to be connected such that the process is taking place in the reactor. Let's start with the second type:

For this, on the one hand, the computed reaction rates need to be connected from the process to the reactor. On the other hand, the process needs the states from the reactor to compute the rates, thus they need to be connected as well. For connecting, the connect function is used, and all connections are collected in a vector.

connect_equations = [
    connect(states(reactor), states(process)),
    connect(rates(process), rates(reactor))
]

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

Now, only the flow from the influent to the reactor remains to be added. Or 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 the vector of equations defined above using the push! function:

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

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

All required connections as well as the individual systems are now defined. Thus, 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, process])

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)

julia> length(equations(model))
19

julia> length(equations(model_simplified))
2

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

Low-Level: Simulating the Model and Plotting

o 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, [])

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

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.00027341797166696594]
 [0.0033436315935539878, 0.0016257321982310197]
 [0.008567256406597243, 0.004071374876262847]
 [0.014882472915433464, 0.007150989121001793]
 [0.020758106178211062, 0.01094217865210492]
 [0.022570828846478947, 0.013492554579288272]
 [0.022264107376017427, 0.01595538018390365]
 [0.020707351289634888, 0.01787757255268276]
 ⋮
 [0.014374184707153242, 0.023153879460625304]
 [0.014102235239505103, 0.023456500242681838]
 [0.013949167411311175, 0.023636166844535723]
 [0.013867334321642855, 0.023735192140278777]
 [0.013829166807753741, 0.02378219788784688]
 [0.013814592970047439, 0.023800875499382475]
 [0.013811696705814775, 0.023806954166967997]
 [0.013814351332833682, 0.0238087203373325]
 [0.013810797855815505, 0.02380886334317415]

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)

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(process).S_O3]

-0.23809523809523808