diff --git a/docs/_toc.yml b/docs/_toc.yml index 784be504d..4b78b0821 100644 --- a/docs/_toc.yml +++ b/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 diff --git a/docs/layout_optimization.md b/docs/layout_optimization.md new file mode 100644 index 000000000..361856579 --- /dev/null +++ b/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) diff --git a/docs/plot_complex_docs.png b/docs/plot_complex_docs.png new file mode 100644 index 000000000..bd4871298 Binary files /dev/null and b/docs/plot_complex_docs.png differ diff --git a/examples/examples_layout_optimization/003_genetic_random_search.py b/examples/examples_layout_optimization/003_genetic_random_search.py new file mode 100644 index 000000000..0955f151e --- /dev/null +++ b/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() diff --git a/floris/flow_visualization.py b/floris/flow_visualization.py index 57ae105a9..b893b172a 100644 --- a/floris/flow_visualization.py +++ b/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]) diff --git a/floris/optimization/layout_optimization/layout_optimization_base.py b/floris/optimization/layout_optimization/layout_optimization_base.py index dd9afaae3..99977c2f5 100644 --- a/floris/optimization/layout_optimization/layout_optimization_base.py +++ b/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 diff --git a/floris/optimization/layout_optimization/layout_optimization_pyoptsparse.py b/floris/optimization/layout_optimization/layout_optimization_pyoptsparse.py index 3a87dff70..794795767 100644 --- a/floris/optimization/layout_optimization/layout_optimization_pyoptsparse.py +++ b/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, diff --git a/floris/optimization/layout_optimization/layout_optimization_random_search.py b/floris/optimization/layout_optimization/layout_optimization_random_search.py new file mode 100644 index 000000000..92e5d7f94 --- /dev/null +++ b/floris/optimization/layout_optimization/layout_optimization_random_search.py @@ -0,0 +1,707 @@ +# Copyright 2022 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 + + +from multiprocessing import Pool +from time import perf_counter as timerpc + +import matplotlib.pyplot as plt +import numpy as np +from numpy import random +from scipy.spatial.distance import cdist, pdist +from shapely.geometry import Point, Polygon + +from floris import FlorisModel +from floris.optimization.yaw_optimization.yaw_optimizer_geometric import ( + YawOptimizationGeometric, +) + +from .layout_optimization_base import LayoutOptimization + + +def _load_local_floris_object( + fmodel_dict, + wind_data=None, +): + # Load local FLORIS object + fmodel = FlorisModel(fmodel_dict) + fmodel.set(wind_data=wind_data) + return fmodel + +def test_min_dist(layout_x, layout_y, min_dist): + coords = np.array([layout_x,layout_y]).T + dist = pdist(coords) + return dist.min() >= min_dist + +def test_point_in_bounds(test_x, test_y, poly_outer): + return poly_outer.contains(Point(test_x, test_y)) + +# Return in MW +def _get_objective( + layout_x, + layout_y, + fmodel, + yaw_angles=None, + use_value=False +): + fmodel.set( + layout_x=layout_x, + layout_y=layout_y, + yaw_angles=yaw_angles + ) + fmodel.run() + + return fmodel.get_farm_AVP() if use_value else fmodel.get_farm_AEP() + +def _gen_dist_based_init( + N, # Number of turbins to place + step_size, #m, courseness of search grid + poly_outer, # Polygon of outer boundary + min_x, + max_x, + min_y, + max_y, + s +): + """ + Generates an initial layout by randomly placing + the first turbine than placing the remaining turbines + as far as possible from the existing turbines. + """ + + # Set random seed + np.random.seed(s) + + # Choose the initial point randomly + init_x = float(random.randint(int(min_x),int(max_x))) + init_y = float(random.randint(int(min_y),int(max_y))) + while not (poly_outer.contains(Point([init_x,init_y]))): + init_x = float(random.randint(int(min_x),int(max_x))) + init_y = float(random.randint(int(min_y),int(max_y))) + + # Intialize the layout arrays + layout_x = np.array([init_x]) + layout_y = np.array([init_y]) + layout = np.array([layout_x, layout_y]).T + + # Now add the remaining points + for i in range(1,N): + + print("Placing turbine {0} of {1}.".format(i, N)) + # Add a new turbine being as far as possible from current + max_dist = 0. + for x in np.arange(min_x, max_x,step_size): + for y in np.arange(min_y, max_y,step_size): + if poly_outer.contains(Point([x,y])): + test_dist = cdist([[x,y]],layout) + min_dist = np.min(test_dist) + if min_dist > max_dist: + max_dist = min_dist + save_x = x + save_y = y + + # Add point to the layout + layout_x = np.append(layout_x,[save_x]) + layout_y = np.append(layout_y,[save_y]) + layout = np.array([layout_x, layout_y]).T + + # Return the layout + return layout_x, layout_y + +class LayoutOptimizationRandomSearch(LayoutOptimization): + def __init__( + self, + fmodel, + boundaries, + min_dist=None, + min_dist_D=None, + distance_pmf=None, + n_individuals=4, + seconds_per_iteration=60., + total_optimization_seconds=600., + interface="multiprocessing", # Options are 'multiprocessing', 'mpi4py', None + max_workers=None, + grid_step_size=100., + relegation_number=1, + enable_geometric_yaw=False, + use_dist_based_init=True, + random_seed=None, + use_value=False, + ): + """ + _summary_ + + Args: + fmodel (_type_): _description_ + boundaries (iterable(float, float)): Pairs of x- and y-coordinates + that represent the boundary's vertices (m). + 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. + min_dist_D (float, optional): The minimum distance to be maintained + between turbines during the optimization, specified as a multiple + of the rotor diameter. + distance_pmf (dict, optional): Probability mass function describing the + length of steps in the random search. Specified as a dictionary with + keys "d" (array of step distances, specified in meters) and "p" + (array of probability of occurrence, should sum to 1). Defaults to + uniform probability between 0.5D and 2D, with some extra mass + to encourage large changes. + n_individuals (int, optional): The number of individuals to use in the + optimization. Defaults to 4. + seconds_per_iteration (float, optional): The number of seconds to + run each step of the optimization for. Defaults to 60. + total_optimization_seconds (float, optional): The total number of + seconds to run the optimization for. Defaults to 600. + interface (str): Parallel computing interface to leverage. Recommended is 'concurrent' + or 'multiprocessing' for local (single-system) use, and 'mpi4py' for high + performance computing on multiple nodes. Defaults to 'multiprocessing'. + max_workers (int): Number of parallel workers, typically equal to the number of cores + you have on your system or HPC. Defaults to None, which will use all + available cores. + grid_step_size (float): The coarseness of the grid used to generate the initial layout. + Defaults to 100. + relegation_number (int): The number of the lowest performing individuals to be replaced + with new individuals generated from the best performing individual. Must + be less than n_individuals / 2. Defaults to 1. + enable_geometric_yaw (bool): Use geometric yaw code to determine approximate wake + steering yaw angles during layout optimization routine. Defaults to False. + use_dist_based_init (bool): Generate initial layouts automatically by placing turbines + as far apart as possible. + random_seed (int or None): Random seed for reproducibility. Defaults to None. + use_value (bool, optional): If True, the layout optimization objective + is to maximize annual value production using the value array in the + FLORIS model's WindData object. If False, the optimization + objective is to maximize AEP. Defaults to False. + """ + # The parallel computing interface to use + if interface == "mpi4py": + import mpi4py.futures as mp + self._PoolExecutor = mp.MPIPoolExecutor + elif interface == "multiprocessing": + import multiprocessing as mp + self._PoolExecutor = mp.Pool + if max_workers is None: + max_workers = mp.cpu_count() + elif interface is None: + if n_individuals > 1 or (max_workers is not None and max_workers > 1): + print( + "Parallelization not possible with interface=None. " + +"Reducing n_individuals to 1 and ignoring max_workers." + ) + self._PoolExecutor = None + max_workers = None + n_individuals = 1 + + # elif interface == "concurrent": + # from concurrent.futures import ProcessPoolExecutor + # self._PoolExecutor = ProcessPoolExecutor + else: + raise ValueError( + f"Interface '{interface}' not recognized. " + "Please use ' 'multiprocessing' or 'mpi4py'." + ) + + # Store the max_workers + self.max_workers = max_workers + + # Store the interface + self.interface = interface + + # Set and store the random seed + self.random_seed = random_seed + + # Confirm the relegation_number is valid + if relegation_number > n_individuals / 2: + raise ValueError("relegation_number must be less than n_individuals / 2.") + self.relegation_number = relegation_number + + # Store the rotor diameter and number of turbines + self.D = fmodel.core.farm.rotor_diameters.max() + if not all(fmodel.core.farm.rotor_diameters == self.D): + self.logger.warning("Using largest rotor diameter for min_dist_D and distance_pmf.") + self.N_turbines = fmodel.n_turbines + + # Make sure not both min_dist and min_dist_D are defined + if min_dist is not None and min_dist_D is not None: + raise ValueError("Only one of min_dist and min_dist_D can be defined.") + + # If min_dist_D is defined, convert to min_dist + if min_dist_D is not None: + min_dist = min_dist_D * self.D + + super().__init__( + fmodel, + boundaries, + min_dist=min_dist, + enable_geometric_yaw=enable_geometric_yaw, + use_value=use_value, + ) + if use_value: + self._obj_name = "value" + self._obj_unit = "" + else: + self._obj_name = "AEP" + self._obj_unit = "[GWh]" + + # Save min_dist_D + self.min_dist_D = self.min_dist / self.D + + # Process and save the step distribution + self._process_dist_pmf(distance_pmf) + + # Store the Core dictionary + self.fmodel_dict = self.fmodel.core.as_dict() + + # Save the grid step size + self.grid_step_size = grid_step_size + + # Save number of individuals + self.n_individuals = n_individuals + + # Store the initial locations + self.x_initial = self.fmodel.layout_x + self.y_initial = self.fmodel.layout_y + + # Store the total optimization seconds + self.total_optimization_seconds = total_optimization_seconds + + # Store the seconds per iteration + self.seconds_per_iteration = seconds_per_iteration + + # Get the initial objective value + self.x = self.x_initial # Required by _get_geoyaw_angles + self.y = self.y_initial # Required by _get_geoyaw_angles + self.objective_initial = _get_objective( + self.x_initial, + self.y_initial, + self.fmodel, + self._get_geoyaw_angles(), + self.use_value, + ) + + # Initialize the objective statistics + self.objective_mean = self.objective_initial + self.objective_median = self.objective_initial + self.objective_max = self.objective_initial + self.objective_min = self.objective_initial + + # Initialize the numpy arrays which will hold the candidate layouts + # these will have dimensions n_individuals x N_turbines + self.x_candidate = np.zeros((self.n_individuals, self.N_turbines)) + self.y_candidate = np.zeros((self.n_individuals, self.N_turbines)) + + # Initialize the array which will hold the objective function values for each candidate + self.objective_candidate = np.zeros(self.n_individuals) + + # Initialize the iteration step + self.iteration_step = -1 + + # Initialize the optimization time + self.opt_time_start = timerpc() + self.opt_time = 0 + + # Generate the initial layouts + if use_dist_based_init: + self._generate_initial_layouts() + else: + print(f'Using supplied initial layout for {self.n_individuals} individuals.') + for i in range(self.n_individuals): + self.x_candidate[i, :] = self.x_initial + self.y_candidate[i, :] = self.y_initial + self.objective_candidate[i] = self.objective_initial + + # Evaluate the initial optimization step + self._evaluate_opt_step() + + # Delete stored x and y to avoid confusion + del self.x, self.y + + def describe(self): + print("Random Layout Optimization") + print(f"Number of turbines to optimize = {self.N_turbines}") + print(f"Minimum distance between turbines = {self.min_dist_D} [D], {self.min_dist} [m]") + print(f"Number of individuals = {self.n_individuals}") + print(f"Seconds per iteration = {self.seconds_per_iteration}") + print(f"Initial {self._obj_name} = {self.objective_initial/1e9:.1f} {self._obj_unit}") + + def _process_dist_pmf(self, dist_pmf): + """ + Check validity of pmf and assign default if none provided. + """ + if dist_pmf is None: + jump_dist = np.min([self.xmax-self.xmin, self.ymax-self.ymin])/2 + jump_prob = 0.05 + + d = np.append(np.linspace(0.0, 2.0*self.D, 99), jump_dist) + p = np.append((1-jump_prob)/len(d)*np.ones(len(d)-1), jump_prob) + p = p / p.sum() + dist_pmf = {"d":d, "p":p} + + # Check correct keys are provided + if not all(k in dist_pmf for k in ("d", "p")): + raise KeyError("distance_pmf must contains keys \"d\" (step distance)"+\ + " and \"p\" (probability of occurrence).") + + # Check entries are in the correct form + if not hasattr(dist_pmf["d"], "__len__") or not hasattr(dist_pmf["d"], "__len__")\ + or len(dist_pmf["d"]) != len(dist_pmf["p"]): + raise TypeError("distance_pmf entries should be numpy arrays or lists"+\ + " of equal length.") + + if not np.isclose(dist_pmf["p"].sum(), 1): + print("Probability mass function does not sum to 1. Normalizing.") + dist_pmf["p"] = np.array(dist_pmf["p"]) / np.array(dist_pmf["p"]).sum() + + self.distance_pmf = dist_pmf + + def _evaluate_opt_step(self): + + # Sort the candidate layouts by objective function value + sorted_indices = np.argsort(self.objective_candidate)[::-1] # Decreasing order + self.objective_candidate = self.objective_candidate[sorted_indices] + self.x_candidate = self.x_candidate[sorted_indices] + self.y_candidate = self.y_candidate[sorted_indices] + + # Update the optimization time + self.opt_time = timerpc() - self.opt_time_start + + # Update the optimizations step + self.iteration_step += 1 + + # Update the objective statistics + self.objective_mean = np.mean(self.objective_candidate) + self.objective_median = np.median(self.objective_candidate) + self.objective_max = np.max(self.objective_candidate) + self.objective_min = np.min(self.objective_candidate) + + # Report the results + increase_mean = ( + 100 * (self.objective_mean - self.objective_initial) / self.objective_initial + ) + increase_median = ( + 100 * (self.objective_median - self.objective_initial) / self.objective_initial + ) + increase_max = 100 * (self.objective_max - self.objective_initial) / self.objective_initial + increase_min = 100 * (self.objective_min - self.objective_initial) / self.objective_initial + print("=======================================") + print(f"Optimization step {self.iteration_step:+.1f}") + print(f"Optimization time = {self.opt_time:+.1f} [s]") + print( + f"Mean {self._obj_name} = {self.objective_mean/1e9:.1f}" + f" {self._obj_unit} ({increase_mean:+.2f}%)" + ) + print( + f"Median {self._obj_name} = {self.objective_median/1e9:.1f}" + f" {self._obj_unit} ({increase_median:+.2f}%)" + ) + print( + f"Max {self._obj_name} = {self.objective_max/1e9:.1f}" + f" {self._obj_unit} ({increase_max:+.2f}%)" + ) + print( + f"Min {self._obj_name} = {self.objective_min/1e9:.1f}" + f" {self._obj_unit} ({increase_min:+.2f}%)" + ) + print("=======================================") + + # Replace the relegation_number worst performing layouts with relegation_number + # best layouts + if self.relegation_number > 0: + self.objective_candidate[-self.relegation_number:] = ( + self.objective_candidate[:self.relegation_number] + ) + self.x_candidate[-self.relegation_number:] = self.x_candidate[:self.relegation_number] + self.y_candidate[-self.relegation_number:] = self.y_candidate[:self.relegation_number] + + + # Private methods + def _generate_initial_layouts(self): + """ + This method generates n_individuals initial layout of turbines. It does + this by calling the _generate_random_layout method within a multiprocessing + pool. + """ + + # Set random seed for initial layout + if self.random_seed is None: + multi_random_seeds = [None]*self.n_individuals + else: + multi_random_seeds = [23 + i for i in range(self.n_individuals)] + # 23 is just an arbitrary choice to ensure different random seeds + # to the evaluation code + + print(f'Generating {self.n_individuals} initial layouts...') + t1 = timerpc() + # Generate the multiargs for parallel execution + multiargs = [ + (self.N_turbines, + self.grid_step_size, + self._boundary_polygon, + self.xmin, + self.xmax, + self.ymin, + self.ymax, + multi_random_seeds[i]) + for i in range(self.n_individuals) + ] + + if self._PoolExecutor: # Parallelized + with self._PoolExecutor(self.max_workers) as p: + # This code is not currently necessary, but leaving in case implement + # concurrent later, based on parallel_computing_interface.py + if (self.interface == "mpi4py") or (self.interface == "multiprocessing"): + out = p.starmap(_gen_dist_based_init, multiargs) + else: # Parallelization not activated + out = [_gen_dist_based_init(*multiargs[0])] + + # Unpack out into the candidate layouts + for i in range(self.n_individuals): + self.x_candidate[i, :] = out[i][0] + self.y_candidate[i, :] = out[i][1] + + # Get the objective function values for each candidate layout + for i in range(self.n_individuals): + self.objective_candidate[i] = _get_objective( + self.x_candidate[i, :], + self.y_candidate[i, :], + self.fmodel, + self._get_geoyaw_angles(), + self.use_value, + ) + + t2 = timerpc() + print(f" Time to generate initial layouts: {t2-t1:.3f} s") + + def _get_initial_and_final_locs(self): + x_initial = self.x_initial + y_initial = self.y_initial + x_opt = self.x_opt + y_opt = self.y_opt + return x_initial, y_initial, x_opt, y_opt + + + # Public methods + + def optimize(self): + """ + Perform the optimization + """ + print(f'Optimizing using {self.n_individuals} individuals.') + opt_start_time = timerpc() + opt_stop_time = opt_start_time + self.total_optimization_seconds + sim_time = 0 + + self.objective_candidate_log = [self.objective_candidate.copy()] + self.num_objective_calls_log = [] + self._num_objective_calls = [0]*self.n_individuals + + while timerpc() < opt_stop_time: + + # Set random seed for the main loop + if self.random_seed is None: + multi_random_seeds = [None]*self.n_individuals + else: + multi_random_seeds = [55 + self.iteration_step + i + for i in range(self.n_individuals)] + # 55 is just an arbitrary choice to ensure different random seeds + # to the initialization code + + # Update the optimization time + sim_time = timerpc() - opt_start_time + print(f'Optimization time: {sim_time:.1f} s / {self.total_optimization_seconds:.1f} s') + + + # Generate the multiargs for parallel execution of single individual optimization + multiargs = [ + (self.seconds_per_iteration, + self.objective_candidate[i], + self.x_candidate[i, :], + self.y_candidate[i, :], + self.fmodel_dict, + self.fmodel.wind_data, + self.min_dist, + self._boundary_polygon, + self.distance_pmf, + self.enable_geometric_yaw, + multi_random_seeds[i], + self.use_value + ) + for i in range(self.n_individuals) + ] + + # Run the single individual optimization in parallel + if self._PoolExecutor: # Parallelized + with self._PoolExecutor(self.max_workers) as p: + out = p.starmap(_single_individual_opt, multiargs) + else: # Parallelization not activated + out = [_single_individual_opt(*multiargs[0])] + + # Unpack the results + for i in range(self.n_individuals): + self.objective_candidate[i] = out[i][0] + self.x_candidate[i, :] = out[i][1] + self.y_candidate[i, :] = out[i][2] + self._num_objective_calls[i] = out[i][3] + self.objective_candidate_log.append(self.objective_candidate) + self.num_objective_calls_log.append(self._num_objective_calls) + + # Evaluate the individuals for this step + self._evaluate_opt_step() + + # Finalize the result + self.objective_final = self.objective_candidate[0] + self.x_opt = self.x_candidate[0, :] + self.y_opt = self.y_candidate[0, :] + + # Print the final result + increase = 100 * (self.objective_final - self.objective_initial) / self.objective_initial + print( + f"Final {self._obj_name} = {self.objective_final/1e9:.1f}" + f" {self._obj_unit} ({increase:+.2f}%)" + ) + + return self.objective_final, self.x_opt, self.y_opt + + + # Helpful visualizations + def plot_distance_pmf(self, ax=None): + """ + Tool to check the used distance pmf. + """ + + if ax is None: + _, ax = plt.subplots(1,1) + + ax.stem(self.distance_pmf["d"], self.distance_pmf["p"], linefmt="k-") + ax.grid(True) + ax.set_xlabel("Step distance [m]") + ax.set_ylabel("Probability") + + return ax + + + +def _single_individual_opt( + seconds_per_iteration, + initial_objective, + layout_x, + layout_y, + fmodel_dict, + wind_data, + min_dist, + poly_outer, + dist_pmf, + enable_geometric_yaw, + s, + use_value +): + # Set random seed + np.random.seed(s) + + # Initialize the optimization time + single_opt_start_time = timerpc() + stop_time = single_opt_start_time + seconds_per_iteration + + num_objective_calls = 0 + + # Get the fmodel + fmodel_ = _load_local_floris_object(fmodel_dict, wind_data) + + # Initialize local variables + num_turbines = len(layout_x) + get_new_point = True # Will always be true, due to hardcoded use_momentum + current_objective = initial_objective + + # Establish geometric yaw optimizer, if desired + if enable_geometric_yaw: + yaw_opt = YawOptimizationGeometric( + fmodel_, + minimum_yaw_angle=-30.0, + maximum_yaw_angle=30.0, + ) + else: # yaw_angles will always be none + yaw_angles = None + + # We have a beta feature to maintain momentum, i.e., if a move improves + # the objective, we try to keep moving in that direction. This is currently + # disabled. + use_momentum = False + + # Loop as long as we've not hit the stop time + while timerpc() < stop_time: + + if not use_momentum: + get_new_point = True + + if get_new_point: #If the last test wasn't successful + + # Randomly select a turbine to nudge + tr = random.randint(0,num_turbines-1) + + # Randomly select a direction to nudge in (uniform direction) + rand_dir = np.random.uniform(low=0.0, high=2*np.pi) + + # Randomly select a distance to travel according to pmf + rand_dist = np.random.choice(dist_pmf["d"], p=dist_pmf["p"]) + + # Get a new test point + test_x = layout_x[tr] + np.cos(rand_dir) * rand_dist + test_y = layout_y[tr] + np.sin(rand_dir) * rand_dist + + # In bounds? + if not test_point_in_bounds(test_x, test_y, poly_outer): + get_new_point = True + continue + + # Make a new layout + original_x = layout_x[tr] + original_y = layout_y[tr] + layout_x[tr] = test_x + layout_y[tr] = test_y + + # Acceptable distances? + if not test_min_dist(layout_x, layout_y,min_dist): + # Revert and continue + layout_x[tr] = original_x + layout_y[tr] = original_y + get_new_point = True + continue + + # Does it improve the objective? + if enable_geometric_yaw: # Select appropriate yaw angles + yaw_opt.fmodel_subset.set(layout_x=layout_x, layout_y=layout_y) + df_opt = yaw_opt.optimize() + yaw_angles = np.vstack(df_opt['yaw_angles_opt']) + + num_objective_calls += 1 + test_objective = _get_objective(layout_x, layout_y, fmodel_, yaw_angles, use_value) + + if test_objective > current_objective: + # Accept the change + current_objective = test_objective + + # If not a random point this cycle and it did improve things + # try not getting a new point + # Feature is currently disabled by use_momentum flag + get_new_point = False + + else: + # Revert the change + layout_x[tr] = original_x + layout_y[tr] = original_y + get_new_point = True + + # Return the best result from this individual + return current_objective, layout_x, layout_y, num_objective_calls diff --git a/floris/optimization/layout_optimization/layout_optimization_scipy.py b/floris/optimization/layout_optimization/layout_optimization_scipy.py index f7ca643b1..ff0e30015 100644 --- a/floris/optimization/layout_optimization/layout_optimization_scipy.py +++ b/floris/optimization/layout_optimization/layout_optimization_scipy.py @@ -5,7 +5,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 LayoutOptimizationScipy(LayoutOptimization): @@ -13,27 +13,6 @@ class LayoutOptimizationScipy(LayoutOptimization): This class provides an interface for optimizing the layout of wind turbines using the Scipy optimization library. The optimization objective is to maximize annual energy production (AEP) or annual value production (AVP). - - - Args: - fmodel (FlorisModel): A FlorisModel object. - boundaries (iterable(float, float)): Pairs of x- and y-coordinates - that represent the boundary's vertices (m). - 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. - 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. - solver (str, optional): Sets the solver used by Scipy. Defaults to 'SLSQP'. - optOptions (dict, optional): Dictionary for setting the - optimization options. Defaults to None. - enable_geometric_yaw (bool, optional): If True, enables geometric yaw - optimization. Defaults to False. - use_value (bool, optional): If True, the layout optimization objective - is to maximize annual value production using the value array in the - FLORIS model's WindData object. If False, the optimization - objective is to maximize AEP. Defaults to False. """ def __init__( self, @@ -46,6 +25,31 @@ def __init__( enable_geometric_yaw=False, use_value=False, ): + """ + Args: + fmodel (FlorisModel): A FlorisModel object. + boundaries (iterable(float, float)): Pairs of x- and y-coordinates + that represent the boundary's vertices (m). + 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. + 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. + solver (str, optional): Sets the solver used by Scipy. Defaults to 'SLSQP'. + optOptions (dict, optional): Dictionary for setting the + optimization options. Defaults to None. + enable_geometric_yaw (bool, optional): If True, enables geometric yaw + optimization. Defaults to False. + use_value (bool, optional): If True, the layout optimization objective + is to maximize annual value production using the value array in the + FLORIS model's WindData object. If False, the optimization + objective is to maximize AEP. Defaults to False. + """ + if list_depth(boundaries) > 1 and hasattr(boundaries[0][0], "__len__"): + raise NotImplementedError( + "LayoutOptimizationScipy is not configured for multiple regions." + ) super().__init__( fmodel, @@ -88,6 +92,10 @@ def __init__( # Private methods def _optimize(self): + + self._num_aep_calls = 0 + self._aep_record = [] + self.residual_plant = minimize( self._obj_func, self.x0, @@ -112,11 +120,16 @@ def _obj_func(self, locs): yaw_angles = self._get_geoyaw_angles() self.fmodel.set_operation(yaw_angles=yaw_angles) self.fmodel.run() + self._num_aep_calls += 1 if self.use_value: - return -1 * self.fmodel.get_farm_AVP() / self.initial_AEP_or_AVP + val = -1 * self.fmodel.get_farm_AVP() / self.initial_AEP_or_AVP + self._aep_record.append(val) + return val else: - return -1 * self.fmodel.get_farm_AEP() / self.initial_AEP_or_AVP + aep = -1 * self.fmodel.get_farm_AEP() / self.initial_AEP_or_AVP + self._aep_record.append(aep) + return aep def _change_coordinates(self, locs): diff --git a/floris/optimization/yaw_optimization/yaw_optimizer_geometric.py b/floris/optimization/yaw_optimization/yaw_optimizer_geometric.py index e78d48c9d..ea68204b4 100644 --- a/floris/optimization/yaw_optimization/yaw_optimizer_geometric.py +++ b/floris/optimization/yaw_optimization/yaw_optimizer_geometric.py @@ -12,6 +12,8 @@ class YawOptimizationGeometric(YawOptimization): :py:class:`floris.optimization.general_library.YawOptimization` that is used to provide a rough estimate of optimal yaw angles based purely on the wind farm geometry. Main use case is for coupled layout and yaw optimization. + + See Stanely et al. (2023) for details: https://wes.copernicus.org/articles/8/1341/2023/ """ def __init__( diff --git a/floris/wind_data.py b/floris/wind_data.py index 926ca9e0e..5745eca54 100644 --- a/floris/wind_data.py +++ b/floris/wind_data.py @@ -4,7 +4,6 @@ from abc import abstractmethod from pathlib import Path -import matplotlib.cm as cm import matplotlib.pyplot as plt import numpy as np import pandas as pd @@ -614,7 +613,7 @@ def plot( wd_step = wd_bins[1] - wd_bins[0] # Get a color array - color_array = cm.get_cmap(color_map, len(ws_bins)) + color_array = plt.get_cmap(color_map, len(ws_bins)) for wd_idx, wd in enumerate(wd_bins): rects = [] @@ -1514,7 +1513,7 @@ def plot( _, ax = plt.subplots(subplot_kw={"polar": True}) # Get a color array - color_array = cm.get_cmap(color_map, len(var_bins)) + color_array = plt.get_cmap(color_map, len(var_bins)) for wd_idx, wd in enumerate(wd_bins): rects = [] diff --git a/tests/layout_optimization_integration_test.py b/tests/layout_optimization_integration_test.py index 0732b969c..18353a8f5 100644 --- a/tests/layout_optimization_integration_test.py +++ b/tests/layout_optimization_integration_test.py @@ -12,6 +12,9 @@ from floris.optimization.layout_optimization.layout_optimization_base import ( LayoutOptimization, ) +from floris.optimization.layout_optimization.layout_optimization_random_search import ( + LayoutOptimizationRandomSearch, +) from floris.optimization.layout_optimization.layout_optimization_scipy import ( LayoutOptimizationScipy, ) @@ -70,3 +73,22 @@ def test_base_class(caplog): LayoutOptimization(fmodel, boundaries, 5) LayoutOptimization(fmodel=fmodel, boundaries=boundaries, min_dist=5) + +def test_LayoutOptimizationRandomSearch(): + fmodel = FlorisModel(configuration=YAML_INPUT) + fmodel.set(layout_x=[0, 500], layout_y = [0, 0]) + + # 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)] + + layout_opt = LayoutOptimizationRandomSearch( + fmodel=fmodel, + boundaries=boundaries, + min_dist_D=5, + seconds_per_iteration=1, + total_optimization_seconds=1, + use_dist_based_init=False, + ) + + # Check that the optimization runs + layout_opt.optimize() diff --git a/tests/reg_tests/random_search_layout_opt_regression_test.py b/tests/reg_tests/random_search_layout_opt_regression_test.py new file mode 100644 index 000000000..f29b40f28 --- /dev/null +++ b/tests/reg_tests/random_search_layout_opt_regression_test.py @@ -0,0 +1,142 @@ + +import numpy as np +import pandas as pd + +from floris import FlorisModel, WindRose +from floris.optimization.layout_optimization.layout_optimization_random_search import ( + LayoutOptimizationRandomSearch, +) +from tests.conftest import ( + assert_results_arrays, +) + + +DEBUG = False +VELOCITY_MODEL = "gauss" +DEFLECTION_MODEL = "gauss" + +locations_baseline_aep = np.array( + [ + [0.0, 571.5416296, 1260.0], + [0.0, 496.57085993, 0.], + ] +) +baseline_aep = 44787987324.21652 + +locations_baseline_value = np.array( + [ + [387.0, 100.0, 200.0, 300.0], + [192.0, 300.0, 100.0, 300.0], + ] +) +baseline_value = 8780876351.32277 + + +def test_random_search_layout_opt(sample_inputs_fixture): + """ + The SciPy optimization method optimizes turbine layout using SciPy's minimize method. This test + compares the optimization results from the SciPy layout optimization for a simple farm with a + simple wind rose to stored baseline results. + """ + sample_inputs_fixture.core["wake"]["model_strings"]["velocity_model"] = VELOCITY_MODEL + sample_inputs_fixture.core["wake"]["model_strings"]["deflection_model"] = DEFLECTION_MODEL + + boundaries = [(0.0, 0.0), (0.0, 1000.0), (1000.0, 1000.0), (1000.0, 0.0), (0.0, 0.0)] + + fmodel = FlorisModel(sample_inputs_fixture.core) + wd_array = np.arange(0, 360.0, 5.0) + ws_array = 8.0 * np.ones_like(wd_array) + + wind_rose = WindRose( + wind_directions=wd_array, + wind_speeds=ws_array, + ti_table=0.1, + ) + D = 126.0 # Rotor diameter for the NREL 5 MW + fmodel.set( + layout_x=[0.0, 5 * D, 10 * D], + layout_y=[0.0, 0.0, 0.0], + wind_data=wind_rose + ) + + layout_opt = LayoutOptimizationRandomSearch( + fmodel=fmodel, + boundaries=boundaries, + min_dist_D=5, + seconds_per_iteration=1, + total_optimization_seconds=1, + use_dist_based_init=False, + random_seed=0, + ) + sol = layout_opt.optimize() + optimized_aep = sol[0] + locations_opt = np.array([sol[1], sol[2]]) + + if DEBUG: + print(locations_opt) + print(optimized_aep) + + assert_results_arrays(locations_opt, locations_baseline_aep) + assert np.isclose(optimized_aep, baseline_aep) + +def test_random_search_layout_opt_value(sample_inputs_fixture): + """ + This test compares the optimization results from the SciPy layout optimization for a simple + farm with a simple wind rose to stored baseline results, optimizing for annual value production + instead of AEP. The value of the energy produced depends on the wind direction, causing the + optimal layout to differ from the case where the objective is maximum AEP. In this case, because + the value is much higher when the wind is from the north or south, the turbines are staggered to + avoid wake interactions for northerly and southerly winds. + """ + sample_inputs_fixture.core["wake"]["model_strings"]["velocity_model"] = VELOCITY_MODEL + sample_inputs_fixture.core["wake"]["model_strings"]["deflection_model"] = DEFLECTION_MODEL + + boundaries = [(0.0, 0.0), (0.0, 400.0), (400.0, 400.0), (400.0, 0.0), (0.0, 0.0)] + + fmodel = FlorisModel(sample_inputs_fixture.core) + + # set wind conditions and values using a WindData object with the default uniform frequency + wd_array = np.arange(0, 360.0, 5.0) + ws_array = np.array([8.0]) + + # Define the value table such that the value of the energy produced is + # significantly higher when the wind direction is close to the north or + # south, and zero when the wind is from the east or west. + value_table = (0.5 + 0.5*np.cos(2*np.radians(wd_array)))**10 + value_table = value_table.reshape((len(wd_array),1)) + + wind_rose = WindRose( + wind_directions=wd_array, + wind_speeds=ws_array, + ti_table=0.1, + value_table=value_table + ) + + # Start with a rectangular 4-turbine array with 2D spacing + D = 126.0 # Rotor diameter for the NREL 5 MW + fmodel.set( + layout_x=200 + np.array([-1 * D, -1 * D, 1 * D, 1 * D]), + layout_y=200 + np.array([-1* D, 1 * D, -1 * D, 1 * D]), + wind_data=wind_rose, + ) + + layout_opt = LayoutOptimizationRandomSearch( + fmodel=fmodel, + boundaries=boundaries, + min_dist_D=5, + seconds_per_iteration=1, + total_optimization_seconds=1, + use_dist_based_init=True, + random_seed=0, + use_value=True, + ) + sol = layout_opt.optimize() + optimized_value = sol[0] + locations_opt = np.array([sol[1], sol[2]]) + + if DEBUG: + print(locations_opt) + print(optimized_value) + + assert_results_arrays(locations_opt, locations_baseline_value) + assert np.isclose(optimized_value, baseline_value)