WRRF/NSS Merged Model

[wobrien3]

Describe the problem Trying to add elements of the Bwaise system into a WRRF model but getting errors with implementation (notably, Excretion & PitLatrine).

To Reproduce Steps to reproduce the behavior, e.g.:

1. Run code below (change su_data_path to your computer's data path)

Screenshots

Environment (please complete the following information):

1. Versions of the following packages:

   • EXPOsan 1.3.2
   • QSDsan 1.3.1
   • BioSTEAM 2.37.6
   • Thermosteam 0.35.2

2. Operating systems: macOS 11.6

Additional context

'''
QSDsan: Quantitative Sustainable Design for sanitation and resource recovery systems

This module is developed by:

    Saumitra Rai <raisaumitra9@gmail.com>

    Joy Zhang <joycheung1994@gmail.com>

This module is under the University of Illinois/NCSA Open Source License.
Please refer to https://github.com/QSD-Group/QSDsan/blob/main/LICENSE.txt
for license details.
'''
import os, pickle, numpy as np, qsdsan as qs
from qsdsan import (
    processes as pc, 
    sanunits as su, 
    WasteStream, 
    System, 
    currency, 
    ImpactItem, 
    PowerUtility, 
    Model, 
    Metric
    )
from qsdsan.utils import (
    AttrFuncSetter,
    AttrSetter,
    DictAttrSetter,
    data_path,
    dct_from_str,
    load_data,
    ospath,
    time_printer,
    get_SRT
    )
from exposan.bsm1 import data_path

#Bwaise imports
from chaospy import distributions as shape
from thermosteam.functional import V_to_rho, rho_to_V 
from exposan.utils import batch_setting_unit_params, add_fugitive_items, run_uncertainty as run
from exposan import bwaise as bw
from exposan.bwaise import (
    _load_components,
    create_system,
    get_alt_salary,
    get_biogas_factor,
    get_decay_k,
    get_LCA_metrics,
    get_TEA_metrics,
    get_recoveries,
    GWP_dct,
    price_dct,
    results_path,
    )

# %%

# =============================================================================
# Parameters and util functions
# =============================================================================

Q = 60000      # influent flowrate [m3/d]
Temp = 273.15+20 # temperature [K]
V_an_1 = 2500   # anoxic zone tank volume [m3] (HRT = 2/2 HRS, M&E)
V_an_2 = 2500 # anoxic zone tank volume [m3] (HRT = 2/2 HRS, M&E)
V_an_3 = 7500 # anoxic zone tank volume [m3] (HRT = 3 HRS, M&E)
V_ae_1 = 20000 # aerated zone tank volume [m3]  (HRT = 8 HRS, M&E)
V_ae_2 = 1875  # aerated zone tank volume [m3]  (HRT = 0.75 HRS, M&E)
# V_ae_2 = 2521  
Q_was = 1700  # Based on trial and error to get SRT between 10-20 days [m3/day]
Q_ras = 45000   # recycle sludge flowrate (75% of influent, M&E) [m3/day] 
biomass_IDs = ('X_H', 'X_AUT')
ammonia_ID = ('S_NH4',)

# aer = pc.DiffusedAeration('Fixed_Aeration', 'S_O', KLa_20=240, SOTE=0.3, V=V_ae,
#                           T_air=Temp, T_water=Temp, d_submergence=4-0.3)
# aer.A = 8.10765
# aer.B = 1750.286
# aer.C = 235.0
# =============================================================================

beijing_inf_kwargs = {
    'concentrations': {
        'S_F': 152,
        'S_A': 15,
        'S_I': 30,
        'X_S': 25,
        'X_I': 0,
        'X_AUT': 0, 
        'X_PAO': 0, 
        'X_H': 0,
        'X_PHA': 10,
        'X_PP': 0,
        'S_NH4': 38,
        'S_PO4': 2.8,
        'S_N2': 0,
        'S_NO3': 0, 
        'S_ALK':7*12,
        },
    'units': ('m3/d', 'mg/L'),
    }

default_asm2d_kwargs = dict(iN_SI=0.01, iN_SF=0.03, iN_XI=0.02, iN_XS=0.04, iN_BM=0.07,
            iP_SI=0.0, iP_SF=0.01, iP_XI=0.01, iP_XS=0.01, iP_BM=0.02,
            iTSS_XI=0.75, iTSS_XS=0.75, iTSS_BM=0.9,
            f_SI=0.0, Y_H=0.625, f_XI_H=0.1,
            Y_PAO=0.625, Y_PO4=0.4, Y_PHA=0.2, f_XI_PAO=0.1,
            Y_A=0.24, f_XI_AUT=0.1,
            K_h=3.0, eta_NO3=0.6, eta_fe=0.4, K_O2=0.2, K_NO3=0.5, K_X=0.1,
            mu_H=6.0, q_fe=3.0, eta_NO3_H=0.8, b_H=0.4, K_O2_H=0.2, K_F=4.0,
            K_fe=4.0, K_A_H=4.0, K_NO3_H=0.5, K_NH4_H=0.05, K_P_H=0.01, K_ALK_H=0.1,
            q_PHA=3.0, q_PP=1.5, mu_PAO=1.0, eta_NO3_PAO=0.6, b_PAO=0.2, b_PP=0.2,
            b_PHA=0.2, K_O2_PAO=0.2, K_NO3_PAO=0.5, K_A_PAO=4.0, K_NH4_PAO=0.05,
            K_PS=0.2, K_P_PAO=0.01, K_ALK_PAO=0.1,
            K_PP=0.01, K_MAX=0.34, K_IPP=0.02, K_PHA=0.01,
            mu_AUT=1.0, b_AUT=0.15, K_O2_AUT=0.5, K_NH4_AUT=1.0, K_ALK_AUT=0.5, K_P_AUT=0.01,
            k_PRE=1.0, k_RED=0.6, K_ALK_PRE=0.5, path=os.path.join(data_path, '_asm2d.tsv'),
            )

default_init_conds = {
        'S_F':5,
        'S_A':2,
        'X_I':1000,
        'X_S':100,
        'X_H':500,
        'X_AUT':100,
        #'X_P':100,
        'S_O2':2,
        'S_NO3':20,
        'S_NH4':2,
        'S_ALK':7*12,
    }

def batch_init(sys, path, sheet):
    df = load_data(path, sheet)
    dct = df.to_dict('index')
    u = sys.flowsheet.unit # unit registry
    for k in sys.units:
        if k.ID.startswith('O'):
            k.set_init_conc(**dct['O'])
        elif k.ID.startswith('A'):
            k.set_init_conc(**dct['A'])           
    c1s = {k:v for k,v in dct['C1_s'].items() if v>0}
    c1x = {k:v for k,v in dct['C1_x'].items() if v>0}
    tss = [v for v in dct['C1_tss'].values() if v>0]
    u.C1.set_init_solubles(**c1s)
    u.C1.set_init_sludge_solids(**c1x)
    u.C1.set_init_TSS(tss)

#%%

# =============================================================================
# Data Sheets
# =============================================================================

join_path = lambda prefix, file_name: os.path.join(prefix, file_name)

su_data_path = '/Users/willobrien/QSDsan/qsdsan/data/sanunit_data'
#su_data_path = ospath.join(data_path, 'sanunit_data')
excretion_data = load_data(join_path(su_data_path, '_excretion.tsv'))
toilet_data = load_data(join_path(su_data_path, '_toilet.tsv'))
pit_latrine_data = load_data(join_path(su_data_path, '_pit_latrine.tsv'))
uddt_data = load_data(join_path(su_data_path, '_uddt.tsv'))
drying_bed_data = load_data(join_path(su_data_path, '_drying_bed.tsv'))
liquid_bed_data = load_data(join_path(su_data_path, '_liquid_treatment_bed.tsv'))
sedimentation_tank_data = load_data(join_path(su_data_path, '_sedimentation_tank.tsv'))
anaerobic_lagoon_data = load_data(join_path(su_data_path, '_anaerobic_lagoon.tsv'))
facultative_lagoon_data = load_data(join_path(su_data_path, '_facultative_lagoon.tsv'))
sludge_separator_data = load_data(join_path(su_data_path, '_sludge_separator.tsv'))
abr_data = load_data(join_path(su_data_path, '_anaerobic_baffled_reactor.tsv'))

#%%

# =============================================================================
# Benchmark Simulation Model No. 1
# =============================================================================

#Changed function name to ensure no issues with _components.py
def create_components():
     asm2d_cmps = pc.create_asm2d_cmps()
     asm2d_cmps.X_S.f_BOD5_COD = 0.54
     # CO2 = qs.Component.from_chemical('S_CO2', search_ID='CO2', particle_size='Dissolved gas', degradability='Undegradable', organic=False)
     # CH4 = qs.Component.from_chemical('S_CH4', search_ID='CH4', particle_size='Dissolved gas', degradability='Undegradable', organic=False)
     # H2 = qs.Component.from_chemical('S_H2', search_ID='H2', particle_size='Dissolved gas', degradability='Undegradable', organic=False)
     # cmps1 = qs.Components.load_default()
     # ash = cmps1.X_Ig_ISS.copy('ash')
     cmps = qs.Components([*asm2d_cmps, 
                           # CO2, CH4, H2, ash
                           ])
     cmps.compile()
     return cmps

def create_system(flowsheet=None, inf_kwargs={}, asm_kwargs={}, init_conds={},
                  aeration_processes=()):
    flowsheet = flowsheet or qs.Flowsheet('bsm1')
    qs.main_flowsheet.set_flowsheet(flowsheet)

    # Components and stream
    cmps = create_components()
    qs.set_thermo(cmps)
    wastewater = WasteStream('wastewater', T=Temp)
    #inf_kwargs = inf_kwargs or default_inf_kwargs
    inf_kwargs = beijing_inf_kwargs 
    wastewater.set_flow_by_concentration(Q, **inf_kwargs)

    effluent = WasteStream('effluent', T=Temp)
    WAS = WasteStream('WAS', T=Temp)
    RAS = WasteStream('RAS', T=Temp)
    flowthrough = WasteStream('flowthrough', T=Temp)
    recycle = WasteStream('recycle', T=Temp)

    #Add Bwaise System Elements
    B1 = su.Excretion('B1', outs=('urine', 'feces'))

    B2 = su.PitLatrine('B2', ins=(B1-0, B1-1,
                                  'toilet_paper', 'flushing_water',
                                  'cleansing_water', 'desiccant'),
                       outs=('mixed_waste', 'leachate', 'A2_CH4', 'A2_N2O'),
                       OPEX_over_CAPEX=0.05,
                       decay_k_COD=3,
                       decay_k_N=3,
                       max_CH4_emission=0.25
                       )
    B3 = su.Trucking('B3', ins=B2-0, outs=('transported', 'conveyance_loss'),
                     load_type='mass', distance=5, distance_unit='km',
                     interval=0.8, interval_unit='yr',
                     loss_ratio=0.02)

    B4 = su.Sedimentation('B4', ins=B3-0,
                          outs=('liq', 'sol', 'A5_CH4', 'A5_N2O'),
                          decay_k_COD=3,
                          decay_k_N=3,
                          max_CH4_emission=0.25)
    B5 = su.DryingBed('B5', ins=B4-1, outs=('dried_sludge', 'evaporated',
                                            'DB1_CH4', 'DB1_N2O'),
                      design_type='unplanted',
                      decay_k_COD=3,
                      decay_k_N=3,
                      max_CH4_emission=0.25)

    M1 = su.Mixer('M1', ins=(B2-2, B4-2, B5-2), outs='CH4')
    #A7 = su.Mixer('A7', ins=(A2-1, A4-2, A5-2), outs=streamA.CH4)
    M1.add_specification(lambda: add_fugitive_items(M1, 'CH4_item'))
    M1.line = 'fugitive CH4 mixer'

    M2 = su.Mixer('M2', ins=(B2-3, B4-3, B5-3), outs='N2O')
    M2.add_specification(lambda: add_fugitive_items(M2, 'N2O_item'))
    M2.line = 'fugitive N2O mixer'

    # Process models
    # if aeration_processes:
    #     aer1, aer2, aer3 = aeration_processes
    # else:
    #     aer1 = aer2 = pc.DiffusedAeration('aer1', 'S_O', KLa=240, DOsat=8.0, V=V_ae)
    #     aer3 = pc.DiffusedAeration('aer3', 'S_O', KLa=84, DOsat=8.0, V=V_ae)
    # # asm_kwargs = asm_kwargs or default_asm_kwargs
    asm2d = pc.ASM2d(iP_SF=0.005, iP_XS=0.005, iP_XI=0.005, iN_BM=0.1, iTSS_XI=0.72)

    A1 = su.CSTR('A1', [wastewater, RAS], V_max=V_an_1,
                  aeration=None, suspended_growth_model=asm2d)

    A2 = su.CSTR('A2', [recycle, A1-0], V_max=V_an_2,
                  aeration=None, suspended_growth_model=asm2d)

    O1 = su.CSTR('O1', A2-0, V_max=V_ae_1, aeration= 2, DO_ID='S_O2', 
                 suspended_growth_model=asm2d)

    S1 = qs.sanunits.Splitter('S1', ins=O1-0, outs=(flowthrough, recycle),
                          split=0.3684, init_with='WasteStream')

    A3 = su.CSTR('A3', flowthrough, V_max=V_an_3,
                  aeration=None, suspended_growth_model=asm2d)

    O2 = su.CSTR('O2', A3-0, V_max=V_ae_2, aeration= 2, 
                 DO_ID='S_O2', suspended_growth_model=asm2d)

    C1 = su.FlatBottomCircularClarifier('C1', O2-0, [effluent, RAS, WAS],
                                        underflow=Q_ras, wastage=Q_was, 
                                        surface_area= 24226, height=6, N_layer=10, 
                                        feed_layer=5, X_threshold=3000, v_max=474, 
                                        v_max_practical=250, rh=5.76e-4, rp=2.86e-3, 
                                        fns=2.28e-3)

    sys = qs.System('metro_ASM2d', path=(B1, B2, B3, B4, B5, M1, M2, A1, A2, O1, 
                                         S1, A3, O2, C1), recycle = [RAS, recycle])

    return sys
#%%
@time_printer

# =============================================================================
# def create_components():
#      asm1_cmps = pc.create_asm1_cmps()
#      CO2 = qs.Component.from_chemical('S_CO2', search_ID='CO2', particle_size='Dissolved gas', degradability='Undegradable', organic=False)
#      CH4 = qs.Component.from_chemical('S_CH4', search_ID='CH4', particle_size='Dissolved gas', degradability='Undegradable', organic=False)
#      H2 = qs.Component.from_chemical('S_H2', search_ID='H2', particle_size='Dissolved gas', degradability='Undegradable', organic=False)
#      cmps1 = qs.Components.load_default()
#      ash = cmps1.X_Ig_ISS.copy('ash')
#      cmps = qs.Components([*asm1_cmps, CO2, CH4, H2, ash])
#      cmps.compile()
#      return cmps
# 
# =============================================================================
def run(t, t_step, method=None, **kwargs):
    sys = create_system()
    # for u in sys.units:
    #     if u.ID in ('O1', 'A1', 'A2', 'A3', 'O2', 'O3', 'O4', 'O5', 'O6', 'O7', 'O8', 'O9'):
    #         u.set_init_conc(**default_init_conds)
    bw_path = os.path.dirname(__file__)
    excel_path = os.path.join(bw_path, "data/initial_conditions_ASM2d.xlsx") 
    batch_init(sys, excel_path, sheet='t=10')

    # RAS = sys.flowsheet.stream.RAS
    # C1 = sys.flowsheet.unit.C1
    # sys.set_dynamic_tracker(RAS, C1)
    sys.set_dynamic_tracker(*sys.products)

    sys.simulate(
        state_reset_hook='reset_cache',
        t_span=(0,t),
        t_eval=np.arange(0, t+t_step, t_step),
        method=method,
        # rtol=1e-2,
        # atol=1e-3,
        # export_state_to=f'results/sol_{t}d_{method}.xlsx',
        **kwargs)
    sys.diagram()
    return sys

if __name__ == '__main__':
    t = 1
    t_step = 1
    method = 'RK45'
    # method = 'RK23' 
    # method = 'DOP853'
    # method = 'Radau'
    # method = 'BDF'
    # method = 'LSODA'
    msg = f'Method {method}'
    print(f'\n{msg}\n{"-"*len(msg)}') # long live OCD!
    print(f'Time span 0-{t}d \n')
    system = run(t, t_step, method=method)

    act_units = [u.ID for u in system.units if isinstance(u, (su.CSTR, su.FlatBottomCircularClarifier))]
    fs = system.flowsheet.stream
    srt = get_SRT(system, biomass_IDs, wastage= [fs.WAS, fs.effluent], active_unit_IDs=act_units)
    print(f'Estimated SRT assuming at steady state is {round(srt, 2)} days')

    #cmps = qs.get_components()
    #cmps = create_components()
    f = qs.main_flowsheet
    unit = system.flowsheet.unit
    fs = system.flowsheet.stream
    #fig, axis = fs.RAS.scope.plot_time_series(( 'S_NH')) 
    # fig, axis = fs.RAS.scope.plot_time_series(('S_S','S_NH')) 
    fig, axis = fs.effluent.scope.plot_time_series(('S_F', 'S_A', 'S_NH4', 'S_NO3')) 
    fig
# =============================================================================
#     fig, axis = unit.C1.scope.plot_time_series(('S_S', 'S_NH')) 
#     fig
# =============================================================================

flowsheet = qs.Flowsheet('bsm1')
qs.main_flowsheet.set_flowsheet(flowsheet)

# Components and stream
cmps = create_components()
qs.set_thermo(cmps)
wastewater = WasteStream('wastewater', T=Temp)
#inf_kwargs = inf_kwargs or default_inf_kwargs
inf_kwargs = beijing_inf_kwargs 
wastewater.set_flow_by_concentration(Q, **inf_kwargs)

effluent = WasteStream('effluent', T=Temp)
WAS = WasteStream('WAS', T=Temp)
RAS = WasteStream('RAS', T=Temp)
flowthrough = WasteStream('flowthrough', T=Temp)
recycle = WasteStream('recycle', T=Temp)
asm2d = pc.ASM2d(iP_SF=0.005, iP_XS=0.005, iP_XI=0.005, iN_BM=0.1, iTSS_XI=0.72)

__all__ = (
    'biomass_IDs',
    'create_system',
    'default_asm2d_kwargs', 'beijing_inf_kwargs', 'default_init_conds',
    'Q', 'Q_ras', 'Q_was', 'Temp', 'V_an_1', 'V_an_2', 'V_an_3', 'V_ae_1', 'V_ae_2', 
    'wastewater', 'inf_kwargs', 'effluent', 'WAS', 'RAS', 'flowthrough', 'recycle','asm2d'
    )