Smooth output of FlorisStandin to prevent oscilations in closed-loop response

hercules/floris_standin.py

@@ -27,7 +27,7 @@
 from pathlib import Path
 
 import numpy as np
-from floris.tools import FlorisInterface
+from floris import FlorisModel
 from floris.turbine_library import build_cosine_loss_turbine_dict
 
 from hercules.amr_wind_standin import AMRWindStandin, read_amr_wind_input
@@ -119,16 +119,36 @@ def construct_floris_from_amr_input(amr_wind_input):
         turb_dict["power_thrust_model"] = "mixed"
 
         # load a default model
-        fi = FlorisInterface(default_floris_dict)
-        fi.set(
+        fmodel = FlorisModel(default_floris_dict)
+        fmodel.set(
             layout_x=layout_x, layout_y=layout_y, turbine_type=[turb_dict] * len(layout_x)
         )
 
-    return fi
+    return fmodel
 
 
 class FlorisStandin(AMRWindStandin):
-    def __init__(self, config_dict, amr_input_file, amr_standin_data_file=None):
+    """
+    FlorisStandin class, which stands in for AMR-Wind. 
+    Arguments:
+    config_dict: dictionary of configuration parameters
+    amr_input_file: path to the AMR-Wind input file
+    amr_standin_data_file [optional]: path to the AMR-Wind standin data file.
+        Defaults to None
+    smoothing_coefficient [optional]: smoothing coefficient for turbine power
+        output. Must be in [0, 1). If 0, no smoothing is applied; if near 1,
+        the output is heavily smoothed. Defaults to 0.5.
+    """
+    def __init__(
+            self,
+            config_dict,
+            amr_input_file,
+            amr_standin_data_file=None,
+            smoothing_coefficient=0.5
+        ):
+        """
+        Constructor for the FlorisStandin class
+        """
         # Ensure outputs folder exists
         Path("outputs").mkdir(parents=True, exist_ok=True)
 
@@ -139,16 +159,22 @@ def __init__(self, config_dict, amr_input_file, amr_standin_data_file=None):
         )
 
         # Construct the floris object
-        self.fi = construct_floris_from_amr_input(amr_input_file)
+        self.fmodel = construct_floris_from_amr_input(amr_input_file)
 
         # Get the number of turbines
-        self.num_turbines = len(self.fi.layout_x)
+        self.num_turbines = len(self.fmodel.layout_x)
 
         # Print the number of turbines
         logger.info("Number of turbines: {}".format(self.num_turbines))
 
         # Initialize storage
         self.yaw_angles_stored = [0.0] * self.num_turbines
+        self.turbine_powers_prev = np.zeros(self.num_turbines)
+
+        # Check and save smoothing coefficient
+        if smoothing_coefficient < 0 or smoothing_coefficient >= 1:
+            raise ValueError("Smoothing coefficient must be in [0, 1).")
+        self.smoothing_coefficient = smoothing_coefficient
 
     def get_step(self, sim_time_s, yaw_angles=None, power_setpoints=None):
         """Retreive or calculate wind speed, direction, and turbine powers
@@ -211,14 +237,22 @@ def get_step(self, sim_time_s, yaw_angles=None, power_setpoints=None):
             power_setpoints = power_setpoints * 1000 # in W
 
         # Set up and solve FLORIS
-        self.fi.set(
+        self.fmodel.set(
             wind_speeds=[amr_wind_speed],
             wind_directions=[amr_wind_direction],
             yaw_angles=yaw_misalignments,
             power_setpoints=power_setpoints
         )
-        self.fi.run()
-        turbine_powers = (self.fi.get_turbine_powers() / 1000).flatten().tolist()  # in kW
+        self.fmodel.run()
+        turbine_powers_floris = (self.fmodel.get_turbine_powers() / 1000).flatten()  # in kW
+
+        # Smooth output
+        turbine_powers = (
+            self.smoothing_coefficient*self.turbine_powers_prev
+            + (1-self.smoothing_coefficient)*turbine_powers_floris
+        )
+        self.turbine_powers_prev = turbine_powers
+        turbine_powers = turbine_powers.tolist()
 
         return (
             amr_wind_speed,

tests/floris_standin_test.py

@@ -2,7 +2,7 @@
 
 import numpy as np
 import pytest
-from floris.tools import FlorisInterface
+from floris import FlorisModel
 from hercules.amr_wind_standin import AMRWindStandin
 from hercules.floris_standin import (
     construct_floris_from_amr_input,
@@ -33,8 +33,8 @@
 
 
 def test_construct_floris_from_amr_input():
-    fi_test = construct_floris_from_amr_input(AMR_INPUT)
-    assert isinstance(fi_test, FlorisInterface)
+    fmodel_test = construct_floris_from_amr_input(AMR_INPUT)
+    assert isinstance(fmodel_test, FlorisModel)
 
 
 def test_FlorisStandin_instantiation():
@@ -46,114 +46,114 @@ def test_FlorisStandin_instantiation():
     assert isinstance(floris_standin, AMRWindStandin)
 
     # Get FLORIS equivalent, match layout and turbines
-    fi_true = FlorisInterface(default_floris_dict)
-    fi_true.set(
-        layout_x=floris_standin.fi.layout_x,
-        layout_y=floris_standin.fi.layout_y,
-        turbine_type=floris_standin.fi.floris.farm.turbine_definitions,
+    fmodel_true = FlorisModel(default_floris_dict)
+    fmodel_true.set(
+        layout_x=floris_standin.fmodel.layout_x,
+        layout_y=floris_standin.fmodel.layout_y,
+        turbine_type=floris_standin.fmodel.core.farm.turbine_definitions,
     )
 
-    assert fi_true.floris.as_dict() == floris_standin.fi.floris.as_dict()
+    assert fmodel_true.core.as_dict() == floris_standin.fmodel.core.as_dict()
 
 
 def test_FlorisStandin_get_step_yaw_angles():
-    floris_standin = FlorisStandin(CONFIG, AMR_INPUT)
+    floris_standin = FlorisStandin(CONFIG, AMR_INPUT, smoothing_coefficient=0.0)
 
     # Get FLORIS equivalent, match layout and turbines
-    fi_true = FlorisInterface(default_floris_dict)
-    fi_true.set(
-        layout_x=floris_standin.fi.layout_x,
-        layout_y=floris_standin.fi.layout_y,
-        turbine_type=floris_standin.fi.floris.farm.turbine_definitions,
+    fmodel_true = FlorisModel(default_floris_dict)
+    fmodel_true.set(
+        layout_x=floris_standin.fmodel.layout_x,
+        layout_y=floris_standin.fmodel.layout_y,
+        turbine_type=floris_standin.fmodel.core.farm.turbine_definitions,
     )
 
     default_wind_direction = 240.0  # Matches default in FlorisStandin
     default_wind_speed = 8.0  # Matches default in FlorisStandin
 
     # Test with None yaw angles
     fs_ws, fs_wd, fs_tp, fs_twd = floris_standin.get_step(5.0)
-    fi_true.set(wind_speeds=[default_wind_speed], wind_directions=[default_wind_direction])
-    fi_true.run()
-    fi_true_tp = fi_true.get_turbine_powers() / 1000 # kW expected
+    fmodel_true.set(wind_speeds=[default_wind_speed], wind_directions=[default_wind_direction])
+    fmodel_true.run()
+    fmodel_true_tp = fmodel_true.get_turbine_powers() / 1000 # kW expected
 
     assert fs_ws == default_wind_speed
     assert fs_wd == default_wind_direction
     assert fs_twd == [default_wind_direction] * 2
-    assert np.allclose(fs_tp, fi_true_tp.flatten().tolist())
+    assert np.allclose(fs_tp, fmodel_true_tp.flatten().tolist())
 
     # Test with any "no value" yaw angles (should apply no yaw angle)
     fs_ws, fs_wd, fs_tp, fs_twd = floris_standin.get_step(5.0, yaw_angles=[-1000, 20])
 
     assert fs_ws == default_wind_speed
     assert fs_wd == default_wind_direction
     assert fs_twd == [default_wind_direction] * 2
-    assert np.allclose(fs_tp, fi_true_tp.flatten().tolist())
+    assert np.allclose(fs_tp, fmodel_true_tp.flatten().tolist())
 
     # Test with aligned turbines
     yaw_angles = [240.0, 240.0]
     fs_ws, fs_wd, fs_tp, fs_twd = floris_standin.get_step(5.0, yaw_angles)
-    fi_true.set(wind_speeds=[default_wind_speed], wind_directions=[default_wind_direction])
-    fi_true.run()
-    fi_true_tp = fi_true.get_turbine_powers() / 1000 # kW expected
+    fmodel_true.set(wind_speeds=[default_wind_speed], wind_directions=[default_wind_direction])
+    fmodel_true.run()
+    fmodel_true_tp = fmodel_true.get_turbine_powers() / 1000 # kW expected
 
-    assert np.allclose(fs_tp, fi_true_tp.flatten().tolist())
+    assert np.allclose(fs_tp, fmodel_true_tp.flatten().tolist())
 
     # Test with misaligned turbines
     yaw_angles = [260.0, 230.0]
     fs_ws, fs_wd, fs_tp, fs_twd = floris_standin.get_step(5.0, yaw_angles)
-    fi_true.set(wind_speeds=[default_wind_speed], wind_directions=[default_wind_direction])
-    fi_true.run()  # Don't expect to work
-    fi_true_tp = fi_true.get_turbine_powers() / 1000
-    assert not np.allclose(fs_tp, fi_true_tp.flatten().tolist())
+    fmodel_true.set(wind_speeds=[default_wind_speed], wind_directions=[default_wind_direction])
+    fmodel_true.run()  # Don't expect to work
+    fmodel_true_tp = fmodel_true.get_turbine_powers() / 1000
+    assert not np.allclose(fs_tp, fmodel_true_tp.flatten().tolist())
 
     # Correct yaw angles
-    fi_true.set(yaw_angles=default_wind_direction - np.array([yaw_angles]))
-    fi_true.run()
-    fi_true_tp = fi_true.get_turbine_powers() / 1000 # kW expected
-    assert np.allclose(fs_tp, fi_true_tp.flatten().tolist())
+    fmodel_true.set(yaw_angles=default_wind_direction - np.array([yaw_angles]))
+    fmodel_true.run()
+    fmodel_true_tp = fmodel_true.get_turbine_powers() / 1000 # kW expected
+    assert np.allclose(fs_tp, fmodel_true_tp.flatten().tolist())
 
     # Test that yaw angles are maintained from the previous step if large misalignments are provided
     yaw_angles = [0.0, 10.0]
     _, _, fs_tp2, _ = floris_standin.get_step(5.0, yaw_angles)
     assert np.allclose(fs_tp, fs_tp2)
     assert np.allclose(
-        default_wind_direction-floris_standin.fi.floris.farm.yaw_angles,
+        default_wind_direction-floris_standin.fmodel.core.farm.yaw_angles,
         [260.0, 230.0]
     )
 
 def test_FlorisStandin_get_step_power_setpoints():
-    floris_standin = FlorisStandin(CONFIG, AMR_INPUT)
+    floris_standin = FlorisStandin(CONFIG, AMR_INPUT, smoothing_coefficient=0.0)
 
     # Get FLORIS equivalent, match layout and turbines
-    fi_true = FlorisInterface(default_floris_dict)
-    fi_true.set(
-        layout_x=floris_standin.fi.layout_x,
-        layout_y=floris_standin.fi.layout_y,
-        turbine_type=floris_standin.fi.floris.farm.turbine_definitions,
+    fmodel_true = FlorisModel(default_floris_dict)
+    fmodel_true.set(
+        layout_x=floris_standin.fmodel.layout_x,
+        layout_y=floris_standin.fmodel.layout_y,
+        turbine_type=floris_standin.fmodel.core.farm.turbine_definitions,
     )
 
     default_wind_direction = 240.0  # Matches default in FlorisStandin
     default_wind_speed = 8.0  # Matches default in FlorisStandin
 
     # Test with power setpoints
     fs_ws, fs_wd, fs_tp, fs_twd = floris_standin.get_step(5.0, power_setpoints=[1e3, 1e3])
-    fi_true.set(wind_speeds=[default_wind_speed], wind_directions=[default_wind_direction])
-    fi_true.run() # don't expect to work
-    fi_true_tp = fi_true.get_turbine_powers() / 1000
-    assert not np.allclose(fs_tp, fi_true_tp.flatten().tolist())
+    fmodel_true.set(wind_speeds=[default_wind_speed], wind_directions=[default_wind_direction])
+    fmodel_true.run() # don't expect to work
+    fmodel_true_tp = fmodel_true.get_turbine_powers() / 1000
+    assert not np.allclose(fs_tp, fmodel_true_tp.flatten().tolist())
 
     # Correct power setpoints
-    fi_true.set(power_setpoints=np.array([[1e6, 1e6]]))
-    fi_true.run()
-    fi_true_tp = fi_true.get_turbine_powers() / 1000 # kW expected
-    assert np.allclose(fs_tp, fi_true_tp.flatten().tolist())
+    fmodel_true.set(power_setpoints=np.array([[1e6, 1e6]]))
+    fmodel_true.run()
+    fmodel_true_tp = fmodel_true.get_turbine_powers() / 1000 # kW expected
+    assert np.allclose(fs_tp, fmodel_true_tp.flatten().tolist())
 
     # Mixed power setpoints
     fs_ws, fs_wd, fs_tp, fs_twd = floris_standin.get_step(5.0, power_setpoints=[None, 1e3])
-    fi_true.set(power_setpoints=np.array([[None, 1e6]]))
-    fi_true.run()
-    fi_true_tp = fi_true.get_turbine_powers() / 1000
-    assert np.allclose(fs_tp, fi_true_tp.flatten().tolist())
+    fmodel_true.set(power_setpoints=np.array([[None, 1e6]]))
+    fmodel_true.run()
+    fmodel_true_tp = fmodel_true.get_turbine_powers() / 1000
+    assert np.allclose(fs_tp, fmodel_true_tp.flatten().tolist())
 
     # Test with invalid combination of yaw angles and power setpoints
     with pytest.raises(ValueError):
@@ -169,17 +169,17 @@ def test_FlorisStandin_get_step_power_setpoints():
     )
     floris_power_setpoints = np.array([power_setpoints])
     floris_power_setpoints[0,1] *= 1e3 
-    fi_true.set(
+    fmodel_true.set(
         yaw_angles=default_wind_direction - np.array([yaw_angles]),
         power_setpoints=floris_power_setpoints
     )
-    fi_true.run()
-    fi_true_tp = fi_true.get_turbine_powers() / 1000
-    assert np.allclose(fs_tp, fi_true_tp.flatten().tolist())
+    fmodel_true.run()
+    fmodel_true_tp = fmodel_true.get_turbine_powers() / 1000
+    assert np.allclose(fs_tp, fmodel_true_tp.flatten().tolist())
 
 
 def test_FlorisStandin_with_standin_data_yaw_angles():
-    floris_standin = FlorisStandin(CONFIG, AMR_INPUT, AMR_EXTERNAL_DATA)
+    floris_standin = FlorisStandin(CONFIG, AMR_INPUT, AMR_EXTERNAL_DATA, smoothing_coefficient=0.0)
 
     yaw_angles_all = [
         [240.0, 240.0],
@@ -231,7 +231,7 @@ def test_FlorisStandin_with_standin_data_yaw_angles():
     assert fs_tp_all[9, :].sum() > fs_tp_all[7, :].sum()
 
 def test_FlorisStandin_with_standin_data_power_setpoints():
-    floris_standin = FlorisStandin(CONFIG, AMR_INPUT, AMR_EXTERNAL_DATA)
+    floris_standin = FlorisStandin(CONFIG, AMR_INPUT, AMR_EXTERNAL_DATA, smoothing_coefficient=0.0)
 
     power_setpoints_all = [
         [None, None],
@@ -276,3 +276,25 @@ def test_FlorisStandin_with_standin_data_power_setpoints():
     assert (fs_tp_all[6, :] == 1e3).all()
     assert fs_tp_all[7, 0] > 1e3
     assert fs_tp_all[7, 1] <= 1e3
+
+def test_FlorisStandin_smoothing_coefficient():
+    floris_standin_no_smoothing = FlorisStandin(CONFIG, AMR_INPUT, smoothing_coefficient=0.0)
+    floris_standin_default_smoothing = FlorisStandin(CONFIG, AMR_INPUT)
+    floris_standin_heavy_smoothing = FlorisStandin(CONFIG, AMR_INPUT, smoothing_coefficient=0.9)
+
+    # Start at zero power
+    floris_standin_no_smoothing.turbine_powers_prev = np.zeros(2)
+    floris_standin_default_smoothing.turbine_powers_prev = np.zeros(2)
+    floris_standin_heavy_smoothing.turbine_powers_prev = np.zeros(2)
+
+    # Step forward
+    fs_tp_no_smoothing = floris_standin_no_smoothing.get_step(1.0)[2]
+    fs_tp_default_smoothing = floris_standin_default_smoothing.get_step(1.0)[2]
+    fs_tp_heavy_smoothing = floris_standin_heavy_smoothing.get_step(1.0)[2]
+
+    # Check smoothing ordering correct
+    assert (np.array(fs_tp_no_smoothing) > np.array(fs_tp_default_smoothing)).all()
+    assert (np.array(fs_tp_default_smoothing) > np.array(fs_tp_heavy_smoothing)).all()
+
+    # Check magnitude is correct
+    assert np.allclose(0.1*np.array(fs_tp_no_smoothing), np.array(fs_tp_heavy_smoothing))