AutoControlCourse/case_demo_functions.py

import os
import sys
# ===== OpenMP冲突兼容设置，避免PyTorch初始化报错 =====
os.environ.setdefault("KMP_DUPLICATE_LIB_OK", "TRUE")
import numpy as np
import matplotlib
matplotlib.use('Agg')
import matplotlib.pyplot as plt

# English-only plot style (no Chinese fonts needed)
plt.rcParams['font.family'] = 'serif'
plt.rcParams['font.serif'] = ['DejaVu Serif', 'Times New Roman']
plt.rcParams['axes.unicode_minus'] = True

MODEL_SRC_PATH = os.path.join(os.path.dirname(os.path.abspath(__file__)), "Model", "src")
if MODEL_SRC_PATH not in sys.path:
    sys.path.insert(0, MODEL_SRC_PATH)

MODEL_DATA_PATH = os.path.join(os.path.dirname(os.path.abspath(__file__)), "Model", "data")


# ============================================================
# 阶段零：模型蒸馏演示 (GPR/CSV → NN)
# ============================================================
def run_distillation_demo(epochs, learning_rate, hidden_size, progress=None):
    """
    运行蒸馏并返回训练曲线 + 对比散点图。
    """
    try:
        from distill_gpr_to_nn import distill_from_csv

        csv_path = os.path.join(MODEL_DATA_PATH, "Cleaned_Engine_Data_Full.csv")
        nn_path = os.path.join(MODEL_DATA_PATH, "engine_nn_proxy.pth")

        epochs = int(np.clip(epochs, 500, 8000))
        learning_rate = float(np.clip(learning_rate, 1e-4, 1e-2))
        hidden_size = int(np.clip(hidden_size, 16, 256))

        if progress is not None:
            progress(0.0, desc="加载数据...")

        def _progress_cb(epoch, total, loss):
            if progress is not None:
                progress(epoch / total, desc=f"训练中 Epoch {epoch}/{total}, Loss={loss:.6f}")

        result = distill_from_csv(
            csv_path, nn_path,
            epochs=epochs, lr=learning_rate, hidden_size=hidden_size,
            verbose=True, progress_callback=_progress_cb
        )

        loss_hist = result['loss_history']
        Y_true = result['Y_train']
        Y_pred = result['Y_nn']

        # ---- Plots ----
        fig, axes = plt.subplots(1, 3, figsize=(16, 5))
        fig.suptitle('GPR / Data  →  NN Knowledge Distillation Results',
                     fontweight='bold', fontsize=13)

        # Loss curve
        axes[0].semilogy(loss_hist, 'b-', lw=1.2)
        axes[0].set_xlabel('Epoch')
        axes[0].set_ylabel('MSE Loss (log scale)')
        axes[0].set_title('Training Loss Curve')
        axes[0].grid(True, linestyle=':', alpha=0.7)

        # Fuel flow parity
        axes[1].scatter(Y_true[:, 0], Y_pred[:, 0], s=8, alpha=0.5, c='tab:orange')
        lim = [0, max(Y_true[:, 0].max(), Y_pred[:, 0].max()) * 1.05]
        axes[1].plot(lim, lim, 'k--', lw=1, alpha=0.7)
        axes[1].set_xlabel('True Fuel Flow (kg/h)')
        axes[1].set_ylabel('NN Predicted Fuel Flow (kg/h)')
        axes[1].set_title(f'Fuel Flow Parity (MAPE={result["rel_error_fuel"]:.2f}%)')
        axes[1].set_xlim(lim); axes[1].set_ylim(lim)
        axes[1].set_aspect('equal')
        axes[1].grid(True, linestyle=':', alpha=0.7)

        # Power parity
        axes[2].scatter(Y_true[:, 1], Y_pred[:, 1], s=8, alpha=0.5, c='tab:blue')
        lim = [0, max(Y_true[:, 1].max(), Y_pred[:, 1].max()) * 1.05]
        axes[2].plot(lim, lim, 'k--', lw=1, alpha=0.7)
        axes[2].set_xlabel('True Power (kW)')
        axes[2].set_ylabel('NN Predicted Power (kW)')
        axes[2].set_title(f'Power Parity (MAPE={result["rel_error_power"]:.2f}%)')
        axes[2].set_xlim(lim); axes[2].set_ylim(lim)
        axes[2].set_aspect('equal')
        axes[2].grid(True, linestyle=':', alpha=0.7)

        fig.tight_layout(rect=[0, 0, 1, 0.94])

        n_params = sum(p.numel() for p in result['nn_model'].parameters())
        summary = (
            f"### Distillation Results\n"
            f"- **NN Architecture**: MLP 3→{hidden_size}→{hidden_size}→2 (Tanh)\n"
            f"- **Parameters**: {n_params:,}\n"
            f"- **Training Samples**: {len(Y_true)}\n"
            f"- **Epochs**: {epochs}, LR: {learning_rate:.1e}\n"
            f"- **Final Loss**: {loss_hist[-1]:.6f}\n"
            f"- **Fuel Flow MAPE**: {result['rel_error_fuel']:.2f}%\n"
            f"- **Power MAPE**: {result['rel_error_power']:.2f}%\n"
            f"- **Model saved** to `Model/data/engine_nn_proxy.pth`"
        )

        return fig, summary

    except Exception as e:
        import traceback
        return None, f"Distillation failed: {e}\n```\n{traceback.format_exc()}\n```"


# ============================================================
# GPR 模型训练/加载 与 可视化
# ============================================================
def run_gpr_training(mode="load", progress=None):
    """
    GPR 模型训练或加载已有模型，并生成可视化图表。

    Parameters
    ----------
    mode : str
        "train" — 从头训练（需 botorch/gpytorch/sklearn）
        "load"  — 加载已有的 .pth 权重文件
    progress : gr.Progress or None
    """
    import csv as csv_mod
    try:
        csv_path = os.path.join(MODEL_DATA_PATH, "Cleaned_Engine_Data_Full.csv")
        gpr_pth  = os.path.join(MODEL_DATA_PATH, "engine_gpr_model.pth")

        # ---------- 读取 CSV 原始数据（不依赖 pandas）----------
        if progress is not None:
            progress(0.05, desc="读取 CSV 数据...")
        with open(csv_path, 'r', encoding='utf-8') as f:
            reader = csv_mod.reader(f)
            header = next(reader)
            rows = [r for r in reader]
        col_idx = {name: i for i, name in enumerate(header)}
        data = np.array([[float(x) for x in r] for r in rows], dtype=np.float64)
        X_cols = ['Altitude_m', 'Mach', 'RPM']
        Y_cols = ['WF_kg_h', 'Power_kW']
        X_raw = data[:, [col_idx[c] for c in X_cols]]
        Y_raw = data[:, [col_idx[c] for c in Y_cols]]
        n_total = len(data)

        # ---------- 尝试导入 GPR 依赖 ----------
        gpr_available = False
        gpr_model = None
        try:
            from engine_gpr_class import EngineGPRModel
            gpr_available = True
        except ImportError:
            gpr_available = False

        if mode == "train":
            if not gpr_available:
                return None, ("### ⚠️ GPR 训练失败\n\n"
                              "缺少依赖包：`botorch`, `gpytorch`, `sklearn`。\n\n"
                              "请执行 `pip install botorch gpytorch scikit-learn` 后重试，"
                              "或选择 **加载已有模型** 模式。")
            if progress is not None:
                progress(0.10, desc="初始化 GPR 模型 ...")
            gpr_model = EngineGPRModel(csv_path=csv_path)
            if progress is not None:
                progress(0.15, desc="训练 GPR（超参数优化中）...")
            gpr_model.train(save_path=gpr_pth)
            if progress is not None:
                progress(0.80, desc="GPR 训练完成，生成可视化 ...")

        elif mode == "load":
            if not gpr_available:
                # --- 无 botorch：仅展示原始数据统计 ---
                if progress is not None:
                    progress(0.30, desc="绘制数据统计图 ...")
                fig = _plot_data_overview(X_raw, Y_raw, X_cols, Y_cols)
                summary = (
                    f"### 📊 数据概览（无 GPR 依赖）\n"
                    f"- **数据集**: Cleaned_Engine_Data_Full.csv\n"
                    f"- **样本数**: {n_total}\n"
                    f"- **输入特征**: {', '.join(X_cols)}\n"
                    f"- **输出目标**: {', '.join(Y_cols)}\n\n"
                    f"> ⚠️ 未安装 `botorch`/`gpytorch`，无法加载 GPR 模型。\n"
                    f"> 请执行 `pip install botorch gpytorch scikit-learn` 后重试。"
                )
                return fig, summary

            if not os.path.exists(gpr_pth):
                return None, ("### ⚠️ 未找到已训练的 GPR 权重文件\n\n"
                              f"路径：`{gpr_pth}`\n\n"
                              "请先选择 **从头训练** 模式。")
            if progress is not None:
                progress(0.10, desc="加载 GPR 模型 ...")
            gpr_model = EngineGPRModel(csv_path=csv_path)
            ok = gpr_model.load_model(pth_path=gpr_pth)
            if not ok:
                return None, "### ⚠️ GPR 模型加载失败，请检查权重文件完整性。"
            if progress is not None:
                progress(0.40, desc="加载完成，生成可视化 ...")

        # ---------- GPR 模型已就绪，生成可视化 ----------
        if progress is not None:
            progress(0.50, desc="GPR 网格预测 ...")

        # 预测网格
        H_range  = np.linspace(X_raw[:, 0].min(), X_raw[:, 0].max(), 40)
        Ma_range = np.linspace(X_raw[:, 1].min(), X_raw[:, 1].max(), 5)
        RPM_range = np.linspace(max(X_raw[:, 2].min(), 1000), X_raw[:, 2].max(), 40)
        H, Ma, RPM = np.meshgrid(H_range, Ma_range, RPM_range, indexing='ij')
        X_grid = np.column_stack([H.ravel(), Ma.ravel(), RPM.ravel()])

        pred_mean, pred_var = gpr_model.predict(X_grid)

        if progress is not None:
            progress(0.75, desc="绘图中 ...")

        # 在训练数据点上的预测精度（过滤极小值，与 NN 训练一致）
        valid_mask = (Y_raw[:, 0] > 0.5) & (Y_raw[:, 1] > 0.5) & (X_raw[:, 2] > 500)
        X_eval = X_raw[valid_mask]
        Y_eval = Y_raw[valid_mask]
        train_pred, _ = gpr_model.predict(X_eval)
        mape_fuel = np.mean(np.abs(train_pred[:, 0] - Y_eval[:, 0]) / np.maximum(Y_eval[:, 0], 1e-6)) * 100
        mape_power = np.mean(np.abs(train_pred[:, 1] - Y_eval[:, 1]) / np.maximum(Y_eval[:, 1], 1e-6)) * 100

        # ---- Plot ----
        fig, axes = plt.subplots(2, 2, figsize=(14, 10))
        fig.suptitle('GPR Model Training / Validation Results', fontweight='bold', fontsize=13)

        # Fuel flow parity
        axes[0, 0].scatter(Y_eval[:, 0], train_pred[:, 0], s=10, alpha=0.5, c='tab:orange')
        lim = [0, max(Y_eval[:, 0].max(), train_pred[:, 0].max()) * 1.05]
        axes[0, 0].plot(lim, lim, 'k--', lw=1)
        axes[0, 0].set_xlabel('True Fuel Flow (kg/h)')
        axes[0, 0].set_ylabel('GPR Predicted Fuel Flow (kg/h)')
        axes[0, 0].set_title(f'Fuel Flow Parity (MAPE={mape_fuel:.2f}%)')
        axes[0, 0].set_xlim(lim); axes[0, 0].set_ylim(lim)
        axes[0, 0].set_aspect('equal'); axes[0, 0].grid(True, ls=':', alpha=0.7)

        # Power parity
        axes[0, 1].scatter(Y_eval[:, 1], train_pred[:, 1], s=10, alpha=0.5, c='tab:blue')
        lim = [0, max(Y_eval[:, 1].max(), train_pred[:, 1].max()) * 1.05]
        axes[0, 1].plot(lim, lim, 'k--', lw=1)
        axes[0, 1].set_xlabel('True Power (kW)')
        axes[0, 1].set_ylabel('GPR Predicted Power (kW)')
        axes[0, 1].set_title(f'Power Parity (MAPE={mape_power:.2f}%)')
        axes[0, 1].set_xlim(lim); axes[0, 1].set_ylim(lim)
        axes[0, 1].set_aspect('equal'); axes[0, 1].grid(True, ls=':', alpha=0.7)

        # Variance heatmap (Mach=0 slice)
        ma0_idx = np.argmin(np.abs(Ma_range - 0.0))
        var_reshaped = pred_var.reshape(len(H_range), len(Ma_range), len(RPM_range), 2)
        wf_var_slice = var_reshaped[:, ma0_idx, :, 0]
        pow_var_slice = var_reshaped[:, ma0_idx, :, 1]

        RPM_g, H_g = np.meshgrid(RPM_range, H_range)
        cf0 = axes[1, 0].contourf(RPM_g, H_g, np.log10(np.maximum(wf_var_slice, 1e-16)),
                                    levels=30, cmap='jet', alpha=0.85)
        axes[1, 0].set_xlabel('RPM')
        axes[1, 0].set_ylabel('Altitude (m)')
        axes[1, 0].set_title('Fuel Flow Variance (log₁₀, Mach=0)')
        fig.colorbar(cf0, ax=axes[1, 0], shrink=0.8)

        cf1 = axes[1, 1].contourf(RPM_g, H_g, np.log10(np.maximum(pow_var_slice, 1e-16)),
                                    levels=30, cmap='jet', alpha=0.85)
        axes[1, 1].set_xlabel('RPM')
        axes[1, 1].set_ylabel('Altitude (m)')
        axes[1, 1].set_title('Power Variance (log₁₀, Mach=0)')
        fig.colorbar(cf1, ax=axes[1, 1], shrink=0.8)

        fig.tight_layout(rect=[0, 0, 1, 0.94])

        mode_label = "从头训练" if mode == "train" else "加载已有模型"
        import torch as _torch
        device_info = "CUDA" if _torch.cuda.is_available() else "CPU"
        summary = (
            f"### GPR 模型结果\n"
            f"- **模式**: {mode_label}\n"
            f"- **计算设备**: {device_info}\n"
            f"- **训练样本**: {n_total}（有效评估样本: {int(valid_mask.sum())}）\n"
            f"- **网格预测点**: {len(X_grid)}\n"
            f"- **Fuel Flow MAPE**: {mape_fuel:.2f}%\n"
            f"- **Power MAPE**: {mape_power:.2f}%\n"
            f"- **模型文件**: `Model/data/engine_gpr_model.pth`"
        )

        if progress is not None:
            progress(1.0, desc="完成")
        return fig, summary

    except Exception as e:
        import traceback
        return None, f"GPR 训练/加载失败: {e}\n```\n{traceback.format_exc()}\n```"


def _plot_data_overview(X_raw, Y_raw, X_cols, Y_cols):
    """当 GPR 依赖不可用时，仅绘制原始数据统计概览。"""
    fig, axes = plt.subplots(2, 2, figsize=(14, 10))
    fig.suptitle('Engine Data Overview (GPR dependencies unavailable)', fontweight='bold', fontsize=13)

    # Altitude vs Fuel Flow
    axes[0, 0].scatter(X_raw[:, 0], Y_raw[:, 0], s=6, alpha=0.4, c='tab:orange')
    axes[0, 0].set_xlabel('Altitude (m)')
    axes[0, 0].set_ylabel('Fuel Flow (kg/h)')
    axes[0, 0].set_title('Altitude vs Fuel Flow')
    axes[0, 0].grid(True, ls=':', alpha=0.7)

    # RPM vs Power
    axes[0, 1].scatter(X_raw[:, 2], Y_raw[:, 1], s=6, alpha=0.4, c='tab:blue')
    axes[0, 1].set_xlabel('RPM')
    axes[0, 1].set_ylabel('Power (kW)')
    axes[0, 1].set_title('RPM vs Power')
    axes[0, 1].grid(True, ls=':', alpha=0.7)

    # RPM vs Fuel Flow colored by Mach
    sc = axes[1, 0].scatter(X_raw[:, 2], Y_raw[:, 0], s=6, alpha=0.4, c=X_raw[:, 1], cmap='viridis')
    axes[1, 0].set_xlabel('RPM')
    axes[1, 0].set_ylabel('Fuel Flow (kg/h)')
    axes[1, 0].set_title('RPM vs Fuel Flow (color=Mach)')
    fig.colorbar(sc, ax=axes[1, 0], shrink=0.8, label='Mach')
    axes[1, 0].grid(True, ls=':', alpha=0.7)

    # Fuel Flow vs Power
    axes[1, 1].scatter(Y_raw[:, 0], Y_raw[:, 1], s=6, alpha=0.4, c='tab:green')
    axes[1, 1].set_xlabel('Fuel Flow (kg/h)')
    axes[1, 1].set_ylabel('Power (kW)')
    axes[1, 1].set_title('Fuel Flow vs Power')
    axes[1, 1].grid(True, ls=':', alpha=0.7)

    fig.tight_layout(rect=[0, 0, 1, 0.94])
    return fig


# ============================================================
# 阶段一：发动机控制器设计 (PID / MPC 可选)
# ============================================================
def run_engine_design(sim_time_s, dt, initial_power_kw, target_power_kw,
                      controller_type,
                      kp, ki, kd, tau_fuel, K_inertia,
                      mpc_horizon, mpc_W_power, mpc_W_dcost, mpc_overshoot_limit,
                      progress=None):
    """发动机控制器阶跃响应仿真"""
    try:
        import torch
        from lightweight_model import EngineNNProxy

        # 参数裁剪
        sim_time_s = float(np.clip(sim_time_s, 5, 120))
        dt = float(np.clip(dt, 0.01, 0.2))
        initial_power_kw = float(np.clip(initial_power_kw, 20, 260))
        target_power_kw = float(np.clip(target_power_kw, 20, 300))
        tau_fuel = float(np.clip(tau_fuel, 0.02, 2.0))
        K_inertia = float(np.clip(K_inertia, 5, 1000))

        nn_pth = os.path.join(MODEL_DATA_PATH, "engine_nn_proxy.pth")
        if not os.path.exists(nn_pth):
            return None, "Error: `engine_nn_proxy.pth` not found. Please run the **Distillation** tab first."

        engine_nn = EngineNNProxy()
        engine_nn.load_state_dict(torch.load(nn_pth, map_location='cpu'))
        engine_nn.eval()

        # 初始稳态转速 — 纯Python二分法（不依赖scipy）
        def _bisect(func, a, b, tol=1e-4, maxiter=50):
            fa, fb = func(a), func(b)
            if fa * fb > 0:
                return a if abs(fa) < abs(fb) else b
            for _ in range(maxiter):
                c = (a + b) / 2.0
                fc = func(c)
                if abs(fc) < tol or (b - a) / 2 < tol:
                    return c
                if fa * fc < 0:
                    b, fb = c, fc
                else:
                    a, fa = c, fc
            return (a + b) / 2.0

        def _solve_rpm(target_p):
            def obj(n):
                inp = torch.tensor([[0.0, 0.0, n]], dtype=torch.float32)
                with torch.no_grad():
                    return engine_nn(inp).numpy()[0, 1] - target_p
            return _bisect(obj, 1000, 58000)

        N_current = _solve_rpm(initial_power_kw)
        with torch.no_grad():
            pred0 = engine_nn(torch.tensor([[0.0, 0.0, N_current]], dtype=torch.float32)).numpy()
        Wf_act = pred0[0, 0]
        Wf_cmd = Wf_act

        use_mpc = (controller_type == "MPC")

        if use_mpc:
            from mpc_controller import TurboShaftMPCController
            mpc_horizon = int(np.clip(mpc_horizon, 3, 30))
            mpc_W_power = float(np.clip(mpc_W_power, 1, 1000))
            mpc_W_dcost = float(np.clip(mpc_W_dcost, 0.01, 50))
            mpc_overshoot_limit = float(np.clip(mpc_overshoot_limit, 0.01, 0.30))
            mpc = TurboShaftMPCController(
                tau_fuel=tau_fuel, K_inertia=K_inertia, dt=dt,
                horizon=mpc_horizon, min_fuel=5.0, max_fuel=400.0,
                overshoot_limit=mpc_overshoot_limit
            )
            mpc.W_power = mpc_W_power
            mpc.W_dcost = mpc_W_dcost
            mpc.reset(initial_output=Wf_cmd, initial_N=N_current)
        else:
            from increPID import IncrementalPIDController
            kp = float(np.clip(kp, 0.01, 30))
            ki = float(np.clip(ki, 0.0, 30))
            kd = float(np.clip(kd, 0.0, 10))
            pid = IncrementalPIDController(
                kp=kp, ki=ki, kd=kd, dt=dt,
                output_min=5.0, output_max=400.0,
                input_scale=300.0, output_scale=400.0
            )
            pid.reset(initial_output=Wf_act)

        time_array = np.arange(0, sim_time_s, dt)
        t_step = sim_time_s * 0.15

        logs = {'N': [], 'Wf_act': [], 'Wf_cmd': [], 'Power': [], 'Power_target': []}

        n_steps = len(time_array)
        for step_i, t in enumerate(time_array):
            if progress is not None and step_i % max(1, n_steps // 20) == 0:
                progress(step_i / n_steps, desc=f"发动机仿真 {step_i}/{n_steps} (t={t:.1f}s)")
            target_p = initial_power_kw if t < t_step else target_power_kw

            delta_N = 5.0
            batch_inp = torch.tensor([
                [0.0, 0.0, N_current],
                [0.0, 0.0, N_current + delta_N]
            ], dtype=torch.float32)
            with torch.no_grad():
                batch_pred = engine_nn(batch_inp).numpy()

            Wf_req = batch_pred[0, 0]
            Power = batch_pred[0, 1]

            if use_mpc:
                k_wf = (batch_pred[1, 0] - batch_pred[0, 0]) / delta_N
                k_p = (batch_pred[1, 1] - batch_pred[0, 1]) / delta_N
                Wf_cmd = mpc.compute(
                    current_N=N_current, current_Wfact=Wf_act,
                    target_power=target_p,
                    precalc_params=(Wf_req, Power, k_wf, k_p)
                )
            else:
                Wf_cmd = pid.compute(setpoint=target_p, measurement=Power)

            dWf = (Wf_cmd - Wf_act) / tau_fuel
            Wf_act += dWf * dt
            dN = K_inertia * (Wf_act - Wf_req)
            N_current += dN * dt

            logs['N'].append(N_current)
            logs['Wf_act'].append(Wf_act)
            logs['Wf_cmd'].append(Wf_cmd)
            logs['Power'].append(Power)
            logs['Power_target'].append(target_p)

        if progress is not None:
            progress(1.0, desc="绘图中...")
        # ---- Performance metrics ----
        power_arr = np.array(logs['Power'])
        step_idx = int(t_step / dt)
        post_step = power_arr[step_idx:]

        tail = max(1, len(post_step) // 10)
        ss_error = np.mean(np.abs(post_step[-tail:] - target_power_kw))

        delta = target_power_kw - initial_power_kw
        overshoot = 0.0
        if abs(delta) > 1:
            if delta > 0:
                overshoot = max(0, (np.max(post_step) - target_power_kw) / delta * 100)
            else:
                overshoot = max(0, (target_power_kw - np.min(post_step)) / abs(delta) * 100)

        rise_time = float('nan')
        if abs(delta) > 1:
            thresh_10 = initial_power_kw + 0.1 * delta
            thresh_90 = initial_power_kw + 0.9 * delta
            t10 = t90 = None
            for i in range(step_idx, len(power_arr)):
                if delta > 0:
                    if t10 is None and power_arr[i] >= thresh_10:
                        t10 = time_array[i] - time_array[step_idx]
                    if t90 is None and power_arr[i] >= thresh_90:
                        t90 = time_array[i] - time_array[step_idx]
                else:
                    if t10 is None and power_arr[i] <= thresh_10:
                        t10 = time_array[i] - time_array[step_idx]
                    if t90 is None and power_arr[i] <= thresh_90:
                        t90 = time_array[i] - time_array[step_idx]
            if t10 is not None and t90 is not None:
                rise_time = t90 - t10

        settling_time = float('nan')
        if abs(delta) > 1:
            band = abs(delta) * 0.02
            for i in range(len(post_step) - 1, -1, -1):
                if abs(post_step[i] - target_power_kw) > band:
                    settling_time = (i + 1) * dt
                    break

        # ---- Plots (English) ----
        ctrl_label = "MPC" if use_mpc else "PID"
        fig, axes = plt.subplots(3, 1, figsize=(12, 9), sharex=True)
        fig.suptitle(f'Engine Controller Step Response ({ctrl_label})',
                     fontweight='bold', fontsize=13)

        axes[0].plot(time_array, logs['Power_target'], 'k--', lw=1.5, label='Target Power')
        axes[0].plot(time_array, logs['Power'], 'r-', lw=1.5, label='Actual Power')
        axes[0].set_ylabel('Power (kW)')
        axes[0].set_title('Power Tracking')
        axes[0].grid(True, linestyle=':'); axes[0].legend()

        axes[1].plot(time_array, logs['N'], 'b-', lw=1.5)
        axes[1].set_ylabel('Speed (RPM)')
        axes[1].set_title('Engine Rotor Speed')
        axes[1].grid(True, linestyle=':')

        axes[2].plot(time_array, logs['Wf_cmd'], 'k--', lw=1.2, label='Fuel Command')
        axes[2].plot(time_array, logs['Wf_act'], 'r-', lw=1.2, label='Actual Fuel')
        axes[2].set_ylabel('Fuel Flow (kg/h)')
        axes[2].set_xlabel('Time (s)')
        axes[2].set_title('Fuel Control Signal')
        axes[2].grid(True, linestyle=':'); axes[2].legend()

        fig.tight_layout(rect=[0, 0, 1, 0.96])

        if use_mpc:
            param_str = (f"Horizon={mpc_horizon}, W_power={mpc_W_power:.1f}, "
                         f"W_Δcost={mpc_W_dcost:.2f}, Overshoot≤{mpc_overshoot_limit*100:.0f}%")
        else:
            param_str = f"Kp={kp:.3f}, Ki={ki:.3f}, Kd={kd:.3f}"

        summary = (
            f"### Engine Controller Results ({ctrl_label})\n"
            f"- **Controller**: {ctrl_label} — {param_str}\n"
            f"- **Power Step**: {initial_power_kw:.0f} → {target_power_kw:.0f} kW\n"
            f"- **Steady-State Error**: {ss_error:.2f} kW\n"
            f"- **Overshoot**: {overshoot:.1f}%\n"
            f"- **Rise Time (10%-90%)**: {rise_time:.3f} s\n"
            f"- **Settling Time (2% band)**: {settling_time:.3f} s\n"
            f"- **Fuel Actuator τ**: {tau_fuel:.2f} s | Rotor Inertia K: {K_inertia:.0f}"
        )
        return fig, summary

    except Exception as e:
        import traceback
        return None, f"Engine simulation failed: {e}\n```\n{traceback.format_exc()}\n```"


# ============================================================
# 阶段二：电机控制器设计 (PID / MPC 可选)
# ============================================================
def run_motor_design(sim_time_s, dt, target_rpm, load_torque,
                     controller_type,
                     kp, ki, kd, J,
                     mpc_W_speed, mpc_W_dcost, mpc_overshoot_limit,
                     progress=None):
    """电机控制器阶跃响应 + 负载扰动仿真"""
    try:
        from motor_sim import MotorSim

        sim_time_s = float(np.clip(sim_time_s, 5, 120))
        dt = float(np.clip(dt, 0.01, 0.2))
        target_rpm = float(np.clip(target_rpm, 200, 6000))
        load_torque = float(np.clip(load_torque, 5, 500))
        J = float(np.clip(J, 0.1, 10.0))

        use_mpc = (controller_type == "MPC")

        if use_mpc:
            mpc_W_speed = float(np.clip(mpc_W_speed, 1, 500))
            mpc_W_dcost = float(np.clip(mpc_W_dcost, 0.01, 50))
            mpc_overshoot_limit = float(np.clip(mpc_overshoot_limit, 0.01, 0.30))
            motor = MotorSim(
                P_rate=300e3, w_rate=575.95, J=J,
                mpc_W_speed=mpc_W_speed, mpc_W_dcost=mpc_W_dcost,
                mpc_overshoot_limit=mpc_overshoot_limit,
            )
        else:
            kp = float(np.clip(kp, 0.01, 80))
            ki = float(np.clip(ki, 0.0, 150))
            kd = float(np.clip(kd, 0.0, 10))
            motor = MotorSim(
                P_rate=300e3, w_rate=575.95, J=J,
                mpc_W_speed=0.0, mpc_W_dcost=0.0,
                mpc_overshoot_limit=0.05,
            )

        time_array = np.arange(0, sim_time_s, dt)
        t_step = sim_time_s * 0.10
        t_load_step = sim_time_s * 0.60
        v_bus = 520.0
        p_supply = 0.0

        if not use_mpc:
            from increPID import IncrementalPIDController
            w_rate = 575.95
            tau_rate = 300e3 / w_rate
            pid_motor = IncrementalPIDController(
                kp=kp, ki=ki, kd=kd, dt=dt,
                output_min=-tau_rate, output_max=tau_rate,
                input_scale=w_rate, output_scale=tau_rate
            )
            pid_motor.reset(0.0)

        logs = {'rpm': [], 'target_rpm': [], 'torque': [],
                'p_bus_req': [], 'p_shaft': [], 'p_loss': [], 'load': []}

        n_steps = len(time_array)
        for step_i, t in enumerate(time_array):
            if progress is not None and step_i % max(1, n_steps // 20) == 0:
                progress(step_i / n_steps, desc=f"电机仿真 {step_i}/{n_steps} (t={t:.1f}s)")
            n_set = 500.0 if t < t_step else target_rpm
            if t < t_step:
                load_t = 20.0
            elif t < t_load_step:
                load_t = load_torque
            else:
                load_t = load_torque * 1.5

            if not use_mpc:
                w_set_rad = n_set * 2 * np.pi / 60.0
                T_cmd = pid_motor.compute(setpoint=w_set_rad, measurement=motor.w_M)
                w_eff = max(abs(motor.w_M), 1.0)
                p_cmd_kw = -T_cmd * w_eff / 1000.0
                state = motor.step(dt=dt, n_setpoint=n_set, p_bus_actual_kw=p_cmd_kw,
                                   v_bus=v_bus, t_load=load_t, t_ext=0.0)
            else:
                state = motor.step(dt=dt, n_setpoint=n_set, p_bus_actual_kw=p_supply,
                                   v_bus=v_bus, t_load=load_t, t_ext=0.0)

            logs['rpm'].append(state['n_rpm'])
            logs['target_rpm'].append(n_set)
            logs['torque'].append(state['t_motor'])
            logs['p_bus_req'].append(state['p_bus_req_kw'])
            logs['p_shaft'].append(state['p_shaft_kw'])
            logs['p_loss'].append(state['p_loss_kw'])
            logs['load'].append(load_t)
            p_supply = state['p_bus_req_kw']

        if progress is not None:
            progress(1.0, desc="绘图中...")
        # ---- Performance ----
        rpm_arr = np.array(logs['rpm'])
        step_idx = int(t_step / dt)
        post_step_rpm = rpm_arr[step_idx:]

        tail = max(1, len(post_step_rpm) // 10)
        ss_error = np.mean(np.abs(post_step_rpm[-tail:] - target_rpm))

        delta_rpm = target_rpm - 500.0
        overshoot = 0.0
        if abs(delta_rpm) > 1 and delta_rpm > 0:
            overshoot = max(0, (np.max(post_step_rpm) - target_rpm) / delta_rpm * 100)

        load_step_idx = int(t_load_step / dt)
        max_dip = 0
        if load_step_idx < len(rpm_arr):
            post_load = rpm_arr[load_step_idx:]
            max_dip = max(0, target_rpm - np.min(post_load)) if len(post_load) > 0 else 0

        # ---- Plots (English) ----
        ctrl_label = "MPC" if use_mpc else "PID"
        fig, axes = plt.subplots(3, 1, figsize=(12, 9), sharex=True)
        fig.suptitle(f'Motor Controller — Step + Load Disturbance ({ctrl_label})',
                     fontweight='bold', fontsize=13)

        axes[0].plot(time_array, logs['target_rpm'], 'k--', lw=1.5, label='Target Speed')
        axes[0].plot(time_array, logs['rpm'], 'b-', lw=1.5, label='Actual Speed')
        axes[0].axvline(t_load_step, color='orange', linestyle=':', lw=1, alpha=0.7, label='Load Disturbance')
        axes[0].set_ylabel('Speed (RPM)')
        axes[0].set_title('Speed Tracking')
        axes[0].grid(True, linestyle=':'); axes[0].legend()

        axes[1].plot(time_array, logs['torque'], 'r-', lw=1.2, label='Motor Torque')
        axes[1].plot(time_array, logs['load'], 'k--', lw=1, alpha=0.6, label='Load Torque')
        axes[1].set_ylabel('Torque (Nm)')
        axes[1].set_title('Torque Response')
        axes[1].grid(True, linestyle=':'); axes[1].legend()

        axes[2].plot(time_array, logs['p_bus_req'], 'g-', lw=1.2, label='Bus Power Request')
        axes[2].plot(time_array, logs['p_shaft'], 'b--', lw=1.2, label='Shaft Power')
        axes[2].plot(time_array, logs['p_loss'], 'r:', lw=1.2, label='Loss Power')
        axes[2].set_ylabel('Power (kW)')
        axes[2].set_xlabel('Time (s)')
        axes[2].set_title('Power Distribution')
        axes[2].grid(True, linestyle=':'); axes[2].legend()

        fig.tight_layout(rect=[0, 0, 1, 0.96])

        if use_mpc:
            param_str = (f"W_speed={mpc_W_speed:.1f}, W_Δcost={mpc_W_dcost:.2f}, "
                         f"Overshoot≤{mpc_overshoot_limit*100:.0f}%")
        else:
            param_str = f"Kp={kp:.3f}, Ki={ki:.3f}, Kd={kd:.3f}"

        summary = (
            f"### Motor Controller Results ({ctrl_label})\n"
            f"- **Controller**: {ctrl_label} — {param_str}\n"
            f"- **Speed Step**: 500 → {target_rpm:.0f} RPM\n"
            f"- **Load Torque**: {load_torque:.0f} Nm → {load_torque*1.5:.0f} Nm\n"
            f"- **Steady-State Error**: {ss_error:.1f} RPM\n"
            f"- **Overshoot**: {overshoot:.1f}%\n"
            f"- **Max Load Dip**: {max_dip:.1f} RPM\n"
            f"- **Inertia J**: {J:.2f} kg·m²"
        )
        return fig, summary

    except Exception as e:
        import traceback
        return None, f"Motor simulation failed: {e}\n```\n{traceback.format_exc()}\n```"


# ============================================================
# 工况配置辅助
# ============================================================
def _profile_points(profile_name):
    if profile_name == "高机动阶跃":
        return [(0., 1600., 70.), (8., 3200., 220.), (20., 2500., 130.),
                (35., 3400., 250.), (50., 1800., 80.)]
    if profile_name == "经济巡航":
        return [(0., 1500., 60.), (15., 2100., 95.), (35., 2300., 105.),
                (55., 2000., 90.)]
    return [(0., 1500., 50.), (10., 3000., 200.), (30., 2800., 150.),
            (50., 1800., 60.)]


def _target_from_profile(t, points):
    rpm, torque = points[0][1], points[0][2]
    for p_t, p_rpm, p_torque in points:
        if t >= p_t:
            rpm, torque = p_rpm, p_torque
        else:
            break
    return rpm, torque


# ============================================================
# 阶段三：能量管理策略设计
# ============================================================
def run_hybrid_demo(sim_time_s, dt, initial_soc_pct, initial_engine_power_kw,
                    profile_name,
                    eng_controller_type, eng_kp, eng_ki, eng_kd,
                    eng_mpc_horizon, eng_mpc_W_power, eng_mpc_W_dcost, eng_mpc_overshoot,
                    mot_controller_type, mot_kp, mot_ki, mot_kd, mot_J,
                    mot_mpc_W_speed, mot_mpc_W_dcost, mot_mpc_overshoot,
                    soc_target_pct, soc_low_pct, soc_high_pct,
                    p_eng_min, p_eng_max, p_charge, k_soc,
                    power_reserve_pct, battery_capacity_kwh,
                    progress=None):
    """混动系统能量管理策略仿真 (规则 + 滞环)"""
    try:
        import torch
        from lightweight_model import EngineNNProxy
        from motor_sim import MotorSim
        from battery_sim import BatterySim

        # 参数裁剪
        sim_time_s = float(np.clip(sim_time_s, 10, 240))
        dt = float(np.clip(dt, 0.01, 0.2))
        initial_soc_pct = float(np.clip(initial_soc_pct, 10, 95))
        initial_engine_power_kw = float(np.clip(initial_engine_power_kw, 20, 260))
        mot_J = float(np.clip(mot_J, 0.1, 10.0))

        soc_target = float(np.clip(soc_target_pct, 20, 80)) / 100.0
        soc_low = float(np.clip(soc_low_pct, 10, 60)) / 100.0
        soc_high = float(np.clip(soc_high_pct, 50, 95)) / 100.0
        if soc_low >= soc_high:
            soc_high = soc_low + 0.1
        p_eng_min = float(np.clip(p_eng_min, 10, 100))
        p_eng_max = float(np.clip(p_eng_max, 100, 350))
        if p_eng_min >= p_eng_max:
            p_eng_min = p_eng_max * 0.1
        p_charge = float(np.clip(p_charge, 50, 300))
        k_soc = float(np.clip(k_soc, 0, 500))
        power_reserve_pct = float(np.clip(power_reserve_pct, 0, 50))
        battery_capacity_kwh = float(np.clip(battery_capacity_kwh, 10, 200))

        # ---- 发动机 ----
        nn_pth = os.path.join(MODEL_DATA_PATH, "engine_nn_proxy.pth")
        if not os.path.exists(nn_pth):
            return None, "Error: engine_nn_proxy.pth not found.", []

        engine_nn = EngineNNProxy()
        engine_nn.load_state_dict(torch.load(nn_pth, map_location='cpu'))
        engine_nn.eval()
        tau_fuel, K_inertia = 0.15, 100.0

        def _bisect(func, a, b, tol=1e-4, maxiter=50):
            fa, fb = func(a), func(b)
            if fa * fb > 0:
                return a if abs(fa) < abs(fb) else b
            for _ in range(maxiter):
                c = (a + b) / 2.0
                fc = func(c)
                if abs(fc) < tol or (b - a) / 2 < tol:
                    return c
                if fa * fc < 0:
                    b, fb = c, fc
                else:
                    a, fa = c, fc
            return (a + b) / 2.0

        def _solve_rpm(target_p):
            def obj(n):
                with torch.no_grad():
                    return engine_nn(torch.tensor([[0.,0.,n]], dtype=torch.float32)).numpy()[0,1] - target_p
            return _bisect(obj, 1000, 58000)

        eng_N = _solve_rpm(initial_engine_power_kw)
        with torch.no_grad():
            pred_init = engine_nn(torch.tensor([[0.,0.,eng_N]], dtype=torch.float32)).numpy()
        eng_Wf_act = pred_init[0, 0]
        eng_Wf_cmd = eng_Wf_act

        eng_use_mpc = (eng_controller_type == "MPC")
        if eng_use_mpc:
            from mpc_controller import TurboShaftMPCController
            eng_mpc = TurboShaftMPCController(
                tau_fuel=tau_fuel, K_inertia=K_inertia, dt=dt,
                horizon=int(np.clip(eng_mpc_horizon, 3, 30)),
                overshoot_limit=float(np.clip(eng_mpc_overshoot, 0.01, 0.30))
            )
            eng_mpc.W_power = float(np.clip(eng_mpc_W_power, 1, 1000))
            eng_mpc.W_dcost = float(np.clip(eng_mpc_W_dcost, 0.01, 50))
            eng_mpc.reset(initial_output=eng_Wf_cmd, initial_N=eng_N)
        else:
            from increPID import IncrementalPIDController
            eng_kp = float(np.clip(eng_kp, 0.01, 30))
            eng_ki = float(np.clip(eng_ki, 0.0, 30))
            eng_kd = float(np.clip(eng_kd, 0.0, 10))
            eng_pid = IncrementalPIDController(
                kp=eng_kp, ki=eng_ki, kd=eng_kd, dt=dt,
                output_min=5.0, output_max=400.0,
                input_scale=300.0, output_scale=400.0
            )
            eng_pid.reset(initial_output=eng_Wf_act)

        # ---- 电机 ----
        mot_use_mpc = (mot_controller_type == "MPC")
        if mot_use_mpc:
            drive_motor = MotorSim(
                P_rate=300e3, w_rate=575.95, J=mot_J,
                mpc_W_speed=float(np.clip(mot_mpc_W_speed, 1, 500)),
                mpc_W_dcost=float(np.clip(mot_mpc_W_dcost, 0.01, 50)),
                mpc_overshoot_limit=float(np.clip(mot_mpc_overshoot, 0.01, 0.30)),
            )
        else:
            drive_motor = MotorSim(
                P_rate=300e3, w_rate=575.95, J=mot_J,
                mpc_W_speed=0., mpc_W_dcost=0., mpc_overshoot_limit=0.05,
            )
            mot_kp = float(np.clip(mot_kp, 0.01, 80))
            mot_ki = float(np.clip(mot_ki, 0.0, 150))
            mot_kd = float(np.clip(mot_kd, 0.0, 10))
            from increPID import IncrementalPIDController
            w_rate = 575.95; tau_rate = 300e3 / w_rate
            mot_pid = IncrementalPIDController(
                kp=mot_kp, ki=mot_ki, kd=mot_kd, dt=dt,
                output_min=-tau_rate, output_max=tau_rate,
                input_scale=w_rate, output_scale=tau_rate
            )
            mot_pid.reset(0.0)

        battery = BatterySim(capacity_kwh=battery_capacity_kwh,
                             initial_soc=initial_soc_pct / 100.0)
        bus_voltage = battery._get_ocv(battery.SOC)
        motor_actual_power_kw = 0.0
        charge_mode = (initial_soc_pct / 100.0) < soc_low
        power_reserve = p_eng_max * power_reserve_pct / 100.0
        points = _profile_points(profile_name)
        time_array = np.arange(0, sim_time_s, dt)

        log = {k: [] for k in [
            'soc', 'bus_voltage', 'prop_speed_rpm', 'target_prop_rpm',
            'target_engine_pwr', 'p_engine_out_kw', 'p_drive_req_kw',
            'p_batt_actual_kw', 'wf_kg_h', 'ems_mode'
        ]}

        n_steps = len(time_array)
        for step_i, t in enumerate(time_array):
            if progress is not None and step_i % max(1, n_steps // 20) == 0:
                progress(step_i / n_steps, desc=f"混动仿真 {step_i}/{n_steps} (t={t:.1f}s)")
            target_rpm, load_torque_t = _target_from_profile(t, points)

            # Motor step
            if not mot_use_mpc:
                w_set_rad = target_rpm * 2 * np.pi / 60.0
                T_cmd = mot_pid.compute(setpoint=w_set_rad, measurement=drive_motor.w_M)
                p_cmd_kw = -T_cmd * max(abs(drive_motor.w_M), 1.0) / 1000.0
                motor_state = drive_motor.step(dt=dt, n_setpoint=target_rpm,
                    p_bus_actual_kw=p_cmd_kw, v_bus=bus_voltage,
                    t_load=load_torque_t, t_ext=0.0)
            else:
                motor_state = drive_motor.step(dt=dt, n_setpoint=target_rpm,
                    p_bus_actual_kw=motor_actual_power_kw, v_bus=bus_voltage,
                    t_load=load_torque_t, t_ext=0.0)

            p_drive_req = motor_state['p_bus_req_kw']
            actual_rpm = motor_state['n_rpm']

            # EMS
            soc = battery.SOC
            if soc < soc_low: charge_mode = True
            elif soc > soc_high: charge_mode = False

            if soc < 0.10:
                target_engine_pwr = p_eng_max; ems_mode_str = "Emergency Charge"
            elif soc > 0.95:
                target_engine_pwr = p_eng_min; ems_mode_str = "Overcharge Prot."
            elif charge_mode:
                target_engine_pwr = p_charge; ems_mode_str = "Charge Mode"
            else:
                soc_error = soc_target - soc
                target_engine_pwr = p_drive_req + power_reserve + soc_error * k_soc
                ems_mode_str = "Power Follow"
            target_engine_pwr = float(np.clip(target_engine_pwr, p_eng_min, p_eng_max))

            # Engine step
            batch_inp = torch.tensor([[0.,0.,eng_N],[0.,0.,eng_N+5.]], dtype=torch.float32)
            with torch.no_grad():
                bp = engine_nn(batch_inp).numpy()
            Wf_req, P_eng_out = bp[0,0], bp[0,1]

            if eng_use_mpc:
                k_wf = (bp[1,0]-bp[0,0])/5.0; k_p = (bp[1,1]-bp[0,1])/5.0
                eng_Wf_cmd = eng_mpc.compute(current_N=eng_N, current_Wfact=eng_Wf_act,
                    target_power=target_engine_pwr, precalc_params=(Wf_req, P_eng_out, k_wf, k_p))
            else:
                eng_Wf_cmd = eng_pid.compute(setpoint=target_engine_pwr, measurement=P_eng_out)

            eng_Wf_act += (eng_Wf_cmd - eng_Wf_act) / tau_fuel * dt
            eng_N += K_inertia * (eng_Wf_act - Wf_req) * dt

            # Battery
            p_batt_req = p_drive_req - P_eng_out
            p_batt_actual, v_bus, i_batt, soc_new = battery.step(dt, p_batt_req)
            bus_voltage = v_bus
            motor_actual_power_kw = P_eng_out + p_batt_actual

            log['soc'].append(soc_new * 100.0)
            log['bus_voltage'].append(v_bus)
            log['prop_speed_rpm'].append(actual_rpm)
            log['target_prop_rpm'].append(target_rpm)
            log['target_engine_pwr'].append(target_engine_pwr)
            log['p_engine_out_kw'].append(P_eng_out)
            log['p_drive_req_kw'].append(p_drive_req)
            log['p_batt_actual_kw'].append(p_batt_actual)
            log['wf_kg_h'].append(eng_Wf_act)
            log['ems_mode'].append(ems_mode_str)

        if progress is not None:
            progress(1.0, desc="绘图中...")
        # ---- Plots (English) ----
        soc_arr = np.array(log['soc'])
        speed_error = np.array(log['target_prop_rpm']) - np.array(log['prop_speed_rpm'])

        fig, axes = plt.subplots(4, 1, figsize=(12, 14), sharex=True)
        fig.suptitle('Hybrid EMS Validation (Rule-Based + Hysteresis)',
                     fontweight='bold', fontsize=13)

        mode_colors = {'Power Follow': '#E3F2FD', 'Charge Mode': '#FFEBEE',
                       'Emergency Charge': '#FFCDD2', 'Overcharge Prot.': '#E8F5E9'}

        axes[0].plot(time_array, log['target_prop_rpm'], 'k--', lw=1.5, label='Target')
        axes[0].plot(time_array, log['prop_speed_rpm'], 'b-', lw=1.5, label='Actual')
        axes[0].set_ylabel('Speed (RPM)'); axes[0].set_title('Propulsion Speed')
        axes[0].grid(True, linestyle=':'); axes[0].legend()

        modes = log['ems_mode']
        i = 0; added = set()
        while i < len(modes):
            m = modes[i]; j = i
            while j < len(modes) and modes[j] == m: j += 1
            col = mode_colors.get(m, '#F5F5F5')
            lbl = m if m not in added else None
            axes[1].axvspan(time_array[i], time_array[min(j-1, len(time_array)-1)],
                            alpha=0.3, color=col, label=lbl)
            if lbl: added.add(m)
            i = j

        axes[1].plot(time_array, log['p_drive_req_kw'], 'k--', lw=1.2, label='Motor Demand')
        axes[1].plot(time_array, log['target_engine_pwr'], color='darkred', ls=':', lw=1, label='Eng Target')
        axes[1].plot(time_array, log['p_engine_out_kw'], 'r-', lw=1.2, label='Eng Output')
        axes[1].plot(time_array, log['p_batt_actual_kw'], 'g-', lw=1.2, label='Battery')
        axes[1].axhline(0, color='gray', lw=0.8)
        axes[1].set_ylabel('Power (kW)'); axes[1].set_title('Power Allocation (bg=EMS mode)')
        axes[1].grid(True, linestyle=':'); axes[1].legend(ncol=3, fontsize=8, loc='upper right')

        axes[2].plot(time_array, log['bus_voltage'], 'm-', lw=1.2, label='Bus Voltage')
        axes[2].set_ylabel('Voltage (V)'); axes[2].set_title('Electrical State & SOC')
        axes[2].grid(True, linestyle=':'); axes[2].legend(loc='upper left')
        ax_soc = axes[2].twinx()
        ax_soc.plot(time_array, log['soc'], 'c-', lw=1.6, label='SOC')
        ax_soc.axhline(soc_low*100, color='r', ls=':', lw=1, alpha=0.7, label=f'Low ({soc_low*100:.0f}%)')
        ax_soc.axhline(soc_high*100, color='g', ls=':', lw=1, alpha=0.7, label=f'High ({soc_high*100:.0f}%)')
        ax_soc.axhline(soc_target*100, color='b', ls='-.', lw=1, alpha=0.5, label=f'Target ({soc_target*100:.0f}%)')
        ax_soc.set_ylabel('SOC (%)'); ax_soc.legend(loc='upper right', fontsize=8)

        axes[3].plot(time_array, log['wf_kg_h'], 'tab:orange', lw=1.2, label='Fuel Flow')
        axes[3].set_ylabel('Fuel (kg/h)'); axes[3].set_xlabel('Time (s)')
        axes[3].set_title('Fuel Consumption')
        axes[3].grid(True, linestyle=':'); axes[3].legend()

        fig.tight_layout(rect=[0, 0, 1, 0.96])

        mode_times = {}
        for m_ in modes: mode_times[m_] = mode_times.get(m_, 0) + dt
        mode_str = ', '.join([f'{k}: {v:.1f}s' for k, v in mode_times.items()])

        fuel_arr = np.array(log['wf_kg_h'])
        engine_pwr_arr = np.array(log['p_engine_out_kw'])
        batt_pwr_arr = np.array(log['p_batt_actual_kw'])

        summary = (
            f"### Hybrid Simulation Summary\n"
            f"- **Duration**: {sim_time_s:.0f}s, dt={dt:.3f}s\n"
            f"- **Controllers**: Engine={eng_controller_type}, Motor={mot_controller_type}\n"
            f"- **Max Speed Error**: {np.max(np.abs(speed_error)):.1f} RPM\n"
            f"- **SOC**: {soc_arr[0]:.1f}% → {soc_arr[-1]:.1f}% "
            f"(min {np.min(soc_arr):.1f}%, max {np.max(soc_arr):.1f}%)\n"
            f"- **Avg Engine Power**: {np.mean(engine_pwr_arr):.1f} kW\n"
            f"- **Avg Battery**: {np.mean(batt_pwr_arr):.1f} kW (+discharge/−charge)\n"
            f"- **Avg Fuel**: {np.mean(fuel_arr):.1f} kg/h\n"
            f"- **EMS Modes**: {mode_str}"
        )

        pick_idx = np.linspace(0, len(time_array)-1, 8, dtype=int)
        table_data = [[
            round(float(time_array[idx]),2),
            round(float(log['target_prop_rpm'][idx]),1),
            round(float(log['prop_speed_rpm'][idx]),1),
            round(float(log['p_engine_out_kw'][idx]),2),
            round(float(log['p_batt_actual_kw'][idx]),2),
            round(float(log['soc'][idx]),2),
            log['ems_mode'][idx],
        ] for idx in pick_idx]

        return fig, summary, table_data

    except Exception as e:
        import traceback
        return None, f"Hybrid simulation failed: {e}\n```\n{traceback.format_exc()}\n```", []