Add yaw optimization with SciPy (#306)

* Move generalized code in SR to base yaw optimization code. Also rename _optimize to optimize for SR optimizer.

* Add pandas df creation to base yaw opt class. Also, improve normalization functionality in preparation for SciPy fmin

* Implement basic SciPy based yaw optimizer

* Add  example to compare SciPy with SR optimizer

* Reintroduce option to verify convergence for yaw optimizer

* Remove cluster_turbines from yaw optimization tools for the time being. Probably not needed / no need to be reimplemented in v3.

* preliminary changes to allow wind speed range in yaw optimizer

* Add compatibility with range of wind speeds to yaw optimizer

* Add range of wind speeds compatibility to SciPy yaw optimizer

* Add yaw optimization example for multiple wind speeds

examples/08_opt_yaw_single_ws.py

@@ -39,7 +39,7 @@
 
 # Initialize optimizer object and run optimization using the Serial-Refine method
 yaw_opt = YawOptimizationSR(fi)#, exploit_layout_symmetry=False)
-df_opt = yaw_opt._optimize()
+df_opt = yaw_opt.optimize()
 
 print("Optimization results:")
 print(df_opt)

examples/09_opt_yaw_multiple_ws.py

@@ -0,0 +1,100 @@
+# Copyright 2022 NREL
+
+# Licensed under the Apache License, Version 2.0 (the "License"); you may not
+# use this file except in compliance with the License. You may obtain a copy of
+# the License at http://www.apache.org/licenses/LICENSE-2.0
+
+# Unless required by applicable law or agreed to in writing, software
+# distributed under the License is distributed on an "AS IS" BASIS, WITHOUT
+# WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. See the
+# License for the specific language governing permissions and limitations under
+# the License.
+
+# See https://floris.readthedocs.io for documentation
+
+import numpy as np
+import matplotlib.pyplot as plt
+from floris.tools import FlorisInterface
+from floris.tools.optimization.yaw_optimization.yaw_optimizer_sr import (
+    YawOptimizationSR,
+)
+
+"""
+This example demonstrates how to perform a yaw optimization for multiple wind directions and multiple wind speeds.
+
+First, we initialize our Floris Interface, and then generate a 3 turbine wind farm. Next, we create the yaw optimization object `yaw_opt` and perform the optimization using the SerialRefine method. Finally, we plot the results.
+"""
+
+# Load the default example floris object
+fi = FlorisInterface("inputs/gch.yaml") # GCH model matched to the default "legacy_gauss" of V2
+# fi = FlorisInterface("inputs/cc.yaml") # New CumulativeCurl model
+
+# Reinitialize as a 3-turbine farm with range of WDs and 1 WS
+D = 126.0 # Rotor diameter for the NREL 5 MW
+fi.reinitialize(
+    layout=[[0.0, 5 * D, 10 * D], [0.0, 0.0, 0.0]],
+    wind_directions=np.arange(0.0, 360.0, 3.0), 
+    wind_speeds=np.arange(2.0, 18.0, 1.0),
+)
+
+# Initialize optimizer object and run optimization using the Serial-Refine method
+yaw_opt = YawOptimizationSR(fi)
+df_opt = yaw_opt.optimize()
+
+print("Optimization results:")
+print(df_opt)
+
+# Split out the turbine results
+for t in range(3):
+    df_opt['t%d' % t] = df_opt.yaw_angles_opt.apply(lambda x: x[t])
+
+# Show the results: optimal yaw angles
+fig, axarr = plt.subplots(nrows=4, ncols=4, sharex=True, sharey=True, figsize=(10, 8))
+jj = 0
+for ii, ws in enumerate(fi.floris.flow_field.wind_speeds):
+    xi = np.remainder(ii, 4)
+    if ((ii > 0) & (xi == 0)):
+        jj += 1
+    ax = axarr[np.remainder(ii, 4)][jj]
+    ids = (df_opt.wind_speed == ws)
+    wd = df_opt.loc[ids, "wind_direction"]
+    for t in range(3):
+        yaw_opt = df_opt.loc[ids, "t{:d}".format(t)]
+        ax.plot(wd, yaw_opt, label='Turbine {:d}'.format(t))
+    ax.set_title("Wind speed = {:.1f} m/s".format(ws), size=10)
+    if ((ii == 0) & (jj == 0)):
+        ax.legend()
+    ax.grid(True)
+    if jj == 0:
+        ax.set_ylabel('Yaw angle (deg)', size=10)
+    if xi == 3:
+        axarr[xi][jj].set_xlabel('Wind Direction (deg)', size=10)
+
+    plt.tight_layout()
+
+# Show the results: baseline and optimized farm power
+fig, axarr = plt.subplots(nrows=4, ncols=4, sharex=True, figsize=(10, 8))
+jj = 0
+for ii, ws in enumerate(fi.floris.flow_field.wind_speeds):
+    xi = np.remainder(ii, 4)
+    if ((ii > 0) & (xi == 0)):
+        jj += 1
+    ax = axarr[np.remainder(ii, 4)][jj]
+    ids = (df_opt.wind_speed == ws)
+    wd = df_opt.loc[ids, "wind_direction"]
+    power_baseline = df_opt.loc[ids, "farm_power_baseline"]
+    power_opt = df_opt.loc[ids, "farm_power_opt"]
+    ax.plot(wd, power_baseline / 1e6, color='k', label='Baseline')
+    ax.plot(wd, power_opt / 1e6, color='r', label='Optimized')
+    ax.set_title("Wind speed = {:.1f} m/s".format(ws), size=10)
+    if ((ii == 0) & (jj == 0)):
+        ax.legend()
+    ax.grid(True)
+    if jj == 0:
+        ax.set_ylabel('Farm Power (MW)', size=10)
+    if xi == 3:
+        axarr[xi][jj].set_xlabel('Wind Direction (deg)', size=10)
+
+    plt.tight_layout()
+
+plt.show()

examples/10_compare_yaw_optimizers.py

@@ -0,0 +1,86 @@
+# Copyright 2022 NREL
+
+# Licensed under the Apache License, Version 2.0 (the "License"); you may not
+# use this file except in compliance with the License. You may obtain a copy of
+# the License at http://www.apache.org/licenses/LICENSE-2.0
+
+# Unless required by applicable law or agreed to in writing, software
+# distributed under the License is distributed on an "AS IS" BASIS, WITHOUT
+# WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. See the
+# License for the specific language governing permissions and limitations under
+# the License.
+
+# See https://floris.readthedocs.io for documentation
+
+from time import perf_counter as timerpc
+import numpy as np
+import matplotlib.pyplot as plt
+from floris.tools import FlorisInterface
+from floris.tools.optimization.yaw_optimization.yaw_optimizer_scipy import (
+    YawOptimizationScipy
+)
+from floris.tools.optimization.yaw_optimization.yaw_optimizer_sr import (
+    YawOptimizationSR
+)
+
+
+"""
+This example compares the SciPy-based yaw optimizer with the new Serial-Refine optimizer.
+
+First, we initialize our Floris Interface, and then generate a 3 turbine wind farm. Next, we create two yaw optimization objects, `yaw_opt_sr` and `yaw_opt_scipy` for the Serial-Refine and SciPy methods, respectively. We then perform the optimization using both methods. Finally, we compare the time it took to find the optimal angles and plot the optimal yaw angles and resulting wind farm powers.
+"""
+
+# Load the default example floris object
+fi = FlorisInterface("inputs/gch.yaml") # GCH model matched to the default "legacy_gauss" of V2
+# fi = FlorisInterface("inputs/cc.yaml") # New CumulativeCurl model
+
+# Reinitialize as a 3-turbine farm with range of WDs and 1 WS
+D = 126.0 # Rotor diameter for the NREL 5 MW
+fi.reinitialize(
+    layout=[[0.0, 5 * D, 10 * D], [0.0, 0.0, 0.0]],
+    wind_directions=np.arange(0.0, 360.0, 3.0), 
+    wind_speeds=[8.0],
+)
+
+print("Performing optimizations with SciPy...")
+start_time = timerpc()
+yaw_opt_scipy = YawOptimizationScipy(fi)
+df_opt_scipy = yaw_opt_scipy.optimize()
+time_scipy = timerpc() - start_time
+
+print("Performing optimizations with Serial Refine...")
+start_time = timerpc()
+yaw_opt_sr = YawOptimizationSR(fi)
+df_opt_sr = yaw_opt_sr.optimize()
+time_sr = timerpc() - start_time
+
+# Print time spent
+print("\n Time spent, Serial Refine: {:.2f} s.".format(time_sr))
+print(" Time spent, SciPy (SLSQP): {:.2f} s.\n".format(time_scipy))
+
+# Split out the turbine results
+yaw_angles_opt_sr = np.vstack(df_opt_sr.yaw_angles_opt)
+yaw_angles_opt_scipy = np.vstack(df_opt_scipy.yaw_angles_opt)
+
+# Yaw results
+for t in range(3):
+    fig, ax = plt.subplots()
+    ax.plot(df_opt_sr.wind_direction, yaw_angles_opt_sr[:, t],color='r',label='Serial Refine')
+    ax.plot(df_opt_scipy.wind_direction, yaw_angles_opt_scipy[:, t],'--', color='g', label='SciPy')
+    ax.grid(True)
+    ax.set_ylabel('Yaw Offset (deg')
+    ax.legend()
+    ax.grid(True)
+    ax.set_title("Turbine {:d}".format(t))
+
+# Power results
+fig, ax = plt.subplots()
+ax.plot(df_opt_sr.wind_direction,df_opt_sr.farm_power_baseline,color='k',label='Baseline')
+ax.plot(df_opt_sr.wind_direction,df_opt_sr.farm_power_opt,color='r',label='Optimized, Serial Refine')
+ax.plot(df_opt_scipy.wind_direction,df_opt_scipy.farm_power_opt,'--', color='g', label='Optimized, SciPy')
+ax.set_ylabel('Wind Farm Power (W)')
+ax.set_xlabel('Wind Direction (deg)')
+ax.legend()
+ax.grid(True)
+
+plt.show()