solph.components._generic_storage.GenericInvestmentStorageBlock._create()

Complexity

Conditions

52

Size

Total Lines	651
Code Lines	354

Duplication

Lines	25
Ratio	3.84 %

Importance

Changes

0

# -*- coding: utf-8 -
    def _create(self, group):
        """Create a storage block for investment modeling"""
        m = self.parent_block()

        # ########################## CHECKS ###################################
        if m.es.periods is not None:
            for n in group:
                error_fixed_absolute_losses = (
                    "For a multi-period investment model, fixed absolute"
                    " losses are not supported. Please remove parameter."
                )
                if n.fixed_losses_absolute[0] != 0:
                    raise ValueError(error_fixed_absolute_losses)
                error_initial_storage_level = (
                    "For a multi-period model, initial_storage_level is"
                    " not supported.\nIt needs to be removed since it"
                    " has no effect.\nstorage_content will be zero,"
                    " until there is some usable storage capacity installed."
                )
                if n.initial_storage_level is not None:
                    raise ValueError(error_initial_storage_level)

        # ########################## SETS #####################################

        self.INVESTSTORAGES = Set(initialize=[n for n in group])

        self.CONVEX_INVESTSTORAGES = Set(
            initialize=[n for n in group if n.investment.nonconvex is False]
        )

        self.NON_CONVEX_INVESTSTORAGES = Set(
            initialize=[n for n in group if n.investment.nonconvex is True]
        )

        self.INVESTSTORAGES_BALANCED = Set(
            initialize=[n for n in group if n.balanced is True]
        )

        self.INVESTSTORAGES_NO_INIT_CONTENT = Set(
            initialize=[n for n in group if n.initial_storage_level is None]
        )

        self.INVESTSTORAGES_INIT_CONTENT = Set(
            initialize=[
                n for n in group if n.initial_storage_level is not None
            ]
        )

        self.INVEST_REL_CAP_IN = Set(
            initialize=[
                n
                for n in group
                if n.invest_relation_input_capacity[0] is not None
            ]
        )

        self.INVEST_REL_CAP_OUT = Set(
            initialize=[
                n
                for n in group
                if n.invest_relation_output_capacity[0] is not None
            ]
        )

        self.INVEST_REL_IN_OUT = Set(
            initialize=[
                n
                for n in group
                if n.invest_relation_input_output[0] is not None
            ]
        )

        # The storage content is a non-negative variable, therefore it makes no
        # sense to create an additional constraint if the lower bound is zero
        # for all time steps.
        self.MIN_INVESTSTORAGES = Set(
            initialize=[
                n
                for n in group
                if sum([n.min_storage_level[t] for t in m.TIMESTEPS]) > 0
            ]
        )

        self.OVERALL_MAXIMUM_INVESTSTORAGES = Set(
            initialize=[
                n for n in group if n.investment.overall_maximum is not None
            ]
        )

        self.OVERALL_MINIMUM_INVESTSTORAGES = Set(
            initialize=[
                n for n in group if n.investment.overall_minimum is not None
            ]
        )

        self.EXISTING_INVESTSTORAGES = Set(
            initialize=[n for n in group if n.investment.existing is not None]
        )

        # ######################### Variables  ################################
        if not m.TSAM_MODE:
            self.storage_content = Var(
                self.INVESTSTORAGES, m.TIMEPOINTS, within=NonNegativeReals
            )
        else:
            self.inter_storage_content = Var(
                self.INVESTSTORAGES, m.CLUSTERS_OFFSET, within=NonNegativeReals
            )
            self.intra_storage_delta = Var(
                self.INVESTSTORAGES, m.TIMEINDEX_TYPICAL_CLUSTER_OFFSET
            )
            # set the initial intra storage content
            # first timestep in intra storage is always zero
            for n in group:
                for p, k in m.TYPICAL_CLUSTERS:
                    self.intra_storage_delta[n, p, k, 0] = 0
                    self.intra_storage_delta[n, p, k, 0].fix()

        def _storage_investvar_bound_rule(_, n, p):
            """
            Rule definition to bound the invested storage capacity `invest`.
            """
            if n in self.CONVEX_INVESTSTORAGES:
                return n.investment.minimum[p], n.investment.maximum[p]
            else:  # n in self.NON_CONVEX_INVESTSTORAGES
                return 0, n.investment.maximum[p]

        self.invest = Var(
            self.INVESTSTORAGES,
            m.PERIODS,
            within=NonNegativeReals,
            bounds=_storage_investvar_bound_rule,
        )

        # Total capacity
        self.total = Var(
            self.INVESTSTORAGES,
            m.PERIODS,
            within=NonNegativeReals,
            initialize=0,
        )

        if m.es.periods is not None:
            # Old capacity to be decommissioned (due to lifetime)
            self.old = Var(
                self.INVESTSTORAGES, m.PERIODS, within=NonNegativeReals
            )

            # Old endogenous capacity to be decommissioned (due to lifetime)
            self.old_end = Var(
                self.INVESTSTORAGES, m.PERIODS, within=NonNegativeReals
            )

            # Old exogenous capacity to be decommissioned (due to lifetime)
            self.old_exo = Var(
                self.INVESTSTORAGES, m.PERIODS, within=NonNegativeReals
            )

        # create status variable for a non-convex investment storage
        self.invest_status = Var(
            self.NON_CONVEX_INVESTSTORAGES, m.PERIODS, within=Binary
        )

        # ######################### CONSTRAINTS ###############################
        i = {n: [i for i in n.inputs][0] for n in group}
        o = {n: [o for o in n.outputs][0] for n in group}

        # Handle unit lifetimes
        def _total_storage_capacity_rule(block):
            """Rule definition for determining total installed
            capacity (taking decommissioning into account)
            """
            for n in self.INVESTSTORAGES:
                for p in m.PERIODS:
                    if p == 0:
                        expr = (
                            self.total[n, p]
                            == self.invest[n, p] + n.investment.existing
                        )
                        self.total_storage_rule.add((n, p), expr)
                    else:
                        expr = (
                            self.total[n, p]
                            == self.invest[n, p]
                            + self.total[n, p - 1]
                            - self.old[n, p]
                        )
                        self.total_storage_rule.add((n, p), expr)

        self.total_storage_rule = Constraint(
            self.INVESTSTORAGES, m.PERIODS, noruleinit=True
        )

        self.total_storage_rule_build = BuildAction(
            rule=_total_storage_capacity_rule
        )

        # multi-period storage implementation for time intervals
        if m.es.periods is not None:

            def _old_storage_capacity_rule_end(block):
                """Rule definition for determining old endogenously installed
                capacity to be decommissioned due to reaching its lifetime.
                Investment and decommissioning periods are linked within
                the constraint. The respective decommissioning period is
                determined for every investment period based on the components
                lifetime and a matrix describing its age of each endogenous
                investment. Decommissioning can only occur at the beginning of
                each period.

                Note
                ----
                For further information on the implementation check
                PR#957 https://github.com/oemof/oemof-solph/pull/957
                """
                for n in self.INVESTSTORAGES:
                    lifetime = n.investment.lifetime
                    if lifetime is None:
                        msg = (
                            "You have to specify a lifetime "
                            "for a Flow going into or out of "
                            "a GenericStorage unit "
                            "in a multi-period model!"
                            f" Value for {n} is missing."
                        )
                        raise ValueError(msg)
                    # get the period matrix describing the temporal distance
                    # between all period combinations.
                    periods_matrix = m.es.periods_matrix

                    # get the index of the minimum value in each row greater
                    # equal than the lifetime. This value equals the
                    # decommissioning period if not zero. The index of this
                    # value represents the investment period. If np.where
                    # condition is not met in any row, min value will be zero
                    decomm_periods = np.argmin(
                        np.where(
                            (periods_matrix >= lifetime),
                            periods_matrix,
                            np.inf,
                        ),
                        axis=1,
                    )

                    # no decommissioning in first period
                    expr = self.old_end[n, 0] == 0
                    self.old_rule_end.add((n, 0), expr)

                    # all periods not in decomm_periods have no decommissioning
                    # zero is excluded
                    for p in m.PERIODS:
                        if p not in decomm_periods and p != 0:
                            expr = self.old_end[n, p] == 0
                            self.old_rule_end.add((n, p), expr)

                    # multiple invests can be decommissioned in the same period
                    # but only sequential ones, thus a bookkeeping is
                    # introduced andconstraints are added to equation one
                    # iteration later.
                    last_decomm_p = np.nan
                    # loop over invest periods (values are decomm_periods)
                    for invest_p, decomm_p in enumerate(decomm_periods):
                        # Add constraint of iteration before
                        # (skipped in first iteration by last_decomm_p = nan)
                        if (decomm_p != last_decomm_p) and (
                            last_decomm_p is not np.nan
                        ):
                            expr = self.old_end[n, last_decomm_p] == expr
                            self.old_rule_end.add((n, last_decomm_p), expr)

                        # no decommissioning if decomm_p is zero
                        if decomm_p == 0:
                            # overwrite decomm_p with zero to avoid
                            # chaining invest periods in next iteration
                            last_decomm_p = 0

                        # if decomm_p is the same as the last one chain invest
                        # period
                        elif decomm_p == last_decomm_p:
                            expr += self.invest[n, invest_p]
                            # overwrite decomm_p
                            last_decomm_p = decomm_p

                        # if decomm_p is not zero, not the same as the last one
                        # and it's not the first period
                        else:
                            expr = self.invest[n, invest_p]
                            # overwrite decomm_p
                            last_decomm_p = decomm_p

                    # Add constraint of very last iteration
                    if last_decomm_p != 0:
                        expr = self.old_end[n, last_decomm_p] == expr
                        self.old_rule_end.add((n, last_decomm_p), expr)

            self.old_rule_end = Constraint(
                self.INVESTSTORAGES, m.PERIODS, noruleinit=True
            )

            self.old_rule_end_build = BuildAction(
                rule=_old_storage_capacity_rule_end
            )

            def _old_storage_capacity_rule_exo(block):
                """Rule definition for determining old exogenously given
                capacity to be decommissioned due to reaching its lifetime
                """
                for n in self.INVESTSTORAGES:
                    age = n.investment.age
                    lifetime = n.investment.lifetime
                    is_decommissioned = False
                    for p in m.PERIODS:
                        # No shutdown in first period
                        if p == 0:
                            expr = self.old_exo[n, p] == 0
                            self.old_rule_exo.add((n, p), expr)
                        elif lifetime - age <= m.es.periods_years[p]:
                            # Track decommissioning status
                            if not is_decommissioned:
                                expr = (
                                    self.old_exo[n, p] == n.investment.existing
                                )
                                is_decommissioned = True
                            else:
                                expr = self.old_exo[n, p] == 0
                            self.old_rule_exo.add((n, p), expr)
                        else:
                            expr = self.old_exo[n, p] == 0
                            self.old_rule_exo.add((n, p), expr)

            self.old_rule_exo = Constraint(
                self.INVESTSTORAGES, m.PERIODS, noruleinit=True
            )

            self.old_rule_exo_build = BuildAction(
                rule=_old_storage_capacity_rule_exo
            )

            def _old_storage_capacity_rule(block):
                """Rule definition for determining (overall) old capacity
                to be decommissioned due to reaching its lifetime
                """
                for n in self.INVESTSTORAGES:
                    for p in m.PERIODS:
                        expr = (
                            self.old[n, p]
                            == self.old_end[n, p] + self.old_exo[n, p]
                        )
                        self.old_rule.add((n, p), expr)

            self.old_rule = Constraint(
                self.INVESTSTORAGES, m.PERIODS, noruleinit=True
            )

            self.old_rule_build = BuildAction(rule=_old_storage_capacity_rule)

            def _initially_empty_rule(_):
                """Ensure storage to be empty initially"""
                for n in self.INVESTSTORAGES:
                    expr = self.storage_content[n, 0] == 0
                    self.initially_empty.add((n, 0), expr)

            if not m.TSAM_MODE:
                # inter and intra initial storage contents are handled above
                self.initially_empty = Constraint(
                    self.INVESTSTORAGES, m.TIMESTEPS, noruleinit=True
                )

                self.initially_empty_build = BuildAction(
                    rule=_initially_empty_rule
                )

        # Standard storage implementation for discrete time points
        else:

            def _inv_storage_init_content_max_rule(block, n):
                """Constraint for a variable initial storage capacity."""
                if not m.TSAM_MODE:
                    lhs = block.storage_content[n, 0]
                else:
                    lhs = block.intra_storage_delta[n, 0, 0, 0]
                return lhs <= n.investment.existing + block.invest[n, 0]

            self.init_content_limit = Constraint(
                self.INVESTSTORAGES_NO_INIT_CONTENT,
                rule=_inv_storage_init_content_max_rule,
            )

            def _inv_storage_init_content_fix_rule(block, n):
                """Constraint for a fixed initial storage capacity."""
                if not m.TSAM_MODE:
                    lhs = block.storage_content[n, 0]
                else:
                    lhs = block.intra_storage_delta[n, 0, 0, 0]
                return lhs == n.initial_storage_level * (
                    n.investment.existing + block.invest[n, 0]
                )

            self.init_content_fix = Constraint(
                self.INVESTSTORAGES_INIT_CONTENT,
                rule=_inv_storage_init_content_fix_rule,
            )

        def _storage_balance_rule(block, n, p, t):
            """
            Rule definition for the storage balance of every storage n and
            every timestep.
            """
            expr = 0
            expr += block.storage_content[n, t + 1]
            expr += (
                -block.storage_content[n, t]
                * (1 - n.loss_rate[t]) ** m.timeincrement[t]
            )
            expr += (
                n.fixed_losses_relative[t]
                * self.total[n, p]
                * m.timeincrement[t]
            )
            expr += n.fixed_losses_absolute[t] * m.timeincrement[t]
            expr += (
                -m.flow[i[n], n, t] * n.inflow_conversion_factor[t]
            ) * m.timeincrement[t]
            expr += (
                m.flow[n, o[n], t] / n.outflow_conversion_factor[t]
            ) * m.timeincrement[t]
            return expr == 0

        def _intra_storage_balance_rule(block, n, p, k, g):

This code seems to be duplicated in your project.

            """
            Rule definition for the storage balance of every storage n and
            every timestep.
            """
            t = m.get_timestep_from_tsam_timestep(p, k, g)
            expr = 0
            expr += block.intra_storage_delta[n, p, k, g + 1]
            expr += (
                -block.intra_storage_delta[n, p, k, g]
                * (1 - n.loss_rate[t]) ** m.timeincrement[t]
            )
            expr += (
                n.fixed_losses_relative[t]
                * self.total[n, p]
                * m.timeincrement[t]
            )
            expr += n.fixed_losses_absolute[t] * m.timeincrement[t]
            expr += (
                -m.flow[i[n], n, t] * n.inflow_conversion_factor[t]
            ) * m.timeincrement[t]
            expr += (
                m.flow[n, o[n], t] / n.outflow_conversion_factor[t]
            ) * m.timeincrement[t]
            return expr == 0

        if not m.TSAM_MODE:
            self.balance = Constraint(
                self.INVESTSTORAGES,
                m.TIMEINDEX,
                rule=_storage_balance_rule,
            )
        else:
            self.intra_balance = Constraint(
                self.INVESTSTORAGES,
                m.TIMEINDEX_TYPICAL_CLUSTER,
                rule=_intra_storage_balance_rule,
            )

        def _inter_storage_balance_rule(block, n, i):
            """
            Rule definition for the storage balance of every storage n and
            every timestep.
            """
            ii = 0
            for p in m.PERIODS:
                ii += len(m.es.tsa_parameters[p]["order"])
                if ii > i:
                    ii -= len(m.es.tsa_parameters[p]["order"])
                    ii = i - ii
                    break

            k = m.es.tsa_parameters[p]["order"][ii]

The variable p does not seem to be defined in case the for loop on line 1810 is not entered. Are you sure this can never be the case?

            t = m.get_timestep_from_tsam_timestep(
                p, k, m.es.tsa_parameters[p]["timesteps"] - 1
            )
            expr = 0
            expr += block.inter_storage_content[n, i + 1]
            expr += -block.inter_storage_content[n, i] * (
                1 - n.loss_rate[t]
            ) ** (m.timeincrement[t] * m.es.tsa_parameters[p]["timesteps"])
            expr += -self.intra_storage_delta[
                n, p, k, m.es.tsa_parameters[p]["timesteps"]
            ]
            return expr == 0

        if m.TSAM_MODE:
            self.inter_balance = Constraint(
                self.INVESTSTORAGES,
                m.CLUSTERS,
                rule=_inter_storage_balance_rule,
            )

        if m.es.periods is None and not m.TSAM_MODE:

            def _balanced_storage_rule(block, n):
                return (
                    block.storage_content[n, m.TIMEPOINTS.at(-1)]
                    == block.storage_content[n, m.TIMEPOINTS.at(1)]
                )

            self.balanced_cstr = Constraint(
                self.INVESTSTORAGES_BALANCED, rule=_balanced_storage_rule
            )

        def _power_coupled(block):
            """
            Rule definition for constraint to connect the input power
            and output power
            """
            for n in self.INVEST_REL_IN_OUT:
                for p in m.PERIODS:
                    expr = (
                        m.InvestmentFlowBlock.total[n, o[n], p]
                    ) * n.invest_relation_input_output[p] == (
                        m.InvestmentFlowBlock.total[i[n], n, p]
                    )
                    self.power_coupled.add((n, p), expr)

        self.power_coupled = Constraint(
            self.INVEST_REL_IN_OUT, m.PERIODS, noruleinit=True
        )

        self.power_coupled_build = BuildAction(rule=_power_coupled)

        def _storage_capacity_inflow_invest_rule(block):
            """
            Rule definition of constraint connecting the inflow
            `InvestmentFlowBlock.invest of storage with invested capacity
            `invest` by nominal_storage_capacity__inflow_ratio
            """
            for n in self.INVEST_REL_CAP_IN:
                for p in m.PERIODS:
                    expr = (
                        m.InvestmentFlowBlock.total[i[n], n, p]
                        == self.total[n, p]
                        * n.invest_relation_input_capacity[p]
                    )
                    self.storage_capacity_inflow.add((n, p), expr)

        self.storage_capacity_inflow = Constraint(
            self.INVEST_REL_CAP_IN, m.PERIODS, noruleinit=True
        )

        self.storage_capacity_inflow_build = BuildAction(
            rule=_storage_capacity_inflow_invest_rule
        )

        def _storage_capacity_outflow_invest_rule(block):
            """
            Rule definition of constraint connecting outflow
            `InvestmentFlowBlock.invest` of storage and invested capacity
            `invest` by nominal_storage_capacity__outflow_ratio
            """
            for n in self.INVEST_REL_CAP_OUT:
                for p in m.PERIODS:
                    expr = (
                        m.InvestmentFlowBlock.total[n, o[n], p]
                        == self.total[n, p]
                        * n.invest_relation_output_capacity[p]
                    )
                    self.storage_capacity_outflow.add((n, p), expr)

        self.storage_capacity_outflow = Constraint(
            self.INVEST_REL_CAP_OUT, m.PERIODS, noruleinit=True
        )

        self.storage_capacity_outflow_build = BuildAction(
            rule=_storage_capacity_outflow_invest_rule
        )

        self._add_storage_limit_constraints()

        def maximum_invest_limit(block, n, p):
            """
            Constraint for the maximal investment in non convex investment
            storage.
            """
            return (
                n.investment.maximum[p] * self.invest_status[n, p]
                - self.invest[n, p]
            ) >= 0

        self.limit_max = Constraint(
            self.NON_CONVEX_INVESTSTORAGES,
            m.PERIODS,
            rule=maximum_invest_limit,
        )

        def smallest_invest(block, n, p):
            """
            Constraint for the minimal investment in non convex investment
            storage if the invest is greater than 0. So the invest variable
            can be either 0 or greater than the minimum.
            """
            return (
                self.invest[n, p]
                - n.investment.minimum[p] * self.invest_status[n, p]
                >= 0
            )

        self.limit_min = Constraint(
            self.NON_CONVEX_INVESTSTORAGES, m.PERIODS, rule=smallest_invest
        )

        if m.es.periods is not None:

            def _overall_storage_maximum_investflow_rule(block):
                """Rule definition for maximum overall investment
                in investment case.
                """
                for n in self.OVERALL_MAXIMUM_INVESTSTORAGES:
                    for p in m.PERIODS:
                        expr = self.total[n, p] <= n.investment.overall_maximum
                        self.overall_storage_maximum.add((n, p), expr)

            self.overall_storage_maximum = Constraint(
                self.OVERALL_MAXIMUM_INVESTSTORAGES, m.PERIODS, noruleinit=True
            )

            self.overall_maximum_build = BuildAction(
                rule=_overall_storage_maximum_investflow_rule
            )

            def _overall_minimum_investflow_rule(block):
                """Rule definition for minimum overall investment
                in investment case.

                Note: This is only applicable for the last period
                """
                for n in self.OVERALL_MINIMUM_INVESTSTORAGES:
                    expr = (
                        n.investment.overall_minimum
                        <= self.total[n, m.PERIODS[-1]]
                    )
                    self.overall_minimum.add(n, expr)

            self.overall_minimum = Constraint(
                self.OVERALL_MINIMUM_INVESTSTORAGES, noruleinit=True
            )

            self.overall_minimum_build = BuildAction(
                rule=_overall_minimum_investflow_rule
            )