From 3263e31a37867fd05fa0189477ef51ed688f8865 Mon Sep 17 00:00:00 2001 From: Paul Date: Tue, 15 Mar 2022 11:15:07 -0600 Subject: [PATCH 01/27] Add function to read in wind rose as csv --- floris/tools/wind_rose.py | 16 ++++++++++++++++ 1 file changed, 16 insertions(+) diff --git a/floris/tools/wind_rose.py b/floris/tools/wind_rose.py index 8f3afb56f..245279085 100644 --- a/floris/tools/wind_rose.py +++ b/floris/tools/wind_rose.py @@ -634,6 +634,22 @@ def make_wind_rose_from_user_data( self.internal_resample_wind_direction(wd=wd) return self.df + + def read_wind_rose_csv( + self, + filename + ): + + #Read in the csv + self.df = pd.read_csv(filename) + + # Renormalize the frequency column + self.df["freq_val"] = self.df["freq_val"] / self.df["freq_val"].sum() + + # Call the resample function in order to set all the internal variables + self.internal_resample_wind_speed(ws=self.df.ws.unique()) + self.internal_resample_wind_direction(wd=self.df.wd.unique()) + def make_wind_rose_from_user_dist( self, From 8505a6834f3d239ba264357ea3206d325e8747fb Mon Sep 17 00:00:00 2001 From: Paul Date: Tue, 15 Mar 2022 11:15:21 -0600 Subject: [PATCH 02/27] add aep function which accepts wind rose object --- floris/tools/floris_interface.py | 102 +++++++++++++++++++++++++++++++ 1 file changed, 102 insertions(+) diff --git a/floris/tools/floris_interface.py b/floris/tools/floris_interface.py index f82ce65be..e3419d6f4 100644 --- a/floris/tools/floris_interface.py +++ b/floris/tools/floris_interface.py @@ -772,6 +772,108 @@ def get_farm_AEP( return aep + def get_farm_AEP_wind_rose_class( + self, + wind_rose, + cut_in_wind_speed=0.001, + cut_out_wind_speed=None, + yaw_angles=None, + no_wake=False, + ) -> float: + """ + Estimate annual energy production (AEP) for distributions of wind speed, wind + direction, frequency of occurrence, and yaw offset. + + Args: + wind_rose (wind_rose): An object of the wind rose class + cut_in_wind_speed (float, optional): Wind speed in m/s below which + any calculations are ignored and the wind farm is known to + produce 0.0 W of power. Note that to prevent problems with the + wake models at negative / zero wind speeds, this variable must + always have a positive value. Defaults to 0.001 [m/s]. + cut_out_wind_speed (float, optional): Wind speed above which the + wind farm is known to produce 0.0 W of power. If None is + specified, will assume that the wind farm does not cut out + at high wind speeds. Defaults to None. + yaw_angles (NDArrayFloat | list[float] | None, optional): + The relative turbine yaw angles in degrees. If None is + specified, will assume that the turbine yaw angles are all + zero degrees for all conditions. Defaults to None. + no_wake: (bool, optional): When *True* updates the turbine + quantities without calculating the wake or adding the wake to + the flow field. This can be useful when quantifying the loss + in AEP due to wakes. Defaults to *False*. + + Returns: + float: + The Annual Energy Production (AEP) for the wind farm in + watt-hours. + """ + + # Hold the starting values of wind speed and direction + wind_speeds = np.array(self.floris.flow_field.wind_speeds, copy=True) + wind_directions = np.array(self.floris.flow_field.wind_directions, copy=True) + + # Now set FLORIS wind speed and wind direction + # over to those values in the wind rose class + wind_speeds_wind_rose = wind_rose.df.ws.unique() + wind_directions_wind_rose = wind_rose.df.wd.unique() + self.reinitialize(wind_speeds=wind_speeds_wind_rose, wind_directions=wind_directions_wind_rose) + + # Build the frequency matrix from wind rose + freq = wind_rose.df.set_index(['wd','ws']).unstack().values + + + # Verify dimensions of the variable "freq" + if not ( + (np.shape(freq)[0] == self.floris.flow_field.n_wind_directions) + & (np.shape(freq)[1] == self.floris.flow_field.n_wind_speeds) + & (len(np.shape(freq)) == 2) + ): + raise UserWarning( + "'freq' should be a two-dimensional array with dimensions" + + " (n_wind_directions, n_wind_speeds)." + ) + + # Check if frequency vector sums to 1.0. If not, raise a warning + if np.abs(np.sum(freq) - 1.0) > 0.001: + self.logger.warning( + "WARNING: The frequency array provided to get_farm_AEP() " + + "does not sum to 1.0. " + ) + + # Copy the full wind speed array from the floris object and initialize + # the the farm_power variable as an empty array. + wind_speeds = np.array(self.floris.flow_field.wind_speeds, copy=True) + farm_power = np.zeros( + (self.floris.flow_field.n_wind_directions, len(wind_speeds)) + ) + + # Determine which wind speeds we must evaluate in floris + conditions_to_evaluate = (wind_speeds >= cut_in_wind_speed) + if cut_out_wind_speed is not None: + conditions_to_evaluate = conditions_to_evaluate & ( + wind_speeds < cut_out_wind_speed + ) + + # Evaluate the conditions in floris + if np.any(conditions_to_evaluate): + wind_speeds_subset = wind_speeds[conditions_to_evaluate] + yaw_angles_subset = None + if yaw_angles is not None: + yaw_angles_subset = yaw_angles[:, conditions_to_evaluate] + self.reinitialize(wind_speeds=wind_speeds_subset) + self.calculate_wake(yaw_angles=yaw_angles_subset, no_wake=no_wake) + farm_power[:, conditions_to_evaluate] = self.get_farm_power() + + # Finally, calculate AEP in GWh + aep = np.sum(np.multiply(freq, farm_power) * 365 * 24) + + # Reset the FLORIS object to the full wind speed and wind direction arrays + self.reinitialize(wind_speeds=wind_speeds, wind_directions=wind_directions) + + return aep + @property def layout_x(self): """ From 3a07bd7688dedf06f2c656a436f9b9ac75e9f9b6 Mon Sep 17 00:00:00 2001 From: Paul Date: Tue, 15 Mar 2022 11:15:37 -0600 Subject: [PATCH 03/27] add an example using wind rose object aep --- examples/07b_calc_aep_from_rose_use_class.py | 74 ++++++++++++++++++++ 1 file changed, 74 insertions(+) create mode 100644 examples/07b_calc_aep_from_rose_use_class.py diff --git a/examples/07b_calc_aep_from_rose_use_class.py b/examples/07b_calc_aep_from_rose_use_class.py new file mode 100644 index 000000000..358fbc19e --- /dev/null +++ b/examples/07b_calc_aep_from_rose_use_class.py @@ -0,0 +1,74 @@ +# 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 + + +import numpy as np +import pandas as pd +import matplotlib.pyplot as plt +from scipy.interpolate import NearestNDInterpolator +from floris.tools import FlorisInterface, WindRose, wind_rose + +""" +This example demonstrates how to calculate the Annual Energy Production (AEP) +of a wind farm using wind rose information stored in a .csv file. + +The wind rose information is first loaded, after which we initialize our Floris +Interface. A 3 turbine farm is generated, and then the turbine wakes and powers +are calculated across all the wind directions. Finally, the farm power is +converted to AEP and reported out. +""" + +# Read in the wind rose using the class +wind_rose = WindRose() +wind_rose.read_wind_rose_csv("inputs/wind_rose.csv") + +# Show the wind rose +wind_rose.plot_wind_rose() + +# Load the FLORIS object +fi = FlorisInterface("inputs/gch.yaml") # GCH model +# fi = FlorisInterface("inputs/cc.yaml") # CumulativeCurl model + +# Assume a three-turbine wind farm with 5D spacing. We reinitialize the +# floris object and assign the layout, wind speed and wind direction arrays. +D = 126.0 # Rotor diameter for the NREL 5 MW +fi.reinitialize( + layout=[[0.0, 5* D, 10 * D], [0.0, 0.0, 0.0]] +) + +# Compute the AEP using the default settings +aep = fi.get_farm_AEP_wind_rose_class(wind_rose=wind_rose) +print("Farm AEP (default options): {:.3f} GWh".format(aep / 1.0e9)) + +# Compute the AEP again while specifying a cut-in and cut-out wind speed. +# The wake calculations are skipped for any wind speed below respectively +# above the cut-in and cut-out wind speed. This can speed up computation and +# prevent unexpected behavior for zero/negative and very high wind speeds. +# In this example, the results should not change between this and the default +# call to 'get_farm_AEP()'. +aep = fi.get_farm_AEP_wind_rose_class( + wind_rose=wind_rose, + cut_in_wind_speed=3.0, # Wakes are not evaluated below this wind speed + cut_out_wind_speed=25.0, # Wakes are not evaluated above this wind speed +) +print("Farm AEP (with cut_in/out specified): {:.3f} GWh".format(aep / 1.0e9)) + +# Finally, we can also compute the AEP while ignoring all wake calculations. +# This can be useful to quantity the annual wake losses in the farm. Such +# calculations can be facilitated by enabling the 'no_wake' handle. +aep_no_wake = fi.get_farm_AEP_wind_rose_class(wind_rose=wind_rose, no_wake=True) +print("Farm AEP (no_wake=True): {:.3f} GWh".format(aep_no_wake / 1.0e9)) + + +plt.show() \ No newline at end of file From 1f040821b0e48dd63e56491013dfd3bf6467d9b3 Mon Sep 17 00:00:00 2001 From: Paul Date: Thu, 7 Apr 2022 12:14:50 -0600 Subject: [PATCH 04/27] update how no_wake is passed --- floris/tools/floris_interface.py | 5 ++++- 1 file changed, 4 insertions(+), 1 deletion(-) diff --git a/floris/tools/floris_interface.py b/floris/tools/floris_interface.py index e13b80121..5baac6f2b 100644 --- a/floris/tools/floris_interface.py +++ b/floris/tools/floris_interface.py @@ -813,7 +813,10 @@ def get_farm_AEP_wind_rose_class( if yaw_angles is not None: yaw_angles_subset = yaw_angles[:, conditions_to_evaluate] self.reinitialize(wind_speeds=wind_speeds_subset) - self.calculate_wake(yaw_angles=yaw_angles_subset, no_wake=no_wake) + if no_wake: + self.calculate_no_wake(yaw_angles=yaw_angles_subset) + else: + self.calculate_wake(yaw_angles=yaw_angles_subset) farm_power[:, conditions_to_evaluate] = self.get_farm_power() # Finally, calculate AEP in GWh From 15f30d5570a00261ae84c04ad759b7cef399a075 Mon Sep 17 00:00:00 2001 From: bayc Date: Wed, 27 Apr 2022 09:42:20 -0600 Subject: [PATCH 05/27] preventing negative sqrts that produce NANs in the cumulative model --- .../wake_velocity/cumulative_gauss_curl.py | 27 ++++++++++--------- 1 file changed, 15 insertions(+), 12 deletions(-) diff --git a/floris/simulation/wake_velocity/cumulative_gauss_curl.py b/floris/simulation/wake_velocity/cumulative_gauss_curl.py index 6a3f219e3..a55d2dd2f 100644 --- a/floris/simulation/wake_velocity/cumulative_gauss_curl.py +++ b/floris/simulation/wake_velocity/cumulative_gauss_curl.py @@ -175,21 +175,24 @@ def function( a2 = 2 ** (4 / n - 2) # based on Blondel model, modified to include cumulative effects - C = a1 - np.sqrt( - a2 - - ( - (n * turbine_Ct[:,:,ii:ii+1]) - * cosd(turbine_yaw) - / ( - 16.0 - * gamma(2 / n) - * np.sign(sigma_n) - * (np.abs(sigma_n) ** (4 / n)) - * (1 - sum_lbda) ** 2 - ) + tmp = a2 - ( + (n * turbine_Ct[:,:,ii:ii+1]) + * cosd(turbine_yaw) + / ( + 16.0 + * gamma(2 / n) + * np.sign(sigma_n) + * (np.abs(sigma_n) ** (4 / n)) + * (1 - sum_lbda) ** 2 ) ) + # for some low wind speeds, tmp can become slightly negative, which causes NANs, + # so replace the slightly negative values with zeros + tmp = tmp * np.array(tmp >= 0) + + C = a1 - np.sqrt(tmp) + C = C * (1 - sum_lbda) Ctmp[ii] = C From 8f6905304adb8a8b3643d1d097f68db2ebd00bdf Mon Sep 17 00:00:00 2001 From: bayc Date: Fri, 29 Apr 2022 14:47:23 -0600 Subject: [PATCH 06/27] add the base layout optimization class --- floris/tools/optimization/__init__.py | 2 +- .../layout_optimization/__init__.py | 0 .../layout_optimization_base.py | 83 +++++++++++++++++++ 3 files changed, 84 insertions(+), 1 deletion(-) create mode 100644 floris/tools/optimization/layout_optimization/__init__.py create mode 100644 floris/tools/optimization/layout_optimization/layout_optimization_base.py diff --git a/floris/tools/optimization/__init__.py b/floris/tools/optimization/__init__.py index f40bb816e..1ce9f032a 100644 --- a/floris/tools/optimization/__init__.py +++ b/floris/tools/optimization/__init__.py @@ -1 +1 @@ -from . import other, scipy, pyoptsparse, yaw_optimization +from . import other, scipy, pyoptsparse, yaw_optimization, layout_optimization diff --git a/floris/tools/optimization/layout_optimization/__init__.py b/floris/tools/optimization/layout_optimization/__init__.py new file mode 100644 index 000000000..e69de29bb diff --git a/floris/tools/optimization/layout_optimization/layout_optimization_base.py b/floris/tools/optimization/layout_optimization/layout_optimization_base.py new file mode 100644 index 000000000..7c71bb86b --- /dev/null +++ b/floris/tools/optimization/layout_optimization/layout_optimization_base.py @@ -0,0 +1,83 @@ +# 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 + +import numpy as np +import matplotlib.pyplot as plt +from shapely.geometry import Polygon, LineString + +from ....logging_manager import LoggerBase + +class LayoutOptimization(LoggerBase): + def __init__(self, fi, boundaries, min_dist=None, freq=None): + self.fi = fi + self.boundaries = boundaries + + self.boundary_polygon = Polygon(self.boundaries) + self.boundary_line = LineString(self.boundaries) + + 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]) + self.x0 = self._norm(self.fi.layout_x, self.xmin, self.xmax) + self.y0 = self._norm(self.fi.layout_y, self.ymin, self.ymax) + + if min_dist is None: + self.min_dist = 2 * self.rotor_diameter + else: + self.min_dist = min_dist + + if freq is None: + self.freq = 1 + else: + self.freq = freq + + self.wdir = self.fi.floris.flow_field.wind_directions + self.wspd = self.fi.floris.flow_field.wind_speeds + self.initial_AEP = np.sum(self.fi.get_farm_power() * self.freq * 8760) + + def __str__(self): + return "layout" + + def _norm(self, val, x1, x2): + return (val - x1) / (x2 - x1) + + def _unnorm(self, val, x1, x2): + return np.array(val) * (x2 - x1) + x1 + + # Public methods required for each optimization class + + def optimize(self): + sol = self._optimize() + return sol + + ########################################################################### + # Properties + ########################################################################### + + @property + def nturbs(self): + """ + This property returns the number of turbines in the FLORIS + object. + + Returns: + nturbs (int): The number of turbines in the FLORIS object. + """ + self._nturbs = self.fi.floris.farm.n_turbines + return self._nturbs + + @property + def rotor_diameter(self): + return self.fi.floris.farm.rotor_diameters[0][0][0] \ No newline at end of file From ff66aba07e01850ed89104327a9d33d38b0e3644 Mon Sep 17 00:00:00 2001 From: bayc Date: Fri, 29 Apr 2022 14:47:41 -0600 Subject: [PATCH 07/27] add the scipy layout optimzation class --- .../layout_optimization_scipy.py | 316 ++++++++++++++++++ 1 file changed, 316 insertions(+) create mode 100644 floris/tools/optimization/layout_optimization/layout_optimization_scipy.py diff --git a/floris/tools/optimization/layout_optimization/layout_optimization_scipy.py b/floris/tools/optimization/layout_optimization/layout_optimization_scipy.py new file mode 100644 index 000000000..c7d641290 --- /dev/null +++ b/floris/tools/optimization/layout_optimization/layout_optimization_scipy.py @@ -0,0 +1,316 @@ +# 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 + +import numpy as np +import matplotlib.pyplot as plt +from scipy.optimize import minimize +from shapely.geometry import Point +from scipy.spatial.distance import cdist + +from .layout_optimization_base import LayoutOptimization + +class LayoutOptimizationScipy(LayoutOptimization): + def __init__( + self, + fi, + boundaries, + freq=None, + x0=None, + bnds=None, + min_dist=None, + solver='SLSQP', + optOptions=None, + ): + super().__init__(fi, boundaries, min_dist=min_dist, freq=freq) + + if optOptions is None: + self.optOptions = {"maxiter": 100, "disp": True, "iprint": 2, "ftol": 1e-9, "eps":0.01} + + self.reinitialize_opt( + boundaries=boundaries, + freq=freq, + x0=x0, + bnds=bnds, + solver=solver, + optOptions=self.optOptions, + ) + + + # Private methods + + def _optimize(self): + self.residual_plant = minimize( + self._obj_func, + self.x0, + method=self.solver, + bounds=self.bnds, + constraints=self.cons, + options=self.optOptions, + ) + + return self.residual_plant.x + + def _obj_func(self, locs): + locs_unnorm = [ + self._unnorm(valx, self.bndx_min, self.bndx_max) + for valx in locs[0 : self.nturbs] + ] + [ + self._unnorm(valy, self.bndy_min, self.bndy_max) + for valy in locs[self.nturbs : 2 * self.nturbs] + ] + self._change_coordinates(locs_unnorm) + self.fi.calculate_wake() + AEP_sum = np.sum(self.fi.get_farm_power() * self.freq * 8760) + return -1 * AEP_sum / self.initial_AEP + + def _change_coordinates(self, locs): + # Parse the layout coordinates + layout_x = locs[0 : self.nturbs] + layout_y = locs[self.nturbs : 2 * self.nturbs] + layout_array = (layout_x, layout_y) + + # Update the turbine map in floris + self.fi.reinitialize(layout=layout_array) + + def _generate_constraints(self): + # grad_constraint1 = grad(self._space_constraint) + # grad_constraint2 = grad(self._distance_from_boundaries) + + tmp1 = { + "type": "ineq", + "fun": lambda x, *args: self._space_constraint(x), + } + tmp2 = { + "type": "ineq", + "fun": lambda x: self._distance_from_boundaries(x), + } + + self.cons = [tmp1, tmp2] + + def _set_opt_bounds(self): + self.bnds = [(0.0, 1.0) for _ in range(2 * self.nturbs)] + + def _space_constraint(self, x_in): + rho=500 + x = [ + self._unnorm(valx, self.bndx_min, self.bndx_max) + for valx in x_in[0 : self.nturbs] + ] + y = [ + self._unnorm(valy, self.bndy_min, self.bndy_max) + for valy in x_in[self.nturbs : 2 * self.nturbs] + ] + # x = np.array(x_in[0 : self.nturbs]) + # y = np.array(x_in[self.nturbs :]) + # Calculate distances between turbines + locs = np.vstack((x, y)).T + distances = cdist(locs, locs) + arange = np.arange(distances.shape[0]) + distances[arange, arange] = 1e10 + dist = np.min(distances, axis=0) + + g = 1 - np.array(dist) / self.min_dist + + # Following code copied from OpenMDAO KSComp(). + # Constraint is satisfied when KS_constraint <= 0 + g_max = np.max(np.atleast_2d(g), axis=-1)[:, np.newaxis] + g_diff = g - g_max + exponents = np.exp(rho * g_diff) + summation = np.sum(exponents, axis=-1)[:, np.newaxis] + KS_constraint = g_max + 1.0 / rho * np.log(summation) + + return -1*KS_constraint[0][0] + + def _distance_from_boundaries(self, x_in): + # x = np.array(x_in[0 : self.nturbs]) + # y = np.array(x_in[self.nturbs :]) + x = [ + self._unnorm(valx, self.bndx_min, self.bndx_max) + for valx in x_in[0 : self.nturbs] + ] + y = [ + self._unnorm(valy, self.bndy_min, self.bndy_max) + for valy in x_in[self.nturbs : 2 * self.nturbs] + ] + boundary_con = np.zeros(self.nturbs) + for i in range(self.nturbs): + loc = Point(x[i], y[i]) + boundary_con[i] = loc.distance(self.boundary_line) + if self.boundary_polygon.contains(loc)==True: + boundary_con[i] *= 1.0 + + return boundary_con + + # Public methods + + def optimize(self): + """ + This method finds the optimized layout of wind turbines for power + production given the provided frequencies of occurance of wind + conditions (wind speed, direction). + + Returns: + opt_locs (iterable): A list of the optimized locations of each + turbine (m). + """ + print("=====================================================") + print("Optimizing turbine layout...") + print("Number of parameters to optimize = ", len(self.x0)) + print("=====================================================") + + opt_locs_norm = self._optimize() + + print("Optimization complete.") + + opt_locs = [ + [ + self._unnorm(valx, self.bndx_min, self.bndx_max) + for valx in opt_locs_norm[0 : self.nturbs] + ], + [ + self._unnorm(valy, self.bndy_min, self.bndy_max) + for valy in opt_locs_norm[self.nturbs : 2 * self.nturbs] + ], + ] + + return opt_locs + + def reinitialize_opt( + self, + boundaries=None, + freq=None, + x0=None, + bnds=None, + solver=None, + optOptions=None, + ): + """ + This method reinitializes any optimization parameters that are + specified. Otherwise, the current parameter values are kept. + + Args: + boundaries (iterable(float, float)): Pairs of x- and y-coordinates + that represent the boundary's vertices (m). + wd (np.array): An array of wind directions (deg). Defaults to None. + ws (np.array): An array of wind speeds (m/s). Defaults to None. + freq (np.array): An array of the frequencies of occurance + correponding to each pair of wind direction and wind speed + values. Defaults to None. + AEP_initial (float): The initial Annual Energy + Production used for normalization in the optimization (Wh). If + not specified, initializes to the AEP of the current Floris + object. Defaults to None. + x0 (iterable, optional): The initial turbine locations, + ordered by x-coordinate and then y-coordiante (ie. [x1, x2, ... + , xn, y1, y2, ..., yn] (m)). If none are provided, x0 + initializes to the current turbine locations. 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 the min. and max. of the + boundaries iterable. 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. + opt_method (str, optional): The optimization method for + scipy.optimize.minize to use. If none is specified, initializes + to 'SLSQP'. Defaults to None. + opt_options (dict, optional): Dicitonary for setting the + optimization options. Defaults to None. + """ + if boundaries is not None: + self.boundaries = boundaries + self.bndx_min = np.min([val[0] for val in boundaries]) + self.bndy_min = np.min([val[1] for val in boundaries]) + self.bndx_max = np.max([val[0] for val in boundaries]) + self.bndy_max = np.max([val[1] for val in boundaries]) + self.boundaries_norm = [ + [ + self._norm(val[0], self.bndx_min, self.bndx_max), + self._norm(val[1], self.bndy_min, self.bndy_max), + ] + for val in self.boundaries + ] + if freq is not None: + self.freq = freq + if x0 is not None: + self.x0 = x0 + else: + self.x0 = [ + self._norm(x, self.bndx_min, self.bndx_max) + for x in self.fi.layout_x + ] + [ + self._norm(y, self.bndy_min, self.bndy_max) + for y in self.fi.layout_y + ] + if bnds is not None: + self.bnds = bnds + else: + self._set_opt_bounds() + if solver is not None: + self.solver = solver + if optOptions is not None: + self.optOptions = optOptions + + self._generate_constraints() + + def plot_layout_opt_results(self): + """ + This method plots the original and new locations of the turbines in a + wind farm after layout optimization. + """ + locsx_old = [ + self._unnorm(valx, self.bndx_min, self.bndx_max) + for valx in self.x0[0 : self.nturbs] + ] + locsy_old = [ + self._unnorm(valy, self.bndy_min, self.bndy_max) + for valy in self.x0[self.nturbs : 2 * self.nturbs] + ] + locsx = [ + self._unnorm(valx, self.bndx_min, self.bndx_max) + for valx in self.residual_plant.x[0 : self.nturbs] + ] + locsy = [ + self._unnorm(valy, self.bndy_min, self.bndy_max) + for valy in self.residual_plant.x[self.nturbs : 2 * self.nturbs] + ] + + plt.figure(figsize=(9, 6)) + fontsize = 16 + plt.plot(locsx_old, locsy_old, "ob") + plt.plot(locsx, locsy, "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"], + 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" + ) + + plt.show() From 21e73c2d48d1bd3ef34cc3a251256da3658869cf Mon Sep 17 00:00:00 2001 From: bayc Date: Fri, 29 Apr 2022 14:47:57 -0600 Subject: [PATCH 08/27] add the pyoptsparse layout optimization class --- .../layout_optimization_pyoptsparse.py | 200 ++++++++++++++++++ 1 file changed, 200 insertions(+) create mode 100644 floris/tools/optimization/layout_optimization/layout_optimization_pyoptsparse.py diff --git a/floris/tools/optimization/layout_optimization/layout_optimization_pyoptsparse.py b/floris/tools/optimization/layout_optimization/layout_optimization_pyoptsparse.py new file mode 100644 index 000000000..63ee31a6f --- /dev/null +++ b/floris/tools/optimization/layout_optimization/layout_optimization_pyoptsparse.py @@ -0,0 +1,200 @@ +# 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 + + +import numpy as np +import matplotlib.pyplot as plt +from shapely.geometry import Point +from scipy.spatial.distance import cdist + +from .layout_optimization_base import LayoutOptimization + +class LayoutOptimizationPyOptSparse(LayoutOptimization): + def __init__( + self, + fi, + boundaries, + min_dist=None, + freq=None, + solver=None, + optOptions=None, + ): + super().__init__(fi, boundaries, min_dist=min_dist, freq=freq) + self._reinitialize(solver=solver, optOptions=optOptions) + + def _reinitialize(self, solver=None, optOptions=None): + try: + import pyoptsparse + except ImportError: + err_msg = ( + "It appears you do not have pyOptSparse installed. " + + "Please refer to https://pyoptsparse.readthedocs.io/ for " + + "guidance on how to properly install the module." + ) + self.logger.error(err_msg, stack_info=True) + raise ImportError(err_msg) + + # Insantiate ptOptSparse optimization object with name and objective function + self.optProb = pyoptsparse.Optimization('layout', self._obj_func) + + self.optProb = self.add_var_group(self.optProb) + self.optProb = self.add_con_group(self.optProb) + self.optProb.addObj("obj") + + if solver is not None: + self.solver = solver + print("Setting up optimization with user's choice of solver: ", self.solver) + else: + self.solver = "SLSQP" + print("Setting up optimization with default solver: SLSQP.") + if optOptions is not None: + self.optOptions = optOptions + else: + if self.solver == "SNOPT": + self.optOptions = {"Major optimality tolerance": 1e-7} + else: + self.optOptions = {} + + exec("self.opt = pyoptsparse." + self.solver + "(options=self.optOptions)") + + def _optimize(self): + if hasattr(self, "_sens"): + self.sol = self.opt(self.optProb, sens=self._sens) + else: + self.sol = self.opt(self.optProb, sens="CDR", storeHistory='hist.hist') + return self.sol + + def _obj_func(self, varDict): + # Parse the variable dictionary + self.parse_opt_vars(varDict) + + # Update turbine map with turbince locations + self.fi.reinitialize(layout=[self.x, self.y]) + self.fi.calculate_wake() + + # Compute the objective function + funcs = {} + funcs["obj"] = ( + -1 * np.sum(self.fi.get_farm_power() * self.freq * 8760) / self.initial_AEP + ) + + # Compute constraints, if any are defined for the optimization + funcs = self.compute_cons(funcs, self.x, self.y) + + fail = False + return funcs, fail + + # Optionally, the user can supply the optimization with gradients + # def _sens(self, varDict, funcs): + # funcsSens = {} + # fail = False + # return funcsSens, fail + + def parse_opt_vars(self, varDict): + self.x = self._unnorm(varDict["x"], self.xmin, self.xmax) + self.y = self._unnorm(varDict["y"], self.ymin, self.ymax) + + def parse_sol_vars(self, sol): + self.x = list(self._unnorm(sol.getDVs()["x"], self.xmin, self.xmax))[0] + self.y = list(self._unnorm(sol.getDVs()["y"], self.ymin, self.ymax))[1] + + def add_var_group(self, optProb): + optProb.addVarGroup( + "x", self.nturbs, type="c", lower=0.0, upper=1.0, value=self.x0 + ) + optProb.addVarGroup( + "y", self.nturbs, type="c", lower=0.0, upper=1.0, value=self.y0 + ) + + return optProb + + def add_con_group(self, optProb): + optProb.addConGroup("boundary_con", self.nturbs, upper=0.0) + optProb.addConGroup("spacing_con", 1, upper=0.0) + + return optProb + + def compute_cons(self, funcs, x, y): + funcs["boundary_con"] = self.distance_from_boundaries(x, y) + funcs["spacing_con"] = self.space_constraint(x, y) + + return funcs + + def space_constraint(self, x, y, rho=500): + # Calculate distances between turbines + locs = np.vstack((x, y)).T + distances = cdist(locs, locs) + arange = np.arange(distances.shape[0]) + distances[arange, arange] = 1e10 + dist = np.min(distances, axis=0) + + g = 1 - np.array(dist) / self.min_dist + + # Following code copied from OpenMDAO KSComp(). + # Constraint is satisfied when KS_constraint <= 0 + g_max = np.max(np.atleast_2d(g), axis=-1)[:, np.newaxis] + g_diff = g - g_max + exponents = np.exp(rho * g_diff) + summation = np.sum(exponents, axis=-1)[:, np.newaxis] + KS_constraint = g_max + 1.0 / rho * np.log(summation) + + return KS_constraint[0][0] + + def distance_from_boundaries(self, x, y): + boundary_con = np.zeros(self.nturbs) + for i in range(self.nturbs): + loc = Point(x[i], y[i]) + boundary_con[i] = loc.distance(self.boundary_line) + if self.boundary_polygon.contains(loc)==True: + boundary_con[i] *= -1.0 + + return boundary_con + + def plot_layout_opt_results(self): + """ + Method to plot the old and new locations of the layout opitimization. + """ + locsx = self._unnorm(self.sol.getDVs()["x"], self.xmin, self.xmax) + locsy = self._unnorm(self.sol.getDVs()["y"], self.ymin, self.ymax) + x0 = self._unnorm(self.x0, self.xmin, self.xmax) + y0 = self._unnorm(self.y0, self.ymin, self.ymax) + + plt.figure(figsize=(9, 6)) + fontsize = 16 + plt.plot(x0, y0, "ob") + plt.plot(locsx, locsy, "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"], + 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" + ) + + plt.show() From 582d12bd4cd661af7e571e9b7bd4e3e5a1c1e3b1 Mon Sep 17 00:00:00 2001 From: Paul Date: Mon, 2 May 2022 12:32:42 -0600 Subject: [PATCH 09/27] Add label option --- floris/tools/wind_rose.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/floris/tools/wind_rose.py b/floris/tools/wind_rose.py index 245279085..e1e1ebe37 100644 --- a/floris/tools/wind_rose.py +++ b/floris/tools/wind_rose.py @@ -1299,7 +1299,7 @@ def indices_for_coord(self, f, lat_index, lon_index): ij = [int(round(x / 2000)) for x in delta] return tuple(reversed(ij)) - def plot_wind_speed_all(self, ax=None): + def plot_wind_speed_all(self, ax=None, label=None): """ This method plots the wind speed frequency distribution of the WindRose object averaged across all wind directions. If no axis is provided, a @@ -1313,7 +1313,7 @@ def plot_wind_speed_all(self, ax=None): _, ax = plt.subplots() df_plot = self.df.groupby("ws").sum() - ax.plot(self.ws, df_plot.freq_val) + ax.plot(self.ws, df_plot.freq_val, label=label) def plot_wind_speed_by_direction(self, dirs, ax=None): """ From 3b36d6ab4e95949af651d25c2b637217156c7b8d Mon Sep 17 00:00:00 2001 From: Paul Date: Mon, 2 May 2022 12:33:03 -0600 Subject: [PATCH 10/27] add x_20 turbine option --- floris/turbine_library/x_20MW.yaml | 176 +++++++++++++++++++++++++++++ 1 file changed, 176 insertions(+) create mode 100644 floris/turbine_library/x_20MW.yaml diff --git a/floris/turbine_library/x_20MW.yaml b/floris/turbine_library/x_20MW.yaml new file mode 100644 index 000000000..436a83b52 --- /dev/null +++ b/floris/turbine_library/x_20MW.yaml @@ -0,0 +1,176 @@ +turbine_type: 'x_20MW' +generator_efficiency: 1.0 +hub_height: 165.0 +pP: 1.88 +pT: 1.88 +rotor_diameter: 252.0 +TSR: 8.0 +power_thrust_table: + power: + - 0.000000 + - 0.000000 + - 0.074000 + - 0.325100 + - 0.376200 + - 0.402700 + - 0.415600 + - 0.423000 + - 0.427400 + - 0.429300 + - 0.429800 + - 0.429800 + - 0.429800 + - 0.429800 + - 0.429800 + - 0.429800 + - 0.429800 + - 0.429800 + - 0.429800 + - 0.429800 + - 0.429800 + - 0.429800 + - 0.429800 + - 0.429800 + - 0.429800 + - 0.429800 + - 0.429800 + - 0.429800 + - 0.429800 + - 0.429800 + - 0.429800 + - 0.429603 + - 0.354604 + - 0.316305 + - 0.281478 + - 0.250068 + - 0.221924 + - 0.196845 + - 0.174592 + - 0.154919 + - 0.137570 + - 0.122300 + - 0.108881 + - 0.097094 + - 0.086747 + - 0.077664 + - 0.069686 + - 0.062677 + - 0.056511 + - 0.051083 + - 0.046299 + - 0.043182 + - 0.033935 + - 0.000000 + - 0.000000 + thrust: + - 0.000000 + - 0.000000 + - 0.770100 + - 0.770100 + - 0.776300 + - 0.782400 + - 0.782000 + - 0.780200 + - 0.777200 + - 0.771900 + - 0.776800 + - 0.776800 + - 0.776800 + - 0.776800 + - 0.776800 + - 0.776800 + - 0.776800 + - 0.776800 + - 0.776800 + - 0.776800 + - 0.776800 + - 0.776800 + - 0.776800 + - 0.776800 + - 0.776800 + - 0.776800 + - 0.776800 + - 0.776800 + - 0.776800 + - 0.767500 + - 0.765100 + - 0.758700 + - 0.505600 + - 0.431000 + - 0.370800 + - 0.320900 + - 0.278800 + - 0.243200 + - 0.212800 + - 0.186800 + - 0.164500 + - 0.145400 + - 0.128900 + - 0.114700 + - 0.102400 + - 0.091800 + - 0.082500 + - 0.074500 + - 0.067500 + - 0.061300 + - 0.055900 + - 0.051200 + - 0.047000 + - 0.000000 + - 0.000000 + wind_speed: + - 0.000000 + - 2.900000 + - 3.000000 + - 4.000000 + - 4.514700 + - 5.000800 + - 5.457400 + - 5.883300 + - 6.277700 + - 6.639700 + - 6.968400 + - 7.263200 + - 7.523400 + - 7.748400 + - 7.937700 + - 8.090900 + - 8.207700 + - 8.287700 + - 8.330800 + - 8.337000 + - 8.367800 + - 8.435600 + - 8.540100 + - 8.681200 + - 8.858500 + - 9.071700 + - 9.320200 + - 9.603500 + - 9.921000 + - 10.272000 + - 10.655700 + - 11.507700 + - 12.267700 + - 12.744100 + - 13.249400 + - 13.782400 + - 14.342000 + - 14.926900 + - 15.535900 + - 16.167500 + - 16.820400 + - 17.493200 + - 18.184200 + - 18.892100 + - 19.615200 + - 20.351900 + - 21.100600 + - 21.859600 + - 22.627300 + - 23.401900 + - 24.181700 + - 24.750000 + - 25.010000 + - 25.020000 + - 50.000000 \ No newline at end of file From 84933206ba386135aab3d44605743763c1036e39 Mon Sep 17 00:00:00 2001 From: bayc Date: Tue, 10 May 2022 10:46:36 -0600 Subject: [PATCH 11/27] Add boundary grid functions for layout optimization --- .../layout_optimization_boundary_grid.py | 632 ++++++++++++++++++ 1 file changed, 632 insertions(+) create mode 100644 floris/tools/optimization/layout_optimization/layout_optimization_boundary_grid.py diff --git a/floris/tools/optimization/layout_optimization/layout_optimization_boundary_grid.py b/floris/tools/optimization/layout_optimization/layout_optimization_boundary_grid.py new file mode 100644 index 000000000..733aaaac3 --- /dev/null +++ b/floris/tools/optimization/layout_optimization/layout_optimization_boundary_grid.py @@ -0,0 +1,632 @@ +# 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 + + +import numpy as np +import matplotlib.pyplot as plt +from shapely.geometry import Point, Polygon, LineString +from scipy.spatial.distance import cdist + +from .layout_optimization_base import LayoutOptimization + +class LayoutOptimizationBoundaryGrid(LayoutOptimization): + def __init__( + self, + fi, + boundaries, + start, + x_spacing, + y_spacing, + shear, + rotation, + center_x, + center_y, + boundary_setback, + n_boundary_turbines=None, + boundary_spacing=None, + ): + self.fi = fi + + self.boundary_x = np.array([val[0] for val in boundaries]) + self.boundary_y = np.array([val[1] for val in boundaries]) + boundary = np.zeros((len(self.boundary_x), 2)) + boundary[:, 0] = self.boundary_x[:] + boundary[:, 1] = self.boundary_y[:] + self.boundary_polygon = Polygon(boundary) + + self.start = start + self.x_spacing = x_spacing + self.y_spacing = y_spacing + self.shear = shear + self.rotation = rotation + self.center_x = center_x + self.center_y = center_y + self.boundary_setback = boundary_setback + self.n_boundary_turbines = n_boundary_turbines + self.boundary_spacing = boundary_spacing + + def _discontinuous_grid( + self, + nrows, + ncols, + farm_width, + farm_height, + shear, + rotation, + center_x, + center_y, + shrink_boundary, + boundary_x, + boundary_y, + eps=1e-3, + ): + """ + Map from grid design variables to turbine x and y locations. Includes integer design variables and the formulation + results in a discontinous design space. + + TODO: shrink_boundary doesn't work well with concave boundaries, or with boundary angles less than 90 deg + + Args: + nrows (Int): number of rows in the grid. + ncols (Int): number of columns in the grid. + farm_width (Float): total grid width (before shear). + farm_height (Float): total grid height. + shear (Float): grid shear (rad). + rotation (Float): rotation about grid center (rad). + center_x (Float): location of grid x center. + center_y (Float): location of grid y center. + shrink_boundary (Float): how much to shrink the boundary that the grid can occupy. + boundary_x (Array(Float)): x boundary points. + boundary_y (Array(Float)): y boundary points. + + Returns: + grid_x (Array(Float)): turbine x locations. + grid_y (Array(Float)): turbine y locations. + """ + # create grid + nrows = int(nrows) + ncols = int(ncols) + xlocs = np.linspace(0.0, farm_width, ncols) + ylocs = np.linspace(0.0, farm_height, nrows) + y_spacing = ylocs[1] - ylocs[0] + nturbs = nrows * ncols + grid_x = np.zeros(nturbs) + grid_y = np.zeros(nturbs) + turb = 0 + for i in range(nrows): + for j in range(ncols): + grid_x[turb] = xlocs[j] + float(i) * y_spacing * np.tan(shear) + grid_y[turb] = ylocs[i] + turb += 1 + + # rotate + grid_x, grid_y = ( + np.cos(rotation) * grid_x - np.sin(rotation) * grid_y, + np.sin(rotation) * grid_x + np.cos(rotation) * grid_y, + ) + + # move center of grid + grid_x = (grid_x - np.mean(grid_x)) + center_x + grid_y = (grid_y - np.mean(grid_y)) + center_y + + # arrange the boundary + + # boundary = np.zeros((len(boundary_x),2)) + # boundary[:,0] = boundary_x[:] + # boundary[:,1] = boundary_y[:] + # poly = Polygon(boundary) + # centroid = poly.centroid + + # boundary[:,0] = (boundary_x[:]-centroid.x)*boundary_mult + centroid.x + # boundary[:,1] = (boundary_y[:]-centroid.y)*boundary_mult + centroid.y + # poly = Polygon(boundary) + + boundary = np.zeros((len(boundary_x), 2)) + boundary[:, 0] = boundary_x[:] + boundary[:, 1] = boundary_y[:] + poly = Polygon(boundary) + + if shrink_boundary != 0.0: + nBounds = len(boundary_x) + for i in range(nBounds): + point = Point(boundary_x[i] + eps, boundary_y[i]) + if poly.contains(point) is True or poly.touches(point) is True: + boundary[i, 0] = boundary_x[i] + shrink_boundary + else: + boundary[i, 0] = boundary_x[i] - shrink_boundary + + point = Point(boundary_x[i], boundary_y[i] + eps) + if poly.contains(point) is True or poly.touches(point) is True: + boundary[i, 1] = boundary_y[i] + shrink_boundary + else: + boundary[i, 1] = boundary_y[i] - shrink_boundary + + poly = Polygon(boundary) + + # get rid of points outside of boundary + index = 0 + for i in range(len(grid_x)): + point = Point(grid_x[index], grid_y[index]) + if poly.contains(point) is False and poly.touches(point) is False: + grid_x = np.delete(grid_x, index) + grid_y = np.delete(grid_y, index) + else: + index += 1 + + return grid_x, grid_y + + def _discrete_grid( + self, + x_spacing, + y_spacing, + shear, + rotation, + center_x, + center_y, + boundary_setback, + boundary_poly + ): + """ + returns grid turbine layout. Assumes the turbines fill the entire plant area + + Args: + x_spacing (Float): grid spacing in the unrotated x direction (m) + y_spacing (Float): grid spacing in the unrotated y direction (m) + shear (Float): grid shear (rad) + rotation (Float): grid rotation (rad) + center_x (Float): the x coordinate of the grid center (m) + center_y (Float): the y coordinate of the grid center (m) + boundary_poly (Polygon): a shapely Polygon of the wind plant boundary + + Returns + return_x (Array(Float)): turbine x locations + return_y (Array(Float)): turbine y locations + """ + + shrunk_poly = boundary_poly.buffer(-boundary_setback) + if shrunk_poly.area <= 0: + return np.array([]), np.array([]) + # create grid + minx, miny, maxx, maxy = shrunk_poly.bounds + width = maxx-minx + height = maxy-miny + + center_point = Point((center_x,center_y)) + poly_to_center = center_point.distance(shrunk_poly.centroid) + + width = np.max([width,poly_to_center]) + height = np.max([height,poly_to_center]) + nrows = int(np.max([width,height])/np.min([x_spacing,y_spacing]))*2 + 1 + ncols = nrows + + xlocs = np.arange(0,ncols)*x_spacing + ylocs = np.arange(0,nrows)*y_spacing + row_number = np.arange(0,nrows) + + d = np.array([i for x in xlocs for i in row_number]) + layout_x = np.array([x for x in xlocs for y in ylocs]) + d*y_spacing*np.tan(shear) + layout_y = np.array([y for x in xlocs for y in ylocs]) + + # rotate + rotate_x = np.cos(rotation)*layout_x - np.sin(rotation)*layout_y + rotate_y = np.sin(rotation)*layout_x + np.cos(rotation)*layout_y + + # move center of grid + rotate_x = (rotate_x - np.mean(rotate_x)) + center_x + rotate_y = (rotate_y - np.mean(rotate_y)) + center_y + + # get rid of points outside of boundary polygon + meets_constraints = np.zeros(len(rotate_x),dtype=bool) + for i in range(len(rotate_x)): + pt = Point(rotate_x[i],rotate_y[i]) + if shrunk_poly.contains(pt) or shrunk_poly.touches(pt): + meets_constraints[i] = True + + # arrange final x,y points + return_x = rotate_x[meets_constraints] + return_y = rotate_y[meets_constraints] + + return return_x, return_y + + def find_lengths(self, x, y, npoints): + length = np.zeros(len(x) - 1) + for i in range(npoints): + length[i] = np.sqrt((x[i + 1] - x[i]) ** 2 + (y[i + 1] - y[i]) ** 2) + return length + + # def _place_boundary_turbines(self, n_boundary_turbs, start, boundary_x, boundary_y): + # """ + # Place turbines equally spaced traversing the perimiter if the wind farm along the boundary + + # Args: + # n_boundary_turbs (Int): number of turbines to be placed on the boundary + # start (Float): where the first turbine should be placed + # boundary_x (Array(Float)): x boundary points + # boundary_y (Array(Float)): y boundary points + + # Returns + # layout_x (Array(Float)): turbine x locations + # layout_y (Array(Float)): turbine y locations + # """ + + # # check if the boundary is closed, correct if not + # if boundary_x[-1] != boundary_x[0] or boundary_y[-1] != boundary_y[0]: + # boundary_x = np.append(boundary_x, boundary_x[0]) + # boundary_y = np.append(boundary_y, boundary_y[0]) + + # # make the boundary + # boundary = np.zeros((len(boundary_x), 2)) + # boundary[:, 0] = boundary_x[:] + # boundary[:, 1] = boundary_y[:] + # poly = Polygon(boundary) + # perimeter = poly.length + + # # get the flattened turbine locations + # spacing = perimeter / float(n_boundary_turbs) + # flattened_locs = np.linspace(start, perimeter + start - spacing, n_boundary_turbs) + + # # set all of the flattened values between 0 and the perimeter + # for i in range(n_boundary_turbs): + # while flattened_locs[i] < 0.0: + # flattened_locs[i] += perimeter + # if flattened_locs[i] > perimeter: + # flattened_locs[i] = flattened_locs[i] % perimeter + + # # place the turbines around the perimeter + # nBounds = len(boundary_x) + # layout_x = np.zeros(n_boundary_turbs) + # layout_y = np.zeros(n_boundary_turbs) + + # lenBound = np.zeros(nBounds - 1) + # for i in range(nBounds - 1): + # lenBound[i] = Point(boundary[i]).distance(Point(boundary[i + 1])) + # for i in range(n_boundary_turbs): + # for j in range(nBounds - 1): + # if flattened_locs[i] < sum(lenBound[0 : j + 1]): + # layout_x[i] = ( + # boundary_x[j] + # + (boundary_x[j + 1] - boundary_x[j]) + # * (flattened_locs[i] - sum(lenBound[0:j])) + # / lenBound[j] + # ) + # layout_y[i] = ( + # boundary_y[j] + # + (boundary_y[j + 1] - boundary_y[j]) + # * (flattened_locs[i] - sum(lenBound[0:j])) + # / lenBound[j] + # ) + # break + + # return layout_x, layout_y + + def _place_boundary_turbines(self, start, boundary_poly, nturbs=None, spacing=None): + xBounds, yBounds = boundary_poly.boundary.coords.xy + print('5: ', nturbs) + print('6: ', spacing) + + if xBounds[-1] != xBounds[0]: + xBounds = np.append(xBounds, xBounds[0]) + yBounds = np.append(yBounds, yBounds[0]) + + nBounds = len(xBounds) + lenBound = self.find_lengths(xBounds, yBounds, len(xBounds) - 1) + circumference = sum(lenBound) + + if nturbs is not None and spacing is None: + # When the number of boundary turbines is specified + nturbs = int(nturbs) + bound_loc = np.linspace( + start, start + circumference - circumference / float(nturbs), nturbs + ) + elif spacing is not None and nturbs is None: + # When the spacing of boundary turbines is specified + nturbs = int(np.floor(circumference / spacing)) + bound_loc = np.linspace( + start, start + circumference - circumference / float(nturbs), nturbs + ) + else: + raise ValueError("Please specify either nturbs or spacing.") + + x = np.zeros(nturbs) + y = np.zeros(nturbs) + + if spacing is None: + # When the number of boundary turbines is specified + for i in range(nturbs): + if bound_loc[i] > circumference: + bound_loc[i] = bound_loc[i] % circumference + while bound_loc[i] < 0.0: + bound_loc[i] += circumference + for i in range(nturbs): + done = False + for j in range(nBounds): + if done == False: + if bound_loc[i] < sum(lenBound[0:j+1]): + point_x = xBounds[j] + (xBounds[j+1]-xBounds[j])*(bound_loc[i]-sum(lenBound[0:j]))/lenBound[j] + point_y = yBounds[j] + (yBounds[j+1]-yBounds[j])*(bound_loc[i]-sum(lenBound[0:j]))/lenBound[j] + done = True + x[i] = point_x + y[i] = point_y + else: + # When the spacing of boundary turbines is specified + additional_space = 0.0 + end_loop = False + for i in range(nturbs): + done = False + for j in range(nBounds): + while done == False: + dist = start + i*spacing + additional_space + if dist < sum(lenBound[0:j+1]): + point_x = xBounds[j] + (xBounds[j+1]-xBounds[j])*(dist -sum(lenBound[0:j]))/lenBound[j] + point_y = yBounds[j] + (yBounds[j+1]-yBounds[j])*(dist -sum(lenBound[0:j]))/lenBound[j] + + # Check if turbine is too close to previous turbine + if i > 0: + # Check if turbine just placed is to close to first turbine + min_dist = cdist([(point_x, point_y)], [(x[0], y[0])]) + if min_dist < spacing: + end_loop = True + ii = i + break + + min_dist = cdist([(point_x, point_y)], [(x[i-1], y[i-1])]) + if min_dist < spacing: + additional_space += 1.0 + else: + done = True + x[i] = point_x + y[i] = point_y + elif i == 0: + # If first turbine, just add initial turbine point + done = True + x[i] = point_x + y[i] = point_y + else: + pass + else: + break + if end_loop == True: + break + if end_loop == True: + x = x[:ii] + y = y[:ii] + break + print('xx: ', x) + print('yy: ', y) + return x, y + + def _place_boundary_turbines_with_specified_spacing(self, spacing, start, boundary_x, boundary_y): + """ + Place turbines equally spaced traversing the perimiter if the wind farm along the boundary + + Args: + n_boundary_turbs (Int): number of turbines to be placed on the boundary + start (Float): where the first turbine should be placed + boundary_x (Array(Float)): x boundary points + boundary_y (Array(Float)): y boundary points + + Returns + layout_x (Array(Float)): turbine x locations + layout_y (Array(Float)): turbine y locations + """ + + # check if the boundary is closed, correct if not + if boundary_x[-1] != boundary_x[0] or boundary_y[-1] != boundary_y[0]: + boundary_x = np.append(boundary_x, boundary_x[0]) + boundary_y = np.append(boundary_y, boundary_y[0]) + + # make the boundary + boundary = np.zeros((len(boundary_x), 2)) + boundary[:, 0] = boundary_x[:] + boundary[:, 1] = boundary_y[:] + poly = Polygon(boundary) + perimeter = poly.length + + # get the flattened turbine locations + n_boundary_turbs = int(perimeter / float(spacing)) + flattened_locs = np.linspace(start, perimeter + start - spacing, n_boundary_turbs) + + # set all of the flattened values between 0 and the perimeter + for i in range(n_boundary_turbs): + while flattened_locs[i] < 0.0: + flattened_locs[i] += perimeter + if flattened_locs[i] > perimeter: + flattened_locs[i] = flattened_locs[i] % perimeter + + # place the turbines around the perimeter + nBounds = len(boundary_x) + layout_x = np.zeros(n_boundary_turbs) + layout_y = np.zeros(n_boundary_turbs) + + lenBound = np.zeros(nBounds - 1) + for i in range(nBounds - 1): + lenBound[i] = Point(boundary[i]).distance(Point(boundary[i + 1])) + for i in range(n_boundary_turbs): + for j in range(nBounds - 1): + if flattened_locs[i] < sum(lenBound[0 : j + 1]): + layout_x[i] = ( + boundary_x[j] + + (boundary_x[j + 1] - boundary_x[j]) + * (flattened_locs[i] - sum(lenBound[0:j])) + / lenBound[j] + ) + layout_y[i] = ( + boundary_y[j] + + (boundary_y[j + 1] - boundary_y[j]) + * (flattened_locs[i] - sum(lenBound[0:j])) + / lenBound[j] + ) + break + + return layout_x, layout_y + + def boundary_grid( + self, + start, + x_spacing, + y_spacing, + shear, + rotation, + center_x, + center_y, + boundary_setback, + n_boundary_turbines=None, + boundary_spacing=None, + ): + """ + Place turbines equally spaced traversing the perimiter if the wind farm along the boundary + + Args: + n_boundary_turbs,start: boundary variables + nrows,ncols,farm_width,farm_height,shear,rotation,center_x,center_y,shrink_boundary,eps: grid variables + boundary_x,boundary_y: boundary points + + Returns + layout_x (Array(Float)): turbine x locations + layout_y (Array(Float)): turbine y locations + """ + + print('3: ', n_boundary_turbines) + print('4: ', boundary_spacing) + boundary_turbines_x, boundary_turbines_y = self._place_boundary_turbines( + start, self.boundary_polygon, nturbs=n_boundary_turbines, spacing=boundary_spacing + ) + # boundary_turbines_x, boundary_turbines_y = self._place_boundary_turbines_with_specified_spacing( + # spacing, start, boundary_x, boundary_y + # ) + + # grid_turbines_x, grid_turbines_y = self._discontinuous_grid( + # nrows, + # ncols, + # farm_width, + # farm_height, + # shear, + # rotation, + # center_x, + # center_y, + # shrink_boundary, + # boundary_x, + # boundary_y, + # eps=eps, + # ) + + grid_turbines_x, grid_turbines_y = self._discrete_grid( + x_spacing, + y_spacing, + shear, + rotation, + center_x, + center_y, + boundary_setback, + self.boundary_polygon, + ) + + layout_x = np.append(boundary_turbines_x, grid_turbines_x) + layout_y = np.append(boundary_turbines_y, grid_turbines_y) + + return layout_x, layout_y + + def reinitialize_bg( + self, + n_boundary_turbines=None, + start=None, + x_spacing=None, + y_spacing=None, + shear=None, + rotation=None, + center_x=None, + center_y=None, + boundary_setback=None, + boundary_x=None, + boundary_y=None, + boundary_spacing=None, + ): + + if n_boundary_turbines is not None: + self.n_boundary_turbines = n_boundary_turbines + if start is not None: + self.start = start + if x_spacing is not None: + self.x_spacing = x_spacing + if y_spacing is not None: + self.y_spacing = y_spacing + if shear is not None: + self.shear = shear + if rotation is not None: + self.rotation = rotation + if center_x is not None: + self.center_x = center_x + if center_y is not None: + self.center_y = center_y + if boundary_setback is not None: + self.boundary_setback = boundary_setback + if boundary_x is not None: + self.boundary_x = boundary_x + if boundary_y is not None: + self.boundary_y = boundary_y + if boundary_spacing is not None: + self.boundary_spacing = boundary_spacing + + def reinitialize_xy(self): + print('1: ', self.n_boundary_turbines) + print('2: ', self.boundary_spacing) + layout_x, layout_y = self.boundary_grid( + self.start, + self.x_spacing, + self.y_spacing, + self.shear, + self.rotation, + self.center_x, + self.center_y, + self.boundary_setback, + self.n_boundary_turbines, + self.boundary_spacing, + ) + + self.fi.reinitialize(layout=(layout_x, layout_y)) + + def plot_layout(self): + plt.figure(figsize=(9, 6)) + fontsize = 16 + + plt.plot(self.fi.layout_x, self.fi.layout_y, "ob") + # plt.plot(locsx, locsy, "or") + + 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.show() + + def space_constraint(self, x, y, min_dist, rho=500): + # Calculate distances between turbines + locs = np.vstack((x, y)).T + distances = cdist(locs, locs) + arange = np.arange(distances.shape[0]) + distances[arange, arange] = 1e10 + dist = np.min(distances, axis=0) + + g = 1 - np.array(dist) / min_dist + + # Following code copied from OpenMDAO KSComp(). + # Constraint is satisfied when KS_constraint <= 0 + g_max = np.max(np.atleast_2d(g), axis=-1)[:, np.newaxis] + g_diff = g - g_max + exponents = np.exp(rho * g_diff) + summation = np.sum(exponents, axis=-1)[:, np.newaxis] + KS_constraint = g_max + 1.0 / rho * np.log(summation) + + return KS_constraint[0][0], dist From 02b425d629aaeb9380b55b0413780c359a6c95a7 Mon Sep 17 00:00:00 2001 From: bayc Date: Thu, 12 May 2022 17:46:22 -0600 Subject: [PATCH 12/27] removing print statements --- .../layout_optimization_boundary_grid.py | 8 -------- 1 file changed, 8 deletions(-) diff --git a/floris/tools/optimization/layout_optimization/layout_optimization_boundary_grid.py b/floris/tools/optimization/layout_optimization/layout_optimization_boundary_grid.py index 733aaaac3..08ed0198d 100644 --- a/floris/tools/optimization/layout_optimization/layout_optimization_boundary_grid.py +++ b/floris/tools/optimization/layout_optimization/layout_optimization_boundary_grid.py @@ -312,8 +312,6 @@ def find_lengths(self, x, y, npoints): def _place_boundary_turbines(self, start, boundary_poly, nturbs=None, spacing=None): xBounds, yBounds = boundary_poly.boundary.coords.xy - print('5: ', nturbs) - print('6: ', spacing) if xBounds[-1] != xBounds[0]: xBounds = np.append(xBounds, xBounds[0]) @@ -402,8 +400,6 @@ def _place_boundary_turbines(self, start, boundary_poly, nturbs=None, spacing=No x = x[:ii] y = y[:ii] break - print('xx: ', x) - print('yy: ', y) return x, y def _place_boundary_turbines_with_specified_spacing(self, spacing, start, boundary_x, boundary_y): @@ -497,8 +493,6 @@ def boundary_grid( layout_y (Array(Float)): turbine y locations """ - print('3: ', n_boundary_turbines) - print('4: ', boundary_spacing) boundary_turbines_x, boundary_turbines_y = self._place_boundary_turbines( start, self.boundary_polygon, nturbs=n_boundary_turbines, spacing=boundary_spacing ) @@ -579,8 +573,6 @@ def reinitialize_bg( self.boundary_spacing = boundary_spacing def reinitialize_xy(self): - print('1: ', self.n_boundary_turbines) - print('2: ', self.boundary_spacing) layout_x, layout_y = self.boundary_grid( self.start, self.x_spacing, From c7004a8e692640a746785d4675dafc363e54e79a Mon Sep 17 00:00:00 2001 From: bayc Date: Thu, 12 May 2022 18:33:10 -0600 Subject: [PATCH 13/27] adding check for second turbine --- .../layout_optimization_boundary_grid.py | 10 +++++++--- 1 file changed, 7 insertions(+), 3 deletions(-) diff --git a/floris/tools/optimization/layout_optimization/layout_optimization_boundary_grid.py b/floris/tools/optimization/layout_optimization/layout_optimization_boundary_grid.py index 08ed0198d..8a483b672 100644 --- a/floris/tools/optimization/layout_optimization/layout_optimization_boundary_grid.py +++ b/floris/tools/optimization/layout_optimization/layout_optimization_boundary_grid.py @@ -374,9 +374,13 @@ def _place_boundary_turbines(self, start, boundary_poly, nturbs=None, spacing=No # Check if turbine just placed is to close to first turbine min_dist = cdist([(point_x, point_y)], [(x[0], y[0])]) if min_dist < spacing: - end_loop = True - ii = i - break + # TODO: make this more robust; pass is needed if 2nd turbine is too close to the first + if i == 1: + pass + else: + end_loop = True + ii = i + break min_dist = cdist([(point_x, point_y)], [(x[i-1], y[i-1])]) if min_dist < spacing: From 95cdf5466ce1d36ebde5811232c9955a7fabe98a Mon Sep 17 00:00:00 2001 From: Paul Fleming Date: Thu, 26 May 2022 10:21:22 -0600 Subject: [PATCH 14/27] add time limit and storeHistory options --- .../layout_optimization_pyoptsparse.py | 11 ++++++++++- 1 file changed, 10 insertions(+), 1 deletion(-) diff --git a/floris/tools/optimization/layout_optimization/layout_optimization_pyoptsparse.py b/floris/tools/optimization/layout_optimization/layout_optimization_pyoptsparse.py index 63ee31a6f..1d1669190 100644 --- a/floris/tools/optimization/layout_optimization/layout_optimization_pyoptsparse.py +++ b/floris/tools/optimization/layout_optimization/layout_optimization_pyoptsparse.py @@ -29,10 +29,15 @@ def __init__( freq=None, solver=None, optOptions=None, + timeLimit=None, + storeHistory='hist.hist' ): super().__init__(fi, boundaries, min_dist=min_dist, freq=freq) self._reinitialize(solver=solver, optOptions=optOptions) + self.storeHistory = storeHistory + self.timeLimit = timeLimit + def _reinitialize(self, solver=None, optOptions=None): try: import pyoptsparse @@ -72,7 +77,11 @@ def _optimize(self): if hasattr(self, "_sens"): self.sol = self.opt(self.optProb, sens=self._sens) else: - self.sol = self.opt(self.optProb, sens="CDR", storeHistory='hist.hist') + if self.timeLimit is not None: + print(self.timeLimit) + self.sol = self.opt(self.optProb, sens="CDR", storeHistory=self.storeHistory, timeLimit=self.timeLimit) + else: + self.sol = self.opt(self.optProb, sens="CDR", storeHistory=self.storeHistory) return self.sol def _obj_func(self, varDict): From 18f746b120390328a9dd5494e609aa3c7c938cb8 Mon Sep 17 00:00:00 2001 From: Paul Date: Tue, 31 May 2022 15:23:21 -0600 Subject: [PATCH 15/27] add hotStart option --- .../layout_optimization_pyoptsparse.py | 9 +++++---- 1 file changed, 5 insertions(+), 4 deletions(-) diff --git a/floris/tools/optimization/layout_optimization/layout_optimization_pyoptsparse.py b/floris/tools/optimization/layout_optimization/layout_optimization_pyoptsparse.py index 1d1669190..701d7fffc 100644 --- a/floris/tools/optimization/layout_optimization/layout_optimization_pyoptsparse.py +++ b/floris/tools/optimization/layout_optimization/layout_optimization_pyoptsparse.py @@ -30,13 +30,15 @@ def __init__( solver=None, optOptions=None, timeLimit=None, - storeHistory='hist.hist' + storeHistory='hist.hist', + hotStart=None ): super().__init__(fi, boundaries, min_dist=min_dist, freq=freq) self._reinitialize(solver=solver, optOptions=optOptions) self.storeHistory = storeHistory self.timeLimit = timeLimit + self.hotStart = hotStart def _reinitialize(self, solver=None, optOptions=None): try: @@ -78,10 +80,9 @@ def _optimize(self): self.sol = self.opt(self.optProb, sens=self._sens) else: if self.timeLimit is not None: - print(self.timeLimit) - self.sol = self.opt(self.optProb, sens="CDR", storeHistory=self.storeHistory, timeLimit=self.timeLimit) + self.sol = self.opt(self.optProb, sens="CDR", storeHistory=self.storeHistory, timeLimit=self.timeLimit, hotStart=self.hotStart) else: - self.sol = self.opt(self.optProb, sens="CDR", storeHistory=self.storeHistory) + self.sol = self.opt(self.optProb, sens="CDR", storeHistory=self.storeHistory, hotStart=self.hotStart) return self.sol def _obj_func(self, varDict): From b113df08664326aaff91966c18ba7183324c9e50 Mon Sep 17 00:00:00 2001 From: Paul Date: Wed, 1 Jun 2022 15:51:11 -0600 Subject: [PATCH 16/27] call get_farm_aep from within the wind rose ver --- floris/tools/floris_interface.py | 59 ++++++-------------------------- 1 file changed, 11 insertions(+), 48 deletions(-) diff --git a/floris/tools/floris_interface.py b/floris/tools/floris_interface.py index 5baac6f2b..d85ff8056 100644 --- a/floris/tools/floris_interface.py +++ b/floris/tools/floris_interface.py @@ -773,60 +773,23 @@ def get_farm_AEP_wind_rose_class( # Build the frequency matrix from wind rose freq = wind_rose.df.set_index(['wd','ws']).unstack().values + # Now compute aep + aep = self.get_farm_AEP( + freq, + cut_in_wind_speed=cut_in_wind_speed, + cut_out_wind_speed=cut_out_wind_speed, + yaw_angles=yaw_angles, + no_wake=no_wake) - # Verify dimensions of the variable "freq" - if not ( - (np.shape(freq)[0] == self.floris.flow_field.n_wind_directions) - & (np.shape(freq)[1] == self.floris.flow_field.n_wind_speeds) - & (len(np.shape(freq)) == 2) - ): - raise UserWarning( - "'freq' should be a two-dimensional array with dimensions" - + " (n_wind_directions, n_wind_speeds)." - ) - # Check if frequency vector sums to 1.0. If not, raise a warning - if np.abs(np.sum(freq) - 1.0) > 0.001: - self.logger.warning( - "WARNING: The frequency array provided to get_farm_AEP() " - + "does not sum to 1.0. " - ) - - # Copy the full wind speed array from the floris object and initialize - # the the farm_power variable as an empty array. - wind_speeds = np.array(self.floris.flow_field.wind_speeds, copy=True) - farm_power = np.zeros( - (self.floris.flow_field.n_wind_directions, len(wind_speeds)) - ) - - # Determine which wind speeds we must evaluate in floris - conditions_to_evaluate = (wind_speeds >= cut_in_wind_speed) - if cut_out_wind_speed is not None: - conditions_to_evaluate = conditions_to_evaluate & ( - wind_speeds < cut_out_wind_speed - ) - - # Evaluate the conditions in floris - if np.any(conditions_to_evaluate): - wind_speeds_subset = wind_speeds[conditions_to_evaluate] - yaw_angles_subset = None - if yaw_angles is not None: - yaw_angles_subset = yaw_angles[:, conditions_to_evaluate] - self.reinitialize(wind_speeds=wind_speeds_subset) - if no_wake: - self.calculate_no_wake(yaw_angles=yaw_angles_subset) - else: - self.calculate_wake(yaw_angles=yaw_angles_subset) - farm_power[:, conditions_to_evaluate] = self.get_farm_power() - - # Finally, calculate AEP in GWh - aep = np.sum(np.multiply(freq, farm_power) * 365 * 24) - - # Reset the FLORIS object to the full wind speed and wind direction arrays + # Reset the FLORIS object to the original wind speed and directions self.reinitialize(wind_speeds=wind_speeds, wind_directions=wind_directions) + return aep + + @property def layout_x(self): """ From 36227391785ef29d7aa25e1223c82ba20ae23a2f Mon Sep 17 00:00:00 2001 From: bayc Date: Fri, 3 Jun 2022 15:31:23 -0600 Subject: [PATCH 17/27] add deflection to the turbopark model --- floris/simulation/solver.py | 43 ++++++++++++++++---- floris/simulation/wake_velocity/turbopark.py | 5 ++- floris/tools/floris_interface.py | 6 --- 3 files changed, 37 insertions(+), 17 deletions(-) diff --git a/floris/simulation/solver.py b/floris/simulation/solver.py index d25d18755..4e126fd06 100644 --- a/floris/simulation/solver.py +++ b/floris/simulation/solver.py @@ -702,6 +702,7 @@ def turbopark_solver(farm: Farm, flow_field: FlowField, grid: TurbineGrid, model w_wake = np.zeros_like(flow_field.w_initial_sorted) shape = (farm.n_turbines,) + np.shape(flow_field.u_initial_sorted) velocity_deficit = np.zeros(shape) + deflection_field = np.zeros_like(flow_field.u_initial_sorted) turbine_turbulence_intensity = flow_field.turbulence_intensity * np.ones((flow_field.n_wind_directions, flow_field.n_wind_speeds, farm.n_turbines, 1, 1)) ambient_turbulence_intensity = flow_field.turbulence_intensity @@ -768,15 +769,38 @@ def turbopark_solver(farm: Farm, flow_field: FlowField, grid: TurbineGrid, model # Model calculations # NOTE: exponential - deflection_field = model_manager.deflection_model.function( - x_i, - y_i, - effective_yaw_i, - turbulence_intensity_i, - ct_i, - rotor_diameter_i, - **deflection_model_args - ) + if (farm.yaw_angles_sorted == 0).all(): + pass + else: + for ii in range(i): + x_ii = np.mean(grid.x_sorted[:, :, ii:ii+1], axis=(3, 4)) + x_ii = x_ii[:, :, :, None, None] + y_ii = np.mean(grid.y_sorted[:, :, ii:ii+1], axis=(3, 4)) + y_ii = y_ii[:, :, :, None, None] + + yaw_ii = farm.yaw_angles_sorted[:, :, ii:ii+1, None, None] + turbulence_intensity_ii = turbine_turbulence_intensity[:, :, ii:ii+1] + ct_ii = Ct( + velocities=flow_field.u_sorted, + yaw_angle=farm.yaw_angles_sorted, + fCt=farm.turbine_fCts, + turbine_type_map=farm.turbine_type_map_sorted, + ix_filter=[ii] + ) + ct_ii = ct_ii[:, :, 0:1, None, None] + rotor_diameter_ii = farm.rotor_diameters_sorted[: ,:, ii:ii+1, None, None] + + deflection_field_ii = model_manager.deflection_model.function( + x_ii, + y_ii, + yaw_ii, + turbulence_intensity_ii, + ct_ii, + rotor_diameter_ii, + **deflection_model_args + ) + + deflection_field[:,:,ii:ii+1,:,:] = deflection_field_ii[:,:,i:i+1,:,:] if model_manager.enable_transverse_velocities: v_wake, w_wake = calculate_transverse_velocity( @@ -815,6 +839,7 @@ def turbopark_solver(farm: Farm, flow_field: FlowField, grid: TurbineGrid, model rotor_diameter_i, farm.rotor_diameters_sorted[:, :, :, None, None], i, + deflection_field, **deficit_model_args ) diff --git a/floris/simulation/wake_velocity/turbopark.py b/floris/simulation/wake_velocity/turbopark.py index 20bf8ffa2..d16382cdf 100644 --- a/floris/simulation/wake_velocity/turbopark.py +++ b/floris/simulation/wake_velocity/turbopark.py @@ -72,6 +72,7 @@ def function( rotor_diameter_i: np.ndarray, rotor_diameters: np.ndarray, i: int, + deflection_field: np.ndarray, # enforces the use of the below as keyword arguments and adherence to the # unpacking of the results from prepare_function() *, @@ -89,8 +90,8 @@ def function( x_dist = (x_i - x) * downstream_mask / rotor_diameters # Radial distance between turbine i and the centerlines of wakes from all real/image turbines - r_dist = np.sqrt((y_i - y) ** 2 + (z_i - z) ** 2) - r_dist_image = np.sqrt((y_i - y) ** 2 + (z_i - (-z)) ** 2) + r_dist = np.sqrt((y_i - (y + deflection_field)) ** 2 + (z_i - z) ** 2) + r_dist_image = np.sqrt((y_i - (y + deflection_field)) ** 2 + (z_i - (-z)) ** 2) Cts[:,:,i:,:,:] = 0.00001 diff --git a/floris/tools/floris_interface.py b/floris/tools/floris_interface.py index 4f3b33723..e7ebb9c3a 100644 --- a/floris/tools/floris_interface.py +++ b/floris/tools/floris_interface.py @@ -113,9 +113,6 @@ def calculate_wake( # TODO decide where to handle this sign issue if (yaw_angles is not None) and not (np.all(yaw_angles==0.)): - if self.floris.wake.model_strings["velocity_model"] == "turbopark": - # TODO: Implement wake steering for the TurbOPark model - raise ValueError("Non-zero yaw angles given and for TurbOPark model; wake steering with this model is not yet implemented.") self.floris.farm.yaw_angles = yaw_angles # Initialize solution space @@ -142,9 +139,6 @@ def calculate_no_wake( # TODO decide where to handle this sign issue if (yaw_angles is not None) and not (np.all(yaw_angles==0.)): - if self.floris.wake.model_strings["velocity_model"] == "turbopark": - # TODO: Implement wake steering for the TurbOPark model - raise ValueError("Non-zero yaw angles given and for TurbOPark model; wake steering with this model is not yet implemented.") self.floris.farm.yaw_angles = yaw_angles # Initialize solution space From 60a5ffcb12714c72f9936e19fbdc39da2cf2f04f Mon Sep 17 00:00:00 2001 From: Paul Date: Fri, 3 Jun 2022 21:39:25 -0600 Subject: [PATCH 18/27] Add a simplified spreading optimization --- .../layout_optimization_pyoptsparse_spread.py | 218 ++++++++++++++++++ 1 file changed, 218 insertions(+) create mode 100644 floris/tools/optimization/layout_optimization/layout_optimization_pyoptsparse_spread.py diff --git a/floris/tools/optimization/layout_optimization/layout_optimization_pyoptsparse_spread.py b/floris/tools/optimization/layout_optimization/layout_optimization_pyoptsparse_spread.py new file mode 100644 index 000000000..772fa0fab --- /dev/null +++ b/floris/tools/optimization/layout_optimization/layout_optimization_pyoptsparse_spread.py @@ -0,0 +1,218 @@ +# 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 + + +import numpy as np +import matplotlib.pyplot as plt +from shapely.geometry import Point +from scipy.spatial.distance import cdist + +from .layout_optimization_base import LayoutOptimization + +class LayoutOptimizationPyOptSparse(LayoutOptimization): + def __init__( + self, + fi, + boundaries, + 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) + self._reinitialize(solver=solver, optOptions=optOptions) + + self.storeHistory = storeHistory + self.timeLimit = timeLimit + self.hotStart = hotStart + + def _reinitialize(self, solver=None, optOptions=None): + try: + import pyoptsparse + except ImportError: + err_msg = ( + "It appears you do not have pyOptSparse installed. " + + "Please refer to https://pyoptsparse.readthedocs.io/ for " + + "guidance on how to properly install the module." + ) + self.logger.error(err_msg, stack_info=True) + raise ImportError(err_msg) + + # Insantiate ptOptSparse optimization object with name and objective function + self.optProb = pyoptsparse.Optimization('layout', self._obj_func) + + self.optProb = self.add_var_group(self.optProb) + self.optProb = self.add_con_group(self.optProb) + self.optProb.addObj("obj") + + if solver is not None: + self.solver = solver + print("Setting up optimization with user's choice of solver: ", self.solver) + else: + self.solver = "SLSQP" + print("Setting up optimization with default solver: SLSQP.") + if optOptions is not None: + self.optOptions = optOptions + else: + if self.solver == "SNOPT": + self.optOptions = {"Major optimality tolerance": 1e-7} + else: + self.optOptions = {} + + exec("self.opt = pyoptsparse." + self.solver + "(options=self.optOptions)") + + def _optimize(self): + if hasattr(self, "_sens"): + self.sol = self.opt(self.optProb, sens=self._sens) + else: + if self.timeLimit is not None: + self.sol = self.opt(self.optProb, sens="CDR", storeHistory=self.storeHistory, timeLimit=self.timeLimit, hotStart=self.hotStart) + else: + self.sol = self.opt(self.optProb, sens="CDR", storeHistory=self.storeHistory, hotStart=self.hotStart) + return self.sol + + def _obj_func(self, varDict): + # Parse the variable dictionary + self.parse_opt_vars(varDict) + + # Update turbine map with turbince locations + # self.fi.reinitialize(layout=[self.x, self.y]) + # self.fi.calculate_wake() + + # Compute the objective function + 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 + funcs = self.compute_cons(funcs, self.x, self.y) + + fail = False + return funcs, fail + + # Optionally, the user can supply the optimization with gradients + # def _sens(self, varDict, funcs): + # funcsSens = {} + # fail = False + # return funcsSens, fail + + def parse_opt_vars(self, varDict): + self.x = self._unnorm(varDict["x"], self.xmin, self.xmax) + self.y = self._unnorm(varDict["y"], self.ymin, self.ymax) + + def parse_sol_vars(self, sol): + self.x = list(self._unnorm(sol.getDVs()["x"], self.xmin, self.xmax))[0] + self.y = list(self._unnorm(sol.getDVs()["y"], self.ymin, self.ymax))[1] + + def add_var_group(self, optProb): + optProb.addVarGroup( + "x", self.nturbs, type="c", lower=0.0, upper=1.0, value=self.x0 + ) + optProb.addVarGroup( + "y", self.nturbs, type="c", lower=0.0, upper=1.0, value=self.y0 + ) + + return optProb + + def add_con_group(self, optProb): + optProb.addConGroup("boundary_con", self.nturbs, upper=0.0) + optProb.addConGroup("spacing_con", 1, upper=0.0) + + return optProb + + def compute_cons(self, funcs, x, y): + funcs["boundary_con"] = self.distance_from_boundaries(x, y) + funcs["spacing_con"] = self.space_constraint(x, y) + + return funcs + + def mean_distance(self, x, y): + + locs = np.vstack((x, y)).T + distances = cdist(locs, locs) + return np.mean(distances) + + + def space_constraint(self, x, y, rho=500): + # Calculate distances between turbines + locs = np.vstack((x, y)).T + distances = cdist(locs, locs) + arange = np.arange(distances.shape[0]) + distances[arange, arange] = 1e10 + dist = np.min(distances, axis=0) + + g = 1 - np.array(dist) / self.min_dist + + # Following code copied from OpenMDAO KSComp(). + # Constraint is satisfied when KS_constraint <= 0 + g_max = np.max(np.atleast_2d(g), axis=-1)[:, np.newaxis] + g_diff = g - g_max + exponents = np.exp(rho * g_diff) + summation = np.sum(exponents, axis=-1)[:, np.newaxis] + KS_constraint = g_max + 1.0 / rho * np.log(summation) + + return KS_constraint[0][0] + + def distance_from_boundaries(self, x, y): + boundary_con = np.zeros(self.nturbs) + for i in range(self.nturbs): + loc = Point(x[i], y[i]) + boundary_con[i] = loc.distance(self.boundary_line) + if self.boundary_polygon.contains(loc)==True: + boundary_con[i] *= -1.0 + + return boundary_con + + def plot_layout_opt_results(self): + """ + Method to plot the old and new locations of the layout opitimization. + """ + locsx = self._unnorm(self.sol.getDVs()["x"], self.xmin, self.xmax) + locsy = self._unnorm(self.sol.getDVs()["y"], self.ymin, self.ymax) + x0 = self._unnorm(self.x0, self.xmin, self.xmax) + y0 = self._unnorm(self.y0, self.ymin, self.ymax) + + plt.figure(figsize=(9, 6)) + fontsize = 16 + plt.plot(x0, y0, "ob") + plt.plot(locsx, locsy, "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"], + 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" + ) + + plt.show() From 894ca22b63a965206d9b42f69f76a8feca77f5ad Mon Sep 17 00:00:00 2001 From: Michael LoCascio Date: Fri, 29 Jul 2022 09:31:57 -0700 Subject: [PATCH 19/27] changes to layout optimization with flowers --- floris/simulation/turbine.py | 4 +++- .../layout_optimization_base.py | 2 +- .../layout_optimization_pyoptsparse.py | 12 ++++++++---- .../optimization/pyoptsparse/optimization.py | 19 +++++++++++++------ 4 files changed, 25 insertions(+), 12 deletions(-) diff --git a/floris/simulation/turbine.py b/floris/simulation/turbine.py index ef52e281a..9dffbb511 100644 --- a/floris/simulation/turbine.py +++ b/floris/simulation/turbine.py @@ -375,7 +375,9 @@ def __attrs_post_init__(self) -> None: inner_power = 0.5 * self.rotor_area * self.fCp_interp(wind_speeds) * self.generator_efficiency * wind_speeds ** 3 self.power_interp = interp1d( wind_speeds, - inner_power + inner_power, + fill_value = 0.0, + bounds_error=False, ) """ diff --git a/floris/tools/optimization/layout_optimization/layout_optimization_base.py b/floris/tools/optimization/layout_optimization/layout_optimization_base.py index 7c71bb86b..f079300e0 100644 --- a/floris/tools/optimization/layout_optimization/layout_optimization_base.py +++ b/floris/tools/optimization/layout_optimization/layout_optimization_base.py @@ -80,4 +80,4 @@ def nturbs(self): @property def rotor_diameter(self): - return self.fi.floris.farm.rotor_diameters[0][0][0] \ No newline at end of file + return self.fi.floris.farm.rotor_diameters_sorted[0][0][0] \ No newline at end of file diff --git a/floris/tools/optimization/layout_optimization/layout_optimization_pyoptsparse.py b/floris/tools/optimization/layout_optimization/layout_optimization_pyoptsparse.py index 701d7fffc..80d8d8399 100644 --- a/floris/tools/optimization/layout_optimization/layout_optimization_pyoptsparse.py +++ b/floris/tools/optimization/layout_optimization/layout_optimization_pyoptsparse.py @@ -25,6 +25,8 @@ def __init__( self, fi, boundaries, + scale_dv=1.0, + scale_con=1.0, min_dist=None, freq=None, solver=None, @@ -34,6 +36,8 @@ def __init__( hotStart=None ): super().__init__(fi, boundaries, min_dist=min_dist, freq=freq) + self.scale_dv = scale_dv + self.scale_con = scale_con self._reinitialize(solver=solver, optOptions=optOptions) self.storeHistory = storeHistory @@ -121,17 +125,17 @@ def parse_sol_vars(self, sol): def add_var_group(self, optProb): optProb.addVarGroup( - "x", self.nturbs, type="c", lower=0.0, upper=1.0, value=self.x0 + "x", self.nturbs, type="c", lower=0.0, upper=1.0, value=self.x0, scale=self.scale_dv ) optProb.addVarGroup( - "y", self.nturbs, type="c", lower=0.0, upper=1.0, value=self.y0 + "y", self.nturbs, type="c", lower=0.0, upper=1.0, value=self.y0, scale=self.scale_dv ) return optProb def add_con_group(self, optProb): - optProb.addConGroup("boundary_con", self.nturbs, upper=0.0) - optProb.addConGroup("spacing_con", 1, upper=0.0) + optProb.addConGroup("boundary_con", self.nturbs, upper=0.0, scale=self.scale_con) + optProb.addConGroup("spacing_con", 1, upper=0.0, scale=self.scale_con) return optProb diff --git a/floris/tools/optimization/pyoptsparse/optimization.py b/floris/tools/optimization/pyoptsparse/optimization.py index 65e6c2a49..ecf884d9c 100644 --- a/floris/tools/optimization/pyoptsparse/optimization.py +++ b/floris/tools/optimization/pyoptsparse/optimization.py @@ -27,7 +27,7 @@ class Optimization(LoggerBase): Optimization: An instantiated Optimization object. """ - def __init__(self, model, solver=None): + def __init__(self, model, solver=None, storeHistory='hist.hist', optOptions=None, timeLimit=None): """ Instantiate Optimization object and its parameters. """ @@ -51,7 +51,11 @@ def __init__(self, model, solver=None): + str(self.solver_choices) ) - self.reinitialize(solver=solver) + self.optOptions = optOptions + self.reinitialize(solver=solver, optOptions=optOptions) + self.storeHistory = storeHistory + self.timeLimit = timeLimit + # Private methods @@ -85,7 +89,7 @@ def _reinitialize(self, solver=None, optOptions=None): if self.solver == "SNOPT": self.optOptions = {"Major optimality tolerance": 1e-7} else: - self.optOptions = {} + self.optOptions = {"IPRINT": 0} exec("self.opt = pyoptsparse." + self.solver + "(options=self.optOptions)") @@ -93,12 +97,15 @@ def _optimize(self): if hasattr(self.model, "_sens"): self.sol = self.opt(self.optProb, sens=self.model._sens) else: - self.sol = self.opt(self.optProb, sens="CDR", storeHistory='hist.hist') + if self.timeLimit is not None: + self.sol = self.opt(self.optProb, sens="CDR", storeHistory=self.storeHistory, timeLimit=self.timeLimit) + else: + self.sol = self.opt(self.optProb, sens="CDR", storeHistory=self.storeHistory) # Public methods - def reinitialize(self, solver=None): - self._reinitialize(solver=solver) + def reinitialize(self, solver=None, optOptions=None): + self._reinitialize(solver=solver, optOptions=optOptions) def optimize(self): self._optimize() From 1d94a3ad9fcaa49fb0df8b018d631102f35f6eb4 Mon Sep 17 00:00:00 2001 From: Michael LoCascio Date: Fri, 29 Jul 2022 09:35:24 -0700 Subject: [PATCH 20/27] testing new branch --- .../layout_optimization/layout_optimization_pyoptsparse.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/floris/tools/optimization/layout_optimization/layout_optimization_pyoptsparse.py b/floris/tools/optimization/layout_optimization/layout_optimization_pyoptsparse.py index 80d8d8399..f4d2c12d2 100644 --- a/floris/tools/optimization/layout_optimization/layout_optimization_pyoptsparse.py +++ b/floris/tools/optimization/layout_optimization/layout_optimization_pyoptsparse.py @@ -12,7 +12,7 @@ # See https://floris.readthedocs.io for documentation - +# Test import numpy as np import matplotlib.pyplot as plt from shapely.geometry import Point From 13730cd1662644dd4c8d8428195d3cd75600f8ea Mon Sep 17 00:00:00 2001 From: Michael LoCascio Date: Fri, 29 Jul 2022 10:48:27 -0600 Subject: [PATCH 21/27] testing changes --- .../layout_optimization/layout_optimization_pyoptsparse.py | 1 - 1 file changed, 1 deletion(-) diff --git a/floris/tools/optimization/layout_optimization/layout_optimization_pyoptsparse.py b/floris/tools/optimization/layout_optimization/layout_optimization_pyoptsparse.py index f4d2c12d2..71e75b248 100644 --- a/floris/tools/optimization/layout_optimization/layout_optimization_pyoptsparse.py +++ b/floris/tools/optimization/layout_optimization/layout_optimization_pyoptsparse.py @@ -12,7 +12,6 @@ # See https://floris.readthedocs.io for documentation -# Test import numpy as np import matplotlib.pyplot as plt from shapely.geometry import Point From 1e3dc3284255733c76ebb11f3b548f2444173758 Mon Sep 17 00:00:00 2001 From: Michael LoCascio Date: Fri, 29 Jul 2022 12:06:38 -0600 Subject: [PATCH 22/27] remote branch testing --- .../layout_optimization/layout_optimization_pyoptsparse.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/floris/tools/optimization/layout_optimization/layout_optimization_pyoptsparse.py b/floris/tools/optimization/layout_optimization/layout_optimization_pyoptsparse.py index 71e75b248..72b6453be 100644 --- a/floris/tools/optimization/layout_optimization/layout_optimization_pyoptsparse.py +++ b/floris/tools/optimization/layout_optimization/layout_optimization_pyoptsparse.py @@ -11,7 +11,7 @@ # the License. # See https://floris.readthedocs.io for documentation - +#Test import numpy as np import matplotlib.pyplot as plt from shapely.geometry import Point From 49af1b6fd1cb9adb2bb544ebb2c83e9bf7eef810 Mon Sep 17 00:00:00 2001 From: Michael LoCascio Date: Fri, 29 Jul 2022 11:07:00 -0700 Subject: [PATCH 23/27] test again --- .../layout_optimization/layout_optimization_pyoptsparse.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/floris/tools/optimization/layout_optimization/layout_optimization_pyoptsparse.py b/floris/tools/optimization/layout_optimization/layout_optimization_pyoptsparse.py index 72b6453be..71e75b248 100644 --- a/floris/tools/optimization/layout_optimization/layout_optimization_pyoptsparse.py +++ b/floris/tools/optimization/layout_optimization/layout_optimization_pyoptsparse.py @@ -11,7 +11,7 @@ # the License. # See https://floris.readthedocs.io for documentation -#Test + import numpy as np import matplotlib.pyplot as plt from shapely.geometry import Point From 1e8cf244ff3d735308cc967b592cba05f61922ac Mon Sep 17 00:00:00 2001 From: Michael LoCascio Date: Mon, 1 Aug 2022 07:47:20 -0700 Subject: [PATCH 24/27] snopt finite difference --- .../layout_optimization/layout_optimization_pyoptsparse.py | 6 +++--- floris/tools/optimization/pyoptsparse/optimization.py | 6 +++--- 2 files changed, 6 insertions(+), 6 deletions(-) diff --git a/floris/tools/optimization/layout_optimization/layout_optimization_pyoptsparse.py b/floris/tools/optimization/layout_optimization/layout_optimization_pyoptsparse.py index 71e75b248..a95a3142f 100644 --- a/floris/tools/optimization/layout_optimization/layout_optimization_pyoptsparse.py +++ b/floris/tools/optimization/layout_optimization/layout_optimization_pyoptsparse.py @@ -82,10 +82,10 @@ def _optimize(self): if hasattr(self, "_sens"): self.sol = self.opt(self.optProb, sens=self._sens) else: - if self.timeLimit is not None: - self.sol = self.opt(self.optProb, sens="CDR", storeHistory=self.storeHistory, timeLimit=self.timeLimit, hotStart=self.hotStart) + if self.timeLimit is not None: #sens="CDR", + self.sol = self.opt(self.optProb, storeHistory=self.storeHistory, timeLimit=self.timeLimit, hotStart=self.hotStart) else: - self.sol = self.opt(self.optProb, sens="CDR", storeHistory=self.storeHistory, hotStart=self.hotStart) + self.sol = self.opt(self.optProb, storeHistory=self.storeHistory, hotStart=self.hotStart) return self.sol def _obj_func(self, varDict): diff --git a/floris/tools/optimization/pyoptsparse/optimization.py b/floris/tools/optimization/pyoptsparse/optimization.py index ecf884d9c..5087e7aa9 100644 --- a/floris/tools/optimization/pyoptsparse/optimization.py +++ b/floris/tools/optimization/pyoptsparse/optimization.py @@ -97,10 +97,10 @@ def _optimize(self): if hasattr(self.model, "_sens"): self.sol = self.opt(self.optProb, sens=self.model._sens) else: - if self.timeLimit is not None: - self.sol = self.opt(self.optProb, sens="CDR", storeHistory=self.storeHistory, timeLimit=self.timeLimit) + if self.timeLimit is not None: #sens="CDR" + self.sol = self.opt(self.optProb, storeHistory=self.storeHistory, timeLimit=self.timeLimit) else: - self.sol = self.opt(self.optProb, sens="CDR", storeHistory=self.storeHistory) + self.sol = self.opt(self.optProb, storeHistory=self.storeHistory) # Public methods From 6b61549f1e61ada61dde67452dbcc51b30560a14 Mon Sep 17 00:00:00 2001 From: Michael LoCascio Date: Mon, 8 Aug 2022 13:34:04 -0600 Subject: [PATCH 25/27] change to layout optimization --- .../layout_optimization/layout_optimization_pyoptsparse.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/floris/tools/optimization/layout_optimization/layout_optimization_pyoptsparse.py b/floris/tools/optimization/layout_optimization/layout_optimization_pyoptsparse.py index 71e75b248..7b68924f8 100644 --- a/floris/tools/optimization/layout_optimization/layout_optimization_pyoptsparse.py +++ b/floris/tools/optimization/layout_optimization/layout_optimization_pyoptsparse.py @@ -133,7 +133,7 @@ def add_var_group(self, optProb): return optProb def add_con_group(self, optProb): - optProb.addConGroup("boundary_con", self.nturbs, upper=0.0, scale=self.scale_con) + optProb.addConGroup("boundary_con", self.nturbs, upper=-1.0, scale=self.scale_con) optProb.addConGroup("spacing_con", 1, upper=0.0, scale=self.scale_con) return optProb From 7da56ae1c7a4dd3f7f6674088eda4dee1dcbc6ad Mon Sep 17 00:00:00 2001 From: Michael LoCascio Date: Wed, 15 Mar 2023 18:05:06 -0700 Subject: [PATCH 26/27] turbine model with constant coefficients --- floris/turbine_library/nrel_5MW_cc.yaml | 164 ++++++++++++++++++++++++ 1 file changed, 164 insertions(+) create mode 100644 floris/turbine_library/nrel_5MW_cc.yaml diff --git a/floris/turbine_library/nrel_5MW_cc.yaml b/floris/turbine_library/nrel_5MW_cc.yaml new file mode 100644 index 000000000..242ed21ae --- /dev/null +++ b/floris/turbine_library/nrel_5MW_cc.yaml @@ -0,0 +1,164 @@ +turbine_type: 'nrel_5MW' +generator_efficiency: 1.0 +hub_height: 90.0 +pP: 1.88 +pT: 1.88 +rotor_diameter: 126.0 +TSR: 8.0 +power_thrust_table: + power: + - 0.0 + - 0.000000 + - 0.000000 + - 0.43 + - 0.43 + - 0.43 + - 0.43 + - 0.43 + - 0.43 + - 0.43 + - 0.43 + - 0.43 + - 0.43 + - 0.43 + - 0.43 + - 0.43 + - 0.43 + - 0.43 + - 0.43 + - 0.43 + - 0.43 + - 0.43 + - 0.43 + - 0.43 + - 0.43 + - 0.43 + - 0.43 + - 0.43 + - 0.43 + - 0.43 + - 0.43 + - 0.43 + - 0.43 + - 0.43 + - 0.43 + - 0.43 + - 0.43 + - 0.43 + - 0.43 + - 0.43 + - 0.43 + - 0.43 + - 0.43 + - 0.43 + - 0.43 + - 0.43 + - 0.43 + - 0.43 + - 0.43 + - 0.0 + - 0.0 + thrust: + - 0.0 + - 0.0 + - 0.0 + - 0.75 + - 0.75 + - 0.75 + - 0.75 + - 0.75 + - 0.75 + - 0.75 + - 0.75 + - 0.75 + - 0.75 + - 0.75 + - 0.75 + - 0.75 + - 0.75 + - 0.75 + - 0.75 + - 0.75 + - 0.75 + - 0.75 + - 0.75 + - 0.75 + - 0.75 + - 0.75 + - 0.75 + - 0.75 + - 0.75 + - 0.75 + - 0.75 + - 0.75 + - 0.75 + - 0.75 + - 0.75 + - 0.75 + - 0.75 + - 0.75 + - 0.75 + - 0.75 + - 0.75 + - 0.75 + - 0.75 + - 0.75 + - 0.75 + - 0.75 + - 0.75 + - 0.75 + - 0.75 + - 0.0 + - 0.0 + wind_speed: + - 0.0 + - 2.0 + - 2.5 + - 3.0 + - 3.5 + - 4.0 + - 4.5 + - 5.0 + - 5.5 + - 6.0 + - 6.5 + - 7.0 + - 7.5 + - 8.0 + - 8.5 + - 9.0 + - 9.5 + - 10.0 + - 10.5 + - 11.0 + - 11.5 + - 12.0 + - 12.5 + - 13.0 + - 13.5 + - 14.0 + - 14.5 + - 15.0 + - 15.5 + - 16.0 + - 16.5 + - 17.0 + - 17.5 + - 18.0 + - 18.5 + - 19.0 + - 19.5 + - 20.0 + - 20.5 + - 21.0 + - 21.5 + - 22.0 + - 22.5 + - 23.0 + - 23.5 + - 24.0 + - 24.5 + - 25.0 + - 25.01 + - 25.02 + - 50.0 From 2183858c53b2cdc8574278944235e9fb4d20f967 Mon Sep 17 00:00:00 2001 From: Michael LoCascio Date: Fri, 17 Mar 2023 09:51:07 -0600 Subject: [PATCH 27/27] layout optimization changes --- .../layout_optimization_pyoptsparse.py | 16 ++++++++-------- 1 file changed, 8 insertions(+), 8 deletions(-) diff --git a/floris/tools/optimization/layout_optimization/layout_optimization_pyoptsparse.py b/floris/tools/optimization/layout_optimization/layout_optimization_pyoptsparse.py index 158e7577d..011d3ba5f 100644 --- a/floris/tools/optimization/layout_optimization/layout_optimization_pyoptsparse.py +++ b/floris/tools/optimization/layout_optimization/layout_optimization_pyoptsparse.py @@ -133,8 +133,8 @@ def add_var_group(self, optProb): return optProb def add_con_group(self, optProb): - optProb.addConGroup("boundary_con", self.nturbs, upper=-1.0, scale=self.scale_con) - optProb.addConGroup("spacing_con", 1, upper=0.0, scale=self.scale_con) + optProb.addConGroup("boundary_con", self.nturbs, upper=0.0, scale=self.scale_con) + optProb.addConGroup("spacing_con", self.nturbs, upper=0.0, scale=self.scale_con) return optProb @@ -156,13 +156,13 @@ def space_constraint(self, x, y, rho=500): # Following code copied from OpenMDAO KSComp(). # Constraint is satisfied when KS_constraint <= 0 - g_max = np.max(np.atleast_2d(g), axis=-1)[:, np.newaxis] - g_diff = g - g_max - exponents = np.exp(rho * g_diff) - summation = np.sum(exponents, axis=-1)[:, np.newaxis] - KS_constraint = g_max + 1.0 / rho * np.log(summation) + # g_max = np.max(np.atleast_2d(g), axis=-1)[:, np.newaxis] + # g_diff = g - g_max + # exponents = np.exp(rho * g_diff) + # summation = np.sum(exponents, axis=-1)[:, np.newaxis] + # KS_constraint = g_max + 1.0 / rho * np.log(summation) - return KS_constraint[0][0] + return g def distance_from_boundaries(self, x, y): boundary_con = np.zeros(self.nturbs)