diff --git a/docs/wind.md b/docs/wind.md index 95e896f8..e83034c4 100644 --- a/docs/wind.md +++ b/docs/wind.md @@ -1,16 +1,38 @@ # Wind Farm Components -## Wind_MesoToPower +Hercules provides four wind farm simulation components that differ in their approach to wake modeling and data sources. The first three components support both simple filter-based turbine models and 1-degree-of-freedom (1-DOF) turbine dynamics, while the fourth component uses SCADA power data directly. -Wind_MesoToPower is a comprehensive wind farm simulator that focuses on meso-scale phenomena by applying a separate wind speed time signal to each turbine model derived from data. It combines FLORIS wake modeling with detailed turbine dynamics for long-term wind farm performance analysis. +## Wind_MesoToPower (Dynamic Wake Model) + +Wind_MesoToPower is a comprehensive wind farm simulator that computes wake effects dynamically at each time step (or at intervals specified by `floris_update_time_s`). It focuses on meso-scale phenomena by applying a separate wind speed time signal to each turbine model derived from data. This model combines FLORIS wake modeling with detailed turbine dynamics for long-term wind farm performance analysis. + +**Use this model when:** +- Turbines have individual power setpoints or non-uniform operation +- Precise wake modeling is required for each control action +- Turbines may be partially derated or individually controlled + +## Wind_MesoToPowerPrecomFloris (Precomputed Wake Model) + +Wind_MesoToPowerPrecomFloris is an optimized variant that pre-computes all FLORIS wake deficits at initialization for improved simulation performance. This approach provides significant speed improvements while conservatively assuming wakes are always based on nominal operation. + +**Use this model when:** +- Not investigating wakes of derated turbines or wake losses can be conservatively estimated. + + +## Wind_MesoToPowerNoAddedWakes (No Wake Modeling) + +Wind_MesoToPowerNoAddedWakes assumes that wake effects are already included in the input wind data and performs no wake modeling during simulation. Model is appropriate for using SCADA data of operational farm since wake losses already included in data. + + +## WindFarmSCADAPower (SCADA Power Data) + +WindFarmSCADAPower uses SCADA power measurements directly rather than computing power from wind speeds and turbine models. This component applies a filter to the SCADA power data to simulate turbine response dynamics and respects power setpoint constraints. -## Wind_MesoToPowerPrecomFloris -Wind_MesoToPowerPrecomFloris is an optimized variant of Wind_MesoToPower that pre-computes FLORIS wake deficits for improved simulation performance. This approach trades some accuracy for significant speed improvements in specific operating scenarios. ## Overview -Both wind farm components integrate FLORIS for wake effects with individual turbine models to simulate wind farm behavior over extended periods. They support both simple filter-based turbine models and 1-degree-of-freedom (1-DOF) turbine dynamics. +The first three wind farm components apply wind speed time signals to turbine models to simulate wind farm behavior over extended periods. They differ only in how wake effects are computed and applied. The WindFarmSCADAPower component uses a fundamentally different approach by using actual SCADA power measurements as input. ### Precomputed FLORIS Approach @@ -45,19 +67,49 @@ Required parameters for Wind_MesoToPowerPrecomFloris: - `floris_update_time_s`: Determines the cadence of wake precomputation. At each cadence tick, the last `floris_update_time_s` seconds are averaged and used to evaluate FLORIS. The computed wake deficits are then applied until the next cadence tick. - `log_channels`: List of output channels to log. See [Logging Configuration](#logging-configuration) section below for details. +### Wind_MesoToPowerNoAddedWakes Specific Parameters + +Required parameters for Wind_MesoToPowerNoAddedWakes: +- `floris_update_time_s`: Required for interface consistency but not used (no FLORIS calculations performed) +- `floris_input_file`: Still required to read turbine power curve and properties +- `log_channels`: List of output channels to log. See [Logging Configuration](#logging-configuration) section below for details. + +### WindFarmSCADAPower Specific Parameters + +Required parameters for WindFarmSCADAPower: +- `scada_filename`: Path to SCADA data file (CSV, pickle, or feather format) +- `turbine_file_name`: Turbine model configuration (for filter parameters) +- `log_channels`: List of output channels to log. See [Logging Configuration](#logging-configuration) section below for details. + +**SCADA File Format:** + +The SCADA file must contain the following columns: +- `time_utc`: Timestamps in UTC (ISO 8601 format or parseable datetime strings) +- `wd_mean`: Mean wind direction in degrees +- `pow_###`: Power output for each turbine (e.g., `pow_000`, `pow_001`, `pow_002`) + +Optional columns: +- `ws_###`: Wind speed for each turbine (e.g., `ws_000`, `ws_001`, `ws_002`) +- `ws_mean`: Mean wind speed (used if individual turbine speeds not provided) +- `ti_###`: Turbulence intensity for each turbine (defaults to 0.08 if not provided) + +The number of turbines and rated power are automatically inferred from the SCADA data. + ## Turbine Models +**Note:** WindFarmSCADAPower uses only the filter model for power smoothing, as power values come directly from SCADA data rather than being computed from wind speeds. + ### Filter Model Simple first-order filter for power output smoothing with configurable time constants. ### 1-DOF Model -Advanced model with rotor dynamics, pitch control, and generator torque control. +Advanced model with rotor dynamics, pitch control, and generator torque control. Not applicable to WindFarmSCADAPower. ## Outputs ### Common Outputs -Both components provide these outputs in the h_dict at each simulation step: +All four components provide these outputs in the h_dict at each simulation step: - `power`: Total wind farm power (kW) - `turbine_powers`: Individual turbine power outputs (array, kW) - `turbine_power_setpoints`: Current power setpoint values (array, kW) @@ -67,6 +119,8 @@ Both components provide these outputs in the h_dict at each simulation step: - `wind_speeds_background`: Per-turbine background wind speeds (array, m/s) - `wind_speeds_withwakes`: Per-turbine with-wakes wind speeds (array, m/s) +**Note for Wind_MesoToPowerNoAddedWakes and WindFarmSCADAPower:** In these models (no wake modeling), `wind_speeds_withwakes` equals `wind_speeds_background` and `wind_speed_mean_withwakes` equals `wind_speed_mean_background`. + ## Logging Configuration The `log_channels` parameter controls which outputs are written to the HDF5 output file. This is a list of channel names. The `power` channel is always logged, even if not explicitly specified. diff --git a/examples/02c_wind_farm_realistic_inflow_direct/README.md b/examples/02c_wind_farm_realistic_inflow_direct/README.md new file mode 100644 index 00000000..a6116a74 --- /dev/null +++ b/examples/02c_wind_farm_realistic_inflow_direct/README.md @@ -0,0 +1,31 @@ +# Example 02c: Wind Farm Realistic Inflow (Direct - No Wake Modeling) + +## Description + +This example demonstrates the `Wind_MesoToPowerNoAddedWakes` wake model, which assumes that wake effects are already included in the input wind data and performs no additional wake modeling. + +The `Wind_MesoToPowerNoAddedWakes` component type uses `wake_model="no_added_wakes"` internally, which means: +- No FLORIS calculations are performed during the simulation (only at initialization to read turbine properties) +- `wind_speeds_withwakes` equals `wind_speeds_background` at all times +- Wake deficits are always zero +- Turbine dynamics (filter model or DOF1 model) still operate normally + +This example automatically generates the necessary input files in the centralized `examples/inputs/` folder when first run. + +## Running + +To run the example, execute the following command in the terminal: + +```bash +python hercules_runscript.py +``` + +## Outputs + +To plot the outputs run the following command in the terminal: + +```bash +python plot_outputs.py +``` + + diff --git a/examples/02c_wind_farm_realistic_inflow_direct/hercules_input.yaml b/examples/02c_wind_farm_realistic_inflow_direct/hercules_input.yaml new file mode 100644 index 00000000..d209b50c --- /dev/null +++ b/examples/02c_wind_farm_realistic_inflow_direct/hercules_input.yaml @@ -0,0 +1,42 @@ +# Input YAML for hercules + +# Name +name: example_02c + +### +# Describe this simulation setup +description: Wind Only Realistic Inflow (Direct - No Wake Modeling) + +dt: 1.0 +starttime_utc: "2024-06-24T16:59:08Z" # Jun 24, 2024 16:59:08 UTC (Zulu time) +endtime_utc: "2024-06-26T16:59:00Z" # ≈48 hours later (Jun 26, 2024 16:59:00 UTC) +verbose: False + +plant: + interconnect_limit: 45000 # kW + +wind_farm: + + component_type: Wind_MesoToPowerNoAddedWakes + floris_input_file: ../inputs/floris_input_large.yaml + wind_input_filename: ../inputs/wind_input_large.ftr + turbine_file_name: ../inputs/turbine_filter_model.yaml + log_file_name: outputs/log_wind_sim.log + log_channels: + - power + - wind_speed_mean_background + - wind_speed_mean_withwakes + - wind_direction_mean + - turbine_powers + - wind_speeds_withwakes + - wind_speeds_background + - turbine_power_setpoints + floris_update_time_s: 300.0 # Not used but kept for interface consistency + + +controller: + + + + + diff --git a/examples/02c_wind_farm_realistic_inflow_direct/hercules_runscript.py b/examples/02c_wind_farm_realistic_inflow_direct/hercules_runscript.py new file mode 100644 index 00000000..8e06f8a7 --- /dev/null +++ b/examples/02c_wind_farm_realistic_inflow_direct/hercules_runscript.py @@ -0,0 +1,53 @@ +import numpy as np +from hercules.hercules_model import HerculesModel +from hercules.utilities_examples import ensure_example_inputs_exist, prepare_output_directory + +prepare_output_directory() + +# Ensure example inputs exist +ensure_example_inputs_exist() + +# Initialize the Hercules model +hmodel = HerculesModel("hercules_input.yaml") + + +# Define a simple controller that sets all power setpoints to full rating +class ControllerFullRating: + """A simple controller that sets all turbines to full rating. + + This controller is appropriate for the direct wake model where + wake effects are already included in the input wind data. + """ + + def __init__(self, h_dict): + """Initialize the controller. + + Args: + h_dict (dict): The hercules input dictionary. + """ + pass + + def step(self, h_dict): + """Execute one control step. + + Args: + h_dict (dict): The hercules input dictionary. + + Returns: + dict: The updated hercules input dictionary. + """ + # Set all turbines to full rating + h_dict["wind_farm"]["turbine_power_setpoints"] = 5000 * np.ones( + h_dict["wind_farm"]["n_turbines"] + ) + + return h_dict + + +# Assign the controller to the Hercules model +hmodel.assign_controller(ControllerFullRating(hmodel.h_dict)) + +# Run the simulation +hmodel.run() + +hmodel.logger.info("Process completed successfully") diff --git a/examples/02c_wind_farm_realistic_inflow_direct/plot_outputs.py b/examples/02c_wind_farm_realistic_inflow_direct/plot_outputs.py new file mode 100644 index 00000000..c63bfcd6 --- /dev/null +++ b/examples/02c_wind_farm_realistic_inflow_direct/plot_outputs.py @@ -0,0 +1,79 @@ +# Plot the outputs of the simulation + +import matplotlib.pyplot as plt +from hercules import HerculesOutput + +# Read the Hercules output file using HerculesOutput +ho = HerculesOutput("outputs/hercules_output.h5") + +# Print metadata information +print("Simulation Metadata:") +ho.print_metadata() +print() + +# Create a shortcut to the dataframe +df = ho.df + +# Limit to the first 4 hours +df = df.iloc[: 3600 * 4] + +# Set number of turbines +turbines_to_plot = [0, 8] + +# Define a consistent color map with 9 +colors = [ + "tab:blue", + "tab:orange", + "tab:green", + "tab:red", + "tab:purple", + "tab:brown", + "tab:pink", + "tab:gray", + "tab:olive", +] + +fig, axarr = plt.subplots(2, 1, sharex=True) + +# Plot the wind speeds +ax = axarr[0] +for t_idx in turbines_to_plot: + ax.plot( + df["time"], + df[f"wind_farm.wind_speeds_background.{t_idx:03}"], + label=f"Wind Speed {t_idx}", + color=colors[t_idx], + ) + +# Note: In direct mode, wind_speeds_withwakes == wind_speeds_background + +# Plot the mean wind speed +ax.plot( + df["time"], + df["wind_farm.wind_speed_mean_background"], + label="Mean Wind Speed", + color="black", + lw=2, +) + +ax.grid(True) +ax.legend() +ax.set_ylabel("Wind Speed [m/s]") +ax.set_title("Direct Wake Model (No Wake Modeling)") + + +# Plot the power +ax = axarr[1] +for t_idx in turbines_to_plot: + ax.plot( + df["time"], + df[f"wind_farm.turbine_powers.{t_idx:03}"], + label=f"Turbine {t_idx}", + color=colors[t_idx], + ) + +ax.grid(True) +ax.legend() +ax.set_xlabel("Time [s]") +ax.set_ylabel("Power [kW]") +plt.show() diff --git a/hercules/hybrid_plant.py b/hercules/hybrid_plant.py index 28e66378..c5cca2a9 100644 --- a/hercules/hybrid_plant.py +++ b/hercules/hybrid_plant.py @@ -4,8 +4,8 @@ from hercules.plant_components.battery_simple import BatterySimple from hercules.plant_components.electrolyzer_plant import ElectrolyzerPlant from hercules.plant_components.solar_pysam_pvwatts import SolarPySAMPVWatts -from hercules.plant_components.wind_meso_to_power import Wind_MesoToPower -from hercules.plant_components.wind_meso_to_power_precom_floris import Wind_MesoToPowerPrecomFloris +from hercules.plant_components.wind_farm import WindFarm +from hercules.plant_components.wind_farm_scada_power import WindFarmSCADAPower from hercules.utilities import get_available_component_names, get_available_generator_names @@ -95,25 +95,39 @@ def get_plant_component(self, component_name, h_dict): Raises: Exception: If the component_type is not recognized. """ - if h_dict[component_name]["component_type"] == "Wind_MesoToPower": - return Wind_MesoToPower(h_dict) - - if h_dict[component_name]["component_type"] == "Wind_MesoToPowerPrecomFloris": - return Wind_MesoToPowerPrecomFloris(h_dict) - - if h_dict[component_name]["component_type"] == "SolarPySAMPVWatts": + component_type = h_dict[component_name]["component_type"] + + # Handle wind farm component types with unified WindFarm class + if component_type in [ + "Wind_MesoToPower", + "Wind_MesoToPowerPrecomFloris", + "Wind_MesoToPowerNoAddedWakes", + ]: + # Map component_type to wake_model + wake_model_map = { + "Wind_MesoToPower": "dynamic", + "Wind_MesoToPowerPrecomFloris": "precomputed", + "Wind_MesoToPowerNoAddedWakes": "no_added_wakes", + } + wake_model = wake_model_map[component_type] + return WindFarm(h_dict, wake_model=wake_model) + + if component_type == "WindFarmSCADAPower": + return WindFarmSCADAPower(h_dict) + + if component_type == "SolarPySAMPVWatts": return SolarPySAMPVWatts(h_dict) - if h_dict[component_name]["component_type"] == "BatteryLithiumIon": + if component_type == "BatteryLithiumIon": return BatteryLithiumIon(h_dict) - if h_dict[component_name]["component_type"] == "BatterySimple": + if component_type == "BatterySimple": return BatterySimple(h_dict) - if h_dict[component_name]["component_type"] == "ElectrolyzerPlant": + if component_type == "ElectrolyzerPlant": return ElectrolyzerPlant(h_dict) - raise Exception("Unknown component_type: ", h_dict[component_name]["component_type"]) + raise Exception("Unknown component_type: ", component_type) def step(self, h_dict): """Execute one simulation step for all plant components. diff --git a/hercules/plant_components/wind_meso_to_power.py b/hercules/plant_components/wind_farm.py similarity index 61% rename from hercules/plant_components/wind_meso_to_power.py rename to hercules/plant_components/wind_farm.py index 755373ed..264aef91 100644 --- a/hercules/plant_components/wind_meso_to_power.py +++ b/hercules/plant_components/wind_farm.py @@ -1,9 +1,10 @@ -# Implements the meso-scale wind model for Hercules. - +# Unified wind farm model for Hercules supporting multiple wake modeling strategies. import numpy as np import pandas as pd -from floris import FlorisModel +from floris import ApproxFlorisModel, FlorisModel +from floris.core import average_velocity +from floris.uncertain_floris_model import map_turbine_powers_uncertain from hercules.plant_components.component_base import ComponentBase from hercules.utilities import ( hercules_float_type, @@ -11,33 +12,69 @@ load_perffile, load_yaml, ) -from scipy.interpolate import interp1d from scipy.optimize import minimize_scalar from scipy.stats import circmean RPM2RADperSec = 2 * np.pi / 60.0 -class Wind_MesoToPower(ComponentBase): - def __init__(self, h_dict): - """Initialize the Wind_MesoToPower class. +class WindFarm(ComponentBase): + """Unified wind farm model with configurable wake modeling strategies. + + This model simulates wind farm performance by applying wind speed time signals + to turbine models. It supports three wake modeling strategies: + + 1. **dynamic**: Real-time FLORIS wake calculations at each time step or interval. + Use when turbines may have individual setpoints or non-uniform operation. - This model focuses on meso-scale wind phenomena by applying a separate wind speed - time signal to each turbine model derived from data. It combines FLORIS wake - modeling with detailed turbine dynamics for wind farm performance analysis. + 2. **precomputed**: Pre-computed FLORIS wake deficits for all conditions. + Use when all turbines operate uniformly (all on, all off, or uniform curtailment). + More efficient but less flexible than dynamic. + + 3. **none**: No wake modeling - wind speeds are used directly. + Use when wake effects are already included in the input data or when + wake modeling is not needed. + + All three strategies support detailed turbine dynamics (filter_model or dof1_model). + """ + + def __init__(self, h_dict, wake_model=None): + """Initialize the WindFarm class. Args: - h_dict (dict): Dictionary containing values for the simulation. + h_dict (dict): Dictionary containing simulation parameters. + wake_model (str, optional): Wake modeling strategy: "dynamic", "precomputed", + or "no_added_wakes". If None, infers from component_type for backward compatibility. + Defaults to None. + + Raises: + ValueError: If wake_model is invalid or required parameters are missing. """ # Store the name of this component self.component_name = "wind_farm" - # Store the type of this component - self.component_type = "Wind_MesoToPower" + # Infer wake_model from component_type if not provided (backward compatibility) + if wake_model is None: + wake_model = self._infer_wake_model_from_component_type(h_dict) + + # Validate wake_model + if wake_model not in ["dynamic", "precomputed", "no_added_wakes"]: + raise ValueError( + f"wake_model must be 'dynamic', 'precomputed', or " + f"'no_added_wakes', got '{wake_model}'" + ) + + self.wake_model = wake_model + + # Store the type of this component (for backward compatibility) + component_type = h_dict[self.component_name].get("component_type", "WindFarm") + self.component_type = component_type # Call the base class init super().__init__(h_dict, self.component_name) + self.logger.info(f"Initializing WindFarm with wake_model='{self.wake_model}'") + # Track the number of FLORIS calculations self.num_floris_calcs = 0 @@ -46,17 +83,16 @@ def __init__(self, h_dict): self.wind_input_filename = h_dict[self.component_name]["wind_input_filename"] self.turbine_file_name = h_dict[self.component_name]["turbine_file_name"] - # Require floris_update_time_s to be in the h_dict + # Require floris_update_time_s for interface consistency if "floris_update_time_s" not in h_dict[self.component_name]: raise ValueError("floris_update_time_s must be in the h_dict") - - # Save the floris update time and make sure it is at least 1 second self.floris_update_time_s = h_dict[self.component_name]["floris_update_time_s"] if self.floris_update_time_s < 1: raise ValueError("FLORIS update time must be at least 1 second") + self.logger.info("Reading in wind input file...") + # Read in the weather file data - # If a csv file is provided, read it in if self.wind_input_filename.endswith(".csv"): df_wi = pd.read_csv(self.wind_input_filename) elif self.wind_input_filename.endswith(".p") | self.wind_input_filename.endswith(".pkl"): @@ -68,6 +104,8 @@ def __init__(self, h_dict): else: raise ValueError("Wind input file must be a .csv or .p, .f or .ftr file") + self.logger.info("Checking wind input file...") + # Make sure the df_wi contains a column called "time_utc" if "time_utc" not in df_wi.columns: raise ValueError("Wind input file must contain a column called 'time_utc'") @@ -83,7 +121,7 @@ def __init__(self, h_dict): df_wi["time_utc"] = pd.to_datetime(df_wi["time_utc"], utc=True) # Ensure time_utc is timezone-aware (UTC) - if not pd.api.types.is_datetime64tz_dtype(df_wi["time_utc"]): + if not isinstance(df_wi["time_utc"].dtype, pd.DatetimeTZDtype): df_wi["time_utc"] = df_wi["time_utc"].dt.tz_localize("UTC") # Get starttime_utc and endtime_utc from h_dict @@ -119,7 +157,7 @@ def __init__(self, h_dict): f"in the wind input file ({max_time})" ) - # Set starttime_utc (zero_time_utc is redundant since time=0 corresponds to starttime_utc) + # Set starttime_utc self.starttime_utc = starttime_utc # Determine the dt implied by the weather file @@ -130,17 +168,264 @@ def __init__(self, h_dict): self.logger.info(f"dt_wi = {self.dt_wi}") self.logger.info(f"dt = {self.dt}") + self.logger.info("Interpolating wind input file...") + # Interpolate df_wi on to the time steps time_steps_all = np.arange(self.starttime, self.endtime, self.dt) df_wi = interpolate_df(df_wi, time_steps_all) - # FLORIS PREPARATION + # INITIALIZE FLORIS BASED ON WAKE MODEL + if self.wake_model == "precomputed": + self._init_floris_precomputed(df_wi) + elif self.wake_model == "dynamic": + self._init_floris_dynamic(df_wi) + else: # wake_model == "no_added_wakes" + self._init_floris_none(df_wi) + + # Common post-FLORIS initialization + self.logger.info("Initializing turbines...") + + # Get the turbine information + self.turbine_dict = load_yaml(self.turbine_file_name) + self.turbine_model_type = self.turbine_dict["turbine_model_type"] + + # Initialize the turbine array + if self.turbine_model_type == "filter_model": + # Use vectorized implementation for improved performance + self.turbine_array = TurbineFilterModelVectorized( + self.turbine_dict, self.dt, self.fmodel, self.wind_speeds_withwakes + ) + self.use_vectorized_turbines = True + elif self.turbine_model_type == "dof1_model": + self.turbine_array = [ + Turbine1dofModel( + self.turbine_dict, self.dt, self.fmodel, self.wind_speeds_withwakes[t_idx] + ) + for t_idx in range(self.n_turbines) + ] + self.use_vectorized_turbines = False + else: + raise Exception("Turbine model type should be either filter_model or dof1_model") + + # Initialize the power array to the initial wind speeds + if self.use_vectorized_turbines: + self.turbine_powers = self.turbine_array.prev_powers.copy() + else: + self.turbine_powers = np.array( + [self.turbine_array[t_idx].prev_power for t_idx in range(self.n_turbines)], + dtype=hercules_float_type, + ) + + # Get the rated power of the turbines + if self.use_vectorized_turbines: + self.rated_turbine_power = self.turbine_array.get_rated_power() + else: + self.rated_turbine_power = self.turbine_array[0].get_rated_power() + + # Get the capacity of the farm + self.capacity = self.n_turbines * self.rated_turbine_power + + # Update the user + self.logger.info( + f"Initialized WindFarm with {self.n_turbines} turbines " + f"(wake_model='{self.wake_model}')" + ) + + def _infer_wake_model_from_component_type(self, h_dict): + """Infer wake_model from component_type for backward compatibility. + + Args: + h_dict (dict): Dictionary containing simulation parameters. + + Returns: + str: Inferred wake_model ("dynamic", "precomputed", or "no_added_wakes"). + """ + component_type = h_dict[self.component_name].get("component_type", "WindFarm") + + if component_type == "Wind_MesoToPower": + return "dynamic" + elif component_type == "Wind_MesoToPowerPrecomFloris": + return "precomputed" + elif component_type == "Wind_MesoToPowerNoAddedWakes": + return "no_added_wakes" + else: + # Default to dynamic for unknown types + return "dynamic" + + def _init_floris_precomputed(self, df_wi): + """Initialize FLORIS with precomputed wake deficits. + + Args: + df_wi (pd.DataFrame): Interpolated wind input dataframe. + """ + self.logger.info("Initializing FLORIS (precomputed mode)...") + + # Initialize the FLORIS model as an ApproxFlorisModel + self.fmodel = ApproxFlorisModel( + self.floris_input_file, + wd_resolution=1.0, + ws_resolution=1.0, + ) + + # Get the layout and number of turbines from FLORIS + self.layout_x = self.fmodel.layout_x + self.layout_y = self.fmodel.layout_y + self.n_turbines = self.fmodel.n_turbines + + self.logger.info("Converting wind input file to numpy matrices...") + + # Convert the wind directions and wind speeds and ti to numpy matrices + if "ws_mean" in df_wi.columns and "ws_000" not in df_wi.columns: + self.ws_mat = np.tile( + df_wi["ws_mean"].values.astype(hercules_float_type)[:, np.newaxis], + (1, self.n_turbines), + ) + else: + self.ws_mat = df_wi[[f"ws_{t_idx:03d}" for t_idx in range(self.n_turbines)]].to_numpy( + dtype=hercules_float_type + ) + + # Compute the turbine averaged wind speeds (axis = 1) using mean + self.ws_mat_mean = np.mean(self.ws_mat, axis=1, dtype=hercules_float_type) + + self.initial_wind_speeds = self.ws_mat[0, :] + self.wind_speed_mean_background = self.ws_mat_mean[0] + + # For now require "wd_mean" to be in the df_wi + if "wd_mean" not in df_wi.columns: + raise ValueError("Wind input file must contain a column called 'wd_mean'") + self.wd_mat_mean = df_wi["wd_mean"].values.astype(hercules_float_type) + + if "ti_000" in df_wi.columns: + self.ti_mat = df_wi[[f"ti_{t_idx:03d}" for t_idx in range(self.n_turbines)]].to_numpy( + dtype=hercules_float_type + ) + + # Compute the turbine averaged turbulence intensities (axis = 1) using mean + self.ti_mat_mean = np.mean(self.ti_mat, axis=1, dtype=hercules_float_type) + + self.initial_tis = self.ti_mat[0, :] + + else: + self.ti_mat_mean = 0.08 * np.ones_like(self.ws_mat_mean, dtype=hercules_float_type) + + # Precompute the wake deficits at the cadence specified by floris_update_time_s + self.logger.info("Precomputing FLORIS wake deficits...") + + # Derived step count + self.floris_update_steps = max(1, int(self.floris_update_time_s / self.dt)) + + # Determine update step cadence and indices to evaluate FLORIS + update_steps = self.floris_update_steps + n_steps = len(self.ws_mat_mean) + eval_indices = np.arange(update_steps - 1, n_steps, update_steps) + # Ensure at least the final time is evaluated + if eval_indices.size == 0: + eval_indices = np.array([n_steps - 1]) + elif eval_indices[-1] != n_steps - 1: + eval_indices = np.append(eval_indices, n_steps - 1) + + # Build right-aligned windowed means for ws, wd, ti at the evaluation indices + def window_mean(arr_1d, idx, win): + start = max(0, idx - win + 1) + return np.mean(arr_1d[start : idx + 1], dtype=hercules_float_type) + + def window_circmean(arr_1d, idx, win): + start = max(0, idx - win + 1) + return circmean(arr_1d[start : idx + 1], high=360.0, low=0.0, nan_policy="omit") + + ws_eval = np.array( + [window_mean(self.ws_mat_mean, i, update_steps) for i in eval_indices], + dtype=hercules_float_type, + ) + wd_eval = np.array( + [window_circmean(self.wd_mat_mean, i, update_steps) for i in eval_indices], + dtype=hercules_float_type, + ) + if np.isscalar(self.ti_mat_mean): + ti_eval = self.ti_mat_mean * np.ones_like(ws_eval, dtype=hercules_float_type) + else: + ti_eval = np.array( + [window_mean(self.ti_mat_mean, i, update_steps) for i in eval_indices], + dtype=hercules_float_type, + ) + + # Evaluate FLORIS at the evaluation cadence + self.fmodel.set( + wind_directions=wd_eval, + wind_speeds=ws_eval, + turbulence_intensities=ti_eval, + ) + self.logger.info("Running FLORIS...") + self.fmodel.run() + self.num_floris_calcs = 1 + self.logger.info("FLORIS run complete") + + # TODO: THIS CODE WILL WORK IN THE FUTURE + # https://github.com/NREL/floris/pull/1135 + # floris_velocities = self.fmodel.turbine_average_velocities + + # For now compute in place here (replace later) + expanded_velocities = average_velocity( + velocities=self.fmodel.fmodel_expanded.core.flow_field.u, + method=self.fmodel.fmodel_expanded.core.grid.average_method, + cubature_weights=self.fmodel.fmodel_expanded.core.grid.cubature_weights, + ) + + floris_velocities = map_turbine_powers_uncertain( + unique_turbine_powers=expanded_velocities, + map_to_expanded_inputs=self.fmodel.map_to_expanded_inputs, + weights=self.fmodel.weights, + n_unexpanded=self.fmodel.n_unexpanded, + n_sample_points=self.fmodel.n_sample_points, + n_turbines=self.fmodel.n_turbines, + ).astype(hercules_float_type) + + # Determine the free_stream velocities as the maximum velocity in each row + free_stream_velocities = np.tile( + np.max(floris_velocities, axis=1)[:, np.newaxis], (1, self.n_turbines) + ).astype(hercules_float_type) + + # Compute wake deficits at evaluation times + floris_wake_deficits_eval = free_stream_velocities - floris_velocities + + # Expand the wake deficits to all time steps by holding constant within each interval + deficits_all = np.zeros_like(self.ws_mat, dtype=hercules_float_type) + # For each block, fill with the corresponding deficits + prev_end = -1 + for block_idx, end_idx in enumerate(eval_indices): + start_idx = prev_end + 1 + prev_end = end_idx + # Use deficits from this evaluation time for the whole block + deficits_all[start_idx : end_idx + 1, :] = floris_wake_deficits_eval[block_idx, :] + + # Compute all the withwakes wind speeds from background minus deficits + self.wind_speeds_withwakes_all = self.ws_mat - deficits_all + + # Initialize the turbine powers to nan + self.turbine_powers = np.zeros(self.n_turbines, dtype=hercules_float_type) * np.nan + + # Get the initial background wind speeds + self.wind_speeds_background = self.ws_mat[0, :] + + # Compute initial withwakes wind speeds + self.wind_speeds_withwakes = self.wind_speeds_withwakes_all[0, :] + + # Get the initial FLORIS wake deficits + self.floris_wake_deficits = self.wind_speeds_background - self.wind_speeds_withwakes + + def _init_floris_dynamic(self, df_wi): + """Initialize FLORIS for dynamic wake calculation. + + Args: + df_wi (pd.DataFrame): Interpolated wind input dataframe. + """ + self.logger.info("Initializing FLORIS (dynamic mode)...") # Initialize the FLORIS model self.fmodel = FlorisModel(self.floris_input_file) - # Change to the simple-derating model turbine - # (Note this could also be done with the mixed model) + # Change to the mixed operation model self.fmodel.set_operation_model("mixed") # Get the layout and number of turbines from FLORIS @@ -162,8 +447,7 @@ def __init__(self, h_dict): # Add an initial non-nan value to be over-written on first step self.turbine_power_setpoints_buffer[0, :] = 1e12 - # Convert the wind directions and wind speeds and ti to simply numpy matrices - # Starting with wind speeds + # Convert the wind directions and wind speeds and ti to numpy matrices if "ws_mean" in df_wi.columns and "ws_000" not in df_wi.columns: self.ws_mat = np.tile( df_wi["ws_mean"].values.astype(hercules_float_type)[:, np.newaxis], @@ -185,25 +469,7 @@ def __init__(self, h_dict): raise ValueError("Wind input file must contain a column called 'wd_mean'") self.wd_mat_mean = df_wi["wd_mean"].values.astype(hercules_float_type) - # OLD APPROACH - # # Now the wind directions - # if "wd_000" in df_wi.columns: - # self.wd_mat = df_wi[ - # [f"wd_{t_idx:03d}" for t_idx in range(self.n_turbines)] - # ].to_numpy() - - # # Compute the turbine-averaged wind directions (axis = 1) using circmean - # self.wd_mat_mean = np.apply_along_axis( - # lambda x: circmean(x, high=360.0, low=0.0, nan_policy="omit"), - # axis=1, - # arr=self.wd_mat, - # ) - - # self.initial_wind_directions = self.wd_mat[0, :] - # elif "wd_mean" in df_wi.columns: - # self.wd_mat_mean = df_wi["wd_mean"].values - - # Compute the initial floris wind direction and wind speed as at the start index + # Compute the initial floris wind direction self.floris_wind_direction = self.wd_mat_mean[0] if "ti_000" in df_wi.columns: @@ -233,56 +499,79 @@ def __init__(self, h_dict): self.turbine_powers = np.zeros(self.n_turbines, dtype=hercules_float_type) * np.nan # Get the initial background wind speeds - # TODO: This is more a debugging thing, not really necessary self.wind_speeds_background = self.ws_mat[0, :] - # # Compute the initial waked wind speeds + # Compute the initial waked wind speeds self.update_wake_deficits(step=0) # Compute withwakes wind speeds self.wind_speeds_withwakes = self.ws_mat[0, :] - self.floris_wake_deficits - # Get the turbine information - self.turbine_dict = load_yaml(self.turbine_file_name) - self.turbine_model_type = self.turbine_dict["turbine_model_type"] + def _init_floris_none(self, df_wi): + """Initialize without wake modeling. - # Initialize the turbine array - if self.turbine_model_type == "filter_model": - # Use vectorized implementation for improved performance - self.turbine_array = TurbineFilterModelVectorized( - self.turbine_dict, self.dt, self.fmodel, self.wind_speeds_withwakes + Args: + df_wi (pd.DataFrame): Interpolated wind input dataframe. + """ + self.logger.info("Initializing FLORIS (no wake modeling)...") + + # Initialize the FLORIS model (still needed for turbine power curve) + self.fmodel = FlorisModel(self.floris_input_file) + + # Get the layout and number of turbines from FLORIS + self.layout_x = self.fmodel.layout_x + self.layout_y = self.fmodel.layout_y + self.n_turbines = self.fmodel.n_turbines + + # floris_update_steps not used but set for consistency + self.floris_update_steps = max(1, int(self.floris_update_time_s / self.dt)) + + # Convert the wind directions and wind speeds and ti to numpy matrices + if "ws_mean" in df_wi.columns and "ws_000" not in df_wi.columns: + self.ws_mat = np.tile( + df_wi["ws_mean"].values.astype(hercules_float_type)[:, np.newaxis], + (1, self.n_turbines), ) - self.use_vectorized_turbines = True - elif self.turbine_model_type == "dof1_model": - self.turbine_array = [ - Turbine1dofModel( - self.turbine_dict, self.dt, self.fmodel, self.wind_speeds_withwakes[t_idx] - ) - for t_idx in range(self.n_turbines) - ] - self.use_vectorized_turbines = False else: - raise Exception("Turbine model type should be either filter_model or dof1_model") + self.ws_mat = df_wi[[f"ws_{t_idx:03d}" for t_idx in range(self.n_turbines)]].to_numpy( + dtype=hercules_float_type + ) - # Initialize the power array to the initial wind speeds - if self.use_vectorized_turbines: - self.turbine_powers = self.turbine_array.prev_powers.copy() - else: - self.turbine_powers = np.array( - [self.turbine_array[t_idx].prev_power for t_idx in range(self.n_turbines)] - ).astype(hercules_float_type) + # Compute the turbine averaged wind speeds (axis = 1) using mean + self.ws_mat_mean = np.mean(self.ws_mat, axis=1, dtype=hercules_float_type) + + self.initial_wind_speeds = self.ws_mat[0, :] + self.wind_speed_mean_background = self.ws_mat_mean[0] + + # For now require "wd_mean" to be in the df_wi + if "wd_mean" not in df_wi.columns: + raise ValueError("Wind input file must contain a column called 'wd_mean'") + self.wd_mat_mean = df_wi["wd_mean"].values.astype(hercules_float_type) + + if "ti_000" in df_wi.columns: + self.ti_mat = df_wi[[f"ti_{t_idx:03d}" for t_idx in range(self.n_turbines)]].to_numpy( + dtype=hercules_float_type + ) + + # Compute the turbine averaged turbulence intensities (axis = 1) using mean + self.ti_mat_mean = np.mean(self.ti_mat, axis=1, dtype=hercules_float_type) + + self.initial_tis = self.ti_mat[0, :] - # Get the rated power of the turbines, for now assume all turbines have the same rated power - if self.use_vectorized_turbines: - self.rated_turbine_power = self.turbine_array.get_rated_power() else: - self.rated_turbine_power = self.turbine_array[0].get_rated_power() + self.ti_mat_mean = 0.08 * np.ones_like(self.ws_mat_mean, dtype=hercules_float_type) - # Get the capacity of the farm - self.capacity = self.n_turbines * self.rated_turbine_power + # No wake deficits + self.floris_wake_deficits = np.zeros(self.n_turbines, dtype=hercules_float_type) - # Update the user - self.logger.info(f"Initialized Wind_MesoToPower with {self.n_turbines} turbines") + # Initialize the turbine powers to nan + self.turbine_powers = np.zeros(self.n_turbines, dtype=hercules_float_type) * np.nan + + # Get the initial background wind speeds + self.wind_speeds_background = self.ws_mat[0, :] + + # No wakes: withwakes == background + self.wind_speeds_withwakes = self.wind_speeds_background.copy() def get_initial_conditions_and_meta_data(self, h_dict): """Add any initial conditions or meta data to the h_dict. @@ -311,7 +600,7 @@ def get_initial_conditions_and_meta_data(self, h_dict): return h_dict def update_wake_deficits(self, step): - """Update the wake deficits in the FLORIS model based on the current simulation step. + """Update the wake deficits in the FLORIS model (dynamic mode only). This method computes the necessary FLORIS inputs (wind direction, wind speed, turbulence intensity, and power_setpoints) over a specified time window. If any of these @@ -325,7 +614,6 @@ def update_wake_deficits(self, step): window_start = max(0, step - self.floris_update_steps + 1) # Compute new values of the floris inputs - # TODO: CONFIRM THE +1 in the slice is right self.floris_wind_direction = circmean( self.wd_mat_mean[window_start : step + 1], high=360.0, low=0.0, nan_policy="omit" ) @@ -356,7 +644,7 @@ def update_wake_deficits(self, step): self.num_floris_calcs += 1 def update_power_setpoints_buffer(self, turbine_power_setpoints): - """Update the power_setpoints buffer with the turbine_power_setpoints values. + """Update the power_setpoints buffer (dynamic mode only). This method stores the given power setpoint values in the current position of the power_setpoints buffer and updates the index to point to the next position in a @@ -379,9 +667,8 @@ def update_power_setpoints_buffer(self, turbine_power_setpoints): def step(self, h_dict): """Execute one simulation step for the wind farm. - Updates wake deficits, computes waked velocities, calculates turbine powers, - and updates the simulation dictionary with results. Handles power_setpoint - signals and optional extra logging outputs. + Updates wake deficits (if applicable), computes waked velocities, calculates + turbine powers, and updates the simulation dictionary with results. Args: h_dict (dict): Dictionary containing current simulation state including @@ -396,22 +683,37 @@ def step(self, h_dict): if self.verbose: self.logger.info(f"step = {step} (of {self.n_steps})") - # Grab the instantaneous turbine power setpoint signal and update the power_setpoints buffer + # Grab the instantaneous turbine power setpoint signal turbine_power_setpoints = h_dict[self.component_name]["turbine_power_setpoints"] - self.update_power_setpoints_buffer(turbine_power_setpoints) - # Get the background wind speeds - # TODO: This is more a debugging thing, not really necessary - self.wind_speeds_background = self.ws_mat[step, :] + # Update wind speeds based on wake model + if self.wake_model == "dynamic": + # Update power setpoints buffer + self.update_power_setpoints_buffer(turbine_power_setpoints) - # Check if it is time to update the withwakes wind speeds - if step % self.floris_update_steps == 0: - self.update_wake_deficits(step) + # Get the background wind speeds + self.wind_speeds_background = self.ws_mat[step, :] - # Compute withwakes wind speeds - self.wind_speeds_withwakes = self.ws_mat[step, :] - self.floris_wake_deficits + # Check if it is time to update the withwakes wind speeds + if step % self.floris_update_steps == 0: + self.update_wake_deficits(step) + + # Compute withwakes wind speeds + self.wind_speeds_withwakes = self.ws_mat[step, :] - self.floris_wake_deficits - # Update the turbine powers + elif self.wake_model == "precomputed": + # Update all the wind speeds + self.wind_speeds_background = self.ws_mat[step, :] + self.wind_speeds_withwakes = self.wind_speeds_withwakes_all[step, :] + self.floris_wake_deficits = self.wind_speeds_background - self.wind_speeds_withwakes + + else: # wake_model == "no_added_wakes" + # No wake modeling - use background speeds directly + self.wind_speeds_background = self.ws_mat[step, :] + self.wind_speeds_withwakes = self.wind_speeds_background.copy() + self.floris_wake_deficits = np.zeros(self.n_turbines, dtype=hercules_float_type) + + # Update the turbine powers (common for all wake models) if self.use_vectorized_turbines: # Vectorized calculation for all turbines at once self.turbine_powers = self.turbine_array.step( @@ -441,113 +743,8 @@ def step(self, h_dict): ) h_dict[self.component_name]["wind_speeds_withwakes"] = self.wind_speeds_withwakes h_dict[self.component_name]["wind_speeds_background"] = self.wind_speeds_background - return h_dict - - -class TurbineFilterModel: - """Simple filter-based wind turbine model for power output simulation. - - This model simulates wind turbine power output using a first-order filter - to smooth the response to changing wind conditions, providing a simplified - representation of turbine dynamics. - - NOTE: This class is now unused and kept for backward compatibility. - The filter_model turbine_model_type now uses TurbineFilterModelVectorized - for improved performance. - """ - - def __init__(self, turbine_dict, dt, fmodel, initial_wind_speed): - """Initialize the turbine filter model. - - Args: - turbine_dict (dict): Dictionary containing turbine configuration, - including filter model parameters and other turbine-specific data. - dt (float): Time step for the simulation in seconds. - fmodel (FlorisModel): FLORIS model of the farm. - initial_wind_speed (float): Initial wind speed in m/s to initialize - the simulation. - """ - # Save the time step - self.dt = dt - - # Save the turbine dict - self.turbine_dict = turbine_dict - # Save the filter time constant - self.filter_time_constant = turbine_dict["filter_model"]["time_constant"] - - # Solve for the filter alpha value given dt and the time constant - self.alpha = 1 - np.exp(-self.dt / self.filter_time_constant) - - # Grab the wind speed power curve from the fmodel and define a simple 1D LUT - turbine_type = fmodel.core.farm.turbine_definitions[0] - wind_speeds = turbine_type["power_thrust_table"]["wind_speed"] - powers = turbine_type["power_thrust_table"]["power"] - self.power_lut = interp1d( - wind_speeds, - powers, - fill_value=0.0, - bounds_error=False, - ) - - # Initialize the previous power to the initial wind speed - self.prev_power = self.power_lut(initial_wind_speed) - - def get_rated_power(self): - """Get the rated power of the turbine. - - Returns: - float: The rated power of the turbine in kW. - """ - return np.max(self.power_lut(np.arange(0, 25, 1.0, dtype=hercules_float_type))) - - def step(self, wind_speed, power_setpoint): - """Simulate a single time step of the wind turbine power output. - - This method calculates the power output of a wind turbine based on the - given wind speed and power_setpoint. The power output is - smoothed using an exponential moving average to simulate the turbine's - response to changing wind conditions. - - Args: - wind_speed (float): The current wind speed in meters per second (m/s). - power_setpoint (float): The maximum allowable power output in kW. - - Returns: - float: The calculated power output of the wind turbine, constrained - by the power_setpoint and smoothed using the exponential moving average. - """ - # Instantaneous power - instant_power = self.power_lut(wind_speed) - - # Limit the current power to not be greater than power_setpoint - instant_power = min(instant_power, power_setpoint) - - # Limit the instant power to be greater than 0 - instant_power = max(instant_power, 0.0) - - # TEMP: why are NaNs occurring? - if np.isnan(instant_power): - print( - f"NaN instant power at wind speed {wind_speed} m/s, " - f"power setpoint {power_setpoint} kW, prev power {self.prev_power} kW" - ) - instant_power = self.prev_power - - # Update the power - power = self.alpha * instant_power + (1 - self.alpha) * self.prev_power - - # Limit the power to not be greater than power_setpoint - power = min(power, power_setpoint) - - # Limit the power to be greater than 0 - power = max(power, 0.0) - - # Update the previous power - self.prev_power = power - - # Return the power - return power + return h_dict class TurbineFilterModelVectorized: @@ -637,8 +834,6 @@ def step(self, wind_speeds, power_setpoints): # Handle NaNs by replacing with previous power values nan_mask = np.isnan(instant_powers) if np.any(nan_mask): - # Log warning for NaN values (but don't print every occurrence for performance) - # Could add logging here if needed instant_powers[nan_mask] = self.prev_powers[nan_mask] # Vectorized exponential filter update @@ -753,7 +948,6 @@ def step(self, wind_speed, power_setpoint): / self.turbine_dict["dof1_model"]["rotor_inertia"] ) omegaf = (1 - self.filteralpha) * omega + self.filteralpha * (self.prev_omegaf) - # print(omegaf-omega) pitch, gentq = self.simplecontroller(wind_speed, omegaf) tsr = float(omegaf * self.rotor_radius / wind_speed) if power_setpoint > 0: @@ -785,9 +979,6 @@ def step(self, wind_speed, power_setpoint): * self.perffuncs["Cq"]([tsr, pitch]) ) - # power = ( - # self.perffuncs["Cp"]([tsr, pitch]) * 0.5 * self.rho * self.rotor_area * wind_speed**3 - # ) power = gentq * omega * self.turbine_dict["dof1_model"]["gearbox_ratio"] self.prev_omega = omega @@ -813,9 +1004,6 @@ def simplecontroller(self, wind_speed, omegaf): tuple: (pitch_angle, generator_torque) where pitch is in radians and generator torque is in N⋅m. """ - # if omega <= self.turbine_dict['dof1_model']['rated_wind_speed']: pitch = 0.0 gentorque = self.turbine_dict["dof1_model"]["controller"]["r2_k_torque"] * omegaf**2 - # else - # raise Exception("Region-3 controller not implemented yet") return pitch, gentorque diff --git a/hercules/plant_components/wind_farm_scada_power.py b/hercules/plant_components/wind_farm_scada_power.py new file mode 100644 index 00000000..3759d3fa --- /dev/null +++ b/hercules/plant_components/wind_farm_scada_power.py @@ -0,0 +1,345 @@ +# Unified wind farm model for Hercules supporting multiple wake modeling strategies. + +import numpy as np +import pandas as pd +from hercules.plant_components.component_base import ComponentBase +from hercules.utilities import ( + hercules_float_type, + interpolate_df, +) + + +class WindFarmSCADAPower(ComponentBase): + """Wind farm model that uses SCADA power data to simulate wind farm performance.""" + + def __init__(self, h_dict, wake_model=None): + """Initialize the WindFarm class. + + Args: + h_dict (dict): Dictionary containing simulation parameters. + wake_model (str, optional): Wake modeling strategy: "dynamic", "precomputed", + or "no_added_wakes". If None, infers from component_type for backward compatibility. + Defaults to None. + + Raises: + ValueError: If wake_model is invalid or required parameters are missing. + """ + # Store the name of this component + self.component_name = "wind_farm" + + self.component_type = "WindFarmSCADAPower" + + # Call the base class init + super().__init__(h_dict, self.component_name) + + self.logger.info("Initializing WindFarmSCADAPower") + + # Track the number of FLORIS calculations + self.num_floris_calcs = 0 + + # Read in the input file names + self.scada_filename = h_dict[self.component_name]["scada_filename"] + self.turbine_file_name = h_dict[self.component_name]["turbine_file_name"] + + self.logger.info("Reading in SCADA power data...") + + # Read in the scada file + if self.scada_filename.endswith(".csv"): + df_scada = pd.read_csv(self.scada_filename) + elif self.scada_filename.endswith(".p") | self.scada_filename.endswith(".pkl"): + df_scada = pd.read_pickle(self.scada_filename) + elif (self.scada_filename.endswith(".f")) | (self.scada_filename.endswith(".ftr")): + df_scada = pd.read_feather(self.scada_filename) + else: + raise ValueError("SCADA file must be a .csv or .p, .f or .ftr file") + + self.logger.info("Checking SCADA file...") + + # Make sure the df_scada contains a column called "time_utc" + if "time_utc" not in df_scada.columns: + raise ValueError("SCADA file must contain a column called 'time_utc'") + + # Convert time_utc to datetime if it's not already + if not pd.api.types.is_datetime64_any_dtype(df_scada["time_utc"]): + # Strip whitespace from time_utc values to handle CSV formatting issues + df_scada["time_utc"] = df_scada["time_utc"].astype(str).str.strip() + try: + df_scada["time_utc"] = pd.to_datetime( + df_scada["time_utc"], format="ISO8601", utc=True + ) + except (ValueError, TypeError): + # If ISO8601 format fails, try parsing without specifying format + df_scada["time_utc"] = pd.to_datetime(df_scada["time_utc"], utc=True) + + # Ensure time_utc is timezone-aware (UTC) + if not isinstance(df_scada["time_utc"].dtype, pd.DatetimeTZDtype): + df_scada["time_utc"] = df_scada["time_utc"].dt.tz_localize("UTC") + + # Get starttime_utc and endtime_utc from h_dict + starttime_utc = h_dict["starttime_utc"] + endtime_utc = h_dict["endtime_utc"] + + # Ensure starttime_utc is timezone-aware (UTC) + if not isinstance(starttime_utc, pd.Timestamp): + starttime_utc = pd.to_datetime(starttime_utc, utc=True) + elif starttime_utc.tz is None: + starttime_utc = starttime_utc.tz_localize("UTC") + + # Ensure endtime_utc is timezone-aware (UTC) + if not isinstance(endtime_utc, pd.Timestamp): + endtime_utc = pd.to_datetime(endtime_utc, utc=True) + elif endtime_utc.tz is None: + endtime_utc = endtime_utc.tz_localize("UTC") + + # Generate time column internally: time = 0 corresponds to starttime_utc + df_scada["time"] = (df_scada["time_utc"] - starttime_utc).dt.total_seconds() + + # Validate that starttime_utc and endtime_utc are within the time_utc range + if df_scada["time_utc"].min() > starttime_utc: + min_time = df_scada["time_utc"].min() + raise ValueError( + f"Start time UTC {starttime_utc} is before the earliest time " + f"in the SCADA file ({min_time})" + ) + if df_scada["time_utc"].max() < endtime_utc: + max_time = df_scada["time_utc"].max() + raise ValueError( + f"End time UTC {endtime_utc} is after the latest time " + f"in the SCADA file ({max_time})" + ) + + # Set starttime_utc + self.starttime_utc = starttime_utc + + # Determine the dt implied by the weather file + self.dt_scada = df_scada["time"][1] - df_scada["time"][0] + + # Log the values + if self.verbose: + self.logger.info(f"dt_scada = {self.dt_scada}") + self.logger.info(f"dt = {self.dt}") + + self.logger.info("Interpolating SCADA file...") + + # Interpolate df_scada on to the time steps + time_steps_all = np.arange(self.starttime, self.endtime, self.dt) + df_scada = interpolate_df(df_scada, time_steps_all) + + # Get a list of power columns and infer number of turbines + self.power_columns = sorted([col for col in df_scada.columns if col.startswith("pow_")]) + + # Infer the number of turbines by the number of power columns + self.n_turbines = len(self.power_columns) + + self.logger.info(f"Inferred number of turbines: {self.n_turbines}") + + # Collect the turbine powers + self.scada_powers = df_scada[self.power_columns].to_numpy(dtype=hercules_float_type) + + # Now get the wind speed and directions + + # Convert the wind directions and wind speeds and ti to numpy matrices + if "ws_mean" in df_scada.columns and "ws_000" not in df_scada.columns: + self.ws_mat = np.tile( + df_scada["ws_mean"].values.astype(hercules_float_type)[:, np.newaxis], + (1, self.n_turbines), + ) + else: + self.ws_mat = df_scada[ + [f"ws_{t_idx:03d}" for t_idx in range(self.n_turbines)] + ].to_numpy(dtype=hercules_float_type) + + # Compute the turbine averaged wind speeds (axis = 1) using mean + self.ws_mat_mean = np.mean(self.ws_mat, axis=1, dtype=hercules_float_type) + + self.initial_wind_speeds = self.ws_mat[0, :] + self.wind_speed_mean_background = self.ws_mat_mean[0] + + # Get the initial background wind speeds + self.wind_speeds_background = self.ws_mat[0, :] + + # No wakes: withwakes == background + self.wind_speeds_withwakes = self.wind_speeds_background.copy() + + # For now require "wd_mean" to be in the df_scada + if "wd_mean" not in df_scada.columns: + raise ValueError("SCADA file must contain a column called 'wd_mean'") + self.wd_mat_mean = df_scada["wd_mean"].values.astype(hercules_float_type) + + if "ti_000" in df_scada.columns: + self.ti_mat = df_scada[ + [f"ti_{t_idx:03d}" for t_idx in range(self.n_turbines)] + ].to_numpy(dtype=hercules_float_type) + + # Compute the turbine averaged turbulence intensities (axis = 1) using mean + self.ti_mat_mean = np.mean(self.ti_mat, axis=1, dtype=hercules_float_type) + + self.initial_tis = self.ti_mat[0, :] + + else: + self.ti_mat_mean = 0.08 * np.ones_like(self.ws_mat_mean, dtype=hercules_float_type) + + # No wake deficits + self.floris_wake_deficits = np.zeros(self.n_turbines, dtype=hercules_float_type) + + # Infer the rated power as the 99 percentile of the power column of 0th turbine + self.rated_turbine_power = np.percentile(df_scada[self.power_columns[0]], 99) + + # Get the capacity of the farm + self.capacity = self.n_turbines * self.rated_turbine_power + + self.logger.info(f"Inferred rated turbine power: {self.rated_turbine_power}") + self.logger.info(f"Inferred capacity: {self.capacity/1E3} MW") + + # Initialize the turbine array + self.turbine_array = TurbineFilterModelVectorizedSCADA(self.dt, self.scada_powers[0, :]) + + # Initialize the turbine powers to the starting row + self.turbine_powers = self.turbine_array.prev_powers.copy() + + def get_initial_conditions_and_meta_data(self, h_dict): + """Add any initial conditions or meta data to the h_dict. + + Meta data is data not explicitly in the input yaml but still useful for other + modules. + + Args: + h_dict (dict): Dictionary containing simulation parameters. + + Returns: + dict: Dictionary containing simulation parameters with initial conditions and meta data. + """ + h_dict["wind_farm"]["n_turbines"] = self.n_turbines + h_dict["wind_farm"]["capacity"] = self.capacity + h_dict["wind_farm"]["rated_turbine_power"] = self.rated_turbine_power + h_dict["wind_farm"]["wind_direction_mean"] = self.wd_mat_mean[0] + h_dict["wind_farm"]["wind_speed_mean_background"] = self.ws_mat_mean[0] + h_dict["wind_farm"]["turbine_powers"] = self.turbine_powers + h_dict["wind_farm"]["power"] = np.sum(self.turbine_powers) + + # Log the start time UTC if available + if hasattr(self, "starttime_utc"): + h_dict["wind_farm"]["starttime_utc"] = self.starttime_utc + + return h_dict + + def step(self, h_dict): + """Execute one simulation step for the wind farm. + + Updates wake deficits (if applicable), computes waked velocities, calculates + turbine powers, and updates the simulation dictionary with results. + + Args: + h_dict (dict): Dictionary containing current simulation state including + step number and power_setpoint values for each turbine. + + Returns: + dict: Updated simulation dictionary with wind farm outputs including + turbine powers, total power, and optional extra outputs. + """ + # Get the current step + step = h_dict["step"] + if self.verbose: + self.logger.info(f"step = {step} (of {self.n_steps})") + + # Grab the instantaneous turbine power setpoint signal + turbine_power_setpoints = h_dict[self.component_name]["turbine_power_setpoints"] + + # Update wind speeds based on wake model + + # No wake modeling - use background speeds directly + self.wind_speeds_background = self.ws_mat[step, :] + self.wind_speeds_withwakes = self.wind_speeds_background.copy() + self.floris_wake_deficits = np.zeros(self.n_turbines, dtype=hercules_float_type) + + # Update the turbine powers (common for all wake models) + # Vectorized calculation for all turbines at once + self.turbine_powers = self.turbine_array.step( + self.scada_powers[step, :], + turbine_power_setpoints, + ) + + # Update instantaneous wind direction and wind speed + self.wind_direction_mean = self.wd_mat_mean[step] + self.wind_speed_mean_background = self.ws_mat_mean[step] + + # Update the h_dict with outputs + h_dict[self.component_name]["power"] = np.sum(self.turbine_powers) + h_dict[self.component_name]["turbine_powers"] = self.turbine_powers + h_dict[self.component_name]["turbine_power_setpoints"] = turbine_power_setpoints + h_dict[self.component_name]["wind_direction_mean"] = self.wind_direction_mean + h_dict[self.component_name]["wind_speed_mean_background"] = self.wind_speed_mean_background + h_dict[self.component_name]["wind_speed_mean_withwakes"] = np.mean( + self.wind_speeds_withwakes, dtype=hercules_float_type + ) + h_dict[self.component_name]["wind_speeds_withwakes"] = self.wind_speeds_withwakes + h_dict[self.component_name]["wind_speeds_background"] = self.wind_speeds_background + + return h_dict + + +class TurbineFilterModelVectorizedSCADA: + """Vectorized filter-based wind turbine model for power output simulation. + + This model does not use a filter so that SCADA data can be used directly. + """ + + def __init__(self, dt, initial_scada_powers): + """Initialize the vectorized turbine filter model. + + Args: + dt (float): Time step for the simulation in seconds. + initial_scada_powers (np.ndarray): Initial wind speeds in m/s for all turbines + to initialize the simulation. + """ + # Save the time step + self.dt = dt + + # Number of turbines + self.n_turbines = len(initial_scada_powers) + + # Initialize the previous powers for all turbines + self.prev_powers = initial_scada_powers.copy() + + def step(self, scada_powers, power_setpoints): + """Simulate a single time step for all wind turbines simultaneously. + + This method calculates the power output of all wind turbines based on the + given wind speeds and power setpoints. The power outputs are + smoothed using an exponential moving average to simulate the turbines' + response to changing wind conditions. + + Args: + scada_powers (np.ndarray): Current SCADA powers for all turbines. + power_setpoints (np.ndarray): Maximum allowable power outputs in kW for all turbines. + + Returns: + np.ndarray: Calculated power outputs of all wind turbines, constrained + by the power setpoints and smoothed using the exponential moving average. + """ + + # Vectorized limiting: current power not greater than power_setpoint + instant_powers = np.minimum(scada_powers, power_setpoints) + + # Vectorized limiting: instant power not less than 0 + instant_powers = np.maximum(instant_powers, 0.0) + + # Handle NaNs by replacing with previous power values + nan_mask = np.isnan(instant_powers) + if np.any(nan_mask): + instant_powers[nan_mask] = self.prev_powers[nan_mask] + + # Simple update without any filtering + powers = instant_powers + + # Vectorized limiting: power not greater than power_setpoint + powers = np.minimum(powers, power_setpoints) + + # Vectorized limiting: power not less than 0 + powers = np.maximum(powers, 0.0) + + # Update the previous powers for all turbines + self.prev_powers = powers.copy() + + # Return the powers + return powers diff --git a/hercules/plant_components/wind_meso_to_power_precom_floris.py b/hercules/plant_components/wind_meso_to_power_precom_floris.py deleted file mode 100644 index 42dd14fe..00000000 --- a/hercules/plant_components/wind_meso_to_power_precom_floris.py +++ /dev/null @@ -1,572 +0,0 @@ -# Implements the meso-scale wind model for Hercules. - - -import numpy as np -import pandas as pd -from floris import ApproxFlorisModel -from floris.core import average_velocity -from floris.uncertain_floris_model import map_turbine_powers_uncertain -from hercules.plant_components.component_base import ComponentBase -from hercules.plant_components.wind_meso_to_power import ( - Turbine1dofModel, - TurbineFilterModelVectorized, -) -from hercules.utilities import ( - hercules_float_type, - interpolate_df, - load_yaml, -) -from scipy.interpolate import interp1d -from scipy.stats import circmean - -RPM2RADperSec = 2 * np.pi / 60.0 - - -class Wind_MesoToPowerPrecomFloris(ComponentBase): - def __init__(self, h_dict): - """Initialize the Wind_MesoToPowerPrecomFloris class. - - This model focuses on meso-scale wind phenomena by applying a separate wind speed - time signal to each turbine model derived from data. It combines FLORIS wake - modeling with detailed turbine dynamics for wind farm performance analysis. - - In contrast to the Wind_MesoToPower class, this class pre-computes the FLORIS wake - deficits for all wind speeds and wind directions. This is done by running FLORIS - once for all wind speeds and wind directions (but not for varying power setpoints). - This is valid - for cases where the wind farm is operating: - - all turbines operating normally - - all turbines off - - following a wind-farm wide derating level - - It is in practice conservative with respect to the wake deficits, but it is more efficient - than running FLORIS for each condition. In cases where turbines are: - - partially derated below the curtailment level - - not uniformly curtailed or some turbines are off - - This is not an appropriate model and the more general Wind_MesoToPower class should be used. - - Args: - h_dict (dict): Dictionary containing values for the simulation. - - Required keys in `h_dict['wind_farm']`: - - `floris_input_file` (str): Path to FLORIS configuration file. - - `wind_input_filename` (str): Path to wind input data file. - - `turbine_file_name` (str): Path to turbine configuration file. - - `floris_update_time_s` (float): Update period in seconds. This value - determines the cadence of the wake precomputation. Wind inputs are - averaged over the most recent `floris_update_time_s` and FLORIS is - evaluated at that interval. The resulting wake deficits are then held - constant until the next FLORIS update. - """ - # Store the name of this component - self.component_name = "wind_farm" - - # Store the type of this component - self.component_type = "Wind_MesoToPowerPrecomFloris" - - # Call the base class init - super().__init__(h_dict, self.component_name) - - self.logger.info("Completed base class init...") - - # Track the number of FLORIS calculations - self.num_floris_calcs = 0 - - self.logger.info("Reading in FLORIS input files...") - - # Read in the input file names - self.floris_input_file = h_dict[self.component_name]["floris_input_file"] - self.wind_input_filename = h_dict[self.component_name]["wind_input_filename"] - self.turbine_file_name = h_dict[self.component_name]["turbine_file_name"] - - # Require floris_update_time_s for interface consistency, though it is unused - if "floris_update_time_s" not in h_dict[self.component_name]: - raise ValueError("floris_update_time_s must be in the h_dict") - self.floris_update_time_s = h_dict[self.component_name]["floris_update_time_s"] - if self.floris_update_time_s < 1: - raise ValueError("FLORIS update time must be at least 1 second") - # Derived step count (not used by this precomputed model, but kept for parity) - self.floris_update_steps = max(1, int(self.floris_update_time_s / self.dt)) - - self.logger.info("Reading in wind input file...") - - # Read in the weather file data - # If a csv file is provided, read it in - if self.wind_input_filename.endswith(".csv"): - df_wi = pd.read_csv(self.wind_input_filename) - elif self.wind_input_filename.endswith(".p") | self.wind_input_filename.endswith(".pkl"): - df_wi = pd.read_pickle(self.wind_input_filename) - elif (self.wind_input_filename.endswith(".f")) | ( - self.wind_input_filename.endswith(".ftr") - ): - df_wi = pd.read_feather(self.wind_input_filename) - else: - raise ValueError("Wind input file must be a .csv or .p, .f or .ftr file") - - self.logger.info("Checking wind input file...") - - # Make sure the df_wi contains a column called "time_utc" - if "time_utc" not in df_wi.columns: - raise ValueError("Wind input file must contain a column called 'time_utc'") - - # Convert time_utc to datetime if it's not already - if not pd.api.types.is_datetime64_any_dtype(df_wi["time_utc"]): - # Strip whitespace from time_utc values to handle CSV formatting issues - df_wi["time_utc"] = df_wi["time_utc"].astype(str).str.strip() - try: - df_wi["time_utc"] = pd.to_datetime(df_wi["time_utc"], format="ISO8601", utc=True) - except (ValueError, TypeError): - # If ISO8601 format fails, try parsing without specifying format - df_wi["time_utc"] = pd.to_datetime(df_wi["time_utc"], utc=True) - - # Ensure time_utc is timezone-aware (UTC) - if not pd.api.types.is_datetime64tz_dtype(df_wi["time_utc"]): - df_wi["time_utc"] = df_wi["time_utc"].dt.tz_localize("UTC") - - # Get starttime_utc and endtime_utc from h_dict - starttime_utc = h_dict["starttime_utc"] - endtime_utc = h_dict["endtime_utc"] - - # Ensure starttime_utc is timezone-aware (UTC) - if not isinstance(starttime_utc, pd.Timestamp): - starttime_utc = pd.to_datetime(starttime_utc, utc=True) - elif starttime_utc.tz is None: - starttime_utc = starttime_utc.tz_localize("UTC") - - # Ensure endtime_utc is timezone-aware (UTC) - if not isinstance(endtime_utc, pd.Timestamp): - endtime_utc = pd.to_datetime(endtime_utc, utc=True) - elif endtime_utc.tz is None: - endtime_utc = endtime_utc.tz_localize("UTC") - - # Generate time column internally: time = 0 corresponds to starttime_utc - df_wi["time"] = (df_wi["time_utc"] - starttime_utc).dt.total_seconds() - - # Validate that starttime_utc and endtime_utc are within the time_utc range - if df_wi["time_utc"].min() > starttime_utc: - min_time = df_wi["time_utc"].min() - raise ValueError( - f"Start time UTC {starttime_utc} is before the earliest time " - f"in the wind input file ({min_time})" - ) - if df_wi["time_utc"].max() < endtime_utc: - max_time = df_wi["time_utc"].max() - raise ValueError( - f"End time UTC {endtime_utc} is after the latest time " - f"in the wind input file ({max_time})" - ) - - # Set starttime_utc (zero_time_utc is redundant since time=0 corresponds to starttime_utc) - self.starttime_utc = starttime_utc - - # Determine the dt implied by the weather file - self.dt_wi = df_wi["time"][1] - df_wi["time"][0] - - # Log the values - if self.verbose: - self.logger.info(f"dt_wi = {self.dt_wi}") - self.logger.info(f"dt = {self.dt}") - - self.logger.info("Interpolating wind input file...") - - # Interpolate df_wi on to the time steps - time_steps_all = np.arange(self.starttime, self.endtime, self.dt) - df_wi = interpolate_df(df_wi, time_steps_all) - - # FLORIS PRECOMPUTATION - - # Initialize the FLORIS model as an ApproxFlorisModel - self.fmodel = ApproxFlorisModel( - self.floris_input_file, - wd_resolution=1.0, - ws_resolution=1.0, - ) - - # Get the layout and number of turbines from FLORIS - self.layout_x = self.fmodel.layout_x - self.layout_y = self.fmodel.layout_y - self.n_turbines = self.fmodel.n_turbines - - self.logger.info("Converting wind input file to numpy matrices...") - - # Convert the wind directions and wind speeds and ti to simply numpy matrices - # Starting with wind speed - # Apply the Hercules float type to the wind speeds - if "ws_mean" in df_wi.columns and "ws_000" not in df_wi.columns: - self.ws_mat = np.tile( - df_wi["ws_mean"].values.astype(hercules_float_type)[:, np.newaxis], - (1, self.n_turbines), - ) - else: - self.ws_mat = df_wi[[f"ws_{t_idx:03d}" for t_idx in range(self.n_turbines)]].to_numpy( - dtype=hercules_float_type - ) - - # Compute the turbine averaged wind speeds (axis = 1) using mean - self.ws_mat_mean = np.mean(self.ws_mat, axis=1, dtype=hercules_float_type) - - self.initial_wind_speeds = self.ws_mat[0, :] - self.wind_speed_mean_background = self.ws_mat_mean[0] - - # For now require "wd_mean" to be in the df_wi - if "wd_mean" not in df_wi.columns: - raise ValueError("Wind input file must contain a column called 'wd_mean'") - self.wd_mat_mean = df_wi["wd_mean"].values.astype(hercules_float_type) - - if "ti_000" in df_wi.columns: - self.ti_mat = df_wi[[f"ti_{t_idx:03d}" for t_idx in range(self.n_turbines)]].to_numpy( - dtype=hercules_float_type - ) - - # Compute the turbine averaged turbulence intensities (axis = 1) using mean - self.ti_mat_mean = np.mean(self.ti_mat, axis=1, dtype=hercules_float_type) - - self.initial_tis = self.ti_mat[0, :] - - else: - self.ti_mat_mean = 0.08 * np.ones_like(self.ws_mat_mean, dtype=hercules_float_type) - - # Precompute the wake deficits at the cadence specified by floris_update_time_s - self.logger.info("Precomputing FLORIS wake deficits...") - - # Determine update step cadence and indices to evaluate FLORIS - update_steps = self.floris_update_steps - n_steps = len(self.ws_mat_mean) - eval_indices = np.arange(update_steps - 1, n_steps, update_steps) - # Ensure at least the final time is evaluated - if eval_indices.size == 0: - eval_indices = np.array([n_steps - 1]) - elif eval_indices[-1] != n_steps - 1: - eval_indices = np.append(eval_indices, n_steps - 1) - - # Build right-aligned windowed means for ws, wd, ti at the evaluation indices - def window_mean(arr_1d, idx, win): - start = max(0, idx - win + 1) - return np.mean(arr_1d[start : idx + 1], dtype=hercules_float_type) - - def window_circmean(arr_1d, idx, win): - start = max(0, idx - win + 1) - return circmean(arr_1d[start : idx + 1], high=360.0, low=0.0, nan_policy="omit") - - ws_eval = np.array( - [window_mean(self.ws_mat_mean, i, update_steps) for i in eval_indices], - dtype=hercules_float_type, - ) - wd_eval = np.array( - [window_circmean(self.wd_mat_mean, i, update_steps) for i in eval_indices], - dtype=hercules_float_type, - ) - if np.isscalar(self.ti_mat_mean): - ti_eval = self.ti_mat_mean * np.ones_like(ws_eval, dtype=hercules_float_type) - else: - ti_eval = np.array( - [window_mean(self.ti_mat_mean, i, update_steps) for i in eval_indices], - dtype=hercules_float_type, - ) - - # Evaluate FLORIS at the evaluation cadence - self.fmodel.set( - wind_directions=wd_eval, - wind_speeds=ws_eval, - turbulence_intensities=ti_eval, - ) - self.logger.info("Running FLORIS...") - self.fmodel.run() - self.num_floris_calcs = 1 - self.logger.info("FLORIS run complete") - - # TODO: THIS CODE WILL WORK IN THE FUTURE - # https://github.com/NREL/floris/pull/1135 - # floris_velocities = ( - # self.fmodel.turbine_average_velocities - # ) # This is a 2D array of shape (len(wind_directions), n_turbines) - - # For now compute in place here (replace later) - expanded_velocities = average_velocity( - velocities=self.fmodel.fmodel_expanded.core.flow_field.u, - method=self.fmodel.fmodel_expanded.core.grid.average_method, - cubature_weights=self.fmodel.fmodel_expanded.core.grid.cubature_weights, - ) - - floris_velocities = map_turbine_powers_uncertain( - unique_turbine_powers=expanded_velocities, - map_to_expanded_inputs=self.fmodel.map_to_expanded_inputs, - weights=self.fmodel.weights, - n_unexpanded=self.fmodel.n_unexpanded, - n_sample_points=self.fmodel.n_sample_points, - n_turbines=self.fmodel.n_turbines, - ).astype(hercules_float_type) - - # Determine the free_stream velocities as the maximum velocity in each row - # of floris velocities. Make sure to keep shape (len(wind_directions), n_turbines) - # by repeating the maximum velocity accross each column for each row - free_stream_velocities = np.tile( - np.max(floris_velocities, axis=1)[:, np.newaxis], (1, self.n_turbines) - ).astype(hercules_float_type) - - # Compute wake deficits at evaluation times - floris_wake_deficits_eval = free_stream_velocities - floris_velocities - - # Expand the wake deficits to all time steps by holding constant within each interval - deficits_all = np.zeros_like(self.ws_mat, dtype=hercules_float_type) - # For each block, fill with the corresponding deficits - prev_end = -1 - for block_idx, end_idx in enumerate(eval_indices): - start_idx = prev_end + 1 - prev_end = end_idx - # Use deficits from this evaluation time for the whole block - deficits_all[start_idx : end_idx + 1, :] = floris_wake_deficits_eval[block_idx, :] - - # Compute all the withwakes wind speeds from background minus deficits - self.wind_speeds_withwakes_all = self.ws_mat - deficits_all - - # Initialize the turbine powers to nan - self.turbine_powers = np.zeros(self.n_turbines, dtype=hercules_float_type) * np.nan - - # Get the initial background wind speeds - self.wind_speeds_background = self.ws_mat[0, :] - - # Compute initial withwakes wind speeds - self.wind_speeds_withwakes = self.wind_speeds_withwakes_all[0, :] - - # Get the initial FLORIS wake deficits - self.floris_wake_deficits = self.wind_speeds_background - self.wind_speeds_withwakes - - # Get the turbine information - self.turbine_dict = load_yaml(self.turbine_file_name) - self.turbine_model_type = self.turbine_dict["turbine_model_type"] - - # Initialize the turbine array - if self.turbine_model_type == "filter_model": - # Use vectorized implementation for improved performance - self.turbine_array = TurbineFilterModelVectorized( - self.turbine_dict, self.dt, self.fmodel, self.wind_speeds_withwakes - ) - self.use_vectorized_turbines = True - elif self.turbine_model_type == "dof1_model": - self.turbine_array = [ - Turbine1dofModel( - self.turbine_dict, self.dt, self.fmodel, self.wind_speeds_withwakes[t_idx] - ) - for t_idx in range(self.n_turbines) - ] - self.use_vectorized_turbines = False - else: - raise Exception("Turbine model type should be either filter_model or dof1_model") - - # Initialize the power array to the initial wind speeds - if self.use_vectorized_turbines: - self.turbine_powers = self.turbine_array.prev_powers.copy() - else: - self.turbine_powers = np.array( - [self.turbine_array[t_idx].prev_power for t_idx in range(self.n_turbines)], - dtype=hercules_float_type, - ) - - # Get the rated power of the turbines, for now assume all turbines have the same rated power - if self.use_vectorized_turbines: - self.rated_turbine_power = self.turbine_array.get_rated_power() - else: - self.rated_turbine_power = self.turbine_array[0].get_rated_power() - - # Get the capacity of the farm - self.capacity = self.n_turbines * self.rated_turbine_power - - # Update the user - self.logger.info( - f"Initialized Wind_MesoToPowerPrecomFloris with {self.n_turbines} turbines" - ) - - def get_initial_conditions_and_meta_data(self, h_dict): - """Add any initial conditions or meta data to the h_dict. - - Meta data is data not explicitly in the input yaml but still useful for other - modules. - - Args: - h_dict (dict): Dictionary containing simulation parameters. - - Returns: - dict: Dictionary containing simulation parameters with initial conditions and meta data. - """ - h_dict["wind_farm"]["n_turbines"] = self.n_turbines - h_dict["wind_farm"]["capacity"] = self.capacity - h_dict["wind_farm"]["rated_turbine_power"] = self.rated_turbine_power - h_dict["wind_farm"]["wind_direction_mean"] = self.wd_mat_mean[0] - h_dict["wind_farm"]["wind_speed_mean_background"] = self.ws_mat_mean[0] - h_dict["wind_farm"]["turbine_powers"] = self.turbine_powers - h_dict["wind_farm"]["power"] = np.sum(self.turbine_powers) - - # Log the start time UTC if available - if hasattr(self, "starttime_utc"): - h_dict["wind_farm"]["starttime_utc"] = self.starttime_utc - - return h_dict - - def step(self, h_dict): - """Execute one simulation step for the wind farm. - - Calculates turbine powers, - and updates the simulation dictionary with results. Handles power_setpoint - signals and optional extra logging outputs. - - Args: - h_dict (dict): Dictionary containing current simulation state including - step number and power_setpoint values for each turbine. - - Returns: - dict: Updated simulation dictionary with wind farm outputs including - turbine powers, total power, and optional extra outputs. - """ - # Get the current step - step = h_dict["step"] - if self.verbose: - self.logger.info(f"step = {step} (of {self.n_steps})") - - # Grab the instantaneous turbine power setpoint signal and update the power_setpoints buffer - turbine_power_setpoints = h_dict[self.component_name]["turbine_power_setpoints"] - - # Update all the wind speeds - self.wind_speeds_background = self.ws_mat[step, :] - self.wind_speeds_withwakes = self.wind_speeds_withwakes_all[step, :] - self.floris_wake_deficits = self.wind_speeds_background - self.wind_speeds_withwakes - - # Update the turbine powers - if self.use_vectorized_turbines: - # Vectorized calculation for all turbines at once - self.turbine_powers = self.turbine_array.step( - self.wind_speeds_withwakes, - turbine_power_setpoints, - ) - else: - # Original loop-based calculation - for t_idx in range(self.n_turbines): - self.turbine_powers[t_idx] = self.turbine_array[t_idx].step( - self.wind_speeds_withwakes[t_idx], - power_setpoint=turbine_power_setpoints[t_idx], - ) - - # Update instantaneous wind direction and wind speed - self.wind_direction_mean = self.wd_mat_mean[step] - self.wind_speed_mean_background = self.ws_mat_mean[step] - - # Update the h_dict with outputs - h_dict[self.component_name]["power"] = np.sum(self.turbine_powers) - h_dict[self.component_name]["turbine_powers"] = self.turbine_powers - h_dict[self.component_name]["turbine_power_setpoints"] = turbine_power_setpoints - h_dict[self.component_name]["wind_direction_mean"] = self.wind_direction_mean - h_dict[self.component_name]["wind_speed_mean_background"] = self.wind_speed_mean_background - h_dict[self.component_name]["wind_speed_mean_withwakes"] = np.mean( - self.wind_speeds_withwakes, dtype=hercules_float_type - ) - h_dict[self.component_name]["wind_speeds_withwakes"] = self.wind_speeds_withwakes - h_dict[self.component_name]["wind_speeds_background"] = self.wind_speeds_background - - return h_dict - - -class TurbineFilterModel: - """Simple filter-based wind turbine model for power output simulation. - - This model simulates wind turbine power output using a first-order filter - to smooth the response to changing wind conditions, providing a simplified - representation of turbine dynamics. - - NOTE: This class is now unused and kept for backward compatibility. - The filter_model turbine_model_type now uses TurbineFilterModelVectorized - for improved performance. - """ - - def __init__(self, turbine_dict, dt, fmodel, initial_wind_speed): - """Initialize the turbine filter model. - - Args: - turbine_dict (dict): Dictionary containing turbine configuration, - including filter model parameters and other turbine-specific data. - dt (float): Time step for the simulation in seconds. - fmodel (FlorisModel): FLORIS model of the farm. - initial_wind_speed (float): Initial wind speed in m/s to initialize - the simulation. - """ - # Save the time step - self.dt = dt - - # Save the turbine dict - self.turbine_dict = turbine_dict - - # Save the filter time constant - self.filter_time_constant = turbine_dict["filter_model"]["time_constant"] - - # Solve for the filter alpha value given dt and the time constant - self.alpha = 1 - np.exp(-self.dt / self.filter_time_constant) - - # Grab the wind speed power curve from the fmodel and define a simple 1D LUT - turbine_type = fmodel.core.farm.turbine_definitions[0] - wind_speeds = turbine_type["power_thrust_table"]["wind_speed"] - powers = turbine_type["power_thrust_table"]["power"] - self.power_lut = interp1d( - wind_speeds, - powers, - fill_value=0.0, - bounds_error=False, - ) - - # Initialize the previous power to the initial wind speed - self.prev_power = self.power_lut(initial_wind_speed) - - def get_rated_power(self): - """Get the rated power of the turbine. - - Returns: - float: The rated power of the turbine in kW. - """ - return np.max(self.power_lut(np.arange(0, 25, 1.0, dtype=hercules_float_type))) - - def step(self, wind_speed, power_setpoint): - """Simulate a single time step of the wind turbine power output. - - This method calculates the power output of a wind turbine based on the - given wind speed and power_setpoint. The power output is - smoothed using an exponential moving average to simulate the turbine's - response to changing wind conditions. - - Args: - wind_speed (float): The current wind speed in meters per second (m/s). - power_setpoint (float): The maximum allowable power output in kW. - - Returns: - float: The calculated power output of the wind turbine, constrained - by the power_setpoint and smoothed using the exponential moving average. - """ - # Instantaneous power - instant_power = self.power_lut(wind_speed) - - # Limit the current power to not be greater than power_setpoint - instant_power = min(instant_power, power_setpoint) - - # Limit the instant power to be greater than 0 - instant_power = max(instant_power, 0.0) - - # TEMP: why are NaNs occurring? - if np.isnan(instant_power): - print( - f"NaN instant power at wind speed {wind_speed} m/s, " - f"power setpoint {power_setpoint} kW, prev power {self.prev_power} kW" - ) - instant_power = self.prev_power - - # Update the power - power = self.alpha * instant_power + (1 - self.alpha) * self.prev_power - - # Limit the power to not be greater than power_setpoint - power = min(power, power_setpoint) - - # Limit the power to be greater than 0 - power = max(power, 0.0) - - # Update the previous power - self.prev_power = power - - # Return the power - return power diff --git a/hercules/utilities.py b/hercules/utilities.py index 6c38223f..a6443644 100644 --- a/hercules/utilities.py +++ b/hercules/utilities.py @@ -50,7 +50,11 @@ def get_available_component_types(): dict: Component names mapped to available simulation types. """ return { - "wind_farm": ["Wind_MesoToPower", "Wind_MesoToPowerPrecomFloris"], + "wind_farm": [ + "Wind_MesoToPower", + "Wind_MesoToPowerPrecomFloris", + "Wind_MesoToPowerNoAddedWakes", + ], "solar_farm": ["SolarPySAMPVWatts"], "battery": ["BatterySimple", "BatteryLithiumIon"], "electrolyzer": ["ElectrolyzerPlant"], diff --git a/tests/test_inputs/scada_input.csv b/tests/test_inputs/scada_input.csv new file mode 100644 index 00000000..b40f4fee --- /dev/null +++ b/tests/test_inputs/scada_input.csv @@ -0,0 +1,13 @@ +time_utc,wd_mean,ws_000,ws_001,ws_002,pow_000,pow_001,pow_002 +2018-05-10 12:31:00,180.5,8.2,8.1,8.3,2500.0,2400.0,2600.0 +2018-05-10 12:31:01,185.2,9.1,9.0,9.2,3200.0,3100.0,3300.0 +2018-05-10 12:31:02,190.8,7.8,7.7,7.9,2200.0,2100.0,2300.0 +2018-05-10 12:31:03,175.3,6.5,6.4,6.6,1500.0,1400.0,1600.0 +2018-05-10 12:31:04,170.1,10.2,10.1,10.3,4200.0,4100.0,4300.0 +2018-05-10 12:31:05,165.7,11.5,11.4,11.6,5000.0,4900.0,5000.0 +2018-05-10 12:31:06,160.4,9.8,9.7,9.9,5000.0,3800.0,4000.0 +2018-05-10 12:31:07,155.9,8.7,8.6,8.8,3000.0,2900.0,3100.0 +2018-05-10 12:31:08,150.2,7.3,7.2,7.4,1900.0,1800.0,2000.0 +2018-05-10 12:31:09,145.6,6.9,6.8,7.0,1700.0,1600.0,1800.0 +2018-05-10 12:31:10,140.3,8.4,8.3,8.5,2700.0,2600.0,2800.0 + diff --git a/tests/wind_farm_direct_test.py b/tests/wind_farm_direct_test.py new file mode 100644 index 00000000..cdd96b26 --- /dev/null +++ b/tests/wind_farm_direct_test.py @@ -0,0 +1,165 @@ +"""Tests for the WindFarm class in direct wake mode (Wind_MesoToPowerNoAddedWakes).""" + +import copy + +import numpy as np +from hercules.plant_components.wind_farm import WindFarm +from hercules.utilities import hercules_float_type + +from tests.test_inputs.h_dict import h_dict_wind + +# Create a base test dictionary for Wind_MesoToPowerNoAddedWakes +h_dict_wind_direct = copy.deepcopy(h_dict_wind) +# Update component type +h_dict_wind_direct["wind_farm"]["component_type"] = "Wind_MesoToPowerNoAddedWakes" + + +def test_wind_farm_direct_initialization(): + """Test that WindFarm initializes correctly with wake_model='no_added_wakes'.""" + wind_sim = WindFarm(h_dict_wind_direct) + + assert wind_sim.component_name == "wind_farm" + assert wind_sim.component_type == "Wind_MesoToPowerNoAddedWakes" + assert wind_sim.wake_model == "no_added_wakes" + assert wind_sim.n_turbines == 3 + assert wind_sim.dt == 1.0 + assert wind_sim.starttime == 0.0 + assert wind_sim.endtime == 10.0 + # No FLORIS calculations in direct mode + assert wind_sim.num_floris_calcs == 0 + assert wind_sim.floris_update_time_s == h_dict_wind_direct["wind_farm"]["floris_update_time_s"] + + +def test_wind_farm_direct_no_wakes(): + """Test that no wake deficits are applied in direct mode.""" + wind_sim = WindFarm(h_dict_wind_direct) + + # Verify initial wake deficits are zero + assert np.all(wind_sim.floris_wake_deficits == 0.0) + + # Verify initial wind speeds with wakes equal background + assert np.allclose(wind_sim.wind_speeds_withwakes, wind_sim.wind_speeds_background) + + +def test_wind_farm_direct_step(): + """Test that the step method works correctly in direct mode.""" + wind_sim = WindFarm(h_dict_wind_direct) + + # Add power setpoint values to the step h_dict + step_h_dict = {"step": 1} + step_h_dict["wind_farm"] = { + "turbine_power_setpoints": np.array([1000.0, 1500.0, 2000.0]), + } + + result = wind_sim.step(step_h_dict) + + # Verify outputs exist + assert "turbine_powers" in result["wind_farm"] + assert "power" in result["wind_farm"] + assert len(result["wind_farm"]["turbine_powers"]) == 3 + assert isinstance(result["wind_farm"]["turbine_powers"], np.ndarray) + assert "power" in result["wind_farm"] + assert isinstance(result["wind_farm"]["power"], (int, float)) + + # Verify no wake deficits applied + assert np.all(wind_sim.floris_wake_deficits == 0.0) + assert np.allclose( + result["wind_farm"]["wind_speeds_withwakes"], + result["wind_farm"]["wind_speeds_background"], + ) + + +def test_wind_farm_direct_no_wake_deficits_over_time(): + """Test that wake deficits remain zero throughout simulation.""" + wind_sim = WindFarm(h_dict_wind_direct) + + # Run multiple steps + for step in range(5): + step_h_dict = {"step": step} + step_h_dict["wind_farm"] = { + "turbine_power_setpoints": np.ones(3, dtype=hercules_float_type) * 5000.0, + } + + result = wind_sim.step(step_h_dict) + + # Verify no wakes at each step + assert np.all(wind_sim.floris_wake_deficits == 0.0) + assert np.allclose( + result["wind_farm"]["wind_speeds_withwakes"], + result["wind_farm"]["wind_speeds_background"], + ) + + +def test_wind_farm_direct_turbine_dynamics(): + """Test that turbine dynamics still work in direct mode.""" + wind_sim = WindFarm(h_dict_wind_direct) + + # Run a step with very low power setpoint + step_h_dict = {"step": 1} + step_h_dict["wind_farm"] = { + "turbine_power_setpoints": np.array([100.0, 100.0, 100.0]), + } + + result = wind_sim.step(step_h_dict) + + # Turbine powers should be limited by setpoint + assert np.all(result["wind_farm"]["turbine_powers"] <= 100.0 + 1e-6) + + +def test_wind_farm_direct_power_setpoint_zero(): + """Test that turbine powers go to zero when setpoint is zero.""" + wind_sim = WindFarm(h_dict_wind_direct) + + # Run multiple steps with zero setpoint to ensure filter settles + for step in range(10): + step_h_dict = {"step": step} + step_h_dict["wind_farm"] = { + "turbine_power_setpoints": np.zeros(3, dtype=hercules_float_type), + } + result = wind_sim.step(step_h_dict) + + # After multiple steps, powers should be effectively zero + assert np.all(result["wind_farm"]["turbine_powers"] < 1.0) + + +def test_wind_farm_direct_initial_conditions(): + """Test that initial conditions are correctly set in h_dict.""" + wind_sim = WindFarm(h_dict_wind_direct) + + initial_h_dict = copy.deepcopy(h_dict_wind_direct) + result_h_dict = wind_sim.get_initial_conditions_and_meta_data(initial_h_dict) + + assert "n_turbines" in result_h_dict["wind_farm"] + assert "capacity" in result_h_dict["wind_farm"] + assert "rated_turbine_power" in result_h_dict["wind_farm"] + assert "wind_direction_mean" in result_h_dict["wind_farm"] + assert "wind_speed_mean_background" in result_h_dict["wind_farm"] + assert "turbine_powers" in result_h_dict["wind_farm"] + assert "power" in result_h_dict["wind_farm"] + + assert result_h_dict["wind_farm"]["n_turbines"] == 3 + assert result_h_dict["wind_farm"]["capacity"] > 0 + assert result_h_dict["wind_farm"]["rated_turbine_power"] > 0 + + +def test_wind_farm_direct_output_consistency(): + """Test that outputs are consistent with no wake modeling.""" + wind_sim = WindFarm(h_dict_wind_direct) + + # Run a step + step_h_dict = {"step": 2} + step_h_dict["wind_farm"] = { + "turbine_power_setpoints": np.ones(3, dtype=hercules_float_type) * 5000.0, + } + + result = wind_sim.step(step_h_dict) + + # Calculate expected mean withwakes speed (should equal background mean) + expected_mean_withwakes = np.mean( + result["wind_farm"]["wind_speeds_background"], dtype=hercules_float_type + ) + + assert np.isclose(result["wind_farm"]["wind_speed_mean_withwakes"], expected_mean_withwakes) + + # Total power should be sum of turbine powers + assert np.isclose(result["wind_farm"]["power"], np.sum(result["wind_farm"]["turbine_powers"])) diff --git a/tests/wind_farm_scada_power_test.py b/tests/wind_farm_scada_power_test.py new file mode 100644 index 00000000..b355068e --- /dev/null +++ b/tests/wind_farm_scada_power_test.py @@ -0,0 +1,450 @@ +"""Tests for the WindFarmSCADAPower class.""" + +import copy +import os +import tempfile + +import numpy as np +import pandas as pd +import pytest +from hercules.plant_components.wind_farm_scada_power import WindFarmSCADAPower +from hercules.utilities import hercules_float_type + +from tests.test_inputs.h_dict import h_dict_wind + +# Create a base test dictionary for WindFarmSCADAPower +h_dict_wind_scada = copy.deepcopy(h_dict_wind) +# Update component type and remove unneeded parameters +h_dict_wind_scada["wind_farm"]["component_type"] = "WindFarmSCADAPower" +h_dict_wind_scada["wind_farm"]["scada_filename"] = "tests/test_inputs/scada_input.csv" +# Keep turbine_file_name for filter model parameters +# Remove FLORIS-specific parameters +del h_dict_wind_scada["wind_farm"]["floris_input_file"] +del h_dict_wind_scada["wind_farm"]["wind_input_filename"] +del h_dict_wind_scada["wind_farm"]["floris_update_time_s"] + + +def test_wind_farm_scada_power_initialization(): + """Test that WindFarmSCADAPower initializes correctly with valid inputs.""" + wind_sim = WindFarmSCADAPower(h_dict_wind_scada) + + assert wind_sim.component_name == "wind_farm" + assert wind_sim.component_type == "WindFarmSCADAPower" + assert wind_sim.n_turbines == 3 + assert wind_sim.dt == 1.0 + assert wind_sim.starttime == 0.0 + assert wind_sim.endtime == 10.0 + # No FLORIS calculations in SCADA power mode + assert wind_sim.num_floris_calcs == 0 + + +def test_wind_farm_scada_power_infers_n_turbines(): + """Test that number of turbines is correctly inferred from power columns.""" + wind_sim = WindFarmSCADAPower(h_dict_wind_scada) + + assert wind_sim.n_turbines == 3 + assert len(wind_sim.power_columns) == 3 + assert wind_sim.power_columns == ["pow_000", "pow_001", "pow_002"] + + +def test_wind_farm_scada_power_infers_rated_power(): + """Test that rated power is correctly inferred from 99th percentile.""" + wind_sim = WindFarmSCADAPower(h_dict_wind_scada) + + # Check that rated power is positive and reasonable + assert wind_sim.rated_turbine_power == 5000.0 + assert wind_sim.capacity == wind_sim.n_turbines * wind_sim.rated_turbine_power + + +def test_wind_farm_scada_power_no_wakes(): + """Test that no wake deficits are applied in SCADA power mode.""" + wind_sim = WindFarmSCADAPower(h_dict_wind_scada) + + # Verify initial wake deficits are zero + assert np.all(wind_sim.floris_wake_deficits == 0.0) + + # Verify initial wind speeds with wakes equal background + assert np.allclose(wind_sim.wind_speeds_withwakes, wind_sim.wind_speeds_background) + + +def test_wind_farm_scada_power_step(): + """Test that the step method works correctly.""" + wind_sim = WindFarmSCADAPower(h_dict_wind_scada) + + # Add power setpoint values to the step h_dict + step_h_dict = {"step": 1} + step_h_dict["wind_farm"] = { + "turbine_power_setpoints": np.array([5000.0, 5000.0, 5000.0]), + } + + result = wind_sim.step(step_h_dict) + + # Verify outputs exist + assert "turbine_powers" in result["wind_farm"] + assert "power" in result["wind_farm"] + assert len(result["wind_farm"]["turbine_powers"]) == 3 + assert isinstance(result["wind_farm"]["turbine_powers"], np.ndarray) + assert "power" in result["wind_farm"] + assert isinstance(result["wind_farm"]["power"], (int, float)) + + # Verify no wake deficits applied + assert np.all(wind_sim.floris_wake_deficits == 0.0) + assert np.allclose( + result["wind_farm"]["wind_speeds_withwakes"], + result["wind_farm"]["wind_speeds_background"], + ) + + # Verify turbine powers + assert np.allclose(result["wind_farm"]["turbine_powers"], [3200.0, 3100.0, 3300.0]) + assert np.isclose(result["wind_farm"]["power"], 3200.0 + 3100.0 + 3300.0) + + +def test_wind_farm_scada_power_power_setpoint_applies(): + """Test that turbine powers are limited by power setpoint when setpoint is low.""" + wind_sim = WindFarmSCADAPower(h_dict_wind_scada) + + # Set very low power setpoint values that should definitely limit power output + # Run multiple steps to let filter settle (within available data range 0-9) + for step in range(wind_sim.n_steps): + step_h_dict = {"step": step} + step_h_dict["wind_farm"] = { + "turbine_power_setpoints": np.array([100.0, 200.0, 300.0]), # Very low setpoints + } + result = wind_sim.step(step_h_dict) + + # Verify that turbine powers are at or below power setpoint limits + turbine_powers = result["wind_farm"]["turbine_powers"] + power_setpoints = [100.0, 200.0, 300.0] + + for i, (power, setpoint) in enumerate(zip(turbine_powers, power_setpoints)): + assert ( + power <= setpoint + 1e-6 + ), f"Turbine {i} power {power} exceeds power setpoint {setpoint}" + + +def test_wind_farm_scada_power_power_setpoint_zero(): + """Test that turbine powers go to zero when setpoint is zero.""" + wind_sim = WindFarmSCADAPower(h_dict_wind_scada) + + # Run multiple steps with zero setpoint to ensure filter settles (within available data range) + for step in range(wind_sim.n_steps): + step_h_dict = {"step": step} + step_h_dict["wind_farm"] = { + "turbine_power_setpoints": np.zeros(3, dtype=hercules_float_type), + } + result = wind_sim.step(step_h_dict) + + # After multiple steps, powers should be effectively zero + assert np.all(result["wind_farm"]["turbine_powers"] < 1.0) + + +def test_wind_farm_scada_power_get_initial_conditions_and_meta_data(): + """Test that get_initial_conditions_and_meta_data adds correct metadata to h_dict.""" + wind_sim = WindFarmSCADAPower(h_dict_wind_scada) + + # Create a copy of the input h_dict to avoid modifying the original + test_h_dict_copy = copy.deepcopy(h_dict_wind_scada) + + # Call the method + result = wind_sim.get_initial_conditions_and_meta_data(test_h_dict_copy) + + # Verify that the method returns the modified h_dict + assert result is test_h_dict_copy + + # Verify that all expected metadata is added to the wind_farm section + assert "n_turbines" in result["wind_farm"] + assert "capacity" in result["wind_farm"] + assert "rated_turbine_power" in result["wind_farm"] + assert "wind_direction_mean" in result["wind_farm"] + assert "wind_speed_mean_background" in result["wind_farm"] + assert "turbine_powers" in result["wind_farm"] + + # Verify the values match the wind_sim attributes + assert result["wind_farm"]["n_turbines"] == wind_sim.n_turbines + assert result["wind_farm"]["capacity"] == wind_sim.capacity + assert result["wind_farm"]["rated_turbine_power"] == wind_sim.rated_turbine_power + assert result["wind_farm"]["wind_direction_mean"] == wind_sim.wd_mat_mean[0] + assert result["wind_farm"]["wind_speed_mean_background"] == wind_sim.ws_mat_mean[0] + + # Verify turbine_powers is a numpy array with correct length + assert isinstance(result["wind_farm"]["turbine_powers"], np.ndarray) + assert len(result["wind_farm"]["turbine_powers"]) == wind_sim.n_turbines + np.testing.assert_array_equal(result["wind_farm"]["turbine_powers"], wind_sim.turbine_powers) + + # Verify that the original h_dict structure is preserved + assert "dt" in result + assert "starttime" in result + assert "endtime" in result + assert "plant" in result + + +def test_wind_farm_scada_power_time_utc_handling(): + """Test that time_utc is correctly parsed and validated.""" + # Create wind input data with time_utc columns + scada_data = { + "time_utc": [ + "2023-01-01T00:00:00Z", + "2023-01-01T00:00:01Z", + "2023-01-01T00:00:02Z", + "2023-01-01T00:00:03Z", + "2023-01-01T00:00:04Z", + ], + "wd_mean": [270.0, 275.0, 280.0, 285.0, 290.0], + "ws_000": [8.0, 9.0, 10.0, 11.0, 12.0], + "ws_001": [8.5, 9.5, 10.5, 11.5, 12.5], + "ws_002": [9.0, 10.0, 11.0, 12.0, 13.0], + "pow_000": [2500.0, 3200.0, 4000.0, 4500.0, 5000.0], + "pow_001": [2400.0, 3100.0, 3900.0, 4400.0, 4900.0], + "pow_002": [2600.0, 3300.0, 4100.0, 4600.0, 5000.0], + } + + # Create temporary file + with tempfile.NamedTemporaryFile(mode="w", suffix=".csv", delete=False) as f: + df = pd.DataFrame(scada_data) + df.to_csv(f.name, index=False) + temp_scada_file = f.name + + try: + # Create test h_dict with the temporary scada file + test_h_dict = copy.deepcopy(h_dict_wind_scada) + test_h_dict["wind_farm"]["scada_filename"] = temp_scada_file + test_h_dict["starttime"] = 0.0 + test_h_dict["endtime"] = 4.0 + test_h_dict["starttime_utc"] = "2023-01-01T00:00:00Z" + test_h_dict["endtime_utc"] = "2023-01-01T00:00:04Z" + test_h_dict["dt"] = 1.0 + + # Initialize wind simulation + wind_sim = WindFarmSCADAPower(test_h_dict) + + # Verify that starttime_utc is set correctly + assert hasattr(wind_sim, "starttime_utc"), "starttime_utc should be set" + + expected_start_time = pd.to_datetime("2023-01-01T00:00:00Z", utc=True) + + # Convert to pandas Timestamp for comparison + actual_start_time = pd.Timestamp(wind_sim.starttime_utc) + + # Compare datetime values + assert actual_start_time.replace(tzinfo=None) == expected_start_time.replace( + tzinfo=None + ), f"starttime_utc mismatch: expected {expected_start_time}, got {actual_start_time}" + + finally: + # Clean up temporary file + if os.path.exists(temp_scada_file): + os.unlink(temp_scada_file) + + +def test_wind_farm_scada_power_time_utc_validation_start_too_early(): + """Test that error is raised when starttime_utc is before earliest SCADA data.""" + # Create SCADA data starting at 2023-01-01T00:00:00Z + scada_data = { + "time_utc": [ + "2023-01-01T00:00:00Z", + "2023-01-01T00:00:01Z", + "2023-01-01T00:00:02Z", + ], + "wd_mean": [270.0, 275.0, 280.0], + "ws_000": [8.0, 9.0, 10.0], + "ws_001": [8.5, 9.5, 10.5], + "ws_002": [9.0, 10.0, 11.0], + "pow_000": [2500.0, 3200.0, 4000.0], + "pow_001": [2400.0, 3100.0, 3900.0], + "pow_002": [2600.0, 3300.0, 4100.0], + } + + with tempfile.NamedTemporaryFile(mode="w", suffix=".csv", delete=False) as f: + df = pd.DataFrame(scada_data) + df.to_csv(f.name, index=False) + temp_scada_file = f.name + + try: + test_h_dict = copy.deepcopy(h_dict_wind_scada) + test_h_dict["wind_farm"]["scada_filename"] = temp_scada_file + test_h_dict["starttime"] = 0.0 + test_h_dict["endtime"] = 2.0 + # Try to start before earliest SCADA data + test_h_dict["starttime_utc"] = "2022-12-31T23:59:59Z" + test_h_dict["endtime_utc"] = "2023-01-01T00:00:02Z" + test_h_dict["dt"] = 1.0 + + with pytest.raises(ValueError, match="Start time UTC .* is before the earliest time"): + WindFarmSCADAPower(test_h_dict) + + finally: + if os.path.exists(temp_scada_file): + os.unlink(temp_scada_file) + + +def test_wind_farm_scada_power_time_utc_validation_end_too_late(): + """Test that error is raised when endtime_utc is after latest SCADA data.""" + # Create SCADA data ending at 2023-01-01T00:00:02Z + scada_data = { + "time_utc": [ + "2023-01-01T00:00:00Z", + "2023-01-01T00:00:01Z", + "2023-01-01T00:00:02Z", + ], + "wd_mean": [270.0, 275.0, 280.0], + "ws_000": [8.0, 9.0, 10.0], + "ws_001": [8.5, 9.5, 10.5], + "ws_002": [9.0, 10.0, 11.0], + "pow_000": [2500.0, 3200.0, 4000.0], + "pow_001": [2400.0, 3100.0, 3900.0], + "pow_002": [2600.0, 3300.0, 4100.0], + } + + with tempfile.NamedTemporaryFile(mode="w", suffix=".csv", delete=False) as f: + df = pd.DataFrame(scada_data) + df.to_csv(f.name, index=False) + temp_scada_file = f.name + + try: + test_h_dict = copy.deepcopy(h_dict_wind_scada) + test_h_dict["wind_farm"]["scada_filename"] = temp_scada_file + test_h_dict["starttime"] = 0.0 + test_h_dict["endtime"] = 5.0 + test_h_dict["starttime_utc"] = "2023-01-01T00:00:00Z" + # Try to end after latest SCADA data + test_h_dict["endtime_utc"] = "2023-01-01T00:00:05Z" + test_h_dict["dt"] = 1.0 + + with pytest.raises(ValueError, match="End time UTC .* is after the latest time"): + WindFarmSCADAPower(test_h_dict) + + finally: + if os.path.exists(temp_scada_file): + os.unlink(temp_scada_file) + + +def test_wind_farm_scada_power_ws_mean_handling(): + """Test that ws_mean is correctly handled when individual speeds are not present.""" + # Create SCADA data with ws_mean but no individual speeds + scada_data = { + "time_utc": [ + "2023-01-01T00:00:00Z", + "2023-01-01T00:00:01Z", + "2023-01-01T00:00:02Z", + "2023-01-01T00:00:03Z", + "2023-01-01T00:00:04Z", + ], + "wd_mean": [270.0, 275.0, 280.0, 285.0, 290.0], + "ws_mean": [8.0, 9.0, 10.0, 11.0, 12.0], + "pow_000": [2500.0, 3200.0, 4000.0, 4500.0, 5000.0], + "pow_001": [2400.0, 3100.0, 3900.0, 4400.0, 4900.0], + "pow_002": [2600.0, 3300.0, 4100.0, 4600.0, 5000.0], + } + + with tempfile.NamedTemporaryFile(mode="w", suffix=".csv", delete=False) as f: + df = pd.DataFrame(scada_data) + df.to_csv(f.name, index=False) + temp_scada_file = f.name + + try: + test_h_dict = copy.deepcopy(h_dict_wind_scada) + test_h_dict["wind_farm"]["scada_filename"] = temp_scada_file + test_h_dict["starttime"] = 0.0 + test_h_dict["endtime"] = 4.0 + test_h_dict["starttime_utc"] = "2023-01-01T00:00:00Z" + test_h_dict["endtime_utc"] = "2023-01-01T00:00:04Z" + test_h_dict["dt"] = 1.0 + + wind_sim = WindFarmSCADAPower(test_h_dict) + + # Verify that ws_mat is properly tiled from ws_mean + assert wind_sim.ws_mat.shape == (4, 3) + # All turbines should have the same wind speed (from ws_mean) + assert (wind_sim.ws_mat[:, 0] == wind_sim.ws_mat[:, 1]).all() + assert (wind_sim.ws_mat[:, 1] == wind_sim.ws_mat[:, 2]).all() + + finally: + if os.path.exists(temp_scada_file): + os.unlink(temp_scada_file) + + +def test_wind_farm_scada_power_output_consistency(): + """Test that outputs are consistent with no wake modeling.""" + wind_sim = WindFarmSCADAPower(h_dict_wind_scada) + + # Run a step + step_h_dict = {"step": 2} + step_h_dict["wind_farm"] = { + "turbine_power_setpoints": np.ones(3, dtype=hercules_float_type) * 5000.0, + } + + result = wind_sim.step(step_h_dict) + + # Calculate expected mean withwakes speed (should equal background mean) + expected_mean_withwakes = np.mean( + result["wind_farm"]["wind_speeds_background"], dtype=hercules_float_type + ) + + assert np.isclose(result["wind_farm"]["wind_speed_mean_withwakes"], expected_mean_withwakes) + + # Total power should be sum of turbine powers + assert np.isclose(result["wind_farm"]["power"], np.sum(result["wind_farm"]["turbine_powers"])) + + +def test_wind_farm_scada_power_multiple_file_formats(): + """Test that SCADA data can be loaded from different file formats.""" + # Test CSV (already tested above, but included for completeness) + wind_sim_csv = WindFarmSCADAPower(h_dict_wind_scada) + assert wind_sim_csv.n_turbines == 3 + + # Test pickle format + current_dir = os.path.dirname(__file__) + df_scada = pd.read_csv(current_dir + "/test_inputs/scada_input.csv") + + # Create temporary pickle file + with tempfile.NamedTemporaryFile(suffix=".pkl", delete=False) as f: + df_scada.to_pickle(f.name) + temp_pickle_file = f.name + + try: + test_h_dict = copy.deepcopy(h_dict_wind_scada) + test_h_dict["wind_farm"]["scada_filename"] = temp_pickle_file + + wind_sim_pkl = WindFarmSCADAPower(test_h_dict) + assert wind_sim_pkl.n_turbines == 3 + + finally: + if os.path.exists(temp_pickle_file): + os.unlink(temp_pickle_file) + + # Test feather format + with tempfile.NamedTemporaryFile(suffix=".ftr", delete=False) as f: + df_scada.to_feather(f.name) + temp_feather_file = f.name + + try: + test_h_dict = copy.deepcopy(h_dict_wind_scada) + test_h_dict["wind_farm"]["scada_filename"] = temp_feather_file + + wind_sim_ftr = WindFarmSCADAPower(test_h_dict) + assert wind_sim_ftr.n_turbines == 3 + + finally: + if os.path.exists(temp_feather_file): + os.unlink(temp_feather_file) + + +def test_wind_farm_scada_power_invalid_file_format(): + """Test that invalid file format raises ValueError.""" + test_h_dict = copy.deepcopy(h_dict_wind_scada) + test_h_dict["wind_farm"]["scada_filename"] = "tests/test_inputs/invalid.txt" + + # Create a dummy file + with tempfile.NamedTemporaryFile(mode="w", suffix=".txt", delete=False) as f: + f.write("dummy") + temp_file = f.name + + try: + test_h_dict["wind_farm"]["scada_filename"] = temp_file + + with pytest.raises(ValueError, match="SCADA file must be a .csv or .p, .f or .ftr file"): + WindFarmSCADAPower(test_h_dict) + + finally: + if os.path.exists(temp_file): + os.unlink(temp_file) diff --git a/tests/wind_meso_to_power_precom_floris_test.py b/tests/wind_meso_to_power_precom_floris_test.py index 2590b6bb..9f15f725 100644 --- a/tests/wind_meso_to_power_precom_floris_test.py +++ b/tests/wind_meso_to_power_precom_floris_test.py @@ -1,4 +1,4 @@ -"""Tests for the Wind_MesoToPowerPrecomFloris class.""" +"""Tests for the WindFarm class in precomputed wake mode (Wind_MesoToPowerPrecomFloris).""" import copy import os @@ -7,9 +7,7 @@ import numpy as np import pandas as pd import pytest -from hercules.plant_components.wind_meso_to_power_precom_floris import ( - Wind_MesoToPowerPrecomFloris, -) +from hercules.plant_components.wind_farm import WindFarm from hercules.utilities import hercules_float_type from tests.test_inputs.h_dict import h_dict_wind @@ -22,7 +20,7 @@ def test_wind_meso_to_power_precom_floris_initialization(): """Test that Wind_MesoToPowerPrecomFloris initializes correctly with valid inputs.""" - wind_sim = Wind_MesoToPowerPrecomFloris(h_dict_wind_precom_floris) + wind_sim = WindFarm(h_dict_wind_precom_floris) assert wind_sim.component_name == "wind_farm" assert wind_sim.component_type == "Wind_MesoToPowerPrecomFloris" @@ -52,7 +50,7 @@ def test_wind_meso_to_power_precom_floris_ws_mean(): # Test that, since individual speed are specified, ws_mean is ignored # Note that h_dict_wind_precom_floris specifies an end time of 10. - wind_sim = Wind_MesoToPowerPrecomFloris(test_h_dict) + wind_sim = WindFarm(test_h_dict) assert ( wind_sim.ws_mat[:, 0] == df_input["ws_000"].to_numpy(dtype=hercules_float_type)[:10] ).all() @@ -67,7 +65,7 @@ def test_wind_meso_to_power_precom_floris_ws_mean(): df_input = df_input.drop(columns=["ws_000", "ws_001", "ws_002"]) df_input.to_csv(current_dir + "/test_inputs/wind_input_temp.csv") - wind_sim = Wind_MesoToPowerPrecomFloris(test_h_dict) + wind_sim = WindFarm(test_h_dict) assert (wind_sim.ws_mat_mean == 10.0).all() assert (wind_sim.ws_mat[:, :] == 10.0).all() @@ -81,7 +79,7 @@ def test_wind_meso_to_power_precom_floris_requires_floris_update_time(): del test_h_dict["wind_farm"]["floris_update_time_s"] with pytest.raises(ValueError, match="floris_update_time_s must be in the h_dict"): - Wind_MesoToPowerPrecomFloris(test_h_dict) + WindFarm(test_h_dict) def test_wind_meso_to_power_precom_floris_invalid_update_time(): @@ -90,12 +88,12 @@ def test_wind_meso_to_power_precom_floris_invalid_update_time(): test_h_dict["wind_farm"]["floris_update_time_s"] = 0.5 with pytest.raises(ValueError, match="FLORIS update time must be at least 1 second"): - Wind_MesoToPowerPrecomFloris(test_h_dict) + WindFarm(test_h_dict) def test_wind_meso_to_power_precom_floris_step(): """Test that the step method updates outputs correctly.""" - wind_sim = Wind_MesoToPowerPrecomFloris(h_dict_wind_precom_floris) + wind_sim = WindFarm(h_dict_wind_precom_floris) # Add power setpoint values to the step h_dict step_h_dict = {"step": 1} @@ -115,7 +113,7 @@ def test_wind_meso_to_power_precom_floris_step(): def test_wind_meso_to_power_precom_floris_power_setpoint_applies(): """Test that turbine powers equal power setpoint when setpoint is very low.""" - wind_sim = Wind_MesoToPowerPrecomFloris(h_dict_wind_precom_floris) + wind_sim = WindFarm(h_dict_wind_precom_floris) # Set very low power setpoint values that should definitely limit power output step_h_dict = {"step": 1} @@ -137,7 +135,7 @@ def test_wind_meso_to_power_precom_floris_power_setpoint_applies(): def test_wind_meso_to_power_precom_floris_get_initial_conditions_and_meta_data(): """Test that get_initial_conditions_and_meta_data adds correct metadata to h_dict.""" - wind_sim = Wind_MesoToPowerPrecomFloris(h_dict_wind_precom_floris) + wind_sim = WindFarm(h_dict_wind_precom_floris) # Create a copy of the input h_dict to avoid modifying the original test_h_dict_copy = copy.deepcopy(h_dict_wind_precom_floris) @@ -177,7 +175,7 @@ def test_wind_meso_to_power_precom_floris_get_initial_conditions_and_meta_data() def test_wind_meso_to_power_precom_floris_precomputed_wake_deficits(): """Test that wake deficits are precomputed and stored correctly.""" - wind_sim = Wind_MesoToPowerPrecomFloris(h_dict_wind_precom_floris) + wind_sim = WindFarm(h_dict_wind_precom_floris) # Verify that precomputed wake wind speeds exist assert hasattr(wind_sim, "wind_speeds_withwakes_all") @@ -231,7 +229,7 @@ def test_wind_meso_to_power_precom_floris_velocities_update_correctly(): test_h_dict["dt"] = 1.0 # Initialize wind simulation - wind_sim = Wind_MesoToPowerPrecomFloris(test_h_dict) + wind_sim = WindFarm(test_h_dict) # Store initial wind speeds initial_background = wind_sim.wind_speeds_background.copy() @@ -303,7 +301,7 @@ def test_wind_meso_to_power_precom_floris_time_utc_reconstruction(): test_h_dict["dt"] = 1.0 # Initialize wind simulation - wind_sim = Wind_MesoToPowerPrecomFloris(test_h_dict) + wind_sim = WindFarm(test_h_dict) # Verify that starttime_utc is set correctly assert hasattr(wind_sim, "starttime_utc"), "starttime_utc should be set" @@ -440,7 +438,7 @@ def test_wind_meso_to_power_precom_floris_time_utc_different_starttime(): test_h_dict["dt"] = 1.0 # Initialize wind simulation - wind_sim = Wind_MesoToPowerPrecomFloris(test_h_dict) + wind_sim = WindFarm(test_h_dict) # Verify that starttime_utc is set correctly assert hasattr(wind_sim, "starttime_utc"), "starttime_utc should be set" diff --git a/tests/wind_meso_to_power_test.py b/tests/wind_meso_to_power_test.py index 89481e10..58c768a4 100644 --- a/tests/wind_meso_to_power_test.py +++ b/tests/wind_meso_to_power_test.py @@ -1,4 +1,4 @@ -"""Tests for the Wind_MesoToPower class.""" +"""Tests for the WindFarm class in dynamic wake mode (Wind_MesoToPower).""" import copy import os @@ -7,15 +7,15 @@ import numpy as np import pandas as pd import pytest -from hercules.plant_components.wind_meso_to_power import TurbineFilterModel, Wind_MesoToPower +from hercules.plant_components.wind_farm import WindFarm from hercules.utilities import hercules_float_type from tests.test_inputs.h_dict import h_dict_wind def test_wind_meso_to_power_initialization(): - """Test that Wind_MesoToPower initializes correctly with valid inputs.""" - wind_sim = Wind_MesoToPower(h_dict_wind) + """Test that WindFarm initializes correctly with valid inputs (dynamic mode).""" + wind_sim = WindFarm(h_dict_wind) assert wind_sim.component_name == "wind_farm" assert wind_sim.component_type == "Wind_MesoToPower" @@ -41,7 +41,7 @@ def test_wind_meso_to_power_precom_floris_ws_mean(): # Test that, since individual speed are specified, ws_mean is ignored # Note that h_dict_wind specifies an end time of 10. - wind_sim = Wind_MesoToPower(test_h_dict) + wind_sim = WindFarm(test_h_dict) assert ( wind_sim.ws_mat[:, 0] == df_input["ws_000"].to_numpy(dtype=hercules_float_type)[:10] ).all() @@ -56,7 +56,7 @@ def test_wind_meso_to_power_precom_floris_ws_mean(): df_input = df_input.drop(columns=["ws_000", "ws_001", "ws_002"]) df_input.to_csv(current_dir + "/test_inputs/wind_input_temp.csv") - wind_sim = Wind_MesoToPower(test_h_dict) + wind_sim = WindFarm(test_h_dict) assert (wind_sim.ws_mat_mean == 10.0).all() assert (wind_sim.ws_mat[:, :] == 10.0).all() @@ -70,7 +70,7 @@ def test_wind_meso_to_power_missing_floris_update_time(): del test_h_dict["wind_farm"]["floris_update_time_s"] with pytest.raises(ValueError, match="floris_update_time_s must be in the h_dict"): - Wind_MesoToPower(test_h_dict) + WindFarm(test_h_dict) def test_wind_meso_to_power_invalid_update_time(): @@ -79,7 +79,7 @@ def test_wind_meso_to_power_invalid_update_time(): test_h_dict["wind_farm"]["floris_update_time_s"] = 0.5 # Less than 1 second with pytest.raises(ValueError, match="FLORIS update time must be at least 1 second"): - Wind_MesoToPower(test_h_dict) + WindFarm(test_h_dict) def test_wind_meso_to_power_step(): @@ -88,7 +88,7 @@ def test_wind_meso_to_power_step(): # Set a shorter update time for testing test_h_dict["wind_farm"]["floris_update_time_s"] = 1.0 - wind_sim = Wind_MesoToPower(test_h_dict) + wind_sim = WindFarm(test_h_dict) # Add power setpoint values to the step h_dict step_h_dict = {"step": 1} @@ -106,66 +106,9 @@ def test_wind_meso_to_power_step(): assert isinstance(result["wind_farm"]["power"], (int, float)) -def test_turbine_filter_model_initialization(): - """Test that TurbineFilterModel initializes correctly.""" - from floris import FlorisModel - - turbine_dict = {"filter_model": {"time_constant": 12.0}} - - # Use actual FLORIS model - fmodel = FlorisModel("tests/test_inputs/floris_input.yaml") - - turbine = TurbineFilterModel(turbine_dict, dt=1.0, fmodel=fmodel, initial_wind_speed=8.0) - - assert turbine.dt == 1.0 - assert turbine.filter_time_constant == 12.0 - assert turbine.alpha > 0.0 - assert turbine.alpha < 1.0 - assert isinstance(turbine.prev_power, (int, float, np.ndarray)) - - -def test_turbine_filter_model_step(): - """Test that TurbineFilterModel step method works correctly.""" - from floris import FlorisModel - - turbine_dict = {"filter_model": {"time_constant": 12.0}} - - # Use actual FLORIS model - fmodel = FlorisModel("tests/test_inputs/floris_input.yaml") - - turbine = TurbineFilterModel(turbine_dict, dt=1.0, fmodel=fmodel, initial_wind_speed=8.0) - - # Test step with different wind speeds - power1 = turbine.step(wind_speed=10.0, power_setpoint=1000.0) - power2 = turbine.step(wind_speed=12.0, power_setpoint=1500.0) - - assert isinstance(power1, (int, float)) - assert isinstance(power2, (int, float)) - assert power1 >= 0.0 - assert power2 >= 0.0 - - -def test_turbine_filter_model_power_setpoint_limit(): - """Test that TurbineFilterModel respects power setpoint limits.""" - from floris import FlorisModel - - turbine_dict = {"filter_model": {"time_constant": 12.0}} - - # Use actual FLORIS model - fmodel = FlorisModel("tests/test_inputs/floris_input.yaml") - - turbine = TurbineFilterModel(turbine_dict, dt=1.0, fmodel=fmodel, initial_wind_speed=8.0) - - # Test with low power setpoint limit - power = turbine.step(wind_speed=15.0, power_setpoint=500.0) - - assert power <= 500.0 - assert power >= 0.0 - - def test_wind_meso_to_power_time_utc_conversion(): """Test that time_utc column is properly converted to datetime.""" - wind_sim = Wind_MesoToPower(h_dict_wind) + wind_sim = WindFarm(h_dict_wind) # Check that time_utc was converted to datetime type # The wind_sim should have successfully processed the CSV with time_utc column @@ -184,7 +127,7 @@ def test_wind_meso_to_power_power_setpoint_too_high(): test_h_dict = copy.deepcopy(h_dict_wind) test_h_dict["wind_farm"]["floris_update_time_s"] = 1.0 - wind_sim = Wind_MesoToPower(test_h_dict) + wind_sim = WindFarm(test_h_dict) # Set very high power setpoint values that should not limit power output step_h_dict = {"step": 1} @@ -207,7 +150,7 @@ def test_wind_meso_to_power_power_setpoint_applies(): test_h_dict = copy.deepcopy(h_dict_wind) test_h_dict["wind_farm"]["floris_update_time_s"] = 1.0 - wind_sim = Wind_MesoToPower(test_h_dict) + wind_sim = WindFarm(test_h_dict) # Set very low power setpoint values that should definitely limit power output step_h_dict = {"step": 1} @@ -229,7 +172,7 @@ def test_wind_meso_to_power_power_setpoint_applies(): def test_wind_meso_to_power_get_initial_conditions_and_meta_data(): """Test that get_initial_conditions_and_meta_data adds correct metadata to h_dict.""" - wind_sim = Wind_MesoToPower(h_dict_wind) + wind_sim = WindFarm(h_dict_wind) # Create a copy of the input h_dict to avoid modifying the original test_h_dict = copy.deepcopy(h_dict_wind) @@ -308,7 +251,7 @@ def test_wind_meso_to_power_regular_floris_updates(): test_h_dict["dt"] = 1.0 # Initialize wind simulation - wind_sim = Wind_MesoToPower(test_h_dict) + wind_sim = WindFarm(test_h_dict) # Run 5 steps with constant power setpoints floris_calc_counts = [] @@ -370,7 +313,7 @@ def test_wind_meso_to_power_power_setpoints_buffer(): test_h_dict["dt"] = 1.0 # Initialize wind simulation - wind_sim = Wind_MesoToPower(test_h_dict) + wind_sim = WindFarm(test_h_dict) # Run steps with varying power setpoints for step in range(5):