Randomized layout optimization (#697)

* Add random optimization class

* Add temporary example

* rename example

* Adding support for arbitrary pmf.

* Plotting for pmf; update example.

* Respecting new examples in develop.

* Update base class to handle multiple regions of layout optimization.

* geoyaw option added.

* Terminology change: particle -> individual.

* Adding functionality for plotting the optimization boundary only.

* Adding grid.

* Adding visualizations.

* Allow random seed setting; non-parallel runs; log status.

* Updating link to Stanleys published paper.

* Enabling geometric yaw option.

* WIP; storing for hotfix.

* Bug fixed: updating layout when called.

* Ruff

* fi_subset needed updating rather than fi.

* WIP commit to change branch.

* runs.

* Store initial aep correctly.

* Minor update to defualt pmf'

* Runs, but freq data not passed correctly.

* Pipe through wind_data for freq information.

* Move example to subdirectory.

* Ruff.

* Working on adding tests; still some work to do on value optimization.

* Enable value optimization.

* Handling TODO items

* UserWarning -> ValueError; logger.warnings

* Simple documentation added.

* White space in docs.

* Remove cm.get_cmap in favor of plt.get_cmap

* Improve documentation and example.

* renable geometric yaw option for v4; fix scipy layoutopt args.

* remove self.x and self.y after initialization.

* Improve progress plots and add to example.

---------

Co-authored-by: misi9170 <[email protected]>

docs/_toc.yml

@@ -18,6 +18,7 @@ parts:
     - file: floating_wind_turbine
     - file: turbine_interaction
     - file: operation_models_user
+    - file: layout_optimization
     - file: input_reference_main
     - file: input_reference_turbine
     - file: examples

docs/layout_optimization.md

@@ -0,0 +1,81 @@
+
+(layout_optimization)=
+# Layout optimization
+
+The FLORIS package provides layout optimization tools to place turbines within a specified
+boundary area to optimize annual energy production (AEP) or wind plant value. Layout
+optimizers accept an instantiated `FlorisModel` and alter the turbine layouts in order to
+improve the objective function value (AEP or value).
+
+## Background
+
+Layout optimization entails placing turbines in a wind farm in a configuration that maximizes
+an objective function, usually the AEP. Turbines are moved to minimize their wake interactions
+in the most dominant wind directions, while respecting the boundaries of the area for turbine
+placement as well as minimum distance requirements between neighboring turbines.
+
+Mathematically, we represent this as a (nonconvex) optimization problem.
+Let $x = \{x_i\}_{i=1,\dots,N}$, $x_i \in \mathbb{R}^2$ represent the set of
+coordinates of turbines within a farm (that is, $x_i$ represents the $(X, Y)$
+location of turbine $i$). Further, let $R \subset \mathbb{R}^2$ be a closed
+region in which to place turbines. Finally, let $d$ represent the minimum
+allowable distance between two turbines. Then, the layout optimization problem
+is expressed as
+
+$$
+\begin{aligned}
+\underset{x}{\text{maximize}} & \:\: f(x)\\
+\text{subject to} & \:\: x \subset R \\
+& \:\: ||x_i - x_j|| \geq d \:\: \forall j = 1,\dots,N, j\neq i
+\end{aligned}
+$$
+
+Here, $||\cdot||$ denotes the Euclidean norm, and $f(x)$ is the cost function to be maximized.
+
+When maximizing the AEP, $f = \sum_w P(w, x)p_W(w)$, where $w$ is the wind condition bin
+(e.g., wind speed, wind direction pair); $P(w, x)$ is the power produced by the wind farm in
+condition $w$ with layout $x$; and $p_W(w)$ is the annual frequency of occurrence of
+condition $w$.
+
+Layout optimizers take iterative approaches to solving the layout optimization problem
+specified above. Optimization routines available in FLORIS are described below.
+
+## Scipy layout optimization
+The `LayoutOptimizationScipy` class is built around `scipy.optimize`s `minimize`
+routine, using the `SLSQP` solver by default. Options for adjusting
+`minimize`'s behavior are exposed to the user with the `optOptions` argument.
+Other options include enabling fast wake steering at each layout optimizer
+iteration with the `enable_geometric_yaw` argument, and switch from AEP
+optimization to value optimization with the `use_value` argument.
+
+## Genetic random search layout optimization
+The `LayoutOptimizationRandomSearch` class is a custom optimizer designed specifically for
+layout optimization via random perturbations of the turbine locations. It is designed to have
+the following features:
+- Robust to complex wind conditions and complex boundaries, including disjoint regions for
+turbine placement
+- Easy to parallelize and wrapped in a genetic algorithm for propagating candidate solutions
+- Simple to set up and tune for non-optimization experts
+- Set up to run cheap constraint checks prior to more expensive objective function evaluations
+to accelerate optimization
+
+The algorithm, described in full in an upcoming paper that will be linked here when it is
+publicly available, moves a random turbine and random distance in a random direction; checks
+that constraints are satisfied; evaluates the objective function (AEP or value); and then
+commits to the move if there is an objective function improvement. The main tuning parameter
+is the probability mass function for the random movement distance, which is a dictionary to be
+passed to the `distance_pmf` argument.
+
+The `distance_pmf` dictionary should contain two keys, each containing a 1D array of equal
+length: `"d"`, which specifies the perturbation distance _in units of the rotor diameter_,
+and `"p"`, which specifies the probability that the corresponding perturbation distance is
+chosen at any iteration of the random search algorithm. The `distance_pmf` can therefore be
+used to encourage or discourage more aggressive or more conservative movements, and to enable
+or disable jumps between disjoint regions for turbine placement.
+
+The figure below shows an example of the optimized layout of a farm using the GRS algorithm, with
+the black dots indicating the initial layout; red dots indicating the final layout; and blue
+shading indicating wind speed heterogeneity (lighter shade is lower wind speed, darker shade is
+higher wind speed). The progress of each of the genetic individuals in the optimization process is
+shown in the right-hand plot.
+![](plot_complex_docs.png)

examples/examples_layout_optimization/003_genetic_random_search.py

@@ -0,0 +1,82 @@
+"""Example: Layout optimization with genetic random search
+This example shows a layout optimization using the genetic random search
+algorithm. It provides options for the users to try different distance
+probability mass functions for the random search perturbations.
+"""
+
+import matplotlib.pyplot as plt
+import numpy as np
+from scipy.stats import gamma
+
+from floris import FlorisModel, WindRose
+from floris.optimization.layout_optimization.layout_optimization_random_search import (
+    LayoutOptimizationRandomSearch,
+)
+
+
+if __name__ == '__main__':
+    # Set up FLORIS
+    fmodel = FlorisModel('../inputs/gch.yaml')
+
+
+    # Setup 72 wind directions with a random 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.0
+    # Shape frequency distribution to match number of wind directions and wind speeds
+    freq = (
+        np.abs(
+            np.sort(
+                np.random.randn(len(wind_directions))
+            )
+        )
+        .reshape( ( len(wind_directions), len(wind_speeds) ) )
+    )
+    freq = freq / freq.sum()
+    fmodel.set(
+        wind_data=WindRose(
+            wind_directions=wind_directions,
+            wind_speeds=wind_speeds,
+            freq_table=freq,
+            ti_table=0.06
+        )
+    )
+
+    # Set the boundaries
+    # 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)]
+
+    # Set turbine locations to 4 turbines in a rectangle
+    D = 126.0 # rotor diameter for the NREL 5MW
+    layout_x = [0, 0, 6 * D, 6 * D]
+    layout_y = [0, 4 * D, 0, 4 * D]
+    fmodel.set(layout_x=layout_x, layout_y=layout_y)
+
+    # Perform the optimization
+    distance_pmf = None
+
+    # Other options that users can try
+    # 1.
+    # distance_pmf = {"d": [100, 1000], "p": [0.8, 0.2]}
+    # 2.
+    # p = gamma.pdf(np.linspace(0, 900, 91), 15, scale=20); p = p/p.sum()
+    # distance_pmf = {"d": np.linspace(100, 1000, 91), "p": p}
+
+    layout_opt = LayoutOptimizationRandomSearch(
+        fmodel,
+        boundaries,
+        min_dist_D=5.,
+        seconds_per_iteration=10,
+        total_optimization_seconds=60.,
+        distance_pmf=distance_pmf
+    )
+    layout_opt.describe()
+    layout_opt.plot_distance_pmf()
+
+    layout_opt.optimize()
+
+    layout_opt.plot_layout_opt_results()
+
+    layout_opt.plot_progress()
+
+    plt.show()

floris/flow_visualization.py

@@ -369,7 +369,7 @@ def plot_rotor_values(
         plot_rotor_values(floris.flow_field.w, findex=0, n_rows=1, ncols=4, show=True)
     """
 
-    cmap = plt.cm.get_cmap(name=cmap)
+    cmap = plt.get_cmap(name=cmap)
 
     if t_range is None:
         t_range = range(values.shape[1])

floris/optimization/layout_optimization/layout_optimization_base.py

@@ -1,7 +1,7 @@
 
 import matplotlib.pyplot as plt
 import numpy as np
-from shapely.geometry import LineString, Polygon
+from shapely.geometry import MultiPolygon, Polygon
 
 from floris import TimeSeries
 from floris.optimization.yaw_optimization.yaw_optimizer_geometric import (
@@ -45,13 +45,28 @@ def __init__(
         self.enable_geometric_yaw = enable_geometric_yaw
         self.use_value = use_value
 
-        self._boundary_polygon = Polygon(self.boundaries)
-        self._boundary_line = LineString(self.boundaries)
+        # Allow boundaries to be set either as a list of corners or as a
+        # nested list of corners (for seperable regions)
+        self.boundaries = boundaries
+        b_depth = list_depth(boundaries)
+
+        boundary_specification_error_msg = (
+            "boundaries should be a list of coordinates (specifed as (x,y) "+\
+            "tuples) or as a list of list of tuples (for seperable regions)."
+        )
 
-        self.xmin = np.min([tup[0] for tup in boundaries])
-        self.xmax = np.max([tup[0] for tup in boundaries])
-        self.ymin = np.min([tup[1] for tup in boundaries])
-        self.ymax = np.max([tup[1] for tup in boundaries])
+        if b_depth == 1:
+            self._boundary_polygon = MultiPolygon([Polygon(self.boundaries)])
+            self._boundary_line = self._boundary_polygon.boundary
+        elif b_depth == 2:
+            if not isinstance(self.boundaries[0][0], tuple):
+                raise TypeError(boundary_specification_error_msg)
+            self._boundary_polygon = MultiPolygon([Polygon(p) for p in self.boundaries])
+            self._boundary_line = self._boundary_polygon.boundary
+        else:
+            raise TypeError(boundary_specification_error_msg)
+
+        self.xmin, self.ymin, self.xmax, self.ymax = self._boundary_polygon.bounds
 
         # If no minimum distance is provided, assume a value of 2 rotor diameters
         if min_dist is None:
@@ -115,36 +130,106 @@ def optimize(self):
         sol = self._optimize()
         return sol
 
-    def plot_layout_opt_results(self):
+    def plot_layout_opt_results(self, plot_boundary_dict={}, ax=None, fontsize=16):
+
         x_initial, y_initial, x_opt, y_opt = self._get_initial_and_final_locs()
 
-        plt.figure(figsize=(9, 6))
-        fontsize = 16
-        plt.plot(x_initial, y_initial, "ob")
-        plt.plot(x_opt, y_opt, "or")
-        # plt.title('Layout Optimization Results', fontsize=fontsize)
-        plt.xlabel("x (m)", fontsize=fontsize)
-        plt.ylabel("y (m)", fontsize=fontsize)
-        plt.axis("equal")
-        plt.grid()
-        plt.tick_params(which="both", labelsize=fontsize)
-        plt.legend(
-            ["Old locations", "New locations"],
+        # Generate axis, if needed
+        if ax is None:
+            fig = plt.figure(figsize=(9,6))
+            ax = fig.add_subplot(111)
+            ax.set_aspect("equal")
+
+        default_plot_boundary_dict = {
+            "color":"None",
+            "alpha":1,
+            "edgecolor":"b",
+            "linewidth":2
+        }
+        plot_boundary_dict = {**default_plot_boundary_dict, **plot_boundary_dict}
+
+        self.plot_layout_opt_boundary(plot_boundary_dict, ax=ax)
+        ax.plot(x_initial, y_initial, "ob", label="Initial locations")
+        ax.plot(x_opt, y_opt, "or", label="New locations")
+        ax.set_xlabel("x (m)", fontsize=fontsize)
+        ax.set_ylabel("y (m)", fontsize=fontsize)
+        ax.grid(True)
+        ax.tick_params(which="both", labelsize=fontsize)
+        ax.legend(
             loc="lower center",
             bbox_to_anchor=(0.5, 1.01),
             ncol=2,
             fontsize=fontsize,
         )
 
-        verts = self.boundaries
-        for i in range(len(verts)):
-            if i == len(verts) - 1:
-                plt.plot([verts[i][0], verts[0][0]], [verts[i][1], verts[0][1]], "b")
-            else:
-                plt.plot(
-                    [verts[i][0], verts[i + 1][0]], [verts[i][1], verts[i + 1][1]], "b"
+        return ax
+
+    def plot_layout_opt_boundary(self, plot_boundary_dict={}, ax=None):
+
+        # Generate axis, if needed
+        if ax is None:
+            fig = plt.figure(figsize=(9,6))
+            ax = fig.add_subplot(111)
+            ax.set_aspect("equal")
+
+        default_plot_boundary_dict = {
+            "color":"k",
+            "alpha":0.1,
+            "edgecolor":None
+        }
+
+        plot_boundary_dict = {**default_plot_boundary_dict, **plot_boundary_dict}
+
+        for line in self._boundary_line.geoms:
+            xy = np.array(line.coords)
+            ax.fill(xy[:,0], xy[:,1], **plot_boundary_dict)
+        ax.grid(True)
+
+        return ax
+
+    def plot_progress(self, ax=None):
+
+        if not hasattr(self, "objective_candidate_log"):
+            raise NotImplementedError(
+                "plot_progress not yet configured for "+self.__class__.__name__
+            )
+
+        if ax is None:
+            _, ax = plt.subplots(1,1)
+
+        objective_log_array = np.array(self.objective_candidate_log)
+
+        if len(objective_log_array.shape) == 1: # Just one AEP candidate per step
+            ax.plot(np.arange(len(objective_log_array)), objective_log_array, color="k")
+        elif len(objective_log_array.shape) == 2: # Multiple AEP candidates per step
+            for i in range(objective_log_array.shape[1]):
+                ax.plot(
+                    np.arange(len(objective_log_array)),
+                    objective_log_array[:,i],
+                    color="lightgray"
                 )
 
+        ax.scatter(
+            np.zeros(objective_log_array.shape[1]),
+            objective_log_array[0,:],
+            color="b",
+            label="Initial"
+        )
+        ax.scatter(
+            objective_log_array.shape[0]-1,
+            objective_log_array[-1,:].max(),
+            color="r",
+            label="Final"
+        )
+
+        # Plot aesthetics
+        ax.grid(True)
+        ax.set_xlabel("Optimization step [-]")
+        ax.set_ylabel("Objective function")
+        ax.legend()
+
+        return ax
+
 
     ###########################################################################
     # Properties
@@ -165,3 +250,11 @@ def nturbs(self):
     @property
     def rotor_diameter(self):
         return self.fmodel.core.farm.rotor_diameters_sorted[0][0]
+
+# Helper functions
+
+def list_depth(x):
+    if isinstance(x, list) and len(x) > 0:
+        return 1 + max(list_depth(item) for item in x)
+    else:
+        return 0

floris/optimization/layout_optimization/layout_optimization_pyoptsparse.py

@@ -4,7 +4,7 @@
 from scipy.spatial.distance import cdist
 from shapely.geometry import Point
 
-from .layout_optimization_base import LayoutOptimization
+from .layout_optimization_base import LayoutOptimization, list_depth
 
 
 class LayoutOptimizationPyOptSparse(LayoutOptimization):
@@ -54,6 +54,10 @@ def __init__(
         enable_geometric_yaw=False,
         use_value=False,
     ):
+        if list_depth(boundaries) > 1 and hasattr(boundaries[0][0], "__len__"):
+            raise NotImplementedError(
+                "LayoutOptimizationPyOptSparse is not configured for multiple regions."
+            )
 
         super().__init__(
             fmodel,