Pareto Optimization

Determine the optimal combination of segments and periods for time series aggregation.

Author: Leander Kotzur

from pathlib import Path

import numpy as np
import pandas as pd
import plotly.express as px
import plotly.io as pio

import tsam
from tsam import ClusterConfig
from tsam.tuning import find_pareto_front

pio.renderers.default = "notebook"

# Ensure results directory exists
RESULTS_DIR = Path("results")
RESULTS_DIR.mkdir(exist_ok=True)

Input data

Read in time series from testdata.csv with pandas

raw = pd.read_csv("testdata.csv", index_col=0)
raw = raw.rename(
    columns={"T": "Temperature", "Load": "Demand", "Wind": "Wind", "GHI": "Solar"}
)
period_duration = 24

Plot the original data

# Original data heatmaps using tsam.unstack_to_periods() with plotly
unstacked = tsam.unstack_to_periods(raw, period_duration=period_duration)
for col in raw.columns:
    px.imshow(
        unstacked[col].values.T,
        labels={"x": "Day", "y": "Hour", "color": col},
        title=f"Original {col}",
        aspect="auto",
    ).show()

Find Pareto-optimal aggregations

Use find_pareto_front() to explore the Pareto-optimal combinations of periods and segments.

pareto_results = find_pareto_front(
    raw,
    period_duration=period_duration,
    timesteps=np.geomspace(5, 8760, 50).astype(int).tolist(),
    cluster=ClusterConfig(method="hierarchical", representation="distribution"),
    n_jobs=-1,
)

Testing configurations (2 workers): 100%|██████████| 141/141 [00:27<00:00,  5.10it/s]

Visualize the Pareto front - the trade-off between compression and accuracy.

pareto_df = pd.DataFrame(
    [
        {
            "timesteps": r.n_clusters * r.n_segments,
            "periods": r.n_clusters,
            "segments": r.n_segments,
            "rmse": r.accuracy.rmse.mean(),
        }
        for r in pareto_results
    ]
)

fig = px.line(
    pareto_df,
    x="timesteps",
    y="rmse",
    markers=True,
    labels={"timesteps": "Timesteps (periods x segments)", "rmse": "RMSE"},
    title="Pareto Front: Compression vs Accuracy",
    hover_data=["periods", "segments"],
    range_y=(0, None),
)
fig.show()

Show the final result

last_result = pareto_results[-1]
print(
    f"Final: {last_result.n_clusters} periods, {last_result.n_segments} segments, RMSE: {last_result.accuracy.rmse.mean():.4f}"
)

Final: 365 periods, 24 segments, RMSE: 0.1535

# Reconstructed data heatmaps
reconstructed = last_result.reconstructed
unstacked_recon = tsam.unstack_to_periods(
    reconstructed, period_duration=period_duration
)
for col in reconstructed.columns:
    px.imshow(
        unstacked_recon[col].values.T,
        labels={"x": "Day", "y": "Hour", "color": col},
        title=f"Reconstructed {col}",
        aspect="auto",
    ).show()

Animated visualization

Animate through all Pareto-optimal aggregations to visualize the trade-off between compression and accuracy.

n_days = len(raw) // period_duration
n_vars = len(raw.columns)

# Get normalization parameters from original data
raw_min = raw.min()
raw_range = raw.max() - raw.min()

frames_data, labels = [], []
for result in reversed(pareto_results):
    p, s = result.n_clusters, result.n_segments
    labels.append(f"{round((1 - s * p / len(raw)) * 100, 1)}% ({p}p x {s}s)")

    # Normalize at DataFrame level, then reshape
    reconstructed = result.reconstructed
    normalized = (reconstructed - raw_min) / raw_range
    data = normalized.values.reshape(n_days, period_duration, n_vars).transpose(2, 1, 0)

    frames_data.append(data.reshape(-1, n_days))

img_stack = np.stack(frames_data)

fig = px.imshow(
    img_stack,
    animation_frame=0,
    color_continuous_scale="RdYlBu_r",
    aspect="auto",
    labels={"x": "Day", "y": "Hour"},
    title="Time Series Aggregation",
)

for i, step in enumerate(fig.layout.sliders[0].steps):
    step["label"] = labels[i]

tickvals = [period_duration * i + period_duration // 2 for i in range(n_vars)]
fig.update_yaxes(tickvals=tickvals, ticktext=list(raw.columns))
fig.update_layout(height=600, coloraxis_showscale=False)
fig.show()

Save results

pareto_df.to_csv(RESULTS_DIR / "paretoOptimalAggregation.csv")
fig.write_html(RESULTS_DIR / "animation.html")