Documentation Collimator Review by Iniyaa
csv files in the BSM2_Develop
The list of csv files available in the model can found by clicking the data source block diagram and the list of all csv files will appear under the properties tab under parameters-> file name
The csv zip folder downloaded from the sharepoint- It had 2 csv files which was identical with different names so I deleted one of them.
Deleted file: reac1_particulate(1).csv
CSV file input:
The Data columns in the input CSV files were not set properly. 'influent_SO_SNO_SNH_SND.csv'- the data columns was set to Temperature (incorrect)
Modified the data column from temperature to SO,SNO,SNH,SND
Similarly, 'influent_XND_SALK_TSS.csv' was modified.
ERRORS
1) len() of unsized object (Flow combiner_1->flowcombiner) For len() error, the code was modified from len() to ().size. For example: "num_flow = len(Q_array)" was changed to "num_flow = len(Q_array).size"
2) Jax error (ActivatedSludge->bioreactor1->Particulate->TSS)
3) PythonScriptError (flowcombiner_3): PythonScriptError in block flowcombiner_3.flowcombiner: Too many indices for array: 1 non-None/Ellipsis indices for dim 0.
Barbara's physical model
Barbara's physical model project was shared and duplicated to "MY PROJECTS" and named it as "BSM2_Physical_IT". The physical model under the project is "BSM2_Physical_Model" We made a copy of it, named it "BSM2_Physical_Model_modified", and made changes to "BSM2_Physical_Model_modified". 'BSM2_Physical_Model_IT_INIT.py' is the init script used in the "BSM2_Physical_Model_modified"
Q, T, Soluble and Particulate- all of them were combined into a single csv file nameD 'influents.csv' which was fed to the 'Q_bypass_PrimaryClarifier_split'
Flow combiner
Changes were made to: FlowCombiner, FlowCombiner_1, flowcombiner_2 and flowcombiner_3 Modified code for the flowcombiner which has 4 inputs (Depending on the inputs to the flow combiner, the value of len_Q and len_T changes):
# Inputs
# Q_array
# T_array
# S_array
# P_array
# Outputs
# Q
# T
# Soluble_conc
# Particulate_conc
import numpy as np
len_Q = 4
len_T = 4
S_types = 7
P_types = 7
Q= sum(Q_array)
Soluble_mass = np.zeros(S_types)
Soluble_conc = np.zeros(S_types)
Particulate_mass = np.zeros(P_types)
Particulate_conc = np.zeros(P_types)
T_total = 0.0
for i in range(len_Q):
T_total = T_total + (T_array[i]*Q_array[i])
T= T_total/Q
for j in range(S_types):
for i in range(len_Q):
Soluble_mass[j] = Soluble_mass[j] + S_array[S_types*i+j]*Q_array[i]
Soluble_conc[j]= Soluble_mass[j]/Q
for j in range(P_types):
for i in range(len_Q):
Particulate_mass[j] = Particulate_mass[j] + P_array[P_types*i+j]*Q_array[i]
Particulate_conc[j]= Particulate_mass[j]/Q
Activated Sludge
Modified the following: ActivatedSludge -> bioreactor_1 -> Particulate -> TSS
# import numpy as np
# import pandas as pd
# Particulate = in_0
# Particulate_out = Particulate
# # Read parameters from csv files
# df = pd.read_csv('Particulate2TSS.csv')
# ratio2TSS = df.iloc[0].to_numpy()
# Particulate2TSS = Particulate[0:5]
# Particulate_out[-1] = sum(ratio2TSS*Particulate2TSS)
# out_0 = Particulate_out
import pandas as pd
import numpy as np
Particulate = in_0
Particulate_out = Particulate
# Read parameters from csv files
df = pd.read_csv('Particulate2TSS.csv')
ratio2TSS = df.iloc[0,1]
#result= ratio2TSS * Particulate_out[1]
# Use the .at method to set the value in Particulate_out
#Particulate_out=Particulate_out.at[5].set(result)
out_0 = Particulate_out
The above code is incorrect refer to original BSM2 documentation.
Anaerobic Digester
Made changes to the ASM2ADM python script:
import numpy as np
import pandas as pd
import math
# List of Activated Sludge (AS) variables indexes
Si_i = 0
Ss_i = 1
So_i = 2
Sno_i = 3
Snh_i = 4
Snd_i = 5
Salk_i = 6
Xi_i = 0
Xs_i = 1
Xbh_i = 2
Xba_i = 3
Xp_i = 4
Xnd_i = 5
TSS_i = 6
# List of Anaerobic Digester (AD) variables indexes
Ssu_adm_i=0
Saa_adm_i=1
Sfa_adm_i=2
Sva_adm_i=3
Sbu_adm_i=4
Spro_adm_i=5
Sac_adm_i=6
Sh2_adm_i=7
Sch4_adm_i=8
Sic_adm_i=9
Sin_adm_i=10
Si_adm_i=11
Xxc_adm_i=12
Xch_adm_i=13
Xpr_adm_i=14
Xli_adm_i=15
Xsu_adm_i=16
Xaa_adm_i=17
Xfa_adm_i=18
Xc4_adm_i=19
Xpro_adm_i=20
Xac_adm_i=21
Xh2_adm_i=22
Xi_adm_i=23
scat_adm_i=24
san_adm_i=25
# Input
Q = in_0
T = in_1
Soluble = in_2
Particulate = in_3
pH_adm = in_4
# Create temporary variables
Soluble_temp1 = np.zeros((7)) # Soluble_temp1 = Soluble
Soluble_temp2 = np.zeros((7)) # Soluble_temp2 = Soluble
Particulate_temp1 = np.zeros((7)) # Particulate_temp1 = Particulate
Particulate_temp2 = np.zeros((7)) # Particulate_temp2 = Particulate
pH_temp1 = np.zeros((7)) # pH_temp1 = pH_adm
pH_temp2 = np.zeros((7)) # pH_temp2 = pH_adm
# Initialize output variables
y_output = np.zeros(26, dtype='double')
# Read parameters from csv files
df = pd.read_csv('INTERFACEPAR.csv')
CODequiv = df['CODequiv'][0]
fnaa = df['fnaa'][0] # Fraction of N in amino acids (Saa) and Xpr
fnxc = df['fnxc'][0] # N content of composite material based on BSM2
fnbac = df['fnbac'][0] # N content of biomass based on BSM1
fxni = df['fxni'][0] # N content of inerts XI and XP
fsni = df['fsni'][0] # N content of SI
fsni_adm = df['fsni_adm'][0] # N content of SI in the AD system
frlixs = df['frlixs'][0] # lipid fraction of non-nitrogenous Xs
frlibac = df['frlibac'][0] # lipid fraction of non-nitrogenous biomass
frxs_adm = df['frxs_adm'][0]
fdegrade_adm = df['fdegrade_adm'][0]
frxs_as = df['frxs_as'][0]
fdegrade_as = df['fdegrade_as'][0]
R = df['R'][0]
T_base = df['T_base'][0]
T_op = df['T_op'][0] # should be an input variable if dynamic temperature control is used
pK_w_base = df['pK_w_base'][0]
pK_a_va_base = df['pK_a_va_base'][0]
pK_a_bu_base = df['pK_a_bu_base'][0]
pK_a_pro_base = df['pK_a_pro_base'][0]
pK_a_ac_base = df['pK_a_ac_base'][0]
pK_a_co2_base = df['pK_a_co2_base'][0]
pK_a_IN_base = df['pK_a_IN_base'][0]
# Calculate parameters
factor = (1.0/T_base - 1.0/T_op)/(100.0*R)
pK_w = pK_w_base - math.log10(math.exp(55900.0*factor))
pK_a_co2 = pK_a_co2_base - math.log10(math.exp(7646.0*factor))
pK_a_IN = pK_a_IN_base - math.log10(math.exp(51965.0*factor))
alfa_va = 1.0/208.0*(-1.0/(1.0 + pow(10, pK_a_va_base - pH_adm)))
alfa_bu = 1.0/160.0*(-1.0/(1.0 + pow(10, pK_a_bu_base - pH_adm)))
alfa_pro = 1.0/112.0*(-1.0/(1.0 + pow(10, pK_a_pro_base - pH_adm)))
alfa_ac = 1.0/64.0*(-1.0/(1.0 + pow(10, pK_a_ac_base - pH_adm)))
alfa_co2 = -1.0/(1.0 + pow(10, pK_a_co2 - pH_adm))
alfa_IN = (pow(10, pK_a_IN - pH_adm))/(1.0 + pow(10, pK_a_IN - pH_adm))
alfa_NH = 1.0/14000.0 # convert mgN/l into kmoleN/m3 */
alfa_alk = -0.001 # convert moleHCO3/m3 into kmoleHCO3/m3 */
alfa_NO = -1.0/14000.0 # convert mgN/l into kmoleN/m3 */
"""
Step 1: The COD demand is caculated from the So and Sno states and
subtracted from the COD pool with an associated loss of Ss, Xs, Xbh,
and Xba (arrange in the order of depletion)
"""
# Let CODdemand be the COD demand of available electron
# acceptors prior to the anaerobic digester, i.e. oxygen and nitrate
CODdemand = Soluble[So_i] + CODequiv*Soluble[Sno_i] # So + CODequiv*Sno
Soluble_temp1[So_i] = 0.0 # So
Soluble_temp1[Sno_i] = 0.0 # Sno
"""
if extreme detail was used then some extra NH4 would be transformed
into N bound in biomass and some biomass would be formed when
removing the CODdemand (based on the yield). But on a total COD balance
approach the below is correct (neglecting the N need for biomass growth)
The COD is reduced in a hierarchical approach in the order:
1) Ss (Soluble[Ss_i]); 2) Xs (Particulate[Xs_i]); 3) Xbh(Particulate[Xbh_i]); 4) Xba (Particulate[Xba_i]).
It is no real improvement to remove SS and add biomass.
The net result is the same.
"""
CODdemand_removal_var = [Soluble[Ss_i], Particulate[Xs_i], Particulate[Xbh_i], Particulate[Xba_i]]
remainCOD = CODdemand
fulfilled_at_i_1 = len(CODdemand_removal_var)
for i in range(len(CODdemand_removal_var)):
# Check if COD demand can be fulfilled by Ss, Xx, Xbh, and Xba
if remainCOD > CODdemand_removal_var[i]:
remainCOD = remainCOD - CODdemand_removal_var[i]
# Break the loop if the COD demand is fulfilled
else:
fulfilled_at_i_1 = i
break
if fulfilled_at_i_1 == 0: # All COD Demand is reduced by use of Ss
Soluble_temp1[Ss_i] = Soluble[Ss_i] - remainCOD
elif fulfilled_at_i_1 == 1: # All COD demand is reduced by use of Ss and Xs
Soluble_temp1[Ss_i] = 0.0
Particulate_temp1[Xs_i] = Particulate[Xs_i] - remainCOD
elif fulfilled_at_i_1 == 2: # All COD demand is reduced by use of Ss, Xs, and Xbh
Soluble_temp1[Ss_i] = 0.0
Particulate_temp1[Xs_i] = 0.0
Particulate_temp1[Xbh_i] = Particulate[Xbh_i] - remainCOD
Soluble_temp1[Snh_i] = Soluble_temp1[Snh_i] + remainCOD*fnbac # related to Snh (ammonia and ammonium)
elif fulfilled_at_i_1 == 3: # All COD demand is reduced by use of Ss, Xs, Xbh, and Xba
Soluble_temp1[Ss_i] = 0.0
Particulate_temp1[Xs_i] = 0.0
Particulate_temp1[Xbh_i] = 0.0
Soluble_temp1[Snh_i] = Soluble_temp1[Snh_i] + Particulate[Xbh_i]*fnbac # when COD demand is reduced by Xbh
Particulate_temp1[Xba_i] = Particulate[Xba_i] - remainCOD
Soluble_temp1[Snh_i] = Soluble_temp1[Snh_i] + remainCOD*fnbac # when COD demand is reduced by Xba
else: # Not all COD demand is reduced by use of Ss, Xs, Xbh, and Xba
"""
If we reached here, there is problem with carbon shortage.
We will need to print an error message
"""
print("There is problem with carbon shortage. Occurs in Step 1 of asm2adm.py")
Soluble_temp1[Ss_i] = 0.0
Particulate_temp1[Xs_i] = 0.0
Particulate_temp1[Xbh_i] = 0.0
Soluble_temp1[Snh_i] += Particulate[Xbh_i]*fnbac
Particulate_temp1[Xba_i] = 0.0
Soluble_temp1[Snh_i] += Particulate[Xba_i]*fnbac
Soluble_temp1[So_i] = remainCOD # So
"""
Step 2: The Soluble organic nitrogen (Snd) is then transformed
along with a calculated amount of Ss (based on the amino acid nitrogen content)
to amino acids (Saa). Any remaining Ss is mapped to Ssu and the remaining Snd is stored
in a temporary variable. Enough Snd must be available for mapping to Saa.
"""
# Snd (organic nitrogen associated with Ss (readily biodegradable substrate))
sorgn = Soluble[Snd_i]/fnaa # Saa COD equivalent to Snd
# Mapping Snd-N in terms of COD to Saa
# If not all Snd-N in termos of COD fits into amino acids (Saa)
if (sorgn >= Soluble_temp1[Ss_i]):
y_output[Saa_adm_i] = Soluble_temp1[Ss_i] # map all Ss COD into Saa
Soluble_temp1[Ss_i] = 0.0 # All Ss is used
Soluble_temp1[Snd_i] += - Soluble_temp1[Ss_i]*fnaa # excess Snd
# If all Snd-N fits into amino acids (Saa)
else:
y_output[Saa_adm_i] = sorgn # map all Snd related COD into Saa
Soluble_temp1[Ss_i] += - sorgn # excess Ss will become sugar in ADM1 (no nitrogen association)
Soluble_temp1[Snd_i] = 0.0 # All Snd is used
"""
Step 3: The remaining Xnd (Particulate[Xnd_i]) and Xs (Particulate[Xs_i]) are mapped to Xpr (proteins). Any Xs left is further
mapped to Xch and Xli (no N contents) according to the lipid content of Xs. The remaining Xnd
is stored in a temporary variable. Enough Xnd must be available for mapping to Xpr.
"""
xorgn = Particulate[Xnd_i]/fnaa # Xpr COD equivalent to Xnd (Particulate[Xnd_i])
if (xorgn >= Particulate_temp1[Xs_i]):
xprtemp = Particulate_temp1[Xs_i] # map all Xs COD into Xpr from step 3
Particulate_temp1[Xs_i] = 0.0 # All Xs is used
xlitemp = 0.0 # no Xs left to map to Xli
xchtemp = 0.0 # no Xs left to map to Xch
Particulate_temp1[Xnd_i] += - Particulate_temp1[Xs_i]*fnaa # excess Xnd
else:
xprtemp = xorgn # map all Xnd related COD to Xpr from step 3
Particulate_temp1[Xnd_i] = 0.0 # All Xnd is used
xlitemp = frlixs*(Particulate_temp1[Xs_i]-xorgn) # Part of Xnd map to Xli
xchtemp = (1- frlixs)*(Particulate_temp1[Xs_i]-xorgn) # Part of Xnd map to Xch
Particulate_temp1[Xs_i] = 0.0 # All Xs is used to map to Xpr, Xli, and Xch
"""
Step 4: Xbh (Particulate[Xbh_i]) and Xba (Particulate[Xba_i]) are first split down in a degradable and a non-degredable part. The degradable fraction is
partially mapped to Xpr (taking into account its nitrogen content and using the remaining Xnd). Any
remnants are mapped to Xch and Xli according to the lipid content of biomass (Xbh and Xba). The remaining Xnd
is stored in a temporary variable.
"""
biomass = Particulate_temp1[Xbh_i] + Particulate_temp1[Xba_i] # Sum of Xbh ([4]) and Xba ([5])
# Split to degradable and non-degredable
biomass_nobio = biomass*(1.0 - frxs_adm)
biomass_bioN = ((biomass*fnbac) - (biomass_nobio*fxni)) # biomass associate with N
biomass_bioCOD = biomass - biomass_nobio
if (biomass_bioN < 0.0):
print('ERROR: not enough biomass N to map the requested inert part')
# Problem: if here we should print 'ERROR: not enough biomass N to map the requested inert part'
if ((biomass_bioN/fnaa) <= biomass_bioCOD):# all biomass N is used
xprtemp2 = biomass_bioN/fnaa
remainCOD = biomass_bioCOD - xprtemp2
# part of remaining Xnd-N to form proteins
if ((Particulate_temp1[Xnd_i]/fnaa) > remainCOD):
xprtemp2 += remainCOD # Xpr form from step 4
Particulate_temp1[Xnd_i] += - remainCOD*fnaa # remaining Xnd
remainCOD = 0.0
# remaining Xnd-N to form protens
else:
xprtemp2 += (Particulate_temp1[Xnd_i]/fnaa)
remainCOD += - (Particulate_temp1[Xnd_i]/fnaa) # remaining COD
Particulate_temp1[Xnd_i] = 0.0
# Remaining COD map to Xli and Xch
xlitemp2 = frlibac*remainCOD
xchtemp2 = (1-frlibac)*remainCOD
else: #all biomass COD is used
xprtemp2 = biomass_bioCOD
Particulate_temp1[Xnd_i] += biomass_bioN - (xprtemp2*fnaa) # remaining N in Xnd
# All Xbh and Xba are used
Particulate_temp1[Xbh_i] = 0.0 # Xbh
Particulate_temp1[Xba_i] = 0.0 # Xba
"""
Step 5a: Soluble inerts (Si) are mapped directly across the interface accounting
for the difference in nitrogen content of Si in ASM1 and ADM1 respectively
(using all remaining Snd, Xnd, and Snh).
N content may be different, take first from Si-N, then Snd-N, then Xnd-N, then Snh.
"""
inertS = 0.0 # Si in ADM1
N_content_variable_5a = [Soluble_temp1[Snd_i], Particulate_temp1[Xnd_i], Soluble_temp1[Snh_i]]
fulfilled_at_i_5a = len( N_content_variable_5a)
if (fsni < fsni_adm): # if N-content in Si in ASM < N-content in Si in ADM
inertS += Soluble_temp1[Si_i]*fsni/fsni_adm
Soluble_temp1[Si_i] += -(Soluble_temp1[Si_i]*fsni/fsni_adm) # remaining Si in ASM that needed to map to ADM
# move to next N-content variables
for i in range(fulfilled_at_i_5a):
if (N_content_variable_5a[i] < Soluble_temp1[Si_i]*fsni_adm): # Check if N-content in each N-variable is enough before moving to the next
inertS += N_content_variable_5a[i]/fsni_adm
Soluble_temp1[Si_i] += - (N_content_variable_5a[i]/fsni_adm)
else: # break the loop if the N-content variable is enough
fulfilled_at_i_5a = i
inertS += Soluble_temp1[Si_i]
break
if fulfilled_at_i_5a == 0: # N in Snd is enough for mapping
Soluble_temp1[Snd_i] += -(Soluble_temp1[Si_i]*fsni_adm)
elif fulfilled_at_i_5a == 1: # N in Xnd is enough for mapping
Particulate_temp1[Xnd_i] += -(Soluble_temp1[Si_i]*fsni_adm)
Soluble_temp1[Snd_i] = 0.0
elif fulfilled_at_i_5a == 2: # N in Snh is enough for mapping
Soluble_temp1[Snh_i] += -(Soluble_temp1[Si_i]*fsni_adm)
Soluble_temp1[Snd_i] = 0.0
Particulate_temp1[Xnd_i] = 0.0
else: # N in Snd, Xnd, and Snh are not enough for mapping
"""
Should print warning: Nitrogen shortage when converting Si
Put remaining Si as monosacharides
"""
Soluble_temp1[Ss_i] += Soluble_temp1[Si_i] # add remaining Si to monosacharides
Soluble_temp1[Snd_i] = 0.0
Particulate_temp1[Xnd_i] = 0.0
Soluble_temp1[Snh_i] = 0.0
else: # N in Si in ASM is enough for mapping
inertS += Soluble_temp1[Si_i]
Soluble_temp1[Snd_i] = Soluble_temp1[Snd_i] + Soluble_temp1[Si_i]*(fsni-fsni_adm) # put remaining N in Si (AS) in Snd
Soluble_temp1[Si_i] = 0.0
"""
Step 5b: Particulate inerts (Xi) (Particulate[0]) and the anaerobically
non-degradable biomass fraction (Xp) (Particulate[4]) are directly mapped to ADM1 Xi
if fdegrade_adm = 0
assumption: they have the same N-Content in ASM1 and ADM1 particulate inerts
"""
inertX = (1-fdegrade_adm)*(Particulate_temp1[Xi_i] + Particulate_temp1[Xp_i]) # ADM1 Xi
"""
if part of Xi and Xp in the ASM can be degraded in the AD (fdegrade_adm > 0)
- When bsm2 was created, they have no knowledge about the contents so they put it
in a composits (Xc)
- They need to keep track of the associated nitrogen
- N content which may be different, take first from Xi and Xp-N, then Xnd-N, then Snd-N, then Snh.
- A similar principle could be used for other states.
"""
xc = 0.0 # initialize composites
xlitemp3 = 0.0
xchtemp3 = 0.0
# Mapping noninertX to Xc
if (fdegrade_adm>0):
# degradable Xi
noninertX = fdegrade_adm*(Particulate_temp1[Xi_i] + Particulate_temp1[Xp_i])
# temporary variables Xnd, Snd, and Snh
N_content_variable_5b = [Particulate_temp1[Xnd_i], Soluble_temp1[Snd_i], Soluble_temp1[Snh_i]]
if (fxni < fnxc): #N-content in Xi and Xp < N-content in Xc
xc = noninertX*fxni/fnxc # convert N- content biodegradable Xi and Xp to Xc
noninertX += -xc # remaining degradable Xi
# move to Xnd-N, then Snd-N, then Snh
fulfilled_at_i_5b = len( N_content_variable_5b)
for i in range(fulfilled_at_i_5b):
if (N_content_variable_5b[i] < noninertX*fnxc): # Check if N-content in each N-variable is enough before moving to the next
xc += N_content_variable_5b[i]/fnxc
noninertX += -N_content_variable_5b[i]/fnxc
else: # break the loop if the N-content variable is enough
fulfilled_at_i_5b = i
xc += noninertX
break
if fulfilled_at_i_5b == 0: # N in Xnd is enough for mapping
Particulate_temp1[Xnd_i] += -(noninertX*fnxc)
noninertX = 0.0
elif fulfilled_at_i_5b == 1: # N in Snd is enough for mapping
Soluble_temp1[Snd_i] += -(noninertX*fnxc)
Particulate_temp1[Xnd_i] = 0.0
noninertX = 0.0
elif fulfilled_at_i_5b == 2: # N in Snh is enough for mapping
Soluble_temp1[Snh_i] += -(noninertX*fnxc)
Particulate_temp1[Xnd_i] = 0.0
Soluble_temp1[Snd_i] = 0.0
noninertX = 0.0
else: # N in Xnd, Snd, Snh are not enough for mapping
"""
Should be a WARNING printout: Nitrogen shortage when converting biodegradable XI&XP
Putting remaining XI&XP as lipids (50%) and carbohydrates (50%)"""
xlitemp3 = 0.5*noninertX
xchtemp3 = 0.5*noninertX
noninertX = 0.0
else: # N-content in Xi and Xp is enough for mapping
xc += noninertX
Particulate_temp1[Xnd_i] += noninertX*(fxni-fnxc) # put the remaining N in as Xnd
noninertX = 0.0
Particulate_temp1[Xi_i] = 0.0
Particulate_temp1[Xp_i] = 0.0
"""
Step 6: All remaining nitrogen compounds (Snd, Xnd, Snh) are directly mapped to
the inorganic nitrogen fraction Sin (y[10])
"""
Sin = (Soluble_temp1[Snd_i] + Particulate_temp1[Xnd_i] + Soluble_temp1[Snh_i])/14000.0 # convert from g N/m3 to kmole N/m3
"""
Final mapping of variables and assures the maintenance of the charge balance.
It uses all charged components at both sides of the interface (also accounting for dissociation
state with respect to the respective pH and T values) and determines the SIC, San and Scat in order to
close the balance
"""
Q_out = Q
T_out = T_op - 273.15 # T_op can be dynamic input. For this case, it was set to 35 degC
y_output[Ssu_adm_i] = Soluble_temp1[Ss_i]/1000.0 # Ssu = monosacharides (kg COD/m3)
y_output[Saa_adm_i] = y_output[Saa_adm_i]/1000.0 # Saa = amino acids (kg COD/m3)
y_output[Sin_adm_i] = Sin # Sin = inorganic nitrogen (kmole N/m3)
y_output[Si_adm_i] = inertS/1000.0 # Si = Soluble inerts (kg COD/m3)
y_output[Xxc_adm_i] = xc/1000.0 # Xc = composites (kg COD/m3)
y_output[Xch_adm_i] = (xchtemp + xchtemp2 + xchtemp3)/1000.0 # Xch = carbohydrates (kg COD/m3)
y_output[Xpr_adm_i] = (xprtemp + xprtemp2)/1000.0 # Xpr = proteins (kg COD/m3)
y_output[Xli_adm_i] = (xlitemp + xlitemp2 + xlitemp3)/1000.0 # Xli = lipids (kg COD/m3)
y_output[Xi_adm_i] = (biomass_nobio + inertX)/1000.0 # Xi = particulate inerts (kg COD/m3)
# Charge balance, output S_IC
y_output[Sic_adm_i] = ((Soluble_temp2[Sno_i]*alfa_NO + Soluble_temp2[Snh_i]*alfa_NH + Soluble_temp2[Salk_i]*alfa_alk)\
- (y_output[Sva_adm_i]*alfa_va+y_output[Sbu_adm_i]*alfa_bu+y_output[Spro_adm_i]*alfa_pro+y_output[Sac_adm_i]*alfa_ac+y_output[Sin_adm_i]*alfa_IN))\
/alfa_co2
# Calculate anions and cations based on full charge balance including H+ and OH-
ScatminusSan = y_output[Sva_adm_i]*alfa_va + y_output[Sbu_adm_i]*alfa_bu + y_output[Spro_adm_i]*alfa_pro + y_output[Sac_adm_i]*alfa_ac\
+ y_output[Sin_adm_i]*alfa_IN + y_output[Sic_adm_i]*alfa_co2 + pow(10, (-pK_w + pH_adm)) - pow(10, -pH_adm)
if (ScatminusSan > 0):
y_output[scat_adm_i] = ScatminusSan
y_output[san_adm_i] = 0.0
else:
y_output[scat_adm_i] = 0.0
y_output[san_adm_i] = -1.0*ScatminusSan
# Input-output mass balance check for COD and N
COD_input_ASM = 0.0 # need to calculate later
COD_output_ADM = (y_output[Ssu_adm_i] + y_output[Saa_adm_i] + y_output[Si_adm_i] + y_output[Xxc_adm_i] + y_output[Xch_adm_i] + y_output[Xpr_adm_i]\
+ y_output[Xli_adm_i] + y_output[Xi_adm_i])*1000
N_input_ASM = 0.0 # need to calculate later
ixn = 0.08 # based on the paper but need to recheck with matlab
N_output_ASM = (y_output[Ssu_adm_i]+y_output[Saa_adm_i])*(fnaa/14.0) + y_output[Sin_adm_i] + (y_output[Si_adm_i]+y_output[Xi_adm_i])*ixn # need to check the unit
# Output
y_Soluble = y_output[0:12]
y_Particulate = y_output[12:24]
y_ions = y_output[24:26]
out_0 = Q_out
out_1 = T_out
out_2 = y_Soluble
out_3 = y_Particulate
out_4 = y_ions
BSM2 model
Downloaded the BSM2 files from github: https://github.com/wwtmodels/Benchmark-Simulation-Models
Comparison of BSM2_Physical_Model_modified with BSM2 [up to ActivatedSludge]
BSM2 [up to ActivatedSludge] model was sent by the Collimator support team. This model was compared with the BSM2_Physical_Model_modified model under BSM2_Physical_IT project
BSM2 [up to ActivatedSludge] model
BSM2 [up to ActivatedSludge] model only has BSM2 upto Activated Sludge (with only one bioreactor in it) :
BSM2 [up to ActivatedSludge]_ ActivatedSludge_flowcombiner: The python script uses numpy array
The input and output of the bioreactor 1 is as same as the BSM2_Physical_Model_modified model.
BSM2 [up to ActivatedSludge]_ ActivatedSludge_flowcombiner_bioreactor_1:
BSM2 [up to ActivatedSludge]_ ActivatedSludge_flowcombiner_bioreactor_1_carbon_flowcombiner:
The carbon_flowcombiner python script was modified
carbonflow = in_0 #m^3/day
carbonconc = in_1 #mg COD/l (gCOD/m^3)
flow = in_2 #m^3/day
temp = in_3 #degree celcius
soluble_var = in_4 #gCOD/m^3
particulate_var = in_5 #gCOD/m^3
comb_flow = carbonflow + flow
out_soluble = np.array(soluble_var*flow/comb_flow)
out_soluble[1] = (soluble_var[1]*flow+ carbonconc*carbonflow)/comb_flow
out_particulate = particulate_var*flow/comb_flow
# outputs
out_0 = comb_flow
out_1 = temp
out_2 = out_soluble
out_3 = out_particulate
BSM2 [up to ActivatedSludge]_ ActivatedSludge_flowcombiner_bioreactor_1_Soluble:
The ODE mass balance block, product and the adder block is replaced with a python script and an intregrator block without needing to mux and demux it.
BSM2 [up to ActivatedSludge]_ ActivatedSludge_flowcombiner_bioreactor_1_Soluble_dot_state:
#init
import numpy as np
idx_So = 2
#end of init
# set reaction term for So
reac_So = (So_sat_temp - state[idx_So])*KLa_temp + reac_in[idx_So]
reac_in_mod = reac_in.copy()
reac_in_mod[idx_So] = reac_So
# compute state derivatives
dot_state = (Soluble_in - state)*Q_per_VOL + reac_in_mod
BSM2 [up to ActivatedSludge]_ ActivatedSludge_flowcombiner_bioreactor_1_Soluble_x_gt_or_equal_zero:
#init
import numpy as np
# List of indexes
So_i = 2
#end of init
state_sol = np.where(Soluble_in<0.0, 0.0, Soluble_in)
if Kla<0.0:
state_sol[So_i] = np.abs(Kla)
else:
# this else block can set state_sol[So_i] back to negative and therefore shouldn't be included (logic from original model)
state_sol[So_i] = Soluble_in[So_i]
# Outputs
state_So_i = state_sol[So_i]
Similarly for Particulates, ODE mass balance blocks were replaced with a python script and an Integrator.
BSM2 [up to ActivatedSludge]_ActivatedSludge_flowcombiner_bioreactor_1_Particulate:
BSM2 [up to ActivatedSludge]_ActivatedSludge_flowcombiner_bioreactor_1_Particulate_TSS_0: The python script was modified.
BSM2_Physical_Model_modified model
BSM2_Physical_Model_modified_ActivatedSludge_flowcombiner: The python script uses jax array
BSM2_Physical_Model_modified_ActivatedSludge_flowcombiner_bioreactor_1:
BSM2_Physical_Model_modified_ActivatedSludge_flowcombiner_bioreactor_1_carbon_flowcombiner:
#init
import numpy as np
# List of indexes
Si_i = 0
Ss_i = 1
So_i = 2
Sno_i = 3
Snh_i = 4
Snd_i = 5
Salk_i = 6
Xi_i = 0
Xs_i = 1
Xbh_i = 2
Xba_i = 3
Xp_i = 4
Xnd_i = 5
TSS_i = 6
#end of init
import numpy as np
carbonflow = in_0 #m^3/day
carbonconc = in_1 #mg COD/l (gCOD/m^3)
flow = in_2 #m^3/day
temp = in_3 #degree celcius
soluble_var = in_4 #gCOD/m^3
particulate_var = in_5 #gCOD/m^3
comb_flow = carbonflow + flow
out_soluble = np.empty(len(soluble_var),dtype='double')
out_particulate=np.empty(len(particulate_var),dtype='double')
for i in range(len(soluble_var)):
if i == 1:
out_soluble[i]= (soluble_var[i]*flow+ carbonconc*carbonflow)/comb_flow
else:
out_soluble[i]= (soluble_var[i]*flow)/comb_flow
for j in range(len(particulate_var)):
out_particulate[j] = (particulate_var[j]*flow)/comb_flow
# print("Carbon flowcombiner: in_4[4]", in_4[4], "out_soluble[4]" , out_soluble[4])
# outputs
out_0 = comb_flow
out_1 = temp
out_2 = out_soluble
out_3 = out_particulate
BSM2_Physical_Model_modified_ActivatedSludge_flowcombiner_bioreactor_1_Soluble:
BSM2_Physical_Model_modified_ActivatedSludge_flowcombiner_bioreactor_1_Soluble_x_gt_or_equal_zero:
Soluble_in = in_0
Kla = in_1
state_sol = np.empty(Soluble_in.size, dtype='double')
for i in range(len(Soluble_in)):
if i == So_i:
if Kla < 0:
state_sol[i] = abs(Kla)
else:
state_sol[i] = Soluble_in[i]
else:
if in_0[i] < 0:
state_sol[i] = 0.0
else:
state_sol[i] = Soluble_in[i]
# Outputs
out_0 = state_sol
out_1 = state_sol[So_i]
BSM2_Physical_Model_modified_ActivatedSludge_flowcombiner_bioreactor_1_Particulate:
