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*

bash_demo.sh

@@ -0,0 +1,26 @@
+# Example bash for running things locally
+# I just run these one at a t time
+
+# A lot of modules and conda stuff
+conda activate hercules
+
+# Set the helics port to use: 
+export HELICS_PORT=32000
+
+#make sure you use the same port number in the amr_input.inp and hercules_input_000.yaml files. 
+rm logsender logcomputer
+
+# Set up the helics broker
+helics_broker -t zmq  -f 2 --loglevel="debug" --local_port=$HELICS_PORT & 
+#helics_broker -f 2 --consoleloglevel=trace --loglevel=debug --local_port=$HELICS_PORT >> loghelics &
+
+# Need to set this to your hercules folder
+# cd /home/pfleming/hercules/hercules
+python3 seas_demo.py 0 >> logsender 2>&1 & # Start the controller center and pass in input file
+
+
+python3 seas_demo.py 1 >> logcomputer 2>&1
+# Now go back to scratch folder and launch the job
+
+# cd /scratch/pfleming/c2c/example_sim_02
+# mpirun -n 72 /home/pfleming/amr-wind/build/amr_wind amr_input.inp >> logamr 

hercules/amr_wind_standin.py

@@ -160,18 +160,24 @@ def run(self):
         # logger.info("** Initial Message Sent: {}".format(message_from_client_array))
 
         # Subscribe to helics messages:
-        incoming_messages = self.helics_connector.get_all_waiting_messages()
-        if incoming_messages != {}:
-            try:
-                message_from_server = list(ast.literal_eval(incoming_messages))
-            except Exception:
+        for _ in range(1):
+            incoming_messages = self.helics_connector.get_all_waiting_messages()
+            if incoming_messages != {}:
+                try:
+                    message_from_server = list(ast.literal_eval(incoming_messages))
+                except Exception:
+                    print("here here")
+                    message_from_server = None
+            else:
+                print("here 2")
                 message_from_server = None
-        else:
-            message_from_server = None
+            #self.sync_time_helics(self.absolute_helics_time)
 
         # Synchronize time bewteen control center and AMRWind
-        self.sync_time_helics(self.absolute_helics_time + self.deltat)
+        self.sync_time_helics(self.absolute_helics_time - self.deltat)
+        #self.sync_time_helics(0)
         sim_time_s = float(self.absolute_helics_time)
+        print("sim_time_s before while loop", sim_time_s)
 
         logger.info("** Initial Received reply: {}".format(message_from_server))
 
@@ -185,19 +191,22 @@ def run(self):
         #while sim_time_s <= (self.endtime - self.starttime):
             # SIMULATE A CALCULATION STEP IN AMR WIND=========================
             sim_time_s = float(self.absolute_helics_time)
+            print("sim_time_s inside while loop", sim_time_s)
             logger.info("Calculating simulation time: %.3f" % sim_time_s)
 
             # 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 +252,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,29 +134,61 @@ 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
+        print("First receive")
+#        for _ in range(2):
+        self.receive_amrwind_data()
+        print(self.main_dict)
+        print("Sending -1s")
         self.send_via_helics("control", str("[-1,-1,-1]"))
         print(" #### Entering main loop #### ")
-        self.sync_time_helics(self.absolute_helics_time + self.deltat)
+        self.sync_time_helics(self.absolute_helics_time - self.deltat)
+        print(self.absolute_helics_time)
         # Initialize the first iteration flag
         self.first_iteration = True
 
         # Run simulation till  endtime
         # while self.absolute_helics_time < self.endtime:
+        idx = 0
         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
             # helics time? Why aren't they the same?
             self.main_dict["time"] = self.absolute_helics_time
+            print(self.main_dict["external_signals"])
             self.main_dict = self.controller.step(self.main_dict)
+            print(self.main_dict['hercules_comms']['amr_wind']['wind_farm_0']['turbine_power_setpoints'])
             self.py_sims.step(self.main_dict)
             self.main_dict["py_sims"] = self.py_sims.get_py_sim_dict()
 
@@ -168,6 +211,9 @@ def run(self):
                 self.save_main_dict_as_text()
                 self.first_iteration = False
 
+            #self.sync_time_helics(self.absolute_helics_time + self.deltat)
+            idx += 1
+
     def receive_amrwind_data(self):
         # Subscribe to helics messages:
         incoming_messages = self.helics_connector.get_all_waiting_messages()
@@ -279,6 +325,8 @@ def log_main_dict(self):
         keys = list(self.main_dict_flat.keys())
         values = list(self.main_dict_flat.values())
 
+        print('OUTPUT FILE', self.output_file)
+
         # 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 +392,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):