Chapter 2: Clustering Methodology

This part exists because a Distribution System Operator (DSO) can no longer monitor its Distributed Energy Resources (DER)-connected customers one by one once their numbers move past a few hundred, the same low-voltage bottleneck Chapter 1 opened on, and clustering is how that population becomes tractable again. Chapter 1 left four questions open, on purpose: which objective family fits a given population, what distance and representation a load curve actually needs, how a discovered grouping can be checked without ground truth, and how background knowledge should constrain a clustering directly. Those were the right questions to raise before clustering a single customer, and the wrong ones to answer in the abstract. This chapter settles all four, so the rest of this part has a fixed recipe to build on instead of working the same ground over again.

A search across 1,284 London households, ten feature sets, five scalers, five dimension-reduction families, and three clustering algorithms once returned a combination that scored a silhouette of 0.994, the best result of several hundred trials. Checking what that combination actually produced, rather than trusting the number, told a different story: 1,282 households in one cluster, and one household each in two more. The “clusters” were not archetypes. One household, MAC000037, records almost no evening consumption at all, giving it a morning_evening_ratio of 1,296 against a population median under 1. The search had not found a customer segment. It had found one broken ratio, and rewarded it handsomely.

Chapter 1 already warned against trusting a single validity index, citing Jain’s retrospective on fifty years of work past k-means [1]. The standard response is to check several indices together, silhouette alongside two others this chapter introduces properly shortly, reasoning that if independently constructed formulas agree, the agreement must mean something. MAC000037 shows why that reasoning fails: all three rated the degenerate split as excellent, because however differently each is built, all three only measure how separated groups are. An outlier, by definition, is perfectly separated. None of them can ask the question that matters: is this a customer archetype, or an artefact of one broken feature? Four questions follow from that gap, answered here in order:

Choosing a pipeline

Chapter 1’s “shape, not magnitude” convention and its discussion of Dynamic Time Warping (DTW), tensor factorisation, and diffusion maps established that distance and representation matter, without settling which combination wins for a given population [2], [3], [4]. Answering that concretely means trying enough combinations that the answer is a finding, not a hand-picked example. The ten feature sets on offer span full weekly and seasonal load-shape profiles, raw magnitude, day-to-day volatility, a linear trend, peak-normalised shape, autocorrelation, calendar-peak timing, and the marginal distribution of demand itself, deliberately mixing whole-profile representations with single, targeted summary statistics rather than ten variations on the same idea. A fixed grid, those ten feature sets times five scalers times a handful of dimension-reduction settings times three clusterers times three values of \(K\), is the obvious way to search this space, and also the wrong tool: a dimension-reduction technique’s continuous hyperparameters, Uniform Manifold Approximation and Projection (UMAP)’s neighbourhood size and minimum distance, t-SNE’s perplexity, cannot be explored by a grid without enumerating every value in advance [5], [6].

This chapter instead runs a proper multi-objective search: Non-dominated Sorting Genetic Algorithm (NSGA)-II, a genetic algorithm for problems with more than one objective at once, evaluated across 400 trials over feature set, scaler, dimension-reduction technique (Principal Component Analysis, Kernel Principal Component Analysis (PCA), t-SNE, UMAP, and PaCMAP, each with its own hyperparameters), clusterer, and \(K \in [2, 12]\) [7], [8], [9]. NSGA-II already draws its own initial population at random before crossover and mutation begin, giving the space coverage a later importance analysis needs without a separate warm-up phase. The search runs on Optuna, an open-source hyperparameter-optimisation framework built around this same trial-and-study abstraction [10]. Three objectives are optimised at once, all maximised: silhouette, a balance-entropy score, and a cheap stability proxy, both explained properly further into this chapter. None of the several hundred resulting trials dominates every other on all three axes at once; the Pareto front below is what survives that comparison.

The search itself is not free: 400 trials against 1,284 households, most of it spent fitting UMAP, t-SNE, and PaCMAP over and over, took about 53 minutes on a single machine. Everything this chapter builds on top of it afterwards, the audit, the cross-dataset check three utilities deep, added only another three, because Optuna had already saved its own study to disk rather than trusting anyone to remember to. That habit paid for itself the first time a downstream analysis needed rerunning and the search did not.

Every point on this front is a real trade-off, not noise: no other trial beats it on separation, balance, and stability simultaneously. Which of the three search dimensions actually drives that trade-off, quantified rather than eyeballed from a chart, comes from functional Analysis of Variance (fANOVA), functional analysis of variance over a random-forest surrogate fitted to the completed trials [11].

Which dimension-reduction technique is chosen accounts for more than half the variance in both silhouette and balance: whether a clustering looks good is mostly a question of geometry, how the reduction step warps local neighbourhood structure, not which features fed it. Feature-set choice is a clear second on those same two measures. Neither comes close to explaining the third bar, though: there, \(K\) alone dominates, more than every other parameter combined. Scaler choice barely moves either silhouette or balance, confirming, with a quantified decomposition behind it rather than a spread of medians, a finding an earlier exploratory pass over this same data already suggested informally. None of this ranks a feature set that was never in the running to begin with, though; ten menus were compared, not the space of every menu that could have been written.

fANOVA is a real, peer-reviewed method, not a tool to trust blindly because it produced a chart [11]. It fits a random-forest surrogate over the completed trials and reads importance off that surrogate, so the estimate is only as reliable as the surrogate’s fit, and this search tried PCA in only 14 of its 400 trials against 117 for t-SNE, meaning fANOVA’s read on PCA rests on far less evidence than its read on t-SNE. Importance is not direction either: it says how much a choice moves an outcome, never which specific value wins or why. A newer evaluator, PED-ANOVA, designed to handle exactly this kind of uneven trial count more robustly, gives a useful cross-check rather than a second opinion to pick and choose from [12]:

parameter fANOVA PED-ANOVA
0 k 0.621 0.913
1 clusterer 0.176 0.031
2 feature 0.104 0.028
3 reducer 0.068 0.020
4 scaler 0.031 0.007

Both agree on the ranking, \(K\) first by a wide margin, then clusterer, feature, reducer, scaler in that order, but PED-ANOVA is even more emphatic about \(K\)’s dominance for stability, 0.913 against fANOVA’s 0.629. Where two independently built estimators agree this strongly on the ranking, that agreement is worth trusting; where they diverge on magnitude, the honest reading is “the effect is real and \(K\)-dominated,” not a specific decimal to defend.

None of this says which feature set, dimension-reduction technique, scaler, or clustering algorithm actually wins, or why. A tempting shortcut is to group this search’s 400 completed trials by each choice and read the group means as that choice’s effect. That shortcut is actively misleading here.

Common Mistake

Grouping a search’s trial results by one parameter and reading the group means as that parameter’s effect. NSGA-II’s search population is not a balanced factorial design, so a category’s mean mixes its good pairings with its bad ones and answers a different question than the one asked. Checked directly against this search: the robust scaler shows the lowest mean cheap-stability of the five scalers tried, and Gaussian mixture the lowest mean of the three clusterers tried, despite robust and Gaussian mixture being the actual settled recipe. A marginal mean answers “how did this choice do on average across everything else that happened to vary alongside it,” not “what happens if only this choice changes.”

The honest way to answer “what happens if only this changes, holding everything else fixed” is a controlled ablation: fix every other choice at the settled recipe’s winning values, shape, robust, PaCMAP, Gaussian mixture, \(K = 2\), and vary exactly one axis at a time. Worth naming plainly first, though: this ablation is a local sensitivity check around the point the global search already found, not a substitute for it. NSGA-II found one good region of a search space it explored efficiently rather than exhaustively; the ablations below explain that region, they do not re-search for a better one. Where one surfaces a real surprise, that is evidence the original search’s space was incomplete, worth reporting as an open question, not license to quietly redefine the settled recipe from a single local test.

technique silhouette balance min_cluster_size
0 PCA 0.997 0.009 1
1 KernelPCA 0.502 0.767 287
2 TSNE 0.431 0.964 500
3 UMAP 0.598 0.788 303
4 PaCMAP 0.784 0.787 302
5 none 0.996 0.009 1

Every dimension-reduction technique was tried here, plus a sixth condition never in the original search space at all: no reduction, clustering directly on the scaled 48-dimensional shape vector. The result answers the manifold-versus-time-series question directly. PCA and no reduction land on the same failure, a silhouette near 0.997 and a smallest cluster of exactly one household, the identical degenerate split MAC000037 produced at the start of this chapter. A single linear projection, or no projection at all, cannot separate a genuine two-way behavioural split from one outlier sitting far from everything else; both read as “perfectly separated” for the same structural reason silhouette itself cannot see population share. Kernel PCA, t-SNE, and UMAP all avoid that failure, evidence that some nonlinearity is what matters, not PaCMAP specifically, but none balances separation and population balance as well PaCMAP does, 0.784 silhouette against 0.787 balance, against UMAP’s weaker 0.598 silhouette or t-SNE’s 0.431. PaCMAP’s design explicitly preserves mid-near point pairs alongside near and far ones, a deliberate balance between local and global structure UMAP and t-SNE’s local-neighbourhood-only objectives do not share [9], plausibly why it tolerates a smooth, continuous load-shape manifold better than either.

feature silhouette balance min_cluster_size
6 shape 0.784 0.787 302
2 profile_seasonal 0.713 0.266 58
0 profile 0.671 0.276 61
1 profile_weekly 0.651 0.272 60
5 trend 0.488 0.998 608
4 volatility 0.476 0.989 563
3 magnitude 0.452 1.000 626
8 calendar_peak 0.431 0.986 552
7 autocorr 0.409 0.887 391
9 distribution 0.393 0.957 486

shape wins clearly among feature sets too, ahead of the three profile-based sets and well ahead of the six single summary statistics. But it does not win everywhere. Reused against AusNet and GoiEner through the same shared ark.cluster.datasets loaders Part 5 Chapter 3’s generalisation check uses, this same fixed-context ablation finds a different winner on each population: autocorr on AusNet, profile_weekly on GoiEner, shape only on the London population the search itself was run against.

AusNet's own winner: autocorr
GoiEner's own winner: profile_weekly
feature silhouette min_cluster_size
7 autocorr 0.545 164
2 profile_seasonal 0.537 154
5 trend 0.532 139
0 profile 0.530 156
8 calendar_peak 0.525 154
1 profile_weekly 0.525 153
3 magnitude 0.447 113
9 distribution 0.431 97
4 volatility 0.417 161
6 shape 0.395 149

Whether shape wins because it is the right representation for behavioural load data in general, or only because London’s population happens to favour it, is exactly the question a single-dataset search cannot answer on its own. What each of the ten feature sets is actually built to ask is worth seeing directly, not just naming:

A natural next question is whether combining feature sets beats the single best one alone. It does not, at least not among the pairings tried here:

combination n_dims silhouette balance
0 shape 3 0.784 0.787
1 shape+volatility 5 0.613 0.738
2 shape+autocorr 8 0.631 0.885
3 shape+calendar_peak 7 0.589 0.854
4 shape+distribution 9 0.525 0.951

Every combination adds real dimensions and loses silhouette against shape alone, shape+distribution the worst of the four at 0.525 against 0.784. More features are not automatically better features; some of what volatility, calendar_peak, and distribution each capture on their own is redundant with, or actively dilutes, the same behavioural signal shape already carries cleanly.

scaler silhouette balance min_cluster_size
2 robust 0.784 0.787 302
4 power 0.769 0.780 297
3 quantile 0.769 0.272 60
1 minmax 0.532 0.988 560
0 standard 0.396 0.986 553

robust and power are close, both well ahead of minmax and standard. That gap has a checkable cause: shape’s feature columns carry a mean absolute skew of 12.6 and roughly a third of a percent of points sitting beyond three real standard deviations, a genuinely skewed, outlier-bearing distribution standard and minmax both let dominate the scaled feature space, and robust’s median-and-interquartile-range construction is built specifically to resist. quantile matches robust and power on silhouette but collapses on balance, 0.272 against 0.787, a reminder that a scaler can preserve separation while still distorting the population’s relative sizes.

clusterer silhouette balance min_cluster_size
0 kmeans 0.784 0.787 302
1 gmm 0.784 0.787 302
2 agg 0.784 0.787 302

The clustering algorithm, by contrast, makes no difference at all here: k-means, Gaussian mixture, and agglomerative clustering return the identical partition, identical silhouette, identical balance, identical smallest cluster. On a population this cleanly separated once shape, robust scaling, and PaCMAP have already done their work, there is only one reasonable two-way split left to find, and all three algorithms find it. This is the local-sensitivity point from earlier made concrete: fANOVA and PED-ANOVA both rank clusterer choice as mattering more than feature or scaler choice across the whole search, where plenty of trials sit on messier, less-separated embeddings; at this specific, already-settled operating point, it matters not at all.

Three algorithms agreeing is suggestive; it is not a test of Chapter 1’s four objective families, since k-means, Gaussian mixture, and agglomerative clustering are all partition- or hierarchy-based. Density-connectivity and graph-based clustering, the other two families Chapter 1 introduced, had never actually been run on this data. Both join the comparison here, at the same fixed feature, scaler, and embedding. DBSCAN’s radius comes from Ester, Kriegel, Sander, and Xu’s k-distance elbow, not a guess. Spectral clustering runs on a nearest-neighbour similarity graph, Ng, Jordan, and Weiss’s normalised cut, at the same \(K = 2\).

shape + robust + PaCMAP, k=2 (settled): eigengap heuristic's own suggested k (independent of the fixed k=2 above) = 11
condition objective_family clusters_found noise_points silhouette balance ari_vs_reference
0 shape + robust + PaCMAP, k=2 (settled) k-means 2 0 0.784 0.787 1.000
1 shape + robust + PaCMAP, k=2 (settled) Gaussian mixture 2 0 0.784 0.787 1.000
2 shape + robust + PaCMAP, k=2 (settled) Agglomerative 2 0 0.784 0.787 1.000
3 shape + robust + PaCMAP, k=2 (settled) DBSCAN (density-connectivity) 2 1 0.784 0.786 0.999
4 shape + robust + PaCMAP, k=2 (settled) Spectral (graph/normalised cut, k=2) 2 0 0.784 0.787 1.000

At this operating point, all five agree almost perfectly: DBSCAN’s two clusters match Gaussian mixture’s on all but a single household (Adjusted Rand Index (ARI) 0.999), and spectral clustering matches it exactly (ARI 1.000). Four objective families, tested for the first time on this data rather than assumed, unanimous. But this only shows that a population this well separated leaves no real disagreement for any family to have. It does not show the families are interchangeable in general. A second, less trivially separated condition below finds real disagreement between them.

None of the ablations above touch the one open question worth naming plainly rather than closing quietly: is \(K = 2\) the only real, resampling-trustworthy answer this search could have given, or did the cheap stability proxy simply never send anything else to the expensive real audit? The composite score’s top-5 finalists, every one of them \(K = 2\) or \(K = 3\), are not the only trials with real, decent balance and cheap stability at higher \(K\). Five candidates at \(K = 4\) through \(K = 10\), pulled from the full 400 completed trials rather than only the non-dominated Pareto front, share that same profile without ever reaching the real audit:

auditing 5 real k=4-10 candidates out of 388 completed trials
trial_number feature scaler reducer k silhouette balance cheap_stability
369 379 volatility power PaCMAP 4 0.465241 0.987722 0.925340
350 360 volatility minmax PaCMAP 10 0.454167 0.995144 0.896487
210 217 shape robust TSNE 6 0.445897 0.979646 0.887984
236 244 trend standard TSNE 5 0.393980 0.997911 0.884165
265 273 shape minmax KernelPCA 4 0.387714 0.860514 0.871961
trial_number feature reducer k prediction_strength min_cluster_stability mean_cluster_stability
0 379 volatility PaCMAP 4 0.893 0.953 0.972
1 360 volatility PaCMAP 10 0.541 0.847 0.904
2 217 shape TSNE 6 0.855 0.960 0.974
3 244 trend TSNE 5 0.505 0.657 0.807
4 273 shape KernelPCA 4 0.876 0.943 0.973

Three of the five pass both Tibshirani and Walther’s and Hennig’s thresholds just as convincingly as the settled recipe does: volatility+PaCMAP at \(K = 4\) (prediction strength 0.893, minimum cluster stability 0.953), shape+t-SNE at \(K = 6\) (0.855, 0.960), and shape+Kernel PCA at \(K = 4\) (0.876, 0.943). The other two do not, volatility+PaCMAP at \(K = 10\) recovering only a weak 0.541 prediction strength despite reasonable stability, trend+t-SNE at \(K = 5\) failing both. This is the honest answer to the open question just asked: \(K = 2\) is a genuine, resampling-validated answer, not the only one. Real four- and six-archetype structure, built from different feature and dimension-reduction pairings than the settled recipe, holds up under the same test just as well. The cheap proxy that drove the broad search never surfaced them for that test at all, not because they failed it, but because it never asked.

What do these richer splits actually look like as load curves, rather than as numbers in a table?

shape + t-SNE, K = 6: k=6, mean pairwise archetype-curve distance = 2.276, min = 0.971
volatility + PaCMAP, K = 4: k=4, mean pairwise archetype-curve distance = 0.509, min = 0.376
settled K=2 recipe: mean pairwise archetype-curve distance = 1.468

shape+t-SNE’s six archetypes are not a subtler version of the same coarse split. Their own mean pairwise distance between archetype curves, 2.276, is considerably larger than the settled \(K = 2\) recipe’s 1.468. That is a genuinely more differentiated set of behavioural rhythms, not a finer subdivision of the same two groups. volatility+PaCMAP’s four archetypes, by contrast, read as far less distinct in raw shape terms, 0.509. They were never built to separate on shape at all. They separate households by how erratic day-to-day demand is, a different question the shape-based view cannot see by construction.

Which representation is chosen does not just change the score. It changes what kind of structure the search is even capable of finding. Shape asks whether two households share a rhythm regardless of size, and finds a coarse two-way split or a sharper six-way one depending on how locally the embedding preserves neighbours. Volatility asks a different question, whether demand is routine or erratic, and finds a real but differently-shaped four-way answer neither shape view ever could.

A second, harder condition tests whether the four objective families still agree once the representation is not this well separated: the strongest real \(K = 4\) candidate above, volatility+PaCMAP, run through the same five families.

volatility + power + PaCMAP, k=4 (best k=4 candidate): eigengap heuristic's own suggested k (independent of the fixed k=4 above) = 11
condition objective_family clusters_found noise_points silhouette balance ari_vs_reference
0 volatility + power + PaCMAP, k=4 (best k=4 can... k-means 4 0 0.465 0.988 1.000
1 volatility + power + PaCMAP, k=4 (best k=4 can... Gaussian mixture 4 0 0.462 0.990 0.927
2 volatility + power + PaCMAP, k=4 (best k=4 can... Agglomerative 4 0 0.445 0.980 0.847
3 volatility + power + PaCMAP, k=4 (best k=4 can... DBSCAN (density-connectivity) 1 0 NaN 1.000 0.000
4 volatility + power + PaCMAP, k=4 (best k=4 can... Spectral (graph/normalised cut, k=4) 4 0 0.461 0.994 0.899

Here the families genuinely disagree. k-means, Gaussian mixture, and spectral clustering broadly agree with each other (ARI 0.85 to 1.00). DBSCAN collapses the entire population into a single cluster instead, ARI 0.000 against every other family. Density-connectivity’s notion of a cluster, a core point with enough neighbours within a fixed radius, simply does not recognise four separate dense regions here, the way a partition or a graph cut is built to assume one exists. This is Chapter 1’s “four ways to define a cluster” claim, demonstrated rather than only stated. Different families see different structure, not as a hypothetical, but concretely, on this book’s data, the moment the operating point stops being trivially easy.

The eigengap heuristic on the graph Laplacian’s spectrum is an entirely independent signal, one Optuna’s search never computed. It suggests \(K = 11\) at both operating points tried, a third, separate piece of evidence pointing toward the same richer-structure conclusion as the resampling audit above. Worth naming honestly as a further open question, not a settled third answer.

The honest next step these findings call for is a second, expanded search. It would score every Pareto-optimal trial’s cheap stability against this same audit, not just the composite ranking’s top five, and treat DBSCAN’s radius and spectral clustering’s eigengap as first-class search dimensions, not a fixed-context afterthought. This chapter’s single search has not already settled everything worth knowing.

How many components to keep

Once a dimension-reduction technique that keeps more than two dimensions is in play, a question Chapter 1 raised but never settled comes back: how many components to keep. This book’s existing convention elsewhere is a fixed rule, keep enough components to reach 90% of cumulative explained variance, a threshold chosen for convenience rather than derived from the data’s structure. Zhu and Ghodsi’s profile-likelihood method offers a principled alternative: model the top \(q\) eigenvalues and the remaining \(p - q\) as two Gaussian groups sharing one pooled variance, and choose the \(q\) that maximises the resulting profile log-likelihood [13], \[ \ell(q) = \sum_{i=1}^{q} \log \mathcal{N}(\lambda_i \mid \mu_1, \sigma^2) + \sum_{i=q+1}^{p} \log \mathcal{N}(\lambda_i \mid \mu_2, \sigma^2), \] where \(\lambda_1 \geq \dots \geq \lambda_p\) are the ordered eigenvalues, \(\mu_1, \mu_2\) are each group’s mean, and \(\sigma^2\) is pooled across both so that the profile likelihood cannot diverge by shrinking one group’s variance to zero. The method finds an automatic change point in the spectrum itself, with no threshold to tune.

Tried on four feature sets, chosen to span native dimensionality from 6 to 264 columns, the two methods disagree on every one, and not narrowly:

feature n_features n_components_variance_threshold n_components_profile_likelihood agree variance_threshold_pct_of_features profile_likelihood_pct_of_features
0 profile 264 14 1 False 0.053030 0.003788
1 profile_weekly 168 13 1 False 0.077381 0.005952
2 profile_seasonal 96 11 1 False 0.114583 0.010417
3 distribution 6 3 1 False 0.500000 0.166667

Profile-likelihood recommends a single component for all four, including distribution, the smallest and least redundant of the four sets. This is not the method failing; it is answering a different question. These load-shape spectra share one dominant direction, essentially how large and how peaky a household’s overall demand is, followed by a long, smoothly decaying tail that a shared-variance Gaussian model correctly treats as one homogeneous group rather than a second, confidently separable signal. Reaching 90% cumulative variance still costs 11 to 14 components on the three profile sets, each contributing a small, gradually diminishing slice of a fixed, somewhat arbitrary target, not a second structure in the data. The practical stakes are not abstract: the best clustering results below come from two-component embeddings, far closer to profile-likelihood’s recommendation than the variance rule’s.

When separation lies

Chapter 1 already gave silhouette, Davies-Bouldin, and Calinski-Harabasz their own formal definitions [14], [15], [16]. All three share the same structural blind spot: each is a function of distances or variances alone, with no term that penalises one group holding almost none of the population.

Figure 12.1: Same formulas, same blind spot. A real, even split scores 0.78; MAC000037’s degenerate split, where two isolated households make up barely one percent of the population between them, scores higher: 0.99. Neither silhouette nor Davies-Bouldin nor Calinski-Harabasz ever looks at how many points are actually in a group.

MAC000037’s cluster satisfies \(b(i) \gg a(i)\) trivially, since nothing else in London is remotely close to it, and every one of the three formulas reports that trivial fact as a triumph: silhouette 0.994, Davies-Bouldin 0.002 (lower is meant to be better), Calinski-Harabasz 3.18 million (higher is meant to be better). Three formulas, three different scales, one household, and all three call it a win.

Common Mistake

Checking multiple internal validity indices and treating their agreement as confirmation. Silhouette, Davies-Bouldin, and Calinski-Harabasz are three different formulas built from the same two ingredients, within-cluster compactness and between-cluster separation, so they agree precisely where all three are equally blind: an isolated outlier, a single household whose own recorded behaviour is a data artefact rather than a genuine pattern, satisfies every separation-based definition of a good cluster at once. Agreement among these three is evidence they were computed from the same clustering, not evidence the clustering found something real.

The fix needs a different kind of evidence, not one more formula built from the same two ingredients. Chapter 1’s balance measure, entropy over the cluster-size distribution normalised to \([0, 1]\), at least notices the shape of the problem [17]. MAC000037’s \([1282, 1, 1]\) split scores a balance of 0.009, about as far from even as three clusters can get. Balance is cheap, needs no resampling, and catches this specific failure, but it is a symptom check, not proof of realness, exactly the limit Chapter 1 already named: a population could, in principle, contain a rare archetype that also happens to be a small cluster. Telling the two apart needs a different kind of evidence entirely.

Validating without labels

Chapter 1’s “Checking stability” section introduced Strehl and Ghosh’s ensemble approach, building a consensus matrix from several independent clustering runs of the same data [18]. That answers “does this algorithm keep finding the same structure,” a useful question, and a different one from “does this structure survive at all if the data itself had come out slightly differently.” This chapter needs the second, complementary notion, because an outlier-driven cluster can be perfectly reproducible across repeated runs on the same fixed dataset while still being entirely an artefact of that one household’s presence.

Figure 12.2: Two ways to ask whether a clustering survives data coming out slightly differently. Prediction strength splits the data once and checks whether one half’s centroids still predict the other half’s pairs correctly. Cluster-wise stability bootstraps the whole dataset and checks how well each original cluster survives being found again. The panel below reads each method’s verdict on this book’s data: prediction strength clears Tibshirani and Walther’s threshold comfortably for London and falls just short of it for AusNet, while cluster-wise stability reaches 1.000 for a genuine archetype and lands at 0.63, inside Hennig’s “patterned, not stable” band, for MAC000037’s isolated cluster.

Chapter 1 already gave both resampling-based methods that answer this second question their own formal definitions: Tibshirani and Walther’s prediction strength [19], and Hennig’s cluster-wise stability just below. Neither, as far as this book’s literature search found, has been applied to smart-meter customer archetyping before. As a brief reminder, prediction strength splits the data into two random halves \(A\) and \(B\), fits centroids on \(A\)’s clustering, and checks whether those centroids still predict \(B\)’s same-cluster pairs correctly. \[ \text{ps}(C) = \frac{1}{n_C(n_C - 1)} \sum_{\substack{i, i' \in C \\ i \neq i'}} \mathbb{1}\!\left[\text{pred}(i) = \text{pred}(i')\right] \] Prediction strength for \(K\) is the minimum of \(\text{ps}(C)\) over \(B\)’s clusters, and the recommendation is to pick the largest \(K\) scoring above 0.8 to 0.9. A genuine limitation is worth stating plainly rather than glossed over: a cluster with fewer than two members contributes no pairs, so a singleton is excluded from its own split’s minimum rather than penalised by it. Prediction strength validates whether \(K\) itself is well supported globally; it is not, by itself, guaranteed to flag one small outlier cluster hiding inside an otherwise well-supported partition.

Key Concept

Prediction strength asks one checkable question: if this population had come out slightly differently, would a clustering fit on half of it still make the same calls on the other half’s pairs? A \(K\) that only holds together on the exact rows it was fit against is not a customer segment, it is a coincidence of that particular sample. The score is reported once per \(K\), a judgement about whether that many groups is safe to trust at all, before any single cluster inside that partition gets scrutinised on its own.

This is the gap Hennig’s cluster-wise stability closes [20], bootstrap-resampling the data, refitting centroids, and checking the Jaccard overlap between each original cluster and its best match among the resample’s predicted clusters. \[ J(C, C') = \frac{|C \cap C'|}{|C \cup C'|}, \qquad \text{stability}(C) = \frac{1}{B}\sum_{b=1}^{B} \max_{C'_b} J(C, C'_b) \] A real archetype’s centroid is recovered reliably even from a partial resample, reclaiming essentially all of its original members every time, stability near 1.0. A cluster built around one extreme household only has a centroid to speak of in resamples where that household happened to be drawn, and is silently absorbed into whichever cluster is nearest otherwise. Hennig’s interpretation bands read stability below 0.5 as dissolved, 0.6 to 0.75 as patterned but not stable, and above 0.75 as stable. Checked against MAC000037’s cluster, stability lands at 0.63, inside the “patterned but not stable” band, not near zero: whenever that household is present in a resample it isolates itself again, so its score settles close to \(1 - 1/e \approx 0.632\), the plain probability that a single point survives one bootstrap draw, a value this book’s test suite for the method reproduces analytically before trusting the result on real data.

Key Concept

Cluster-wise stability asks the question prediction strength cannot: not whether \(K\) as a whole is safe to trust, but whether this specific cluster survives being found again. Bootstrap-resample the population, re-cluster, and check how much of the original cluster comes back. A real archetype’s centroid has hundreds of households anchoring it, so it comes back almost every time. A cluster built around one household only exists in the resamples lucky enough to draw that one household, and Hennig’s interpretation bands turn that difference into a number a practitioner can act on rather than eyeball.

A trust-gated composite score

Chapter 1’s closing warning, that no single validity index is universally correct, needs a combination rule, not just a longer list of indices to eyeball together [1]. The naive fix, averaging separation and trust into one additive score, does not work: MAC000037’s combination scores a near-perfect separation ensemble and a trust factor of 0.009, and a flat average lets the separation term partially buy the trust deficit back, the kind of compensation a workable decision rule for this problem cannot allow.

Figure 12.3: Two illustrative candidates, one with excellent separation and almost no trust, one with moderate separation and real trust. An additive average lets the first buy back its own trust deficit and wins outright, 0.51 against 0.47. A multiplicative gate, separation times trust, correctly prefers the trustworthy candidate instead, 0.23 against 0.02.

Outside clustering, this is a well-studied problem: multi-criteria decision analysis calls it a noncompensatory decision rule. An alternative that fails badly enough on one criterion gets excluded, or vetoed, no matter how well it performs on everything else; it never gets the chance to trade one off against the rest [21]. Roy’s outranking theory, the foundation of the Elimination and Choice Expressing Reality (ELECTRE) family of methods, formalises this with a veto threshold: a candidate whose performance on one criterion falls far enough behind is excluded from consideration outright [22]. The composite score below is a direct, continuous version of that same principle, not an additive vote: \[ \text{composite}(x) = S(x) \times T(x), \] where \(S(x)\) is a diversity-weighted rank ensemble over silhouette, Calinski-Harabasz, and Davies-Bouldin, and \(T(x) = \min\bigl(\text{balance}(x), \text{stability}(x)\bigr)\) is a trust factor built from whichever of balance or stability is weaker for that candidate, a Rawlsian “weakest link” reading rather than an average of the two. The diversity weighting itself follows the same logic Chapter 1’s choice of index already implied: metrics that are highly correlated with each other across the candidate pool contribute redundant information and are down-weighted, so a near-unanimous cluster of agreeing separation metrics cannot outvote whichever one is actually saying something different. Multiplication, not addition, is the load-bearing choice: \(S(x)\) near its maximum cannot rescue a \(T(x)\) near zero, since the product of a large number and a number near zero is itself near zero, the same veto behaviour Roy’s theory formalises and a flat weighted sum structurally cannot provide.

Key Concept

The composite score folds three separation metrics and two trust signals into one number a search can rank candidates by, without letting either side hide the other. High separation cannot buy back low trust, and high trust alone is not enough if nothing about the split is actually separated. A candidate only reaches the top of this ranking by being both well separated and safe to believe, the same standard a human reviewer would apply by eye if there were time to check every trial by hand.

Ranking by raw silhouette alone would shortlist MAC000037’s degenerate trial, and two others just like it, for further, expensive scrutiny, the exact bias this whole approach exists to avoid. Ranking by composite score instead, on the same Pareto-optimal trials, drops every outlier-isolation case to the bottom of the list. Five balanced candidates rise instead, each with a real multi-hundred-household minority cluster, forming a shortlist worth the expensive, well-powered resampling audit: ten-split prediction strength and hundred-bootstrap cluster stability. Those checks run once per shortlisted candidate, not on every trial, since the cheap versions already did the work of finding a fair shortlist to spend that cost on.

Every trial silhouette alone would send to the expensive audit has a smallest cluster of at most a few dozen households, three of the five a literal single point; every trial the composite score sends instead has a smallest cluster in the hundreds. The same Pareto-optimal trials, read through two different rules, produce two shortlists that share not one candidate.

The gate has a real seam, though, worth naming rather than papering over. Taking the minimum of balance and stability means a genuinely rare archetype, an early-adopter segment that only ever amounts to a few percent of a population, gets vetoed by the same low balance reading a fake one would, no matter how stable it turns out to be. Nothing in these three populations ever needed the distinction, since neither London’s nor AusNet’s nor GoiEner’s real split was ever small. A population with a real minority segment worth finding on purpose would need trust read off stability alone.

Choosing K

The fANOVA decomposition earlier in this chapter already answers this, once the question is asked the right way: \(K\) dominates the variance in this search’s cheap stability signal, more than every other parameter combined, and explains almost nothing of silhouette or balance. Chapter 1’s “Choosing K” section warned against the elbow method on distortion, since distortion falls monotonically by construction and offers no signal about when to stop [1]. This finding sharpens that warning into something a practitioner can act on: whichever dimension-reduction technique and feature set are chosen mostly determine how separated a clustering looks, but \(K\) itself mostly determines whether that clustering is trustworthy under resampling, a different axis a separation-only view, however many indices it checks, cannot see.

Several real answers, not one

The point of this whole search was never to crown a single winner. It was to understand how the choices this chapter opened with, feature set, scaler, dimension-reduction technique, clusterer, actually shape what structure a search finds. The honest answer this section states plainly: several genuinely different, independently validated answers, not one. The audit’s verdict on the five composite-score finalists is still worth stating first. It remains the single strongest result on its own terms.

trial_number feature scaler reducer clusterer k silhouette calinski_harabasz davies_bouldin balance prediction_strength prediction_strength_std min_cluster_stability mean_cluster_stability separation_score trust_factor composite_score
0 343 shape robust PaCMAP gmm 2 0.784203 7204.961816 0.272857 0.786971 1.000000 0.000000 1.000000 1.000000 1.000000 0.786971 0.786971
1 399 volatility minmax PaCMAP kmeans 2 0.608491 3311.485168 0.558540 0.988487 0.967819 0.021060 0.985515 0.987053 0.718519 0.985515 0.708111
2 298 volatility minmax TSNE kmeans 2 0.588867 3138.718700 0.590206 0.995794 0.978515 0.019673 0.989180 0.989938 0.518519 0.989180 0.512908
3 150 shape minmax PaCMAP kmeans 3 0.570801 4432.848835 0.603476 0.998915 0.905078 0.032414 0.965516 0.976168 0.444444 0.965516 0.429118
4 251 shape minmax PaCMAP kmeans 2 0.575562 2929.900183 0.596423 0.998207 0.915511 0.042446 0.980830 0.981724 0.318519 0.980830 0.312412

shape features (peak-normalised daily demand, the same convention established for a different reason in Chapter 1), a robust scaler, PaCMAP, and Gaussian mixture clustering at \(K = 2\) wins the composite ranking outright: a prediction strength of 1.000 with a standard deviation of 0.000 across all ten held-out splits, every split agreeing completely, and a minimum bootstrap cluster stability of 1.000, both clusters recovered in every one of the hundred resamples.

cluster sizes: [982 302]

The two clusters hold 982 and 302 of London’s 1,284 households, a sharper evening-peaking rhythm against a flatter, lower one. But the earlier ablation section already found real, differently-shaped structure this composite ranking never surfaced. It never surfaced them only because its own cheap stability proxy never sent them to the expensive real audit. volatility+PaCMAP at \(K = 4\), shape+t-SNE at \(K = 6\), and shape+Kernel PCA at \(K = 4\) all pass the same real resampling thresholds just as convincingly. Each reveals a different kind of behavioural structure: a routine-versus-erratic split rather than a shape split, or a sharper six-way shape split whose own archetypes are more differentiated from each other than the settled two-way split is. None of this makes \(K = 2\) wrong. It makes \(K = 2\) one well-supported answer among several, the one this search’s cheap proxy happened to rank first, not the only structure this population actually contains.

A settled recipe that only works on the dataset it was tuned against is not settled at all. So the single strongest candidate, the composite ranking’s \(K = 2\) recipe, was reconstructed with its identical feature set, scaler, dimension-reduction technique and its own winning hyperparameters, clusterer, and \(K\). It was then applied unchanged to AusNet’s 342 households and an 800-household sample from GoiEner, a second, independently metered Spanish utility. The richer \(K = 4\) to \(K = 6\) alternatives above are real, validated findings on London specifically. They have not yet been checked the same way against AusNet or GoiEner, an honest scope boundary worth stating rather than implying a generality this chapter has not actually earned yet.

dataset n_households silhouette balance prediction_strength min_cluster_stability mean_cluster_stability cluster_sizes
0 London 1284 0.784203 0.786971 1.000000 1.000000 1.000000 [982, 302]
1 AusNet 342 0.429859 0.977683 0.749357 0.844909 0.859937 [201, 141]
2 GoiEner 800 0.482046 0.999116 0.860307 0.913367 0.915853 [414, 386]

Every population clears Hennig’s 0.75 stability bar on its weakest cluster, AusNet and GoiEner comfortably so, and every split is balanced rather than one dominant group and stragglers. AusNet’s prediction strength, 0.749, falls just under Tibshirani and Walther’s suggested 0.8 to 0.9 range, worth reporting as it stands rather than rounding favourably; with only 342 households to begin with, AusNet has the least data of the three to support a resampling-based check, so a slightly weaker reading here is an expected consequence of sample size, not a sign the recipe itself fails to travel. A number this close to a published threshold is a call for more splits, not a verdict either way. Silhouette is the one number that does not transfer cleanly, London’s 0.784 well above AusNet’s 0.430 or GoiEner’s 0.482, and balance explains why: AusNet and GoiEner’s splits land close to even, 59/41 and 52/48, while London’s is more skewed at 76/24, and a more even split leaves less silhouette headroom between two group centroids relative to their own internal spread by construction. The same recipe finds a stable, non-degenerate two-way behavioural split in each of the three independent populations tried, differently proportioned in every case, which is what a settled methodology has to show before the next chapter can treat it as decided rather than still being tested.

Why bother

It would have been faster to stop at the very first search this chapter ran, publish the 0.994-silhouette combination, and move on. Every later chapter in this part would then have inherited MAC000037’s cluster as “the” London archetype, and the same failure mode would have waited quietly inside AusNet’s and GoiEner’s searches too, since nothing about a single broken feature is unique to one dataset.

For anyone applying clustering to a new population, the trust-gated composite score and the two resampling checks behind it are not specific to smart-meter data. Nothing in prediction strength, cluster-wise stability, or a multiplicative trust gate assumes load curves, households, or kilowatt-hours. The same recipe applies wherever a clustering has to be trusted without ground truth, anywhere a real outlier risks quietly winning the search.

For the next chapter, this one hands over a validated evaluation method, not a fixed recipe to copy unchanged. Whatever number of archetypes Chapter 3 settles on for AusNet, for its own operational reasons rather than London’s, it does not have to guess whether that choice is safe. It can check it against the exact same real audit built here, prediction strength and cluster stability, and report the honest result, whichever way it lands.

For anyone who has to defend a clustering result to someone else, a regulator, a colleague, a reviewer, the honest answer to “how do you know this is real” stops being “the silhouette was high.” It becomes “prediction strength was 1.000 across ten independent splits, and every cluster survived a hundred bootstrap resamples,” a claim that can be checked, not just asserted.

Where this leads

None of the four questions this chapter opened with had a knowable answer in advance, and pretending otherwise, publishing whichever combination scored best on the first metric checked, is how MAC000037 becomes an archetype instead of a footnote. What settles here is not one fixed number of archetypes. It is a validated method for finding trustworthy structure at all: a resampling-validated dimensionality choice, a trust-gated composite score standing in for a single validity index. Applied honestly, that method surfaced several real, differently-shaped answers, not only the first one checked. The clearest of them, load-shape features, PaCMAP, and Gaussian mixture clustering at \(K = 2\), is the one checked once against a second internal signal and then against two independent utilities. It is not the only recipe the next chapter has to use. It is the standard the next chapter’s choices get held to. Chapter 3 builds AusNet archetypes with its own pipeline, compares them against a deep-clustering alternative, calibrates confidence in each customer’s membership, and checks against an independently labelled population that Chapter 1 already flagged as a genuine limit worth reporting honestly rather than smoothing over.

References