-
Notifications
You must be signed in to change notification settings - Fork 24
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
Changes from all commits
15fd0d0
69d0a7c
def5c82
bd6828d
f08063b
36e8fca
7396fe1
File filter
Filter by extension
Conversations
Jump to
Diff view
Diff view
There are no files selected for viewing
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 | ||
} |
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 | ||
|
||
|
@@ -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) | ||
There was a problem hiding this comment. Choose a reason for hiding this commentThe reason will be displayed to describe this comment to others. Learn more. Maybe similar to the comment above about There was a problem hiding this comment. Choose a reason for hiding this commentThe reason will be displayed to describe this comment to others. Learn more. Yes, I think in our case this would be equivalent There was a problem hiding this comment. Choose a reason for hiding this commentThe reason will be displayed to describe this comment to others. Learn more. I've now added some comments describing the wind speeds |
||
|
||
# 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( | ||
|
There was a problem hiding this comment.
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?There was a problem hiding this comment.
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?
There was a problem hiding this comment.
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.