question: D3 Part 2 Q2

battery_local_control.py

@@ -1,110 +0,0 @@
-from util import pos, clamp, soc_scaler, cast_to_cm
-import parameters
-from time import time, sleep
-import zmq
-import logging
-import syslab
-logging.basicConfig(level=logging.INFO)
-
-# Controller Parameters
-Kp = 0.4  # P factor for our controller.
-Ki = 0.8  # I factor for our controller.
-
-### Variables
-# Target that we are trying to reach at the grid connection.
-pcc_target = 0.0
-
-# Controller variables
-x_battery = 1  # Default splitting factor
-
-# Info on battery state
-battery_setpoint = 0.0  # Default setpoint
-battery_soc = 0.5  # Default SOC (= state of charge)
-
-### Communication
-# Make a context that we use to set up sockets
-context = zmq.Context()
-
-# Set up a socket we can use to publish our soc on
-soc_out_socket = context.socket(zmq.PUB)
-soc_out_socket.bind(f"tcp://*:{parameters.BATTERY_SOC_PORT}")
-
-# Set up a socket to subscribe to our splitting factor
-splitting_in_socket = context.socket(zmq.SUB)
-splitting_in_socket.connect(f"tcp://{parameters.SUPERVISOR_IP}:{parameters.SUPERVISOR_PORT}")
-
-# Ensure we only see message on the battery's splitting factor
-splitting_in_socket.subscribe(parameters.TOPIC_BATTERY_SPLITTING)
-
-### Unit connections
-# TODO step 1.2: Set up connection to control the battery and reconstruct the pcc (remember that vswitchboard is still not working)
-gaia = syslab.WindTurbine("vgaia1")
-dumpload = syslab.Dumpload("vload1")
-mobload1 = syslab.Dumpload("vmobload1")
-pv319 = syslab.Photovoltaics("vpv319")
-batt = syslab.Battery('vbatt1')
-
-### Import your controller class
-# TODO step 1.2: Import the controller class from "simlab_controller_d5_batt.py" or copy/paste it here and pick reasonable controller parameters
-# Note: The controller is identical to Day 5 with the exception of incorporating the splitting factor
-from simlab_controller_d5_batt import PIDController
-
-pid = PIDController(Kp=Kp, Ki=Ki, Kd=0.0,
-                    u_min=parameters.MIN_BATTERY_P,
-                    u_max=parameters.MAX_BATTERY_P,
-                    Iterm=0.0)
-
-
-# Put everything in a "try" block so clean-up is easy
-try:
-    while True:
-        try:
-            # Try to connect to the supervisor.
-            # If we have a connection to the supervisor, get our requested splitting factor.
-            # If none have come in, continue with previous splitting factor.
-            # We put this in a while-loop to ensure we empty the queue each time,
-            # so we always have the latest value.
-            while True:
-                # Receive the latest splitting factor
-                incoming_str = splitting_in_socket.recv_string(flags=zmq.NOBLOCK)
-                # The incoming string will look like "batt_split;0.781",
-                # so we split it up and take the last bit forward.
-                x_battery = float(incoming_str.split(" ")[-1])
-                logging.info(f"New splitting factor: {x_battery}")
-        except zmq.Again as e:
-            # No (more) new messages, move along
-            pass
-
-        # Poll the grid connection to get the current grid exchange.
-        pcc_p = cast_to_cm(-(
-                gaia.getActivePower().value 
-              + batt.getActivePower().value 
-              + pv319.getACActivePower().value
-              - dumpload.getActivePower().value 
-              - mobload1.getActivePower().value 
-        ))
-        # Check our own state of charge
-        battery_soc = batt.getSOC()  # TODO step 1.2: Read the true battery SOC
-
-        # Calculate new requests using PID controller
-        battery_setpoint = pid.update(pcc_p.value, x_batt=x_battery)
-
-        # Ensure we don't exceed our bounds for the battery
-        battery_setpoint = clamp(parameters.MIN_BATTERY_P, battery_setpoint, parameters.MAX_BATTERY_P)
-
-        # Send the new setpoint to the battery
-        # TODO step 1.2: Send the new setpoint to the battery
-        batt.setActivePower(battery_setpoint)
-        logging.info(f"Sent setpoint: {battery_setpoint}")
-
-        # Publish our current state of charge for the supervisory controller to see
-        soc_out_socket.send_string(f"{parameters.TOPIC_BATTERY_SOC} {battery_soc:.06f}")
-
-        # Loop once more in a second
-        sleep(1)
-finally:
-    # Clean up by closing our sockets.
-    # TODO step 1.2: Set the setpoint of the battery to zero after use
-    batt.setActivePower(0.)
-    splitting_in_socket.close()
-    soc_out_socket.close()

data/measurements/dummy_file

@@ -0,0 +1 @@
+Makes folder visible for git.

data/setpoints/dummy_file

@@ -0,0 +1 @@
+Makes folder visible for git.

data/setpoints/setpoints_1718013923.json

@@ -0,0 +1,4 @@
+{"unit": "dumpload_sp", "value": 0.0, "time": 1718013928.6948164}
+{"unit": "dumpload_sp", "value": 20.0, "time": 1718013959.6977577}
+{"unit": "dumpload_sp", "value": 10.0, "time": 1718013989.6989655}
+{"unit": "dumpload_sp", "value": 0.0, "time": 1718014018.716735}

demo_datalogger.py

@@ -0,0 +1,48 @@
+#!./venv/bin/python3
+import syslab
+from json import dump
+from time import sleep, time
+
+# Defines location and name of the measurements file
+LOG_FILE = f'data/measurements/measurements_{time():.00f}.json'
+print(f"Logging to file {LOG_FILE}")
+
+# Set up a connection to the switchboards
+sb3192 = syslab.SwitchBoard('319-2')
+sb3193 = syslab.SwitchBoard('319-3')
+sb33012 = syslab.SwitchBoard("330-12")
+sb1172 = syslab.SwitchBoard('117-2')
+
+# Convenience function to
+def take_measurements():
+    measurements = {
+        "pcc_p": sb3192.getActivePower('Grid'),
+        "pcc_q": sb3192.getReactivePower('Grid'),
+        "pv319_p": sb3192.getActivePower('PV'),
+        "pv319_q": sb3192.getReactivePower('PV'),
+        "dumpload_p": sb3192.getActivePower('Dumpload'),
+        "dumpload_q": sb3192.getReactivePower('Dumpload'),
+        "gaia_p": sb33012.getActivePower('Gaia'),
+        "gaia_q": sb33012.getReactivePower('Gaia'),
+        "pv330_p": sb33012.getActivePower('PV_1'),
+        "pv330_q": sb33012.getReactivePower('PV_1'),
+        "b2b_p": sb3193.getActivePower('ABB_Sec'),
+        "b2b_q": sb3193.getReactivePower('ABB_Sec'),
+        "battery_p": sb1172.getActivePower('Battery'),
+        "battery_q": sb1172.getReactivePower('Battery'),
+    }
+    return [{'unit': k, 'value': meas.value, 'time': meas.timestampMicros/1e6} for k, meas in measurements.items()]
+
+
+while True:
+    measurement = take_measurements()
+
+    # Open the output file in "append" mode which adds lines to the end
+    with open(LOG_FILE, 'a') as file:
+        for m in measurement:
+            # Convert the dictionary m to a json string and put it
+            # in the file.
+            dump(m, file)
+            # Write a newline for each measurement to make loading easier
+            file.write('\n')
+    sleep(1)

demo_plotter.py

@@ -0,0 +1,73 @@
+#!./venv/bin/python3
+import pandas as pd
+import numpy as np
+import json
+import matplotlib.pyplot as plt
+import os
+from datetime import timedelta
+
+## Read the measurements data file ##
+DATA_MEAS_DIR = 'data/measurements'
+# Always plot latest datafile - replace [-1] with another index if you want to plot a specific file.
+MEAS_LOG_FILE = sorted(os.listdir(DATA_MEAS_DIR))[-1]
+
+# Store each dictionary of the measurements json in a list
+with open(os.path.join(DATA_MEAS_DIR, MEAS_LOG_FILE)) as f:
+    meas_data = [json.loads(line) for line in f]
+
+# Use setpoint logger (only necessary for part two of the exercise "collecting fresh data")
+use_setpoint_log = False
+
+
+## Read the setpoints data file ##
+if use_setpoint_log:
+    DATA_SP_DIR = 'data/setpoints'
+    # Always plot latest datafile
+    SP_LOG_FILE = sorted(os.listdir(DATA_SP_DIR))[-1]
+
+    # Store each dictionary of the setpoints json in a list
+    with open(os.path.join(DATA_SP_DIR, SP_LOG_FILE)) as f:
+        sp_data = [json.loads(line) for line in f]
+
+    # Merge measurements and setpoints in one list
+    data = meas_data + sp_data
+
+else:
+    data = meas_data
+
+################################################################################
+################## Part 2 ######################################################
+################################################################################
+
+def overshoot(df, T1, T2):
+    yT1, yT2 = df[T1], df[T2]
+    over = 1 / (yT1 - yT2) * np.max(yT2 - df)
+    return over
+
+SETPOINT_UNIX = 1718013959.6977577
+SETPOINT_TS = pd.to_datetime(SETPOINT_UNIX, unit='s')
+WINDOW = pd.to_datetime(SETPOINT_UNIX+25, unit='s')
+
+## The controller is reasonably fast at reacting to changes; the sum of in and
+## out is at zero roughly 5-10 seconds after a change.
+
+# Construct a dataframe and pivot it to obtain a dataframe with a column per unit, and a row per timestamp.
+df = pd.DataFrame.from_records(data)
+df['time'] = pd.to_datetime(df['time'], unit='s')
+df_pivot = df.pivot_table(values='value', columns='unit', index='time')
+df_resampled = df_pivot.resample('0.1s').mean()
+df_resampled.interpolate(method='linear', inplace=True)
+df_resampled = pd.DataFrame(df_resampled)
+
+# Plot the data. Note, that the data will mostly not be plotted with lines.
+plt.ion()  # Turn interactive mode on
+plt.figure()
+ax1, ax2 = plt.subplot(211), plt.subplot(212)
+df_resampled[[c for c in df_resampled.columns if '_p' in c]].plot(marker='.', ax=ax1, linewidth=3)
+ax2.plot(df_resampled['pcc_p'][SETPOINT_TS:WINDOW], marker='.', linewidth=3, label='pcc_p')
+ax2.plot(df_resampled['dumpload_p'][SETPOINT_TS:WINDOW], marker='.', linewidth=3, label='dumpload')
+plt.legend()
+
+# print(overshoot(df_resampled['pcc_p'][SETPOINT_TS:WINDOW], SETPOINT_TS, WINDOW))
+
+plt.show(block=True)

mobileload_local_control.py

@@ -1,109 +0,0 @@
-from util import pos, clamp, cast_to_cm
-import parameters
-from time import time, sleep
-import zmq
-import syslab
-import logging
-
-# Ensure that we see "info" statements below
-logging.basicConfig(level=logging.INFO)
-
-# # Parameters
-Kp = 0.4  # P factor for our controller.
-Ki = 0.8  # I factor for our controller.
-
-# # Variables
-# Target that we are trying to reach at the grid connection.
-pcc_target = 0.0
-
-# Controller variables
-x_load = 1  # Default splitting factor
-
-# Info on battery state
-mobileload_setpoint = 0.0  # Default setpoint
-battery_soc = 0.5  # Default SOC (= state of charge)
-
-### Communication
-# Handle the sockets we need
-# Make a context that we use to set up sockets
-context = zmq.Context()
-
-# Set up a socket to subscribe to our splitting factor
-splitting_in_socket = context.socket(zmq.SUB)
-splitting_in_socket.connect(f"tcp://{parameters.SUPERVISOR_IP}:{parameters.SUPERVISOR_PORT}")
-
-# Ensure we only see message on the load's splitting factor
-splitting_in_socket.subscribe(parameters.TOPIC_LOAD_SPLITTING)
-
-### Unit connections
-# TODO step 1.2: Set up connection to control the battery and reconstruct the pcc (remember that vswitchboard is still not working)
-gaia = syslab.WindTurbine("vgaia1")
-dumpload = syslab.Dumpload("vload1")
-mobload1 = syslab.Dumpload("vmobload1")
-pv319 = syslab.Photovoltaics("vpv319")
-batt = syslab.Battery('vbatt1')
-
-### Import your controller class
-# TODO step 1.2: Import the controller class from "simlab_controller_d5_load.py" or copy/paste it here and pick reasonable controller parameters
-# Note: The controller is identical to Day 5 with the exception of incorporating the splitting factor
-from simlab_controller_d5_load import PIDController
-
-pid = PIDController(Kp=Kp, Ki=Ki, Kd=0.0,
-                    u_min=parameters.MIN_LOAD_P,
-                    u_max=parameters.MAX_LOAD_P,
-                    Iterm=0.0)
-
-
-# Ensure the mobile load is on before we start (The sleeps wait for the load to respond.)
-while not mobload1.isLoadOn().value:
-    print("Starting load")
-    mobload1.startLoad()
-    sleep(2)
-
-
-try:
-    while True:
-        try:
-            # Try to connect the the supervisor.
-            # If we have a connection to the supervisor, get our requested splitting factor.
-            # If none have come in, continue with previous splitting factor.
-            # We put this in a while-loop to ensure we empty the queue each time.
-            while True:
-                # Receive the latest splitting factor
-                incoming_str = splitting_in_socket.recv_string(flags=zmq.NOBLOCK)
-                # The incoming string will look like "load_split;0.781",
-                # so we split it up and take the last part.
-                x_load = float(incoming_str.split(" ")[-1])
-                logging.info(f"New splitting factor: {x_load}")
-        except zmq.Again as e:
-            # No (more) new messages, move along
-            pass
-
-        # Poll the grid connection to get the current grid exchange.
-        pcc_p = cast_to_cm(-(
-                gaia.getActivePower().value 
-              + batt.getActivePower().value 
-              + pv319.getACActivePower().value
-              - dumpload.getActivePower().value 
-              - mobload1.getActivePower().value 
-        ))
-
-        # Calculate new requests using PID controller
-        mobileload_setpoint = pid.update(pcc_p.value, x_load=x_load)
-
-        # Ensure we don't exceed our bounds for the load
-        mobileload_setpoint = clamp(parameters.MIN_LOAD_P, mobileload_setpoint, parameters.MAX_LOAD_P)
-
-        # Send the new setpoint to the load
-        # TODO step 1.2: Send the new setpoint to the load
-        mobload1.setPowerSetPoint(mobileload_setpoint)
-        logging.info(f"Sent setpoint: {mobileload_setpoint}")
-
-        # Loop once more in a second
-        sleep(1)
-finally:
-    # Clean up by closing our socket.
-    # TODO step 1.2: Set the setpoint of the mobile load to zero and shut it down after use
-    mobload1.setPowerSetPoint(0.)
-    mobload1.stopLoad()
-    splitting_in_socket.close()

parameters.py

@@ -1,21 +0,0 @@
-
-# Parameters
-
-# Location of supervisor
-SUPERVISOR_IP = 'localhost'
-SUPERVISOR_PORT = 6001
-
-# Location of battery
-BATTERY_SOC_IP = 'localhost'
-BATTERY_SOC_PORT = 6002
-
-# Unit parameters
-MAX_BATTERY_P = 15
-MIN_BATTERY_P = -15
-MAX_LOAD_P = 15.0
-MIN_LOAD_P = 0
-
-# Topics for pub/sub
-TOPIC_BATTERY_SPLITTING = "batt_split"
-TOPIC_BATTERY_SOC = "batt_soc"
-TOPIC_LOAD_SPLITTING = "load_split"

plotter.py

@@ -0,0 +1,28 @@
+#!./venv/bin/python3
+import pandas as pd
+import json
+import matplotlib.pyplot as plt
+import os
+
+DATA_MEAS_DIR = 'data/measurements'
+SPECIFIC_FILE = ''
+MEAS_LOG_FILE = sorted(os.listdir(DATA_MEAS_DIR))[-1] if not SPECIFIC_FILE else SPECIFIC_FILE
+
+with open(os.path.join(DATA_MEAS_DIR, MEAS_LOG_FILE)) as f:
+    meas_data = [json.loads(line) for line in f]
+data = meas_data
+
+df = pd.DataFrame.from_records(data)
+df['time'] = pd.to_datetime(df['time'], unit='s')
+df_pivot = df.pivot_table(values='value', columns='unit', index='time')
+df_resampled = df_pivot.resample('s').mean()
+df_resampled.interpolate(method='linear', inplace=True)
+df_resampled = pd.DataFrame(df_resampled)
+
+# Plot the data. Note, that the data will mostly not be plotted with lines.
+plt.ion()  # Turn interactive mode on
+plt.figure()
+ax1, ax2 = plt.subplot(211), plt.subplot(212)
+df_resampled[[c for c in df_resampled.columns if '_p' in c]].plot(marker='.', ax=ax1, linewidth=3)
+df_resampled[[c for c in df_resampled.columns if '_q' in c]].plot(marker='.', ax=ax2, linewidth=3)
+plt.show(block=True)

readme.md

@@ -1,37 +0,0 @@
-# Power-sharing controller 
-
-This is an example of a power sharing controller for the Microgrid case.
-
-The case consists of the following units:
-- A battery
-- A mobile load, representing controllable consumption that gives us money, e.g. an electrolyzer which produces hydrogen.
-- A wind turbine
-- A solar panel
-- A dump load, which is controlled externally from the controller, representing load that must be served, for instance the local neighbourhood.
-
-The controller uses the battery and mobile loads to attempt to keep the active power at 0 on the grid connection (PCC).
-Further, we try to keep the battery's state of charge somewhat away from being full or empty, to give the most flexibility if we are called on to deliver services.
-The service delivery is not implemented here, and so we just try to keep the battery in a certain range of state of charge.
-
-If we need more power in the system, i.e. we are importing from the grid, the battery will be asked to provide the missing amount, up to its maximum.
-
-If we have too much power, both the battery and mobile load will act together to consume the excess.
-How much will be consumed by each unit is determined by a *splitting factor*:
-- If the battery's state of charge is less than an `soc_lower`, the battery will consume all excess power.
-- If the battery's state of charge is above an `soc_upper`, the mobile load will consume all excess power.
-- If the battery's state of charge is *in between* these two values, the excess will be split between the two in a manner which linearly depends on the battery's state of charge; the higher the state of charge, the more power is given to the mobile load.
-
-See the included "d7_controller_splitting.pdf" for a graph of the splitting factor.
-Essentially, we are trying to keep the battery charged to above `soc_lower` plus some buffer, and then use the rest of the power for production.
-
-In addition to changing the splitting factor, the boundaries `soc_lower` and `soc_upper` change depending on the currently available renewable energy:
-- If there is a lot of renewable energy available, the range from `soc_lower` to `soc_upper` gets wider. We are less worried about getting the battery charged, so we use more power for the dump load.
-- If there is only a small amount of renewable energy available, the range from `soc_lower` to `soc_upper` gets narrower. We need to ensure the battery gets charged up so we are ready to handle a high load or a request for energy from elsewhere.
-
-## Included scripts
-
-- `battery_local_control.py`: Connects to the battery and acts on the splitting factor delivered by the supervisor. Makes the battery state of charge available to the supervisor.
-- `mobileload_local_control.py`: Connects to the battery and acts on the splitting factor delivered by the supervisor.
-- `supervisory_controller.py`: Is in charge of calculating the splitting factor depending on the available renewable production.
-
-Each script can be started independently, and will use defaults or re-use last known values if the other scripts are down.

simlab_controller_d5_batt.py

@@ -1,178 +0,0 @@
-import syslab
-from time import sleep, time
-from dataclasses import dataclass
-from util import clamp, cast_to_cm
-import sys
-import zmq
-import logging
-
-# Make a context that we use to set up sockets
-context = zmq.Context()
-
-# Set up a socket to subscribe to the tester
-reference_in_socket = context.socket(zmq.SUB)
-reference_in_socket.connect(f"tcp://localhost:5000")
-
-# Ensure we only see message on the setpoint reference
-reference_in_socket.subscribe("setpoint_reference")
-
-# Used to log the incoming reference setpoints
-logging.basicConfig(level=logging.INFO)
-
-@dataclass
-class PIDController:
-    Kp: float = 0.0  # P factor
-    Ki: float = 0.0  # I factor
-    Kd: float = 0.0  # D factor
-    r: float = 0.0  # Reference which we want the measured y to be
-    Pterm: float = 0.0  # P part
-    Iterm: float = 0.0  # Integral part
-    Dterm: float = 0.0  # Differential part
-    previous_error: float = 0.0  # used to calculate differential part
-    previous_Iterm: float = 0.0  # used to calculate integral part
-    current_time: float = None  # used to calculate differential and integral part
-    previous_time: float = None  # used to calculate differential and integral part
-    u_max: float = None  # used to calculate controller input saturation
-    u_min: float = None  # used to calculate controller input saturation
-    u_p: float = 0.0
-    u: float = 0.0
-
-    def __post_init__(self):
-        self.previous_time = self.previous_time if self.previous_time is not None else time()
-        self.current_time = self.current_time if self.current_time is not None else time()
-
-    def set_reference(self, new_reference):
-        self.r = new_reference
-
-    def update(self, new_y_meas, x_batt, current_time=None):
-
-        # Ensure we have the time elapsed since our last value; we'll use that for the integral and differential parts
-        self.current_time = current_time if current_time is not None else time()
-        delta_time = self.current_time - self.previous_time
-
-        # Filter the incoming value
-        y_2 = self._filter_input(new_y_meas)
-
-        # Calculate the error to the setpoint
-        error = self._calculate_error(self.r, y_2)
-
-        # Apply splitting factor  # TODO step 1.2: Familiarize with the usage of the splitting factor
-        error = error - (1-x_batt)*max(error, 0)
-
-        # Calculate the PID terms
-        Pterm = self._calc_Pterm(error)
-        Iterm = self._calc_Iterm(error, delta_time)
-        Dterm = self._calc_Dterm(error, delta_time)
-
-        # Calculate our raw response
-        self.u_p = Pterm + Iterm + Dterm
-
-        # Filter our response
-        self.u = self._filter_output(self.u_p)
-        u_new = clamp(self.u_min, self.u, self.u_max)
-
-        # Remember values for next update
-        if self.u == u_new:
-            self.previous_Iterm = Iterm
-        else:  # in saturation - don't sum up Iterm
-            self.u = u_new
-
-        self.previous_time = self.current_time
-        self.previous_error = error
-
-        # Return the filtered response
-        return self.u
-
-    def _filter_input(self, new_y_hat):
-        # optional code to filter the measurement signal
-        return new_y_hat
-
-    def _filter_output(self, u_p):
-        # optional code to filter the input / actuation signal
-        return u_p
-
-    def _calculate_error(self, y_hat_2):
-        # TODO step 1.2: calculate and return the error between reference and measurement
-        return self.r - y_hat_2
-
-    def _calc_Pterm(self, error):
-        # TODO step 1.2: calculate the proportional term based on the error and self.Kp
-        return error * self.Kp
-
-    def _calc_Iterm(self, error, delta_time):
-        # TODO step 1.2:  calculate the integral term based on error, last error and deltaT, and self.Ki
-        return 0.0
-
-    def _calc_Dterm(self, error, delta_time):
-        # calculate the proportional term based on error, last error and deltaT
-        return self.Kd * 0.0
-
-if __name__ == "__main__":
-    # "main" method if this class is run as a script
-    BATT_MIN_P, BATT_MAX_P, = -15, 15
-    #LOAD_MIN_P, LOAD_MAX_P = 0, 20
-    pid = PIDController(Kp=0.4, Ki=0.8, Kd=0.0,  # 0.5 , 0.2 / 0.4 , 0.1, Kp 0.7
-                        u_min=BATT_MIN_P,
-                        u_max=BATT_MAX_P)
-
-    # Establish connection to the switchboard to get the pcc power measurement (that's the usual practice, however, vswitchboard currently not working)
-    # sb = syslab.SwitchBoard('vswitchboard')
-
-    # Instead, establish connection to all units to replace switchboard measurements
-    gaia = syslab.WindTurbine("vgaia1")
-    dumpload = syslab.Dumpload("vload1")
-    mobload1 = syslab.Dumpload("vmobload1")
-    pv319 = syslab.Photovoltaics("vpv319")
-    batt = syslab.Battery('vbatt1')
-
-    # Initiate reference signal
-    r_old = 0.0
-    r = 0.0
-
-    try:
-        while True:
-
-            # We loop over the broadcast queue to empty it and only keep the latest value
-            try:
-                while True:
-                    # Receive reference setpoint
-                    incoming_str = reference_in_socket.recv_string(flags=zmq.NOBLOCK)
-                    # The incoming string will look like "setpoint_reference;0.547",
-                    # so we split it up and take the last bit forward.
-                    r = float(incoming_str.split(" ")[-1])
-                    #logging.info(f"New reference setpoint: {r}")
-
-            # An error will indicate that the queue is empty so we move along.
-            except zmq.Again as e:
-                pass
-
-            # Overwrite setpoint reference r if different from previous value
-            pid.set_reference(r) if r != r_old else None
-
-            # Store reference for comparison in next loop
-            r_old = r
-
-            # To ensure accurate
-            start_time = time()
-            # Get measurements of the current power exchange at the grid connection (=Point of Common Coupling (PCC))
-            #pcc_p = sb.getActivePower(0)
-            pcc_p = cast_to_cm(-(gaia.getActivePower().value - dumpload.getActivePower().value - mobload1.getActivePower().value + batt.getActivePower().value + pv319.getACActivePower().value))
-
-            # Calculate new requests using PID controller
-            p_request_total = pid.update(pcc_p.value)
-
-            # Ensure requests do not exceed limits of battery (extra safety, should be ensured by controller)
-            batt_p = clamp(BATT_MIN_P, p_request_total, BATT_MAX_P)
-
-            # # Send new setpoints
-            batt.setActivePower(batt_p)
-
-            #print(f"Measured: P:{pcc_p.value:.02f}, want to set: P {batt_p:.02f}")
-            print(f"Setpoint:{r:.02f} Measured: P:{pcc_p.value:.02f}, want to set: P {batt_p:.02f}")
-
-            # Run the loop again once at most a second has passed.
-            sleep(1)
-    finally:
-        # When closing, set all setpoints to 0.
-        batt.setActivePower(0.0)
-        reference_in_socket.close()

simlab_controller_d5_load.py

@@ -1,178 +0,0 @@
-import syslab
-from time import sleep, time
-from dataclasses import dataclass
-from util import clamp, cast_to_cm
-import sys
-import zmq
-import logging
-
-# Make a context that we use to set up sockets
-context = zmq.Context()
-
-# Set up a socket to subscribe to the tester
-reference_in_socket = context.socket(zmq.SUB)
-reference_in_socket.connect(f"tcp://localhost:5000")
-
-# Ensure we only see message on the setpoint reference
-reference_in_socket.subscribe("setpoint_reference")
-
-# Used to log the incoming reference setpoints
-logging.basicConfig(level=logging.INFO)
-
-@dataclass
-class PIDController:
-    Kp: float = 0.0  # P factor
-    Ki: float = 0.0  # I factor
-    Kd: float = 0.0  # D factor
-    r: float = 0.0  # Reference which we want the measured y to be
-    Pterm: float = 0.0  # P part
-    Iterm: float = 0.0  # Integral part
-    Dterm: float = 0.0  # Differential part
-    previous_error: float = 0.0  # used to calculate differential part
-    previous_Iterm: float = 0.0  # used to calculate integral part
-    current_time: float = None  # used to calculate differential and integral part
-    previous_time: float = None  # used to calculate differential and integral part
-    u_max: float = None  # used to calculate controller input saturation
-    u_min: float = None  # used to calculate controller input saturation
-    u_p: float = 0.0
-    u: float = 0.0
-
-    def __post_init__(self):
-        self.previous_time = self.previous_time if self.previous_time is not None else time()
-        self.current_time = self.current_time if self.current_time is not None else time()
-
-    def set_reference(self, new_reference):
-        self.r = new_reference
-
-    def update(self, new_y_meas, x_load, current_time=None):
-
-        # Ensure we have the time elapsed since our last value; we'll use that for the integral and differential parts
-        self.current_time = current_time if current_time is not None else time()
-        delta_time = self.current_time - self.previous_time
-
-        # Filter the incoming value
-        y_2 = self._filter_input(new_y_meas)
-
-        # Calculate the error to the setpoint
-        error = self._calculate_error(self.r, y_2)
-
-        # Apply splitting factor  # TODO step 1.2: Familiarize with the usage of the splitting factor
-        error = x_load*max(error, 0)
-
-        # Calculate the PID terms
-        Pterm = self._calc_Pterm(error)
-        Iterm = self._calc_Iterm(error, delta_time)
-        Dterm = self._calc_Dterm(error, delta_time)
-
-        # Calculate our raw response
-        self.u_p = Pterm + Iterm + Dterm
-
-        # Filter our response
-        self.u = self._filter_output(self.u_p)
-        u_new = clamp(self.u_min, self.u, self.u_max)
-
-        # Remember values for next update
-        if self.u == u_new:
-            self.previous_Iterm = Iterm
-        else:  # in saturation - don't sum up Iterm
-            self.u = u_new
-
-        self.previous_time = self.current_time
-        self.previous_error = error
-
-        # Return the filtered response
-        return self.u
-
-    def _filter_input(self, new_y_hat):
-        # optional code to filter the measurement signal
-        return new_y_hat
-
-    def _filter_output(self, u_p):
-        # optional code to filter the input / actuation signal
-        return u_p
-
-    def _calculate_error(self, y_hat_2):
-        # TODO step 1.2: calculate and return the error between reference and measurement
-        return self.r - y_hat_2
-
-    def _calc_Pterm(self, error):
-        # TODO step 1.2: calculate the proportional term based on the error and self.Kp
-        return error * self.Kp
-
-    def _calc_Iterm(self, error, delta_time):
-        # TODO step 1.2:  calculate the integral term based on error, last error and deltaT, and self.Ki
-        return 0.0
-
-    def _calc_Dterm(self, error, delta_time):
-        # calculate the proportional term based on error, last error and deltaT
-        return self.Kd * 0.0
-
-if __name__ == "__main__":
-    # "main" method if this class is run as a script
-    BATT_MIN_P, BATT_MAX_P, = -15, 15
-    #LOAD_MIN_P, LOAD_MAX_P = 0, 20
-    pid = PIDController(Kp=0.4, Ki=0.8, Kd=0.0,  # 0.5 , 0.2 / 0.4 , 0.1, Kp 0.7
-                        u_min=BATT_MIN_P,
-                        u_max=BATT_MAX_P)
-
-    # Establish connection to the switchboard to get the pcc power measurement (that's the usual practice, however, vswitchboard currently not working)
-    # sb = syslab.SwitchBoard('vswitchboard')
-
-    # Instead, establish connection to all units to replace switchboard measurements
-    gaia = syslab.WindTurbine("vgaia1")
-    dumpload = syslab.Dumpload("vload1")
-    mobload1 = syslab.Dumpload("vmobload1")
-    pv319 = syslab.Photovoltaics("vpv319")
-    batt = syslab.Battery('vbatt1')
-
-    # Initiate reference signal
-    r_old = 0.0
-    r = 0.0
-
-    try:
-        while True:
-
-            # We loop over the broadcast queue to empty it and only keep the latest value
-            try:
-                while True:
-                    # Receive reference setpoint
-                    incoming_str = reference_in_socket.recv_string(flags=zmq.NOBLOCK)
-                    # The incoming string will look like "setpoint_reference;0.547",
-                    # so we split it up and take the last bit forward.
-                    r = float(incoming_str.split(" ")[-1])
-                    #logging.info(f"New reference setpoint: {r}")
-
-            # An error will indicate that the queue is empty so we move along.
-            except zmq.Again as e:
-                pass
-
-            # Overwrite setpoint reference r if different from previous value
-            pid.set_reference(r) if r != r_old else None
-
-            # Store reference for comparison in next loop
-            r_old = r
-
-            # To ensure accurate
-            start_time = time()
-            # Get measurements of the current power exchange at the grid connection (=Point of Common Coupling (PCC))
-            #pcc_p = sb.getActivePower(0)
-            pcc_p = cast_to_cm(-(gaia.getActivePower().value - dumpload.getActivePower().value - mobload1.getActivePower().value + batt.getActivePower().value + pv319.getACActivePower().value))
-
-            # Calculate new requests using PID controller
-            p_request_total = pid.update(pcc_p.value)
-
-            # Ensure requests do not exceed limits of battery (extra safety, should be ensured by controller)
-            batt_p = clamp(BATT_MIN_P, p_request_total, BATT_MAX_P)
-
-            # # Send new setpoints
-            batt.setActivePower(batt_p)
-
-            #print(f"Measured: P:{pcc_p.value:.02f}, want to set: P {batt_p:.02f}")
-            print(f"Setpoint:{r:.02f} Measured: P:{pcc_p.value:.02f}, want to set: P {batt_p:.02f}")
-
-            # Run the loop again once at most a second has passed.
-            sleep(1)
-    finally:
-        # When closing, set all setpoints to 0.
-        batt.setActivePower(0.0)
-        reference_in_socket.close()

supervisory_controller.py

@@ -1,93 +0,0 @@
-from util import pos, clamp
-import parameters
-from time import time, sleep
-import zmq
-import syslab
-from syslab.comm.LogUtils import setup_event_logger
-logging = setup_event_logger()
-
-
-### Parameters
-# If soc is below this limit, divert all power to the battery.
-base_soc_lower_limit = 0.2
-# If soc is above this limit, divert all power to the mobile load.
-base_soc_upper_limit = 0.8
-# This much renewable production shifts our soc_lower_limit down by 0.1 and soc_upper_limit up by 0.1.
-# I.e. if there is a lot of renewable production, the range over which we split between battery and
-# load widens.
-base_renewable_shift = 10.0
-
-
-# # Variables
-# Controller variables
-x_battery = 0.5  # Default splitting factor
-x_load = 0.5  # Default splitting factor
-# Info on battery state
-battery_soc = 0.5  # Default SOC (= state of charge)
-
-soc_lower_limit = 0.2
-soc_upper_limit = 0.8
-
-# # Communication
-# Handle the sockets we need
-# Make a context that we use to set up sockets
-context = zmq.Context()
-
-# Set up a socket we can use to broadcast splitting factors on
-splitting_out_socket = context.socket(zmq.PUB)
-splitting_out_socket.bind(f"tcp://*:{parameters.SUPERVISOR_PORT}")
-
-# Set up a socket to subscribe to the battery's soc
-soc_in_socket = context.socket(zmq.SUB)
-soc_in_socket.connect(f"tcp://{parameters.BATTERY_SOC_IP}:{parameters.BATTERY_SOC_PORT}")
-
-# Ensure we only see message on the battery's soc
-soc_in_socket.subscribe(parameters.TOPIC_BATTERY_SOC)
-
-### Unit connections
-# TODO step 1.3: Set up connection to the units
-
-try:
-    while True:
-        try:
-            # Try to connect the the battery.
-            # If we have a connection to the battery, get its current soc.
-            # If none have come in, continue with previous soc.
-            # We put this in a while-loop to ensure we empty the queue each time,
-            # so we always have the latest value.
-            while True:
-                # Receive the latest splitting factor
-                incoming_str = soc_in_socket.recv_string(flags=zmq.NOBLOCK)
-                # The incoming string will look like "batt_split;0.781",
-                # so we split it up and take the last bit forward.
-                battery_soc = float(incoming_str.split(" ")[-1])
-                logging.info(f"New battery soc: {battery_soc}")
-        except zmq.Again as e:
-            # No (more) new messages, move along
-            pass
-
-        # Poll the grid connection to get the current production of renewables.
-        wind_p = 4.0  # TODO Task 1.3: Change into syslab connection. Remember that the measurements need to be collected directly from the units since vswitchboard is not working
-        solar_p = 7.0
-        logging.info(f"Current renewable production: {wind_p + solar_p} kW.")
-
-        soc_lower_limit = base_soc_lower_limit - (wind_p + solar_p) / base_renewable_shift / 10.0
-        soc_upper_limit = base_soc_upper_limit + (wind_p + solar_p) / base_renewable_shift / 10.0
-
-        # Calculate a new splitting factor for the battery.
-        x_battery = clamp(0, (soc_upper_limit - battery_soc)/(soc_upper_limit - soc_lower_limit), 1)
-
-        # If anything is not taken by the battery, put it in the load.
-        x_load = 1 - x_battery
-
-        # Publish the new splitting factors.
-        splitting_out_socket.send_string(f"{parameters.TOPIC_BATTERY_SPLITTING} {x_battery:.06f}")
-        splitting_out_socket.send_string(f"{parameters.TOPIC_LOAD_SPLITTING} {x_load:.06f}")
-        logging.info(f"Battery split: {x_battery:.03f}; Load split: {x_load:.03f}")
-
-        # Loop once more in a second
-        sleep(1)
-finally:
-    # Clean up by closing our sockets.
-    soc_in_socket.close()
-    splitting_out_socket.close()

syslab-python/.gitignore

@@ -0,0 +1,109 @@
+# Byte-compiled / optimized / DLL files
+__pycache__/
+*.py[cod]
+*$py.class
+
+# C extensions
+*.so
+
+# Distribution / packaging
+.Python
+build/
+develop-eggs/
+dist/
+downloads/
+eggs/
+.eggs/
+lib/
+lib64/
+parts/
+sdist/
+var/
+wheels/
+*.egg-info/
+.installed.cfg
+*.egg
+MANIFEST
+
+# PyInstaller
+#  Usually these files are written by a python script from a template
+#  before PyInstaller builds the exe, so as to inject date/other infos into it.
+*.manifest
+*.spec
+
+# Installer logs
+pip-log.txt
+pip-delete-this-directory.txt
+
+# Unit test / coverage reports
+htmlcov/
+.tox/
+.coverage
+.coverage.*
+.cache
+nosetests.xml
+coverage.xml
+*.cover
+.hypothesis/
+.pytest_cache/
+
+# Translations
+*.mo
+*.pot
+
+# Django stuff:
+*.log
+local_settings.py
+db.sqlite3
+
+# Flask stuff:
+instance/
+.webassets-cache
+
+# Scrapy stuff:
+.scrapy
+
+# Sphinx documentation
+docs/_build/
+
+# PyBuilder
+target/
+
+# Jupyter Notebook
+.ipynb_checkpoints
+
+# pyenv
+.python-version
+
+# celery beat schedule file
+celerybeat-schedule
+
+# SageMath parsed files
+*.sage.py
+
+# Environments
+.env
+.venv
+env/
+venv/
+ENV/
+env.bak/
+venv.bak/
+
+# Spyder project settings
+.spyderproject
+.spyproject
+
+# Rope project settings
+.ropeproject
+
+# mkdocs documentation
+/site
+
+# mypy
+.mypy_cache/
+
+\.DS_Store
+
+# Pycharm
+.idea/

syslab-python/README.md

@@ -0,0 +1,27 @@
+# SYSLAB Python Interface
+
+This project provides a Python interface to several SYSLAB components.
+
+## Required software
+
+- Python (>=3.7)
+- python-requests (>=2.18)
+- python-bs4 (>= 4.7)
+
+## Installation
+
+To install the package in your local path, run
+
+```shell
+    python setup.py install
+```
+
+Alternatively, copy the `syslab` folder into your project directory.
+
+# Contributors
+
+- Anders Thavlov: Initial implementation
+- Oliver Gehrke
+- Daniel Esteban Morales Bondy
+- Tue Vissing Jensen
+- Federico Zarelli

syslab-python/example.py

@@ -0,0 +1,11 @@
+from syslab import SwitchBoard
+
+name = '319-2'
+SB_connection = SwitchBoard(name)
+print("Let's look at what is going on in the switchboard {}.".format(name))
+
+for bay in range(SB_connection.getNumBays()):
+	print(SB_connection.getBayName(bay),' : ',SB_connection.getActivePower(bay))
+
+
+

syslab-python/notes.md

@@ -0,0 +1,99 @@
+# Notes on SOAP implementation
+
+    - Units should reflect SOAP methods onto their own namespace (1)
+    - CompositeMeasurement should convert to/from SOAP
+    - Type checking via client.get\_type('ns0:compositeMeasurement').elements
+    - 
+
+# Notes about this module - to be discussed
+
+SYSLAB\_Unit.py has a whole bunch of static methods - should these be split into a util library instead?
+Generally, many places where static methods are used for things that should perhaps just be functions...
+
+The following files are empty:
+    - BattOpMode.py
+    - FlowBatteryState.py
+    - GaiaWindTurbine.py
+
+
+To check the methods available, use, e.g.:
+http://syslab-33.syslab.dk:8080/typebased_WebService_HeatSubstation/HeatSwitchboardWebService/716-h1/resourceNames
+
+To figure out this URL, look at software.xml for the corresponding machine, and use this template:
+
+http://(machineName).syslab.dk:(port)/(interfaceName)/(shortServerName)/(unitname)/resourceNames
+
+| field | corresponds to | notes |
+| ----- | -------------- | ----- |
+| machineName | N/A | Look this up on the wiki |
+| port | N/A | Dynamically allocated, starting at 8080 - good luck! |
+| interfaceName | typeBasedWebService, interfaceName | |
+| shortServerName | typeBasedWebService, serverClass | Remove the "Server" at the end |
+| unitname | dataLogger, unit | Also defined as "name" in hardware.xml |
+
+
+-------------------------
+
+
+SYSLAB COMMON
+Broadcast event logger:
+
+Transcode to python:
+https://git.elektro.dtu.dk/syslab/syslab-common/-/blob/master/src/main/java/risoe/syslab/comm/broadcast/BroadcastLogSender.java
+
+broadcast log sender 
+:: transcode the "send()" method to python. It byte-encodes the message for UDP. 
+Java: 
+send(String origin, byte[] origIP, long timestamp, int ploadType, String message,
+                    int level, int flags, String[] tags)
+------------
+Python:
+----------.
+def send(origin, origIP, timestamp, ploadType, message, level, flags, tags):
+    ploadbytes = message[:min(1024, len(message))].encode()
+    origbytes = origin[:min(32, len(origin))].encode()
+    tagbytes = tagsToBytes(tags, 256)
+    pktlen = 2 + 2 + 1 + len(origbytes) + 4 + 2 + 2 + 8 + 1 + 2 + len(ploadbytes) + len(tagbytes)
+    buf = bytearray(pktlen)
+    buf[0] = BroadcastLogConstants.BROADCASTLOG_PKTID >> 8
+    buf[1] = BroadcastLogConstants.BROADCASTLOG_PKTID & 0xff
+    buf[2] = (pktlen >> 8) & 0xff
+    buf[3] = pktlen & 0xff
+    buf[4] = len(origbytes)
+    buf[5:5+len(origbytes)] = origbytes
+    writePtr = 5 + len(origbytes)
+    buf[writePtr:writePtr+4] = origIP
+    writePtr += 4
+    buf[writePtr] = (level >> 8) & 0xff
+    buf[writePtr+1] = level & 0xff
+    buf[writePtr+2] = (flags >> 8) & 0xff
+    buf[writePtr+3] = flags & 0xff
+    for i in range(8):
+        buf[writePtr+7-i] = timestamp & 0xff
+        timestamp >>= 8
+    writePtr += 8
+    buf[writePtr] = ploadType & 0xff
+    buf[writePtr+1] = (len(ploadbytes) >> 8) & 0xff
+    buf[writePtr+2] = len(ploadbytes) & 0xff
+    buf[writePtr+3:writePtr+3+len(ploadbytes)] = ploadbytes
+    writePtr += len(ploadbytes)
+    buf[writePtr:writePtr+len(tagbytes)] = tagbytes
+    pack = n
+    pack = DatagramPacket(buf, len(buf), InetAddress.getByName("localhost"), 4445)
+    sock.send(pack)
+
+
+------------
+
+broadcast log receiver
++ needs a logger
+listener ist interface for receiver
+
+gui wall (SYSLAB Userspacce)
+broadcast log displet
+https://git.elektro.dtu.dk/syslab/syslab-userspace/-/blob/master/src/main/java/risoe/syslab/gui/wall/displets/BroadcastLogDisplet.java
+... maybe extend with simple log file writer. 
+
+
+
+

syslab-python/requirements.txt

@@ -0,0 +1,2 @@
+requests >= 2.18
+beautifulsoup4 >= 3.7

syslab-python/setup.cfg

@@ -0,0 +1,10 @@
+[flake8]
+ignore =
+max-line-length = 79
+max-complexity = 11
+
+[pytest]
+addopts = --doctest-glob="*.rst"
+
+[wheel]
+universal = True

syslab-python/setup.py

@@ -0,0 +1,36 @@
+from setuptools import setup, find_packages
+
+
+setup(
+    name='syslab',
+    version='0.3.0',
+    author='Tue Vissing Jensen',
+    author_email='tvjens at elektro.dtu.dk',
+    description=('SYSLAB webservice client library.'),
+    long_description=(''),
+    url='https://www.syslab.dk',
+    install_requires=[
+        'requests>=2.18',
+        'beautifulsoup4>=3.7',
+    ],
+    packages=find_packages(exclude=['tests*']),
+    include_package_data=True,
+    entry_points={
+        'console_scripts': [
+        ],
+    },
+    classifiers=[
+        'Development Status :: 3 - Alpha',
+        'Environment :: Console',
+        'Intended Audience :: Science/Research',
+        'License :: Other/Proprietary License',
+        'Natural Language :: English',
+        'Operating System :: OS Independent',
+        'Programming Language :: Python',
+        'Programming Language :: Python :: 3',
+        'Programming Language :: Python :: 3.7',
+        'Programming Language :: Python :: 3.8',
+        'Topic :: Scientific/Engineering',
+        'Topic :: Software Development :: Libraries :: Python Modules',
+    ],
+)

syslab-python/syslab/core/SyslabUnit.py

@@ -231,7 +231,7 @@ def parse_resources(url):
     for resource in resources:
         try:
             # TODO - handle java arrays, ie, []
-            m = re.match("(\w+)\[?\]? (\w+)(\([\w ,]*\))?", resource).groups()
+            m = re.match(r"(\w+)\[?\]? (\w+)(\([\w ,]*\))?", resource).groups()
         except AttributeError as e:
             print(e)
             continue

util.py

@@ -1,42 +1,5 @@
-import numpy as np
-from syslab.core.datatypes import CompositeMeasurement
-from typing import Union
-from time import time
-
-
-def pos(x):
-    """
-    :param x: input
-    :return: x if x > 0 else 0
-    """
-    return max(x, 0)
-
-
 def clamp(a, x, b):
     """
         Restrict x to lie in the range [a, b]
     """
     return max(a, min(x, b))
-
-def cast_to_cm(m: Union[CompositeMeasurement, float]):
-    if type(m) == float:
-        request = CompositeMeasurement(m, timestampMicros=time()*1e6, timePrecision=1000)
-    elif type(m) == CompositeMeasurement:
-        request = m
-    else:
-        raise TypeError(f"Unknown request type: {type(m)}")
-    return request
-
-def soc_scaler(soc, min_soc=0.0, max_soc=1.0):
-    """
-    Scale an soc to emulate having a smaller battery.
-    If the actual battery has a capacity of 14 kWh, the rescaled battery will have a capacity
-    of 14*(max_soc - min_soc) kWh.
-    Does not check for negative soc.
-    :param soc: current actual soc
-    :param min_soc: actual soc that should correspond to an soc of 0.0
-    :param max_soc: actual soc that should correspond to an soc of 1.0
-    :return: rescaled soc
-    """
-
-    return (soc - min_soc)/(max_soc - min_soc)