Logical flow

osier/__init__.py

@@ -10,7 +10,9 @@
 
 
 from .technology import *
+from .models.model import *
 from .models.dispatch import *
+from .models.logic_dispatch import *
 from .models.capacity_expansion import *
 from .models.deap_runner import *
 from .utils import *

osier/models/capacity_expansion.py

@@ -7,12 +7,17 @@
 import functools
 import time
 
-from osier import DispatchModel
+from osier import DispatchModel, LogicDispatchModel
 
 from pymoo.core.problem import ElementwiseProblem
 
 LARGE_NUMBER = 1e20
 
+dispatch_models = {'lp': DispatchModel,
+                   'hierarchical': LogicDispatchModel,
+                   'logical': LogicDispatchModel,
+                   'rule_based': LogicDispatchModel}
+
 
 class CapacityExpansion(ElementwiseProblem):
     """
@@ -98,6 +103,7 @@ def __init__(self,
                  allow_blackout=False,
                  verbosity=50,
                  solver='cbc',
+                 model_engine='lp',
                  **kwargs):
         self.technology_list = deepcopy(technology_list)
         self.demand = demand
@@ -109,6 +115,7 @@ def __init__(self,
         self.curtailment = curtailment
         self.allow_blackout = allow_blackout
         self.verbosity = verbosity
+        self.model_engine = model_engine.lower()
         self.solver = solver
 
         if isinstance(demand, unyt_array):
@@ -186,13 +193,14 @@ def _evaluate(self, x, out, *args, **kwargs):
             - wind_gen \
             - solar_gen
 
-        model = DispatchModel(technology_list=self.dispatchable_techs,
-                              net_demand=net_demand,
-                              power_units=self.power_units,
-                              curtailment=self.curtailment,
-                              allow_blackout=self.allow_blackout,
-                              solver=self.solver,
-                              verbosity=self.verbosity)
+        dispatch_model = dispatch_models[self.model_engine]
+        model = dispatch_model(technology_list=self.dispatchable_techs,
+                               net_demand=net_demand,
+                               power_units=self.power_units,
+                               curtailment=self.curtailment,
+                               allow_blackout=self.allow_blackout,
+                               solver=self.solver,
+                               verbosity=self.verbosity)
         model.solve()
 
         if model.results is not None:

osier/models/dispatch.py

@@ -108,7 +108,7 @@ class DispatchModel():
         Can be overridden by specifying a unit with the value.
     solver : str
         Indicates which solver to use. May require separate installation.
-        Accepts: ['cplex', 'cbc']. Other solvers will be added in the future.
+        Accepts: ['cplex', 'cbc', 'appsi_highs']. Other solvers will be added in the future.
     lower_bound : float
         The minimum amount of energy each technology can produce per time
         period. Default is 0.0.
@@ -598,10 +598,10 @@ def solve(self, solver=None):
         else:
             optimizer = po.SolverFactory(self.solver)
 
-        results = optimizer.solve(self.model, tee=self.verbose)
         try:
+            optimizer.solve(self.model, tee=self.verbose)
             self.objective = self.model.objective()
-        except ValueError:
+        except (ValueError, RuntimeError):
             if self.verbosity <= 30:
                 warnings.warn(
                     f"Infeasible or no solution. Objective set to {LARGE_NUMBER}")

osier/models/logic_dispatch.py

@@ -0,0 +1,100 @@
+from osier import OsierModel
+from osier.utils import get_tech_names
+from copy import deepcopy
+from unyt import unyt_array, MW
+import pandas as pd
+import numpy as np
+import warnings
+
+LARGE_NUMBER = 1e20
+
+
+class LogicDispatchModel(OsierModel):
+
+    def __init__(self,
+                 technology_list,
+                 net_demand,
+                 allow_blackout=False,
+                 curtailment=True,
+                 verbosity=50,
+                 *args, **kwargs):
+        super().__init__(technology_list=technology_list,
+                         net_demand=net_demand,
+                         *args, **kwargs)
+        self.technology_list = technology_list
+        self.technology_list.sort()
+        self.original_order = get_tech_names(self.technology_list)
+        self.cost_history = np.zeros(len(net_demand))
+        self.covered_demand = None
+        self.objective = None
+        self.results = None
+        self.verbosity = verbosity
+        self.allow_blackout = allow_blackout
+        self.curtailment = curtailment
+
+    def _reset_all(self):
+        for t in self.technology_list:
+            t.reset_history()
+        return
+
+    def _format_results(self):
+        data = {}
+        for t in self.technology_list:
+            data[f"{t.technology_name}"] = unyt_array(
+                t.power_history).to_ndarray()
+            if t.technology_type == 'storage':
+                data[f"{t.technology_name}_level"] = unyt_array(
+                    t.storage_history).to_ndarray()
+                data[f"{t.technology_name}_charge"] = unyt_array(
+                    t.charge_history).to_ndarray()
+        data["Curtailment"] = np.array(
+            [v if v <= 0 else 0 for v in self.covered_demand])
+        data["Shortfall"] = np.array(
+            [v if v > 0 else 0 for v in self.covered_demand])
+        self.results = pd.DataFrame(data)
+        return
+
+    def _calculate_objective(self):
+        self.objective = sum(np.array(t.power_history).sum()
+                             * t.variable_cost.to_value()
+                             for t in self.technology_list)
+        return
+
+    def solve(self):
+        """
+        This function executes the model solve with a rule-based approach.
+        Net demand is copied, then the technology histories are reset.
+        """
+        self.covered_demand = self.net_demand.copy()
+        self._reset_all()
+        try:
+            for i, v in enumerate(self.covered_demand):
+                for t in self.technology_list:
+                    power_out = t.power_output(v, time_delta=self.time_delta)
+                    v -= power_out
+
+                self.covered_demand[i] = v
+                if not self.allow_blackout and (v > 0):
+                    if self.verbosity <= 20:
+                        print('solve failed -- unmet demand')
+                    raise ValueError
+
+                if not self.curtailment and (v < 0):
+                    if self.verbosity <= 20:
+                        print(
+                            ('solve failed -- '
+                            'too much overproduction '
+                            '(no curtailment allowed)'))
+                    raise ValueError
+
+            self._format_results()
+            self._calculate_objective()
+        except ValueError:
+            if self.verbosity <= 30:
+                warnings.warn(
+                    (f"Infeasible or no solution." 
+                     f"Objective set to {LARGE_NUMBER}")
+                    )
+            self.objective = LARGE_NUMBER
+
+        return

osier/models/model.py

@@ -0,0 +1,80 @@
+import pandas as pd
+import numpy as np
+from unyt import unyt_array, hr, MW, kW, GW
+import itertools as it
+from osier import Technology
+from osier.technology import _validate_quantity, _validate_unit
+from osier.utils import synchronize_units
+import warnings
+import logging
+
+_freq_opts = {'D': 'day',
+              'H': 'hour',
+              'S': 'second',
+              'T': 'minute'}
+
+
+class OsierModel():
+    """
+    A class for instantiating energy models in Osier.
+    """
+
+    def __init__(self,
+                 technology_list,
+                 net_demand,
+                 time_delta=None,
+                 power_units=MW,
+                 **kwargs):
+        self.net_demand = net_demand
+        self.time_delta = time_delta
+        self.results = None
+        self.objective = None
+
+        if isinstance(net_demand, unyt_array):
+            self.power_units = net_demand.units
+        elif isinstance(net_demand, (np.ndarray, list)):
+            self.power_units = power_units
+            self.net_demand = np.array(self.net_demand) * self.power_units
+        elif isinstance(net_demand, pd.core.series.Series):
+            self.power_units = power_units
+            self.net_demand = np.array(self.net_demand) * self.power_units
+        else:
+            self.power_units = power_units
+
+        self.technology_list = synchronize_units(
+            technology_list,
+            unit_power=self.power_units,
+            unit_time=self.time_delta.units)
+
+    @property
+    def time_delta(self):
+        return self._time_delta
+
+    @time_delta.setter
+    def time_delta(self, value):
+        if value:
+            valid_quantity = _validate_quantity(value, dimension='time')
+            self._time_delta = valid_quantity
+        else:
+            if isinstance(self.net_demand, pd.DataFrame):
+                try:
+                    freq_list = list(self.net_demand.index.inferred_freq)
+                    freq_key = freq_list[-1]
+                    try:
+                        value = float(freq_list[0])
+                    except ValueError:
+                        warnings.warn((f"Could not convert value "
+                                       f"{freq_list[0]} to float. "
+                                       "Setting to 1.0."),
+                                      UserWarning)
+                        value = 1.0
+                    self._time_delta = _validate_quantity(
+                        f"{value} {_freq_opts[freq_key]}", dimension='time')
+                except KeyError:
+                    warnings.warn(
+                        (f"Could not infer time delta with freq {freq_key} "
+                         "from pandas dataframe. Setting delta to 1 hour."),
+                        UserWarning)
+                    self._time_delta = 1 * hr
+            else:
+                self._time_delta = 1 * hr

osier/tech_library.py

@@ -96,7 +96,7 @@
                             fuel_cost=0*to_MDOLLARS,
                             storage_duration=4,
                             efficiency=0.85,
-                            initial_storage=0,
+                            initial_storage=0.0*MW*hr,
                             capacity_credit=0.5,
                             lifecycle_co2_rate=3.3e-5*co2_eq_units,
                             land_intensity=0.0,