Add support to shut off turbines

examples/41_test_disable_turbines.py

@@ -0,0 +1,97 @@
+# Copyright 2023 NREL
+
+# Licensed under the Apache License, Version 2.0 (the "License"); you may not
+# use this file except in compliance with the License. You may obtain a copy of
+# the License at http://www.apache.org/licenses/LICENSE-2.0
+
+# Unless required by applicable law or agreed to in writing, software
+# distributed under the License is distributed on an "AS IS" BASIS, WITHOUT
+# WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. See the
+# License for the specific language governing permissions and limitations under
+# the License.
+
+# See https://floris.readthedocs.io for documentation
+
+# Example adapted from https://github.com/NREL/floris/pull/693 contributed by Elie Kadoche
+
+
+import matplotlib.pyplot as plt
+import numpy as np
+import yaml
+
+from floris.tools import FlorisInterface
+
+
+"""
+This example demonstrates the ability of FLORIS to shut down some turbines
+during a simulation.
+"""
+
+# Initialize the FLORIS interface
+fi = FlorisInterface("inputs/gch.yaml")
+
+# Change to the mixed model turbine
+with open(
+    str(
+        fi.floris.as_dict()["farm"]["turbine_library_path"]
+        / (fi.floris.as_dict()["farm"]["turbine_type"][0] + ".yaml")
+    )
+) as t:
+    turbine_type = yaml.safe_load(t)
+turbine_type["power_thrust_model"] = "mixed"
+fi.reinitialize(turbine_type=[turbine_type])
+
+# Consider a wind farm of 3 aligned wind turbines
+layout = np.array([[0.0, 0.0], [500.0, 0.0], [1000.0, 0.0]])
+
+# Run the computations for 2 identical wind data
+# (n_findex = 2)
+wind_directions = np.array([270.0, 270.0])
+wind_speeds = np.array([8.0, 8.0])
+
+# Shut down the first 2 turbines for the second findex
+# 2 findex x 3 turbines
+disable_turbines = np.array([[False, False, False], [True, True, False]])
+
+# Simulation
+# ------------------------------------------
+
+# Reinitialize flow field
+fi.reinitialize(
+    layout_x=layout[:, 0],
+    layout_y=layout[:, 1],
+    wind_directions=wind_directions,
+    wind_speeds=wind_speeds,
+)
+
+# # Compute wakes
+fi.calculate_wake(disable_turbines=disable_turbines)
+
+# Results
+# ------------------------------------------
+
+# Get powers and effective wind speeds
+turbine_powers = fi.get_turbine_powers()
+turbine_powers = np.round(turbine_powers * 1e-3, decimals=2)
+effective_wind_speeds = fi.turbine_average_velocities
+
+
+# Plot the results
+fig, axarr = plt.subplots(2, 1, sharex=True)
+
+# Plot the power
+ax = axarr[0]
+ax.plot(["T0", "T1", "T2"], turbine_powers[0, :], "ks-", label="All on")
+ax.plot(["T0", "T1", "T2"], turbine_powers[1, :], "ro-", label="T0 & T1 disabled")
+ax.set_ylabel("Power (kW)")
+ax.grid(True)
+ax.legend()
+
+ax = axarr[1]
+ax.plot(["T0", "T1", "T2"], effective_wind_speeds[0, :], "ks-", label="All on")
+ax.plot(["T0", "T1", "T2"], effective_wind_speeds[1, :], "ro-", label="T0 & T1 disabled")
+ax.set_ylabel("Effective wind speeds (m/s)")
+ax.grid(True)
+ax.legend()
+
+plt.show()

floris/tools/floris_interface.py

@@ -31,7 +31,11 @@
 )
 from floris.tools.cut_plane import CutPlane
 from floris.tools.wind_data import WindDataBase
-from floris.type_dec import floris_array_converter, NDArrayFloat
+from floris.type_dec import (
+    floris_array_converter,
+    NDArrayBool,
+    NDArrayFloat,
+)
 
 
 class FlorisInterface(LoggingManager):
@@ -122,6 +126,7 @@ def calculate_wake(
         yaw_angles: NDArrayFloat | list[float] | None = None,
         # tilt_angles: NDArrayFloat | list[float] | None = None,
         power_setpoints: NDArrayFloat | list[float] | list[float, None] | None = None,
+        disable_turbines: NDArrayBool | list[bool] | None = None,
     ) -> None:
         """
         Wrapper to the :py:meth:`~.Farm.set_yaw_angles` and
@@ -133,6 +138,9 @@ def calculate_wake(
             power_setpoints (NDArrayFloat | list[float] | None, optional): Turbine power setpoints.
                 May be specified with some float values and some None values; power maximization
                 will be assumed for any None value. Defaults to None.
+            disable_turbines (NDArrayBool | list[bool] | None, optional): NDArray with dimensions
+                n_findex x n_turbines.  True values indicate the turbine is disabled at that findex
+                and the power setpoint at that position is set to 0.  Defaults to None
         """
 
         if yaw_angles is None:
@@ -160,6 +168,33 @@ def calculate_wake(
         ] = POWER_SETPOINT_DEFAULT
         power_setpoints = floris_array_converter(power_setpoints)
 
+        # Check for turbines to disable
+        if disable_turbines is not None:
+
+            # Force to numpy array
+            disable_turbines = np.array(disable_turbines)
+
+            # Must have first dimension = n_findex
+            if disable_turbines.shape[0] != self.floris.flow_field.n_findex:
+                raise ValueError(
+                    f"disable_turbines has a size of {disable_turbines.shape[0]} "
+                    f"in the 0th dimension, must be equal to "
+                    f"n_findex={self.floris.flow_field.n_findex}"
+                )
+
+            # Must have first dimension = n_turbines
+            if disable_turbines.shape[1] != self.floris.farm.n_turbines:
+                raise ValueError(
+                    f"disable_turbines has a size of {disable_turbines.shape[1]} "
+                    f"in the 1th dimension, must be equal to "
+                    f"n_turbines={self.floris.farm.n_turbines}"
+                )
+
+            # Set power_setpoints and yaw_angles to 0 in all locations where
+            # disable_turbines is True
+            yaw_angles[disable_turbines] = 0.0
+            power_setpoints[disable_turbines] = 0.001 # Not zero to avoid numerical problems
+
         self.floris.farm.power_setpoints = power_setpoints
 
         # # TODO is this required?
@@ -179,6 +214,8 @@ def calculate_wake(
     def calculate_no_wake(
         self,
         yaw_angles: NDArrayFloat | list[float] | None = None,
+        power_setpoints: NDArrayFloat | list[float] | list[float, None] | None = None,
+        disable_turbines: NDArrayBool | list[bool] | None = None,
     ) -> None:
         """
         This function is similar to `calculate_wake()` except
@@ -201,6 +238,51 @@ def calculate_no_wake(
             )
         self.floris.farm.yaw_angles = yaw_angles
 
+        if power_setpoints is None:
+            power_setpoints = POWER_SETPOINT_DEFAULT * np.ones(
+                (
+                    self.floris.flow_field.n_findex,
+                    self.floris.farm.n_turbines,
+                )
+            )
+        else:
+            power_setpoints = np.array(power_setpoints)
+
+        # Convert any None values to the default power setpoint
+        power_setpoints[
+            power_setpoints == np.full(power_setpoints.shape, None)
+        ] = POWER_SETPOINT_DEFAULT
+        power_setpoints = floris_array_converter(power_setpoints)
+
+        # Check for turbines to disable
+        if disable_turbines is not None:
+
+            # Force to numpy array
+            # disable_turbines = np.array(disable_turbines)
+
+            # Must have first dimension = n_findex
+            if disable_turbines.shape[0] != self.floris.flow_field.n_findex:
+                raise ValueError(
+                    f"disable_turbines has a size of {disable_turbines.shape[0]} "
+                    f"in the 0th dimension, must be equal to "
+                    f"n_findex={self.floris.flow_field.n_findex}"
+                )
+
+            # Must have first dimension = n_turbines
+            if disable_turbines.shape[1] != self.floris.farm.n_turbines:
+                raise ValueError(
+                    f"disable_turbines has a size of {disable_turbines.shape[1]} "
+                    f"in the 1th dimension, must be equal to "
+                    f"n_turbines={self.floris.farm.n_turbines}"
+                )
+
+            # Set power_setpoints and yaw_angles to 0 in all locations where
+            # disable_turbines is True
+            yaw_angles[disable_turbines] = 0.0
+            power_setpoints[disable_turbines] = 0.001 # Not zero to avoid numerical problems
+
+        self.floris.farm.power_setpoints = power_setpoints
+
         # Initialize solution space
         self.floris.initialize_domain()
 

tests/floris_interface_test.py

@@ -2,6 +2,7 @@
 
 import numpy as np
 import pytest
+import yaml
 
 from floris.simulation.turbine.operation_models import POWER_SETPOINT_DEFAULT
 from floris.tools.floris_interface import FlorisInterface
@@ -15,8 +16,6 @@ def test_read_yaml():
     fi = FlorisInterface(configuration=YAML_INPUT)
     assert isinstance(fi, FlorisInterface)
 
-
-
 def test_calculate_wake():
     """
     In FLORIS v3.2, running calculate_wake twice incorrectly set the yaw angles when the first time
@@ -143,6 +142,79 @@ def test_get_farm_power():
     farm_power_from_turbine = turbine_powers.sum(axis=1)
     np.testing.assert_almost_equal(farm_power_from_turbine, farm_powers)
 
+def test_disable_turbines():
+
+    fi = FlorisInterface(configuration=YAML_INPUT)
+
+    # Set to mixed turbine model
+    with open(
+        str(
+            fi.floris.as_dict()["farm"]["turbine_library_path"]
+            / (fi.floris.as_dict()["farm"]["turbine_type"][0] + ".yaml")
+        )
+    ) as t:
+        turbine_type = yaml.safe_load(t)
+    turbine_type["power_thrust_model"] = "mixed"
+    fi.reinitialize(turbine_type=[turbine_type])
+
+    # Init to n-findex = 2, n_turbines = 3
+    fi.reinitialize(
+        wind_speeds=np.array([8.,8.,]),
+        wind_directions=np.array([270.,270.]),
+        layout_x = [0,1000,2000],
+        layout_y=[0,0,0]
+    )
+
+    # Confirm that passing in a disable value with wrong n_findex raises error
+    with pytest.raises(ValueError):
+        fi.calculate_wake(disable_turbines=np.zeros((10, 3), dtype=bool))
+
+    # Confirm that passing in a disable value with wrong n_turbines raises error
+    with pytest.raises(ValueError):
+        fi.calculate_wake(disable_turbines=np.zeros((2, 10), dtype=bool))
+
+    # Confirm that if all turbines are disabled, power is near 0 for all turbines
+    fi.calculate_wake(disable_turbines=np.ones((2, 3), dtype=bool))
+    turbines_powers = fi.get_turbine_powers()
+    np.testing.assert_allclose(turbines_powers,0,atol=0.1)
+
+    # Confirm the same for calculate_no_wake
+    fi.calculate_no_wake(disable_turbines=np.ones((2, 3), dtype=bool))
+    turbines_powers = fi.get_turbine_powers()
+    np.testing.assert_allclose(turbines_powers,0,atol=0.1)
+
+    # Confirm that if all disabled values set to false, equivalent to running normally
+    fi.calculate_wake()
+    turbines_powers_normal = fi.get_turbine_powers()
+    fi.calculate_wake(disable_turbines=np.zeros((2, 3), dtype=bool))
+    turbines_powers_false_disable = fi.get_turbine_powers()
+    np.testing.assert_allclose(turbines_powers_normal,turbines_powers_false_disable,atol=0.1)
+
+    # Confirm the same for calculate_no_wake
+    fi.calculate_no_wake()
+    turbines_powers_normal = fi.get_turbine_powers()
+    fi.calculate_no_wake(disable_turbines=np.zeros((2, 3), dtype=bool))
+    turbines_powers_false_disable = fi.get_turbine_powers()
+    np.testing.assert_allclose(turbines_powers_normal,turbines_powers_false_disable,atol=0.1)
+
+    # Confirm the shutting off the middle turbine is like removing from the layout
+    # In terms of impact on third turbine
+    disable_turbines = np.zeros((2, 3), dtype=bool)
+    disable_turbines[:,1] = [True, True]
+    fi.calculate_wake(disable_turbines=disable_turbines)
+    power_with_middle_disabled = fi.get_turbine_powers()
+
+    fi.reinitialize(layout_x = [0,2000],layout_y = [0, 0])
+    fi.calculate_wake()
+    power_with_middle_removed = fi.get_turbine_powers()
+
+    np.testing.assert_almost_equal(power_with_middle_disabled[0,2], power_with_middle_removed[0,1])
+    np.testing.assert_almost_equal(power_with_middle_disabled[1,2], power_with_middle_removed[1,1])
+
+    # Check that yaw angles are correctly set when turbines are disabled
+    fi.reinitialize(layout_x = [0,1000,2000],layout_y = [0,0,0])
+    fi.calculate_wake(disable_turbines=disable_turbines, yaw_angles=np.ones((2, 3)))
+    assert (fi.floris.farm.yaw_angles == np.array([[1.0, 0.0, 1.0], [1.0, 0.0, 1.0]])).all()
 
 def test_get_farm_aep():
     fi = FlorisInterface(configuration=YAML_INPUT)