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```", []