"""Wind turbine performance models."""
from __future__ import annotations
import math
from abc import ABC, abstractmethod
from dataclasses import KW_ONLY, dataclass
from functools import cache, cached_property
from pathlib import Path
from typing import TYPE_CHECKING, Any, overload
import numpy as np
from ._repr import _FieldRepr
from ._validation import (
_validate_positive,
_validate_real,
_validate_real_array,
_validate_unit_interval,
)
if TYPE_CHECKING:
from collections.abc import Callable
from numpy.typing import ArrayLike
__all__ = [
"Controlled",
"Idealized",
"Tabulated",
"Turbine",
]
_BETZ_LIMIT = 16.0 / 27.0
# Controlled: off-design power loss exp(-decay |λ/λ_opt - 1| ** exponent).
_CP_LOSS_DECAY = 0.8
_CP_LOSS_EXPONENT = 3.0
# Buhl (2005, https://doi.org/10.2172/15016819) replaces momentum thrust above
# this induction.
_BUHL_INDUCTION = 0.4
_BUHL_THRUST = 0.96
# Wind speed about which thrust_slope tilts the design thrust, in m/s.
_THRUST_SLOPE_SPEED = 8.0
# Larger intensities inflate the apparent C_P well above cp_max.
_MAX_TURBULENCE_SMOOTHING = 0.1
# Exported curves add knots at 0 m/s, just outside cut-in and cut-out and at
# twice cut-out, so the zero a table reads outside them keeps idle and parked
# thrust. The gap stays representable after rounding.
_IDLE_EXPORT_SPEED = 0.0
_PARKED_EXPORT_MARGIN = 2.0
_STATE_GAP = 1.0e-6
# Sea-level standard atmosphere, assumed for windIO Cp curves, which carry no
# air density.
_AIR_DENSITY = 1.225
# Speeds evaluated per block, so the curve temporaries stay in cache.
_BLOCK_SIZE = 32768
# Rated-speed scan points evaluated per block, from the idealized rated speed
# up to the first crossing.
_RATED_SCAN_BLOCK = 256
def _validate_cp(name: str, value: float) -> float:
value = _validate_positive(name, value)
if value > _BETZ_LIMIT:
raise ValueError(f"{name} must not exceed the Betz limit")
return value
def _validate_thrust(name: str, value: float) -> float:
value = _validate_positive(name, value)
if value > 1.0:
raise ValueError(f"{name} must not exceed 1")
return value
def _as_speed_array(wind_speed) -> tuple[np.ndarray, bool]:
values = _validate_real_array("wind_speed", wind_speed)
scalar = values.ndim == 0
values = np.atleast_1d(values)
if not np.all(np.isfinite(values)):
raise ValueError("wind_speed must contain only finite values")
if np.any(values < 0.0):
raise ValueError("wind_speed must be non-negative")
return values, scalar
def _maybe_scalar(values: np.ndarray, scalar: bool) -> float | np.ndarray:
return float(values[0]) if scalar else values
def _curve_speeds(cutin: float, rated: float, cutout: float, n_points: int):
"""Export speeds, ``n_points`` from cut-in to cut-out plus the state knots.
``rated`` adds a knot when it is reached, that is when it is not ``nan``.
"""
if n_points < 2:
raise ValueError("n_points must be at least 2")
# The idle knot sits just below cut-in, and halfway to zero for a cut-in
# closer to zero than the gap.
idle = max(cutin - _STATE_GAP, 0.5 * cutin)
speeds = np.concatenate(
[
np.asarray([_IDLE_EXPORT_SPEED, idle], dtype=float),
np.linspace(cutin, cutout, int(n_points)),
np.asarray(
[cutout + _STATE_GAP, cutout * _PARKED_EXPORT_MARGIN], dtype=float
),
np.asarray([] if np.isnan(rated) else [rated], dtype=float),
]
)
return np.unique(np.round(speeds, 12))
def _axial_induction(cp_values):
"""Invert ``C_P = 4 a (1 - a)^2`` on its lower branch, ``a <= 1/3``."""
cp_values = np.asarray(cp_values, dtype=float)
normalized = np.clip(cp_values / _BETZ_LIMIT, 0.0, 1.0)
argument = np.clip(1.0 - 2.0 * normalized, -1.0, 1.0)
return (2.0 / 3.0) * (1.0 - np.cos(np.arccos(argument) / 3.0))
def _ct_from_cp(cp_values, cp_rated: float, ct_rated: float):
"""Thrust scaled from ``ct_rated`` by ``a (1 - a)`` of the induction."""
induction = _axial_induction(cp_values)
induction_rated = float(_axial_induction(cp_rated))
denominator = induction_rated * (1.0 - induction_rated)
return ct_rated * induction * (1.0 - induction) / denominator
def _momentum_cp(induction):
return 4.0 * induction * (1.0 - induction) ** 2
def _momentum_ct(induction):
"""Momentum thrust of an actuator disc, with Buhl's correction above a = 0.4."""
momentum = 4.0 * induction * (1.0 - induction)
buhl = 8.0 / 9.0 - (4.0 / 9.0) * induction + (14.0 / 9.0) * induction**2
return np.where(induction > _BUHL_INDUCTION, buhl, momentum)
def _induction_from_ct(ct):
"""Invert :func:`_momentum_ct` on its lower branch."""
ct = np.asarray(ct, dtype=float)
momentum = 0.5 * (1.0 - np.sqrt(np.clip(1.0 - ct, 0.0, None)))
buhl = (4.0 + np.sqrt(np.clip(16.0 - 56.0 * (8.0 - 9.0 * ct), 0.0, None))) / 28.0
return np.where(ct > _BUHL_THRUST, buhl, momentum)
def _scalar_momentum_cp(induction: float) -> float:
""":func:`_momentum_cp` of one float, rounded as the array version."""
complement = 1.0 - induction
return 4.0 * induction * (complement * complement)
def _scalar_induction_from_ct(ct: float) -> float:
""":func:`_induction_from_ct` of one float."""
if ct > _BUHL_THRUST:
return (4.0 + math.sqrt(max(16.0 - 56.0 * (8.0 - 9.0 * ct), 0.0))) / 28.0
return 0.5 * (1.0 - math.sqrt(max(1.0 - ct, 0.0)))
def _blockwise(function, values: np.ndarray, size: int) -> np.ndarray:
"""Apply the elementwise *function* to *values* in blocks of *size*."""
if len(values) <= size:
return function(values)
result = np.empty_like(values)
for start in range(0, len(values), size):
result[start : start + size] = function(values[start : start + size])
return result
@cache
def _smoothing_rule() -> tuple[np.ndarray, np.ndarray]:
"""Nine-point Gauss-Hermite nodes and weights for a standard normal."""
nodes, weights = np.polynomial.hermite_e.hermegauss(9)
return nodes, weights / weights.sum()
[docs]
class Turbine(ABC):
"""Wind turbine with power and thrust-coefficient curves.
The base of every OOWL turbine model. Wind speeds are in m/s, power in
watts and heights and diameters in meters.
"""
name: str | None
hub_height: float
rotor_diameter: float
rated_power: float
cutin_wind_speed: float
cutout_wind_speed: float
@property
@abstractmethod
def rated_wind_speed(self) -> float:
"""Wind speed at which rated power is reached, ``nan`` if never."""
@overload
def power(self, wind_speed: float) -> float: ...
@overload
def power(self, wind_speed: ArrayLike) -> np.ndarray: ...
[docs]
@abstractmethod
def power(self, wind_speed: ArrayLike) -> float | np.ndarray:
"""Power in watts at one or more wind speeds."""
@overload
def ct(self, wind_speed: float) -> float: ...
@overload
def ct(self, wind_speed: ArrayLike) -> np.ndarray: ...
[docs]
@abstractmethod
def ct(self, wind_speed: ArrayLike) -> float | np.ndarray:
"""Thrust coefficient at one or more wind speeds."""
[docs]
def save(self, path: str | Path, **kwargs: Any) -> Path:
"""Write this turbine to *path* as a windIO ``turbine`` document.
Wraps :func:`oowl.save` and returns the written path. ``name`` defaults
to the turbine's name.
"""
from ._windio import save
return save(self, path, **kwargs)
def _to_windio_turbine(self, name: str | None = None, n_points: int = 50) -> dict:
"""Return the windIO turbine mapping, with curves sampled at ``n_points``."""
speeds = _curve_speeds(
self.cutin_wind_speed,
self.rated_wind_speed,
self.cutout_wind_speed,
n_points,
)
return self._windio_mapping(name, speeds, self.power(speeds), self.ct(speeds))
def _windio_mapping(self, name, speeds, power_values, ct_values) -> dict:
speeds = np.asarray(speeds, dtype=float)
return {
"name": name or self.name or type(self).__name__,
"performance": {
"power_curve": {
"power_values": np.asarray(power_values, dtype=float).tolist(),
"power_wind_speeds": speeds.tolist(),
},
"Ct_curve": {
"Ct_values": np.asarray(ct_values, dtype=float).tolist(),
"Ct_wind_speeds": speeds.tolist(),
},
},
"hub_height": float(self.hub_height),
"rotor_diameter": float(self.rotor_diameter),
}
class _AnalyticalTurbine(_FieldRepr, Turbine):
"""Validation and helpers shared by the analytical models.
Concrete models supply ``cp``, ``ct``, ``power``, ``rated_wind_speed`` and
the operating thrust at cut-out, the parked default.
"""
rated_power: float
cp_max: float
air_density: float
generator_efficiency: float
ct_idle: float | None
ct_parked: float | None
_terminal_operating_ct: float
@overload
def cp(self, wind_speed: float) -> float: ...
@overload
def cp(self, wind_speed: ArrayLike) -> np.ndarray: ...
@abstractmethod
def cp(self, wind_speed: ArrayLike) -> float | np.ndarray:
"""Power coefficient at one or more wind speeds."""
def _validate_rotor(self) -> None:
self.rated_power = _validate_positive("rated_power", self.rated_power)
self.rotor_diameter = _validate_positive("rotor_diameter", self.rotor_diameter)
self.hub_height = _validate_positive("hub_height", self.hub_height)
if self.hub_height < self.rotor_diameter / 2.0:
raise ValueError("hub_height must keep the rotor above ground")
self.cp_max = _validate_cp("cp_max", self.cp_max)
def _validate_operation(self) -> None:
if self.ct_idle is not None:
self.ct_idle = _validate_unit_interval("ct_idle", self.ct_idle)
if self.ct_parked is not None:
self.ct_parked = _validate_unit_interval("ct_parked", self.ct_parked)
self.air_density = _validate_positive("air_density", self.air_density)
self.generator_efficiency = _validate_positive(
"generator_efficiency", self.generator_efficiency
)
if self.generator_efficiency > 1.0:
raise ValueError("generator_efficiency must not exceed 1")
self.cutin_wind_speed = _validate_positive(
"cutin_wind_speed", self.cutin_wind_speed
)
self.cutout_wind_speed = _validate_positive(
"cutout_wind_speed", self.cutout_wind_speed
)
if self.cutout_wind_speed <= self.cutin_wind_speed:
raise ValueError("cutout_wind_speed must exceed cutin_wind_speed")
@property
def rotor_area(self) -> float:
"""Rotor swept area in square meters."""
return float(np.pi * (self.rotor_diameter / 2.0) ** 2)
def _operating(self, speeds: np.ndarray) -> np.ndarray:
return (speeds >= self.cutin_wind_speed) & (speeds <= self.cutout_wind_speed)
def _power_scale(self) -> float:
"""Watts per unit ``C_P`` and cubic meter per cubic second of wind."""
return self.generator_efficiency * 0.5 * self.air_density * self.rotor_area
def _idle_and_parked(self, values: np.ndarray, speeds: np.ndarray) -> None:
"""Set ``ct_idle`` below cut-in and ``ct_parked`` above cut-out, in place.
``values`` arrives zero outside the operating band, the idle default.
"""
if self.ct_idle is not None:
values[speeds < self.cutin_wind_speed] = self.ct_idle
above_cutout = speeds > self.cutout_wind_speed
if np.any(above_cutout):
values[above_cutout] = (
self._terminal_operating_ct
if self.ct_parked is None
else self.ct_parked
)
[docs]
@dataclass(repr=False)
class Idealized(_AnalyticalTurbine):
"""Idealized turbine with a constant below-rated ``C_P``.
Follows the generalized turbine relations of Sørensen (2023,
https://doi.org/10.5194/wes-8-1017-2023), including the idealized cut-in
discontinuity. The rated wind speed is derived, not supplied.
Parameters
----------
rated_power
Electrical rated power in watts.
rotor_diameter, hub_height
Rotor size and height in meters.
cp_max
Power coefficient below rated, at most the Betz limit 16/27.
ct_rated
Thrust coefficient below rated.
air_density
Air density in kg/m³.
generator_efficiency
Share of the rotor power delivered as electrical power.
cutin_wind_speed, cutout_wind_speed
Operating band in m/s.
ct_idle, ct_parked
Thrust below cut-in, default 0, and above cut-out, default the
operating value at cut-out.
name
Label used in plots and windIO documents.
"""
rated_power: float
rotor_diameter: float
hub_height: float
cp_max: float
_: KW_ONLY
ct_rated: float = 0.8
air_density: float = 1.225
generator_efficiency: float = 1.0
cutin_wind_speed: float = 4.0
cutout_wind_speed: float = 25.0
ct_idle: float | None = None
ct_parked: float | None = None
name: str | None = None
def __post_init__(self):
self._validate_rotor()
self.ct_rated = _validate_unit_interval("ct_rated", self.ct_rated)
self._validate_operation()
if not self.cutin_wind_speed < self.rated_wind_speed < self.cutout_wind_speed:
raise ValueError(
"derived rated_wind_speed must lie between cut-in and cut-out"
)
@overload
def ct(self, wind_speed: float) -> float: ...
@overload
def ct(self, wind_speed: ArrayLike) -> np.ndarray: ...
[docs]
def ct(self, wind_speed: ArrayLike) -> float | np.ndarray:
"""Return the thrust coefficient in [0, 1] at ``wind_speed`` in m/s.
``ct_rated`` below rated. Above it, the thrust follows the momentum
induction of the falling ``C_P``. Below cut-in it is ``ct_idle``
(default 0), above cut-out ``ct_parked``, which defaults to the
operating value at cut-out.
"""
speeds, scalar = _as_speed_array(wind_speed)
values = np.zeros_like(speeds)
operating = self._operating(speeds)
values[operating & (speeds <= self.rated_wind_speed)] = self.ct_rated
above_rated = operating & (speeds > self.rated_wind_speed)
if np.any(above_rated):
values[above_rated] = _ct_from_cp(
self.cp(speeds[above_rated]), self.cp_max, self.ct_rated
)
self._idle_and_parked(values, speeds)
return _maybe_scalar(values, scalar)
@cached_property
def _terminal_operating_ct(self) -> float:
"""Operating C_T at cut-out, the parked fallback."""
cp_cutout = float(self.cp(self.cutout_wind_speed))
return float(_ct_from_cp(cp_cutout, self.cp_max, self.ct_rated))
@property
def rated_wind_speed(self) -> float:
"""Wind speed in m/s at which rated power is reached.
Derived from ``rated_power``, ``rotor_diameter``, ``cp_max``,
``generator_efficiency`` and ``air_density``.
"""
cube = self.rated_power / (self._power_scale() * self.cp_max)
return float(cube ** (1.0 / 3.0))
@overload
def cp(self, wind_speed: float) -> float: ...
@overload
def cp(self, wind_speed: ArrayLike) -> np.ndarray: ...
[docs]
def cp(self, wind_speed: ArrayLike) -> float | np.ndarray:
"""Return the power coefficient at ``wind_speed`` in m/s.
Scalar in, scalar out. Zero outside the cut-in to cut-out band.
"""
speeds, scalar = _as_speed_array(wind_speed)
values = np.zeros_like(speeds)
operating = self._operating(speeds)
values[operating & (speeds <= self.rated_wind_speed)] = self.cp_max
above_rated = operating & (speeds > self.rated_wind_speed)
values[above_rated] = (
self.cp_max * (self.rated_wind_speed / speeds[above_rated]) ** 3
)
return _maybe_scalar(values, scalar)
@overload
def power(self, wind_speed: float) -> float: ...
@overload
def power(self, wind_speed: ArrayLike) -> np.ndarray: ...
[docs]
def power(self, wind_speed: ArrayLike) -> float | np.ndarray:
"""Return the electrical power in watts at ``wind_speed`` in m/s."""
speeds, scalar = _as_speed_array(wind_speed)
values = self._power_scale() * np.asarray(self.cp(speeds)) * speeds**3
return _maybe_scalar(values, scalar)
[docs]
@dataclass(repr=False)
class Controlled(_AnalyticalTurbine):
"""Variable-speed turbine whose thrust and power follow one axial induction.
The rotor tracks ``optimal_tip_speed_ratio``, held between
``minimum_rotor_speed`` and ``maximum_tip_speed``. The induction follows
the realized tip-speed ratio as ``a = a_opt (λ / λ_opt) ** induction_exponent``,
where ``a_opt`` gives ``ct_optimal``. Thrust is the momentum value
``4 a (1 - a)``, with Buhl's correction above ``a = 0.4``. Power follows the
same induction, and its power coefficient never exceeds ``cp_max`` unless
``turbulence_smoothing`` averages it. Above rated the blades pitch
to hold rated power, and the induction falls with the power coefficient.
At the minimum rotor speed, near cut-in, the rotor over-induces and the
thrust rises above ``ct_optimal``. At the tip-speed limit it falls before
rated. ``induction_exponent=0`` gives a flat thrust below rated.
Rotor speed, tip-speed ratio and induction are diagnostics only and are
not exported to windIO.
Parameters
----------
rated_power
Electrical rated power in watts.
rotor_diameter, hub_height
Rotor size and height in meters.
cp_max
Power coefficient at the optimal tip-speed ratio.
optimal_tip_speed_ratio
Tip-speed ratio the controller tracks.
minimum_rotor_speed
Lowest rotor speed in rpm.
maximum_tip_speed
Highest blade tip speed in m/s.
ct_optimal
Thrust coefficient at the optimal tip-speed ratio.
induction_exponent
How strongly the induction follows the tip-speed ratio.
ct_peak_shaving
Thrust cap, or ``None``. The lower induction also lowers the power.
thrust_slope
Change of the design thrust per m/s about 8 m/s below rated, as in
fine-pitch schedules. Above rated the blades pitch instead.
turbulence_smoothing
Turbulence intensity, at most 0.1, over which power and thrust are
averaged. It rounds the knee at rated as in measured curves and keeps
cut-in and cut-out sharp. As in measured curves, the apparent ``C_P``
then exceeds ``cp_max``, by about 3% at 0.1 and up to about 13% just
above cut-in, where the average sees no idling speeds.
ct_idle, ct_parked
Thrust below cut-in, default 0, and above cut-out, default the
operating value at cut-out.
air_density
Air density in kg/m³.
generator_efficiency
Share of the rotor power delivered as electrical power.
cutin_wind_speed, cutout_wind_speed
Operating band in m/s.
name
Label used in plots and windIO documents.
"""
rated_power: float
rotor_diameter: float
hub_height: float
cp_max: float
optimal_tip_speed_ratio: float
minimum_rotor_speed: float
maximum_tip_speed: float
_: KW_ONLY
ct_optimal: float = 0.8
induction_exponent: float = 0.35
ct_peak_shaving: float | None = None
thrust_slope: float = 0.0
turbulence_smoothing: float = 0.0
air_density: float = 1.225
generator_efficiency: float = 0.96
cutin_wind_speed: float = 3.0
cutout_wind_speed: float = 25.0
ct_idle: float | None = None
ct_parked: float | None = None
name: str | None = None
def __post_init__(self):
self._validate_rotor()
self.optimal_tip_speed_ratio = _validate_positive(
"optimal_tip_speed_ratio", self.optimal_tip_speed_ratio
)
self.minimum_rotor_speed = _validate_positive(
"minimum_rotor_speed", self.minimum_rotor_speed
)
self.maximum_tip_speed = _validate_positive(
"maximum_tip_speed", self.maximum_tip_speed
)
if self._omega_min >= self._omega_max:
raise ValueError("minimum_rotor_speed must stay below maximum_tip_speed")
self.ct_optimal = _validate_thrust("ct_optimal", self.ct_optimal)
if self.ct_peak_shaving is not None:
self.ct_peak_shaving = _validate_thrust(
"ct_peak_shaving", self.ct_peak_shaving
)
self.induction_exponent = _validate_real(
"induction_exponent", self.induction_exponent
)
if not np.isfinite(self.induction_exponent) or self.induction_exponent < 0.0:
raise ValueError("induction_exponent must be non-negative and finite")
self.thrust_slope = _validate_real("thrust_slope", self.thrust_slope)
if not np.isfinite(self.thrust_slope):
raise ValueError("thrust_slope must be finite")
self.turbulence_smoothing = _validate_real(
"turbulence_smoothing", self.turbulence_smoothing
)
if not 0.0 <= self.turbulence_smoothing <= _MAX_TURBULENCE_SMOOTHING:
raise ValueError(
f"turbulence_smoothing must be in [0, {_MAX_TURBULENCE_SMOOTHING}]"
)
self._validate_operation()
band = np.asarray([self.cutin_wind_speed, self.cutout_wind_speed])
design_ct = self._design_ct(band)
if np.any(design_ct <= 0.0) or np.any(design_ct > 1.0):
raise ValueError(
"thrust_slope must keep the design thrust in (0, 1] from cut-in "
"to cut-out"
)
if self._blade_efficiency > 1.0:
raise ValueError("cp_max is unreachable at ct_optimal")
if not self.cutin_wind_speed < self.rated_wind_speed < self.cutout_wind_speed:
raise ValueError(
"controlled turbine must reach rated power after cut-in and before "
"cut-out"
)
# Above rated the pitched induction is at most 1/3 and the thrust at
# most 8/9, so only speeds up to rated can push the thrust above 1.
scan = np.linspace(self.cutin_wind_speed, self.cutout_wind_speed, 1025)
tracked = scan[scan <= self.rated_wind_speed]
if np.max(_momentum_ct(self._tracked_induction(tracked, tracked)[2])) > 1.0:
raise ValueError(
"thrust coefficient exceeds 1, lower induction_exponent or raise "
"minimum_rotor_speed"
)
@property
def rotor_radius(self) -> float:
"""Rotor radius in meters."""
return self.rotor_diameter / 2.0
@property
def _omega_min(self) -> float:
return self.minimum_rotor_speed * 2.0 * np.pi / 60.0
@property
def _omega_max(self) -> float:
return self.maximum_tip_speed / self.rotor_radius
@overload
def rotor_speed(self, wind_speed: float) -> float: ...
@overload
def rotor_speed(self, wind_speed: ArrayLike) -> np.ndarray: ...
[docs]
def rotor_speed(self, wind_speed: ArrayLike) -> float | np.ndarray:
"""Return the rotor speed in rpm, zero outside the operating band."""
speeds, scalar = _as_speed_array(wind_speed)
return _maybe_scalar(self._omega(speeds) * 60.0 / (2.0 * np.pi), scalar)
@overload
def tip_speed_ratio(self, wind_speed: float) -> float: ...
@overload
def tip_speed_ratio(self, wind_speed: ArrayLike) -> np.ndarray: ...
[docs]
def tip_speed_ratio(self, wind_speed: ArrayLike) -> float | np.ndarray:
"""Return the realized tip-speed ratio, zero outside the operating band."""
speeds, scalar = _as_speed_array(wind_speed)
return _maybe_scalar(self._tip_speed_ratio(speeds), scalar)
@overload
def axial_induction(self, wind_speed: float) -> float: ...
@overload
def axial_induction(self, wind_speed: ArrayLike) -> np.ndarray: ...
[docs]
def axial_induction(self, wind_speed: ArrayLike) -> float | np.ndarray:
"""Return the rotor-averaged axial induction, before turbulence smoothing."""
speeds, scalar = _as_speed_array(wind_speed)
return _maybe_scalar(self._induction(speeds), scalar)
def _band_omega(self, speeds: np.ndarray) -> np.ndarray:
"""Rotor speed in rad/s, for speeds inside the operating band."""
return np.clip(
self.optimal_tip_speed_ratio * speeds / self.rotor_radius,
self._omega_min,
self._omega_max,
)
def _omega(self, speeds: np.ndarray) -> np.ndarray:
return np.where(self._operating(speeds), self._band_omega(speeds), 0.0)
def _tip_speed_ratio(self, speeds: np.ndarray) -> np.ndarray:
return np.divide(
self._omega(speeds) * self.rotor_radius,
speeds,
out=np.zeros_like(speeds),
where=speeds > 0.0,
)
def _design_ct(self, speeds: np.ndarray) -> np.ndarray:
return self.ct_optimal + self.thrust_slope * (speeds - _THRUST_SLOPE_SPEED)
@cached_property
def _optimal_momentum_cp(self) -> float:
"""Momentum power coefficient of the induction that gives ``ct_optimal``."""
return float(_momentum_cp(_induction_from_ct(np.asarray(self.ct_optimal))))
@cached_property
def _blade_efficiency(self) -> float:
"""Share of the momentum power converted at the design point."""
return float(self.cp_max / self._optimal_momentum_cp)
@cached_property
def _induction_cap(self) -> float:
if self.ct_peak_shaving is None:
return float("inf")
return float(_induction_from_ct(np.asarray(self.ct_peak_shaving)))
@cached_property
def _optimal_induction(self) -> np.ndarray:
"""Design induction without ``thrust_slope``, as a one-element array."""
return _induction_from_ct(np.asarray([self.ct_optimal]))
def _tracked_induction(
self, speeds: np.ndarray, design_speeds: np.ndarray | None = None
) -> tuple[np.ndarray, np.ndarray, np.ndarray]:
"""Tip-speed ratio over its optimum, design induction and induction.
For speeds inside the operating band, before pitch regulation. The
design thrust is read at ``design_speeds``, by default ``speeds`` held
at rated, where pitch regulation takes over from the ``thrust_slope``
schedule.
"""
ratio = (
self._band_omega(speeds)
* self.rotor_radius
/ speeds
/ self.optimal_tip_speed_ratio
)
if self.thrust_slope == 0.0:
design = self._optimal_induction
else:
if design_speeds is None:
design_speeds = np.minimum(speeds, self.rated_wind_speed)
design = _induction_from_ct(self._design_ct(design_speeds))
induction = np.minimum(
design * ratio**self.induction_exponent, self._induction_cap
)
return ratio, design, induction
def _tracking_cp(
self, speeds: np.ndarray, design_speeds: np.ndarray | None = None
) -> np.ndarray:
"""Power coefficient before pitch regulation, inside the operating band.
Induction beyond the design value adds thrust but no power, and a
``thrust_slope`` that raises the design thrust adds no power either, so
``C_P`` never exceeds ``cp_max``. ``design_speeds`` is as in
:meth:`_tracked_induction`.
"""
ratio, design, induction = self._tracked_induction(speeds, design_speeds)
momentum = np.minimum(
np.minimum(_momentum_cp(induction), _momentum_cp(design)),
self._optimal_momentum_cp,
)
loss = np.exp(-_CP_LOSS_DECAY * np.abs(ratio - 1.0) ** _CP_LOSS_EXPONENT)
return self._blade_efficiency * momentum * loss
def _scalar_tracking(self) -> Callable[[float], tuple[float, float]]:
"""Return the induction and ``C_P`` of :meth:`_tracking_cp` as a function.
The function takes one in-band speed and reads the design thrust there.
It repeats the array arithmetic in scalar form, so a bisection step
takes about a microsecond.
"""
tip_speed_ratio = self.optimal_tip_speed_ratio
radius = self.rotor_radius
omega_min, omega_max = self._omega_min, self._omega_max
ct_optimal, thrust_slope = self.ct_optimal, self.thrust_slope
exponent, cap = self.induction_exponent, self._induction_cap
optimal_momentum = self._optimal_momentum_cp
efficiency = self._blade_efficiency
def tracking(speed: float) -> tuple[float, float]:
omega = min(max(tip_speed_ratio * speed / radius, omega_min), omega_max)
ratio = omega * radius / speed / tip_speed_ratio
design = _scalar_induction_from_ct(
ct_optimal + thrust_slope * (speed - _THRUST_SLOPE_SPEED)
)
# NumPy takes the square root for an exponent of 0.5, which pow
# can round differently.
induction = min(design * float(np.power(ratio, exponent)), cap)
momentum = min(
_scalar_momentum_cp(induction),
_scalar_momentum_cp(design),
optimal_momentum,
)
loss = math.exp(-_CP_LOSS_DECAY * abs(ratio - 1.0) ** _CP_LOSS_EXPONENT)
return induction, efficiency * momentum * loss
return tracking
@cached_property
def rated_wind_speed(self) -> float:
"""Wind speed in m/s at which rated power is first reached.
``nan`` when rated power is never reached, which the constructor
rejects.
"""
# The first crossing is found on the unregulated curve, so the design
# thrust is read at each scanned speed rather than held at rated.
scan = np.linspace(self.cutin_wind_speed, self.cutout_wind_speed, 4097)
scale = self._power_scale()
# C_P never exceeds cp_max, so no speed below the idealized rated speed
# reaches rated power. The margin covers rounding.
floor = (self.rated_power / (scale * self.cp_max)) ** (1.0 / 3.0)
start = int(np.searchsorted(scan, floor * (1.0 - 1e-9)))
for block_start in range(start, len(scan), _RATED_SCAN_BLOCK):
block = scan[block_start : block_start + _RATED_SCAN_BLOCK]
power = scale * self._tracking_cp(block, block) * block**3
crossings = np.flatnonzero(power >= self.rated_power)
if len(crossings):
index = block_start + int(crossings[0])
break
else:
return float("nan")
if index == 0:
return float(scan[0])
lower, upper = float(scan[index - 1]), float(scan[index])
tracking = self._scalar_tracking()
for _ in range(60):
midpoint = 0.5 * (lower + upper)
# Adjacent floats, so further steps change nothing.
if midpoint in (lower, upper):
break
if scale * tracking(midpoint)[1] * midpoint**3 >= self.rated_power:
upper = midpoint
else:
lower = midpoint
return 0.5 * (lower + upper)
@cached_property
def _pitch_efficiency(self) -> float:
"""Share of the momentum power converted at rated, held while pitching."""
induction, cp = self._scalar_tracking()(self.rated_wind_speed)
return cp / _momentum_cp(induction)
def _induction(self, speeds: np.ndarray) -> np.ndarray:
induction = np.zeros_like(speeds)
operating = self._operating(speeds)
induction[operating] = _blockwise(
self._band_induction, speeds[operating], _BLOCK_SIZE
)
return induction
def _band_induction(self, speeds: np.ndarray) -> np.ndarray:
"""Induction inside the operating band, pitched above rated."""
above = speeds > self.rated_wind_speed
below = ~above
induction = np.empty_like(speeds)
# Up to rated the design thrust is read at the speed itself.
tracked = speeds[below]
induction[below] = self._tracked_induction(tracked, tracked)[2]
pitched = speeds[above]
cp = self.rated_power / (self._power_scale() * pitched**3)
induction[above] = np.minimum(
_axial_induction(cp / self._pitch_efficiency), self._induction_cap
)
return induction
def _sharp_power(self, speeds: np.ndarray) -> np.ndarray:
power = self._power_scale() * self._tracking_cp(speeds) * speeds**3
return np.minimum(power, self.rated_power)
def _sharp_ct(self, speeds: np.ndarray) -> np.ndarray:
return _momentum_ct(self._band_induction(speeds))
def _smoothed(self, curve, speeds: np.ndarray) -> np.ndarray:
"""Average *curve* over turbulence, for speeds inside the operating band.
Speeds pass through *curve* in blocks, which keeps temporaries small.
"""
if self.turbulence_smoothing == 0.0:
return _blockwise(curve, speeds, _BLOCK_SIZE)
nodes, weights = _smoothing_rule()
factors = 1.0 + self.turbulence_smoothing * nodes
def average(block: np.ndarray) -> np.ndarray:
shifted = np.clip(
block[:, None] * factors,
self.cutin_wind_speed,
self.cutout_wind_speed,
)
values = curve(shifted.ravel()).reshape(shifted.shape)
# Summing node by node gives each speed the same value in any
# batch, which a BLAS matrix product does not.
total = values[:, 0] * weights[0]
for node in range(1, len(weights)):
total += values[:, node] * weights[node]
return total
return _blockwise(average, speeds, _BLOCK_SIZE // len(factors))
def _operating_curve(self, curve, wind_speed, fill: float = 0.0):
speeds, scalar = _as_speed_array(wind_speed)
values = np.full_like(speeds, fill)
operating = self._operating(speeds)
if np.any(operating):
values[operating] = self._smoothed(curve, speeds[operating])
return values, speeds, scalar
@overload
def power(self, wind_speed: float) -> float: ...
@overload
def power(self, wind_speed: ArrayLike) -> np.ndarray: ...
[docs]
def power(self, wind_speed: ArrayLike) -> float | np.ndarray:
"""Return the electrical power in watts at ``wind_speed`` in m/s."""
values, _, scalar = self._operating_curve(self._sharp_power, wind_speed)
return _maybe_scalar(values, scalar)
@overload
def cp(self, wind_speed: float) -> float: ...
@overload
def cp(self, wind_speed: ArrayLike) -> np.ndarray: ...
[docs]
def cp(self, wind_speed: ArrayLike) -> float | np.ndarray:
"""Return the power coefficient implied by :meth:`power`.
Scalar in, scalar out. Zero outside the cut-in to cut-out band.
"""
values, speeds, scalar = self._operating_curve(self._sharp_power, wind_speed)
operating = self._operating(speeds)
values[operating] /= self._power_scale() * speeds[operating] ** 3
return _maybe_scalar(values, scalar)
@overload
def ct(self, wind_speed: float) -> float: ...
@overload
def ct(self, wind_speed: ArrayLike) -> np.ndarray: ...
[docs]
def ct(self, wind_speed: ArrayLike) -> float | np.ndarray:
"""Return the thrust coefficient at ``wind_speed`` in m/s.
Below cut-in it is ``ct_idle`` (default 0). Above cut-out it is
``ct_parked``, which defaults to the operating value at cut-out.
"""
values, speeds, scalar = self._operating_curve(self._sharp_ct, wind_speed)
self._idle_and_parked(values, speeds)
return _maybe_scalar(values, scalar)
@cached_property
def _terminal_operating_ct(self) -> float:
"""Operating C_T at cut-out, the parked fallback."""
cutout = np.asarray([self.cutout_wind_speed])
return float(self._smoothed(self._sharp_ct, cutout)[0])
[docs]
@dataclass(init=False, repr=False)
class Tabulated(_FieldRepr, Turbine):
"""Turbine defined by tabulated power and thrust-coefficient curves.
Curves are linearly interpolated inside the tabulated range and are zero
outside it, so the first and last speeds act as cut-in and cut-out.
Parameters
----------
rotor_diameter, hub_height
Rotor size and height in meters.
wind_speeds
Strictly increasing, non-negative wind speeds in m/s.
power_values
Electrical power in watts at each speed.
ct_values
Thrust coefficient in [0, 1] at each speed.
name
Label used in plots and windIO documents.
rated_power
Rated power in watts. Defaults to the largest tabulated power.
"""
rotor_diameter: float
hub_height: float
wind_speeds: np.ndarray
power_values: np.ndarray
ct_values: np.ndarray
name: str | None = None
rated_power: float | None = None
def __init__(
self,
rotor_diameter: float,
hub_height: float,
wind_speeds: ArrayLike,
power_values: ArrayLike,
ct_values: ArrayLike,
*,
name: str | None = None,
rated_power: float | None = None,
) -> None:
self.rotor_diameter = _validate_positive("rotor_diameter", rotor_diameter)
self.hub_height = _validate_positive("hub_height", hub_height)
if self.hub_height < self.rotor_diameter / 2.0:
raise ValueError("hub_height must keep the rotor above ground")
# Copies, so later changes to the caller's arrays cannot bypass validation.
self.wind_speeds = np.array(wind_speeds, dtype=float)
self.power_values = np.array(power_values, dtype=float)
self.ct_values = np.array(ct_values, dtype=float)
if (
self.wind_speeds.ndim != 1
or self.power_values.shape != self.wind_speeds.shape
or self.ct_values.shape != self.wind_speeds.shape
or len(self.wind_speeds) < 2
):
raise ValueError("tabulated curves must be equal-length 1-D arrays")
if not (
np.all(np.isfinite(self.wind_speeds))
and np.all(np.isfinite(self.power_values))
and np.all(np.isfinite(self.ct_values))
):
raise ValueError("tabulated curves must contain only finite values")
if np.any(np.diff(self.wind_speeds) <= 0):
raise ValueError("wind_speeds must be strictly increasing")
if np.any(self.wind_speeds < 0) or np.any(self.power_values < 0):
raise ValueError("wind speeds and power values must be non-negative")
if np.any((self.ct_values < 0) | (self.ct_values > 1)):
raise ValueError("ct_values must be in [0, 1]")
self.name = name
if rated_power is None:
rated_power = float(np.max(self.power_values))
self.rated_power = _validate_positive("rated_power", rated_power)
def __eq__(self, other: object) -> bool:
"""Whether ``other`` is a tabulated turbine with the same curves and sizes."""
if other.__class__ is not self.__class__:
return NotImplemented
return (
self.rotor_diameter == other.rotor_diameter
and self.hub_height == other.hub_height
and self.name == other.name
and self.rated_power == other.rated_power
and np.array_equal(self.wind_speeds, other.wind_speeds)
and np.array_equal(self.power_values, other.power_values)
and np.array_equal(self.ct_values, other.ct_values)
)
@property
def cutin_wind_speed(self) -> float:
"""Lowest tabulated wind speed in m/s."""
return float(self.wind_speeds[0])
@property
def cutout_wind_speed(self) -> float:
"""Highest tabulated wind speed in m/s."""
return float(self.wind_speeds[-1])
@property
def rated_wind_speed(self) -> float:
"""Wind speed where the first upward segment reaches ``rated_power``.
Solved linearly on the tabulation. ``nan`` if never reached.
"""
rated_power = self.rated_power
assert rated_power is not None
for index in range(1, len(self.wind_speeds)):
lower_power = float(self.power_values[index - 1])
upper_power = float(self.power_values[index])
if lower_power < rated_power <= upper_power:
fraction = (rated_power - lower_power) / (upper_power - lower_power)
lower_speed = float(self.wind_speeds[index - 1])
upper_speed = float(self.wind_speeds[index])
return lower_speed + fraction * (upper_speed - lower_speed)
if self.power_values[0] >= rated_power:
return float(self.wind_speeds[0])
return float("nan")
def _interp(self, wind_speed, values):
speeds, scalar = _as_speed_array(wind_speed)
result = np.zeros_like(speeds)
operating = (speeds >= self.cutin_wind_speed) & (
speeds <= self.cutout_wind_speed
)
result[operating] = np.interp(speeds[operating], self.wind_speeds, values)
return _maybe_scalar(result, scalar)
@overload
def power(self, wind_speed: float) -> float: ...
@overload
def power(self, wind_speed: ArrayLike) -> np.ndarray: ...
[docs]
def power(self, wind_speed: ArrayLike) -> float | np.ndarray:
"""Return the electrical power in watts at ``wind_speed`` in m/s."""
return self._interp(wind_speed, self.power_values)
@overload
def ct(self, wind_speed: float) -> float: ...
@overload
def ct(self, wind_speed: ArrayLike) -> np.ndarray: ...
[docs]
def ct(self, wind_speed: ArrayLike) -> float | np.ndarray:
"""Return the thrust coefficient in [0, 1] at ``wind_speed`` in m/s."""
return self._interp(wind_speed, self.ct_values)
def _to_windio_turbine(self, name: str | None = None, n_points: int = 50) -> dict:
"""Return the windIO turbine mapping of the tabulated curves.
``n_points`` is ignored because the tabulated speeds are exported as given.
"""
return self._windio_mapping(
name, self.wind_speeds, self.power_values, self.ct_values
)
def _compress_turbines(turbines, n_turbines: int):
"""Return unique types and position-aligned indices in first-use order."""
if turbines is None:
return (), np.empty(0, dtype=int)
values = (turbines,) * n_turbines if isinstance(turbines, Turbine) else turbines
if len(values) != n_turbines:
raise ValueError("number of turbines must match the number of positions")
unique = []
indices = []
identities = {}
for turbine in values:
if not isinstance(turbine, Turbine):
raise TypeError("turbines must contain only Turbine objects")
key = id(turbine)
if key not in identities:
identities[key] = len(unique)
unique.append(turbine)
indices.append(identities[key])
return tuple(unique), np.asarray(indices, dtype=int)
def _from_windio_turbines(turbine_types, type_indices):
"""One shared turbine, or one per position when the types differ."""
if not isinstance(turbine_types, dict):
raise ValueError("windIO turbine_types must be a mapping")
try:
keys = sorted(int(key) for key in turbine_types)
except (TypeError, ValueError):
raise ValueError("windIO turbine_types keys must be integers") from None
types = {}
for key in keys:
entry = turbine_types.get(str(key), turbine_types.get(key))
if not isinstance(entry, dict):
raise ValueError(f"windIO turbine type {key} must be a mapping")
for field in ("performance", "rotor_diameter", "hub_height"):
if field not in entry:
raise ValueError(f"windIO turbine type {key} is missing {field!r}")
types[key] = _turbine_from_windio_performance(entry)
message = f"windIO layout turbine_types must refer to defined types {keys}"
# int() would read 0.7 as type 0 and True as type 1.
if any(isinstance(index, (bool, float)) for index in type_indices):
raise ValueError(message)
try:
values = tuple(types[int(index)] for index in type_indices)
except (KeyError, TypeError, ValueError):
raise ValueError(message) from None
return values[0] if values and all(item is values[0] for item in values) else values
def _operating_grid(cutin: float, rated: float, cutout: float, ct_speeds) -> np.ndarray:
"""Half-meter-per-second operating grid that keeps the tabulated knots."""
speeds = np.concatenate(
[
np.arange(cutin, cutout, 0.5),
np.asarray(ct_speeds, dtype=float),
np.asarray([cutin, cutout, rated], dtype=float),
]
)
speeds = speeds[(speeds >= cutin) & (speeds <= cutout)]
return np.unique(np.round(speeds, 6))
def _windio_curve(curve, curve_key: str, speeds_key: str, values_key: str):
"""Return a windIO performance curve as equal-length float arrays."""
from ._windio import _number_list
if not isinstance(curve, dict):
raise ValueError(f"windIO {curve_key} must be a mapping")
for key in (speeds_key, values_key):
if key not in curve:
raise ValueError(f"windIO {curve_key} is missing {key!r}")
speeds = _number_list(curve[speeds_key], f"{curve_key} {speeds_key}")
values = _number_list(curve[values_key], f"{curve_key} {values_key}")
if len(values) != len(speeds):
raise ValueError(
f"windIO {curve_key} needs one of {values_key} per {speeds_key}"
)
# Interpolation reads these speeds, and np.interp does not check them.
if not np.all(np.isfinite(speeds)) or np.any(np.diff(speeds) <= 0):
raise ValueError(
f"windIO {curve_key} {speeds_key} must be finite and strictly increasing"
)
return speeds, values
def _turbine_from_windio_performance(entry) -> Tabulated:
"""Build a :class:`Tabulated` from any windIO performance form.
A ``Cp`` curve gives ``P = 0.5 rho A C_P u^3`` at sea-level density. The
rated form rises cubically to rated power. Thrust is interpolated from the
``Ct`` curve. Power and ``Ct`` curves on different speeds are linearly
interpolated onto the union of both, each reading zero outside its own
speeds as a :class:`Tabulated` curve does.
"""
from ._windio import _number
performance = entry["performance"]
if not isinstance(performance, dict):
raise ValueError("windIO turbine performance must be a mapping")
ct = performance.get("Ct_curve")
if not isinstance(ct, dict):
raise ValueError("windIO turbine performance must contain a Ct curve")
rotor_diameter = _number(entry["rotor_diameter"], "turbine rotor_diameter")
shared = {
"rotor_diameter": rotor_diameter,
"hub_height": _number(entry["hub_height"], "turbine hub_height"),
"name": entry.get("name"),
}
ct_speeds, ct_values = _windio_curve(ct, "Ct_curve", "Ct_wind_speeds", "Ct_values")
power = performance.get("power_curve")
if power is not None:
speeds, power_values = _windio_curve(
power, "power_curve", "power_wind_speeds", "power_values"
)
if not np.array_equal(speeds, ct_speeds):
power_speeds = speeds
speeds = np.union1d(power_speeds, ct_speeds)
power_values = np.interp(
speeds, power_speeds, power_values, left=0.0, right=0.0
)
ct_values = np.interp(speeds, ct_speeds, ct_values, left=0.0, right=0.0)
return Tabulated(
wind_speeds=speeds,
power_values=power_values,
ct_values=ct_values,
**shared,
)
area = np.pi * (rotor_diameter / 2.0) ** 2
efficiency = _number(
performance.get("generator_efficiency", 1.0), "generator_efficiency"
)
cp_curve = performance.get("Cp_curve")
if cp_curve is not None:
cp_speeds, cp_values = _windio_curve(
cp_curve, "Cp_curve", "Cp_wind_speeds", "Cp_values"
)
speeds = np.union1d(cp_speeds, ct_speeds)
cp = np.interp(speeds, cp_speeds, cp_values, left=0.0, right=0.0)
return Tabulated(
wind_speeds=speeds,
power_values=0.5 * _AIR_DENSITY * area * cp * speeds**3 * efficiency,
ct_values=np.interp(speeds, ct_speeds, ct_values, left=0.0, right=0.0),
**shared,
)
required = (
"rated_power",
"rated_wind_speed",
"cutin_wind_speed",
"cutout_wind_speed",
)
if any(name not in performance for name in required):
raise ValueError(
"windIO turbine performance needs a power curve, a Cp curve, or "
"rated power with its operating interval"
)
rated_power, rated_speed, cutin, cutout = (
_number(performance[name], name) for name in required
)
speeds = _operating_grid(cutin, rated_speed, cutout, ct_speeds)
below = speeds < rated_speed
power_values = np.full(speeds.shape, rated_power)
power_values[below] = rated_power * (speeds[below] / rated_speed) ** 3
return Tabulated(
wind_speeds=speeds,
power_values=power_values,
ct_values=np.interp(speeds, ct_speeds, ct_values, left=0.0, right=0.0),
rated_power=rated_power,
**shared,
)
