From 724f452a97a2a848a842307df63e4ef430f04a4c Mon Sep 17 00:00:00 2001 From: misi9170 <39596329+misi9170@users.noreply.github.com> Date: Fri, 22 Mar 2024 17:09:48 -0400 Subject: [PATCH] Save WindData onto FlorisModel and simplify post-`run()` calls (#849) * Add functions to reshape turbine and farm power to wd x ws * Add tests of new wind rose functions * wind_data saved onto FlorisModel; functions partially built. * Update tests; rename _wind_data to wind_data since optimizers may need access. * 07 example updated (waked and no_wake match previous output). * Remove wind_data need from layout optimizers. * Bugfix; copy did not bring over wind_data." * Update examples 13, 15 * Removing unneeded methods. * Group hidden and outer get_turbine_powers methods. * Ruff and isort. * Add getter for wind_data. * Rename converters to be more explicit about output. * Updating tests. * Copy up docstring from hidden version. * Fix a couple more examples. * Fix scaling on uniform frequency. * ruff and isort. * Log warnings when freq not provided and test. * Update uncertain model for new paradigm * Add test of wind rose setting * bufgix * change _expanded suffix to _rose to avoid confusion in UncertainFlorisModel. * to_ methods specify class rather than suggested instantiation. --------- Co-authored-by: Paul --- examples/07_calc_aep_from_rose.py | 23 +- .../13_optimize_yaw_with_neighboring_farm.py | 3 + examples/15_optimize_layout.py | 6 +- .../16c_optimize_layout_with_heterogeneity.py | 13 +- examples/29_floating_vs_fixedbottom_farm.py | 1 + examples/34_wind_data.py | 10 +- floris/floris_model.py | 316 +++++++++------- .../layout_optimization_base.py | 25 +- .../layout_optimization_pyoptsparse.py | 11 +- .../layout_optimization_pyoptsparse_spread.py | 3 +- .../layout_optimization_scipy.py | 15 +- floris/uncertain_floris_model.py | 355 ++++++++++-------- floris/wind_data.py | 23 +- tests/floris_model_integration_test.py | 156 +++++--- tests/layout_optimization_integration_test.py | 43 ++- .../parallel_floris_model_integration_test.py | 3 + ...uncertain_floris_model_integration_test.py | 49 ++- tests/wind_data_integration_test.py | 18 +- 18 files changed, 625 insertions(+), 448 deletions(-) diff --git a/examples/07_calc_aep_from_rose.py b/examples/07_calc_aep_from_rose.py index cc2de88d4..135a4c119 100644 --- a/examples/07_calc_aep_from_rose.py +++ b/examples/07_calc_aep_from_rose.py @@ -55,26 +55,15 @@ wind_speeds=wind_speeds, turbulence_intensities=turbulence_intensities, ) +fmodel.run() # Compute the AEP using the default settings aep = fmodel.get_farm_AEP(freq=freq) -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 = fmodel.get_farm_AEP( - freq=freq, - 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)) +print("Farm AEP: {:.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 = fmodel.get_farm_AEP(freq, no_wake=True) -print("Farm AEP (no_wake=True): {:.3f} GWh".format(aep_no_wake / 1.0e9)) +# calculations can be facilitated by first running with run_no_wake(). +fmodel.run_no_wake() +aep_no_wake = fmodel.get_farm_AEP(freq=freq) +print("Farm AEP (no wakes): {:.3f} GWh".format(aep_no_wake / 1.0e9)) diff --git a/examples/13_optimize_yaw_with_neighboring_farm.py b/examples/13_optimize_yaw_with_neighboring_farm.py index 18d5e1b26..300748341 100644 --- a/examples/13_optimize_yaw_with_neighboring_farm.py +++ b/examples/13_optimize_yaw_with_neighboring_farm.py @@ -202,6 +202,7 @@ def yaw_opt_interpolant(wd, ws): print(" ") print("===========================================================") print("Calculating baseline annual energy production (AEP)...") + fmodel_aep.run() aep_bl_subset = 1.0e-9 * fmodel_aep.get_farm_AEP( freq=freq_windrose, turbine_weights=turbine_weights @@ -247,11 +248,13 @@ def yaw_opt_interpolant(wd, ws): print("===========================================================") print("Calculating annual energy production with wake steering (AEP)...") fmodel_aep.set(yaw_angles=yaw_angles_opt_nonb_AEP) + fmodel_aep.run() aep_opt_subset_nonb = 1.0e-9 * fmodel_aep.get_farm_AEP( freq=freq_windrose, turbine_weights=turbine_weights, ) fmodel_aep.set(yaw_angles=yaw_angles_opt_AEP) + fmodel_aep.run() aep_opt_subset = 1.0e-9 * fmodel_aep.get_farm_AEP( freq=freq_windrose, turbine_weights=turbine_weights, diff --git a/examples/15_optimize_layout.py b/examples/15_optimize_layout.py index 071a62b87..df0f1d460 100644 --- a/examples/15_optimize_layout.py +++ b/examples/15_optimize_layout.py @@ -54,7 +54,7 @@ fmodel.set(layout_x=layout_x, layout_y=layout_y) # Setup the optimization problem -layout_opt = LayoutOptimizationScipy(fmodel, boundaries, wind_data=wind_rose) +layout_opt = LayoutOptimizationScipy(fmodel, boundaries) # Run the optimization sol = layout_opt.optimize() @@ -62,10 +62,10 @@ # Get the resulting improvement in AEP print('... calcuating improvement in AEP') fmodel.run() -base_aep = fmodel.get_farm_AEP_with_wind_data(wind_data=wind_rose) / 1e6 +base_aep = fmodel.get_farm_AEP() / 1e6 fmodel.set(layout_x=sol[0], layout_y=sol[1]) fmodel.run() -opt_aep = fmodel.get_farm_AEP_with_wind_data(wind_data=wind_rose) / 1e6 +opt_aep = fmodel.get_farm_AEP() / 1e6 percent_gain = 100 * (opt_aep - base_aep) / base_aep diff --git a/examples/16c_optimize_layout_with_heterogeneity.py b/examples/16c_optimize_layout_with_heterogeneity.py index 616b60e68..069511cd8 100644 --- a/examples/16c_optimize_layout_with_heterogeneity.py +++ b/examples/16c_optimize_layout_with_heterogeneity.py @@ -87,7 +87,6 @@ layout_opt = LayoutOptimizationScipy( fmodel, boundaries, - wind_data=wind_rose, min_dist=2*D, optOptions={"maxiter":maxiter} ) @@ -100,10 +99,10 @@ print('... calcuating improvement in AEP') fmodel.run() -base_aep = fmodel.get_farm_AEP_with_wind_data(wind_data=wind_rose) / 1e6 +base_aep = fmodel.get_farm_AEP() / 1e6 fmodel.set(layout_x=sol[0], layout_y=sol[1]) fmodel.run() -opt_aep = fmodel.get_farm_AEP_with_wind_data(wind_data=wind_rose) / 1e6 +opt_aep = fmodel.get_farm_AEP() / 1e6 percent_gain = 100 * (opt_aep - base_aep) / base_aep @@ -128,7 +127,6 @@ layout_opt = LayoutOptimizationScipy( fmodel, boundaries, - wind_data=wind_rose, min_dist=2*D, enable_geometric_yaw=True, optOptions={"maxiter":maxiter} @@ -142,12 +140,11 @@ print('... calcuating improvement in AEP') fmodel.set(yaw_angles=np.zeros_like(layout_opt.yaw_angles)) -base_aep = fmodel.get_farm_AEP_with_wind_data(wind_data=wind_rose) / 1e6 +fmodel.run() +base_aep = fmodel.get_farm_AEP() / 1e6 fmodel.set(layout_x=sol[0], layout_y=sol[1], yaw_angles=layout_opt.yaw_angles) fmodel.run() -opt_aep = fmodel.get_farm_AEP_with_wind_data( - wind_data=wind_rose -) / 1e6 +opt_aep = fmodel.get_farm_AEP() / 1e6 percent_gain = 100 * (opt_aep - base_aep) / base_aep diff --git a/examples/29_floating_vs_fixedbottom_farm.py b/examples/29_floating_vs_fixedbottom_farm.py index e04ac3f98..ef9745621 100644 --- a/examples/29_floating_vs_fixedbottom_farm.py +++ b/examples/29_floating_vs_fixedbottom_farm.py @@ -124,6 +124,7 @@ wind_speeds= ws_grid.flatten(), turbulence_intensities=0.06 * np.ones_like(wd_grid.flatten()) ) + fmodel.run() # Compute the AEP aep_fixed = fmodel_fixed.get_farm_AEP(freq=freq) diff --git a/examples/34_wind_data.py b/examples/34_wind_data.py index 3a4d56fe5..0d17e7924 100644 --- a/examples/34_wind_data.py +++ b/examples/34_wind_data.py @@ -39,7 +39,7 @@ time_series = TimeSeries(wd_array, ws_array, turbulence_intensities=ti_array) # Now build the wind rose -wind_rose = time_series.to_wind_rose() +wind_rose = time_series.to_WindRose() # Plot the wind rose fig, ax = plt.subplots(subplot_kw={"polar": True}) @@ -47,7 +47,7 @@ fig.suptitle("WindRose Plot") # Now build a wind rose with turbulence intensity -wind_ti_rose = time_series.to_wind_ti_rose() +wind_ti_rose = time_series.to_WindTIRose() # Plot the wind rose with TI fig, axs = plt.subplots(2, 1, figsize=(6,8), subplot_kw={"polar": True}) @@ -78,9 +78,9 @@ wind_rose_power = fmodel_wind_rose.get_farm_power() wind_ti_rose_power = fmodel_wind_ti_rose.get_farm_power() -time_series_aep = fmodel_time_series.get_farm_AEP_with_wind_data(time_series) -wind_rose_aep = fmodel_wind_rose.get_farm_AEP_with_wind_data(wind_rose) -wind_ti_rose_aep = fmodel_wind_ti_rose.get_farm_AEP_with_wind_data(wind_ti_rose) +time_series_aep = fmodel_time_series.get_farm_AEP() +wind_rose_aep = fmodel_wind_rose.get_farm_AEP() +wind_ti_rose_aep = fmodel_wind_ti_rose.get_farm_AEP() print(f"AEP from TimeSeries {time_series_aep / 1e9:.2f} GWh") print(f"AEP from WindRose {wind_rose_aep / 1e9:.2f} GWh") diff --git a/floris/floris_model.py b/floris/floris_model.py index 2b0f6cb9a..548f2e9f6 100644 --- a/floris/floris_model.py +++ b/floris/floris_model.py @@ -35,7 +35,12 @@ nested_set, print_nested_dict, ) -from floris.wind_data import WindDataBase +from floris.wind_data import ( + TimeSeries, + WindDataBase, + WindRose, + WindTIRose, +) class FlorisModel(LoggingManager): @@ -104,6 +109,8 @@ def __init__(self, configuration: dict | str | Path): ) raise ValueError("turbine_grid_points must be less than or equal to 3.") + # Initialize stored wind_data object to None + self._wind_data = None ### Methods for setting and running the FlorisModel @@ -159,30 +166,31 @@ def _reinitialize( flow_field_dict = floris_dict["flow_field"] farm_dict = floris_dict["farm"] - # Make the given changes - - # First check if wind data is not None, - # if not, get wind speeds, wind direction and - # turbulence intensity using the unpack_for_reinitialize - # method - if wind_data is not None: - if ( - (wind_directions is not None) - or (wind_speeds is not None) - or (turbulence_intensities is not None) - or (heterogenous_inflow_config is not None) - ): + # + if ( + (wind_directions is not None) + or (wind_speeds is not None) + or (turbulence_intensities is not None) + or (heterogenous_inflow_config is not None) + ): + if wind_data is not None: raise ValueError( "If wind_data is passed to reinitialize, then do not pass wind_directions, " "wind_speeds, turbulence_intensities or " "heterogenous_inflow_config as this is redundant" ) - ( - wind_directions, - wind_speeds, - turbulence_intensities, - heterogenous_inflow_config, - ) = wind_data.unpack_for_reinitialize() + elif self.wind_data is not None: + self.logger.warning("Deleting stored wind_data information.") + self._wind_data = None + if wind_data is not None: + # Unpack wind data for reinitialization and save wind_data for use in output + ( + wind_directions, + wind_speeds, + turbulence_intensities, + heterogenous_inflow_config, + ) = wind_data.unpack_for_reinitialize() + self._wind_data = wind_data ## FlowField if wind_speeds is not None: @@ -403,7 +411,7 @@ def run_no_wake(self) -> None: ### Methods for extracting turbine performance after running - def get_turbine_powers(self) -> NDArrayFloat: + def _get_turbine_powers(self) -> NDArrayFloat: """Calculates the power at each turbine in the wind farm. Returns: @@ -413,8 +421,7 @@ def get_turbine_powers(self) -> NDArrayFloat: # Confirm calculate wake has been run if self.core.state is not State.USED: raise RuntimeError( - "Can't run function `FlorisModel.get_turbine_powers` without " - "first running `FlorisModel.run`." + "Can't compute turbine powers without first running `FlorisModel.run()`." ) # Check for negative velocities, which could indicate bad model # parameters or turbines very closely spaced. @@ -436,7 +443,44 @@ def get_turbine_powers(self) -> NDArrayFloat: ) return turbine_powers - def get_farm_power( + + def get_turbine_powers(self): + """ + Calculates the power at each turbine in the wind farm. + + Returns: + NDArrayFloat: Powers at each turbine. + """ + turbine_powers = self._get_turbine_powers() + + if self.wind_data is not None: + if type(self.wind_data) is WindRose: + turbine_powers_rose = np.full( + (len(self.wind_data.wd_flat), self.core.farm.n_turbines), + np.nan + ) + turbine_powers_rose[self.wind_data.non_zero_freq_mask, :] = turbine_powers + turbine_powers = turbine_powers_rose.reshape( + len(self.wind_data.wind_directions), + len(self.wind_data.wind_speeds), + self.core.farm.n_turbines + ) + elif type(self.wind_data) is WindTIRose: + turbine_powers_rose = np.full( + (len(self.wind_data.wd_flat), self.core.farm.n_turbines), + np.nan + ) + turbine_powers_rose[self.wind_data.non_zero_freq_mask, :] = turbine_powers + turbine_powers = turbine_powers_rose.reshape( + len(self.wind_data.wind_directions), + len(self.wind_data.wind_speeds), + len(self.wind_data.turbulence_intensities), + self.core.farm.n_turbines + ) + + return turbine_powers + + def _get_farm_power( self, turbine_weights=None, use_turbulence_correction=False, @@ -462,9 +506,9 @@ def get_farm_power( objective function. If None, this is an array with all values 1.0 and with shape equal to (n_findex, n_turbines). Defaults to None. - use_turbulence_correction: (bool, optional): When *True* uses a + use_turbulence_correction: (bool, optional): When True uses a turbulence parameter to adjust power output calculations. - Defaults to *False*. + Defaults to False. Not currently implemented. Returns: float: Sum of wind turbine powers in W. @@ -475,12 +519,16 @@ def get_farm_power( # TODO: Uncomment out the following two lines once the above are resolved # for turbine in self.core.farm.turbines: # turbine.use_turbulence_correction = use_turbulence_correction + if use_turbulence_correction: + raise NotImplementedError( + "Turbulence correction is not yet implemented in the power calculation." + ) - # Confirm calculate wake has been run + # Confirm run() has been run if self.core.state is not State.USED: raise RuntimeError( - "Can't run function `FlorisModel.get_turbine_powers` without " - "first running `FlorisModel.calculate_wake`." + "Can't run function `FlorisModel.get_farm_power` without " + "first running `FlorisModel.run`." ) if turbine_weights is None: @@ -499,38 +547,82 @@ def get_farm_power( ) # Calculate all turbine powers and apply weights - turbine_powers = self.get_turbine_powers() + turbine_powers = self._get_turbine_powers() turbine_powers = np.multiply(turbine_weights, turbine_powers) return np.sum(turbine_powers, axis=1) - def get_farm_AEP( + def get_farm_power( self, - freq, - cut_in_wind_speed=0.001, - cut_out_wind_speed=None, turbine_weights=None, - no_wake=False, + use_turbulence_correction=False, + ): + """ + Report wind plant power from instance of floris. Optionally includes + uncertainty in wind direction and yaw position when determining power. + Uncertainty is included by computing the mean wind farm power for a + distribution of wind direction and yaw position deviations from the + original wind direction and yaw angles. + + Args: + turbine_weights (NDArrayFloat | list[float] | None, optional): + weighing terms that allow the user to emphasize power at + particular turbines and/or completely ignore the power + from other turbines. This is useful when, for example, you are + modeling multiple wind farms in a single floris object. If you + only want to calculate the power production for one of those + farms and include the wake effects of the neighboring farms, + you can set the turbine_weights for the neighboring farms' + turbines to 0.0. The array of turbine powers from floris + is multiplied with this array in the calculation of the + objective function. If None, this is an array with all values + 1.0 and with shape equal to (n_findex, n_turbines). + Defaults to None. + use_turbulence_correction: (bool, optional): When True uses a + turbulence parameter to adjust power output calculations. + Defaults to False. Not currently implemented. + + Returns: + float: Sum of wind turbine powers in W. + """ + farm_power = self._get_farm_power(turbine_weights, use_turbulence_correction) + + if self.wind_data is not None: + if type(self.wind_data) is WindRose: + farm_power_rose = np.full(len(self.wind_data.wd_flat), np.nan) + farm_power_rose[self.wind_data.non_zero_freq_mask] = farm_power + farm_power = farm_power_rose.reshape( + len(self.wind_data.wind_directions), + len(self.wind_data.wind_speeds) + ) + elif type(self.wind_data) is WindTIRose: + farm_power_rose = np.full(len(self.wind_data.wd_flat), np.nan) + farm_power_rose[self.wind_data.non_zero_freq_mask] = farm_power + farm_power = farm_power_rose.reshape( + len(self.wind_data.wind_directions), + len(self.wind_data.wind_speeds), + len(self.wind_data.turbulence_intensities) + ) + + return farm_power + + def get_expected_farm_power( + self, + freq=None, + turbine_weights=None, ) -> float: """ - Estimate annual energy production (AEP) for distributions of wind speed, wind - direction, frequency of occurrence, and yaw offset. + Compute the expected (mean) power of the wind farm. Args: freq (NDArrayFloat): NumPy array with shape (n_findex) with the frequencies of each wind direction and wind speed combination. These frequencies should typically sum up to 1.0 and are used to weigh the wind farm power for every - condition in calculating the wind farm's AEP. - 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. + condition in calculating the wind farm's AEP. Defaults to None. + If None and a WindData object was supplied, the WindData object's + frequencies will be used. Otherwise, uniform frequencies are assumed + (i.e., a simple mean over the findices is computed). turbine_weights (NDArrayFloat | list[float] | None, optional): weighing terms that allow the user to emphasize power at particular turbines and/or completely ignore the power @@ -544,98 +636,36 @@ def get_farm_AEP( objective function. If None, this is an array with all values 1.0 and with shape equal to (n_findex, n_turbines). 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. """ - # Verify dimensions of the variable "freq" - if np.shape(freq)[0] != self.core.flow_field.n_findex: - raise UserWarning( - "'freq' should be a one-dimensional array with dimensions (n_findex). " - f"Given shape is {np.shape(freq)}" - ) + farm_power = self._get_farm_power(turbine_weights=turbine_weights) - # 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.core.flow_field.wind_speeds, copy=True) - wind_directions = np.array(self.core.flow_field.wind_directions, copy=True) - turbulence_intensities = np.array(self.core.flow_field.turbulence_intensities, copy=True) - farm_power = np.zeros(self.core.flow_field.n_findex) - - # Determine which wind speeds we must evaluate - 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] - wind_directions_subset = wind_directions[conditions_to_evaluate] - turbulence_intensities_subset = turbulence_intensities[conditions_to_evaluate] - self.set( - wind_speeds=wind_speeds_subset, - wind_directions=wind_directions_subset, - turbulence_intensities=turbulence_intensities_subset, - ) - if no_wake: - self.run_no_wake() + if freq is None: + if self.wind_data is None: + freq = np.array([1.0/self.core.flow_field.n_findex]) else: - self.run() - farm_power[conditions_to_evaluate] = self.get_farm_power( - turbine_weights=turbine_weights - ) + freq = self.wind_data.unpack_freq() - # Finally, calculate AEP in GWh - aep = np.sum(np.multiply(freq, farm_power) * 365 * 24) + return np.nansum(np.multiply(freq, farm_power)) - # Reset the FLORIS object to the full wind speed array - self.set( - wind_speeds=wind_speeds, - wind_directions=wind_directions, - turbulence_intensities=turbulence_intensities - ) - - return aep - - def get_farm_AEP_with_wind_data( + def get_farm_AEP( self, - wind_data, - cut_in_wind_speed=0.001, - cut_out_wind_speed=None, + freq=None, turbine_weights=None, - no_wake=False, + hours_per_year=8760, ) -> float: """ Estimate annual energy production (AEP) for distributions of wind speed, wind direction, frequency of occurrence, and yaw offset. Args: - wind_data: (type(WindDataBase)): TimeSeries or WindRose object containing - the wind conditions over which to calculate the AEP. Should match the wind_data - object passed to reinitialize(). - 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. + freq (NDArrayFloat): NumPy array with shape (n_findex) + with the frequencies of each wind direction and + wind speed combination. These frequencies should typically sum + up to 1.0 and are used to weigh the wind farm power for every + condition in calculating the wind farm's AEP. Defaults to None. + If None and a WindData object was supplied, the WindData object's + frequencies will be used. Otherwise, uniform frequencies are assumed. turbine_weights (NDArrayFloat | list[float] | None, optional): weighing terms that allow the user to emphasize power at particular turbines and/or completely ignore the power @@ -649,31 +679,27 @@ def get_farm_AEP_with_wind_data( objective function. If None, this is an array with all values 1.0 and with shape equal to (n_findex, n_turbines). 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*. + hours_per_year (float, optional): Number of hours in a year. Defaults to 365 * 24. Returns: float: The Annual Energy Production (AEP) for the wind farm in watt-hours. """ + if ( + freq is None + and not isinstance(self.wind_data, WindRose) + and not isinstance(self.wind_data, WindTIRose) + ): + self.logger.warning( + "Computing AEP with uniform frequencies. Results results may not reflect annual " + "operation." + ) - # Verify the wind_data object matches FLORIS' initialization - if wind_data.n_findex != self.core.flow_field.n_findex: - raise ValueError("WindData object and floris do not have same findex") - - # Get freq directly from wind_data - freq = wind_data.unpack_freq() - - return self.get_farm_AEP( - freq, - cut_in_wind_speed=cut_in_wind_speed, - cut_out_wind_speed=cut_out_wind_speed, - turbine_weights=turbine_weights, - no_wake=no_wake, - ) + return self.get_expected_farm_power( + freq=freq, + turbine_weights=turbine_weights + ) * hours_per_year def get_turbine_ais(self) -> NDArrayFloat: turbine_ais = axial_induction( @@ -1376,6 +1402,10 @@ def turbine_average_velocities(self) -> NDArrayFloat: cubature_weights=self.core.grid.cubature_weights, ) + @property + def wind_data(self): + return self._wind_data + ### v3 functions that are removed - raise an error if used diff --git a/floris/optimization/layout_optimization/layout_optimization_base.py b/floris/optimization/layout_optimization/layout_optimization_base.py index c8e192d1a..d52e6b1f2 100644 --- a/floris/optimization/layout_optimization/layout_optimization_base.py +++ b/floris/optimization/layout_optimization/layout_optimization_base.py @@ -21,16 +21,15 @@ class LayoutOptimization(LoggingManager): fmodel (FlorisModel): A FlorisModel object. boundaries (iterable(float, float)): Pairs of x- and y-coordinates that represent the boundary's vertices (m). - wind_data (TimeSeries | WindRose): A TimeSeries or WindRose object - values. min_dist (float, optional): The minimum distance to be maintained between turbines during the optimization (m). If not specified, initializes to 2 rotor diameters. Defaults to None. enable_geometric_yaw (bool, optional): If True, enables geometric yaw optimization. Defaults to False. """ - def __init__(self, fmodel, boundaries, wind_data, min_dist=None, enable_geometric_yaw=False): - self.fmodel = fmodel.copy() + def __init__(self, fmodel, boundaries, min_dist=None, enable_geometric_yaw=False): + self.fmodel = fmodel.copy() # Does not copy over the wind_data object + self.fmodel.set(wind_data=fmodel.wind_data) self.boundaries = boundaries self.enable_geometric_yaw = enable_geometric_yaw @@ -49,12 +48,15 @@ def __init__(self, fmodel, boundaries, wind_data, min_dist=None, enable_geometri self.min_dist = min_dist # Check that wind_data is a WindDataBase object - if (not isinstance(wind_data, WindDataBase)): - raise ValueError( - "wind_data entry is not an object of WindDataBase" - " (eg TimeSeries, WindRose, WindTIRose)" + if (not isinstance(self.fmodel.wind_data, WindDataBase)): + # NOTE: it is no longer strictly necessary that fmodel use + # a WindData object, but it is still recommended. + self.logger.warning( + "Running layout optimization without a WindData object (e.g. TimeSeries, WindRose, " + "WindTIRose). We suggest that the user set the wind conditions on the FlorisModel " + " using the wind_data keyword argument for layout optimizations to capture " + "frequencies accurately." ) - self.wind_data = wind_data # Establish geometric yaw class if self.enable_geometric_yaw: @@ -63,8 +65,9 @@ def __init__(self, fmodel, boundaries, wind_data, min_dist=None, enable_geometri minimum_yaw_angle=-30.0, maximum_yaw_angle=30.0, ) - - self.initial_AEP = fmodel.get_farm_AEP_with_wind_data(self.wind_data) + # TODO: is this being used? + fmodel.run() + self.initial_AEP = fmodel.get_farm_AEP() def __str__(self): return "layout" diff --git a/floris/optimization/layout_optimization/layout_optimization_pyoptsparse.py b/floris/optimization/layout_optimization/layout_optimization_pyoptsparse.py index 9d26bc616..959b152a3 100644 --- a/floris/optimization/layout_optimization/layout_optimization_pyoptsparse.py +++ b/floris/optimization/layout_optimization/layout_optimization_pyoptsparse.py @@ -12,7 +12,6 @@ def __init__( self, fmodel, boundaries, - wind_data, min_dist=None, solver=None, optOptions=None, @@ -21,7 +20,7 @@ def __init__( hotStart=None, enable_geometric_yaw=False, ): - super().__init__(fmodel, boundaries, wind_data=wind_data, min_dist=min_dist, + super().__init__(fmodel, boundaries, min_dist=min_dist, enable_geometric_yaw=enable_geometric_yaw) self.x0 = self._norm(self.fmodel.layout_x, self.xmin, self.xmax) @@ -95,15 +94,11 @@ def _obj_func(self, varDict): yaw_angles = self._get_geoyaw_angles() # Update turbine map with turbine locations and yaw angles self.fmodel.set(layout_x=self.x, layout_y=self.y, yaw_angles=yaw_angles) + self.fmodel.run() # Compute the objective function funcs = {} - funcs["obj"] = ( - - -1 * self.fmodel.get_farm_AEP_with_wind_data(self.wind_data) - / self.initial_AEP - - ) + funcs["obj"] = -1 * self.fmodel.get_farm_AEP() / self.initial_AEP # Compute constraints, if any are defined for the optimization funcs = self.compute_cons(funcs, self.x, self.y) diff --git a/floris/optimization/layout_optimization/layout_optimization_pyoptsparse_spread.py b/floris/optimization/layout_optimization/layout_optimization_pyoptsparse_spread.py index aa8d9f54e..ac568d4de 100644 --- a/floris/optimization/layout_optimization/layout_optimization_pyoptsparse_spread.py +++ b/floris/optimization/layout_optimization/layout_optimization_pyoptsparse_spread.py @@ -12,7 +12,6 @@ def __init__( self, fmodel, boundaries, - wind_data, min_dist=None, solver=None, optOptions=None, @@ -20,7 +19,7 @@ def __init__( storeHistory='hist.hist', hotStart=None ): - super().__init__(fmodel, boundaries, wind_data=wind_data, min_dist=min_dist) + super().__init__(fmodel, boundaries, min_dist=min_dist) self._reinitialize(solver=solver, optOptions=optOptions) self.storeHistory = storeHistory diff --git a/floris/optimization/layout_optimization/layout_optimization_scipy.py b/floris/optimization/layout_optimization/layout_optimization_scipy.py index 23c866071..ff3048cae 100644 --- a/floris/optimization/layout_optimization/layout_optimization_scipy.py +++ b/floris/optimization/layout_optimization/layout_optimization_scipy.py @@ -13,7 +13,6 @@ def __init__( self, fmodel, boundaries, - wind_data, bnds=None, min_dist=None, solver='SLSQP', @@ -27,8 +26,6 @@ def __init__( fmodel (FlorisModel): A FlorisModel object. boundaries (iterable(float, float)): Pairs of x- and y-coordinates that represent the boundary's vertices (m). - wind_data (TimeSeries | WindRose): A TimeSeries or WindRose object - values. If None, equal weight is given to each pair of wind conditions bnds (iterable, optional): Bounds for the optimization variables (pairs of min/max values for each variable (m)). If none are specified, they are set to 0 and 1. Defaults to None. @@ -39,8 +36,12 @@ def __init__( optOptions (dict, optional): Dicitonary for setting the optimization options. Defaults to None. """ - super().__init__(fmodel, boundaries, min_dist=min_dist, wind_data=wind_data, - enable_geometric_yaw=enable_geometric_yaw) + super().__init__( + fmodel, + boundaries, + min_dist=min_dist, + enable_geometric_yaw=enable_geometric_yaw + ) self.boundaries_norm = [ [ @@ -98,9 +99,9 @@ def _obj_func(self, locs): # Compute turbine yaw angles using PJ's geometric code (if enabled) yaw_angles = self._get_geoyaw_angles() self.fmodel.set(yaw_angles=yaw_angles) + self.fmodel.run() - return (-1 * self.fmodel.get_farm_AEP_with_wind_data(self.wind_data) / - self.initial_AEP) + return -1 * self.fmodel.get_farm_AEP() / self.initial_AEP def _change_coordinates(self, locs): diff --git a/floris/uncertain_floris_model.py b/floris/uncertain_floris_model.py index 2cfc85b0b..2242f4075 100644 --- a/floris/uncertain_floris_model.py +++ b/floris/uncertain_floris_model.py @@ -12,7 +12,12 @@ NDArrayFloat, ) from floris.utilities import wrap_180 -from floris.wind_data import WindDataBase +from floris.wind_data import ( + TimeSeries, + WindDataBase, + WindRose, + WindTIRose, +) class UncertainFlorisModel(LoggingManager): @@ -99,21 +104,6 @@ def __init__( # Instantiate the expanded FlorisModel # self.core_interface = FlorisModel(configuration) - def copy(self): - """Create an independent copy of the current UncertainFlorisModel object""" - return UncertainFlorisModel( - self.fmodel_unexpanded.core.as_dict(), - wd_resolution=self.wd_resolution, - ws_resolution=self.ws_resolution, - ti_resolution=self.ti_resolution, - yaw_resolution=self.yaw_resolution, - power_setpoint_resolution=self.power_setpoint_resolution, - wd_std=self.wd_std, - wd_sample_points=self.wd_sample_points, - fix_yaw_to_nominal_direction=self.fix_yaw_to_nominal_direction, - verbose=self.verbose, - ) - def set( self, **kwargs, @@ -207,20 +197,6 @@ def _set_uncertain( ], ) - def run(self): - """ - Run the simulation in the underlying FlorisModel object. - """ - - self.fmodel_expanded.run() - - def run_no_wake(self): - """ - Run the simulation in the underlying FlorisModel object without wakes. - """ - - self.fmodel_expanded.run_no_wake() - def reset_operation(self): """ Reset the operation of the underlying FlorisModel object. @@ -235,7 +211,21 @@ def reset_operation(self): # Calling set_uncertain again to reset the expanded FlorisModel self._set_uncertain() - def get_turbine_powers(self): + def run(self): + """ + Run the simulation in the underlying FlorisModel object. + """ + + self.fmodel_expanded.run() + + def run_no_wake(self): + """ + Run the simulation in the underlying FlorisModel object without wakes. + """ + + self.fmodel_expanded.run_no_wake() + + def _get_turbine_powers(self): """Calculates the power at each turbine in the wind farm. This method calculates the power at each turbine in the wind farm, considering @@ -248,7 +238,7 @@ def get_turbine_powers(self): # Pass to off-class function result = map_turbine_powers_uncertain( - unique_turbine_powers=self.fmodel_expanded.get_turbine_powers(), + unique_turbine_powers=self.fmodel_expanded._get_turbine_powers(), map_to_expanded_inputs=self.map_to_expanded_inputs, weights=self.weights, n_unexpanded=self.n_unexpanded, @@ -258,9 +248,58 @@ def get_turbine_powers(self): return result - def get_farm_power( + def get_turbine_powers(self): + """ + Calculate the power at each turbine in the wind farm. If WindRose or + WindTIRose is passed in, result is reshaped to match + + Returns: + NDArrayFloat: An array containing the powers at each turbine for each findex. + """ + + turbine_powers = self._get_turbine_powers() + + if self.fmodel_unexpanded.wind_data is not None: + if type(self.fmodel_unexpanded.wind_data) is WindRose: + turbine_powers_rose = np.full( + ( + len(self.fmodel_unexpanded.wind_data.wd_flat), + self.fmodel_unexpanded.core.farm.n_turbines, + ), + np.nan, + ) + turbine_powers_rose[ + self.fmodel_unexpanded.wind_data.non_zero_freq_mask, : + ] = turbine_powers + turbine_powers = turbine_powers_rose.reshape( + len(self.fmodel_unexpanded.wind_data.wind_directions), + len(self.fmodel_unexpanded.wind_data.wind_speeds), + self.fmodel_unexpanded.core.farm.n_turbines, + ) + elif type(self.fmodel_unexpanded.wind_data) is WindTIRose: + turbine_powers_rose = np.full( + ( + len(self.fmodel_unexpanded.wind_data.wd_flat), + self.fmodel_unexpanded.core.farm.n_turbines, + ), + np.nan, + ) + turbine_powers_rose[ + self.fmodel_unexpanded.wind_data.non_zero_freq_mask, : + ] = turbine_powers + turbine_powers = turbine_powers_rose.reshape( + len(self.fmodel_unexpanded.wind_data.wind_directions), + len(self.fmodel_unexpanded.wind_data.wind_speeds), + len(self.fmodel_unexpanded.wind_data.turbulence_intensities), + self.fmodel_unexpanded.core.farm.n_turbines, + ) + + return turbine_powers + + def _get_farm_power( self, turbine_weights=None, + use_turbulence_correction=False, ): """ Report wind plant power from instance of floris with uncertainty. @@ -279,10 +318,23 @@ def get_farm_power( objective function. If None, this is an array with all values 1.0 and with shape equal to (n_findex, n_turbines). Defaults to None. + use_turbulence_correction: (bool, optional): When True uses a + turbulence parameter to adjust power output calculations. + Defaults to False. Not currently implemented. Returns: float: Sum of wind turbine powers in W. """ + # TODO: Turbulence correction used in the power calculation, but may not be in + # the model yet + # TODO: Turbines need a switch for using turbulence correction + # TODO: Uncomment out the following two lines once the above are resolved + # for turbine in self.core.farm.turbines: + # turbine.use_turbulence_correction = use_turbulence_correction + if use_turbulence_correction: + raise NotImplementedError( + "Turbulence correction is not yet implemented in the power calculation." + ) if turbine_weights is None: # Default to equal weighing of all turbines when turbine_weights is None @@ -300,38 +352,86 @@ def get_farm_power( ) # Calculate all turbine powers and apply weights - turbine_powers = self.get_turbine_powers() + turbine_powers = self._get_turbine_powers() turbine_powers = np.multiply(turbine_weights, turbine_powers) return np.sum(turbine_powers, axis=1) - def get_farm_AEP( + def get_farm_power( self, - freq, - cut_in_wind_speed=0.001, - cut_out_wind_speed=None, turbine_weights=None, - no_wake=False, + use_turbulence_correction=False, + ): + """ + Report wind plant power from instance of floris. Optionally includes + uncertainty in wind direction and yaw position when determining power. + Uncertainty is included by computing the mean wind farm power for a + distribution of wind direction and yaw position deviations from the + original wind direction and yaw angles. + + Args: + turbine_weights (NDArrayFloat | list[float] | None, optional): + weighing terms that allow the user to emphasize power at + particular turbines and/or completely ignore the power + from other turbines. This is useful when, for example, you are + modeling multiple wind farms in a single floris object. If you + only want to calculate the power production for one of those + farms and include the wake effects of the neighboring farms, + you can set the turbine_weights for the neighboring farms' + turbines to 0.0. The array of turbine powers from floris + is multiplied with this array in the calculation of the + objective function. If None, this is an array with all values + 1.0 and with shape equal to (n_findex, n_turbines). + Defaults to None. + use_turbulence_correction: (bool, optional): When True uses a + turbulence parameter to adjust power output calculations. + Defaults to False. Not currently implemented. + + Returns: + float: Sum of wind turbine powers in W. + """ + farm_power = self._get_farm_power(turbine_weights, use_turbulence_correction) + + if self.fmodel_unexpanded.wind_data is not None: + if type(self.fmodel_unexpanded.wind_data) is WindRose: + farm_power_rose = np.full(len(self.fmodel_unexpanded.wind_data.wd_flat), np.nan) + farm_power_rose[ + self.fmodel_unexpanded.wind_data.non_zero_freq_mask + ] = farm_power + farm_power = farm_power_rose.reshape( + len(self.fmodel_unexpanded.wind_data.wind_directions), + len(self.fmodel_unexpanded.wind_data.wind_speeds), + ) + elif type(self.fmodel_unexpanded.wind_data) is WindTIRose: + farm_power_rose = np.full(len(self.fmodel_unexpanded.wind_data.wd_flat), np.nan) + farm_power_rose[ + self.fmodel_unexpanded.wind_data.non_zero_freq_mask + ] = farm_power + farm_power = farm_power_rose.reshape( + len(self.fmodel_unexpanded.wind_data.wind_directions), + len(self.fmodel_unexpanded.wind_data.wind_speeds), + len(self.fmodel_unexpanded.wind_data.turbulence_intensities), + ) + + return farm_power + + def get_expected_farm_power( + self, + freq=None, + turbine_weights=None, ) -> float: """ - Estimate annual energy production (AEP) for distributions of wind speed, wind - direction, frequency of occurrence, and yaw offset. + Compute the expected (mean) power of the wind farm. Args: freq (NDArrayFloat): NumPy array with shape (n_findex) with the frequencies of each wind direction and wind speed combination. These frequencies should typically sum up to 1.0 and are used to weigh the wind farm power for every - condition in calculating the wind farm's AEP. - 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. + condition in calculating the wind farm's AEP. Defaults to None. + If None and a WindData object was supplied, the WindData object's + frequencies will be used. Otherwise, uniform frequencies are assumed + (i.e., a simple mean over the findices is computed). turbine_weights (NDArrayFloat | list[float] | None, optional): weighing terms that allow the user to emphasize power at particular turbines and/or completely ignore the power @@ -345,92 +445,36 @@ def get_farm_AEP( objective function. If None, this is an array with all values 1.0 and with shape equal to (n_findex, n_turbines). 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. """ - # Verify dimensions of the variable "freq" - if np.shape(freq)[0] != self.n_unexpanded: - raise UserWarning( - "'freq' should be a one-dimensional array with dimensions (self.n_unexpanded). " - f"Given shape is {np.shape(freq)}" - ) + farm_power = self._get_farm_power(turbine_weights=turbine_weights) - # 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_directions = np.array(self.wind_directions_unexpanded, copy=True) - wind_speeds = np.array(self.wind_speeds_unexpanded, copy=True) - farm_power = np.zeros_like(wind_directions) - - # Determine which wind speeds we must evaluate - 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] - wind_directions_subset = wind_directions[conditions_to_evaluate] - self.set( - wind_speeds=wind_speeds_subset, - wind_directions=wind_directions_subset, - ) - - if no_wake: - self.run_no_wake() + if freq is None: + if self.fmodel_unexpanded.wind_data is None: + freq = np.array([1.0 / self.core.flow_field.n_findex]) else: - self.run() - farm_power[conditions_to_evaluate] = self.get_farm_power( - turbine_weights=turbine_weights - ) + freq = self.fmodel_unexpanded.wind_data.unpack_freq() - # Finally, calculate AEP in GWh - aep = np.sum(np.multiply(freq, farm_power) * 365 * 24) + return np.nansum(np.multiply(freq, farm_power)) - # Reset the FLORIS object to the full wind speed array - self.set(wind_speeds=wind_speeds, wind_directions=wind_directions) - - return aep - - def get_farm_AEP_with_wind_data( + def get_farm_AEP( self, - wind_data, - cut_in_wind_speed=0.001, - cut_out_wind_speed=None, + freq=None, turbine_weights=None, - no_wake=False, + hours_per_year=8760, ) -> float: """ Estimate annual energy production (AEP) for distributions of wind speed, wind direction, frequency of occurrence, and yaw offset. Args: - wind_data: (type(WindDataBase)): TimeSeries or WindRose object containing - the wind conditions over which to calculate the AEP. Should match the wind_data - object passed to reinitialize(). - 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. + freq (NDArrayFloat): NumPy array with shape (n_findex) + with the frequencies of each wind direction and + wind speed combination. These frequencies should typically sum + up to 1.0 and are used to weigh the wind farm power for every + condition in calculating the wind farm's AEP. Defaults to None. + If None and a WindData object was supplied, the WindData object's + frequencies will be used. Otherwise, uniform frequencies are assumed. turbine_weights (NDArrayFloat | list[float] | None, optional): weighing terms that allow the user to emphasize power at particular turbines and/or completely ignore the power @@ -444,51 +488,28 @@ def get_farm_AEP_with_wind_data( objective function. If None, this is an array with all values 1.0 and with shape equal to (n_findex, n_turbines). 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*. + hours_per_year (float, optional): Number of hours in a year. Defaults to 365 * 24. Returns: float: The Annual Energy Production (AEP) for the wind farm in watt-hours. """ + if ( + freq is None + and not isinstance(self.fmodel_unexpanded.wind_data, WindRose) + and not isinstance(self.fmodel_unexpanded.wind_data, WindTIRose) + ): + self.logger.warning( + "Computing AEP with uniform frequencies. Results results may not reflect annual " + "operation." + ) - # Verify the wind_data object matches FLORIS' initialization - if wind_data.n_findex != self.n_unexpanded: - raise ValueError("WindData object findex not length n_unexpanded") - - # Get freq directly from wind_data - freq = wind_data.unpack_freq() - - return self.get_farm_AEP( - freq, - cut_in_wind_speed=cut_in_wind_speed, - cut_out_wind_speed=cut_out_wind_speed, - turbine_weights=turbine_weights, - no_wake=no_wake, + return ( + self.get_expected_farm_power(freq=freq, turbine_weights=turbine_weights) + * hours_per_year ) - # def copy(self): - # """Create an independent copy of the current UncertainFlorisModel object""" - # return UncertainFlorisModel( - # self.fmodel_unexpanded.core.as_dict(), - # wd_resolution=self.wd_resolution, - # ws_resolution=self.ws_resolution, - # ti_resolution=self.ti_resolution, - # yaw_resolution=self.yaw_resolution, - # power_setpoint_resolution=self.power_setpoint_resolution, - # wd_std=self.wd_std, - # wd_sample_points=self.wd_sample_points, - # verbose=self.verbose, - # ) - - # @property - # def core(self): - # """Return core of underlying expanded FlorisModel object""" - # return self.fmodel_expanded.core - def _get_rounded_inputs( self, input_array, @@ -617,7 +638,6 @@ def _expand_wind_directions( # If fix_yaw_to_nominal_direction is True, set the yaw angle to relative # to the nominal wind direction if fix_yaw_to_nominal_direction: - # Wrap between -180 and 180 output_array[start_idx:end_idx, 3 : 3 + n_turbines] = wrap_180( output_array[start_idx:end_idx, 3 : 3 + n_turbines] + wd_sample_points[i] @@ -670,6 +690,21 @@ def _get_weights(self, wd_std, wd_sample_points): return weights + def copy(self): + """Create an independent copy of the current UncertainFlorisModel object""" + return UncertainFlorisModel( + self.fmodel_unexpanded.core.as_dict(), + wd_resolution=self.wd_resolution, + ws_resolution=self.ws_resolution, + ti_resolution=self.ti_resolution, + yaw_resolution=self.yaw_resolution, + power_setpoint_resolution=self.power_setpoint_resolution, + wd_std=self.wd_std, + wd_sample_points=self.wd_sample_points, + fix_yaw_to_nominal_direction=self.fix_yaw_to_nominal_direction, + verbose=self.verbose, + ) + @property def layout_x(self): """ diff --git a/floris/wind_data.py b/floris/wind_data.py index 2ecac6fac..2b8952e9f 100644 --- a/floris/wind_data.py +++ b/floris/wind_data.py @@ -49,16 +49,7 @@ def unpack_for_reinitialize(self): def unpack_freq(self): """Unpack frequency weighting""" - ( - _, - _, - _, - freq_table_unpack, - _, - _, - ) = self.unpack() - - return freq_table_unpack + return self.unpack()[3] def check_heterogenous_inflow_config_by_wd(self, heterogenous_inflow_config_by_wd): """ @@ -398,7 +389,7 @@ def resample_wind_rose(self, wd_step=None, ws_step=None): ) # Now build a new wind rose using the new steps - return time_series.to_wind_rose( + return time_series.to_WindRose( wd_step=wd_step, ws_step=ws_step, bin_weights=self.freq_table_flat ) @@ -619,7 +610,7 @@ def read_csv_long(file_path: str, time_series = TimeSeries(wind_directions, wind_speeds, turbulence_intensities) # Now build a new wind rose using the new steps - return time_series.to_wind_rose( + return time_series.to_WindRose( wd_step=wd_step, ws_step=ws_step, bin_weights=freq_values ) @@ -851,7 +842,7 @@ def resample_wind_rose(self, wd_step=None, ws_step=None, ti_step=None): ) # Now build a new wind rose using the new steps - return time_series.to_wind_ti_rose( + return time_series.to_WindTIRose( wd_step=wd_step, ws_step=ws_step, ti_step=ti_step, bin_weights=self.freq_table_flat ) @@ -1058,7 +1049,7 @@ def read_csv_long(file_path: str, time_series = TimeSeries(wind_directions, wind_speeds, turbulence_intensities) # Now build a new wind rose using the new steps - return time_series.to_wind_ti_rose( + return time_series.to_WindTIRose( wd_step=wd_step, ws_step=ws_step, ti_step=ti_step,bin_weights=freq_values ) @@ -1285,7 +1276,7 @@ def iref_func(wind_directions, wind_speeds): self.assign_ti_using_wd_ws_function(iref_func) - def to_wind_rose( + def to_WindRose( self, wd_step=2.0, ws_step=1.0, wd_edges=None, ws_edges=None, bin_weights=None ): """ @@ -1425,7 +1416,7 @@ def to_wind_rose( self.heterogenous_inflow_config_by_wd, ) - def to_wind_ti_rose( + def to_WindTIRose( self, wd_step=2.0, ws_step=1.0, diff --git a/tests/floris_model_integration_test.py b/tests/floris_model_integration_test.py index 3bb210cda..ae5f07558 100644 --- a/tests/floris_model_integration_test.py +++ b/tests/floris_model_integration_test.py @@ -1,10 +1,11 @@ +import logging from pathlib import Path import numpy as np import pytest import yaml -from floris import FlorisModel +from floris import FlorisModel, WindRose from floris.core.turbine.operation_models import POWER_SETPOINT_DEFAULT @@ -341,7 +342,7 @@ def test_disable_turbines(): fmodel.run() assert (fmodel.core.farm.yaw_angles == np.array([[1.0, 0.0, 1.0], [1.0, 0.0, 1.0]])).all() -def test_get_farm_aep(): +def test_get_farm_aep(caplog): fmodel = FlorisModel(configuration=YAML_INPUT) wind_speeds = np.array([8.0, 8.0, 8.0]) @@ -369,56 +370,25 @@ def test_get_farm_aep(): freq = np.ones(n_findex) freq = freq / np.sum(freq) - farm_aep = fmodel.get_farm_AEP(freq=freq) - - aep = np.sum(np.multiply(freq, farm_powers) * 365 * 24) + # Check warning raised if freq not passed; no warning if freq passed + with caplog.at_level(logging.WARNING): + fmodel.get_farm_AEP() + assert caplog.text != "" # Checking not empty + caplog.clear() + with caplog.at_level(logging.WARNING): + fmodel.get_farm_AEP(freq=freq) + assert caplog.text == "" # Checking empty - # In this case farm_aep should match farm powers - np.testing.assert_allclose(farm_aep, aep) - -def test_get_farm_aep_with_conditions(): - fmodel = FlorisModel(configuration=YAML_INPUT) - - wind_speeds = np.array([5.0, 8.0, 8.0, 8.0, 20.0]) - wind_directions = np.array([270.0, 270.0, 270.0, 270.0, 270.0]) - turbulence_intensities = np.array([0.06, 0.06, 0.06, 0.06, 0.06]) - n_findex = len(wind_directions) - - layout_x = np.array([0, 0]) - layout_y = np.array([0, 1000]) - # n_turbines = len(layout_x) - - fmodel.set( - wind_speeds=wind_speeds, - wind_directions=wind_directions, - turbulence_intensities=turbulence_intensities, - layout_x=layout_x, - layout_y=layout_y, - ) - - fmodel.run() - - farm_powers = fmodel.get_farm_power() - - # Start with uniform frequency - freq = np.ones(n_findex) - freq = freq / np.sum(freq) - - # Get farm AEP with conditions on minimun and max wind speed - # which exclude the first and last findex - farm_aep = fmodel.get_farm_AEP(freq=freq, cut_in_wind_speed=6.0, cut_out_wind_speed=15.0) + farm_aep = fmodel.get_farm_AEP(freq=freq) - # In this case the aep should be computed assuming 0 power - # for the 0th and last findex - farm_powers[0] = 0 - farm_powers[-1] = 0 aep = np.sum(np.multiply(freq, farm_powers) * 365 * 24) # In this case farm_aep should match farm powers np.testing.assert_allclose(farm_aep, aep) - #Confirm n_findex reset after the operation - assert n_findex == fmodel.core.flow_field.n_findex + # Also check get_expected_farm_power + expected_farm_power = fmodel.get_expected_farm_power(freq=freq) + np.testing.assert_allclose(expected_farm_power, aep / (365 * 24)) def test_set_ti(): fmodel = FlorisModel(configuration=YAML_INPUT) @@ -494,6 +464,102 @@ def test_calculate_planes(): with pytest.raises(ValueError): fmodel.calculate_cross_plane(500.0, ws=[wind_speeds[0]], wd=[wind_directions[0]]) +def test_get_turbine_powers_with_WindRose(): + fmodel = FlorisModel(configuration=YAML_INPUT) + + wind_speeds = np.array([8.0, 10.0, 12.0, 8.0, 10.0, 12.0]) + wind_directions = np.array([270.0, 270.0, 270.0, 280.0, 280.0, 280.0]) + turbulence_intensities = 0.06 * np.ones_like(wind_speeds) + + fmodel.set( + wind_speeds=wind_speeds, + wind_directions=wind_directions, + turbulence_intensities=turbulence_intensities, + layout_x=[0, 1000, 2000, 3000], + layout_y=[0, 0, 0, 0] + ) + fmodel.run() + turbine_powers_simple = fmodel.get_turbine_powers() + + # Now declare a WindRose with 2 wind directions and 3 wind speeds + # uniform TI and frequency + wind_rose = WindRose( + wind_directions=np.unique(wind_directions), + wind_speeds=np.unique(wind_speeds), + ti_table=0.06 + ) + + # Set this wind rose, run + fmodel.set(wind_data=wind_rose) + fmodel.run() + + # Get the turbine powers in the wind rose + turbine_powers_windrose = fmodel.get_turbine_powers() + + # Turbine power should have shape (n_wind_directions, n_wind_speeds, n_turbines) + assert turbine_powers_windrose.shape == (2, 3, 4) + assert np.allclose(turbine_powers_simple.reshape(2, 3, 4), turbine_powers_windrose) + assert np.allclose(turbine_powers_simple, turbine_powers_windrose.reshape(2*3, 4)) + + # Test that if certain combinations in the wind rose have 0 frequency, the power in + # those locations is nan + wind_rose = WindRose( + wind_directions = np.array([270.0, 280.0]), + wind_speeds = np.array([8.0, 10.0, 12.0]), + ti_table=0.06, + freq_table=np.array([[0.25, 0.25, 0.0], [0.0, 0.0, 0.5]]) + ) + fmodel.set(wind_data=wind_rose) + fmodel.run() + turbine_powers = fmodel.get_turbine_powers() + assert np.isnan(turbine_powers[0, 2, 0]) + +def test_get_powers_with_wind_data(): + fmodel = FlorisModel(configuration=YAML_INPUT) + + wind_speeds = np.array([8.0, 10.0, 12.0, 8.0, 10.0, 12.0]) + wind_directions = np.array([270.0, 270.0, 270.0, 280.0, 280.0, 280.0]) + turbulence_intensities = 0.06 * np.ones_like(wind_speeds) + + fmodel.set( + wind_speeds=wind_speeds, + wind_directions=wind_directions, + turbulence_intensities=turbulence_intensities, + layout_x=[0, 1000, 2000, 3000], + layout_y=[0, 0, 0, 0] + ) + fmodel.run() + farm_power_simple = fmodel.get_farm_power() + + # Now declare a WindRose with 2 wind directions and 3 wind speeds + # uniform TI and frequency + wind_rose = WindRose( + wind_directions=np.unique(wind_directions), + wind_speeds=np.unique(wind_speeds), + ti_table=0.06 + ) + + # Set this wind rose, run + fmodel.set(wind_data=wind_rose) + fmodel.run() + + farm_power_windrose = fmodel.get_farm_power() + + # Check dimensions and that the farm power is the sum of the turbine powers + assert farm_power_windrose.shape == (2, 3) + assert np.allclose(farm_power_windrose, fmodel.get_turbine_powers().sum(axis=2)) + + # Check that simple and windrose powers are consistent + assert np.allclose(farm_power_simple.reshape(2, 3), farm_power_windrose) + assert np.allclose(farm_power_simple, farm_power_windrose.flatten()) + + # Test that if the last turbine's weight is set to 0, the farm power is the same as the + # sum of the first 3 turbines + turbine_weights = np.array([1.0, 1.0, 1.0, 0.0]) + farm_power_weighted = fmodel.get_farm_power(turbine_weights=turbine_weights) + + assert np.allclose(farm_power_weighted, fmodel.get_turbine_powers()[:,:,:-1].sum(axis=2)) + def test_get_and_set_param(): fmodel = FlorisModel(configuration=YAML_INPUT) diff --git a/tests/layout_optimization_integration_test.py b/tests/layout_optimization_integration_test.py index dafd5e0d6..0732b969c 100644 --- a/tests/layout_optimization_integration_test.py +++ b/tests/layout_optimization_integration_test.py @@ -1,3 +1,4 @@ +import logging from pathlib import Path import numpy as np @@ -21,7 +22,7 @@ YAML_INPUT = TEST_DATA / "input_full.yaml" -def test_base_class(): +def test_base_class(caplog): # Get a test fi fmodel = FlorisModel(configuration=YAML_INPUT) @@ -32,24 +33,40 @@ def test_base_class(): # (this should fail) freq = np.ones((5, 5)) freq = freq / freq.sum() - with pytest.raises(ValueError): - LayoutOptimization(fmodel, boundaries, freq, 5) - # Passing as a keyword freq to wind_data should also fail - with pytest.raises(ValueError): - LayoutOptimization(fmodel=fmodel, boundaries=boundaries, wind_data=freq, min_dist=5,) + # Check that warning is raised if fmodel does not contain wind_data + with caplog.at_level(logging.WARNING): + LayoutOptimization(fmodel, boundaries, 5) + assert caplog.text != "" # Checking not empty + + caplog.clear() + with caplog.at_level(logging.WARNING): + LayoutOptimization(fmodel=fmodel, boundaries=boundaries, min_dist=5,) + assert caplog.text != "" # Checking not empty time_series = TimeSeries( wind_directions=fmodel.core.flow_field.wind_directions, wind_speeds=fmodel.core.flow_field.wind_speeds, turbulence_intensities=fmodel.core.flow_field.turbulence_intensities, ) - wind_rose = time_series.to_wind_rose() + fmodel.set(wind_data=time_series) + + caplog.clear() + with caplog.at_level(logging.WARNING): + LayoutOptimization(fmodel, boundaries, 5) + assert caplog.text != "" # Not empty, because get_farm_AEP called on TimeSeries + + # Passing without keyword arguments should work, or with keyword arguments + LayoutOptimization(fmodel, boundaries, 5) + LayoutOptimization(fmodel=fmodel, boundaries=boundaries, min_dist=5) + + # Check with WindRose on fmodel + fmodel.set(wind_data=time_series.to_WindRose()) - # Passing wind_data objects in the 3rd position should not fail - LayoutOptimization(fmodel, boundaries, time_series, 5) - LayoutOptimization(fmodel, boundaries, wind_rose, 5) + caplog.clear() + with caplog.at_level(logging.WARNING): + LayoutOptimization(fmodel, boundaries, 5) + assert caplog.text == "" # Empty - # Passing wind_data objects by keyword should not fail - LayoutOptimization(fmodel=fmodel, boundaries=boundaries, wind_data=time_series, min_dist=5) - LayoutOptimization(fmodel=fmodel, boundaries=boundaries, wind_data=wind_rose, min_dist=5) + LayoutOptimization(fmodel, boundaries, 5) + LayoutOptimization(fmodel=fmodel, boundaries=boundaries, min_dist=5) diff --git a/tests/parallel_floris_model_integration_test.py b/tests/parallel_floris_model_integration_test.py index e5d603adf..4b4d5aeec 100644 --- a/tests/parallel_floris_model_integration_test.py +++ b/tests/parallel_floris_model_integration_test.py @@ -60,6 +60,8 @@ def test_parallel_get_AEP(sample_inputs_fixture): fmodel = FlorisModel(sample_inputs_fixture.core) pfmodel_input = copy.deepcopy(fmodel) + + fmodel.run() serial_farm_AEP = fmodel.get_farm_AEP(freq=freq) pfmodel = ParallelFlorisModel( @@ -122,6 +124,7 @@ def test_parallel_uncertain_get_AEP(sample_inputs_fixture): wd_std=3 ) pfmodel_input = copy.deepcopy(ufmodel) + ufmodel.run() serial_farm_AEP = ufmodel.get_farm_AEP(freq=freq) pfmodel = ParallelFlorisModel( diff --git a/tests/uncertain_floris_model_integration_test.py b/tests/uncertain_floris_model_integration_test.py index c6bfb0f8e..42ac9ec8a 100644 --- a/tests/uncertain_floris_model_integration_test.py +++ b/tests/uncertain_floris_model_integration_test.py @@ -6,7 +6,7 @@ from floris import FlorisModel from floris.core.turbine.operation_models import POWER_SETPOINT_DEFAULT -from floris.uncertain_floris_model import UncertainFlorisModel +from floris.uncertain_floris_model import UncertainFlorisModel, WindRose TEST_DATA = Path(__file__).resolve().parent / "data" @@ -215,3 +215,50 @@ def test_uncertain_floris_model_setpoints(): unc_powers = ufmodel.get_turbine_powers()[:, 1].flatten() np.testing.assert_allclose(np.sum(nom_powers * weights), unc_powers) + + +def test_get_powers_with_wind_data(): + ufmodel = UncertainFlorisModel(configuration=YAML_INPUT) + + wind_speeds = np.array([8.0, 10.0, 12.0, 8.0, 10.0, 12.0]) + wind_directions = np.array([270.0, 270.0, 270.0, 280.0, 280.0, 280.0]) + turbulence_intensities = 0.06 * np.ones_like(wind_speeds) + + ufmodel.set( + wind_speeds=wind_speeds, + wind_directions=wind_directions, + turbulence_intensities=turbulence_intensities, + layout_x=[0, 1000, 2000, 3000], + layout_y=[0, 0, 0, 0] + ) + ufmodel.run() + farm_power_simple = ufmodel.get_farm_power() + + # Now declare a WindRose with 2 wind directions and 3 wind speeds + # uniform TI and frequency + wind_rose = WindRose( + wind_directions=np.unique(wind_directions), + wind_speeds=np.unique(wind_speeds), + ti_table=0.06 + ) + + # Set this wind rose, run + ufmodel.set(wind_data=wind_rose) + ufmodel.run() + + farm_power_windrose = ufmodel.get_farm_power() + + # Check dimensions and that the farm power is the sum of the turbine powers + assert farm_power_windrose.shape == (2, 3) + assert np.allclose(farm_power_windrose, ufmodel.get_turbine_powers().sum(axis=2)) + + # Check that simple and windrose powers are consistent + assert np.allclose(farm_power_simple.reshape(2, 3), farm_power_windrose) + assert np.allclose(farm_power_simple, farm_power_windrose.flatten()) + + # Test that if the last turbine's weight is set to 0, the farm power is the same as the + # sum of the first 3 turbines + turbine_weights = np.array([1.0, 1.0, 1.0, 0.0]) + farm_power_weighted = ufmodel.get_farm_power(turbine_weights=turbine_weights) + + assert np.allclose(farm_power_weighted, ufmodel.get_turbine_powers()[:,:,:-1].sum(axis=2)) diff --git a/tests/wind_data_integration_test.py b/tests/wind_data_integration_test.py index 778c35403..c6398a1fa 100644 --- a/tests/wind_data_integration_test.py +++ b/tests/wind_data_integration_test.py @@ -218,12 +218,12 @@ def test_wrap_wind_directions_near_360(): assert np.allclose(wd_wrapped, expected_result) -def test_time_series_to_wind_rose(): +def test_time_series_to_WindRose(): # Test just 1 wind speed wind_directions = np.array([259.8, 260.2, 264.3]) wind_speeds = np.array([5.0, 5.0, 5.1]) time_series = TimeSeries(wind_directions, wind_speeds, 0.06) - wind_rose = time_series.to_wind_rose(wd_step=2.0, ws_step=1.0) + wind_rose = time_series.to_WindRose(wd_step=2.0, ws_step=1.0) # The wind directions should be 260, 262 and 264 because they're binned # to the nearest 2 deg increment @@ -243,7 +243,7 @@ def test_time_series_to_wind_rose(): wind_directions = np.array([259.8, 260.2, 264.3]) wind_speeds = np.array([5.0, 5.0, 6.1]) time_series = TimeSeries(wind_directions, wind_speeds, 0.06) - wind_rose = time_series.to_wind_rose(wd_step=2.0, ws_step=1.0) + wind_rose = time_series.to_WindRose(wd_step=2.0, ws_step=1.0) # The wind directions should be 260, 262 and 264 assert np.allclose(wind_rose.wind_directions, [260, 262, 264]) @@ -267,11 +267,11 @@ def test_time_series_to_wind_rose(): assert np.allclose(ti_table[~np.isnan(ti_table)], 0.06) -def test_time_series_to_wind_rose_wrapping(): +def test_time_series_to_WindRose_wrapping(): wind_directions = np.arange(0.0, 360.0, 0.25) wind_speeds = 8.0 * np.ones_like(wind_directions) time_series = TimeSeries(wind_directions, wind_speeds, 0.06) - wind_rose = time_series.to_wind_rose(wd_step=2.0, ws_step=1.0) + wind_rose = time_series.to_WindRose(wd_step=2.0, ws_step=1.0) # Expert for the first bin in this case to be 0, and the final to be 358 # and both to have equal numbers of points @@ -280,7 +280,7 @@ def test_time_series_to_wind_rose_wrapping(): np.testing.assert_almost_equal(wind_rose.freq_table[0, 0], wind_rose.freq_table[-1, 0]) -def test_time_series_to_wind_rose_with_ti(): +def test_time_series_to_WindRose_with_ti(): wind_directions = np.array([259.8, 260.2, 260.3, 260.1]) wind_speeds = np.array([5.0, 5.0, 5.1, 7.2]) turbulence_intensities = np.array([0.5, 1.0, 1.5, 2.0]) @@ -289,7 +289,7 @@ def test_time_series_to_wind_rose_with_ti(): wind_speeds, turbulence_intensities=turbulence_intensities, ) - wind_rose = time_series.to_wind_rose(wd_step=2.0, ws_step=1.0) + wind_rose = time_series.to_WindRose(wd_step=2.0, ws_step=1.0) # Turbulence intensity should average to 1 in the 5 m/s bin and 2 in the 7 m/s bin ti_table = wind_rose.ti_table @@ -460,7 +460,7 @@ def test_wind_ti_rose_resample(): ) -def test_time_series_to_wind_ti_rose(): +def test_time_series_to_WindTIRose(): wind_directions = np.array([259.8, 260.2, 260.3, 260.1]) wind_speeds = np.array([5.0, 5.0, 5.1, 7.2]) turbulence_intensities = np.array([0.05, 0.1, 0.15, 0.2]) @@ -469,7 +469,7 @@ def test_time_series_to_wind_ti_rose(): wind_speeds, turbulence_intensities=turbulence_intensities, ) - wind_rose = time_series.to_wind_ti_rose(wd_step=2.0, ws_step=1.0, ti_step=0.1) + wind_rose = time_series.to_WindTIRose(wd_step=2.0, ws_step=1.0, ti_step=0.1) # The binning should result in turbulence intensity bins of 0.1 and 0.2 tis_windrose = wind_rose.turbulence_intensities