bsm2_base ¶

def step(self, i: int, *args, **kwargs):
    # --8<-- [end:step_8]
    """Simulates one time step of the BSM2 model.

    Parameters
    ----------
    i : int
        Index of the current time step [-].
    """

    # --8<-- [start:step_9]
    # timestep = timesteps[i]
    step: float = self.simtime[i]
    stepsize: float = self.timesteps[i]

    self.reactor1.kla, self.reactor2.kla, self.reactor3.kla, self.reactor4.kla, self.reactor5.kla = self.klas
    # --8<-- [end:step_9]

    # --8<-- [start:step_10]
    # get influent data that is smaller than and closest to current time step
    y_in_timestep = self.y_in[np.where(self.data_time <= step)[0][-1], :]

    iqi = self.performance.iqi(y_in_timestep)[0]
    self.iqi_all[i] = iqi
    # --8<-- [end:step_10]

    # --8<-- [start:step_11]
    yp_in_c, y_in_bp = self.input_splitter.output(y_in_timestep, (0.0, 0.0), float(reginit.QBYPASS))
    y_plant_bp, y_in_as_c = self.bypass_plant.output(y_in_bp, (1 - reginit.QBYPASSPLANT, reginit.QBYPASSPLANT))

    yp_in = self.combiner_primclar_pre.output(yp_in_c, self.yst_sp_p, self.yt_sp_p)
    self.yp_uf, self.yp_of, self.yp_internal = self.primclar.output(stepsize, step, yp_in)
    y_c_as_bp = self.combiner_primclar_post.output(self.yp_of[:21], y_in_as_c)
    # --8<-- [end:step_11]

    # --8<-- [start:step_12]
    y_bp_as, y_as_bp_c_eff = self.bypass_reactor.output(y_c_as_bp, (1 - reginit.QBYPASSAS, reginit.QBYPASSAS))
    y_in1 = self.combiner_reactor.output(self.ys_r, y_bp_as, self.yst_sp_as, self.yt_sp_as, self.y_out5_r)
    self.y_out1 = self.reactor1.output(stepsize, step, y_in1)
    self.y_out2 = self.reactor2.output(stepsize, step, self.y_out1)
    self.y_out3 = self.reactor3.output(stepsize, step, self.y_out2)
    self.y_out4 = self.reactor4.output(stepsize, step, self.y_out3)
    self.y_out5 = self.reactor5.output(stepsize, step, self.y_out4)
    ys_in, self.y_out5_r = self.splitter_reactor.output(
        self.y_out5, (max(self.y_out5[14] - self.qintr, 0.0), float(self.qintr))
    )
    self.ys_r, self.ys_was, self.ys_of, _, self.ys_tss_internal = self.settler.output(stepsize, step, ys_in)
    # --8<-- [end:step_12]

    # --8<-- [start:step_13]
    self.y_eff = self.combiner_effluent.output(y_plant_bp, y_as_bp_c_eff, self.ys_of)

    eqi = self.performance.eqi(self.ys_of, y_plant_bp, y_as_bp_c_eff)[0]
    self.eqi_all[i] = eqi

    self.yt_uf, yt_of = self.thickener.output(self.ys_was)
    self.yt_sp_p, self.yt_sp_as = self.splitter_thickener.output(
        yt_of, (1 - reginit.QTHICKENER2AS, reginit.QTHICKENER2AS)
    )
    # --8<-- [end:step_13]

    # --8<-- [start:step_14]
    self.yd_in = self.combiner_adm1.output(self.yt_uf, self.yp_uf)
    self.yi_out2, self.yd_out, _ = self.adm1_reactor.output(stepsize, step, self.yd_in, reginit.T_OP)

    self.ydw_s, ydw_r = self.dewatering.output(self.yi_out2)

    self.yst_out, self.yst_vol = self.storage.output(stepsize, step, ydw_r, reginit.QSTORAGE)
    self.yst_sp_p, self.yst_sp_as = self.splitter_storage.output(
        self.yst_out, (1 - reginit.QSTORAGE2AS, reginit.QSTORAGE2AS)
    )
    # --8<-- [end:step_14]

    # --8<-- [start:step_15]
    vol = np.array(
        [
            self.reactor1.volume,
            self.reactor2.volume,
            self.reactor3.volume,
            self.reactor4.volume,
            self.reactor5.volume,
            self.adm1_reactor.volume_liq,
        ]
    )
    sosat = np.array([asm1init.SOSAT1, asm1init.SOSAT2, asm1init.SOSAT3, asm1init.SOSAT4, asm1init.SOSAT5])
    self.ae = self.performance.aerationenergy_step(self.klas, vol[0:5], sosat)
    flows = np.array([self.qintr, asm1init.QR, asm1init.QW, self.yp_uf[14], self.yt_uf[14], self.ydw_s[14]])
    self.pe = self.performance.pumpingenergy_step(flows, pp_init.PP_PAR[10:16])
    self.me = self.performance.mixingenergy_step(self.klas, vol, pp_init.PP_PAR[16])
    # --8<-- [end:step_15]

    # --8<-- [start:step_16]
    tss_mass = self.performance.tss_mass_bsm2(
        self.yp_of,
        self.yp_uf,
        self.yp_internal,
        self.y_out1,
        self.y_out2,
        self.y_out3,
        self.y_out4,
        self.y_out5,
        self.ys_tss_internal,
        self.yd_out,
        self.yst_out,
        self.yst_vol,
    )
    # --8<-- [end:step_16]

    # --8<-- [start:step_17]
    ydw_s_tss_flow = self.performance.tss_flow(self.ydw_s)
    y_eff_tss_flow = self.performance.tss_flow(self.y_eff)
    carb = reginit.CARB1 + reginit.CARB2 + reginit.CARB3 + reginit.CARB4 + reginit.CARB5
    added_carbon_mass = self.performance.added_carbon_mass(carb, reginit.CARBONSOURCECONC)
    self.heat_demand = self.performance.heat_demand_step(self.yd_in, reginit.T_OP)[0]
    ch4_prod, h2_prod, co2_prod, q_gas = self.performance.gas_production(self.yd_out, reginit.T_OP)
    # --8<-- [end:step_17]

    # --8<-- [start:step_18]
    # This calculates an approximate oci value for each time step,
    # neglecting changes in the tss mass inside the whole plant
    self.oci_all[i] = self.performance.oci(
        self.pe * 24,
        self.ae * 24,
        self.me * 24,
        ydw_s_tss_flow,
        added_carbon_mass,
        self.heat_demand * 24,
        ch4_prod,
    )
    # --8<-- [end:step_18]

    # --8<-- [start:step_19]
    # These values are used to calculate the exact performance values at the end of the simulation
    self.perf_factors_all[i, :12] = [
        self.pe * 24,
        self.ae * 24,
        self.me * 24,
        ydw_s_tss_flow,
        y_eff_tss_flow,
        tss_mass,
        added_carbon_mass,
        self.heat_demand * 24,
        ch4_prod,
        h2_prod,
        co2_prod,
        q_gas,
    ]
    # --8<-- [end:step_19]

    # --8<-- [start:step_20]
    self.y_in_all[i] = y_in_timestep
    self.y_eff_all[i] = self.y_eff
    self.y_in_bp_all[i] = y_in_bp
    self.to_primary_all[i] = yp_in_c
    self.prim_in_all[i] = yp_in
    self.qpass_plant_all[i] = y_plant_bp
    self.qpassplant_to_as_all[i] = y_in_as_c
    self.qpassAS_all[i] = y_as_bp_c_eff
    self.to_as_all[i] = y_bp_as
    self.feed_settler_all[i] = ys_in
    self.qthick2AS_all[i] = self.yt_sp_as
    self.qthick2prim_all[i] = self.yt_sp_p
    self.qstorage2AS_all[i] = self.yst_sp_as
    self.qstorage2prim_all[i] = self.yst_sp_p
    self.sludge_all[i] = self.ydw_s

    # data for calculation of final oci
    self.violation_all[i] = self.performance.violation_step(self.y_eff[SNH], 4)[0]
    self.y_out1_all[i] = self.y_out1
    self.y_out2_all[i] = self.y_out2
    self.y_out3_all[i] = self.y_out3
    self.y_out4_all[i] = self.y_out4
    self.y_out5_all[i] = self.y_out5
    self.ys_r_all[i] = self.ys_r
    self.ys_was_all[i] = self.ys_was
    self.ys_of_all[i] = self.ys_of
    self.ys_tss_internal_all[i] = self.ys_tss_internal
    self.yp_uf_all[i] = self.yp_uf
    self.yp_of_all[i] = self.yp_of
    self.yt_uf_all[i] = self.yt_uf
    self.yd_out_all[i] = self.yd_out
    self.yi_out2_all[i] = self.yi_out2
    self.yst_out_all[i] = self.yst_out
    self.yst_vol_all[i] = self.yst_vol, step
    self.ydw_s_all[i] = self.ydw_s
    self.yp_internal_all[i] = self.yp_internal

def stabilize(self, atol: float = 1e-3):
    """Stabilizes the plant.

    Parameters
    ----------
    atol : float (optional)
        Absolute tolerance for the stabilization. <br>
        Default is 1e-3.

    Returns
    -------
    stable : bool
        Returns `True` if plant is stabilized after iterations.
    """
    check_vars = [
        'y_eff',
        'y_out1',
        'y_out2',
        'y_out3',
        'y_out4',
        'y_out5',
        'yt_sp_as',
        'yt_sp_p',
        'yst_sp_as',
        'yst_sp_p',
        'ydw_s',
        'y_out1',
        'y_out2',
        'y_out3',
        'y_out4',
        'y_out5',
        'ys_of',
        'yp_uf',
    ]
    stable = super()._stabilize(check_vars=check_vars, atol=atol)
    return stable

def simulate(self, *, plot=True, export=True):
    """Simulates the plant.

    Parameters
    ----------
    plot : bool (optional)
        If `True`, the data is plotted. <br>
        Default is `True`.
    export : bool (optional)
        If `True`, the data is exported. <br>
        Default is `True`.
    """

    for i, _ in enumerate(tqdm(self.simtime)):
        self.step(i)

        if i == 0:
            self.evaluator.add_new_data('iqi', 'iqi')
            self.evaluator.add_new_data('eqi', 'eqi')
            self.evaluator.add_new_data('oci', 'oci')
            self.evaluator.add_new_data('oci_final', 'oci_final')
        if self.evaltime[0] <= self.simtime[i] <= self.evaltime[1]:
            self.evaluator.update_data('iqi', self.iqi_all[i], self.simtime[i])
            self.evaluator.update_data('eqi', self.eqi_all[i], self.simtime[i])
            self.evaluator.update_data('oci', self.oci_all[i], self.simtime[i])

    self.oci_final = self.get_final_performance()[-1]
    self.evaluator.update_data('oci_final', self.oci_final, self.evaltime[1])
    self.finish_evaluation(plot=plot)

def finish_evaluation(self, *, plot=True):
    """Finishes the evaluation of the plant.

    Parameters
    ----------
    plot : bool (optional)
        If `True`, the data is plotted. <br>
        Default is `True`.
    """

    if plot:
        self.evaluator.plot_data()
    self.evaluator.export_data()

def get_violations(self, comp: tuple = ('SNH',), lim: tuple = (4,)):
    """Returns the violations of the given components within the evaluation interval.

    Parameters
    ----------
    comp : tuple(str) (optional)
        List of components to check for violations. <br>
        Default is ('SNH').
    lim : tuple(float) (optional)
        List of limits for the components. <br>
        Default is (4).

    Returns
    -------
    violations : dict{'comp': float}
        Dictionary with the components as keys and the violation time [d].
    """

    comp_dict = {
        'SI': 0,
        'SS': 1,
        'XI': 2,
        'XS': 3,
        'XBH': 4,
        'XBA': 5,
        'XP': 6,
        'SO': 7,
        'SNO': 8,
        'SNH': 9,
        'SND': 10,
        'XND': 11,
        'SALK': 12,
        'TSS': 13,
        'Q': 14,
        'TEMP': 15,
        'SD1': 16,
        'SD2': 17,
        'SD3': 18,
        'XD4': 19,
        'XD5': 20,
    }
    if len(comp) != len(lim):
        raise ValueError('The length of comp and lim has to be the same.')
    for c in comp:
        if c not in comp_dict:
            err = f'The component {c} is not in the list of components.'
            raise ValueError(err)
    comps = [comp_dict[c] for c in comp]

    violations = {}
    for i in range(len(comp)):
        comp_eff = self.y_eff_all[self.eval_idx[0] : self.eval_idx[1], comps[i]]
        violations[comp[i]] = np.sum(self.performance.violation_step(comp_eff, lim[i])) / 60 / 24
    return violations

def get_final_performance(self):
    """Returns the final performance values for evaluation period.

    Returns
    -------
    iqi_eval : float
        Final iqi value [kg ⋅ d⁻¹].
    eqi_eval : float
        Final eqi value [kg ⋅ d⁻¹].
    tot_sludge_prod : float
        Total sludge production [kg ⋅ d⁻¹].
    tot_tss_mass : float
        Total tss mass [kg ⋅ d⁻¹].
    carb_mass : float
        Carbon mass [kg(COD) ⋅ d⁻¹].
    ch4_prod : float
        Methane production [kg(CH4) ⋅ d⁻¹].
    h2_prod : float
        Hydrogen production [kg(H2) ⋅ d⁻¹].
    co2_prod : float
        Carbon dioxide production [kg(CO2) ⋅ d⁻¹].
    q_gas : float
        Total gas production [m³ ⋅ d⁻¹].
    heat_demand : float
        Heat demand [kWh ⋅ d⁻¹].
    mixingenergy : float
        Mixing energy [kWh ⋅ d⁻¹].
    pumpingenergy : float
        Pumping energy [kWh ⋅ d⁻¹].
    aerationenergy : float
        Aeration energy [kWh ⋅ d⁻¹].
    oci_eval : float
        Final operational cost index value [-].
    """

    # calculate the final performance values

    num_eval_timesteps = self.eval_idx[1] - self.eval_idx[0]

    iqi_eval = np.sum(self.iqi_all[self.eval_idx[0] : self.eval_idx[1]]) / num_eval_timesteps
    eqi_eval = np.sum(self.eqi_all[self.eval_idx[0] : self.eval_idx[1]]) / num_eval_timesteps

    oci_factors_eval = self.perf_factors_all[self.eval_idx[0] : self.eval_idx[1]]
    pumpingenergy = np.sum(oci_factors_eval[:, 0]) / num_eval_timesteps
    aerationenergy = np.sum(oci_factors_eval[:, 1]) / num_eval_timesteps
    mixingenergy = np.sum(oci_factors_eval[:, 2]) / num_eval_timesteps
    tot_tss_mass = np.sum(oci_factors_eval[:, 3]) / num_eval_timesteps + (
        oci_factors_eval[-1, 5] - oci_factors_eval[0, 5]
    ) / (self.evaltime[-1] - self.evaltime[0])
    tot_sludge_prod = (np.sum(oci_factors_eval[:, 4]) + np.sum(oci_factors_eval[:, 3])) / num_eval_timesteps + (
        oci_factors_eval[-1, 5] - oci_factors_eval[0, 5]
    ) / (self.evaltime[-1] - self.evaltime[0])
    carb_mass = np.sum(oci_factors_eval[:, 6]) / num_eval_timesteps
    heat_demand = np.sum(oci_factors_eval[:, 7]) / num_eval_timesteps
    ch4_prod = np.sum(oci_factors_eval[:, 8]) / num_eval_timesteps
    h2_prod = np.sum(oci_factors_eval[:, 9]) / num_eval_timesteps
    co2_prod = np.sum(oci_factors_eval[:, 10]) / num_eval_timesteps
    q_gas = np.sum(oci_factors_eval[:, 11]) / num_eval_timesteps

    oci_eval = self.performance.oci(
        pumpingenergy, aerationenergy, mixingenergy, tot_tss_mass, carb_mass, heat_demand, ch4_prod
    )

    return (
        iqi_eval,
        eqi_eval,
        tot_sludge_prod,
        tot_tss_mass,
        carb_mass,
        ch4_prod,
        h2_prod,
        co2_prod,
        q_gas,
        heat_demand,
        mixingenergy,
        pumpingenergy,
        aerationenergy,
        oci_eval,
    )