Allow yaw optimization with disabled turbines

examples/examples_control_optimization/008_optimize_yaw_with_disabled_turbines.py

@@ -0,0 +1,48 @@
+"""Example: Optimizing yaw angles with disabled turbines
+
+This example demonstrates how to optimize yaw angles in FLORIS, when some turbines are disabled.
+The example optimization is run using both YawOptimizerSR and YawOptimizerGeometric, the two
+yaw optimizers that support disabling turbines.
+"""
+
+import numpy as np
+
+from floris import FlorisModel
+from floris.optimization.yaw_optimization.yaw_optimizer_geometric import YawOptimizationGeometric
+from floris.optimization.yaw_optimization.yaw_optimizer_sr import YawOptimizationSR
+
+
+# Load a 3-turbine model
+fmodel = FlorisModel("../inputs/gch.yaml")
+
+# Set wind conditions to be the same for two cases
+fmodel.set(wind_directions=[270.]*2, wind_speeds=[8.]*2, turbulence_intensities=[.06]*2)
+
+# First run the case where all turbines are active and print results
+yaw_opt = YawOptimizationSR(fmodel)
+df_opt = yaw_opt.optimize()
+print("Serial Refine optimized yaw angles (all turbines active) [deg]:\n", df_opt.yaw_angles_opt)
+
+yaw_opt = YawOptimizationGeometric(fmodel)
+df_opt = yaw_opt.optimize()
+print("\nGeometric optimized yaw angles (all turbines active) [deg]:\n", df_opt.yaw_angles_opt)
+
+# Disable turbines (different pattern for each of the two cases)
+# First case: disable the middle turbine
+# Second case: disable the front turbine
+fmodel.set_operation_model('mixed')
+fmodel.set(disable_turbines=np.array([[False, True, False], [True, False, False]]))
+
+# Rerun optimizations and print results
+yaw_opt = YawOptimizationSR(fmodel)
+df_opt = yaw_opt.optimize()
+print(
+    "\nSerial Refine optimized yaw angles (some turbines disabled) [deg]:\n",
+    df_opt.yaw_angles_opt
+)
+# Note that disabled turbines may not have a zero yaw angle, because a disabled turbine's
+# yaw angle does not affect the total power output.
+
+yaw_opt = YawOptimizationGeometric(fmodel)
+df_opt = yaw_opt.optimize()
+print("\nGeometric optimized yaw angles (some turbines disabled) [deg]:\n", df_opt.yaw_angles_opt)

floris/core/turbine/operation_models.py

@@ -382,22 +382,22 @@ def axial_induction(
 @define
 class MixedOperationTurbine(BaseOperationModel):
 
+    @staticmethod
     def power(
         yaw_angles: NDArrayFloat,
         power_setpoints: NDArrayFloat,
         **kwargs
     ):
-        # Yaw angles mask all yaw_angles not equal to zero
-        yaw_angles_mask = yaw_angles != 0.0
-        power_setpoints_mask = power_setpoints < POWER_SETPOINT_DEFAULT
-        neither_mask = np.logical_not(yaw_angles_mask) & np.logical_not(power_setpoints_mask)
-
-        if (power_setpoints_mask & yaw_angles_mask).any():
-            raise ValueError((
-                "Power setpoints and yaw angles are incompatible."
-                "If yaw_angles entry is nonzero, power_setpoints must be greater than"
-                " or equal to {0}.".format(POWER_SETPOINT_DEFAULT)
-            ))
+        (
+            yaw_angles,
+            power_setpoints,
+            yaw_angles_mask,
+            power_setpoints_mask,
+            neither_mask
+        ) = MixedOperationTurbine._handle_mixed_operation_setpoints(
+            yaw_angles=yaw_angles,
+            power_setpoints=power_setpoints
+        )
 
         powers = np.zeros_like(power_setpoints)
         powers[yaw_angles_mask] += CosineLossTurbine.power(
@@ -414,21 +414,22 @@ def power(
 
         return powers
 
+    @staticmethod
     def thrust_coefficient(
         yaw_angles: NDArrayFloat,
         power_setpoints: NDArrayFloat,
         **kwargs
     ):
-        yaw_angles_mask = yaw_angles != 0.0
-        power_setpoints_mask = power_setpoints < POWER_SETPOINT_DEFAULT
-        neither_mask = np.logical_not(yaw_angles_mask) & np.logical_not(power_setpoints_mask)
-
-        if (power_setpoints_mask & yaw_angles_mask).any():
-            raise ValueError((
-                "Power setpoints and yaw angles are incompatible."
-                "If yaw_angles entry is nonzero, power_setpoints must be greater than"
-                " or equal to {0}.".format(POWER_SETPOINT_DEFAULT)
-            ))
+        (
+            yaw_angles,
+            power_setpoints,
+            yaw_angles_mask,
+            power_setpoints_mask,
+            neither_mask
+        ) = MixedOperationTurbine._handle_mixed_operation_setpoints(
+            yaw_angles=yaw_angles,
+            power_setpoints=power_setpoints
+        )
 
         thrust_coefficients = np.zeros_like(power_setpoints)
         thrust_coefficients[yaw_angles_mask] += CosineLossTurbine.thrust_coefficient(
@@ -445,21 +446,22 @@ def thrust_coefficient(
 
         return thrust_coefficients
 
+    @staticmethod
     def axial_induction(
         yaw_angles: NDArrayFloat,
         power_setpoints: NDArrayFloat,
         **kwargs
     ):
-        yaw_angles_mask = yaw_angles != 0.0
-        power_setpoints_mask = power_setpoints < POWER_SETPOINT_DEFAULT
-        neither_mask = np.logical_not(yaw_angles_mask) & np.logical_not(power_setpoints_mask)
-
-        if (power_setpoints_mask & yaw_angles_mask).any():
-            raise ValueError((
-                "Power setpoints and yaw angles are incompatible."
-                "If yaw_angles entry is nonzero, power_setpoints must be greater than"
-                " or equal to {0}.".format(POWER_SETPOINT_DEFAULT)
-            ))
+        (
+            yaw_angles,
+            power_setpoints,
+            yaw_angles_mask,
+            power_setpoints_mask,
+            neither_mask
+        ) = MixedOperationTurbine._handle_mixed_operation_setpoints(
+            yaw_angles=yaw_angles,
+            power_setpoints=power_setpoints
+        )
 
         axial_inductions = np.zeros_like(power_setpoints)
         axial_inductions[yaw_angles_mask] += CosineLossTurbine.axial_induction(
@@ -476,6 +478,34 @@ def axial_induction(
 
         return axial_inductions
 
+    @staticmethod
+    def _handle_mixed_operation_setpoints(
+        yaw_angles: NDArrayFloat,
+        power_setpoints: NDArrayFloat,
+    ):
+        """
+        Check for incompatible yaw angles and power setpoints and raise an error if found.
+        Return masks and updated setpoints.
+        """
+        # If any turbines are disabled, set their yaw angles to zero
+        yaw_angles[power_setpoints <= POWER_SETPOINT_DISABLED] = 0.0
+
+        # Create masks for whether yaw angles and power setpoints are set
+        yaw_angles_mask = yaw_angles != 0.0
+        power_setpoints_mask = power_setpoints < POWER_SETPOINT_DEFAULT
+        neither_mask = np.logical_not(yaw_angles_mask) & np.logical_not(power_setpoints_mask)
+
+        # Check for incompatibility and raise error if found.
+        if (power_setpoints_mask & yaw_angles_mask).any():
+            raise ValueError((
+                "Power setpoints and yaw angles are incompatible."
+                "If yaw_angles entry is nonzero, power_setpoints must be greater than"
+                " or equal to {0}.".format(POWER_SETPOINT_DEFAULT)
+            ))
+
+        # Return updated setpoints as well as masks
+        return yaw_angles, power_setpoints, yaw_angles_mask, power_setpoints_mask, neither_mask
+
 @define
 class AWCTurbine(BaseOperationModel):
     """

floris/floris_model.py

@@ -453,10 +453,7 @@ def set(
         # previous setting
         if not (_yaw_angles == 0).all():
             self.core.farm.set_yaw_angles(_yaw_angles)
-        if not (
-            (_power_setpoints == POWER_SETPOINT_DEFAULT)
-            | (_power_setpoints == POWER_SETPOINT_DISABLED)
-        ).all():
+        if not (_power_setpoints == POWER_SETPOINT_DEFAULT).all():
             self.core.farm.set_power_setpoints(_power_setpoints)
         if _awc_modes is not None:
             self.core.farm.set_awc_modes(_awc_modes)

floris/optimization/yaw_optimization/yaw_optimization_base.py

@@ -99,7 +99,7 @@ def __init__(
         """
 
         # Save turbine object to self
-        self.fmodel = fmodel.copy()
+        self.fmodel = copy.deepcopy(fmodel)
         self.nturbs = len(self.fmodel.layout_x)
 
         # # Check floris options
@@ -224,7 +224,7 @@ def _reduce_control_problem(self):
         self.turbs_to_opt = (self.maximum_yaw_angle - self.minimum_yaw_angle >= 0.001)
 
         # Initialize subset variables as full set
-        self.fmodel_subset = self.fmodel.copy()
+        self.fmodel_subset = copy.deepcopy(self.fmodel)
         n_findex_subset = copy.deepcopy(self.fmodel.core.flow_field.n_findex)
         minimum_yaw_angle_subset = copy.deepcopy(self.minimum_yaw_angle)
         maximum_yaw_angle_subset = copy.deepcopy(self.maximum_yaw_angle)
@@ -301,6 +301,7 @@ def _calculate_farm_power(
             ti_array=None,
             turbine_weights=None,
             heterogeneous_speed_multipliers=None,
+            power_setpoints=None,
         ):
         """
         Calculate the wind farm power production assuming the predefined
@@ -353,6 +354,7 @@ def _calculate_farm_power(
             wind_speeds=ws_array,
             turbulence_intensities=ti_array,
             yaw_angles=yaw_angles,
+            power_setpoints=power_setpoints,
         )
         fmodel_subset.run()
         turbine_power = fmodel_subset.get_turbine_powers()

floris/optimization/yaw_optimization/yaw_optimizer_geometric.py

@@ -1,6 +1,7 @@
 
 import numpy as np
 
+from floris.core.turbine.operation_models import POWER_SETPOINT_DISABLED
 from floris.utilities import rotate_coordinates_rel_west
 
 from .yaw_optimization_base import YawOptimization
@@ -46,10 +47,11 @@ def optimize(self):
         # Loop through every WD individually. WS ignored!
         wd_array = self.fmodel_subset.core.flow_field.wind_directions
 
+        active_turbines = self.fmodel_subset.core.farm.power_setpoints > POWER_SETPOINT_DISABLED
         for nwdi, wd in enumerate(wd_array):
-            self._yaw_angles_opt_subset[nwdi, :] = geometric_yaw(
-                self.fmodel_subset.layout_x,
-                self.fmodel_subset.layout_y,
+            self._yaw_angles_opt_subset[nwdi, active_turbines[nwdi]] = geometric_yaw(
+                self.fmodel_subset.layout_x[active_turbines[nwdi]],
+                self.fmodel_subset.layout_y[active_turbines[nwdi]],
                 wd,
                 self.fmodel.core.farm.turbine_definitions[0]["rotor_diameter"],
                 top_left_yaw_upper=self.maximum_yaw_angle[0, 0],

floris/optimization/yaw_optimization/yaw_optimizer_scipy.py

@@ -30,6 +30,12 @@ def __init__(
         Instantiate YawOptimizationScipy object with a FlorisModel object
         and assign parameter values.
         """
+        valid_op_models = ["cosine-loss"]
+        if fmodel.get_operation_model() not in valid_op_models:
+            raise ValueError(
+                "YawOptimizationScipy is currently limited to the following operation models: "
+                + ", ".join(valid_op_models)
+            )
         if opt_options is None:
             # Default SciPy parameters
             opt_options = {

floris/optimization/yaw_optimization/yaw_optimizer_sr.py

@@ -93,6 +93,7 @@ def _calc_powers_with_memory(self, yaw_angles_subset, use_memory=True):
         wd_array_subset = self.fmodel_subset.core.flow_field.wind_directions
         ws_array_subset = self.fmodel_subset.core.flow_field.wind_speeds
         ti_array_subset = self.fmodel_subset.core.flow_field.turbulence_intensities
+        power_setpoints_subset = self.fmodel_subset.core.farm.power_setpoints
         turbine_weights_subset = self._turbine_weights_subset
 
         # Reformat yaw_angles_subset, if necessary
@@ -108,6 +109,7 @@ def _calc_powers_with_memory(self, yaw_angles_subset, use_memory=True):
             wd_array_subset = np.tile(wd_array_subset, Ny)
             ws_array_subset = np.tile(ws_array_subset, Ny)
             ti_array_subset = np.tile(ti_array_subset, Ny)
+            power_setpoints_subset = np.tile(power_setpoints_subset, (Ny, 1))
             turbine_weights_subset = np.tile(turbine_weights_subset, (Ny, 1))
 
         # Initialize empty matrix for floris farm power outputs
@@ -143,7 +145,8 @@ def _calc_powers_with_memory(self, yaw_angles_subset, use_memory=True):
                 ti_array=ti_array_subset[~idx],
                 turbine_weights=turbine_weights_subset[~idx, :],
                 yaw_angles=yaw_angles_subset[~idx, :],
-                heterogeneous_speed_multipliers=het_sm
+                heterogeneous_speed_multipliers=het_sm,
+                power_setpoints=power_setpoints_subset[~idx, :],
             )
             self.time_spent_in_floris += (timerpc() - start_time)