Add de-rating op model

examples/40_test_derating.py

@@ -0,0 +1,116 @@
+# Copyright 2024 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
+
+import matplotlib.pyplot as plt
+import numpy as np
+import yaml
+
+from floris.tools import FlorisInterface
+
+
+"""
+Example to test out derating of turbines and mixed derating and yawing. Will be refined before
+release. TODO: Demonstrate shutting off turbines also, once developed.
+"""
+
+# Grab model of FLORIS and update to deratable turbines
+fi = FlorisInterface("inputs/gch.yaml")
+
+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"] = "simple-derating"
+
+# Convert to a simple two turbine layout with derating turbines
+fi.reinitialize(layout_x=[0, 1000.0], layout_y=[0.0, 0.0], turbine_type=[turbine_type])
+
+# Set the wind directions and speeds to be constant over n_findex = N time steps
+N = 50
+fi.reinitialize(wind_directions=270 * np.ones(N), wind_speeds=10.0 * np.ones(N))
+fi.calculate_wake()
+turbine_powers_orig = fi.get_turbine_powers()
+
+# Add derating
+power_setpoints = np.tile(np.linspace(1, 6e6, N), 2).reshape(2, N).T
+fi.calculate_wake(power_setpoints=power_setpoints)
+turbine_powers_derated = fi.get_turbine_powers()
+
+# Compute available power at downstream turbine
+power_setpoints_2 = np.array([np.linspace(1, 6e6, N), np.full(N, None)]).T
+fi.calculate_wake(power_setpoints=power_setpoints_2)
+turbine_powers_avail_ds = fi.get_turbine_powers()[:,1]
+
+# Plot the results
+fig, ax = plt.subplots(1, 1)
+ax.plot(power_setpoints[:, 0]/1000, turbine_powers_derated[:, 0]/1000, color="C0", label="Upstream")
+ax.plot(
+    power_setpoints[:, 1]/1000,
+    turbine_powers_derated[:, 1]/1000,
+    color="C1",
+    label="Downstream"
+)
+ax.plot(
+    power_setpoints[:, 0]/1000,
+    turbine_powers_orig[:, 0]/1000,
+    color="C0",
+    linestyle="dotted",
+    label="Upstream available"
+)
+ax.plot(
+    power_setpoints[:, 1]/1000,
+    turbine_powers_avail_ds/1000,
+    color="C1",
+    linestyle="dotted", label="Downstream available"
+)
+ax.plot(
+    power_setpoints[:, 1]/1000,
+    np.ones(N)*np.max(turbine_type["power_thrust_table"]["power"]),
+    color="k",
+    linestyle="dashed",
+    label="Rated power"
+)
+ax.grid()
+ax.legend()
+ax.set_xlim([0, 6e3])
+ax.set_xlabel("Power setpoint (kW)")
+ax.set_ylabel("Power produced (kW)")
+
+# Second example showing mixed model use.
+turbine_type["power_thrust_model"] = "mixed"
+yaw_angles = np.array([
+    [0.0, 0.0],
+    [0.0, 0.0],
+    [20.0, 10.0],
+    [0.0, 10.0],
+    [20.0, 0.0]
+])
+power_setpoints = np.array([
+    [None, None],
+    [2e6, 1e6],
+    [None, None],
+    [2e6, None,],
+    [None, 1e6]
+])
+fi.reinitialize(
+    wind_directions=270 * np.ones(len(yaw_angles)),
+    wind_speeds=10.0 * np.ones(len(yaw_angles)),
+    turbine_type=[turbine_type]*2
+)
+fi.calculate_wake(yaw_angles=yaw_angles, power_setpoints=power_setpoints)
+turbine_powers = fi.get_turbine_powers()
+print(turbine_powers)
+
+plt.show()

floris/simulation/farm.py

@@ -32,6 +32,7 @@
     Turbine,
 )
 from floris.simulation.rotor_velocity import compute_tilt_angles_for_floating_turbines_map
+from floris.simulation.turbine.operation_models import POWER_SETPOINT_DEFAULT
 from floris.type_dec import (
     convert_to_path,
     floris_array_converter,
@@ -92,6 +93,9 @@ class Farm(BaseClass):
     tilt_angles: NDArrayFloat = field(init=False)
     tilt_angles_sorted: NDArrayFloat = field(init=False)
 
+    power_setpoints: NDArrayFloat = field(init=False)
+    power_setpoints_sorted: NDArrayFloat = field(init=False)
+
     hub_heights: NDArrayFloat = field(init=False)
     hub_heights_sorted: NDArrayFloat = field(init=False, factory=list)
 
@@ -233,6 +237,11 @@ def initialize(self, sorted_indices):
             sorted_indices[:, :, 0, 0],
             axis=1,
         )
+        self.power_setpoints_sorted = np.take_along_axis(
+            self.power_setpoints,
+            sorted_indices[:, :, 0, 0],
+            axis=1,
+        )
         self.state = State.INITIALIZED
 
     def construct_hub_heights(self):
@@ -341,6 +350,10 @@ def set_tilt_to_ref_tilt(self, n_findex: int):
             * self.ref_tilts
         )
 
+    def set_power_setpoints(self, n_findex: int):
+        self.power_setpoints = POWER_SETPOINT_DEFAULT * np.ones((n_findex, self.n_turbines))
+        self.power_setpoints_sorted = POWER_SETPOINT_DEFAULT * np.ones((n_findex, self.n_turbines))
+
     def calculate_tilt_for_eff_velocities(self, rotor_effective_velocities):
         tilt_angles = compute_tilt_angles_for_floating_turbines_map(
             self.turbine_type_map_sorted,

floris/simulation/floris.py

@@ -98,6 +98,7 @@ def __attrs_post_init__(self) -> None:
         self.farm.construct_turbine_correct_cp_ct_for_tilt()
         self.farm.set_yaw_angles(self.flow_field.n_findex)
         self.farm.set_tilt_to_ref_tilt(self.flow_field.n_findex)
+        self.farm.set_power_setpoints(self.flow_field.n_findex)
 
         if self.solver["type"] == "turbine_grid":
             self.grid = TurbineGrid(

floris/simulation/solver.py

@@ -101,8 +101,10 @@ def sequential_solver(
 
         ct_i = thrust_coefficient(
             velocities=flow_field.u_sorted,
+            air_density=flow_field.air_density,
             yaw_angles=farm.yaw_angles_sorted,
             tilt_angles=farm.tilt_angles_sorted,
+            power_setpoints=farm.power_setpoints_sorted,
             thrust_coefficient_functions=farm.turbine_thrust_coefficient_functions,
             tilt_interps=farm.turbine_tilt_interps,
             correct_cp_ct_for_tilt=farm.correct_cp_ct_for_tilt_sorted,
@@ -118,8 +120,10 @@ def sequential_solver(
         ct_i = ct_i[:, 0:1, None, None]
         axial_induction_i = axial_induction(
             velocities=flow_field.u_sorted,
+            air_density=flow_field.air_density,
             yaw_angles=farm.yaw_angles_sorted,
             tilt_angles=farm.tilt_angles_sorted,
+            power_setpoints=farm.power_setpoints_sorted,
             axial_induction_functions=farm.turbine_axial_induction_functions,
             tilt_interps=farm.turbine_tilt_interps,
             correct_cp_ct_for_tilt=farm.correct_cp_ct_for_tilt_sorted,
@@ -330,8 +334,10 @@ def full_flow_sequential_solver(
 
         ct_i = thrust_coefficient(
             velocities=turbine_grid_flow_field.u_sorted,
+            air_density=turbine_grid_flow_field.air_density,
             yaw_angles=turbine_grid_farm.yaw_angles_sorted,
             tilt_angles=turbine_grid_farm.tilt_angles_sorted,
+            power_setpoints=turbine_grid_farm.power_setpoints_sorted,
             thrust_coefficient_functions=turbine_grid_farm.turbine_thrust_coefficient_functions,
             tilt_interps=turbine_grid_farm.turbine_tilt_interps,
             correct_cp_ct_for_tilt=turbine_grid_farm.correct_cp_ct_for_tilt_sorted,
@@ -344,8 +350,10 @@ def full_flow_sequential_solver(
         ct_i = ct_i[:, 0:1, None, None]
         axial_induction_i = axial_induction(
             velocities=turbine_grid_flow_field.u_sorted,
+            air_density=turbine_grid_flow_field.air_density,
             yaw_angles=turbine_grid_farm.yaw_angles_sorted,
             tilt_angles=turbine_grid_farm.tilt_angles_sorted,
+            power_setpoints=turbine_grid_farm.power_setpoints_sorted,
             axial_induction_functions=turbine_grid_farm.turbine_axial_induction_functions,
             tilt_interps=turbine_grid_farm.turbine_tilt_interps,
             correct_cp_ct_for_tilt=turbine_grid_farm.correct_cp_ct_for_tilt_sorted,
@@ -495,8 +503,10 @@ def cc_solver(
         turb_avg_vels = average_velocity(turb_inflow_field)
         turb_Cts = thrust_coefficient(
             turb_avg_vels,
+            flow_field.air_density,
             farm.yaw_angles_sorted,
             farm.tilt_angles_sorted,
+            farm.power_setpoints_sorted,
             farm.turbine_thrust_coefficient_functions,
             tilt_interps=farm.turbine_tilt_interps,
             correct_cp_ct_for_tilt=farm.correct_cp_ct_for_tilt_sorted,
@@ -508,8 +518,10 @@ def cc_solver(
         turb_Cts = turb_Cts[:, :, None, None]
         turb_aIs = axial_induction(
             turb_avg_vels,
+            flow_field.air_density,
             farm.yaw_angles_sorted,
             farm.tilt_angles_sorted,
+            farm.power_setpoints_sorted,
             farm.turbine_axial_induction_functions,
             tilt_interps=farm.turbine_tilt_interps,
             correct_cp_ct_for_tilt=farm.correct_cp_ct_for_tilt_sorted,
@@ -526,8 +538,10 @@ def cc_solver(
 
         axial_induction_i = axial_induction(
             velocities=flow_field.u_sorted,
+            air_density=flow_field.air_density,
             yaw_angles=farm.yaw_angles_sorted,
             tilt_angles=farm.tilt_angles_sorted,
+            power_setpoints=farm.power_setpoints_sorted,
             axial_induction_functions=farm.turbine_axial_induction_functions,
             tilt_interps=farm.turbine_tilt_interps,
             correct_cp_ct_for_tilt=farm.correct_cp_ct_for_tilt_sorted,
@@ -737,8 +751,10 @@ def full_flow_cc_solver(
         turb_avg_vels = average_velocity(turbine_grid_flow_field.u_sorted)
         turb_Cts = thrust_coefficient(
             velocities=turb_avg_vels,
+            air_density=flow_field_grid.air_density,
             yaw_angles=turbine_grid_farm.yaw_angles_sorted,
             tilt_angles=turbine_grid_farm.tilt_angles_sorted,
+            power_setpoints=turbine_grid_farm.power_setpoints_sorted,
             thrust_coefficient_functions=turbine_grid_farm.turbine_thrust_coefficient_functions,
             tilt_interps=turbine_grid_farm.turbine_tilt_interps,
             correct_cp_ct_for_tilt=turbine_grid_farm.correct_cp_ct_for_tilt_sorted,
@@ -751,8 +767,10 @@ def full_flow_cc_solver(
 
         axial_induction_i = axial_induction(
             velocities=turbine_grid_flow_field.u_sorted,
+            air_density=turbine_grid_flow_field.air_density,
             yaw_angles=turbine_grid_farm.yaw_angles_sorted,
             tilt_angles=turbine_grid_farm.tilt_angles_sorted,
+            power_setpoints=turbine_grid_farm.power_setpoints_sorted,
             axial_induction_functions=turbine_grid_farm.turbine_axial_induction_functions,
             tilt_interps=turbine_grid_farm.turbine_tilt_interps,
             correct_cp_ct_for_tilt=turbine_grid_farm.correct_cp_ct_for_tilt_sorted,
@@ -891,8 +909,10 @@ def turbopark_solver(
 
         Cts = thrust_coefficient(
             velocities=flow_field.u_sorted,
+            air_density=flow_field.air_density,
             yaw_angles=farm.yaw_angles_sorted,
             tilt_angles=farm.tilt_angles_sorted,
+            power_setpoints=farm.power_setpoints_sorted,
             thrust_coefficient_functions=farm.turbine_thrust_coefficient_functions,
             tilt_interps=farm.turbine_tilt_interps,
             correct_cp_ct_for_tilt=farm.correct_cp_ct_for_tilt_sorted,
@@ -904,8 +924,10 @@ def turbopark_solver(
 
         ct_i = thrust_coefficient(
             velocities=flow_field.u_sorted,
+            air_density=flow_field.air_density,
             yaw_angles=farm.yaw_angles_sorted,
             tilt_angles=farm.tilt_angles_sorted,
+            power_setpoints=farm.power_setpoints_sorted,
             thrust_coefficient_functions=farm.turbine_thrust_coefficient_functions,
             tilt_interps=farm.turbine_tilt_interps,
             correct_cp_ct_for_tilt=farm.correct_cp_ct_for_tilt_sorted,
@@ -920,8 +942,10 @@ def turbopark_solver(
         ct_i = ct_i[:, 0:1, None, None]
         axial_induction_i = axial_induction(
             velocities=flow_field.u_sorted,
+            air_density=flow_field.air_density,
             yaw_angles=farm.yaw_angles_sorted,
             tilt_angles=farm.tilt_angles_sorted,
+            power_setpoints=farm.power_setpoints_sorted,
             axial_induction_functions=farm.turbine_axial_induction_functions,
             tilt_interps=farm.turbine_tilt_interps,
             correct_cp_ct_for_tilt=farm.correct_cp_ct_for_tilt_sorted,
@@ -978,8 +1002,10 @@ def turbopark_solver(
                 turbulence_intensity_ii = turbine_turbulence_intensity[:, ii:ii+1]
                 ct_ii = thrust_coefficient(
                     velocities=flow_field.u_sorted,
+                    air_density=flow_field.air_density,
                     yaw_angles=farm.yaw_angles_sorted,
                     tilt_angles=farm.tilt_angles_sorted,
+                    power_setpoints=farm.power_setpoints_sorted,
                     thrust_coefficient_functions=farm.turbine_thrust_coefficient_functions,
                     tilt_interps=farm.turbine_tilt_interps,
                     correct_cp_ct_for_tilt=farm.correct_cp_ct_for_tilt_sorted,
@@ -1174,8 +1200,10 @@ def empirical_gauss_solver(
 
         ct_i = thrust_coefficient(
             velocities=flow_field.u_sorted,
+            air_density=flow_field.air_density,
             yaw_angles=farm.yaw_angles_sorted,
             tilt_angles=farm.tilt_angles_sorted,
+            power_setpoints=farm.power_setpoints_sorted,
             thrust_coefficient_functions=farm.turbine_thrust_coefficient_functions,
             tilt_interps=farm.turbine_tilt_interps,
             correct_cp_ct_for_tilt=farm.correct_cp_ct_for_tilt_sorted,
@@ -1190,8 +1218,10 @@ def empirical_gauss_solver(
         ct_i = ct_i[:, 0:1, None, None]
         axial_induction_i = axial_induction(
             velocities=flow_field.u_sorted,
+            air_density=flow_field.air_density,
             yaw_angles=farm.yaw_angles_sorted,
             tilt_angles=farm.tilt_angles_sorted,
+            power_setpoints=farm.power_setpoints_sorted,
             axial_induction_functions=farm.turbine_axial_induction_functions,
             tilt_interps=farm.turbine_tilt_interps,
             correct_cp_ct_for_tilt=farm.correct_cp_ct_for_tilt_sorted,
@@ -1375,8 +1405,10 @@ def full_flow_empirical_gauss_solver(
 
         ct_i = thrust_coefficient(
             velocities=turbine_grid_flow_field.u_sorted,
+            air_density=turbine_grid_flow_field.air_density,
             yaw_angles=turbine_grid_farm.yaw_angles_sorted,
             tilt_angles=turbine_grid_farm.tilt_angles_sorted,
+            power_setpoints=turbine_grid_farm.power_setpoints_sorted,
             thrust_coefficient_functions=turbine_grid_farm.turbine_thrust_coefficient_functions,
             tilt_interps=turbine_grid_farm.turbine_tilt_interps,
             correct_cp_ct_for_tilt=turbine_grid_farm.correct_cp_ct_for_tilt_sorted,
@@ -1389,8 +1421,10 @@ def full_flow_empirical_gauss_solver(
         ct_i = ct_i[:, 0:1, None, None]
         axial_induction_i = axial_induction(
             velocities=turbine_grid_flow_field.u_sorted,
+            air_density=turbine_grid_flow_field.air_density,
             yaw_angles=turbine_grid_farm.yaw_angles_sorted,
             tilt_angles=turbine_grid_farm.tilt_angles_sorted,
+            power_setpoints=turbine_grid_farm.power_setpoints_sorted,
             axial_induction_functions=turbine_grid_farm.turbine_axial_induction_functions,
             tilt_interps=turbine_grid_farm.turbine_tilt_interps,
             correct_cp_ct_for_tilt=turbine_grid_farm.correct_cp_ct_for_tilt_sorted,

floris/simulation/turbine/__init__.py

@@ -14,5 +14,7 @@
 
 from floris.simulation.turbine.operation_models import (
     CosineLossTurbine,
+    MixedOperationTurbine,
+    SimpleDeratingTurbine,
     SimpleTurbine,
 )