Feature/load heuristic

examples/examples_control_types/001b_derating_for_load.py

@@ -0,0 +1,90 @@
+"""Example of de-rating for load
+
+TBD.
+"""
+
+import matplotlib.pyplot as plt
+import numpy as np
+
+from floris import FlorisModel
+
+
+fmodel = FlorisModel("../inputs/gch.yaml")
+
+# Change to the simple-derating model turbine
+# (Note this could also be done with the mixed model)
+fmodel.set_operation_model("simple-derating")
+
+# Convert to a simple two turbine layout with derating turbines
+fmodel.set(layout_x=[0, 126.0 * 5], layout_y=[0.0, 0.0])
+
+
+# Set the wind directions and speeds to be constant over n_findex = N time steps
+N = 50
+fmodel.set(
+    wind_directions=270 * np.ones(N),
+    wind_speeds=7.0 * np.ones(N),
+    turbulence_intensities=0.06 * np.ones(N),
+)
+fmodel.run()
+turbine_powers_orig = fmodel.get_turbine_powers()
+
+# Get the nominal operating power of the upstream turbine
+nom_pow = turbine_powers_orig[0,0]
+
+# # Derate the front turbine from full power to 75 %
+power_setpoints_front = np.linspace(nom_pow,nom_pow *.75,N)
+full_rating = np.ones_like(power_setpoints_front) * 5E6
+
+# Only derate the front turbine
+power_setpoints = np.column_stack([power_setpoints_front, full_rating])
+fmodel.set(power_setpoints=power_setpoints)
+fmodel.run()
+turbine_powers_derated = fmodel.get_turbine_powers()
+
+# Compute the mean power
+power_mean = np.mean(turbine_powers_derated,axis=1)
+
+# Define de-rating as percent reduction
+de_rating = 100 * (nom_pow - power_setpoints_front) / nom_pow
+
+# Grab the load heuristic
+load_h = fmodel._get_turbine_load_h()
+
+# Compute mean load
+load_h_mean = np.mean(load_h,axis=1)
+
+# Plot the results
+fig, axarr = plt.subplots(2,3,sharex=True, sharey='row',figsize=(12,8))
+
+ax = axarr[0,0]
+ax.plot(de_rating, turbine_powers_derated[:,0]/turbine_powers_derated[0,0],'k')
+ax.set_title('Turbine 0')
+ax.set_ylabel('Power (ratio to nomimal)')
+
+ax = axarr[0,1]
+ax.plot(de_rating, turbine_powers_derated[:,1]/turbine_powers_derated[0,1],'k')
+ax.set_title('Turbine 1')
+
+ax = axarr[0,2]
+ax.plot(de_rating, power_mean/power_mean[0],'k')
+ax.set_title('Mean')
+
+ax = axarr[1,0]
+ax.plot(de_rating, load_h[:,0]/load_h[0,0],'k')
+ax.set_ylabel('Load H (ratio to nomimal)')
+ax.set_xlabel('Front De-Rating %')
+
+ax = axarr[1,1]
+ax.plot(de_rating, load_h[:,1]/load_h[0,1],'k')
+ax.set_xlabel('Front De-Rating %')
+
+ax = axarr[1,2]
+ax.plot(de_rating, load_h_mean/load_h_mean[0],'k')
+ax.set_xlabel('Front De-Rating %')
+
+for ax in axarr.flatten():
+    ax.grid(True)
+
+
+plt.show()

floris/floris_model.py

@@ -36,6 +36,7 @@
     nested_get,
     nested_set,
     print_nested_dict,
+    rotate_coordinates_rel_west,
 )
 from floris.wind_data import (
     TimeSeries,
@@ -503,8 +504,112 @@ def run_no_wake(self) -> None:
         self.core.finalize()
 
 
-    ### Methods for extracting turbine performance after running
+    #TODO: This could probably be moved elsewhere
+    #TODO: load_h is a stand-in term for "load heuristic" until something better comes to mind
+    def _get_turbine_load_h(self) -> NDArrayFloat:
+
+        #TODO: Probably these should be function input parameters
+        wake_slope = 0.3
+        max_dist_D = 10.0
+        weight_ws = 1.0
+        weight_ct = 1.0
+        weight_ti = 1.0
+
+        # This one I'm not sure, maybe it will be ambient ti, maybe a 1D TD load_ambient_ti
+        load_ambient_ti = 0.10
+
+        # Store for convenience
+        n_turbines = self.n_turbines
+        n_findex = self.n_findex
+
+        # TODO for now assume the first one
+        D = self.core.farm.rotor_diameters[0,0]
+
+        # Get the wind direction rotated farm layouts with wind flowing left to right
+        coord_rotate = rotate_coordinates_rel_west(self.wind_directions,
+                                                    self.core.grid.turbine_coordinates)
+        x_rot = coord_rotate[0]
+        y_rot = coord_rotate[1]
+
+        # Get all turbine thrust coefficients
+        cts = self.get_turbine_thrust_coefficients()
+
+        # Get ws, was originally thinking to use the actual rotor wind speeds (commented first line)
+        # but think more consistent with load methods to use ambient here since those ambient flows
+        # will surround the wake somehow for load inducing purposes
+        # ws = np.mean(self.core.flow_field.u, axis=(2,3))
+        ws = self.wind_speeds
+        # Make the same for every turbine
+        ws = np.tile(ws[:, np.newaxis], (1, n_turbines))
+
+        # Initialize the load_ti to the load_ambient_ti
+        load_ti = np.ones((n_findex, n_turbines)) * load_ambient_ti
+
+        # For each findex and turbine, compute the "load turbulence"
+        for findex in range(n_findex):
+
+            # t_i is the index of the "current" turbine whose "load_ti" we want to compute
+            for t_i in range(n_turbines):
+                nearest_t_dist = 1 + np.max(x_rot) - np.min(x_rot) # Initialize to the max + 1
+
+                # t_j is the index of the turbine who is proposed to be waking t_i
+                for t_j in range(n_turbines):
+
+                    # If the same turbine skip
+                    if t_i == t_j:
+                        continue
+
+                    # Compute the distance from proposed waking to turbine to the current turbine
+                    # with the distance positive if the current turbine is downstream
+                    x_diff = x_rot[findex, t_i] - x_rot[findex, t_j]
+
+                    # If the proposed waking turbine is downstream skip
+                    if x_diff<=0:
+                        continue
+
+                    # If the distance between the turbines is more than the max skip
+                    if x_diff >= max_dist_D * D:
+                        continue
 
+                    # If this is not closer than previouly found wakes skip
+                    if x_diff >= nearest_t_dist:
+                        continue
+
+                    # Compute the cross-stream distance
+                    y_diff_abs = np.abs(y_rot[findex, t_i] - y_rot[findex, t_j])
+
+                    # If the t_i turbine is above the wake cone of the
+                    # proposed waking turbine skip
+                    if y_diff_abs > D + x_diff * wake_slope:
+                        continue
+
+                    # Set nearest_t_dist to the current value
+                    nearest_t_dist = x_diff
+
+                    # Get the c_t of the proposed waking turbine
+                    c_t_j = cts[findex, t_j]
+
+                    # Get the ambient ws (note all turbines same value
+                    #  but might as well use the current turbine)
+                    ambient_ws = ws[findex,t_i]
+
+                    # Compute the added turbulence
+                    ti_add = ambient_ws / (1.5 + 0.8 * (x_diff/D) / np.sqrt(c_t_j) )
+
+                    # Update the ti_term
+                    load_ti[findex, t_i] = (
+
+                        np.sqrt(ti_add**2 + (load_ambient_ti * ambient_ws)**2) / ambient_ws
+
+                    )
+
+
+
+        # Compute load_h
+        return weight_ws * ws + weight_ct * cts + weight_ti * load_ti
+
+
+    ### Methods for extracting turbine performance after running
     def _get_turbine_powers(self) -> NDArrayFloat:
         """Calculates the power at each turbine in the wind farm.
 

tests/floris_model_integration_test.py

@@ -193,6 +193,37 @@ def test_run_no_wake():
     power_no_wake = fmodel.get_turbine_powers()
     assert len(np.unique(power_no_wake)) == 1
 
+def test_get_turbine_load_h():
+
+    fmodel = FlorisModel(configuration=YAML_INPUT)
+
+    wind_speeds = np.array([8.0, 8.0, 8.0])
+    wind_directions = np.array([270.0, 270.0, 280.0])
+    turbulence_intensities = np.array([0.06, 0.06, 0.06])
+    n_findex = len(wind_directions)
+
+    layout_x = np.array([0, 1000])
+    layout_y = np.array([0, 1])
+    n_turbines = len(layout_x)
+
+    fmodel.set(
+        wind_speeds=wind_speeds,
+        wind_directions=wind_directions,
+        turbulence_intensities=turbulence_intensities,
+        layout_x=layout_x,
+        layout_y=layout_y,
+    )
+
+    fmodel.run()
+
+    load_h = fmodel._get_turbine_load_h()
+
+    assert load_h.shape[0] == n_findex
+    assert load_h.shape[1] == n_turbines
+    assert load_h[0, 0] == load_h[1, 0]
+
+
+
 def test_get_turbine_powers():
     # Get turbine powers should return n_findex x n_turbine powers
     # Apply the same wind speed and direction multiple times and confirm all equal