Account for remaining values in multi-period models

docs/whatsnew/v0-5-2.rst

@@ -13,6 +13,11 @@ Documentation
 Bug fixes
 #########
 
+* Fix handling of investment annuities and fixed costs for multi-period models:
+  Limit to costs that occur within the optimization horizon to prevent a
+  bias towards investments happening earlier in the optimization horizon.
+* Fix bugs in multi-period documentation.
+
 Testing
 #######
 
@@ -25,3 +30,6 @@ Contributors
 ############
 
 * Patrik Schönfeldt
+* Johannes Kochems
+* Julian Endres
+* Hendrik Huyskens

src/oemof/solph/_energy_system.py

@@ -172,8 +172,10 @@ def __init__(
             )
             warnings.warn(msg, debugging.SuspiciousUsageWarning)
         self.periods = periods
-        self._extract_periods_years()
-        self._extract_periods_matrix()
+        if self.periods is not None:
+            self._extract_periods_years()
+            self._extract_periods_matrix()
+            self._extract_end_year_of_optimization()
 
     def _extract_periods_years(self):
         """Map years in optimization to respective period based on time indices
@@ -183,13 +185,12 @@ def _extract_periods_years(self):
         start of the optimization run and starting with 0.
         """
         periods_years = [0]
-        if self.periods is not None:
-            start_year = self.periods[0].min().year
-            for k, v in enumerate(self.periods):
-                if k >= 1:
-                    periods_years.append(v.min().year - start_year)
+        start_year = self.periods[0].min().year
+        for k, v in enumerate(self.periods):
+            if k >= 1:
+                periods_years.append(v.min().year - start_year)
 
-            self.periods_years = periods_years
+        self.periods_years = periods_years
 
     def _extract_periods_matrix(self):
         """Determines a matrix describing the temporal distance to each period.
@@ -200,13 +201,38 @@ def _extract_periods_matrix(self):
         between each investment period to each decommissioning period.
         """
         periods_matrix = []
-        if self.periods is not None:
-            period_years = np.array(self.periods_years)
-            for v in period_years:
-                row = period_years - v
-                row = np.where(row < 0, 0, row)
-                periods_matrix.append(row)
-            self.periods_matrix = np.array(periods_matrix)
+        period_years = np.array(self.periods_years)
+        for v in period_years:
+            row = period_years - v
+            row = np.where(row < 0, 0, row)
+            periods_matrix.append(row)
+        self.periods_matrix = np.array(periods_matrix)
+
+    def _extract_end_year_of_optimization(self):
+        """Extract the end of the optimization in years"""
+        duration_last_period = self.get_period_duration(-1)
+        self.end_year_of_optimization = (
+            self.periods_years[-1] + duration_last_period
+        )
+
+    def get_period_duration(self, period):
+        """Get duration of a period in full years
+
+        Parameters
+        ----------
+        period : int
+            Period for which the duration in years shall be obtained
+
+        Returns
+        -------
+        int
+            Duration of the period
+        """
+        return (
+            self.periods[period].max().year
+            - self.periods[period].min().year
+            + 1
+        )
 
 
 def create_time_index(

src/oemof/solph/components/_generic_storage.py

@@ -362,7 +362,7 @@ class GenericStorageBlock(ScalarBlock):
 
     Set storage_content of last time step to one at t=0 if balanced == True
         .. math::
-            E(t_{last}) = &E(-1)
+            E(t_{last}) = E(-1)
 
     Storage balance :attr:`om.Storage.balance[n, t]`
         .. math:: E(t) = &E(t-1) \cdot
@@ -375,8 +375,8 @@ class GenericStorageBlock(ScalarBlock):
     Connect the invest variables of the input and the output flow.
         .. math::
           InvestmentFlowBlock.invest(source(n), n, p) + existing = \\
-          (InvestmentFlowBlock.invest(n, target(n), p) + existing) * \\
-          invest\_relation\_input\_output(n) \\
+          (InvestmentFlowBlock.invest(n, target(n), p) + existing) \\
+          * invest\_relation\_input\_output(n) \\
           \forall n \in \textrm{INVEST\_REL\_IN\_OUT} \\
           \forall p \in \textrm{PERIODS}
 
@@ -429,7 +429,7 @@ class GenericStorageBlock(ScalarBlock):
 
     * :attr: `storage_costs` not 0
 
-        ..math::
+        .. math::
             \sum_{t \in \textrm{TIMESTEPS}} c_{storage}(t) \cdot E(t)
 
 
@@ -438,11 +438,14 @@ class GenericStorageBlock(ScalarBlock):
     * :attr:`fixed_costs` not None
 
         .. math::
-            \sum_{p \in \textrm{PERIODS}} E_{nom}
-            \cdot c_{fixed}(p) \cdot DF^{-p}
+            \displaystyle \sum_{pp=0}^{year_{max}} E_{nom}
+            \cdot c_{fixed}(pp) \cdot DF^{-pp}
 
     whereby:
-    :math:`DF=(1+dr)` is the discount factor with discount rate :math:`dr`
+
+    * :math:`DF=(1+dr)` is the discount factor with discount rate :math:`dr`
+    * :math:`year_{max}` denotes the last year of the optimization
+      horizon, i.e. at the end of the last period.
 
     """  # noqa: E501
 
@@ -599,12 +602,12 @@ def _objective_expression(self):
         if m.es.periods is not None:
             for n in self.STORAGES:
                 if n.fixed_costs[0] is not None:
-                    for p in m.PERIODS:
-                        fixed_costs += (
-                            n.nominal_storage_capacity
-                            * n.fixed_costs[p]
-                            * ((1 + m.discount_rate) ** -p)
-                        )
+                    fixed_costs += sum(
+                        n.nominal_storage_capacity
+                        * n.fixed_costs[pp]
+                        * ((1 + m.discount_rate) ** (-pp))
+                        for pp in range(m.es.end_year_of_optimization)
+                    )
         self.fixed_costs = Expression(expr=fixed_costs)
 
         storage_costs = 0
@@ -620,8 +623,9 @@ def _objective_expression(self):
                     )
 
         self.storage_costs = Expression(expr=storage_costs)
+        self.costs = Expression(expr=storage_costs + fixed_costs)
 
-        return self.fixed_costs + self.storage_costs
+        return self.costs
 
 
 class GenericInvestmentStorageBlock(ScalarBlock):
@@ -883,31 +887,36 @@ class GenericInvestmentStorageBlock(ScalarBlock):
                 E_{invest}(0) \cdot c_{invest,var}(0)
                 + c_{invest,fix}(0) \cdot b_{invest}(0)\\
 
+    Whereby 0 denotes the 0th (investment) period since
+    in a standard model, there is only this one period.
+
     *Multi-period model*
 
         * :attr:`nonconvex = False`
 
             .. math::
                 &
-                E_{invest}(p) \cdot A(c_{invest,var}(p), l, ir) \cdot l
-                \cdot DF^{-p}\\
+                E_{invest}(p) \cdot A(c_{invest,var}(p), l, ir)
+                \cdot \frac {1}{ANF(d, ir)} \cdot DF^{-p}\\
                 &
                 \forall p \in \textrm{PERIODS}
 
         * :attr:`nonconvex = True`
 
             .. math::
                 &
-                E_{invest}(p) \cdot A(c_{invest,var}(p), l, ir) \cdot l
-                \cdot DF^{-p} +  c_{invest,fix}(p) \cdot b_{invest}(p)\\
+                (E_{invest}(p) \cdot A(c_{invest,var}(p), l, ir)
+                \cdot \frac {1}{ANF(d, ir)}\\
+                &
+                +  c_{invest,fix}(p) \cdot b_{invest}(p)) \cdot DF^{-p} \\
                 &
                 \forall p \in \textrm{PERIODS}
 
         * :attr:`fixed_costs` not None for investments
 
             .. math::
                 &
-                \sum_{pp=year(p)}^{year(p)+l}
+                \sum_{pp=year(p)}^{limit_{end}}
                 E_{invest}(p) \cdot c_{fixed}(pp) \cdot DF^{-pp})
                 \cdot DF^{-p}\\
                 &
@@ -916,27 +925,45 @@ class GenericInvestmentStorageBlock(ScalarBlock):
         * :attr:`fixed_costs` not None for existing capacity
 
             .. math::
-                \sum_{pp=0}^{l-a} E_{exist} \cdot c_{fixed}(pp)
+                \sum_{pp=0}^{limit_{exo}} E_{exist} \cdot c_{fixed}(pp)
                 \cdot DF^{-pp}
 
 
-        whereby:
-
-        * :math:`A(c_{invest,var}(p), l, ir)` A is the annuity for
-          investment expenses :math:`c_{invest,var}(p)` lifetime :math:`l` and
-          interest rate :math:`ir`
-        * :math:`DF=(1+dr)` is the discount factor with discount rate math:`dr`
+    whereby:
 
-    The annuity hereby is:
+    * :math:`A(c_{invest,var}(p), l, ir)` A is the annuity for
+      investment expenses :math:`c_{invest,var}(p)`, lifetime :math:`l`
+      and interest rate :math:`ir`.
+    * :math:`ANF(d, ir)` is the annuity factor for duration :math:`d`
+      and interest rate :math:`ir`.
+    * :math:`d=min\{year_{max} - year(p), l\}` defines the
+      number of years within the optimization horizon that investment
+      annuities are accounted for.
+    * :math:`year(p)` denotes the start year of period :math:`p`.
+    * :math:`year_{max}` denotes the last year of the optimization
+      horizon, i.e. at the end of the last period.
+    * :math:`limit_{end}=min\{year_{max}, year(p) + l\}` is used as an
+      upper bound to ensure fixed costs for endogenous investments
+      to occur within the optimization horizon.
+    * :math:`limit_{exo}=min\{year_{max}, l - a\}` is used as an
+      upper bound to ensure fixed costs for existing capacities to occur
+      within the optimization horizon. :math:`a` is the initial age
+      of an asset.
+    * :math:`DF=(1+dr)` is the discount factor.
+
+    The annuity / annuity factor hereby is:
 
         .. math::
-
+            &
             A(c_{invest,var}(p), l, ir) = c_{invest,var}(p) \cdot
-                \frac {(1+i)^l \cdot i} {(1+i)^l - 1} \cdot
+                \frac {(1+i)^l \cdot i} {(1+i)^l - 1}\\
+            &\\
+            &
+            ANF(d, ir)=\frac {(1+ir)^d \cdot ir} {(1+ir)^d - 1}
 
-    It is retrieved, using oemof.tools.economics annuity function. The
-    interest rate is defined as a weighted average costs of capital (wacc) and
-    assumed constant over time.
+    They are retrieved, using oemof.tools.economics annuity function. The
+    interest rate :math:`i` for the annuity is defined as weighted
+    average costs of capital (wacc) and assumed constant over time.
 
     The overall summed cost expressions for all *InvestmentFlowBlock* objects
     can be accessed by
@@ -1712,7 +1739,6 @@ def _objective_expression(self):
                 "social planner point of view and does not reflect "
                 "microeconomic interest requirements."
             )
-
             for n in self.CONVEX_INVESTSTORAGES:
                 lifetime = n.investment.lifetime
                 interest = n.investment.interest_rate
@@ -1728,12 +1754,16 @@ def _objective_expression(self):
                         n=lifetime,
                         wacc=interest,
                     )
-                    investment_costs_increment = (
-                        self.invest[n, p]
-                        * annuity
-                        * lifetime
-                        * ((1 + m.discount_rate) ** (-m.es.periods_years[p]))
+                    duration = min(
+                        m.es.end_year_of_optimization - m.es.periods_years[p],
+                        lifetime,
+                    )
+                    present_value_factor = 1 / economics.annuity(
+                        capex=1, n=duration, wacc=interest
                     )
+                    investment_costs_increment = (
+                        self.invest[n, p] * annuity * present_value_factor
+                    ) * (1 + m.discount_rate) ** (-m.es.periods_years[p])
                     investment_costs += investment_costs_increment
                     period_investment_costs[p] += investment_costs_increment
 
@@ -1752,36 +1782,50 @@ def _objective_expression(self):
                         n=lifetime,
                         wacc=interest,
                     )
+                    duration = min(
+                        m.es.end_year_of_optimization - m.es.periods_years[p],
+                        lifetime,
+                    )
+                    present_value_factor = 1 / economics.annuity(
+                        capex=1, n=duration, wacc=interest
+                    )
                     investment_costs_increment = (
-                        self.invest[n, p] * annuity * lifetime
+                        self.invest[n, p] * annuity * present_value_factor
                         + self.invest_status[n, p] * n.investment.offset[p]
-                    ) * ((1 + m.discount_rate) ** (-m.es.periods_years[p]))
+                    ) * (1 + m.discount_rate) ** (-m.es.periods_years[p])
                     investment_costs += investment_costs_increment
                     period_investment_costs[p] += investment_costs_increment
 
             for n in self.INVESTSTORAGES:
                 if n.investment.fixed_costs[0] is not None:
                     lifetime = n.investment.lifetime
                     for p in m.PERIODS:
+                        range_limit = min(
+                            m.es.end_year_of_optimization,
+                            m.es.periods_years[p] + lifetime,
+                        )
                         fixed_costs += sum(
                             self.invest[n, p]
                             * n.investment.fixed_costs[pp]
-                            * ((1 + m.discount_rate) ** (-pp))
+                            * (1 + m.discount_rate) ** (-pp)
                             for pp in range(
                                 m.es.periods_years[p],
-                                m.es.periods_years[p] + lifetime,
+                                range_limit,
                             )
-                        ) * ((1 + m.discount_rate) ** (-m.es.periods_years[p]))
+                        ) * (1 + m.discount_rate) ** (-m.es.periods_years[p])
 
             for n in self.EXISTING_INVESTSTORAGES:
                 if n.investment.fixed_costs[0] is not None:
                     lifetime = n.investment.lifetime
                     age = n.investment.age
+                    range_limit = min(
+                        m.es.end_year_of_optimization, lifetime - age
+                    )
                     fixed_costs += sum(
                         n.investment.existing
                         * n.investment.fixed_costs[pp]
-                        * ((1 + m.discount_rate) ** (-pp))
-                        for pp in range(0, lifetime - age)
+                        * (1 + m.discount_rate) ** (-pp)
+                        for pp in range(range_limit)
                     )
 
         self.investment_costs = Expression(expr=investment_costs)