# -*- coding: utf-8 -*- """ ============================================================================= ĐÁ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) # 1. Nạp dữ liệu 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" # Lọc mask đồng nhất 14,385 mẫu 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 # Chuẩn bị 5 seeds độc lập 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) # Tính toán metrics metrics = evaluate_metrics(y_eval, pred_eval, f"Residual KAN (Seed {s})") metrics["seed"] = s metrics["training_time_sec"] = round(dur, 2) # Đánh giá cải thiện so với Persistence Naive 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) # Thống kê tổng hợp Mean ± Std 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) # Lưu bảng CSV df_seeds.to_csv(reports_dir / "kan_multi_seed_stability.csv", index=False, encoding="utf-8-sig") # Tính Ensemble trung bình của 5 seeds 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) # Kiểm định thống kê Contiguous MBB trên mô hình Ensemble 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) # 4. Vẽ biểu đồ thẩm mỹ cao thể hiện độ ổn định qua các seeds fig, axes = plt.subplots(1, 2, figsize=(14, 5)) # Đồ thị MAE qua các Seed 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") # Đồ thị RMSE qua các Seed 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()