Merge pull request #921 from NREL/develop

FLORIS v4.1

.gitignore

@@ -40,3 +40,9 @@ examples/SLSQP.out
 
 # Log files
 *.log
+
+# Temp files in convert test
+tests/v3_to_v4_convert_test/convert_turbine_v3_to_v4.py
+tests/v3_to_v4_convert_test/convert_floris_input_v3_to_v4.py
+tests/v3_to_v4_convert_test/gch_v4.yaml
+tests/v3_to_v4_convert_test/nrel_5MW_v3_v4.yaml

README.md

@@ -3,7 +3,7 @@
 FLORIS is a controls-focused wind farm simulation software incorporating
 steady-state engineering wake models into a performance-focused Python
 framework. It has been in active development at NREL since 2013 and the latest
-release is [FLORIS v4.0.1](https://github.com/NREL/floris/releases/latest).
+release is [FLORIS v4.1](https://github.com/NREL/floris/releases/latest).
 Online documentation is available at https://nrel.github.io/floris.
 
 The software is in active development and engagement with the development team
@@ -79,7 +79,7 @@ PACKAGE CONTENTS
     wind_data
 
 VERSION
-    4
+    4.1
 
 FILE
     ~/floris/floris/__init__.py

docs/_toc.yml

@@ -14,9 +14,11 @@ parts:
     - file: intro_concepts
     - file: advanced_concepts
     - file: wind_data_user
+    - file: heterogeneous_map
     - file: floating_wind_turbine
     - file: turbine_interaction
     - file: operation_models_user
+    - file: layout_optimization
     - file: input_reference_main
     - file: input_reference_turbine
     - file: examples

docs/empirical_gauss_model.md

@@ -172,7 +172,7 @@ The effect of AWC is represented by updating the
 wake-induced mixing term as follows:
 
 $$ \text{WIM}_j = \sum_{i \in T^{\text{up}}(j)} \frac{A_{ij} a_i} {(x_j - x_i)/D_i} +
-\frac{\beta_{j}^{p}{d}$$
+\frac{\beta_{j}^{p}}{d}$$
 
 where $\beta_{j}$ is the AWC amplitude of turbine $j$, and the exponent $p$ and
 denominator $d$ are tuning parameters that can be set in the `emgauss.yaml` file with

docs/layout_optimization.md

@@ -0,0 +1,81 @@
+
+(layout_optimization)=
+# Layout optimization
+
+The FLORIS package provides layout optimization tools to place turbines within a specified
+boundary area to optimize annual energy production (AEP) or wind plant value. Layout
+optimizers accept an instantiated `FlorisModel` and alter the turbine layouts in order to
+improve the objective function value (AEP or value).
+
+## Background
+
+Layout optimization entails placing turbines in a wind farm in a configuration that maximizes
+an objective function, usually the AEP. Turbines are moved to minimize their wake interactions
+in the most dominant wind directions, while respecting the boundaries of the area for turbine
+placement as well as minimum distance requirements between neighboring turbines.
+
+Mathematically, we represent this as a (nonconvex) optimization problem.
+Let $x = \{x_i\}_{i=1,\dots,N}$, $x_i \in \mathbb{R}^2$ represent the set of
+coordinates of turbines within a farm (that is, $x_i$ represents the $(X, Y)$
+location of turbine $i$). Further, let $R \subset \mathbb{R}^2$ be a closed
+region in which to place turbines. Finally, let $d$ represent the minimum
+allowable distance between two turbines. Then, the layout optimization problem
+is expressed as
+
+$$
+\begin{aligned}
+\underset{x}{\text{maximize}} & \:\: f(x)\\
+\text{subject to} & \:\: x \subset R \\
+& \:\: ||x_i - x_j|| \geq d \:\: \forall j = 1,\dots,N, j\neq i
+\end{aligned}
+$$
+
+Here, $||\cdot||$ denotes the Euclidean norm, and $f(x)$ is the cost function to be maximized.
+
+When maximizing the AEP, $f = \sum_w P(w, x)p_W(w)$, where $w$ is the wind condition bin
+(e.g., wind speed, wind direction pair); $P(w, x)$ is the power produced by the wind farm in
+condition $w$ with layout $x$; and $p_W(w)$ is the annual frequency of occurrence of
+condition $w$.
+
+Layout optimizers take iterative approaches to solving the layout optimization problem
+specified above. Optimization routines available in FLORIS are described below.
+
+## Scipy layout optimization
+The `LayoutOptimizationScipy` class is built around `scipy.optimize`s `minimize`
+routine, using the `SLSQP` solver by default. Options for adjusting
+`minimize`'s behavior are exposed to the user with the `optOptions` argument.
+Other options include enabling fast wake steering at each layout optimizer
+iteration with the `enable_geometric_yaw` argument, and switch from AEP
+optimization to value optimization with the `use_value` argument.
+
+## Genetic random search layout optimization
+The `LayoutOptimizationRandomSearch` class is a custom optimizer designed specifically for
+layout optimization via random perturbations of the turbine locations. It is designed to have
+the following features:
+- Robust to complex wind conditions and complex boundaries, including disjoint regions for
+turbine placement
+- Easy to parallelize and wrapped in a genetic algorithm for propagating candidate solutions
+- Simple to set up and tune for non-optimization experts
+- Set up to run cheap constraint checks prior to more expensive objective function evaluations
+to accelerate optimization
+
+The algorithm, described in full in an upcoming paper that will be linked here when it is
+publicly available, moves a random turbine and random distance in a random direction; checks
+that constraints are satisfied; evaluates the objective function (AEP or value); and then
+commits to the move if there is an objective function improvement. The main tuning parameter
+is the probability mass function for the random movement distance, which is a dictionary to be
+passed to the `distance_pmf` argument.
+
+The `distance_pmf` dictionary should contain two keys, each containing a 1D array of equal
+length: `"d"`, which specifies the perturbation distance _in units of the rotor diameter_,
+and `"p"`, which specifies the probability that the corresponding perturbation distance is
+chosen at any iteration of the random search algorithm. The `distance_pmf` can therefore be
+used to encourage or discourage more aggressive or more conservative movements, and to enable
+or disable jumps between disjoint regions for turbine placement.
+
+The figure below shows an example of the optimized layout of a farm using the GRS algorithm, with
+the black dots indicating the initial layout; red dots indicating the final layout; and blue
+shading indicating wind speed heterogeneity (lighter shade is lower wind speed, darker shade is
+higher wind speed). The progress of each of the genetic individuals in the optimization process is
+shown in the right-hand plot.
+![](plot_complex_docs.png)

docs/operation_models_user.ipynb

@@ -380,6 +380,92 @@
    "cell_type": "markdown",
    "id": "25f9c86c",
    "metadata": {},
+   "source": [
+    "### Peak shaving model\n",
+    "\n",
+    "User-level name: `\"peak-shaving\"`\n",
+    "\n",
+    "Underlying class: `PeakShavingTurbine`\n",
+    "\n",
+    "Required data on `power_thrust_table`:\n",
+    "- `ref_air_density` (scalar)\n",
+    "- `ref_tilt` (scalar)\n",
+    "- `wind_speed` (list)\n",
+    "- `power` (list)\n",
+    "- `thrust_coefficient` (list)\n",
+    "- `peak_shaving_fraction` (scalar)\n",
+    "- `peak_shaving_TI_threshold` (scalar)\n",
+    "\n",
+    "The `\"peak-shaving\"` operation model allows users to implement peak shaving, where the thrust\n",
+    "of the wind turbine is reduced from the nominal curve near rated to reduce unwanted structural\n",
+    "loading. Peak shaving here is implemented here by reducing the thrust by a fixed fraction from\n",
+    "the peak thrust on the nominal thrust curve, as specified by `peak_shaving_fraction`.This only\n",
+    "affects wind speeds near the peak in the thrust\n",
+    "curve (usually near rated wind speed), as thrust values away from the peak will be below the\n",
+    "fraction regardless. Further, peak shaving is only applied if the turbulence intensity experienced\n",
+    "by the turbine meets the `peak_shaving_TI_threshold`. To apply peak shaving in all wind conditions,\n",
+    "`peak_shaving_TI_threshold` may be set to zero.\n",
+    "\n",
+    "When the turbine is peak shaving to reduce thrust, the power output is updated accordingly. Letting\n",
+    "$C_{T}$ represent the thrust coefficient when peak shaving (at given wind speed), and $C_{T}'$\n",
+    "represent the thrust coefficient that the turbine would be operating at under nominal control, then\n",
+    "the power $P$ due to peak shaving (compared to the power $P'$ available under nominal control) is \n",
+    "computed (based on actuator disk theory) as\n",
+    "\n",
+    "$$ P = \\frac{C_T (1 - a)}{C_T' (1 - a')} P'$$\n",
+    "\n",
+    "where $a$ (respectively, $a'$) is the axial induction factor corresponding to $C_T$\n",
+    "(respectively, $C_T'$), computed using the usual relationship from actuator disk theory,\n",
+    "i.e. the lesser solution to $C_T=4a(1-a)$.\n"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "1eff05f3",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "# Set up the FlorisModel\n",
+    "fmodel = FlorisModel(\"../examples/inputs/gch.yaml\")\n",
+    "fmodel.set(\n",
+    "    layout_x=[0.0], layout_y=[0.0],\n",
+    "    wind_data=TimeSeries(\n",
+    "        wind_speeds=wind_speeds,\n",
+    "        wind_directions=np.ones(100) * 270.0,\n",
+    "        turbulence_intensities=0.2 # Higher than threshold value of 0.1\n",
+    "    )\n",
+    ")\n",
+    "fmodel.reset_operation()\n",
+    "fmodel.set_operation_model(\"simple\")\n",
+    "fmodel.run()\n",
+    "powers_base = fmodel.get_turbine_powers()/1000\n",
+    "thrust_coefficients_base = fmodel.get_turbine_thrust_coefficients()\n",
+    "fmodel.set_operation_model(\"peak-shaving\")\n",
+    "fmodel.run()\n",
+    "powers_peak_shaving = fmodel.get_turbine_powers()/1000\n",
+    "thrust_coefficients_peak_shaving = fmodel.get_turbine_thrust_coefficients()\n",
+    "\n",
+    "fig, ax = plt.subplots(2,1,sharex=True)\n",
+    "ax[0].plot(wind_speeds, thrust_coefficients_base, label=\"Without peak shaving\", color=\"black\")\n",
+    "ax[0].plot(wind_speeds, thrust_coefficients_peak_shaving, label=\"With peak shaving\", color=\"C0\")\n",
+    "ax[1].plot(wind_speeds, powers_base, label=\"Without peak shaving\", color=\"black\")\n",
+    "ax[1].plot(wind_speeds, powers_peak_shaving, label=\"With peak shaving\", color=\"C0\")\n",
+    "\n",
+    "ax[1].grid()\n",
+    "ax[0].grid()\n",
+    "ax[0].legend()\n",
+    "ax[0].set_ylabel(\"Thrust coefficient [-]\")\n",
+    "ax[1].set_xlabel(\"Wind speed [m/s]\")\n",
+    "ax[1].set_ylabel(\"Power [kW]\")"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "92912bf7",
+   "metadata": {},
+   "outputs": [],
    "source": []
   }
  ],

examples/002_visualizations.py

@@ -69,8 +69,10 @@
 # we show the horizontal plane at hub height, further examples are provided within
 # the examples_visualizations folder
 
-# For flow visualizations, the FlorisModel must be set to run a single condition
-# (n_findex = 1)
+# For flow visualizations, the FlorisModel may be set to run a single condition
+# (n_findex = 1). Otherwise, the user may set multiple conditions and then use
+# the findex_for_viz keyword argument to calculate_horizontal_plane to specify which
+# flow condition to visualize.
 fmodel.set(wind_speeds=[8.0], wind_directions=[290.0], turbulence_intensities=[0.06])
 horizontal_plane = fmodel.calculate_horizontal_plane(
     x_resolution=200,

examples/examples_control_types/005_peak_shaving.py

@@ -0,0 +1,124 @@
+"""Example of using the peak-shaving turbine operation model.
+
+This example demonstrates how to use the peak-shaving operation model in FLORIS.
+The peak-shaving operation model allows the user to a thrust reduction near rated wind speed to
+reduce loads on the turbine. The power is reduced accordingly, and wind turbine wakes
+are shallower due to the reduced thrust.
+
+"""
+
+import matplotlib.pyplot as plt
+import numpy as np
+
+from floris import FlorisModel, TimeSeries
+
+
+fmodel = FlorisModel("../inputs/gch.yaml")
+fmodel.set(layout_x=[0, 1000.0], layout_y=[0.0, 0.0])
+wind_speeds = np.linspace(0, 30, 100)
+fmodel.set(
+    wind_data=TimeSeries(
+        wind_directions=270 * np.ones_like(wind_speeds),
+        wind_speeds=wind_speeds,
+        turbulence_intensities=0.10, # High enough to engage peak shaving
+    )
+)
+
+# Start with "normal" operation under the simple turbine operation model
+fmodel.set_operation_model("simple")
+fmodel.run()
+powers_base = fmodel.get_turbine_powers()/1000
+thrust_coefficients_base = fmodel.get_turbine_thrust_coefficients()
+
+# Switch to the peak-shaving operation model
+fmodel.set_operation_model("peak-shaving")
+fmodel.run()
+powers_peak_shaving = fmodel.get_turbine_powers()/1000
+thrust_coefficients_peak_shaving = fmodel.get_turbine_thrust_coefficients()
+
+# Compare the power and thrust coefficients of the upstream turbine
+fig, ax = plt.subplots(2,1,sharex=True)
+ax[0].plot(
+    wind_speeds,
+    thrust_coefficients_base[:,0],
+    label="Without peak shaving",
+    color="black"
+)
+ax[0].plot(
+    wind_speeds,
+    thrust_coefficients_peak_shaving[:,0],
+    label="With peak shaving",
+    color="C0"
+)
+ax[1].plot(
+    wind_speeds,
+    powers_base[:,0],
+    label="Without peak shaving",
+    color="black"
+)
+ax[1].plot(
+    wind_speeds,
+    powers_peak_shaving[:,0],
+    label="With peak shaving",
+    color="C0"
+)
+ax[1].grid()
+ax[0].grid()
+ax[0].legend()
+ax[0].set_ylabel("Thrust coefficient [-]")
+ax[1].set_xlabel("Wind speed [m/s]")
+ax[1].set_ylabel("Power [kW]")
+
+# Look at the total power across the two turbines for each case
+fig, ax = plt.subplots(2,1,sharex=True,sharey=True)
+ax[0].fill_between(
+    wind_speeds,
+    0,
+    powers_base[:, 0]/1e6,
+    color='C0',
+    label='Turbine 1'
+)
+ax[0].fill_between(
+    wind_speeds,
+    powers_base[:, 0]/1e6,
+    powers_base[:, :2].sum(axis=1)/1e6,
+    color='C1',
+    label='Turbine 2'
+    )
+ax[0].plot(
+    wind_speeds,
+    powers_base[:,:2].sum(axis=1)/1e6,
+    color='k',
+    label='Farm'
+)
+ax[1].fill_between(
+    wind_speeds,
+    0,
+    powers_peak_shaving[:, 0]/1e6,
+    color='C0',
+    label='Turbine 1'
+)
+ax[1].fill_between(
+    wind_speeds,
+    powers_peak_shaving[:, 0]/1e6,
+    powers_peak_shaving[:, :2].sum(axis=1)/1e6,
+    color='C1',
+    label='Turbine 2'
+    )
+ax[1].plot(
+    wind_speeds,
+    powers_peak_shaving[:,:2].sum(axis=1)/1e6,
+    color='k',
+    label='Farm'
+)
+ax[0].legend()
+ax[0].set_title("Without peak shaving")
+ax[1].set_title("With peak shaving")
+ax[0].set_ylabel("Power [MW]")
+ax[1].set_ylabel("Power [MW]")
+ax[0].grid()
+ax[1].grid()
+
+ax[1].set_xlabel("Free stream wind speed [m/s]")
+
+plt.show()

examples/examples_heterogeneous/001_heterogeneous_inflow_single.py

@@ -17,7 +17,7 @@
 import numpy as np
 
 from floris import FlorisModel, TimeSeries
-from floris.flow_visualization import visualize_cut_plane
+from floris.flow_visualization import visualize_heterogeneous_cut_plane
 from floris.layout_visualization import plot_turbine_labels
 
 
@@ -65,15 +65,20 @@
     x_resolution=200, y_resolution=100, height=90.0
 )
 
-# Plot the horizontal plane
+# Plot the horizontal plane using the visualize_heterogeneous_cut_plane.
+# Note that this function is not very different than the standard
+# visualize_cut_plane except that it accepts the fmodel object in order to
+# visualize the boundary of the heterogeneous inflow region.
 fig, ax = plt.subplots()
-visualize_cut_plane(
+visualize_heterogeneous_cut_plane(
     horizontal_plane,
+    fmodel=fmodel,
     ax=ax,
     title="Horizontal plane at hub height",
     color_bar=True,
     label_contours=True,
 )
 plot_turbine_labels(fmodel, ax)
+ax.legend()
 
 plt.show()