"""Provides various utility functions that don't fit within a common theme."""
from __future__ import annotations
import math
from string import digits, punctuation
from typing import TYPE_CHECKING, Callable
from functools import cache
import numpy as np
import pandas as pd
# Don't repeat the most common inputs that occur when there are no state changes, but
# the result is needed again for logging
[docs]
@cache
def _mean(*args) -> float:
"""Multiplies two numbers. Used for a reduce operation.
Parameters
----------
args : int | float
The values to compute the mean over
Returns
-------
float
The average of the values provided
"""
return np.mean(args)
[docs]
def create_variable_from_string(string: str) -> str:
"""Creates a valid Python variable style string from a passed string.
Parameters
----------
string : str
The string to convert into a Python-friendly variable name.
Returns
-------
str
A Python-valid variable name.
Examples
--------
>>> example_string = "*Electrical!*_ _System$*_"
>>> print(create_variable_from_string(example_string))
'electrical_system'
"""
new_string = (
string.lstrip(punctuation + digits)
.translate(str.maketrans("", "", punctuation.replace("_", " ")))
.lower()
)
return new_string
[docs]
def IEC_power_curve(
windspeed_column: np.ndarray | pd.Series,
power_column: np.ndarray | pd.Series,
bin_width: float = 0.5,
windspeed_start: float = 0.0,
windspeed_end: float = 30.0,
) -> Callable:
"""Direct copyfrom OpenOA:
https://github.com/NREL/OpenOA/blob/main/operational_analysis/toolkits/power_curve/functions.py#L16-L57
Use IEC 61400-12-1-2 method for creating wind-speed binned power curve.
Parameters
----------
windspeed_column : np.ndarray | pandas.Series
The power curve's windspeed values, in m/s.
power_column : np.ndarray | pandas.Series
The power curve's output power values, in kW.
bin_width : float
Width of windspeed bin, default is 0.5 m/s according to standard, by default
0.5.
windspeed_start : float
Left edge of first windspeed bin, where all proceeding values will be 0.0,
by default 0.0.
windspeed_end : float
Right edge of last windspeed bin, where all following values will be 0.0, by
default 30.0.
Returns
-------
Callable
Python function of the power curve, of type (Array[float] -> Array[float]),
that maps input windspeed value(s) to ouptut power value(s).
"""
if not isinstance(windspeed_column, pd.Series):
windspeed_column = pd.Series(windspeed_column)
if not isinstance(power_column, pd.Series):
power_column = pd.Series(power_column)
# Set up evenly spaced bins of fixed width, with np.inf for any value over the max
n_bins = int(np.ceil((windspeed_end - windspeed_start) / bin_width)) + 1
bins = np.append(np.linspace(windspeed_start, windspeed_end, n_bins), [np.inf])
# Initialize an array which will hold the mean values of each bin
P_bin = np.ones(len(bins) - 1) * np.nan
# Compute the mean of each bin and set corresponding P_bin
for ibin in range(0, len(bins) - 1):
indices = (windspeed_column >= bins[ibin]) & (windspeed_column < bins[ibin + 1])
P_bin[ibin] = power_column.loc[indices].mean()
# Linearly interpolate any missing bins
P_bin = pd.Series(data=P_bin).interpolate(method="linear").bfill().values
# Create a closure over the computed bins which computes the power curve value for
# arbitrary array-like input
def pc_iec(x):
P = np.zeros(np.shape(x))
for i in range(0, len(bins) - 1):
idx = np.where((x >= bins[i]) & (x < bins[i + 1]))
P[idx] = P_bin[i]
cutoff_idx = (x < windspeed_start) | (x > windspeed_end)
P[cutoff_idx] = 0.0
return P
return pc_iec
[docs]
def calculate_stack_current(
power: np.ndarray | pd.Series,
p1: int | float,
p2: int | float,
p3: int | float,
p4: int | float,
p5: int | float,
) -> np.ndarray | pd.Series:
"""Convert power produced, in kW to current, in amperes (I) using the efficiency
curve from the
`H2Integrage PEM electrolysis model <https://github.com/NREL/H2Integrate/blob/main/h2integrate/simulation/technologies/hydrogen/electrolysis/PEM_H2_LT_electrolyzer_Clusters.py>`_.
.. math::
i_{stack} = p1 * power^{3} + p2 * power^{2} + p3 * power + p4 * power^{1/2} + p5
Parameters
----------
power : np.ndarray | pd.Series
Power produced, in kW.
p1 : int | float
First coefficient.
p2 : int | float
Second coefficient.
p3 : int | float
Third coefficient
p4 : int | float
Fourth coefficient.
p5 : int | float
Intercept.
Returns
-------
np.ndarray | pd.Series
Converted current from the farm's power production.
"""
return p1 * power**3 + p2 * power**2 + p3 * power + p4 * np.sqrt(power) + p5
[docs]
def calculate_hydrogen_production(
rated_capacity: float,
FE: float = 0.9999999,
n_cells: int = 135,
turndown_ratio: float = 0.1,
*,
efficiency_rate: float | None = None,
p1: int | float | None = None,
p2: int | float | None = None,
p3: int | float | None = None,
p4: int | float | None = None,
p5: int | float | None = None,
) -> Callable:
"""Create the hydrogen production curve for an electrolyzer. One of
:py:attr:`efficiency_rate` or a complete set of polynomial values (:py:attr:`p1`,
:py:attr:`p2`, :py:attr:`p3`, :py:attr:`p4`, and :py:attr:`p5`) for the power
to current conversion must be provided.
kWh are converted to current using the following efficiency formulation when the
polynomial inputs are all provided, and the :py:attr:`efficiency_rate` is nonzero.
.. math::
i_{stack} = p1 * power^{3} + p2 * power^{2} + p3 * power + p4 * power^{1/2} + p5
Parameters
----------
rated_capacity : float
Rated maximum input power of the electrolyzer, in kW.
FE : float, optional
Faradic efficiency, by default 0.9999999.
n_cells : int, optional
Number of cells per 1 MW stack, by default 135.
turndown_ratio : float, optional
Minimum input power as a ratio of the rated capacity, by default 0.1.
efficiency_rate : float, optional
Energy efficiency in kWh per kg H2. Required if ``p1`` through ``p5 are not
provided.
p1 : int | float | None
First coefficient of the efficiency polynomial curve. Required if
:py:attr:`efficiency_rate` is not provided, by default None.
p2 : int | float | None
Second coefficient of the efficiency polynomial curve. Required if
:py:attr:`efficiency_rate` is not provided, by default None.
p3 : int | float | None
Third coefficient of the efficiency polynomial curve. Required if
:py:attr:`efficiency_rate` is not provided, by default None.
p4 : int | float | None
Fourth coefficient of the efficiency polynomial curve. Required if
:py:attr:`efficiency_rate` is not provided, by default None.
p5 : int | float | None
Intercept of the efficiency polynomial curve. Required if
:py:attr:`efficiency_rate` is not provided, by default None.
Returns
-------
Callable
Function to calculate H₂ production (kg/hr) from power input (kW).
Raises
------
ValueError
Raised when :py:attr:`turndown_ratio` is outside the range [0, 1].
ValueError
Raised when :py:attr:`efficiency_rate` is provided and less than or equal to 0.
ValueError
Raised when neither the :py:attr:`efficiency_rate` nor polynomial values (
:py:attr:`p1` through :py:attr:`p5`) are provided.
ValueError
Raised both the :py:attr:`efficiency_rate` and polynomial values (
:py:attr:`p1` through :py:attr:`p5`) are provided.
"""
if turndown_ratio > 1 or turndown_ratio < 0:
msg = f"`turndown_ratio` must be between 0 and 1, not: {turndown_ratio}"
raise ValueError(msg)
use_efficiency = False
if efficiency_rate is not None:
if efficiency_rate > 0:
use_efficiency = True
else:
raise ValueError("`efficiency_rate` must be a positive number.")
poly_names = ("p1", "p2", "p3", "p4", "p5")
poly_vals = (p1, p2, p3, p4, p5)
invalid_polynomials = [n for n, v in zip(poly_names, poly_vals) if v is None]
use_polynomial = len(invalid_polynomials) == 0
if not use_efficiency and not use_polynomial:
msg = (
"One of `efficiency_rate` or all polynomial values (`p1`, `p2`, `p3`,"
f" `p4`, `p5`) must be provided. Received {efficiency_rate=},"
f" {p1=}, {p2=}, {p3=}, {p4=}, {p5=}."
)
raise ValueError(msg)
if use_efficiency and use_polynomial:
msg = (
"Only one of `efficiency_rate` or all polynomial values (`p1`, `p2`, `p3`,"
f" `p4`, `p5`) must be provided. Received {efficiency_rate=},"
f" {p1=}, {p2=}, {p3=}, {p4=}, {p5=}."
)
raise ValueError(msg)
F = 96485.34 # Faraday's Constant (C/mol)
SECONDS = 3600 # seconds per hour
molar_mass_h2 = 2.016 # molecular weight of H2 (grams/mol)
min_power = rated_capacity * turndown_ratio
def h2(power):
power = np.where(power > rated_capacity, rated_capacity, power)
h2_rate_kg_hr = np.zeros_like(power)
if use_efficiency:
h2_rate_kg_hr = power / efficiency_rate
else:
if TYPE_CHECKING:
assert isinstance(p1, int | float)
assert isinstance(p2, int | float)
assert isinstance(p3, int | float)
assert isinstance(p4, int | float)
assert isinstance(p5, int | float)
i = calculate_stack_current(power, p1, p2, p3, p4, p5)
h2_rate_g_s = i * n_cells * FE * molar_mass_h2 / (2 * F)
h2_rate_kg_hr = h2_rate_g_s * SECONDS / 1e3
h2_rate_kg_hr[h2_rate_kg_hr < 0] = 0
h2_rate_kg_hr[power < min_power] = 0
return h2_rate_kg_hr
return h2
[docs]
def calculate_windfarm_operational_level(
operations: pd.DataFrame,
turbine_id: np.ndarray | list[str],
turbine_weights: pd.DataFrame,
substation_turbine_map: dict[str, dict[str, np.ndarray]],
) -> pd.DataFrame:
"""Calculates the overall wind farm operational level, accounting for substation
downtime by multiplying the sum of all downstream turbine operational levels by
the substation's operational level.
Parameters
----------
operations : pd.DataFrame
The turbine and substation operational level DataFrame.
turbine_id : np.ndarray | list[str]
The turbine ids that match :py:attr:`Windfarm.turbine_id`.
turbine_weights : pd.DataFrame
The turbine weights, coming from :py:attr:`Windfarm.turbine_weights`.
substation_turbine_map : dict[str, dict[str, np.ndarray]]
The :py:attr:`Windfarm.substation_turbine_map`.
Notes
-----
This is a crude cap on the operations, and so a smarter way of capping
the availability should be added in the future.
Returns
-------
pd.DataFrame
The aggregate wind farm operational level.
"""
turbines = turbine_weights[turbine_id].values * operations[turbine_id]
total = np.sum(
[
np.array(
[[math.fsum(row)] for _, row in turbines[val["turbines"]].iterrows()]
).reshape(-1, 1)
for sub, val in substation_turbine_map.items()
],
axis=0,
)
return total
