← Overview
Partitional clustering: k-means and k-medoids¶
Partitional clustering assigns every period to exactly one cluster and describes each cluster by a single representative period. k-means and k-medoids share the same goal — both minimize the total within-cluster distance $J$ introduced in section 2 — and differ only in how they reach it. k-medoids restricts every center to a real period, which turns the search into a finite problem an MILP (Mixed-Integer Linear Program) solver can settle exactly. k-means allows a synthetic center and reaches it approximately, via Lloyd's heuristic.
Each algorithm also carries a standard representation — the kind of center it produces, and the reason for its name: k-means averages its members, k-medoids picks a real member. That pairing is a convention of the classical algorithms, not a constraint: tsam lets you combine any clustering with any representation (see Representation).
| Method | Standard representation | Representatives chosen to… | How it is solved |
|---|---|---|---|
| k-means | mean |
minimize the total within-cluster distance $J$ (approximately) | Lloyd's heuristic → local optimum |
| k-medoids | medoid |
minimize that same $J$, centers restricted to real periods (exactly) | MILP → global optimum (needs a solver) |
This notebook is about how each method formulates and solves that problem: what goes in, what happens, and what comes out. For applying them to your own data, see the how-to guides.
1 What comes in: the preprocessed period matrix¶
Clustering starts from the preprocessed data: the normalized,
unstacked period matrix $D$. Each of the six rows is one day, flattened into eight
coordinates — two attributes (solar, load) × four timesteps.
import numpy as np
import pandas as pd
import plotly.express as px
import plotly.io as pio
pio.renderers.default = "notebook_connected"
# Preprocessed period matrix D (normalized + unstacked) from 01_preprocessing.
tiny_period_df = pd.read_csv(
"../../../data/tiny_periods.csv", header=[0, 1], index_col=0
)
tiny_period_array = tiny_period_df.values # shape (6, 8): six periods, eight features
N_PERIODS = tiny_period_array.shape[0]
N_ATTRS, N_TIMESTEPS = 2, 4
tiny_period_df.round(4)
| solar | load | |||||||
|---|---|---|---|---|---|---|---|---|
| TimeStep | 0 | 1 | 2 | 3 | 0 | 1 | 2 | 3 |
| PeriodNum | ||||||||
| 0 | 0.0 | 1.000 | 0.750 | 0.0 | 0.0000 | 0.0000 | 0.1429 | 0.2857 |
| 1 | 0.0 | 0.875 | 0.875 | 0.0 | 0.0000 | 0.0000 | 0.1429 | 0.1429 |
| 2 | 0.0 | 0.375 | 0.250 | 0.0 | 0.1429 | 0.2857 | 0.2857 | 0.4286 |
| 3 | 0.0 | 0.250 | 0.375 | 0.0 | 0.1429 | 0.1429 | 0.4286 | 0.2857 |
| 4 | 0.0 | 0.125 | 0.125 | 0.0 | 0.4286 | 0.4286 | 0.5714 | 0.5714 |
| 5 | 0.0 | 0.125 | 0.000 | 0.0 | 0.4286 | 0.5714 | 0.7143 | 1.0000 |
2 The objective: the distance that gets minimized in k-means and k-medoids¶
The shared building block of k-means and k-medoids is the distance between a period and a center. Writing $x_{p,a,t}$ for the value of attribute $a$ at timestep $t$ in period $p$:
$$ \text{dist}(x_p, c_k) = \sqrt{\sum_{a=1}^{N_a} \sum_{t=1}^{N_t} (x_{p,a,t} - c_{k,a,t})^2} $$
The double sum walks every coordinate of the period vector:
- $a = 1 \dots N_a$ indexes the attributes — here $N_a = 2$ (
solar,load); - $t = 1 \dots N_t$ indexes the timesteps within a period — here $N_t = 4$ (the four 6-hourly steps of a day).
so it runs over all $N_a \cdot N_t = 2 \times 4 = 8$ coordinates — exactly the eight columns of $D$.
The subtraction inside the sum is element-wise: coordinate $(a,t)$ of the period is only
ever compared against the same coordinate $(a,t)$ of the center — solar@t1 of period $p$
against solar@t1 of center $k$, load@t3 against load@t3, and so on. Nothing is ever
compared across attributes or across timesteps. So $x_p$ and $c_k$ are points of the same
shape, and a center is directly comparable to a period whether it is a synthetic mean or a
real day.
This is also why the column order matters: the pairing is positional, so $x_p$ and $c_k$
must carry their coordinates in the same order for solar@t1 to meet solar@t1. tsam keeps
that order fixed for you (it changed in v4 — see
column order).
From this single distance, the quality of a whole clustering is the total within-cluster distance $J$ — every period summed against the center of the cluster it lands in:
$$ J = \sum_{k=1}^{N_k} \sum_{p \in \mathbb{C}_k} \text{dist}(x_p, c_k)^2 $$
| Symbol | Meaning |
|---|---|
| $x_p$ | period $p$ — one row of $D$ (D_arr[p]), a point in 8-D space |
| $\mathbb{C}_k$ | cluster $k$: the set of periods assigned to group $k$ |
| $c_k$ | the center of cluster $k$ |
| $N_k$ | number of clusters ( = n_clusters) |
$J$ is the yardstick k-means and k-medoids both drive down. The following cell shows how the distance is calculated, using day 0 as the period and day 2 as a stand-in center.
def euclidean_dist(x_p, c_k):
"""Euclidean distance between a period and a cluster center.
x_p, c_k : np.ndarray, shape (N_a * N_t,) == (8,)
A period and a center — same shape and coordinate order, one entry per
(attribute a, timestep t) pair.
Returns the scalar sqrt(sum((x_p - c_k) ** 2)); zero only for identical vectors.
"""
return float(np.sqrt(np.sum((x_p - c_k) ** 2)))
# A concrete center pair: period 0 vs a stand-in center (period 2).
x_p = tiny_period_array[0] # day 0 (a sunny day)
c_k = tiny_period_array[2] # day 2 (an overcast day), used here as a stand-in center
# Label every coordinate with the (a, t) it belongs to — the 8 columns of D.
terms = pd.DataFrame(
{
"a (attribute)": tiny_period_df.columns.get_level_values(0),
"t (timestep)": tiny_period_df.columns.get_level_values(1).astype(int),
"x_p": x_p,
"c_k": c_k,
"gap = x_p - c_k": x_p - c_k,
"gap**2": (x_p - c_k) ** 2,
}
)
print(terms.round(3).to_string(index=False))
# One term of the double sum picked out, e.g. a = solar, t = 1 (coordinate index 1):
print(
"\nExample single term a=solar, t=1:",
f"(x = {x_p[1]:.3f} - c = {c_k[1]:.3f})**2 = {(x_p[1] - c_k[1]) ** 2:.3f}",
)
print(
"sum over all",
N_ATTRS * N_TIMESTEPS,
"terms =",
round(float(((x_p - c_k) ** 2).sum()), 4),
)
print("dist(x_p, c_k) = sqrt(sum) =", round(euclidean_dist(x_p, c_k), 4))
a (attribute) t (timestep) x_p c_k gap = x_p - c_k gap**2
solar 0 0.000 0.000 0.000 0.000
solar 1 1.000 0.375 0.625 0.391
solar 2 0.750 0.250 0.500 0.250
solar 3 0.000 0.000 0.000 0.000
load 0 0.000 0.143 -0.143 0.020
load 1 0.000 0.286 -0.286 0.082
load 2 0.143 0.286 -0.143 0.020
load 3 0.286 0.429 -0.143 0.020
Example single term a=solar, t=1: (x = 1.000 - c = 0.375)**2 = 0.391
sum over all 8 terms = 0.7835
dist(x_p, c_k) = sqrt(sum) = 0.8851
3 Approximate solution — k-means (Lloyd's algorithm)¶
Mechanism: Lloyd's algorithm does not solve the objective exactly; it converges to a local optimum by alternating two cheap steps until assignments stop changing:
- Initialize $k$ centers $c_1, \dots, c_k$ (e.g. k-means++).
- Assignment step: put each period with its nearest center, $\text{cluster}(p) = \arg\min_k \text{dist}(x_p, c_k)$.
- Update step: move each center to the mean of its members, $c_k = \frac{1}{|\mathbb{C}_k|} \sum_{p \in \mathbb{C}_k} x_p$.
- Repeat from step 2.
Its default representative is therefore the cluster mean — a synthetic centroid that need not match any real day.
From-scratch Lloyd iteration on the tiny series¶
With the distance function in hand, the assignment step is just "call euclidean_dist for
every period against every center and take the nearest", and the update step is the
centroid mean. Tracing it on the six periods to find three typical days ($k=3$), with days 0, 2, 4 as the initial cluster centers:
# Reuse euclidean_dist(x_p, c_k) defined above.
n_clusters = 3
# Deterministic initialization: pick periods 0, 2, 4 as initial centers
centers = tiny_period_array[[0, 2, 4]].copy().astype(float)
print("Initial centers (periods 0, 2, 4):")
for k, center in enumerate(centers):
print(f" c{k} = {center.round(3)}")
for iteration in range(6):
# --- Step 1: assign each period to its nearest center ---
distances = np.empty((N_PERIODS, n_clusters))
for period_idx in range(N_PERIODS):
for k in range(n_clusters):
distances[period_idx, k] = euclidean_dist(
tiny_period_array[period_idx], centers[k]
)
assignments = np.empty(N_PERIODS, dtype=int)
for period_idx in range(N_PERIODS):
assignments[period_idx] = np.argmin(distances[period_idx])
# --- Step 2: recompute each center as the mean of its assigned periods ---
new_centers = np.empty_like(centers)
for k in range(n_clusters):
cluster_members = []
for p in range(N_PERIODS):
if assignments[p] == k:
cluster_members.append(tiny_period_array[p])
if len(cluster_members) > 0:
new_centers[k] = np.mean(cluster_members, axis=0)
else:
new_centers[k] = centers[k] # keep old center if no period was assigned
converged = np.allclose(centers, new_centers)
centers = new_centers
print(f"\nIteration {iteration + 1}: assignments = {assignments}")
if converged:
print(" Converged.")
break
print("\nFinal cluster assignments:")
for period_idx in range(N_PERIODS):
print(f" day_{period_idx} -> cluster {assignments[period_idx]}")
Initial centers (periods 0, 2, 4): c0 = [0. 1. 0.75 0. 0. 0. 0.143 0.286] c1 = [0. 0.375 0.25 0. 0.143 0.286 0.286 0.429] c2 = [0. 0.125 0.125 0. 0.429 0.429 0.571 0.571] Iteration 1: assignments = [0 0 1 1 2 2] Iteration 2: assignments = [0 0 1 1 2 2] Converged. Final cluster assignments: day_0 -> cluster 0 day_1 -> cluster 0 day_2 -> cluster 1 day_3 -> cluster 1 day_4 -> cluster 2 day_5 -> cluster 2
What comes out: the cluster assignment. Lloyd's loop returns one number per period — the
array printed above. Entry p is the cluster that day p landed in, in period order, so
[0 0 1 1 2 2] reads "days 0–1 together, days 2–3 together, days 4–5 together" — exactly the
three shape-pairs of the dataset (sunny / overcast / cloudy).
That array is the output of clustering. Everything the rest of the pipeline does — choosing
a representation per group, rescaling
the results, counting how often each typical period occurs — reads only this grouping, plus
the periods themselves. The centers Lloyd computed along the way are scaffolding: k-means
keeps them as its mean representative, but any other representation recomputes its own from
the members.
One property matters downstream: the cluster labels are arbitrary names. Which group is called 0 and which is called 2 depends on where the centers were initialized, so a different method — or a different starting point — can produce the same grouping under different numbers. Only the partition is meaningful: which days share a cluster. Never write code that depends on a particular label.
4 Exact solution — k-medoids (MILP)¶
Mechanism: where Lloyd's algorithm only approximates the minimum of $J$, k-medoids solves it exactly. Restricting every center to be an actual period turns the search into a finite combinatorial problem — the classic $p$-median / facility-location problem (known in spatial planning as the Hess model) — which can be written as a Mixed-Integer Linear Program (MILP) and handed to a solver (here HiGHS) for a provably globally optimal clustering.
The model¶
The only data the solver needs is the matrix of pairwise distances $d_{i,j}$ between periods. A single family of binary variables encodes the whole clustering:
$$ z_{i,j} = \begin{cases} 1 & \text{period } j \text{ is assigned to center } i \\ 0 & \text{otherwise} \end{cases} $$
The diagonal $z_{i,i}=1$ marks period $i$ as one of the open centers (a medoid). The program minimizes the total within-cluster distance,
$$ \min_{z}\ \sum_{i}\sum_{j} d_{i,j}\, z_{i,j} \;=\; J, $$
subject to three constraints — each one line of _setup_k_medoids in k_medoids_exact.py:
| # | constraint | meaning | rule in _setup_k_medoids |
|---|---|---|---|
| 1 | $\sum_i z_{i,j} = 1\ \ \forall j$ | every period is assigned to exactly one center | candToClusterRule |
| 2 | $\sum_i z_{i,i} = k$ | exactly $k$ centers are opened | noClustersRule |
| 3 | $z_{i,j} \le z_{i,i}\ \ \forall i,j$ | a period may only be assigned to an open center | clusterRelationRule |
What $i$, $j$, $z$, $d$ mean on our six days¶
Everything above is concrete on the tiny dataset:
- $i$ and $j$ both run over the six days $0, 1, \dots, 5$ — every day is at once a candidate center (index $i$) and a period to be assigned (index $j$).
- $d_{i,j}$ is simply the distance between day $i$ and day $j$ — one cell of the $6\times6$
matrix below. For example $d_{0,2}\approx 0.89$ is exactly the
dist(day 0, day 2) = 0.8851worked out in section 2. - $z_{i,j}=1$ reads "day $j$ is represented by day $i$"; the diagonal $z_{i,i}=1$ flags day $i$ as one of the three chosen medoids.
So the optimizer never looks at solar/load profiles — it works from the distance matrix alone, choosing which 3 days act as centers and attaching every day to the cheapest one. A caveat lives in the formulation itself: there are $n^2$ binary variables for $n$ periods, so the exact route is practical for hundreds of periods, not tens of thousands.
The solver's only input: the distance matrix¶
The exact-MILP k-medoids never sees the day profiles at all — its solver receives only
the matrix of pairwise distances $d_{i,j}$ (M.d inside _setup_k_medoids). Below we build
that full $6\times6$ matrix for the tiny dataset and hand it straight to tsam's own model
builder and solver, then read the optimal medoids and assignment back out.
# Build the 6x6 pairwise distance matrix — exactly the input `M.d` the MILP receives.
Dmat = np.array(
[
[
euclidean_dist(tiny_period_array[i], tiny_period_array[j])
for j in range(N_PERIODS)
]
for i in range(N_PERIODS)
]
)
# Hand it to tsam's *own* k-medoids model builder + solver. No day profiles involved.
from tsam.algorithms.k_medoids_exact import _setup_k_medoids, _solve_given_pyomo_model
model = _setup_k_medoids(Dmat, n_clusters=3)
r_x, r_y, r_obj = _solve_given_pyomo_model(model, solver="highs")
medoids = [
p for p, opened in enumerate(r_y) if opened == 1
] # the periods with z_ii = 1
assignment = r_x.argmax(axis=1) # period j -> its medoid i
print("opened medoids (z_ii = 1): ", medoids)
print("each period -> its medoid: ", assignment.tolist())
print("objective sum d_ij * z_ij = ", round(r_obj, 4))
opened medoids (z_ii = 1): [0, 3, 5] each period -> its medoid: [0, 0, 3, 3, 5, 5] objective sum d_ij * z_ij = 1.0214
labels = [f"day_{p}" for p in range(N_PERIODS)]
fig = px.imshow(
Dmat,
x=labels,
y=labels,
text_auto=".2f",
color_continuous_scale="Blues",
labels={"x": "candidate center", "y": "period", "color": "distance"},
title="Pairwise distance matrix — the exact input to the k-medoids optimizer",
)
# Outline the cell each period contributes: period j -> the medoid the MILP assigned it to.
for j in range(N_PERIODS):
m = int(assignment[j])
fig.add_shape(
type="rect",
x0=m - 0.5,
x1=m + 0.5,
y0=j - 0.5,
y1=j + 0.5,
line={"color": "#EF553B", "width": 3},
)
fig.update_layout(width=560, height=520)
fig.show()
obj = sum(Dmat[j, int(assignment[j])] for j in range(N_PERIODS))
print(f"Objective = sum of outlined cells (day -> its medoid) = {obj:.3f}")
Objective = sum of outlined cells (day -> its medoid) = 1.021
Both axes are indexed by the same day_0 … day_5 from the legend, so a cell (i, j) is the
distance between two specific days. The block structure is the whole story: days 0–1 (sunny),
2–3 (overcast) and 4–5 (cloudy) sit cheaply close (dark), while crossing between blocks is
expensive (light). The optimizer opens one center per block ($z_{i,i}=1$, the medoids) and
assigns every period to it ($z_{i,j}=1$, the outlined cells) — the explicit form of the
abstract "$p \in \mathbb{C}_k$" membership from section 2. Here the opened medoids [0, 3, 5]
are exactly one day from each shape-block: a sunny day (0), an overcast day (3) and the
cloudy/extreme-load day (5). It minimizes $\sum_{i,j} d_{i,j}\, z_{i,j}$, the sum of the
outlined cells: here $\approx 1.02$ (the diagonal day → itself terms are zero).
Notice the optimum can be degenerate: inside a symmetric pair either day is an equally good medoid, so the solver's particular pick is an arbitrary tie-break — but the partition into pairs is unique and globally optimal.
Up next:
- Agglomerative clustering — hierarchical Ward and contiguous Ward, which reach a partition by merging bottom-up instead of refining iteratively
- Extremal-prototype selection — k-maxoids, which maximizes spread between representatives instead of minimizing $J$
- Comparing clustering methods — watch the methods reach different partitions on the same data
See also:
- Representation — choosing what the cluster center looks like, independently of the method that formed the cluster
- Extreme periods — preserving peaks through clustering