from collections.abc import Sequence
from typing import Any
import foxes.constants as FC
import foxes.variables as FV
import numpy as np
from foxes_opt.core.farm_constraint import FarmConstraint
from foxes_opt.core.farm_opt_problem import FarmOptProblem
[docs]
class MinDistConstraint(FarmConstraint):
"""
Turbines must keep at least a minimal
spatial distance.
"""
def __init__(
self,
problem: FarmOptProblem,
min_dist: float,
min_dist_unit: str = "m",
name: str = "dist",
sel_turbines: list[int] | None = None,
infer_vars: bool = True,
check_only_selected: bool = False,
**kwargs: Any,
) -> None:
"""
Parameters
----------
problem
The underlying optimization problem
min_dist
The minimal distance
min_dist_unit
The minimal distance unit, either m or D
name
The name of the constraint
sel_turbines
The selected turbines
infer_vars
Whether to infer the float variables automatically
check_only_selected
Whether to check distances only among selected turbines
kwargs
Additional parameters for `iwopy.Constraint`
"""
self.min_dist: float = min_dist
self.min_dist_unit: str = min_dist_unit
self.check_only_selected = check_only_selected
selt = problem.sel_turbines if sel_turbines is None else sel_turbines
vrs = []
for ti in selt:
vrs += [problem.tvar(FV.X, ti), problem.tvar(FV.Y, ti)]
vnames_float: list[str] | None = vrs
if infer_vars:
vnames_float = None
super().__init__(
problem, name, sel_turbines, vnames_float=vnames_float, **kwargs
)
[docs]
def initialize(self, verbosity: int = 0) -> None:
"""
Initialize the constaint.
Parameters
----------
verbosity
The verbosity level, 0 = silent
"""
N = self.farm.n_turbines
i2t: list[list[int]] = [] # i --> (ti, tj)
self._t2i: np.ndarray[tuple[int, int], np.dtype[np.int_]] = np.full(
[N, N], -1
) # (ti, tj) --> i
i = 0
for ti in self.sel_turbines:
other_turbines = self.sel_turbines if self.check_only_selected else range(N)
for tj in other_turbines:
if ti != tj and self._t2i[ti, tj] < 0:
i2t.append([ti, tj])
self._t2i[ti, tj] = i
self._t2i[tj, ti] = i
i += 1
self._i2t: np.ndarray[tuple[int, int], np.dtype[np.int_]] = np.asarray(
i2t, dtype=int
)
self._cnames: list[str] = [f"{self.name}_{ti}_{tj}" for ti, tj in self._i2t]
super().initialize(verbosity)
[docs]
def n_components(self) -> int:
"""
Returns the number of components of the
function.
Returns
-------
value
The number of components.
"""
return len(self._i2t)
[docs]
def active_components(
self, vars_float: np.ndarray, lower: np.ndarray, upper: np.ndarray
) -> np.ndarray:
"""Select pairs that can approach the spacing limit within bounds.
Parameters
----------
vars_float
Initial values in this constraint's float-variable order.
lower, upper
Enforced coordinate bounds in the same order, in metres.
Unselected turbine positions remain fixed.
Returns
-------
components
Stable pair-component indices whose coordinate boxes can
approach the nominal minimum distance. Omitted pairs are
separated throughout the box. Component metadata, full
values and feasibility checks remain unchanged. Unbounded
coordinates are handled conservatively. Unsupported variable
parameterizations, variable diameters and state-dependent
coordinates retain all pairs.
Raises
------
ValueError
If the minimum-distance unit is neither m nor D.
"""
coordinates = [
np.asarray(turbine.xy, dtype=np.float64) for turbine in self.farm.turbines
]
if any(
coordinate.shape != (2,) or not np.all(np.isfinite(coordinate))
for coordinate in coordinates
):
return np.arange(self.n_components(), dtype=int)
positions = np.asarray(coordinates)
point_lower, point_upper = positions.copy(), positions.copy()
for variable_index, name in enumerate(self.var_names_float):
try:
coordinate, turbine_index = self.problem.parse_tvar(name)
except ValueError:
return np.arange(self.n_components(), dtype=int)
if coordinate not in (FV.X, FV.Y):
return np.arange(self.n_components(), dtype=int)
axis = 0 if coordinate == FV.X else 1
point_lower[turbine_index, axis] = lower[variable_index]
point_upper[turbine_index, axis] = upper[variable_index]
pairs = self._i2t
gap = np.maximum(
np.maximum(
point_lower[pairs[:, 0]] - point_upper[pairs[:, 1]],
point_lower[pairs[:, 1]] - point_upper[pairs[:, 0]],
),
0.0,
)
separation = np.hypot(gap[:, 0], gap[:, 1])
if self.min_dist_unit == "m":
required = np.full(len(pairs), self.min_dist)
elif self.min_dist_unit == "D":
diameters = self.farm.get_rotor_diameters(self.problem.algo)
required = self.min_dist * np.maximum(
diameters[pairs[:, 0]], diameters[pairs[:, 1]]
)
else:
raise ValueError(
f"Constraint '{self.name}': Unknown min_dist_unit '{self.min_dist_unit}'"
)
reserve = 32.0 * np.finfo(np.float64).eps * np.maximum(required, 1.0)
return np.flatnonzero(separation <= required + reserve)
[docs]
def vardeps_float(self) -> np.ndarray[tuple[int, int], np.dtype[Any]]:
"""
Gets the dependencies of all components
on the function float variables
Returns
-------
deps
The dependencies of components on function
variables, shape
"""
turbines = np.asarray(self.problem.sel_turbines, dtype=int)
deps = np.zeros((self.n_components(), 2 * len(turbines)), dtype=bool)
turbine_to_var = np.full(self.farm.n_turbines, -1, dtype=int)
turbine_to_var[turbines] = np.arange(len(turbines))
for pair_column in range(2):
var_indices = turbine_to_var[self._i2t[:, pair_column]]
rows = np.flatnonzero(var_indices >= 0)
columns = 2 * var_indices[rows]
deps[rows, columns] = True
deps[rows, columns + 1] = True
return deps
[docs]
def ana_deriv(
self,
vars_int: np.ndarray,
vars_float: np.ndarray,
var: int,
components: Sequence[int] | np.ndarray | None = None,
) -> np.ndarray:
"""
Calculate the analytical derivative for one float variable.
Parameters
----------
vars_int
The integer variable values
vars_float
The float variable values
var
The float variable index
components
The selected components, or None for all
Returns
-------
deriv
The derivative values, shape: (n_sel_components,)
"""
cidx = (
np.arange(self.n_components(), dtype=int)
if components is None
else np.asarray(components, dtype=int)
)
deriv = np.zeros(len(cidx), dtype=np.float64)
try:
vname, ti = self.problem.parse_tvar(self.var_names_float[var])
except (IndexError, ValueError):
return np.full(len(cidx), np.nan, dtype=np.float64)
if vname not in (FV.X, FV.Y):
return np.full(len(cidx), np.nan, dtype=np.float64)
pairs = self._i2t[cidx]
affected = np.any(pairs == ti, axis=1)
if not np.any(affected):
return deriv
xy = np.asarray(
[np.asarray(t.xy).reshape(-1, 2)[0] for t in self.farm.turbines],
dtype=np.float64,
)
for name, value in zip(self.var_names_float, vars_float):
try:
xyname, xyti = self.problem.parse_tvar(name)
except (IndexError, ValueError):
continue
if xyname in (FV.X, FV.Y):
xy[xyti, 0 if xyname == FV.X else 1] = value
apairs = pairs[affected]
delta = xy[apairs[:, 0]] - xy[apairs[:, 1]]
distance = np.linalg.norm(delta, axis=1)
coord = 0 if vname == FV.X else 1
values = np.full(len(distance), np.nan, dtype=np.float64)
nz = distance > 0.0
sign = np.where(apairs[:, 0] == ti, -1.0, 1.0)
values[nz] = sign[nz] * delta[nz, coord] / distance[nz]
deriv[affected] = values
return deriv
[docs]
def calc_individual(
self,
vars_int: np.ndarray,
vars_float: np.ndarray,
problem_results: Any,
components: list[int] | None = None,
) -> np.ndarray:
"""
Calculate values for a single individual of the
underlying problem.
Parameters
----------
vars_int
The integer variable values, shape: (n_vars_int,)
vars_float
The float variable values, shape: (n_vars_float,)
problem_results
The results of the variable application
to the problem
components
The selected components or None for all
Returns
-------
values
The component values, shape: (n_sel_components,)
"""
xy = np.stack(
[problem_results[FV.X].to_numpy(), problem_results[FV.Y].to_numpy()],
axis=-1,
)
if not np.all(np.abs(np.min(xy, axis=0) - np.max(xy, axis=0)) < 1e-13):
raise ValueError(f"Constraint '{self.name}': Require state independet XY")
xy = xy[0]
cidx = np.arange(self.n_components(), dtype=int)
if components is not None and len(components) < self.n_components():
cidx = np.asarray(components, dtype=int)
a: np.ndarray[tuple[int, ...], np.dtype[Any]] = np.take_along_axis(
xy, self._i2t[cidx][:, 0, None], axis=0
)
b: np.ndarray[tuple[int, ...], np.dtype[Any]] = np.take_along_axis(
xy, self._i2t[cidx][:, 1, None], axis=0
)
d = np.linalg.norm(a - b, axis=-1)
if self.min_dist_unit == "m":
mind: Any = self.min_dist
elif self.min_dist_unit == "D":
D = problem_results[FV.D].to_numpy()
if not np.all(np.abs(np.min(D, axis=0) - np.max(D, axis=0)) < 1e-13):
raise ValueError(
f"Constraint '{self.name}': Require state independet D"
)
D = D[0]
Da = np.take_along_axis(D, self._i2t[cidx][:, 0], axis=0)
Db = np.take_along_axis(D, self._i2t[cidx][:, 1], axis=0)
mind = self.min_dist * np.maximum(Da, Db)
return mind - d
[docs]
def calc_population(
self,
vars_int: np.ndarray,
vars_float: np.ndarray,
problem_results: Any,
components: list[int] | None = None,
) -> np.ndarray:
"""
Calculate values for all individuals of a population.
Parameters
----------
vars_int
The integer variable values, shape: (n_pop, n_vars_int)
vars_float
The float variable values, shape: (n_pop, n_vars_float)
problem_results
The results of the variable application
to the problem
components
The selected components or None for all
Returns
-------
values
The component values, shape: (n_pop, n_sel_components)
"""
n_pop = problem_results["n_pop"].values
n_states = problem_results["n_org_states"].values
n_turbines = problem_results.sizes[FC.TURBINE]
xy = np.stack(
[problem_results[FV.X].to_numpy(), problem_results[FV.Y].to_numpy()],
axis=-1,
)
xy = xy.reshape(n_pop, n_states, n_turbines, 2)
if not np.all(np.abs(np.min(xy, axis=1) - np.max(xy, axis=1)) < 1e-13):
raise ValueError(f"Constraint '{self.name}': Require state independet XY")
xy = xy[:, 0]
cidx = np.arange(self.n_components(), dtype=int)
if components is not None and len(components) < self.n_components():
cidx = np.asarray(components, dtype=int)
a: np.ndarray[tuple[int, ...], np.dtype[Any]] = np.take_along_axis(
xy, self._i2t[cidx][None, :, 0, None], axis=1
)
b: np.ndarray[tuple[int, ...], np.dtype[Any]] = np.take_along_axis(
xy, self._i2t[cidx][None, :, 1, None], axis=1
)
d = np.linalg.norm(a - b, axis=-1)
if self.min_dist_unit == "m":
mind: Any = self.min_dist
elif self.min_dist_unit == "D":
D = problem_results[FV.D].to_numpy().reshape(n_pop, n_states, n_turbines)
if not np.all(np.abs(np.min(D, axis=1) - np.max(D, axis=1)) < 1e-13):
raise ValueError(
f"Constraint '{self.name}': Require state independet D"
)
D = D[:, 0]
Da = np.take_along_axis(D, self._i2t[cidx][None, :, 0], axis=1)
Db = np.take_along_axis(D, self._i2t[cidx][None, :, 1], axis=1)
mind = self.min_dist * np.maximum(Da, Db)
return mind - d
[docs]
class MinDistLocalConstraint(MinDistConstraint):
"""
Turbines must keep at least a minimal spatial distance,
checking only pairs that can violate the distance in a local
move problem.
"""
def __init__(
self,
problem: FarmOptProblem,
min_dist: float,
min_dist_unit: str = "m",
name: str = "dist",
sel_turbines: list[int] | None = None,
infer_vars: bool = True,
check_only_selected: bool = False,
max_move: float | None = None,
**kwargs: Any,
) -> None:
"""
Parameters
----------
problem
The underlying local move optimization problem
min_dist
The minimal distance
min_dist_unit
The minimal distance unit, either m or D
name
The name of the constraint
sel_turbines
The selected turbines
infer_vars
Whether to infer the float variables automatically
check_only_selected
Whether to check distances only among selected turbines
max_move
The maximal local displacement per selected turbine. If None,
infer from problem.radius or problem.square_length.
kwargs
Additional parameters for `iwopy.Constraint`
"""
self.max_move = max_move
super().__init__(
problem,
min_dist,
min_dist_unit,
name,
sel_turbines,
infer_vars,
check_only_selected,
**kwargs,
)
[docs]
def initialize(self, verbosity: int = 0) -> None:
"""
Initialize the constaint.
Parameters
----------
verbosity
The verbosity level, 0 = silent
"""
if self.max_move is None:
if hasattr(self.problem, "radius"):
max_move = self.problem.radius
elif hasattr(self.problem, "square_length"):
max_move = np.sqrt(2) * self.problem.square_length / 2
else:
raise ValueError(
f"Constraint '{self.name}': Failed to infer local movement limit"
)
else:
max_move = self.max_move
N = self.farm.n_turbines
xy0 = np.asarray([t.xy for t in self.farm.turbines], dtype=float)
move = np.zeros(N, dtype=float)
move[self.problem.sel_turbines] = max_move
if self.min_dist_unit == "m":
mind = np.full((N, N), self.min_dist, dtype=float)
elif self.min_dist_unit == "D":
D = self.farm.get_rotor_diameters(self.problem.algo)
mind = self.min_dist * np.maximum(D[:, None], D[None, :])
else:
raise ValueError(
f"Constraint '{self.name}': Unknown min_dist_unit '{self.min_dist_unit}'"
)
i2t: list[list[int]] = [] # i --> (ti, tj)
self._t2i: np.ndarray[tuple[int, int], np.dtype[np.int_]] = np.full(
[N, N], -1
) # (ti, tj) --> i
i = 0
for ti in self.sel_turbines:
other_turbines = self.sel_turbines if self.check_only_selected else range(N)
for tj in other_turbines:
if ti != tj and self._t2i[ti, tj] < 0:
dist0 = np.linalg.norm(xy0[ti] - xy0[tj])
if dist0 <= mind[ti, tj] + move[ti] + move[tj]:
i2t.append([ti, tj])
self._t2i[ti, tj] = i
self._t2i[tj, ti] = i
i += 1
self._i2t: np.ndarray[tuple[int, int], np.dtype[np.int_]] = np.asarray(
i2t, dtype=int
).reshape((-1, 2))
self._cnames: list[str] = [f"{self.name}_{ti}_{tj}" for ti, tj in self._i2t]
FarmConstraint.initialize(self, verbosity)