Allow for fixing investments

docs/usage.rst

@@ -1279,6 +1279,33 @@ Besides the `invest` variable, new variables are introduced as well. These are:
  monthly periods, but you would need to be very careful in parameterizing your energy system and your model and also,
  this would mean monthly discounting (if applicable) as well as specifying your plants lifetimes in months.
 
+Repeated solving with fixed investment results
+----------------------------------------------
+
+You can rerun a given optimization model with fixed investments.
+Let's assume you have set up a model instance called ``model``
+containing investments and solved it by calling:
+
+.. code-block:: python
+
+ model = solph.Model(energysystem)
+ model.solve()
+
+You can now fix the investment results from the previous solving and
+only solve the dispatch-related part via a call of the method
+``fix_investments`` as well as a second solving of your model:
+
+.. code-block:: python
+
+ model.fix_investments()
+ model.solve()
+
+This will now take the investment variables as parameters. Thus, their
+values from the previous solving are parameters in the second solving. Only
+the dispatch-related variables will be optimized. This two-stage approach
+for instance may allow you to interpret dual values of balancing constraints
+which you cannot sensibly interpret in the presence of investment variables.
+
 Modelling cellular energy systems and modularizing energy system models
 -----------------------------------------------------------------------
 

docs/whatsnew/v0-5-2.rst

@@ -5,13 +5,16 @@ API changes
 ###########
 
 * New bool attribute `use_remaining_value` of `oemof.solph.EnergySystem`
+* New method `fix_investments` of `oemof.solph.Model`
 
 New features
 ############
 
 * Allow for evaluating differences in the remaining vs. the original value
  for multi-period investments.
 * Allow to define minimum up- and down-time per time step
+* Allow for fixing investment variables to the values of a previous
+ optimization run, enabling a repeated solve with fixed investments.
 
 Documentation
 #############

src/oemof/solph/_models.py

@@ -317,6 +317,9 @@ class Model(BaseModel):
  discount_rate : float or None
  The rate used for discounting in a multi-period model.
  A 2% discount rate needs to be defined as 0.02.
+ rounding_precision : int or None
+ Number of decimal digits to apply for rounding when fixing variables
+ in a repeated model solve.
 
  Note
  ----
@@ -383,6 +386,8 @@ def __init__(self, energysystem, discount_rate=None, **kwargs):
  self._set_discount_rate_with_warning()
  else:
  pass
+ self._fix_investments = False
+ self.rounding_precision = None
  super().__init__(energysystem, **kwargs)
 
  def _set_discount_rate_with_warning(self):
@@ -520,3 +525,31 @@ def _add_parent_block_variables(self):
  if (o, i) in self.UNIDIRECTIONAL_FLOWS:
  for p, t in self.TIMEINDEX:
  self.flow[o, i, p, t].setlb(0)
+
+ def fix_investments(self, rounding_precision=None):
+ """Fix investment results of an already solved model
+
+ Parameters
+ ----------
+ rounding_precision : int or None
+ If not None, round investments to given number of decimal digits
+ """
+ if self.solver_results is None:
+ msg = (
+ "Cannot fix investments as model has not yet been solved!\n"
+ "You have to first solve your model and then call method "
+ "`fix_investments()` on your model instance."
+ )
+ raise ValueError(msg)
+ self._fix_investments = True
+ self.rounding_precision = rounding_precision
+ if hasattr(self, "InvestmentFlowBlock"):
+ self.InvestmentFlowBlock.fix_investments_results()
+ if hasattr(self, "GenericInvestmentStorageBlock"):
+ self.GenericInvestmentStorageBlock.fix_investments_results()
+ if hasattr(self, "SinkDSMOemofInvestmentBlock"):
+ self.SinkDSMOemofInvestmentBlock.fix_investments_results()
+ if hasattr(self, "SinkDSMDIWInvestmentBlock"):
+ self.SinkDSMDIWInvestmentBlock.fix_investments_results()
+ if hasattr(self, "SinkDSMDLRInvestmentBlock"):
+ self.SinkDSMDLRInvestmentBlock.fix_investments_results()

src/oemof/solph/components/_generic_storage.py

@@ -35,6 +35,7 @@
 from pyomo.environ import NonNegativeReals
 from pyomo.environ import Set
 from pyomo.environ import Var
+from pyomo.environ import value
 
 from oemof.solph._helpers import check_node_object_for_missing_attribute
 from oemof.solph._options import Investment
@@ -1968,3 +1969,66 @@ def _evaluate_remaining_value_difference(
  return 0
  else:
  return 0
+
+ def fix_investments_results(self):
+ """Fix investments if `_fix_investments` is set to True for model"""
+ m = self.parent_block()
+ for n in self.INVESTSTORAGES:
+ for p in m.PERIODS:
+ if m.rounding_precision is not None:
+ self.invest[n, p] = self._round_and_ensure_bounds(
+ self.invest[n, p],
+ lb=n.investment.minimum[p],
+ ub=n.investment.maximum[p],
+ )
+ ub_total = (
+ n.investment.overall_maximum
+ if n.investment.overall_maximum is not None
+ else float("+inf")
+ )
+ self.total[n, p] = self._round_and_ensure_bounds(
+ self.total[n, p],
+ ub=ub_total,
+ )
+ self.invest[n, p].fix()
+ self.total[n, p].fix()
+ if m.es.periods is not None:
+ if m.rounding_precision is not None:
+ self.old[n, p] = self._round_and_ensure_bounds(
+ self.old[n, p],
+ )
+ self.old_end[n, p] = self._round_and_ensure_bounds(
+ self.old_end[n, p],
+ )
+ self.old_exo[n, p] = self._round_and_ensure_bounds(
+ self.old_exo[n, p],
+ )
+ self.old[n, p].fix()
+ self.old_end[n, p].fix()
+ self.old_exo[n, p].fix()
+
+ def _round_and_ensure_bounds(self, var, lb=0, ub=float("+inf")):
+ """Round given investments and ensure within bounds
+
+ Parameters
+ ----------
+ var : pyomo.core.IndexedVar
+ variable to round
+
+ lb : float
+ Lower bound to be ensured
+
+ ub : float
+ Upper bound to be ensured
+ """
+ m = self.parent_block()
+ return min(
+ ub,
+ max(
+ lb,
+ round(
+ value(var),
+ m.rounding_precision,
+ ),
+ ),
+ )

src/oemof/solph/components/experimental/_sink_dsm.py

@@ -35,6 +35,7 @@
 from pyomo.environ import NonNegativeReals
 from pyomo.environ import Set
 from pyomo.environ import Var
+from pyomo.environ import value
 
 from oemof.solph._options import Investment
 from oemof.solph._plumbing import sequence
@@ -1536,6 +1537,69 @@ def _evaluate_remaining_value_difference(
  else:
  return 0
 
+ def fix_investments_results(self):
+ """Fix investments if `_fix_investments` is set to True for model"""
+ m = self.parent_block()
+ for g in self.investdsm:
+ for p in m.PERIODS:
+ if m.rounding_precision is not None:
+ self.invest[g, p] = self._round_and_ensure_bounds(
+ self.invest[g, p],
+ lb=g.investment.minimum[p],
+ ub=g.investment.maximum[p],
+ )
+ ub_total = (
+ g.investment.overall_maximum
+ if g.investment.overall_maximum is not None
+ else float("+inf")
+ )
+ self.total[g, p] = self._round_and_ensure_bounds(
+ self.total[g, p],
+ ub=ub_total,
+ )
+ self.invest[g, p].fix()
+ self.total[g, p].fix()
+ if m.es.periods is not None:
+ if m.rounding_precision is not None:
+ self.old[g, p] = self._round_and_ensure_bounds(
+ self.old[g, p],
+ )
+ self.old_end[g, p] = self._round_and_ensure_bounds(
+ self.old_end[g, p],
+ )
+ self.old_exo[g, p] = self._round_and_ensure_bounds(
+ self.old_exo[g, p],
+ )
+ self.old[g, p].fix()
+ self.old_end[g, p].fix()
+ self.old_exo[g, p].fix()
+
+ def _round_and_ensure_bounds(self, var, lb=0, ub=float("+inf")):
+ """Round given investments and ensure within bounds
+
+ Parameters
+ ----------
+ var : pyomo.core.IndexedVar
+ variable to round
+
+ lb : float
+ Lower bound to be ensured
+
+ ub : float
+ Upper bound to be ensured
+ """
+ m = self.parent_block()
+ return min(
+ ub,
+ max(
+ lb,
+ round(
+ value(var),
+ m.rounding_precision,
+ ),
+ ),
+ )
+
 
 class SinkDSMDIWBlock(ScalarBlock):
  r"""Constraints for SinkDSM with "DIW" approach
@@ -3390,6 +3454,69 @@ def _evaluate_remaining_value_difference(
  else:
  return 0
 
+ def fix_investments_results(self):
+ """Fix investments if `_fix_investments` is set to True for model"""
+ m = self.parent_block()
+ for g in self.investdsm:
+ for p in m.PERIODS:
+ if m.rounding_precision is not None:
+ self.invest[g, p] = self._round_and_ensure_bounds(
+ self.invest[g, p],
+ lb=g.investment.minimum[p],
+ ub=g.investment.maximum[p],
+ )
+ ub_total = (
+ g.investment.overall_maximum
+ if g.investment.overall_maximum is not None
+ else float("+inf")
+ )
+ self.total[g, p] = self._round_and_ensure_bounds(
+ self.total[g, p],
+ ub=ub_total,
+ )
+ self.invest[g, p].fix()
+ self.total[g, p].fix()
+ if m.es.periods is not None:
+ if m.rounding_precision is not None:
+ self.old[g, p] = self._round_and_ensure_bounds(
+ self.old[g, p],
+ )
+ self.old_end[g, p] = self._round_and_ensure_bounds(
+ self.old_end[g, p],
+ )
+ self.old_exo[g, p] = self._round_and_ensure_bounds(
+ self.old_exo[g, p],
+ )
+ self.old[g, p].fix()
+ self.old_end[g, p].fix()
+ self.old_exo[g, p].fix()
+
+ def _round_and_ensure_bounds(self, var, lb=0, ub=float("+inf")):
+ """Round given investments and ensure within bounds
+
+ Parameters
+ ----------
+ var : pyomo.core.IndexedVar
+ variable to round
+
+ lb : float
+ Lower bound to be ensured
+
+ ub : float
+ Upper bound to be ensured
+ """
+ m = self.parent_block()
+ return min(
+ ub,
+ max(
+ lb,
+ round(
+ value(var),
+ m.rounding_precision,
+ ),
+ ),
+ )
+
 
 class SinkDSMDLRBlock(ScalarBlock):
  r"""Constraints for SinkDSM with "DLR" approach
@@ -5888,3 +6015,66 @@ def _evaluate_remaining_value_difference(
  return 0
  else:
  return 0
+
+ def fix_investments_results(self):
+ """Fix investments if `_fix_investments` is set to True for model"""
+ m = self.parent_block()
+ for g in self.INVESTDR:
+ for p in m.PERIODS:
+ if m.rounding_precision is not None:
+ self.invest[g, p] = self._round_and_ensure_bounds(
+ self.invest[g, p],
+ lb=g.investment.minimum[p],
+ ub=g.investment.maximum[p],
+ )
+ ub_total = (
+ g.investment.overall_maximum
+ if g.investment.overall_maximum is not None
+ else float("+inf")
+ )
+ self.total[g, p] = self._round_and_ensure_bounds(
+ self.total[g, p],
+ ub=ub_total,
+ )
+ self.invest[g, p].fix()
+ self.total[g, p].fix()
+ if m.es.periods is not None:
+ if m.rounding_precision is not None:
+ self.old[g, p] = self._round_and_ensure_bounds(
+ self.old[g, p],
+ )
+ self.old_end[g, p] = self._round_and_ensure_bounds(
+ self.old_end[g, p],
+ )
+ self.old_exo[g, p] = self._round_and_ensure_bounds(
+ self.old_exo[g, p],
+ )
+ self.old[g, p].fix()
+ self.old_end[g, p].fix()
+ self.old_exo[g, p].fix()
+
+ def _round_and_ensure_bounds(self, var, lb=0, ub=float("+inf")):
+ """Round given investments and ensure within bounds
+
+ Parameters
+ ----------
+ var : pyomo.core.IndexedVar
+ variable to round
+
+ lb : float
+ Lower bound to be ensured
+
+ ub : float
+ Upper bound to be ensured
+ """
+ m = self.parent_block()
+ return min(
+ ub,
+ max(
+ lb,
+ round(
+ value(var),
+ m.rounding_precision,
+ ),
+ ),
+ )