Drive layout optimizations using WindData objects (#822)

examples/15_optimize_layout.py

@@ -4,7 +4,7 @@
 import matplotlib.pyplot as plt
 import numpy as np
 
-from floris.tools import FlorisInterface
+from floris.tools import FlorisInterface, WindRose
 from floris.tools.optimization.layout_optimization.layout_optimization_scipy import (
     LayoutOptimizationScipy,
 )
@@ -24,15 +24,22 @@
 file_dir = os.path.dirname(os.path.abspath(__file__))
 fi = FlorisInterface('inputs/gch.yaml')
 
-# Setup 72 wind directions with a random wind speed and frequency distribution
+# Setup 72 wind directions with a 1 wind speed and frequency distribution
 wind_directions = np.arange(0, 360.0, 5.0)
-np.random.seed(1)
-wind_speeds = 8.0 + np.random.randn(1) * 0.5 * np.ones_like(wind_directions)
+wind_speeds = np.array([8.0])
+
 # Shape frequency distribution to match number of wind directions and wind speeds
-freq = (np.abs(np.sort(np.random.randn(len(wind_directions)))))
-freq = freq / freq.sum()
+freq_table = np.zeros((len(wind_directions), len(wind_speeds)))
+np.random.seed(1)
+freq_table[:,0] = (np.abs(np.sort(np.random.randn(len(wind_directions)))))
+freq_table = freq_table / freq_table.sum()
 
-fi.set(wind_directions=wind_directions, wind_speeds=wind_speeds)
+# Establish a TimeSeries object
+wind_rose = WindRose(wind_directions=wind_directions,
+                     wind_speeds=wind_speeds,
+                     freq_table=freq_table)
+
+fi.set(wind_data=wind_rose)
 
 # The boundaries for the turbines, specified as vertices
 boundaries = [(0.0, 0.0), (0.0, 1000.0), (1000.0, 1000.0), (1000.0, 0.0), (0.0, 0.0)]
@@ -44,18 +51,19 @@
 fi.set(layout_x=layout_x, layout_y=layout_y)
 
 # Setup the optimization problem
-layout_opt = LayoutOptimizationScipy(fi, boundaries, freq=freq)
+layout_opt = LayoutOptimizationScipy(fi, boundaries, wind_data=wind_rose)
 
 # Run the optimization
 sol = layout_opt.optimize()
 
 # Get the resulting improvement in AEP
 print('... calcuating improvement in AEP')
 fi.run()
-base_aep = fi.get_farm_AEP(freq=freq) / 1e6
+base_aep = fi.get_farm_AEP_with_wind_data(wind_data=wind_rose) / 1e6
 fi.set(layout_x=sol[0], layout_y=sol[1])
 fi.run()
-opt_aep = fi.get_farm_AEP(freq=freq) / 1e6
+opt_aep = fi.get_farm_AEP_with_wind_data(wind_data=wind_rose)  / 1e6
+
 percent_gain = 100 * (opt_aep - base_aep) / base_aep
 
 # Print and plot the results

examples/16c_optimize_layout_with_heterogeneity.py

@@ -4,7 +4,7 @@
 import matplotlib.pyplot as plt
 import numpy as np
 
-from floris.tools import FlorisInterface
+from floris.tools import FlorisInterface, WindRose
 from floris.tools.optimization.layout_optimization.layout_optimization_scipy import (
     LayoutOptimizationScipy,
 )
@@ -28,12 +28,17 @@
 
 # Setup 2 wind directions (due east and due west)
 # and 1 wind speed with uniform probability
-wind_directions = [270., 90.]
+wind_directions = np.array([270., 90.])
 n_wds = len(wind_directions)
-wind_speeds = [8.0] * np.ones_like(wind_directions)
+wind_speeds = np.array([8.0])
 # Shape frequency distribution to match number of wind directions and wind speeds
-freq = np.ones((len(wind_directions), len(wind_speeds)))
-freq = freq / freq.sum()
+freq_table = np.ones((len(wind_directions), len(wind_speeds)))
+freq_table = freq_table / freq_table.sum()
+
+# Establish a TimeSeries object
+wind_rose = WindRose(wind_directions=wind_directions,
+                     wind_speeds=wind_speeds,
+                     freq_table=freq_table)
 
 # The boundaries for the turbines, specified as vertices
 D = 126.0 # rotor diameter for the NREL 5MW
@@ -66,8 +71,7 @@
 fi.set(
     layout_x=layout_x,
     layout_y=layout_y,
-    wind_directions=wind_directions,
-    wind_speeds=wind_speeds,
+    wind_data=wind_rose,
     heterogenous_inflow_config=heterogenous_inflow_config
 )
 
@@ -76,7 +80,7 @@
 layout_opt = LayoutOptimizationScipy(
     fi,
     boundaries,
-    freq=freq,
+    wind_data=wind_rose,
     min_dist=2*D,
     optOptions={"maxiter":maxiter}
 )
@@ -87,11 +91,13 @@
 
 # Get the resulting improvement in AEP
 print('... calcuating improvement in AEP')
+
 fi.run()
-base_aep = fi.get_farm_AEP(freq=freq) / 1e6
+base_aep = fi.get_farm_AEP_with_wind_data(wind_data=wind_rose) / 1e6
 fi.set(layout_x=sol[0], layout_y=sol[1])
 fi.run()
-opt_aep = fi.get_farm_AEP(freq=freq) / 1e6
+opt_aep = fi.get_farm_AEP_with_wind_data(wind_data=wind_rose) / 1e6
+
 percent_gain = 100 * (opt_aep - base_aep) / base_aep
 
 # Print and plot the results
@@ -115,7 +121,7 @@
 layout_opt = LayoutOptimizationScipy(
     fi,
     boundaries,
-    freq=freq,
+    wind_data=wind_rose,
     min_dist=2*D,
     enable_geometric_yaw=True,
     optOptions={"maxiter":maxiter}
@@ -127,10 +133,15 @@
 
 # Get the resulting improvement in AEP
 print('... calcuating improvement in AEP')
+
 fi.set(yaw_angles=np.zeros_like(layout_opt.yaw_angles))
-base_aep = fi.get_farm_AEP(freq=freq) / 1e6
+base_aep = fi.get_farm_AEP_with_wind_data(wind_data=wind_rose) / 1e6
 fi.set(layout_x=sol[0], layout_y=sol[1], yaw_angles=layout_opt.yaw_angles)
-opt_aep = fi.get_farm_AEP(freq=freq) / 1e6
+fi.run()
+opt_aep = fi.get_farm_AEP_with_wind_data(
+    wind_data=wind_rose
+) / 1e6
+
 percent_gain = 100 * (opt_aep - base_aep) / base_aep
 
 # Print and plot the results

floris/tools/__init__.py

@@ -40,7 +40,6 @@
 # from floris.tools import (
 #     cut_plane,
 #     floris_interface,
-#     interface_utilities,
 #     layout_visualization,
 #     optimization,
 #     plotting,

floris/tools/optimization/layout_optimization/layout_optimization_base.py

@@ -3,15 +3,33 @@
 import numpy as np
 from shapely.geometry import LineString, Polygon
 
+from floris.tools import TimeSeries
 from floris.tools.optimization.yaw_optimization.yaw_optimizer_geometric import (
     YawOptimizationGeometric,
 )
+from floris.tools.wind_data import WindDataBase
 
 from ....logging_manager import LoggingManager
 
 
 class LayoutOptimization(LoggingManager):
-    def __init__(self, fi, boundaries, min_dist=None, freq=None, enable_geometric_yaw=False):
+    """
+    Base class for layout optimization. This class should not be used directly
+    but should be subclassed by a specific optimization method.
+
+    Args:
+        fi (FlorisInterface): A FlorisInterface object.
+        boundaries (iterable(float, float)): Pairs of x- and y-coordinates
+            that represent the boundary's vertices (m).
+        wind_data (TimeSeries | WindRose): A TimeSeries or WindRose object
+            values.
+        min_dist (float, optional): The minimum distance to be maintained
+            between turbines during the optimization (m). If not specified,
+            initializes to 2 rotor diameters. Defaults to None.
+        enable_geometric_yaw (bool, optional): If True, enables geometric yaw
+            optimization. Defaults to False.
+    """
+    def __init__(self, fi, boundaries, wind_data, min_dist=None, enable_geometric_yaw=False):
         self.fi = fi.copy()
         self.boundaries = boundaries
         self.enable_geometric_yaw = enable_geometric_yaw
@@ -30,12 +48,13 @@ def __init__(self, fi, boundaries, min_dist=None, freq=None, enable_geometric_ya
         else:
             self.min_dist = min_dist
 
-        # If freq is not provided, give equal weight to all wind conditions
-        if freq is None:
-            self.freq = np.ones((self.fi.floris.flow_field.n_findex,))
-            self.freq = self.freq / self.freq.sum()
-        else:
-            self.freq = freq
+        # Check that wind_data is a WindDataBase object
+        if (not isinstance(wind_data, WindDataBase)):
+            raise ValueError(
+                "wind_data entry is not an object of WindDataBase"
+                " (eg TimeSeries, WindRose, WindTIRose)"
+            )
+        self.wind_data = wind_data
 
         # Establish geometric yaw class
         if self.enable_geometric_yaw:
@@ -45,7 +64,7 @@ def __init__(self, fi, boundaries, min_dist=None, freq=None, enable_geometric_ya
                 maximum_yaw_angle=30.0,
             )
 
-        self.initial_AEP = fi.get_farm_AEP(self.freq)
+        self.initial_AEP = fi.get_farm_AEP_with_wind_data(self.wind_data)
 
     def __str__(self):
         return "layout"

floris/tools/optimization/layout_optimization/layout_optimization_pyoptsparse.py

@@ -12,16 +12,16 @@ def __init__(
         self,
         fi,
         boundaries,
+        wind_data,
         min_dist=None,
-        freq=None,
         solver=None,
         optOptions=None,
         timeLimit=None,
         storeHistory='hist.hist',
         hotStart=None,
         enable_geometric_yaw=False,
     ):
-        super().__init__(fi, boundaries, min_dist=min_dist, freq=freq,
+        super().__init__(fi, boundaries, wind_data=wind_data, min_dist=min_dist,
                          enable_geometric_yaw=enable_geometric_yaw)
 
         self.x0 = self._norm(self.fi.layout_x, self.xmin, self.xmax)
@@ -99,7 +99,10 @@ def _obj_func(self, varDict):
         # Compute the objective function
         funcs = {}
         funcs["obj"] = (
-            -1 * self.fi.get_farm_AEP(self.freq) / self.initial_AEP
+
+            -1 * self.fi.get_farm_AEP_with_wind_data(self.wind_data)
+              / self.initial_AEP
+
         )
 
         # Compute constraints, if any are defined for the optimization

floris/tools/optimization/layout_optimization/layout_optimization_pyoptsparse_spread.py

@@ -12,15 +12,15 @@ def __init__(
         self,
         fi,
         boundaries,
+        wind_data,
         min_dist=None,
-        freq=None,
         solver=None,
         optOptions=None,
         timeLimit=None,
         storeHistory='hist.hist',
         hotStart=None
     ):
-        super().__init__(fi, boundaries, min_dist=min_dist, freq=freq)
+        super().__init__(fi, boundaries,  wind_data=wind_data, min_dist=min_dist)
         self._reinitialize(solver=solver, optOptions=optOptions)
 
         self.storeHistory = storeHistory
@@ -95,7 +95,6 @@ def _obj_func(self, varDict):
         funcs = {}
         funcs["obj"] = (
             -1 * self.mean_distance(self.x, self.y)
-           #  -1 * np.sum(self.fi.get_farm_power() * self.freq * 8760) / self.initial_AEP
         )
 
         # Compute constraints, if any are defined for the optimization

floris/tools/optimization/layout_optimization/layout_optimization_scipy.py

@@ -13,7 +13,7 @@ def __init__(
         self,
         fi,
         boundaries,
-        freq=None,
+        wind_data,
         bnds=None,
         min_dist=None,
         solver='SLSQP',
@@ -27,10 +27,8 @@ def __init__(
             fi (_type_): _description_
             boundaries (iterable(float, float)): Pairs of x- and y-coordinates
                 that represent the boundary's vertices (m).
-            freq (np.array): An array of the frequencies of occurance
-                correponding to each pair of wind direction and wind speed
+            wind_data (TimeSeries | WindRose): A TimeSeries or WindRose object
                 values. If None, equal weight is given to each pair of wind conditions
-                Defaults to None.
             bnds (iterable, optional): Bounds for the optimization
                 variables (pairs of min/max values for each variable (m)). If
                 none are specified, they are set to 0 and 1. Defaults to None.
@@ -41,7 +39,7 @@ def __init__(
             optOptions (dict, optional): Dicitonary for setting the
                 optimization options. Defaults to None.
         """
-        super().__init__(fi, boundaries, min_dist=min_dist, freq=freq,
+        super().__init__(fi, boundaries, min_dist=min_dist, wind_data=wind_data,
                     enable_geometric_yaw=enable_geometric_yaw)
 
         self.boundaries_norm = [
@@ -100,7 +98,10 @@ def _obj_func(self, locs):
         # Compute turbine yaw angles using PJ's geometric code (if enabled)
         yaw_angles = self._get_geoyaw_angles()
         self.fi.set(yaw_angles=yaw_angles)
-        return -1 * self.fi.get_farm_AEP(self.freq) /self.initial_AEP
+
+        return (-1 * self.fi.get_farm_AEP_with_wind_data(self.wind_data) /
+                self.initial_AEP)
+
 
     def _change_coordinates(self, locs):
         # Parse the layout coordinates

tests/layout_optimization_integration_test.py

@@ -0,0 +1,55 @@
+from pathlib import Path
+
+import numpy as np
+import pytest
+
+from floris.tools import (
+    TimeSeries,
+    WindRose,
+)
+from floris.tools.floris_interface import FlorisInterface
+from floris.tools.optimization.layout_optimization.layout_optimization_base import (
+    LayoutOptimization,
+)
+from floris.tools.optimization.layout_optimization.layout_optimization_scipy import (
+    LayoutOptimizationScipy,
+)
+from floris.tools.wind_data import WindDataBase
+
+
+TEST_DATA = Path(__file__).resolve().parent / "data"
+YAML_INPUT = TEST_DATA / "input_full.yaml"
+
+
+def test_base_class():
+    # Get a test fi
+    fi = FlorisInterface(configuration=YAML_INPUT)
+
+    # Set up a sample boundary
+    boundaries = [(0.0, 0.0), (0.0, 1000.0), (1000.0, 1000.0), (1000.0, 0.0), (0.0, 0.0)]
+
+    # Now initiate layout optimization with a frequency matrix passed in the 3rd position
+    # (this should fail)
+    freq = np.ones((5, 5))
+    freq = freq / freq.sum()
+    with pytest.raises(ValueError):
+        LayoutOptimization(fi, boundaries, freq, 5)
+
+    # Passing as a keyword freq to wind_data should also fail
+    with pytest.raises(ValueError):
+        LayoutOptimization(fi=fi, boundaries=boundaries, wind_data=freq, min_dist=5,)
+
+    time_series = TimeSeries(
+        wind_directions=fi.floris.flow_field.wind_directions,
+        wind_speeds=fi.floris.flow_field.wind_speeds,
+        turbulence_intensities=fi.floris.flow_field.turbulence_intensities,
+    )
+    wind_rose = time_series.to_wind_rose()
+
+    # Passing wind_data objects in the 3rd position should not fail
+    LayoutOptimization(fi, boundaries, time_series, 5)
+    LayoutOptimization(fi, boundaries, wind_rose, 5)
+
+    # Passing wind_data objects by keyword should not fail
+    LayoutOptimization(fi=fi, boundaries=boundaries, wind_data=time_series, min_dist=5)
+    LayoutOptimization(fi=fi, boundaries=boundaries, wind_data=wind_rose, min_dist=5)