Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Add FLORIS based wind-speed estimator to FLASC #228

Merged
merged 7 commits into from
Nov 26, 2024
Merged
Show file tree
Hide file tree
Changes from all commits
Commits
File filter

Filter by extension

Filter by extension

Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
Original file line number Diff line number Diff line change
@@ -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
}
92 changes: 91 additions & 1 deletion flasc/utilities/floris_tools.py
Original file line number Diff line number Diff line change
@@ -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

Expand Down Expand Up @@ -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.
Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

How strictly do we enforce that df_scada's power should be specified in kW? Presumably, if the power is not in kW, the approach will fail?

Copy link
Collaborator Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

That's true, I think we only recommend it. I guess it will fail very badly at least, but we could think about adding some tests that power is in kW. Any values in excess of 20,000 would be a clue?

Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Yeah, I'm trying to think whether any other procedures in FLASC assume that power is in kW. Possibly not, since we are often dividing a power by a reference power, so the units cancel out. But maybe there are other places on floris_tools.py where we make this assumption.

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)
Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Maybe similar to the comment above about ws3, it seems strange to define a gain value of -1 for the ws3 entry and then clip any interpolated value back to 0? Would that be equivalent to setting ws3 closer to ws2 and then setting the last value of ws_est_gain_y to 0?

Copy link
Collaborator Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Yes, I think in our case this would be equivalent

Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I've now added some comments describing the wind speeds ws0, ws1, ws2, ws3.


# 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(
Expand Down
41 changes: 41 additions & 0 deletions tests/floris_tools_test.py
Original file line number Diff line number Diff line change
Expand Up @@ -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,
)
Expand Down Expand Up @@ -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])
Loading