Process

Click the badge below to try this tutorial interactively in your browser:

You can also run this tutorial inGoogle Colab. It takes a one-time setup per session: follow theColab instructions.

• Prepared by:

  • Joy Zhang

• Learning objectives. After this tutorial, you will be able to:

  • Define a Process with rate expressions and stoichiometry

  • Group processes into a Processes collection

  • Use processes inside a dynamic unit to drive state evolution

• Prerequisites: 6. System

• Covered topics:

  • 1. Understanding Process/Processes/process_models

  • 2. Making a simple biokinetic model

  • 3. Advanced features

Companion video. A walkthrough of this tutorial is available on YouTube (Part 1, Part 2), presented by Joy Zhang. Recorded against QSDsan v1.2.5. The concepts still apply, but if the code on screen differs from this notebook, follow the notebook.

Setup

Import QSDsan and confirm the installed version.

import qsdsan as qs
print(f'This tutorial was made with qsdsan v{qs.__version__}.')

This tutorial was made with qsdsan v1.5.3.

1. Understanding Process/Processes/process_models

This tutorial focuses on the Process module in qsdsan, which provides the modeling capacity for dynamic conversions of waste streams in a SanUnit. With it, simulations can describe not only a system’s mass and energy flows at steady state but also how those flows change over time.

Similar to the Component module, the Process module includes three main classes: Process, Processes, and CompiledProcesses.

The Process class is used to model individual dynamic processes with user-defined stoichiometry and reaction kinetics. The Processes class is simply a collection of multiple Process objects. Compilation of Processes yields CompiledProcesses, which provides additional attributes of the collection. It is often used to implement (pseudo-)mechanistic models with parallel biological and/or physiochemical processes in dynamic simulations of treatment systems.

from IPython.display import Image
Image(url='https://lucid.app/publicSegments/view/2c231fa2-6065-46b9-83af-a790ce38b6c0/image.png', width=600)

In Section 1, we will focus on understanding the structures of these classes by going over some existing processes in qsdsan. Next section will guide you through the creation of your own Processes subclass.

1.1. Processes vs. process_models

You might remember SanUnit vs. unit_operations from the SanUnit tutorial. Processes and process_models have the same kind of distinction: Processes is a class, while process_models is the package that holds the pre-built process models.

# If you check
qs.Process
# or
# qs.Processes
# or
# qs.CompiledProcesses

qsdsan._process.Process

This means Process/Processes/CompiledProcesses is a class in the _process module (i.e., the codes defining these classes are in the _process.py file within the qsdsan directory).

# If you check
qs.process_models

<module 'qsdsan.process_models' from 'C:\\Users\\Yalin\\Documents\\Coding\\QSDsan-platform\\QSDsan\\qsdsan\\process_models\\__init__.py'>

This shows that process_models is actually a package (i.e., a directory) rather than a single class. (Why? See here.) It bundles the pre-built process models shipped with qsdsan.

Note: qsdsan.processes is a legacy alias for qsdsan.process_models (just as qsdsan.sanunits aliases qsdsan.unit_operations). Both names point to the same package, but process_models is the current, preferred name.

For the full inventory of available models and helper functions, see the Process Models API reference. Briefly, what qsdsan.process_models.__all__ exposes falls into three groups:

• CompiledProcesses subclasses such as ASM1, ASM2d, mASM2d, ADM1, ADM1p, and PM2 (parallel multi-process biokinetic models),

• helper functions named create_<modelname>_cmps that build a matching CompiledComponents set, plus rate-function utilities like pH_inhibit and Hill_inhibit, and

• standalone classes like Decay and KineticReaction for simple processes whose state variables can be solved analytically.

Submodule scripts inside the package carry a single leading underscore by convention (e.g., _asm1 for the ASM1 implementation) but can still be accessed directly:

# For example, `_asm1` is the submodule that defines the ASM1 model.
qs.process_models._asm1

<module 'qsdsan.process_models._asm1' from 'C:\\Users\\Yalin\\Documents\\Coding\\QSDsan-platform\\QSDsan\\qsdsan\\process_models\\_asm1.py'>

For the rest of this section, we will use activated sludge model no. 1 (ASM1) as an example to learn about the core classes in the _process module.

ASM1

ASM1 is a mathematical model for the activated sludge processes in wastewater treatment systems. It uses 13 components to characterize wastewater and activated sludge, including biomass, organic substrates (COD), organic and inorganic nitrogens (N), dissolved oxygen (DO), alkalinity etc.

ASM1 includes 8 transformation processes of the components (e.g., hydrolysis, biomass growth and decay). The stoichiometry of each process considers conservations of COD, N, and electric charges. The process rates are dependent on environmental conditions and state variables (i.e., relevant component concentrations).

# Before we get to the subsections, let's get ready by loading ASM1-related objects in qsdsan
from qsdsan.process_models import create_asm1_cmps, ASM1
cmps = create_asm1_cmps()
cmps.show()

CompiledComponents([
    S_I,   S_S,  X_I,  X_S,
    X_BH,  X_BA, X_P,  S_O,
    S_NO,  S_NH, S_ND, X_ND,
    S_ALK, S_N2, H2O,
])

This means the create_asm1_cmps function returns a CompiledComponents object corresponding to the 13 components in ASM1 besides water.

Note: Classes in the Process module cannot be used independent of the Component module. Therefore, it is common practice to create corresponding CompiledComponents objects and set thermo before creating Process instances.

# By default, the thermo is set with this `CompiledComponents` object upon its creation.
# We can verify that by calling the `get_thermo` function
qs.get_thermo()

Thermo(
    chemicals=CompiledComponents([
        S_I,   S_S,  X_I,  X_S,
        X_BH,  X_BA, X_P,  S_O,
        S_NO,  S_NH, S_ND, X_ND,
        S_ALK, S_N2, H2O,
    ]),
    mixture=IdealMixture(...
        include_excess_energies=False
    ),
    Gamma=DortmundActivityCoefficients,
    Phi=IdealFugacityCoefficients,
    PCF=MockPoyintingCorrectionFactors
)

# Next we need to create an instance of the ASM1 model
# Without getting into the details of ASM1, we will leave all parameters at their default values
asm1 = ASM1()
asm1.show()

ASM1([aero_growth_hetero, anox_growth_hetero, aero_growth_auto, decay_hetero, decay_auto, ammonification, hydrolysis, hydrolysis_N])

This shows that asm1 is composed of the 8 transformation processes.

1.2. Process

The Process class is developed to model individual dynamic processes. Dynamic means that the rate of reaction is finite (i.e., the reaction is not instantaneous) and often changes with the state of the system. Each of the eight transformation process in ASM1 is described by a Process object.

# Let's take the 0th process in `asm1` as an example to learn about `Process`.
# Each process in a `CompiledProcesses` object is stored as a named attribute,
# so we can grab one with normal attribute access (equivalent to `asm1.tuple[0]`):
p1 = asm1.aero_growth_hetero
p1.show()

Process: aero_growth_hetero
[stoichiometry]      S_S: -1.0/Y_H
                     X_BH: 1.00
                     S_O: 1.0*(Y_H - 1.0)/Y_H
                     S_NH: -0.0800
                     S_ALK: -0.0686
[reference]          X_BH
[rate equation]      S_NH*S_O*S_S*X_BH*mu_H/((K_N...
[parameters]         Y_H: 0.67
                     Y_A: 0.24
                     f_P: 0.08
                     mu_H: 4
                     K_S: 10
                     K_O_H: 0.2
                     K_NO: 0.5
                     b_H: 0.3
                     mu_A: 0.5
                     K_NH: 1
                     K_O_A: 0.4
                     b_A: 0.05
                     eta_g: 0.8
                     k_a: 0.05
                     k_h: 3
                     K_X: 0.1
                     eta_h: 0.8
[dynamic parameters]

Python Aside: dot access and getattr (click to expand)

obj.attr and getattr(obj, 'attr') retrieve the same value; the difference is that the second form takes the attribute name as a string, so it’s the one to reach for when the name is computed at runtime:

process_id = 'aero_growth_hetero'
asm1.aero_growth_hetero    # direct attribute access
getattr(asm1, process_id)  # equivalent, lookup by string

That matters when you want to iterate over a list of process IDs or pass an attribute name through a function. For everything else, the dot form reads more naturally.

p1 is an instance of Process describing the growth of heterotrophic biomass under aerobic conditions. Its show method prints p1’s defining characteristics, all stored as attributes on the object: stoichiometry, reference, rate equation, parameters, and dynamic parameters.

# For example, we can retrieve information on the stoichiometry of this process
p1.stoichiometry

{'S_S': -1.49253731343284,
 'X_BH': 1.00000000000000,
 'S_O': -0.492537313432836,
 'S_NH': -0.0800000000000000,
 'S_ALK': -0.0685997415522571}

Unlike what’s printed by the show method, a dictionary of evaluated stoichiometric coefficients is returned. The coefficients indicate that to produce 1g COD’s worth of X_BH, approximately 1.49g COD’s worth of S_S, 0.08g N’s worth of S_NH, and 0.49g S_O are consumed. The values of these stoichiometric coefficients are dependent on the parameter values in the formulas. We’ll see how stoichiometry will change in response to different parameter values later.

Note: The stoichiometry must be defined and interpreted in consistency with the components’ units of measure. You can check the units by calling <Component>.measured_as.

Note that X_BH has a stoichiometric coefficient of 1, which means X_BH is a “reference component” (i.e., the process rate calculated by the rate equation is also X_BH’s rate of production/consumption).

# This information can also be accessed by calling the `ref_component` property.
p1.ref_component

'X_BH'

# Another defining characteristics of a process is its rate equation, which is stored as a
# property of the `Process` object
p1.rate_equation

$\displaystyle \frac{4.0 S_{NH} S_{O} S_{S} X_{BH}}{\left(S_{NH} + 1.0\right) \left(S_{O} + 0.2\right) \left(S_{S} + 10.0\right)}$

The output above is the rate equation with parameter values already substituted in, which is convenient for evaluation but hides the parameter symbols. The four symbols that remain (\(S_{NH}\), \(S_O\), \(S_S\), \(X_{BH}\)) are the state variables this process’s rate depends on.

How do we know this function takes concentrations and not, say, mass flowrates? It comes down to the units of the parameters. To see the un-substituted form (with parameters still as symbols), Process exposes the symbolic equation via the _rate_equation attribute. The leading underscore marks it as a private/internal attribute (not part of the everyday public API), but it’s the cleanest way to surface which parameters a specific process’s rate depends on:

# The un-substituted, symbolic form of the rate equation
p1._rate_equation

$\displaystyle \frac{S_{NH} S_{O} S_{S} X_{BH} \mu_{H}}{\left(K_{NH} + S_{NH}\right) \left(K_{O H} + S_{O}\right) \left(K_{S} + S_{S}\right)}$

Now we can see \(\mu_H\), \(K_S\), \(K_{NH}\), and \(K_{O,H}\) are this process’s rate parameters. Their units aren’t stored on the Process object itself; they’re part of the model’s specification, documented in the ASM1 class docstring. We can pull just the relevant parameter lines:

# Each parameter is documented as a two-line block in `ASM1.__doc__`:
#   <name> : float, optional
#       <description>, in [<unit>]. The default is ...
# We slice out the blocks for the four parameters in p1's rate equation.
lines = ASM1.__doc__.splitlines()
for i, line in enumerate(lines):
    for p in ('mu_H ', 'K_S ', 'K_NH ', 'K_O_H '):
        if line.strip().startswith(p):
            print(line.strip())
            print('  ' + lines[i + 1].strip())
            print()
            break

mu_H : float, optional
  Heterotrophic maximum specific growth rate, in [d^(-1)]. The default

K_S : float, optional
  Readily biodegradable substrate half saturation coefficient, in [g COD/m^3].

K_O_H : float, optional
  Oxygen half saturation coefficient for heterotrophic growth, in [g O2/m^3].

K_NH : float, optional
  Ammonium (nutrient) half saturation coefficient, in [g N/m^3]. The default

Since \(\mu_H\) is in \(d^{-1}\) and \(K_S\), \(K_{NH}\), \(K_{O,H}\) are concentrations in \(g/m^3\), the rate equation produces a process rate in \(d^{-1}\) (the time derivative of \(X_{BH}\) concentration) only when \(S_{NH}\), \(S_O\), \(S_S\), and \(X_{BH}\) are also in \(g/m^3\). This is part of ASM1’s mathematical specification: biokinetic models conventionally define state variables as component concentrations.

Now that we know what units the inputs need, let’s evaluate the rate equation: substitute numerical values for the state variables (\(S_{NH}\), \(S_O\), \(S_S\), \(X_{BH}\)) and compute a single number, the process rate at that operating point. Parameter values are already baked into the formula, so we only need to supply state values.

Rate functions in qsdsan expect a 1-D array of length len(cmps) + 1: the first len(cmps) entries are component concentrations in the order given by cmps.IDs, and the last entry is the volumetric flowrate \(Q\). (\(Q\) doesn’t appear in the ASM1 process rates, but it’s part of the state because reactor-scale dynamics need it.) For a more meaningful evaluation, let’s build a snapshot of a BSM1-style aerobic reactor with realistic mixed-liquor concentrations:

import numpy as np

# Snapshot of a BSM1-style aerobic reactor. Concentrations follow the units
# documented by the ASM1 paper: COD-based components in g COD/m^3, nitrogen
# components in g N/m^3, dissolved oxygen in g O2/m^3, alkalinity in mol/m^3.
state = np.array([
    30.0,    # S_I    soluble inert
    1.5,     # S_S    readily biodegradable substrate (low: mostly consumed)
    51.0,    # X_I    particulate inert
    50.0,    # X_S    slowly biodegradable substrate
    2500.0,  # X_BH   heterotrophic biomass
    150.0,   # X_BA   autotrophic biomass
    600.0,   # X_P    inert decay products
    2.0,     # S_O    dissolved oxygen (aerated)
    8.0,     # S_NO   nitrate + nitrite
    2.0,     # S_NH   ammonium (residual)
    0.7,     # S_ND   soluble organic N
    5.0,     # X_ND   particulate organic N
    7.0,     # S_ALK  alkalinity
    15.0,    # S_N2   dinitrogen
    1000.0,  # H2O    bulk
    18446.0, # Q      volumetric flowrate [m^3/d] (BSM1 dry-weather)
])
state

array([3.000e+01, 1.500e+00, 5.100e+01, 5.000e+01, 2.500e+03, 1.500e+02,
       6.000e+02, 2.000e+00, 8.000e+00, 2.000e+00, 7.000e-01, 5.000e+00,
       7.000e+00, 1.500e+01, 1.000e+03, 1.845e+04])

# Note that the `rate_equation` attribute only stores the formula.
# The evaluation of process rate is done through the `rate_function` attribute, which is rendered
# automatically from the rate equation.
p1.rate_function(state)

790.5138339920948

At this operating point, the rate function returns about \(790~g\,COD/m^3/d\): the rate at which heterotrophic biomass (\(X_{BH}\)) is being produced per unit reactor volume per day. As a sanity check, with \(X_{BH} \approx 2500~g\,COD/m^3\) this corresponds to a specific growth rate of roughly \(790/2500 \approx 0.32~d^{-1}\), well below the \(\mu_H = 4~d^{-1}\) ceiling because \(S_S\) is the limiting substrate (its Monod term \(S_S/(K_S+S_S) = 1.5/11.5 \approx 0.13\) is by far the smallest of the three).

Note: Both stoichiometry and rate_equation are properties of the Process object. Changing the formulas of stoichiometric coefficients or rate equation are prohibited after initiation of the object.

As you may remember from the ASM1 documentation and the asm1.show() printout above, parameters and dynamic parameters are two more important attributes of a Process object. Both represent intrinsic characteristics of the process and affect the evaluation of stoichiometry and process rate, but they differ in how they’re evaluated:

• Parameters are static numerical values that don’t change during simulation (e.g., \(\mu_H = 4~d^{-1}\)). All of ASM1’s parameters are of this kind.

• Dynamic parameters are functions of state variables or time that get re-evaluated as the system state evolves. They’re useful when a stoichiometric or kinetic coefficient depends on the current operating point but can’t be folded neatly into the rate equation itself. ASM1 has none; we’ll create one in Section 3.

Let’s first look at the static parameters.

# For `asm1` and the individual processes within `asm1`, parameters are stored as a dictionary
p1.parameters

{'Y_H': 0.67,
 'Y_A': 0.24,
 'f_P': 0.08,
 'mu_H': 4.0,
 'K_S': 10.0,
 'K_O_H': 0.2,
 'K_NO': 0.5,
 'b_H': 0.3,
 'mu_A': 0.5,
 'K_NH': 1.0,
 'K_O_A': 0.4,
 'b_A': 0.05,
 'eta_g': 0.8,
 'k_a': 0.05,
 'k_h': 3.0,
 'K_X': 0.1,
 'eta_h': 0.8}

Note that p1 and asm1 share the dictionary of parameters, although not all parameters included in this dictionary are used in the stoichiometry or process rate calculations for p1. This is to make sure the parameter values are consistent across different processes in asm1.

asm1.parameters
p1.parameters is asm1.parameters

True

# If you update parameter values for `p1`, the same parameters will be updated accordingly for
# `asm1` and any other processes in `asm1`.
p1.set_parameters(mu_H=6.0, K_S=8.0, Y_H=0.8)
asm1.parameters

{'Y_H': 0.8,
 'Y_A': 0.24,
 'f_P': 0.08,
 'mu_H': 6.0,
 'K_S': 8.0,
 'K_O_H': 0.2,
 'K_NO': 0.5,
 'b_H': 0.3,
 'mu_A': 0.5,
 'K_NH': 1.0,
 'K_O_A': 0.4,
 'b_A': 0.05,
 'eta_g': 0.8,
 'k_a': 0.05,
 'k_h': 3.0,
 'K_X': 0.1,
 'eta_h': 0.8}

# If you evaluate process rate or stoichiometry again with the same input,
# you should now expect different output.
# p1.stoichiometry
p1.rate_function(state)

1435.4066985645932

Now the dynamic parameters. Since all of ASM1’s parameters are static, p1’s dynamic_parameters dict should be empty. We can verify this directly:

# `dynamic_parameters` is the public dict of `DynamicParameter` objects attached
# to this process; empty when all parameters are static (as in ASM1).
p1.dynamic_parameters

{}

In addition to the attributes above, conserved_for is another useful attribute, especially for biochemical processes. It dictates the conservation principles followed by the stoichiometric coefficients of the the process.

p1.conserved_for

('COD', 'charge', 'N')

We can see p1 is subject to conservations of COD, N, and electric charge, as defined by ASM1. This information is useful for automatic derivation of unknown stoichiometric coefficients upon initiation of the Process object. If numerical values are provided for all stoichiometric coefficients when creating the Process object, this information can then be used by check_conservation to check if these values satisfy the conservation principles.

# No return indicates that all materials in `conserved_for` are conserved.
# Otherwise, a warning or a `RuntimeError` will be raised.
p1.check_conservation()

C:\Users\Yalin\Documents\Coding\QSDsan-platform\QSDsan\qsdsan\_process.py:518: UserWarning: The following materials aren't strictly conserved by the stoichiometric coefficients. A positive value means the material is created, a negative value means the material is destroyed:
 charge: -5.20417042793042E-18
  warn("The following materials aren't strictly conserved by the "

1.3. Processes and CompiledProcesses

It is more often than not you’re working with more than one dynamic process, which often makes it impossible to derive an analytical solution to the dynamics or the steady state of the system. This is when Processes and CompiledProcesses come in handy.

As is mentioned above, Processes is simply a collection of multiple Process objects. Compiling a Processes object yields a CompiledProcesses object. Note that CompiledProcesses is also a subclass of Processes, i.e., it has all of Processes’s attributes and more.

We know that ASM1 is a subclass of CompiledProcesses, which makes asm1 an instance of CompiledProcesses.

What’s actually inside a CompiledProcesses object versus a plain Processes? compile() doesn’t add new information — it derives extra structures (the shared parameter dictionaries, the stoichiometry matrix, the compiled rate function) from the constituent Process objects and stores them on the collection. We can see the difference by rebuilding asm1 as a plain Processes from its tuple of Process objects, then re-compiling and comparing __dict__:

# Rebuild a plain `Processes` from the same `Process` objects that make up `asm1`
asm1 = qs.Processes(asm1.tuple)
asm1.show()

Processes([aero_growth_hetero, anox_growth_hetero, aero_growth_auto, decay_hetero, decay_auto, ammonification, hydrolysis, hydrolysis_N])

# A plain `Processes` only stores individual `Process` objects as attributes
asm1.__dict__

{'aero_growth_hetero': <qsdsan._process.Process at 0x245893ea450>,
 'anox_growth_hetero': <qsdsan._process.Process at 0x2458b1ff170>,
 'aero_growth_auto': <qsdsan._process.Process at 0x24589dcdfa0>,
 'decay_hetero': <qsdsan._process.Process at 0x24589e5c800>,
 'decay_auto': <qsdsan._process.Process at 0x245899e44d0>,
 'ammonification': <qsdsan._process.Process at 0x2458b830ce0>,
 'hydrolysis': <qsdsan._process.Process at 0x2458b8e6f30>,
 'hydrolysis_N': <qsdsan._process.Process at 0x2458b806300>}

A plain Processes object only stores individual Process objects as attributes.

# Re-compile to get the `CompiledProcesses` form back
asm1.compile()
asm1.__dict__

{'aero_growth_hetero': <qsdsan._process.Process at 0x245893ea450>,
 'anox_growth_hetero': <qsdsan._process.Process at 0x2458b1ff170>,
 'aero_growth_auto': <qsdsan._process.Process at 0x24589dcdfa0>,
 'decay_hetero': <qsdsan._process.Process at 0x24589e5c800>,
 'decay_auto': <qsdsan._process.Process at 0x245899e44d0>,
 'ammonification': <qsdsan._process.Process at 0x2458b830ce0>,
 'hydrolysis': <qsdsan._process.Process at 0x2458b8e6f30>,
 'hydrolysis_N': <qsdsan._process.Process at 0x2458b806300>,
 'tuple': (<qsdsan._process.Process at 0x245893ea450>,
  <qsdsan._process.Process at 0x2458b1ff170>,
  <qsdsan._process.Process at 0x24589dcdfa0>,
  <qsdsan._process.Process at 0x24589e5c800>,
  <qsdsan._process.Process at 0x245899e44d0>,
  <qsdsan._process.Process at 0x2458b830ce0>,
  <qsdsan._process.Process at 0x2458b8e6f30>,
  <qsdsan._process.Process at 0x2458b806300>),
 'size': 8,
 'IDs': ('aero_growth_hetero',
  'anox_growth_hetero',
  'aero_growth_auto',
  'decay_hetero',
  'decay_auto',
  'ammonification',
  'hydrolysis',
  'hydrolysis_N'),
 '_index': {'aero_growth_hetero': 0,
  'anox_growth_hetero': 1,
  'aero_growth_auto': 2,
  'decay_hetero': 3,
  'decay_auto': 4,
  'ammonification': 5,
  'hydrolysis': 6,
  'hydrolysis_N': 7},
 '_components': CompiledComponents([S_I, S_S, X_I, X_S, X_BH, X_BA, X_P, S_O, S_NO, S_NH, S_ND, X_ND, S_ALK, S_N2, H2O]),
 '_parameters': {'Y_H': 0.8,
  'Y_A': 0.24,
  'f_P': 0.08,
  'mu_H': 6.0,
  'K_S': 8.0,
  'K_O_H': 0.2,
  'K_NO': 0.5,
  'b_H': 0.3,
  'mu_A': 0.5,
  'K_NH': 1.0,
  'K_O_A': 0.4,
  'b_A': 0.05,
  'eta_g': 0.8,
  'k_a': 0.05,
  'k_h': 3.0,
  'K_X': 0.1,
  'eta_h': 0.8},
 '_dyn_params': {},
 '_stoichiometry': [[0,
   -1.25000000000000,
   0,
   0,
   1.00000000000000,
   0,
   0,
   -0.250000000000000,
   0,
   -0.0800000000000000,
   0,
   0,
   -0.0685997415522571,
   0,
   0],
  [0,
   -1.25000000000000,
   0,
   0,
   1.00000000000000,
   0,
   0,
   0,
   -0.0875000000000000,
   -0.0800000000000001,
   0,
   0,
   0.00643122577052393,
   0.0875000000000000,
   0],
  [0,
   0,
   0,
   0,
   0,
   1.00000000000000,
   0,
   -18.0476190476190,
   4.16666666666667,
   -4.24666666666667,
   0,
   0,
   -7.21440615324570,
   0,
   0],
  [0,
   0,
   0,
   0.920000000000000,
   -1.00000000000000,
   0,
   0.0800000000000000,
   0,
   0,
   0,
   0,
   0.0752000000000000,
   0,
   0,
   0],
  [0,
   0,
   0,
   0.920000000000000,
   0,
   -1.00000000000000,
   0.0800000000000000,
   0,
   0,
   0,
   0,
   0.0752000000000000,
   0,
   0,
   0],
  [0,
   0,
   0,
   0,
   0,
   0,
   0,
   0,
   0,
   1.00000000000000,
   -1.00000000000000,
   0,
   0.857496769403214,
   0,
   0],
  [0, 1.0, 0, -1.0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],
  [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1.0, -1.0, 0, 0, 0]],
 '_stoichio_lambdified': None,
 '_rate_equations': (S_NH*S_O*S_S*X_BH*mu_H/((K_NH + S_NH)*(K_O_H + S_O)*(K_S + S_S)),
  K_O_H*S_NH*S_NO*S_S*X_BH*eta_g*mu_H/((K_NH + S_NH)*(K_NO + S_NO)*(K_O_H + S_O)*(K_S + S_S)),
  S_NH*S_O*X_BA*mu_A/((K_NH + S_NH)*(K_O_A + S_O)),
  X_BH*b_H,
  X_BA*b_A,
  S_ND*X_BH*k_a,
  X_BH*X_S*k_h*(K_O_H*S_NO*eta_h/((K_NO + S_NO)*(K_O_H + S_O)) + S_O/(K_O_H + S_O))/(K_X*X_BH + X_S),
  X_BH*X_ND*k_h*(K_O_H*S_NO*eta_h/((K_NO + S_NO)*(K_O_H + S_O)) + S_O/(K_O_H + S_O))/(K_X*X_BH + X_S)),
 '_production_rates': [0,
  -1.25*K_O_H*S_NH*S_NO*S_S*X_BH*eta_g*mu_H/((K_NH + S_NH)*(K_NO + S_NO)*(K_O_H + S_O)*(K_S + S_S)) - 1.25*S_NH*S_O*S_S*X_BH*mu_H/((K_NH + S_NH)*(K_O_H + S_O)*(K_S + S_S)) + 1.0*X_BH*X_S*k_h*(K_O_H*S_NO*eta_h/((K_NO + S_NO)*(K_O_H + S_O)) + S_O/(K_O_H + S_O))/(K_X*X_BH + X_S),
  0,
  0.92*X_BA*b_A - 1.0*X_BH*X_S*k_h*(K_O_H*S_NO*eta_h/((K_NO + S_NO)*(K_O_H + S_O)) + S_O/(K_O_H + S_O))/(K_X*X_BH + X_S) + 0.92*X_BH*b_H,
  1.0*K_O_H*S_NH*S_NO*S_S*X_BH*eta_g*mu_H/((K_NH + S_NH)*(K_NO + S_NO)*(K_O_H + S_O)*(K_S + S_S)) + 1.0*S_NH*S_O*S_S*X_BH*mu_H/((K_NH + S_NH)*(K_O_H + S_O)*(K_S + S_S)) - 1.0*X_BH*b_H,
  1.0*S_NH*S_O*X_BA*mu_A/((K_NH + S_NH)*(K_O_A + S_O)) - 1.0*X_BA*b_A,
  0.08*X_BA*b_A + 0.08*X_BH*b_H,
  -0.25*S_NH*S_O*S_S*X_BH*mu_H/((K_NH + S_NH)*(K_O_H + S_O)*(K_S + S_S)) - 18.047619047619*S_NH*S_O*X_BA*mu_A/((K_NH + S_NH)*(K_O_A + S_O)),
  -0.0875*K_O_H*S_NH*S_NO*S_S*X_BH*eta_g*mu_H/((K_NH + S_NH)*(K_NO + S_NO)*(K_O_H + S_O)*(K_S + S_S)) + 4.16666666666667*S_NH*S_O*X_BA*mu_A/((K_NH + S_NH)*(K_O_A + S_O)),
  -0.0800000000000001*K_O_H*S_NH*S_NO*S_S*X_BH*eta_g*mu_H/((K_NH + S_NH)*(K_NO + S_NO)*(K_O_H + S_O)*(K_S + S_S)) + 1.0*S_ND*X_BH*k_a - 0.08*S_NH*S_O*S_S*X_BH*mu_H/((K_NH + S_NH)*(K_O_H + S_O)*(K_S + S_S)) - 4.24666666666667*S_NH*S_O*X_BA*mu_A/((K_NH + S_NH)*(K_O_A + S_O)),
  -1.0*S_ND*X_BH*k_a + 1.0*X_BH*X_ND*k_h*(K_O_H*S_NO*eta_h/((K_NO + S_NO)*(K_O_H + S_O)) + S_O/(K_O_H + S_O))/(K_X*X_BH + X_S),
  0.0752*X_BA*b_A - 1.0*X_BH*X_ND*k_h*(K_O_H*S_NO*eta_h/((K_NO + S_NO)*(K_O_H + S_O)) + S_O/(K_O_H + S_O))/(K_X*X_BH + X_S) + 0.0752*X_BH*b_H,
  0.00643122577052393*K_O_H*S_NH*S_NO*S_S*X_BH*eta_g*mu_H/((K_NH + S_NH)*(K_NO + S_NO)*(K_O_H + S_O)*(K_S + S_S)) + 0.857496769403214*S_ND*X_BH*k_a - 0.0685997415522571*S_NH*S_O*S_S*X_BH*mu_H/((K_NH + S_NH)*(K_O_H + S_O)*(K_S + S_S)) - 7.2144061532457*S_NH*S_O*X_BA*mu_A/((K_NH + S_NH)*(K_O_A + S_O)),
  0.0875*K_O_H*S_NH*S_NO*S_S*X_BH*eta_g*mu_H/((K_NH + S_NH)*(K_NO + S_NO)*(K_O_H + S_O)*(K_S + S_S)),
  0],
 '_rate_function': None}

A CompiledProcesses object has many more attributes that streamline the use of dynamic process models. We didn’t supply any extra information to compile(): every new attribute is derived from the constituent Process objects.

Note: a ``Process`` belongs to one ``CompiledProcesses`` at a time. Unlike Component, which can be reused freely across different CompiledComponents sets, compile() rewires each constituent process’s parameter dictionary to point at the new collection’s. If you put the same Process object into a second Processes and compile that, the process is then bound to the second collection, and set_parameters on the first no longer propagates to it. In practice, build a fresh Process (or copy it with Process.copy) when you want it in a different CompiledProcesses.

The .stoichiometry attribute of a CompiledProcesses returns a table with processes as rows, components as columns, and the stoichiometric coefficients in the cells. This layout is known as the Petersen (Gujer) matrix. It is widely used in biokinetic models such as Activated Sludge Model (ASM1/2/3) and the Anaerobic Digestion Model No. 1 (Tutorial 12). Calling the matrix \(\mathbf{A}\) (rows = processes, columns = components), the per-component net production or consumption rate is \(\mathbf{r} = \mathbf{A}^{T} \boldsymbol{\rho}\), with \(\boldsymbol{\rho}\) the vector of process rates — exactly the form CompiledProcesses.rate_function evaluates below.

# `.stoichiometry` returns the Petersen matrix described above
asm1.stoichiometry

	S_I	S_S	X_I	X_S	X_BH	...	S_ND	X_ND	S_ALK	S_N2	H2O
aero_growth_hetero	0	-1.25	0	0	1	...	0	0	-0.0686	0	0
anox_growth_hetero	0	-1.25	0	0	1	...	0	0	0.00643	0.0875	0
aero_growth_auto	0	0	0	0	0	...	0	0	-7.21	0	0
decay_hetero	0	0	0	0.92	-1	...	0	0.0752	0	0	0
decay_auto	0	0	0	0.92	0	...	0	0.0752	0	0	0
ammonification	0	0	0	0	0	...	-1	0	0.857	0	0
hydrolysis	0	1	0	-1	0	...	0	0	0	0	0
hydrolysis_N	0	0	0	0	0	...	1	-1	0	0	0

8 rows × 15 columns

# Similarly for the rate equations
asm1.rate_equations

	rate_equation
aero_growth_hetero	6.0*S_NH*S_O*S_S*X_BH/((S_NH + ...
anox_growth_hetero	0.96*S_NH*S_NO*S_S*X_BH/((S_NH ...
aero_growth_auto	0.5*S_NH*S_O*X_BA/((S_NH + 1.0)...
decay_hetero	0.3*X_BH
decay_auto	0.05*X_BA
ammonification	0.05*S_ND*X_BH
hydrolysis	3.0*X_BH*X_S*(0.16*S_NO/((S_NO ...
hydrolysis_N	3.0*X_BH*X_ND*(0.16*S_NO/((S_NO...

# More importantly, the `rate_function` attribute of a `CompiledProcesses` now outputs an array
# with each element corresponding orderly to the individual processes
asm1.rate_function(state)

array([1435.407,  108.078,   41.667,  750.   ,    7.5  ,   87.5  ,
       1221.925,  122.193])

With this, the evaluation of all process rates can be done with one command and the output data type is convenient for matrix operations.

To calculate the overall rate of production for a component, we need to first estimate its rates of production (positive means being produced, negative means being consumed) by each process and then sum them up. In mathematical forms, let’s say \(a_{ij}\) represents the stoichiometric coefficient of component \(j\) in process \(i\) (i.e., value on the \(i\) th row and \(j\) th column of the stoichiometry matrix). Denote \(\rho_i\) as process \(i\)’s reaction rate, now the overall production rate of component \(j\) should be

\[r_j = \sum_i{a_{ij}\cdot\rho_i}\]

In matrix notation, this calculation can be neatly described as

\[\mathbf{r} = \mathbf{A^T} \mathbf{\rho}\]

where \(\mathbf{A}\) is the stoichiometry matrix and \(\mathbf{\rho}\) is the array of process rates.

# This matrix operation is already streamlined for `CompiledProcesses` objects, you can see the
# mathematical form of the rates of production as a function of component concentrations.
asm1.production_rates

	rate_of_production
S_I	0
S_S	-1.2*S_NH*S_NO*S_S*X_BH/((S_NH ...
X_I	0
X_S	0.046*X_BA - 3.0*X_BH*X_S*(0.16...
X_BH	0.96*S_NH*S_NO*S_S*X_BH/((S_NH ...
X_BA	0.5*S_NH*S_O*X_BA/((S_NH + 1.0)...
X_P	0.004*X_BA + 0.024*X_BH
S_O	-1.5*S_NH*S_O*S_S*X_BH/((S_NH +...
S_NO	-0.084*S_NH*S_NO*S_S*X_BH/((S_N...
S_NH	0.05*S_ND*X_BH - 0.076800000000...
S_ND	-0.05*S_ND*X_BH + 3.0*X_BH*X_ND...
X_ND	0.00376*X_BA - 3.0*X_BH*X_ND*(0...
S_ALK	0.0428748384701607*S_ND*X_BH + ...
S_N2	0.084*S_NH*S_NO*S_S*X_BH/((S_NH...
H2O	0

# To evaluate the rates of production for all components, all you need to
# do is to call the `production_rates_eval` method.
asm1.production_rates_eval(state)

array([    0.   ,  -707.43 ,     0.   ,  -525.025,   793.484,    34.167,
          60.6  , -1110.836,   164.154,  -212.923,    34.693,   -65.229,
        -323.343,     9.457,     0.   ])

With these convenient evaluation methods, asm1 can now be used in SanUnit to model system dynamics (e.g., reactor effluent concentrations over time). See the Dynamic Simulation tutorial for how to integrate process models with SanUnit subclasses and enable dynamic simulations.

Congratulations! You’ve understood the structure and the core capacities of the _process module. 🎉

2. Making a simple biokinetic model

In this section, we will get some hands-on experience with the _process module. For demonstration purpose, we will implement a simple biokinetic model with the classes we learned in the previous section. It’s even better if you have your own mathematical model of dynamic processes in mind.

Let’s get started with defining the mathematical model.

The biokinetic model

Let’s consider a simplified activated sludge system with 5 components besides water:

• particulate organic substrate (\(X_S\), measured as COD)

• soluble organic substrate (\(S_S\), measured as COD)

• dissolved oxygen (\(O_2\), measured as itself)

• dissolved CO2 (\(CO_2\), measured as C)

• biomass (\(X_B\), measured as COD)

This model uses 3 dynamic processes to describe the transformations among the components:

1. Hydrolysis (i.e., \(X_S\) transforming into \(S_S\))

2. Growth of biomass (i.e., \(X_B\) utilizing \(S_S\) and \(O_2\) to grow)

3. Decay of biomass (i.e., \(X_B\) dying and transforming into \(X_S\))

To mathematically characterize these processes, we need to define the stoichiometry and the process rate equation of each process. Let’s say, for process 1 and 3:

\[X_S \Rightarrow S_S,~~\rho_1 = k_{hyd}X_S\]

\[X_B \Rightarrow X_S,~~\rho_3 = b\cdot X_B\]

where \(k_{hyd}\) is the hydrolysis rate constant, \(b\) is the biomass decay rate constant, both in \(d^{-1}\cdot m^3\cdot g^{-1}\).

For biomass growth, assuming the yield of biomass is \(y_B\) [gCOD \(X_B\)/gCOD \(S_S\)], then the stoichiometry can be defined as:

\[\frac{1}{y_B}S_S + a_1O_2 \Rightarrow X_B + a_2CO_2\]

where \(a_1\) and \(a_2\) need to be determined by conservations of COD and carbon. Similar to other activated sludge models, we can use a Monod-type function to represent biomass growth rate, i.e.,

\[\rho_2 = \mu_B\cdot \frac{S_S}{K_S+S_S}X_B\]

2.1. Defining the components

For convenience, we will find similar components that have been defined in the default set. Note that we will need to consider conservations of COD and C in the stoichiometry of process 2. Therefore, we should make sure the Component objects we pick have appropriate values for COD content (i_COD) and carbon content (i_C).

# Load the default set of components
cmps_all = qs.Components.load_default()

# Copy select components from the default set and give them new IDs
X_S = cmps_all.X_B_Subst.copy('X_S')
S_S = cmps_all.S_F.copy('S_S')
O2 = cmps_all.S_O2.copy('O2')
CO2 = cmps_all.S_CO3.copy('CO2')
X_B = cmps_all.X_OHO.copy('X_B')

# Compile them into a new set. Don't forget the water!
cmps_bkm = qs.Components([X_S, S_S, O2, CO2, X_B, cmps_all.H2O])
cmps_bkm.compile(ignore_inaccurate_molar_weight=True)
cmps_bkm.show()

CompiledComponents([
    X_S, S_S, O2, CO2,
    X_B, H2O,
])

Note: You might’ve noticed H2O isn’t an active component in any process but is always included in the CompiledComponents set. This is because water is assumed to be the bulk liquid and is needed for the convergence of total volumetric flows in simulation regardless of the process model.

# Now we can check if their `measured_as` attributes are correctly set
(X_S.measured_as == 'COD',
 CO2.measured_as == 'C',
 O2.measured_as is None)

(True, True, True)

# Then you can check relevant `i_` properties of the components
# For example, O2 should have a negative COD content, or more specifically, -1 gCOD/gO2
O2.i_COD

-1.0

Tip: It is not always/just i_COD and i_C that are important. For example, if a model involves nitrification/denitrification processes, it is critical that the components’ i_N values are set correctly. It ultimately boils down to what conservation principles need to be considered.

# A final preparation step is to set thermo with the components
qs.set_thermo(cmps_bkm)

It is common practice to wrap up the creation of model-specific CompiledComponents objects in a function for convenient reuse. The create_<ModelName>_cmps functions we’ve seen in the previous section is a good example.

2.2. Defining the processes

With the components ready, we can now define individual dynamic processes in the biokinetic model as Process objects.

For process 1 and process 3, the stoichiometric coefficients are known and the rate equations take simple forms with just one parameter.

# To define `Process` objects for these two processes, we just need to input their
# stoichiometries and rate equations

# Process 1: hydrolysis
pc1 = qs.Process(ID='hydrolysis',
                reaction='X_S -> S_S',
                ref_component='X_S',
                rate_equation='k_hyd*X_S',
                conserved_for=(),
                parameters=('k_hyd',))

Note that the information on stoichiometry is reflected in the reaction argument. ref_component is set to be X_S because its stoichiometric coefficient is -1.

No conservation principle is used in the initiation of pc1 (i.e., conserved_for is an empty tuple). This is fine, because there’s no unknown stoichiometric coefficient that relies on conservation principles to solve. However, this doesn’t mean the process isn’t subject to any conservation principle. The 1 to 1 stoichiometry of this process reflects the conservation of COD, as both \(X_S\) and \(S_S\) are measured as COD.

# We can verify this with the `check_conservation` method
# To do that, `conserved_for` needs to be set first
pc1.conserved_for = ('COD',)
pc1.check_conservation()

No return means pc1 conserves COD.

Note: conserved_for is treated as an iterable in the check_conservation method. So, when there’s only one conservation principle, we should make sure the entire string is recognized as one element in iteration.

# For example, if we do this, we will get an error message because the function interprets it as
# checking conservations of 'C', 'O', and 'D', respectively
pc1.conserved_for = ('COD') # note the missing comma, which makes it a string instead of a tuple
pc1.check_conservation()

---------------------------------------------------------------------------
AttributeError                            Traceback (most recent call last)
File ~\Documents\Coding\QSDsan-platform\QSDsan\qsdsan\_process.py:608, in Process.conserved_for(self, materials)
    607 for mat in materials:
--> 608     try: get(Component, f'i_{mat}')
    609     except: raise ValueError(f'Components do not have i_{mat} attribute.')

AttributeError: type object 'Component' has no attribute 'i_O'

During handling of the above exception, another exception occurred:

ValueError                                Traceback (most recent call last)
Cell In[38], line 3
      1 # For example, if we do this, we will get an error message because the function interprets it as
      2 # checking conservations of 'C', 'O', and 'D', respectively
----> 3 pc1.conserved_for = ('COD') # note the missing comma, which makes it a string instead of a tuple
      4 pc1.check_conservation()

File ~\Documents\Coding\QSDsan-platform\QSDsan\qsdsan\_process.py:609, in Process.conserved_for(self, materials)
    607 for mat in materials:
    608     try: get(Component, f'i_{mat}')
--> 609     except: raise ValueError(f'Components do not have i_{mat} attribute.')
    610 self._conserved_for = materials

ValueError: Components do not have i_O attribute.

Note that the initiation also asks for parameters (k_hyd in this case). This is to help identify any symbolic parameter in the stoichiometry and the rate equation.

# Upon initiation, the parameters are stored as symbols. We still need to set values to them
# before we can evaluate process rate.
pc1.show()

Process: hydrolysis
[stoichiometry]      X_S: -1
                     S_S: 1
[reference]          X_S
[rate equation]      X_S*k_hyd
[parameters]         k_hyd: k_hyd
[dynamic parameters]

# At this point, the initiation of process 3 should be quite straightforward.
# Here shows an alternative way to input stoichiometry
pc3 = qs.Process('decay',
                reaction={'X_B':-1, 'X_S':1},
                ref_component='X_B',
                rate_equation='b*X_B',
                conserved_for=('COD',),
                parameters=('b',))

# Note that when initiating a process with unknown stoichiometric coefficients,
# we can use `?` as placeholders in the reaction string.
pc2 = qs.Process('growth',
                reaction='[1/y_B]S_S + [?]O2 -> X_B + [?]CO2',
                ref_component='X_B',
                rate_equation='mu_B * S_S/(K_S+S_S) * X_B',
                conserved_for=('COD', 'C'),
                parameters=('y_B', 'mu_B', 'K_S'))
pc2.show()

Process: growth
[stoichiometry]      S_S: -1/y_B
                     O2: (y_B - 1.0)/y_B
                     CO2: 0.002*(160.0 - 183.0*y_B)/y_B
                     X_B: 1.00
[reference]          X_B
[rate equation]      S_S*X_B*mu_B/(K_S + S_S)
[parameters]         y_B: y_B
                     mu_B: mu_B
                     K_S: K_S
[dynamic parameters]

Because there’re two unknown stoichiometric coefficients, we need two conservation principles (i.e., two equations) to determine their values. It is convenient to solve for stoichiometric coefficients marked as ? by specifying conserved_for in the initiation of Process objects.

2.3. Compiling the processes

# Now the final step is to compile the individual processes into a biokinetic model
bkm = qs.Processes([pc1, pc2, pc3])
bkm.compile()
bkm.show()

CompiledProcesses([hydrolysis, growth, decay])

Warning: Within a single CompiledProcesses, parameters across all of its constituent Process objects are merged into a single shared dictionary on compile. If two processes in the same collection happen to declare the same symbol with different intended meanings, the second one silently overwrites the first. Use distinct parameter symbols within a collection (e.g., K_S_growth vs. K_S_decay, not K_S in both). However, reusing the same symbol in separate CompiledProcesses collections is fine: each collection has its own parameter dictionary.

# All process parameters are merged into one shared dictionary on compile
bkm.parameters

{'k_hyd': k_hyd, 'y_B': y_B, 'mu_B': mu_B, 'K_S': K_S, 'b': b}

# After setting parameter values, the model will be ready
bkm.set_parameters(k_hyd=3.0, y_B=0.8, mu_B=4.0, K_S=9.0, b=0.4)
bkm.stoichiometry

	X_S	S_S	O2	CO2	X_B	H2O
hydrolysis	-1	1	0	0	0	0
growth	0	-1.25	-0.25	0.034	1	0
decay	1	0	0	0	-1	0

bkm.production_rates

	rate_of_production
X_S	0.4*X_B - 3.0*X_S
S_S	-5.0*S_S*X_B/(S_S + 9.0) + 3.0*X_S
O2	-1.0*S_S*X_B/(S_S + 9.0)
CO2	0.136*S_S*X_B/(S_S + 9.0)
X_B	4.0*S_S*X_B/(S_S + 9.0) - 0.4*X_B
H2O	0

2.4. Batch mode for creating processes

It’d be very tedious to create individual Process objects one by one and compile them, especially for models with a large number of parallel processes and/or components. Therefore, the _process module includes convenient classmethods to create processes in batch.

# Stoichiometry and rate equations are commonly stored in tabular form
# (.csv, .tsv, or Excel). `Processes.load_from_file` consumes such a table
# directly. Here's the example for our biokinetic model:
from qsdsan.utils import load_data
df_bkm = load_data('assets/_bkm.csv', index_col=0)
df_bkm

	X_S	S_S	O2	CO2	X_B	H2O	Unnamed: 7
hydrolysis	-1	1	0	0	0	0	k_hyd*X_S
growth	0	(-1)/y_B	?	?	1	0	mu_B*S_S/(K_S + S_S)*X_B
decay	1	0	0	0	-1	0	b*X_B

This table contains the stoichiometry matrix and the array of rate equations. If the table is organized with process IDs as row indices and component IDs as column names, we can use it as direct input to the Processes.load_from_file classmethod to create processes in batch.

# `conserved_for` is required: pass a tuple to apply the same rules to every
# process. (We'll see two more forms below.)
bkm_batch = qs.Processes.load_from_file(
    path='assets/_bkm.csv',
    conserved_for=('COD', 'C'),
    parameters=('k_hyd', 'y_B', 'mu_B', 'K_S', 'b'),
    compile=True,
)
bkm_batch.show()

CompiledProcesses([hydrolysis, growth, decay])

# You can see `bkm_batch` is equivalent to the `bkm` we created above
bkm_batch.stoichiometry
# bkm_batch.rate_equations
# bkm_batch.parameters

	X_S	S_S	O2	CO2	X_B	H2O
hydrolysis	-1	1	0	0	0	0
growth	0	-1.0/y_B	1.0*(y_B - 1.0)/y_B	0.002*(160.0 - 183.0*y_B)/y_B	1.00000000000000	0
decay	1	0	0	0	-1	0

# The reference component of each process is inferred from its stoichiometry
bkm_batch.decay.ref_component

'X_S'

Notice that we passed a single tuple conserved_for=('COD', 'C') to load_from_file above. Unlike the manual construction in section 2.2 (where each Process could carry its own conserved_for), the tuple form of load_from_file applies the same conservation rules to every process in the batch. That’s convenient but may not be accurate: the decay process, for instance, doesn’t actually conserve carbon in this model. We can see the consequence by running check_conservation on the decay process:

# Decay's stoichiometry doesn't balance carbon, so check_conservation
# raises a RuntimeError when called against the uniformly-applied
# ("COD", "C") rules.
bkm_batch.decay.check_conservation()

---------------------------------------------------------------------------
RuntimeError                              Traceback (most recent call last)
Cell In[49], line 4
      1 # Decay's stoichiometry doesn't balance carbon, so check_conservation
      2 # raises a RuntimeError when called against the uniformly-applied
      3 # ("COD", "C") rules.
----> 4 bkm_batch.decay.check_conservation()

File ~\Documents\Coding\QSDsan-platform\QSDsan\qsdsan\_process.py:505, in Process.check_conservation(self, rtol, atol)
    503         materials = self._conserved_for
    504         unconserved = [(materials[i], ic_dot_v[i]) for i, conserved in enumerate(conserved_arr) if not conserved]
--> 505         raise RuntimeError("The following materials are unconserved by the "
    506                            "stoichiometric coefficients. A positive value "
    507                            "means the material is created, a negative value "
    508                            "means the material is destroyed:\n "
    509                            + "\n ".join([f"{material}: {value:.2f}" for material, value in unconserved]))
    510 elif isa(v, list):
    511     if ic is None:

RuntimeError: The following materials are unconserved by the stoichiometric coefficients. A positive value means the material is created, a negative value means the material is destroyed:
 C: -0.05

The error above tells us decay conserves COD but not carbon under the current i_C values. Actually, we wanted decay checked only for COD (matching the manual pc3 in section 2.2). To recover per-process rules without leaving the batch workflow, pass conserved_for as a dict keyed by process ID:

# `i_C` values for each component,
# if you want to manually check the `i_C` values for all components in the model
dict(zip(cmps_bkm.IDs, cmps_bkm.i_C))

{'X_S': 0.32, 'S_S': 0.32, 'O2': 0.0, 'CO2': 1.0, 'X_B': 0.366, 'H2O': 0.0}

# Dict form: per-process conservation rules. Process IDs absent from the dict
# fall back to the default ('COD', 'N', 'P', 'charge').
bkm_batch = qs.Processes.load_from_file(
    path='assets/_bkm.csv',
    conserved_for={
        'hydrolysis': (),         # known 1-to-1 COD swap, no further check needed
        'growth': ('COD', 'C'),   # solves the two `?` coefficients
        'decay': ('COD',),        # carbon isn't tracked through decay here
    },
    parameters=('k_hyd', 'y_B', 'mu_B', 'K_S', 'b'),
    compile=True,
)
# decay now passes its own (COD-only) conservation check
bkm_batch.decay.check_conservation()

No warning this time: each process is checked against the rules we actually want.

Embedding rules in the file

When the conservation rules are an intrinsic part of the model rather than something to tweak per call, you can carry them inside the data file itself. Add a column named conserved_for whose entries are comma-separated material names (empty cell = no conservation). assets/_bkm_with_conserved.csv is the same biokinetic model with that column appended:

df_bkm_cf = load_data('assets/_bkm_with_conserved.csv', index_col=0)
df_bkm_cf

	X_S	S_S	O2	CO2	X_B	H2O	conserved_for	Unnamed: 8
hydrolysis	-1	1	0	0	0	0	NaN	k_hyd*X_S
growth	0	(-1)/y_B	?	?	1	0	COD,C	mu_B*S_S/(K_S + S_S)*X_B
decay	1	0	0	0	-1	0	COD	b*X_B

# Passing `conserved_for=None` tells `load_from_file` to use the column
# from the file. Processes without a value in the column get the default tuple.
bkm_batch = qs.Processes.load_from_file(
    path='assets/_bkm_with_conserved.csv',
    conserved_for=None,
    parameters=('k_hyd', 'y_B', 'mu_B', 'K_S', 'b'),
    compile=True,
)
# decay was assigned ('COD',) directly by the file column
bkm_batch.decay.conserved_for, bkm_batch.decay.check_conservation()

(('COD',), None)

Same model, same per-process rules as the dict form, but the rules now live alongside the stoichiometry and rate equations. Use whichever form fits your workflow: the tuple form for quick uniform application, the dict form for per-call overrides, and the file column when the rules belong with the model definition.

3. Advanced features

To provide a flexible and user-friendly platform for process modeling, many advanced features are included in the _process module to streamline the implementation of process models with various mathematical complexity.

Let’s again use the biokinetic model above as an example to learn about these features.

3.1. DynamicParameter

What if the yield of biomass growth is dependent on the abundance of soluble organic substrate? Assume this relation can be described as

\[y_B = y_{max}\cdot\sqrt{\frac{S_S}{X_S + S_S}}\]

In this case, we will need to define \(y_B\) as a DynamicParameter object because its value changes with system state during simulation, and we will use the decorators to do so.

@pc2.dynamic_parameter(symbol='y_B', params={'y_max':0.8})
def yB_eval(state_arr, params):
    X_S, S_S = state_arr[:2]  # assuming the first 6 elements are component concentrations
    y_max = params['y_max']
    y_B = y_max * (S_S/(X_S+S_S))**(1/2)
    return y_B

Three things are happening in the cell above:

• pc2.dynamic_parameter(...) is a decorator factory: it collects the arguments below and returns the actual decorator that wraps yB_eval. (For a refresher on decorators and decorator factories, see the Python Aside in the SanUnit (advanced) tutorial.)

• symbol='y_B' names the dynamic parameter. It must match the symbol used in the process’s stoichiometry or rate equation (here, y_B from the [1/y_B]S_S + ... stoichiometry in pc2).

• params={'y_max': 0.8} holds fixed values that the evaluation function needs but that do not change during simulation. They are handed back to the function on every call.

The decorated function yB_eval accepts up to two positional arguments (it raises ValueError if you declare more):

• state_arr: a numpy array of current state-variable values (component concentrations followed by \(Q\), in the order of cmps.IDs). Slice it to pick out the components your formula needs.

• params: the same dict that was passed to the decorator. Look up params['y_max'] (and any other fixed inputs).

It returns a single number — the current value of y_B — which QSDsan splices into the process’s stoichiometry and rate equation each time params_eval (or production_rates_eval) runs.

# At this point, the value for `y_B` is not updated yet, since we haven't evaluated it
# with input of component concentrations
bkm.parameters

{'k_hyd': 3.0, 'y_B': 0.8, 'mu_B': 4.0, 'K_S': 9.0, 'b': 0.4}

# But the list of dynamic parameters have been updated
pc2.show()

Process: growth
[stoichiometry]      S_S: -1/y_B
                     O2: (y_B - 1.0)/y_B
                     CO2: 0.002*(160.0 - 183.0*y_B)/y_B
                     X_B: 1.00
[reference]          X_B
[rate equation]      S_S*X_B*mu_B/(K_S + S_S)
[parameters]         k_hyd: 3
                     y_B: 0.8
                     mu_B: 4
                     K_S: 9
                     b: 0.4
[dynamic parameters] <DynamicParameter: y_B>

# `y_B` is now a `DynamicParameter` object stored in the `dynamic_parameters` dict
type(pc2.dynamic_parameters['y_B'])

qsdsan._process.DynamicParameter

# Assuming component concentrations are all 1.
state_bkm = np.ones(6)

# This function evaluates the dynamic parameters
bkm.params_eval(state_bkm)
bkm.parameters

{'k_hyd': 3.0, 'y_B': 0.5656854249492381, 'mu_B': 4.0, 'K_S': 9.0, 'b': 0.4}

# Then the stoichiometry is updated accordingly
bkm.stoichiometry

	X_S	S_S	O2	CO2	X_B	H2O
hydrolysis	-1	1	0	0	0	0
growth	0	-1.77	-0.768	0.2	1	0
decay	1	0	0	0	-1	0

Note that in simulations, you don’t need to call the params_eval method explicitly unless you would like to capture the dynamics of parameter values. The evaluation of dynamic parameters is performed within production_rates_eval (i.e., the calculation of overall production rates of the components).

3.2. Kinetics and MultiKinetics

Kinetics and MultiKinetics are classes that allow defining process rates as functions.

Kinetics is a subclass of DynamicParameter, because process rate in essence is also a function of state variables and additional static parameters.

# The use of `Kinetics` to define process rate is also similar
# For example, for the decay process
@pc3.kinetics(parameters={'b':0.4})
def rho3_eval(state_arr, params):
    X_B = state_arr[4]
    b = params['b']
    return b*X_B

pc3.rate_function

<Kinetics: decay>

This may seem trivial, but it’s helpful when the form of process rate gets more complex. For example, the rate of substrate uptake in anaerobic digestion model no.1 is dependent on pH, which is dynamically determined by a linear equation system and thus cannot be expressed neatly as a function of state variables.

MultiKinetics serves a similar purpose but calculates the rate of multiple processes at once. It is designed to be used with CompiledProcesses objects. This is convenient for models with similar process rate algorithms across multiple processes in the sense that vectorizing the shared algorithm can often improve computational speed in simulation. Another application of MultiKinetics is models sharing intermediate variables (e.g., pH) across multiple process rate calculation.

Let’s assume the reaction rates of all three processes in the biokinetic model above are now dependent on the abundance of oxygen. Here are the new rate equations:

\[\rho_1 = k_{hyd}X_S\cdot \sqrt{\frac{O_2}{O_2+K_{O2}}}\]

\[\rho_2 = \mu_BX_B\cdot \frac{S_S}{K_S+S_S} \sqrt{\frac{O_2}{O_2+K_{O2}}}\]

\[\rho_3 = b\cdot X_B\cdot (1.5- \sqrt{\frac{O_2}{O_2+K_{O2}}})\]

How do we use MultiKinetics to implement this change?

# Define the function
def rhos_eval(state_arr, params):
    X_S, S_S, O2, CO2, X_B = state_arr[:5]
    I_O2 = (O2/(O2+params['K_O2']))**(1/2)
    rho1 = params['k_hyd']*X_S*I_O2
    rho2 = params['mu_B']*X_B*I_O2
    rho3 = params['b']*X_B*(1.5-I_O2)
    return np.array([rho1, rho2, rho3])

# Set it as the rate function of the `CompiledProcesses` object
bkm.set_rate_function(rhos_eval)

# Assign parameter values
bkm.rate_function.set_params(K_O2=2.0, **bkm.parameters)

# The evaluation of the rate function is the same as before
bkm.rate_function(state_bkm)

array([1.732, 2.309, 0.369])

---

↑ Back to top

↑ Back to top