Add FLORIS based wind-speed estimator to FLASC

examples_artificial_data/04_floris_tuning/01_wind_speed_estimator.ipynb

@@ -0,0 +1,180 @@
+{
+ "cells": [
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "# Wind speed estimation\n",
+    "\n",
+    "This is small example demonstrates the usage of `estimate_ws_with_floris` function to estimate the wind speed at a given point using the FLORIS model. "
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "In this example the wind speed is estimated from artificial power data generated using the FLORIS model.  The data includes a 1 m/s bias to the original wind speed which is corrected via the estimator"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "from pathlib import Path\n",
+    "\n",
+    "import matplotlib.pyplot as plt\n",
+    "import numpy as np\n",
+    "import seaborn as sns\n",
+    "from floris import FlorisModel, TimeSeries\n",
+    "\n",
+    "from flasc import FlascDataFrame\n",
+    "from flasc.utilities.floris_tools import estimate_ws_with_floris"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "### Establish biased wind speed data"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "file_path = Path.cwd()\n",
+    "fm_path = file_path / \"../floris_input_artificial/gch.yaml\"\n",
+    "fm = FlorisModel(fm_path)\n",
+    "\n",
+    "# Set to 1 turbine layout with a linear sweep over wind speeds\n",
+    "N = 25\n",
+    "wind_speeds = np.linspace(0.01, 20.0, N)\n",
+    "time_series = TimeSeries(\n",
+    "    wind_speeds=wind_speeds, wind_directions=270.0, turbulence_intensities=0.06\n",
+    ")\n",
+    "fm.set(layout_x=[0], layout_y=[0], wind_data=time_series)\n",
+    "fm.run()\n",
+    "\n",
+    "# Construct df_scada from the FLORIS output\n",
+    "df_scada = FlascDataFrame(\n",
+    "    {\n",
+    "        \"time\": np.arange(0, N),\n",
+    "        \"pow_000\": fm.get_turbine_powers().squeeze() / 1000.0,\n",
+    "        \"ws_000\": wind_speeds + 1.0,  # Add 1m/s bias\n",
+    "    }\n",
+    ")\n",
+    "print(df_scada.head())"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "### Run wind speed estimation procedure"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "df_scada = estimate_ws_with_floris(df_scada, fm)\n",
+    "print(df_scada.head())"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "### Calculate power with wind speed and estimated wind speed"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "# Calculate with biased measurement\n",
+    "time_series = TimeSeries(\n",
+    "    wind_speeds=df_scada.ws_000.values, wind_directions=270.0, turbulence_intensities=0.06\n",
+    ")\n",
+    "fm.set(wind_data=time_series)\n",
+    "fm.run()\n",
+    "power_from_original_ws = fm.get_turbine_powers().squeeze() / 1000.0\n",
+    "\n",
+    "# Calculate with estimated wind speed\n",
+    "time_series = TimeSeries(\n",
+    "    wind_speeds=df_scada.ws_est_000.values, wind_directions=270.0, turbulence_intensities=0.06\n",
+    ")\n",
+    "fm.set(wind_data=time_series)\n",
+    "fm.run()\n",
+    "power_from_estimated_ws = fm.get_turbine_powers().squeeze() / 1000.0\n",
+    "\n",
+    "# Compute the error of each relative to measured power\n",
+    "original_ws_error = df_scada.pow_000.values - power_from_original_ws\n",
+    "estimated_ws_error = df_scada.pow_000.values - power_from_estimated_ws\n",
+    "\n",
+    "# Plot the error against the measured power\n",
+    "fig, ax = plt.subplots()\n",
+    "sns.scatterplot(x=df_scada.pow_000, y=original_ws_error, ax=ax, label=\"Original WS\", s=100)\n",
+    "sns.scatterplot(x=df_scada.pow_000, y=estimated_ws_error, ax=ax, label=\"Estimated WS\")\n",
+    "ax.set_xlabel(\"Measured Power [kW]\")\n",
+    "ax.set_ylabel(\"Error [kW]\")\n",
+    "ax.set_title(\"Error vs Measured Power\")\n",
+    "ax.grid(True)\n",
+    "plt.show()\n",
+    "\n",
+    "# Plot estimated vs real wind speed\n",
+    "fig = plt.figure()\n",
+    "gs = fig.add_gridspec(2, 1, height_ratios=[2, 1])\n",
+    "\n",
+    "ax = fig.add_subplot(gs[0])\n",
+    "ax.plot(wind_speeds, wind_speeds, color=\"black\", linestyle=\"--\", label=\"Truth\")\n",
+    "sns.scatterplot(x=wind_speeds, y=df_scada.ws_000, ax=ax, label=\"Original WS\", s=100)\n",
+    "sns.scatterplot(x=wind_speeds, y=df_scada.ws_est_000, ax=ax, label=\"Estimated WS\")\n",
+    "ax.grid(True)\n",
+    "ax.set_ylabel(\"Meas.and Est. Wind Speed [m/s]\")\n",
+    "\n",
+    "ax = fig.add_subplot(gs[1])\n",
+    "ax.plot(wind_speeds, df_scada.ws_est_gain_000, color=\"black\")\n",
+    "ax.set_xlabel(\"True Wind Speed [m/s]\")\n",
+    "ax.set_ylabel(\"Est. gain [-]\")\n",
+    "ax.grid(True)"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": []
+  }
+ ],
+ "metadata": {
+  "kernelspec": {
+   "display_name": ".venv",
+   "language": "python",
+   "name": "python3"
+  },
+  "language_info": {
+   "codemirror_mode": {
+    "name": "ipython",
+    "version": 3
+   },
+   "file_extension": ".py",
+   "mimetype": "text/x-python",
+   "name": "python",
+   "nbconvert_exporter": "python",
+   "pygments_lexer": "ipython3",
+   "version": "3.10.4"
+  }
+ },
+ "nbformat": 4,
+ "nbformat_minor": 2
+}

flasc/utilities/floris_tools.py

@@ -1,16 +1,20 @@
 """Utility functions that use FlorisModels."""
 
+from __future__ import annotations
+
 import copy
 from time import perf_counter as timerpc
+from typing import Union
 
 import matplotlib.pyplot as plt
 import numpy as np
 import pandas as pd
-from floris import TimeSeries, WindTIRose
+from floris import FlorisModel, TimeSeries, WindTIRose
 from floris.utilities import wrap_360
 from scipy import interpolate
 from scipy.stats import norm
 
+from flasc import FlascDataFrame
 from flasc.logging_manager import LoggingManager
 from flasc.utilities import utilities as fsut
 
@@ -503,6 +507,92 @@ def add_gaussian_blending_to_floris_approx_table(df_fi_approx, wd_std=3.0, pdf_c
     return df_fi_approx_gauss
 
 
+def estimate_ws_with_floris(
+    df_scada: Union[pd.DataFrame, FlascDataFrame],
+    fm: FlorisModel,
+    verbose: bool = False,
+) -> Union[pd.DataFrame, FlascDataFrame]:
+    """Estimate the wind speed at the turbine locations using a FLORIS model.
+
+    This function estimates the wind speed at the turbine locations using a FLORIS model.
+    The approach follows the example from the RES wind-up method `add_ws_est_one_ttype`
+    (https://github.com/resgroup/wind-up/blob/main/wind_up/ws_est.py) by Alex Clerc.
+    However, in this implementation, FLORIS provides the power curves directly rather
+    than their being learned from data.  In this way, the estimated wind speed is
+    the speed which would cause FLORIS to predict a power matched to the SCADA data.
+    This will be especially useful in model fitting to avoid any error on un-waked turbines.
+
+    Args:
+        df_scada (pd.DataFrame | FlascDataFrame): Pandas DataFrame with the SCADA data.
+        fm (FlorisModel): FLORIS model object.
+        verbose (bool, optional): Print warnings and information to the console. Defaults to False.
+
+    Returns:
+        pd.DataFrame | FlascDataFrame: Pandas DataFrame with the estimated wind speed
+            columns added.
+    """
+    # To allow that turbine types might not be the same, loop over all turbines
+    for ti in range(fm.n_turbines):
+        if verbose:
+            logger.info(f"Estimating wind speed for turbine {ti} of {fm.n_turbines}.")
+
+        # Get the power curve for this turbine from FLORIS
+        ws_pc = fm.core.farm.turbine_definitions[ti]["power_thrust_table"]["wind_speed"]
+        pow_pc = fm.core.farm.turbine_definitions[ti]["power_thrust_table"]["power"]
+
+        # Take the diff of the power curve
+        pow_pc_diff = np.diff(pow_pc)
+
+        # If first entry is not positive, remove it from ws_pc and pow_pc, loop until not true
+        while pow_pc_diff[0] <= 0:
+            ws_pc = ws_pc[1:]
+            pow_pc = pow_pc[1:]
+            pow_pc_diff = np.diff(pow_pc)
+
+        # Find the first point where the diff is not positive and remove it and all following points
+        if np.any(pow_pc_diff <= 0):
+            idx = np.where(pow_pc_diff <= 0)[0][0]
+            ws_pc = ws_pc[: idx + 1]
+            pow_pc = pow_pc[: idx + 1]
+            pow_pc_diff = np.diff(pow_pc)
+
+        # Identify certain quantities for establishing the estimator and blend schedule
+        rated_power = np.max(pow_pc)
+
+        # Following the gain scheduling approach of `add_ws_est_one_ttype`
+        # define the gain via four wind speeds, however, more deterministically since
+        # driven by FLORIS
+        # ws0 (zero wind speed) is the wind speed at which partial power dependence begins
+        # ws1 (10% of rated) is the wind speed at which full power dependence begins
+        # ws2 (90% of rated) is the wind speed at which full power dependence ends
+        # ws3 (90% of rated + 0.5 m/s) is the wind speed at which partial power dependence ends
+        ws0 = 0
+        ws1 = float(np.interp(rated_power * 0.1, pow_pc, ws_pc))
+        ws2 = float(np.interp(rated_power * 0.9, pow_pc, ws_pc))
+        ws3 = ws2 + 0.5
+
+        # For starters make simple gain schedule
+        ws_est_gain_x = [ws0, ws1, ws2, ws3]
+        ws_est_gain_y = [0, 1, 1, 0]
+
+        ws_est_gain = np.interp(df_scada["ws_{:03d}".format(ti)], ws_est_gain_x, ws_est_gain_y)
+        ws_est_gain = np.clip(ws_est_gain, 0, 1)
+
+        # Now estimate wind speed from power
+        ws_est_from_pow = np.interp(df_scada["pow_{:03d}".format(ti)], pow_pc, ws_pc)
+        ws_est = (
+            ws_est_gain * ws_est_from_pow + (1.0 - ws_est_gain) * df_scada["ws_{:03d}".format(ti)]
+        )
+
+        # Add the estimated wind speed to the dataframe
+        df_scada["ws_est_{:03d}".format(ti)] = ws_est
+
+        # Store the gain as well, at least for now
+        df_scada["ws_est_gain_{:03d}".format(ti)] = ws_est_gain
+
+    return df_scada
+
+
 # TODO Is this function in the right module?
 # TODO Should include itself have a default?
 def get_turbs_in_radius(

tests/floris_tools_test.py

@@ -7,6 +7,7 @@
 from flasc.utilities.floris_tools import (
     add_gaussian_blending_to_floris_approx_table,
     calc_floris_approx_table,
+    estimate_ws_with_floris,
     get_dependent_turbines_by_wd,
     interpolate_floris_from_df_approx,
 )
@@ -151,3 +152,43 @@ def test_get_dependent_turbines_by_wd(self):
         # Test the limit_number
         dep = get_dependent_turbines_by_wd(fi, 2, np.array([226]), limit_number=1)
         self.assertEqual(dep[0], [1])
+
+    def test_estimate_ws_with_floris(self):
+        # Load FLORIS object
+        fm, _ = load_floris()
+
+        # Set as two turbine layout
+        fm.set(layout_x=[0.0, 0.0], layout_y=[0.0, 500.0])
+
+        # Create a sample SCADA dataframe
+        df_scada = pd.DataFrame(
+            {
+                "ws_000": [0.5, 8.5, 20.0],
+                "ws_001": [8.0, 8.5, 9.0],
+                "pow_000": [100.0, 100.0, 100.0],
+                "pow_001": [1000.0, 1000.0, 1000.0],
+            }
+        )
+
+        # Estimate wind speed using FLORIS
+        df_estimated = estimate_ws_with_floris(df_scada, fm)
+
+        # Check if the estimated wind speed columns are added
+        self.assertTrue("ws_est_000" in df_estimated.columns)
+        self.assertTrue("ws_est_001" in df_estimated.columns)
+
+        # Check if the estimated wind speed gain columns are added
+        self.assertTrue("ws_est_gain_000" in df_estimated.columns)
+        self.assertTrue("ws_est_gain_001" in df_estimated.columns)
+
+        # Check that the third element of ws_est_000 are
+        # unchanged from ws_000
+        self.assertTrue(
+            np.all(df_estimated["ws_est_000"].values[[2]] == df_scada["ws_000"].values[[2]])
+        )
+
+        # Check the the middle element of ws_est_000 is less than from ws_000
+        self.assertTrue(df_estimated["ws_est_000"].values[1] < df_scada["ws_000"].values[1])
+
+        # Check that estimated middle value for turbine 1 is greater that that of turbine 0
+        self.assertTrue(df_estimated["ws_est_001"].values[1] > df_estimated["ws_est_000"].values[1])