commit 2757c8afd4486c0e53f6d5c06612acf7c2d9f4b2 Author: Hongru Date: Tue Dec 30 19:10:49 2025 +0800 Initial commit diff --git a/.gitignore b/.gitignore new file mode 100644 index 0000000..74ec90c --- /dev/null +++ b/.gitignore @@ -0,0 +1,126 @@ +# ======================================================= +# Python 通用忽略规则 +# ======================================================= + +# 字节码编译文件 +__pycache__/ +*.py[cod] +*$py.class + +# C 扩展 +*.so + +# 分发与打包 +.Python +build/ +develop-eggs/ +dist/ +downloads/ +eggs/ +.eggs/ +lib/ +lib64/ +parts/ +sdist/ +var/ +wheels/ +share/python-wheels/ +*.egg-info/ +.installed.cfg +*.egg +MANIFEST + +# PyInstaller +*.manifest +*.spec + +# 安装日志 +pip-log.txt +pip-delete-this-directory.txt + +# 单元测试 / 覆盖率 +htmlcov/ +.tox/ +.nox/ +.coverage +.coverage.* +.cache +nosetests.xml +coverage.xml +*.cover +*.py,cover +.hypothesis/ +.pytest_cache/ +cover/ + +# Jupyter Notebook +.ipynb_checkpoints + +# 环境配置 +.env +.venv +env/ +venv/ +ENV/ +env.bak/ +venv.bak/ + +# ======================================================= +# IDE 与 编辑器 +# ======================================================= +.vscode/ +.idea/ +*.swp +*.swo + +# ======================================================= +# 操作系统生成文件 +# ======================================================= +Thumbs.db +ehthumbs.db +Desktop.ini +.DS_Store + +# ======================================================= +# 强化学习 (RL) 与 深度学习特定忽略规则 +# ======================================================= + +# 训练日志与 TensorBoard +logs/ +runs/ +tensorboard/ +wandb/ +mlruns/ + +# 模型权重与检查点 (通常很大，不建议提交) +models/ +checkpoints/ +saved_models/ +*.pth +*.pt +*.ckpt +*.h5 +*.onnx +*.pkl +*.pb + +# 数据集与回放缓冲区 +data/ +datasets/ +replay_buffer/ +*.npz + +# 仿真输出数据 +results/ +*.csv +*.xlsx + +# ======================================================= +# 项目特定 +# ======================================================= + +# *.dll + +# 临时文件 +*.tmp +*.bak diff --git a/engine_env/DLL_Usage_Guide.md b/engine_env/DLL_Usage_Guide.md new file mode 100644 index 0000000..a37e8fb --- /dev/null +++ b/engine_env/DLL_Usage_Guide.md @@ -0,0 +1,227 @@ +# libEngine.dll 调用指南 + +本文档详细说明了如何通过 C++ 和 Python 调用 `libEngine.dll` 进行发动机仿真。 + +## 1. 数据结构定义 + +无论使用哪种语言，都需要严格遵循以下数据结构定义（内存布局）。 + +### 1.1 Engine_Identity (身份/配置参数) + +| 字段名 | 类型 (C++) | 类型 (Python ctypes) | 说明 | +| :------------------------ | :--------- | :------------------- | :------------------------- | +| `Identity_Ok` | `int` | `c_int` | 身份验证位，通常设为 1 | +| `Cof_Eff_Low_Identity` | `double` | `c_double` | 低压部件效率系数，默认 1.0 | +| `Cof_Eff_High_Identity` | `double` | `c_double` | 高压部件效率系数，默认 1.0 | + +### 1.2 EngInPut (输入控制量) + +| 字段名 | 类型 (C++) | 类型 (Python ctypes) | 说明 | +| :----------------- | :--------- | :------------------- | :--------------------- | +| `Altp` | `double` | `c_double` | 飞行高度 (m) | +| `Ma0` | `double` | `c_double` | 飞行马赫数 | +| `dT0` | `double` | `c_double` | 与标准大气温差 (K) | +| `StepTime` | `double` | `c_double` | 仿真步长 (s) | +| `Wf` | `double` | `c_double` | 主燃烧室供油量 (kg/h) | +| `Wf_After` | `double` | `c_double` | 加力供油量 (kg/h) | +| `Angle_FanVane` | `double` | `c_double` | 风扇导流叶片角度 | +| `Angle_CompVane` | `double` | `c_double` | 高压压气机导流叶片角度 | +| `A8` | `double` | `c_double` | 喷口临界面积 (m^2) | +| `AddPower` | `double` | `c_double` | 附加功率/起动功率 (W) | + +### 1.3 EngOutPut (输出状态量) + +| 字段名 | 类型 (C++) | 类型 (Python ctypes) | 说明 | +| :---------- | :--------- | :------------------- | :-------------------------- | +| `NL` | `double` | `c_double` | 风扇相对转速 (%) | +| `NH` | `double` | `c_double` | 高压相对转速 (%) | +| `T1t` | `double` | `c_double` | 进气温度 (K) | +| `P1t` | `double` | `c_double` | 进气总压 (kPa) | +| `P1s` | `double` | `c_double` | 风扇进口静压 (kPa) | +| `P3s` | `double` | `c_double` | 高压压气机后静压 (kPa) | +| `T5t` | `double` | `c_double` | 涡轮后温度 (K) | +| `P5t` | `double` | `c_double` | 涡轮后压力 (kPa) | +| `Wf_Main` | `double` | `c_double` | 实际燃烧主燃油流量 (kg/h) | +| `T5t_Gas` | `double` | `c_double` | 涡轮后气体温度 (无惯性) (K) | + +--- + +## 2. C++ 调用方法 + +在 C++ 中，通常使用 `LoadLibrary` 和 `GetProcAddress` 进行显式调用（动态加载）。 + +### 函数原型 + +```cpp +// 假设使用 stdcall 调用约定 +typedef struct EngOutPut (__stdcall *CreateEngFunc)(int **pEngine, struct Engine_Identity m_Identity); +typedef struct EngOutPut (__stdcall *EngStepGoFunc)(int *pEngine, struct EngInPut m_EngInPut); +typedef struct EngOutPut (__stdcall *DestroyEngFunc)(int *pEngine); +``` + +### 示例代码 + +```cpp +#include +#include + +// 定义结构体 (需与上述定义一致) +struct Engine_Identity { + int Identity_Ok; + double Cof_Eff_Low_Identity; + double Cof_Eff_High_Identity; +}; + +struct EngInPut { + double Altp, Ma0, dT0, StepTime; + double Wf, Wf_After, Angle_FanVane, Angle_CompVane, A8, AddPower; +}; + +struct EngOutPut { + double NL, NH, T1t, P1t, P1s, P3s, T5t, P5t, Wf_Main, T5t_Gas; +}; + +// 定义函数指针类型 +typedef EngOutPut (__stdcall *CreateEngFunc)(int**, Engine_Identity); +typedef EngOutPut (__stdcall *EngStepGoFunc)(int*, EngInPut); +typedef EngOutPut (__stdcall *DestroyEngFunc)(int*); + +int main() { + // 1. 加载 DLL + HMODULE hDll = LoadLibrary("libEngine.dll"); + if (!hDll) { + std::cerr << "无法加载 DLL" << std::endl; + return 1; + } + + // 2. 获取函数地址 + CreateEngFunc CreateEng = (CreateEngFunc)GetProcAddress(hDll, "CreateEng"); + EngStepGoFunc EngStepGo = (EngStepGoFunc)GetProcAddress(hDll, "EngStepGo"); + DestroyEngFunc DestroyEng = (DestroyEngFunc)GetProcAddress(hDll, "DestroyEng"); + + if (!CreateEng || !EngStepGo || !DestroyEng) { + std::cerr << "无法获取函数地址" << std::endl; + FreeLibrary(hDll); + return 1; + } + + // 3. 创建发动机实例 + int* hEngine = nullptr; // 句柄指针 + Engine_Identity identity = {1, 1.0, 1.0}; + + // 注意：CreateEng 需要传入指针的地址 (&hEngine) + CreateEng(&hEngine, identity); + + if (!hEngine) { + std::cerr << "发动机创建失败" << std::endl; + FreeLibrary(hDll); + return 1; + } + + // 4. 仿真循环 + EngInPut input = {0}; + input.StepTime = 0.02; + input.Wf = 100.0; + input.A8 = 0.1; + input.AddPower = 90000.0; // 起动功率 + + for (int i = 0; i < 100; ++i) { + // 执行单步 + EngOutPut output = EngStepGo(hEngine, input); + + std::cout << "Step: " << i + << " NH: " << output.NH + << " T5t: " << output.T5t << std::endl; + + // 简单的起动逻辑示例 + if (output.NH > 0.25) { + input.AddPower = 0.0; + input.Wf = 300.0; + } + } + + // 5. 销毁与释放 + DestroyEng(hEngine); + FreeLibrary(hDll); + + return 0; +} +``` + +--- + +## 3. Python 调用方法 + +Python 中推荐使用 `ctypes` 库进行调用。 + +### 核心要点 + +1. 使用 `ctypes.Structure` 定义对应的 C 结构体。 +2. 使用 `WinDLL` 加载 DLL（因为是 `__stdcall` 调用约定）。 +3. 配置 `argtypes` 和 `restype` 以确保参数传递正确。 + +### 示例代码 (基于封装好的类) + +```python +import ctypes +from ctypes import * +import os + +# --- 结构体定义 (略，见 core_model.py) --- + +class EngineSim: + def __init__(self, dll_path): + self.lib = WinDLL(dll_path) + + # 配置函数原型 + # CreateEng: 传入 int** (POINTER(POINTER(c_int))) + self.lib.CreateEng.argtypes = [POINTER(POINTER(c_int)), Engine_Identity] + self.lib.CreateEng.restype = EngOutPut + + # EngStepGo: 传入 int* (POINTER(c_int)) + self.lib.EngStepGo.argtypes = [POINTER(c_int), EngInPut] + self.lib.EngStepGo.restype = EngOutPut + + # DestroyEng: 传入 int* + self.lib.DestroyEng.argtypes = [POINTER(c_int)] + self.lib.DestroyEng.restype = EngOutPut + + self.h_engine = None + + def create(self): + self.h_engine = POINTER(c_int)() # 创建一个空指针用于接收句柄 + identity = Engine_Identity(1, 1.0, 1.0) + # 传入指针的引用 byref + self.lib.CreateEng(byref(self.h_engine), identity) + + def step(self, input_data): + return self.lib.EngStepGo(self.h_engine, input_data) + + def close(self): + if self.h_engine: + self.lib.DestroyEng(self.h_engine) + self.h_engine = None + +# --- 使用 --- +if __name__ == "__main__": + dll_path = os.path.join(os.path.dirname(__file__), "libEngine.dll") + sim = EngineSim(dll_path) + sim.create() + + inp = EngInPut() + inp.StepTime = 0.02 + inp.Wf = 100.0 + inp.AddPower = 90000.0 + inp.A8 = 0.1 + + for i in range(100): + out = sim.step(inp) + print(f"NH: {out.NH:.2f}, T5t: {out.T5t:.2f}") + + sim.close() +``` + +### 注意事项 + +* **指针传递**: `CreateEng` 在 C++ 中接收 `int**`，在 Python `ctypes` 中对应 `byref(h_engine)`，其中 `h_engine` 是 `POINTER(c_int)()`。 +* **调用约定**: 必须使用 `WinDLL` 而不是 `CDLL`，除非 DLL 编译时使用的是 `cdecl`。根据现有代码推断为 `stdcall`。 diff --git a/engine_env/core_model.py b/engine_env/core_model.py new file mode 100644 index 0000000..67959ea --- /dev/null +++ b/engine_env/core_model.py @@ -0,0 +1,128 @@ +import ctypes +import os +from ctypes import * + +# ========================================== +# 1. 结构体定义 +# ========================================== + +class Engine_Identity(Structure): + _fields_ = [ + ("Identity_Ok", c_int), + ("Cof_Eff_Low_Identity", c_double), + ("Cof_Eff_High_Identity", c_double) + ] + +class EngInPut(Structure): + _fields_ = [ + ("Altp", c_double), # 飞行高度 m + ("Ma0", c_double), # 马赫数 + ("dT0", c_double), # 温差 K + ("StepTime", c_double), # 步长 + ("Wf", c_double), # 主燃油流量 kg/h + ("Wf_After", c_double), # 加力燃油 kg/h + ("Angle_FanVane", c_double), # 风扇导叶 + ("Angle_CompVane", c_double), # 高压导叶 + ("A8", c_double), # 喉道面积 m^2 + ("AddPower", c_double) # 起动力矩 W + ] + +class EngOutPut(Structure): + _fields_ = [ + ("NL", c_double), # 风扇转速 (0,1) + ("NH", c_double), # 高压转速 (0,1) + ("T1t", c_double), # 进口温度 K + ("P1t", c_double), # 进口压力 kPa + ("P1s", c_double), # 进口静压 kPa + ("P3s", c_double), # 压气机出口压力 kPa + ("T5t", c_double), # 涡轮出口温度 K + ("P5t", c_double), # 涡轮出口压力 kPa + ("Wf_Main", c_double), # 实际主燃油流量 kg/h + ("T5t_Gas", c_double) # 涡轮出口实际气体温度 K + ] + +# ========================================== +# 2. DLL 封装类 +# ========================================== + +class AeroEngineDLL: + def __init__(self, dll_path=None): + if dll_path is None: + dll_path = os.path.join(os.path.dirname(__file__), "libEngine.dll") + + self.lib = WinDLL(dll_path) # 使用 WinDLL 因为是 stdcall + + # 配置函数原型 + self.lib.CreateEng.argtypes = [POINTER(POINTER(c_int)), Engine_Identity] + self.lib.CreateEng.restype = EngOutPut + + self.lib.EngStepGo.argtypes = [POINTER(c_int), EngInPut] + self.lib.EngStepGo.restype = EngOutPut + + self.lib.DestroyEng.argtypes = [POINTER(c_int)] + self.lib.DestroyEng.restype = EngOutPut + + self.h_engine = None + self.dt = 0.02 # 默认仿真步长 20ms + + def reset(self): + # 如果已经存在实例，先销毁 + if self.h_engine: + self.lib.DestroyEng(self.h_engine) + + # 创建新实例 + self.h_engine = POINTER(c_int)() + identity = Engine_Identity(1, 1.0, 1.0) + self.lib.CreateEng(byref(self.h_engine), identity) + + # 初始化一个默认的输入状态 + self.current_input = EngInPut() + self.current_input.Altp = 0.0 + self.current_input.Ma0 = 0.0 + self.current_input.dT0 = 0.0 + self.current_input.StepTime = self.dt + self.current_input.Wf = 100.0 # 初始燃油 + self.current_input.A8 = 0.1 # 初始 A8 + self.current_input.AddPower = 0.0 + + # 跑一步获取初始状态 + out = self.lib.EngStepGo(self.h_engine, self.current_input) + return out + + def step(self, action_dict): + """ + action_dict: 包含具体物理值的字典 + {'Wf': 200.0, 'A8': 0.12, ...} + """ + # 更新输入结构体 + if 'Wf' in action_dict: self.current_input.Wf = action_dict['Wf'] + if 'A8' in action_dict: self.current_input.A8 = action_dict['A8'] + if 'FanVane' in action_dict: self.current_input.Angle_FanVane = action_dict['FanVane'] + if 'CompVane' in action_dict: self.current_input.Angle_CompVane = action_dict['CompVane'] + if 'AddPower' in action_dict: self.current_input.AddPower = action_dict['AddPower'] + + # 确保步长正确 + self.current_input.StepTime = self.dt + + # 调用 DLL + out = self.lib.EngStepGo(self.h_engine, self.current_input) + return out + + def close(self): + if self.h_engine: + try: + self.lib.DestroyEng(self.h_engine) + except OSError: + # 忽略销毁时的访问冲突，可能是 DLL 内部状态问题 + pass + self.h_engine = None + + # 尝试释放 DLL + # 注意：这通常不是必须的，但在某些情况下（如反复加载/卸载或 DLL 内部有全局状态）可能有帮助 + if hasattr(self, 'lib'): + try: + import _ctypes + _ctypes.FreeLibrary(self.lib._handle) + except: + pass + del self.lib diff --git a/engine_env/demo.py b/engine_env/demo.py new file mode 100644 index 0000000..d8f4378 --- /dev/null +++ b/engine_env/demo.py @@ -0,0 +1,189 @@ +import matplotlib.pyplot as plt +import matplotlib.ticker as ticker +import numpy as np +from core_model import AeroEngineDLL +import os + +# ========================================== +# 范围测试示例 +# ========================================== + +def run_simulation(): + print("=== 开始发动机控制量范围测试 (优化版) ===") + + # 1. 初始化模型 + dll_path = os.path.join(os.path.dirname(__file__), "libEngine.dll") + if not os.path.exists(dll_path): + dll_path = os.path.join(os.path.dirname(__file__), "..", "libEngine.dll") + + sim = AeroEngineDLL(dll_path) + out = sim.reset() + + # 初始状态 + current_wf = 100.0 + current_add_power = 90000.0 + current_a8 = 0.1 + current_fan_vane = 0.0 + current_comp_vane = 0.0 + + # 数据记录 + data = { + 'time': [], 'nh': [], 'nl': [], 't5': [], 'p3': [], 'wf': [], + 'a8': [], 'fan_vane': [], 'comp_vane': [] + } + + # 辅助记录函数 + def record(t, o, wf, a8, fv, cv): + data['time'].append(t) + data['nh'].append(o.NH) + data['nl'].append(o.NL) + data['t5'].append(o.T5t) + data['p3'].append(o.P3s) + data['wf'].append(wf) + data['a8'].append(a8) + data['fan_vane'].append(fv) + data['comp_vane'].append(cv) + + record(0.0, out, current_wf, current_a8, current_fan_vane, current_comp_vane) + + total_time = 0.0 + dt = 0.02 + + # ========================================== + # Phase 1: 起动与稳态 (0 - 20s) + # ========================================== + print("Phase 1: 起动与稳态 (0-20s)...") + while total_time < 20.0: + if out.NH > 0.25 and current_add_power > 0: + current_add_power = 0.0 + current_wf = 2000.0 # 提高稳态燃油至 2000 (原 1000) 以测试高工况下的导叶效果 + + action = { + 'Wf': current_wf, 'AddPower': current_add_power, 'A8': current_a8, + 'FanVane': current_fan_vane, 'CompVane': current_comp_vane + } + out = sim.step(action) + total_time += dt + record(total_time, out, current_wf, current_a8, current_fan_vane, current_comp_vane) + + # ========================================== + # Phase 2: A8 扫描 (Step Response) + # A8: 0.01 -> 1.0 + # ========================================== + print("Phase 2: 测试 A8 阶跃变化 (0.01 -> 1.0)...") + steps_per_point = 50 # 1秒 per point + a8_targets = np.linspace(0.01, 1.0, 40) + + for val in a8_targets: + current_a8 = val + for _ in range(steps_per_point): + action = {'Wf': current_wf, 'AddPower': 0.0, 'A8': current_a8, 'FanVane': current_fan_vane, 'CompVane': current_comp_vane} + out = sim.step(action) + total_time += dt + record(total_time, out, current_wf, current_a8, current_fan_vane, current_comp_vane) + + # 恢复 A8 并稳定 + current_a8 = 0.1 + print("Restabilizing A8...") + for _ in range(100): # 2秒 + action = {'Wf': current_wf, 'AddPower': 0.0, 'A8': current_a8, 'FanVane': current_fan_vane, 'CompVane': current_comp_vane} + out = sim.step(action) + total_time += dt + record(total_time, out, current_wf, current_a8, current_fan_vane, current_comp_vane) + + # ========================================== + # Phase 3: Fan Vane 扫描 (Step Response) + # Range: -20 -> +20 度 + # ========================================== + print("Phase 3: 测试 Fan Vane 阶跃变化 (-20 -> 20)...") + vane_targets = np.linspace(-20.0, 20.0, 40) + + for val in vane_targets: + current_fan_vane = val + for _ in range(steps_per_point): + action = {'Wf': current_wf, 'AddPower': 0.0, 'A8': current_a8, 'FanVane': current_fan_vane, 'CompVane': current_comp_vane} + out = sim.step(action) + total_time += dt + record(total_time, out, current_wf, current_a8, current_fan_vane, current_comp_vane) + + # 恢复 Fan Vane 并稳定 + current_fan_vane = 0.0 + print("Restabilizing Fan Vane...") + for _ in range(100): + action = {'Wf': current_wf, 'AddPower': 0.0, 'A8': current_a8, 'FanVane': current_fan_vane, 'CompVane': current_comp_vane} + out = sim.step(action) + total_time += dt + record(total_time, out, current_wf, current_a8, current_fan_vane, current_comp_vane) + + # ========================================== + # Phase 4: Comp Vane 扫描 (Step Response) + # Range: -20 -> +20 度 + # ========================================== + print("Phase 4: 测试 Comp Vane 阶跃变化 (-20 -> 20)...") + + for val in vane_targets: + current_comp_vane = val + for _ in range(steps_per_point): + action = {'Wf': current_wf, 'AddPower': 0.0, 'A8': current_a8, 'FanVane': current_fan_vane, 'CompVane': current_comp_vane} + out = sim.step(action) + total_time += dt + record(total_time, out, current_wf, current_a8, current_fan_vane, current_comp_vane) + + sim.close() + print("测试结束，正在绘图...") + plot_range_test(data) + +def plot_range_test(data): + t = data['time'] + + plt.rcParams['font.family'] = 'sans-serif' + fig, axes = plt.subplots(5, 1, figsize=(12, 16), sharex=True) + + # 1. Speed (NH & NL) + axes[0].plot(t, data['nh'], color='blue', label='NH') + axes[0].plot(t, data['nl'], color='cyan', linestyle='--', label='NL') + axes[0].set_ylabel('Speed (%)') + axes[0].grid(True, linestyle='--', alpha=0.7) + axes[0].legend(loc='upper right') + axes[0].set_title('Engine Response to Control Inputs') + + # 2. T5 + axes[1].plot(t, data['t5'], color='red', label='T5') + axes[1].set_ylabel('T5 (K)') + axes[1].grid(True, linestyle='--', alpha=0.7) + axes[1].legend(loc='upper right') + + # 3. P3 + axes[2].plot(t, data['p3'], color='orange', label='P3') + axes[2].set_ylabel('P3 (kPa)') + axes[2].grid(True, linestyle='--', alpha=0.7) + axes[2].legend(loc='upper right') + + # 4. A8 Input + axes[3].plot(t, data['a8'], color='green', label='A8 Input') + axes[3].set_ylabel('A8 (m^2)') + axes[3].grid(True, linestyle='--', alpha=0.7) + axes[3].legend(loc='upper right') + axes[3].yaxis.set_major_locator(ticker.MaxNLocator(nbins=5)) + axes[3].minorticks_on() + axes[3].grid(which='minor', linestyle=':', alpha=0.4) + + # 5. Vane Inputs + axes[4].plot(t, data['fan_vane'], color='purple', label='Fan Vane') + axes[4].plot(t, data['comp_vane'], color='brown', label='Comp Vane') + axes[4].set_ylabel('Vane Angle (deg)') + axes[4].set_xlabel('Time (s)') + axes[4].grid(True, linestyle='--', alpha=0.7) + axes[4].legend(loc='upper right') + axes[4].yaxis.set_major_locator(ticker.MaxNLocator(nbins=5)) + axes[4].minorticks_on() + axes[4].grid(which='minor', linestyle=':', alpha=0.4) + + plt.tight_layout() + save_path = os.path.join(os.path.dirname(__file__), 'range_test_result.png') + plt.savefig(save_path) + print(f"结果图已保存至: {save_path}") + plt.show() + +if __name__ == "__main__": + run_simulation() diff --git a/engine_env/engine_env.py b/engine_env/engine_env.py new file mode 100644 index 0000000..3fb7d45 --- /dev/null +++ b/engine_env/engine_env.py @@ -0,0 +1,277 @@ +import gymnasium as gym +from gymnasium import spaces +import numpy as np + +from engine_env.core_model import AeroEngineDLL +from engine_env.schedule import ControlSchedule + +class AeroEngineGymEnv(gym.Env): + """ + 航空发动机控制环境 (Aero-Engine Control Environment) + + 特点: + 1. 外部通过 set_pla() 控制任务输入。 + 2. 包含起动阶段和稳态控制阶段的复杂约束判定。 + 3. 奖励函数包含调节时间、超调量、起动功、平滑性等工程指标。 + + Observation Space (8维): + [NH, NL, T5, P3, Target_NH, Target_NL, Error_NH, Error_NL] + + Action Space (5维, 范围 [-1, 1]): + [FanVane, CompVane, AddPower, Wf, A8] + """ + metadata = {"render_modes": ["human"]} + + def __init__(self): + # 1. 加载核心组件 + self.engine = AeroEngineDLL() + self.schedule = ControlSchedule() + + # 2. 仿真参数 + self.current_step = 0 + self.dt = 0.02 # 仿真步长 20ms + self.max_steps = 1000 # 最大步数 (20秒) + + # 3. 核心状态变量 + self.pla = 0.0 # 油门杆角度 (由外部控制) + + # 4. 定义动作空间 (5个控制量) + # 顺序: [风扇导叶, 高压导叶, 起动功率, 燃油Wf, A8] + self.action_space = spaces.Box(low=-1.0, high=1.0, shape=(5,), dtype=np.float32) + + # 5. 定义观测空间 (8维) + self.observation_space = spaces.Box(low=-np.inf, high=np.inf, shape=(8,), dtype=np.float32) + + # 6. 物理范围映射 (Physical Limits) + self.limits = { + 'Wf': (50.0, 5000.0), # kg/h (燃油流量) + 'A8': (0.2, 0.6), # m^2 (喷管面积) + 'FanVane': (-20.0, 20.0), # deg (风扇导叶) + 'CompVane': (-20.0, 20.0), # deg (压气机导叶) + 'AddPower': (0.0, 100000.0) # W (起动电机功率) + } + + # 内部缓存变量 + self.t_nh, self.t_nl = 0.0, 0.0 + self.limit_t5, self.limit_p3 = 1200.0, 2000.0 + + # 计时器 (用于约束判定) + self.timer_startup = 0.0 + self.timer_settle_5 = 0.0 + self.timer_settle_1 = 0.0 + + # 平滑性权重 (Smoothness Weights) + # 燃油(Wf)和喷管(A8)给高权重，AddPower允许突变 + self.smooth_weights = np.array([1.0, 1.0, 0.5, 5.0, 2.0], dtype=np.float32) + + # 初始化上一帧动作缓存 + self.last_action = np.zeros(5, dtype=np.float32) + + def reset(self, seed=None, options=None): + """ + 重置环境到初始状态 + options: {'pla': float} 可指定初始 PLA + """ + super().reset(seed=seed) + self.current_step = 0 + + # 重置所有计时器 + self.timer_startup = 0.0 + self.timer_settle_5 = 0.0 + self.timer_settle_1 = 0.0 + + # 重置动作缓存 + self.last_action = np.zeros(5, dtype=np.float32) + + # 设定初始 PLA + self.pla = options.get('pla', 0.0) if options else 0.0 + + # 复位底层模型 + out = self.engine.reset() + + # 更新初始目标 + self.t_nh, self.t_nl, self.limit_t5, self.limit_p3 = self.schedule.get_targets(self.pla) + + return self._get_obs(out), {"PLA": self.pla} + + def set_pla(self, pla_value): + """ + 【外部接口】设置当前的 PLA (油门杆角度) + """ + self.pla = np.clip(pla_value, 0.0, 110.0) + + def step(self, action): + """ + 环境交互核心函数 + action: AI 输出的归一化动作 [-1, 1] + """ + self.current_step += 1 + + # ======================================================= + # 1. 查表：更新控制目标和限制 (基于当前的 self.pla) + # ======================================================= + self.t_nh, self.t_nl, self.limit_t5, self.limit_p3 = self.schedule.get_targets(self.pla) + + # ======================================================= + # 2. 动作反归一化 (AI [-1, 1] -> Physical Value) + # ======================================================= + phys_action = { + 'FanVane': self._denormalize(action[0], 'FanVane'), + 'CompVane': self._denormalize(action[1], 'CompVane'), + 'AddPower': self._denormalize(action[2], 'AddPower'), + 'Wf': self._denormalize(action[3], 'Wf'), + 'A8': self._denormalize(action[4], 'A8'), + } + + # ======================================================= + # 3. 执行物理仿真 + # ======================================================= + out = self.engine.step(phys_action) + + # ======================================================= + # 4. 计算奖励 (Reward Function) + # ======================================================= + # 传入 normalized action (action) 用于计算平滑度 + total_reward = self._calculate_official_reward(out, phys_action, self.pla, action) + + # ======================================================= + # 5. 终止条件判定 (Terminated Check) + # ======================================================= + terminated, term_penalty = self._check_terminated(out, self.pla) + total_reward += term_penalty + + truncated = (self.current_step >= self.max_steps) + + # ======================================================= + # 6. 组装返回信息 + # ======================================================= + obs = self._get_obs(out) + + info = { + "PLA": self.pla, + "Real_NH": out.NH, + "Target_NH": self.t_nh, + "Real_T5": out.T5t, + "Action_Physical": phys_action + } + + return obs, total_reward, terminated, truncated, info + + def _check_terminated(self, out, current_pla): + """ + 终止判定 + """ + # ======================================================= + # 全局硬约束 - 只有炸机才重开 + # ======================================================= + + # 1. 严重超温 + if out.T5t > 1400.0: + # print(f"[Terminated] T5 High: {out.T5t:.1f}") + return True, -1000.0 + + # 2. 严重超转 + + if out.NH > 1.10: + # print(f"[Terminated] NH High: {out.NH:.3f}") + return True, -1000.0 + + if out.NL > 1.10: + # print(f"[Terminated] NL High: {out.NL:.3f}") + return True, -1000.0 + + + return False, 0.0 + + def _calculate_official_reward(self, out, phys_action, current_pla, normalized_action): + """ + 官方奖励计算函数 (Official Reward Function) + """ + step_reward = 0.0 + + # 提取目标和当前值 + t_nh, t_nl = self.t_nh, self.t_nl + nh, nl = out.NH, out.NL + + # 计算误差 + err_nh = t_nh - nh + err_nl = t_nl - nl + abs_err_nh = abs(err_nh) + + # ======================================================= + # 阶段 A: 起动阶段 - PLA < 15 + # ======================================================= + if current_pla < 15.0: + # 1. 【起动时间】 + if nh < 0.60: + step_reward -= 0.5 # 时间惩罚基数 + + # 2. 【起动功】 + # 逻辑：使用的电功率越大，扣分越多 + power_penalty = (phys_action['AddPower'] / 100000.0) * 0.2 + step_reward -= power_penalty + + # 3. 跟踪引导 (主要看 NH) + step_reward -= (err_nh ** 2) * 20.0 + + # ======================================================= + # 阶段 B: 正常工作阶段- PLA >= 15 + # ======================================================= + else: + # 1. 【调节时间】 + if abs_err_nh > 0.02: # 2% 误差带 + step_reward -= 0.5 # 调节时间惩罚 + + # 2. 【超调量】 + if nh > t_nh: + overshoot_amount = nh - t_nh + step_reward -= overshoot_amount * 50.0 # 重罚超调 + + # 3. 基础跟踪 (主要看 NL 推力, 兼顾 NH) + step_reward -= (err_nl ** 2) * 20.0 + step_reward -= (err_nh ** 2) * 10.0 + + # ======================================================= + # 全局安全约束 + # ======================================================= + + # 1. 【涡轮后温度】 + if out.T5t > self.limit_t5: + over_temp = out.T5t - self.limit_t5 + step_reward -= over_temp * 1 # 每超 1K 扣 1分 + + # 2. 【压气机后压力】 + if out.P3s > self.limit_p3: + over_pres = out.P3s - self.limit_p3 + step_reward -= over_pres * 0.5 # 每超 1kPa 扣 0.5分 + + # ======================================================= + # 3. 动作平滑性 (Action Smoothness) + # ======================================================= + # 计算动作变化率: delta = current - last + delta_action = normalized_action - self.last_action + + # 计算加权平方和 (L2 Norm) + smoothness_cost = np.sum(self.smooth_weights * np.square(delta_action)) + + # 乘以系数 (调节平滑性在总分中的占比) + step_reward -= smoothness_cost * 0.1 + + # 更新上一帧动作缓存 + self.last_action = normalized_action.copy() + + return step_reward + + def _get_obs(self, out): + state = np.array([out.NH, out.NL, out.T5t, out.P3s], dtype=np.float32) + targets = np.array([self.t_nh, self.t_nl], dtype=np.float32) + errors = targets - state[:2] + return np.concatenate([state, targets, errors]) + + def _denormalize(self, val, key): + min_v, max_v = self.limits[key] + return (val + 1.0) / 2.0 * (max_v - min_v) + min_v + + def close(self): + if self.engine: + self.engine.close() diff --git a/engine_env/libEngine.dll b/engine_env/libEngine.dll new file mode 100644 index 0000000..1da2190 Binary files /dev/null and b/engine_env/libEngine.dll differ diff --git a/engine_env/range_test_result.png b/engine_env/range_test_result.png new file mode 100644 index 0000000..0915761 Binary files /dev/null and b/engine_env/range_test_result.png differ diff --git a/engine_env/schedule.py b/engine_env/schedule.py new file mode 100644 index 0000000..59359b0 --- /dev/null +++ b/engine_env/schedule.py @@ -0,0 +1,95 @@ +import numpy as np +from scipy.interpolate import interp1d + +class ControlSchedule: + """ + 航空发动机控制计划 (Control Schedule) + 输入: PLA (角度 0~110) + 输出: Target NH, Target NL (0.0~1.0), Limit T5 + """ + def __init__(self): + + # ========================================== + # 1. 定义控制计划数据点 + # ========================================== + + # --- NH (高压转速) 计划 --- + # 依据图片1：在 PLA=80 时达到 1.0 + # 0 -> 14.99: 0 + # 15.0 : 0.7 + # 80.0 : 1.0 + self._pla_nh = np.array([0.0, 14.99, 15.0, 80.0, 110.0]) + self._val_nh = np.array([0.0, 0.0, 0.70, 1.00, 1.00]) + + # --- NL (低压转速) 计划 --- + # 依据图片2：在 PLA=90 时达到 1.0 + # 0 -> 14.99: 0 + # 15.0 : 0.5 + # 90.0 : 1.0 + self._pla_nl = np.array([0.0, 14.99, 15.0, 90.0, 110.0]) + self._val_nl = np.array([0.0, 0.0, 0.50, 1.00, 1.00]) + + # --- 限制值 (T5, P3) 计划 --- + # 保持原有逻辑：在 PLA=100 时达到最大限制 + self._pla_lim = np.array([0.0, 14.99, 15.0, 100.0, 110.0]) + + # Limit T5 (温度限制) [K] + # 15度时1200K, 100度时1500K + self._val_t5 = np.array([1200.0, 1200.0, 1200.0, 1500.0, 1500.0]) + + # Limit P3 (压力限制) [kPa] + # 15度时2000kPa, 100度时2500kPa + self._val_p3 = np.array([800.0, 800.0, 2000.0, 2500.0, 2500.0]) + + + # ========================================== + # 2. 构建插值函数 + # ========================================== + self.f_nh = interp1d(self._pla_nh, self._val_nh, kind='linear', fill_value="extrapolate") + self.f_nl = interp1d(self._pla_nl, self._val_nl, kind='linear', fill_value="extrapolate") + self.f_t5 = interp1d(self._pla_lim, self._val_t5, kind='linear', fill_value="extrapolate") + self.f_p3 = interp1d(self._pla_lim, self._val_p3, kind='linear', fill_value="extrapolate") + + def get_targets(self, pla): + """ + 输入: pla (角度, 0.0 ~ 110.0) + """ + # 稍微做个限幅，防止超出定义域太多 + pla = np.clip(pla, 0.0, 110.0) + + t_nh = float(self.f_nh(pla)) + t_nl = float(self.f_nl(pla)) + l_t5 = float(self.f_t5(pla)) + l_p3 = float(self.f_p3(pla)) + + return t_nh, t_nl, l_t5, l_p3 + +# ========================================== +# 自测代码 +# ========================================== +if __name__ == "__main__": + import matplotlib.pyplot as plt + + sch = ControlSchedule() + + # 测试关键点 + test_plas = [0, 10, 14.9, 15.0, 15.1, 57.5, 100, 105] + print(f"{'PLA':<10} {'NH':<10} {'NL':<10}") + print("-" * 30) + for p in test_plas: + nh, nl, _, _ = sch.get_targets(p) + print(f"{p:<10} {nh:<10.4f} {nl:<10.4f}") + + # 画图确认阶跃形状 + x = np.linspace(0, 110, 500) + y_nh = [sch.get_targets(i)[0] for i in x] + + plt.figure() + plt.plot(x, y_nh, label='Target NH') + plt.axvline(15, color='r', linestyle='--', alpha=0.5, label='Idle Point (15 deg)') + plt.title("PLA to Engine Speed Schedule") + plt.xlabel("PLA (Degree)") + plt.ylabel("Normalized Speed") + plt.legend() + plt.grid(True) + plt.show() diff --git a/test_pid.py b/test_pid.py new file mode 100644 index 0000000..9fa3a10 --- /dev/null +++ b/test_pid.py @@ -0,0 +1,188 @@ +import numpy as np +import matplotlib.pyplot as plt +from engine_env.engine_env import AeroEngineGymEnv + +# ========================================== +# 1. 增量式 PID 类 (保持不变) +# ========================================== +class IncrementalPIDController: + def __init__(self, kp, ki, kd, output_min=-1.0, output_max=1.0): + self.kp = kp + self.ki = ki + self.kd = kd + + self.min_val = output_min + self.max_val = output_max + self.error_prev = 0.0 + self.error_prev2 = 0.0 + self.current_output = 0.0 + + def set_current_output(self, value): + self.current_output = np.clip(value, self.min_val, self.max_val) + self.error_prev = 0.0 + self.error_prev2 = 0.0 + + def update(self, error, dt): + delta_p = self.kp * (error - self.error_prev) + delta_i = self.ki * error * dt + delta_d = self.kd * (error - 2*self.error_prev + self.error_prev2) / dt + + delta_u = delta_p + delta_i + delta_d + self.current_output += delta_u + self.current_output = np.clip(self.current_output, self.min_val, self.max_val) + + self.error_prev2 = self.error_prev + self.error_prev = error + return self.current_output + +# ========================================== +# 2. 主测试逻辑 +# ========================================== +def run_pid_test(): + env = AeroEngineGymEnv() + + # ----------------------------------------------------- + # 【关键修改】定义正确的 PID 配对 + # ----------------------------------------------------- + + # Loop 1: 目标 NH -> 控制 Wf (燃油) + # 逻辑: NH 低 -> 加油 (正反馈) + pid_nh_wf = IncrementalPIDController(kp=10.0, ki=6.0, kd=0.1) + + # Loop 2: 目标 NL -> 控制 A8 (喷管) + pid_nl_a8 = IncrementalPIDController(kp=-4.0, ki=3.0, kd=0.05) + + total_steps = 3000 # 60秒 + obs, info = env.reset(options={'pla': 0.0}) + + # 状态标志位 + pid_initialized = False + + history = { + 'time': [], 'PLA': [], 'NH': [], 'Target_NH': [], 'NL': [], 'Target_NL': [], + 'Wf_Action': [], 'A8_Action': [], 'Reward': [] + } + + print("开始修正后的 PID 仿真 (NH->Wf, NL->A8)...") + + for step in range(total_steps): + t = step * env.dt + + # --- A. 任务剖面 --- + if t < 2.0: current_pla = 0.0 + elif t < 20.0: current_pla = 15.0 + else: current_pla = 60.0 + + env.set_pla(current_pla) + + # --- B. 获取状态 --- + real_nh, real_nl = obs[0], obs[1] + target_nh, target_nl = obs[4], obs[5] + + # 计算误差 + error_nh = target_nh - real_nh + error_nl = target_nl - real_nl + + action = np.zeros(5) + + # ========================================================= + # C. 分阶段控制 + # ========================================================= + + # --- 阶段 1: 起动 (Open Loop) --- + if current_pla >= 15.0 and real_nh < 0.60: + action[2] = 0.8 # 起动机 + action[3] = -0.6 # 点火油量 + action[4] = 1.0 # 喷管全开 + pid_initialized = False + + # --- 阶段 2: 闭环控制 (Closed Loop) --- + elif current_pla >= 15.0 and real_nh >= 0.60: + + # 无扰切换 + if not pid_initialized: + print(f"[Switch] t={t:.2f}s, 切入闭环。PID1: NH->Fuel(-0.6), PID2: NL->A8(1.0)") + # NH控制燃油，继承 -0.6 + pid_nh_wf.set_current_output(-0.6) + # NL控制喷管，继承 1.0 + pid_nl_a8.set_current_output(1.0) + pid_initialized = True + + action[2] = -1.0 + + # 【核心修正】 PID 计算 + # Loop 1: 用 NH 的误差算 Wf + wf_out = pid_nh_wf.update(error_nh, env.dt) + + # Loop 2: 用 NL 的误差算 A8 + a8_out = pid_nl_a8.update(error_nl, env.dt) + + action[3] = wf_out + action[4] = a8_out + + # --- 阶段 3: 停机 --- + else: + action = np.array([0,0,-1,-1,1]) + pid_initialized = False + + # --- D. 执行 --- + obs, reward, terminated, truncated, info = env.step(action) + + # 记录 + history['time'].append(t) + history['PLA'].append(current_pla) + history['NH'].append(real_nh) + history['Target_NH'].append(target_nh) + history['NL'].append(real_nl) + history['Target_NL'].append(target_nl) + history['Wf_Action'].append(action[3]) + history['A8_Action'].append(action[4]) + history['Reward'].append(reward) + + if terminated: + print(f"Terminated at {t:.2f}s") + break + + env.close() + plot_results(history) + +def plot_results(history): + t = history['time'] + plt.figure(figsize=(12, 8)) + + # 1. Speed + plt.subplot(2, 2, 1) + plt.plot(t, history['Target_NH'], 'r--', alpha=0.6) + plt.plot(t, history['NH'], 'r', label='NH (Controlled by Fuel)') + plt.plot(t, history['Target_NL'], 'b--', alpha=0.6) + plt.plot(t, history['NL'], 'b', label='NL (Controlled by A8)') + plt.legend() + plt.title('Rotor Speeds') + plt.grid(True) + + # 2. Action + plt.subplot(2, 2, 2) + plt.plot(t, history['Wf_Action'], 'orange', label='Fuel (Wf)') + plt.plot(t, history['A8_Action'], 'green', label='Nozzle (A8)') + plt.legend() + plt.title('Control Actions') + plt.ylim(-1.1, 1.1) + plt.grid(True) + + # 3. PLA + plt.subplot(2, 2, 3) + plt.plot(t, history['PLA'], 'k') + plt.title('PLA') + plt.grid(True) + + # 4. Reward + plt.subplot(2, 2, 4) + plt.plot(t, history['Reward']) + plt.title('Reward') + plt.grid(True) + + plt.tight_layout() + plt.show() + +if __name__ == "__main__": + run_pid_test()