Machine learning for BMS analysis and optimisation

In [29], the problem of estimating sparse graphs, graphs with only few edges, by a lasso penalty to the inverse covariance matrix is considered. Let us consider the case where X₁, ..., X_n are independent and identically distributed N_p (μ, Σ) and being the estimated precision matrix, which is the inverse matrix from the covariance matrix, denoted Ω ≡ Σ⁻¹. This function solves the following optimisation problem:

$$\begin{eqnarray}&&\mathrm{Objective}:\quad {\hat{{\rm{\Omega }}}}_{\lambda }={{argmin}}_{{\rm{\Omega }}\in {S}_{+}^{P}}\{{Tr}(S{\rm{\Omega }})-\mathrm{log}{\det }({\rm{\Omega }})+| | {\rm{\Omega }}| {| }_{1}\}\end{eqnarray} \tag{ 1 }$$

where tr is the trace, λ > 0 and the penalisation parameter $| | {\rm{\Omega }}| {| }_{1}={\sum }_{i,j}| {{\rm{\Omega }}}_{i,j}|$ is the L₁-Norm of Σ (Lasso regularisation parameter). In order to provide useful information, this problem proposes to maximise the penalized log likelihood with respect to Ω, so the nodes are not fully connected and the connections kept on each cluster have useful information concerning the relationships of time series on each cluster. This happens because Lasso regularisation parameter shrinks the less important features coefficient to zero, removing less meaningful coefficients. The algorithm employed to solve this problem is the GLasso algorithm, which is explained in [29], where they consider the problem of estimating sparse graphs by a lasso penalty, which is the penalty applied to non-zero coefficients by the sum of their absolute values, applied to the inverse covariance matrix.

With this estimate of the inverse of the correlation matrix, we have the partial independence relationship. If two features are independent conditionally on the others, the corresponding coefficient in the inverse of the covariance matrix would be zero, as it learns independence relations from the data, instead of being a distance measure itself between time series.

According to the authors of the implementation package in 'scikit-learn' in [30], the search for the optimal penalization parameter is done on an iteratively refined grid: first the cross-validated scores on a grid are computed, then a new refined grid is centred around the maximum, and so on. One of the challenges here is that the solvers can fail to converge to a well-conditioned estimate. The corresponding values of alpha then come out as missing values, but the optimum may be close to these missing values.

Agglomerative Hierarchical Clustering (AHC), has a long history, especially in taxonomy or classificatory systems, and phylogenetics [31, 32]. Further studies generalised this algorithm, [33], and have further developed and improved [34, 35].

Base on the definition given in [36], the goal of hierarchical clustering is to create a sequence of nested partitions or clusters, which can be conveniently visualised via a tree or hierarchy of clusters, also called the cluster dendrogram. In AHC, it starts with each of the n points in a separate cluster, and then merging the two closest clusters until all points are members of the same cluster. Algorithm 1 shows this procedure, with ${\bf{D}}=\{{{\bf{x}}}_{1},...,{{\bf{x}}}_{n}\}$, where ${{\bf{x}}}_{i}\in {{\mathbb{R}}}^{d}$, a clustering ζ = {C₁, ..., C_k } is a partition of ${\bf{D}}$.

Algorithm 1. Agglomerative Hierarchical Clustering (D, k), reproduced with permission from [36]:

1 $\zeta \leftarrow \{{C}_{i}={{\bf{x}}}_{i}\}\ | {{\bf{x}}}_{i}\in {\bf{D}}|$;               // Each time series is in a separate cluster initially

2 ${\rm{\Delta }}\leftarrow \{\delta ({{\bf{x}}}_{i},{{\bf{x}}}_{j}):\ {{\bf{x}}}_{i},{{\bf{x}}}_{j}\in {\bf{D}}\}$;                                           // Compute matrix with distances

3 repeat

4     Find the closest pair of clusters Ci , Cj ∈ ζ ;

5     ${C}_{{ij}}\leftarrow {C}_{i}\cup {C}_{j}$;                                                 // Merge the clusters

6     $\zeta \leftarrow (\zeta \ \setminus \ \{{C}_{i},{C}_{j}\ \})\cup \{{C}_{{ij}}\ \}$;                                              // Update the clustering

7     Update distance matrix Δ to reflect new clustering;

8 until $| \zeta | =k$

Different distance measures between time series can be used for clustering: pearson's correlation coefficient, dynamic time warping and integrated periodogram distance.

2.2.1. Pearson's correlation similarity

Let $x={[{x}_{1}{x}_{2}...{x}_{L}]}^{T}$ and $y={[{y}_{1}{y}_{2}...{y}_{L}]}^{T}$ be two zero-mean real-valued random vectors of length L. As described in [37], the Pearson's correlation coefficient between x and y is

$$\begin{eqnarray}&&{\rho }^{2}(x,y)=\displaystyle \frac{{E}^{2}({x}^{T}y)}{E({x}^{T}x)E({y}^{T}y)}.\end{eqnarray} \tag{ 2 }$$

With E being the expected value. According to the 'tsclust' package documentation in [38], which is the one used for the purpose of this study, two different measures of dissimilarity between two time series based on the estimated Pearson's correlation can be computed. These can be ${d}_{1}=\sqrt{2(1-\rho )}$ or ${d}_{2}=\sqrt{{\left(\tfrac{1-\rho }{1+\rho }\right)}^{\beta }}$, where β specifies the regulation of the convergence.

2.2.2. Dynamic Time Warping

DTW has the basic idea behind that the sequences are extended by repeating elements and the distance is calculated between extended sequences. Therefore, DTW can handle input sequences with different lengths [39].

Let $x={[{x}_{1}{x}_{2}...{x}_{r}]}^{T}$ and $y={[{y}_{1}{y}_{2}...{y}_{s}]}^{T}$ be two time series, where lengths r and s are not necessarily equal. Let M be an r × s matrix with the (i, j) element containing the squared Euclidean distance between two points x_i and y_j . Each element (i, j) in M corresponds to the alignment between two points x_i and y_j . Let $W={w}_{1},{w}_{2},...,{w}_{k}$ be a warping path, where the kth element w_k  = (i_k , j_k ). Then $\max r,s\leqslant K\lt r+s-1$ with the warping paths having the following restrictions: monotonicity, continuity and boundary conditions. There are exponentially many paths that satisfy these conditions, being the optimal path the one which minimizes the warping cost [40]:

$$\begin{eqnarray}&&{DTW}(x,y)=\displaystyle \frac{\sqrt{{\displaystyle \sum }_{l=1}^{K}{w}_{l}}}{K}=\displaystyle \frac{\sqrt{{\displaystyle \sum }_{l=1}^{K}d({x}_{{i}_{l}},{y}_{{j}_{l}})}}{K},\end{eqnarray} \tag{ 3 }$$

where (i_l , j_l ) = w_l for l = 1, 2, ..., K. Then the optimal path can be found through dynamic programming according to $\gamma (i,j)=d({x}_{i},{y}_{i})+\min \gamma (i-1,j),\gamma (i-1,j-1),\gamma (i,j-1)$, where γ(i, j) is the dynamic time warping distance between the sub sequences x₁, x₂, ..., x_i and y₁, y₂, ..., y_j .

2.2.3. Integrated Periodogram Distance

The distance based on the normalised periodogram was introduced in [21]. Let ${P}_{x}({w}_{j})=(1/n)| {\sum }_{t=1}^{n}{x}_{t}{e}^{-{{itw}}_{j}}{| }^{2}$ and ${P}_{y}({w}_{j})=(1/n)| {\sum }_{t=1}^{n}{y}_{t}{e}^{-{{itw}}_{j}}{| }^{2}$ be periodograms of time series x and y, respectively, at frequencies ${w}_{j}=2\pi j/n,j=1,\ldots ,[n/2]$ in the range 0 to π, [n/2] being the largest integer less or equal to n/2. We are interested only on its correlation structure, so it is better to use normalized periodogram defined by ${NP}({w}_{j})=P({w}_{j})/\hat{\gamma }$, where $\hat{\gamma }$ is the sample variance of the time series. Also, since the variance of the periodogram ordinates is proportional to the spectrum value at the corresponding frequencies, logarithms can be taken and therefore, the distance between x and y can be defined by

$$\begin{eqnarray}&&{d}_{{LNP}}(x,y)=\sqrt{\displaystyle \sum _{j=1}^{[n/2]}{\left[\mathrm{log}({{NP}}_{x}({w}_{j}))-\mathrm{log}({{NP}}_{y}({w}_{j}))\right]}^{2}}.\end{eqnarray} \tag{ 4 }$$

Knowing that the periodogram has the equivalent representation $P({w}_{j})=2\left[{\hat{\gamma }}_{0}+{\displaystyle \sum }_{k=1}^{n-1}{\hat{\gamma }}_{k}{\cos }({w}_{j}k)\right]$, where ${\hat{\gamma }}_{k}$ is the sample autocovariance function (defined with more detail in [41]) which, according to the authors, leads to

$$\begin{eqnarray}&&{d}_{{NP}}(x,y)=\sqrt{16\left[\displaystyle \frac{1}{4}{\left({\hat{\rho }}_{1,x}-{\hat{\rho }}_{1,y}\right)}^{2}+...+\displaystyle \frac{1}{4}{\left({\hat{\rho }}_{n-1,x}-{\hat{\rho }}_{n-1,y}\right)}^{2}\right]}=2\sqrt{n}\sqrt{\displaystyle \sum _{k=1}^{n-1}{\left({\hat{\rho }}_{k,x}-{\hat{\rho }}_{k,y}\right)}^{2}}.\end{eqnarray} \tag{ 5 }$$

For the implementation of the AHC methodologies (Pearson's correlation similarity, dynamic time warping and integrated periodogram distance), R package TSclust has been used [38]. For lasso clustering, package Scikit-Learn from Python version 3.7 has been used [42].

The manufacturing part is comprised of several spaces, every space dedicated to a part of the process, as this is an automotive manufacturing plant. In terms of cooling, heating and ventilation systems, the building has 3 chillers, 6 boilers, 6 fan coil units and 16 multi-speed fan AHUs. As the focus of this paper are AHU systems, the actual linkage between AHU-space and type can be seen below in table 1. Some of the AHUs have been discarded due to the lack of sensors in their respective areas of influence, which makes them irrelevant to this study.

We based our experiments on ground truth table 1. The data has been extracted directly from the manufacturer for every AHU individually [43], technical specifications include also dimensions of boxes and specs. of motors, fans, thermal wheel (if applicable), coils, etc. The data set comprises BMS sensor records. Every sensor is tagged with a specific name describing the company and followed by building location and subset (Ventilation, Metering, Cooling, Heating, Globals, Terminals and Lighting).

The building used for this case study is a Rolls Royce plant located in Rotherham, with 1639 BMS points in total. The period chosen for this study is June to July 2018 (30 days). The facility comprises two main areas: the production plant on the ground floor and the offices in the upper floor.

Time series data provided by sensor points is not often reliable and there are data gaps to be filled, units to be removed, and anomalies, which are bits of data misplaced between points because of extraction. The used data points are:

• AHU Supply Air Temperature (SAT) points:  These points measure the temperature of the air supplied to the area. There is one of these points per AHU.
• AHU Return Air Temperature (RAT) points: These points measure the temperature of the air extracted from the space prior to re-circulation or disposal.
• Room temperature points: Sensors are located in each room to measure temperature. There are between 1 and 4 sensors located on each room, and consists of averages of the sensors.
• Fan Coil Unit (FCU) room temperature: Output temperature of FCUs. This piece of equipment is connected (or close to) an AHU, and we define the relationships in one of the experiments below.

The office spaces contains the meeting rooms and common spaces such as restaurant, canteen, kitchen, changing room and open spaces. In this experiment we establish correlations-based relationships of these two types of equipment. This correlations study would illustrate the way in which different equipment types are linked, so the main influence areas are clearly defined. A schematic is shown in figure 3. In this case study building, FCUs feed air to specific meeting rooms in the office above the manufacturing facilities.

Zoom In Zoom Out Reset image size

Different meeting rooms receive the following generic names: Rhenium, Tugsten, Nickel, Cobalt, Tantalum and the gym. The main goal of this experiment is to link the respective rooms, each being fed by a different FCU, with their closest AHU unit's area of influence, being the AHU's the same ones as in experiments 1 and 2.

In order to test if the correlations are reliable, we need to ensure their significance. If the link between room sensors and SAT is established, the internal difference in temperature is an important factor to take into consideration to define this significance. All spaces in the facility may generate heat internally (kitchen, manufacturing process, people, computer equipment, etc.) which could make difficult identification of spaces with their respective AHUs. This is why we measure this internal disruption in the first place, and then gradually add complexity to the clustering.

For this purpose, the difference between supply and return air temperature has been taken into account,

$$\begin{eqnarray}&&{\rm{\Delta }}T=| {SAT}-{RAT}| \ {(}^{\circ }C).\end{eqnarray} \tag{ 6 }$$

Table 2 represents the mean, variance, minimum and maximum value of the element-wise difference between supply and return air temperatures for each AHU.

Table 2 shows that some AHUs, such as 6, 7A, 7B, 9 and 10 present a much higher mean value of their difference in supply and return temperatures. However it can be seen that the variance in some of them is not very high, meaning that the high mean difference in temperature is stable. Probably one area has been constantly influenced by other areas. On the other hand we can find areas with a relatively low mean with a high variance, meaning that the physical space presents temperature disruptions very often. The reason to choose this value is that the clustering algorithms fail in linking AHUs with sensors above an approximate difference in temperature of 8 degrees Celsius in experiment 2.

These values are to be used as a measure on how much the defined AHU linkage can be trusted, therefore we decide to test the following hypothesis: we define a limit for the mean difference in temperature as 7.5, above which the cluster may not be reliable as the internal heat disruption may be too high so we are not able to ensure that the cluster is comprising these points is correct.

For the first experiment, three AHUs have been selected. We chose, according to table 2, the ones with a relatively low mean and variance thus, the ones with a lower internal difference in temperature. Also it is necessary to validate how good the solution is, so the site engineers are consulted on the room names to which the AHUs are supplying air to. The chosen AHUs and corresponding rooms are:

• AHU05 supplying air to the 'NPI room'
• AHU11 supplying air to the 'Visual inspection', 'Manual inspection' and 'X Ray' rooms
• AHU16 supplying air to the 'Shell drying' room

We apply the methodologies discussed to this first experiment and later we compare our results with the real connections. In a properly clustered group, the AHU's SAT should be in the same group together with their corresponding space temperature sensors. The dendrograms correspond to the AHC methodology with the three distance metrics. AHC can be represented as a dendrogram because the algorithm progressively separates the clusters based on the different distance metrics until all the points form a separate cluster. For this reason, the lasso clustering methodology is represented separately in table 3, as this algorithm does not respond to this type of representation because the search for the optimal penalisation parameter is done in an iteratively refined grid.

For the experiments, we assume that we know already that the BMS points belong to the different AHU systems and room temperature sensors. What we assume not to know beforehand are which of the 13 AHUs is associated with which of the 11 physical spaces. We know the parent-child relationship once the AHU has been associated with its corresponding physical space.

In the dendrograms, a distance between clusters has been chosen to the best convenience. The distance chosen determines which branches below that distance form a cluster. The same distances are used later in the experiment with the whole building evaluation. An ideal cluster should contain the AHU sensors together with their respective temperature sensors within the same. Figures 4–6 show the clustering results of the different distance metrics in the form of a dendrogram.

Zoom In Zoom Out Reset image size

The fact that the AHUs' SAT belongs in the same cluster as their respective room temperature sensors means in the case of AHC that, in the lowest levels of the dendrogram, the time series are closer (in terms of the chosen distance) to each other, and form groups that are more distant to each other as the branches go up. In the case of graphical lasso, it means that the elements corresponding to the estimated inverse of the covariance matrix are zero between groups of time series data or clusters, thus forming groups with the time series data that present the most similar number of common features.

Figure 4 shows that correlations-based AHC clusters AHU05 and AHU11 properly with their corresponding temperature sensors. AHU16 is in a cluster together with one of its sensors, however the other sensor is excluded, thus forming a separate cluster. Figure 5 shows the same clustering methodology but based on DTW. On this, it can be seen that the clusters do not respond to the logic of the physical connections, independently of the distance chosen to form the clusters. Similarly in figure 6 the clusters do not respond to the ideal behavior either. In table 3 we show the results of lasso clustering. In this case, the clusters respond to the most ideal behavior. The three AHUs are properly linked to their respective temperature sensors within the same cluster.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Results have been summarised in table 4. The percentages are obtained based on the rate of room temperature sensors contained within the cluster with their corresponding AHU each. For instance, correlations AHC clusters elements properly within AHUs 05 and 11 but missed one of the two elements in cluster containing AHU16. The one showing the best results is the graphical lasso clustering, properly grouping AHUs with their respective areas of influence. This performance is followed by correlations clustering.

Now that the methodologies that perform best have been proved, we are proceeding with them both for the next experiment.

In experiment 2, we use the two best performing methodologies discussed in section 4.2, Pearson's distance-based AHC and lasso clustering, with all AHUs and all temperature sensors within the building. We can see the dendrogram as a result of Pearson's correlation distance-based AHC figure A1, and the results for lasso clustering in table A1 in the appendix. In the results of the Pearson's correlation distance-based AHC shown in figure A1. Correct clusters are considered by using the same distance to analyse the clusters as in experiment 1. As an example for correct cluster, figure A1 shows that AHU01 is contained within the same branch levels as the four lab rooms, as expected. Similarly, AHU07A and AHU07B are clustered with the Foundry area. Other areas such as Finish rooms, Wax rooms and the canteen, they are in different clusters with respect to their AHUs. Interestingly, we observe that the related rooms are usually found together in the lowest level of the branches. Another observation is that AHUs that share common spaces are also clustered together, such as AHU07A and AHU07B. In the results of lasso clustering in table A1, as in experiment 1, we used two columns to separate AHUs and room temperature sensors within the same cluster. Some of the clusters only group room temperature sensors, but no AHU SAT is present, as is the case of clusters 2, 5 and 9. On the other hand, clusters with only AHUs and no temperature sensors are observed in cluster 11. AHU06 is discarded for the lack of sensors in the room.

As defined in section 4.1, we set a limit value of 7.5 for the difference between SAT and RAT. Below in table 5 we summarised the results obtained. The table describes the result, the average difference in temperature, confirmation/denial of linkage successful in the appropriate category and the issue associated with the AHU if applicable.

• Issue (a): Open plan space. Heat exchange occurring between nearby rooms.
• Issue (b): The AHU is enabled and demand is 100% cooling near constantly but is having no effect on setpoint (the temperature value set) or is not enough to cool to the set-point. This implies the AHU is not mechanically sound or capable to meet requirements but onsite investigation would confirm this.

Table 5 shows the results of this experiment. The performance of both methodologies is very similar, except that lasso fails to cluster wax rooms with their respective temperature sensors. In general, we can say that correlations based AHC performs better than graphical lasso in terms of number of sensors belonging to the correct cluster. 70% reveals that this methodology discovers underlying physical patterns between AHUs and their respective spaces. Also, it can be observed that the AHUs with a mean ΔT above the set-up limit fail to predict the matches between AHUs and physical spaces in general, although exceptions can be seen in graphical lasso and AHU03. AHUs 07 A and B are in the same branch as their respective sensors in the AHC technique. AHUs 2 and 12 have a very high variance (as shown in table 2), which does not seem to affect the performance of the algorithms in the case of AHU 2. For AHU 12, both methods fail to predict its only sensor with the canteen. When looking at its variance, it seems to be higher than other AHUs, as the canteen is crowded mainly during lunch time. This, together with the fact that its only sensor may be misplaced, could explain this issue.

We wanted to test the hypothesis of the difference in temperature having a determinant effect on obtaining these relationships. In the case of AHU07 A & B and despite being an open space, correlations-based AHC has been able to properly identify all 8 sensors in this more challenging case. Therefore we conclude this hypothesis cannot be confirmed or the experiment is insufficient.

The main goal of this experiment is to link the respective rooms, each being fed by a different FCU, with their closest AHU unit. The strength of these links has been defined by the Person's correlation coefficient, with high showing a strong, very likely correlation, medium showing a moderate strength and low a weak link. These relationships are shown per FCU, and the links of these with all (if any) of the correlated AHUs. This is summarised in table 6.

We conclude from correlations that the most influential AHUs to the FCUs are AHU 01, 02 and 13. Without having information about the equipment links, we have some degree of confidence to discard low-flagged category AHUs. After validation with the engineers, we found out that AHU02 is the one that has the highest influence over these areas. Therefore, we uncovered some other potential dependencies from the supply air temperatures of AHUs 01 and 13 to some of the rooms (Nickel, Tugsten and Gym for AHU01 and Cobalt, Rhenium and Tantalum for AHU13).