From f6b1c1814b7724897178d4e3058f82dbf17bae57 Mon Sep 17 00:00:00 2001
From: "Hammond, Rob" <13874373+RHammond2@users.noreply.github.com>
Date: Mon, 10 Aug 2026 14:19:47 -0700
Subject: [PATCH 01/24] add initial prototype for weibull reliability

---
 h2integrate/core/model_baseclasses.py |  52 +++++++++++++
 h2integrate/reliability/__init__.py   |   0
 h2integrate/reliability/models.py     | 105 ++++++++++++++++++++++++++
 3 files changed, 157 insertions(+)
 create mode 100644 h2integrate/reliability/__init__.py
 create mode 100644 h2integrate/reliability/models.py

diff --git a/h2integrate/core/model_baseclasses.py b/h2integrate/core/model_baseclasses.py
index 4f1160f53..c6f9e2f4e 100644
--- a/h2integrate/core/model_baseclasses.py
+++ b/h2integrate/core/model_baseclasses.py
@@ -150,6 +150,25 @@ def compute(self, inputs, outputs, discrete_inputs, discrete_outputs):
         raise NotImplementedError("This method should be implemented in a subclass.")
 
 
+class ReliabilityModelBaseClass(om.ExplicitComponent):
+    def initialize(self):
+        self.options.declare("driver_config", types=dict)
+        self.options.declare("plant_config", types=dict)
+        self.options.declare("tech_config", types=dict)
+
+    def setup(self):
+        self.n_timesteps = self.options["plant_config"]["plant"]["simulation"]["n_timesteps"]
+        self.availability = np.ones(self.n_timesteps, dtype=float)
+
+        self.add_output(
+            f"{self.commodity}_availability",
+            val=self.availability,
+            shape=self.n_timesteps,
+            units="unitless",
+            desc="Production-based availability for a commodity.",
+        )
+
+
 @define(kw_only=True)
 class CostModelBaseConfig(BaseConfig):
     cost_year: int = field(converter=int)
@@ -545,3 +564,36 @@ def compute(self, inputs, outputs, discrete_inputs, discrete_outputs):
         # self.cache_outputs(inputs, outputs, discrete_inputs, discrete_outputs)
 
         raise NotImplementedError("This method should be implemented in a subclass.")
+
+    def apply_curtailment(self, outputs):
+        """Apply curtailment to ``{commodity}_out`` based on ``{commodity}_command_value``.
+
+        Copies the current ``{commodity}_out`` into ``uncurtailed_{commodity}_out``,
+        then clips ``{commodity}_out`` to ``min(uncurtailed, command_value)`` element-wise.
+
+        Only operates when the model has ``_control_classifier == "flexible"``.
+        Should be called at the end of each flexible model's ``compute()`` method
+        after the raw production has been written to ``outputs[f"{commodity}_out"]``.
+        """
+        if "system_level_control" in self.options["plant_config"]:
+            if getattr(self, "_control_classifier", None) != "flexible":
+                return
+
+            commodity_out_key = f"{self.commodity}_out"
+            uncurtailed_key = f"uncurtailed_{self.commodity}_out"
+            command_value_key = f"{self.commodity}_command_value"
+
+            uncurtailed = np.array(outputs[commodity_out_key])
+            outputs[uncurtailed_key] = uncurtailed
+
+            command_value = self._inputs[command_value_key]
+            outputs[commodity_out_key] = np.minimum(uncurtailed, command_value)
+
+    def compute(self, inputs, outputs, discrete_inputs, discrete_outputs):  # noqa: F811
+        """
+        Computation for the OM component.
+
+        For a template class this is not implement and raises an error.
+        """
+
+        raise NotImplementedError("This method should be implemented in a subclass.")
diff --git a/h2integrate/reliability/__init__.py b/h2integrate/reliability/__init__.py
new file mode 100644
index 000000000..e69de29bb
diff --git a/h2integrate/reliability/models.py b/h2integrate/reliability/models.py
new file mode 100644
index 000000000..e2d2d3ee8
--- /dev/null
+++ b/h2integrate/reliability/models.py
@@ -0,0 +1,105 @@
+import numpy as np
+from attrs import field, define, validators
+
+from h2integrate.core.utilities import BaseConfig, merge_shared_inputs
+from h2integrate.core.model_baseclasses import ReliabilityModelBaseClass
+
+
+@define(kw_only=True)
+class WeibullReliabilityConfig(BaseConfig):
+    r"""Basic reliability model for operating/not operating statuses.
+
+    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:
+        - 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: float = field(validator=validators.gt(1))
+    rng: np.random._generator.Generator = field(
+        default=np.random.default_rng(),
+        init=False,
+        validator=validators.instance_of(np.random._generator.Generator),
+    )
+
+
+class WeibullReliabilityModel(ReliabilityModelBaseClass):
+    """
+    Performance model for natural gas power plants.
+
+    This model calculates electricity output from natural gas input based on
+    the plant's heat rate. It can be used for both natural gas combustion
+    turbines (NGCT) and natural gas combined cycle (NGCC) plants by providing
+    appropriate heat rate values.
+
+    The model implements the relationship:
+        electricity_out = natural_gas_in / heat_rate
+
+    Inputs:
+        system_capacity (float): Natural gas plant rated capacity in MW
+        natural_gas_in (array): Natural gas input energy in MMBtu/h
+        heat_rate_mmbtu_per_mwh (float): Plant heat rate in MMBtu/MWh
+        electricity_command_value (array): Electricity command value in MW for each timestep
+
+    Outputs:
+        electricity_out (array): Electricity output in MW for each timestep
+        natural_gas_consumed (array): Natural gas consumed in MMBtu/h
+
+    """
+
+    def initialize(self):
+        super().initialize()
+        self.commodity = "electricity"
+
+    def setup(self):
+        super().setup()
+
+        self.config = WeibullReliabilityConfig.from_dict(
+            merge_shared_inputs(self.options["tech_config"]["model_inputs"], "reliability"),
+            additional_cls_name=self.__class__.__name__,
+        )
+
+        self.create_downtime_events()
+
+    def create_downtime_events(self):
+        """Creates a ``time_to_failure`` and ``downtime_per_event`` array based on the distributions
+        described in ``WeibullReliabilityConfig``.
+        """
+        self.time_to_failures = np.floor(
+            self.config.rng.weibul(self.config.shape, size=12) * self.config.scale * 8760
+        ).astype(int)
+        downtime_per_event = np.floor(np.rng.normal(loc=self.config.downtime, size=12)).astype(int)
+        self.downtime_per_event = np.where(downtime_per_event >= 1, downtime_per_event, 1)
+
+    def calculate_availability(self):
+        """Determine the timing and duration of outages for a single year of simulation time."""
+        accumulated = 0
+        while accumulated < 8760:
+            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 > 8760:
+                break
+
+            start = accumulated + event
+            end = start + duration
+            self.availability[start:end] = 0
+            accumulated = start
+            if not self.time_to_failures:
+                self.create_downtime_events()

From 19e3e1cc4ee5b8d4698801315046267cc9535498 Mon Sep 17 00:00:00 2001
From: "Hammond, Rob" <13874373+RHammond2@users.noreply.github.com>
Date: Tue, 11 Aug 2026 09:52:36 -0700
Subject: [PATCH 02/24] simplify relibaility model non openmdao model

---
 .../natural_gas/natural_gas_cc_ct.py          | 11 +++
 h2integrate/core/model_baseclasses.py         | 19 -----
 h2integrate/reliability/models.py             | 69 +++++++------------
 3 files changed, 35 insertions(+), 64 deletions(-)

diff --git a/h2integrate/converters/natural_gas/natural_gas_cc_ct.py b/h2integrate/converters/natural_gas/natural_gas_cc_ct.py
index fde5b8590..b1d960eb5 100644
--- a/h2integrate/converters/natural_gas/natural_gas_cc_ct.py
+++ b/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(
diff --git a/h2integrate/core/model_baseclasses.py b/h2integrate/core/model_baseclasses.py
index c6f9e2f4e..9964c286f 100644
--- a/h2integrate/core/model_baseclasses.py
+++ b/h2integrate/core/model_baseclasses.py
@@ -150,25 +150,6 @@ def compute(self, inputs, outputs, discrete_inputs, discrete_outputs):
         raise NotImplementedError("This method should be implemented in a subclass.")
 
 
-class ReliabilityModelBaseClass(om.ExplicitComponent):
-    def initialize(self):
-        self.options.declare("driver_config", types=dict)
-        self.options.declare("plant_config", types=dict)
-        self.options.declare("tech_config", types=dict)
-
-    def setup(self):
-        self.n_timesteps = self.options["plant_config"]["plant"]["simulation"]["n_timesteps"]
-        self.availability = np.ones(self.n_timesteps, dtype=float)
-
-        self.add_output(
-            f"{self.commodity}_availability",
-            val=self.availability,
-            shape=self.n_timesteps,
-            units="unitless",
-            desc="Production-based availability for a commodity.",
-        )
-
-
 @define(kw_only=True)
 class CostModelBaseConfig(BaseConfig):
     cost_year: int = field(converter=int)
diff --git a/h2integrate/reliability/models.py b/h2integrate/reliability/models.py
index e2d2d3ee8..4ac973cd1 100644
--- a/h2integrate/reliability/models.py
+++ b/h2integrate/reliability/models.py
@@ -1,14 +1,19 @@
 import numpy as np
 from attrs import field, define, validators
 
-from h2integrate.core.utilities import BaseConfig, merge_shared_inputs
-from h2integrate.core.model_baseclasses import ReliabilityModelBaseClass
+from h2integrate.core.utilities import BaseConfig
+
+
+N_TIMESTEPS = 8760
 
 
 @define(kw_only=True)
-class WeibullReliabilityConfig(BaseConfig):
+class WeibullReliabilityModel(BaseConfig):
     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
@@ -24,6 +29,7 @@ class WeibullReliabilityConfig(BaseConfig):
         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
     """
 
@@ -35,66 +41,39 @@ class WeibullReliabilityConfig(BaseConfig):
         init=False,
         validator=validators.instance_of(np.random._generator.Generator),
     )
+    availability: np.ndarray = field(
+        default=np.ones(N_TIMESTEPS), init=False, validator=validators.instance_of(np.ndarray)
+    )
+    downtime_per_event: np.ndarray = field(
+        default=np.zeros(N_TIMESTEPS), init=False, validator=validators.instance_of(np.ndarray)
+    )
 
-
-class WeibullReliabilityModel(ReliabilityModelBaseClass):
-    """
-    Performance model for natural gas power plants.
-
-    This model calculates electricity output from natural gas input based on
-    the plant's heat rate. It can be used for both natural gas combustion
-    turbines (NGCT) and natural gas combined cycle (NGCC) plants by providing
-    appropriate heat rate values.
-
-    The model implements the relationship:
-        electricity_out = natural_gas_in / heat_rate
-
-    Inputs:
-        system_capacity (float): Natural gas plant rated capacity in MW
-        natural_gas_in (array): Natural gas input energy in MMBtu/h
-        heat_rate_mmbtu_per_mwh (float): Plant heat rate in MMBtu/MWh
-        electricity_command_value (array): Electricity command value in MW for each timestep
-
-    Outputs:
-        electricity_out (array): Electricity output in MW for each timestep
-        natural_gas_consumed (array): Natural gas consumed in MMBtu/h
-
-    """
-
-    def initialize(self):
-        super().initialize()
-        self.commodity = "electricity"
-
-    def setup(self):
-        super().setup()
-
-        self.config = WeibullReliabilityConfig.from_dict(
-            merge_shared_inputs(self.options["tech_config"]["model_inputs"], "reliability"),
-            additional_cls_name=self.__class__.__name__,
-        )
-
+    def __attrs_post_init__(self):
         self.create_downtime_events()
+        self.calculate_availability()
 
     def create_downtime_events(self):
         """Creates a ``time_to_failure`` and ``downtime_per_event`` array based on the distributions
         described in ``WeibullReliabilityConfig``.
         """
-        self.time_to_failures = np.floor(
-            self.config.rng.weibul(self.config.shape, size=12) * self.config.scale * 8760
+        # 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 = np.ceil(
+            self.config.rng.weibul(self.config.shape, size=30) * self.config.scale * N_TIMESTEPS
         ).astype(int)
-        downtime_per_event = np.floor(np.rng.normal(loc=self.config.downtime, size=12)).astype(int)
+        downtime_per_event = np.ceil(np.rng.normal(loc=self.config.downtime, size=30)).astype(int)
         self.downtime_per_event = np.where(downtime_per_event >= 1, downtime_per_event, 1)
 
     def calculate_availability(self):
         """Determine the timing and duration of outages for a single year of simulation time."""
         accumulated = 0
-        while accumulated < 8760:
+        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 > 8760:
+            if event + accumulated > N_TIMESTEPS:
                 break
 
             start = accumulated + event

From fccb5bd1fecbe665908cbfb62c2c2f8798dc8e22 Mon Sep 17 00:00:00 2001
From: "Hammond, Rob" <13874373+RHammond2@users.noreply.github.com>
Date: Tue, 11 Aug 2026 12:21:45 -0700
Subject: [PATCH 03/24] add model retrieval and remove extrasfrom baseclasses

---
 h2integrate/core/model_baseclasses.py | 33 ---------------------------
 h2integrate/reliability/models.py     | 10 ++++++++
 2 files changed, 10 insertions(+), 33 deletions(-)

diff --git a/h2integrate/core/model_baseclasses.py b/h2integrate/core/model_baseclasses.py
index 9964c286f..4f1160f53 100644
--- a/h2integrate/core/model_baseclasses.py
+++ b/h2integrate/core/model_baseclasses.py
@@ -545,36 +545,3 @@ def compute(self, inputs, outputs, discrete_inputs, discrete_outputs):
         # self.cache_outputs(inputs, outputs, discrete_inputs, discrete_outputs)
 
         raise NotImplementedError("This method should be implemented in a subclass.")
-
-    def apply_curtailment(self, outputs):
-        """Apply curtailment to ``{commodity}_out`` based on ``{commodity}_command_value``.
-
-        Copies the current ``{commodity}_out`` into ``uncurtailed_{commodity}_out``,
-        then clips ``{commodity}_out`` to ``min(uncurtailed, command_value)`` element-wise.
-
-        Only operates when the model has ``_control_classifier == "flexible"``.
-        Should be called at the end of each flexible model's ``compute()`` method
-        after the raw production has been written to ``outputs[f"{commodity}_out"]``.
-        """
-        if "system_level_control" in self.options["plant_config"]:
-            if getattr(self, "_control_classifier", None) != "flexible":
-                return
-
-            commodity_out_key = f"{self.commodity}_out"
-            uncurtailed_key = f"uncurtailed_{self.commodity}_out"
-            command_value_key = f"{self.commodity}_command_value"
-
-            uncurtailed = np.array(outputs[commodity_out_key])
-            outputs[uncurtailed_key] = uncurtailed
-
-            command_value = self._inputs[command_value_key]
-            outputs[commodity_out_key] = np.minimum(uncurtailed, command_value)
-
-    def compute(self, inputs, outputs, discrete_inputs, discrete_outputs):  # noqa: F811
-        """
-        Computation for the OM component.
-
-        For a template class this is not implement and raises an error.
-        """
-
-        raise NotImplementedError("This method should be implemented in a subclass.")
diff --git a/h2integrate/reliability/models.py b/h2integrate/reliability/models.py
index 4ac973cd1..cf8ca4dba 100644
--- a/h2integrate/reliability/models.py
+++ b/h2integrate/reliability/models.py
@@ -7,6 +7,16 @@
 N_TIMESTEPS = 8760
 
 
+def create_reliability(config: dict):
+    """Retrieves and initializes a matching reliability model."""
+    name = config.pop("reliability")
+    match name:
+        case "WeibullReliabilityModel":
+            return WeibullReliabilityModel.from_dict(config)
+        case _:
+            raise NotImplementedError(f"{name} is not a valid model name")
+
+
 @define(kw_only=True)
 class WeibullReliabilityModel(BaseConfig):
     r"""Basic reliability model for operating/not operating statuses.

From 58bdff9c3629a12435f636c8ebef12bda99e577f Mon Sep 17 00:00:00 2001
From: "Hammond, Rob" <13874373+RHammond2@users.noreply.github.com>
Date: Mon, 17 Aug 2026 12:50:09 -0700
Subject: [PATCH 04/24] make downtime and reliability separate models

---
 h2integrate/reliability/models.py | 106 +++++++++++++++++++++++++-----
 1 file changed, 88 insertions(+), 18 deletions(-)

diff --git a/h2integrate/reliability/models.py b/h2integrate/reliability/models.py
index cf8ca4dba..6c4bf0050 100644
--- a/h2integrate/reliability/models.py
+++ b/h2integrate/reliability/models.py
@@ -1,9 +1,16 @@
+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
 
 
@@ -11,14 +18,49 @@ def create_reliability(config: dict):
     """Retrieves and initializes a matching reliability model."""
     name = config.pop("reliability")
     match name:
-        case "WeibullReliabilityModel":
-            return WeibullReliabilityModel.from_dict(config)
+        case "WeibullReliability":
+            return WeibullReliability.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 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 WeibullReliabilityModel(BaseConfig):
+class WeibullReliability(BaseConfig):
     r"""Basic reliability model for operating/not operating statuses.
 
     Assumes a full operational shutdown with zero ramping of production for an hourly, 1 year
@@ -45,34 +87,62 @@ class WeibullReliabilityModel(BaseConfig):
 
     scale: float = field(validator=validators.instance_of(float))
     shape: float = field(validator=validators.instance_of(float))
-    downtime: float = field(validator=validators.gt(1))
-    rng: np.random._generator.Generator = field(
-        default=np.random.default_rng(),
-        init=False,
-        validator=validators.instance_of(np.random._generator.Generator),
-    )
+    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)
     )
-    downtime_per_event: np.ndarray = field(
-        default=np.zeros(N_TIMESTEPS), init=False, validator=validators.instance_of(np.ndarray)
-    )
 
     def __attrs_post_init__(self):
         self.create_downtime_events()
         self.calculate_availability()
 
     def create_downtime_events(self):
-        """Creates a ``time_to_failure`` and ``downtime_per_event`` array based on the distributions
-        described in ``WeibullReliabilityConfig``.
-        """
+        """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 = np.ceil(
-            self.config.rng.weibul(self.config.shape, size=30) * self.config.scale * N_TIMESTEPS
+            self.rng.weibul(self.shape, size=100) * self.scale * N_TIMESTEPS
         ).astype(int)
-        downtime_per_event = np.ceil(np.rng.normal(loc=self.config.downtime, size=30)).astype(int)
-        self.downtime_per_event = np.where(downtime_per_event >= 1, downtime_per_event, 1)
+        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 = start
+            if not self.time_to_failures:
+                self.create_downtime_events()
+
+
+@define(kw_only=True)
+class FixedIntervalReliability(BaseConfig):
+    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 create_downtime_events(self):
+        """Creates a ``time_to_failure`` and ``downtime_per_event``."""
+        self.time_to_failures = ...  # TODO
+        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."""

From 4283ab591519753608d3d533b9ef7b45116e775c Mon Sep 17 00:00:00 2001
From: "Hammond, Rob" <13874373+RHammond2@users.noreply.github.com>
Date: Tue, 18 Aug 2026 09:39:18 -0700
Subject: [PATCH 05/24] update sampling and fix interval availability

---
 h2integrate/reliability/models.py | 80 +++++++++++++++++++++++--------
 1 file changed, 61 insertions(+), 19 deletions(-)

diff --git a/h2integrate/reliability/models.py b/h2integrate/reliability/models.py
index 6c4bf0050..850a4dcee 100644
--- a/h2integrate/reliability/models.py
+++ b/h2integrate/reliability/models.py
@@ -14,12 +14,26 @@
 N_TIMESTEPS = 8760
 
 
-def create_reliability(config: dict):
+def create_failure_model(config: dict):
     """Retrieves and initializes a matching reliability model."""
-    name = config.pop("reliability")
+    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")
 
@@ -39,6 +53,12 @@ class BaseDowntime(ABC, BaseConfig):
     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.
@@ -60,7 +80,7 @@ def sample_downtime(self) -> np.ndarray:
 
 
 @define(kw_only=True)
-class WeibullReliability(BaseConfig):
+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
@@ -97,13 +117,14 @@ 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 = np.ceil(
-            self.rng.weibul(self.shape, size=100) * self.scale * N_TIMESTEPS
-        ).astype(int)
+        self.time_to_failures = self.sample_events()
         self.downtime_per_event = self.downtime.sample_downtime()
 
     def calculate_availability(self):
@@ -121,13 +142,25 @@ def calculate_availability(self):
             start = accumulated + event
             end = start + duration
             self.availability[start:end] = 0
-            accumulated = start
+            accumulated = end
             if not self.time_to_failures:
                 self.create_downtime_events()
 
 
 @define(kw_only=True)
-class FixedIntervalReliability(BaseConfig):
+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))
@@ -139,26 +172,35 @@ 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.time_to_failures = ...  # TODO
+        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."""
-        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:],
+        )
+        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:],
             )
-            if event + accumulated > N_TIMESTEPS:
-                break
-
-            start = accumulated + event
+            # TODO: handle start > 8760
             end = start + duration
+            end = np.where(end > 8760, 8760, end)
             self.availability[start:end] = 0
-            accumulated = start
-            if not self.time_to_failures:
-                self.create_downtime_events()
+            accumulated = end
+            if not self.downtime_per_event:
+                self.downtime.sample_downtime()

From f5e61da7dfa150c5b616906f688bb9cd5d71ad80 Mon Sep 17 00:00:00 2001
From: "Hammond, Rob" <13874373+RHammond2@users.noreply.github.com>
Date: Wed, 19 Aug 2026 08:37:58 -0700
Subject: [PATCH 06/24] update downtime for matrix and generation helpers

---
 h2integrate/reliability/models.py | 50 +++++++++++++++++++++++++------
 1 file changed, 41 insertions(+), 9 deletions(-)

diff --git a/h2integrate/reliability/models.py b/h2integrate/reliability/models.py
index 850a4dcee..a60484edf 100644
--- a/h2integrate/reliability/models.py
+++ b/h2integrate/reliability/models.py
@@ -3,6 +3,7 @@
 
 import numpy as np
 from attrs import field, define, validators
+from numpy.typing import ArrayLike
 
 from h2integrate.core.utilities import BaseConfig
 
@@ -16,12 +17,13 @@
 
 def create_failure_model(config: dict):
     """Retrieves and initializes a matching reliability model."""
-    name = config.pop("failure_model")
+    name = config["failure_model"]
+    fail_config = config["failure_parameters"]
     match name:
         case "WeibullReliability":
-            return WeibullReliability.from_dict(config)
+            return WeibullReliability.from_dict()
         case "FixedIntervalReliability":
-            return FixedIntervalReliability.from_dict(config)
+            return FixedIntervalReliability.from_dict(fail_config)
         case _:
             raise NotImplementedError(f"{name} is not a valid model name")
 
@@ -29,11 +31,12 @@ def create_failure_model(config: dict):
 def create_maintenance_model(config: dict):
     """Retrieves and initializes a matching reliability model."""
     name = config.pop("maintenance_model")
+    maintenance_config = config["maintenance_parameters"]
     match name:
         case "WeibullReliability":
-            return WeibullReliability.from_dict(config)
+            return WeibullReliability.from_dict(maintenance_config)
         case "FixedIntervalReliability":
-            return FixedIntervalReliability.from_dict(config)
+            return FixedIntervalReliability.from_dict(maintenance_config)
         case _:
             raise NotImplementedError(f"{name} is not a valid model name")
 
@@ -59,6 +62,14 @@ class BaseReliability(ABC, BaseConfig):
     def sample_events(self) -> np.ndarray: ...
 
 
+def float_array_converter(val: int | float | ArrayLike):
+    return np.ndarray(val).astype(float).reshape(-1, 1)
+
+
+def int_array_converter(val: int | float | ArrayLike):
+    return np.ndarray(val).astype(float).reshape(-1, 1)
+
+
 @define(kw_only=True)
 class LogNormalDowntime(BaseDowntime):
     """Basic log-normal downtime model for generating the length of downtime for a given event.
@@ -70,13 +81,34 @@ class LogNormalDowntime(BaseDowntime):
             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)))
+    mean: float = field(
+        converter=float_array_converter,
+        validator=(validators.instance_of(np.ndarray), validators.ge(0)),
+    )
+    sigma: float = field(
+        converter=float_array_converter,
+        validator=(validators.instance_of(np.ndarray), validators.ge(0)),
+    )
     n_components: int = field(default=1, validator=(validators.instance_of(int), validators.ge(1)))
 
+    @sigma.validator
+    def validate_shape(self, attribute, value):
+        """Validates that :py:attr:`mean` and :py:attr:`sigma` are the same size."""
+        if self.mean.shape != value.shape:
+            msg = (
+                "Inputs to 'mean' and 'sigma' must be the same length. Received"
+                f" 'mean': {self.mean.size}, 'sigma': {value.size}"
+            )
+            raise ValueError(msg)
+
+    def __attrs_post_init__(self):
+        if self.mean.size == 1 and self.n_components > 1:
+            broadcaster = np.ones((self.n_components, 1))
+            self.mean *= broadcaster
+            self.sigma *= broadcaster
+
     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)
+        return rng.lognormal(self.mean, self.sigma, size=(self.mean.shape[0], 100))
 
 
 @define(kw_only=True)

From dce308e585e874f0f64e9e66897c7b7dd7af2244 Mon Sep 17 00:00:00 2001
From: "Hammond, Rob" <13874373+RHammond2@users.noreply.github.com>
Date: Wed, 19 Aug 2026 10:21:13 -0700
Subject: [PATCH 07/24] make weibull model definiton multi-component-compatible

---
 h2integrate/reliability/models.py | 30 ++++++++++++++++++++----------
 1 file changed, 20 insertions(+), 10 deletions(-)

diff --git a/h2integrate/reliability/models.py b/h2integrate/reliability/models.py
index a60484edf..155f7de01 100644
--- a/h2integrate/reliability/models.py
+++ b/h2integrate/reliability/models.py
@@ -75,10 +75,10 @@ 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 | array-like): Average length of downtime per event, in hours.
+        sigma (float | array-like): Standard deviation of the distribution(s), in hours.
+        n_components (int): Number of identical components to sample to avoid defining an array of
+            mean and sigma values when they are the same. Defaults to 1.
     """
 
     mean: float = field(
@@ -137,20 +137,29 @@ class WeibullReliability(BaseReliability):
         - 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))
+    scale: float = field(
+        converter=float_array_converter,
+        validator=validators.instance_of(np.ndarray),
+    )
+    shape: float = field(
+        converter=float_array_converter,
+        validator=validators.instance_of(np.ndarray),
+    )
+    n_components: int = field(default=1, validator=(validators.instance_of(int), validators.ge(1)))
     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)
-    )
+    availability: np.ndarray = field(init=False, validator=validators.instance_of(np.ndarray))
 
     def __attrs_post_init__(self):
+        self.availability = np.ones(self.scale.size, N_TIMESTEPS, dtype=float)
         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)
+        """Samples 100 events for each simulated component, rounding up to the nearest timestep."""
+        return np.ceil(
+            self.rng.weibul(self.shape, size=(self.shape, 100)) * self.scale * N_TIMESTEPS
+        ).astype(int)
 
     def create_downtime_events(self):
         """Creates a ``time_to_failure`` and ``downtime_per_event``."""
@@ -161,6 +170,7 @@ def create_downtime_events(self):
 
     def calculate_availability(self):
         """Determine the timing and duration of outages for a single year of simulation time."""
+        # TODO: convert to matrix compatible variation, not just single component version
         accumulated = 0
         while accumulated < N_TIMESTEPS:
             event, self.time_to_failures = self.time_to_failures[0], self.time_to_failures[1:]

From 2463e78f4032545883da71c0e3900ce74f6372fc Mon Sep 17 00:00:00 2001
From: "Hammond, Rob" <13874373+RHammond2@users.noreply.github.com>
Date: Wed, 19 Aug 2026 12:48:18 -0700
Subject: [PATCH 08/24] getting setup partially working

---
 h2integrate/reliability/models.py | 23 ++++++++++++-----------
 1 file changed, 12 insertions(+), 11 deletions(-)

diff --git a/h2integrate/reliability/models.py b/h2integrate/reliability/models.py
index 155f7de01..ebf79986d 100644
--- a/h2integrate/reliability/models.py
+++ b/h2integrate/reliability/models.py
@@ -21,7 +21,7 @@ def create_failure_model(config: dict):
     fail_config = config["failure_parameters"]
     match name:
         case "WeibullReliability":
-            return WeibullReliability.from_dict()
+            return WeibullReliability.from_dict(fail_config)
         case "FixedIntervalReliability":
             return FixedIntervalReliability.from_dict(fail_config)
         case _:
@@ -41,7 +41,9 @@ def create_maintenance_model(config: dict):
             raise NotImplementedError(f"{name} is not a valid model name")
 
 
-def generate_downtime_model(config: dict):
+def generate_downtime_model(config: dict | int):
+    if not isinstance(config, dict):
+        return config
     name = config.pop("model")
     match name:
         case "LogNormalDowntime":
@@ -63,11 +65,11 @@ def sample_events(self) -> np.ndarray: ...
 
 
 def float_array_converter(val: int | float | ArrayLike):
-    return np.ndarray(val).astype(float).reshape(-1, 1)
+    return np.array(val).astype(float).reshape(-1, 1)
 
 
 def int_array_converter(val: int | float | ArrayLike):
-    return np.ndarray(val).astype(float).reshape(-1, 1)
+    return np.array(val).astype(float).reshape(-1, 1)
 
 
 @define(kw_only=True)
@@ -134,7 +136,6 @@ class WeibullReliability(BaseReliability):
 
     TODO:
         - how to pass n_timesteps through from plant?
-        - stabilize random generator/determine how to manage random seeding across library
     """
 
     scale: float = field(
@@ -146,31 +147,31 @@ class WeibullReliability(BaseReliability):
         validator=validators.instance_of(np.ndarray),
     )
     n_components: int = field(default=1, validator=(validators.instance_of(int), validators.ge(1)))
-    downtime: Any = field(converter=generate_downtime_model)
+    downtime: int | BaseDowntime = field(converter=generate_downtime_model)
+    time_to_failures: np.ndarray = field(init=False, validator=validators.instance_of(np.ndarray))
     downtime_per_event: np.ndarray = field(init=False, validator=validators.instance_of(np.ndarray))
     availability: np.ndarray = field(init=False, validator=validators.instance_of(np.ndarray))
 
     def __attrs_post_init__(self):
-        self.availability = np.ones(self.scale.size, N_TIMESTEPS, dtype=float)
+        self.availability = np.ones((self.scale.size, N_TIMESTEPS), dtype=float)
         self.create_downtime_events()
-        self.calculate_availability()
+        # self.calculate_availability()
 
     def sample_events(self):
         """Samples 100 events for each simulated component, rounding up to the nearest timestep."""
         return np.ceil(
-            self.rng.weibul(self.shape, size=(self.shape, 100)) * self.scale * N_TIMESTEPS
+            rng.weibull(self.shape, size=(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."""
         # TODO: convert to matrix compatible variation, not just single component version
+        # only modify rows that haven't reached 8760 yet, then
         accumulated = 0
         while accumulated < N_TIMESTEPS:
             event, self.time_to_failures = self.time_to_failures[0], self.time_to_failures[1:]

From efaf66541d44dc03512de25b8508a35e6335d6cc Mon Sep 17 00:00:00 2001
From: "Hammond, Rob" <13874373+RHammond2@users.noreply.github.com>
Date: Wed, 19 Aug 2026 14:36:27 -0700
Subject: [PATCH 09/24] update docstrings and add fixed downtime

---
 h2integrate/reliability/models.py | 22 ++++++++++++++++++++++
 1 file changed, 22 insertions(+)

diff --git a/h2integrate/reliability/models.py b/h2integrate/reliability/models.py
index ebf79986d..34b4988c5 100644
--- a/h2integrate/reliability/models.py
+++ b/h2integrate/reliability/models.py
@@ -72,6 +72,27 @@ def int_array_converter(val: int | float | ArrayLike):
     return np.array(val).astype(float).reshape(-1, 1)
 
 
+@define(kw_only=True)
+class FixedDowntime(BaseDowntime):
+    """Basic log-normal downtime model for generating the length of downtime for a given event.
+
+    Args:
+        hours (int | array-like): Length of downtime per event, in hours. Must be at least 1 hour.
+        n_components (int): Number of identical components to sample to avoid defining an array of
+            :py:attr:`hours` values when they are the same. Defaults to 1.
+    """
+
+    hours: int | ArrayLike = field(
+        converter=int_array_converter,
+        validator=(validators.instance_of(np.ndarray), validators.ge(1)),
+    )
+    n_components: int = field(default=1, validator=(validators.instance_of(int), validators.ge(1)))
+
+    def sample_downtime(self):
+        """Return an array of 100 :py:attr:`hours`."""
+        return np.ones((1, 100)) * self.hours
+
+
 @define(kw_only=True)
 class LogNormalDowntime(BaseDowntime):
     """Basic log-normal downtime model for generating the length of downtime for a given event.
@@ -110,6 +131,7 @@ def __attrs_post_init__(self):
             self.sigma *= broadcaster
 
     def sample_downtime(self) -> np.ndarray:
+        """Return an array of 100 samples of each lognormal distribution."""
         return rng.lognormal(self.mean, self.sigma, size=(self.mean.shape[0], 100))
 
 

From e4bfcb8a5f25f8e9930dd8bdac0516ac1935bb58 Mon Sep 17 00:00:00 2001
From: "Hammond, Rob" <13874373+RHammond2@users.noreply.github.com>
Date: Wed, 19 Aug 2026 16:20:51 -0700
Subject: [PATCH 10/24] fix small bugs and ensure matrix works for weibull

---
 h2integrate/reliability/models.py | 56 ++++++++++++++++++++-----------
 1 file changed, 36 insertions(+), 20 deletions(-)

diff --git a/h2integrate/reliability/models.py b/h2integrate/reliability/models.py
index 34b4988c5..6051b1a62 100644
--- a/h2integrate/reliability/models.py
+++ b/h2integrate/reliability/models.py
@@ -30,7 +30,7 @@ def create_failure_model(config: dict):
 
 def create_maintenance_model(config: dict):
     """Retrieves and initializes a matching reliability model."""
-    name = config.pop("maintenance_model")
+    name = config["maintenance_model"]
     maintenance_config = config["maintenance_parameters"]
     match name:
         case "WeibullReliability":
@@ -43,11 +43,14 @@ def create_maintenance_model(config: dict):
 
 def generate_downtime_model(config: dict | int):
     if not isinstance(config, dict):
-        return config
-    name = config.pop("model")
+        return FixedDowntime(hours=config)
+    name = config["model"]
+    parameters = {k: v for k, v in config.items() if k != "model"}
     match name:
         case "LogNormalDowntime":
-            return LogNormalDowntime.from_dict(config)
+            return LogNormalDowntime.from_dict(parameters)
+        case "FixedDowntime":
+            return FixedDowntime.from_dict(parameters)
         case _:
             raise NotImplementedError(f"{name} is not a valid model name")
 
@@ -69,7 +72,7 @@ def float_array_converter(val: int | float | ArrayLike):
 
 
 def int_array_converter(val: int | float | ArrayLike):
-    return np.array(val).astype(float).reshape(-1, 1)
+    return np.array(val).astype(int).reshape(-1, 1)
 
 
 @define(kw_only=True)
@@ -84,13 +87,24 @@ class FixedDowntime(BaseDowntime):
 
     hours: int | ArrayLike = field(
         converter=int_array_converter,
-        validator=(validators.instance_of(np.ndarray), validators.ge(1)),
+        validator=validators.instance_of(np.ndarray),
     )
     n_components: int = field(default=1, validator=(validators.instance_of(int), validators.ge(1)))
 
+    @hours.validator
+    def hours_validator(self, attribute, value):
+        """Validates that all values of :py:attr:`hours` are greater than or equal to 1."""
+        if not np.all(value > 1):
+            raise ValueError("All values passed to 'hours' must be greater than or equal to 1.")
+
+    def __attrs_post_init__(self):
+        if self.hours.size == 1 and self.n_components > 1:
+            self.hours = np.broadcast_to(self.hours, (self.n_components, 1))
+        # ...
+
     def sample_downtime(self):
         """Return an array of 100 :py:attr:`hours`."""
-        return np.ones((1, 100)) * self.hours
+        return np.ones((1, 100), dtype=int) * self.hours
 
 
 @define(kw_only=True)
@@ -158,6 +172,7 @@ class WeibullReliability(BaseReliability):
 
     TODO:
         - how to pass n_timesteps through from plant?
+        - burn-in
     """
 
     scale: float = field(
@@ -192,24 +207,25 @@ def create_downtime_events(self):
 
     def calculate_availability(self):
         """Determine the timing and duration of outages for a single year of simulation time."""
-        # TODO: convert to matrix compatible variation, not just single component version
-        # only modify rows that haven't reached 8760 yet, then
-        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:
+        accumulated = np.zeros_like(self.shape, dtype=int)
+        while any(accumulated < N_TIMESTEPS):
+            if not self.time_to_failures.size == 0:
+                self.create_downtime_events()
+            event = self.time_to_failures[:, 0].reshape(-1, 1)
+            self.time_to_failures = self.time_to_failures[:, 1:]
+            duration = self.downtime_per_event[:, 0].reshape(-1, 1)
+            self.downtime_per_event = self.downtime_per_event[:, 1:]
+
+            if all(event + accumulated) > N_TIMESTEPS:
                 break
 
             start = accumulated + event
             end = start + duration
-            self.availability[start:end] = 0
+            for i, (s, e) in enumerate(zip(start.flatten(), end.flatten())):
+                if s < N_TIMESTEPS:
+                    e = min(N_TIMESTEPS, e)
+                    self.availability[i, s:e] = 0
             accumulated = end
-            if not self.time_to_failures:
-                self.create_downtime_events()
 
 
 @define(kw_only=True)

From 6388656f39538515e00c9384ddac2b0a69ade4ce Mon Sep 17 00:00:00 2001
From: "Hammond, Rob" <13874373+RHammond2@users.noreply.github.com>
Date: Wed, 19 Aug 2026 16:30:20 -0700
Subject: [PATCH 11/24] system availability

---
 h2integrate/reliability/models.py | 5 +++++
 1 file changed, 5 insertions(+)

diff --git a/h2integrate/reliability/models.py b/h2integrate/reliability/models.py
index 6051b1a62..02f7edfd4 100644
--- a/h2integrate/reliability/models.py
+++ b/h2integrate/reliability/models.py
@@ -188,6 +188,9 @@ class WeibullReliability(BaseReliability):
     time_to_failures: np.ndarray = field(init=False, validator=validators.instance_of(np.ndarray))
     downtime_per_event: np.ndarray = field(init=False, validator=validators.instance_of(np.ndarray))
     availability: np.ndarray = field(init=False, validator=validators.instance_of(np.ndarray))
+    system_availability: np.ndarray = field(
+        init=False, validator=validators.instance_of(np.ndarray)
+    )
 
     def __attrs_post_init__(self):
         self.availability = np.ones((self.scale.size, N_TIMESTEPS), dtype=float)
@@ -227,6 +230,8 @@ def calculate_availability(self):
                     self.availability[i, s:e] = 0
             accumulated = end
 
+        self.system_availability = np.min(self.availability, axis=0)
+
 
 @define(kw_only=True)
 class FixedIntervalReliability(BaseReliability):

From 081386b83e5b7aaa53d7bdf2bd746ad8396adb9f Mon Sep 17 00:00:00 2001
From: "Hammond, Rob" <13874373+RHammond2@users.noreply.github.com>
Date: Wed, 19 Aug 2026 16:42:20 -0700
Subject: [PATCH 12/24] move calculate_availability to base and simplify
 broadcasting

---
 h2integrate/reliability/models.py | 89 +++++++++++++------------------
 1 file changed, 36 insertions(+), 53 deletions(-)

diff --git a/h2integrate/reliability/models.py b/h2integrate/reliability/models.py
index 02f7edfd4..84a5ac496 100644
--- a/h2integrate/reliability/models.py
+++ b/h2integrate/reliability/models.py
@@ -62,9 +62,33 @@ def sample_downtime(self) -> np.ndarray: ...
 
 
 @define(kw_only=True)
-class BaseReliability(ABC, BaseConfig):
-    @abstractmethod
-    def sample_events(self) -> np.ndarray: ...
+class BaseReliability(BaseConfig):
+    def sample_events(self) -> np.ndarray:
+        raise NotImplementedError("Failed to successfully subclass, please implement me.")
+
+    def calculate_availability(self):
+        """Determine the timing and duration of outages for a single year of simulation time."""
+        accumulated = np.zeros_like(self.shape, dtype=int)
+        while any(accumulated < N_TIMESTEPS):
+            if not self.time_to_failures.size == 0:
+                self.create_downtime_events()
+            event = self.time_to_failures[:, 0].reshape(-1, 1)
+            self.time_to_failures = self.time_to_failures[:, 1:]
+            duration = self.downtime_per_event[:, 0].reshape(-1, 1)
+            self.downtime_per_event = self.downtime_per_event[:, 1:]
+
+            if all(event + accumulated) > N_TIMESTEPS:
+                break
+
+            start = accumulated + event
+            end = start + duration
+            for i, (s, e) in enumerate(zip(start.flatten(), end.flatten())):
+                if s < N_TIMESTEPS:
+                    e = min(N_TIMESTEPS, e)
+                    self.availability[i, s:e] = 0
+            accumulated = end
+
+        self.system_availability = np.min(self.availability, axis=0)
 
 
 def float_array_converter(val: int | float | ArrayLike):
@@ -140,9 +164,9 @@ def validate_shape(self, attribute, value):
 
     def __attrs_post_init__(self):
         if self.mean.size == 1 and self.n_components > 1:
-            broadcaster = np.ones((self.n_components, 1))
-            self.mean *= broadcaster
-            self.sigma *= broadcaster
+            shape = (self.n_components, 1)
+            self.mean = np.broadcast_to(self.mean, shape)
+            self.sigma = np.broadcast_to(self.sigma, shape)
 
     def sample_downtime(self) -> np.ndarray:
         """Return an array of 100 samples of each lognormal distribution."""
@@ -193,9 +217,14 @@ class WeibullReliability(BaseReliability):
     )
 
     def __attrs_post_init__(self):
+        if self.scale.size == 1 and self.n_components > 1:
+            shape = (self.n_components, 1)
+            self.scale = np.broadcast_to(self.scale, shape)
+            self.shape = np.broadcast_to(self.shape, shape)
+
         self.availability = np.ones((self.scale.size, N_TIMESTEPS), dtype=float)
         self.create_downtime_events()
-        # self.calculate_availability()
+        self.calculate_availability()
 
     def sample_events(self):
         """Samples 100 events for each simulated component, rounding up to the nearest timestep."""
@@ -208,30 +237,6 @@ def create_downtime_events(self):
         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 = np.zeros_like(self.shape, dtype=int)
-        while any(accumulated < N_TIMESTEPS):
-            if not self.time_to_failures.size == 0:
-                self.create_downtime_events()
-            event = self.time_to_failures[:, 0].reshape(-1, 1)
-            self.time_to_failures = self.time_to_failures[:, 1:]
-            duration = self.downtime_per_event[:, 0].reshape(-1, 1)
-            self.downtime_per_event = self.downtime_per_event[:, 1:]
-
-            if all(event + accumulated) > N_TIMESTEPS:
-                break
-
-            start = accumulated + event
-            end = start + duration
-            for i, (s, e) in enumerate(zip(start.flatten(), end.flatten())):
-                if s < N_TIMESTEPS:
-                    e = min(N_TIMESTEPS, e)
-                    self.availability[i, s:e] = 0
-            accumulated = end
-
-        self.system_availability = np.min(self.availability, axis=0)
-
 
 @define(kw_only=True)
 class FixedIntervalReliability(BaseReliability):
@@ -268,25 +273,3 @@ 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()

From e5cbdd9fe26aae00b07f145838d8dac0f45137ac Mon Sep 17 00:00:00 2001
From: "Hammond, Rob" <13874373+RHammond2@users.noreply.github.com>
Date: Thu, 20 Aug 2026 09:54:39 -0700
Subject: [PATCH 13/24] convert fixexd interval to componentizable model and
 update n_components in every model

---
 h2integrate/reliability/models.py | 73 +++++++++++++++++++++++--------
 1 file changed, 54 insertions(+), 19 deletions(-)

diff --git a/h2integrate/reliability/models.py b/h2integrate/reliability/models.py
index 84a5ac496..1a28a6476 100644
--- a/h2integrate/reliability/models.py
+++ b/h2integrate/reliability/models.py
@@ -1,5 +1,4 @@
 from abc import ABC, abstractmethod
-from typing import Any
 
 import numpy as np
 from attrs import field, define, validators
@@ -68,7 +67,7 @@ def sample_events(self) -> np.ndarray:
 
     def calculate_availability(self):
         """Determine the timing and duration of outages for a single year of simulation time."""
-        accumulated = np.zeros_like(self.shape, dtype=int)
+        accumulated = np.zeros((self.n_components, 1), dtype=int)
         while any(accumulated < N_TIMESTEPS):
             if not self.time_to_failures.size == 0:
                 self.create_downtime_events()
@@ -163,10 +162,13 @@ def validate_shape(self, attribute, value):
             raise ValueError(msg)
 
     def __attrs_post_init__(self):
-        if self.mean.size == 1 and self.n_components > 1:
-            shape = (self.n_components, 1)
-            self.mean = np.broadcast_to(self.mean, shape)
-            self.sigma = np.broadcast_to(self.sigma, shape)
+        if self.mean.size == 1:
+            if self.n_components > 1:
+                shape = (self.n_components, 1)
+                self.mean = np.broadcast_to(self.mean, shape)
+                self.sigma = np.broadcast_to(self.sigma, shape)
+            else:
+                self.n_components = 1
 
     def sample_downtime(self) -> np.ndarray:
         """Return an array of 100 samples of each lognormal distribution."""
@@ -209,6 +211,7 @@ class WeibullReliability(BaseReliability):
     )
     n_components: int = field(default=1, validator=(validators.instance_of(int), validators.ge(1)))
     downtime: int | BaseDowntime = field(converter=generate_downtime_model)
+
     time_to_failures: np.ndarray = field(init=False, validator=validators.instance_of(np.ndarray))
     downtime_per_event: np.ndarray = field(init=False, validator=validators.instance_of(np.ndarray))
     availability: np.ndarray = field(init=False, validator=validators.instance_of(np.ndarray))
@@ -217,10 +220,13 @@ class WeibullReliability(BaseReliability):
     )
 
     def __attrs_post_init__(self):
-        if self.scale.size == 1 and self.n_components > 1:
-            shape = (self.n_components, 1)
-            self.scale = np.broadcast_to(self.scale, shape)
-            self.shape = np.broadcast_to(self.shape, shape)
+        if self.scale.size == 1:
+            if self.n_components > 1:
+                shape = (self.n_components, 1)
+                self.scale = np.broadcast_to(self.scale, shape)
+                self.shape = np.broadcast_to(self.shape, shape)
+            else:
+                self.n_components = 1
 
         self.availability = np.ones((self.scale.size, N_TIMESTEPS), dtype=float)
         self.create_downtime_events()
@@ -249,27 +255,56 @@ class FixedIntervalReliability(BaseReliability):
             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)
+    frequency: float = field(
+        converter=float_array_converter,
+        validator=validators.instance_of(np.ndarray),
+    )
+    downtime: int | BaseDowntime = 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)
+    n_components: int = field(default=1, validator=(validators.instance_of(int), validators.ge(1)))
+
+    time_to_failures: np.ndarray = field(init=False, validator=validators.instance_of(np.ndarray))
+    availability: np.ndarray = field(init=False, validator=validators.instance_of(np.ndarray))
+    system_availability: np.ndarray = field(
+        init=False, validator=validators.instance_of(np.ndarray)
     )
 
     def __attrs_post_init__(self):
+        if self.frequency.size == 1:
+            if self.n_components > 1:
+                shape = (self.n_components, 1)
+                self.frequency = np.broadcast_to(self.frequency, shape)
+            else:
+                self.n_components = 1
+
+        self.availability = np.ones((self.frequency.size, N_TIMESTEPS), dtype=float)
         self.create_downtime_events()
         self.calculate_availability()
 
+    @frequency.validator
+    def validate_frequency(self, attribute, value):
+        """Validates all passed values to :py:attr:`frequency` are greater than 0."""
+        if any(value <= 0):
+            raise ValueError("All values of 'frequency' must be greater than 0.")
+
     def sample_events(self):
-        frequency = np.array([0.1, 20, 30]).reshape(-1, 1)
-        interval = np.ceil(8760 / frequency).astype(int)
+        """Creates the time to next failure array for each event's modality with the first event
+        occurring randomly either in the first year or first failure interval, whichever
+        comes first.
+
+        Returns:
+            time_to_failures (np.ndarray): An array of the next 100 events' time to next failure.
+        """
+        interval = np.ceil(8760 / (1 / self.frequency)).astype(int)
         first_occurrence = rng.integers(0, np.where(interval > 8760, 8760, interval))
-        return first_occurrence, interval
+        time_to_failures = np.hstack(
+            (first_occurrence, np.broadcast_to(interval, (interval.size, 99)))
+        )
+        return time_to_failures
 
     def create_downtime_events(self):
         """Creates a ``time_to_failure`` and ``downtime_per_event``."""
-        self.first_occurrence, self.interval = self.sample_events()
+        self.time_to_failures = self.sample_events()
         self.downtime_per_event = self.downtime.sample_downtime()

From 4e326bbab1ee8ac9d7245276a60bcd595f0632db Mon Sep 17 00:00:00 2001
From: "Hammond, Rob" <13874373+RHammond2@users.noreply.github.com>
Date: Thu, 20 Aug 2026 09:58:08 -0700
Subject: [PATCH 14/24] update fixed downtime

---
 h2integrate/reliability/models.py | 9 ++++++---
 1 file changed, 6 insertions(+), 3 deletions(-)

diff --git a/h2integrate/reliability/models.py b/h2integrate/reliability/models.py
index 1a28a6476..1f8cceeb8 100644
--- a/h2integrate/reliability/models.py
+++ b/h2integrate/reliability/models.py
@@ -121,9 +121,12 @@ def hours_validator(self, attribute, value):
             raise ValueError("All values passed to 'hours' must be greater than or equal to 1.")
 
     def __attrs_post_init__(self):
-        if self.hours.size == 1 and self.n_components > 1:
-            self.hours = np.broadcast_to(self.hours, (self.n_components, 1))
-        # ...
+        if self.hours.size == 1:
+            if self.n_components > 1:
+                shape = (self.n_components, 1)
+                self.hours = np.broadcast_to(self.hours, shape)
+            else:
+                self.n_components = 1
 
     def sample_downtime(self):
         """Return an array of 100 :py:attr:`hours`."""

From b3544be25e5540096ce61095fd4a40dc915a6bb8 Mon Sep 17 00:00:00 2001
From: "Hammond, Rob" <13874373+RHammond2@users.noreply.github.com>
Date: Thu, 20 Aug 2026 10:11:22 -0700
Subject: [PATCH 15/24] update lognormal downtime model

---
 h2integrate/reliability/models.py | 15 ++++++++++++---
 1 file changed, 12 insertions(+), 3 deletions(-)

diff --git a/h2integrate/reliability/models.py b/h2integrate/reliability/models.py
index 1f8cceeb8..c25f776ec 100644
--- a/h2integrate/reliability/models.py
+++ b/h2integrate/reliability/models.py
@@ -146,14 +146,21 @@ class LogNormalDowntime(BaseDowntime):
 
     mean: float = field(
         converter=float_array_converter,
-        validator=(validators.instance_of(np.ndarray), validators.ge(0)),
+        validator=(validators.instance_of(np.ndarray)),
     )
     sigma: float = field(
         converter=float_array_converter,
-        validator=(validators.instance_of(np.ndarray), validators.ge(0)),
+        validator=(validators.instance_of(np.ndarray)),
     )
     n_components: int = field(default=1, validator=(validators.instance_of(int), validators.ge(1)))
 
+    @mean.validator
+    @sigma.validator
+    def validate_ge_zero(self, attribute, value):
+        """Validates the :py:attr:`sigma` and :py:attr:`mean` values are at least 0."""
+        if any(value < 0):
+            raise ValueError(f"'{attribute.name}' must have all values of at least 0.")
+
     @sigma.validator
     def validate_shape(self, attribute, value):
         """Validates that :py:attr:`mean` and :py:attr:`sigma` are the same size."""
@@ -175,7 +182,9 @@ def __attrs_post_init__(self):
 
     def sample_downtime(self) -> np.ndarray:
         """Return an array of 100 samples of each lognormal distribution."""
-        return rng.lognormal(self.mean, self.sigma, size=(self.mean.shape[0], 100))
+        return np.ceil(rng.lognormal(self.mean, self.sigma, size=(self.mean.shape[0], 100))).astype(
+            int
+        )
 
 
 @define(kw_only=True)

From c809904d0bf97032ed2caef152ff6ff34ba29f15 Mon Sep 17 00:00:00 2001
From: "Hammond, Rob" <13874373+RHammond2@users.noreply.github.com>
Date: Thu, 20 Aug 2026 11:00:10 -0700
Subject: [PATCH 16/24] add burn-in and fix bugs in n_components setting

---
 h2integrate/reliability/models.py | 119 ++++++++++++++++--------------
 1 file changed, 64 insertions(+), 55 deletions(-)

diff --git a/h2integrate/reliability/models.py b/h2integrate/reliability/models.py
index c25f776ec..9317ed649 100644
--- a/h2integrate/reliability/models.py
+++ b/h2integrate/reliability/models.py
@@ -62,13 +62,51 @@ def sample_downtime(self) -> np.ndarray: ...
 
 @define(kw_only=True)
 class BaseReliability(BaseConfig):
+    """Base reliability class responsible for common definitions and functionality.
+
+    Args:
+        burn_in (float): Number of years into the simulation to use as the starting point of the
+            the simulation's availability record. Defaults to 0.
+        n_components (int): Number of identical components to sample to avoid defining an array of
+            mean and sigma values when they are the same. After initialization this value changes
+            to align with the number of components being simulated, regardless of identicality.
+            Defaults to 1.
+        downtime (dict): Configuration for a downtime model.
+
+    Attributes:
+        time_to_failures (np.ndarray): Array of the hours to the next failure for each modeled
+            component. Generated by each subclass' :py:meth:`sample_events`.
+        downtime_per_event (np.ndarray): Number of hours of downtime corresponding to each downtime
+            event in :py:attr:`time_to_failures`.
+        availability (np.ndarray): Operational ratio of each modeled component for every time step
+            of the simulation.
+        system_availability (np.ndarray): Minimum operational ratio of all components for every time
+            time step of the simulation.
+    """
+
+    burn_in: float = field(default=0, converter=float, validator=validators.ge(0))
+    n_components: int = field(default=1, validator=(validators.instance_of(int), validators.ge(1)))
+    downtime: dict | BaseDowntime = field(converter=generate_downtime_model)
+
+    time_to_failures: np.ndarray = field(init=False, validator=validators.instance_of(np.ndarray))
+    downtime_per_event: np.ndarray = field(init=False, validator=validators.instance_of(np.ndarray))
+    availability: np.ndarray = field(init=False, validator=validators.instance_of(np.ndarray))
+    system_availability: np.ndarray = field(
+        init=False, validator=validators.instance_of(np.ndarray)
+    )
+
     def sample_events(self) -> np.ndarray:
         raise NotImplementedError("Failed to successfully subclass, please implement me.")
 
     def calculate_availability(self):
         """Determine the timing and duration of outages for a single year of simulation time."""
+        burn_in_time = np.ceil(self.burn_in * N_TIMESTEPS).astype(int)
+        simulation_end = burn_in_time + N_TIMESTEPS
+
         accumulated = np.zeros((self.n_components, 1), dtype=int)
-        while any(accumulated < N_TIMESTEPS):
+        availability = np.ones((self.n_components, simulation_end), dtype=float)
+
+        while any(accumulated < simulation_end):
             if not self.time_to_failures.size == 0:
                 self.create_downtime_events()
             event = self.time_to_failures[:, 0].reshape(-1, 1)
@@ -76,17 +114,18 @@ def calculate_availability(self):
             duration = self.downtime_per_event[:, 0].reshape(-1, 1)
             self.downtime_per_event = self.downtime_per_event[:, 1:]
 
-            if all(event + accumulated) > N_TIMESTEPS:
+            if all(event + accumulated) > simulation_end:
                 break
 
             start = accumulated + event
             end = start + duration
             for i, (s, e) in enumerate(zip(start.flatten(), end.flatten())):
-                if s < N_TIMESTEPS:
-                    e = min(N_TIMESTEPS, e)
-                    self.availability[i, s:e] = 0
+                if s < simulation_end:
+                    e = min(simulation_end, e)
+                    availability[i, s:e] = 0
             accumulated = end
 
+        self.availability = availability[:, burn_in_time:simulation_end]
         self.system_availability = np.min(self.availability, axis=0)
 
 
@@ -121,12 +160,10 @@ def hours_validator(self, attribute, value):
             raise ValueError("All values passed to 'hours' must be greater than or equal to 1.")
 
     def __attrs_post_init__(self):
-        if self.hours.size == 1:
-            if self.n_components > 1:
-                shape = (self.n_components, 1)
-                self.hours = np.broadcast_to(self.hours, shape)
-            else:
-                self.n_components = 1
+        if self.hours.size == 1 and self.n_components > 1:
+            shape = (self.n_components, 1)
+            self.hours = np.broadcast_to(self.hours, shape)
+        self.n_components = self.hours.size
 
     def sample_downtime(self):
         """Return an array of 100 :py:attr:`hours`."""
@@ -140,8 +177,6 @@ class LogNormalDowntime(BaseDowntime):
     Args:
         mean (float | array-like): Average length of downtime per event, in hours.
         sigma (float | array-like): Standard deviation of the distribution(s), in hours.
-        n_components (int): Number of identical components to sample to avoid defining an array of
-            mean and sigma values when they are the same. Defaults to 1.
     """
 
     mean: float = field(
@@ -172,13 +207,11 @@ def validate_shape(self, attribute, value):
             raise ValueError(msg)
 
     def __attrs_post_init__(self):
-        if self.mean.size == 1:
-            if self.n_components > 1:
-                shape = (self.n_components, 1)
-                self.mean = np.broadcast_to(self.mean, shape)
-                self.sigma = np.broadcast_to(self.sigma, shape)
-            else:
-                self.n_components = 1
+        if self.mean.size == 1 and self.n_components > 1:
+            shape = (self.n_components, 1)
+            self.mean = np.broadcast_to(self.mean, shape)
+            self.sigma = np.broadcast_to(self.sigma, shape)
+        self.n_components = self.mean.size
 
     def sample_downtime(self) -> np.ndarray:
         """Return an array of 100 samples of each lognormal distribution."""
@@ -221,26 +254,14 @@ class WeibullReliability(BaseReliability):
         converter=float_array_converter,
         validator=validators.instance_of(np.ndarray),
     )
-    n_components: int = field(default=1, validator=(validators.instance_of(int), validators.ge(1)))
-    downtime: int | BaseDowntime = field(converter=generate_downtime_model)
-
-    time_to_failures: np.ndarray = field(init=False, validator=validators.instance_of(np.ndarray))
-    downtime_per_event: np.ndarray = field(init=False, validator=validators.instance_of(np.ndarray))
-    availability: np.ndarray = field(init=False, validator=validators.instance_of(np.ndarray))
-    system_availability: np.ndarray = field(
-        init=False, validator=validators.instance_of(np.ndarray)
-    )
 
     def __attrs_post_init__(self):
-        if self.scale.size == 1:
-            if self.n_components > 1:
-                shape = (self.n_components, 1)
-                self.scale = np.broadcast_to(self.scale, shape)
-                self.shape = np.broadcast_to(self.shape, shape)
-            else:
-                self.n_components = 1
-
-        self.availability = np.ones((self.scale.size, N_TIMESTEPS), dtype=float)
+        if self.scale.size == 1 and self.n_components > 1:
+            shape = (self.n_components, 1)
+            self.scale = np.broadcast_to(self.scale, shape)
+            self.shape = np.broadcast_to(self.shape, shape)
+        self.n_components = self.scale.size
+
         self.create_downtime_events()
         self.calculate_availability()
 
@@ -273,25 +294,13 @@ class FixedIntervalReliability(BaseReliability):
         converter=float_array_converter,
         validator=validators.instance_of(np.ndarray),
     )
-    downtime: int | BaseDowntime = field(converter=generate_downtime_model)
-    downtime_per_event: np.ndarray = field(init=False, validator=validators.instance_of(np.ndarray))
-    n_components: int = field(default=1, validator=(validators.instance_of(int), validators.ge(1)))
-
-    time_to_failures: np.ndarray = field(init=False, validator=validators.instance_of(np.ndarray))
-    availability: np.ndarray = field(init=False, validator=validators.instance_of(np.ndarray))
-    system_availability: np.ndarray = field(
-        init=False, validator=validators.instance_of(np.ndarray)
-    )
 
     def __attrs_post_init__(self):
-        if self.frequency.size == 1:
-            if self.n_components > 1:
-                shape = (self.n_components, 1)
-                self.frequency = np.broadcast_to(self.frequency, shape)
-            else:
-                self.n_components = 1
-
-        self.availability = np.ones((self.frequency.size, N_TIMESTEPS), dtype=float)
+        if self.frequency.size == 1 and self.n_components > 1:
+            shape = (self.n_components, 1)
+            self.frequency = np.broadcast_to(self.frequency, shape)
+        self.n_components = self.frequency.size
+
         self.create_downtime_events()
         self.calculate_availability()
 

From 2bde5e350dc606dc6511883bdc992e827b6ef60b Mon Sep 17 00:00:00 2001
From: "Hammond, Rob" <13874373+RHammond2@users.noreply.github.com>
Date: Thu, 20 Aug 2026 11:23:18 -0700
Subject: [PATCH 17/24] update use_reliability

---
 h2integrate/converters/natural_gas/natural_gas_cc_ct.py | 5 ++---
 1 file changed, 2 insertions(+), 3 deletions(-)

diff --git a/h2integrate/converters/natural_gas/natural_gas_cc_ct.py b/h2integrate/converters/natural_gas/natural_gas_cc_ct.py
index b1d960eb5..1b493d5c0 100644
--- a/h2integrate/converters/natural_gas/natural_gas_cc_ct.py
+++ b/h2integrate/converters/natural_gas/natural_gas_cc_ct.py
@@ -69,7 +69,6 @@ def initialize(self):
         self.commodity_rate_units = "MW"
         self.commodity_amount_units = "MW*h"
         self.reliability_model = None
-        self.use_reliability = False
 
     def setup(self):
         super().setup()
@@ -78,12 +77,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"]:
+        if use_reliability := "reliability" in self.options["tech_config"]["model_inputs"]:
             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
+        self.use_reliability = use_reliability
 
         # Add natural gas consumed output
         self.add_output(

From 33f80a33fa494af0fcc13562775bced33094ecd4 Mon Sep 17 00:00:00 2001
From: "Hammond, Rob" <13874373+RHammond2@users.noreply.github.com>
Date: Thu, 20 Aug 2026 15:43:23 -0700
Subject: [PATCH 18/24] add validator for array variation of validators.ge

---
 h2integrate/reliability/models.py | 27 +++++++++++----------------
 1 file changed, 11 insertions(+), 16 deletions(-)

diff --git a/h2integrate/reliability/models.py b/h2integrate/reliability/models.py
index 9317ed649..5088a391b 100644
--- a/h2integrate/reliability/models.py
+++ b/h2integrate/reliability/models.py
@@ -14,6 +14,14 @@
 N_TIMESTEPS = 8760
 
 
+def array_ge(val):
+    """Validates that all values of an array are greater than or equal to :py:attr:`val`."""
+
+    def validator(instance, attribute, value):
+        if any(value < val):
+            raise ValueError(f"'{attribute.name}' must have all values of at least 0.")
+
+
 def create_failure_model(config: dict):
     """Retrieves and initializes a matching reliability model."""
     name = config["failure_model"]
@@ -149,16 +157,10 @@ class FixedDowntime(BaseDowntime):
 
     hours: int | ArrayLike = field(
         converter=int_array_converter,
-        validator=validators.instance_of(np.ndarray),
+        validator=(validators.instance_of(np.ndarray), array_ge(1)),
     )
     n_components: int = field(default=1, validator=(validators.instance_of(int), validators.ge(1)))
 
-    @hours.validator
-    def hours_validator(self, attribute, value):
-        """Validates that all values of :py:attr:`hours` are greater than or equal to 1."""
-        if not np.all(value > 1):
-            raise ValueError("All values passed to 'hours' must be greater than or equal to 1.")
-
     def __attrs_post_init__(self):
         if self.hours.size == 1 and self.n_components > 1:
             shape = (self.n_components, 1)
@@ -181,21 +183,14 @@ class LogNormalDowntime(BaseDowntime):
 
     mean: float = field(
         converter=float_array_converter,
-        validator=(validators.instance_of(np.ndarray)),
+        validator=(validators.instance_of(np.ndarray), array_ge(0)),
     )
     sigma: float = field(
         converter=float_array_converter,
-        validator=(validators.instance_of(np.ndarray)),
+        validator=(validators.instance_of(np.ndarray), array_ge(0)),
     )
     n_components: int = field(default=1, validator=(validators.instance_of(int), validators.ge(1)))
 
-    @mean.validator
-    @sigma.validator
-    def validate_ge_zero(self, attribute, value):
-        """Validates the :py:attr:`sigma` and :py:attr:`mean` values are at least 0."""
-        if any(value < 0):
-            raise ValueError(f"'{attribute.name}' must have all values of at least 0.")
-
     @sigma.validator
     def validate_shape(self, attribute, value):
         """Validates that :py:attr:`mean` and :py:attr:`sigma` are the same size."""

From 8568a2c13765610b767e487b6f9e1e4f49649b10 Mon Sep 17 00:00:00 2001
From: "Hammond, Rob" <13874373+RHammond2@users.noreply.github.com>
Date: Thu, 20 Aug 2026 16:01:02 -0700
Subject: [PATCH 19/24] standardize dimensionality update

---
 h2integrate/reliability/models.py | 49 +++++++++++++++++++------------
 1 file changed, 31 insertions(+), 18 deletions(-)

diff --git a/h2integrate/reliability/models.py b/h2integrate/reliability/models.py
index 5088a391b..9218783a5 100644
--- a/h2integrate/reliability/models.py
+++ b/h2integrate/reliability/models.py
@@ -62,6 +62,29 @@ def generate_downtime_model(config: dict | int):
             raise NotImplementedError(f"{name} is not a valid model name")
 
 
+def update_dimensions(n_components: int, *args: np.ndarray):
+    """Update the dimensionality of a series of arrays and the value of :py:attr:`n_components`
+    to match. If the size of an array passed :py:attr:`args` is 1 and :py:attr:`n_components`
+    is greater than 1, all arrays passed to :py:attr:`args` will be broadcast to an array shaped
+    (:py:attr:`n_components`, 1). If the arrays are already larger than 1, then
+    :py:attr:`n_components will be updated to the size of the arrays.
+
+    Args:
+        n_components (int): Number of components in the model
+        args (np.ndarray): NumPy array of attribute values used to create a model.
+
+    Returns:
+        n_components: The value as passed or updated to match the size of arrays in :py:attr:`args`.
+        args: The arrays as passed or the arrays reshaped to shape (:py:attr:`n_components`, 1).
+    """
+    if args[0].size == 1 and n_components > 1:
+        shape = (n_components, 1)
+        for arg in args:
+            arg = np.broadcast_to(arg, shape)
+    n_components = args[0].size
+    return n_components, *args
+
+
 @define(kw_only=True)
 class BaseDowntime(ABC, BaseConfig):
     @abstractmethod
@@ -162,10 +185,7 @@ class FixedDowntime(BaseDowntime):
     n_components: int = field(default=1, validator=(validators.instance_of(int), validators.ge(1)))
 
     def __attrs_post_init__(self):
-        if self.hours.size == 1 and self.n_components > 1:
-            shape = (self.n_components, 1)
-            self.hours = np.broadcast_to(self.hours, shape)
-        self.n_components = self.hours.size
+        self.n_components, self.hours = update_dimensions(self.n_components, self.hours)
 
     def sample_downtime(self):
         """Return an array of 100 :py:attr:`hours`."""
@@ -202,11 +222,9 @@ def validate_shape(self, attribute, value):
             raise ValueError(msg)
 
     def __attrs_post_init__(self):
-        if self.mean.size == 1 and self.n_components > 1:
-            shape = (self.n_components, 1)
-            self.mean = np.broadcast_to(self.mean, shape)
-            self.sigma = np.broadcast_to(self.sigma, shape)
-        self.n_components = self.mean.size
+        self.n_components, self.mean, self.sigma = update_dimensions(
+            self.n_components, self.mean, self.sigma
+        )
 
     def sample_downtime(self) -> np.ndarray:
         """Return an array of 100 samples of each lognormal distribution."""
@@ -251,11 +269,9 @@ class WeibullReliability(BaseReliability):
     )
 
     def __attrs_post_init__(self):
-        if self.scale.size == 1 and self.n_components > 1:
-            shape = (self.n_components, 1)
-            self.scale = np.broadcast_to(self.scale, shape)
-            self.shape = np.broadcast_to(self.shape, shape)
-        self.n_components = self.scale.size
+        self.n_components, self.scale, self.shape = update_dimensions(
+            self.n_components, self.scale, self.shape
+        )
 
         self.create_downtime_events()
         self.calculate_availability()
@@ -291,10 +307,7 @@ class FixedIntervalReliability(BaseReliability):
     )
 
     def __attrs_post_init__(self):
-        if self.frequency.size == 1 and self.n_components > 1:
-            shape = (self.n_components, 1)
-            self.frequency = np.broadcast_to(self.frequency, shape)
-        self.n_components = self.frequency.size
+        self.n_components, self.frequency = update_dimensions(self.n_components, self.frequency)
 
         self.create_downtime_events()
         self.calculate_availability()

From 7e9d4730b1cb6928cc0bd0b28267da2f70529dde Mon Sep 17 00:00:00 2001
From: "Hammond, Rob" <13874373+RHammond2@users.noreply.github.com>
Date: Thu, 20 Aug 2026 16:07:36 -0700
Subject: [PATCH 20/24] minor reorg

---
 h2integrate/reliability/models.py | 34 ++++++++++++++++---------------
 1 file changed, 18 insertions(+), 16 deletions(-)

diff --git a/h2integrate/reliability/models.py b/h2integrate/reliability/models.py
index 9218783a5..964c8de3a 100644
--- a/h2integrate/reliability/models.py
+++ b/h2integrate/reliability/models.py
@@ -14,14 +14,6 @@
 N_TIMESTEPS = 8760
 
 
-def array_ge(val):
-    """Validates that all values of an array are greater than or equal to :py:attr:`val`."""
-
-    def validator(instance, attribute, value):
-        if any(value < val):
-            raise ValueError(f"'{attribute.name}' must have all values of at least 0.")
-
-
 def create_failure_model(config: dict):
     """Retrieves and initializes a matching reliability model."""
     name = config["failure_model"]
@@ -85,6 +77,24 @@ def update_dimensions(n_components: int, *args: np.ndarray):
     return n_components, *args
 
 
+def float_array_converter(val: int | float | ArrayLike):
+    return np.array(val).astype(float).reshape(-1, 1)
+
+
+def int_array_converter(val: int | float | ArrayLike):
+    return np.array(val).astype(int).reshape(-1, 1)
+
+
+# NOTE: yes, I see the irony in creating a custom converter after removing all the attrs duplicates
+# NOTE: yes, I'll also be removing the above comment, but I'm sure someone will appreciate it
+def array_ge(val):
+    """Validates that all values of an array are greater than or equal to :py:attr:`val`."""
+
+    def validator(instance, attribute, value):
+        if any(value < val):
+            raise ValueError(f"'{attribute.name}' must have all values of at least 0.")
+
+
 @define(kw_only=True)
 class BaseDowntime(ABC, BaseConfig):
     @abstractmethod
@@ -160,14 +170,6 @@ def calculate_availability(self):
         self.system_availability = np.min(self.availability, axis=0)
 
 
-def float_array_converter(val: int | float | ArrayLike):
-    return np.array(val).astype(float).reshape(-1, 1)
-
-
-def int_array_converter(val: int | float | ArrayLike):
-    return np.array(val).astype(int).reshape(-1, 1)
-
-
 @define(kw_only=True)
 class FixedDowntime(BaseDowntime):
     """Basic log-normal downtime model for generating the length of downtime for a given event.

From b4a0ff50811b9edc37bc8ec0cf6472a4eb830881 Mon Sep 17 00:00:00 2001
From: "Hammond, Rob" <13874373+RHammond2@users.noreply.github.com>
Date: Thu, 20 Aug 2026 16:11:29 -0700
Subject: [PATCH 21/24] add missing return

---
 h2integrate/reliability/models.py | 2 ++
 1 file changed, 2 insertions(+)

diff --git a/h2integrate/reliability/models.py b/h2integrate/reliability/models.py
index 964c8de3a..0d33817a1 100644
--- a/h2integrate/reliability/models.py
+++ b/h2integrate/reliability/models.py
@@ -94,6 +94,8 @@ def validator(instance, attribute, value):
         if any(value < val):
             raise ValueError(f"'{attribute.name}' must have all values of at least 0.")
 
+    return validator
+
 
 @define(kw_only=True)
 class BaseDowntime(ABC, BaseConfig):

From 40ac1c99fb31ddc7df79e72b057cdb0b1387de86 Mon Sep 17 00:00:00 2001
From: "Hammond, Rob" <13874373+RHammond2@users.noreply.github.com>
Date: Thu, 20 Aug 2026 16:17:53 -0700
Subject: [PATCH 22/24] add shape checking validator

---
 h2integrate/reliability/models.py | 21 +++++++++++++++++++--
 1 file changed, 19 insertions(+), 2 deletions(-)

diff --git a/h2integrate/reliability/models.py b/h2integrate/reliability/models.py
index 0d33817a1..744517de8 100644
--- a/h2integrate/reliability/models.py
+++ b/h2integrate/reliability/models.py
@@ -97,6 +97,23 @@ def validator(instance, attribute, value):
     return validator
 
 
+def match_shape(other):
+    """Validates the shape of an array matches the shape of :py:attr:`other`."""
+
+    def validator(instance, attribute, value):
+        other_arr = getattr(instance, other, None)
+        if other_arr is None:
+            raise ValueError(f"'{other}' does not exist.")
+        if other_arr.shape != value.shape:
+            msg = (
+                f"Shape of '{attribute.name}' {value.shape} does not match"
+                f" '{other}' {other_arr.shape}."
+            )
+            raise ValueError(msg)
+
+    return validator
+
+
 @define(kw_only=True)
 class BaseDowntime(ABC, BaseConfig):
     @abstractmethod
@@ -211,7 +228,7 @@ class LogNormalDowntime(BaseDowntime):
     )
     sigma: float = field(
         converter=float_array_converter,
-        validator=(validators.instance_of(np.ndarray), array_ge(0)),
+        validator=(validators.instance_of(np.ndarray), match_shape("mean"), array_ge(0)),
     )
     n_components: int = field(default=1, validator=(validators.instance_of(int), validators.ge(1)))
 
@@ -269,7 +286,7 @@ class WeibullReliability(BaseReliability):
     )
     shape: float = field(
         converter=float_array_converter,
-        validator=validators.instance_of(np.ndarray),
+        validator=(validators.instance_of(np.ndarray), match_shape("scale")),
     )
 
     def __attrs_post_init__(self):

From 823ea1995e9e95d30ec27eb9a6a11b096e5ca42f Mon Sep 17 00:00:00 2001
From: "Hammond, Rob" <13874373+RHammond2@users.noreply.github.com>
Date: Thu, 20 Aug 2026 16:18:37 -0700
Subject: [PATCH 23/24] remove extra shape check

---
 h2integrate/reliability/models.py | 10 ----------
 1 file changed, 10 deletions(-)

diff --git a/h2integrate/reliability/models.py b/h2integrate/reliability/models.py
index 744517de8..a5e91ef26 100644
--- a/h2integrate/reliability/models.py
+++ b/h2integrate/reliability/models.py
@@ -232,16 +232,6 @@ class LogNormalDowntime(BaseDowntime):
     )
     n_components: int = field(default=1, validator=(validators.instance_of(int), validators.ge(1)))
 
-    @sigma.validator
-    def validate_shape(self, attribute, value):
-        """Validates that :py:attr:`mean` and :py:attr:`sigma` are the same size."""
-        if self.mean.shape != value.shape:
-            msg = (
-                "Inputs to 'mean' and 'sigma' must be the same length. Received"
-                f" 'mean': {self.mean.size}, 'sigma': {value.size}"
-            )
-            raise ValueError(msg)
-
     def __attrs_post_init__(self):
         self.n_components, self.mean, self.sigma = update_dimensions(
             self.n_components, self.mean, self.sigma

From 6850587612e8b6375a36b3d9e4d51d1b5ba62b0d Mon Sep 17 00:00:00 2001
From: "Hammond, Rob" <13874373+RHammond2@users.noreply.github.com>
Date: Fri, 21 Aug 2026 09:01:49 -0700
Subject: [PATCH 24/24] move downtime event definition to base class and add
 array_gt

---
 h2integrate/reliability/models.py | 37 +++++++++++++++----------------
 1 file changed, 18 insertions(+), 19 deletions(-)

diff --git a/h2integrate/reliability/models.py b/h2integrate/reliability/models.py
index a5e91ef26..6d7eb32a5 100644
--- a/h2integrate/reliability/models.py
+++ b/h2integrate/reliability/models.py
@@ -92,7 +92,17 @@ def array_ge(val):
 
     def validator(instance, attribute, value):
         if any(value < val):
-            raise ValueError(f"'{attribute.name}' must have all values of at least 0.")
+            raise ValueError(f"'{attribute.name}' must have all values of at least {val}.")
+
+    return validator
+
+
+def array_gt(val):
+    """Validates that all values of an array are greater than or equal to :py:attr:`val`."""
+
+    def validator(instance, attribute, value):
+        if any(value <= val):
+            raise ValueError(f"'{attribute.name}' must have all values greater than {val}.")
 
     return validator
 
@@ -121,7 +131,7 @@ def sample_downtime(self) -> np.ndarray: ...
 
 
 @define(kw_only=True)
-class BaseReliability(BaseConfig):
+class BaseReliability(ABC, BaseConfig):
     """Base reliability class responsible for common definitions and functionality.
 
     Args:
@@ -158,6 +168,11 @@ class BaseReliability(BaseConfig):
     def sample_events(self) -> np.ndarray:
         raise NotImplementedError("Failed to successfully subclass, please implement me.")
 
+    def create_downtime_events(self):
+        """Creates the ``time_to_failure`` and ``downtime_per_event`` arrays."""
+        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."""
         burn_in_time = np.ceil(self.burn_in * N_TIMESTEPS).astype(int)
@@ -293,11 +308,6 @@ def sample_events(self):
             rng.weibull(self.shape, size=(self.shape.size, 100)) * self.scale * N_TIMESTEPS
         ).astype(int)
 
-    def create_downtime_events(self):
-        """Creates a ``time_to_failure`` and ``downtime_per_event``."""
-        self.time_to_failures = self.sample_events()
-        self.downtime_per_event = self.downtime.sample_downtime()
-
 
 @define(kw_only=True)
 class FixedIntervalReliability(BaseReliability):
@@ -314,7 +324,7 @@ class FixedIntervalReliability(BaseReliability):
 
     frequency: float = field(
         converter=float_array_converter,
-        validator=validators.instance_of(np.ndarray),
+        validator=(validators.instance_of(np.ndarray), array_gt(0)),
     )
 
     def __attrs_post_init__(self):
@@ -323,12 +333,6 @@ def __attrs_post_init__(self):
         self.create_downtime_events()
         self.calculate_availability()
 
-    @frequency.validator
-    def validate_frequency(self, attribute, value):
-        """Validates all passed values to :py:attr:`frequency` are greater than 0."""
-        if any(value <= 0):
-            raise ValueError("All values of 'frequency' must be greater than 0.")
-
     def sample_events(self):
         """Creates the time to next failure array for each event's modality with the first event
         occurring randomly either in the first year or first failure interval, whichever
@@ -343,8 +347,3 @@ def sample_events(self):
             (first_occurrence, np.broadcast_to(interval, (interval.size, 99)))
         )
         return time_to_failures
-
-    def create_downtime_events(self):
-        """Creates a ``time_to_failure`` and ``downtime_per_event``."""
-        self.time_to_failures = self.sample_events()
-        self.downtime_per_event = self.downtime.sample_downtime()