How to generate loadshape clusters from AMI data.

Daily loadshape identification is a common problem in utility data analysis. Before the advent of AMI metering, end-use metering surveys and simulation were often used to develop diurnal loadshapes for different seasons and weekday types. These surveys usually involved relative modest numbers of customers, whose customer class was usually well-known, and the loadshapes were often one of primary deliverables. Consequently generating the range of loadshapes observed in the surveys was not typically a problem for utility analysts.

With the advent of automated metering infrastructure (AMI), the volume of data increased significantly. Unfortunately, the analysis methods used for the end-use load surveys are not always helpful because the AMI data is not labeled. To address this problem, unsupervised machine learning methods can be used to identify the principal loadshapes found in the AMI data.

This post describes the methodology use by the OpenFIDO loadshape pipeline. The use-case presented here is based on a data set from a utility containing 1 year of AMI data from 2015 for 973 customers. The original dataset has been sanitized to preserve privacy and reduce the size of the data repository. In addition the subhourly data has been rolled up to hourly data.

The methodology used to create the loadshape data is based on an unsupervised learning algorithm known as k-means clustering which looks for $k$ loadshapes in the data. This method partitions $N$ observations into $k$ clusters in which each observation belongs to the cluster with the nearest mean. Each cluster mean serves as a prototype of the cluster. This partition produces a data space divided into so-called “Voronoi” cells. The clusters minimize the loadshape variances (i.e., the Euclidean distances).

The analysis performed here applies the loadshape clustering method step by step, visualizing the intermediate results for inspection, and comparing the final result with an alternate t-distributed stochastic neighbor embedded ($t$-SNE) partitioning method. The $t$-SNE algorithm is a dimensionality reduction method that can be used to visualize high-dimensional data as a two or three-dimensional map that does not share any common elements with the $k$-means method being validated.

The entire process is divided into seven steps, briefly summarized as follows:

1. Load the AMI data - The original data was downloaded from a utility AMI system and stored in a CSV file. This file was cleaned of identifiable information and aggregated to hourly interval energy values before being stored for analysis.

2. Identify hour types - The hours of the year are grouped into a total of 192 hour types for each season, day type, and hour that comprise the loadshapes to be found in the data.

3. Aggregate hour types - The hourly meter data is aggregated by hour type to find the mean value for each meter during each hour type.

4. Cluster the meters - The meters are clustered using $k$-means into $K$ loadshapes and each meter is assigned to the loadshape with the smallest mean-squared difference.

5. Rescale the loadshapes - Each meter is assigned a loadshape scale and offset to ensure that the loadshape has the same variance as the meter which is assigned to it.

6. Generate the loadshape weights - Each loadshape is assigned a weight based on the fraction of the total annual energy associated with that loadshape.

7. Validate the results - The results are validated using a $t$-SNE partition method.

Load the AMI data

The first task is to load the AMI data from the sanitized AMI data repository. In this case we are only interested in the timestamps with timezones, the meter ids, and the power values. All other data in the repository is ignored.

import pandas as pd
import datetime as dt
pd.set_option("max_columns",11)
pd.set_option("max_rows",9)

data = pd.read_csv("ami.csv.gz",
    converters = {
        "timestamp" : dt.datetime.fromisoformat,
        "meter_id" : int,
        "tz" : int,
        "power" : float,
    })

data

The data set contains nearly 9 million records, with meter id, timestamps, timezone, and power measurements. The timestamps are provided in UTC and the timezone provides the offset of the local time from the UTC timestamp. These will be necessary to determine time-dependent information such as the season, day of week, and local hour of day.

Identify hour types

The timestamps are then analysed to identify the hour types that are used to generate the loadshapes. The hour type is based on the season (winter=0, spring=1, summer=2, and fall=3) and the day type (weekday=0 and weekend=1). These, along with the 24 hours of the day, result in 192 hour types, which are associated with each meter observation.

hour = (data["timestamp"].apply(lambda x:x.hour-1)+data["tz"]).astype(int).mod(24)
weekend = data["timestamp"].apply(lambda x:(int(x.weekday()/5)))
season = data["timestamp"].apply(lambda x:x.quarter-1)
data["hourtype"] = season*48 + weekend*24 + hour

data

Aggregate hour types

The next task is to calculate the average load for each hour type of each meter. The result is then pivoted to create one column for each hour type and one row for each meter.

groups = data.groupby(["meter_id","hourtype"]).mean().reset_index()
groups = groups.pivot(index="meter_id",columns="hourtype",values="power")
groups.fillna(method="ffill",inplace=True)

groups

Note that some meters may not have sufficient data to generate loads shapes. The following identifies which meters do not have loadshapes associated with them.

missing = data.iloc[list(set(data["meter_id"].unique()).difference(set(groups.index.values)))]
#missing.to_csv("missing.csv")
list(missing["meter_id"].unique())

[]

We can visualize the average load for each meter by hour type, where hours 0-23 are winter weekdays, hours 24-47 are winter weekends, hours 48-71 are spring weekdays, etc.

groups.T.plot(figsize=(15,8), legend=False, color='blue', alpha=0.05,grid=True,xlim=[0,191],xticks=range(0,192,24));

Clustering the meters

The next task is to find the clusters that partition this dataset among loadshapes with minimal Euclidean distances to the cluster centroid. We use the $k$-means algorithm to identify 10 clusters on the normalized hour type profiles and search for the smallest number of clusters necessary to ensure every group has more than 1 loadshape in it.

from sklearn.cluster import KMeans
from sklearn.preprocessing import MinMaxScaler
rescaled = MinMaxScaler().fit_transform(groups.values.copy())
ok = False
n_clusters = 10
while not ok and n_clusters > 1:
    n_clusters -= 1
    groups = data.groupby(["meter_id","hourtype"]).mean().reset_index()
    groups = groups.pivot(index="meter_id",columns="hourtype",values="power")
    groups.fillna(method="ffill",inplace=True)    
    kmeans = KMeans(n_clusters=n_clusters)
    group_found = kmeans.fit_predict(rescaled)
    group_data = pd.Series(group_found, name="loadshape")
    groups.set_index(group_data, append=True, inplace=True)
    ok = True
    for group in groups.index.get_level_values(level=1).unique():
        if len(groups.xs(group, level=1)) == 1:
            ok = False

A histogram of the loadshapes found shows that only a few loadshapes are common and the remaining ones are idiosyncratic.

import matplotlib.pyplot as plt
plt.figure(figsize=(10,10))
ax = groups.reset_index().groupby("loadshape")["meter_id"].count().plot.pie()
ax.set_title("Loadshape customer count")
ax.set_ylabel("");

The loadshape data is obtained by aggregating the median load and relabeling the columns to improve readability.

loadshapes = groups.groupby("loadshape").mean()
seasons = ["win","spr","sum","fal"]
weekdays = ["wd","we"]
loadshapes.columns = [f"{seasons[season]}_{weekdays[weekend]}_{hour}h" for season in range(4) for weekend in range(2) for hour in range(24)]
loadshapes

The final results are plotted with the original data to visualize how the loadshapes are representative of each cluster.

rows = int(loadshapes.shape[0]/2+0.5)
fig,ax = plt.subplots(rows,2, figsize=(18,5*rows))
for group in groups.index.get_level_values(level=1).unique():
    subset = groups.xs(group, level=1)
    ax = plt.subplot(rows,2,group+1)
    plt.plot(subset.T,alpha=0.05,ls='-',color="blue")
    plt.plot(subset.median(), alpha=1.0, ls='-',color="black")
    ax.set_title(f"Loadshape {group} (N={len(subset)})")
    ax.set_ylabel('Power [kW]')
    ax.set_xlabel('\nHour type')
    ax.grid()
    y0 = ax.get_ylim()[0]
    for s in range(4):
        for w in range(2):
            label = plt.text(s*48+w*24+12,y0,["Weekday","Weekend"][w],axes=ax)
            label.set_ha("center")
            label.set_va("top")
        label = plt.text(s*48+24,y0,["\nWinter","\nSpring","\nSummer","\nFall"][s],axes=ax)    
        label.set_ha("center")
        label.set_va("top")
    plt.xticks(range(0,192,24),[],axes=ax)
    plt.xlim([0,191])            

Rescaling the loadshapes

When the identified loadshapes are used to generate a load model, the mean loadshape power (as well as the total energy) will always match the original meters’ means and total energy because the methodology uses the mean power to identify the cluster centroids. However, the methodology does not address the different variances among the meters and these will not match the loadshapes centroids’ variances.

To address this, we compute the standard deviation for each meter power reading and its corresponding loadshape value.

values = groups.reset_index().set_index(["meter_id","loadshape"]).melt(ignore_index=False).reset_index().set_index(["meter_id","loadshape","hourtype"])
meters = data.drop(["tz","timestamp"],axis=1).set_index(["meter_id","hourtype"])
mapping = meters.join(values)
mapping.groupby("meter_id").std()

To rescale the loadshape so that each meter’s loadshape has the same variance we use a linear function for each meter $m$

where the scale $scale_m = { \sigma^2_m \over \sigma^2_l}$ with $\sigma^2_m$ is the variance of the meter data and $\sigma^2_l$ is the variance corresponding loadshape data, and $offset_m = (1-scale_m)~\mu_m$ with $\mu_m$ being the mean of the meter data.

The histograms of the scales and offsets verifies the reasonableness of the range of values found.

fig,ax = plt.subplots(1,2,figsize=(18,5))
scale = pd.DataFrame(mapping.groupby("meter_id").std().power/mapping.groupby("meter_id").std().value,columns=["scale"])
scale.hist(figsize=(10,5),density=True,ax=ax[0]);
offset = pd.DataFrame((1-scale).scale*mapping.groupby("meter_id").mean()["power"],columns=["offset"])
offset.hist(figsize=(10,5),density=True,ax=ax[1]);

The scales and offsets are now joined with the original meter data to compute the rescaled values.

mapping = mapping.join(scale)
mapping = mapping.join(offset)
mapping["rescaled"] = mapping.value*mapping.scale + mapping.offset
mapping

The final rescaled loadshape values have the same variance, mean power, and total energy use as the original AMI meter data:

final=mapping.drop(["scale","offset","value"],axis=1)
pd.DataFrame(list(map(lambda x:(x.power-x.rescaled)*(x.power-x.rescaled).max(),[final.std(),final.mean(),final.sum()])),index=["std","mean","sum"],columns=["MaxSE"]).T

Generate the loadshape weights

The loadshape weights can be generated by computing the total energy in each loadshape meter group and dividing it by the total energy over all loadshapes.

loadshapes = groups.drop(range(192),axis=1)
total = data.drop("tz",axis=1).set_index(["meter_id","hourtype"]).join(loadshapes).groupby("loadshape").sum()
total.columns = ["energy"]
(total/total.sum()).plot.pie(y="energy",figsize=(8,8),normalize=True);
(total/total.sum()).round(3).T

Validation

The final result is validated using a $t$-distributed stochastic neighbor embedding ($t$-SNE) algorithm to reduce the dimensionality of the original data set to two dimensions. The color of each point in the map corresponds to each loadshape. The grouping of similar loadshapes shows how the cluster methods produce similar results.

from sklearn.manifold import TSNE
import matplotlib.colors
import numpy as np
tsne = TSNE()
results_tsne = tsne.fit_transform(rescaled)
color_list = list(zip(np.arange(0,1,0.1),np.arange(0,1,0.1),np.arange(1,0,-0.1)))
cmap = matplotlib.colors.LinearSegmentedColormap.from_list(groups.index.get_level_values(level=1).unique(), color_list)
plt.figure(figsize=(10,10))
plt.scatter(results_tsne[:,0], results_tsne[:,1],
    c=groups.index.get_level_values(level=1),
    cmap=cmap, 
    alpha=0.6, 
    );

This provides an simple example from which most analysts and researchers can develop their own loadshapes from AMI data. In some cases it may be necessary to do some preprocessing to the AMI data to normalize timestamps, aggregate subhourly data up to hourly day, or removed extraneous columns. These tasks are relatively simple and straightforward, and are common enough that the reader should have no problem implementing these for themselves.