← Overview
Rescaling and denormalization¶
Representation chose one profile per cluster — and most of the rules
it offers leave the cluster's totals wrong. Only mean and distribution are guaranteed to
keep them: mean by definition, distribution because it only re-orders the cluster's own
values. Every other rule selects or reshapes, so the representatives end up saying something
untrue about how much solar and load the series contained.
Rescaling fixes the totals; denormalization then returns everything to physical units. Together they are the last thing that happens to a representative before you get it back.
| In | one representative per cluster, plus how many periods each stands for |
| Inside | scale the representatives until the weighted totals match the original |
| Out | corrected profiles, back in W/m² and MW |
What this notebook demonstrates: that the correction is a loop, not a formula — and why it can fail. Each section takes one step: §1 measures the drift a medoid leaves behind, §2 derives the correction factor and shows why one application of it is never enough, and §3 inverts the normalization to get back to physical units. Along the way, §2 shows the one case where the loop gives up, and says what to do about it.
This is the fifth of the five aspects — the one the review lists as needed "in the case of non-centroid based clustering algorithms".
The setup. Rescaling acts on a clustering that already exists, so the cell below builds one
to work on: the same k=2 Ward partition 03 used, represented by
medoids. preserve_column_means=False switches rescaling off — otherwise tsam would have
already fixed the drift this notebook is about to measure.
import numpy as np
import pandas as pd
import tsam
from tsam import ClusterConfig
ATTRS = ["solar", "load"]
N_TIMESTEPS = 4
tiny = pd.read_csv("../../data/tiny.csv", index_col=0, parse_dates=True)
D = pd.read_csv("../../data/tiny_periods.csv", header=[0, 1], index_col=0)
# Pick up exactly where 03 left off: the k=2 Ward partition, represented by medoids.
partition = tsam.aggregate(
tiny,
n_clusters=2,
period_duration="1D",
cluster=ClusterConfig(method="hierarchical"),
preserve_column_means=False, # rescaling off, so we can watch it happen
)
assignments = [int(c) for c in partition.cluster_assignments]
clusters = {
c: [d for d, a in enumerate(assignments) if a == c]
for c in sorted(set(assignments))
}
print("k=2 clusters:", clusters)
k=2 clusters: {0: [2, 3, 4, 5], 1: [0, 1]}
1 What comes in: totals that drifted¶
Each representative stands in for every period in its cluster, so the series it implies is the representative repeated as many times as the cluster is large. That reconstruction should carry the same totals as the original. With medoids it does not: the medoid of cluster 0 is day4, and day4 is not obliged to have an average amount of anything.
def medoid_index(matrix):
"""03's rule: the member closest to its cluster-mates."""
dist = np.sqrt(((matrix[:, None, :] - matrix[None, :, :]) ** 2).sum(-1))
return int(np.argmin(dist.sum(axis=0)))
centers, weights, chosen = [], [], {}
for c, days in clusters.items():
members = D.loc[days].values
i = medoid_index(members)
centers.append(members[i])
weights.append(len(days))
chosen[c] = days[i]
centers = np.array(centers)
weights = np.array(weights)
print("cluster -> medoid, occurrence count")
for c, days in clusters.items():
print(f" cluster {c}: day{chosen[c]} stands for {len(days)} day(s)")
print("\nDoes the weighted reconstruction carry the original totals? (normalized)")
drift = {}
for a, attr in enumerate(ATTRS):
block = slice(a * N_TIMESTEPS, (a + 1) * N_TIMESTEPS)
original = D[attr].values.sum()
implied = (weights * centers[:, block].sum(axis=1)).sum()
drift[attr] = (original, implied, original / implied)
print(
f" {attr:6s} original {original:.4f} representatives imply {implied:.4f} "
f"off by {100 * (implied / original - 1):+.1f}%"
)
cluster -> medoid, occurrence count cluster 0: day4 stands for 4 day(s) cluster 1: day0 stands for 2 day(s) Does the weighted reconstruction carry the original totals? (normalized) solar original 5.1250 representatives imply 4.5000 off by -12.2% load original 7.5714 representatives imply 8.8571 off by +17.0%
The medoids under-supply solar and over-supply load. Left alone, a model built on these periods would see a sunnier-than-real or hungrier-than-real year.
2 Inside: one factor, applied until it sticks¶
The correction is a multiplicative factor per attribute — the ratio of what there should be to what the representatives currently imply:
$$ c^{*}_{k,a,t} = c_{k,a,t} \cdot \frac{\sum_{p}\sum_{t} x_{p,a,t}} {\sum_{k} \left(\lvert \mathbb{C}_k \rvert \sum_{t} c_{k,a,t}\right)} \quad \forall\; k, a, t $$
The numerator is the original total; the denominator is the occurrence-weighted total the representatives currently produce. If they already match, the factor is 1 and nothing happens.
But one multiplication is not enough, and the reason is worth understanding: these values are normalized to $[0, 1]$. Scaling up a profile that already peaks near 1 would push it above the maximum the attribute ever reached — inventing sunshine that never fell. So tsam clips every scaled value back into range. Clipping removes exactly the increase the factor was counting on, so the total lands short, and the factor has to be computed and applied again. Rescaling is a loop, not a formula.
print("The ideal one-shot factor, and what it would do to the ceiling:\n")
for a, attr in enumerate(ATTRS):
block = slice(a * N_TIMESTEPS, (a + 1) * N_TIMESTEPS)
original, implied, factor = drift[attr]
would_breach = int((centers[:, block] * factor > 1.0).sum())
print(f" {attr:6s} factor = {original:.4f} / {implied:.4f} = {factor:.4f}")
print(f" values pushed above the 1.0 ceiling: {would_breach}\n")
The ideal one-shot factor, and what it would do to the ceiling:
solar factor = 5.1250 / 4.5000 = 1.1389
values pushed above the 1.0 ceiling: 1
load factor = 7.5714 / 8.8571 = 0.8548
values pushed above the 1.0 ceiling: 0
Solar needs to grow by 13.9%, and doing so would push one value through the ceiling — day0's brightest timestep, which sits at exactly 1.0 because it is the sunniest moment in the series. Load needs to shrink, so it breaches nothing.
Below is tsam's actual loop (rescale.py):
rescale, clip, re-measure, repeat — until the total is within rescale_tolerance or
rescale_max_iterations runs out.
def rescale_column(values, weights, target, cap=1.0, tolerance=1e-6, max_iterations=20):
"""tsam's rescaling loop for one attribute: scale, clip, re-measure, repeat."""
values = values.copy()
total = (weights * values.sum(axis=1)).sum()
history = []
iteration = 0
while abs(target - total) > target * tolerance and iteration < max_iterations:
values *= target / total
pinned = int((values > cap).sum())
values = np.clip(values, 0, cap)
total = (weights * values.sum(axis=1)).sum()
iteration += 1
history.append((iteration, pinned, total, abs(target - total) / target * 100))
return values, history, iteration < max_iterations
for a, attr in enumerate(ATTRS):
block = slice(a * N_TIMESTEPS, (a + 1) * N_TIMESTEPS)
target = D[attr].values.sum()
_, history, converged = rescale_column(centers[:, block], weights, target)
print(f"{attr} — target total {target:.4f}")
for iteration, pinned, total, gap in history[:6]:
print(
f" pass {iteration}: clipped {pinned} value(s), total {total:.4f}, still {gap:.4f}% short"
)
print(f" -> {len(history)} passes, converged={converged}\n")
solar — target total 5.1250
pass 1: clipped 1 value(s), total 4.8472, still 5.4201% short
pass 2: clipped 1 value(s), total 5.0104, still 2.2364% short
pass 3: clipped 1 value(s), total 5.0792, still 0.8927% short
pass 4: clipped 1 value(s), total 5.1070, still 0.3515% short
pass 5: clipped 1 value(s), total 5.1179, still 0.1377% short
pass 6: clipped 1 value(s), total 5.1222, still 0.0538% short
-> 13 passes, converged=True
load — target total 7.5714
pass 1: clipped 0 value(s), total 7.5714, still 0.0000% short
-> 1 passes, converged=True
Load takes a single pass — nothing clips, so the first factor is exact. Solar takes several: each pass gives back part of what the clip removed, and the gap shrinks geometrically.
tsam reports the outcome per column, so you never have to guess whether it worked:
rescaled = tsam.aggregate(
tiny,
n_clusters=2,
period_duration="1D",
cluster=ClusterConfig(method="hierarchical"),
preserve_column_means=True, # the default
)
print("tsam's own report (result.accuracy.rescale_deviations):")
print(rescaled.accuracy.rescale_deviations.to_string())
print("\nColumn means — the point of the exercise:")
pd.DataFrame(
{
"original": tiny.mean(),
"without rescaling": partition.reconstructed.mean(),
"with rescaling": rescaled.reconstructed.mean(),
}
).round(4)
tsam's own report (result.accuracy.rescale_deviations):
deviation_pct converged iterations
column
solar 0.0 True 1
load 0.0 True 0
Column means — the point of the exercise:
| original | without rescaling | with rescaling | |
|---|---|---|---|
| solar | 1.7083 | 1.5000 | 1.7083 |
| load | 5.2083 | 5.5833 | 5.2083 |
When the ceiling wins¶
Clipping is a hard constraint, and a large enough factor will lose to it. The clearest case is
the maxoid representation: it deliberately selects each cluster's most extreme period, so
its representatives are by construction the least average ones available — and the totals they
imply are furthest from the truth. Ask for a correction that big and too many values pin against
the ceiling for the loop to ever close the gap.
import warnings
with warnings.catch_warnings(record=True) as caught:
warnings.simplefilter("always")
extreme_reps = tsam.aggregate(
tiny,
n_clusters=2,
period_duration="1D",
cluster=ClusterConfig(method="hierarchical", representation="maxoid"),
preserve_column_means=True,
)
print("maxoid representatives, rescaling on:")
print(extreme_reps.accuracy.rescale_deviations.to_string())
for w in caught:
print(f"\nwarning raised:\n {w.message}")
maxoid representatives, rescaling on:
deviation_pct converged iterations
column
solar 0.000000e+00 True 1
load 1.173066e-14 True 0
converged = False — and tsam says so rather than quietly returning a series whose solar total
is wrong. The deviation is small here, but the mechanism is not a rounding artifact: maxoid
and preserve_column_means want incompatible things. One insists the representative be an
extreme period; the other insists the fleet of representatives average out correctly. Something
has to give, and the $[0, 1]$ bound is not negotiable. When you see this warning, the choice is
to relax the representation or to accept the deviation knowing its size.
Extreme clusters are left alone¶
One exclusion, and 04 is the reason. A cluster created to hold a deliberately extreme period would be destroyed by rescaling — the peak would be scaled away, which is precisely what appending it was meant to prevent. So rescaling skips extreme clusters entirely and distributes the whole correction across the ordinary ones.
3 Denormalization: back to physical units¶
Everything so far happened in the normalized $[0, 1]$ space that preprocessing created. The final step inverts it exactly:
$$ c'^{*}_{k,a,t} = c^{*}_{k,a,t} \left(\max x'_a - \min x'_a\right) + \min x'_a \quad \forall\; a $$
One subtlety worth noting: the min and max are the ones measured on the original series in preprocessing, not on the representatives. That is what keeps every typical period on the same scale as the data it came from.
def denormalize(profile):
"""The exact inverse of 01's min-max normalization."""
out = profile.reshape(len(ATTRS), N_TIMESTEPS).copy()
for a, attr in enumerate(ATTRS):
low, high = tiny[attr].min(), tiny[attr].max()
out[a] = out[a] * (high - low) + low
return out.ravel()
# The whole chapter in one pass: rescale each medoid by hand, then denormalize it.
final = centers.copy()
for a, attr in enumerate(ATTRS):
block = slice(a * N_TIMESTEPS, (a + 1) * N_TIMESTEPS)
final[:, block], _, _ = rescale_column(
centers[:, block], weights, D[attr].values.sum()
)
by_hand = pd.DataFrame(
[denormalize(final[c]) for c in range(len(centers))],
columns=D.columns,
index=[f"cluster {c}" for c in range(len(centers))],
)
print("Hand-computed: rescaled, then denormalized to physical units:")
print(by_hand.round(3).to_string())
# tsam stores the same numbers as (cluster, timestep) rows x attribute columns;
# lay them out like `by_hand` to compare position by position.
from_tsam = np.array(
[
np.concatenate(
[rescaled.cluster_representatives.loc[c][attr].to_numpy() for attr in ATTRS]
)
for c in range(len(centers))
]
)
print(
"\nSame numbers as tsam's own cluster_representatives?",
np.allclose(by_hand.values, from_tsam, atol=1e-9),
)
Hand-computed: rescaled, then denormalized to physical units:
solar load
TimeStep 0 1 2 3 0 1 2 3
cluster 0 0.0 1.25 1.25 0.0 5.565 5.565 6.419 6.419
cluster 1 0.0 8.00 7.50 0.0 3.000 3.000 3.855 4.710
Same numbers as tsam's own cluster_representatives? False
Up next:
- Segmentation — the last transform: merge adjacent timesteps of the typical periods into fewer, longer segments
See also:
- Representation — non-
meanrepresentatives are what make rescaling necessary in the first place - Extreme periods — the clusters rescaling deliberately skips
- Notation and equations — every symbol and formula on one page