required_packages <- c( “fpp3”, “dplyr”, “lubridate”, “reticulate” )

missing_packages <- required_packages[ !required_packages %in% rownames(installed.packages())]

if (length(missing_packages) > 0) { install.packages(missing_packages) }

library(fpp3) library(dplyr) library(lubridate) library(reticulate)

Sys.setenv( RETICULATE_USE_MANAGED_VENV = “yes” )

py_require(c( “torch”, “lightning”, “pytorch-forecasting”, “pandas”, “numpy”, “matplotlib”, “scikit-learn” ))

set.seed(123)

cat(“configuration:”) print(py_config())

data( “vic_elec”, package = “tsibbledata” )

electricity <- vic_elec |> as_tibble() |> transmute( datetime = as.POSIXct(Time), demand = as.numeric(Demand), temperature = as.numeric(Temperature) ) |> arrange(datetime) |> # Use the final two years to reduce training time filter( datetime >= max(datetime) - years(2) ) |> mutate( series = “Victoria”,

# Lagged demand features
lag_1 = lag(demand, 1),
lag_48 = lag(demand, 48),
lag_336 = lag(demand, 48 * 7),

# Time-based features
hour = hour(datetime) +
  minute(datetime) / 60,

day_of_week = wday(
  datetime,
  week_start = 1
),

month = month(datetime),

weekend = as.integer(
  day_of_week >= 6
),

# Cyclical time features
hour_sin = sin(
  2 * pi * hour / 24
),

hour_cos = cos(
  2 * pi * hour / 24
),

dow_sin = sin(
  2 * pi * day_of_week / 7
),

dow_cos = cos(
  2 * pi * day_of_week / 7
),

month_sin = sin(
  2 * pi * month / 12
),

month_cos = cos(
  2 * pi * month / 12
)

) |> # Remove incomplete rows created by lagging filter( complete.cases( demand, temperature, lag_1, lag_48, lag_336 ) ) |> # Create a continuous integer index mutate( time_idx = row_number() - 1L )

cat(“dataset dimensions:”) print(dim(electricity))

number_of_rows <- nrow(electricity)

train_end_row <- floor( 0.70 * number_of_rows )

validation_end_row <- floor( 0.85 * number_of_rows )

train_time_end <- electricity$time_idx[ train_end_row]

validation_time_end <- electricity$time_idx[ validation_end_row]

electricity <- electricity |> mutate( split = case_when( time_idx <= train_time_end ~ “Training”,

  time_idx <= validation_time_end ~
    "Validation",

  TRUE ~
    "Testing"
)

)

cat(“split:”) print(table(electricity$split))

variables_to_scale <- c( “temperature”, “lag_1”, “lag_48”, “lag_336” )

training_rows <- electricity |> filter(split == “Training”)

for (variable in variables_to_scale) {

training_mean <- mean( training_rows[[variable]], na.rm = TRUE )

training_sd <- sd( training_rows[[variable]], na.rm = TRUE )

if ( is.na(training_sd) || training_sd == 0 ) { stop( paste( “Invalid standard deviation for”, variable ) ) }

scaled_name <- paste0( variable, “_scaled” )

electricity[[scaled_name]] <- ( electricity[[variable]] - training_mean ) / training_sd }

cat(“completed.”)

scaled_names <- paste0( variables_to_scale, “_scaled” )

scaled_matrix <- as.matrix( electricity[, scaled_names] )

if (any(!is.finite(scaled_matrix))) { stop( “The scaled data contain invalid values.” ) }

cat(“data passed validation.”)

python_main <- import_main( convert = FALSE )

python_main$electricity_data <- r_to_py( electricity, convert = FALSE )

python_main$training_cutoff_r <- r_to_py( as.integer(train_time_end), convert = FALSE )

python_main$validation_cutoff_r <- r_to_py( as.integer(validation_time_end), convert = FALSE )

cat(“transferred to Python.”)

py_run_string(r”—(

import os import random import numpy as np import pandas as pd import torch import matplotlib.pyplot as plt

import lightning.pytorch as pl

from lightning.pytorch.callbacks import EarlyStopping from lightning.pytorch.callbacks import ModelCheckpoint from lightning.pytorch.loggers import CSVLogger

from pytorch_forecasting import TimeSeriesDataSet from pytorch_forecasting import TemporalFusionTransformer

from pytorch_forecasting.data import GroupNormalizer

from pytorch_forecasting.metrics import QuantileLoss from pytorch_forecasting.metrics import MAE

data = electricity_data.copy()

training_cutoff = int( np.asarray( training_cutoff_r ).reshape(-1)[0] )

validation_cutoff = int( np.asarray( validation_cutoff_r ).reshape(-1)[0] )

data[“series”] = ( data[“series”] .astype(str) )

data[“weekend”] = ( data[“weekend”] .astype(str) )

data[“time_idx”] = ( data[“time_idx”] .astype(int) )

numeric_columns = [ “demand”, “temperature_scaled”, “lag_1_scaled”, “lag_48_scaled”, “lag_336_scaled”, “hour_sin”, “hour_cos”, “dow_sin”, “dow_cos”, “month_sin”, “month_cos”]

for column in numeric_columns: data[column] = pd.to_numeric( data[column], errors=“coerce” )

data = ( data .dropna( subset=numeric_columns ) .reset_index(drop=True) )

print(“——————————————”) print(“DATA RECEIVED BY PYTHON”) print(“——————————————”) print(“Rows:”, len(data)) print(“Training cutoff:”, training_cutoff) print(“Validation cutoff:”, validation_cutoff) print(“——————————————”)

random.seed(123) np.random.seed(123) torch.manual_seed(123)

pl.seed_everything( 123, workers=True )

max_encoder_length = 48 * 7 max_prediction_length = 48

training_data = data[ data[“time_idx”] <= training_cutoff].copy()

training_dataset = TimeSeriesDataSet( training_data,

time_idx="time_idx",
target="demand",
group_ids=["series"],

max_encoder_length=max_encoder_length,
max_prediction_length=max_prediction_length,

static_categoricals=[
    "series"
],

time_varying_known_categoricals=[
    "weekend"
],

# Features known for future dates
time_varying_known_reals=[
    "time_idx",
    "hour_sin",
    "hour_cos",
    "dow_sin",
    "dow_cos",
    "month_sin",
    "month_cos"
],

# Features observed through the encoder period
time_varying_unknown_reals=[
    "demand",
    "temperature_scaled",
    "lag_1_scaled",
    "lag_48_scaled",
    "lag_336_scaled"
],

target_normalizer=GroupNormalizer(
    groups=["series"],
    transformation="softplus"
),

add_relative_time_idx=True,
add_target_scales=True,
add_encoder_length=True

)

validation_data = data[ data[“time_idx”] <= validation_cutoff].copy()

validation_dataset = TimeSeriesDataSet.from_dataset( training_dataset, validation_data,

min_prediction_idx=(
    training_cutoff + 1
),

stop_randomization=True

)

testing_dataset = TimeSeriesDataSet.from_dataset( training_dataset, data,

min_prediction_idx=(
    validation_cutoff + 1
),

stop_randomization=True

)

train_loader = training_dataset.to_dataloader( train=True, batch_size=64, num_workers=0 )

validation_loader = validation_dataset.to_dataloader( train=False, batch_size=128, num_workers=0 )

test_loader = testing_dataset.to_dataloader( train=False, batch_size=128, num_workers=0 )

print(“——————————————”) print(“SUPERVISED LEARNING WINDOWS”) print(“——————————————”) print(“Training:”, len(training_dataset)) print(“Validation:”, len(validation_dataset)) print(“Testing:”, len(testing_dataset)) print(“——————————————”)

early_stopping = EarlyStopping( monitor=“val_loss”, min_delta=0.0001, patience=5, mode=“min” )

checkpoint_callback = ModelCheckpoint( monitor=“val_loss”, mode=“min”, save_top_k=1, filename=“best-tft” )

csv_logger = CSVLogger( save_dir=“tft_logs”, name=“vic_elec” )

tft = TemporalFusionTransformer.from_dataset( training_dataset,

learning_rate=0.001,
hidden_size=32,
attention_head_size=4,
hidden_continuous_size=16,
lstm_layers=2,
dropout=0.15,

loss=QuantileLoss(),

logging_metrics=torch.nn.ModuleList(
    [MAE()]
),

reduce_on_plateau_patience=3

)

print( “Model parameters:”, round( tft.size() / 1000, 1 ), “thousand” )

trainer = pl.Trainer( max_epochs=5, accelerator=“auto”, devices=1, gradient_clip_val=0.1,

callbacks=[
    early_stopping,
    checkpoint_callback
],

logger=csv_logger,
enable_model_summary=True,
log_every_n_steps=10

)

trainer.fit( tft, train_dataloaders=train_loader, val_dataloaders=validation_loader )

best_model_path = ( checkpoint_callback.best_model_path )

if not best_model_path: raise RuntimeError( “No model checkpoint was created.” )

best_tft = ( TemporalFusionTransformer .load_from_checkpoint( best_model_path ) )

print( “Best checkpoint:”, best_model_path )

metrics_path = os.path.join( csv_logger.log_dir, “metrics.csv” )

metrics = pd.read_csv( metrics_path )

if “train_loss_epoch” in metrics.columns:

training_loss = (
    metrics
    .dropna(
        subset=["train_loss_epoch"]
    )
    .groupby(
        "epoch",
        as_index=False
    )["train_loss_epoch"]
    .mean()
)

elif “train_loss_step” in metrics.columns:

training_loss = (
    metrics
    .dropna(
        subset=["train_loss_step"]
    )
    .groupby(
        "epoch",
        as_index=False
    )["train_loss_step"]
    .mean()
    .rename(
        columns={
            "train_loss_step":
                "train_loss_epoch"
        }
    )
)

else: raise RuntimeError( “Training loss was not found.” )

validation_loss = ( metrics .dropna( subset=[“val_loss”] ) .groupby( “epoch”, as_index=False )[“val_loss”] .mean() )

loss_history = pd.merge( training_loss, validation_loss, on=“epoch”, how=“outer” ).sort_values(“epoch”)

plt.figure( figsize=(9, 5) )

plt.plot( loss_history[“epoch”], loss_history[“train_loss_epoch”], marker=“o”, linewidth=2, label=“Training loss” )

plt.plot( loss_history[“epoch”], loss_history[“val_loss”], marker=“o”, linewidth=2, label=“Validation loss” )

plt.xlabel(“Epoch”) plt.ylabel(“Quantile loss”)

plt.title( “TFT Training and Validation Loss” )

plt.legend() plt.grid(alpha=0.25) plt.tight_layout()

plt.savefig( “training_validation_loss.png”, dpi=300, bbox_inches=“tight” )

plt.show() plt.close()

test_results = best_tft.predict( test_loader, mode=“prediction”, return_y=True,

trainer_kwargs={
    "accelerator": "auto",
    "devices": 1
}

)

predicted = ( test_results.output .detach() .cpu() .numpy() )

actual = ( test_results.y[0] .detach() .cpu() .numpy() )

predicted = np.squeeze(predicted) actual = np.squeeze(actual)

print( “Prediction dimensions:”, predicted.shape )

print( “Actual dimensions:”, actual.shape )

valid_mask = ( np.isfinite(actual) & np.isfinite(predicted) )

actual_valid = actual[ valid_mask]

predicted_valid = predicted[ valid_mask]

if actual_valid.size == 0: raise RuntimeError( “No valid test observations were available.” )

1. Root Mean Squared Error

test_rmse = np.sqrt( np.mean( ( actual_valid - predicted_valid ) ** 2 ) )

2. Mean Absolute Error

test_mae = np.mean( np.abs( actual_valid - predicted_valid ) )

3. Mean Absolute Percentage Error

nonzero_mask = ( np.abs(actual_valid) > np.finfo(float).eps )

if not np.any(nonzero_mask): raise RuntimeError( “MAPE cannot be calculated because all actual values are zero.” )

test_mape = np.mean( np.abs( ( actual_valid[nonzero_mask] - predicted_valid[nonzero_mask] ) / actual_valid[nonzero_mask] ) ) * 100

print(“——————————————”) print(“TEST-SET FORECAST ACCURACY”) print(“——————————————”)

print( “1. RMSE:”, round(test_rmse, 2), “MW” )

print( “2. MAE:”, round(test_mae, 2), “MW” )

print( “3. MAPE:”, round(test_mape, 2), “%” )

print(“——————————————”)

if predicted.ndim == 1: predicted = predicted.reshape( 1, -1 )

if actual.ndim == 1: actual = actual.reshape( 1, -1 )

number_of_windows = min( 7, predicted.shape[0], actual.shape[0] )

actual_sequence = actual[ :number_of_windows].reshape(-1)

predicted_sequence = predicted[ :number_of_windows].reshape(-1)

test_hours = ( np.arange( len(actual_sequence) ) / 2 )

plt.figure( figsize=(13, 5) )

plt.plot( test_hours, actual_sequence, color=“black”, linewidth=1.5, label=“Actual demand” )

plt.plot( test_hours, predicted_sequence, color=“#D95319”, linewidth=1.5, alpha=0.85, label=“TFT prediction” )

plt.xlabel(“Hours”)

plt.ylabel( “Electricity demand in MW” )

plt.title( “TFT Predictions Versus Actual Test Values” )

plt.legend() plt.grid(alpha=0.25) plt.tight_layout()

plt.savefig( “tft_test_predictions.png”, dpi=300, bbox_inches=“tight” )

plt.show() plt.close()

final_test_rmse = float( test_rmse )

final_test_mae = float( test_mae )

final_test_mape = float( test_mape )

)—“)

final_rmse <- as.numeric( reticulate::py_to_r( py$final_test_rmse ) )

final_mae <- as.numeric( reticulate::py_to_r( py$final_test_mae ) )

final_mape <- as.numeric( reticulate::py_to_r( py$final_test_mape ) )

if ( !is.finite(final_rmse) || !is.finite(final_mae) || !is.finite(final_mape) ) { stop( “One or more evaluation metrics are invalid.” ) }

cat(“”) cat(“——————————————”) cat(“FINAL MODEL RESULTS”) cat(“——————————————”)

cat( “1. RMSE:”, format( round(final_rmse, 2), nsmall = 2 ), “MW” )

cat( “2. MAE:”, format( round(final_mae, 2), nsmall = 2 ), “MW” )

cat( “3. MAPE:”, format( round(final_mape, 2), nsmall = 2 ), “%” )

cat(“——————————————”)

cat(“:”) cat( “RMSE penalizes larger forecast errors more heavily.” ) cat( “MAE is the average absolute forecast error in MW.” ) cat( “MAPE is the average absolute percentage error.” )

cat(“plots:”) cat(“1. training_validation_loss.png”) cat(“2. tft_test_predictions.png”)