📚 时域分析常用公式

📐 二阶系统标准形式

标准传递函数：

G(s) =

ω_n²

s² + 2ζω_ns + ω_n²

📊 关键参数

参数	含义	影响
ζ (zeta)	阻尼比	控制超调量和振荡
ω_n	无阻尼自然频率 (rad/s)	决定响应速度

🎯 阻尼比分类

ζ < 0: 不稳定系统
ζ = 0: 无阻尼振荡
0 < ζ < 1: 欠阻尼（有振荡）⭐ 最常见
ζ = 1: 临界阻尼（无超调）
ζ > 1: 过阻尼（响应慢）

📊 时域性能指标（欠阻尼系统）

1️⃣ 上升时间 (Rise Time, t_r)

响应从 10% 上升到 90% 终值所需时间

t_r ≈

1.8

ω_n

2️⃣ 峰值时间 (Peak Time, t_p)

响应达到第一个峰值的时间

t_p =

π

ω_n√(1−ζ²)

=

π

ω_d

其中 ω_d = ω_n√(1−ζ²) 为阻尼振荡频率

3️⃣ 超调量 (Overshoot, σ%)

响应超过稳态值的最大百分比

σ% = e^{−πζ/√(1−ζ²)} × 100%

💡 仅与 ζ 有关！

ζ = 0.5	σ% ≈ 16%
ζ = 0.707	σ% ≈ 4.3% ⭐ 最佳

4️⃣ 调节时间 (Settling Time, t_s)

响应达到并保持在稳态值 ±2%（或±5%）范围内的时间

t_s ≈

4

ζω_n

(2% 误差带)

t_s ≈

3

ζω_n

(5% 误差带)

🎯 一阶系统特性

标准传递函数：

G(s) =

K

τs + 1

阶跃响应：

y(t) = K(1 − e^−t/τ)

💡 关键时间点：

t = τ: 达到稳态值的 63.2%
t = 3τ: 达到稳态值的 95%
t = 4τ: 达到稳态值的 98.2%
t = 5τ: 达到稳态值的 99.3%

📚 频域分析常用方法

🎯 增益裕度 (Gain Margin, GM)

定义：在相角为 -180° 时，系统增益可增加的最大倍数

GM_dB = −20log₁₀|G(jω_pc)|

ω_pc - 相角交越频率（相角 = -180°）

判断准则

GM > 0 dB	✅ 系统稳定
GM = 0 dB	⚠️ 临界稳定
GM < 0 dB	❌ 系统不稳定

💡 工程要求：通常要求 GM ≥ 6 dB（约2倍增益余量）

🎯 相角裕度 (Phase Margin, PM)

定义：在增益为 1 (0dB) 时，系统相角与 -180° 的差值

PM = 180° + ∠G(jω_gc)

ω_gc - 增益交越频率（幅值 = 1 或 0dB）

判断准则

PM > 0°	✅ 系统稳定
PM = 0°	⚠️ 临界稳定
PM < 0°	❌ 系统不稳定

💡 工程要求：通常要求 PM ∈ [30°, 60°]

• PM ≈ 45° ~ 60°: 良好阻尼特性

• PM 越大，超调量越小

与时域性能的关系

ζ ≈

PM

100

(PM 以度为单位)

PM	ζ	超调量
30°	≈ 0.3	≈ 37%
45°	≈ 0.45	≈ 20%
60°	≈ 0.6	≈ 10% ⭐

📊 Bode 图与 Nyquist 图

Bode 图（伯德图）

组成：幅频特性图 + 相频特性图

横轴：频率 ω (对数刻度)
纵轴（幅频）：20log|G(jω)| (dB)
纵轴（相频）：∠G(jω) (度)

💡 优点：

便于绘制（渐近线近似）
直观读取稳定裕度
串联系统可图形叠加

Nyquist 图（奈奎斯特图）

定义：开环频率特性在复平面上的轨迹

🎯 奈奎斯特稳定判据：

Z = P − N

• Z: 闭环右半平面极点数
• P: 开环右半平面极点数
• N: 曲线逆时针包围 (−1, j0) 的圈数

稳定条件：Z = 0

🔧 典型环节频率特性

1️⃣ 比例环节 K

• 幅频：20logK (dB) - 水平线
• 相频：0° - 水平线

2️⃣ 积分环节 1/s

• 幅频：-20dB/dec 斜率
• 相频：-90° - 水平线

3️⃣ 惯性环节 1/(Ts+1)

• 转折频率：ω = 1/T
• 幅频：低频 0dB，高频 -20dB/dec
• 相频：0° → -90°

📚 根轨迹分析要点

🎯 根轨迹基本概念

定义：当开环增益 K 从 0 → ∞ 变化时，闭环特征方程根（极点）在 s 平面上的运动轨迹

特征方程

1 + KG(s)H(s) = 0

主要作用

✅ 直观显示增益对极点位置的影响
✅ 判断系统稳定性
✅ 选择合适的增益值
✅ 设计控制器参数

📐 根轨迹绘制基本条件

1️⃣ 幅值条件（充要条件）

|KG(s₀)H(s₀)| = 1

💡 用途：确定增益 K 的值

2️⃣ 相角条件（充要条件）

∠G(s₀)H(s₀) = (2k+1) × 180°

其中 k = 0, ±1, ±2, ...

💡 用途：判断 s 平面某点是否在根轨迹上

🌟 根轨迹基本性质

规则	内容
起点	K=0，位于开环极点
终点	K→∞，位于开环零点或无穷远
分支数	max(n, m)，n=极点数，m=零点数
实轴段	右侧零极点总数为奇数的区段
渐近线数	n − m 条
渐近线角度	φ_a = (2k+1)×180°/(n−m)
渐近线交点	σ_a = (Σp − Σz)/(n−m)

🛡️ s 平面稳定性区域

区域	条件	稳定性
左半平面	Re(s) < 0	✅ 稳定
虚轴	Re(s) = 0	⚠️ 临界稳定
右半平面	Re(s) > 0	❌ 不稳定

💡 虚轴（x=0）即为稳定性边界

根轨迹穿越虚轴点对应临界增益 K_c

📏 阻尼比等值线

从原点出发的射线代表恒定阻尼比 ζ 的轨迹：

θ = arccos(ζ)

阻尼比 ζ	角度 θ	系统特性
0.5	60°	欠阻尼，较大超调
0.707	45°	⭐ 最佳阻尼

💡 工程应用：通过根轨迹与阻尼比等值线的交点，选择满足性能要求的增益 K

🔧 发动机控制器设计理论

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

标准PID传递函数：

u(t) = K_p e(t) + K_i ∫e(τ)dτ + K_d de(t) dt

增量式（速度型）PID：

Δu(k) = K_p[e(k) − e(k−1)] + K_i e(k) + K_d[e(k) − 2e(k−1) + e(k−2)]

💡 增量式优势：

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

📊 各参数作用与调节规律

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

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

设 K_i=0, K_d=0，逐渐增大 K_p 直到系统出现等幅振荡
记录此时的临界增益 K_u 和振荡周期 T_u
按经验公式设定：K_p=0.6K_u, K_i=2K_p/T_u, K_d=K_pT_u/8
根据实际响应微调各参数

🧠 MPC — 模型预测控制

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

🔑 优化目标函数

min_U J = Σ_k=1^H [ W_p (P_k − P_ref)² + W_Δu (Δu_k)² ]

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

🔑 MPC 核心概念

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

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

N_k+1 = N_k + Δt · K_inertia · (W_f,act − W_f,req)
W_f,act,k+1 = (1 − Δt/τ) · W_f,act,k + (Δt/τ) · u_k

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

⚖️ PID vs MPC 对比：

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

✈️ 涡轴发动机模型

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

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

💡 仿真流程：

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

📏 性能评价指标：

稳态误差：功率达到稳态后与目标的平均偏差
超调量 (σ%)：输出超过目标值的最大百分比
上升时间 (t_r)：功率从10%变化到90%所需时间
调节时间 (t_s)：功率进入并持续保持在 ±2% 误差带内的时间

📐 PID 控制器 — Laplace 域传递函数

连续域 PID 传递函数：

C(s) = K_p ( 1 + 1 T_i s + T_d s )

其中 T_i = K_p/K_i 为积分时间常数， T_d = K_d/K_p 为微分时间常数。

Ziegler-Nichols 整定公式表

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

控制器类型	K_p	T_i	T_d
P	0.5 K_u	—	—
PI	0.45 K_u	T_u / 1.2	—
PID	0.6 K_u	T_u / 2	T_u / 8

等价转换：K_i = K_p/T_i，K_d = K_p·T_d。 Z-N 公式通常产生约 25% 超调，实际应用中需进一步微调。

🧮 MPC 矩阵化优化问题

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

离散状态空间模型

x_k+1 = Ax_k + Bu_k
y_k = Cx_k

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

预测矩阵展开

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

Y = Ψx₀ + ΘU

Y = [y₁, ..., y_H]^T：预测输出序列
U = [u₀, ..., u_H−1]^T：控制输入序列
Ψ = [CA, CA², ..., CA^H]^T：自由响应矩阵
Θ：Toeplitz 控制响应矩阵

带约束 QP 问题

min_U (Y − R)^T Q (Y − R) + ΔU^T R_Δ ΔU

s.t. u_min ≤ u_k ≤ u_max, |y_k − r_k| ≤ ε_overshoot

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

⚡ 电机控制器设计理论

⚙️ PMSM 永磁同步电机模型

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

机械动力学方程

J dω dt = −(T_motor + T_load + T_ext + Bω)

其中 T_motor 为电磁转矩（电机产生），T_load 为负载转矩， T_ext 为外部扰动，Bω 为粘性摩擦力矩。

电磁转矩表达式

T_e = 1.5 · n_p · ψ_f · i_q

对表贴式PMSM (L_d = L_q)，转矩仅与q轴电流 i_q 成正比。

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

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

🎯 转速控制环路

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

ω_ref → [转速控制器] → T_cmd → [电流环 + SVPWM] → PMSM → ω_actual

PID 模式：

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

MPC 模式：

利用电机离散动力学模型预测 H 步后的转速轨迹
优化 H 步的转矩序列，使转速跟踪误差和控制增量同时最小化
支持超调硬约束：转速偏差不得超过设定百分比

min_{T_0..H} Σ [ W_ω(ω_k − ω_ref)² + W_Δ(T_k − T_k−1)² ]

离散状态方程：ω_k+1 = (1−dt·B/J)·ω_k − (dt/J)·T_cmd − (dt/J)·(T_ext+T_load)

🔍 负载扰动测试设计

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

📊 抗扰性能指标

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

⚙️ 转动惯量 J 的影响：

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

💡 MPC 调参建议：

如果扰动后转速恢复慢→增大 W_speed 或减小 W_Δcost
如果转速振荡→增大 W_Δcost 或减小 W_speed
如果超调过大→降低超调限制百分比
负载很大时，确认电机额定功率是否足以提供所需转矩

📐 PMSM dq 坐标系电压方程

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

u_d = R_s i_d + L_d di_d dt − ω_e L_q i_q
u_q = R_s i_q + L_q di_q dt + ω_e (L_d i_d + ψ_f)

ω_e = n_p · ω_m：电气角速度
−ω_eL_qi_q：d 轴耦合项（q 轴电流产生的交叉耦合）
+ω_eψ_f：永磁体反电动势（ωψ_f = E_back-EMF）
表贴式 PMSM 中 L_d = L_q，无磁阻转矩分量

磁场定向控制 (FOC) 原理

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

T_e = 1.5 · n_p · ψ_f · i_q

外环（转速环）→ 给出转矩参考 → 换算为 i_q 参考
内环（电流环）→ PI 控制 i_d=0, i_q=目标值
逆 Park 变换 → SVPWM 产生三相驱动电压

坐标变换

变换	作用	公式要点
Clarke (ABC→αβ)	三相→两相静止坐标	i_α = i_a，i_β = (i_a + 2i_b)/√3
Park (αβ→dq)	静止→旋转坐标	i_d = i_αcosθ + i_βsinθ

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

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

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

📌 核心定义

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

f(x) ~ GP( m(x), k(x, x') )

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

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

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

🔬 核函数 (Kernel / Covariance Function) 详解

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

1️⃣ RBF 核 (径向基函数核 / 平方指数核)

k_RBF(x, x') = σ_f² · exp ( − ‖x − x'‖² 2ℓ² )

σ_f² (输出尺度)：控制函数值的整体变化幅度
ℓ (长度尺度)：控制平滑度——ℓ 越大，函数变化越缓慢；ℓ 越小，函数变化越剧烈
RBF 核假设函数无限可微，适合建模光滑连续的物理过程
当 ‖x − x'‖ ≫ ℓ 时，k → 0（远距离点不相关）

2️⃣ Matérn 核（本项目使用 Matérn-5/2）

k_5/2(r) = σ_f² (1 + √5 · r ℓ + 5r² 3ℓ² ) exp(−√5 · r ℓ )

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

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

3️⃣ ARD (自动相关性判定) 机制

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

r_ARD² = Σ_j=1^d (x_j − x'_j)² ℓ_j²

若 ℓ_j 很大 → 函数对第 j 维输入不敏感（自动"忽略"该维度）
若 ℓ_j 很小 → 函数在第 j 维变化剧烈（高度相关）
ARD 通过数据自动学习每维的重要性，实现隐式特征选择
本项目：3 个输入 [高度, 马赫数, RPM] → 3 个独立的 ℓ_j

📐 GPR 后验推断公式（预测均值与方差）

训练数据：给定 N 个观测 D = {(x_i, y_i)}_i=1^N，假设 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)⁻¹ (y − m)

预测方差： σ_*² = k_** − k_*^T (K + σ_n²I)⁻¹ k_*

其中：

K ∈ ℝ^N×N：训练数据之间的核矩阵，K_ij = k(x_i, x_j)
k_* ∈ ℝ^N：测试点与所有训练点之间的核向量
k_**：测试点自身的先验方差 k(x_*, x_*)
(K + σ_n²I)⁻¹：核心计算瓶颈，复杂度 O(N³)

📊 方差的物理含义：

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

⚙️ 超参数优化 — 最大化边际似然

GPR 的超参数 θ = {σ_f, ℓ₁, ..., ℓ_d, σ_n} 通过 最大化对数边际似然 (Log Marginal Likelihood) 自动确定：

log p(y|X, θ) = − 1 2 y^TK_y⁻¹y − 1 2 log|K_y| − N 2 log(2π)

其中 K_y = K + σ_n²I

📐 三项分解的物理意义：

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

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

🔄 数据预处理流水线

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

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

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

GPR 在 Z-Score 空间中训练和预测
NN（蒸馏学生模型）使用 GPR 相同的 scaler 参数
这确保了两个模型的输入输出在相同数值空间中进行比较

📊 模型评价指标与计算复杂度

评价指标公式

MAPE (Mean Absolute Percentage Error)：

MAPE = 1 N Σ_i=1^N |y_i − ŷ_i| y_i × 100%

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

R² (决定系数)：

R² = 1 − Σ(y_i − ŷ_i)² Σ(y_i − ȳ)²

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 神经网络 — 知识蒸馏与 MLP 详解

📖 知识蒸馏 (Knowledge Distillation) 原理

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

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

蒸馏流程

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

教师 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) ∈ ℝ^d_l×d_l−1：第 l 层的权重矩阵
b^(l) ∈ ℝ^d_l：第 l 层的偏置向量
σ(·)：激活函数（本项目使用 Tanh）
a⁽⁰⁾ = x（输入向量）

本网络的完整前向计算

h₁ = tanh(W₁ · x̂ + b₁)     // 隐藏层 1: ℝ³ → ℝ⁶⁴
h₂ = tanh(W₂ · h₁ + b₂)     // 隐藏层 2: ℝ⁶⁴ → ℝ⁶⁴
ŷ = W₃ · h₂ + b₃             // 输出层: ℝ⁶⁴ → ℝ² (无激活)

参数量计算

层	形状	权重数	偏置数	合计
Linear-1	3 → 64	192	64	256
Linear-2	64 → 64	4096	64	4160
Linear-3	64 → 2	128	2	130
总参数量				4,546

📐 激活函数详解

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

Tanh (双曲正切) — 本项目所用

tanh(x) = e^x − e^−x e^x + e^−x

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

常用激活函数对比

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

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

📉 损失函数与反向传播

损失函数 (Loss Function)

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

L = 1 N Σ_i=1^N ‖ŷ_i − y_i‖²

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

反向传播算法 (Backpropagation)

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

输出层梯度： ∂L ∂W₃ = ∂L ∂ŷ · h₂^T
隐藏层梯度： ∂L ∂W_l = δ^(l) · a^(l−1)T 其中 δ^(l) = (W_l+1^T δ^(l+1)) ⊙ σ'(z^(l))

⊙ 表示逐元素乘法 (Hadamard product)
σ' 为激活函数的导数：tanh'(x) = 1 − tanh²(x)
反向传播的梯度从后向前"流动"，自动计算所有参数的偏导数
PyTorch 的 autograd 引擎自动构建计算图并执行反向传播

🔧 优化器与学习率调度

Adam 优化器

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

m_t = β₁ m_t−1 + (1−β₁) g_t   // 一阶矩（梯度均值）
v_t = β₂ v_t−1 + (1−β₂) g_t²   // 二阶矩（梯度方差）
m̂_t = m_t / (1−β₁^t)                  // 偏差修正
v̂_t = v_t / (1−β₂^t)
θ_t+1 = θ_t − η · m̂_t / (√v̂_t + ε)

β₁=0.9, β₂=0.999：动量和二阶矩的衰减率（默认值）
ε=10⁻⁸：防止除零的小常数
η：初始学习率（本项目约 1e-3 ~ 5e-3）
Adam 为每个参数维护独立的自适应学习率

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

η_t = η_min + η_max − η_min 2 ( 1 + cos ( πt T ) )

学习率从 η_max 平滑下降到 η_min，呈余弦形状
训练初期学习率大 → 快速探索参数空间
训练后期学习率小 → 精细化调整，不会跳出最优解
比阶梯式衰减更平滑，收敛更稳定

🔄 数据处理流水线与超参数建议

完整数据流水线

原始输入：[高度(m), 马赫数, 转速(RPM)]
输入 Z-Score：x̂ = (x − μ_x) / σ_x
NN 前向传播：ŷ_norm = MLP(x̂)
反 Z-Score：ŷ_log = ŷ_norm · σ_y + μ_y
反 Log1p：ŷ = expm1(ŷ_log) = e^ŷ_log − 1
物理输出：[燃油流量(kg/h), 功率(kW)]

超参数调优建议

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

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

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

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

📊 SOC 迟滞控制 (Hysteresis Control)

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

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

运行模式详解

模式	发动机功率	触发条件	设计目的
紧急充电	最大功率 P_eng,max	SOC < 10%	防止电池深放电损坏
充电模式	恒定充电功率 P_charge	SOC < SOC_low（进入）	高效率恒定工况点运行
功率跟随	P_demand + Reserve + K_SOC·ΔSOC	SOC > SOC_high（退出充电）	实时匹配负载需求
过充保护	最低功率 P_eng,min	SOC > 95%	防止电池过充损坏

💡 迟滞控制原理：

上阈值 SOC_high（如 0.7）：SOC 超过此值时从充电切换到功率跟随
下阈值 SOC_low（如 0.3）：SOC 降到此值时从功率跟随切换到充电
死区宽度 = SOC_high − SOC_low，越宽则切换频率越低，但 SOC 波动范围越大
死区内保持上一步的运行模式不变 → 有效抑制"抖动"(chattering)

🔌 串联混合动力架构

[涡轴发动机] → 发电机 → 直流母线 ← → [锂电池组]
直流母线 → 逆变器 → [PMSM 驱动电机] → 螺旋桨

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

架构特点

特性	串联构型	并联构型（对比）
发动机-负载耦合	完全解耦	机械直连
发动机工况	可固定在最优效率点	随负载波动
EMS 复杂度	较低（功率分配）	较高（扭矩耦合）
能量传递效率	经两次电能转换（略低）	机械直驱（略高）

⚙️ 关键参数与调节指南

参数	作用	调大效果	调小效果
P_charge	充电模式发动机功率	SOC 恢复更快	充电更温和，效率稍高
K_SOC	功率跟随模式 SOC 反馈增益	SOC 偏差纠正更强	发动机功率更平稳
功率裕度 %	发动机额外功率余量	抗瞬态能力增强	油耗更低
电池容量	能量缓冲区大小	SOC 波动更小	系统更轻、响应更快

🎯 仿真观察要点：

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

💡 调参建议：

先用默认参数跑一次，观察 SOC 整体走势
若 SOC 持续下降 → 增大 P_charge 或 K_SOC
若模式切换过于频繁 → 加大 SOC_high − SOC_low 死区宽度
若发动机功率波动大 → 降低 K_SOC，增大功率裕度

🔋 电池 SOC 动态方程与功率平衡

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

dSOC dt = − P_bat Q_bat

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

系统功率平衡方程

P_demand(t) = P_engine(t) + P_battery(t)

P_demand > P_engine：电池补充差额（放电），SOC 下降
P_demand < P_engine：多余功率给电池充电，SOC 上升
EMS 的核心任务就是合理分配 P_engine，使 SOC 维持在安全范围

功率跟随模式公式

P_eng,cmd = P_demand × (1 + Reserve%) + K_SOC × (SOC_target − SOC)
P_eng,cmd = clamp(P_eng,cmd, P_eng,min, P_eng,max)

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

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

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

min_{P_eng} J = ṁ_fuel(P_eng) + s(t) · P_bat η_bat · Q_LHV

s(t)：等效因子，将电能消耗折算为燃油消耗
Q_LHV：燃料低热值 (MJ/kg)
η_bat：电池充放电效率
当 s(t) 自适应调整以维持 SOC 时，ECMS 近似全局最优

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