Adding new Processes

This section explains how to add new processes to the framework.

There are two ways to specify new processes:

• The first one is the most general one and allows (or requires) to specify all equations directly in code. It is discussed in the section Adding a new process using Equations.
• The second one provides a way to specify it using the stoichiometric matrix notation known from the ASM-type models found in (Henze et al., 2006). This notation is also known as Gujer- or Petersen-Matrix. This variant is discussed in the section Adding a new process using Matrices.

About Processes

First, before we dive into how to define new processes, we will have a look into what we consider to be a Process and what possible elements it comprises:

Thus, a process is independent from a reactor or similar in which it might be simulated and does not contain any information about flows and balances of it. Instead, for usage within a simulation, it is intended to be connected to a Reactor via ports for the process rates (the port called rates) and for the current states (the port called states). This has the advantage, that it allows for simple combination of different processes (e.g. an ASM with an aeration model) by adding the process rates before connecting to the reactor.

For it's computations, each Process has the following properties (some are optional):

• A name
• The rates port
• The states port
• equations describing the relation to compute the rates
• t: The independent variable for time
• parameters for the equations (can be empty)
• Internal variables used for the equations, but not part of neither rates nor states (optional)
• Additional inputs for the computations (e.g. Oxygen Flow for the aeration) (optional)

Adding a new process using Equations

The most general way to add new processes works by directly specifying the equations as Julia code.

The advantage of this way is that it gives the full flexibility of Julia, on the down side, it requires you to write the Julia code for it.

But now to this approach. To use it there are four required steps:

1. Creating connectors for the input states and for the output rates
2. Generating symbols for the parameters and eventual internal variables and additional inputs
3. Specifying the process equations
4. Generating a Process struct with all this information

To be able to reuse it and to make the usage in building a system more clear, this should be done in a (properly-named) function which takes the name and values for the parameters as inputs and returns the resulting Process struct.

In the following this process will be shown on the example of the simple aeration model with oxygen transfer rate

\[OTR = k_La*(S_{O2}^{max} - S_{O2})\]

where $S_{O2}^{max}$ is the saturation concentration of oxygen, $S_{O2}$ the current oxygen concentration and $k_La$ the volumetric mass transfer rate.

As the Process should only return the computed process rate for further use (without any integration or similar), the output of this system is $S_{O2, rate} = OTR$.

This system has one parameter ($S_{O2}^{max}$) and one external input ($k_La$). Further it gets one state ($S_{O2}$) and should return one process rate (for $S_{O2}$).

Now in the following chapters you will find a description how to generate such a system.

The function header

First, the function header has to be built. It should take the name of the system and the parameter values. Further, we give it the initial value for the $k_La$ input:

@component function Aeration(; name, k_La_init = 240, S_O2_max = 8)

end 

Here we used keyword arguments for the name, the initial $k_La$ and the parameter $S_{O2}^{max}$ as this allows to provide them by name instead of position. This is convenient especially with many parameters. Further, default values are provided for the initial $k_La$ (240) and the parameter $S_{O2}^{max}$ (8) for convenience for the user. However, the name does not have a default and thus has to be provided. This is like that, as no two systems should have the same name, thus it is better not to set a default here. Further, the $@component$ macro has been used on the function to mark it as component for the underlying ModelingToolkit library.

Creating the connectors

Then, we need to create the connectors to get the current $S_{O2}$ state as input and providing the rate as output. To do so, we use the StateInputPort and ReactionOutputPort functions. First we start with the state and call the StateInputPort function:

@named states = StateInputPort(Dict([
# State name => particle size (whether particulate colloidal or soluble)
  "S_O2" => (particulate = false, colloidal = false, soluble = true),
  ]))

Here we provide a Dict as input, with keys being the name of the state (here S_O2) and values being a NamedTuple with contents describing wheter this state is particulate, colloidal or soluble. All properties have to be provided for every single state. Actually, this function would also need a name but instead of explicitly specifying it, we used the convenient @named macro here, which automatically assignes the variable name (here states) as name.

Then, we continue with the rates output port, where we can take all info directly from the states system:

@named rates = ReactionOutputPort(unknowns(states)) 

Here, we used again the @named macro for the same purpose as above, but the specifying of the variables is much easier: With the call unknowns(states) we get directly all variables from the StateInputPort defined above and the ReactionOutputPort then extracts all needed information directly from them.

Generating symbols for the parameters and creating the connectors for the additional inputs

To generate the symbols for the parameters we can use the @parameters macro:

parameters = @parameters S_O2_max=S_O2_max

Here, a parameter called S_O2_max is created and the current value of the variable S_O2_max is set as its default value. Further the parameter is from now on available under the variable S_O2_max (overwriting the default value, as it is already set and not needed anymore after here, this is not a problem) and additionally the variable parameters now is a vector with all parameters specified in this line.

  parameters = @parameters param1=param1_value param2=param2_value

Then, a port for the input of the mass transfer rate $k_La$ needs to be generated. For this we use the RealInput from the ModelingToolkitStandardLibrary:

@named k_La = ModelingToolkitStandardLibrary.Blocks.RealInput(nin=1, guess=k_La_init)

Here again we need to provide a name (here through the @named macro), the number of inputs specified as nin and the guess specifies the initial value. See the documentation there for detailed info on the arguments.

Specifying the equations

Now, as we have all ports and parameters, we can specify the equations:

eqs = [
    rates.S_O2 ~ k_La.u * (S_O2_max - states.S_O2)
]

Here, rates.S_O2 takes the S_O2 variable of the rates port. As this is the one we want to specify, it is on the left hand side of the equation (the ~ designates the equal sign = and is used as the equal sign is reserved for assignements). Then, the left hand side of the equation is set to the equation: k_La.u is the input from the port (the .u is the default for getting the variable of the port here, see the docs of the ModelingToolkitStandardLibrary), then we have the parameter S_O2_max defined above and the state input from states.S_O2. If multiple equations are needed, simply add the second one to this array.

Generating the Process

Finally, we can assemble the Process struct using everything above:

Process(eqs, t, parameters, rates, states, [k_La]; name=name)

Hereby the arguments are:

• The vector of equations
• The time (independent variable). The library exports this as t.
• The vector of parameters (if none specify as [])
• The rates port
• The states port
• The vector of additional inputs (here only one, but needs to be a vector, thus [k_La])
• The name as keyword argument, here just forwarding the name provided to the function

The resulting function

Putting everything together, we get the following function:

@component function Aeration(; name, k_La_init = 240, S_O2_max = 8)

    # Generate the ports
    @named states = StateInputPort(Dict([
           "S_O2" => (particulate = false, colloidal = false, soluble = true),  # oxygen
           ]))
    @named rates = ReactionOutputPort(unknowns(states))

    # Generate the Parameters and additional inputs
    # parameters
    parameters = @parameters S_O2_max=S_O2_max
    # additional inputs
    @named k_La = ModelingToolkitStandardLibrary.Blocks.RealInput(nin=1, guess=k_La_init)

    # Specifying the Equations
    eqs = [
        rates.S_O2 ~ k_La.u * (S_O2_max - states.S_O2)
    ]

    # Generating the Process and returning it
    Process(eqs, t, parameters, rates, states, [k_La]; name=name)
end 

Now, this function can then be used as any other process that is already defined in the package.

Adding a new process using Matrices

This way allows to define models by specifying the stoichiometric matrix of it, also known as Gujer- or Petersen matrix. This is for example the default way the ASM-type models are provided (Henze et al., 2006). This is convenient for all models that can be represented as such (which is the case for many process models in the water engineering world), but is also restricted to those models.

To do so, 6 csv-files are neccessary to define (in the filenames, replace MODELNAME by the name of the model):

• MODELNAME.md: Short description of the model
• MODELNAME_states.csv: Listing the states and their properties
• MODELNAME_processrates.csv: Listing the process rates, their properties and equations
• MODELNAME_matrix.csv: The stoichiometric matrix
• MODELNAME_compositionmatrix.csv: The composition matrix
• MODELNAME_parameters.csv: The model parameters and their properties

In the section Explanation of the files, these files will be explained in detail.

Once all these files are defined, one can read in the model in different ways:

• By providing the paths to all single files

  Process(states_path, parameters_path, processrates_path, matrix_path, compositionmatrix_path)

• By providing the path of the folder where the files are located and the MODELNAME. In this case the files have to be named as in the list above. In case that some of the files are not in the folder or do not adhere to the naming specification, optional keyword arguments are available to specify the paths for single files. The corresponding names are states_path, parameters_path, processrates_path, matrix_path or compositionmatrix_path.

  Process(folder_path, "MODELNAME")

• If the model is officially added to the package, one can simply specify the model name:

  Process("MODELNAME")

In all cases, additional keyword arguments to set the parameter values (names according to the model specification) can be provided.

Explanation of the files

In the following, the neccessary files and their contents are desribed in detail. For the filenames, which are neccessary to be as menitoned here for the cases where the individual paths are not provided, MODELNAME is to be replaced by the actual name of the model. As an example in the following, the implemented ASM1 is used.

Model Description

Filename: MODELNAME.md

Content: A short description of the model and references for it. This file is not strictly required for executing the model, but still a good idea to add to know what the model is about.

If it is for inclusion in the official package, we require this file. Further, in this case, you can use the following:

• For citing papers or similar, consider adding the citation to the refs.bib in the docs folder. Then the citations referenced as [CITE_KEY](@cite) and to generate the list of references, use:

  ```@bibliography
  Pages = []
  Canonical = false
  
  CITE_KEY1
  CITE_KEY2
  ```

• For automatically inserting a table with the state, parameter or processrate properties or the stoichiometric or composition matrix you can use:

  ```@eval
  Main.state_descriptions[:MODELNAME] # For the state property table
  Main.processrate_descriptions[:MODELNAME] # For the processrate property table
  Main.stoichmat[:MODELNAME] # For the Stoichiometric matrix as table
  Main.compositionmat[:MODELNAME] # For the composition matrix as table
  Main.parameter_descriptions[:MODELNAME] # For the parameter property table
  ```

Example for ASM1:

# Activated Sludge Model nr. 1 (ASM1)
The original source of the activated sludge model is [Henze:1987](@cite). This model has been tested and extended since, so that we decided to further include the information from [Hauduc:2010](@cite) and [Corominas:2010](@cite). The goal was to provide these models in a machine readable form with a mass balance that sums up to zero.

!!! note "Using this Model"
    Create a process with this model using

    ```julia
    Process("ASM1")
    ```
    
    Add keyword arguments to overwrite the default values of the parameters. E.g.
    
    ```julia
    Process("ASM1"; Y_OHO=1)
    ```

    This works for all parameters below and just add multiples for every to overwrite.

!!! warning "Initial State"
    When using this Model, one needs to set the initial state/initial condition **in the corresponding reactor**. This is because the default values are all `0`, and the some of the process equations have a division by `0` if all states are `0`!

## States
```@eval
Main.state_descriptions[:ASM1]
```

## Process Rates
```@eval
Main.processrate_descriptions[:ASM1]
```

## Stoichiometric Matrix
```@eval
Main.stoichmat[:ASM1]
```

## Composition Matrix
```@eval
Main.compositionmat[:ASM1]
```

## Parameters
```@eval
Main.parameter_descriptions[:ASM1]
```

## References
```@bibliography
Pages = []
Canonical = false

Hauduc:2010
Henze:1987
Corominas:2010
```

Model States

Filename: MODELNAME_states.csv

Contents: A csv table with a list of all states and the properties. It has the following columns:

• name is the name of the state that the model will use.
• description is a string that describes the state.
• particle_size describes the size of the particles described by the state. particulate, colloidalor soluble are possible. If multiple apply they are separated by comma (,)

Example for ASM1:

name;description;particle_size
S_B;Soluble biodegradable organics;soluble
S_U;Soluble nondegradable organics;soluble
S_O2;Dissolved oxygen;soluble
XC_B;Particulate and colloidal biodegradable organics;particulate,colloidal
X_UInf;Particulate nonbiodegradable organics from the influent;particulate
X_UE;Particulate nonbiodegradable endogenous products;particulate
S_NHx;Ammonia (NH4 + NH3);soluble
S_NOx;Nitrate and nitrite (NO3 + NO2) (considered to be NO3 only for stoichiometry);soluble
XC_BN;Particulate and colloidal biodegradable organic N;particulate,colloidal
S_BN;Soluble biodegradable organic N;soluble
X_OHO;Ordinary heterotrophic organisms;particulate
X_ANO;Autotrophic nitrifying organisms (NH4+ to NO3-);particulate
S_Alk;Alkalinity (HCO3-);soluble
S_N2;Dissolved nitrogen (gas, N2);soluble

Model Processrates

Filename: MODELNAME_processrates.csv

Contents: A csv table with a list of all processrates, the properties and equations. It has the following columns:

• name is the name of the process that the model will use.
• description is a string that describes the process rates.
• equation is the equation that describes the process rate taken from Hauduc at al. (2010), which originates from Henze et al. (1987)

Example for ASM1:

name;description;equation
g_hO2;Aerobic growth of heterotrophs;m_OHOMax*(S_B/(K_SBOHO+S_B))*(S_O2/(K_O2OHO+S_O2))*(S_NHx/(K_NHxOHO+S_NHx))*X_OHO
g_hAn;Anoxic growth of heterotrophs;m_OHOMax*(S_B/(K_SBOHO+S_B))*(K_O2OHO/(K_O2OHO+S_O2))*(S_NOx/(K_NOxOHO+S_NOx))*(S_NHx/(K_NHxOHO+S_NHx))*n_mOHOAx*X_OHO
g_aO2;Aerobic growth of autotrophs;m_ANOMax*(S_NHx/(K_NHxANO+S_NHx))*(S_O2/(K_O2ANO+S_O2))*X_ANO
d_h;Decay of heterotrophs;b_OHO*X_OHO
d_a;Decay of autotrophs;b_ANO*X_ANO
am_N;Ammonification of soluble organic Nitrogen;q_am*S_BN*X_OHO
ho;Hydrolysis of entrapped organics;q_XCBSBhyd*((XC_B/X_OHO)/(K_XCBhyd+XC_B/X_OHO))*((S_O2/(K_O2OHO+S_O2))+n_qhydAx*(K_O2OHO/(K_O2OHO+S_O2))*(S_NOx/(K_NOxOHO+S_NOx)))*X_OHO
ho_N;Hydrolysis of entrapped organic nitrogen;q_XCBSBhyd*(XC_BN/XC_B)*((XC_B/X_OHO)/(K_XCBhyd+XC_B/X_OHO))*((S_O2/(K_O2OHO+S_O2))+n_qhydAx*(K_O2OHO/(K_O2OHO+S_O2))*(S_NOx/(K_NOxOHO+S_NOx)))*X_OHO

Model Stoichiometric Matrix

Filename: MODELNAME_matrix.csv

Contents: The stoichiometric matrix defining the stoichiometry of the model. Contains the names of all states (has to match the states description above) and the process rate names (also has to match the names in the corresponding file). Further all used parameters have to match the ones specified in the parameters file. Zeros in the stoichiometry can be ommitted.

Example for ASM1:

process\state;S_U;S_B;X_UInf;XC_B;X_OHO;X_ANO;X_UE;S_O2;S_NOx;S_NHx;S_BN;XC_BN;S_Alk;S_N2
g_hO2;;-1/Y_OHO;;;1;;;-(1-Y_OHO)/Y_OHO;;-i_NXBio;;;-i_NXBio*i_ChargeSNHx;
g_hAn;;-1/Y_OHO;;;1;;;;-(1-Y_OHO)/(i_NO3N2*Y_OHO);-i_NXBio;;;-(1-Y_OHO)/(i_NO3N2*Y_OHO)*i_ChargeSNOx-i_NXBio*i_ChargeSNHx;(1-Y_OHO)/(i_NO3N2*Y_OHO)
g_aO2;;;;;;1;;-(-i_CODNO3-Y_ANO)/Y_ANO;1/Y_ANO;-i_NXBio-1/Y_ANO;;;-(i_NXBio+1/Y_ANO)*i_ChargeSNHx+(1/Y_ANO)*i_ChargeSNOx;
d_h;;;;1-f_XUBiolys;-1;;f_XUBiolys;;;;;i_NXBio-f_XUBiolys*i_NXUE;;
d_a;;;;1-f_XUBiolys;;-1;f_XUBiolys;;;;;i_NXBio-f_XUBiolys*i_NXUE;;
am_N;;;;;;;;;;1;-1;;i_ChargeSNHx;
ho;;1;;-1;;;;;;;;;;
ho_N;;;;;;;;;;;1;-1;;

Model Composition Matrix

Filename: MODELNAME_compositionmatrix.csv

Contents: The composition matrix is included to be able to check the models mass balance symbolically. !!!note The composition matrix is currently only used to check the model and not to infer parameters.

Example for ASM1:

composition\state;S_U;S_B;X_UInf;XC_B;X_OHO;X_ANO;X_UE;S_O2;S_NOx;S_NHx;S_BN;XC_BN;S_Alk;S_N2
COD;1;1;1;1;1;1;1;-1;i_CODNO3;0;0;0;0;i_CODN2
N;0;0;0;0;i_NXBio;i_NXBio;i_NXUE;0;1;1;1;1;0;1
Charge;0;0;0;0;0;0;0;0;i_ChargeSNOx;i_ChargeSNHx;0;0;-1;0

Model Parameters

Filename: MODELNAME_parameters.csv

Contents: A csv table with a list of all parameters and their properties. It has the following columns:

• name The name of the parameter that the model will use. We took it from Hauduc et al. (2010) and changed it into a machine readable form.
• description A string that describes the parameter.
• unit_description A description of the unit containing for example if only grams of Nitrogen or COD are included, i.e. what is actually measured.
• value A standard value for the parameter. This might need to be calibrated.
• expression Some parameters are calculated from other parameter. They have an expression. A parameter should either have an expression or a value, not both.
• temperature The temperature for which the typical parameter value is chosen.

Example for ASM1:

name;description;latex;unit;unit_description;value;expression;temperature;type
Y_OHO;Yield for XOHO growth;\gamma_\text{h};g/g;gXOHO/gXCB;0.67;;20;stoichiometric
f_XUBiolys;Fraction of XU generated in biomass decay;\zeta_\text{Xnb,l};g/g;gramXU/gramXBio1;0.08;;20;stoichiometric
Y_ANO;Yield of XANO growth per SNO3;\gamma_\text{a};g/g;gramXAUT/gramSNO3;0.24;;20;stoichiometric
i_NXBio;N content of biomass (XOHO, XPAO, XANO);\iota_\text{NXb};g/g;gramN/gramXBio;0.086;;20;stoichiometric
i_NXUE;N content of products from biomass;\iota_\text{NXnb};g/g;gramN/gramXUE;0.06;;20;stoichiometric
i_NO3N2;Conversion factor for NO3 reduction to N2;\iota_\text{NO3,N2};g/g;gramCOD/gramN;;(COD_N-(COD_N+3*COD_O+COD_neg))/(M_N);20;stoichiometric
i_CODNO3;Conversion factor for NO3 in COD;\iota_\text{COD,NO3};g/g;gramCOD/gramN;;(COD_N+3*COD_O+COD_neg)/M_N;20;stoichiometric
i_CODN2;Conversion factor for N2 in COD;\iota_\text{COD,N2};g/g;gramCOD/gramN;;(2*COD_N)/(2*M_N);20;stoichiometric
i_ChargeSNHx;Conversion factor for NHx in charge;\iota_\text{cSNHx};mol/g;Charge/gramN;;1/M_N;20;stoichiometric
i_ChargeSNOx;Conversion factor for NO3 in charge;\iota_\text{cSNOx};mol/g;Charge/gramN;;-1/M_N;20;stoichiometric
q_XCBSBhyd;Maximum specific hydrolysis rate of particulate and soluble biodegradable organics;\lambda_\text{hyd,b};g/g;gramXCB/gramXOHO/d;3;;20;kinetic
K_XCBhyd;Saturation coefficient for XB/XOHO;\kappa_\text{Xhyd};g/g;gramXCB/gramXOHO;0.03;;20;kinetic
n_qhydAx;Correction factor for hydrolysis under anoxic conditions;\eta_\text{hyd,an};-;-;0.4;;20;kinetic
m_OHOMax;Maximum growth rate of XOHO;\mu_\text{max,Xh};1/day;per day;6;;20;kinetic
n_mOHOAx;Reduction factor for anoxic growth of XOHO;\eta_\text{an,h};-;-;0.8;;20;kinetic
K_SBOHO;Half-saturation coefficient for SB;\kappa_\text{SB};g/m^3;gramSB/m^3;20;;20;kinetic
b_OHO;Decay rate for XOHO;\beta_\text{h};1/day;per day;0.62;;20;kinetic
K_O2OHO;Half-saturation coefficient for SO2 XOHO;\kappa_\text{O2,h};g/m^3;gramSO2/m^3;0.2;;20;kinetic
K_NOxOHO;Half-saturation coefficient for SNOx XOHO;\kappa_\text{NOx,h};g/m^3;gramSNOx/m^3;0.5;;20;kinetic
K_NHxOHO;Half-saturation coefficient  for NH4*;\kappa_\text{NHx,h};g/m^3;gramSNHx/m^3;0.05;;20;kinetic
m_ANOMax;Maximum growth rate of XANO;\mu_\text{max,a};1/day;per day;0.8;;20;kinetic
b_ANO;Decay rate for XANO;\beta_\text{a};1/day;per day;0.15;;20;kinetic
q_am;Rate constant for ammonification;\lambda_\text{am};m^3/g/day;m^3/gramXCB,N/day;0.08;;20;kinetic
K_O2ANO;Half-saturation coefficient for SO2 for XANO;\kappa_\text{O2,a};g/m^3;gramSO2/m^3;0.4;;20;kinetic
K_NHxANO;Half-saturation coefficient for SNHx for XANO;\kappa_\text{NHx,a};g/m^3;gramSNHx/m^3;1;;20;kinetic
COD_neg; Theoretical COD of negative charge; \text{COD}_{-};g/mol;gramCOD/mol;8;;;factor
COD_pos; Theoretical COD of positive charge; \text{COD}_{+};g/mol;gramCOD/mol;-8;;;factor
COD_C; Theoretical COD of molar carbon; \text{COD}_\text{C};g/mol;gramCOD/mol;32;;;factor
COD_N; Theoretical COD of molar nitrogen; \text{COD}_\text{N};g/mol;gramCOD/mol;-24;;;factor
COD_H; Theoretical COD of molar hydrogen; \text{COD}_\text{H};g/mol;gramCOD/mol;8;;;factor
COD_O; Theoretical COD of molar oxygen; \text{COD}_\text{O};g/mol;gramCOD/mol;-16;;;factor
COD_S; Theoretical COD of molar sulphur; \text{COD}_\text{S};g/mol;gramCOD/mol;48;;;factor
COD_P; Theoretical COD of molar phosphorus; \text{COD}_\text{P};g/mol;gramCOD/mol;40;;;factor
COD_Fe; Theoretical COD of molar iron; \text{COD}_\text{Fe};g/mol;gramCOD/mol;24;;;factor
M_N; atomic molar mass of nitrogen; M_\text{N};g/mol;gram/mol;14;;;mass