Skip to content

Sg/floris example fix #23

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.

Sign up for GitHub

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

Merged
skygering merged 2 commits into from
Aug 14, 2026
Merged

Sg/floris example fix #23

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
2 commits
Select commit Hold shift + click to select a range
ca79a1e
Clean up example 8
skygering Aug 14, 2026
f22aa60
Update lock file
skygering Aug 14, 2026
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
2 changes: 1 addition & 1 deletion MITRotor/ReferenceTurbines/ROSCO_IEA15MW.yaml
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -31,7 +31,7 @@ controller_params:
F_LPFType: 2 # {1: first-order low-pass filter, 2: second-order low-pass filter}, [rad/s] (currently filters generator speed and pitch control signals)
F_NotchType: 0 # Notch on the measured generator speed {0: disable, 1: enable}
IPC_ControlMode: 0 # Turn Individual Pitch Control (IPC) for fatigue load reductions (pitch contribution) {0: off, 1: 1P reductions, 2: 1P+2P reductions}
VS_ControlMode: 1 # Generator torque control mode in above rated conditions {0: constant torque, 1: constant power, 2: TSR tracking PI control}
VS_ControlMode: 3 # Generator torque control mode in above rated conditions {0: constant torque, 1: constant power, 2: TSR tracking PI control}
PC_ControlMode: 1 # Blade pitch control mode {0: No pitch, fix to fine pitch, 1: active PI blade pitch control}
Y_ControlMode: 0 # Yaw control mode {0: no yaw control, 1: yaw rate control, 2: yaw-by-IPC}
SS_Mode: 1 # Setpoint Smoother mode {0: no setpoint smoothing, 1: introduce setpoint smoothing}
Expand Down
110 changes: 26 additions & 84 deletions examples/example_08_floris_rosco.py
View file
Open in desktop
Original file line number Diff line number Diff line change
@@ -1,12 +1,8 @@
# Python modules
import os
import matplotlib.pyplot as plt
import tempfile
from pathlib import Path
from ruamel.yaml import YAML
import numpy as np
import pandas as pd
import polars as pl
from scipy.interpolate import interp1d
import time

Expand All @@ -22,28 +18,6 @@

figdir = Path("fig")

def change_control_param(
param_key, param_value, template_yaml, bem,
regenerate = False, save_control_file = "control.csv", yaw_grid_deg = np.linspace(-25.0,25.0,50),
):
yaml = YAML()
with tempfile.TemporaryDirectory() as tmpdir:
temp_yaml = Path(tmpdir) / "rosco_temp.yaml"

with open(template_yaml) as f:
data = yaml.load(f)

data["controller_params"][param_key] = param_value

with open(temp_yaml, "w") as f:
yaml.dump(data, f)

pitch_interp, tsr_interp, rated_rotor_speed = get_rosco_control_interps(
temp_yaml, bem,
regenerate = regenerate, save_control_file = save_control_file, yaw_grid_deg = yaw_grid_deg
)
return pitch_interp, tsr_interp, rated_rotor_speed

def get_turbine_power_coefficent(fturbine, fmodel, wind_speeds):
rotor_area = np.pi * fturbine.bem_model.rotor.R**2
floris_power = np.squeeze(fmodel.get_turbine_powers())
Expand Down Expand Up @@ -78,7 +52,7 @@ def main():
cache_file = cache_dir / "lut.csv"
lut_model = UnifiedMomentumLUT(
cache_fn=cache_file,
regenerate=False,
regenerate=True,
LUT_Cts=np.linspace(-0.5,1.5,40),
LUT_yaws=yaws_degs, # Note here that we are only allowing for yaws <25 degrees! Update if needed!
)
Expand All @@ -87,38 +61,33 @@ def main():

# Load IEA15MW ROSCO parameters to make controllers
start = time.time()
rosco_yaml_VS_1 = "MITRotor/ReferenceTurbines/ROSCO_IEA15MW.yaml"
rosco_yaml_VS_3 = "MITRotor/ReferenceTurbines/ROSCO_IEA15MW.yaml"

# takes ~25min with 8 workers
rosco_VS_1_pitch_interp, rosco_VS_1_tsr_interp, rosco_VS_1_rated_rotorspeed = get_rosco_control_interps(
rosco_yaml_VS_1, bem,
regenerate = False, save_control_file = "control_vs_1.csv",
yaw_grid_deg=yaws_degs
rosco_VS_3_pitch_interp, rosco_VS_3_tsr_interp, rosco_VS_3_rated_rotorspeed = get_rosco_control_interps(
rosco_yaml_VS_3, bem,
regenerate = True, save_control_file = "control_VS_3.csv",
yaw_grid_deg=yaws_degs,
n_jobs = 8 # NOTE: making more than one control scheme with n_jobs > 1 within the same file tends to cause crashes from within ROSCO
# due to memory issues. Make each control scheme in its own python file to avoid.
)

end = time.time()
print(f"Time to make control CSV: {end - start}")
# takes ~25min with 8 workers
rosco_VS_3_pitch_interp, rosco_VS_3_tsr_interp, rosco_VS_3_rated_rotorspeed = change_control_param(
"VS_ControlMode", 3, rosco_yaml_VS_1, bem,
regenerate = True, save_control_file = "control_vs_3.csv",
yaw_grid_deg=yaws_degs
)

# Plot IEA15MW control from ROSCO paper, ROSOC control with VS_Control_Mode = 1, and ROSOC control withVS_Control_Mode = 3
# Plot IEA15MW control from ROSCO paper, ROSOC control with VS_Control_Mode = 3, and ROSOC control withVS_Control_Mode = 3
pitch_abbas = abbas_pitch_interp(wind_speeds)
pitch_ps0 = np.rad2deg(query_controls(rosco_VS_3_pitch_interp, wind_speeds, 0.0))
pitch_ps3 = np.rad2deg(query_controls(rosco_VS_1_pitch_interp, wind_speeds, 0.0))
pitch_vs1 = np.rad2deg(query_controls(rosco_VS_3_pitch_interp, wind_speeds, 0.0))

tsr_abbas = abbas_tsr_interp(wind_speeds)
tsr_ps0 = query_controls(rosco_VS_3_tsr_interp, wind_speeds, 0.0)
tsr_ps3 = query_controls(rosco_VS_1_tsr_interp, wind_speeds, 0.0)
tsr_vs1 = query_controls(rosco_VS_3_tsr_interp, wind_speeds, 0.0)

# Create figure
fig, ax = plt.subplots(1, 2, figsize = (10,4), sharey = True, constrained_layout=True)
fig.suptitle(fr"Setpoint Trajectories")
# Plot pitch
ax[0].plot(wind_speeds, pitch_ps0, label="ROSCO:VS_Control_Mode = 3", lw=3)
ax[0].plot(wind_speeds, pitch_ps3, label="ROSCO: VS_Control_Mode = 1", lw=3, linestyle = "dashed")
ax[0].plot(wind_speeds, pitch_vs1, label="ROSCO: VS_Control_Mode = 3", lw=3, linestyle = "dashed")
ax[0].plot(wind_speeds, pitch_abbas, label="Abbas et al. Fig 2", lw=3, linestyle = "dotted")

ax[0].set_xlabel("Wind Speed [m/s]")
Expand All @@ -127,8 +96,7 @@ def main():
ax[0].grid(True)

# Plot TSR
ax[1].plot(wind_speeds, tsr_ps0, label="MITRotor+FLORIS+ROSCO:VS_Control_Mode = 3", lw=3)
ax[1].plot(wind_speeds, tsr_ps3, label="MITRotor+FLORIS+ROSCO: VS_Control_Mode = 1", lw=3, linestyle = "dashed")
ax[1].plot(wind_speeds, tsr_vs1, label="MITRotor+FLORIS+ROSCO: VS_Control_Mode = 3", lw=3, linestyle = "dashed")
ax[1].plot(wind_speeds, tsr_abbas, label="Abbas et al. ROSCO Control", lw=3, linestyle = "dotted")

ax[1].set_xlabel("Wind Speed [m/s]")
Expand Down Expand Up @@ -158,14 +126,6 @@ def main():
tsr_interp = rosco_VS_3_tsr_interp,
rated_rotor_speed = rosco_VS_3_rated_rotorspeed,
)

floris_VS_1_turbine = MITRotorTurbine(
bem_model = bem,
pitch_interp = rosco_VS_1_pitch_interp,
pitch_rad = True,
tsr_interp = rosco_VS_1_tsr_interp,
rated_rotor_speed = rosco_VS_1_rated_rotorspeed,
)
yaw = 0.0 # degrees
yaw_angles = [[yaw] for _ in np.arange(len(wind_speeds))]
fmodel_default = FlorisModel("defaults")
Expand All @@ -182,20 +142,14 @@ def main():
Ct_abbas = fmodel_abbas.get_turbine_thrust_coefficients()
Cp_abbas, P_abbas = get_turbine_power_coefficent(floris_abbas_turbine, fmodel_abbas, wind_speeds)


fmodel_VS_3 = FlorisModel("defaults")
fmodel_VS_3.set(layout_x = [0.0], layout_y = [0.0], wind_data = time_series, yaw_angles = yaw_angles)
fmodel_VS_3.set_operation_model(floris_VS_3_turbine)
fmodel_VS_3.run()
fmodel_VS_3.run()
Ct_VS_3 = fmodel_VS_3.get_turbine_thrust_coefficients()
Cp_VS_3, P_VS_3 = get_turbine_power_coefficent(floris_VS_3_turbine, fmodel_VS_3, wind_speeds)

fmodel_VS_1 = FlorisModel("defaults")
fmodel_VS_1.set(layout_x = [0.0], layout_y = [0.0], wind_data = time_series, yaw_angles = yaw_angles)
fmodel_VS_1.set_operation_model(floris_VS_1_turbine)
fmodel_VS_1.run()
Ct_VS_1 = fmodel_VS_1.get_turbine_thrust_coefficients()
Cp_VS_1, P_VS_1 = get_turbine_power_coefficent(floris_VS_1_turbine, fmodel_VS_1, wind_speeds)

# plot CT values against one another and against IEA15MW from figure 3.1-C (https://docs.nlr.gov/docs/fy20osti/75698.pdf)
fig, (ax1, ax2, ax3) = plt.subplots(ncols = 3, sharey = False, figsize = (16,4), constrained_layout=True)
fig.suptitle(fr"$C_T$ and $C_P$ with Yaw = {yaw} for Different Control Strategies")
Expand All @@ -207,11 +161,7 @@ def main():
wind_speeds, Ct_abbas, label="FLORIS+MITROTOR Abbas et al. Control",
lw=3, linestyle = "solid", zorder = 1,
)
ax1.plot(
wind_speeds, Ct_VS_3, label="FLORIS+MITROTOR+ROSCO:VS_Control_Mode = 3",
lw=3, linestyle = "dashed", zorder = 1,
)
ax1.plot(wind_speeds, Ct_VS_1, label="FLORIS+MITROTOR+ROSCO: VS_Control_Mode = 1",
ax1.plot(wind_speeds, Ct_VS_3, label="FLORIS+MITROTOR+ROSCO: VS_Control_Mode = 3",
lw=3, linestyle = "dotted", zorder = 1,
)
ax1.set_xlabel("Wind Speed [m/s]")
Expand All @@ -227,11 +177,7 @@ def main():
wind_speeds, Cp_abbas, label="Abbas et al. ROSCO Control",
lw=3, linestyle = "solid", zorder = 1,
)
ax2.plot(
wind_speeds, Cp_VS_3, label="FLORIS + MITROTOR + ROSCO:VS_Control_Mode = 3",
lw=3, linestyle = "dashed", zorder = 1,
)
ax2.plot(wind_speeds, Cp_VS_1, label="FLORIS + MITROTOR + ROSCO: VS_Control_Mode = 1",
ax2.plot(wind_speeds, Cp_VS_3, label="FLORIS + MITROTOR + ROSCO: VS_Control_Mode = 3",
lw=3, linestyle = "dotted", zorder = 1,
)
ax2.set_xlabel("Wind Speed [m/s]")
Expand All @@ -247,11 +193,7 @@ def main():
wind_speeds, P_abbas, label="Abbas et al. ROSCO Control",
lw=3, linestyle = "solid", zorder = 1,
)
ax3.plot(
wind_speeds, P_VS_3, label="FLORIS + MITROTOR + ROSCO:VS_Control_Mode = 3",
lw=3, linestyle = "dashed", zorder = 1,
)
ax3.plot(wind_speeds, P_VS_1, label="FLORIS + MITROTOR + ROSCO: VS_Control_Mode = 1",
ax3.plot(wind_speeds, P_VS_3, label="FLORIS + MITROTOR + ROSCO: VS_Control_Mode = 3",
lw=3, linestyle = "dotted", zorder = 1,
)
ax3.axhline(y = 15, lw=3, linestyle = "dotted", zorder = 0, color = "k", label = "15MW")
Expand All @@ -266,15 +208,15 @@ def main():
# plot the difference in CT/CP curves for different yaw values
yaw_list = [0.0, 10.0, 20.0] # degrees
fig, (ax1, ax2, ax3) = plt.subplots(ncols = 3, sharey = False, figsize = (10,4), constrained_layout=True)
fig.suptitle(fr"$C_T$ and $C_P$ under Yaw for FLORIS + MITROTOR + ROSCO: VS_Control_Mode = 1")
fig.suptitle(fr"$C_T$ and $C_P$ under Yaw for FLORIS + MITROTOR + ROSCO: VS_Control_Mode = 3")
for (i, yaw) in enumerate(yaw_list):
yaw_angles = [[yaw] for _ in np.arange(len(wind_speeds))]
fmodel_VS_1 = FlorisModel("defaults")
fmodel_VS_1.set(layout_x = [0.0], layout_y = [0.0], wind_data = time_series, yaw_angles = yaw_angles)
fmodel_VS_1.set_operation_model(floris_VS_1_turbine)
fmodel_VS_1.run()
ct = fmodel_VS_1.get_turbine_thrust_coefficients()
cp, p = get_turbine_power_coefficent(floris_VS_1_turbine, fmodel_VS_1, wind_speeds)
fmodel_VS_3 = FlorisModel("defaults")
fmodel_VS_3.set(layout_x = [0.0], layout_y = [0.0], wind_data = time_series, yaw_angles = yaw_angles)
fmodel_VS_3.set_operation_model(floris_VS_3_turbine)
fmodel_VS_3.run()
ct = fmodel_VS_3.get_turbine_thrust_coefficients()
cp, p = get_turbine_power_coefficent(floris_VS_3_turbine, fmodel_VS_3, wind_speeds)

ax1.plot(wind_speeds, ct, label=fr"Yaw = ${yaw}^\circ$",
lw=3, linestyle = "solid", zorder = 1,
Expand Down
Loading
Loading