Reliability Prototype

h2integrate/converters/natural_gas/natural_gas_cc_ct.py

@@ -3,6 +3,7 @@
 
 from h2integrate.core.utilities import BaseConfig, merge_shared_inputs
 from h2integrate.core.validators import gt_zero, gte_zero
+from h2integrate.reliability.models import WeibullReliabilityModel
 from h2integrate.core.model_baseclasses import (
     CostModelBaseClass,
     CostModelBaseConfig,
@@ -67,6 +68,8 @@ def initialize(self):
         self.commodity = "electricity"
         self.commodity_rate_units = "MW"
         self.commodity_amount_units = "MW*h"
+        self.reliability_model = None
+        self.use_reliability = False
 
     def setup(self):
         super().setup()
@@ -75,6 +78,12 @@ def setup(self):
             merge_shared_inputs(self.options["tech_config"]["model_inputs"], "performance"),
             additional_cls_name=self.__class__.__name__,
         )
+        if self.options["tech_config"]["model_inputs"]["reliability"]:
+            self.reliability_model = WeibullReliabilityModel.from_dict(
+                merge_shared_inputs(self.options["tech_config"]["model_inputs"], "reliability"),
+                additional_cls_name=self.__class__.__name__,
+            )
+            self.use_reliability = True
 
         # Add natural gas consumed output
         self.add_output(
@@ -154,6 +163,8 @@ def compute(self, inputs, outputs):
             inputs["electricity_command_value"],
         )
         natural_gas_demand = electricity_command_value * heat_rate_mmbtu_per_mwh
+        if self.use_reliability:
+            natural_gas_demand * self.reliability_model.availability
 
         # available feedstock, saturated at maximum system feedstock consumption
         natural_gas_available = np.where(

h2integrate/reliability/models.py

@@ -0,0 +1,206 @@
+from abc import ABC, abstractmethod
+from typing import Any
+
+import numpy as np
+from attrs import field, define, validators
+
+from h2integrate.core.utilities import BaseConfig
+
+
+# generated from np.random.SeedSequence().entropy
+rng = np.random.default_rng(279299947538423226929715083173412195503)
+
+
+N_TIMESTEPS = 8760
+
+
+def create_failure_model(config: dict):
+    """Retrieves and initializes a matching reliability model."""
+    name = config.pop("failure_model")
+    match name:
+        case "WeibullReliability":
+            return WeibullReliability.from_dict(config)
+        case "FixedIntervalReliability":
+            return FixedIntervalReliability.from_dict(config)
+        case _:
+            raise NotImplementedError(f"{name} is not a valid model name")
+
+
+def create_maintenance_model(config: dict):
+    """Retrieves and initializes a matching reliability model."""
+    name = config.pop("maintenance_model")
+    match name:
+        case "WeibullReliability":
+            return WeibullReliability.from_dict(config)
+        case "FixedIntervalReliability":
+            return FixedIntervalReliability.from_dict(config)
+        case _:
+            raise NotImplementedError(f"{name} is not a valid model name")
+
+
+def generate_downtime_model(config: dict):
+    name = config.pop("model")
+    match name:
+        case "LogNormalDowntime":
+            return LogNormalDowntime.from_dict(config)
+        case _:
+            raise NotImplementedError(f"{name} is not a valid model name")
+
+
+@define(kw_only=True)
+class BaseDowntime(ABC, BaseConfig):
+    @abstractmethod
+    def sample_downtime(self) -> np.ndarray: ...
+
+
+@define(kw_only=True)
+class BaseReliability(ABC, BaseConfig):
+    @abstractmethod
+    def sample_events(self) -> np.ndarray: ...
+
+
+@define(kw_only=True)
+class LogNormalDowntime(BaseDowntime):
+    """Basic log-normal downtime model for generating the length of downtime for a given event.
+
+    Args:
+        mean (float): Average length of downtime per event, in hours.
+        sigma (float): Standard deviation of the distribution(s), in hours.
+        n_components (int): Number of identical components to sample. Primarily for convenience.
+            Defaults to 1.
+    """
+
+    mean: float = field(validator=(validators.instance_of(float), validators.ge(0)))
+    sigma: float = field(validator=(validators.instance_of(float), validators.ge(0)))
+    n_components: int = field(default=1, validator=(validators.instance_of(int), validators.ge(1)))
+
+    def sample_downtime(self) -> np.ndarray:
+        size = (self.n_components, 100) if isinstance(self.mean, int | float) else 100
+        return rng.lognormal(self.mean, self.sigma, size=size)
+
+
+@define(kw_only=True)
+class WeibullReliability(BaseReliability):
+    r"""Basic reliability model for operating/not operating statuses.
+
+    Assumes a full operational shutdown with zero ramping of production for an hourly, 1 year
+    simulation.
+
+    Args:
+        scale (float): Also referred to as :math:`\lambda` or :math:`\alpha`. Determines
+            the scale of distribution, and is equivalent to the mean time
+            between failure in years (MTBF), or 1 / annual failure rate.
+        shape (float): Also referred to as ``k`` or :math:`\beta`. A value less than 1
+            corresponds to a decreasing hazard rate over time (break-in period failures);
+            a value greater than 1 corresponds to an increasing hazard rate over time (
+            aging/wear-out failures); and a value of 1 corresponds to a constant hazard
+            rate over time (exponential distribution).
+        downtime (float): Average amount of downtime per failure.
+
+    Attributes:
+        rng (np.random._generator.Generator): NumPy random generator object.
+
+    TODO:
+        - how to pass n_timesteps through from plant?
+        - stabilize random generator/determine how to manage random seeding across library
+    """
+
+    scale: float = field(validator=validators.instance_of(float))
+    shape: float = field(validator=validators.instance_of(float))
+    downtime: Any = field(converter=generate_downtime_model)
+    downtime_per_event: np.ndarray = field(init=False, validator=validators.instance_of(np.ndarray))
+    availability: np.ndarray = field(
+        default=np.ones(N_TIMESTEPS), init=False, validator=validators.instance_of(np.ndarray)
+    )
+
+    def __attrs_post_init__(self):
+        self.create_downtime_events()
+        self.calculate_availability()
+
+    def sample_events(self):
+        return np.ceil(self.rng.weibul(self.shape, size=100) * self.scale * N_TIMESTEPS).astype(int)
+
+    def create_downtime_events(self):
+        """Creates a ``time_to_failure`` and ``downtime_per_event``."""
+        # NOTE: Arrays are default length 30 to ensure enough events are created for a 1-year
+        # simulation without burdening the memory usage.
+        self.time_to_failures = self.sample_events()
+        self.downtime_per_event = self.downtime.sample_downtime()
+
+    def calculate_availability(self):
+        """Determine the timing and duration of outages for a single year of simulation time."""
+        accumulated = 0
+        while accumulated < N_TIMESTEPS:
+            event, self.time_to_failures = self.time_to_failures[0], self.time_to_failures[1:]
+            duration, self.downtime_per_event = (
+                self.downtime_per_event[0],
+                self.downtime_per_event[1:],
+            )
+            if event + accumulated > N_TIMESTEPS:
+                break
+
+            start = accumulated + event
+            end = start + duration
+            self.availability[start:end] = 0
+            accumulated = end
+            if not self.time_to_failures:
+                self.create_downtime_events()
+
+
+@define(kw_only=True)
+class FixedIntervalReliability(BaseReliability):
+    """Basic fixed interval downtime reliability model.
+
+    Args:
+        frequency (int | float | array-like): The annual frequency of events, e.g., 4 is equivalent
+            to a quarterly downtime event and 0.25 is equivalent to an every 4 years downtime event.
+            For all events the timing of the first event will be sampled within the first year or
+            interval period to offset events from being based on January 1st in an 8760.
+        downtime (int | float | dict): Either fixed length of each downtime, in hours, or a
+            configuration dictionary for a downtime length model.
+
+    """
+
+    frequency: float = field(validator=(validators.instance_of((float, int)), validators.gt(0)))
+    downtime: Any = field(converter=generate_downtime_model)
+    downtime_per_event: np.ndarray = field(init=False, validator=validators.instance_of(np.ndarray))
+    availability: np.ndarray = field(
+        default=np.ones(N_TIMESTEPS), init=False, validator=validators.instance_of(np.ndarray)
+    )
+
+    def __attrs_post_init__(self):
+        self.create_downtime_events()
+        self.calculate_availability()
+
+    def sample_events(self):
+        frequency = np.array([0.1, 20, 30]).reshape(-1, 1)
+        interval = np.ceil(8760 / frequency).astype(int)
+        first_occurrence = rng.integers(0, np.where(interval > 8760, 8760, interval))
+        return first_occurrence, interval
+
+    def create_downtime_events(self):
+        """Creates a ``time_to_failure`` and ``downtime_per_event``."""
+        self.first_occurrence, self.interval = self.sample_events()
+        self.downtime_per_event = self.downtime.sample_downtime()
+
+    def calculate_availability(self):
+        """Determine the timing and duration of outages for a single year of simulation time."""
+        duration, self.downtime_per_event = (
+            self.downtime_per_event[0],
+            self.downtime_per_event[1:],
+        )
+        end = self.first_occurrence + duration
+        self.availability[self.first_occurrence : end] = 0
+        accumulated = end
+        while (start := accumulated + self.interval) < N_TIMESTEPS:
+            duration, self.downtime_per_event = (
+                self.downtime_per_event[0],
+                self.downtime_per_event[1:],
+            )
+            # TODO: handle start > 8760
+            end = start + duration
+            end = np.where(end > 8760, 8760, end)
+            self.availability[start:end] = 0
+            accumulated = end
+            if not self.downtime_per_event:
+                self.downtime.sample_downtime()