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

Low Level: Types of Systems

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

• Reactor: Defines the differential equations for the mass balance in the reactor. 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 direct 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 now 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.

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 downside to this is, that the building of models is slightly more complicated as the process and the reactor have to be defined individually, but this should be feasible. For example, for a Completely Stirred Tank Reactor (CSTR) with the Activated Sludge Model Nr. 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.

In the future, we plan to develop a high-level interface which does this automatically.

Low Level: 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 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 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 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, which 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, 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), 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, 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.

Low Level: 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:

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.

Then, 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 ththe 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.

And making it 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: First 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 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} \]

Having this, 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} \]

Now all required connections are there as well as the individual systems. 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, 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 = structural_simplify(model)

julia> length(equations(model))
19

julia> length(equations(model_simplified))
2

Then, the simulation and plotting step can be started.

Low Level: 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, [])

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

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.004071374876262847, 0.008567256406597243]
 [0.007150989121001792, 0.014882472915433466]
 [0.010942178652104926, 0.020758106178211083]
 [0.013492554579288274, 0.022570828846478944]
 [0.015955380183903647, 0.022264107376017413]
 [0.01787757255268275, 0.020707351289634864]
 ⋮
 [0.023153879460625297, 0.01437418470715324]
 [0.02345650024268183, 0.014102235239505101]
 [0.02363616684453572, 0.013949167411311183]
 [0.023735192140278777, 0.013867334321642869]
 [0.023782197887846877, 0.013829166807753712]
 [0.023800875499382472, 0.01381459297004733]
 [0.023806954166967983, 0.01381169670581443]
 [0.023808720337332497, 0.013814351332832831]
 [0.023808863343174145, 0.013810797855815073]

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

-0.23809523809523808