Example: Layout optimization with heterogeneous inflow

"""Example: Layout optimization with heterogeneous inflow
This example shows a layout optimization using the geometric yaw option. It
combines elements of layout optimization and heterogeneous
inflow for demonstrative purposes.

Heterogeneity in the inflow provides the necessary driver for coupled yaw
and layout optimization to be worthwhile. First, a layout optimization is
run without coupled yaw optimization; then a coupled optimization is run to
show the benefits of coupled optimization when flows are heterogeneous.
"""


import os

import matplotlib.pyplot as plt
import numpy as np

from floris import FlorisModel, WindRose
from floris.optimization.layout_optimization.layout_optimization_scipy import (
    LayoutOptimizationScipy,
)


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

# Setup 2 wind directions (due east and due west)
# and 1 wind speed with uniform probability
wind_directions = np.array([90.0, 270.0])
n_wds = len(wind_directions)
wind_speeds = np.array([8.0])

# Shape frequency distribution to match number of wind directions and wind speeds
freq_table = np.ones((len(wind_directions), len(wind_speeds)))
freq_table = freq_table / freq_table.sum()


# The boundaries for the turbines, specified as vertices
D = 126.0  # rotor diameter for the NREL 5MW
size_D = 12
boundaries = [(0.0, 0.0), (size_D * D, 0.0), (size_D * D, 0.1), (0.0, 0.1), (0.0, 0.0)]

# Set turbine locations to 4 turbines at corners of the rectangle
# (optimal without flow heterogeneity)
layout_x = [0.1, 0.3 * size_D * D, 0.6 * size_D * D]
layout_y = [0, 0, 0]

# Generate exaggerated heterogeneous inflow (same for all wind directions)
speed_multipliers = np.repeat(np.array([0.5, 1.0, 0.5, 1.0])[None, :], n_wds, axis=0)
x_locs = [0, size_D * D, 0, size_D * D]
y_locs = [-D, -D, D, D]

# Create the configuration dictionary to be used for the heterogeneous inflow.
heterogeneous_inflow_config_by_wd = {
    "speed_multipliers": speed_multipliers,
    "wind_directions": wind_directions,
    "x": x_locs,
    "y": y_locs,
}

# Establish a WindRose object
wind_rose = WindRose(
    wind_directions=wind_directions,
    wind_speeds=wind_speeds,
    freq_table=freq_table,
    ti_table=0.06,
    heterogeneous_inflow_config_by_wd=heterogeneous_inflow_config_by_wd,
)


fmodel.set(
    layout_x=layout_x,
    layout_y=layout_y,
    wind_data=wind_rose,
)

# Setup and solve the layout optimization problem without heterogeneity
maxiter = 100
layout_opt = LayoutOptimizationScipy(
    fmodel, boundaries, min_dist=2 * D, optOptions={"maxiter": maxiter}
)

# Run the optimization
np.random.seed(0)
sol = layout_opt.optimize()

# Get the resulting improvement in AEP
print("... calcuating improvement in AEP")

fmodel.run()
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() / 1e6

percent_gain = 100 * (opt_aep - base_aep) / base_aep

# Print and plot the results
print(f"Optimal layout: {sol}")
print(
    f"Optimal layout improves AEP by {percent_gain:.1f}% "
    f"from {base_aep:.1f} MWh to {opt_aep:.1f} MWh"
)
layout_opt.plot_layout_opt_results()
ax = plt.gca()
fig = plt.gcf()
sm = ax.tricontourf(x_locs, y_locs, speed_multipliers[0], cmap="coolwarm")
fig.colorbar(sm, ax=ax, label="Speed multiplier")
ax.legend(["_Optimization boundary", "Initial layout", "Optimized layout" ])
ax.set_title("Geometric yaw disabled")


# Rerun the layout optimization with geometric yaw enabled
print("\nReoptimizing with geometric yaw enabled.")
fmodel.set(layout_x=layout_x, layout_y=layout_y)
layout_opt = LayoutOptimizationScipy(
    fmodel, boundaries, min_dist=2 * D, enable_geometric_yaw=True, optOptions={"maxiter": maxiter}
)

# Run the optimization
np.random.seed(0)
sol = layout_opt.optimize()

# Get the resulting improvement in AEP
print("... calcuating improvement in AEP")

fmodel.set(yaw_angles=np.zeros_like(layout_opt.yaw_angles))
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() / 1e6

percent_gain = 100 * (opt_aep - base_aep) / base_aep

# Print and plot the results
print(f"Optimal layout: {sol}")
print(
    f"Optimal layout improves AEP by {percent_gain:.1f}% "
    f"from {base_aep:.1f} MWh to {opt_aep:.1f} MWh"
)
layout_opt.plot_layout_opt_results()
ax = plt.gca()
fig = plt.gcf()
sm = ax.tricontourf(x_locs, y_locs, speed_multipliers[0], cmap="coolwarm")
fig.colorbar(sm, ax=ax, label="Speed multiplier")
ax.legend(["_Optimization boundary", "Initial layout", "Optimized layout"])
ax.set_title("Geometric yaw enabled")

print(
    "Turbine geometric yaw angles for wind direction {0:.2f}".format(wind_directions[1])
    + " and wind speed {0:.2f} m/s:".format(wind_speeds[0]),
    f"{layout_opt.yaw_angles[1, :]}",
)

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

=====================================================
Optimizing turbine layout...
Number of parameters to optimize =  6
=====================================================

  NIT    FC           OBJFUN            GNORM
    1     8    -2.212410E+00     2.629932E+00

    2    15    -2.240060E+00     3.547760E+00

floris.floris_model.FlorisModel WARNING Some velocities at the rotor are negative.

floris.floris_model.FlorisModel WARNING Some velocities at the rotor are negative.

    3    22    -3.765389E+00     3.855406E+00

floris.floris_model.FlorisModel WARNING Some velocities at the rotor are negative.

floris.floris_model.FlorisModel WARNING Some velocities at the rotor are negative.

    5    29    -3.637278E+00     1.814021E+02

floris.floris_model.FlorisModel WARNING Some velocities at the rotor are negative.

floris.floris_model.FlorisModel WARNING Some velocities at the rotor are negative.

floris.floris_model.FlorisModel WARNING Some velocities at the rotor are negative.

floris.floris_model.FlorisModel WARNING Some velocities at the rotor are negative.

    6    46    -3.637278E+00     1.814021E+02

floris.floris_model.FlorisModel WARNING Some velocities at the rotor are negative.

floris.floris_model.FlorisModel WARNING Some velocities at the rotor are negative.

floris.floris_model.FlorisModel WARNING Some velocities at the rotor are negative.

floris.floris_model.FlorisModel WARNING Some velocities at the rotor are negative.

    7    63    -3.637233E+00     1.814021E+02

floris.floris_model.FlorisModel WARNING Some velocities at the rotor are negative.

floris.floris_model.FlorisModel WARNING Some velocities at the rotor are negative.

    8    80    -3.748231E+00     1.814021E+02

floris.floris_model.FlorisModel WARNING Some velocities at the rotor are negative.

floris.floris_model.FlorisModel WARNING Some velocities at the rotor are negative.

floris.floris_model.FlorisModel WARNING Some velocities at the rotor are negative.

    9    97    -3.795270E+00     1.813920E+02

floris.floris_model.FlorisModel WARNING Some velocities at the rotor are negative.

floris.floris_model.FlorisModel WARNING Some velocities at the rotor are negative.

floris.floris_model.FlorisModel WARNING Some velocities at the rotor are negative.

floris.floris_model.FlorisModel WARNING Some velocities at the rotor are negative.

floris.floris_model.FlorisModel WARNING Some velocities at the rotor are negative.

   10   114    -7.436518E+00     1.813941E+02

   12   121    -7.456515E-01     6.222679E+02

floris.floris_model.FlorisModel WARNING Some velocities at the rotor are negative.

floris.floris_model.FlorisModel WARNING Some velocities at the rotor are negative.

floris.floris_model.FlorisModel WARNING Some velocities at the rotor are negative.

   13   138    -7.456521E-01     6.214802E+02

floris.floris_model.FlorisModel WARNING Some velocities at the rotor are negative.

floris.floris_model.FlorisModel WARNING Some velocities at the rotor are negative.

floris.floris_model.FlorisModel WARNING Some velocities at the rotor are negative.

floris.logging_manager.LoggingManager WARNING The calculated flow field contains points outside of the the user-defined heterogeneous inflow bounds. For these points, the interpolated value has been filled with the freestream wind speed. If this is not the desired behavior, the user will need to expand the heterogeneous inflow bounds to fully cover the calculated flow field area.

   14   155    -7.456514E-01     6.207013E+02

floris.floris_model.FlorisModel WARNING Some velocities at the rotor are negative.

floris.floris_model.FlorisModel WARNING Some velocities at the rotor are negative.

floris.floris_model.FlorisModel WARNING Some velocities at the rotor are negative.

floris.logging_manager.LoggingManager WARNING The calculated flow field contains points outside of the the user-defined heterogeneous inflow bounds. For these points, the interpolated value has been filled with the freestream wind speed. If this is not the desired behavior, the user will need to expand the heterogeneous inflow bounds to fully cover the calculated flow field area.

   15   172    -7.456514E-01     6.199276E+02

floris.floris_model.FlorisModel WARNING Some velocities at the rotor are negative.

floris.floris_model.FlorisModel WARNING Some velocities at the rotor are negative.

floris.floris_model.FlorisModel WARNING Some velocities at the rotor are negative.

floris.logging_manager.LoggingManager WARNING The calculated flow field contains points outside of the the user-defined heterogeneous inflow bounds. For these points, the interpolated value has been filled with the freestream wind speed. If this is not the desired behavior, the user will need to expand the heterogeneous inflow bounds to fully cover the calculated flow field area.

   16   189    -7.456514E-01     6.191569E+02

floris.floris_model.FlorisModel WARNING Some velocities at the rotor are negative.

floris.floris_model.FlorisModel WARNING Some velocities at the rotor are negative.

floris.floris_model.FlorisModel WARNING Some velocities at the rotor are negative.

   17   206    -7.436518E+00     6.183878E+02

floris.floris_model.FlorisModel WARNING Some velocities at the rotor are negative.

floris.floris_model.FlorisModel WARNING Some velocities at the rotor are negative.

floris.floris_model.FlorisModel WARNING Some velocities at the rotor are negative.

floris.logging_manager.LoggingManager WARNING The calculated flow field contains points outside of the the user-defined heterogeneous inflow bounds. For these points, the interpolated value has been filled with the freestream wind speed. If this is not the desired behavior, the user will need to expand the heterogeneous inflow bounds to fully cover the calculated flow field area.

   18   223    -7.456514E-01     6.183916E+02

floris.floris_model.FlorisModel WARNING Some velocities at the rotor are negative.

floris.floris_model.FlorisModel WARNING Some velocities at the rotor are negative.

floris.floris_model.FlorisModel WARNING Some velocities at the rotor are negative.

   19   240    -7.436518E+00     6.176159E+02

floris.floris_model.FlorisModel WARNING Some velocities at the rotor are negative.

floris.floris_model.FlorisModel WARNING Some velocities at the rotor are negative.

floris.floris_model.FlorisModel WARNING Some velocities at the rotor are negative.

   20   256    -7.379377E+00     6.176204E+02
Inequality constraints incompatible    (Exit mode 4)
            Current function value: -7.379377455018026
            Iterations: 20
            Function evaluations: 256
            Gradient evaluations: 18
Optimization complete.
... calcuating improvement in AEP
Optimal layout: [[np.float64(1504.18580302199), np.float64(1504.1858030205085), np.float64(1504.1858030200362)], [np.float64(0.0005166359525263718), np.float64(0.0005166359727981101), np.float64(0.0005166359733195485)]]
Optimal layout improves AEP by 637.9% from 6198.3 MWh to 45739.8 MWh

Reoptimizing with geometric yaw enabled.
=====================================================
Optimizing turbine layout...
Number of parameters to optimize =  6
=====================================================

  NIT    FC           OBJFUN            GNORM
    1     8    -2.365430E+00     2.961907E+00

    2    15    -2.508748E+00     3.815898E+00

    3    22    -4.348248E+00     3.896698E+00

    4    29    -2.382590E+00     2.495638E+02

    5    36    -2.487043E+00     3.904223E+00

    6    43    -3.678105E+00     8.572127E+00

    7    52    -2.595233E+00     4.406446E+00

    8    59    -2.598515E+00     4.410221E+00

    9    66    -2.599785E+00     4.435161E+00

   10    73    -2.598476E+00     4.770662E+00

   11    90    -2.599786E+00     4.463337E+00

   13    97    -2.598885E+00     4.463341E+00

   14   114    -2.386362E+00     4.770614E+00

   15   131    -2.599782E+00     4.462964E+00

   16   138    -2.599790E+00     4.462966E+00

   17   145    -2.599830E+00     4.462962E+00

   18   152    -2.600029E+00     4.462943E+00

   19   159    -2.601020E+00     4.462849E+00

   20   166    -2.601019E+00     4.764651E+00

   21   183    -2.601025E+00     4.462385E+00

   22   189    -2.601025E+00     4.462390E+00
Optimization terminated successfully    (Exit mode 0)
            Current function value: -2.6010250276521445
            Iterations: 22
            Function evaluations: 189
            Gradient evaluations: 21
Optimization complete.
... calcuating improvement in AEP
Optimal layout: [[np.float64(649.8718401364441), np.float64(1259.6506736338688), np.float64(1511.9999999993222)], [np.float64(2.429451220684333e-18), np.float64(2.004487092314508e-18), np.float64(0.09999999973226106)]]
Optimal layout improves AEP by 160.1% from 6198.3 MWh to 16122.0 MWh
Turbine geometric yaw angles for wind direction 270.00 and wind speed 8.00 m/s: [24.18704639 27.59438202  0.        ]