Activated Sludge System with ASM1

In this example, a small activated sludge system with the activated sludge system nr. 1 (ASM1) will be built. An activated sludge system is a simple setup with a single CSTR:

For building and simulating this in BioChemicalTreatment.jl it has to be noted first, that this actually needs a series of systems, some of which are invisible in the plot above. For example the unifying of the inflow and return sludge in the beginning needs an own system, and the splitting of the wastage stream needs even two: A FlowSeparator to split the flow and then a FlowPump to specify the flow rate for the split. Introducing all of those yields a new diagram:

Now with all these required systems, the building of the model can be started.

Plug and Play Example

using BioChemicalTreatment # Reactors etc.
using ModelingToolkit # Modeling Framework
using DifferentialEquations # Solve the equations
using Plots # Plotting

# Keeping track
sys = []
eqns = Equation[]

# Default Process
@set_process asm1 = Process("ASM1")

# The reactor
@named cstr = CSTR(1000; initial_states=[1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1]) # 1000 m^3 is the assumed volume, we set the initial states to all 1 as the equations do not work with the default (all 0)
push!(sys, cstr)

# The other systems
@named influent = Influent([7.0, 69.5, 6.95, 0.0, 31.56, 0.0, 0.0, 30.0, 202.32, 10.59, 0.0, 28.17, 0.0, 51.2], flowrate = 18446) # This is the stabilization influent from BSM1
push!(sys, influent)

# The unifier with two inflows
@named unifier = FlowUnifier()
push!(sys, unifier)

# The ideal settler
@named settler = IdealClarifier()
push!(sys, settler)

# The settler pump
@named settler_pump = FlowPump(; flowrate=1000) # Flowrate 1000m^3/d
push!(sys, settler_pump)

# The wastage split into 2 flows
@named wastage_split = FlowSeparator()
push!(sys, wastage_split)

# The wastage pump
@named wastage_pump = FlowPump(; flowrate=250) # Flowrate 250 m^3/d wasted
push!(sys, wastage_pump)

# Connect the flows
append!(eqns, [
	# Influent -> unifier
	connect(outflows(influent, 1), inflows(unifier, 1)),
	# Unifier -> cstr
	connect(outflows(unifier, 1), inflows(cstr, 1)),
	# CSTR -> settler
	connect(outflows(cstr, 1), inflows(settler, 1)),
	# Settler -> settler pump
	connect(outflows(settler, :sludge), inflows(settler_pump, 1)), # Use the named access to get the sludge outflow of the settler, is more clear
	# Settler pump -> wastage split
	connect(outflows(settler_pump, 1), inflows(wastage_split, 1)),
	# wastage split back -> unifier
	connect(outflows(wastage_split, 1), inflows(unifier, 2)), # Note the 2 for the unifier: We do not want to connect the same inflow twice!
	# Wastage split -> wastage pump
	connect(outflows(wastage_split, 2), inflows(wastage_pump, 1)), # Also here the 2 for the split: We take the other outflow not going to the unifier!
])

# Assembly of the full system
@named model = ODESystem(eqns, t, systems = sys)
model_simplified = structural_simplify(model)

# Simulation and plotting
## Simulate the system (for 1 day)
# First create a ODEProblem and then solve it
prob = ODEProblem(model_simplified, [], (0, 1); initializealg=NoInit()) # No initialization scheme
sol = solve(prob)

## Plot the ouput, only the effluent of the settler
plot(
	sol,
	idxs=getproperty.(Ref(outflows(settler, :effluent)), propertynames(outflows(settler, :effluent))[2:end-1]), # settler effluent
	legend=:topright,
	title = "Activated Sludge System Effluent"
	)

Step-by-Step Explanations

Setup

First, the needed packages are imported and constants defined:

using BioChemicalTreatment # Reactors etc.
using ModelingToolkit # Modeling Framework
using DifferentialEquations # Solve the equations
using Plots # Plotting

Further, to allow for easy bookkeeping, empty vectors for the systems and all equations are created. Those are not neccessary to make here, but in the end for building the model those lists are required. Thus it is easier to keep these vectors growing as the systems and equations are defined which helps not to forget about some of them.

sys = ODESystem[]
eqns = Equation[]

System Creation

Then, all needed systems are created, including the "hidden" systems, which are usually not shown in a flow diagram. The second flow diagram above shows all of them, the "hidden" ones using a dashed line.

As usual, first the main process for the model is set as default, as this defines the tracked components in the flow, which can then be used automatically. In this example, the ASM1 model defined as matrix is used to show how one of the matrix-defined models can be loaded. This is done using:

@set_process asm1 = Process("ASM1")

Hereby, the default parameters are used. If others are desired, one could provide them directly as keyword arguments to this function or they could as well be overwritten later (for overriding the parameters, see this section).

Note

This way of loading matrix-defined models is equal for all included in the library.

Try it out by using ASM3 instead of ASM1 to work with the ASM3 instead. In this case, it is as well needed to adapt the influent to have values for the correct states.

Continuing, the reactor is created. As in the other examples, this one needs to know its volume and the states. It is then directly added to the system vector.

@named cstr = CSTR(1000; initial_states=[1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1]) # 1000 m^3 is the assumed volume, we set the initial states to all 1 as the equations do not work with the default (all 0)
push!(sys, cstr)
nameof.(sys) # Show all added systems

# output

1-element Vector{Symbol}:
 :cstr

Afterwards, all other systems can be created and added to the vector of systems. A constant influent is assumed (see the BSM1 example for a time varying influent) as well as some default values for the pumped flows.

# First the inflow
@named influent = Influent([7.0, 69.5, 6.95, 0.0, 31.56, 0.0, 0.0, 30.0, 202.32, 10.59, 0.0, 28.17, 0.0, 51.2], flowrate = 18446) # This is the stabilization influent from BSM1
push!(sys, influent)

# The unifier with two inflows
@named unifier = FlowUnifier() # defaults to two inflows, otherwise add more using `n_in`
push!(sys, unifier)

# The ideal settler
@named settler = IdealClarifier()
push!(sys, settler)

# The settler pump
@named settler_pump = FlowPump(; flowrate=1000) # Flowrate 1000m^3/d
push!(sys, settler_pump)

# The wastage split into 2 flows
@named wastage_split = FlowSeparator() # Defaults to two outflows
push!(sys, wastage_split)

# The wastage pump
@named wastage_pump = FlowPump(; flowrate=250) # Flowrate 250 m^3/d wasted
push!(sys, wastage_pump)
nameof.(sys) # Show all added systems

# output

7-element Vector{Symbol}:
 :cstr
 :influent
 :unifier
 :settler
 :settler_pump
 :wastage_split
 :wastage_pump

Now, all systems are created, and the flows can be connected (The other connections for the process to the reactor have already been added)

Connection of the Flows

Then, the connection equations for the systems are created and added to the vector.

append!(eqns, [
	# Influent to the unifier
	connect(outflows(influent, 1), inflows(unifier, 1)),
	# Unifier to the cstr
	connect(outflows(unifier, 1), inflows(cstr, 1)),
	# CSTR to the settler
	connect(outflows(cstr, 1), inflows(settler, 1)),
	# Settler to the settler pump
	connect(outflows(settler, :sludge), inflows(settler_pump, 1)), # Use the named access to get the sludge outflow of the settler, is more clear
	# Settler pump to wastage split
	connect(outflows(settler_pump, 1), inflows(wastage_split, 1)),
	# wastage split back to the unifier
	connect(outflows(wastage_split, 1), inflows(unifier, 2)), # Note the 2 for the unifier: We do not want to connect the same inflow twice!
	# Wastage split to wastage pump
	connect(outflows(wastage_split, 2), inflows(wastage_pump, 1)), # Also here the 2 for the split: We take the other outflow not going to the unifier!
])

The `inflows`and `outflows` functions

The inflows and outflows functions, when applied without any other argument, return a vector of the inflow and outflow connectors. However, also a second argument can be given:

Integer: Only the port at this position is returned. This is used e.g. for the first connection of the influent to the unifier
Symbol(starting with :): Only the port with the name according to the symbol is returned. This is used for the connection from the settler to the pump

Depending on the situation, one or the other is more convenient, but both can be used for every case. For some, the name is just not too helpful and then the numeric variant is easier to read, but for other systems the name might be helpful. For example, for getting the second outflow from the wastage separator, the following expressions are equivalent:

julia> outflows(wastage_split)[2]
Model wastage_split₊outflow_2:
Equations (15):
  15 connecting: see equations(expand_connections(wastage_split₊outflow_2))
Unknowns (15): see unknowns(wastage_split₊outflow_2)
  q(t): Flow Rate
  S_Alk(t): S_Alk concentration
  S_B(t): S_B concentration
  S_BN(t): S_BN concentration
  ⋮

julia> outflows(wastage_split, 2)
Model wastage_split₊outflow_2:
Equations (15):
  15 connecting: see equations(expand_connections(wastage_split₊outflow_2))
Unknowns (15): see unknowns(wastage_split₊outflow_2)
  q(t): Flow Rate
  S_Alk(t): S_Alk concentration
  S_B(t): S_B concentration
  S_BN(t): S_BN concentration
  ⋮

julia> outflows(wastage_split, :outflow_2)
Model wastage_split₊outflow_2:
Equations (15):
  15 connecting: see equations(expand_connections(wastage_split₊outflow_2))
Unknowns (15): see unknowns(wastage_split₊outflow_2)
  q(t): Flow Rate
  S_Alk(t): S_Alk concentration
  S_B(t): S_B concentration
  S_BN(t): S_BN concentration
  ⋮

Here, the named access seems overly verbose. However, for getting the sludge outflow from the settler, the same three variants are:

julia> outflows(settler)[2]
Model settler₊sludge:
Equations (15):
  15 connecting: see equations(expand_connections(settler₊sludge))
Unknowns (15): see unknowns(settler₊sludge)
  q(t): Flow Rate
  S_Alk(t): S_Alk concentration
  S_B(t): S_B concentration
  S_BN(t): S_BN concentration
  ⋮

julia> outflows(settler, 2)
Model settler₊sludge:
Equations (15):
  15 connecting: see equations(expand_connections(settler₊sludge))
Unknowns (15): see unknowns(settler₊sludge)
  q(t): Flow Rate
  S_Alk(t): S_Alk concentration
  S_B(t): S_B concentration
  S_BN(t): S_BN concentration
  ⋮

julia> outflows(settler, :sludge)
Model settler₊sludge:
Equations (15):
  15 connecting: see equations(expand_connections(settler₊sludge))
Unknowns (15): see unknowns(settler₊sludge)
  q(t): Flow Rate
  S_Alk(t): S_Alk concentration
  S_B(t): S_B concentration
  S_BN(t): S_BN concentration
  ⋮

Here instead, the variant with the name seems more clear than the number.

Now, with all systems and connection equations ready, the overall model can be build:

@named model = ODESystem(eqns, t, systems = sys)

And simplify it to get the minimum set of equations:

model_simplified = structural_simplify(model)

Simulation and Plotting

Finally, we simulate and plot the result.

## Simulate the system (for 1 day)
# First create a ODEProblem and then solve it
prob = ODEProblem(model_simplified, [], (0, 1); initializealg=NoInit()) # No initialization scheme for the DAE. See note below
sol = solve(prob)

## Plot the ouput, only the effluent of the settler
plot(
	sol,
	idxs=getproperty.(Ref(outflows(settler, :effluent)), propertynames(outflows(settler, :effluent))[2:end-1]),
	legend=:topright,
	title = "Activated Sludge System Effluent"
	)

About the `initializealg` argument to the problem

For the solution of the system, the algorithm generally needs initial conditions for all states. Due to the model transformations (the call to structural_simplify), the initial conditions needed might be different from the ones of the original problem. That's where the initialization algorithm jumps in and computes these newly needed initial conditions.

This initialization algorithm is set using the initializealg keyword argument when building the ODEProblem. Sometimes, this argument might not be neccessary (mostly for systems without any recirculation), but for others it will be needed. There are multiple options for this algorithm, see [https://docs.sciml.ai/DiffEqDocs/stable/solvers/dae_solve/#Initialization-Schemes] for reference.