Zhluková analýza (hierarchické klastrovanie) – spotreba elektriny

Úvod

Cieľom je ukázať hierarchickú zhlukovú analýzu na dátach o spotrebe elektriny. Databáza obsahuje merania v 10-minútovom kroku (teplota, vlhkosť, vietor, difúzne žiarenie a spotreba v troch zónach).
Keďže ide o časový rad s veľkým počtom riadkov, zhlukovanie na úrovni jednotlivých 10-minútových meraní by bolo neprehľadné. Preto budeme zhlukovať dni: každý deň opíšeme priemernými hodnotami premenných a následne zhlukujeme podobné dni.

---

library(dplyr)
library(knitr)
library(kableExtra)

Načítanie a príprava dát

rm(list = ls())

# 1) Načítanie dát
data_raw <- read.csv("powerconsumption.csv", stringsAsFactors = FALSE)

# 2) Prevod dátumu a času
# Formát v súbore vyzerá napr. "1/1/2017 0:00"
data_raw$Datetime <- as.POSIXct(data_raw$Datetime, format = "%m/%d/%Y %H:%M")

# 3) Kontrola základných informácií
str(data_raw)

'data.frame':   52416 obs. of  9 variables:
 $ Datetime              : POSIXct, format: "2017-01-01 00:00:00" "2017-01-01 00:10:00" ...
 $ Temperature           : num  6.56 6.41 6.31 6.12 5.92 ...
 $ Humidity              : num  73.8 74.5 74.5 75 75.7 76.9 77.7 78.2 78.1 77.3 ...
 $ WindSpeed             : num  0.083 0.083 0.08 0.083 0.081 0.081 0.08 0.085 0.081 0.082 ...
 $ GeneralDiffuseFlows   : num  0.051 0.07 0.062 0.091 0.048 0.059 0.048 0.055 0.066 0.062 ...
 $ DiffuseFlows          : num  0.119 0.085 0.1 0.096 0.085 0.108 0.096 0.093 0.141 0.111 ...
 $ PowerConsumption_Zone1: num  34056 29815 29128 28229 27336 ...
 $ PowerConsumption_Zone2: num  16129 19375 19007 18361 17872 ...
 $ PowerConsumption_Zone3: num  20241 20131 19668 18899 18442 ...

Pre ďalšie kroky vytvoríme denný dataset (priemery za deň). Spotrebu spojíme aj do jednej premennej TotalPower (súčet zón).

data_raw <- data_raw %>%
  mutate(
    Date = as.Date(Datetime),
    TotalPower = PowerConsumption_Zone1 + PowerConsumption_Zone2 + PowerConsumption_Zone3
  )

# Denné priemery
data_day <- data_raw %>%
  group_by(Date) %>%
  summarise(
    across(
      .cols = c(Temperature, Humidity, WindSpeed, GeneralDiffuseFlows, DiffuseFlows,
                PowerConsumption_Zone1, PowerConsumption_Zone2, PowerConsumption_Zone3, TotalPower),
      .fns = ~mean(.x, na.rm = TRUE)
    ),
    .groups = "drop"
  )

# Nastavenie riadkových mien pre prehľadnejšie výstupy
data_day_df <- as.data.frame(data_day)
data_day_df <- data_day_df %>% 
  filter(!is.na(Date))

# Rýchla ukážka (len prvých 10 dní)
kable(head(data_day, 10), caption = "Tab. 1: Ukážka denných priemerov (prvých 10 dní)") %>%
  kable_styling(full_width = FALSE)

Tab. 1: Ukážka denných priemerov (prvých 10 dní)

Date	Temperature	Humidity	WindSpeed	GeneralDiffuseFlows	DiffuseFlows	PowerConsumption_Zone1	PowerConsumption_Zone2	PowerConsumption_Zone3	TotalPower
2017-01-01	9.675299	68.51931	0.3151458	121.39077	25.99392	28465.23	17737.79	17868.80	64071.82
2017-01-02	12.476875	71.45632	0.0765625	120.40449	27.22741	28869.49	19557.73	17820.76	66247.98
2017-01-03	12.100000	74.98167	0.0767153	120.68601	28.57466	30562.45	20057.27	17620.80	68240.52
2017-01-04	10.509479	75.45979	0.0824167	122.95932	28.82722	30689.83	20102.08	17673.69	68465.60
2017-01-05	10.866444	71.04049	0.0838958	118.74986	29.74144	30802.91	20033.94	17664.18	68501.03
2017-01-06	11.685208	77.25778	0.1459514	63.70045	59.33356	30712.24	20139.77	17481.08	68333.09
2017-01-07	12.018819	76.21708	0.0829306	88.99876	48.84979	30612.70	19252.43	17985.94	67851.08
2017-01-08	14.724653	69.10924	0.0813056	71.08221	46.38776	29559.92	17221.12	17810.60	64591.64
2017-01-09	15.490139	68.33236	0.0741736	110.43367	55.82734	30540.25	19102.63	17417.79	67060.68
2017-01-10	14.019583	75.67840	0.0782917	112.49842	47.99520	31863.29	19405.27	17425.34	68693.90

Škálovanie premenných

Pri zhlukovaní používame vzdialenosti. Aby premenné s veľkými jednotkami (napr. spotreba) nedomninovali všetko ostatné, použijeme z-škálovanie:

\[ z = \frac{x-\mu}{\sigma} \]

data_complete <- data_day_df %>%
  dplyr::select(-Date) %>%   # výslovne dplyr::select
  na.omit()

data_scaled <- scale(data_complete)

Boxploty škálovaných premenných

num_vars <- as.data.frame(data_scaled)
num_plots <- ncol(num_vars)

par(mfrow = c(ceiling(sqrt(num_plots)), ceiling(num_plots / ceiling(sqrt(num_plots)))))
par(mar = c(4, 4, 2, 1))

for (col in names(num_vars)) {
  boxplot(
    num_vars[[col]],
    main = col,
    col = "lightblue",
    horizontal = TRUE
  )
}

mtext("Boxploty škálovaných denných premenných", outer = TRUE, cex = 1.2, font = 2)

Korelácie premenných

Ak by boli niektoré premenné extrémne korelované, jedna z dvojice by zbytočne „zdvojovala“ informáciu. Pozrieme sa na korelačnú maticu.

cor_mat <- cor(data_scaled, use = "pairwise.complete.obs")
cor_mat <- round(cor_mat, 2)

kable(cor_mat, caption = "Tab. 2: Korelačná matica (škálované denné dáta)") %>%
  kable_styling(full_width = FALSE)

Tab. 2: Korelačná matica (škálované denné dáta)

	Temperature	Humidity	WindSpeed	GeneralDiffuseFlows	DiffuseFlows	PowerConsumption_Zone1	PowerConsumption_Zone2	PowerConsumption_Zone3	TotalPower
Temperature	1.00	-0.31	0.56	0.64	0.08	0.73	0.42	0.61	0.73
Humidity	-0.31	1.00	-0.17	-0.39	-0.04	-0.21	-0.21	-0.26	-0.29
WindSpeed	0.56	-0.17	1.00	0.37	-0.08	0.47	0.30	0.41	0.49
GeneralDiffuseFlows	0.64	-0.39	0.37	1.00	0.35	0.54	0.11	0.58	0.55
DiffuseFlows	0.08	-0.04	-0.08	0.35	1.00	0.11	-0.25	0.22	0.08
PowerConsumption_Zone1	0.73	-0.21	0.47	0.54	0.11	1.00	0.39	0.71	0.87
PowerConsumption_Zone2	0.42	-0.21	0.30	0.11	-0.25	0.39	1.00	0.23	0.59
PowerConsumption_Zone3	0.61	-0.26	0.41	0.58	0.22	0.71	0.23	1.00	0.90
TotalPower	0.73	-0.29	0.49	0.55	0.08	0.87	0.59	0.90	1.00

Matica vzdialeností (ukážka)

Vzdialenosť dvoch dní počítame euklidovsky zo všetkých použitých premenných. Celá matica je veľká (máme veľa dní), preto ukážeme len prvých 12 × 12.

dist_obj <- dist(data_scaled, method = "euclidean")

# Ukážka prvých 12 dní
dist_mat_small <- as.matrix(dist_obj)[1:12, 1:12]
dist_mat_small <- round(dist_mat_small, 2)

kable(dist_mat_small, caption = "Tab. 3: Ukážka matice vzdialeností (prvých 12 dní)") %>%
  kable_styling(full_width = FALSE)

Tab. 3: Ukážka matice vzdialeností (prvých 12 dní)

1	2	3	4	5	6	7	8	9	10	11	12
0.00	0.95	1.47	1.46	1.38	1.92	1.52	1.34	1.69	1.94	1.54	2.30
0.95	0.00	0.77	0.92	0.86	1.41	1.05	1.29	1.19	1.37	1.25	1.76
1.47	0.77	0.00	0.32	0.43	1.06	0.72	1.58	1.19	0.84	1.17	1.29
1.46	0.92	0.32	0.00	0.40	1.08	0.79	1.75	1.41	0.99	1.30	1.39
1.38	0.86	0.43	0.40	0.00	1.15	0.84	1.61	1.21	0.99	1.05	1.37
1.92	1.41	1.06	1.08	1.15	0.00	0.54	1.59	1.27	0.94	1.21	0.87
1.52	1.05	0.72	0.79	0.84	0.54	0.00	1.25	1.02	0.69	0.94	0.92
1.34	1.29	1.58	1.75	1.61	1.59	1.25	0.00	1.00	1.49	1.03	1.74
1.69	1.19	1.19	1.41	1.21	1.27	1.02	1.00	0.00	0.91	0.51	1.07
1.94	1.37	0.84	0.99	0.99	0.94	0.69	1.49	0.91	0.00	0.92	0.62
1.54	1.25	1.17	1.30	1.05	1.21	0.94	1.03	0.51	0.92	0.00	1.07
2.30	1.76	1.29	1.39	1.37	0.87	0.92	1.74	1.07	0.62	1.07	0.00

Hierarchické zhlukovanie (Wardova metóda)

Wardova metóda je aglomeratívna: začína s jednočlennými zhlukmi a spája ich tak, aby bol nárast vnútornej variability čo najmenší. Prakticky často vytvára kompaktné a interpretovateľné zhluky.

hc <- hclust(dist_obj, method = "ward.D2")

plot(
  hc,
  labels = FALSE,
  main = "Hierarchické zhlukovanie dní (Ward.D2)",
  xlab = "Dni (bez popisov pre čitateľnosť)",
  sub = ""
)

k <- 3
rect.hclust(hc, k = k, border = "red")

Príslušnosť dní do klastrov

cluster_membership <- cutree(hc, k = k)

clusters_df <- data.frame(
  Date = rownames(data_complete),
  Cluster = factor(cluster_membership)
)

kable(head(clusters_df, 20), caption = "Tab. 4: Prvých 20 dní a ich klaster") %>%
  kable_styling(full_width = FALSE)

Tab. 4: Prvých 20 dní a ich klaster

Date	Cluster
1	1
2	1
3	1
4	1
5	1
6	1
7	1
8	1
9	1
10	1
11	1
12	1
13	1
14	1
15	1
16	1
17	1
18	1
19	1
20	1

Centroidy (priemery premenných v klastroch)

Tu už pracujeme s pôvodnými (neškálovanými) dennými priemermi, aby boli hodnoty interpretovateľné.

data_day_for_desc <- as.data.frame(data_complete)
data_day_for_desc$Cluster <- factor(cluster_membership)

centroids <- data_day_for_desc %>%
  group_by(Cluster) %>%
  summarise(across(where(is.numeric), ~mean(.x, na.rm = TRUE)))

kable(centroids, digits = 2, caption = "Tab. 5: Centroidy – priemerné hodnoty premenných v klastroch") %>%
  kable_styling(full_width = FALSE)

Tab. 5: Centroidy – priemerné hodnoty premenných v klastroch

Cluster	Temperature	Humidity	WindSpeed	GeneralDiffuseFlows	DiffuseFlows	PowerConsumption_Zone1	PowerConsumption_Zone2	PowerConsumption_Zone3	TotalPower
1	14.19	71.19	0.76	113.64	65.35	30272.39	20701.00	15304.36	66277.75
2	20.79	67.95	2.17	225.55	87.75	32957.74	19918.76	16876.61	69753.11
3	26.31	61.40	4.59	265.71	71.84	36322.49	24400.04	26459.26	87181.79

Záver

Zhluková analýza rozdelila dni do niekoľkých typov podľa kombinácie počasia (teplota, vlhkosť, vietor, difúzne žiarenie) a spotreby elektriny (zónovo aj spolu). Takéto zhluky sa dajú interpretovať napríklad ako „chladné dni s vyššou spotrebou“, „teplé dni s nižšou spotrebou“, prípadne dni so špecifickým priebehom žiarenia a odlišnou spotrebou. V praxi môže byť výsledok využiteľný na segmentáciu dní pre plánovanie výroby, nákupu energie alebo detekciu atypických období.