diff --git a/.github/workflows/typos.toml b/.github/workflows/typos.toml index eb4df77585..0254951c21 100644 --- a/.github/workflows/typos.toml +++ b/.github/workflows/typos.toml @@ -37,6 +37,9 @@ equil = "equil" astroid = "astroid" # delt == delta delt = "delt" +# FOM for fixed operating and maintenance for DISPATCHES +FOM = "FOM" +fom = "fom" # Minimal Infeasible System mis = "mis" MIS = "MIS" diff --git a/docs/reference_guides/apps/grid_integration/index.rst b/docs/reference_guides/apps/grid_integration/index.rst index 08645e3b36..99714c9434 100644 --- a/docs/reference_guides/apps/grid_integration/index.rst +++ b/docs/reference_guides/apps/grid_integration/index.rst @@ -20,6 +20,7 @@ wholesale electricity markets. For more information, please look at the introduc Bidder Tracker Coordinator + multiperiod/index Cite us ^^^^^^^ diff --git a/docs/reference_guides/apps/grid_integration/multiperiod/Price_Taker.rst b/docs/reference_guides/apps/grid_integration/multiperiod/Price_Taker.rst new file mode 100644 index 0000000000..6ae60840eb --- /dev/null +++ b/docs/reference_guides/apps/grid_integration/multiperiod/Price_Taker.rst @@ -0,0 +1,55 @@ +Price Taker +=========== + +Price takers are entities that must accept market prices since they lack the market share +to directly influence the market price. Likewise, it is assumed that a price taker's resource or energy +system is small enough such that it does not significantly impact the market. When coupled with multi-period modeling, +the price-taker model is able to synthesize grid-centric modeling with steady-state process-centric modeling, as +depicted in figure below. + +.. |pricetaker| image:: images/pricetaker.png + :width: 1200 + :alt: Alternative text + :align: middle + +|pricetaker| + +The following equations represent the multi-period price taker model, where :math:`d` are design decisions, +:math:`u` are operating decisions, :math:`δ` are power decisions, :math:`s` are scenarios (timeseries/representative days), +:math:`w` is weight/frequency, :math:`R` is revenue, :math:`π` is price data, +:math:`C` is capital and operating costs, :math:`g` is the process model, and :math:`h` is the temporal constraint. + + .. math:: + + max_{d,u, δ} = \sum_{s ∈ S} \sum_{t ∈ T} w_{s}[R(d,u_{s,t},δ_{s,t},π_{s,t}) - C(d,u_{s,t},δ_{s,t})] + + .. math:: + + g(d,u_{s,t},δ_{s,t}) = 0; ∀_{s} ∈ S, t ∈ T + + .. math:: + + h(d,u_{s,t},δ_{s,t},u_{s,t+1},δ_{s,t+1}) = 0; ∀_{s} ∈ S, t ∈ T + + +The price taker multi-period modeling workflow involves the integration of multiple software platforms into the IDAES optimization model +and can be broken down into two distinct functions, as shown in the figure below. In part 1, simulated or historical +ISO (Independent System Operator) data is used to generate locational marginal price (LMP) +signals, and production cost models (PCMs) are used to compute and optimize the time-varying dispatch schedules for each +resource based on their respective bid curves. Advanced data analytics (RAVEN) reinterpret the LMP signals and PCM +as stochastic realizations of the LMPs in the form of representative days (or simply the full-year price signals). +In part 2, PRESCIENT uses a variety of input parameters (design capacity, minimum power output, ramp rate, minimum up/down time, marginal cost, no load cost, and startup profile) +to generate data for the market surrogates. Meanwhile, IDAES uses the double loop simulation to integrate detailed +process models (b, ii) into the daily (a, c) and hourly (i, iii) grid operations workflow. + +.. |hybrid_energy_system| image:: images/hybrid_energy_system.png + :width: 1200 + :alt: Alternative text + :align: middle + +|hybrid_energy_system| + +.. module:: idaes.apps.grid_integration.pricetaker.price_taker_model + +.. autoclass:: PriceTakerModel + :members: \ No newline at end of file diff --git a/docs/reference_guides/apps/grid_integration/multiperiod/images/hybrid_energy_system.png b/docs/reference_guides/apps/grid_integration/multiperiod/images/hybrid_energy_system.png new file mode 100644 index 0000000000..119e862d94 Binary files /dev/null and b/docs/reference_guides/apps/grid_integration/multiperiod/images/hybrid_energy_system.png differ diff --git a/docs/reference_guides/apps/grid_integration/multiperiod/images/pricetaker.png b/docs/reference_guides/apps/grid_integration/multiperiod/images/pricetaker.png new file mode 100644 index 0000000000..67abeb1e9a Binary files /dev/null and b/docs/reference_guides/apps/grid_integration/multiperiod/images/pricetaker.png differ diff --git a/docs/reference_guides/apps/grid_integration/multiperiod/index.rst b/docs/reference_guides/apps/grid_integration/multiperiod/index.rst new file mode 100644 index 0000000000..2d2b972b90 --- /dev/null +++ b/docs/reference_guides/apps/grid_integration/multiperiod/index.rst @@ -0,0 +1,13 @@ +Multi-Period Modeling +===================== + +Multi-period modeling can be used to simplify dynamic optimization problems +by linking steady-state models over a time horizon with rate-limiting constraints. +Therefore, sets of constraints at one temporal index will affect decisions taken +at a different moment in time. These interactions can be used to model the relationship +between the energy systems and wholesale electricity markets. + +.. toctree:: + :maxdepth: 2 + + Price_Taker diff --git a/idaes/apps/grid_integration/__init__.py b/idaes/apps/grid_integration/__init__.py index 5d8c6f802e..28a2f92fbc 100644 --- a/idaes/apps/grid_integration/__init__.py +++ b/idaes/apps/grid_integration/__init__.py @@ -15,3 +15,8 @@ from .coordinator import DoubleLoopCoordinator from .forecaster import PlaceHolderForecaster from .multiperiod.multiperiod import MultiPeriodModel +from .pricetaker.price_taker_model import PriceTakerModel +from .pricetaker.design_and_operation_models import ( + DesignModel, + OperationModel, +) diff --git a/idaes/apps/grid_integration/pricetaker/__init__.py b/idaes/apps/grid_integration/pricetaker/__init__.py new file mode 100644 index 0000000000..e69de29bb2 diff --git a/idaes/apps/grid_integration/pricetaker/clustering.py b/idaes/apps/grid_integration/pricetaker/clustering.py new file mode 100644 index 0000000000..2e1cc2758e --- /dev/null +++ b/idaes/apps/grid_integration/pricetaker/clustering.py @@ -0,0 +1,252 @@ +################################################################################# +# The Institute for the Design of Advanced Energy Systems Integrated Platform +# Framework (IDAES IP) was produced under the DOE Institute for the +# Design of Advanced Energy Systems (IDAES). +# +# Copyright (c) 2018-2023 by the software owners: The Regents of the +# University of California, through Lawrence Berkeley National Laboratory, +# National Technology & Engineering Solutions of Sandia, LLC, Carnegie Mellon +# University, West Virginia University Research Corporation, et al. +# All rights reserved. Please see the files COPYRIGHT.md and LICENSE.md +# for full copyright and license information. +################################################################################# + +""" +Contains functions for clustering the price signal into representative days/periods. +""" + +from typing import Union +import matplotlib.pyplot as plt +import numpy as np +import pandas as pd +from pyomo.common.dependencies import attempt_import +import idaes.logger as idaeslog + +sklearn, sklearn_avail = attempt_import("sklearn") + +if sklearn_avail: + from sklearn.cluster import KMeans + from sklearn.metrics import silhouette_score + +_logger = idaeslog.getLogger(__name__) + + +def generate_daily_data(raw_data: list, horizon_length: int): + """ + Function used to generate the daily data in a usable format + from the raw data provided. + + Args: + raw_data : list, + Columnar data for a given LMP signal + + horizon_length : int, + Length of each representative day/period + + Returns: + daily_data: pd.DataFrame, + Correctly formatted daily LMP data for later use + """ + # Converting the data to a list (Needed if the input is not a list, e.g., pd.DataFrame) + raw_data = list(raw_data) + + if horizon_length > len(raw_data): + raise ValueError( + f"Horizon length {horizon_length} exceeds the price signal length of {len(raw_data)}" + ) + + elements_ignored = len(raw_data) % horizon_length + if elements_ignored: + _logger.warning( + f"Length of the price signal is not an integer multiple of horizon_length.\n" + f"\tIgnoring the last {elements_ignored} elements in the price signal." + ) + + daily_data = { + j: raw_data[((j - 1) * horizon_length) : (j * horizon_length)] + for j in range(1, (len(raw_data) // horizon_length) + 1) + } + + # DataFrame arranges the data for each day as a column. Since most clustering techniques + # require different samples as rows, we take the transpose before returning the data. + return pd.DataFrame(daily_data).transpose() + + +def cluster_lmp_data( + raw_data: list, horizon_length: int, n_clusters: int, seed: int = 42 +): + """ + Clusters the given price signal into n_clusters using the k-means clustering + technique. + + Args: + raw_data: list, + Columnar data for a given LMP signal. + + horizon_length: int, + Length of each cluster (representative day/period) + + n_clusters: int, + Number of clusters desired for the data (representative days). + + seed: int, + Seed value for initializing random number generator within Kmeans + + Returns: + lmp_data_clusters: dict + A dictionary of representative day LMP data, indices are indexed + by integers starting at 1. Example: :: + + {1: {1: 4, 2: 3, 3: 5}, + 2: {1: 1, 2: 7, 3: 3}} + + + weights: dict + A dictionary of weights for each representative day, indexed the + same way as lmp_data. Example: :: + + {1: 45, 2: 56} + """ + # reconfiguring raw data + daily_data = generate_daily_data(raw_data, horizon_length) + + # KMeans clustering with the optimal number of clusters + kmeans = KMeans(n_clusters=n_clusters, random_state=seed).fit(daily_data) + centroids = kmeans.cluster_centers_ + labels = kmeans.labels_ + + # Set any centroid values that are < 1e-4 to 0 to avoid noise + centroids = centroids * (abs(centroids) >= 1e-4) + + # Create dicts for lmp data and the weight of each cluster + # By default, the data is of type numpy.int or numpy.float. + # Converting the data to python int/float. Otherwise, Pyomo complains! + num_days, num_time_periods = centroids.shape + lmp_data_clusters = { + d + 1: {t + 1: float(centroids[d, t]) for t in range(num_time_periods)} + for d in range(num_days) + } + weights = {d + 1: int(sum(labels == d)) for d in range(num_days)} + + return lmp_data_clusters, weights + + +def get_optimal_num_clusters( + samples: Union[pd.DataFrame, np.array], + kmin: int = 2, + kmax: int = 30, + method: str = "silhouette", + generate_elbow_plot: bool = False, + seed: int = 42, +): + """ + Determines the appropriate number of clusters needed for a + given price signal. + + Args: + samples: pd.DataFrame | np.array, + Set of points with rows containing samples and columns + containing features + + kmin : int, + Minimum number of clusters + + kmax : int, + Maximum number of clusters + + method : str, + Method for obtaining the optimal number of clusters. + Supported methods are elbow and silhouette + + generate_elbow_plot : bool, + If True, generates an elbow plot for inertia as a function of + number of clusters + + seed : int, + Seed value for random number generator + + Returns: + n_clusters: int, + The optimal number of clusters for the given data + """ + if kmin >= kmax: + raise ValueError(f"kmin must be less than kmax, but {kmin} >= {kmax}") + + # For silhouette method, kmin must be 2. So, we require kmin >= 2 in general + if kmin <= 1: + raise ValueError(f"kmin must be > 1. Received kmin = {kmin}") + + if kmax >= len(samples): + raise ValueError(f"kmax must be < len(samples). Received kmax = {kmax}") + + k_values = list(range(kmin, kmax + 1)) + inertia_values = [] + mean_silhouette = [] + + for k in k_values: + kmeans = KMeans(n_clusters=k, random_state=seed).fit(samples) + inertia_values.append(kmeans.inertia_) + + if method == "silhouette": + # Calculate the average silhouette score, if the chosen method + # is silhouette + mean_silhouette.append(silhouette_score(samples, kmeans.labels_)) + + # Identify the optimal number of clusters + if method == "elbow": + raise NotImplementedError( + "elbow method is not supported currently for finding the optimal " + "number of clusters." + ) + + elif method == "silhouette": + max_index = mean_silhouette.index(max(mean_silhouette)) + n_clusters = k_values[max_index] + + else: + raise ValueError( + f"Unrecognized method {method} for optimal number of clusters." + f"\tSupported methods include elbow and silhouette." + ) + + _logger.info(f"Optimal number of clusters is: {n_clusters}") + + if n_clusters + 2 >= kmax: + _logger.warning( + f"Optimal number of clusters is close to kmax: {kmax}. " + f"Consider increasing kmax." + ) + + if generate_elbow_plot: + plt.plot(k_values, inertia_values) + plt.axvline(x=n_clusters, color="red", linestyle="--", label="Elbow") + plt.xlabel("Number of clusters") + plt.ylabel("Inertia") + plt.title("Elbow Method") + plt.xlim(kmin, kmax) + plt.show() + + return n_clusters + + +def locate_elbow(x: list, y: list): + """ + Identifies the elbow/knee for the input/output data + + Args: + x : list + List of independent variables + y : list + List of dependent variables + + Returns: + opt_x : float + Optimal x at which curvature changes significantly + """ + # The implementation is based on + # Ville Satopaa, Jeannie Albrecht, David Irwin, Barath Raghavan + # Finding a “Kneedle” in a Haystack: + # Detecting Knee Points in System Behavior + # https://raghavan.usc.edu/papers/kneedle-simplex11.pdf + + return len(np.array(x)) + len(np.array(y)) diff --git a/idaes/apps/grid_integration/pricetaker/design_and_operation_models.py b/idaes/apps/grid_integration/pricetaker/design_and_operation_models.py new file mode 100644 index 0000000000..83b8d85acf --- /dev/null +++ b/idaes/apps/grid_integration/pricetaker/design_and_operation_models.py @@ -0,0 +1,226 @@ +################################################################################# +# The Institute for the Design of Advanced Energy Systems Integrated Platform +# Framework (IDAES IP) was produced under the DOE Institute for the +# Design of Advanced Energy Systems (IDAES). +# +# Copyright (c) 2018-2023 by the software owners: The Regents of the +# University of California, through Lawrence Berkeley National Laboratory, +# National Technology & Engineering Solutions of Sandia, LLC, Carnegie Mellon +# University, West Virginia University Research Corporation, et al. +# All rights reserved. Please see the files COPYRIGHT.md and LICENSE.md +# for full copyright and license information. +################################################################################# + +from pyomo.environ import Binary, Param, Var +from pyomo.common.config import Bool, ConfigDict, ConfigValue +from idaes.core.base.process_base import declare_process_block_class +from idaes.core.base.process_base import ProcessBlockData +import idaes.logger as idaeslog + +_logger = idaeslog.getLogger(__name__) + + +# pylint: disable = attribute-defined-outside-init, too-many-ancestors +# pylint: disable = invalid-name, logging-fstring-interpolation +@declare_process_block_class("DesignModel") +class DesignModelData(ProcessBlockData): + """ + Builds the 'design model' for a unit/process. + + Args: + model_func: Function that builds the design model + model_args: Dictionary containing the arguments needed for model_func + + The class defines `install_unit` binary variable that takes the value 1 + if the unit is built/installed, and 0 otherwise. + + Function model_func must declare all the necessary design variables, + relations among design variables, capital cost correlations, and fixed O&M + cost correlations. The function must also define attributes `capex` for + capital cost, and `fom` for fixed O&M cost. If not defined, these attributes + will be set to zero. + + Example Usage: + + .. code-block:: python + + def my_design_model(m, p_min, p_max, cost): + m.power = Var() + m.min_capacity = Constraint( + expr=p_min * m.install_unit <= m.power + ) + m.max_capacity = Constraint( + expr=m.power <= p_max * m.install_unit + ) + + # capex and fom must either be a constant, or Var, or Expression + m.capex = Expression(expr=cost["capex"] * m.power) + m.fom = Expression(expr=cost["fom"] * m.power) + + m = ConcreteModel() + m.unit_1 = DesignModel( + model_func=my_design_model, + model_args={ + "p_min": 150, "p_max": 600, "cost": {"capex": 10, "fom": 1}, + }, + ) + """ + + CONFIG = ConfigDict() + CONFIG.declare( + "model_func", + ConfigValue( + doc="Function that builds the design model", + ), + ) + CONFIG.declare( + "model_args", + ConfigValue( + default={}, + doc="Dictionary containing arguments needed for model_func", + ), + ) + + def build(self): + super().build() + + self.install_unit = Var( + within=Binary, + doc="Binary: 1, if the unit is installed, 0 otherwise", + ) + + if self.config.model_func is None: + # Function that builds the design model is not specified + _logger.warning( + "The function that builds the design model is not specified." + "model_func must declare all the necessary design variables," + "relations among design variables, capital cost correlations," + "and fixed operating and maintenance cost correlations." + ) + return + + # Call the function that builds the design model + self.config.model_func(self, **self.config.model_args) + + # Check if capital and fixed O&M costs are defined + if not hasattr(self, "capex"): + _logger.warning( + "'capex' attribute is not set for the design model " + f"{self.name}. Setting the capital cost of the unit to zero." + ) + self.capex = 0 + + if not hasattr(self, "fom"): + _logger.warning( + "'fom' attribute is not set for the design model " + f"{self.name}. Setting the fixed O&M cost of the unit to zero." + ) + self.fom = 0 + + +@declare_process_block_class("OperationModel") +class OperationModelData(ProcessBlockData): + """ + Builds the 'operation model' for a unit/process. + + Args: + model_func: Function that builds the operation model + model_args: Dictionary containing the arguments needed for model_func + + The class defines `op_mode`, `startup`, and `shutdown` binary variables + to track the operation, startup, and shutdown of the unit/process. + + Function model_func must declare all the necessary operation variables, + relations among operation variables, and variable O&M cost correlations. + + Example Usage: + + .. code-block:: python + + def my_operation_model(m, design_blk): + m.power = Var() + m.fuel_flow = Var() + ... + + m = ConcreteModel() + m.unit_1 = DesignModel( + model_func=my_design_model, + model_args={ + "p_min": 150, "p_max": 600, "cost": {"capex": 10, "fom": 1}, + }, + ) + m.op_unit_1 = OperationModel( + model_func=my_operation_model, model_args={"design_blk": m.unit_1}, + ) + """ + + CONFIG = ConfigDict() + CONFIG.declare( + "model_func", + ConfigValue( + doc="Function that builds the design model", + ), + ) + CONFIG.declare( + "model_args", + ConfigValue( + default={}, + doc="Dictionary containing arguments needed for model_func", + ), + ) + CONFIG.declare( + "declare_op_vars", + ConfigValue( + default=True, + domain=Bool, + doc="Should op_mode, startup, shutdown vars be defined?", + ), + ) + CONFIG.declare( + "declare_lmp_param", + ConfigValue( + default=True, + domain=Bool, + doc="Should LMP data automatically be appended to the model?", + ), + ) + + # noinspection PyAttributeOutsideInit + def build(self): + super().build() + + # Build the model + if self.config.declare_op_vars: + self.op_mode = Var( + within=Binary, + doc="Binary: 1 if the unit is operating, 0 otherwise", + ) + + self.startup = Var( + within=Binary, + doc="Binary: 1 if the startup is initiated, 0 otherwise", + ) + + self.shutdown = Var( + within=Binary, + doc="Binary: 1 if the shutdown is initiated, 0 otherwise", + ) + + if self.config.declare_lmp_param: + self.LMP = Param( + initialize=1, + mutable=True, + doc="Time-varying locational marginal prices (LMPs) [in $/MWh]", + ) + + if self.config.model_func is None: + _logger.warning( + "The function that builds the operation model is not specified." + "model_func must declare all the necessary operation variables," + "relations among operation variables, and variable" + "operating and maintenance cost correlations." + ) + return + + # Call the function that builds the operation model + self.config.model_func(self, **self.config.model_args) diff --git a/idaes/apps/grid_integration/pricetaker/price_taker_model.py b/idaes/apps/grid_integration/pricetaker/price_taker_model.py new file mode 100644 index 0000000000..d3a6b3b056 --- /dev/null +++ b/idaes/apps/grid_integration/pricetaker/price_taker_model.py @@ -0,0 +1,950 @@ +################################################################################# +# The Institute for the Design of Advanced Energy Systems Integrated Platform +# Framework (IDAES IP) was produced under the DOE Institute for the +# Design of Advanced Energy Systems (IDAES). +# +# Copyright (c) 2018-2023 by the software owners: The Regents of the +# University of California, through Lawrence Berkeley National Laboratory, +# National Technology & Engineering Solutions of Sandia, LLC, Carnegie Mellon +# University, West Virginia University Research Corporation, et al. +# All rights reserved. Please see the files COPYRIGHT.md and LICENSE.md +# for full copyright and license information. +################################################################################# + +from typing import Optional, Union, Callable, List +import pandas as pd +from pyomo.environ import ( + ConcreteModel, + Block, + Var, + Param, + RangeSet, + Objective, + Constraint, + NonNegativeReals, + Expression, + maximize, +) +from pyomo.common.config import ( + ConfigDict, + ConfigValue, + NonNegativeInt, + NonNegativeFloat, + ListOf, + PositiveInt, +) + +from idaes.apps.grid_integration.pricetaker.design_and_operation_models import ( + DesignModelData, +) +from idaes.apps.grid_integration.pricetaker.clustering import ( + generate_daily_data, + cluster_lmp_data, + get_optimal_num_clusters, +) +from idaes.apps.grid_integration.pricetaker.unit_commitment import ( + UnitCommitmentData, + capacity_limits, + ramping_limits, + startup_shutdown_constraints, +) +from idaes.core.util.config import ConfigurationError, is_in_range +import idaes.logger as idaeslog + +_logger = idaeslog.getLogger(__name__) + + +# Defining a ConfigDict to simplify the domain validation of +# arguments needed for all the methods of the PriceTakerModel class +CONFIG = ConfigDict() + +# List of arguments needed for the `append_lmp_data` method +CONFIG.declare( + "seed", + ConfigValue( + default=20, + domain=NonNegativeInt, + doc="Seed value for clustering techniques", + ), +) +CONFIG.declare( + "horizon_length", + ConfigValue( + domain=PositiveInt, + doc="Length of each representative day/period", + ), +) +CONFIG.declare( + "num_clusters_range", + ConfigValue( + default=(5, 30), + domain=ListOf(int, PositiveInt), + doc="Range of number of clusters for generating elbow plot", + ), +) +CONFIG.declare( + "num_clusters", + ConfigValue( + domain=PositiveInt, + doc="Number of clusters/representaive days/periods", + ), +) +CONFIG.declare( + "lmp_data", + ConfigValue( + domain=ListOf(float), + doc="List containing locational marginal prices (LMP)", + ), +) + +# List of arguments needed for adding startup/shutdown constraints +CONFIG.declare( + "up_time", + ConfigValue( + domain=PositiveInt, + doc="Minimum uptime [in hours]", + ), +) +CONFIG.declare( + "down_time", + ConfigValue( + domain=PositiveInt, + doc="Minimum downtime [in hours]", + ), +) + +# List of arguments for NPV calculation +CONFIG.declare( + "lifetime", + ConfigValue( + domain=PositiveInt, + doc="Total lifetime of the system [in years]", + ), +) +CONFIG.declare( + "discount_rate", + ConfigValue( + domain=is_in_range(0, 1), + doc="Discount rate for annualization [fraction]", + ), +) +CONFIG.declare( + "corporate_tax_rate", + ConfigValue( + domain=is_in_range(0, 1), + doc="Effective corporate tax rate [fraction]", + ), +) +CONFIG.declare( + "annualization_factor", + ConfigValue( + domain=is_in_range(0, 1), + doc="Capital cost annualization factor [fraction]", + ), +) +CONFIG.declare( + "cash_inflow_scale_factor", + ConfigValue( + domain=NonNegativeFloat, + doc="Scaling factor for net cash inflow calculations", + ), +) + + +# pylint: disable = too-many-ancestors, too-many-instance-attributes +class PriceTakerModel(ConcreteModel): + """Builds a price-taker model for a given system""" + + def __init__(self, *args, **kwds): + super().__init__(*args, **kwds) + self._config = CONFIG() + self._has_hourly_cashflows = False + self._has_overall_cashflows = False + self._op_blk_uc_data = {} # Save UC data for op. blocks + self._op_blk_uptime_downtime = {} + self._linking_constraint_counter = 1 + + @property + def num_representative_days(self): + """Returns the number of representative days""" + return self._config.num_clusters + + @num_representative_days.setter + def num_representative_days(self, value): + """Setter for the num_representative_days property""" + if self._config.num_clusters is not None: + raise ConfigurationError( + f"num_representative_days is already defined as " + f"{self.num_representative_days} and it cannot be overwritten." + f"\n\tInstantiate a new PriceTakerModel object." + ) + + self._config.num_clusters = value + + @property + def horizon_length(self): + """Returns the length of each representative day""" + if self._config.horizon_length is not None: + return self._config.horizon_length + + _logger.warning( + "Attribute horizon_length is not specified. Using 24 hours " + "as the horizon length." + ) + return 24 + + @horizon_length.setter + def horizon_length(self, value): + """Setter for the horizon_length property""" + if self._config.horizon_length is not None: + raise ConfigurationError( + f"horizon_length is already defined as {self.horizon_length} and " + f"it cannot be overwritten.\n\tInstantiate a new PriceTakerModel object." + ) + + self._config.horizon_length = value + + def _assert_lmp_data_exists(self): + """Raise an error if LMP data does not exist""" + if self._config.lmp_data is None: + raise ConfigurationError( + "LMP data is missing. Please append the LMP data using the " + "`append_lmp_data` method." + ) + + def append_lmp_data( + self, + lmp_data: Union[list, pd.DataFrame, pd.Series], + num_representative_days: Optional[int] = None, + horizon_length: Optional[int] = None, + seed: Optional[int] = 42, + ): + """ + Appends the locational marginal price (LMP) data to the + PriceTakerModel object. If desired, the method can cluster the data + into representative periods before appending the data. + + Args: + lmp_data: Union[list, tuple, pd.DataFrame] + List of locational marginal prices + + num_representative_days: Optional[int], default=None + number of clusters or representative periods. + + horizon_length: Optional[int], default=None + Length of each representative period + + seed: Optional[int], default=42 + Seed value for k-means clustering technique + """ + + # Perform domain validation (ConfigDict performs the validation) + if self._config.lmp_data is None: + assert len(lmp_data) >= 2 # Do not remove this check! + self._config.lmp_data = lmp_data + self._config.num_clusters = num_representative_days + self._config.horizon_length = horizon_length + self._config.seed = seed + + else: + # We do not allow modification of lmp data after it is defined + raise ConfigurationError( + "Attempted to overwrite the LMP data. Instantiate a " + "new PriceTakerModel object to change the LMP data." + ) + + def get_optimal_representative_days( + self, + kmin: int = 4, + kmax: int = 30, + method: str = "silhouette", + generate_elbow_plot: bool = True, + ): + """ + Returns the optimal number of representative days for the + given price signal. + """ + self._assert_lmp_data_exists() + # Domain validation of kmin and kmax + self._config.num_clusters_range = (kmin, kmax) + + daily_data = generate_daily_data(self._config.lmp_data, self.horizon_length) + + return get_optimal_num_clusters( + daily_data, kmin, kmax, method, generate_elbow_plot, self._config.seed + ) + + def _assert_mp_model_exists(self): + """Raise an error if the multiperiod model does not exist""" + if not hasattr(self, "period"): + raise ConfigurationError( + "Unable to find the multiperiod model. Please use the " + "build_multiperiod_model method to construct one." + ) + + def build_multiperiod_model( + self, + flowsheet_func: Callable, + flowsheet_options: Optional[dict] = None, + ): + """ + Builds the multiperiod model. + + Args: + flowsheet_func : Callable + A function that returns an instance of the flowsheet + + flowsheet_options : dict, + Optional arguments needed for `flowsheet_func` + + """ + self._assert_lmp_data_exists() + if hasattr(self, "period"): + # Object may contain a multiperiod model. so raise an error + raise ConfigurationError( + "A multiperiod model might already exist, as the object has " + "`period` attribute." + ) + + lmp_data = self._config.lmp_data + + # pylint: disable=attribute-defined-outside-init + if self.num_representative_days is not None: + # Need to use representative days + rep_days_lmp, rep_days_weights = cluster_lmp_data( + raw_data=lmp_data, + horizon_length=self.horizon_length, + n_clusters=self.num_representative_days, + seed=self._config.seed, + ) + + else: + # Use full year price signal + self.num_representative_days = 1 + self.horizon_length = len(lmp_data) + rep_days_lmp = {1: {t + 1: val for t, val in enumerate(lmp_data)}} + rep_days_weights = {1: 1} + + self.set_days = RangeSet(self.num_representative_days) + self.set_time = RangeSet(self.horizon_length) + self.rep_days_lmp = rep_days_lmp + self.rep_days_weights = rep_days_weights + + flowsheet_blk = ConcreteModel() + if flowsheet_options is None: + flowsheet_options = {} + flowsheet_func(flowsheet_blk, **flowsheet_options) + + # Based on some earlier testing, it is faster to clone the model + # than to build it from scratch. So, instead of using the `rule` + # argument, we are cloning the model and then transferring the + # attributes to the `period` block. + self.period = Block(self.set_days, self.set_time) + for d, t in self.period: + self.period[d, t].transfer_attributes_from(flowsheet_blk.clone()) + + # If the flowsheet model has LMP defined, update it. + if hasattr(self.period[d, t], "LMP"): + self.period[d, t].LMP = self.rep_days_lmp[d][t] + + # Iterate through model to append LMP data if it's been defined + for blk in self.period[d, t].component_data_objects(Block): + if hasattr(blk, "LMP"): + blk.LMP = self.rep_days_lmp[d][t] + + def _get_operation_blocks( + self, + blk_name: str, + attribute_list: list, + ): + """ + Returns a dictionary of operational blocks named 'blk_name'. + In addition, it also checks the existence of the operational + blocks, and the existence of specified attributes. + """ + # Ensure that the multiperiod model exists + self._assert_mp_model_exists() + + # pylint: disable=not-an-iterable + op_blocks = { + d: {t: self.period[d, t].find_component(blk_name) for t in self.set_time} + for d in self.set_days + } + + # NOTE: It is sufficient to perform checks only for one block, because + # the rest of them are clones. + blk = op_blocks[1][1] # This object always exists + + # First, check for the existence of the operational block + if blk is None: + raise AttributeError(f"Operational block {blk_name} does not exist.") + + # Next, check for the existence of attributes. + for attribute in attribute_list: + if not hasattr(blk, attribute): + raise AttributeError( + f"Required attribute {attribute} is not found in " + f"the operational block {blk_name}." + ) + + return op_blocks + + def _get_operation_vars(self, var_name: str): + """ + Returns a dictionary of pointers to the var_name variable located in each flowsheet + instance. If the variable is not present, then an error is raised. + """ + # Ensure that the multiperiod model exists + self._assert_mp_model_exists() + + # pylint: disable=not-an-iterable + op_vars = { + d: {t: self.period[d, t].find_component(var_name) for t in self.set_time} + for d in self.set_days + } + + # NOTE: It is sufficient to perform checks only for one variable + if op_vars[1][1] is None: + raise AttributeError( + f"Variable {var_name} does not exist in the multiperiod model." + ) + + return op_vars + + def add_linking_constraints(self, previous_time_var: str, current_time_var: str): + """ + Adds constraints to relate variables across two consecutive time periods. This method is + usually needed if the system has storage. Using this method, the holdup at the end of the + previous time period can be equated to the holdup at the beginning of the current time + period. + + Args: + previous_time_var : str, + Name of the operational variable at the end of the previous time step + + current_time_var : str, + Name of the operational variable at the beginning of the current time step + + """ + old_time_var = self._get_operation_vars(previous_time_var) + new_time_var = self._get_operation_vars(current_time_var) + + def _rule_linking_constraints(_, d, t): + if t == 1: + return Constraint.Skip + return old_time_var[d][t - 1] == new_time_var[d][t] + + setattr( + self, + "variable_linking_constraints_" + str(self._linking_constraint_counter), + Constraint(self.set_days, self.set_time, rule=_rule_linking_constraints), + ) + _logger.info( + f"Linking constraints are added to the model at " + f"variable_linking_constraints_{self._linking_constraint_counter}" + ) + self._linking_constraint_counter += 1 # Increase the counter for the new pair + + def _retrieve_uc_data(self, op_blk: str, commodity: str, **kwargs): + """ + Retrieves unit commitment data for the given operation + block if it exists. Otherwise, it creates a new object and + returns the object containing the unit commitment data + """ + if (op_blk, commodity) in self._op_blk_uc_data: + self._op_blk_uc_data[op_blk, commodity].update(**kwargs) + + else: + self._op_blk_uc_data[op_blk, commodity] = UnitCommitmentData( + blk_name=op_blk, commodity_name=commodity, **kwargs + ) + + return self._op_blk_uc_data[op_blk, commodity] + + def add_capacity_limits( + self, + op_block_name: str, + commodity: str, + capacity: Union[float, Param, Var, Expression], + op_range_lb: float, + ): + """ + Adds capacity limit constraints of the form: + op_range_lb * capacity * op_mode(t) <= commodity(t) <= capacity * op_mode(t) + ex: 0.3 * P_max * op_mode(t) <= P(t) <= P_max * op_mode(t), + where P(t) is power at time t and op_mode(t) = 1 if the system is operating + at time t; and op_mode(t) = 0, otherwise. + + Args: + op_block_name: str, + Name of the operation model block, ex: ("fs.ngcc") + + commodity: str, + Name of the commodity on the model the capacity constraints + will be applied to, ex: ("total_power") + + capacity: float, or Pyomo Var, Param, or Expression + Maximum capacity on the commodity, ex: (650.0, or m.min_power_capacity) + + op_range_lb: float, + Ratio of the capacity at minimum stable operation to the maximum capacity. + Must be a number in the interval [0, 1] + """ + # _get_operation_blocks method ensures that the operation block exists, and + # the commodity variable also exists. + op_blocks = self._get_operation_blocks( + blk_name=op_block_name, attribute_list=["op_mode", commodity] + ) + # This method performs all the necessary data checks + uc_data = self._retrieve_uc_data( + op_blk=op_block_name, + commodity=commodity, + capacity=capacity, + op_range_lb=op_range_lb, + ) + + # Create a block for storing capacity limit constraints + cap_limit_blk_name = op_block_name.split(".")[-1] + f"_{commodity}_limits" + if hasattr(self, cap_limit_blk_name): + raise ConfigurationError( + f"Attempting to overwrite capacity limits for {commodity} in {op_block_name}." + ) + + setattr(self, cap_limit_blk_name, Block(self.set_days)) + cap_limit_blk = getattr(self, cap_limit_blk_name) + + # pylint: disable = not-an-iterable + for d in self.set_days: + capacity_limits( + blk=cap_limit_blk[d], + op_blocks=op_blocks[d], + uc_data=uc_data, + set_time=self.set_time, + ) + + # Logger info for where constraint is located on the model + _logger.info( + f"Created capacity limit constraints for commodity {commodity} in " + f"operation block {op_block_name} at {cap_limit_blk.name}" + ) + + def add_ramping_limits( + self, + op_block_name: str, + commodity: str, + capacity: Union[float, Param, Var, Expression], + startup_rate: float, + shutdown_rate: float, + rampup_rate: float, + rampdown_rate: float, + ): + """ + Adds ramping constraints of the form: + ramping_var[t] - ramping_var[t-1] <= + startup_rate * capacity * startup[t] + rampup_rate * capacity * op_mode[t-1]; + ramping_var[t-1] - ramping_var[t] <= + shutdown_rate * capacity * shutdown[t] + rampdown_rate * capacity * op_mode[t] + + Args: + op_block_name: str, + Name of the operation model block, e.g., ("fs.ngcc") + + commodity: str, + Name of the variable that the ramping constraints will be applied to, + e.g., "power" + + capacity: float, or Pyomo Var, or Param, or Expression + String of the name of the entity on the model the ramping constraints + will be applied to, ex: ("total_power") + + startup_rate: float, + Fraction of the maximum capacity that variable ramping_var can + increase during startup (between 0 and 1) + + shutdown_rate: float, + Fraction of the maximum capacity that variable ramping_var can + decrease during shutdown (between 0 and 1) + + rampup_rate: float, + Fraction of the maximum capacity that variable ramping_var can + increase during operation (between 0 and 1) + + rampdown_rate: float, + Fraction of the maximum capacity that variable ramping_var can + decrease during operation (between 0 and 1) + """ + # Get operational blocks + op_blocks = self._get_operation_blocks( + blk_name=op_block_name, + attribute_list=["op_mode", "startup", "shutdown", commodity], + ) + + # Perform all the necessary datachecks + uc_data = self._retrieve_uc_data( + op_blk=op_block_name, + commodity=commodity, + capacity=capacity, + startup_rate=startup_rate, + shutdown_rate=shutdown_rate, + rampup_rate=rampup_rate, + rampdown_rate=rampdown_rate, + ) + uc_data.assert_ramping_args_present() + + # Creating the pyomo block + ramp_blk_name = op_block_name.split(".")[-1] + f"_{commodity}_ramping" + if hasattr(self, ramp_blk_name): + raise ConfigurationError( + f"Attempting to overwrite ramping limits for {commodity} in {op_block_name}." + ) + + setattr(self, ramp_blk_name, Block(self.set_days)) + ramp_blk = getattr(self, ramp_blk_name) + + # pylint: disable=not-an-iterable + for d in self.set_days: + ramping_limits( + blk=ramp_blk[d], + op_blocks=op_blocks[d], + uc_data=uc_data, + set_time=self.set_time, + ) + + # Logger info for where constraint is located on the model + _logger.info( + f"Created ramping constraints for variable {commodity} " + f"on operational block {op_block_name} at {ramp_blk.name}" + ) + + def add_startup_shutdown( + self, + op_block_name: str, + des_block_name: Optional[str] = None, + up_time: int = 1, + down_time: int = 1, + ): + """ + Adds minimum uptime/downtime constraints for a given unit/process + + Args: + op_block_name: str, + Name of the operation model block, e.g., "fs.ngcc" + + des_block_name: str, default=None, + Name of the design model block for the operation block + op_block_name. This argument is specified if the design is + being optimized simultaneously, e.g., "ngcc_design" + + up_time: int, default=1, + Total uptime (must be >= 1), e.g., 4 time periods + Uptime must include the minimum uptime and the time required + for shutdown. + + down_time: int, default=1, + Total downtime (must be >= 1), e.g., 4 time periods + Downtime must include the minimum downtime and the time + required for startup + """ + # Check up_time and down_time for validity + self.config.up_time = up_time + self.config.down_time = down_time + + op_blocks = self._get_operation_blocks( + blk_name=op_block_name, + attribute_list=["op_mode", "startup", "shutdown"], + ) + if des_block_name is not None: + install_unit = self.find_component(des_block_name + ".install_unit") + + if install_unit is None: + raise AttributeError( + f"Binary variable associated with unit installation is not found " + f"in {des_block_name}. \n\tDo not specify des_block_name argument if " + f"installation of the unit is not a decision variable." + ) + else: + install_unit = 1 + + start_shut_blk_name = op_block_name.split(".")[-1] + "_startup_shutdown" + if hasattr(self, start_shut_blk_name): + raise ConfigurationError( + f"Attempting to overwrite startup/shutdown constraints " + f"for operation block {op_block_name}." + ) + + setattr(self, start_shut_blk_name, Block(self.set_days)) + start_shut_blk = getattr(self, start_shut_blk_name) + + # pylint: disable=not-an-iterable + for d in self.set_days: + startup_shutdown_constraints( + blk=start_shut_blk[d], + op_blocks=op_blocks[d], + install_unit=install_unit, + up_time=up_time, + down_time=down_time, + set_time=self.set_time, + ) + + # Save the uptime and downtime data for reference + self._op_blk_uptime_downtime[op_block_name] = { + "up_time": up_time, + "down_time": down_time, + } + + # Logger info for where constraint is located on the model + _logger.info( + f"Created startup/shutdown constraints for operation model " + f" {op_block_name} at {start_shut_blk.name}." + ) + + def add_hourly_cashflows( + self, + revenue_streams: Optional[List] = None, + operational_costs: Optional[List] = None, + ): + """ + Adds an expression for the net cash inflow for each flowsheet + instance. Operational costs in each flowsheet instance may include + 'non_fuel_vom' (non-fuel variable operating costs), 'fuel_cost' + (cost of fuel), and 'carbon_price' (cost associated with producing + carbon; i.e., a carbon tax). Revenue streams in each flowsheet instance + may include the revenue from the wholesale electricity market, revenue + from alternative products (e.g., hydrogen), etc. + + The net cash inflow is calculated as: + + Sum(revenue streams) - Sum(costs) + + for each time period. + + Args: + revenue_streams: List of strings representing the names of the + revenue streams, default: None + + Example: :: + + ['elec_revenue', 'H2_revenue', ] + + costs: List of strings representing the names of the + costs associated with operating at a time period. + default: None + Example: :: + + ['hourly_fixed_cost', 'electricity_cost',] + """ + # Ensure that multiperiod model exists + self._assert_mp_model_exists() + + if operational_costs is None: + _logger.warning( + "Argument operational_costs is not specified, so the total " + "operational cost will be set to 0." + ) + operational_costs = [] + + if revenue_streams is None: + _logger.warning( + "Argument revenue_streams is not specified, so the total " + "revenue will be set to 0." + ) + revenue_streams = [] + + for p in self.period: + # Set net profit contribution expressions to 0 + total_cost_expr = 0 + total_revenue_expr = 0 + + for blk in self.period[p].component_data_objects(Block): + + # Add costs for each block. If more than one block, may have + # costs that exist on one block and not on another. (i.e., coproduction) + for cost in operational_costs: + curr_cost = blk.find_component(cost) + total_cost_expr += curr_cost if curr_cost is not None else 0 + + # Add revenue streams for each block. If more than one block, may have + # revenues that exist on one block and not on another. (i.e., coproduction) + for revenue in revenue_streams: + curr_rev = blk.find_component(revenue) + total_revenue_expr += curr_rev if curr_rev is not None else 0 + + # Account for flowsheet-level operational_costs + for cost in operational_costs: + curr_cost = self.period[p].find_component(cost) + total_cost_expr += curr_cost if curr_cost is not None else 0 + + # Account for flowsheet-level revenue_streams + for revenue in revenue_streams: + curr_rev = self.period[p].find_component(revenue) + total_revenue_expr += curr_rev if curr_rev is not None else 0 + + # Set total cost expression + self.period[p].total_hourly_cost = Expression(expr=total_cost_expr) + + # Set total revenue expression + self.period[p].total_hourly_revenue = Expression(expr=total_revenue_expr) + + # Set net cash inflow expression + self.period[p].net_hourly_cash_inflow = Expression( + expr=self.period[p].total_hourly_revenue + - self.period[p].total_hourly_cost + ) + + # Logger info for where constraint is located on the model + self._has_hourly_cashflows = True + _logger.info( + "Created hourly cashflow expressions at:\n\t period[d, t].total_hourly_cost, " + "period[d, t].total_hourly_revenue, period[d, t].net_hourly_cash_inflow." + ) + + def add_overall_cashflows( + self, + lifetime: int = 30, + discount_rate: float = 0.08, + corporate_tax_rate: float = 0.2, + annualization_factor: Optional[float] = None, + cash_inflow_scale_factor: Optional[float] = 1.0, + ): + """ + Builds overall cashflow expressions. + + Args: + lifetime: int, default=30, + Lifetime of the unit/process [in years] + + discount_rate: float, default=0.08, + Discount rate [fraction] for calculating the current value + of the cashflow. It is also used to compute the annualization + factor. Must be between 0 and 1. + + corporate_tax_rate: float, default=0.2, + Fractional value of corporate tax used in NPV calculations. + Must be between 0 and 1. + + annualization_factor: float, default=None, + Annualization factor + + cash_inflow_scale_factor: float, default=1.0, + Scaling factor for the sum of net hourly cashflows. + If the specified price signal is for the full year, then the + value of this argument must be 1.0. However, if the specified + price signal is for a short period, say 1 month, then an + appropriate scaling factor can be specified to compute the value + for a year (12, for the case of 1 month's price signal). + + """ + # Ensure that multiperiod model exists + self._assert_mp_model_exists() + + # Domain validation of input arguments + self.config.lifetime = lifetime + self.config.discount_rate = discount_rate + self.config.corporate_tax = corporate_tax_rate + self.config.annualization_factor = annualization_factor + self.config.cash_inflow_scale_factor = cash_inflow_scale_factor + + if not self._has_hourly_cashflows: + raise ConfigurationError( + "Hourly cashflows are not added to the model. Please run " + "add_hourly_cashflows method before calling the " + "add_overall_cashflows method." + ) + + capex_expr = 0 + fom_expr = 0 + + count_des_blks = 0 + for blk in self.component_data_objects(Block): + if isinstance(blk, DesignModelData): + count_des_blks += 1 + capex_expr += blk.capex + fom_expr += blk.fom + + if count_des_blks == 0: + _logger.warning( + "No design blocks were found, so the overall capital cost " + "(capex) and the fixed O&M cost (fom) are set to 0." + ) + + # pylint: disable=attribute-defined-outside-init + self.cashflows = Block() + cf = self.cashflows + cf.capex = Var(within=NonNegativeReals, doc="Total CAPEX") + cf.capex_calculation = Constraint(expr=cf.capex == capex_expr) + + cf.fom = Var(within=NonNegativeReals, doc="Yearly Fixed O&M Cost") + cf.fom_calculation = Constraint(expr=cf.fom == fom_expr) + + cf.depreciation = Var(within=NonNegativeReals, doc="Yearly depreciation") + cf.depreciation_calculation = Constraint( + expr=cf.depreciation == cf.capex / lifetime + ) + + cf.net_cash_inflow = Var(doc="Net cash inflow") + cf.net_cash_inflow_calculation = Constraint( + expr=cf.net_cash_inflow + == cash_inflow_scale_factor + * sum( + self.rep_days_weights[d] * self.period[d, t].net_hourly_cash_inflow + for d, t in self.period + ) + ) + + cf.corporate_tax = Var(within=NonNegativeReals, doc="Corporate tax") + cf.corporate_tax_calculation = Constraint( + expr=cf.corporate_tax + >= corporate_tax_rate * (cf.net_cash_inflow - cf.fom - cf.depreciation) + ) + + cf.net_profit = Var(doc="Net profit after taxes") + cf.net_profit_calculation = Constraint( + expr=cf.net_profit == cf.net_cash_inflow - cf.fom - cf.corporate_tax + ) + + if annualization_factor is None: + # If the annualization factor is not specified + annualization_factor = discount_rate / ( + 1 - (1 + discount_rate) ** (-lifetime) + ) + + cf.lifetime_npv = Expression( + expr=(1 / annualization_factor) * cf.net_profit - cf.capex + ) + cf.npv = Expression( + expr=cf.net_profit - annualization_factor * cf.capex, + ) + + self._has_overall_cashflows = True + _logger.info(f"Overall cashflows are added to the block {cf.name}") + + def add_objective_function(self, objective_type="npv"): + """ + Appends the objective function to the model. + + Args: + objective_type: str, default="npv", + Supported objective functions are + annualized net present value ("npv"), + net profit ("net_profit"), and + lifetime net present value ("lifetime_npv"). + """ + # pylint: disable = attribute-defined-outside-init + if not self._has_overall_cashflows: + raise ConfigurationError( + "Overall cashflows are not appended. Please run the " + "add_overall_cashflows method." + ) + + try: + self.obj = Objective( + expr=getattr(self.cashflows, objective_type), + sense=maximize, + ) + + except AttributeError as msg: + raise ConfigurationError( + f"{objective_type} is not a supported objective function." + f"Please specify either npv, or lifetime_npv, or net_profit " + f"as the objective_type." + ) from msg diff --git a/idaes/apps/grid_integration/pricetaker/tests/__init__.py b/idaes/apps/grid_integration/pricetaker/tests/__init__.py new file mode 100644 index 0000000000..e69de29bb2 diff --git a/idaes/apps/grid_integration/pricetaker/tests/test_clustering.py b/idaes/apps/grid_integration/pricetaker/tests/test_clustering.py new file mode 100644 index 0000000000..8b711e52ae --- /dev/null +++ b/idaes/apps/grid_integration/pricetaker/tests/test_clustering.py @@ -0,0 +1,281 @@ +################################################################################# +# The Institute for the Design of Advanced Energy Systems Integrated Platform +# Framework (IDAES IP) was produced under the DOE Institute for the +# Design of Advanced Energy Systems (IDAES). +# +# Copyright (c) 2018-2023 by the software owners: The Regents of the +# University of California, through Lawrence Berkeley National Laboratory, +# National Technology & Engineering Solutions of Sandia, LLC, Carnegie Mellon +# University, West Virginia University Research Corporation, et al. +# All rights reserved. Please see the files COPYRIGHT.md and LICENSE.md +# for full copyright and license information. +################################################################################# + +from pathlib import Path +import matplotlib.pyplot as plt +import pandas as pd +import pytest + +from idaes.apps.grid_integration.pricetaker.clustering import ( + generate_daily_data, + get_optimal_num_clusters, + cluster_lmp_data, + sklearn_avail, + locate_elbow, +) +import idaes.logger as idaeslog + +pytest.importorskip("sklearn", reason="sklearn not available") + + +@pytest.fixture(name="sample_data") +def sample_data_fixture(): + """Returns a sample price signal for testing""" + file_path = Path(__file__).parent / "FLECCS_princeton.csv" + data = pd.read_csv(file_path) + return data + + +@pytest.fixture(name="dummy_data") +def dummy_data_fixture(): + """Returns dummy data for clustering""" + # Test Data Set: [0,0], [0,1], [1,1], [1, 0], [10, 0], [10, 1], + # [11, 1], [11, 0], [5,10], [5, 11], [6, 11], [6, 10] + # Centroids for this set must be [0.5, 0.5], [10.5, 0.5], [5.5, 10.5] + data = [ + 0, + 0, + 0, + 1, + 1, + 1, + 1, + 0, + 10, + 0, + 10, + 1, + 11, + 1, + 11, + 0, + 5, + 10, + 5, + 11, + 6, + 11, + 6, + 10, + ] + return data + + +@pytest.mark.unit +def test_daily_data_logger_message1(): + """Tests ValueError in generate_daily_data function""" + with pytest.raises( + ValueError, + match="Horizon length 24 exceeds the price signal length of 4", + ): + generate_daily_data(raw_data=[1, 2, 3, 4], horizon_length=24) + + +@pytest.mark.unit +def test_daily_data_logger_message2(caplog): + """Tests non-integer-multiple length warning in generate_daily_data function""" + + with caplog.at_level(idaeslog.WARNING): + daily_data = generate_daily_data(raw_data=list(range(8)), horizon_length=3) + assert ( + "Length of the price signal is not an integer multiple of horizon_length.\n" + "\tIgnoring the last 2 elements in the price signal." + ) in caplog.text + + assert isinstance(daily_data, pd.DataFrame) + assert daily_data.shape == (2, 3) + + +@pytest.mark.unit +def test_generate_daily_data(): + """Tests the generate_daily_data function""" + + raw_data = list(range(9)) + daily_data = generate_daily_data(raw_data, horizon_length=3) + assert daily_data.shape == (3, 3) + assert daily_data.equals( + pd.DataFrame({1: [0, 1, 2], 2: [3, 4, 5], 3: [6, 7, 8]}).transpose() + ) + + +@pytest.mark.unit +@pytest.mark.skipif( + not sklearn_avail, reason="sklearn (optional dependency) not available" +) +def test_cluster_lmp_data(dummy_data): + """Tests the cluster_lmp_data function""" + + data_clusters, weights = cluster_lmp_data( + raw_data=dummy_data, horizon_length=2, n_clusters=3, seed=42 + ) + assert weights == {1: 4, 2: 4, 3: 4} + assert data_clusters == { + 1: {1: 10.5, 2: 0.5}, + 2: {1: 5.5, 2: 10.5}, + 3: {1: 0.5, 2: 0.5}, + } + + +@pytest.mark.unit +@pytest.mark.skipif( + not sklearn_avail, reason="sklearn (optional dependency) not available" +) +def test_optimal_num_clusters_elbow(dummy_data): + """Tests the get_optimal_num_clusters function""" + samples = generate_daily_data(dummy_data, horizon_length=2) + with pytest.raises( + NotImplementedError, + match=( + "elbow method is not supported currently for finding the optimal " + "number of clusters." + ), + ): + get_optimal_num_clusters( + samples=samples, + kmin=2, + kmax=7, + method="elbow", + generate_elbow_plot=True, + seed=20, + ) + + +@pytest.mark.unit +@pytest.mark.skipif( + not sklearn_avail, reason="sklearn (optional dependency) not available" +) +def test_optimal_num_clusters_silhouette(dummy_data): + """Tests the get_optimal_num_clusters function""" + samples = generate_daily_data(dummy_data, horizon_length=2) + n_clusters = get_optimal_num_clusters( + samples=samples, + kmin=2, + kmax=7, + method="silhouette", + generate_elbow_plot=True, + ) + + assert n_clusters == 3 + + # Test that a figure was created + assert plt.gcf() is not None + # Test that axes were created + assert plt.gca() is not None + # Test that the plot has data + # assert plt.gca().has_data() + + plt.close("all") + + +@pytest.mark.unit +@pytest.mark.skipif( + not sklearn_avail, reason="sklearn (optional dependency) not available" +) +def test_optimal_clusters_logger_message1(dummy_data): + """Tests the error message associated with kmin >= kmax""" + + samples = generate_daily_data(dummy_data, 2) + with pytest.raises( + ValueError, + match="kmin must be less than kmax, but 15 >= 10", + ): + get_optimal_num_clusters( + samples=samples, + kmin=15, + kmax=10, + ) + + +@pytest.mark.unit +@pytest.mark.skipif( + not sklearn_avail, reason="sklearn (optional dependency) not available" +) +def test_optimal_clusters_logger_message2(dummy_data): + """Tests the error message associated with kmin < 2""" + + samples = generate_daily_data(dummy_data, 2) + with pytest.raises( + ValueError, + match="kmin must be > 1. Received kmin = 1", + ): + get_optimal_num_clusters( + samples=samples, + kmin=1, + kmax=10, + ) + + +@pytest.mark.unit +@pytest.mark.skipif( + not sklearn_avail, reason="sklearn (optional dependency) not available" +) +def test_optimal_clusters_logger_message3(dummy_data): + """Tests the error message associated with kmax >= len(samples)""" + + samples = generate_daily_data(dummy_data, 2) + with pytest.raises( + ValueError, + match="kmax must be < len\\(samples\\). Received kmax = 12", + ): + get_optimal_num_clusters( + samples=samples, + kmin=2, + kmax=12, + ) + + +@pytest.mark.unit +@pytest.mark.skipif( + not sklearn_avail, reason="sklearn (optional dependency) not available" +) +def test_optimal_clusters_logger_message4(dummy_data): + """Tests the error message associated with unrecognized method""" + + samples = generate_daily_data(dummy_data, 2) + with pytest.raises( + ValueError, + match=( + "Unrecognized method foo for optimal number of clusters." + "\tSupported methods include elbow and silhouette." + ), + ): + get_optimal_num_clusters( + samples=samples, + kmin=2, + kmax=5, + method="foo", + ) + + +@pytest.mark.unit +@pytest.mark.skipif( + not sklearn_avail, reason="sklearn (optional dependency) not available" +) +def test_optimal_clusters_logger_message5(dummy_data, caplog): + """Tests the logger warning associated with n_clusters close to kmax""" + + samples = generate_daily_data(dummy_data, 2) + with caplog.at_level(idaeslog.WARNING): + get_optimal_num_clusters(samples, kmin=2, kmax=5, seed=20) + assert ( + "Optimal number of clusters is close to kmax: 5. Consider increasing kmax." + in caplog.text + ) + + +@pytest.mark.unit +def test_locate_elbow(): + """Tests the locate_elbow function""" + opt_x = locate_elbow([1, 2, 3], [4, 5, 6]) + + assert opt_x == 6 diff --git a/idaes/apps/grid_integration/pricetaker/tests/test_design_and_operation_models.py b/idaes/apps/grid_integration/pricetaker/tests/test_design_and_operation_models.py new file mode 100644 index 0000000000..6dd50459bd --- /dev/null +++ b/idaes/apps/grid_integration/pricetaker/tests/test_design_and_operation_models.py @@ -0,0 +1,154 @@ +################################################################################# +# The Institute for the Design of Advanced Energy Systems Integrated Platform +# Framework (IDAES IP) was produced under the DOE Institute for the +# Design of Advanced Energy Systems (IDAES). +# +# Copyright (c) 2018-2023 by the software owners: The Regents of the +# University of California, through Lawrence Berkeley National Laboratory, +# National Technology & Engineering Solutions of Sandia, LLC, Carnegie Mellon +# University, West Virginia University Research Corporation, et al. +# All rights reserved. Please see the files COPYRIGHT.md and LICENSE.md +# for full copyright and license information. +################################################################################# + +import pytest +from pyomo.environ import ConcreteModel, Constraint, Expression, Var +from idaes.apps.grid_integration.pricetaker.design_and_operation_models import ( + DesignModel, + OperationModel, +) + + +@pytest.mark.unit +def test_design_model_class(caplog): + """Tests the DesignModel class""" + + def _my_design_model(m, p_min, p_max, cost): + m.power = Var() + m.min_capacity = Constraint(expr=p_min * m.install_unit <= m.power) + m.max_capacity = Constraint(expr=m.power <= p_max * m.install_unit) + + # capex and fom must either be a constant, or Var, or Expression + m.capex = Expression(expr=cost["capex"] * m.power) + m.fom = Expression(expr=cost["fom"] * m.power) + + blk = ConcreteModel() + blk.unit_1 = DesignModel( + model_func=_my_design_model, + model_args={ + "p_min": 150, + "p_max": 600, + "cost": {"capex": 10, "fom": 1}, + }, + ) + + # Begin checks + for attr in ["install_unit", "power", "min_capacity", "max_capacity"]: + assert hasattr(blk.unit_1, attr) + + assert isinstance(blk.unit_1.capex, Expression) + assert isinstance(blk.unit_1.fom, Expression) + assert "Setting the capital cost of the unit to zero." not in caplog.text + assert "Setting the fixed O&M cost of the unit to zero." not in caplog.text + + # Test empty design model + blk.unit_2 = DesignModel() + assert hasattr(blk.unit_2, "install_unit") + for attr in ["power", "min_capacity", "max_capacity", "capex", "fom"]: + assert not hasattr(blk.unit_2, attr) + + assert "Setting the capital cost of the unit to zero." not in caplog.text + assert "Setting the fixed O&M cost of the unit to zero." not in caplog.text + + # Test warning messages + def dummy_func(_): + pass + + blk.unit_3 = DesignModel(model_func=dummy_func) + assert hasattr(blk.unit_3, "install_unit") + assert ( + "'capex' attribute is not set for the design model " + "unit_3. Setting the capital cost of the unit to zero." + ) in caplog.text + assert ( + "'fom' attribute is not set for the design model " + "unit_3. Setting the fixed O&M cost of the unit to zero." + ) in caplog.text + assert blk.unit_3.capex == 0 + assert blk.unit_3.fom == 0 + + +@pytest.mark.unit +def test_design_model_class_logger_message1(caplog): + caplog.clear() + + blk = ConcreteModel() + blk.unit_1 = DesignModel( + model_args={ + "p_min": 150, + "p_max": 600, + "cost": {"capex": 10, "fom": 1}, + }, + ) + + assert ( + "The function that builds the design model is not specified." + "model_func must declare all the necessary design variables," + "relations among design variables, capital cost correlations," + "and fixed operating and maintenance cost correlations." in caplog.text + ) + + +@pytest.mark.unit +def test_operation_model_class(): + """Tests the OperationModel class""" + + def des_model(m): + m.power = Var() + m.capex = 0 + m.fom = 0 + + def op_model(m, des_blk): + m.power = Var() + m.scaled_power = Expression( + expr=m.power / des_blk.power, + doc="Instantaneous power scaled with the design power", + ) + + blk = ConcreteModel() + blk.unit_1_design = DesignModel(model_func=des_model) + blk.unit_1_op = OperationModel( + model_func=op_model, model_args={"des_blk": blk.unit_1_design} + ) + + # Begin checks + for attr in ["op_mode", "startup", "shutdown", "power", "LMP"]: + assert hasattr(blk.unit_1_op, attr) + + assert isinstance(blk.unit_1_op.scaled_power, Expression) + + # Empty operation model + blk.unit_2_op = OperationModel(declare_op_vars=False, declare_lmp_param=False) + for attr in ["op_mode", "startup", "shutdown", "power", "LMP"]: + assert not hasattr(blk.unit_2_op, attr) + + +@pytest.mark.unit +def test_operation_model_class_logger_message1(caplog): + caplog.clear() + + def des_model(m): + m.power = Var() + m.capex = 0 + m.fom = 0 + + blk = ConcreteModel() + blk.unit_1_design = DesignModel(model_func=des_model) + blk.unit_1_op = OperationModel(model_args={"des_blk": blk.unit_1_design}) + + assert ( + "The function that builds the operation model is not specified." + "model_func must declare all the necessary operation variables," + "relations among operation variables, and variable" + "operating and maintenance cost correlations." in caplog.text + ) diff --git a/idaes/apps/grid_integration/pricetaker/tests/test_price_taker_model.py b/idaes/apps/grid_integration/pricetaker/tests/test_price_taker_model.py new file mode 100644 index 0000000000..ac1dfc26de --- /dev/null +++ b/idaes/apps/grid_integration/pricetaker/tests/test_price_taker_model.py @@ -0,0 +1,1024 @@ +################################################################################# +# The Institute for the Design of Advanced Energy Systems Integrated Platform +# Framework (IDAES IP) was produced under the DOE Institute for the +# Design of Advanced Energy Systems (IDAES). +# +# Copyright (c) 2018-2023 by the software owners: The Regents of the +# University of California, through Lawrence Berkeley National Laboratory, +# National Technology & Engineering Solutions of Sandia, LLC, Carnegie Mellon +# University, West Virginia University Research Corporation, et al. +# All rights reserved. Please see the files COPYRIGHT.md and LICENSE.md +# for full copyright and license information. +################################################################################# + +import pandas as pd +import pyomo.environ as pyo +import pytest + +import idaes.logger as idaeslog +from idaes.apps.grid_integration import DesignModel, OperationModel +from idaes.apps.grid_integration.pricetaker.price_taker_model import PriceTakerModel +from idaes.apps.grid_integration.pricetaker.unit_commitment import UnitCommitmentData + +# pylint: disable = unused-import +from idaes.apps.grid_integration.pricetaker.tests.test_clustering import ( + dummy_data_fixture, + sklearn_avail, +) +from idaes.core.util.exceptions import ConfigurationError + +pytest.importorskip("sklearn", reason="sklearn not available") + + +def simple_flowsheet_func(m): + """Dummy flowsheet function for testing""" + m.x = pyo.Var() + m.y = pyo.Var() + m.con1 = pyo.Constraint(expr=m.x + m.y == 1) + m.LMP = pyo.Param(initialize=0.0, mutable=True) + + m.blk = pyo.Block() + m.blk.power = pyo.Var() + m.blk.op_mode = pyo.Var(within=pyo.Binary) + m.blk.LMP = pyo.Param(initialize=0.0, mutable=True) + + +def foo_design_model(m, max_power, min_power): + """Dummy design model""" + m.PMAX = pyo.Param(initialize=max_power, mutable=False) + m.PMIN = pyo.Param(initialize=min_power, mutable=False) + + m.capacity = pyo.Var() + + # DesignModel automatically declares install_unit variables + m.low_capacity_limit = pyo.Constraint(expr=m.PMIN * m.install_unit <= m.capacity) + m.up_capacity_limit = pyo.Constraint(expr=m.capacity <= m.PMAX * m.install_unit) + + # Capital investment cost ($) + m.capex = 1000.0 + # Fixed operating and investment ($ / year) + m.fom = 30.0 + + +def foo_operation_model(m, design_blk): + """Dummy operation model""" + # Operation Variables + m.power = pyo.Var(domain=pyo.NonNegativeReals, bounds=(0, design_blk.PMAX)) + + # NOTE: OperationModel automatically declares op_mode, + # startup and shutdown binary variables, and LMP Parameter + # Surrogate cost constraint (combined fuel_cost and non_fuel_vom) + m.fuel_cost = pyo.Expression(expr=23 * m.power + 50 * m.op_mode) + m.elec_revenue = pyo.Expression(expr=m.power * m.LMP) + m.su_sd_costs = pyo.Expression(expr=5 * m.startup + 4 * m.shutdown) + + +def build_foo_flowsheet(m, design_blk): + """Dummy flowsheet model for testing""" + m.op_blk = OperationModel( + model_func=foo_operation_model, + model_args={"design_blk": design_blk}, + ) + + m.fixed_rev = pyo.Expression(expr=100) + m.fixed_cost = pyo.Expression(expr=50) + + +@pytest.mark.unit +def test_num_representative_days(): + """Tests the num_representative_days property""" + m = PriceTakerModel() + + # By default, num_representative_days should be None + assert m.num_representative_days is None + + # Assign number of representative days + m.num_representative_days = 20 + assert m.num_representative_days == 20 + + # Test overwrite error + with pytest.raises( + ConfigurationError, + match=( + "num_representative_days is already defined as 20 " + "and it cannot be overwritten." + "\n\tInstantiate a new PriceTakerModel object." + ), + ): + m.num_representative_days = 25 + + # Ensure that the value remained the same + assert m.num_representative_days == 20 + + # Test error if an invalid value is specified + m = PriceTakerModel() + with pytest.raises(ValueError): + m.num_representative_days = -0.5 + + +@pytest.mark.unit +def test_horizon_length(caplog): + """Tests the horizon_length property""" + m = PriceTakerModel() + + # Test default value warning is not printed, it it is assigned + with caplog.at_level(idaeslog.WARNING): + m.horizon_length = 48 + assert m.horizon_length == 48 + assert ( + "Attribute horizon_length is not specified. Using 24 hours " + "as the horizon length." + ) not in caplog.text + + # Test overwrite error + with pytest.raises( + ConfigurationError, + match=( + "horizon_length is already defined as 48 and it cannot be " + "overwritten.\n\tInstantiate a new PriceTakerModel object." + ), + ): + m.horizon_length = 72 + + # Ensure that the value remained the same + assert m.horizon_length == 48 + + # Test the default value warning + m = PriceTakerModel() + with caplog.at_level(idaeslog.WARNING): + assert m.horizon_length == 24 + assert ( + "Attribute horizon_length is not specified. Using 24 hours " + "as the horizon length." + ) in caplog.text + + # Test error message associated with an infeasible value + with pytest.raises(ValueError): + m.horizon_length = -4 + + +@pytest.mark.unit +def test_append_lmp_data(): + "Tests the append_lmp_data method" + m = PriceTakerModel() + data = [-4, 5.0, 2, 50, -3.2] + + m.append_lmp_data( + lmp_data=pd.Series(data), + num_representative_days=3, + horizon_length=2, + seed=30, + ) + + # pylint: disable = protected-access + assert m._config.lmp_data == data + assert m.num_representative_days == 3 + assert m.horizon_length == 2 + assert m._config.seed == 30 + + # Test the error associated with overwriting LMP data + with pytest.raises( + ConfigurationError, + match=( + "Attempted to overwrite the LMP data. Instantiate a " + "new PriceTakerModel object to change the LMP data." + ), + ): + m.append_lmp_data(lmp_data=data + [4, 7, 8]) + + +@pytest.mark.unit +@pytest.mark.skipif( + not sklearn_avail, reason="sklearn (optional dependency) not available" +) +def test_get_optimal_representative_days(dummy_data): + """Tests the get_optimal_representative_days method""" + m = PriceTakerModel() + + # Test lmp data not available error + with pytest.raises( + ConfigurationError, + match=( + "LMP data is missing. Please append the LMP data using the " + "`append_lmp_data` method." + ), + ): + m.get_optimal_representative_days() + + # Append LMP data + m.append_lmp_data(lmp_data=dummy_data, horizon_length=2) + + # Test invalid kmin and kmax error + with pytest.raises(ValueError): + m.get_optimal_representative_days(kmin=-5, kmax=-3.5) + + # Test if the method works + assert 3 == m.get_optimal_representative_days( + kmin=2, kmax=6, method="silhouette", generate_elbow_plot=False + ) + + +@pytest.mark.unit +def test_mp_model_full_year(dummy_data): + """ + Tests the build_multiperiod_model using the full year price signal + """ + m = PriceTakerModel() + + # Test lmp data not available error + with pytest.raises( + ConfigurationError, + match=( + "LMP data is missing. Please append the LMP data using the " + "`append_lmp_data` method." + ), + ): + m.build_multiperiod_model(flowsheet_func=simple_flowsheet_func) + + # Build the multiperiod model + m.append_lmp_data(lmp_data=dummy_data) + m.build_multiperiod_model(flowsheet_func=simple_flowsheet_func) + + # Test if the model construction is successful + assert m.horizon_length == len(dummy_data) + assert m.num_representative_days == 1 + assert len(m.period) == len(dummy_data) + assert isinstance(m.set_days, pyo.RangeSet) + assert isinstance(m.set_time, pyo.RangeSet) + assert len(m.set_days) == 1 + assert len(m.set_time) == len(dummy_data) + assert m.rep_days_weights == {1: 1} + assert m.rep_days_lmp == {1: {t + 1: val for t, val in enumerate(dummy_data)}} + for d, t in m.period: + assert isinstance(m.period[d, t].x, pyo.Var) + assert isinstance(m.period[d, t].y, pyo.Var) + assert isinstance(m.period[d, t].con1, pyo.Constraint) + + # Test if the LMP data is updated at flowsheet level + assert m.period[d, t].LMP.value == pytest.approx(dummy_data[t - 1]) + + # Test if the LMP data is updated for each operational block + assert m.period[d, t].blk.LMP.value == pytest.approx(dummy_data[t - 1]) + + # Test model overwrite error + with pytest.raises( + ConfigurationError, + match=( + "A multiperiod model might already exist, as the object has " + "`period` attribute." + ), + ): + m.build_multiperiod_model(flowsheet_func=simple_flowsheet_func) + + +@pytest.mark.unit +@pytest.mark.skipif( + not sklearn_avail, reason="sklearn (optional dependency) not available" +) +def test_mp_model_rep_days(dummy_data): + """ + Tests the build_multiperiod_model using representative days + """ + m = PriceTakerModel() + m.append_lmp_data(lmp_data=dummy_data, horizon_length=2, num_representative_days=3) + m.build_multiperiod_model(flowsheet_func=simple_flowsheet_func) + + # Test if the model construction is successful + assert m.horizon_length == 2 + assert m.num_representative_days == 3 + assert len(m.period) == 2 * 3 + assert len(m.set_days) == 3 + assert len(m.set_time) == 2 + assert m.rep_days_weights == {1: 4, 2: 4, 3: 4} + + rep_lmp_data = { + 1: {1: 10.5, 2: 0.5}, + 2: {1: 5.5, 2: 10.5}, + 3: {1: 0.5, 2: 0.5}, + } + assert m.rep_days_lmp == rep_lmp_data + for d, t in m.period: + assert isinstance(m.period[d, t].x, pyo.Var) + assert isinstance(m.period[d, t].y, pyo.Var) + assert isinstance(m.period[d, t].con1, pyo.Constraint) + + # Test if the LMP data is updated at flowsheet level + assert m.period[d, t].LMP.value == pytest.approx(rep_lmp_data[d][t]) + + # Test if the LMP data is updated for each operational block + assert m.period[d, t].blk.LMP.value == pytest.approx(rep_lmp_data[d][t]) + + +@pytest.mark.unit +def test_get_operation_blocks(dummy_data): + """Tests the _get_operation_blocks method""" + m = PriceTakerModel() + m.append_lmp_data(lmp_data=dummy_data) + + # Test multiperiod model does not exist error + with pytest.raises( + ConfigurationError, + match=( + "Unable to find the multiperiod model. Please use the " + "build_multiperiod_model method to construct one." + ), + ): + # pylint: disable = protected-access + m._get_operation_blocks(blk_name="blk", attribute_list=["op_mode"]) + + # Build multiperiod model + m.build_multiperiod_model(flowsheet_func=simple_flowsheet_func) + + # pylint: disable = protected-access + op_blocks = m._get_operation_blocks(blk_name="blk", attribute_list=["op_mode"]) + + assert len(op_blocks) == 1 + assert len(op_blocks[1]) == 24 + + for d, t in m.period: + assert op_blocks[d][t] is m.period[d, t].blk + + # Test operational block not found error + with pytest.raises( + AttributeError, + match="Operational block foo_blk does not exist.", + ): + m._get_operation_blocks(blk_name="foo_blk", attribute_list=["op_mode"]) + + # Test missing attribute error + with pytest.raises( + AttributeError, + match=( + "Required attribute startup is not found in " "the operational block blk." + ), + ): + m._get_operation_blocks(blk_name="blk", attribute_list=["op_mode", "startup"]) + + +@pytest.mark.unit +@pytest.mark.skipif( + not sklearn_avail, reason="sklearn (optional dependency) not available" +) +def test_get_operation_blocks_rep_days(dummy_data): + """Tests the get_operation_blocks method with representative days""" + m = PriceTakerModel() + m.append_lmp_data(lmp_data=dummy_data, horizon_length=2, num_representative_days=3) + m.build_multiperiod_model(flowsheet_func=simple_flowsheet_func) + + # pylint: disable = protected-access + op_blocks = m._get_operation_blocks(blk_name="blk", attribute_list=["op_mode"]) + + assert len(op_blocks) == 3 # 3 representative days + for d in m.rep_days_weights: + assert len(op_blocks[d]) == 2 # horizon length = 2 + + for d, t in m.period: + assert op_blocks[d][t] is m.period[d, t].blk + + +@pytest.mark.unit +def test_get_operation_vars(dummy_data): + """Tests the _get_operation_vars method""" + m = PriceTakerModel() + m.append_lmp_data(lmp_data=dummy_data) + + # Test multiperiod model does not exist error + with pytest.raises( + ConfigurationError, + match=( + "Unable to find the multiperiod model. Please use the " + "build_multiperiod_model method to construct one." + ), + ): + # pylint: disable = protected-access + m._get_operation_vars(var_name="blk.power") + + # Build multiperiod model + m.build_multiperiod_model(flowsheet_func=simple_flowsheet_func) + + # pylint: disable = protected-access + op_vars = m._get_operation_vars(var_name="blk.power") + + assert len(op_vars) == 1 + assert len(op_vars[1]) == 24 + + for d, t in m.period: + assert op_vars[d][t] is m.period[d, t].blk.power + + # Test operational block not found error + with pytest.raises( + AttributeError, + match="Variable foo_blk.foo does not exist in the multiperiod model.", + ): + m._get_operation_vars(var_name="foo_blk.foo") + + +@pytest.mark.unit +@pytest.mark.skipif( + not sklearn_avail, reason="sklearn (optional dependency) not available" +) +def test_get_operation_vars_rep_days(dummy_data): + """Tests the get_operation_vars method with representative days""" + m = PriceTakerModel() + m.append_lmp_data(lmp_data=dummy_data, horizon_length=2, num_representative_days=3) + m.build_multiperiod_model(flowsheet_func=simple_flowsheet_func) + + # pylint: disable = protected-access + op_vars = m._get_operation_vars(var_name="blk.power") + + assert len(op_vars) == 3 # 3 representative days + for d in m.rep_days_weights: + assert len(op_vars[d]) == 2 # horizon length = 2 + + for d, t in m.period: + assert op_vars[d][t] is m.period[d, t].blk.power + + +@pytest.mark.unit +def test_add_linking_constraints(dummy_data): + """Tests the add_linking_constraints method""" + m = PriceTakerModel() + m.append_lmp_data(lmp_data=dummy_data) + m.build_multiperiod_model(flowsheet_func=simple_flowsheet_func) + m.add_linking_constraints( + previous_time_var="blk.power", current_time_var="blk.op_mode" + ) + + assert len(m.variable_linking_constraints_1) == 24 - 1 + assert str(m.variable_linking_constraints_1[1, 10].expr) == ( + "period[1,9].blk.power == period[1,10].blk.op_mode" + ) + assert str(m.variable_linking_constraints_1[1, 24].expr) == ( + "period[1,23].blk.power == period[1,24].blk.op_mode" + ) + + # Add another set of linking constraints + m.add_linking_constraints(previous_time_var="x", current_time_var="y") + assert len(m.variable_linking_constraints_2) == 24 - 1 + assert str(m.variable_linking_constraints_2[1, 10].expr) == ( + "period[1,9].x == period[1,10].y" + ) + assert str(m.variable_linking_constraints_2[1, 24].expr) == ( + "period[1,23].x == period[1,24].y" + ) + + +@pytest.mark.unit +@pytest.mark.skipif( + not sklearn_avail, reason="sklearn (optional dependency) not available" +) +def test_add_linking_constraints_rep_days(dummy_data): + """Tests the add_linking_constraints method with representative days""" + m = PriceTakerModel() + m.append_lmp_data(lmp_data=dummy_data, horizon_length=12, num_representative_days=2) + m.build_multiperiod_model(flowsheet_func=simple_flowsheet_func) + m.add_linking_constraints( + previous_time_var="blk.power", current_time_var="blk.op_mode" + ) + + assert len(m.variable_linking_constraints_1) == 2 * (12 - 1) + assert str(m.variable_linking_constraints_1[1, 10].expr) == ( + "period[1,9].blk.power == period[1,10].blk.op_mode" + ) + assert str(m.variable_linking_constraints_1[2, 12].expr) == ( + "period[2,11].blk.power == period[2,12].blk.op_mode" + ) + + +@pytest.mark.unit +def test_retrieve_uc_data(): + """Tests the _retrieve_uc_data method""" + m = PriceTakerModel() + + # pylint: disable = protected-access + assert len(m._op_blk_uc_data) == 0 + uc_data = m._retrieve_uc_data(op_blk="blk_1", commodity="power", op_range_lb=0.2) + + assert len(m._op_blk_uc_data) == 1 + assert uc_data is m._op_blk_uc_data["blk_1", "power"] + assert isinstance(uc_data, UnitCommitmentData) + + # Test updating the data + uc_data_2 = m._retrieve_uc_data(op_blk="blk_1", commodity="power", startup_rate=0.2) + assert len(m._op_blk_uc_data) == 1 + # It should return the existing object. + assert uc_data_2 is uc_data + assert uc_data_2.config.startup_rate == pytest.approx(0.2) + + # Test invalid data error + with pytest.raises( + ConfigurationError, + match=( + "For commidity power in operational " + "block blk_1, \n\tthe shutdown rate is less than " + "the minimum stable operation value." + ), + ): + m._retrieve_uc_data(op_blk="blk_1", commodity="power", shutdown_rate=0.1) + + # Test addition of different commodity to the same operational block + uc_data_3 = m._retrieve_uc_data( + op_blk="blk_1", commodity="ng_flow", startup_rate=0.2 + ) + assert len(m._op_blk_uc_data) == 2 + assert ("blk_1", "ng_flow") in m._op_blk_uc_data + assert uc_data_3 is not uc_data + assert uc_data_3 is not uc_data_2 + + # Test addition of different operational block + uc_data_4 = m._retrieve_uc_data( + op_blk="blk_new", commodity="ng_flow", startup_rate=0.2 + ) + assert len(m._op_blk_uc_data) == 3 + assert ("blk_new", "ng_flow") in m._op_blk_uc_data + assert uc_data_4 is not uc_data + assert uc_data_4 is not uc_data_2 + assert uc_data_4 is not uc_data_3 + + assert list(m._op_blk_uc_data.keys()) == [ + ("blk_1", "power"), + ("blk_1", "ng_flow"), + ("blk_new", "ng_flow"), + ] + + +@pytest.mark.unit +def test_add_capacity_limits(dummy_data): + """Tests the add_capacity_limits method""" + + m = PriceTakerModel() + m.append_lmp_data(lmp_data=dummy_data) + method_args = { + "op_block_name": "blk", + "commodity": "power", + "capacity": 100, + "op_range_lb": 0.3, + } + + # Build the multiperiod model + m.build_multiperiod_model(flowsheet_func=simple_flowsheet_func) + + # Test operation block does not exist error + with pytest.raises( + AttributeError, + match="Operational block foo_blk does not exist.", + ): + method_args_2 = method_args.copy() + method_args_2["op_block_name"] = "foo_blk" + m.add_capacity_limits(**method_args_2) + + # Test if capacity limits are added correctly + m.add_capacity_limits(**method_args) + + assert len(m.blk_power_limits) == 1 + assert len(m.blk_power_limits[1].capacity_low_limit_con) == 24 + assert len(m.blk_power_limits[1].capacity_high_limit_con) == 24 + + # Test overwriting error + with pytest.raises( + ConfigurationError, + match="Attempting to overwrite capacity limits for power in blk.", + ): + m.add_capacity_limits(**method_args) + + +@pytest.mark.unit +@pytest.mark.skipif( + not sklearn_avail, reason="sklearn (optional dependency) not available" +) +def test_add_capacity_limits_rep_days(dummy_data): + """Tests the add_capacity_limits method with representative days""" + m = PriceTakerModel() + m.append_lmp_data(lmp_data=dummy_data, horizon_length=2, num_representative_days=3) + m.build_multiperiod_model(flowsheet_func=simple_flowsheet_func) + m.add_capacity_limits( + op_block_name="blk", commodity="power", capacity=100, op_range_lb=0.3 + ) + + # Test if capacity limits are added correctly + assert len(m.blk_power_limits) == 3 # 3 representative days + # Following must be valid because the horizon length is 2 + for d in [1, 2, 3]: + assert len(m.blk_power_limits[d].capacity_low_limit_con) == 2 + assert len(m.blk_power_limits[d].capacity_high_limit_con) == 2 + + +@pytest.mark.unit +def test_add_ramping_limits(dummy_data): + """Tests the add_ramping_limits method""" + + m = PriceTakerModel() + m.append_lmp_data(lmp_data=dummy_data) + m.design_blk = DesignModel( + model_func=foo_design_model, model_args={"max_power": 400, "min_power": 300} + ) + method_args = { + "op_block_name": "op_blk", + "commodity": "power", + "capacity": m.design_blk.capacity, + "startup_rate": 0.3, + "shutdown_rate": 0.4, + "rampup_rate": 0.2, + "rampdown_rate": 0.2, + } + + # Build the multiperiod model + m.build_multiperiod_model( + flowsheet_func=build_foo_flowsheet, + flowsheet_options={"design_blk": m.design_blk}, + ) + + # Test necessary attributes does not exist error + with pytest.raises( + AttributeError, + match="Operational block foo_blk does not exist.", + ): + method_args_2 = method_args.copy() + method_args_2["op_block_name"] = "foo_blk" + m.add_ramping_limits(**method_args_2) + + # Test if ramping limits are added correctly + m.add_ramping_limits(**method_args) + + # pylint: disable = protected-access + uc_data = m._op_blk_uc_data["op_blk", "power"].config + assert uc_data.startup_rate == pytest.approx(method_args["startup_rate"]) + assert uc_data.shutdown_rate == pytest.approx(method_args["shutdown_rate"]) + assert uc_data.rampup_rate == pytest.approx(method_args["rampup_rate"]) + assert uc_data.rampdown_rate == pytest.approx(method_args["rampdown_rate"]) + assert uc_data.capacity is m.design_blk.capacity + + assert len(m.op_blk_power_ramping) == 1 + assert len(m.op_blk_power_ramping[1].ramp_up_con) == 24 - 1 + assert len(m.op_blk_power_ramping[1].ramp_down_con) == 24 - 1 + + # Test overwriting error + with pytest.raises( + ConfigurationError, + match="Attempting to overwrite ramping limits for power in op_blk.", + ): + m.add_ramping_limits(**method_args) + + +@pytest.mark.unit +@pytest.mark.skipif( + not sklearn_avail, reason="sklearn (optional dependency) not available" +) +def test_add_ramping_limits_rep_days(dummy_data): + """Tests the add_ramping_limits method with representative days""" + m = PriceTakerModel() + m.append_lmp_data(lmp_data=dummy_data, horizon_length=12, num_representative_days=2) + m.design_blk = DesignModel( + model_func=foo_design_model, model_args={"max_power": 400, "min_power": 300} + ) + m.build_multiperiod_model( + flowsheet_func=build_foo_flowsheet, + flowsheet_options={"design_blk": m.design_blk}, + ) + m.add_ramping_limits( + op_block_name="op_blk", + commodity="power", + capacity=m.design_blk.capacity, + startup_rate=0.3, + shutdown_rate=0.4, + rampup_rate=0.2, + rampdown_rate=0.2, + ) + assert len(m.op_blk_power_ramping) == 2 # 2 representative days + # Following must be valid because the horizon length is 12 + for d in [1, 2]: + assert len(m.op_blk_power_ramping[d].ramp_up_con) == 12 - 1 + assert len(m.op_blk_power_ramping[d].ramp_down_con) == 12 - 1 + + +@pytest.mark.unit +def test_start_shut_no_install(dummy_data): + """Tests the add_startup_shutdown method without install_unit variable.""" + + m = PriceTakerModel() + m.append_lmp_data(lmp_data=dummy_data) + m.design_blk = DesignModel( + model_func=foo_design_model, model_args={"max_power": 400, "min_power": 300} + ) + # Build the multiperiod model + m.build_multiperiod_model( + flowsheet_func=build_foo_flowsheet, + flowsheet_options={"design_blk": m.design_blk}, + ) + + # Test if startup and shutdown constraints are added correctly + m.add_startup_shutdown(op_block_name="op_blk", up_time=3, down_time=4) + + assert hasattr(m, "op_blk_startup_shutdown") + assert len(m.op_blk_startup_shutdown) == 1 + assert len(m.op_blk_startup_shutdown[1].binary_relationship_con) == 24 - 1 + assert len(m.op_blk_startup_shutdown[1].minimum_up_time_con) == 24 - (3 - 1) + assert len(m.op_blk_startup_shutdown[1].minimum_down_time_con) == 24 - (4 - 1) + + con3_1_10 = ( + "period[1,7].op_blk.shutdown + period[1,8].op_blk.shutdown + " + "period[1,9].op_blk.shutdown + period[1,10].op_blk.shutdown <= " + "1 - period[1,10].op_blk.op_mode" + ) + assert con3_1_10 == str(m.op_blk_startup_shutdown[1].minimum_down_time_con[10].expr) + + # Test overwriting error + with pytest.raises( + ConfigurationError, + match=( + "Attempting to overwrite startup/shutdown constraints " + "for operation block op_blk." + ), + ): + m.add_startup_shutdown(op_block_name="op_blk", up_time=5, down_time=3) + + +@pytest.mark.unit +def test_start_shut_with_install(dummy_data): + """Tests the add_startup_shutdown method with install_unit variable.""" + + m = PriceTakerModel() + m.append_lmp_data(lmp_data=dummy_data) + m.design_blk = DesignModel( + model_func=foo_design_model, model_args={"max_power": 400, "min_power": 300} + ) + # Build the multiperiod model + m.build_multiperiod_model( + flowsheet_func=build_foo_flowsheet, + flowsheet_options={"design_blk": m.design_blk}, + ) + + # Test the error associated with the absence of install_unit variable + method_args = { + "op_block_name": "op_blk", + "des_block_name": "design_blk", + "up_time": 3, + "down_time": 4, + } + m.design_blk.del_component(m.design_blk.install_unit) + assert m.find_component("design_blk.install_unit") is None + + with pytest.raises( + AttributeError, + match=( + "Binary variable associated with unit installation is not found " + "in design_blk. \n\tDo not specify des_block_name argument if " + "installation of the unit is not a decision variable." + ), + ): + m.add_startup_shutdown(**method_args) + + # Construct the constraints with install_unit variable + m.design_blk.install_unit = pyo.Var(within=pyo.Binary) + assert m.find_component("design_blk.install_unit") is not None + m.add_startup_shutdown(**method_args) + + con3_1_10 = ( + "period[1,7].op_blk.shutdown + period[1,8].op_blk.shutdown + " + "period[1,9].op_blk.shutdown + period[1,10].op_blk.shutdown <= " + "design_blk.install_unit - period[1,10].op_blk.op_mode" + ) + assert con3_1_10 == str(m.op_blk_startup_shutdown[1].minimum_down_time_con[10].expr) + + +@pytest.mark.unit +@pytest.mark.skipif( + not sklearn_avail, reason="sklearn (optional dependency) not available" +) +def test_start_shut_with_rep_days(dummy_data): + """ + Tests the add_startup_shutdown method with install_unit variable, + and with representative days + """ + + m = PriceTakerModel() + m.append_lmp_data(lmp_data=dummy_data, horizon_length=12, num_representative_days=2) + m.design_blk = DesignModel( + model_func=foo_design_model, model_args={"max_power": 400, "min_power": 300} + ) + m.build_multiperiod_model( + flowsheet_func=build_foo_flowsheet, + flowsheet_options={"design_blk": m.design_blk}, + ) + m.add_startup_shutdown( + op_block_name="op_blk", + des_block_name="design_blk", + up_time=3, + down_time=4, + ) + + assert len(m.op_blk_startup_shutdown) == 2 # 2 representative days + # Following must be valid because the horizon length is 12 + for d in [1, 2]: + assert len(m.op_blk_startup_shutdown[d].binary_relationship_con) == 12 - 1 + assert len(m.op_blk_startup_shutdown[d].minimum_up_time_con) == 12 - (3 - 1) + assert len(m.op_blk_startup_shutdown[d].minimum_down_time_con) == 12 - (4 - 1) + + con3_1_10 = ( + "period[1,7].op_blk.shutdown + period[1,8].op_blk.shutdown + " + "period[1,9].op_blk.shutdown + period[1,10].op_blk.shutdown <= " + "design_blk.install_unit - period[1,10].op_blk.op_mode" + ) + con3_2_10 = ( + "period[2,7].op_blk.shutdown + period[2,8].op_blk.shutdown + " + "period[2,9].op_blk.shutdown + period[2,10].op_blk.shutdown <= " + "design_blk.install_unit - period[2,10].op_blk.op_mode" + ) + assert con3_1_10 == str(m.op_blk_startup_shutdown[1].minimum_down_time_con[10].expr) + assert con3_2_10 == str(m.op_blk_startup_shutdown[2].minimum_down_time_con[10].expr) + + +@pytest.mark.unit +def test_add_hourly_cashflows_warnings(dummy_data, caplog): + """Tests the add_hourly_cashflows method with empty args""" + + m = PriceTakerModel() + m.append_lmp_data(lmp_data=dummy_data) + m.design_blk = DesignModel( + model_func=foo_design_model, model_args={"max_power": 400, "min_power": 300} + ) + # Build the multiperiod model + m.build_multiperiod_model( + flowsheet_func=build_foo_flowsheet, + flowsheet_options={"design_blk": m.design_blk}, + ) + + with caplog.at_level(idaeslog.WARNING): + m.add_hourly_cashflows() + + assert ( + "Argument operational_costs is not specified, so the total " + "operational cost will be set to 0." + ) in caplog.text + + assert ( + "Argument revenue_streams is not specified, so the total " + "revenue will be set to 0." + ) in caplog.text + + for d, t in m.period: + assert str(m.period[d, t].total_hourly_cost.expr) == "0.0" + assert str(m.period[d, t].total_hourly_revenue.expr) == "0.0" + assert str(m.period[d, t].net_hourly_cash_inflow.expr) == "0 - 0" + + +@pytest.mark.unit +def test_add_hourly_cashflows(dummy_data, caplog): + """Tests the add_hourly_cashflows method""" + + m = PriceTakerModel() + m.append_lmp_data(lmp_data=dummy_data) + m.design_blk = DesignModel( + model_func=foo_design_model, model_args={"max_power": 400, "min_power": 300} + ) + # Build the multiperiod model + m.build_multiperiod_model( + flowsheet_func=build_foo_flowsheet, + flowsheet_options={"design_blk": m.design_blk}, + ) + + with caplog.at_level(idaeslog.WARNING): + m.add_hourly_cashflows( + operational_costs=["fuel_cost", "su_sd_costs", "fixed_cost"], + revenue_streams=["elec_revenue", "fixed_rev"], + ) + + assert ( + "Argument operational_costs is not specified, so the total " + "operational cost will be set to 0." + ) not in caplog.text + + assert ( + "Argument revenue_streams is not specified, so the total " + "revenue will be set to 0." + ) not in caplog.text + + for d, t in m.period: + assert str(m.period[d, t].total_hourly_cost.expr) == ( + f"23*period[{d},{t}].op_blk.power + 50*period[{d},{t}].op_blk.op_mode" + f" + (5*period[{d},{t}].op_blk.startup + 4*period[{d},{t}].op_blk.shutdown)" + f" + 50" + ) + assert str(m.period[d, t].total_hourly_revenue.expr) == ( + f"period[{d},{t}].op_blk.LMP*period[{d},{t}].op_blk.power + 100" + ) + + +@pytest.mark.unit +def test_add_overall_cashflows_warnings(dummy_data, caplog): + """Tests the warning/error messages in add_overall_cashflows""" + m = PriceTakerModel() + m.append_lmp_data(lmp_data=dummy_data) + m.build_multiperiod_model(flowsheet_func=simple_flowsheet_func) + + with pytest.raises( + ConfigurationError, + match=( + "Hourly cashflows are not added to the model. Please run " + "add_hourly_cashflows method before calling the " + "add_overall_cashflows method." + ), + ): + m.add_overall_cashflows() + + with caplog.at_level(idaeslog.WARNING): + m.add_hourly_cashflows() + m.add_overall_cashflows() + + assert ( + "Argument operational_costs is not specified, so the total " + "operational cost will be set to 0." + ) in caplog.text + + assert ( + "Argument revenue_streams is not specified, so the total " + "revenue will be set to 0." + ) in caplog.text + + assert ( + "No design blocks were found, so the overall capital cost " + "(capex) and the fixed O&M cost (fom) are set to 0." + ) in caplog.text + + cf = m.cashflows + assert str(cf.capex_calculation.expr) == "cashflows.capex == 0" + assert str(cf.fom_calculation.expr) == "cashflows.fom == 0" + assert str(cf.depreciation_calculation.expr) == ( + "cashflows.depreciation == 0.03333333333333333*cashflows.capex" + ) + assert hasattr(cf, "net_cash_inflow") + assert hasattr(cf, "net_cash_inflow_calculation") + assert str(cf.corporate_tax_calculation.expr) == ( + "0.2*(cashflows.net_cash_inflow - cashflows.fom - cashflows.depreciation)" + " <= cashflows.corporate_tax" + ) + + assert str(cf.net_profit_calculation.expr) == ( + "cashflows.net_profit == cashflows.net_cash_inflow - cashflows.fom" + " - cashflows.corporate_tax" + ) + assert str(cf.npv.expr) == ( + "cashflows.net_profit - 0.08882743338727227*cashflows.capex" + ) + assert str(cf.lifetime_npv.expr) == ( + "11.257783343127485*cashflows.net_profit - cashflows.capex" + ) + + +@pytest.mark.unit +def test_add_overall_cashflows(dummy_data): + """Tests the add_overall_cashflows method""" + + m = PriceTakerModel() + m.append_lmp_data(lmp_data=dummy_data) + m.design_blk = DesignModel( + model_func=foo_design_model, model_args={"max_power": 400, "min_power": 300} + ) + # Build the multiperiod model + m.build_multiperiod_model( + flowsheet_func=build_foo_flowsheet, + flowsheet_options={"design_blk": m.design_blk}, + ) + m.add_hourly_cashflows( + operational_costs=["fuel_cost", "su_sd_costs", "fixed_cost"], + revenue_streams=["elec_revenue", "fixed_rev"], + ) + m.add_overall_cashflows() + cf = m.cashflows + assert str(cf.capex_calculation.expr) == "cashflows.capex == 1000.0" + assert str(cf.fom_calculation.expr) == "cashflows.fom == 30.0" + + +@pytest.mark.unit +def test_add_objective_function(dummy_data): + """Tests the add_objective_function method""" + m = PriceTakerModel() + m.append_lmp_data(lmp_data=dummy_data) + m.build_multiperiod_model(flowsheet_func=simple_flowsheet_func) + + with pytest.raises( + ConfigurationError, + match=( + "Overall cashflows are not appended. Please run the " + "add_overall_cashflows method." + ), + ): + m.add_objective_function(objective_type="npv") + + # Test the error associated with unsupported objective function + m.add_hourly_cashflows() + m.add_overall_cashflows() + + with pytest.raises( + ConfigurationError, + match=( + "foo is not a supported objective function." + "Please specify either npv, or lifetime_npv, or net_profit " + "as the objective_type." + ), + ): + m.add_objective_function(objective_type="foo") + + # Test if the objective function is added + m.add_objective_function(objective_type="net_profit") + assert hasattr(m, "obj") diff --git a/idaes/apps/grid_integration/pricetaker/tests/test_unit_commitment.py b/idaes/apps/grid_integration/pricetaker/tests/test_unit_commitment.py new file mode 100644 index 0000000000..90a08a5af2 --- /dev/null +++ b/idaes/apps/grid_integration/pricetaker/tests/test_unit_commitment.py @@ -0,0 +1,308 @@ +################################################################################# +# The Institute for the Design of Advanced Energy Systems Integrated Platform +# Framework (IDAES IP) was produced under the DOE Institute for the +# Design of Advanced Energy Systems (IDAES). +# +# Copyright (c) 2018-2023 by the software owners: The Regents of the +# University of California, through Lawrence Berkeley National Laboratory, +# National Technology & Engineering Solutions of Sandia, LLC, Carnegie Mellon +# University, West Virginia University Research Corporation, et al. +# All rights reserved. Please see the files COPYRIGHT.md and LICENSE.md +# for full copyright and license information. +################################################################################# + +import pyomo.environ as pyo +import pytest + +import idaes.apps.grid_integration.pricetaker.unit_commitment as uc + + +@pytest.mark.unit +def test_unit_commitment_data(): + """Tests the UnitCommitmentData class""" + m = uc.UnitCommitmentData( + blk_name="ngcc", + commodity_name="power", + op_range_lb=0.3, + ) + assert m.config.startup_rate is None + assert m.config.shutdown_rate is None + assert m.config.rampup_rate is None + assert m.config.rampdown_rate is None + assert m.config.op_range_lb == pytest.approx(0.3) + assert m.blk_name == "ngcc" + assert m.commodity_name == "power" + + # Test the update method + mdl = pyo.ConcreteModel() + mdl.capacity = pyo.Var() + + m.update( + startup_rate=0.4, + rampup_rate=0.4, + capacity=mdl.capacity, + ) + assert m.config.startup_rate == pytest.approx(0.4) + assert m.config.rampup_rate == pytest.approx(0.4) + assert isinstance(m.config.capacity, pyo.Var) + assert m.config.capacity.value is None + + # Test assertion error associated with invalid startup_rate + with pytest.raises( + uc.ConfigurationError, + match=( + "For commidity power in operational " + "block ngcc, \n\tthe startup rate is less than " + "the minimum stable operation value." + ), + ): + m.update(startup_rate=0.2) + + # Test assertion error associated with invalid shutdown_rate + with pytest.raises( + uc.ConfigurationError, + match=( + "For commidity power in operational " + "block ngcc, \n\tthe shutdown rate is less than " + "the minimum stable operation value." + ), + ): + m.update(startup_rate=0.4, shutdown_rate=0.2) + + # Test assertion error associated with incomplete data + with pytest.raises( + uc.ConfigurationError, + match=( + "Necessary arguments needed for the ramping constraints " "are missing." + ), + ): + # rampdown_rate value is missing, so it should throw an error. + m.assert_ramping_args_present() + + +@pytest.mark.unit +def test_startup_shutdown_constraints(): + """Tests the startup_shutdown_constraints function""" + m = pyo.ConcreteModel() + m.set_time = pyo.RangeSet(10) + + @m.Block(m.set_time) + def op_blk(b, _): + b.op_mode = pyo.Var(within=pyo.Binary) + b.startup = pyo.Var(within=pyo.Binary) + b.shutdown = pyo.Var(within=pyo.Binary) + + m.startup_shutdown = pyo.Block() + + # Append startup and shutdown constraints + uc.startup_shutdown_constraints( + blk=m.startup_shutdown, + op_blocks=m.op_blk, + install_unit=1, + up_time=3, + down_time=4, + set_time=m.set_time, + ) + + ss = m.startup_shutdown + + # Tests both the existence of the constraints and + # the number of constraints + # pylint: disable=no-member + assert len(ss.binary_relationship_con) == 9 + assert len(ss.minimum_up_time_con) == 8 + assert len(ss.minimum_down_time_con) == 7 + + # Check constraint is added for correct indices + assert 1 not in ss.binary_relationship_con + assert 1 not in ss.minimum_up_time_con + assert 2 not in ss.minimum_up_time_con + assert 1 not in ss.minimum_down_time_con + assert 2 not in ss.minimum_down_time_con + assert 3 not in ss.minimum_down_time_con + + # Check if the constraints are implemented correctly + con1_2 = ( + "op_blk[2].op_mode - op_blk[1].op_mode == " + "op_blk[2].startup - op_blk[2].shutdown" + ) + con1_10 = ( + "op_blk[10].op_mode - op_blk[9].op_mode == " + "op_blk[10].startup - op_blk[10].shutdown" + ) + assert con1_2 == str(ss.binary_relationship_con[2].expr) + assert con1_10 == str(ss.binary_relationship_con[10].expr) + + con2_3 = ( + "op_blk[1].startup + op_blk[2].startup + " + "op_blk[3].startup <= op_blk[3].op_mode" + ) + con2_10 = ( + "op_blk[8].startup + op_blk[9].startup + " + "op_blk[10].startup <= op_blk[10].op_mode" + ) + assert con2_3 == str(ss.minimum_up_time_con[3].expr) + assert con2_10 == str(ss.minimum_up_time_con[10].expr) + + con3_4 = ( + "op_blk[1].shutdown + op_blk[2].shutdown + " + "op_blk[3].shutdown + op_blk[4].shutdown <= " + "1 - op_blk[4].op_mode" + ) + con3_10 = ( + "op_blk[7].shutdown + op_blk[8].shutdown + " + "op_blk[9].shutdown + op_blk[10].shutdown <= " + "1 - op_blk[10].op_mode" + ) + assert con3_4 == str(ss.minimum_down_time_con[4].expr) + assert con3_10 == str(ss.minimum_down_time_con[10].expr) + + +@pytest.mark.unit +def test_capacity_limits_with_float(): + """ + Tests the capacity_limits function with a constant + for capacity variable + """ + m = pyo.ConcreteModel() + m.set_time = pyo.RangeSet(10) + + @m.Block(m.set_time) + def op_blk(b, _): + b.power = pyo.Var(within=pyo.NonNegativeReals) + b.op_mode = pyo.Var(within=pyo.Binary) + + m.cap_limits = pyo.Block() + uc_data = uc.UnitCommitmentData( + blk_name="op_blk", commodity_name="power", op_range_lb=0.3, capacity=650.0 + ) + uc.capacity_limits( + blk=m.cap_limits, + op_blocks=m.op_blk, + uc_data=uc_data, + set_time=m.set_time, + ) + + cl = m.cap_limits + # Check the existence (and the number) of constraints + # pylint: disable=no-member + assert len(cl.capacity_low_limit_con) == 10 + assert len(cl.capacity_high_limit_con) == 10 + + # Check if the constraint is implemented correctly + con1_1 = "195.0*op_blk[1].op_mode <= op_blk[1].power" + con1_9 = "195.0*op_blk[9].op_mode <= op_blk[9].power" + assert con1_1 == str(cl.capacity_low_limit_con[1].expr) + assert con1_9 == str(cl.capacity_low_limit_con[9].expr) + + con2_1 = "op_blk[1].power <= 650.0*op_blk[1].op_mode" + con2_9 = "op_blk[9].power <= 650.0*op_blk[9].op_mode" + assert con2_1 == str(cl.capacity_high_limit_con[1].expr) + assert con2_9 == str(cl.capacity_high_limit_con[9].expr) + + +@pytest.mark.unit +def test_capacity_limits_with_var(): + """ + Tests the capacity_limits function with a Pyomo Var + for capacity variable + """ + m = pyo.ConcreteModel() + m.set_time = pyo.RangeSet(10) + m.capacity = pyo.Var(bounds=(150, 600)) + + @m.Block(m.set_time) + def op_blk(b, _): + b.power = pyo.Var(within=pyo.NonNegativeReals) + b.op_mode = pyo.Var(within=pyo.Binary) + + m.cap_limits = pyo.Block() + uc_data = uc.UnitCommitmentData( + blk_name="op_blk", commodity_name="power", op_range_lb=0.3, capacity=m.capacity + ) + uc.capacity_limits( + blk=m.cap_limits, + op_blocks=m.op_blk, + uc_data=uc_data, + set_time=m.set_time, + ) + + cl = m.cap_limits + # Check if the constraint is implemented correctly + # pylint: disable=no-member + con1_1 = "0.3*capacity*op_blk[1].op_mode <= op_blk[1].power" + con1_9 = "0.3*capacity*op_blk[9].op_mode <= op_blk[9].power" + assert con1_1 == str(cl.capacity_low_limit_con[1].expr) + assert con1_9 == str(cl.capacity_low_limit_con[9].expr) + + con2_1 = "op_blk[1].power <= capacity*op_blk[1].op_mode" + con2_9 = "op_blk[9].power <= capacity*op_blk[9].op_mode" + assert con2_1 == str(cl.capacity_high_limit_con[1].expr) + assert con2_9 == str(cl.capacity_high_limit_con[9].expr) + + +@pytest.mark.unit +def test_ramping_limits(): + """Tests the ramping_limits constraint""" + m = pyo.ConcreteModel() + m.set_time = pyo.RangeSet(10) + m.capacity = pyo.Var(bounds=(150, 600)) + + @m.Block(m.set_time) + def op_blk(b, _): + b.power = pyo.Var(within=pyo.NonNegativeReals) + b.op_mode = pyo.Var(within=pyo.Binary) + b.startup = pyo.Var(within=pyo.Binary) + b.shutdown = pyo.Var(within=pyo.Binary) + + m.ramp_limits = pyo.Block() + uc_data = uc.UnitCommitmentData( + blk_name="op_blk", + commodity_name="power", + capacity=m.capacity, + startup_rate=0.3, + shutdown_rate=0.3, + rampup_rate=0.2, + rampdown_rate=0.15, + ) + uc.ramping_limits( + blk=m.ramp_limits, + op_blocks=m.op_blk, + uc_data=uc_data, + set_time=m.set_time, + ) + + rl = m.ramp_limits + # Check the existence (and the number) of constraints + # pylint: disable=no-member + assert 1 not in rl.ramp_up_con + assert 1 not in rl.ramp_down_con + assert len(rl.ramp_up_con) == 9 + assert len(rl.ramp_down_con) == 9 + + # Check if the constraint is implemented correctly + con1_2 = ( + "op_blk[2].power - op_blk[1].power <= " + "0.3*capacity*op_blk[2].startup + " + "0.2*capacity*op_blk[1].op_mode" + ) + con1_10 = ( + "op_blk[10].power - op_blk[9].power <= " + "0.3*capacity*op_blk[10].startup + " + "0.2*capacity*op_blk[9].op_mode" + ) + assert con1_2 == str(rl.ramp_up_con[2].expr) + assert con1_10 == str(rl.ramp_up_con[10].expr) + + con2_2 = ( + "op_blk[1].power - op_blk[2].power <= " + "0.3*capacity*op_blk[2].shutdown + " + "0.15*capacity*op_blk[2].op_mode" + ) + con2_10 = ( + "op_blk[9].power - op_blk[10].power <= " + "0.3*capacity*op_blk[10].shutdown + " + "0.15*capacity*op_blk[10].op_mode" + ) + assert con2_2 == str(rl.ramp_down_con[2].expr) + assert con2_10 == str(rl.ramp_down_con[10].expr) diff --git a/idaes/apps/grid_integration/pricetaker/unit_commitment.py b/idaes/apps/grid_integration/pricetaker/unit_commitment.py new file mode 100644 index 0000000000..a63c257e76 --- /dev/null +++ b/idaes/apps/grid_integration/pricetaker/unit_commitment.py @@ -0,0 +1,241 @@ +################################################################################# +# The Institute for the Design of Advanced Energy Systems Integrated Platform +# Framework (IDAES IP) was produced under the DOE Institute for the +# Design of Advanced Energy Systems (IDAES). +# +# Copyright (c) 2018-2023 by the software owners: The Regents of the +# University of California, through Lawrence Berkeley National Laboratory, +# National Technology & Engineering Solutions of Sandia, LLC, Carnegie Mellon +# University, West Virginia University Research Corporation, et al. +# All rights reserved. Please see the files COPYRIGHT.md and LICENSE.md +# for full copyright and license information. +################################################################################# + +""" +This module contains functions that build unit commitment-type +constraints: startup/shutdown, uptime/downtime constraints, +capacity limit constraints, and ramping constraints. +""" + +from typing import Union +from pyomo.common.config import ConfigDict, ConfigValue +from pyomo.environ import Block, Constraint, Var, RangeSet +from idaes.core.util.config import ConfigurationError, is_in_range + + +class UnitCommitmentData: + """Dataclass to store startup, shutdown, and ramp rates""" + + def __init__(self, blk_name: str, commodity_name: str, **kwargs): + self.blk_name = blk_name + self.commodity_name = commodity_name + self.config = self._get_config() + self.update(**kwargs) + + def _get_config(self): + config = ConfigDict() + config.declare( + "startup_rate", + ConfigValue( + domain=is_in_range(0, 1), + doc="Startup rate as a fraction of design capacity", + ), + ) + config.declare( + "shutdown_rate", + ConfigValue( + domain=is_in_range(0, 1), + doc="Shutdown rate as a fraction of design capacity", + ), + ) + config.declare( + "rampup_rate", + ConfigValue( + domain=is_in_range(0, 1), + doc="Rampup rate as a fraction of design capacity", + ), + ) + config.declare( + "rampdown_rate", + ConfigValue( + domain=is_in_range(0, 1), + doc="Rampdown rate as a fraction of design capacity", + ), + ) + config.declare( + "op_range_lb", + ConfigValue( + domain=is_in_range(0, 1), + doc="Minimum stable operation range as a fraction of design capacity", + ), + ) + config.declare( + "capacity", + ConfigValue( + doc=( + "Parameter/variable denoting the maximum " + "capacity of the commodity" + ), + ), + ) + + return config + + def update(self, **kwargs): + """Updates the attribute values""" + self.config.set_value(kwargs) + self.assert_startup_rate_validity() + self.assert_shutdown_rate_validity() + + def assert_startup_rate_validity(self): + """Raises an error if startup rate value is not valid.""" + if None in (self.config.op_range_lb, self.config.startup_rate): + # One of the values is not provided, so skip the check + return + + if self.config.startup_rate < self.config.op_range_lb: + raise ConfigurationError( + f"For commidity {self.commodity_name} in operational " + f"block {self.blk_name}, \n\tthe startup rate is less than " + f"the minimum stable operation value." + ) + + def assert_shutdown_rate_validity(self): + """Raises an error if shutdown rate value is not valid""" + if None in (self.config.op_range_lb, self.config.shutdown_rate): + # One of the values is not provided, so skip the check + return + + if self.config.shutdown_rate < self.config.op_range_lb: + raise ConfigurationError( + f"For commidity {self.commodity_name} in operational " + f"block {self.blk_name}, \n\tthe shutdown rate is less than " + f"the minimum stable operation value." + ) + + def assert_ramping_args_present(self): + """ + Raises an error if any of the arguments needed for the ramping + constraints is missing. + """ + cf = self.config + if ( + None + in ( + cf.startup_rate, + cf.shutdown_rate, + cf.rampup_rate, + cf.rampdown_rate, + ) + or cf.capacity is None + ): + raise ConfigurationError( + "Necessary arguments needed for the ramping constraints are missing." + ) + + +def startup_shutdown_constraints( + blk: Block, + op_blocks: dict, + install_unit: Union[int, Var], + up_time: int, + down_time: int, + set_time: RangeSet, +): + """ + Appends startup and shutdown constraints for a given unit/process + """ + + @blk.Constraint(set_time) + def binary_relationship_con(_, t): + if t == 1: + return Constraint.Skip + + return ( + op_blocks[t].op_mode - op_blocks[t - 1].op_mode + == op_blocks[t].startup - op_blocks[t].shutdown + ) + + @blk.Constraint(set_time) + def minimum_up_time_con(_, t): + if t < up_time: + return Constraint.Skip + + return ( + sum(op_blocks[i].startup for i in range(t - up_time + 1, t + 1)) + <= op_blocks[t].op_mode + ) + + @blk.Constraint(set_time) + def minimum_down_time_con(_, t): + if t < down_time: + return Constraint.Skip + + return ( + sum(op_blocks[i].shutdown for i in range(t - down_time + 1, t + 1)) + <= install_unit - op_blocks[t].op_mode + ) + + +def capacity_limits( + blk: Block, + op_blocks: dict, + uc_data: UnitCommitmentData, + set_time: RangeSet, +): + """ + Appends capacity limit constraints + """ + commodity = uc_data.commodity_name + limits = ( + uc_data.config.op_range_lb * uc_data.config.capacity, + uc_data.config.capacity, + ) + commodity = {t: getattr(blk, commodity) for t, blk in op_blocks.items()} + + @blk.Constraint(set_time) + def capacity_low_limit_con(_, t): + return limits[0] * op_blocks[t].op_mode <= commodity[t] + + @blk.Constraint(set_time) + def capacity_high_limit_con(_, t): + return commodity[t] <= limits[1] * op_blocks[t].op_mode + + +def ramping_limits( + blk: Block, + op_blocks: dict, + uc_data: UnitCommitmentData, + set_time: RangeSet, +): + """ + Appends ramping constraints + """ + _ramping_var = uc_data.commodity_name + startup_rate = uc_data.config.startup_rate * uc_data.config.capacity + shutdown_rate = uc_data.config.shutdown_rate * uc_data.config.capacity + rampup_rate = uc_data.config.rampup_rate * uc_data.config.capacity + rampdown_rate = uc_data.config.rampdown_rate * uc_data.config.capacity + ramping_var = {t: getattr(blk, _ramping_var) for t, blk in op_blocks.items()} + + @blk.Constraint(set_time) + def ramp_up_con(_, t): + if t == 1: + return Constraint.Skip + + return ( + ramping_var[t] - ramping_var[t - 1] + <= startup_rate * op_blocks[t].startup + + rampup_rate * op_blocks[t - 1].op_mode + ) + + @blk.Constraint(set_time) + def ramp_down_con(_, t): + if t == 1: + return Constraint.Skip + + return ( + ramping_var[t - 1] - ramping_var[t] + <= shutdown_rate * op_blocks[t].shutdown + + rampdown_rate * op_blocks[t].op_mode + ) diff --git a/idaes/core/util/config.py b/idaes/core/util/config.py index 3f77e58e10..e9beaf2198 100644 --- a/idaes/core/util/config.py +++ b/idaes/core/util/config.py @@ -192,3 +192,22 @@ def DefaultBool(arg): return arg else: return Bool(arg) + + +def is_in_range(lb, ub): + """ + Domain validator for 1D compact sets. + + Args: + lb: float, lower bound + ub: float, upper bound + """ + + def _in_range(val): + if lb <= val <= ub: + return val + raise ConfigurationError( + f"Value {val} lies outside the admissible range {[lb, ub]}" + ) + + return _in_range diff --git a/idaes/core/util/tests/test_config.py b/idaes/core/util/tests/test_config.py index b7b67446a3..7c7330330d 100644 --- a/idaes/core/util/tests/test_config.py +++ b/idaes/core/util/tests/test_config.py @@ -36,6 +36,7 @@ is_transformation_method, is_transformation_scheme, DefaultBool, + is_in_range, ) from idaes.core.util.exceptions import ConfigurationError @@ -219,3 +220,20 @@ def test_DefaultBool(): assert not DefaultBool(False) with pytest.raises(ValueError): DefaultBool("foo") + + +@pytest.mark.unit +def test_is_in_range(): + assert is_in_range(0, 1)(0) == 0 + assert is_in_range(0, 1)(1) == 1 + assert is_in_range(0, 1)(0.5) == 0.5 + with pytest.raises( + ConfigurationError, + match="Value -1 lies outside the admissible range \\[0, 1\\]", + ): + is_in_range(0, 1)(-1) + + with pytest.raises( + ConfigurationError, match="Value 2 lies outside the admissible range \\[0, 1\\]" + ): + is_in_range(0, 1)(2)