Example: Optimize yaw and compare AEP

This example demonstrates how to perform a yaw optimization and evaluate the performanceover a full wind rose.The script performs the following steps: 1. Load a wind rose from a csv file 2. Calculates the optimal yaw angles for a wind speed of 8 m/s across the directions 3. Applies the optimal yaw angles to the wind rose and calculates the AEP

from time import perf_counter as timerpc

import matplotlib.pyplot as plt
import numpy as np

from floris import (
    FlorisModel,
    TimeSeries,
    WindRose,
)
from floris.optimization.yaw_optimization.yaw_optimizer_sr import YawOptimizationSR


# Load the wind rose from csv
wind_rose = WindRose.read_csv_long(
    "../inputs/wind_rose.csv", wd_col="wd", ws_col="ws", freq_col="freq_val", ti_col_or_value=0.06
)

# Load FLORIS
fmodel = FlorisModel("../inputs/gch.yaml")

# Specify wind farm layout and update in the floris object
N = 2  # number of turbines per row and per column
X, Y = np.meshgrid(
    5.0 * fmodel.core.farm.rotor_diameters_sorted[0][0] * np.arange(0, N, 1),
    5.0 * fmodel.core.farm.rotor_diameters_sorted[0][0] * np.arange(0, N, 1),
)
fmodel.set(layout_x=X.flatten(), layout_y=Y.flatten())

# Get the number of turbines
n_turbines = len(fmodel.layout_x)

# Optimize the yaw angles.  This could be done for every wind direction and wind speed
# but in practice it is much faster to optimize only for one speed and infer the rest
# using a rule of thumb
time_series = TimeSeries(
    wind_directions=wind_rose.wind_directions, wind_speeds=8.0, turbulence_intensities=0.06
)
fmodel.set(wind_data=time_series)

# Get the optimal angles
start_time = timerpc()
yaw_opt = YawOptimizationSR(
    fmodel=fmodel,
    minimum_yaw_angle=0.0,  # Allowable yaw angles lower bound
    maximum_yaw_angle=20.0,  # Allowable yaw angles upper bound
    Ny_passes=[5, 4],
    exclude_downstream_turbines=True,
)
df_opt = yaw_opt.optimize()
end_time = timerpc()
t_tot = end_time - start_time
print("Optimization finished in {:.2f} seconds.".format(t_tot))


# Calculate the AEP in the baseline case
fmodel.set(wind_data=wind_rose)
fmodel.run()
farm_power_baseline = fmodel.get_farm_power()
aep_baseline = fmodel.get_farm_AEP()


# Now need to apply the optimal yaw angles to the wind rose to get the optimized AEP
# do this by applying a rule of thumb where the optimal yaw is applied between 6 and 12 m/s
# and ramped down to 0 above and below this range

# Grab wind speeds and wind directions from the fmodel.  Note that we do this because the
# yaw angles will need to be n_findex long, and accounting for the fact that some wind
# directions and wind speeds may not be present in the wind rose (0 frequency) and aren't
# included in the fmodel
wind_directions = fmodel.wind_directions
wind_speeds = fmodel.wind_speeds
n_findex = fmodel.n_findex


# Now define how the optimal yaw angles for 8 m/s are applied over the other wind speeds
yaw_angles_opt = np.vstack(df_opt["yaw_angles_opt"])
yaw_angles_wind_rose = np.zeros((n_findex, n_turbines))
for i in range(n_findex):
    wind_speed = wind_speeds[i]
    wind_direction = wind_directions[i]

    # Interpolate the optimal yaw angles for this wind direction from df_opt
    id_opt = df_opt["wind_direction"] == wind_direction
    yaw_opt_full = np.array(df_opt.loc[id_opt, "yaw_angles_opt"])[0]

    # Now decide what to do for different wind speeds
    if (wind_speed < 4.0) | (wind_speed > 14.0):
        yaw_opt = np.zeros(n_turbines)  # do nothing for very low/high speeds
    elif wind_speed < 6.0:
        yaw_opt = yaw_opt_full * (6.0 - wind_speed) / 2.0  # Linear ramp up
    elif wind_speed > 12.0:
        yaw_opt = yaw_opt_full * (14.0 - wind_speed) / 2.0  # Linear ramp down
    else:
        yaw_opt = yaw_opt_full  # Apply full offsets between 6.0 and 12.0 m/s

    # Save to collective array
    yaw_angles_wind_rose[i, :] = yaw_opt


# Now apply the optimal yaw angles and get the AEP
fmodel.set(yaw_angles=yaw_angles_wind_rose)
fmodel.run()
aep_opt = fmodel.get_farm_AEP()
aep_uplift = 100.0 * (aep_opt / aep_baseline - 1)
farm_power_opt = fmodel.get_farm_power()

print("Baseline AEP: {:.2f} GWh.".format(aep_baseline/1E9))
print("Optimal AEP: {:.2f} GWh.".format(aep_opt/1E9))
print("Relative AEP uplift by wake steering: {:.3f} %.".format(aep_uplift))

# Use farm_power_baseline, farm_power_opt and wind_data to make a heat map of uplift by
# wind direction and wind speed
wind_directions = wind_rose.wind_directions
wind_speeds = wind_rose.wind_speeds
relative_gain = farm_power_opt - farm_power_baseline

# Plot the heatmap with wind speeds on x, wind directions on y and relative gain as the color
fig, ax = plt.subplots(figsize=(10, 12))
cax = ax.imshow(relative_gain, cmap='viridis', aspect='auto')
fig.colorbar(cax, ax=ax, label="Relative gain (%)")

ax.set_yticks(np.arange(len(wind_directions)))
ax.set_yticklabels(wind_directions)
ax.set_xticks(np.arange(len(wind_speeds)))
ax.set_xticklabels(wind_speeds)
ax.set_ylabel("Wind direction (deg)")
ax.set_xlabel("Wind speed (m/s)")

# Reduce x and y tick font size
for tick in ax.yaxis.get_major_ticks():
    tick.label1.set_fontsize(8)

for tick in ax.xaxis.get_major_ticks():
    tick.label1.set_fontsize(8)

# Set y ticks to be horizontal
for tick in ax.get_yticklabels():
    tick.set_rotation(0)

ax.set_title("Uplift in farm power by wind direction and wind speed", fontsize=12)

plt.tight_layout()
plt.show()
import warnings
warnings.filterwarnings('ignore')

[Serial Refine] Processing pass=0, turbine_depth=0 (0.0%)
[Serial Refine] Processing pass=0, turbine_depth=1 (12.5%)

[Serial Refine] Processing pass=0, turbine_depth=2 (25.0%)
[Serial Refine] Processing pass=0, turbine_depth=3 (37.5%)

[Serial Refine] Processing pass=1, turbine_depth=0 (50.0%)
[Serial Refine] Processing pass=1, turbine_depth=1 (62.5%)

[Serial Refine] Processing pass=1, turbine_depth=2 (75.0%)
[Serial Refine] Processing pass=1, turbine_depth=3 (87.5%)

Optimization finished in 1.56 seconds.

/home/runner/work/floris/floris/floris/core/wake_deflection/gauss.py:328: RuntimeWarning: invalid value encountered in divide
  val = 2 * (avg_v - v_core) / (v_top + v_bottom)
/home/runner/work/floris/floris/floris/core/wake_deflection/gauss.py:163: RuntimeWarning: invalid value encountered in divide
  C0 = 1 - u0 / freestream_velocity
/home/runner/work/floris/floris/floris/core/wake_deflection/gauss.py:498: RuntimeWarning: invalid value encountered in divide
  I_total = np.sqrt((2 / 3) * k_total) / average_u_i
/home/runner/work/floris/floris/floris/core/wake_velocity/gauss.py:80: RuntimeWarning: invalid value encountered in divide
  sigma_z0 = rotor_diameter_i * 0.5 * np.sqrt(uR / (u_initial + u0))

Baseline AEP: 75.48 GWh.
Optimal AEP: 76.18 GWh.
Relative AEP uplift by wake steering: 0.926 %.