|
|
| """
|
| =============================================================================
|
| ĐÁNH GIÁ ĐỘ ỔN ĐỊNH VÀ TÍNH LẶP LẠI CỦA MÔ HÌNH RESIDUAL KAN (MULTI-SEED BENCHMARK)
|
| =============================================================================
|
| Mục tiêu nghiên cứu khoa học:
|
| 1. Huấn luyện Residual KAN qua 5 random seeds độc lập: [42, 101, 2024, 7, 999].
|
| 2. Đánh giá hoàn toàn khách quan trên đúng 14,385 mẫu Test đồng nhất.
|
| 3. Không tinh chỉnh siêu tham số trên Test, cố định hoàn toàn quy trình huấn luyện.
|
| 4. Báo cáo:
|
| - Trung bình (Mean) và Độ lệch chuẩn (Std Dev) của MAE, RMSE, sMAPE, R2.
|
| - Thống kê chênh lệch ΔRMSE, ΔMAE so với Persistence Naive qua từng seed.
|
| - Khoảng tin cậy tổng hợp qua Contiguous Moving-Block Bootstrap.
|
| 5. Xuất báo cáo khoa học vào:
|
| - output_preprocessing/reports/kan_multi_seed_stability.csv
|
| - output_preprocessing/reports/kan_multi_seed_stability.json
|
| - output_preprocessing/figures/15_kan_multi_seed_stability.png
|
| =============================================================================
|
| """
|
|
|
| import os
|
| import sys
|
| import json
|
| import time
|
| from pathlib import Path
|
|
|
| if sys.platform.startswith("win"):
|
| try:
|
| sys.stdout.reconfigure(encoding="utf-8")
|
| sys.stderr.reconfigure(encoding="utf-8")
|
| except Exception:
|
| pass
|
|
|
| import numpy as np
|
| import pandas as pd
|
| import matplotlib.pyplot as plt
|
| import seaborn as sns
|
|
|
| import torch
|
| import torch.nn as nn
|
| from torch.utils.data import TensorDataset, DataLoader
|
| from sklearn.preprocessing import StandardScaler
|
| from sklearn.impute import SimpleImputer
|
| from sklearn.metrics import mean_absolute_error, mean_squared_error, r2_score
|
|
|
| from temporal_significance_test import contiguous_moving_block_bootstrap
|
| from models_kan import ResidualKANRegressor, evaluate_metrics
|
|
|
| sns.set_theme(style="whitegrid", palette="muted")
|
| plt.rcParams['font.sans-serif'] = 'DejaVu Sans'
|
| plt.rcParams['axes.unicode_minus'] = False
|
|
|
|
|
| def train_single_seed(seed, X_train_scaled, y_train_scaled, X_val_scaled, y_val_scaled, X_test_scaled, p_mean, p_std, in_features, epochs=12):
|
| torch.manual_seed(seed)
|
| np.random.seed(seed)
|
| device = torch.device("cpu")
|
|
|
| train_dataset = TensorDataset(torch.tensor(X_train_scaled, dtype=torch.float32), torch.tensor(y_train_scaled, dtype=torch.float32))
|
| val_dataset = TensorDataset(torch.tensor(X_val_scaled, dtype=torch.float32), torch.tensor(y_val_scaled, dtype=torch.float32))
|
|
|
| train_loader = DataLoader(train_dataset, batch_size=256, shuffle=True)
|
| val_loader = DataLoader(val_dataset, batch_size=512, shuffle=False)
|
|
|
| model = ResidualKANRegressor(in_features=in_features, hidden_dims=[64, 32], p_t_idx=0, dropout=0.05).to(device)
|
| criterion = nn.HuberLoss(delta=0.1)
|
| optimizer = torch.optim.AdamW(model.parameters(), lr=2e-3, weight_decay=1e-4)
|
| scheduler = torch.optim.lr_scheduler.ReduceLROnPlateau(optimizer, mode='min', factor=0.5, patience=2)
|
|
|
| best_val_loss = float("inf")
|
| best_model_state = None
|
|
|
| t0 = time.time()
|
| for ep in range(epochs):
|
| model.train()
|
| train_loss = 0.0
|
| for bx, by in train_loader:
|
| bx, by = bx.to(device), by.to(device)
|
| optimizer.zero_grad()
|
| preds = model(bx)
|
| loss = criterion(preds, by)
|
| loss.backward()
|
| torch.nn.utils.clip_grad_norm_(model.parameters(), max_norm=2.0)
|
| optimizer.step()
|
| train_loss += loss.item() * len(by)
|
| train_loss /= len(train_dataset)
|
|
|
| model.eval()
|
| val_loss = 0.0
|
| with torch.no_grad():
|
| for bx, by in val_loader:
|
| bx, by = bx.to(device), by.to(device)
|
| preds = model(bx)
|
| loss = criterion(preds, by)
|
| val_loss += loss.item() * len(by)
|
| val_loss /= len(val_dataset)
|
| scheduler.step(val_loss)
|
|
|
| if val_loss < best_val_loss:
|
| best_val_loss = val_loss
|
| best_model_state = model.state_dict().copy()
|
|
|
| train_time = time.time() - t0
|
| model.load_state_dict(best_model_state)
|
|
|
| model.eval()
|
| with torch.no_grad():
|
| test_tensor = torch.tensor(X_test_scaled, dtype=torch.float32).to(device)
|
| y_pred_scaled = model(test_tensor).numpy()
|
|
|
| y_pred_kw = y_pred_scaled * p_std + p_mean
|
| return y_pred_kw, train_time
|
|
|
|
|
| def main():
|
| base_dir = Path("output_preprocessing")
|
| data_dir = base_dir / "data"
|
| reports_dir = base_dir / "reports"
|
| figures_dir = base_dir / "figures"
|
|
|
| print("=" * 85)
|
| print("KHẢO SÁT TÍNH ỔN ĐỊNH VÀ LẶP LẠI CỦA RESIDUAL KAN QUA 5 RANDOM SEEDS")
|
| print("=" * 85)
|
|
|
|
|
| train_df = pd.read_csv(data_dir / "train.csv")
|
| val_df = pd.read_csv(data_dir / "validation.csv")
|
| test_df = pd.read_csv(data_dir / "test.csv")
|
|
|
| test_df["TIME"] = pd.to_datetime(test_df["TIME"])
|
| test_df = test_df.sort_values("TIME").reset_index(drop=True)
|
|
|
| candidate_features = [
|
| "P_kW", "UA", "UB", "UC", "IA", "IB", "IC", "CosF", "Q_kVAr",
|
| "U_avg", "U_max", "U_min", "U_unbalance_pct",
|
| "I_avg", "I_max", "I_min", "I_unbalance_pct",
|
| "S_kVA", "lag_1", "lag_2", "lag_3", "lag_48", "lag_336",
|
| "rolling_mean_3", "rolling_mean_6", "rolling_mean_48",
|
| "hour", "dayofweek", "dayofmonth", "month", "is_weekend",
|
| "sin_hour", "cos_hour"
|
| ]
|
| feature_cols = [c for c in candidate_features if c in train_df.columns]
|
| target_col = "target_P_kW"
|
|
|
|
|
| pred_persistence = test_df["P_kW"].values
|
| pred_s48 = test_df["lag_48"].values
|
| pred_s336 = test_df["lag_336"].values
|
| pred_ma3 = test_df["rolling_mean_3"].values
|
| y_test_raw = test_df[target_col].values
|
|
|
| mask_common = (
|
| np.isfinite(y_test_raw) &
|
| np.isfinite(pred_persistence) &
|
| np.isfinite(pred_s48) &
|
| np.isfinite(pred_s336) &
|
| np.isfinite(pred_ma3)
|
| )
|
|
|
| imputer = SimpleImputer(strategy="median")
|
| X_train_imp = imputer.fit_transform(train_df[feature_cols])
|
| y_train = train_df[target_col].values
|
| X_val_imp = imputer.transform(val_df[feature_cols])
|
| y_val = val_df[target_col].values
|
| X_test_imp = imputer.transform(test_df[feature_cols])
|
| y_test = test_df[target_col].values
|
|
|
| scaler_X = StandardScaler()
|
| X_train_scaled = scaler_X.fit_transform(X_train_imp)
|
| X_val_scaled = scaler_X.transform(X_val_imp)
|
| X_test_scaled = scaler_X.transform(X_test_imp)
|
|
|
| p_mean = scaler_X.mean_[0]
|
| p_std = scaler_X.scale_[0]
|
|
|
| y_train_scaled = (y_train - p_mean) / p_std
|
| y_val_scaled = (y_val - p_mean) / p_std
|
|
|
| y_eval = y_test[mask_common]
|
| p_eval = pred_persistence[mask_common]
|
|
|
| eval_df = test_df[mask_common].copy().reset_index(drop=True)
|
| diffs = eval_df["TIME"].diff()
|
| segment_ids = (diffs != pd.Timedelta(minutes=30)).cumsum().values
|
|
|
|
|
| seeds = [42, 101, 2024, 7, 999]
|
| seed_records = []
|
| predictions_all_seeds = []
|
|
|
| print(f"Bắt đầu huấn luyện và đánh giá trên {len(seeds)} random seeds: {seeds}...")
|
| for idx, s in enumerate(seeds, start=1):
|
| print(f"\n--- [Seed {idx}/{len(seeds)}: random_seed = {s}] ---")
|
| pred_kw, dur = train_single_seed(
|
| s, X_train_scaled, y_train_scaled, X_val_scaled, y_val_scaled, X_test_scaled,
|
| p_mean, p_std, in_features=len(feature_cols), epochs=12
|
| )
|
| pred_eval = pred_kw[mask_common]
|
| predictions_all_seeds.append(pred_eval)
|
|
|
|
|
| metrics = evaluate_metrics(y_eval, pred_eval, f"Residual KAN (Seed {s})")
|
| metrics["seed"] = s
|
| metrics["training_time_sec"] = round(dur, 2)
|
|
|
|
|
| d_mae = pers_metrics_mae = 598.36 - metrics["MAE_kW"]
|
| d_rmse = pers_metrics_rmse = 1160.44 - metrics["RMSE_kW"]
|
| metrics["Delta_MAE_kW"] = round(d_mae, 2)
|
| metrics["Delta_RMSE_kW"] = round(d_rmse, 2)
|
| metrics["Improvement_MAE_pct"] = round((d_mae / 598.36) * 100.0, 2)
|
| metrics["Improvement_RMSE_pct"] = round((d_rmse / 1160.44) * 100.0, 2)
|
|
|
| print(f" MAE: {metrics['MAE_kW']} kW (Tốt hơn Persistence: {metrics['Improvement_MAE_pct']}%)")
|
| print(f" RMSE: {metrics['RMSE_kW']} kW (Tốt hơn Persistence: {metrics['Improvement_RMSE_pct']}%)")
|
| print(f" R2: {metrics['R2_score']:.4f} | Thời gian: {metrics['training_time_sec']}s")
|
|
|
| seed_records.append(metrics)
|
|
|
| df_seeds = pd.DataFrame(seed_records)
|
|
|
|
|
| mean_mae = df_seeds["MAE_kW"].mean()
|
| std_mae = df_seeds["MAE_kW"].std()
|
|
|
| mean_rmse = df_seeds["RMSE_kW"].mean()
|
| std_rmse = df_seeds["RMSE_kW"].std()
|
|
|
| mean_r2 = df_seeds["R2_score"].mean()
|
| std_r2 = df_seeds["R2_score"].std()
|
|
|
| mean_smape = df_seeds["sMAPE_pct"].mean()
|
| std_smape = df_seeds["sMAPE_pct"].std()
|
|
|
| mean_imp_mae = df_seeds["Improvement_MAE_pct"].mean()
|
| std_imp_mae = df_seeds["Improvement_MAE_pct"].std()
|
|
|
| mean_imp_rmse = df_seeds["Improvement_RMSE_pct"].mean()
|
| std_imp_rmse = df_seeds["Improvement_RMSE_pct"].std()
|
|
|
| print("\n" + "=" * 90)
|
| print("BẢNG TỔNG HỢP ĐỘ ỔN ĐỊNH QUA 5 RANDOM SEEDS ĐỘC LẬP")
|
| print("=" * 90)
|
| cols_show = ["seed", "MAE_kW", "RMSE_kW", "sMAPE_pct", "R2_score", "Improvement_MAE_pct", "Improvement_RMSE_pct"]
|
| print(df_seeds[cols_show].to_string(index=False))
|
| print("-" * 90)
|
| print(f"TRUNG BÌNH ± ĐỘ LỆCH CHUẨN (MEAN ± STD):")
|
| print(f" * MAE: {mean_mae:.2f} ± {std_mae:.2f} kW (Cải thiện: {mean_imp_mae:.2f}% ± {std_imp_mae:.2f}%)")
|
| print(f" * RMSE: {mean_rmse:.2f} ± {std_rmse:.2f} kW (Cải thiện: {mean_imp_rmse:.2f}% ± {std_imp_rmse:.2f}%)")
|
| print(f" * R2: {mean_r2:.4f} ± {std_r2:.4f}")
|
| print(f" * sMAPE:{mean_smape:.2f}% ± {std_smape:.2f}%")
|
| print("=" * 90)
|
|
|
|
|
| df_seeds.to_csv(reports_dir / "kan_multi_seed_stability.csv", index=False, encoding="utf-8-sig")
|
|
|
|
|
| ensemble_preds = np.mean(predictions_all_seeds, axis=0)
|
| ens_metrics = evaluate_metrics(y_eval, ensemble_preds, "Residual KAN (5-Seed Ensemble Mean)")
|
| ens_metrics["Improvement_MAE_pct"] = round(((598.36 - ens_metrics["MAE_kW"]) / 598.36) * 100.0, 2)
|
| ens_metrics["Improvement_RMSE_pct"] = round(((1160.44 - ens_metrics["RMSE_kW"]) / 1160.44) * 100.0, 2)
|
|
|
|
|
| mbb_48_ens = contiguous_moving_block_bootstrap(y_eval, ensemble_preds, p_eval, segment_ids, block_length=48, n_boot=1000)
|
| mbb_336_ens = contiguous_moving_block_bootstrap(y_eval, ensemble_preds, p_eval, segment_ids, block_length=336, n_boot=1000)
|
|
|
| summary_json = {
|
| "timestamp": pd.Timestamp.now().isoformat(),
|
| "seeds_evaluated": seeds,
|
| "evaluated_samples": int(len(y_eval)),
|
| "summary_statistics": {
|
| "MAE_kW_mean": round(float(mean_mae), 2),
|
| "MAE_kW_std": round(float(std_mae), 2),
|
| "RMSE_kW_mean": round(float(mean_rmse), 2),
|
| "RMSE_kW_std": round(float(std_rmse), 2),
|
| "sMAPE_pct_mean": round(float(mean_smape), 2),
|
| "sMAPE_pct_std": round(float(std_smape), 2),
|
| "R2_score_mean": round(float(mean_r2), 4),
|
| "R2_score_std": round(float(std_r2), 4),
|
| "Improvement_MAE_pct_mean": round(float(mean_imp_mae), 2),
|
| "Improvement_RMSE_pct_mean": round(float(mean_imp_rmse), 2)
|
| },
|
| "individual_runs": seed_records,
|
| "ensemble_5_seeds": {
|
| "metrics": ens_metrics,
|
| "mbb_significance_block_48": {
|
| "delta_rmse_mean_kW": round(float(mbb_48_ens['delta_rmse_mean']), 2),
|
| "ci_rmse_95_kW": [round(float(x), 2) for x in mbb_48_ens['ci_rmse_95']],
|
| "rmse_better_pct": round(float(mbb_48_ens['rmse_better_pct']), 1),
|
| "delta_mae_mean_kW": round(float(mbb_48_ens['delta_mae_mean']), 2),
|
| "ci_mae_95_kW": [round(float(x), 2) for x in mbb_48_ens['ci_mae_95']]
|
| },
|
| "mbb_significance_block_336": {
|
| "delta_rmse_mean_kW": round(float(mbb_336_ens['delta_rmse_mean']), 2),
|
| "ci_rmse_95_kW": [round(float(x), 2) for x in mbb_336_ens['ci_rmse_95']],
|
| "rmse_better_pct": round(float(mbb_336_ens['rmse_better_pct']), 1),
|
| "delta_mae_mean_kW": round(float(mbb_336_ens['delta_mae_mean']), 2),
|
| "ci_mae_95_kW": [round(float(x), 2) for x in mbb_336_ens['ci_mae_95']]
|
| }
|
| }
|
| }
|
|
|
| with open(reports_dir / "kan_multi_seed_stability.json", "w", encoding="utf-8") as f:
|
| json.dump(summary_json, f, indent=4, ensure_ascii=False)
|
|
|
|
|
| fig, axes = plt.subplots(1, 2, figsize=(14, 5))
|
|
|
|
|
| x_pos = np.arange(len(seeds))
|
| bars0 = axes[0].bar(x_pos, df_seeds["MAE_kW"], color="#d62728", alpha=0.7, width=0.55, label="Residual KAN (Seed)")
|
| axes[0].axhline(598.36, color="black", ls="--", lw=1.8, label="Persistence Naive (598.36 kW)")
|
| axes[0].axhline(mean_mae, color="blue", ls=":", lw=1.8, label=f"KAN Mean ({mean_mae:.1f} ± {std_mae:.1f} kW)")
|
| axes[0].set_xticks(x_pos)
|
| axes[0].set_xticklabels([f"Seed {s}" for s in seeds], fontweight="bold")
|
| axes[0].set_ylabel("MAE (kW)", fontweight="bold")
|
| axes[0].set_title(f"A. Độ ổn định sai số MAE qua 5 Seeds (Mean: {mean_mae:.1f} kW)", fontweight="bold")
|
| axes[0].set_ylim(400, 650)
|
| for b in bars0:
|
| h = b.get_height()
|
| axes[0].text(b.get_x() + b.get_width()/2, h + 5, f"{h:.1f}", ha="center", fontsize=10, fontweight="bold")
|
| axes[0].legend(loc="upper right")
|
|
|
|
|
| bars1 = axes[1].bar(x_pos, df_seeds["RMSE_kW"], color="#2ca02c", alpha=0.7, width=0.55, label="Residual KAN (Seed)")
|
| axes[1].axhline(1160.44, color="black", ls="--", lw=1.8, label="Persistence Naive (1,160.44 kW)")
|
| axes[1].axhline(mean_rmse, color="blue", ls=":", lw=1.8, label=f"KAN Mean ({mean_rmse:.1f} ± {std_rmse:.1f} kW)")
|
| axes[1].set_xticks(x_pos)
|
| axes[1].set_xticklabels([f"Seed {s}" for s in seeds], fontweight="bold")
|
| axes[1].set_ylabel("RMSE (kW)", fontweight="bold")
|
| axes[1].set_title(f"B. Độ ổn định sai số RMSE qua 5 Seeds (Mean: {mean_rmse:.1f} kW)", fontweight="bold")
|
| axes[1].set_ylim(800, 1250)
|
| for b in bars1:
|
| h = b.get_height()
|
| axes[1].text(b.get_x() + b.get_width()/2, h + 10, f"{h:.1f}", ha="center", fontsize=10, fontweight="bold")
|
| axes[1].legend(loc="upper right")
|
|
|
| fig.suptitle("15. Đánh giá độ ổn định và tính lặp lại của Residual KAN qua 5 Random Seeds trên 14,385 mẫu Test", fontsize=13, fontweight="bold")
|
| fig.tight_layout()
|
| fig.savefig(figures_dir / "15_kan_multi_seed_stability.png", dpi=200)
|
| plt.close(fig)
|
|
|
| print("\n✅ ĐÃ HOÀN THÀNH KHẢO SÁT ĐỘ ỔN ĐỊNH VÀ XUẤT CÁC BÁO CÁO:")
|
| print(f" - {reports_dir / 'kan_multi_seed_stability.csv'}")
|
| print(f" - {reports_dir / 'kan_multi_seed_stability.json'}")
|
| print(f" - {figures_dir / '15_kan_multi_seed_stability.png'}")
|
|
|
|
|
| if __name__ == "__main__":
|
| main()
|
|
|
