Allow yaw optimization with disabled turbines (#1027)

* Update MixedOperationModel to better handle disabling turbines, and remove disabled turbines from default power_setpoints on FlorisModel.

* Update yaw optimization to pass through power_setpoints.

* Raise warning for incompatible op model in scipy optimizer

* Enable geomtric yaw optimization with disabled turbines.

* Add example of optimizing with disabled turbines.

* Update example to include all active case.

* Add tests for yaw opt with turbines disabled, as well as basic serial refine operation.

* Add similar test for Geometric yaw (only farm power comparisons removed).

* Formatting fixes.

* Limit yaw angles to zero for disabled turbines.

* Update comment on disabled turbines' yaw angles.

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 are assigned a zero yaw angle, but their yaw angle is arbitrary as it
+# 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

@@ -5,6 +5,7 @@
 import numpy as np
 import pandas as pd
 
+from floris.core.turbine.operation_models import POWER_SETPOINT_DISABLED
 from floris.logging_manager import LoggingManager
 
 from .yaw_optimization_tools import derive_downstream_turbines
@@ -99,7 +100,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
@@ -131,6 +132,11 @@ def __init__(
         self.minimum_yaw_angle = self._unpack_variable(minimum_yaw_angle)
         self.maximum_yaw_angle = self._unpack_variable(maximum_yaw_angle)
 
+        # Limit yaw angles to zero for disabled turbines
+        active_turbines = fmodel.core.farm.power_setpoints > POWER_SETPOINT_DISABLED
+        self.minimum_yaw_angle[~active_turbines] = 0.0
+        self.maximum_yaw_angle[~active_turbines] = 0.0
+
         # Set initial condition for optimization
         if x0 is not None:
             self.x0 = self._unpack_variable(x0)
@@ -224,7 +230,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 +307,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 +360,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)
 

tests/floris_model_integration_test.py

@@ -10,7 +10,7 @@
     TimeSeries,
     WindRose,
 )
-from floris.core.turbine.operation_models import POWER_SETPOINT_DEFAULT
+from floris.core.turbine.operation_models import POWER_SETPOINT_DEFAULT, POWER_SETPOINT_DISABLED
 
 
 TEST_DATA = Path(__file__).resolve().parent / "data"
@@ -44,7 +44,24 @@ def test_assign_setpoints():
 
     # power_setpoints and disable_turbines (disable_turbines overrides power_setpoints)
     fmodel.set(power_setpoints=[[1e6, 2e6]], disable_turbines=[[True, False]])
-    assert np.allclose(fmodel.core.farm.power_setpoints, np.array([[0.001, 2e6]]))
+    assert np.allclose(fmodel.core.farm.power_setpoints, np.array([[POWER_SETPOINT_DISABLED, 2e6]]))
+
+    # Setting sequentially is equivalent to setting together
+    fmodel.reset_operation()
+    fmodel.set(disable_turbines=[[True, False]])
+    fmodel.set(yaw_angles=[[0, 30]])
+    assert np.allclose(
+        fmodel.core.farm.power_setpoints,
+        np.array([[POWER_SETPOINT_DISABLED, POWER_SETPOINT_DEFAULT]])
+    )
+    assert np.allclose(fmodel.core.farm.yaw_angles, np.array([[0, 30]]))
+
+    fmodel.set(disable_turbines=[[True, False]], yaw_angles=[[0, 30]])
+    assert np.allclose(
+        fmodel.core.farm.power_setpoints,
+        np.array([[POWER_SETPOINT_DISABLED, POWER_SETPOINT_DEFAULT]])
+    )
+    assert np.allclose(fmodel.core.farm.yaw_angles, np.array([[0, 30]]))
 
 def test_set_run():
     """