Hongru/AutoControlCourse
1
0
Fork 0
Code Issues Pull Requests Actions Packages Projects Releases Wiki Activity

增加混电算例功能

Browse Source
This commit is contained in:
Hongru
2026-04-07 19:23:12 +08:00
parent 57ba7b5fee
commit cf9a1a2654
16 changed files with 4338 additions and 838 deletions
Download Patch File Download Diff File Expand all files Collapse all files
Model/data/engine_gpr_model.pth
BIN
View File
Binary file not shown.
Model/data/engine_nn_proxy.pth
BIN
View File
Binary file not shown.
Model/src/distill_gpr_to_nn.py
+275
View File
@@ -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}")
Model/src/engine_dynamic_sim.py
+187 -244
View File
@@ -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')
Model/src/engine_gpr_class.py
+21 -7
View File
@@ -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)
Model/src/lightweight_model.py
+70
View File
@@ -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
Model/src/motor_sim.py
+177 -365
View File
@@ -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 {
Model/src/mpc_controller.py
+298
View File
@@ -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
Model/src/series_hybrid_sim.py
+9 -19
View File
@@ -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]