PartialDischargeCorrelator and PartialDischargeForecaster

CITATION.cff

@@ -0,0 +1,22 @@
+cff-version: 1.2.0
+message: "If you use this software, please cite it as below."
+authors:
+  - family-names: van Dinter
+    given-names: Raymon
+    orcid:  https://orcid.org/0000-0002-1811-8803
+  - family-names: Ekmekci
+    given-names: Görkem
+  - family-names: Netten
+    given-names: Gerdtinus
+  - family-names: Rieken
+    given-names: Sander
+  - family-names: Tekinderdogan
+    given-names: Bedir
+    orcid:  https://orcid.org/0000-0002-8538-7261
+  - family-names: Catal
+    given-names: Cagatay
+    orcid:  https://orcid.org/0000-0003-0959-2930
+title: "A code repository for predictive maintenance on cable joints."
+version: v0.1
+doi: -
+date-released: 2023-04-25

alliander_predictive_maintenance/cognition/visualizer.py

@@ -0,0 +1,98 @@
+import matplotlib.pyplot as plt
+import seaborn as sns
+
+from alliander_predictive_maintenance.conversion.partial_discharge_correlator.partial_discharge_weather_correlator_results import \
+    PartialDischargeWeatherCorrelatorResults
+from alliander_predictive_maintenance.conversion.partial_discharge_correlator.weather_categories import WeatherCategories
+from alliander_predictive_maintenance.cyber.simulation_model.partial_discharge_forecaster_model_results import PartialDischargeForecasterModelResults
+
+
+class Visualizer:
+    """ A class for visualizing results of several models. """
+
+    def plot_forecasted_partial_discharge(self, results: PartialDischargeForecasterModelResults, axes: plt.Axes,
+                                          title: str, ylabel: str, xlabel: str):
+        """ Plot PartialDischargeForecasterModelResults
+
+        :param results: PartialDischargeForecasterModelResults
+        :param axes: axes for subplots
+        :param title: title of the subplot
+        :param ylabel: y label of the subplot
+        :param xlabel: x label of the subplot
+        """
+        self.__plot_partial_discharge_forecaster_model_results(results, "ground truth", "predictions", "red", axes)
+        self.__set_title_and_labels(axes, title, ylabel, xlabel)
+
+    def plot_forecasted_train_partial_discharge(self, train_results: PartialDischargeForecasterModelResults,
+                                                test_results: PartialDischargeForecasterModelResults,
+                                                axes: plt.Axes, title: str, ylabel: str, xlabel: str):
+        """ Plot PartialDischargeForecasterModelResults for after fitting the model
+
+        :param train_results: PartialDischargeForecasterModelResults for the train set
+        :param test_results: PartialDischargeForecasterModelResults for the test set
+        :param axes: axes for subplots
+        :param title: title of the subplot
+        :param ylabel: y label of the subplot
+        :param xlabel: x label of the subplot
+        """
+        self.__plot_partial_discharge_forecaster_model_results(train_results, "train", "train predictions", "green",
+                                                               axes)
+        self.__plot_partial_discharge_forecaster_model_results(test_results, "test", "test predictions", "red", axes)
+        self.__set_title_and_labels(axes, title, ylabel, xlabel)
+
+    def __plot_partial_discharge_forecaster_model_results(self, results: PartialDischargeForecasterModelResults,
+                                                          ground_truth_label: str, predict_label: str, color: str,
+                                                          axes: plt.Axes):
+        """ Plot partial discharge forecaster model results
+
+        :param results: PartialDischargeForecasterModelResults
+        :param ground_truth_label: string for legend
+        :param predict_label: string for legend
+        :param color: color of predicted values
+        :param axes: axes for subplot
+        """
+        axes.plot(results.true_partial_discharge, label=ground_truth_label)
+        axes.plot(results.true_partial_discharge.index, results.partial_discharge, label=predict_label, c=color)
+
+    def __set_title_and_labels(self, axes: plt.Axes, title: str, ylabel: str, xlabel: str):
+        """ Set title and labels of subplots
+
+        :param axes: axes for subplots
+        :param title: title of the subplot
+        :param ylabel: y label of the subplot
+        :param xlabel: x label of the subplot
+        """
+        axes.set_title(title)
+        axes.set_ylabel(ylabel)
+        axes.set_xlabel(xlabel)
+        axes.legend()
+
+    def plot_partial_discharge_weather_correlation(self, correlator_results: PartialDischargeWeatherCorrelatorResults,
+                                                   axes: plt.Axes, title: str, ylabel: str, xlabel: str,
+                                                   correlation_coefficient_threshold: float):
+        """ Plot model results of the PartialDischargeWeatherCorrelator
+
+        :param correlator_results: PartialDischargeWeatherCorrelatorResults
+        :param axes: axes for subplots
+        :param title: title of the subplot
+        :param ylabel: y label of the subplot
+        :param xlabel: x label of the subplot
+        :param correlation_coefficient_threshold: threshold for Pearson correlation coefficient
+        """
+        for weather_category in WeatherCategories:
+            correlating_features = correlator_results.partial_discharge_correlating_weather_features(
+                weather_category, correlation_coefficient_threshold)
+            if len(correlating_features) > 0:
+                print(f"Correlates with {weather_category.name}: {correlating_features[-1]:.2f} for feature "
+                      f"{correlating_features.index[-1][1]}")
+
+        # plot a heatmap of pd features vs weather features
+        corr = correlator_results.correlations
+        pd_columns = [column for column in corr.columns if 'partial_discharge' in column]
+        corr_filtered = corr[~corr.index.isin(pd_columns)]
+        corr_filtered = corr_filtered[pd_columns]
+        significant_features = corr_filtered.apply(lambda row: any(row.abs() > correlation_coefficient_threshold),
+                                                   axis=1)
+        corr_filtered = corr_filtered.loc[significant_features[significant_features].index]
+        sns.heatmap(corr_filtered, annot=True, cmap=plt.cm.PuBu, ax=axes)
+        self.__set_title_and_labels(axes, title, ylabel, xlabel)

alliander_predictive_maintenance/connection/readers/circuit/circuit_factory.py

@@ -1,4 +1,3 @@
-from pathlib import Path
 from typing import Dict
 
 import pandas as pd
@@ -17,12 +16,14 @@ class CircuitFactory:
     def __init__(self,
                  config: Dict,
                  circuit_coordinates_reader: ICircuitCoordinatesReader,
-                 weather_retriever: ICircuitWeatherRetriever,
+                 cds_weather_retriever: ICircuitWeatherRetriever,
+                 knmi_weather_retriever: ICircuitWeatherRetriever,
                  csv_partial_discharge_storage: ICircuitPartialDischargeReader,
                  csv_circuit_config_reader: ICircuitConfigReader):
         self.__config = config
         self.__circuit_coordinates_reader = circuit_coordinates_reader
-        self.__weather_retriever = weather_retriever
+        self.__cds_weather_retriever = cds_weather_retriever
+        self.__knmi_weather_retriever = knmi_weather_retriever
         self.__csv_partial_discharge_storage = csv_partial_discharge_storage
         self.__csv_circuit_config_reader = csv_circuit_config_reader
         self._circuit_weather = None
@@ -39,9 +40,11 @@ def create_circuit(self, circuit_id: int) -> Circuit:
         partial_discharge_data = self.__load_partial_discharge_data(str(circuit_id))
         time_window = TimeWindow(partial_discharge_data[ICircuitPartialDischargeReader.DATETIME_COLUMN].min(),
                                  partial_discharge_data[ICircuitPartialDischargeReader.DATETIME_COLUMN].max())
-        circuit_weather = self.__weather_retriever.get_weather(circuit_coordinate, time_window)
+        circuit_knmi_weather = self.__knmi_weather_retriever.get_weather(circuit_coordinate, time_window)
+        circuit_cds_weather = self.__cds_weather_retriever.get_weather(circuit_coordinate, time_window)
         return Circuit(circuit_id=circuit_id,
-                       weather=circuit_weather,
+                       cds_weather=circuit_cds_weather,
+                       knmi_weather=circuit_knmi_weather,
                        circuit_coordinate=circuit_coordinate,
                        partial_discharge=partial_discharge_data,
                        time_window=time_window,

alliander_predictive_maintenance/constants.py

@@ -18,4 +18,7 @@
 HTTP_ERROR_COULD_NOT_GET_DATA = "Could not get data from the API"
 INVALID_JOINT_LOCATION = "Joint location {location} is invalid with circuit length {circuit_length}"
 INVALID_TIME_WINDOW = "TimeWindow {time_window} is out of range for {min}, {max}"
+WEATHER_DATA_MANDATORY = "Weather data does not contain all mandatory columns {data} not in {mandatory}"
+WEATHER_DATA_DATETIME_INDEX = "Weather data index is not in DateTime format"
+PARTIAL_DISCHARGE_DATA_DATETIME_INDEX = "Partial discharge data index is not in DateTime format"
 NOTIMPLEMENTEDERROR_AWS = "Alliander S3 environment should load {data} data here"

alliander_predictive_maintenance/conversion/data_types/circuit.py

@@ -12,10 +12,11 @@
 
 class Circuit:
     """ An object representing a Circuit consisting of joints """
-    def __init__(self, circuit_id: int, weather: pd.DataFrame, circuit_coordinate: CircuitCoordinate,
+    def __init__(self, circuit_id: int, cds_weather: pd.DataFrame, knmi_weather: pd.DataFrame, circuit_coordinate: CircuitCoordinate,
                  partial_discharge: pd.DataFrame, time_window: TimeWindow, circuit_length: float):
         self.__circuit_id = circuit_id
-        self.__weather = weather
+        self.__cds_weather = cds_weather
+        self.__knmi_weather = knmi_weather
         self.__circuit_coordinate = circuit_coordinate
         self.__partial_discharge = partial_discharge
         self.__time_window = time_window
@@ -24,10 +25,18 @@ def __init__(self, circuit_id: int, weather: pd.DataFrame, circuit_coordinate: C
     @property
     def circuit_id(self):
         return self.__circuit_id
+
+    @property
+    def circuit_length(self):
+        return self.__circuit_length
+
+    @property
+    def cds_weather(self) -> pd.DataFrame:
+        return self.__cds_weather
 
     @property
-    def weather(self) -> pd.DataFrame:
-        return self.__weather
+    def knmi_weather(self) -> pd.DataFrame:
+        return self.__knmi_weather
 
     @property
     def time_window(self):
@@ -46,10 +55,10 @@ def create_joint(self, location: float, time_window: TimeWindow,
             raise ValueError(INVALID_JOINT_LOCATION.format(location=location, circuit_length=self.__circuit_length))
         partial_discharge = self.__get_partial_discharge_at_location(location=location,
                                                                      resampling_strategy=resampling_strategy)
-        if time_window.start_date < partial_discharge.index.min() or time_window.end_date > partial_discharge.index.max():
+        if time_window.start_date > time_window.end_date:
             raise ValueError(INVALID_TIME_WINDOW.format(time_window=time_window,
-                                                        min=partial_discharge.index.min(),
-                                                        max=partial_discharge.index.max()))
+                                                        min=time_window.start_date,
+                                                        max=time_window.end_date))
         partial_discharge = partial_discharge[time_window.start_date:time_window.end_date]
         return Joint(location, partial_discharge)
 

alliander_predictive_maintenance/conversion/data_types/predictive_maintenance_joint_data.py

@@ -0,0 +1,14 @@
+import pandas as pd
+from dataclasses import dataclass
+from typing import Union, Optional
+
+
+@dataclass
+class PredictiveMaintenanceJointData:
+    """ Data needed for predictive maintenance on a cable joint.
+    weather: historical weather data. partial_discharge: historical partial discharge data."""
+    cds_weather: pd.DataFrame
+    knmi_weather: Optional[pd.DataFrame]
+    partial_discharge: Union[pd.DataFrame, pd.Series]
+    circuit_id: int
+    location_in_meters: float

alliander_predictive_maintenance/conversion/partial_discharge_correlator/partial_discharge_weather_correlator.py

@@ -0,0 +1,41 @@
+import numpy as np
+
+from alliander_predictive_maintenance.conversion.data_types.predictive_maintenance_joint_data import PredictiveMaintenanceJointData
+from alliander_predictive_maintenance.conversion.partial_discharge_correlator.partial_discharge_weather_correlator_results import \
+    PartialDischargeWeatherCorrelatorResults
+
+
+class PartialDischargeWeatherCorrelator:
+    """ A class for calculating correlations between partial discharge and weather data. """
+    PARTIAL_DISCHARGE_COLUMN = "partial_discharge"
+
+    def correlate(self, predictive_maintenance_joint_data: PredictiveMaintenanceJointData,
+                  rolling_days: int) -> PartialDischargeWeatherCorrelatorResults:
+        """ Calculate Pearson Correlation Coefficient for Partial Discharge data and weather data
+
+        :param predictive_maintenance_joint_data: PredictiveMaintenanceJointData object
+        :param rolling_days: number of days to use for a moving average
+        :return: PartialDischargeWeatherCorrelatorResults
+        """
+        predictive_maintenance_joint_data.knmi_weather.drop(['lat', 'lon'], axis=1, inplace=True, errors='ignore')
+        predictive_maintenance_joint_data.cds_weather.drop(['lat', 'lon', 'is_permanent_data', 'mean_wave_direction'],
+                                                           axis=1, inplace=True,
+                                                           errors='ignore')
+        data_frame = predictive_maintenance_joint_data.cds_weather.set_index("time").join(
+            predictive_maintenance_joint_data.partial_discharge.rename(self.PARTIAL_DISCHARGE_COLUMN))
+
+        data_frame_resampled = data_frame.resample("1d").agg([np.sum, np.median, np.mean, np.min, np.max])
+        data_frame_resampled['partial_discharge_cumsum'] = data_frame_resampled.partial_discharge['sum'].cumsum()
+        data_frame_resampled['partial_discharge_cumsum_gradient'] = np.gradient(
+            data_frame_resampled['partial_discharge_cumsum'])
+
+        data_frame_resampled.columns = ['_'.join(col).strip('_') for col in data_frame_resampled.columns]
+        data_frame_resampled = data_frame_resampled.join(
+            predictive_maintenance_joint_data.knmi_weather.set_index("time"))
+
+        data_frame_rolling = data_frame_resampled.rolling(rolling_days).mean()
+
+        corr_matrix = data_frame_rolling.corr()
+        corr_matrix.fillna(0, inplace=True)
+
+        return PartialDischargeWeatherCorrelatorResults(corr_matrix)

alliander_predictive_maintenance/conversion/partial_discharge_correlator/partial_discharge_weather_correlator_results.py

@@ -0,0 +1,38 @@
+import pandas as pd
+from dataclasses import dataclass
+
+from alliander_predictive_maintenance.conversion.partial_discharge_correlator.weather_categories import WeatherCategories
+
+
+@dataclass
+class PartialDischargeWeatherCorrelatorResults:
+    """ A data class that contains partial discharge and weather correlation results. """
+    correlations: pd.Series
+
+    @property
+    def sorted_correlation_series(self) -> pd.Series:
+        """ A property returning a sorted series of correlations.
+        The first level of MultiIndex is partial_discharge-related, the second level is weather-related."""
+        corr_matrix_unstacked = self.correlations.unstack()
+        sorted_corr_series = corr_matrix_unstacked.sort_values(key=abs)
+
+        pd_columns = [column for column in self.correlations.columns if 'partial_discharge' in column]
+
+        sorted_corr_series = sorted_corr_series[sorted_corr_series.index.get_level_values(0) !=
+                                                sorted_corr_series.index.get_level_values(1)]
+        sorted_corr_series = sorted_corr_series[sorted_corr_series.index.get_level_values(0).isin(pd_columns)]
+        sorted_corr_series = sorted_corr_series[~sorted_corr_series.index.get_level_values(1).isin(pd_columns)]
+        return sorted_corr_series
+
+    def partial_discharge_correlating_weather_features(self, weather_category: WeatherCategories,
+                                                       correlation_coefficient_threshold: float) -> pd.Series:
+        """ Find weather features in a weather category that correlate with partial discharge
+        :param weather_category: WeatherCategories attribute
+        :param correlation_coefficient_threshold: Pearson Correlation Coefficient threshold
+        :return: series of correlations over the threshold for a specific weather category
+        """
+        sorted_corr_series = self.sorted_correlation_series.copy()
+        sorted_corr_series = sorted_corr_series[sorted_corr_series > correlation_coefficient_threshold]
+        sorted_corr_series = sorted_corr_series[
+            sorted_corr_series.index.get_level_values(1).isin(weather_category.value)]
+        return sorted_corr_series

alliander_predictive_maintenance/conversion/partial_discharge_correlator/weather_categories.py

@@ -0,0 +1,25 @@
+from enum import Enum
+
+
+class WeatherCategories(Enum):
+    """ Categories of Weather data from KNMI Daggegevens and CDS Era5sl"""
+
+    WIND = ['wind_direction', 'FHVEC', 'wind_speed', 'wind_speed_max', 'wind_speed_max_hour', 'wind_speed_min',
+            'wind_speed_min_hour', 'wind_gust_max', 'wind_gust_max_hour', '100m_u_component_of_wind',
+            '100m_v_component_of_wind', '10m_u_component_of_wind', '10m_v_component_of_wind']
+    TEMPERATURE = ['temperature', 'temperature_min', 'temperature_min_hour', 'temperature_max', 'temperature_max_hour',
+                   'T10N', 'T10NH', '2m_temperature']
+    SOIL_TEMPERATURE = ['soil_temperature_level_1', 'soil_temperature_level_2', 'soil_temperature_level_3',
+                        'soil_temperature_level_4']
+    PRECIPITATION = ['precipitation_duration', 'precipitation', 'precipitation_max', 'precipitation_max_hour',
+                     'total_precipitation']
+    HUMIDITY = ['humidity', 'humidity_max', 'humidity_max_hour', 'humidity_min', 'humidity_min_hour']
+    SOLAR = ['sunlight_duration', 'percentage_of_max_possible_sunlight_duration', 'surface_solar_radiation_downwards',
+             'surface_solar_radiation_downward_clear_sky', 'global_radiation']
+    AIR_PRESSURE = ['air_pressure', 'air_pressure_max', 'air_pressure_max_hour', 'air_pressure_min',
+                    'air_pressure_min_hour']
+    VISION = ['VVN', 'VVNH', 'VVX', 'VVXH']
+    WATER_IN_SOIL = ['volumetric_soil_water_layer_1', 'volumetric_soil_water_layer_2', 'volumetric_soil_water_layer_3',
+                     'volumetric_soil_water_layer_4']
+    OTHER = ['EV24', 'cloud_cover', '2m_dewpoint_temperature', 'mean_sea_level_pressure', 'mean_wave_period',
+             'surface_pressure', 'sea_surface_temperature', 'significant_height_of_combined_wind_waves_and_swell']