|
|
| """
|
| =============================================================================
|
| MÔ HÌNH KOLMOGOROV-ARNOLD NETWORK (KAN) CHO DỰ BÁO PHỤ TẢI ĐIỆN 30 PHÚT
|
| BẢN NÂNG CẤP: CHUẨN HÓA MỤC TIÊU (TARGET SCALING) + RESIDUAL LEARNING
|
| =============================================================================
|
| Nguyên lý kỹ thuật:
|
| 1. Target Scaling: target_P_kW có giá trị trung bình 8,240 kW. Nếu không chuẩn hóa y,
|
| mạng nơ-ron sâu với hàm kích hoạt spline/silu biên độ [-1, 1] phải học hệ số khuếch
|
| đại khổng lồ qua hàng trăm epochs mới tới được dải 8,000.
|
| 2. Residual Learning: Mạng KAN học sai số delta: y_pred = P_t + KAN(X)
|
| (Dự báo giá trị gia số tải thay vì dự báo toàn phần).
|
| 3. Đánh giá trên CHÍNH XÁC 14,385 mẫu Test đồng nhất.
|
| 4. Kiểm định Contiguous Moving-Block Bootstrap (Block 48 & Block 336).
|
| =============================================================================
|
| """
|
|
|
| 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
|
| import torch.nn.functional as F
|
| 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
|
|
|
| sns.set_theme(style="whitegrid", palette="muted")
|
| plt.rcParams['font.sans-serif'] = 'DejaVu Sans'
|
| plt.rcParams['axes.unicode_minus'] = False
|
|
|
|
|
| class KANLinear(nn.Module):
|
| def __init__(self, in_features, out_features, grid_size=5, spline_order=3, scale_base=1.0, scale_spline=1.0):
|
| super().__init__()
|
| self.in_features = in_features
|
| self.out_features = out_features
|
| self.grid_size = grid_size
|
| self.spline_order = spline_order
|
|
|
| h = 2.0 / grid_size
|
| grid = torch.linspace(-1.0 - spline_order * h, 1.0 + spline_order * h, grid_size + 2 * spline_order + 1)
|
| self.register_buffer("grid", grid)
|
|
|
| self.base_weight = nn.Parameter(torch.empty(out_features, in_features))
|
| self.spline_weight = nn.Parameter(torch.empty(out_features, in_features * (grid_size + spline_order)))
|
|
|
| self.scale_base = scale_base
|
| self.scale_spline = scale_spline
|
| self.reset_parameters()
|
|
|
| def reset_parameters(self):
|
| nn.init.kaiming_uniform_(self.base_weight, a=np.sqrt(5))
|
| nn.init.normal_(self.spline_weight, mean=0.0, std=0.02)
|
|
|
| def b_splines(self, x):
|
| grid = self.grid
|
| x = x.unsqueeze(-1)
|
| bases = ((x >= grid[:-1]) & (x < grid[1:])).float()
|
| for k in range(1, self.spline_order + 1):
|
| bases = (
|
| (x - grid[: -(k + 1)]) / (grid[k:-1] - grid[: -(k + 1)]) * bases[:, :, :-1]
|
| ) + (
|
| (grid[k + 1 :] - x) / (grid[k + 1 :] - grid[1:(-k)]) * bases[:, :, 1:]
|
| )
|
| return bases.contiguous()
|
|
|
| def forward(self, x):
|
| base_output = F.linear(F.silu(x), self.base_weight)
|
| x_clamped = torch.tanh(x)
|
| spline_basis = self.b_splines(x_clamped).view(x.size(0), -1)
|
| spline_output = F.linear(spline_basis, self.spline_weight)
|
| return self.scale_base * base_output + self.scale_spline * spline_output
|
|
|
|
|
| class ResidualKANRegressor(nn.Module):
|
| """
|
| Mạng KAN học phần bù thặng dư (Residual KAN):
|
| y_hat_scaled = P_t_scaled + KAN(X)
|
| """
|
| def __init__(self, in_features, hidden_dims=[64, 32], p_t_idx=0, dropout=0.05):
|
| super().__init__()
|
| self.p_t_idx = p_t_idx
|
| layers = []
|
| dims = [in_features] + hidden_dims + [1]
|
| for i in range(len(dims) - 1):
|
| layers.append(KANLinear(dims[i], dims[i+1], grid_size=5, spline_order=3))
|
| if i < len(dims) - 2:
|
| layers.append(nn.LayerNorm(dims[i+1]))
|
| layers.append(nn.Dropout(dropout))
|
| self.kan_body = nn.Sequential(*layers)
|
|
|
| def forward(self, x):
|
|
|
| delta = self.kan_body(x).squeeze(-1)
|
| base = x[:, self.p_t_idx]
|
| return base + delta
|
|
|
|
|
| def evaluate_metrics(y_true, y_pred, model_name="Model"):
|
| mae = mean_absolute_error(y_true, y_pred)
|
| rmse = np.sqrt(mean_squared_error(y_true, y_pred))
|
| smape = np.mean(2.0 * np.abs(y_pred - y_true) / (np.abs(y_true) + np.abs(y_pred) + 1e-5)) * 100.0
|
| r2 = r2_score(y_true, y_pred)
|
| return {
|
| "model_name": model_name,
|
| "MAE_kW": round(float(mae), 2),
|
| "RMSE_kW": round(float(rmse), 2),
|
| "MAE_MW": round(float(mae / 1000.0), 4),
|
| "RMSE_MW": round(float(rmse / 1000.0), 4),
|
| "sMAPE_pct": round(float(smape), 2),
|
| "R2_score": round(float(r2), 4)
|
| }
|
|
|
|
|
| def main():
|
| base_dir = Path("output_preprocessing")
|
| data_dir = base_dir / "data"
|
| reports_dir = base_dir / "reports"
|
| figures_dir = base_dir / "figures"
|
|
|
| reports_dir.mkdir(parents=True, exist_ok=True)
|
| figures_dir.mkdir(parents=True, exist_ok=True)
|
|
|
| print("=" * 80)
|
| print("HUẤN LUYỆN RESIDUAL KOLMOGOROV-ARNOLD NETWORK (RESIDUAL KAN)")
|
| print("=" * 80)
|
|
|
|
|
| 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
|
|
|
| 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)
|
|
|
|
|
| torch.manual_seed(42)
|
| np.random.seed(42)
|
| device = torch.device("cpu")
|
|
|
| model = ResidualKANRegressor(in_features=len(feature_cols), 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)
|
|
|
| epochs = 12
|
| best_val_loss = float("inf")
|
| best_model_state = None
|
|
|
| print(f"[3/5] Bắt đầu huấn luyện Residual KAN ({epochs} epochs)...")
|
| 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()
|
|
|
| print(f" Epoch {ep+1:02d}/{epochs} | Scaled Train Loss: {train_loss:.5f} | Val Loss: {val_loss:.5f}")
|
|
|
| train_time = time.time() - t0
|
| model.load_state_dict(best_model_state)
|
|
|
|
|
| print("\n[4/5] Đánh giá trên 14,385 mẫu Test đồng nhất...")
|
| 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_all = y_pred_scaled * p_std + p_mean
|
|
|
| y_eval = y_test[mask_common]
|
| p_eval = pred_persistence[mask_common]
|
| kan_eval = y_pred_all[mask_common]
|
|
|
| kan_metrics = evaluate_metrics(y_eval, kan_eval, "Residual Kolmogorov-Arnold Network (KAN)")
|
| kan_metrics["training_time_sec"] = round(train_time, 2)
|
| kan_metrics["samples_evaluated"] = int(len(y_eval))
|
|
|
| pers_metrics = evaluate_metrics(y_eval, p_eval, "Persistence Naive (P_t)")
|
|
|
| print("\n" + "=" * 90)
|
| print("BẢNG KẾT QUẢ ĐỐI CHIẾU TRÊN 14,385 MẪU TEST:")
|
| print("=" * 90)
|
| df_res = pd.DataFrame([pers_metrics, kan_metrics])
|
| print(df_res[["model_name", "MAE_kW", "RMSE_kW", "MAE_MW", "RMSE_MW", "sMAPE_pct", "R2_score"]].to_string(index=False))
|
| print("=" * 90)
|
|
|
|
|
| print("\n[5/5] Kiểm định thống kê Contiguous Moving-Block Bootstrap...")
|
| 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
|
|
|
| mbb_48 = contiguous_moving_block_bootstrap(y_eval, kan_eval, p_eval, segment_ids, block_length=48, n_boot=1000)
|
| mbb_336 = contiguous_moving_block_bootstrap(y_eval, kan_eval, p_eval, segment_ids, block_length=336, n_boot=1000)
|
|
|
| print(f"\nResidual KAN vs Persistence [Block 48 (24h)]:")
|
| print(f" ΔRMSE trung bình: {mbb_48['delta_rmse_mean']:+.2f} kW | 95% CI: [{mbb_48['ci_rmse_95'][0]:.2f}; {mbb_48['ci_rmse_95'][1]:.2f}] kW")
|
| print(f" Tỷ lệ thắng RMSE: {mbb_48['rmse_better_pct']:.1f}%")
|
| print(f" ΔMAE trung bình: {mbb_48['delta_mae_mean']:+.2f} kW | 95% CI: [{mbb_48['ci_mae_95'][0]:.2f}; {mbb_48['ci_mae_95'][1]:.2f}] kW")
|
|
|
| print(f"\nResidual KAN vs Persistence [Block 336 (7 ngày)]:")
|
| print(f" ΔRMSE trung bình: {mbb_336['delta_rmse_mean']:+.2f} kW | 95% CI: [{mbb_336['ci_rmse_95'][0]:.2f}; {mbb_336['ci_rmse_95'][1]:.2f}] kW")
|
| print(f" Tỷ lệ thắng RMSE: {mbb_336['rmse_better_pct']:.1f}%")
|
| print(f" ΔMAE trung bình: {mbb_336['delta_mae_mean']:+.2f} kW | 95% CI: [{mbb_336['ci_mae_95'][0]:.2f}; {mbb_336['ci_mae_95'][1]:.2f}] kW")
|
|
|
|
|
| kan_summary = {
|
| "timestamp": pd.Timestamp.now().isoformat(),
|
| "architecture": "Residual KAN (y_hat = P_t + KAN_delta) with B-Splines (order=3, grid=5)",
|
| "hidden_dims": [64, 32],
|
| "input_features": len(feature_cols),
|
| "evaluated_samples": len(y_eval),
|
| "metrics": kan_metrics,
|
| "mbb_significance_test": {
|
| "block_48_steps_24h": {
|
| "delta_rmse_mean_kW": round(float(mbb_48['delta_rmse_mean']), 2),
|
| "ci_rmse_95_kW": [round(float(x), 2) for x in mbb_48['ci_rmse_95']],
|
| "rmse_better_pct": round(float(mbb_48['rmse_better_pct']), 1),
|
| "delta_mae_mean_kW": round(float(mbb_48['delta_mae_mean']), 2),
|
| "ci_mae_95_kW": [round(float(x), 2) for x in mbb_48['ci_mae_95']],
|
| "is_significant_at_005": bool(mbb_48['ci_rmse_95'][0] > 0)
|
| },
|
| "block_336_steps_7d": {
|
| "delta_rmse_mean_kW": round(float(mbb_336['delta_rmse_mean']), 2),
|
| "ci_rmse_95_kW": [round(float(x), 2) for x in mbb_336['ci_rmse_95']],
|
| "rmse_better_pct": round(float(mbb_336['rmse_better_pct']), 1),
|
| "delta_mae_mean_kW": round(float(mbb_336['delta_mae_mean']), 2),
|
| "ci_mae_95_kW": [round(float(x), 2) for x in mbb_336['ci_mae_95']],
|
| "is_significant_at_005": bool(mbb_336['ci_rmse_95'][0] > 0)
|
| }
|
| }
|
| }
|
|
|
| with open(reports_dir / "kan_model_results.json", "w", encoding="utf-8") as f:
|
| json.dump(kan_summary, f, indent=4, ensure_ascii=False)
|
|
|
|
|
| fig, axes = plt.subplots(2, 1, figsize=(16, 10))
|
| zoom_slice = slice(-336, None)
|
| time_zoom = test_df["TIME"].iloc[zoom_slice]
|
|
|
| axes[0].plot(time_zoom, y_test[zoom_slice], label="Thực tế (Ground Truth)", color="black", lw=2.0)
|
| axes[0].plot(time_zoom, pred_persistence[zoom_slice], label="Persistence Naive", color="#7f7f7f", ls="--", lw=1.2, alpha=0.7)
|
| axes[0].plot(time_zoom, y_pred_all[zoom_slice], label=f"Residual KAN (MAE: {kan_metrics['MAE_kW']} kW | RMSE: {kan_metrics['RMSE_kW']} kW)", color="#d62728", lw=1.6, alpha=0.9)
|
| axes[0].set_title("14A. Đối chiếu kết quả dự báo phụ tải của Mạng Residual KAN trên tập Test độc lập (Khung 7 ngày)", fontsize=13, fontweight="bold")
|
| axes[0].set_xlabel("Thời gian")
|
| axes[0].set_ylabel("Công suất P (kW)")
|
| axes[0].legend(loc="upper right", frameon=True)
|
|
|
| residuals_kan = y_eval - kan_eval
|
| sns.histplot(residuals_kan, bins=60, kde=True, ax=axes[1], color="#d62728", alpha=0.6)
|
| axes[1].axvline(0, color="black", ls="--", lw=1.5, label="Sai số = 0 kW")
|
| axes[1].axvline(np.mean(residuals_kan), color="blue", ls=":", lw=1.5, label=f"Sai số trung bình: {np.mean(residuals_kan):.2f} kW")
|
| axes[1].set_title(f"14B. Phân phối sai số phần dư của Mạng Residual KAN (MAE: {kan_metrics['MAE_kW']} kW | R2: {kan_metrics['R2_score']:.4f})", fontsize=13, fontweight="bold")
|
| axes[1].set_xlabel("Sai số phần dư = P_true - P_pred (kW)")
|
| axes[1].legend(loc="upper right", frameon=True)
|
|
|
| fig.tight_layout()
|
| fig.savefig(figures_dir / "14_kan_predictions_and_residuals.png", dpi=200)
|
| plt.close(fig)
|
|
|
| print("\n✅ Đã lưu kết quả Residual KAN vào:")
|
| print(f" - {reports_dir / 'kan_model_results.json'}")
|
| print(f" - {figures_dir / '14_kan_predictions_and_residuals.png'}")
|
|
|
|
|
| if __name__ == "__main__":
|
| main()
|
|
|
