# Example: Optimize Layout#

"""Example: Optimize Layout
This example shows a simple layout optimization using the python module Scipy, optimizing for both
annual energy production (AEP) and annual value production (AVP).

First, a 4 turbine array is optimized such that the layout of the turbine produces the
highest AEP based on the given wind resource. The turbines
are constrained to a square boundary and a random wind resource is supplied. The results
of the optimization show that the turbines are pushed to near the outer corners of the boundary,
which, given the generally uniform wind rose, makes sense in order to maximize the energy
production by minimizing wake interactions.

Next, with the same boundary, the same 4 turbine array is optimized to maximize AVP instead of AEP,
using the value table defined in the WindRose object, where value represents the value of the
energy produced for a given wind condition (e.g., the price of electricity). In this example, value
is defined to be significantly higher for northerly and southerly wind directions, and zero when
the wind is from the east or west. Because the value is much higher when the wind is from the north
or south, the turbines are spaced apart roughly evenly in the x direction while being relatively
close in the y direction to avoid wake interactions for northerly and southerly winds. Although the
layout results in large wake losses when the wind is from the east or west, these losses do not
significantly impact the objective function because of the low value for those wind directions.
"""

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,
)

# Define scipy optimization parameters
opt_options = {
"maxiter": 20,
"disp": True,
"iprint": 2,
"ftol": 1e-12,
"eps": 0.05,
}

# Initialize the FLORIS interface fi
fmodel = FlorisModel('../inputs/gch.yaml')

# Setup 72 wind directions with a 1 wind speed and frequency distribution
wind_directions = np.arange(0, 360.0, 5.0)
wind_speeds = np.array([8.0])

# Shape random frequency distribution to match number of wind directions and wind speeds
freq_table = np.zeros((len(wind_directions), len(wind_speeds)))
np.random.seed(1)
freq_table[:,0] = (np.abs(np.sort(np.random.randn(len(wind_directions)))))
freq_table = freq_table / freq_table.sum()

# Define the value table such that the value of the energy produced is
# significantly higher when the wind direction is close to the north or
# south, and zero when the wind is from the east or west. Here, value is
# given a mean value of 25 USD/MWh.
value_table = 25*value_table/np.mean(value_table)
value_table = value_table.reshape((len(wind_directions),1))

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

fmodel.set(wind_data=wind_rose)

# The boundaries for the turbines, specified as vertices
boundaries = [(0.0, 0.0), (0.0, 1000.0), (1000.0, 1000.0), (1000.0, 0.0), (0.0, 0.0)]

# Set turbine locations to 4 turbines in a rectangle
D = 126.0 # rotor diameter for the NREL 5MW
layout_x = [0, 0, 6 * D, 6 * D]
layout_y = [0, 4 * D, 0, 4 * D]
fmodel.set(layout_x=layout_x, layout_y=layout_y)

# Setup the optimization problem to maximize AEP instead of value
layout_opt = LayoutOptimizationScipy(fmodel, boundaries, optOptions=opt_options)

# Run the optimization
sol = layout_opt.optimize()

# Get the resulting improvement in AEP
print('... calculating 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()

# reset to the original layout
fmodel.set(layout_x=layout_x, layout_y=layout_y)

# Now set up the optimization problem to maximize annual value production (AVP)
# using the value table provided in the WindRose object.
layout_opt = LayoutOptimizationScipy(fmodel, boundaries, optOptions=opt_options, use_value=True)

# Run the optimization
sol = layout_opt.optimize()

# Get the resulting improvement in AVP
print('... calculating improvement in annual value production (AVP)')
fmodel.run()
base_avp = fmodel.get_farm_AVP() / 1e6
fmodel.set(layout_x=sol[0], layout_y=sol[1])
fmodel.run()
opt_avp = fmodel.get_farm_AVP() / 1e6

percent_gain = 100 * (opt_avp - base_avp) / base_avp

# Print and plot the results
print(f'Optimal layout: {sol}')
print(
f'Optimal layout improves AVP by {percent_gain:.1f}% '
f'from {base_avp:.1f} dollars to {opt_avp:.1f} dollars'
)
layout_opt.plot_layout_opt_results()

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

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

  NIT    FC           OBJFUN            GNORM
1    10    -1.014190E+00     1.590506E-01

    2    19    -1.033878E+00     1.487938E-01

    3    28    -1.036191E+00     9.303427E-02

    4    37    -1.042592E+00     7.892545E-02

    5    46    -1.043000E+00     6.375930E-02

    6    55    -1.044561E+00     6.598999E-02

    7    64    -1.044745E+00     6.467957E-02

    8    73    -1.044633E+00     6.652021E-02

    9    83    -1.044734E+00     6.623317E-02

   10    94    -1.044523E+00     6.610063E-02

   11   106    -1.044640E+00     6.602156E-02

   12   123    -1.044430E+00     6.624512E-02

   13   142    -1.044429E+00     6.623653E-02

   14   159    -1.044483E+00     6.605804E-02

   15   178    -1.044059E+00     6.607263E-02

   16   197    -1.044841E+00     6.575627E-02

   17   206    -1.044839E+00     6.601919E-02

   18   218    -1.044838E+00     6.605415E-02

Optimization terminated successfully    (Exit mode 0)
Current function value: -1.0448379988083718
Iterations: 18
Function evaluations: 228
Optimization complete.
... calculating improvement in AEP
Optimal layout: [[2.7996863251832383e-18, 126.15002186203341, 794.1056557846636, 1000.0], [1.192172230131564e-14, 999.9698994858535, 7.039764341504759e-24, 1000.0]]
Optimal layout improves AEP by 4.5% from 56299.2 MWh to 58823.5 MWh

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

  NIT    FC           OBJFUN            GNORM
1    10    -1.177117E+00     7.500783E-01

    2    19    -1.159173E+00     5.490016E-01

    3    29    -1.219949E+00     3.415379E-01

    4    38    -1.228863E+00     2.197285E-01

    5    47    -1.237300E+00     2.748039E-01

    6    56    -1.251563E+00     2.651905E-01

    7    65    -1.252709E+00     1.048426E-01

    8    74    -1.252392E+00     4.818055E-02

    9    83    -1.252482E+00     7.784796E-02

   10    93    -1.252233E+00     7.733065E-02

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

   11   112    -1.246641E+00     7.733127E-02

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

   12   131    -1.246692E+00     7.733173E-02

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

   13   150    -1.137010E+00     7.733210E-02

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

   14   169    -1.160673E+00     7.733258E-02

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

   15   188    -1.245428E+00     7.733298E-02

   16   207    -1.251980E+00     7.733323E-02

   17   226    -1.252082E+00     7.733350E-02

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

   18   245    -1.243072E+00     7.733381E-02

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

   19   264    -1.213373E+00     7.733405E-02

   20   282    -1.252434E+00     7.733430E-02
Iteration limit reached    (Exit mode 9)
Current function value: -1.252433893445541
Iterations: 20
Function evaluations: 282