diff --git a/.gradio/certificate.pem b/.gradio/certificate.pem new file mode 100644 index 0000000..b85c803 --- /dev/null +++ b/.gradio/certificate.pem @@ -0,0 +1,31 @@ +-----BEGIN CERTIFICATE----- +MIIFazCCA1OgAwIBAgIRAIIQz7DSQONZRGPgu2OCiwAwDQYJKoZIhvcNAQELBQAw +TzELMAkGA1UEBhMCVVMxKTAnBgNVBAoTIEludGVybmV0IFNlY3VyaXR5IFJlc2Vh +cmNoIEdyb3VwMRUwEwYDVQQDEwxJU1JHIFJvb3QgWDEwHhcNMTUwNjA0MTEwNDM4 +WhcNMzUwNjA0MTEwNDM4WjBPMQswCQYDVQQGEwJVUzEpMCcGA1UEChMgSW50ZXJu +ZXQgU2VjdXJpdHkgUmVzZWFyY2ggR3JvdXAxFTATBgNVBAMTDElTUkcgUm9vdCBY +MTCCAiIwDQYJKoZIhvcNAQEBBQADggIPADCCAgoCggIBAK3oJHP0FDfzm54rVygc +h77ct984kIxuPOZXoHj3dcKi/vVqbvYATyjb3miGbESTtrFj/RQSa78f0uoxmyF+ +0TM8ukj13Xnfs7j/EvEhmkvBioZxaUpmZmyPfjxwv60pIgbz5MDmgK7iS4+3mX6U +A5/TR5d8mUgjU+g4rk8Kb4Mu0UlXjIB0ttov0DiNewNwIRt18jA8+o+u3dpjq+sW +T8KOEUt+zwvo/7V3LvSye0rgTBIlDHCNAymg4VMk7BPZ7hm/ELNKjD+Jo2FR3qyH +B5T0Y3HsLuJvW5iB4YlcNHlsdu87kGJ55tukmi8mxdAQ4Q7e2RCOFvu396j3x+UC +B5iPNgiV5+I3lg02dZ77DnKxHZu8A/lJBdiB3QW0KtZB6awBdpUKD9jf1b0SHzUv +KBds0pjBqAlkd25HN7rOrFleaJ1/ctaJxQZBKT5ZPt0m9STJEadao0xAH0ahmbWn +OlFuhjuefXKnEgV4We0+UXgVCwOPjdAvBbI+e0ocS3MFEvzG6uBQE3xDk3SzynTn +jh8BCNAw1FtxNrQHusEwMFxIt4I7mKZ9YIqioymCzLq9gwQbooMDQaHWBfEbwrbw +qHyGO0aoSCqI3Haadr8faqU9GY/rOPNk3sgrDQoo//fb4hVC1CLQJ13hef4Y53CI +rU7m2Ys6xt0nUW7/vGT1M0NPAgMBAAGjQjBAMA4GA1UdDwEB/wQEAwIBBjAPBgNV +HRMBAf8EBTADAQH/MB0GA1UdDgQWBBR5tFnme7bl5AFzgAiIyBpY9umbbjANBgkq +hkiG9w0BAQsFAAOCAgEAVR9YqbyyqFDQDLHYGmkgJykIrGF1XIpu+ILlaS/V9lZL +ubhzEFnTIZd+50xx+7LSYK05qAvqFyFWhfFQDlnrzuBZ6brJFe+GnY+EgPbk6ZGQ +3BebYhtF8GaV0nxvwuo77x/Py9auJ/GpsMiu/X1+mvoiBOv/2X/qkSsisRcOj/KK +NFtY2PwByVS5uCbMiogziUwthDyC3+6WVwW6LLv3xLfHTjuCvjHIInNzktHCgKQ5 +ORAzI4JMPJ+GslWYHb4phowim57iaztXOoJwTdwJx4nLCgdNbOhdjsnvzqvHu7Ur +TkXWStAmzOVyyghqpZXjFaH3pO3JLF+l+/+sKAIuvtd7u+Nxe5AW0wdeRlN8NwdC +jNPElpzVmbUq4JUagEiuTDkHzsxHpFKVK7q4+63SM1N95R1NbdWhscdCb+ZAJzVc +oyi3B43njTOQ5yOf+1CceWxG1bQVs5ZufpsMljq4Ui0/1lvh+wjChP4kqKOJ2qxq +4RgqsahDYVvTH9w7jXbyLeiNdd8XM2w9U/t7y0Ff/9yi0GE44Za4rF2LN9d11TPA +mRGunUHBcnWEvgJBQl9nJEiU0Zsnvgc/ubhPgXRR4Xq37Z0j4r7g1SgEEzwxA57d +emyPxgcYxn/eR44/KJ4EBs+lVDR3veyJm+kXQ99b21/+jh5Xos1AnX5iItreGCc= +-----END CERTIFICATE----- diff --git a/Model/data/engine_gpr_model.pth b/Model/data/engine_gpr_model.pth index e9dc524..e9238e2 100644 Binary files a/Model/data/engine_gpr_model.pth and b/Model/data/engine_gpr_model.pth differ diff --git a/Model/data/engine_nn_proxy.pth b/Model/data/engine_nn_proxy.pth new file mode 100644 index 0000000..d580d19 Binary files /dev/null and b/Model/data/engine_nn_proxy.pth differ diff --git a/Model/src/distill_gpr_to_nn.py b/Model/src/distill_gpr_to_nn.py new file mode 100644 index 0000000..1e04202 --- /dev/null +++ b/Model/src/distill_gpr_to_nn.py @@ -0,0 +1,275 @@ +""" +GPR → NN 知识蒸馏脚本 +将训练好的高斯过程回归 (GPR) 代理模型的知识蒸馏到轻量级神经网络 (MLP) 中。 + +支持两种模式: + 1. 直接从 CSV 原始数据训练 NN (快速模式, 无需 GPR 依赖) + 2. 从 GPR 教师模型蒸馏 (完整模式, 需 botorch/gpytorch) + +用法: + python distill_gpr_to_nn.py # 从Model目录运行 + python Model/src/distill_gpr_to_nn.py # 从项目根目录运行 + +也可在代码中调用: + from distill_gpr_to_nn import distill, distill_from_csv + result = distill_from_csv(csv_path, nn_save_path, ...) + result = distill(gpr_csv_path, gpr_pth_path, nn_save_path, ...) +""" + +import os +import sys +import csv as csv_mod +import numpy as np +import torch +import torch.nn as nn + +# 确保 src 目录在 path 中 +_this_dir = os.path.dirname(os.path.abspath(__file__)) +if _this_dir not in sys.path: + sys.path.insert(0, _this_dir) + +from lightweight_model import EngineNNProxy + + +# ============================================================ +# 工具函数: 手写 StandardScaler (避免 sklearn import 导致的 +# numpy/scipy 递归问题) +# ============================================================ +class SimpleScaler: + """轻量 Z-Score 归一化, 兼容 EngineNNProxy.set_normalization_params""" + def __init__(self): + self.mean_ = None + self.scale_ = None + + def fit(self, X): + X = np.asarray(X, dtype=np.float64) + self.mean_ = X.mean(axis=0) + self.scale_ = X.std(axis=0) + self.scale_[self.scale_ < 1e-12] = 1.0 + return self + + def transform(self, X): + return (np.asarray(X, dtype=np.float64) - self.mean_) / self.scale_ + + def fit_transform(self, X): + self.fit(X) + return self.transform(X) + + def inverse_transform(self, X): + return np.asarray(X, dtype=np.float64) * self.scale_ + self.mean_ + + +# ============================================================ +# 模式一: 直接从 CSV 训练 (快速, 无需 sklearn/botorch) +# ============================================================ +def distill_from_csv(csv_path, nn_save_path, + epochs=3000, lr=1e-3, batch_size=256, + hidden_size=64, verbose=True, progress_callback=None): + """ + 直接从 CSV 原始发动机数据训练 NN 代理模型。 + + Returns + ------- + dict with keys: loss_history, nn_model, X_train, Y_train, Y_nn, + scaler_X, scaler_Y, rel_error_fuel, rel_error_power + """ + # 1. 读取 CSV + if verbose: + print(f"-> Loading CSV: {csv_path}") + 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]] + + # 过滤零/非物理值 + valid = (Y_raw[:, 0] > 0.5) & (Y_raw[:, 1] > 0.5) & (X_raw[:, 2] > 500) + X_phys = X_raw[valid].astype(np.float32) + Y_phys = Y_raw[valid].astype(np.float32) + if verbose: + print(f"-> Valid samples: {len(X_phys)} / {len(data)}") + + # 2. 归一化 + scaler_X = SimpleScaler() + scaler_X.fit(X_phys) + + Y_log = np.log1p(Y_phys.astype(np.float64)) + scaler_Y = SimpleScaler() + scaler_Y.fit(Y_log) + + return _train_nn(X_phys, Y_phys, scaler_X, scaler_Y, nn_save_path, + epochs, lr, batch_size, hidden_size, verbose, + progress_callback=progress_callback) + + +# ============================================================ +# 模式二: 从 GPR 教师模型蒸馏 (需 botorch/gpytorch/sklearn) +# ============================================================ +def _generate_teacher_data(gpr_csv_path, gpr_pth_path, + n_altitude=25, n_mach=6, n_rpm=30): + """ + 加载 GPR 教师模型，在输入空间的密集网格上生成标注数据。 + """ + from engine_gpr_class import EngineGPRModel + + gpr = EngineGPRModel(csv_path=gpr_csv_path) + success = gpr.load_model(pth_path=gpr_pth_path) + if not success: + raise RuntimeError(f"Failed to load GPR model from {gpr_pth_path}") + + df = gpr.df + H_range = np.linspace(df['Altitude_m'].min(), df['Altitude_m'].max(), n_altitude) + Ma_range = np.linspace(df['Mach'].min(), df['Mach'].max(), n_mach) + RPM_range = np.linspace(max(df['RPM'].min(), 1000), df['RPM'].max(), n_rpm) + + H, Ma, RPM = np.meshgrid(H_range, Ma_range, RPM_range, indexing='ij') + X_grid = np.column_stack([H.ravel(), Ma.ravel(), RPM.ravel()]) + + print(f"-> Querying GPR teacher on {len(X_grid)} grid points ...") + Y_pred, _ = gpr.predict(X_grid) # [WF_kg_h, Power_kW] + + valid = (Y_pred[:, 0] > 0.5) & (Y_pred[:, 1] > 0.5) + X_valid = X_grid[valid].astype(np.float32) + Y_valid = Y_pred[valid].astype(np.float32) + print(f"-> Valid samples: {len(X_valid)} / {len(X_grid)}") + + scaler_X = SimpleScaler() + scaler_X.fit(X_valid) + + Y_log = np.log1p(Y_valid.astype(np.float64)) + scaler_Y = SimpleScaler() + scaler_Y.fit(Y_log) + + return X_valid, Y_valid, scaler_X, scaler_Y + + +def distill(gpr_csv_path, gpr_pth_path, nn_save_path, + n_altitude=25, n_mach=6, n_rpm=30, + epochs=3000, lr=1e-3, batch_size=512, + hidden_size=64, verbose=True): + """ + 从 GPR 教师模型蒸馏到 NN。需安装 botorch / gpytorch。 + """ + X_phys, Y_phys, scaler_X, scaler_Y = _generate_teacher_data( + gpr_csv_path, gpr_pth_path, n_altitude, n_mach, n_rpm + ) + return _train_nn(X_phys, Y_phys, scaler_X, scaler_Y, nn_save_path, + epochs, lr, batch_size, hidden_size, verbose) + + +# ============================================================ +# 通用 NN 训练核心 +# ============================================================ +def _train_nn(X_phys, Y_phys, scaler_X, scaler_Y, nn_save_path, + epochs=3000, lr=1e-3, batch_size=256, + hidden_size=64, verbose=True, progress_callback=None): + """训练 NN 并保存模型。""" + device = torch.device("cuda" if torch.cuda.is_available() else "cpu") + if verbose: + print(f"-> Training device: {device}") + + # 创建模型并移至设备 + nn_model = EngineNNProxy(hidden_size=hidden_size) + nn_model.set_normalization_params(scaler_X, scaler_Y) + nn_model = nn_model.to(device) + + # 准备归一化数据 + X_norm = scaler_X.transform(X_phys).astype(np.float32) + Y_log = np.log1p(Y_phys.astype(np.float64)) + Y_norm = scaler_Y.transform(Y_log).astype(np.float32) + + X_t = torch.tensor(X_norm, dtype=torch.float32) + Y_t = torch.tensor(Y_norm, dtype=torch.float32) + + dataset = torch.utils.data.TensorDataset(X_t, Y_t) + loader = torch.utils.data.DataLoader(dataset, batch_size=batch_size, shuffle=True) + + # 训练 + optimizer = torch.optim.Adam(nn_model.net.parameters(), lr=lr) + scheduler = torch.optim.lr_scheduler.CosineAnnealingLR( + optimizer, T_max=epochs, eta_min=lr * 0.01) + loss_fn = nn.MSELoss() + + loss_history = [] + nn_model.train() + + n_params = sum(p.numel() for p in nn_model.parameters()) + if verbose: + print(f"-> Training NN ({hidden_size}x{hidden_size}, {n_params} params) " + f"for {epochs} epochs on {len(X_phys)} samples ...") + + for epoch in range(epochs): + epoch_loss = 0.0 + n_batches = 0 + for xb, yb in loader: + xb, yb = xb.to(device), yb.to(device) + pred = nn_model.net(xb) + loss = loss_fn(pred, yb) + optimizer.zero_grad() + loss.backward() + optimizer.step() + epoch_loss += loss.item() + n_batches += 1 + + scheduler.step() + avg_loss = epoch_loss / max(n_batches, 1) + loss_history.append(avg_loss) + + if verbose and (epoch + 1) % 500 == 0: + print(f" Epoch {epoch+1:>5d}/{epochs} Loss: {avg_loss:.6f}") + if progress_callback is not None and (epoch + 1) % 50 == 0: + progress_callback(epoch + 1, epochs, avg_loss) + + nn_model.eval() + + # 验证精度（在 CPU 上做推理，避免显存占用） + nn_model_cpu = nn_model.cpu() + with torch.no_grad(): + X_phys_t = torch.tensor(X_phys, dtype=torch.float32) + Y_nn = nn_model_cpu(X_phys_t).numpy() + + rel_err_fuel = np.mean( + np.abs(Y_nn[:, 0] - Y_phys[:, 0]) / np.maximum(Y_phys[:, 0], 1e-6)) * 100 + rel_err_power = np.mean( + np.abs(Y_nn[:, 1] - Y_phys[:, 1]) / np.maximum(Y_phys[:, 1], 1e-6)) * 100 + + if verbose: + print(f"-> Distillation complete!") + print(f" Fuel Flow MAPE: {rel_err_fuel:.2f}%") + print(f" Power MAPE: {rel_err_power:.2f}%") + + # 保存（始终保存 CPU 版，推理时无需 GPU 环境） + if nn_save_path: + os.makedirs(os.path.dirname(os.path.abspath(nn_save_path)), exist_ok=True) + torch.save(nn_model_cpu.state_dict(), nn_save_path) + if verbose: + print(f"-> Saved NN model to {nn_save_path}") + + return { + 'loss_history': loss_history, + 'nn_model': nn_model_cpu, + 'X_train': X_phys, + 'Y_train': Y_phys, + 'Y_nn': Y_nn, + 'scaler_X': scaler_X, + 'scaler_Y': scaler_Y, + 'rel_error_fuel': rel_err_fuel, + 'rel_error_power': rel_err_power, + } + + +if __name__ == "__main__": + model_root = os.path.join(os.path.dirname(os.path.abspath(__file__)), "..") + csv_path = os.path.join(model_root, "data", "Cleaned_Engine_Data_Full.csv") + nn_path = os.path.join(model_root, "data", "engine_nn_proxy.pth") + + # 使用 CSV 直接训练模式 (无需 sklearn/botorch) + result = distill_from_csv(csv_path, nn_path, epochs=3000) + print(f"\nFinal training loss: {result['loss_history'][-1]:.6f}") diff --git a/Model/src/engine_dynamic_sim.py b/Model/src/engine_dynamic_sim.py index d2e52a7..3fd28db 100644 --- a/Model/src/engine_dynamic_sim.py +++ b/Model/src/engine_dynamic_sim.py @@ -1,100 +1,72 @@ import sys import os + sys.path.insert(0, os.path.dirname(os.path.dirname(os.path.abspath(__file__)))) import numpy as np +import torch import matplotlib.pyplot as plt -from src.engine_gpr_class import EngineGPRModel -from src.increPID import IncrementalPIDController +from tqdm import tqdm +from src.lightweight_model import EngineNNProxy +from src.mpc_controller import TurboShaftMPCController class TurboshaftDynamicSim: """ 涡轴发动机动态仿真类 - 包含一阶燃油执行机构和转子动力学模型 + + 包含功能: + 1. 基于 NN Proxy (轻量化神经网络) 的稳态代理模型 + 2. 一阶燃油执行机构动态 (First-order Actuator Dynamics) + 3. 转子动力学积分 (Rotor Dynamics) + 4. MPC (模型预测控制) 用于功率跟随 """ - def __init__(self, gpr_csv_path="data/Cleaned_Engine_Data_Full.csv", gpr_pth_path="data/engine_gpr_model.pth", + def __init__(self, nn_pth_path="data/engine_nn_proxy.pth", tau_fuel=0.15, K_inertia=100.0, - kp=4.652, ki=7.078, kd=0.222, min_fuel=10.0, max_fuel=600.0, - verbose=False): - """ - 初始化仿真环境和代理模型 - :param gpr_csv_path: GPR模型使用的数据集路径 - :param gpr_pth_path: 预训练的模型参数路径 - :param tau_fuel: 燃油执行机构时间常数 (s) - :param K_inertia: 转子惯性增益 (RPM / (kg/h)) - :param kp: PID比例系数 (用于功率跟随控制，归一化域) - :param ki: PID积分系数 (归一化域) - :param kd: PID微分系数 (归一化域) - :param min_fuel: 燃油流量下限 (kg/h) - :param max_fuel: 燃油流量上限 (kg/h) - :param verbose: 是否打印详细信息 - """ - import os - self.verbose = verbose - - if verbose: - print("-> 正在加载 GPR 稳态代理模型...") - self.engine_model = EngineGPRModel(gpr_csv_path) - - # 检查模型文件是否存在，不存在则训练 - if not os.path.exists(gpr_pth_path): - if verbose: - print(f"-> 模型文件 {gpr_pth_path} 不存在，正在训练模型...") - self.engine_model.train(save_path=gpr_pth_path) - else: - success = self.engine_model.load_model(gpr_pth_path) - if not success: - if verbose: - print(f"-> 模型加载失败，正在重新训练...") - self.engine_model.train(save_path=gpr_pth_path) + mpc_horizon=15, mpc_dt=0.02, min_fuel=10.0, max_fuel=400.0, + mpc_overshoot_limit=0.05): + print("-> 正在加载 NN 稳态代理模型...") + self.engine_model = EngineNNProxy() + + if not os.path.exists(nn_pth_path): + raise FileNotFoundError(f"Missing NN model file: {nn_pth_path}. Please run scripts/distill_gpr_to_nn.py first.") + + self.engine_model.load_state_dict(torch.load(nn_pth_path, map_location='cpu')) + self.engine_model.eval() + print("-> NN Model Loaded (CPU Mode for Sim Loop)") - # 动态参数 self.tau_fuel = tau_fuel self.K_inertia = K_inertia - - # 归一化基准 (用于PID计算) self.max_fuel = max_fuel - self.max_power_ref = 300.0 # 300kw为功率归一化基准，根据实际数据调整 + self.max_power_ref = 300.0 - # 控制器初始化 (使用归一化参数，通过 scaling 自动处理) - self.pid = IncrementalPIDController( - kp=kp, ki=ki, kd=kd, dt=0.01, - output_min=min_fuel, output_max=max_fuel, - input_scale=self.max_power_ref, output_scale=self.max_fuel + self.mpc = TurboShaftMPCController( + tau_fuel=tau_fuel, K_inertia=K_inertia, dt=mpc_dt, horizon=mpc_horizon, + min_fuel=min_fuel, max_fuel=max_fuel, overshoot_limit=mpc_overshoot_limit ) - # 环境与状态变量 self.H_env = 0.0 self.Ma_env = 0.0 self.N_current = 0.0 self.Wf_act_current = 0.0 self.power_generated = 0.0 - - # 控制指令 self.Wf_cmd = 0.0 def _solve_steady_rpm(self, target_val, target_type='power'): - """ - 内部方法：反解稳态转速 (RPM) - :param target_val: 目标值 (Power[kW] 或 Fuel[kg/h]) - :param target_type: 'power' 或 'fuel' - :return: 对应的稳态转速 - """ + """数值反解给定功率/燃油下的稳态转速""" from scipy.optimize import brentq - - # 搜索范围 [RPM_min, RPM_max]，根据经验或数据范围设定 low_bound, high_bound = 0.0, 60000.0 def objective(n): - current_input = np.array([[self.H_env, self.Ma_env, n]]) - pred_mean, _ = self.engine_model.predict(current_input) + # 构造输入: [H, Ma, N] + current_input = torch.tensor([[self.H_env, self.Ma_env, n]], dtype=torch.float32) + with torch.no_grad(): + pred_mean = self.engine_model(current_input).numpy() - # index 0: Fuel Flow (kg/h), index 1: Power (kW) + # 输出: [0]=Fuel, [1]=Power val = pred_mean[0, 1] if target_type == 'power' else pred_mean[0, 0] return val - target_val try: - # 简单的边界检查，防止报错 f_low = objective(low_bound) f_high = objective(high_bound) if f_low * f_high > 0: @@ -108,118 +80,104 @@ class TurboshaftDynamicSim: return (low_bound + high_bound) / 2.0 def set_steady_state_by_power(self, H_env, Ma_env, Power_target): - """ - 通过目标功率初始化稳态 - :param Power_target: 目标轴功率 (kW) - :return: 对应的稳态转速 (RPM) - """ + """设定初始稳态工况点""" self.H_env = H_env self.Ma_env = Ma_env - - # 反解转速 self.N_current = self._solve_steady_rpm(Power_target, target_type='power') - # 计算该转速下的稳态燃油 - current_input = np.array([[self.H_env, self.Ma_env, self.N_current]]) - pred_mean, _ = self.engine_model.predict(current_input) + # 计算该稳态下的燃油消耗 + current_input = torch.tensor([[self.H_env, self.Ma_env, self.N_current]], dtype=torch.float32) + with torch.no_grad(): + pred_mean = self.engine_model(current_input).numpy() self.power_generated = pred_mean[0, 1] self.Wf_act_current = pred_mean[0, 0] self.Wf_cmd = self.Wf_act_current - print(f"-> 稳态(Power)已配置: H={H_env}, Ma={Ma_env}, Target_P={Power_target:.1f} kW => N={self.N_current:.1f} RPM, Wf={self.Wf_cmd:.2f} kg/h") if self.verbose else None - - # 重置PID控制器至当前稳态输出 (自动处理归一化) - self.pid.reset(initial_output=self.Wf_cmd) - + print(f"-> 稳态(Power)已配置: H={H_env}, Ma={Ma_env}, Target_P={Power_target:.1f} kW => N={self.N_current:.1f} RPM, Wf={self.Wf_cmd:.2f} kg/h") + self.mpc.reset(initial_output=self.Wf_cmd, initial_N=self.N_current) return self.N_current - def set_steady_state_by_fuel(self, H_env, Ma_env, Wf_target): - """ - 通过目标燃油流量初始化稳态 - :param Wf_target: 目标燃油流量 (kg/h) - :return: 对应的稳态转速 (RPM) - """ - self.H_env = H_env - self.Ma_env = Ma_env - - # 反解转速 - self.N_current = self._solve_steady_rpm(Wf_target, target_type='fuel') - - # 确认该状态下的功率 - current_input = np.array([[self.H_env, self.Ma_env, self.N_current]]) - pred_mean, _ = self.engine_model.predict(current_input) - - self.power_generated = pred_mean[0, 1] - self.Wf_act_current = pred_mean[0, 0] # 理论上应该非常接近 Wf_target - self.Wf_cmd = self.Wf_act_current - - print(f"-> 稳态(Fuel)已配置: H={H_env}, Ma={Ma_env}, Target_Wf={Wf_target:.1f} kg/h => N={self.N_current:.1f} RPM, Power={self.power_generated:.1f} kW") if self.verbose else None - - # 重置PID控制器至当前稳态输出 (自动处理归一化) - self.pid.reset(initial_output=self.Wf_cmd) - - return self.N_current - - def set_fuel_command(self, Wf_cmd): - """ - 更改燃油流量指令 - :param Wf_cmd: 目标燃油流量指令 - """ - self.Wf_cmd = Wf_cmd - - def set_flight_condition(self, H_env=None, Ma_env=None): - """ - 在运行过程中更改当前的飞行条件 (高度和马赫数) - :param H_env: 新的飞行高度 (m)。如果为 None，则保持不变。 - :param Ma_env: 新的飞行马赫数。如果为 None，则保持不变。 - """ - if H_env is not None: - self.H_env = H_env - if Ma_env is not None: - self.Ma_env = Ma_env - - def compute_control_law(self, dt, target_power): - """ - 计算控制律 (PID控制: Power -> Wf) - :param dt: 控制周期 (s) - :param target_power: 期望功率 (kW) - :return: 计算出的燃油指令 - """ - self.pid.dt = dt - # 1. 计算控制增量 (PID内部会自动处理归一化) - self.Wf_cmd = self.pid.compute(setpoint=target_power, measurement=self.power_generated) + def compute_control_law(self, dt, target_power, precalc_params=None): + """调用 MPC 更新控制指令""" + self.Wf_cmd = self.mpc.compute( + current_N=self.N_current, + current_Wfact=self.Wf_act_current, + target_power=target_power, + precalc_params=precalc_params + ) return self.Wf_cmd def step(self, dt, target_power=None): """ - 执行单步动态仿真 - :param dt: 积分步长 (s) - :param target_power: 目标轴功率 (kW)，若不为 None 则执行一次PID控制 - :return: (当前转速, 实际供油量, 当前需要的平衡供油量, 当前功率) + 执行单步动态仿真 (High-Performance Optimized) + + 加速策略: + - 聚合 GPR 预测请求: 将 MPC 所需的梯度计算点与当前物理状态点合并为一个 Batch (Size=2) + - 减少 GPU I/O 次数: 从每步 3 次减少为 1 次 """ - # 0. 闭环控制计算 - if target_power is not None: - self.compute_control_law(dt, target_power) + + # --- 0. 统一 GPU 批次预测 (Batch Prediction) --- + # 构造输入: [Row 0: 当前状态点, Row 1: 用于梯度计算的微扰点] + delta_N = 5.0 - # 1. 燃油执行机构动态 (一阶惯性) + # 判断模型类型: NN 模型使用 forward，GPR 模型使用 predict + if hasattr(self.engine_model, 'predict') and not hasattr(self.engine_model, 'forward'): + # GPR 模型 + inputs = np.array([ + [self.H_env, self.Ma_env, self.N_current], + [self.H_env, self.Ma_env, self.N_current + delta_N] + ]) + pred_mean, _ = self.engine_model.predict(inputs) + else: + # NN 模型 + import torch + inputs = torch.tensor([ + [self.H_env, self.Ma_env, self.N_current], + [self.H_env, self.Ma_env, self.N_current + delta_N] + ], dtype=torch.float32) + with torch.no_grad(): + pred_mean = self.engine_model(inputs).numpy() + + # 核心加速点：一次 GPU 调用获取所有信息 + # pred_mean 形如 [[Wf0, Pow0], [Wf1, Pow1]] + + # 提取结果 + Wf_req_current = pred_mean[0, 0] + Power_current = pred_mean[0, 1] + + Wf_req_pert = pred_mean[1, 0] + Power_pert = pred_mean[1, 1] + + # --- 1. 闭环控制计算 --- + if target_power is not None: + # 在 Python 端快速计算梯度，避免在 MPC 内部再次调用模型 + k_wf = (Wf_req_pert - Wf_req_current) / delta_N + k_p = (Power_pert - Power_current) / delta_N + + # 使用预计算好的参数，MPC 内部将不再调用 engine_model.predict + params = (Wf_req_current, Power_current, k_wf, k_p) + self.compute_control_law(dt, target_power, precalc_params=params) + + # --- 2. 燃油执行机构动态 (一阶惯性) --- dWf_act_dt = (self.Wf_cmd - self.Wf_act_current) / self.tau_fuel Wf_act_next = self.Wf_act_current + dWf_act_dt * dt - # 2. 调用GPR代理模型计算当前转速下的阻力矩(需求燃油) - current_input = np.array([[self.H_env, self.Ma_env, self.N_current]]) - pred_mean, _ = self.engine_model.predict(current_input) - - # GPR 输出: [0]: Fuel Flow, [1]: Shaft Power - Wf_req_current = pred_mean[0, 0] #维持当前转速所需的稳态燃油 + # --- 3. 调用NN代理模型推算当前气动热力参数 --- + current_input = torch.tensor([[self.H_env, self.Ma_env, self.N_current]], dtype=torch.float32) + with torch.no_grad(): + pred_mean = self.engine_model(current_input).numpy() + + # [0]: Fuel Flow (kg/h), [1]: Power (kW) + Wf_req_current = pred_mean[0, 0] Power_current = pred_mean[0, 1] - self.power_generated = Power_current # 更新当前功率状态 + self.power_generated = Power_current - # 3. 转子动力学积分 (燃料差额 -> 转速加速度) + # --- 4. 转子动力学积分 (供油盈余 -> 加速) --- dN_dt = self.K_inertia * (self.Wf_act_current - Wf_req_current) N_next = self.N_current + dN_dt * dt - # 4. 状态更新 + # 更新状态 self.Wf_act_current = Wf_act_next self.N_current = N_next @@ -228,7 +186,8 @@ class TurboshaftDynamicSim: if __name__ == "__main__": # ========================================== - # 【测试示例】利用类运行功率闭环控制仿真 + # 涡轴发动机动态响应测试脚本 + # 模拟复杂剖面: 包含阶跃、正弦、斜坡指令及变高度/马赫数干扰 # ========================================== import matplotlib matplotlib.use('Agg') @@ -238,121 +197,105 @@ if __name__ == "__main__": plt.rcParams['font.serif'] = ['DejaVu Serif', 'Times New Roman'] plt.rcParams['axes.unicode_minus'] = True - # 初始化仿真 (优化后的PID参数) - # kp: 比例系数 - 增大以加快响应速度 - # ki: 积分系数 - 适中以消除稳态误差 - # kd: 微分系数 - 设为0避免控制振荡 - sim = TurboshaftDynamicSim() + sim = TurboshaftDynamicSim(mpc_dt=0.02) - # 1. 设置初始稳态 (通过目标功率设定) - P_initial_target = 100.0 # kW + # 初始状态 + P_initial_target = 100.0 sim.set_steady_state_by_power(H_env=0.0, Ma_env=0.0, Power_target=P_initial_target) - # 仿真参数 dt = 0.02 - t_end = 30.0 # 增加仿真时间以展示更多指令变化 + t_end = 60.0 time_array = np.arange(0, t_end, dt) - # 数据记录 - N_log = [] - Wf_act_log = [] - Wf_cmd_log = [] - Power_log = [] - Power_target_log = [] + N_log, Wf_act_log, Wf_cmd_log, Power_log, Power_target_log = [], [], [], [], [] + H_env_log, Ma_env_log = [], [] - # 定义多段功率指令 (时间[s], 目标功率[kW]) - power_profile = [ - (0.0, 100.0), # 初始稳态 - (3.0, 200.0), # 阶跃上升 - (8.0, 150.0), # 阶跃下降 - (13.0, 250.0), # 阶跃上升至高功率 - (18.0, 100.0), # 快速下降 - (23.0, 180.0), # 再次上升 - ] + print("-> 开始极限工况仿真测试 (大动态指令 + 连续外界干扰)...") - # 斜坡指令测试：从18kW开始以一定速率上升 - ramp_start_time = 25.0 - ramp_rate = 10.0 # kW/s - - # 设定飞行条件改变的时间 - flight_change_time = 28.0 - - def get_target_power(t, profile, ramp_start, ramp_rate, default_power): - """根据时间获取当前目标功率""" - for i, (time, _) in enumerate(profile): - if t < time: - return profile[i-1][1] if i > 0 else default_power - # 如果在斜坡区间 - if t >= ramp_start: - last_static_power = profile[-1][1] - ramp_power = last_static_power + ramp_rate * (t - ramp_start) - return min(ramp_power, 300.0) # 限制最大300kW - return profile[-1][1] - - print("-> 开始仿真步进 (闭环控制)...") - print("-> 功率指令配置文件:") - for time, power in power_profile: - print(f" t={time:.1f}s -> {power:.0f} kW") - print(f" t={ramp_start_time:.1f}s -> 斜坡上升 (速率 {ramp_rate} kW/s)") - print() - - for t in time_array: - # 确定当前的目标功率 - current_target_P = get_target_power(t, power_profile, ramp_start_time, ramp_rate, P_initial_target) + for t in tqdm(time_array, desc="Simulating"): + # --- 1. 生成复杂功率指令 (Setpoint) --- + if t < 10.0: + # 阶跃测试 + target_p = 100.0 if t < 3 else 220.0 if t < 7 else 80.0 + elif t < 25.0: + # 正弦跟踪 (0.2Hz) + target_p = 140.0 + 50.0 * np.sin(2 * np.pi * 0.2 * (t - 10.0)) + elif t < 40.0: + # 锯齿波测试 + cycle = (t - 25.0) % 5.0 + target_p = 80.0 + (80.0 / 5.0) * cycle + else: + # 极限大范围跳变 + target_p = 170.0 if t < 48.0 else 40.0 - # 触发飞行条件变化 (确保仅触发一次以免重复打印) - if flight_change_time <= t < flight_change_time + dt: - print(f"\n[!] 时间 t={t:.2f}s, 模拟飞行条件跃变...") - sim.set_flight_condition(H_env=3000.0, Ma_env=0.2) + # --- 2. 生成环境扰动 (Disturbance) --- + if t < 15.0: + h_env, ma_env = 0.0, 0.0 + elif t < 30.0: + # 爬升阶段: 0->3000m + h_env = 0.0 + (3000.0 / 15.0) * (t - 15.0) + ma_env = 0.0 + (0.2 / 15.0) * (t - 15.0) + elif t < 45.0: + h_env, ma_env = 3000.0, 0.2 + else: + # 突发机动 + h_env, ma_env = 500.0, 0.4 + + sim.H_env = h_env + sim.Ma_env = ma_env - # 执行闭环仿真步进 - N_cur, Wf_act_cur, Wf_req, Power_cur = sim.step(dt, target_power=current_target_P) + # 执行单步仿真 + N_cur, Wf_act_cur, Wf_req, Power_cur = sim.step(dt, target_power=target_p) - # 记录数据 N_log.append(N_cur) Wf_act_log.append(Wf_act_cur) Wf_cmd_log.append(sim.Wf_cmd) Power_log.append(Power_cur) - Power_target_log.append(current_target_P) + Power_target_log.append(target_p) + H_env_log.append(h_env) + Ma_env_log.append(ma_env) - # 绘图 - fig, (ax1, ax2, ax3) = plt.subplots(3, 1, figsize=(12, 10)) - fig.suptitle(f'Turboshaft Engine Dynamic Simulation\nPID Parameters: Kp={sim.pid.kp}, Ki={sim.pid.ki}, Kd={sim.pid.kd}', - fontsize=12, fontweight='bold') + # ========================================== + # 绘图逻辑 + # ========================================== + fig, axes = plt.subplots(4, 1, figsize=(14, 12), gridspec_kw={'height_ratios': [2.5, 2, 2, 1.5]}) + fig.suptitle('Turboshaft Engine: MPC Extreme Stress Test\nSine Tracking + Continuous Disturbances', + fontsize=14, fontweight='bold') - # Plot 1: Power Tracking - ax1.plot(time_array, Power_target_log, 'k--', linewidth=2, label='Target Power') - ax1.plot(time_array, Power_log, 'g-', linewidth=2, label='Actual Power') - ax1.set_ylabel('Shaft Power [kW]') - ax1.set_title('Closed-Loop Control: Power Tracking', fontweight='bold') - ax1.grid(True, linestyle=':', alpha=0.7) - ax1.legend() + # Plot 1: Power Tracking (核心表现) + axes[0].plot(time_array, Power_target_log, 'k--', linewidth=2, label='Target Power (Command)') + axes[0].plot(time_array, Power_log, 'g-', linewidth=2, label='Actual Power (MPC)') + axes[0].set_ylabel('Shaft Power [kW]', fontweight='bold') + axes[0].set_title('Performance: Complex Trajectory Tracking', fontweight='bold') + axes[0].grid(True, linestyle=':', alpha=0.7) + axes[0].legend(loc='upper right') # Plot 2: Rotor Speed - ax2.plot(time_array, N_log, 'b-', linewidth=2, label='Engine Speed (N)') - ax2.set_ylabel('Rotor Speed [RPM]') - ax2.set_title('Engine Response: Rotor Speed', fontweight='bold') - ax2.grid(True, linestyle=':', alpha=0.7) - ax2.legend() + axes[1].plot(time_array, N_log, 'b-', linewidth=2, label='Engine Speed (N)') + axes[1].set_ylabel('Rotor Speed [RPM]', fontweight='bold') + axes[1].set_title('State: Rotor Speed Response', fontweight='bold') + axes[1].grid(True, linestyle=':', alpha=0.7) + axes[1].legend(loc='upper right') # Plot 3: Fuel Flow (Control Input) - ax3.plot(time_array, Wf_cmd_log, 'k--', linewidth=1.5, label='Fuel Command') - ax3.plot(time_array, Wf_act_log, 'r-', linewidth=2, label='Actual Fuel') - ax3.set_xlabel('Time [s]') - ax3.set_ylabel('Fuel Flow [kg/h]') - ax3.set_title('Control Effort: Fuel Flow', fontweight='bold') - ax3.grid(True, linestyle=':', alpha=0.7) - ax3.legend() + axes[2].plot(time_array, Wf_cmd_log, 'r--', linewidth=1.5, label='Fuel Command (MPC Output)') + axes[2].plot(time_array, Wf_act_log, 'm-', linewidth=2, label='Actual Fuel Actuator') + axes[2].set_ylabel('Fuel Flow [kg/h]', fontweight='bold') + axes[2].set_title('Control Effort: Actuator Dynamics', fontweight='bold') + axes[2].grid(True, linestyle=':', alpha=0.7) + axes[2].legend(loc='upper right') + + # Plot 4: Environmental Disturbances + ax4_1 = axes[3] + ax4_2 = ax4_1.twinx() + ax4_1.plot(time_array, H_env_log, 'c-', linewidth=2, label='Altitude (m)') + ax4_2.plot(time_array, Ma_env_log, 'y-', linewidth=2, label='Mach Number') + ax4_1.set_xlabel('Time [s]', fontweight='bold') + ax4_1.set_ylabel('Altitude [m]', color='c', fontweight='bold') + ax4_2.set_ylabel('Mach', color='y', fontweight='bold') + axes[3].set_title('Disturbances: Flight Conditions', fontweight='bold') + axes[3].grid(True, linestyle=':', alpha=0.7) plt.tight_layout(rect=[0, 0, 1, 0.96]) - plt.savefig('figures/engine_dynamic_sim_plot.png') - print('Plot saved to figures/engine_dynamic_sim_plot.png') - - # 保存数据到 .dat 文件 - dat_file = 'data/pid_tuning_data.dat' - with open(dat_file, 'w') as f: - f.write('# Time[s]\tTarget_Power[kW]\tActual_Power[kW]\tRotor_Speed[RPM]\tFuel_Command[kg/h]\tActual_Fuel[kg/h]\n') - for i in range(len(time_array)): - f.write(f'{time_array[i]:.4f}\t{Power_target_log[i]:.2f}\t{Power_log[i]:.2f}\t{N_log[i]:.2f}\t{Wf_cmd_log[i]:.2f}\t{Wf_act_log[i]:.2f}\n') - print(f'Data saved to {dat_file}') - # plt.show() + plt.savefig('figures/engine_mpc_stress_test.png', dpi=200) + print('-> 仿真完成！极限制图已保存至: figures/engine_mpc_stress_test.png') diff --git a/Model/src/engine_gpr_class.py b/Model/src/engine_gpr_class.py index ab18054..b8e52dc 100644 --- a/Model/src/engine_gpr_class.py +++ b/Model/src/engine_gpr_class.py @@ -34,10 +34,13 @@ class EngineGPRModel: self.scaler_X = StandardScaler() self.scaler_Y = StandardScaler() self.model = None + # 自动检测并使用 GPU + self.device = torch.device("cuda" if torch.cuda.is_available() else "cpu") + print(f"Using device: {self.device}") def _prepare_data(self): """数据预处理：读取、筛选、Log变换、标准化""" - self.df = pd.read_csv(self.csv_path) + self.df = pd.read_csv(self.csv_path, encoding='utf-8') X_df = self.df[['Altitude_m', 'Mach', 'RPM']] # 剔除无效特征 @@ -56,7 +59,8 @@ class EngineGPRModel: X_scaled = self.scaler_X.fit_transform(X_numpy) Y_scaled = self.scaler_Y.fit_transform(Y_numpy) - return torch.tensor(X_scaled, dtype=torch.double), torch.tensor(Y_scaled, dtype=torch.double) + # 转换为 Tensor 并移动到 GPU (如果可用) + return torch.tensor(X_scaled, dtype=torch.double).to(self.device), torch.tensor(Y_scaled, dtype=torch.double).to(self.device) def _init_model(self, train_X, train_Y): """内部方法：统一初始化模型结构（包含 Mean 和 Prior 设置）""" @@ -83,6 +87,7 @@ class EngineGPRModel: # 使用统一初始化方法 self.model = self._init_model(train_X, train_Y) + self.model.to(self.device) # 仅在训练开始前设定初始值，引导优化方向 if hasattr(self.model.covar_module, 'base_kernel'): @@ -109,10 +114,12 @@ class EngineGPRModel: train_X, train_Y = self._prepare_data() # 必须使用完全相同的结构初始化，否则 load_state_dict 会报错 self.model = self._init_model(train_X, train_Y) + self.model.to(self.device) # 使用 strict=False 忽略 Prior 缓冲区的差异（例如 _transformed_loc 等内部参数） # 这些参数通常不影响模型预测，只影响后续继续训练时的约束 - self.model.load_state_dict(torch.load(pth_path), strict=True) + state_dict = torch.load(pth_path, map_location=self.device) + self.model.load_state_dict(state_dict, strict=True) self.model.eval() print("-> Model loaded successfully.") return True @@ -125,12 +132,19 @@ class EngineGPRModel: if self.model is None: raise ValueError("Model not initialized.") self.model.eval() - test_X_scaled = torch.tensor(self.scaler_X.transform(test_X_real), dtype=torch.double) - + + # Determine device from model parameters + try: + device = next(self.model.parameters()).device + except StopIteration: + device = torch.device('cpu') + + test_X_scaled = torch.tensor(self.scaler_X.transform(test_X_real), dtype=torch.double, device=device) + with torch.no_grad(): posterior = self.model.posterior(test_X_scaled) - mu_scaled = posterior.mean.numpy() - var_scaled = posterior.variance.numpy() + mu_scaled = posterior.mean.detach().cpu().numpy() + var_scaled = posterior.variance.detach().cpu().numpy() # 反归一化 mu_log_real = self.scaler_Y.inverse_transform(mu_scaled) diff --git a/Model/src/lightweight_model.py b/Model/src/lightweight_model.py new file mode 100644 index 0000000..84cd510 --- /dev/null +++ b/Model/src/lightweight_model.py @@ -0,0 +1,70 @@ +import torch +import torch.nn as nn +import numpy as np + +class EngineNNProxy(nn.Module): + """ + 轻量级神经网络代理模型，用于替代笨重的 GPR 模型。 + 结构: 简单的 MLP (多层感知机) + 输入: [Altitude, Mach, RPM] (未归一化) + 输出: [FuelFlow, Power] (未归一化) + """ + def __init__(self, hidden_size=64): + super(EngineNNProxy, self).__init__() + + # 定义网络结构 + # 对应 distill_gpr_to_nn.py 中的索引访问: + # 0: Linear + # 1: ReLU / Tanh + # 2: Linear + # 3: ReLU / Tanh + # 4: Linear (Output) + self.net = nn.Sequential( + nn.Linear(3, hidden_size), # 0 + nn.Tanh(), # 1: Tanh 通常比 ReLU 更适合平滑的物理函数拟合 + nn.Linear(hidden_size, hidden_size), # 2 + nn.Tanh(), # 3 + nn.Linear(hidden_size, 2) # 4: 输出 2 个物理量 (Fuel, Power) + ) + + # 归一化参数 (注册为 buffer 以便随模型保存) + # 初始化为默认值，防止未调用 set_normalization_params 时报错 + self.register_buffer('x_mean', torch.zeros(3)) + self.register_buffer('x_scale', torch.ones(3)) + self.register_buffer('y_mean', torch.zeros(2)) + self.register_buffer('y_scale', torch.ones(2)) + + def set_normalization_params(self, scaler_X, scaler_Y): + """ + 从 sklearn StandardScaler 中提取参数 + scaler_X: 用于输入的归一化器 + scaler_Y: 用于输出 (Log1p Space) 的归一化器 + """ + if scaler_X is not None: + self.x_mean.copy_(torch.tensor(scaler_X.mean_, dtype=torch.float32)) + self.x_scale.copy_(torch.tensor(scaler_X.scale_, dtype=torch.float32)) + + if scaler_Y is not None: + self.y_mean.copy_(torch.tensor(scaler_Y.mean_, dtype=torch.float32)) + self.y_scale.copy_(torch.tensor(scaler_Y.scale_, dtype=torch.float32)) + + def forward(self, x): + """ + 前向传播: 物理输入 -> 物理输出 + 包含: 归一化 -> NN推理 -> 反归一化 -> expm1 + """ + # 1. 输入归一化 (Z-Score) + # 确保输入 x 与 buffer 在同一设备 + x = x.to(self.x_mean.device) + x_norm = (x - self.x_mean) / self.x_scale + + # 2. 神经网络推理 (预测 Log Normalized Z-Score) + y_norm_pred = self.net(x_norm) + + # 3. 输出反归一化 (Z-Score Inverse) + y_log1p_pred = y_norm_pred * self.y_scale + self.y_mean + + # 4. 指数还原 (Inverse Log1p) + y_phys_pred = torch.expm1(y_log1p_pred) + + return y_phys_pred diff --git a/Model/src/motor_sim.py b/Model/src/motor_sim.py index 122716e..9fdfaa1 100644 --- a/Model/src/motor_sim.py +++ b/Model/src/motor_sim.py @@ -9,18 +9,16 @@ sys.path.insert(0, os.path.dirname(os.path.dirname(os.path.abspath(__file__)))) 模块提供了永磁同步电机 (PMSM) 的离散时间仿真实现，涵盖了闭环转速控制、 物理级的电磁转矩估算及端电压/损耗模型的动力学计算。 -系统符号学约定: +系统符号约定: - 轴系转矩 (T_motor): 正值表示吸收轴系功率 (发电机/负载响应)，负值表示向轴系输出功率 (驱动动力)。 - 直流母线功率 (P_bus): 正值表示从母线汲取有功功率，负值表示向母线回馈有功功率。 """ -from typing import Any, Dict -import sys -import os -sys.path.insert(0, os.path.dirname(os.path.dirname(os.path.abspath(__file__)))) +from typing import Any, Dict, Optional import numpy as np -from src.increPID import IncrementalPIDController +from tqdm import tqdm +from src.mpc_controller import MotorMPCController class MotorSim: @@ -40,12 +38,7 @@ class MotorSim: eta_mot: float = 0.95, eta_gen: float = 0.93, B_visc: float = 3e-4, - k_p_w: float = 13.440362, - k_i_w: float = 42.997816, - k_d_w: float = 0.484660, - k_d_w_error: float = 0.05, p_bus_slew_rate_kw: float = 2000.0, - t_cmd_slew_rate_nm_s: float = 0.0, tau_i: float = 0.004, k_mod: float = 1 / np.sqrt(3), P_const_loss: float = 250.0, @@ -53,9 +46,11 @@ class MotorSim: k_inv: float = 0.015, tau_n_ref: float = 0.35, tau_t_cmd: float = 0.12, - speed_priority_band_rpm: float = 180.0, tau_v_bus: float = 0.08, - i_s_max: float | None = None, + i_s_max: Optional[float] = None, + mpc_W_speed: float = 0.0, + mpc_W_dcost: float = 0.0, + mpc_overshoot_limit: float = 0.05, ): """ 初始化电机动力学仿真器参数。 @@ -72,12 +67,7 @@ class MotorSim: :param eta_mot: 稳态驱动效率参考值 :param eta_gen: 稳态发电效率参考值 :param B_visc: 轴系黏性摩擦阻尼系数 (N·m·s) - :param k_p_w: 速度控制环比例增益 - :param k_i_w: 速度控制环积分增益 - :param k_d_w: 速度控制环微分控制增益(闭环阻尼) - :param k_d_w_error: 目标速度微分前馈增益 :param p_bus_slew_rate_kw: 允许的最大母线请求功率变化率 (kW/s) - :param t_cmd_slew_rate_nm_s: 允许的最大转矩指令变化率 (Nm/s) :param tau_i: 闭环电流/电磁转矩等效一阶延迟时间常数 (s) :param k_mod: DC-AC 变换器电压利用系数 :param P_const_loss: 独立于工况的系统常量损耗 (W) @@ -85,9 +75,11 @@ class MotorSim: :param k_inv: DC-AC 逆变环节损耗系数 :param tau_n_ref: 目标转速给定指令一阶低通滤波时间常数 (s) :param tau_t_cmd: 转矩输出指令一阶低通滤波时间常数 (s) - :param speed_priority_band_rpm: 目标跟随允许容差基准带 (RPM) :param tau_v_bus: DC侧动态响应一阶滤波时间常数 (s) :param i_s_max: 定子相电流约束阈值有效值 (A) + :param mpc_W_speed: MPC转速控制权重 + :param mpc_W_dcost: MPC控制增量权重 + :param mpc_overshoot_limit: MPC超调量硬约束 (0.05 = 5%) """ self._validate_parameters( n_p, @@ -116,12 +108,7 @@ class MotorSim: self.eta_mot = eta_mot self.eta_gen = eta_gen self.B_visc = B_visc - self.k_p_w = k_p_w - self.k_i_w = k_i_w - self.k_d_w = k_d_w - self.k_d_w_error = k_d_w_error self.p_bus_slew_rate_kw = p_bus_slew_rate_kw - self.t_cmd_slew_rate_nm_s = max(0.0, float(t_cmd_slew_rate_nm_s)) self.tau_i = tau_i self.k_mod = k_mod self.P_const_loss = P_const_loss @@ -129,9 +116,11 @@ class MotorSim: self.k_inv = k_inv self.tau_n_ref = max(1e-4, float(tau_n_ref)) self.tau_t_cmd = max(1e-4, float(tau_t_cmd)) - self.speed_priority_band = max(5.0, float(speed_priority_band_rpm)) * 2.0 * np.pi / 60.0 self.tau_v_bus = max(1e-4, float(tau_v_bus)) + # 超调量约束参数 (用于转矩变化率限制) + self._overshoot_limit = mpc_overshoot_limit + self.tau_rate = self.P_rate / self.w_rate if i_s_max is None: base_current = self.tau_rate / (1.5 * self.n_p * max(abs(self.psi_f), 1e-6)) @@ -160,28 +149,20 @@ class MotorSim: self.u_an, self.u_bn, self.u_cn = 0.0, 0.0, 0.0 self.duty_a, self.duty_b, self.duty_c = 0.5, 0.5, 0.5 self._sim_t = 0.0 - self._w_error_int = 0.0 - self._w_error_prev = 0.0 - self._dw_error_f = 0.0 - self._p_bus_req_prev = 0.0 self._w_set_f = 0.0 self._t_ref_f = 0.0 self._v_bus_f = float(self.u_dc) self._t_load_f = 0.0 self._t_ext_f = 0.0 - # 系统机械约束约束配置 + # 系统机械约束配置 self._dw_max = 800.0 - # 增量式PID控制器 (用于转速环功率控制) - # output_scale: 额定功率300kW，input_scale: 额定转速575.95 rad/s - self.speed_pid = IncrementalPIDController( - kp=k_p_w, ki=k_i_w, kd=k_d_w, - dt=0.02, # 将在step中动态更新 - output_min=-self.P_rate/1000.0 * 1.2, # 允许放电 - output_max=self.P_rate/1000.0 * 1.2, # 允许充电 - input_scale=self.w_rate, - output_scale=self.P_rate/1000.0 + # MPC控制器 (用于转速环功率控制) + self.mpc = MotorMPCController( + J=J, B_visc=B_visc, dt=0.02, horizon=10, + W_speed=mpc_W_speed, W_dcost=mpc_W_dcost, + overshoot_limit=mpc_overshoot_limit ) self._log = { @@ -208,7 +189,7 @@ class MotorSim: } # 滤波系数缓存 (用于优化: 避免每步重复计算 exp(-dt/tau)) - self._last_dt: float | None = None + self._last_dt: Optional[float] = None self._a_v: float = 0.0 # 母线电压滤波系数 self._a_n: float = 0.0 # 转速参考滤波系数 self._a_t_load: float = 0.0 # 负载转矩滤波系数 @@ -220,8 +201,6 @@ class MotorSim: self._EPSILON: float = 1e-6 # 数值 epsilon，用于避免除零 self._TAU_LOAD_FILTER: float = 0.05 # 负载转矩滤波时间常数 (s) self._TORQUE_OVERRATE: float = 1.35 # 转矩过载系数 - self._W_ERROR_DEADBAND: float = 0.05 # 转速误差死区 (rad/s) - self._DW_ERROR_FILTER: float = 0.15 # 微分误差滤波系数 # 日志记录配置：可选的最大日志长度限制 (None 表示无限制) self._log_maxlen: int | None = None # 可设置为如 100000 来限制内存使用 @@ -380,10 +359,11 @@ class MotorSim: w_abs = abs(w_M) w_e = self.n_p * w_abs - if w_e < self._EPSILON: + if w_e < 1e-6: return float(max(0.0, abs(self._torque_from_currents(0.0, self.i_s_max)))) - # 优化: 减少网格点数从40到12，使用更智能的搜索策略 + # 优化: 减少网格点数从40到12，在保证精度的同时提高效率 + # 使用更智能的搜索策略：从估计的MTPA点开始搜索 i_d_estimate = max(-self.i_s_max, min(0.0, -0.5 * (self.psi_f / (self.L_q - self.L_d)) if self.L_q != self.L_d else 0.0)) i_d_grid = np.array([ i_d_estimate - 0.3 * self.i_s_max, @@ -407,20 +387,14 @@ class MotorSim: tau_candidate = abs(self._torque_from_currents(i_d, i_q)) if tau_candidate > tau_best: tau_best = tau_candidate - if i_d < -0.01 * self.i_s_max: + if i_d < -0.01 * self.i_s_max: # 继续搜索更负的i_d continue tau_p_max = self.P_rate / max(w_abs, self._W_MIN_RAD) return float(max(0.0, min(tau_best, tau_p_max, self._TORQUE_OVERRATE * self.tau_rate))) - def _apply_power_slew(self, p_bus_req_kw: float, dt: float) -> float: - """对母线请求功率施加斜率限制。""" - if self.p_bus_slew_rate_kw <= 0: - self._p_bus_req_prev = float(p_bus_req_kw) - return float(p_bus_req_kw) - - delta = self.p_bus_slew_rate_kw * dt - self._p_bus_req_prev = float(p_bus_req_kw) + def _apply_power_slew(self, p_bus_req_kw: float) -> float: + """对母线请求功率施加斜率限制（当前实现为 passthrough）。""" return float(p_bus_req_kw) def _power_to_torque_ref(self, p_bus_w: float, w_M: float) -> float: @@ -434,13 +408,13 @@ class MotorSim: if abs(p_bus_w) < 1e-9: return 0.0 - w_eff = max(abs(w_M), 50.0) + w_eff = max(abs(w_M), self._W_MIN_RAD) if p_bus_w >= 0: # 取电模式：从母线取有功，经过内部热量/电磁损耗后，输出转动推力(负转矩) return -(p_bus_w * self.eta_mot) / w_eff else: # 发电模式：电机受到负负载强推，需要输出极大的正机械阻力(查表吸收电能) - return -p_bus_w / (w_eff * max(self.eta_gen, 1e-6)) + return -p_bus_w / (w_eff * max(self.eta_gen, self._EPSILON)) def _current_ref_from_torque(self, t_motor_ref: float, w_M: float, v_bus: float) -> tuple[float, float]: """ @@ -468,12 +442,14 @@ class MotorSim: iq_sign = 1.0 if i_q_ref >= 0 else -1.0 iq_abs_target = abs(i_q_ref) - # 优化: 减少网格点数从80到20 - for i_d in np.linspace(-self.i_s_max, 0.0, 20): + # 优化: 减少网格点数从80到20，使用更高效的搜索 + i_d_candidates = np.linspace(-self.i_s_max, 0.0, 20) + for i_d in i_d_candidates: iq_max = self._max_iq_given_id(i_d, w_e, v_bus) if iq_max >= iq_abs_target: return float(i_d), float(iq_sign * iq_abs_target) + # 如果没找到，使用最大弱磁点 i_d_fw = -self.i_s_max i_q_fw = iq_sign * self._max_iq_given_id(i_d_fw, w_e, v_bus) return float(i_d_fw), float(i_q_fw) @@ -560,10 +536,10 @@ class MotorSim: """ if dt <= 0: raise ValueError('dt must be positive') - if v_bus <= self._EPSILON: + if v_bus <= 1e-6: raise ValueError('v_bus must be positive') - # 更新滤波系数缓存 + # 更新滤波系数缓存 (仅当 dt 变化时重新计算) self._update_filter_coefficients(dt) self._v_bus_f += (float(v_bus) - self._v_bus_f) * self._a_v @@ -576,51 +552,29 @@ class MotorSim: self._w_set_f += (w_set_cmd - self._w_set_f) * self._a_n w_set = self._w_set_f - w_error = w_set - self.w_M - if abs(w_error) < self._W_ERROR_DEADBAND: - w_error = 0.0 - - dw_error = (w_error - self._w_error_prev) / max(dt, self._EPSILON) - self._w_error_prev = w_error - self._dw_error_f = self._dw_error_f * (1 - self._DW_ERROR_FILTER) + dw_error * self._DW_ERROR_FILTER - - self._w_error_int += w_error * dt p_lim_kw = (t_lim * max(self.w_M, self._W_MIN_RAD)) / 1000.0 / max(self.eta_mot, self._EPSILON) - int_limit = 1.2 * p_lim_kw / max(self.k_i_w, self._EPSILON) if self.k_i_w > 0 else 0.0 - self._w_error_int = np.clip(self._w_error_int, -int_limit, int_limit) self._t_load_f += (float(t_load) - self._t_load_f) * self._a_t_load self._t_ext_f += (float(t_ext) - self._t_ext_f) * self._a_t_load - # 阻力前馈折算为补偿所需电功率(kW) - t_ff = self._t_load_f + self._t_ext_f + self.B_visc * self.w_M - p_mech_ff_kw = (t_ff * max(self.w_M, 1e-3)) / 1000.0 - if p_mech_ff_kw > 0: - p_elec_ff_kw = p_mech_ff_kw / max(self.eta_mot, self._EPSILON) - else: - p_elec_ff_kw = p_mech_ff_kw * self.eta_gen - - # 阻力前馈折算为补偿所需电功率(kW) - t_ff = self._t_load_f + self._t_ext_f + self.B_visc * self.w_M - p_mech_ff_kw = (t_ff * max(self.w_M, 1e-3)) / 1000.0 - if p_mech_ff_kw > 0: - p_elec_ff_kw = p_mech_ff_kw / max(self.eta_mot, 1e-6) - else: - p_elec_ff_kw = p_mech_ff_kw * self.eta_gen - # ------------------------------------------------------------- - # A) 需求生成层: 以转速环计算下一拍应向电池索取的 P_bus_req_kw + # A) 需求生成层: 以转速环 MPC 计算下一拍应向电池索取的 P_bus_req_kw # ------------------------------------------------------------- - # 使用增量式PID计算控制量 - self.speed_pid.dt = dt - p_pid = self.speed_pid.compute(setpoint=w_set, measurement=self.w_M) + # 使用 MPC 控制器计算转矩指令，然后转换为功率请求 + t_motor_cmd = self.mpc.compute( + current_w=self.w_M, + target_w=w_set, + t_load=self._t_load_f, + t_ext=self._t_ext_f, + t_lim_upper=t_lim, + t_lim_lower=-t_lim + ) - # 添加微分前馈和阻力前馈 - p_d_ff = self.k_d_w_error * self._dw_error_f - - p_cmd_raw = p_pid + p_d_ff + p_elec_ff_kw + # 将 MPC 转矩指令转换为功率请求 (kW) + # 正转矩(发电) -> 吸收功率，负转矩(驱动) -> 输出功率 + p_cmd_raw = -t_motor_cmd * max(self.w_M, self._W_MIN_RAD) / 1000.0 p_bus_req_kw = float(np.clip(p_cmd_raw, -p_lim_kw, p_lim_kw)) - p_bus_req_kw = self._apply_power_slew(p_bus_req_kw, dt) + p_bus_req_kw = self._apply_power_slew(p_bus_req_kw) self.p_bus_req_kw = p_bus_req_kw @@ -630,6 +584,47 @@ class MotorSim: t_motor_ref = self._power_to_torque_ref(p_bus_actual_kw * 1000.0, self.w_M) t_motor_ref = float(np.clip(t_motor_ref, -t_lim, t_lim)) + # ------------------------------------------------------------- + # 超调量硬约束: 基于转矩变化率限制 + # 当接近目标转速时，限制转矩变化以防止超调 + # ------------------------------------------------------------- + w_target = w_set + w_max_allowed = w_target * (1 + self._overshoot_limit) # 105% + w_min_allowed = w_target * (1 - self._overshoot_limit) # 95% + + # 记录上一时刻的转矩参考 + t_motor_ref_prev = getattr(self, '_t_motor_ref_prev', 0.0) + + if w_target > self._W_MIN_RAD: + # 升速时: 如果当前转速超过允许最大值，强制减速 + if w_set > self.w_M and self.w_M > w_max_allowed: + t_motor_ref = -t_lim # 最大制动 + + # 降速时: 如果当前转速低于允许最小值，强制加速 + elif w_set < self.w_M and self.w_M < w_min_allowed: + t_motor_ref = t_lim # 最大驱动 + + # 接近目标时 (95%-105% 区间): 限制转矩变化率 + elif self.w_M > w_target * 0.90: + # 计算允许的最大转矩变化 + # 使转速不超过 w_max_allowed 的最大加速度 + max_allowed_dw = (w_max_allowed - self.w_M) / dt + max_allowed_T = -self.J * max_allowed_dw - t_load - t_ext - self.B_visc * self.w_M + + # 限制正转矩 (驱动) 不超过计算值 + if t_motor_ref < 0: # 驱动扭矩 + t_motor_ref = max(t_motor_ref, max_allowed_T) + + # 使转速不低于 w_min_allowed 的最小加速度 + min_allowed_dw = (w_min_allowed - self.w_M) / dt + min_allowed_T = -self.J * min_allowed_dw - t_load - t_ext - self.B_visc * self.w_M + + # 限制负转矩 (制动) 不超过计算值 + if t_motor_ref > 0: # 制动扭矩 + t_motor_ref = min(t_motor_ref, min_allowed_T) + + self._t_motor_ref_prev = t_motor_ref + self._t_ref_f += (t_motor_ref - self._t_ref_f) * self._a_t t_motor_ref_filtered = float(self._t_ref_f) @@ -760,12 +755,8 @@ class MotorSim: self.u_an, self.u_bn, self.u_cn = 0.0, 0.0, 0.0 self.duty_a, self.duty_b, self.duty_c = 0.5, 0.5, 0.5 self._sim_t = 0.0 - # 增量式PID控制器重置 - self.speed_pid.reset(initial_output=0.0) - self._w_error_int = 0.0 - self._w_error_prev = 0.0 - self._dw_error_f = 0.0 - self._p_bus_req_prev = 0.0 + # MPC控制器重置 + self.mpc.reset(initial_w=self.w_M) self._w_set_f = 0.0 self._t_ref_f = 0.0 self._v_bus_f = float(self.u_dc) @@ -774,231 +765,6 @@ class MotorSim: self._t_ext_f = 0.0 -def run_motor_test() -> Dict[str, Any]: - """基于接口的独立电机测试 - 复杂测试用例""" - import matplotlib.pyplot as plt - - print('=' * 60) - print('Motor Dynamic Test - Complex Profile') - print('=' * 60) - - # 使用优化后的PID参数 - motor = MotorSim( - P_rate=300e3, - w_rate=575.95, - J=0.8, - k_p_w=13.440362, - k_i_w=42.997816, - k_d_w=0.484660, - P_const_loss=300.0, - k_fe=4e-4, - ) - - dt = 0.02 - t_end = 100.0 - time_array = np.arange(0.0, t_end, dt) - - # 存储关键诊断数据 - torque_balance = [] - power_balance = [] - - p_actual_kw = 0.0 - for t in time_array: - # 复杂测试工况 - 多段转速变化 + 多种扰动 - # 目标转速曲线 - if t < 3.0: - n_setpoint = 0.0 - elif t < 8.0: - n_setpoint = 2000.0 - elif t < 15.0: - n_setpoint = 2800.0 - elif t < 22.0: - n_setpoint = 2200.0 - elif t < 32.0: - n_setpoint = 3500.0 - elif t < 42.0: - n_setpoint = 3800.0 - elif t < 52.0: - n_setpoint = 2800.0 - elif t < 62.0: - n_setpoint = 3000.0 - elif t < 72.0: - n_setpoint = 4000.0 - elif t < 82.0: - n_setpoint = 4200.0 - elif t < 92.0: - n_setpoint = 3200.0 - else: - n_setpoint = 3000.0 - - # 外部转矩扰动 - if t < 5.0: - t_ext = 25.0 - elif t < 12.0: - t_ext = 10.0 - elif t < 18.0: - t_ext = -30.0 - elif t < 25.0: - t_ext = 15.0 - elif t < 35.0: - t_ext = -50.0 - elif t < 45.0: - t_ext = 8.0 - elif t < 55.0: - t_ext = -20.0 - elif t < 65.0: - t_ext = 25.0 - elif t < 75.0: - t_ext = -60.0 - elif t < 85.0: - t_ext = 20.0 - else: - t_ext = 0.0 - - # 负载转矩 - 变化因子 - if t < 10.0: - t_load_factor = 1.0 - elif t < 22.0: - t_load_factor = 1.5 - elif t < 32.0: - t_load_factor = 0.8 - elif t < 42.0: - t_load_factor = 1.8 - elif t < 52.0: - t_load_factor = 1.0 - elif t < 62.0: - t_load_factor = 1.4 - elif t < 72.0: - t_load_factor = 0.9 - elif t < 82.0: - t_load_factor = 1.6 - else: - t_load_factor = 1.1 - - base_t_load = 16.0 + 0.020 * motor.w_M + 1.2e-5 * motor.w_M**2 - t_load = base_t_load * t_load_factor - - # 母线电压波动 - v_bus = 520.0 + 15.0 * np.sin(2 * np.pi * t / 8.0) + 5.0 * np.sin(2 * np.pi * t / 3.0) - - state = motor.step(dt, n_setpoint, p_actual_kw, v_bus, t_load, t_ext) - p_actual_kw = state['p_bus_req_kw'] - - torque_balance.append(state['torque_balance_error']) - power_balance.append(state['power_balance_error']) - - state = motor.get_state() - p_bus_kw = np.array(motor._log['p_bus_kw']) - p_req_kw = np.array(motor._log['p_bus_req_kw']) - p_shaft_kw = np.array(motor._log['p_shaft_kw']) - p_loss_kw = np.array(motor._log['p_loss_kw']) - n_log = np.array(motor._log['n_rpm']) - t_motor_log = np.array(motor._log['t_motor']) - t_load_log = np.array(motor._log['t_load']) - t_ext_log = np.array(motor._log['t_ext']) - v_bus_log = np.array(motor._log['v_bus']) - v_motor_log = np.array(motor._log['v_motor']) - i_bus_log = np.array(motor._log['i_bus_a']) - - # 计算性能指标 - n_setpoints = [] - for t in time_array: - if t < 3.0: - n_setpoints.append(0.0) - elif t < 8.0: - n_setpoints.append(2000.0) - elif t < 15.0: - n_setpoints.append(2800.0) - elif t < 22.0: - n_setpoints.append(2200.0) - elif t < 32.0: - n_setpoints.append(3500.0) - elif t < 42.0: - n_setpoints.append(3800.0) - elif t < 52.0: - n_setpoints.append(2800.0) - elif t < 62.0: - n_setpoints.append(3000.0) - elif t < 72.0: - n_setpoints.append(4000.0) - elif t < 82.0: - n_setpoints.append(4200.0) - elif t < 92.0: - n_setpoints.append(3200.0) - else: - n_setpoints.append(3000.0) - n_setpoints = np.array(n_setpoints) - - error = n_setpoints - n_log - ise = np.sum(error**2) * dt - iae = np.sum(np.abs(error)) * dt - - print(f"Final: speed={state['n_rpm']:.0f} RPM ({state['w_rad_s']:.1f} rad/s)") - print(f" t_motor={state['t_motor']:.2f} Nm, p_bus={state['p_bus_kw']:.2f} kW, v_motor={state['v_motor']:.1f} V") - print(f"Performance - ISE: {ise:.2f}, IAE: {iae:.2f}") - - # 诊断结果 - torque_balance = np.array(torque_balance) - power_balance = np.array(power_balance) - print(f"Diagnostic - Torque balance error (mean): {np.mean(np.abs(torque_balance)):.4f} Nm") - print(f"Diagnostic - Power balance error (mean): {np.mean(np.abs(power_balance)):.4f} W") - - fig, axes = plt.subplots(5, 1, figsize=(14, 14), sharex=True) - - axes[0].plot(time_array, n_setpoints, 'r--', lw=1.0, alpha=0.7, label='Setpoint') - axes[0].plot(time_array, n_log, 'b-', lw=1.6, label='Actual') - axes[0].set_ylabel('Speed [RPM]') - axes[0].set_title('Motor Speed - Complex Profile Test (100s)') - axes[0].grid(True, linestyle=':') - axes[0].legend() - - axes[1].plot(time_array, t_motor_log, 'r-', lw=1.6, label='Motor Torque') - axes[1].plot(time_array, t_load_log, 'k--', lw=1.0, label='Load Torque') - axes[1].plot(time_array, t_ext_log, color='tab:purple', lw=1.0, label='External Torque') - axes[1].set_ylabel('Torque [Nm]') - axes[1].set_title('Torque Balance') - axes[1].axhline(0.0, color='k', lw=0.6) - axes[1].grid(True, linestyle=':') - axes[1].legend() - - axes[2].plot(time_array, p_req_kw, color='tab:gray', lw=1.2, label='Requested Bus Power') - axes[2].plot(time_array, p_bus_kw, color='m', lw=1.4, label='Actual Bus Power') - axes[2].plot(time_array, p_shaft_kw, 'g-', lw=1.2, label='Shaft Power') - axes[2].plot(time_array, p_loss_kw, color='tab:orange', lw=1.0, label='Loss Power') - axes[2].set_ylabel('Power [kW]') - axes[2].set_title('Bus / Shaft / Loss Power') - axes[2].axhline(0.0, color='k', lw=0.6) - axes[2].grid(True, linestyle=':') - axes[2].legend() - - axes[3].plot(time_array, v_bus_log, color='tab:cyan', lw=1.4, label='Bus Voltage') - axes[3].plot(time_array, v_motor_log, color='tab:red', lw=1.2, label='Motor Potential') - axes[3].set_ylabel('Voltage [V]') - axes[3].set_title('Bus Voltage and Motor Equivalent Potential') - axes[3].grid(True, linestyle=':') - axes[3].legend() - - axes[4].plot(time_array, i_bus_log, color='tab:brown', lw=1.4, label='Bus Current') - axes[4].set_ylabel('Current [A]') - axes[4].set_xlabel('Time [s]') - axes[4].set_title('Bus Current') - axes[4].grid(True, linestyle=':') - axes[4].legend() - - plt.tight_layout() - plt.savefig("figures/test_complex_motor_output.png") - - return { - 'final_state': state, - 'peak_speed_rpm': float(np.max(n_log)), - 'peak_bus_current_a': float(np.max(i_bus_log)), - 'torque_balance_error_mean': float(np.mean(np.abs(torque_balance))), - 'power_balance_error_mean': float(np.mean(np.abs(power_balance))), - 'ise': ise, - 'iae': iae, - } - - def run_motor_battery_coupled_test() -> Dict[str, Any]: """执行混合动力电机与电池协同闭环抗扰动验证例程 - 复杂测试用例""" import matplotlib @@ -1019,12 +785,11 @@ def run_motor_battery_coupled_test() -> Dict[str, Any]: P_rate=300e3, w_rate=575.95, J=0.8, - k_p_w=13.440362, - k_i_w=42.997816, - k_d_w=0.484660, P_const_loss=300.0, k_fe=4e-4, u_dc=battery.V_t, + mpc_W_speed=100, + mpc_W_dcost=12, ) dt = 0.02 @@ -1041,36 +806,60 @@ def run_motor_battery_coupled_test() -> Dict[str, Any]: t_ext_log = [] t_load_log = [] + # 螺旋桨负载模型参数 + k_drag = 1.5e-6 # 气动阻力系数 (与转速平方成正比) + J_prop = 0.15 # 螺旋桨等效转动惯量 (kg·m²) + T_min = 10.0 # 最小阻力矩 (Nm) + p_actual_kw = 0.0 - for t in time_array: + n_prev = 0.0 # 上一时刻转速 (RPM) + w_prev = 0.0 # 上一时刻角速度 (rad/s) + for t in tqdm(time_array, desc="Coupled Test"): # 复杂测试工况 - 多段转速变化 + 多种扰动 p_ext_elec = 0.0 # 外部电功率扰动 (kW) - # 目标转速曲线 - 复杂多段变化 + # 目标转速曲线 - 包含阶跃、斜坡、正弦波等复杂信号 if t < 3.0: n_setpoint = 0.0 elif t < 8.0: - n_setpoint = 2000.0 + # 斜坡上升：从0到2000 RPM + n_setpoint = 2000.0 * (t - 3.0) / 5.0 elif t < 15.0: + # 阶跃 n_setpoint = 2800.0 elif t < 22.0: - n_setpoint = 2200.0 + # 正弦波动：2800 ± 600 RPM，周期10s + n_setpoint = 2800.0 + 600.0 * np.sin(2 * np.pi * (t - 15.0) / 10.0) elif t < 32.0: - n_setpoint = 3500.0 + # 斜坡上升 + 正弦波动 + base = 2200.0 + 1300.0 * (t - 22.0) / 10.0 # 斜坡上升 + n_setpoint = base + 300.0 * np.sin(2 * np.pi * (t - 22.0) / 8.0) elif t < 42.0: + # 阶跃到3800 n_setpoint = 3800.0 elif t < 52.0: - n_setpoint = 2800.0 + # 斜坡下降 + 正弦波动 + base = 3800.0 - 1000.0 * (t - 42.0) / 10.0 # 斜坡下降 + n_setpoint = base + 200.0 * np.sin(2 * np.pi * (t - 42.0) / 6.0) elif t < 62.0: - n_setpoint = 3000.0 + # 梯形波：先上升再保持 + if t < 55.0: + n_setpoint = 2800.0 + 400.0 * (t - 52.0) / 3.0 + else: + n_setpoint = 3200.0 elif t < 72.0: - n_setpoint = 4000.0 + # 正弦波动：3000 ± 1000 RPM + n_setpoint = 3000.0 + 1000.0 * np.sin(2 * np.pi * (t - 62.0) / 10.0) elif t < 82.0: - n_setpoint = 4200.0 + # 斜坡上升 + n_setpoint = 3000.0 + 600.0 * (t - 72.0) / 10.0 elif t < 92.0: - n_setpoint = 3200.0 + # 衰减正弦波 + amplitude = 400.0 * np.exp(-(t - 82.0) / 5.0) + n_setpoint = 3200.0 + amplitude * np.sin(2 * np.pi * (t - 82.0) / 8.0) else: - n_setpoint = 3000.0 + # 斜坡下降到3000 + n_setpoint = 3200.0 - 200.0 * (t - 92.0) / 8.0 # 外部转矩扰动 - 模拟发动机并联/涡轴输出变化 if t < 5.0: @@ -1096,38 +885,61 @@ def run_motor_battery_coupled_test() -> Dict[str, Any]: else: t_ext = 0.0 - # 负载转矩 - 模拟风速/气压变化 - if t < 10.0: - t_load = 20.0 - elif t < 22.0: - t_load = 80.0 - elif t < 32.0: - t_load = 40.0 - elif t < 42.0: - t_load = 120.0 - elif t < 52.0: - t_load = 50.0 - elif t < 62.0: - t_load = 90.0 - elif t < 72.0: - t_load = 60.0 - elif t < 82.0: - t_load = 130.0 - else: - t_load = 70.0 + # 螺旋桨负载模型：气动阻力矩 + 惯性负载 + # 气动阻力: T = k_drag * n² (与转速平方成正比) + # 惯性负载: T = J_prop * alpha (与加速度成正比，模拟螺旋桨惯性) + n_current = n_setpoint # 使用目标转速计算负载 + w_current = n_current * 2 * np.pi / 60.0 # rad/s - # 外部电功率扰动 - if 35.0 < t < 40.0: + # 计算角加速度 (rad/s²) + alpha = (w_current - w_prev) / dt if dt > 0 else 0.0 + + # 气动阻力矩 (与转速平方成正比) + t_aero = k_drag * n_current ** 2 + T_min + # 惯性负载 (与加速度成正比) + t_inertia = J_prop * alpha + + t_load = t_aero + t_inertia + + # 更新上一时刻转速 + w_prev = w_current + + # 外部电功率扰动 (发动机/燃气轮机发电功率注入，正值表示充电，负值表示负载) + if 5.0 < t < 10.0: + p_ext_elec = 40.0 # 涡轴发电机注入 40kW + elif 10.0 < t < 15.0: + p_ext_elec = 50.0 # 涡轴发电机注入 50kW + elif 20.0 < t < 23.0: + p_ext_elec = 60.0 # 涡轴发电机注入 60kW + elif 25.0 < t < 30.0: + p_ext_elec = 70.0 # 涡轴发电机注入 70kW + elif 32.0 < t < 35.0: + p_ext_elec = 55.0 # 涡轴发电机注入 55kW + elif 35.0 < t < 40.0: p_ext_elec = 80.0 # 涡轴发电机注入 80kW + elif 45.0 < t < 48.0: + p_ext_elec = 90.0 # 涡轴发电机注入 90kW (短时高功率充电) + elif 50.0 < t < 53.0: + p_ext_elec = 45.0 # 涡轴发电机注入 45kW elif 55.0 < t < 60.0: - p_ext_elec = -100.0 # 其他设备用电 100kW + p_ext_elec = -100.0 # 其他设备用电 100kW (负载) + elif 62.0 < t < 65.0: + p_ext_elec = 50.0 # 涡轴发电机注入 50kW + elif 65.0 < t < 70.0: + p_ext_elec = 65.0 # 涡轴发电机注入 65kW + elif 72.0 < t < 75.0: + p_ext_elec = 55.0 # 涡轴发电机注入 55kW elif 75.0 < t < 80.0: p_ext_elec = 60.0 # 涡轴发电机注入 60kW + elif 82.0 < t < 85.0: + p_ext_elec = 70.0 # 涡轴发电机注入 70kW + elif 88.0 < t < 95.0: + p_ext_elec = 75.0 # 涡轴发电机注入 75kW # 母线电压波动 v_bus = battery.V_t + 15.0 * np.sin(2 * np.pi * t / 8.0) + 5.0 * np.sin(2 * np.pi * t / 3.0) - # 延迟一拍的实际电能参与当前动态步，进而PID在内环打出对下一拍的请求 + # 延迟一拍的实际电能参与当前动态步，进而MPC在内环打出对下一拍的请求 state = motor.step(dt, n_setpoint, p_actual_kw, v_bus, t_load, t_ext) p_req_kw = state['p_bus_req_kw'] @@ -1174,7 +986,7 @@ def run_motor_battery_coupled_test() -> Dict[str, Any]: axes[0].legend() # 2. 转矩响应 - axes[1].plot(time_array, t_motor_log, 'r-', lw=1.6, label='Motor Torque (PID Output)') + axes[1].plot(time_array, t_motor_log, 'r-', lw=1.6, label='Motor Torque (MPC Output)') axes[1].plot(time_array, np.array(t_load_log), 'k--', lw=1.2, label='Propeller Load Torque') axes[1].plot(time_array, np.array(t_ext_log), color='tab:purple', lw=1.2, label='External Torque') axes[1].axhline(0.0, color='k', lw=0.6) @@ -1210,7 +1022,7 @@ def run_motor_battery_coupled_test() -> Dict[str, Any]: axes[4].legend() plt.tight_layout() - plt.savefig("figures/test_complex_output.png") + plt.savefig("figures/test_complex_output.png", dpi=300, bbox_inches='tight') print("Plot saved to figures/test_complex_output.png") return { diff --git a/Model/src/mpc_controller.py b/Model/src/mpc_controller.py new file mode 100644 index 0000000..8118740 --- /dev/null +++ b/Model/src/mpc_controller.py @@ -0,0 +1,298 @@ +import numpy as np + + +# ============================================================================== +# Pure-Python projected gradient descent for box-constrained optimization +# (替代 scipy.optimize.minimize SLSQP，避免 Fortran ABI 兼容性问题) +# ============================================================================== +def _minimize_box(objective, x0, bounds, lr=1.0, max_iter=80, ftol=1e-6): + """ + Projected gradient descent with Armijo backtracking line search. + For small-to-medium MPC horizons (H ≤ 30) this is fast enough. + """ + n = len(x0) + lb = np.array([b[0] for b in bounds], dtype=np.float64) + ub = np.array([b[1] for b in bounds], dtype=np.float64) + x = np.clip(np.array(x0, dtype=np.float64), lb, ub) + + eps = 1e-5 # finite-difference step + f_prev = objective(x) + + for _it in range(max_iter): + # Approximate gradient via central differences + grad = np.empty(n, dtype=np.float64) + for i in range(n): + x_p = x.copy(); x_p[i] += eps + x_m = x.copy(); x_m[i] -= eps + grad[i] = (objective(x_p) - objective(x_m)) / (2 * eps) + + # Backtracking line search (Armijo condition) + step = lr + for _ in range(12): + x_new = np.clip(x - step * grad, lb, ub) + f_new = objective(x_new) + if f_new < f_prev - 1e-4 * step * np.dot(grad, x - x_new): + break + step *= 0.5 + else: + x_new = np.clip(x - step * grad, lb, ub) + f_new = objective(x_new) + + if abs(f_prev - f_new) < ftol: + x = x_new + break + x = x_new + f_prev = f_new + + class _Result: + pass + res = _Result() + res.x = x + res.fun = f_prev + res.success = True + return res + +# ============================================================================== +# 1. 涡轴发动机 MPC 控制器 (基于 Scipy SLSQP 高速求解) +# ============================================================================== +class TurboShaftMPCController: + def __init__(self, tau_fuel, K_inertia, dt=0.02, horizon=15, + min_fuel=10.0, max_fuel=600.0, overshoot_limit=0.05): + """ + 初始化高速 MPC 控制器 (基于 Scipy SLSQP) + 替代 GEKKO 以消除文件 I/O 开销，提升单步推理速度 50x 以上。 + + :param overshoot_limit: 超调量硬约束 (0.05 = 5%) + """ + self.dt = dt + self.H = horizon + self.tau_fuel = tau_fuel + self.K_inertia = K_inertia + self.min_fuel = min_fuel + self.max_fuel = max_fuel + + # 权重参数 (与 GEKKO 版本保持一致) + self.W_power = 200.0 # Power Setpoint Weight + self.W_dcost = 1.5 # Delta Control Weight (DCOST) + + # 超调量约束 + self.overshoot_limit = overshoot_limit # 5% 硬约束 + + # 超调惩罚权重 (使用很大的值使约束"硬"化) + self.W_overshoot_penalty = 1e6 + + # 状态缓存 + self.last_u = 100.0 + self.last_N = 0.0 + + def reset(self, initial_output, initial_N): + self.last_u = initial_output + self.last_N = initial_N + + def compute(self, current_N, current_Wfact, target_power, engine_model=None, H_env=0, Ma_env=0, precalc_params=None, **kwargs): + """ + 计算 MPC 控制律 + :param precalc_params: (必须) 元组 (Wf_req_0, Power_0, k_wf, k_p)。 + FastMPC 必须配合 Batch Prediction 使用。 + """ + if precalc_params is not None: + Wf_req_0, Power_0, k_wf, k_p = precalc_params + else: + raise ValueError("FastMPC 必须配合 Batch Prediction 使用 (提供 precalc_params)") + + N0 = current_N + Wf_act0 = current_Wfact + N_ref = N0 # 用于线性化的参考点 + + # --- 构建线性预测模型 --- + # State x = [N, Wf_act] + # x_{k+1} = A x_k + B u_k + d + # N_{k+1} = N_k + dt * K * (Wf_act_k - (Wf_req_0 + k_wf*(N_k - N0))) + # = (1 - dt*K*k_wf) N_k + (dt*K) Wf_act_k + dt*K*(k_wf*N0 - Wf_req_0) + # Wf_act_{k+1} = Wf_act_k + dt * (u_k - Wf_act_k) / tau + # = (1 - dt/tau) Wf_act_k + (dt/tau) u_k + + dt = self.dt + K = self.K_inertia + tau = self.tau_fuel + + A = np.array([ + [1 - dt * K * k_wf, dt * K], + [0, 1 - dt / tau] + ]) + B = np.array([0, dt / tau]) + d = np.array([dt * K * (k_wf * N_ref - Wf_req_0), 0]) + + # Power Output: P = P0 + k_p * (N - N0) + # = k_p * N + (P0 - k_p * N0) + C_p = k_p + D_p = Power_0 - k_p * N_ref + + # --- 超调量硬约束 --- + # 升功率时: P_k <= target_power * (1 + overshoot_limit) + # 降功率时: P_k >= target_power * (1 - overshoot_limit) + P_max = target_power * (1 + self.overshoot_limit) + P_min = target_power * (1 - self.overshoot_limit) + + # 判断是升功率还是降功率 + is_ramping_up = target_power > Power_0 + + # --- 优化目标函数 --- + # Variables: U = [u_0, ..., u_{H-1}] + # x_0 is fixed. + # Cost = sum_{k=1 to H} W_power * (P_k - P_target)^2 + sum_{k=0 to H-1} W_dcost * (u_k - u_{k-1})^2 + # + W_overshoot_penalty * (违反约束的惩罚) + + def objective(U): + cost = 0.0 + x_k = np.array([N0, Wf_act0]) + u_prev = self.last_u # Use last commanded u for the first delta + + for k in range(self.H): + u_k = U[k] + + # Dynamics Step + x_k = A @ x_k + B * u_k + d + + # Output + P_k = C_p * x_k[0] + D_p + + # Cost Accumulation + cost += self.W_power * (P_k - target_power) ** 2 + cost += self.W_dcost * (u_k - u_prev) ** 2 + + # 超调惩罚: 如果违反约束，添加巨大惩罚 + if is_ramping_up: + if P_k > P_max: + cost += self.W_overshoot_penalty * (P_k - P_max) ** 2 + else: + if P_k < P_min: + cost += self.W_overshoot_penalty * (P_min - P_k) ** 2 + + u_prev = u_k + + return cost + + # --- 求解 --- + # 初始猜测: 保持上一次的输入 + U0 = np.full(self.H, float(self.last_u), dtype=np.float64) + + # 约束: lb <= u <= ub + bounds = [(self.min_fuel, self.max_fuel) for _ in range(self.H)] + + # 使用纯Python投影梯度下降求解带约束优化 (替代 SLSQP) + res = _minimize_box(objective, U0, bounds, lr=5.0, max_iter=60, ftol=1e-3) + + # 更新状态 + u_opt = res.x[0] + self.last_u = u_opt + + return u_opt + + +# ============================================================================== +# 2. 驱动电机 MPC 控制器 (基于 Scipy SLSQP 高速求解) +# ============================================================================== +class MotorMPCController: + def __init__(self, J, B_visc, dt=0.02, horizon=10, W_speed=100.0, W_dcost=1.0, overshoot_limit=0.05): + """ + 初始化高速电机 MPC 控制器 (基于 Scipy SLSQP) + 替代 GEKKO 以消除文件 I/O 开销，提升单步推理速度 50x 以上。 + + :param overshoot_limit: 超调量硬约束 (0.05 = 5%) + """ + self.dt = dt + self.H = horizon + self.J = J + self.B_visc = B_visc + + # 权重参数 (可调) + self.W_speed = W_speed # 转速跟踪权重 + self.W_dcost = W_dcost # 控制变化惩罚 + + # 超调量约束 + self.overshoot_limit = overshoot_limit # 5% 硬约束 + + # 超调惩罚权重 (使用很大的值使约束"硬"化) + self.W_overshoot_penalty = 1e6 + + # 状态缓存 + self.last_T_cmd = 0.0 + + def reset(self, initial_w=0.0): + """重新初始化状态""" + self.last_T_cmd = 0.0 + + def compute(self, current_w, target_w, t_load, t_ext, t_lim_upper, t_lim_lower): + """ + 执行单步 MPC 优化计算 + :param current_w: 当前实际转速 (rad/s) + :param target_w: 目标转速 (rad/s) + :param t_load: 当前气动负载转矩 (Nm) + :param t_ext: 当前外部轴系转矩 (Nm) + :param t_lim_upper: 当前母线电压下，电机能发出的【正向转矩上限】(发电极限) + :param t_lim_lower: 当前母线电压下，电机能发出的【负向转矩下限】(驱动极限) + :return: 最优转矩指令 (Nm) + """ + w0 = current_w + + dt = self.dt + J = self.J + B = self.B_visc + + # --- 离散线性模型 --- + # w_{k+1} = w_k + dt * (-(T_cmd + T_ext + T_load + B*w_k) / J) + # = (1 - dt*B/J) * w_k - dt/J * T_cmd - dt/J * (T_ext + T_load) + a = 1 - dt * B / J + b = -dt / J + d = -dt / J * (t_ext + t_load) + + # --- 超调量硬约束 --- + # 升速时: w_k <= target_w * (1 + overshoot_limit) + # 降速时: w_k >= target_w * (1 - overshoot_limit) + w_max = target_w * (1 + self.overshoot_limit) + w_min = target_w * (1 - self.overshoot_limit) + + # 判断是升速还是降速 + is_ramping_up = target_w > w0 + + # --- 优化目标函数 --- + # Cost = sum_{k=1 to H} W_speed * (w_k - target_w)^2 + sum_{k=0 to H-1} W_dcost * (T_k - T_{k-1})^2 + # + W_overshoot_penalty * (违反约束的惩罚) + + def objective(U): + cost = 0.0 + w_k = w0 + T_prev = self.last_T_cmd + + for k in range(self.H): + T_k = U[k] + w_k = a * w_k + b * T_k + d # 动力学迭代 + cost += self.W_speed * (w_k - target_w) ** 2 + cost += self.W_dcost * (T_k - T_prev) ** 2 + + # 超调惩罚: 如果违反约束，添加巨大惩罚 + if is_ramping_up: + if w_k > w_max: + cost += self.W_overshoot_penalty * (w_k - w_max) ** 2 + else: + if w_k < w_min: + cost += self.W_overshoot_penalty * (w_min - w_k) ** 2 + + T_prev = T_k + + return cost + + # 初始猜测 + U0 = np.full(self.H, float(self.last_T_cmd), dtype=np.float64) + + # 约束: 转矩边界 + bounds = [(t_lim_lower, t_lim_upper) for _ in range(self.H)] + + # 求解 (使用纯Python投影梯度下降，替代 SLSQP) + res = _minimize_box(objective, U0, bounds, lr=50.0, max_iter=120, ftol=1e-5) + + T_opt = res.x[0] + self.last_T_cmd = T_opt + + return T_opt diff --git a/Model/src/series_hybrid_sim.py b/Model/src/series_hybrid_sim.py index 7071cd1..17ce25d 100644 --- a/Model/src/series_hybrid_sim.py +++ b/Model/src/series_hybrid_sim.py @@ -1,6 +1,6 @@ -import os import numpy as np import matplotlib.pyplot as plt +from tqdm import tqdm from engine_dynamic_sim import TurboshaftDynamicSim from motor_sim import MotorSim @@ -10,25 +10,15 @@ class SeriesHybridSystem: """ 串联式混合动力系统总成 """ - def __init__(self): - # ===== 新增：按Model目录定位GPR数据与权重，避免从主项目调用时路径错误 ===== - model_root = os.path.dirname(os.path.dirname(os.path.abspath(__file__))) - gpr_csv_path = os.path.join(model_root, "data", "Cleaned_Engine_Data_Full.csv") - gpr_pth_path = os.path.join(model_root, "data", "engine_gpr_model.pth") - self.genset = TurboshaftDynamicSim( - gpr_csv_path=gpr_csv_path, - gpr_pth_path=gpr_pth_path, - kp=0.5, - ki=0.5, - kd=0.0 - ) + def __init__(self, mpc_overshoot_limit: float = 0.05): + self.genset = TurboshaftDynamicSim(mpc_overshoot_limit=mpc_overshoot_limit) self.drive_motor = MotorSim( P_rate=300e3, w_rate=575.95, - k_p_w=5.0, # 增大比例增益 - k_i_w=2.0, # 增大积分增益 - k_d_w=0.5, # 适度微分 J=1.0, + mpc_W_speed=40, + mpc_W_dcost=5, + mpc_overshoot_limit=mpc_overshoot_limit, # 超调量硬约束 (5%) ) self.battery = BatterySim(capacity_kwh=50.0, initial_soc=0.6) @@ -172,7 +162,7 @@ if __name__ == "__main__": 'p_batt_actual_kw', 'wf_kg_h']} print("-> 开始全系统闭环步进仿真 (总时长 3 分钟)...") - for t in time_array: + for t in tqdm(time_array, desc='Simulating', unit='step'): # 3分钟测试剖面 if t < 15.0: target_rpm, load_torque = 1500.0, 50.0 # 地面滑行 @@ -204,7 +194,7 @@ if __name__ == "__main__": # 2. 功率分配流向 axes[1].plot(time_array, log['p_drive_req_kw'], 'k--', lw=1.5, label='Drive Motor Request') - axes[1].plot(time_array, log['p_engine_out_kw'], 'r-', lw=1.5, label='Engine Output (APU)') + axes[1].plot(time_array, log['p_engine_out_kw'], 'r-', lw=1.5, label='Engine Output (Turbo Shaft)') axes[1].plot(time_array, log['p_batt_actual_kw'], 'g-', lw=1.5, label='Battery Output') axes[1].set_ylabel('Power [kW]') axes[1].set_title('System Power Flow (Energy Management)') @@ -249,7 +239,7 @@ if __name__ == "__main__": os.makedirs(data_dir) dat_path = os.path.join(data_dir, 'series_hybrid_data.dat') - with open(dat_path, 'w') as f: + with open(dat_path, 'w', encoding='utf-8') as f: f.write('# Time(s)\tTarget_RPM\tActual_RPM\tError_RPM\tSOC\tEngine_Power\tDrive_Req\tBatt_Power\n') for i in range(len(time_array)): err = log['target_prop_rpm'][i] - log['prop_speed_rpm'][i] diff --git a/app.py b/app.py index 976d8ef..2506725 100644 --- a/app.py +++ b/app.py @@ -1,3 +1,4 @@ +import os import gradio as gr import time from functools import partial @@ -10,8 +11,8 @@ from analysis_functions import ( frequency_domain_analysis, root_locus_analysis ) -# ===== 新增：算例演示模块函数导入 ===== -from case_demo_functions import run_case_demo +# ===== 算例演示模块函数导入（四阶段设计）===== +from case_demo_functions import run_distillation_demo, run_gpr_training, run_engine_design, run_motor_design, run_hybrid_demo from chatbot import chat_with_ai from user_stats import get_online_status_html, update_user_activity from ui_components import ( @@ -23,8 +24,98 @@ from ui_components import ( create_chatbot_tab ) + +# ===== 系统资源监控 ===== +def get_system_monitor_html(): + """获取 CPU / 内存 / GPU 使用率的 HTML 小组件""" + try: + import psutil + cpu_pct = psutil.cpu_percent(interval=0) + mem = psutil.virtual_memory() + mem_pct = mem.percent + mem_used_gb = mem.used / (1024 ** 3) + mem_total_gb = mem.total / (1024 ** 3) + except ImportError: + return "
psutil 未安装，无法监控系统资源
" + + # GPU 信息 — 优先用 nvidia-smi（不依赖 PyTorch CUDA 版本），再用 torch.cuda 兜底 + gpu_html = "" + try: + import subprocess as _sp + _r = _sp.run( + ['nvidia-smi', '--query-gpu=name,memory.used,memory.total,utilization.gpu', + '--format=csv,noheader,nounits'], + capture_output=True, text=True, timeout=3 + ) + if _r.returncode == 0 and _r.stdout.strip(): + _parts = [p.strip() for p in _r.stdout.strip().split('\n')[0].split(',')] + _gpu_mem_used = float(_parts[1]) / 1024 # MiB → GiB + _gpu_mem_total = float(_parts[2]) / 1024 + _gpu_util = float(_parts[3]) + _gc = '#ff6b6b' if _gpu_util > 80 else '#ffd93d' if _gpu_util > 50 else '#6bcb77' + gpu_html = ( + f"
" + f"🎮 GPU" + f"
" + f"
" + f"
" + f"{_gpu_util:.0f}%   {_gpu_mem_used:.1f}/{_gpu_mem_total:.0f} GB" + f"
" + ) + else: + raise RuntimeError("nvidia-smi no output") + except Exception: + try: + import torch as _torch + if _torch.cuda.is_available(): + _mem_alloc = _torch.cuda.memory_allocated(0) / (1024 ** 3) + _mem_total = _torch.cuda.get_device_properties(0).total_memory / (1024 ** 3) + _util = _mem_alloc / max(_mem_total, 0.01) * 100 + _gc = '#ff6b6b' if _util > 80 else '#ffd93d' if _util > 50 else '#6bcb77' + gpu_html = ( + f"
" + f"🎮 GPU" + f"
" + f"
" + f"
" + f"{_mem_alloc:.1f}/{_mem_total:.0f} GB" + f"
" + ) + else: + gpu_html = "
🎮 GPU N/A
" + except Exception: + gpu_html = "
🎮 GPU N/A
" + + cpu_color = '#ff6b6b' if cpu_pct > 80 else '#ffd93d' if cpu_pct > 50 else '#6bcb77' + mem_color = '#ff6b6b' if mem_pct > 80 else '#ffd93d' if mem_pct > 50 else '#6bcb77' + + html = ( + f"
" + # CPU + f"
" + f"🖥️ CPU" + f"
" + f"
" + f"
" + f"{cpu_pct:.0f}%" + f"
" + # Memory + f"
" + f"💾 RAM" + f"
" + f"
" + f"
" + f"{mem_used_gb:.1f}/{mem_total_gb:.0f} GB ({mem_pct:.0f}%)" + f"
" + # GPU + f"{gpu_html}" + f"
" + ) + return html + # 加载外部CSS文件 -with open("assets/styles.css", "r", encoding="utf-8") as f: +with open(os.path.join(os.path.dirname(__file__), "assets", "styles.css"), "r", encoding="utf-8") as f: custom_css = f.read() # --- 主应用界面 --- @@ -35,24 +126,9 @@ with gr.Blocks(title="自动控制理论学习网站 - AI+数智平台", css=cus # 创建头部信息和在线计数器 online_counter = create_header() - - # 创建共享的输入组件 - with gr.Row(): - with gr.Column(scale=1): - with gr.Group(): - gr.HTML("
📊 通用系统参数
") - num_input = gr.Textbox( - label="传递函数分子系数 (Numerator)", - value="1", - placeholder="例如: 1 或 1,2,3", - info="💡 用逗号分隔多个系数，从最高次项到常数项" - ) - den_input = gr.Textbox( - label="传递函数分母系数 (Denominator)", - value="1,6,11,6", - placeholder="例如: 1,2,1", - info="💡 分母阶数通常高于或等于分子阶数" - ) + + # 系统资源监控（始终可见） + system_monitor = gr.HTML(value=get_system_monitor_html, elem_id="system-monitor") # 创建功能选项卡 with gr.Tabs() as tabs: @@ -80,13 +156,13 @@ with gr.Blocks(title="自动控制理论学习网站 - AI+数智平台", css=cus # --- 时域分析事件 --- time_domain_ui["confirm_button"].click( fn=display_transfer_function, - inputs=[num_input, den_input], + inputs=[time_domain_ui["num_input"], time_domain_ui["den_input"]], outputs=[time_domain_ui["tf_display"]] ).then(lambda: get_online_status_html(), outputs=online_counter) time_domain_ui["analyze_button"].click( fn=time_domain_analysis, - inputs=[num_input, den_input], + inputs=[time_domain_ui["num_input"], time_domain_ui["den_input"]], outputs=[time_domain_ui["output_plot"], time_domain_ui["output_metrics"]] ).then(lambda: get_online_status_html(), outputs=online_counter) @@ -96,7 +172,7 @@ with gr.Blocks(title="自动控制理论学习网站 - AI+数智平台", css=cus fig, metrics, tf_latex, stability = frequency_domain_analysis(num, den, k) return fig, metrics, tf_latex, stability, k, get_online_status_html() - freq_inputs = [num_input, den_input, freq_domain_ui["log_k_slider"]] + freq_inputs = [freq_domain_ui["num_input"], freq_domain_ui["den_input"], freq_domain_ui["log_k_slider"]] freq_outputs = [ freq_domain_ui["plot_output"], freq_domain_ui["metrics_display"], @@ -116,7 +192,7 @@ with gr.Blocks(title="自动控制理论学习网站 - AI+数智平台", css=cus fig, poles, k_val = root_locus_analysis(num, den, log_k) return fig, poles, k_val, get_online_status_html() - rl_inputs = [root_locus_ui["log_k_slider"], num_input, den_input] + rl_inputs = [root_locus_ui["log_k_slider"], root_locus_ui["num_input"], root_locus_ui["den_input"]] rl_outputs = [ root_locus_ui["plot_output"], root_locus_ui["poles_display"], @@ -129,59 +205,170 @@ with gr.Blocks(title="自动控制理论学习网站 - AI+数智平台", css=cus outputs=rl_outputs ) - # 当输入框变化时，也更新频域和根轨迹（如果它们是当前可见的） - def update_all_on_tf_change(num, den, log_k_freq, log_k_rl): - # 更新频域 - k_freq = 10**log_k_freq - fig_freq, metrics, tf_latex, stability = frequency_domain_analysis(num, den, k_freq) - - # 更新根轨迹 - fig_rl, poles, k_val_rl = root_locus_analysis(num, den, log_k_rl) + # 频域：当传递函数输入框变化时自动更新 + freq_domain_ui["num_input"].change( + fn=update_frequency_analysis_wrapper, + inputs=freq_inputs, outputs=freq_outputs + ) + freq_domain_ui["den_input"].change( + fn=update_frequency_analysis_wrapper, + inputs=freq_inputs, outputs=freq_outputs + ) - return ( - fig_freq, metrics, tf_latex, stability, k_freq, - fig_rl, poles, k_val_rl, - get_online_status_html() - ) + # 根轨迹：当传递函数输入框变化时自动更新 + root_locus_ui["num_input"].change( + fn=update_rl_view_wrapper, + inputs=rl_inputs, outputs=rl_outputs + ) + root_locus_ui["den_input"].change( + fn=update_rl_view_wrapper, + inputs=rl_inputs, outputs=rl_outputs + ) - tf_change_inputs = [num_input, den_input, freq_domain_ui["log_k_slider"], root_locus_ui["log_k_slider"]] - tf_change_outputs = freq_outputs[:-1] + rl_outputs[:-1] + [online_counter] - - num_input.change(fn=update_all_on_tf_change, inputs=tf_change_inputs, outputs=tf_change_outputs) - den_input.change(fn=update_all_on_tf_change, inputs=tf_change_inputs, outputs=tf_change_outputs) + # ===== 算例演示事件绑定（四阶段）===== - # ===== 新增：算例演示事件包装器 ===== - def run_case_demo_wrapper(sim_time, dt, initial_soc, initial_engine_power, profile, rpm_scale, load_scale, sid): + # --- 阶段零-A：GPR 模型训练 --- + def run_gpr_wrapper(mode, sid, progress=gr.Progress(track_tqdm=True)): update_user_activity(sid) - fig, summary, table_data = run_case_demo( - sim_time_s=sim_time, - dt=dt, - initial_soc_pct=initial_soc, - initial_engine_power_kw=initial_engine_power, - profile_name=profile, - rpm_scale=rpm_scale, - load_scale=load_scale + fig, summary = run_gpr_training(mode=mode, progress=progress) + return fig, summary, get_online_status_html() + + case_demo_ui["gpr_run_button"].click( + fn=run_gpr_wrapper, + inputs=[case_demo_ui["gpr_mode"], session_id], + outputs=[case_demo_ui["gpr_plot"], case_demo_ui["gpr_summary"], online_counter] + ) + + # --- 阶段零-B：NN 模型训练（蒸馏）--- + def run_distillation_wrapper(epochs, lr, hidden, sid, progress=gr.Progress(track_tqdm=True)): + update_user_activity(sid) + fig, summary = run_distillation_demo(epochs, lr, hidden, progress=progress) + return fig, summary, get_online_status_html() + + case_demo_ui["distill_run_button"].click( + fn=run_distillation_wrapper, + inputs=[ + case_demo_ui["distill_epochs"], case_demo_ui["distill_lr"], + case_demo_ui["distill_hidden"], + session_id + ], + outputs=[case_demo_ui["distill_plot"], case_demo_ui["distill_summary"], online_counter] + ) + + # --- 阶段一：发动机控制器设计 --- + def run_engine_design_wrapper(sim_time, dt, init_power, target_power, + controller_type, + kp, ki, kd, tau_fuel, K_inertia, + mpc_horizon, mpc_W_power, mpc_W_dcost, mpc_overshoot, + sid, progress=gr.Progress(track_tqdm=True)): + update_user_activity(sid) + fig, summary = run_engine_design( + sim_time, dt, init_power, target_power, + controller_type, + kp, ki, kd, tau_fuel, K_inertia, + mpc_horizon, mpc_W_power, mpc_W_dcost, mpc_overshoot / 100.0, + progress=progress + ) + return fig, summary, get_online_status_html() + + case_demo_ui["eng_run_button"].click( + fn=run_engine_design_wrapper, + inputs=[ + case_demo_ui["eng_sim_time"], case_demo_ui["eng_dt"], + case_demo_ui["eng_init_power"], case_demo_ui["eng_target_power"], + case_demo_ui["eng_controller_type"], + case_demo_ui["eng_kp"], case_demo_ui["eng_ki"], case_demo_ui["eng_kd"], + case_demo_ui["eng_tau_fuel"], case_demo_ui["eng_K_inertia"], + case_demo_ui["eng_mpc_horizon"], case_demo_ui["eng_mpc_W_power"], + case_demo_ui["eng_mpc_W_dcost"], case_demo_ui["eng_mpc_overshoot"], + session_id + ], + outputs=[case_demo_ui["eng_plot"], case_demo_ui["eng_summary"], online_counter] + ) + + # --- 阶段二：电机控制器设计 --- + def run_motor_design_wrapper(sim_time, dt, target_rpm, load_torque, + controller_type, + kp, ki, kd, J, + mpc_W_speed, mpc_W_dcost, mpc_overshoot, + sid, progress=gr.Progress(track_tqdm=True)): + update_user_activity(sid) + fig, summary = run_motor_design( + sim_time, dt, target_rpm, load_torque, + controller_type, + kp, ki, kd, J, + mpc_W_speed, mpc_W_dcost, mpc_overshoot / 100.0, + progress=progress + ) + return fig, summary, get_online_status_html() + + case_demo_ui["mot_run_button"].click( + fn=run_motor_design_wrapper, + inputs=[ + case_demo_ui["mot_sim_time"], case_demo_ui["mot_dt"], + case_demo_ui["mot_target_rpm"], case_demo_ui["mot_load_torque"], + case_demo_ui["mot_controller_type"], + case_demo_ui["mot_kp"], case_demo_ui["mot_ki"], case_demo_ui["mot_kd"], + case_demo_ui["mot_J"], + case_demo_ui["mot_mpc_W_speed"], case_demo_ui["mot_mpc_W_dcost"], + case_demo_ui["mot_mpc_overshoot"], + session_id + ], + outputs=[case_demo_ui["mot_plot"], case_demo_ui["mot_summary"], online_counter] + ) + + # --- 阶段三：能量管理策略设计（自动引用前两阶段控制器参数）--- + def run_hybrid_demo_wrapper(sim_time, dt, initial_soc, initial_engine_power, + profile, + eng_ctrl_type, eng_kp, eng_ki, eng_kd, + eng_mpc_horizon, eng_mpc_W_power, eng_mpc_W_dcost, eng_mpc_overshoot, + mot_ctrl_type, mot_kp, mot_ki, mot_kd, mot_J, + mot_mpc_W_speed, mot_mpc_W_dcost, mot_mpc_overshoot, + soc_target, soc_low, soc_high, + p_eng_min, p_eng_max, p_charge, k_soc, + power_reserve, battery_capacity, sid, + progress=gr.Progress(track_tqdm=True)): + update_user_activity(sid) + fig, summary, table_data = run_hybrid_demo( + sim_time, dt, initial_soc, initial_engine_power, profile, + eng_ctrl_type, eng_kp, eng_ki, eng_kd, + eng_mpc_horizon, eng_mpc_W_power, eng_mpc_W_dcost, eng_mpc_overshoot / 100.0, + mot_ctrl_type, mot_kp, mot_ki, mot_kd, mot_J, + mot_mpc_W_speed, mot_mpc_W_dcost, mot_mpc_overshoot / 100.0, + soc_target, soc_low, soc_high, + p_eng_min, p_eng_max, p_charge, k_soc, + power_reserve, battery_capacity, + progress=progress ) return fig, summary, table_data, get_online_status_html() - # ===== 新增：算例演示按钮事件绑定 ===== - case_demo_ui["run_button"].click( - fn=run_case_demo_wrapper, + case_demo_ui["hybrid_run_button"].click( + fn=run_hybrid_demo_wrapper, inputs=[ - case_demo_ui["sim_time"], - case_demo_ui["dt"], - case_demo_ui["initial_soc"], - case_demo_ui["initial_engine_power"], + case_demo_ui["sim_time"], case_demo_ui["dt"], + case_demo_ui["initial_soc"], case_demo_ui["initial_engine_power"], case_demo_ui["profile"], - case_demo_ui["rpm_scale"], - case_demo_ui["load_scale"], + # 发动机控制器参数 + case_demo_ui["eng_controller_type"], + case_demo_ui["eng_kp"], case_demo_ui["eng_ki"], case_demo_ui["eng_kd"], + case_demo_ui["eng_mpc_horizon"], case_demo_ui["eng_mpc_W_power"], + case_demo_ui["eng_mpc_W_dcost"], case_demo_ui["eng_mpc_overshoot"], + # 电机控制器参数 + case_demo_ui["mot_controller_type"], + case_demo_ui["mot_kp"], case_demo_ui["mot_ki"], case_demo_ui["mot_kd"], + case_demo_ui["mot_J"], + case_demo_ui["mot_mpc_W_speed"], case_demo_ui["mot_mpc_W_dcost"], + case_demo_ui["mot_mpc_overshoot"], + # 能量管理策略参数 + case_demo_ui["soc_target"], case_demo_ui["soc_low"], case_demo_ui["soc_high"], + case_demo_ui["p_eng_min"], case_demo_ui["p_eng_max"], + case_demo_ui["p_charge"], case_demo_ui["k_soc"], + case_demo_ui["power_reserve"], case_demo_ui["battery_capacity"], session_id ], outputs=[ - case_demo_ui["plot"], - case_demo_ui["summary"], - case_demo_ui["table"], - online_counter + case_demo_ui["hybrid_plot"], case_demo_ui["hybrid_summary"], + case_demo_ui["hybrid_table"], online_counter ] ) @@ -215,14 +402,13 @@ with gr.Blocks(title="自动控制理论学习网站 - AI+数智平台", css=cus ) # --- 页面加载和定时器事件 --- - def on_page_load(sid): - update_user_activity(sid) - return get_online_status_html() - - demo.load(fn=on_page_load, inputs=[session_id], outputs=[online_counter]) + demo.load(fn=lambda: get_online_status_html(), outputs=[online_counter]) gr.Timer(10).tick(fn=get_online_status_html, outputs=online_counter) + # 系统资源监控定时刷新（每 3 秒） + gr.Timer(3).tick(fn=get_system_monitor_html, outputs=system_monitor) + if __name__ == "__main__": demo.queue().launch( diff --git a/assets/knowledge_cards_html.py b/assets/knowledge_cards_html.py index eba761e..cc19ef3 100644 --- a/assets/knowledge_cards_html.py +++ b/assets/knowledge_cards_html.py @@ -1,5 +1,5 @@ """ -知识卡片内容 - 包含时域、频域、根轨迹分析的公式和方法 +知识卡片内容 - 包含时域、频域、根轨迹、控制器设计、模型蒸馏、能量管理分析的公式和方法 使用纯HTML格式，无需LaTeX渲染库 """ @@ -497,3 +497,1676 @@ ROOT_LOCUS_KNOWLEDGE = """ """ + + +# ============================================================ +# 发动机控制器设计知识卡片 (PID + MPC) — 中文 +# ============================================================ +ENGINE_CONTROL_KNOWLEDGE = """ +
+

+ 🔧 发动机控制器设计理论 +

+ +
+ + 🎛️ PID 控制器 — 增量式形式 + +
+

标准PID传递函数：

+
+
+ u(t) = Kp e(t) + + Ki ∫e(τ)dτ + + Kd + + de(t) + dt + +
+
+

增量式（速度型）PID：

+
+
+ Δu(k) = Kp[e(k) − e(k−1)] + + Ki e(k) + + Kd[e(k) − 2e(k−1) + e(k−2)] +
+
+
+

💡 增量式优势：

+
    +
  • 输出为控制量的增量 Δu，而非绝对值，避免积分饱和（Windup）
    • +
  • 手动/自动切换时输出无跳变，适合工程实际
    • +
  • 计算量恒定，不随运行时间增加
    • +
  • 即使计算失误，影响仅限于当前步的增量
    • +
+
+ +

📊 各参数作用与调节规律

+ + + + + + + + + + + + + + + + + + + + + + + + + +
参数增大效果过大时调节建议
Kp (比例)加快响应速度，减小稳态误差产生振荡甚至不稳定先调大至出现振荡，再回退70%
Ki (积分)彻底消除稳态误差超调加重，积分饱和从小值开始缓慢增大
Kd (微分)抑制超调和振荡，改善动态品质放大噪声，响应迟钝通常取较小值，有噪声时可设为0
+ +
+

🎯 经典调参步骤（Ziegler-Nichols 启发式）：

+
    +
  1. 设 Ki=0, Kd=0，逐渐增大 Kp 直到系统出现等幅振荡
    2. +
  2. 记录此时的临界增益 Ku 和振荡周期 Tu
    3. +
  3. 按经验公式设定：Kp=0.6Ku, Ki=2Kp/Tu, Kd=KpTu/8
    4. +
  4. 根据实际响应微调各参数
    5. +
+
+
+
+ +
+ + 🧠 MPC — 模型预测控制 + +
+

MPC（Model Predictive Control）利用系统的内部动态模型预测未来行为， + 在滚动时域上优化控制序列，仅执行第一步，然后重新优化——这就是"滚动优化"原理。

+ +

🔑 优化目标函数

+
+
+ minU   J = Σk=1H + [ Wp (Pk − Pref)² + + WΔu (Δuk)² ] +
+
+

+ 其中 H 为预测时域长度，Wp 为功率跟踪权重，WΔu 为控制增量惩罚权重。 +

+ +

🔑 MPC 核心概念

+ + + + + + + + + + + + + + + + + + + + + + + + + + +
概念说明调节效果
预测时域 (H)优化器向前看的步数H越大前瞻性越强，但计算量增加
Wpower功率跟踪权重越大跟踪越紧密，但控制信号波动可能增大
WΔcost控制平滑度权重越大控制越平滑，减少执行器磨损
超调限制输出偏差的硬约束（如 ≤5%）越小越安全但可能减缓响应速度
+ +
+

🔬 内部线性化模型：本仿真中，涡轴发动机被线性化为以下离散状态方程：

+
+ Nk+1 = Nk + Δt · Kinertia · (Wf,act − Wf,req)
+ Wf,act,k+1 = (1 − Δt/τ) · Wf,act,k + (Δt/τ) · uk +
+

+ τ 为燃油执行机构时间常数，Kinertia 为转子惯性增益。Wf,req 和功率 P 由 NN 代理模型实时提供。 +

+
+ +
+

⚖️ PID vs MPC 对比：

+ + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + +
特性PIDMPC
原理基于误差的反馈校正（事后调节）基于模型的前馈预测（提前规划）
约束处理需额外加限幅/抗饱和机制在优化中自然处理约束
计算量极低，微秒级较高，需在线优化
模型依赖无需模型需要准确的系统模型
适用场景稳定对象，简单跟踪多变量、有约束的复杂系统
+
+
+
+ +
+ + ✈️ 涡轴发动机模型 + +
+

涡轴发动机由三个耦合子系统建模：

+ + + + + + + + + + + + + + + + + + + + + +
子系统动态特性关键参数
NN 代理模型映射 (高度, 马赫数, 转速) → (燃油流量, 功率)；基于 GPR/数据训练输入3维, 输出2维
燃油执行机构一阶惯性环节: dWf/dt = (u − Wf) / τfuelτ 越大响应越慢
转子动力学dN/dt = Kinertia · (Wf,actual − Wf,required)K 越大加速越快
+ +
+

💡 仿真流程：

+
    +
  1. 给定当前转速 N，NN模型预测所需燃油流量 Wf,req 和当前功率 P
    2. +
  2. 控制器（PID或MPC）根据功率误差计算燃油命令 Wf,cmd
    3. +
  3. 执行机构将命令滤波为实际燃油流量 Wf,act
    4. +
  4. Wf,act − Wf,req 产生净力矩驱动转子加速或减速
    5. +
  5. 转速变化又改变NN模型的输入，形成闭环
    6. +
+
+ +
+

📏 性能评价指标：

+
    +
  • 稳态误差：功率达到稳态后与目标的平均偏差
    • +
  • 超调量 (σ%)：输出超过目标值的最大百分比
    • +
  • 上升时间 (tr)：功率从10%变化到90%所需时间
    • +
  • 调节时间 (ts)：功率进入并持续保持在 ±2% 误差带内的时间
    • +
+
+
+
+ +
+ + 📐 PID 控制器 — Laplace 域传递函数 + +
+

连续域 PID 传递函数：

+
+
+ C(s) = Kp + ( + 1 + + + 1 + Ti s + + + Td s + ) +
+
+

其中 Ti = Kp/Ki 为积分时间常数， + Td = Kd/Kp 为微分时间常数。

+ +

Ziegler-Nichols 整定公式表

+

基于极限灵敏度(临界增益)法——设 Ku 为临界增益，Tu 为临界周期：

+ + + + + + + + + + + + + + + + + + + + + + + + + +
控制器类型KpTiTd
P0.5 Ku——
PI0.45 KuTu / 1.2—
PID0.6 KuTu / 2Tu / 8
+

+ 等价转换：Ki = Kp/Ti，Kd = Kp·Td。 + Z-N 公式通常产生约 25% 超调，实际应用中需进一步微调。 +

+
+
+ +
+ + 🧮 MPC 矩阵化优化问题 + +
+

MPC 将有限时域优化问题转化为二次规划 (QP) 或非线性规划问题：

+ +

离散状态空间模型

+
+
+ xk+1 = Axk + Buk
+ yk = Cxk +
+
+

+ 本系统中：x = [N, Wf,act]T（转速和实际燃油流量）， + u = Wf,cmd（燃油命令），y = P（功率输出）。 +

+ +

预测矩阵展开

+

利用递推关系，将 H 步预测写成矩阵形式：

+
+
+ Y = Ψx0 + ΘU +
+
+
    +
  • Y = [y1, ..., yH]T：预测输出序列
    • +
  • U = [u0, ..., uH−1]T：控制输入序列
    • +
  • Ψ = [CA, CA², ..., CAH]T：自由响应矩阵
    • +
  • Θ：Toeplitz 控制响应矩阵
    • +
+ +

带约束 QP 问题

+
+
+ minU (Y − R)T Q (Y − R) + ΔUT RΔ ΔU

+ s.t.   umin ≤ uk ≤ umax,   + |yk − rk| ≤ εovershoot +
+
+

+ Q 和 RΔ 分别对应界面中的 Wpower 和 WΔcost 权重滑块。 + 本系统使用 SciPy SLSQP 求解器进行约束非线性优化。 +

+
+
+
+""" +MOTOR_CONTROL_KNOWLEDGE = """ +
+

+ ⚡ 电机控制器设计理论 +

+ +
+ + ⚙️ PMSM 永磁同步电机模型 + +
+

永磁同步电机 (PMSM) 是电动推进系统的核心驱动组件，具有功率密度高、效率高、控制精度好等优点。

+ +

机械动力学方程

+
+
+ J + + dω + dt + + = −(Tmotor + Tload + Text + Bω) +
+
+

+ 其中 Tmotor 为电磁转矩（电机产生），Tload 为负载转矩， + Text 为外部扰动，Bω 为粘性摩擦力矩。 +

+ +

电磁转矩表达式

+
+
+ Te = 1.5 · np · ψf · iq +
+
+

+ 对表贴式PMSM (Ld = Lq)，转矩仅与q轴电流 iq 成正比。 +

+ + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + +
符号含义典型值对系统的影响
J转动惯量 (kg·m²)0.5 – 5.0越大响应越慢但越平稳
B粘性摩擦系数 (N·m·s)3×10⁻⁴提供自然阻尼
np极对数3影响电气频率和转矩常数
ψf永磁体磁链 (Wb)0.15决定转矩能力和反电动势
Prate额定功率 (kW)300确定最大可用转矩
+ +
+

💡 电压约束与弱磁： + 当转速升高时，反电动势 e = np·ψf·ω 增大。 + 当反电动势接近母线电压 Vdc 时，电机进入弱磁区， + 需注入负的 id 电流来削弱磁链，此时可用转矩下降。 +

+
+
+
+ +
+ + 🎯 转速控制环路 + +
+

典型电机控制采用双环级联结构：外环为转速环，内环为电流环。

+
+ + ωref → [转速控制器] → Tcmd → [电流环 + SVPWM] → PMSM → ωactual + +
+ +
+

PID 模式：

+
    +
  • 根据转速误差 e = ωref − ω 计算转矩命令 T = Kp·e + Ki·∫e + Kd·ė
    • +
  • 转矩命令转换为功率请求 P = T·ω 送入逆变器
    • +
  • 简单可靠，但无法提前预知负载变化
    • +
+
+
+

MPC 模式：

+
    +
  • 利用电机离散动力学模型预测 H 步后的转速轨迹
    • +
  • 优化 H 步的转矩序列，使转速跟踪误差和控制增量同时最小化
    • +
  • 支持超调硬约束：转速偏差不得超过设定百分比
    • +
+
+ minT0..H Σ [ Wω(ωk − ωref)² + WΔ(Tk − Tk−1)² ] +
+

+ 离散状态方程：ωk+1 = (1−dt·B/J)·ωk − (dt/J)·Tcmd − (dt/J)·(Text+Tload) +

+
+
+
+ +
+ + 🔍 负载扰动测试设计 + +
+

本仿真在运行至 60% 时刻时自动施加 50% 负载阶跃扰动， + 用于检验控制器的抗扰性能。

+ +

📊 抗扰性能指标

+ + + + + + + + + + + + + + + + + + + + + +
指标定义理想范围
转速跌落 (Speed Dip)扰动后转速最大瞬时下降量< 5% 额定转速
恢复时间转速回到 ±2% 设定值的时间< 3 秒
稳态误差扰动后新稳态与目标的偏差< 1% 或指定精度
+ +
+

⚙️ 转动惯量 J 的影响：

+
    +
  • J 较大：响应较慢但转速跌落小，系统惯性大不易被扰动"推走"
    • +
  • J 较小：响应快但对扰动更敏感，转速波动大
    • +
  • 实际工程中，可通过飞轮或增大转子质量来调整等效惯量
    • +
+
+ +
+

💡 MPC 调参建议：

+
    +
  • 如果扰动后转速恢复慢→增大 Wspeed 或减小 WΔcost
    • +
  • 如果转速振荡→增大 WΔcost 或减小 Wspeed
    • +
  • 如果超调过大→降低超调限制百分比
    • +
  • 负载很大时，确认电机额定功率是否足以提供所需转矩
    • +
+
+
+
+ +
+ + 📐 PMSM dq 坐标系电压方程 + +
+

通过 Park 变换将三相 ABC 坐标系转换为旋转 dq 坐标系，PMSM 的电压方程简化为：

+
+
+ ud = Rs id + Ld + + did + dt + + − ωe Lq iq
+ uq = Rs iq + Lq + + diq + dt + + + ωe (Ld id + ψf) +
+
+
    +
  • ωe = np · ωm：电气角速度
    • +
  • −ωeLqiq：d 轴耦合项（q 轴电流产生的交叉耦合）
    • +
  • +ωeψf：永磁体反电动势（ωψf = Eback-EMF）
    • +
  • 表贴式 PMSM 中 Ld = Lq，无磁阻转矩分量
    • +
+ +

磁场定向控制 (FOC) 原理

+
+

FOC 核心思想：通过控制 id = 0（表贴式），使电磁转矩与 iq 线性正比：

+
+ Te = 1.5 · np · ψf · iq +
+
    +
  • 外环（转速环）→ 给出转矩参考 → 换算为 iq 参考
    • +
  • 内环（电流环）→ PI 控制 id=0, iq=目标值
    • +
  • 逆 Park 变换 → SVPWM 产生三相驱动电压
    • +
+
+ +

坐标变换

+ + + + + + + + + + + + + + + + +
变换作用公式要点
Clarke (ABC→αβ)三相→两相静止坐标iα = ia，iβ = (ia + 2ib)/√3
Park (αβ→dq)静止→旋转坐标id = iαcosθ + iβsinθ
+
+
+
+""" +GPR_KNOWLEDGE = """ +
+

+ 📈 GPR 高斯过程回归 — 理论与公式详解 +

+ +
+ + 📖 什么是高斯过程回归 (Gaussian Process Regression)? + +
+

高斯过程 (Gaussian Process, GP) 是一种定义在函数空间上的随机过程。 + 直觉上，GP 将"先验分布"从有限维的参数空间推广到无穷维的函数空间——它不是对参数建模， + 而是直接对函数本身施加概率分布。

+ +
+

📌 核心定义

+

一个高斯过程是指：对于任意有限输入点集合 {x1, ..., xn}， + 其对应的函数值 [f(x1), ..., f(xn)] 服从联合多元高斯分布。

+
+
+ f(x) ~ GP( m(x),  k(x, x') ) +
+
+
    +
  • 均值函数 m(x) = E[f(x)]：描述函数的先验趋势
    • +
  • 协方差函数（核函数） k(x, x') = Cov[f(x), f(x')]： + 刻画不同输入点处函数值之间的相关性
    • +
+
+ +

与参数化回归的本质区别：

+ + + + + + + + + + + + + + + + + + + + + + + + + + +
方面参数化方法（如线性回归）GPR（非参数贝叶斯）
模型形式f(x) = wTφ(x)，参数 w 有限直接在函数空间上建模
复杂度由模型结构（特征维度）固定随数据量自适应增长
不确定性通常只提供点估计天然提供预测方差（置信区间）
小样本容易过拟合或欠拟合先验正则化，小样本表现好
+
+
+ +
+ + 🔬 核函数 (Kernel / Covariance Function) 详解 + +
+

核函数是 GPR 的"灵魂"——它完全决定了 GP 的先验性质（平滑度、周期性、长程相关性等）。 + 核函数 k(x, x') 衡量两个输入点对应函数值的相关程度。

+ +
+ 1️⃣ RBF 核 (径向基函数核 / 平方指数核) +
+
+ kRBF(x, x') = σf2 · exp + ( + − + + ‖x − x'‖² + 2ℓ² + + ) +
+
+
    +
  • σf2 (输出尺度)：控制函数值的整体变化幅度
    • +
  • ℓ (长度尺度)：控制平滑度——ℓ 越大，函数变化越缓慢；ℓ 越小，函数变化越剧烈
    • +
  • RBF 核假设函数无限可微，适合建模光滑连续的物理过程
    • +
  • 当 ‖x − x'‖ ≫ ℓ 时，k → 0（远距离点不相关）
    • +
+
+ +
+ 2️⃣ Matérn 核（本项目使用 Matérn-5/2） +
+
+ k5/2(r) = σf2 + (1 + √5 · + + r + ℓ + + + + 5r² + 3ℓ² + + ) + exp(−√5 · + + r + ℓ + + ) +
+
+

其中 r = ‖x − x'‖ 为欧氏距离

+
    +
  • Matérn 核是 RBF 的推广，通过参数 ν 控制可微性：ν=1/2 → 不可微; ν=∞ → RBF
    • +
  • ν=5/2：函数两次可微，比 RBF 更灵活，是工程建模的常用选择
    • +
  • 相比 RBF 能更好地捕捉工程数据中的局部不规则性
    • +
+
+ +
+ 3️⃣ ARD (自动相关性判定) 机制 +

对于多维输入 x = [x1, ..., xd]，ARD 核为每个输入维度分配独立的长度尺度：

+
+
+ rARD2 = Σj=1d + + (xj − x'j)² + ℓj² + +
+
+
    +
  • 若 ℓj 很大 → 函数对第 j 维输入不敏感（自动"忽略"该维度）
    • +
  • 若 ℓj 很小 → 函数在第 j 维变化剧烈（高度相关）
    • +
  • ARD 通过数据自动学习每维的重要性，实现隐式特征选择
    • +
  • 本项目：3 个输入 [高度, 马赫数, RPM] → 3 个独立的 ℓj
    • +
+
+
+
+ +
+ + 📐 GPR 后验推断公式（预测均值与方差） + +
+

训练数据：给定 N 个观测 D = {(xi, yi)}i=1N， + 假设 y = f(x) + ε，其中 ε ~ N(0, σn²) 为观测噪声。

+ +

联合先验分布：

+
+
+ [ + + y + f* + + ] + ~ N + ( + + m + m* + + , + + [K + σn²I   K*] + [K*T       K**] + + ) +
+
+ +
+

🎯 后验预测分布（关键公式）

+

对新输入 x*，后验 f* | D, x* 仍为高斯分布：

+
+
+ 预测均值：  + μ* = m(x*) + k*T (K + σn²I)−1 (y − m) +
+
+ 预测方差：  + σ*² = k** − k*T (K + σn²I)−1 k* +
+
+

其中：

+
    +
  • K ∈ ℝN×N：训练数据之间的核矩阵，Kij = k(xi, xj)
    • +
  • k* ∈ ℝN：测试点与所有训练点之间的核向量
    • +
  • k**：测试点自身的先验方差 k(x*, x*)
    • +
  • (K + σn²I)−1：核心计算瓶颈，复杂度 O(N³)
    • +
+
+ +
+

📊 方差的物理含义：

+
    +
  • σ*² 小（接近训练数据区域）→ 预测可信度高，GPR "有信心"
    • +
  • σ*² 大（远离训练数据区域）→ 预测不确定性高，是外推区域
    • +
  • 方差热力图可以直观展示模型"知道什么"和"不知道什么"
    • +
  • 这是 GPR 优于 NN 的核心优势——自动量化认知不确定性
    • +
+
+
+
+ +
+ + ⚙️ 超参数优化 — 最大化边际似然 + +
+

GPR 的超参数 θ = {σf, ℓ1, ..., ℓd, σn} 通过 + 最大化对数边际似然 (Log Marginal Likelihood) 自动确定：

+ +
+
+ log p(y|X, θ) = − + + 1 + 2 + + yTKy−1y + − + + 1 + 2 + + log|Ky| + − + + N + 2 + + log(2π) +
+
+

其中 Ky = K + σn²I

+ +
+

📐 三项分解的物理意义：

+ + + + + + + + + + + + + + + + + + + + + +
项数学形式含义
数据拟合项−½ yTKy−1y模型对数据的拟合程度（越大越好拟合）
复杂度惩罚项−½ log|Ky|Occam 剃刀：自动惩罚过于复杂的模型
归一化常数−(N/2) log(2π)与超参数无关，优化时可忽略
+

+ 优化过程在"拟合数据"和"模型简单性"之间自动取得平衡，天然防止过拟合。 + 本项目使用 BoTorch/GPyTorch 框架，通过 L-BFGS 优化器求解。 +

+
+
+
+ +
+ + 🔄 数据预处理流水线 + +
+

发动机数据（燃油流量、功率）的数值范围跨越多个数量级，直接建模会导致核函数 + 无法有效捕捉小值区域的变化。采用两步预处理：

+ +
+
    +
  1. + Log1p 变换： ylog = log(1 + y)
    + → 压缩大值，拉伸小值，使数据分布更均匀。log1p(0)=0，避免 log(0) 的数值问题。 +
    2. +
  2. + Z-Score 标准化： ŷ = (ylog − μ) / σ
    + → 标准化为均值 0、标准差 1，使核函数的超参数在同一尺度上可比。 +
    3. +
  3. + 反变换（预测时）：
    + ylog = ŷ · σ + μ   →   y = expm1(ylog) = eylog − 1
    + → 从标准化空间恢复到物理空间 +
    4. +
+
+ +
+

💡 相同的预处理同时应用于 GPR 和 NN：

+
    +
  • GPR 在 Z-Score 空间中训练和预测
    • +
  • NN（蒸馏学生模型）使用 GPR 相同的 scaler 参数
    • +
  • 这确保了两个模型的输入输出在相同数值空间中进行比较
    • +
+
+
+
+ +
+ + 📊 模型评价指标与计算复杂度 + +
+

评价指标公式

+
+

MAPE (Mean Absolute Percentage Error)：

+
+ MAPE = + + 1 + N + + Σi=1N + + |yi − ŷi| + yi + + × 100% +
+

注意：MAPE 计算排除 y ≈ 0 的数据点（避免除零导致的虚高误差）

+
+
+

R² (决定系数)：

+
+ R² = 1 − + + Σ(yi − ŷi)² + Σ(yi − ȳ)² + +
+

R² = 1 → 完美预测; R² = 0 → 与均值预测一样差; R² < 0 → 更差

+
+ +

计算复杂度分析

+ + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + +
操作时间复杂度空间复杂度瓶颈
训练（Cholesky 分解）O(N³)O(N²)核矩阵 K 的 Cholesky 分解
单点预测O(N)O(N)需与所有训练点计算核值
批量预测 (M 点)O(MN + N²)O(MN)大 M 时成为瓶颈
超参数优化O(N³) × 迭代次数O(N²)每次迭代都需 Cholesky
+

+ ⚠️ 正是 O(N³) 的训练复杂度和 O(N) 的单点预测成本，使得 GPR 不适合 MPC 实时控制 + （每个控制步需评估模型数十次），因此需要将 GPR 的知识蒸馏到 O(1) 推理的 NN 中。 +

+
+
+ +
+ + ✈️ GPR 在发动机代理建模中的应用 + +
+

在本项目中，GPR 被用作涡轴发动机的高精度代理模型 (Surrogate Model)， + 替代复杂的热力学仿真代码。

+ +
+

建模映射关系：

+
+ [高度(m), 马赫数, RPM]  →  GPR  →  [燃油流量(kg/h), 功率(kW)] +
+ + + + + + + + + + + + + + + + + + + + + +
变量含义典型范围
Altitude_m飞行高度0 ~ 10000 m
Mach飞行马赫数0 ~ 0.8
RPM燃气涡轮转速1000 ~ 30000+
+
+ +
+

💡 为什么选择 GPR 作为代理模型？

+
    +
  • 发动机热力学仿真每次调用耗时秒级，而 GPR 训练后毫秒级预测
    • +
  • 发动机测试数据成本高昂，GPR 在小样本下表现优异
    • +
  • GPR 的方差输出可识别数据稀疏区域，指导后续试验设计
    • +
  • 多输出 GPR (Multi-output GP) 可同时建模燃油流量和功率
    • +
+
+
+
+
+""" + + +# ============================================================ +# NN 神经网络知识蒸馏知识卡片 — 中文 +# ============================================================ +NN_KNOWLEDGE = """ +
+

+ 🧠 NN 神经网络 — 知识蒸馏与 MLP 详解 +

+ +
+ + 📖 知识蒸馏 (Knowledge Distillation) 原理 + +
+

知识蒸馏是一种模型压缩技术， + 将大型"教师"模型（如 GPR）的知识迁移到轻量级"学生"模型（如 MLP 神经网络）中， + 在几乎不损失精度的前提下大幅提升推理速度。

+ +
+
+ [GPR 教师模型] +  → 密集网格预测 →  + [NN 学生模型] +
+
+ +

蒸馏流程

+
+
    +
  1. 训练 GPR 教师：用少量真实发动机数据训练 GPR 模型
    2. +
  2. 生成伪数据：在输入空间的密集网格上用 GPR 预测，生成大量"虚拟标签"
    3. +
  3. 训练 NN 学生：用 GPR 的预测结果作为训练标签，训练 MLP 网络
    4. +
  4. 部署学生网络：在 MPC 控制器中使用轻量 NN 替代 GPR
    5. +
+
+ +

教师 vs 学生对比

+ + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + +
指标GPR (教师)NN (学生)
推理速度~10 ms（核矩阵运算）~0.1 ms（矩阵乘法）
内存占用O(N²) — 存储全部训练数据O(W) — 约 4.5K 参数
不确定性预测均值 + 后验方差仅点估计
MPC 适用性推理太慢，不适合实时✅ 满足实时控制需求
批量预测大批量很慢天然并行，GPU 加速
+
+
+ +
+ + 🏗️ MLP 多层感知机架构详解 + +
+

学生网络采用 3 层全连接 MLP (Multi-Layer Perceptron) 架构：

+
+
+ 输入 (3) + → Linear(64) + Tanh + → Linear(64) + Tanh + → Linear(2) + → 输出 (2) +
+
+ +

前向传播数学公式

+

每一层的计算可以表示为仿射变换 + 非线性激活：

+
+
+ 第 l 层：  + z(l) = W(l) a(l−1) + b(l) +
+
+ 激活后：  + a(l) = σ(z(l)) +
+
+
    +
  • W(l) ∈ ℝdl×dl−1：第 l 层的权重矩阵
    • +
  • b(l) ∈ ℝdl：第 l 层的偏置向量
    • +
  • σ(·)：激活函数（本项目使用 Tanh）
    • +
  • a(0) = x（输入向量）
    • +
+ +

本网络的完整前向计算

+
+
+ h1 = tanh(W1 · x̂ + b1)     // 隐藏层 1: ℝ³ → ℝ⁶⁴
+ h2 = tanh(W2 · h1 + b2)     // 隐藏层 2: ℝ⁶⁴ → ℝ⁶⁴
+ ŷ = W3 · h2 + b3             // 输出层: ℝ⁶⁴ → ℝ² (无激活) +
+
+ +

参数量计算

+ + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + +
层形状权重数偏置数合计
Linear-13 → 6419264256
Linear-264 → 644096644160
Linear-364 → 21282130
总参数量4,546
+
+
+ +
+ + 📐 激活函数详解 + +
+

激活函数为神经网络引入非线性——没有激活函数，无论多少层的网络都等价于单层线性变换。

+ +
+ Tanh (双曲正切) — 本项目所用 +
+
+ tanh(x) = + + ex − e−x + ex + e−x + +
+
+
    +
  • 输出范围：(−1, +1)，以零为中心
    • +
  • 导数：tanh'(x) = 1 − tanh²(x)，最大值为 1（x=0 处）
    • +
  • 优点：输出零中心，有助于加速收敛
    • +
  • 缺点：|x| 较大时梯度接近 0（梯度饱和），deep网络可能出现梯度消失
    • +
+
+ +

常用激活函数对比

+ + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + +
函数公式范围特点
Sigmoid1 / (1 + e−x)(0, 1)输出非零中心，双端饱和
Tanh(ex−e−x)/(ex+e−x)(−1, 1)零中心，双端饱和
ReLUmax(0, x)[0, ∞)计算高效，但有"死神经元"问题
GELUx · Φ(x)连续Transformer 常用，性能好
+

对于本项目的回归任务（光滑连续的发动机特性映射），Tanh 是良好的选择， + 因为数据经过 Z-Score 标准化后值域在 [−3, 3] 左右，刚好在 Tanh 的线性区间内。

+
+
+ +
+ + 📉 损失函数与反向传播 + +
+

损失函数 (Loss Function)

+

本项目在标准化 log 空间中使用 MSE (Mean Squared Error) 损失：

+
+
+ L = + + 1 + N + + Σi=1N ‖ŷi − yi‖² +
+
+

+ 其中 ŷ 和 y 均为标准化 log 空间中的值（先 log1p 变换，再 Z-Score）。 + 在标准化空间中训练使得燃油流量和功率两个输出的损失贡献相当。 +

+ +

反向传播算法 (Backpropagation)

+

反向传播基于链式法则 (Chain Rule)，从输出层到输入层逐层计算梯度：

+
+
+ 输出层梯度：  + + ∂L + ∂W3 + + = + + ∂L + ∂ŷ + + · h2T +
+ 隐藏层梯度：  + + ∂L + ∂Wl + + = δ(l) · a(l−1)T +   其中  δ(l) = (Wl+1T δ(l+1)) ⊙ σ'(z(l)) +
+
+
    +
  • ⊙ 表示逐元素乘法 (Hadamard product)
    • +
  • σ' 为激活函数的导数：tanh'(x) = 1 − tanh²(x)
    • +
  • 反向传播的梯度从后向前"流动"，自动计算所有参数的偏导数
    • +
  • PyTorch 的 autograd 引擎自动构建计算图并执行反向传播
    • +
+
+
+ +
+ + 🔧 优化器与学习率调度 + +
+

Adam 优化器

+

Adam (Adaptive Moment Estimation) 结合了动量 (Momentum) 和 RMSProp 的优点：

+
+
+ mt = β1 mt−1 + (1−β1) gt +   // 一阶矩（梯度均值）
+ vt = β2 vt−1 + (1−β2) gt² +   // 二阶矩（梯度方差）
+ m̂t = mt / (1−β1t) +                  + // 偏差修正
+ v̂t = vt / (1−β2t)
+ θt+1 = θt − η · m̂t / (√v̂t + ε) +
+
+
    +
  • β1=0.9, β2=0.999：动量和二阶矩的衰减率（默认值）
    • +
  • ε=10⁻⁸：防止除零的小常数
    • +
  • η：初始学习率（本项目约 1e-3 ~ 5e-3）
    • +
  • Adam 为每个参数维护独立的自适应学习率
    • +
+ +

余弦退火学习率调度 (Cosine Annealing)

+
+
+ ηt = ηmin + + + ηmax − ηmin + 2 + + ( + 1 + cos + ( + + πt + T + + ) + ) +
+
+
    +
  • 学习率从 ηmax 平滑下降到 ηmin，呈余弦形状
    • +
  • 训练初期学习率大 → 快速探索参数空间
    • +
  • 训练后期学习率小 → 精细化调整，不会跳出最优解
    • +
  • 比阶梯式衰减更平滑，收敛更稳定
    • +
+
+
+ +
+ + 🔄 数据处理流水线与超参数建议 + +
+

完整数据流水线

+
+
    +
  1. 原始输入：[高度(m), 马赫数, 转速(RPM)]
    2. +
  2. 输入 Z-Score：x̂ = (x − μx) / σx
    3. +
  3. NN 前向传播：ŷnorm = MLP(x̂)
    4. +
  4. 反 Z-Score：ŷlog = ŷnorm · σy + μy
    5. +
  5. 反 Log1p：ŷ = expm1(ŷlog) = eŷlog − 1
    6. +
  6. 物理输出：[燃油流量(kg/h), 功率(kW)]
    7. +
+
+ +

超参数调优建议

+ + + + + + + + + + + + + + + + + + + + + + + + + +
超参数推荐范围增大效果减小效果
隐藏层宽度32 ~ 128拟合能力增强，速度稍降推理更快，但可能欠拟合
训练轮数2000 ~ 5000精度更高，但耗时增加训练快，但可能未收敛
学习率1e-3 ~ 5e-3收敛更快，但容易振荡收敛更稳，但可能太慢
+ +
+

🎯 观察训练 Loss 曲线的注意事项：

+
    +
  • Loss 持续下降：正常训练中，一切良好
    • +
  • Loss 振荡不降：学习率可能过大，尝试降低
    • +
  • Loss 很快收敛到平台：模型容量不足（增大隐藏层宽度）或学习率太小
    • +
  • Parity Plot 偏离对角线：模型预测有系统偏差，可能需要更多训练轮数
    • +
  • 典型目标：MAPE < 1%，即 NN 与 GPR 的预测差异小于 1%
    • +
+
+
+
+
+""" + + +# ============================================================ +# 能量管理策略知识卡片 — 中文 +# ============================================================ +EMS_KNOWLEDGE = """ +
+

+ 🔋 串联混动能量管理策略 (EMS) +

+ +
+ + 📊 SOC 迟滞控制 (Hysteresis Control) + +
+

EMS 采用迟滞（滞回）控制策略，在"充电模式"与"功率跟随模式" + 之间引入死区，避免因 SOC 在阈值附近波动而导致频繁切换：

+ +
+
+ SOC < SOClow → 充电模式（发动机以固定高功率运行）
+ SOC > SOChigh → 功率跟随模式（发动机跟踪需求功率）
+ SOClow ≤ SOC ≤ SOChigh → 保持当前模式（迟滞死区） +
+
+ +

运行模式详解

+ + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + +
模式发动机功率触发条件设计目的
紧急充电最大功率 Peng,maxSOC < 10%防止电池深放电损坏
充电模式恒定充电功率 PchargeSOC < SOClow（进入）高效率恒定工况点运行
功率跟随Pdemand + Reserve + KSOC·ΔSOCSOC > SOChigh（退出充电）实时匹配负载需求
过充保护最低功率 Peng,minSOC > 95%防止电池过充损坏
+ +
+

💡 迟滞控制原理：

+
    +
  • 上阈值 SOChigh（如 0.7）：SOC 超过此值时从充电切换到功率跟随
    • +
  • 下阈值 SOClow（如 0.3）：SOC 降到此值时从功率跟随切换到充电
    • +
  • 死区宽度 = SOChigh − SOClow，越宽则切换频率越低，但 SOC 波动范围越大
    • +
  • 死区内保持上一步的运行模式不变 → 有效抑制"抖动"(chattering)
    • +
+
+
+
+ +
+ + 🔌 串联混合动力架构 + +
+
+ + [涡轴发动机] → 发电机 → 直流母线 ← → [锂电池组]
+ 直流母线 → 逆变器 → [PMSM 驱动电机] → 螺旋桨 + +
+

发动机与推进轴机械解耦：所有动力通过电气母线传输， + 允许发动机始终运行在最优效率点，不受负载功率变化的直接影响。

+ +

架构特点

+ + + + + + + + + + + + + + + + + + + + + + + + + + +
特性串联构型并联构型（对比）
发动机-负载耦合完全解耦机械直连
发动机工况可固定在最优效率点随负载波动
EMS 复杂度较低（功率分配）较高（扭矩耦合）
能量传递效率经两次电能转换（略低）机械直驱（略高）
+
+
+ +
+ + ⚙️ 关键参数与调节指南 + +
+ + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + +
参数作用调大效果调小效果
Pcharge充电模式发动机功率SOC 恢复更快充电更温和，效率稍高
KSOC功率跟随模式 SOC 反馈增益SOC 偏差纠正更强发动机功率更平稳
功率裕度 %发动机额外功率余量抗瞬态能力增强油耗更低
电池容量能量缓冲区大小SOC 波动更小系统更轻、响应更快
+ +
+

🎯 仿真观察要点：

+
    +
  • SOC 是否维持在安全范围（10% ~ 95%）内
    • +
  • 充电-跟随模式切换是否过于频繁（理想为低频切换）
    • +
  • 发动机功率曲线是否平稳（频繁剧烈波动 → 油耗升高、寿命下降）
    • +
  • 最终 SOC 与初始 SOC 之差（SOC 平衡性）
    • +
  • 电池充放电电流是否超过安全限值
    • +
+
+ +
+

💡 调参建议：

+
    +
  • 先用默认参数跑一次，观察 SOC 整体走势
    • +
  • 若 SOC 持续下降 → 增大 Pcharge 或 KSOC
    • +
  • 若模式切换过于频繁 → 加大 SOChigh − SOClow 死区宽度
    • +
  • 若发动机功率波动大 → 降低 KSOC，增大功率裕度
    • +
+
+
+
+ +
+ + 🔋 电池 SOC 动态方程与功率平衡 + +
+

SOC 状态方程（库仑计数法）

+
+
+ + dSOC + dt + + = − + + Pbat + Qbat + +
+
+

+ Pbat 为电池充放电功率（放电为正），Qbat 为电池总容量 (kWh)。 + 离散化：SOCk+1 = SOCk − (Pbat,k · Δt) / Qbat +

+ +

系统功率平衡方程

+
+
+ Pdemand(t) = Pengine(t) + Pbattery(t) +
+
+
    +
  • Pdemand > Pengine：电池补充差额（放电），SOC 下降
    • +
  • Pdemand < Pengine：多余功率给电池充电，SOC 上升
    • +
  • EMS 的核心任务就是合理分配 Pengine，使 SOC 维持在安全范围
    • +
+ +

功率跟随模式公式

+
+
+ Peng,cmd = Pdemand × (1 + Reserve%) + KSOC × (SOCtarget − SOC)
+ Peng,cmd = clamp(Peng,cmd, Peng,min, Peng,max) +
+
+

+ 第二项为 SOC 反馈补偿：当 SOC 低于目标时增加发动机功率以充电； + 当 SOC 高于目标时减少发动机功率以节油。 +

+
+
+ +
+ + 📚 进阶：ECMS 等效消耗最小化策略 + +
+

ECMS (Equivalent Consumption Minimization Strategy) 是一种更先进的实时优化 EMS。 + 它将电池充放电等效为燃油消耗，从而将双能源分配问题转化为单目标瞬时优化：

+ +
+
+ minPeng   + J = ṁfuel(Peng) + s(t) · + + Pbat + ηbat · QLHV + +
+
+
    +
  • s(t)：等效因子，将电能消耗折算为燃油消耗
    • +
  • QLHV：燃料低热值 (MJ/kg)
    • +
  • ηbat：电池充放电效率
    • +
  • 当 s(t) 自适应调整以维持 SOC 时，ECMS 近似全局最优
    • +
+ +
+

💡 本平台使用的是基于规则的滞环策略（简单可靠）， + ECMS 和 DP (动态规划) 等高级方法可作为后续扩展方向。

+
+
+
+
+""" diff --git a/case_demo_functions.py b/case_demo_functions.py index 7301670..9a09917 100644 --- a/case_demo_functions.py +++ b/case_demo_functions.py @@ -1,145 +1,1049 @@ import os import sys -# ===== 新增：OpenMP冲突兼容设置，避免PyTorch初始化报错 ===== +# ===== 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.0, 1600.0, 70.0), - (8.0, 3200.0, 220.0), - (20.0, 2500.0, 130.0), - (35.0, 3400.0, 250.0), - (50.0, 1800.0, 80.0), - ] + return [(0., 1600., 70.), (8., 3200., 220.), (20., 2500., 130.), + (35., 3400., 250.), (50., 1800., 80.)] if profile_name == "经济巡航": - return [ - (0.0, 1500.0, 60.0), - (15.0, 2100.0, 95.0), - (35.0, 2300.0, 105.0), - (55.0, 2000.0, 90.0), - ] - return [ - (0.0, 1500.0, 50.0), - (10.0, 3000.0, 200.0), - (30.0, 2800.0, 150.0), - (50.0, 1800.0, 60.0), - ] + 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_scale, load_scale): +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 * rpm_scale, torque * load_scale + return rpm, torque -def run_case_demo(sim_time_s, dt, initial_soc_pct, initial_engine_power_kw, profile_name, rpm_scale, load_scale): +# ============================================================ +# 阶段三：能量管理策略设计 +# ============================================================ +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: - try: - # ===== 新增：懒加载混动系统模型，便于捕获缺失依赖 ===== - from series_hybrid_sim import SeriesHybridSystem - except ModuleNotFoundError as e: - if getattr(e, "name", "") == "torch": - return None, "算例仿真失败：缺少依赖 torch，请先在当前环境安装 PyTorch。", [] - return None, f"算例仿真失败：缺少依赖 {e.name}。", [] + 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.0, 240.0)) + # 参数裁剪 + 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.0, 95.0)) - initial_engine_power_kw = float(np.clip(initial_engine_power_kw, 20.0, 260.0)) - rpm_scale = float(np.clip(rpm_scale, 0.5, 1.6)) - load_scale = float(np.clip(load_scale, 0.5, 1.6)) + 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) - system = SeriesHybridSystem() - system.battery.SOC = initial_soc_pct / 100.0 - system.bus_voltage = system.battery._get_ocv(system.battery.SOC) - system.genset.set_steady_state_by_power(H_env=0.0, Ma_env=0.0, Power_target=initial_engine_power_kw) + time_array = np.arange(0, sim_time_s, dt) - time_array = np.arange(0.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" + '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' ]} - for t in time_array: - target_rpm, load_torque = _target_from_profile(t, points, rpm_scale, load_scale) - res = system.step(dt, target_rpm, load_torque) - res["target_prop_rpm"] = target_rpm - for k in log: - log[k].append(res[k]) + 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) - speed_error = np.array(log["target_prop_rpm"]) - np.array(log["prop_speed_rpm"]) - soc_arr = np.array(log["soc"]) - 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"]) + # 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) - fig, axes = plt.subplots(3, 1, figsize=(12, 10), sharex=True) - axes[0].plot(time_array, log["target_prop_rpm"], "k--", lw=1.5, label="目标转速") - axes[0].plot(time_array, log["prop_speed_rpm"], "b-", lw=1.5, label="实际转速") - axes[0].set_ylabel("RPM") - axes[0].set_title("推进轴转速响应") - axes[0].grid(True, linestyle=":") - axes[0].legend() + p_drive_req = motor_state['p_bus_req_kw'] + actual_rpm = motor_state['n_rpm'] - axes[1].plot(time_array, log["p_drive_req_kw"], "k--", lw=1.2, label="电机需求") - axes[1].plot(time_array, log["p_engine_out_kw"], "r-", lw=1.2, label="发动机输出") - axes[1].plot(time_array, log["p_batt_actual_kw"], "g-", lw=1.2, label="电池功率") - axes[1].axhline(0, color="gray", lw=1) - axes[1].set_ylabel("kW") - axes[1].set_title("功率分配") - axes[1].grid(True, linestyle=":") - axes[1].legend() + # EMS + soc = battery.SOC + if soc < soc_low: charge_mode = True + elif soc > soc_high: charge_mode = False - axes[2].plot(time_array, log["bus_voltage"], "m-", lw=1.2, label="母线电压") - axes[2].set_ylabel("V") - axes[2].set_xlabel("时间 (s)") - axes[2].set_title("电气状态") - axes[2].grid(True, linestyle=":") + 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.set_ylabel("SOC (%)") + 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) - fig.tight_layout() + 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"### 算例结果解读\n" - f"- 仿真时长：{sim_time_s:.1f} s，步长：{dt:.3f} s\n" - f"- 最大转速误差：{np.max(np.abs(speed_error)):.1f} RPM\n" - f"- SOC 变化：{soc_arr[0]:.2f}% → {soc_arr[-1]:.2f}%（最小 {np.min(soc_arr):.2f}%）\n" - f"- 平均发动机输出：{np.mean(engine_pwr_arr):.2f} kW\n" - f"- 平均电池功率：{np.mean(batt_pwr_arr):.2f} kW\n" - f"- 平均燃油流量：{np.mean(fuel_arr):.2f} kg/h" + 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 = [] - for idx in pick_idx: - table_data.append([ - 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), - ]) + 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: - return None, f"算例仿真失败：{e}", [] + import traceback + return None, f"Hybrid simulation failed: {e}\n```\n{traceback.format_exc()}\n```", [] diff --git a/data_usage/usage_stats.json b/data_usage/usage_stats.json index 47c26cf..9e8f16e 100644 --- a/data_usage/usage_stats.json +++ b/data_usage/usage_stats.json @@ -1,4 +1,4 @@ { - "total_users": 22, - "last_saved_at": 1775496284.1673155 + "total_users": 37, + "last_saved_at": 1775560914.9575145 } \ No newline at end of file diff --git a/requirements.txt b/requirements.txt index a36ea6f..d0d6784 100644 --- a/requirements.txt +++ b/requirements.txt @@ -1,7 +1,7 @@ # requirements.txt -# ===== 新增：PyTorch CPU 轮子下载源 ===== ---extra-index-url https://download.pytorch.org/whl/cpu +# ===== PyTorch：GPU (CUDA 12.1) 版本，如无 NVIDIA GPU 可改为 cpu ===== +--extra-index-url https://download.pytorch.org/whl/cu121 gradio==4.44.1 gradio-client==1.3.0 @@ -14,10 +14,11 @@ matplotlib==3.9.4 aiohttp==3.13.5 pillow==10.4.0 # ===== 新增：混动模型（Model）运行依赖 ===== -torch==2.4.1 +torch==2.4.1+cu121 botorch==0.14.0 gpytorch==1.14 pyro-ppl==1.9.1 pandas==2.3.3 -scipy==1.15.3 +scipy>=1.10 scikit-learn==1.7.1 +psutil>=5.9 diff --git a/ui_components.py b/ui_components.py index 25af302..f4d7a96 100644 --- a/ui_components.py +++ b/ui_components.py @@ -1,8 +1,13 @@ -import gradio as gr +import gradio as gr from assets.knowledge_cards_html import ( TIME_DOMAIN_KNOWLEDGE, FREQUENCY_DOMAIN_KNOWLEDGE, - ROOT_LOCUS_KNOWLEDGE + ROOT_LOCUS_KNOWLEDGE, + ENGINE_CONTROL_KNOWLEDGE, + MOTOR_CONTROL_KNOWLEDGE, + GPR_KNOWLEDGE, + NN_KNOWLEDGE, + EMS_KNOWLEDGE, ) def create_header(): @@ -66,6 +71,20 @@ def create_time_domain_tab(): ui_dict = {} with gr.Row(): with gr.Column(scale=1): + with gr.Group(): + gr.HTML("
📊 传递函数设定
") + ui_dict["num_input"] = gr.Textbox( + label="分子系数 (Numerator)", + value="1", + placeholder="例如: 1 或 1,2,3", + info="💡 用逗号分隔，从最高次项到常数项" + ) + ui_dict["den_input"] = gr.Textbox( + label="分母系数 (Denominator)", + value="1,6,11,6", + placeholder="例如: 1,2,1", + info="💡 分母阶数通常 ≥ 分子阶数" + ) with gr.Group(): gr.HTML("
🔧 系统模型
") ui_dict["tf_display"] = gr.Markdown(label="当前传递函数", elem_classes="output-display") @@ -116,6 +135,20 @@ def create_frequency_domain_tab(): ui_dict = {} with gr.Row(): with gr.Column(scale=1): + with gr.Group(): + gr.HTML("
📊 传递函数设定
") + ui_dict["num_input"] = gr.Textbox( + label="分子系数 (Numerator)", + value="1", + placeholder="例如: 1 或 1,2,3", + info="💡 用逗号分隔，从最高次项到常数项" + ) + ui_dict["den_input"] = gr.Textbox( + label="分母系数 (Denominator)", + value="1,6,11,6", + placeholder="例如: 1,2,1", + info="💡 分母阶数通常 ≥ 分子阶数" + ) with gr.Group(): gr.HTML("
🎚️ 调整系统增益
") ui_dict["log_k_slider"] = gr.Slider(minimum=-4, maximum=4, value=1, step=0.01, label="对数增益 log₁₀(K)", info="💡 拖动滑块查看实时变化") @@ -127,6 +160,8 @@ def create_frequency_domain_tab(): gr.HTML("
📊 稳定裕度分析
") ui_dict["metrics_display"] = gr.Textbox(label="Stability Margins", lines=4, interactive=False, elem_classes="output-metrics") ui_dict["stability_display"] = gr.Markdown(elem_classes="stability-result") + with gr.Row(): + ui_dict["analyze_button"] = gr.Button("🚀 开始分析", variant="primary", scale=1, elem_classes="primary-btn") with gr.Column(scale=2): ui_dict["plot_output"] = gr.Plot(label="频域响应图", elem_classes="plot-container") # 知识卡片 @@ -163,6 +198,20 @@ def create_root_locus_tab(): ui_dict = {} with gr.Row(): with gr.Column(scale=1): + with gr.Group(): + gr.HTML("
📊 传递函数设定
") + ui_dict["num_input"] = gr.Textbox( + label="分子系数 (Numerator)", + value="1", + placeholder="例如: 1 或 1,2,3", + info="💡 用逗号分隔，从最高次项到常数项" + ) + ui_dict["den_input"] = gr.Textbox( + label="分母系数 (Denominator)", + value="1,6,11,6", + placeholder="例如: 1,2,1", + info="💡 分母阶数通常 ≥ 分子阶数" + ) with gr.Group(): gr.HTML("
🎚️ 调整系统增益
") ui_dict["log_k_slider"] = gr.Slider(minimum=-4, maximum=4, value=1, step=0.01, label="对数增益 log₁₀(K)", info="💡 拖动滑块观察极点移动") @@ -170,6 +219,8 @@ def create_root_locus_tab(): with gr.Group(): gr.HTML("
📍 闭环极点位置
") ui_dict["poles_display"] = gr.Textbox(label="Closed-Loop Pole Locations", lines=6, interactive=False, elem_classes="output-metrics") + with gr.Row(): + ui_dict["analyze_button"] = gr.Button("🚀 开始分析", variant="primary", scale=1, elem_classes="primary-btn") with gr.Column(scale=2): ui_dict["plot_output"] = gr.Plot(label="根轨迹图") gr.HTML(f""" @@ -201,34 +252,285 @@ def create_root_locus_tab(): return ui_dict def create_case_demo_tab(): - """创建算例演示选项卡的UI组件""" + """创建算例演示选项卡 — 四阶段交互设计（蒸馏→发动机→电机→能量管理）""" ui_dict = {} - with gr.Row(): - with gr.Column(scale=1): - with gr.Group(): - gr.HTML("
🧪 算例参数设置
") - ui_dict["profile"] = gr.Dropdown( - choices=["起飞-巡航-降落", "高机动阶跃", "经济巡航"], - value="起飞-巡航-降落", - label="工况模板" - ) - ui_dict["sim_time"] = gr.Slider(minimum=20, maximum=180, value=60, step=5, label="仿真时长 (s)") - ui_dict["dt"] = gr.Dropdown(choices=[0.02, 0.05, 0.1], value=0.02, label="仿真步长 (s)") - ui_dict["initial_soc"] = gr.Slider(minimum=20, maximum=90, value=60, step=1, label="初始SOC (%)") - ui_dict["initial_engine_power"] = gr.Slider(minimum=20, maximum=250, value=50, step=5, label="初始发动机功率 (kW)") - ui_dict["rpm_scale"] = gr.Slider(minimum=0.6, maximum=1.4, value=1.0, step=0.05, label="目标转速缩放系数") - ui_dict["load_scale"] = gr.Slider(minimum=0.6, maximum=1.4, value=1.0, step=0.05, label="负载转矩缩放系数") - ui_dict["run_button"] = gr.Button("🚀 运行混动算例", variant="primary", elem_classes="primary-btn") - with gr.Group(): - gr.HTML("
📝 结果解读
") - ui_dict["summary"] = gr.Markdown() - with gr.Column(scale=2): - ui_dict["plot"] = gr.Plot(label="混动系统响应图") - ui_dict["table"] = gr.Dataframe( - headers=["时间(s)", "目标转速", "实际转速", "发动机功率(kW)", "电池功率(kW)", "SOC(%)"], - label="关键时刻数据", - interactive=False + + # === MathJax re-render helper (reused across tabs) === + def _mathjax_script(div_id): + return f""" + """ + + with gr.Tabs(): + # ========== 阶段零：模型训练（GPR + NN 两个子标签页）========== + with gr.TabItem("🧬 模型训练", id="distill_tab"): + gr.HTML("""
+ 阶段零：模型训练包含两步——先训练/验证 GPR 高斯过程代理模型， + 再将其知识蒸馏为轻量 NN（MLP）用于后续实时控制仿真。
""") + with gr.Tabs(): + # ----- 子标签页 A: GPR 模型训练 ----- + with gr.TabItem("📈 GPR 模型训练", id="gpr_sub_tab"): + with gr.Row(): + with gr.Column(scale=1): + with gr.Group(): + gr.HTML("
🔬 GPR 训练设置
") + ui_dict["gpr_mode"] = gr.Radio( + choices=["load", "train"], value="load", + label="运行模式", + info="load: 加载已有权重 | train: 从头训练（需 botorch）") + ui_dict["gpr_run_button"] = gr.Button( + "🚀 运行 GPR 训练 / 加载", variant="primary", + elem_classes="primary-btn") + with gr.Group(): + gr.HTML("
📝 GPR 结果
") + ui_dict["gpr_summary"] = gr.Markdown() + with gr.Column(scale=2): + ui_dict["gpr_plot"] = gr.Plot(label="GPR 模型结果") + gr.HTML(f""" +
+ {GPR_KNOWLEDGE} +
+ {_mathjax_script('gpr-knowledge')} + """) + + # ----- 子标签页 B: NN 模型训练 ----- + with gr.TabItem("🧠 NN 模型训练", id="nn_sub_tab"): + with gr.Row(): + with gr.Column(scale=1): + with gr.Group(): + gr.HTML("
🧪 NN 训练参数
") + ui_dict["distill_epochs"] = gr.Slider(minimum=500, maximum=8000, value=3000, step=100, + label="训练轮数 (Epochs)", info="越多越精确，但耗时更长") + ui_dict["distill_lr"] = gr.Slider(minimum=1e-4, maximum=1e-2, value=3e-3, step=1e-4, + label="学习率 (LR)", info="推荐 1e-3 ~ 5e-3") + ui_dict["distill_hidden"] = gr.Slider(minimum=16, maximum=256, value=64, step=16, + label="隐藏层宽度", info="MLP每层神经元数") + ui_dict["distill_run_button"] = gr.Button("🚀 开始 NN 训练", variant="primary", + elem_classes="primary-btn") + with gr.Group(): + gr.HTML("
📝 NN 训练结果
") + ui_dict["distill_summary"] = gr.Markdown() + with gr.Column(scale=2): + ui_dict["distill_plot"] = gr.Plot(label="NN 训练结果 (Loss + Parity)") + gr.HTML(f""" +
+ {NN_KNOWLEDGE} +
+ {_mathjax_script('nn-knowledge')} + """) + + # ========== 阶段一：发动机控制器设计 ========== + with gr.TabItem("🔧 发动机控制器设计", id="engine_tab"): + gr.HTML("""
+ 阶段一：选择 PID 或 MPC 控制器，调整参数，运行阶跃响应测试，观察功率跟踪性能。
""") + with gr.Row(): + with gr.Column(scale=1): + with gr.Group(): + gr.HTML("
🎯 控制器选择
") + ui_dict["eng_controller_type"] = gr.Radio( + choices=["PID", "MPC"], value="PID", label="控制器类型", + info="PID: 经典三参数 | MPC: 模型预测控制") + with gr.Group(visible=True) as eng_pid_group: + gr.HTML("
🎛️ PID 参数
") + ui_dict["eng_kp"] = gr.Slider(minimum=0.1, maximum=20, value=4.652, step=0.01, + label="比例增益 Kp", info="增大加快响应，过大导致振荡") + ui_dict["eng_ki"] = gr.Slider(minimum=0.0, maximum=20, value=7.078, step=0.01, + label="积分增益 Ki", info="消除稳态误差，过大导致超调") + ui_dict["eng_kd"] = gr.Slider(minimum=0.0, maximum=5, value=0.222, step=0.001, + label="微分增益 Kd", info="抑制振荡，改善动态特性") + ui_dict["eng_pid_group"] = eng_pid_group + with gr.Group(visible=False) as eng_mpc_group: + gr.HTML("
🎛️ MPC 参数
") + ui_dict["eng_mpc_horizon"] = gr.Slider(minimum=3, maximum=30, value=15, step=1, + label="预测时域 (Horizon)", info="MPC前看步数") + ui_dict["eng_mpc_W_power"] = gr.Slider(minimum=1, maximum=500, value=200, step=1, + label="功率跟踪权重 W_power") + ui_dict["eng_mpc_W_dcost"] = gr.Slider(minimum=0.01, maximum=20, value=1.5, step=0.01, + label="控制增量权重 W_Δcost") + ui_dict["eng_mpc_overshoot"] = gr.Slider(minimum=1, maximum=30, value=5, step=1, + label="超调限制 (%)") + ui_dict["eng_mpc_group"] = eng_mpc_group + with gr.Group(): + gr.HTML("
⚙️ 发动机模型参数
") + ui_dict["eng_tau_fuel"] = gr.Slider(minimum=0.05, maximum=1.0, value=0.15, step=0.01, + label="燃油执行机构时间常数 τ (s)") + ui_dict["eng_K_inertia"] = gr.Slider(minimum=10, maximum=500, value=100, step=5, + label="转子惯性增益 K") + with gr.Group(): + gr.HTML("
🧪 仿真设置
") + ui_dict["eng_sim_time"] = gr.Slider(minimum=5, maximum=60, value=30, step=1, + label="仿真时长 (s)") + ui_dict["eng_dt"] = gr.Dropdown(choices=[0.02, 0.05, 0.1], value=0.02, + label="仿真步长 (s)") + ui_dict["eng_init_power"] = gr.Slider(minimum=20, maximum=250, value=100, step=5, + label="初始功率 (kW)") + ui_dict["eng_target_power"] = gr.Slider(minimum=20, maximum=300, value=200, step=5, + label="目标功率 (kW)") + ui_dict["eng_run_button"] = gr.Button("🚀 运行发动机仿真", variant="primary", + elem_classes="primary-btn") + with gr.Group(): + gr.HTML("
📝 设计结果
") + ui_dict["eng_summary"] = gr.Markdown() + with gr.Column(scale=2): + ui_dict["eng_plot"] = gr.Plot(label="发动机控制器阶跃响应") + gr.HTML(f""" +
+ {ENGINE_CONTROL_KNOWLEDGE} +
+ {_mathjax_script('engine-knowledge')} + """) + + # Radio toggle PID/MPC visibility + ui_dict["eng_controller_type"].change( + fn=lambda ct: (gr.update(visible=(ct == "PID")), gr.update(visible=(ct == "MPC"))), + inputs=[ui_dict["eng_controller_type"]], + outputs=[eng_pid_group, eng_mpc_group], ) + + # ========== 阶段二：电机控制器设计 ========== + with gr.TabItem("⚡ 电机控制器设计", id="motor_tab"): + gr.HTML("""
+ 阶段二：选择 PID 或 MPC 控制器，运行转速跟踪 + 负载扰动测试。 + 在仿真60%时刻自动施加50%负载扰动，检验抗扰能力。
""") + with gr.Row(): + with gr.Column(scale=1): + with gr.Group(): + gr.HTML("
🎯 控制器选择
") + ui_dict["mot_controller_type"] = gr.Radio( + choices=["PID", "MPC"], value="PID", label="控制器类型", + info="PID: 经典三参数 | MPC: 模型预测控制") + with gr.Group(visible=True) as mot_pid_group: + gr.HTML("
🎛️ PID 参数
") + ui_dict["mot_kp"] = gr.Slider(minimum=0.1, maximum=50, value=5.0, step=0.1, + label="比例增益 Kp", info="增大加快转速响应") + ui_dict["mot_ki"] = gr.Slider(minimum=0.0, maximum=100, value=2.0, step=0.1, + label="积分增益 Ki", info="消除转速稳态偏差") + ui_dict["mot_kd"] = gr.Slider(minimum=0.0, maximum=5, value=0.5, step=0.01, + label="微分增益 Kd", info="抑制转速振荡") + ui_dict["mot_pid_group"] = mot_pid_group + with gr.Group(visible=False) as mot_mpc_group: + gr.HTML("
🎛️ MPC 参数
") + ui_dict["mot_mpc_W_speed"] = gr.Slider(minimum=1, maximum=500, value=200, step=1, + label="转速跟踪权重 W_speed") + ui_dict["mot_mpc_W_dcost"] = gr.Slider(minimum=0.01, maximum=20, value=0.3, step=0.01, + label="控制增量权重 W_Δcost") + ui_dict["mot_mpc_overshoot"] = gr.Slider(minimum=1, maximum=30, value=5, step=1, + label="超调限制 (%)") + ui_dict["mot_mpc_group"] = mot_mpc_group + with gr.Group(): + gr.HTML("
⚙️ 电机模型参数
") + ui_dict["mot_J"] = gr.Slider(minimum=0.1, maximum=5.0, value=1.0, step=0.1, + label="转动惯量 J (kg·m²)", info="越大响应越慢但越平稳") + with gr.Group(): + gr.HTML("
🧪 仿真设置
") + ui_dict["mot_sim_time"] = gr.Slider(minimum=5, maximum=60, value=30, step=1, + label="仿真时长 (s)") + ui_dict["mot_dt"] = gr.Dropdown(choices=[0.02, 0.05, 0.1], value=0.02, + label="仿真步长 (s)") + ui_dict["mot_target_rpm"] = gr.Slider(minimum=500, maximum=5000, value=2000, step=50, + label="目标转速 (RPM)") + ui_dict["mot_load_torque"] = gr.Slider(minimum=10, maximum=400, value=80, step=5, + label="负载转矩 (Nm)") + ui_dict["mot_run_button"] = gr.Button("🚀 运行电机仿真", variant="primary", + elem_classes="primary-btn") + with gr.Group(): + gr.HTML("
📝 设计结果
") + ui_dict["mot_summary"] = gr.Markdown() + with gr.Column(scale=2): + ui_dict["mot_plot"] = gr.Plot(label="电机控制器阶跃响应") + gr.HTML(f""" +
+ {MOTOR_CONTROL_KNOWLEDGE} +
+ {_mathjax_script('motor-knowledge')} + """) + + # Radio toggle PID/MPC visibility + ui_dict["mot_controller_type"].change( + fn=lambda ct: (gr.update(visible=(ct == "PID")), gr.update(visible=(ct == "MPC"))), + inputs=[ui_dict["mot_controller_type"]], + outputs=[mot_pid_group, mot_mpc_group], + ) + + # ========== 阶段三：能量管理策略设计 ========== + with gr.TabItem("🔋 能量管理策略设计", id="ems_tab"): + gr.HTML("""
+ 阶段三：设计基于规则的能量管理策略（自动引用前两阶段的控制器参数）。
+ 策略原理：SOC < 下限阈值 → 进入充电模式; + SOC > 上限阈值 → 退出充电，进入功率跟随模式。 + 下限~上限之间为滞环区间，防止模式频繁切换。
""") + with gr.Row(): + with gr.Column(scale=1): + with gr.Group(): + gr.HTML("
📊 SOC规则参数（滞环控制）
") + ui_dict["soc_target"] = gr.Slider(minimum=30, maximum=80, value=60, step=1, + label="SOC目标值 (%)", info="功率跟随模式下的SOC补偿基准") + ui_dict["soc_low"] = gr.Slider(minimum=15, maximum=50, value=30, step=1, + label="SOC下限阈值 (%)", info="低于此值→进入充电模式") + ui_dict["soc_high"] = gr.Slider(minimum=50, maximum=90, value=70, step=1, + label="SOC上限阈值 (%)", info="高于此值→退出充电模式") + with gr.Group(): + gr.HTML("
⚡ 功率规则参数
") + ui_dict["p_eng_min"] = gr.Slider(minimum=10, maximum=100, value=20, step=5, + label="发动机最小功率 (kW)") + ui_dict["p_eng_max"] = gr.Slider(minimum=100, maximum=350, value=300, step=10, + label="发动机最大功率 (kW)") + ui_dict["p_charge"] = gr.Slider(minimum=50, maximum=300, value=200, step=10, + label="充电模式发动机功率 (kW)", info="进入充电模式后发动机固定输出") + ui_dict["k_soc"] = gr.Slider(minimum=0, maximum=200, value=50, step=5, + label="SOC补偿增益 (kW/ΔSOC)", info="跟随模式下对SOC偏差的修正力度") + ui_dict["power_reserve"] = gr.Slider(minimum=0, maximum=50, value=10, step=1, + label="动态功率储备 (%)", info="发动机额外预留功率百分比") + with gr.Group(): + gr.HTML("
🔧 系统与仿真参数
") + ui_dict["battery_capacity"] = gr.Slider(minimum=10, maximum=200, value=50, step=5, + label="电池容量 (kWh)") + ui_dict["initial_soc"] = gr.Slider(minimum=10, maximum=95, value=60, step=1, + label="初始SOC (%)") + ui_dict["initial_engine_power"] = gr.Slider(minimum=20, maximum=250, value=50, step=5, + label="初始发动机功率 (kW)") + ui_dict["profile"] = gr.Dropdown( + choices=["起飞-巡航-降落", "高机动阶跃", "经济巡航"], + value="起飞-巡航-降落", label="工况模板") + ui_dict["sim_time"] = gr.Slider(minimum=20, maximum=180, value=60, step=5, + label="仿真时长 (s)") + ui_dict["dt"] = gr.Dropdown(choices=[0.02, 0.05, 0.1], value=0.02, + label="仿真步长 (s)") + ui_dict["hybrid_run_button"] = gr.Button("🚀 运行混动系统仿真", variant="primary", + size="lg", elem_classes="primary-btn") + with gr.Group(): + gr.HTML("
📝 结果摘要
") + ui_dict["hybrid_summary"] = gr.Markdown() + with gr.Column(scale=2): + ui_dict["hybrid_plot"] = gr.Plot(label="混动系统响应图") + ui_dict["hybrid_table"] = gr.Dataframe( + headers=["时间(s)", "目标转速", "实际转速", "发动机功率(kW)", + "电池功率(kW)", "SOC(%)", "EMS模式"], + label="关键时刻数据", interactive=False + ) + gr.HTML(f""" +
+ {EMS_KNOWLEDGE} +
+ {_mathjax_script('ems-knowledge')} + """) return ui_dict def create_chatbot_tab(): @@ -260,3 +562,4 @@ def create_chatbot_tab(): label="💡 试试这些问题：" ) return ui_dict +