Enable passing power setpoints to wind simulators

.gitignore

@@ -35,6 +35,7 @@ hercules/local_amr_wind_demo/sample_copy.nc
 #Ignore csv files
 *.csv
 !tests/test_inputs/amr_standin_data.csv
+!tests/test_inputs/external_data.csv
 
 # Some output files to ignore
 t_00*

hercules/amr_wind_standin.py

@@ -189,15 +189,17 @@ def run(self):
 
             # Compute the turbine power using a simple formula
             if self.message_from_server is not None:
-                yaw_angles = self.message_from_server[-self.num_turbines :]
+                yaw_angles = self.message_from_server[-2*self.num_turbines:-self.num_turbines]
+                power_setpoints = self.message_from_server[-self.num_turbines:]
             else:
                 yaw_angles = None
+                power_setpoints = None
             (
                 amr_wind_speed,
                 amr_wind_direction,
                 turbine_powers,
                 turbine_wind_directions,
-            ) = self.get_step(sim_time_s, yaw_angles)
+            ) = self.get_step(sim_time_s, yaw_angles, power_setpoints)
 
             # ================================================================
             # Communicate with control center
@@ -243,7 +245,7 @@ def run(self):
 
     # TODO cleanup code to move publish and subscribe here.
 
-    def get_step(self, sim_time_s, yaw_angles=None):
+    def get_step(self, sim_time_s, yaw_angles=None, power_setpoints=None):
         """Retreive or calculate wind speed, direction, and turbine powers
 
         Input:

hercules/emulator.py

@@ -4,6 +4,7 @@
 import sys
 
 import numpy as np
+import pandas as pd
 from SEAS.federate_agent import FederateAgent
 
 LOGFILE = str(dt.datetime.now()).replace(":", "_").replace(" ", "_").replace(".", "_")
@@ -18,9 +19,12 @@ def __init__(self, controller, py_sims, input_dict):
         self.main_dict_flat = {}
 
         # Initialize the output file
-        self.output_file = "hercules_output.csv"
+        if "output_file" in input_dict:
+            self.output_file = input_dict["output_file"]
+        else:
+            self.output_file = "hercules_output.csv"
 
-        # Save timt step
+        # Save time step
         self.dt = input_dict["dt"]
 
         # Initialize components
@@ -39,6 +43,13 @@ def __init__(self, controller, py_sims, input_dict):
         self.hercules_helics_dict = self.hercules_comms_dict["helics"]
         self.helics_config_dict = self.hercules_comms_dict["helics"]["config"]
 
+        # Read in any external data
+        self.external_data_all = {}
+        if "external_data_file" in input_dict:
+            self._read_external_data_file(input_dict["external_data_file"])
+            self.external_signals = {}
+            self.main_dict["external_signals"] = {}
+
         # Write the time step into helics config dict
         self.helics_config_dict["helics"]["deltat"] = self.dt
 
@@ -79,7 +90,7 @@ def __init__(self, controller, py_sims, input_dict):
             )
 
         # TODO For now, need to assume for simplicity there is one and only
-        # one AMR_Wind simualtion
+        # one AMR_Wind simulation
         self.num_turbines = self.amr_wind_dict[self.amr_wind_names[0]]["num_turbines"]
         self.rotor_diameter = self.amr_wind_dict[self.amr_wind_names[0]]["rotor_diameter"]
         self.turbine_locations = self.amr_wind_dict[self.amr_wind_names[0]]["turbine_locations"]
@@ -123,11 +134,29 @@ def __init__(self, controller, py_sims, input_dict):
         # list(self.pub.values())[0].publish(str("[-1,-1,-1]"))
         # self.logger.info(" #### Entering main loop #### ")
 
+    def _read_external_data_file(self, filename):
+        # Read in the external data file
+        df_ext = pd.read_csv(filename)
+        if "time" not in df_ext.columns:
+            raise ValueError("External data file must have a 'time' column")
+
+        # Interpolate the external data according to time
+        times = np.arange(
+            self.helics_config_dict["starttime"],
+            self.helics_config_dict["stoptime"],
+            self.dt
+        )
+        self.external_data_all["time"] = times
+        for c in df_ext.columns:
+            if c != "time":
+                self.external_data_all[c] = np.interp(times, df_ext.time, df_ext[c])
+
     def run(self):
         # TODO In future code that doesnt insist on AMRWInd can make this optional
         print("... waiting for initial connection from AMRWind")
         # Send initial connection signal to AMRWind
         # publish on topic: control
+        self.receive_amrwind_data()
         self.send_via_helics("control", str("[-1,-1,-1]"))
         print(" #### Entering main loop #### ")
         self.sync_time_helics(self.absolute_helics_time + self.deltat)
@@ -137,9 +166,15 @@ def run(self):
         # Run simulation till  endtime
         # while self.absolute_helics_time < self.endtime:
         while self.absolute_helics_time < (self.endtime - self.starttime + 1):
+            print(self.absolute_helics_time)
             # Loop till we reach simulation startime.
             # if self.absolute_helics_time < self.starttime:
             #     continue
+            # Get any external data
+            for k in self.external_data_all:
+                self.main_dict["external_signals"][k] = self.external_data_all[k][
+                    self.external_data_all["time"] == self.absolute_helics_time
+                ][0]
 
             # Update controller and py sims
             # TODO: Should 'time' in the main dict be AMR-wind time or
@@ -279,6 +314,7 @@ def log_main_dict(self):
         keys = list(self.main_dict_flat.keys())
         values = list(self.main_dict_flat.values())
 
+
         # If this is first iteration, write the keys as csv header
         if self.first_iteration:
             with open(self.output_file, "w") as filex:
@@ -344,10 +380,23 @@ def process_periodic_publication(self):
         else:  # set yaw_angles based on self.wind_direction
             yaw_angles = [self.wind_direction] * self.num_turbines
 
+        if (
+            "turbine_power_setpoints"
+            in self.main_dict["hercules_comms"]["amr_wind"][self.amr_wind_names[0]]
+        ):
+            power_setpoints = self.main_dict["hercules_comms"]["amr_wind"][self.amr_wind_names[0]][
+                "turbine_power_setpoints"
+            ]
+        else:  # set yaw_angles based on self.wind_direction
+            power_setpoints = [None] * self.num_turbines
+
         # Send timing and yaw information to AMRWind via helics
         # publish on topic: control
+        # TODO: power_setpoints part will not work with AMRWind proper.
         tmp = np.array(
-            [self.absolute_helics_time, self.wind_speed, self.wind_direction] + yaw_angles
+            [self.absolute_helics_time, self.wind_speed, self.wind_direction]
+            + yaw_angles
+            + power_setpoints
         ).tolist()
 
         self.send_via_helics("control", str(tmp))

hercules/floris_standin.py

@@ -115,6 +115,7 @@ def construct_floris_from_amr_input(amr_wind_input):
             ref_air_density=ref_air_density,
             generator_efficiency=1.0,
         )
+        turb_dict["power_thrust_model"] = "mixed"
 
         # load a default model
         fi = FlorisInterface(default_floris_dict)
@@ -145,16 +146,18 @@ def __init__(self, config_dict, amr_input_file, amr_standin_data_file=None):
         # Initialize storage
         self.yaw_angles_stored = [0.0] * self.num_turbines
 
-    def get_step(self, sim_time_s, yaw_angles=None):
+    def get_step(self, sim_time_s, yaw_angles=None, power_setpoints=None):
         """Retreive or calculate wind speed, direction, and turbine powers
 
         Input:
         sim_time_s: simulation time step
+        yaw_angles: absolute yaw positions for each turbine (deg). Defaults to None.
+        power_setpoints: power setpoints for each turbine (W). Defaults to None.
 
         Output:
-        amr_wind_speed: wind speed at current time step
-        amr_wind_direction: wind direction at current time step
-        turbine_powers: turbine powers at current time step
+        amr_wind_speed: wind speed at current time step [m/s]
+        amr_wind_direction: wind direction at current time step [deg]
+        turbine_powers: turbine powers at current time step [kW]
         """
 
         if hasattr(self, "standin_data"):
@@ -178,7 +181,10 @@ def get_step(self, sim_time_s, yaw_angles=None):
         # Compute the turbine power using FLORIS
         self.fi.reinitialize(wind_speeds=[amr_wind_speed], wind_directions=[amr_wind_direction])
 
-        if yaw_angles is not None:
+        if yaw_angles is None or (np.array(yaw_angles) == -1000).any():
+            # Note: -1000 is the "no value" flag for yaw_angles (NaNs not handled well)
+            yaw_misalignments = None
+        else:
             yaw_misalignments = (amr_wind_direction - np.array(yaw_angles))[
                 None, :
             ]  # TODO: remove 2
@@ -196,10 +202,15 @@ def get_step(self, sim_time_s, yaw_angles=None):
                 yaw_misalignments = (amr_wind_direction - np.array(self.yaw_angles_stored))[None, :]
             else:  # Reasonable yaw angles, save in case bad angles received later
                 self.yaw_angles_stored = yaw_angles
-        else:
-            yaw_misalignments = yaw_angles
-        self.fi.calculate_wake(yaw_angles=yaw_misalignments)
-        # This converts the output power from Floris (in Watts) to kW (standard for Hercules)
+
+        if power_setpoints is not None:
+            power_setpoints = np.array(power_setpoints)[None,:]
+            # Set invalid power setpoints to a large value
+            power_setpoints[power_setpoints == np.full(power_setpoints.shape, None)] = 1e9
+            power_setpoints[power_setpoints < 0] = 1e9
+            # Note conversion from Watts (used in Floris) and back to kW (used in Hercules)
+            power_setpoints = power_setpoints * 1000 # in W
+        self.fi.calculate_wake(yaw_angles=yaw_misalignments, power_setpoints=power_setpoints)
         turbine_powers = (self.fi.get_turbine_powers() / 1000).flatten().tolist()  # in kW
 
         return (
@@ -216,7 +227,7 @@ def process_periodic_endpoint(self):
         pass
 
     def process_periodic_publication(self):
-        # Periodically publish data to the surrpogate
+        # Periodically publish data to the surrogate
         pass
 
     def process_subscription_messages(self, msg):

tests/emulator_test.py

@@ -0,0 +1,82 @@
+from hercules.controller_standin import ControllerStandin
+from hercules.emulator import Emulator
+from hercules.py_sims import PySims
+
+# Copy data from hercules_input.yaml into input_dict
+test_input_dict = {
+    "name": "test_input_dict",
+    "description": "Test input dictionary for Hercules",
+    "dt": 1,
+    "hercules_comms": {
+        "amr_wind": {
+            "test_farm": {
+                "type": "amr_wind",
+                # requires running pytest from the top level directory
+                "amr_wind_input_file": "tests/test_inputs/amr_input_florisstandin.inp",
+            }
+        },
+        "helics": {
+            "config": {
+                "name": "hercules",
+                "use_dash_frontend": False,
+                "KAFKA": False,
+                "KAFKA_topics": "EMUV1py",
+                "helics": { 
+                    # deltat: 1 # This will be assigned in software
+                    "subscription_topics": ["status"],
+                    "publication_topics": ["control"],
+                    "endpoints": [],
+                    "helicsport" : 32000,
+                },
+                "publication_interval": 1,
+                "endpoint_interval": 1,
+                "starttime": 0,
+                "stoptime": 100,
+                "Agent": "ControlCenter",
+            },
+        },
+    },
+    "py_sims": {
+        "test_solar": {
+            "py_sim_type": "SimpleSolar",
+            "capacity": 50,
+            "efficiency": 0.5,
+            "initial_conditions": {
+                "power": 10.0,
+                "irradiance": 1000,
+            },
+        }
+    },
+    "controller": {
+        "num_turbines": 2,
+        "initial_conditions": {
+            "yaw": [270.0, 270.0],
+        },
+    },
+}
+
+def test_Emulator_instantiation():
+
+
+    controller = ControllerStandin(test_input_dict)
+    py_sims = PySims(test_input_dict)
+
+    emulator = Emulator(controller, py_sims, test_input_dict)
+
+    # Check default settings
+    assert emulator.output_file == "hercules_output.csv"
+    assert emulator.external_data_all == {}
+
+    test_input_dict_2 = test_input_dict.copy()
+    test_input_dict_2["external_data_file"] = "tests/test_inputs/external_data.csv"
+    test_input_dict_2["output_file"] = "test_output.csv"
+    test_input_dict_2["dt"] = 0.5
+
+    emulator = Emulator(controller, py_sims, test_input_dict_2)
+
+    assert emulator.external_data_all["power_reference"][0] == 1000
+    assert emulator.external_data_all["power_reference"][1] == 1500
+    assert emulator.external_data_all["power_reference"][2] == 2000
+    assert emulator.external_data_all["power_reference"][-1] == 3000
+
+    assert emulator.output_file == "test_output.csv"