diff --git a/brayton_cycle/__init__.py b/brayton_cycle/__init__.py index 2e26aab..fda55e3 100644 --- a/brayton_cycle/__init__.py +++ b/brayton_cycle/__init__.py @@ -9,7 +9,16 @@ from .components import ( Turbine, ) from .cycles import BraytonCycle +from .optimization import ( + optimize_rc_fixed_param, + optimize_rc_param, + plot_optimization_landscape, + plot_sweep_optimization_results, + scan_rc_efficiency, + sweep_and_optimize_rc, +) from .properties import CO2PropertyCalculator +from .sensitivity import evaluate_rc_efficiency __all__ = [ "BraytonCycle", @@ -20,4 +29,11 @@ __all__ = [ "Heater", "Recuperator", "Turbine", + "evaluate_rc_efficiency", + "optimize_rc_fixed_param", + "optimize_rc_param", + "plot_optimization_landscape", + "plot_sweep_optimization_results", + "scan_rc_efficiency", + "sweep_and_optimize_rc", ] diff --git a/brayton_cycle/cycles.py b/brayton_cycle/cycles.py index 42c6b51..43e5f9c 100644 --- a/brayton_cycle/cycles.py +++ b/brayton_cycle/cycles.py @@ -1,9 +1,5 @@ # -*- coding: utf-8 -*- -"""Cycle definitions, efficiency calculations, optimization, and plotting.""" - -from scipy.optimize import minimize_scalar -import numpy as np -import matplotlib.pyplot as plt +"""Cycle definitions and efficiency calculations.""" from .components import ( Compressor, @@ -16,11 +12,11 @@ from .components import ( from .properties import CO2PropertyCalculator -class BraytonCycle(): - """循环计算""" - def __init__(self, name, refprop_path = "C:/Program Files (x86)/REFPROP 10.0+/REFPROP"): +class BraytonCycle: + """CO2 Brayton cycle calculator.""" + + def __init__(self, name, refprop_path="C:/Program Files (x86)/REFPROP 10.0+/REFPROP"): self.name = name - self.property_calculator = None self.compressor = None self.turbine = None self.recuperator = None @@ -31,445 +27,218 @@ class BraytonCycle(): self.property_calculator = CO2PropertyCalculator(self.refprop_path) def cycle_eff_calculator(self): - if type(self.compressor) == list: - Wc = (self.compressor[0].variables['Wc'] + - self.compressor[1].variables['Wc']) + if isinstance(self.compressor, list): + Wc = sum(compressor.variables["Wc"] for compressor in self.compressor) else: - Wc = self.compressor.variables['Wc'] - Wt = self.turbine.variables['Wt'] - Q_input = self.heater.variables['Q_in'] - cycle_eff = (Wt - Wc) / Q_input - return cycle_eff - - def SC(self, T_low, T_high, p_low, p_high, param = None): - if param == None: - param = { - 'compressor_eff': 0.98, - 'turbine_eff': 0.95 - } - # 定义循环组件 - self.heater = Heater(name = "Main heater") - self.condenser = Condenser(name = "Main condenser") - self.compressor = Compressor(name = "Main compressor", eff = param['compressor_eff']) - self.turbine = Turbine(name = "Turbine", eff = param['turbine_eff']) - - # 计算循环参数 + Wc = self.compressor.variables["Wc"] + + Wt = self.turbine.variables["Wt"] + Q_input = self.heater.variables["Q_in"] + return (Wt - Wc) / Q_input + + def SC(self, T_low, T_high, p_low, p_high, param=None): + """Simple Brayton cycle.""" + defaults = { + "compressor_eff": 0.98, + "turbine_eff": 0.95, + } + cycle_params = defaults if param is None else {**defaults, **param} + + self.heater = Heater(name="Main heater") + self.condenser = Condenser(name="Main condenser") + self.compressor = Compressor( + name="Main compressor", + eff=cycle_params["compressor_eff"], + ) + self.turbine = Turbine(name="Turbine", eff=cycle_params["turbine_eff"]) + self.compressor.calculator(p_low, T_low, p_high, self.property_calculator) self.turbine.calculator(p_high, T_high, p_low, self.property_calculator) - self.heater.calculator(self.compressor.variables['outlet_state'], self.turbine.variables['inlet_state']) - self.condenser.calculator(self.turbine.variables['outlet_state'], self.compressor.variables['inlet_state']) - - # 计算循环效率 - cycle_eff = self.cycle_eff_calculator() - return cycle_eff - + self.heater.calculator( + self.compressor.variables["outlet_state"], + self.turbine.variables["inlet_state"], + ) + self.condenser.calculator( + self.turbine.variables["outlet_state"], + self.compressor.variables["inlet_state"], + ) + + return self.cycle_eff_calculator() + def SRC(self, T_low, T_high, p_low, p_high, param=None, ploss=0.0): - if param is None or 'recuperator_eff' not in param: - param = { - 'compressor_eff': 0.98, - 'turbine_eff': 0.95, - 'recuperator_eff': 0.85 - } - - # 定义循环组件 - self.heater = Heater(name = "Main heater") - self.condenser = Condenser(name = "Main condenser") - self.compressor = Compressor(name = "Main compressor", eff = param['compressor_eff']) - self.turbine = Turbine(name = "Turbine", eff = param['turbine_eff']) - self.recuperator = Recuperator(name = "Recuperator", eff = param['recuperator_eff']) - - # 计算循环参数 + """Simple recuperated Brayton cycle.""" + defaults = { + "compressor_eff": 0.98, + "turbine_eff": 0.95, + "recuperator_eff": 0.85, + } + cycle_params = defaults if param is None else {**defaults, **param} + + self.heater = Heater(name="Main heater") + self.condenser = Condenser(name="Main condenser") + self.compressor = Compressor( + name="Main compressor", + eff=cycle_params["compressor_eff"], + ) + self.turbine = Turbine(name="Turbine", eff=cycle_params["turbine_eff"]) + self.recuperator = Recuperator( + name="Recuperator", + eff=cycle_params["recuperator_eff"], + ) + self.compressor.calculator(p_low, T_low, p_high, self.property_calculator) self.turbine.calculator(p_high, T_high, p_low, self.property_calculator) - self.recuperator.calculator(self.compressor.variables['outlet_state'], - self.turbine.variables['outlet_state'], 0, - ploss, self.property_calculator) - self.condenser.calculator(self.recuperator.variables['hot_outlet_state'], - self.compressor.variables['inlet_state']) - self.heater.calculator(self.recuperator.variables['cold_outlet_state'], - self.turbine.variables['inlet_state']) - - # 计算循环效率 - cycle_eff = self.cycle_eff_calculator() - return cycle_eff - - def RC(self, T_low, T_high, p_low, p_high, ploss, param = None): - if (param is None or - 'recompressor_eff' not in param or - 'highT_recuperator_eff' not in param or - 'x' not in param): - param = { - 'compressor_eff': 0.98, - 'turbine_eff': 0.95, - 'recuperator_eff': 0.85, - 'recompressor_eff': 0.98, - 'highT_recuperator_eff': 0.85, - 'x': 0.9 - } - self.x = param['x'] - # 定义循环组件 - self.heater = Heater(name = "Main heater") - self.condenser = Condenser(name = "Main_condenser") - self.turbine = Turbine(name = "Main Turbine", - eff = param['turbine_eff']) - main_compressor = Compressor(name = "Main compressor", - eff = param['compressor_eff']) - recompressor = Compressor(name = "Recompressor", - eff = param['recompressor_eff']) - lrecuperator = Recuperator(name = "Low Temperature recuprerator", - eff = param['recuperator_eff'], - x = self.x) - hrecuperator = Recuperator(name = "High Temperature recuperator", - eff = param['highT_recuperator_eff']) - self.concentrator = Concentrator(name = "Concentrator") + self.recuperator.calculator( + self.compressor.variables["outlet_state"], + self.turbine.variables["outlet_state"], + 0, + ploss, + self.property_calculator, + ) + self.condenser.calculator( + self.recuperator.variables["hot_outlet_state"], + self.compressor.variables["inlet_state"], + ) + self.heater.calculator( + self.recuperator.variables["cold_outlet_state"], + self.turbine.variables["inlet_state"], + ) + + return self.cycle_eff_calculator() + + def RC(self, T_low, T_high, p_low, p_high, ploss, param=None): + """Recompression Brayton cycle.""" + defaults = { + "compressor_eff": 0.98, + "turbine_eff": 0.95, + "recuperator_eff": 0.85, + "recompressor_eff": 0.98, + "highT_recuperator_eff": 0.85, + "x": 0.9, + } + cycle_params = defaults if param is None else {**defaults, **param} + self.x = cycle_params["x"] + + self.heater = Heater(name="Main heater") + self.condenser = Condenser(name="Main_condenser") + self.turbine = Turbine( + name="Main Turbine", + eff=cycle_params["turbine_eff"], + ) + main_compressor = Compressor( + name="Main compressor", + eff=cycle_params["compressor_eff"], + ) + recompressor = Compressor( + name="Recompressor", + eff=cycle_params["recompressor_eff"], + ) + lrecuperator = Recuperator( + name="Low Temperature recuprerator", + eff=cycle_params["recuperator_eff"], + x=self.x, + ) + hrecuperator = Recuperator( + name="High Temperature recuperator", + eff=cycle_params["highT_recuperator_eff"], + ) + self.concentrator = Concentrator(name="Concentrator") self.recuperator = [] self.compressor = [] - - - # 计算组件进出口参数 - # 先算压缩机和涡轮 + main_compressor.calculator(p_low, T_low, p_high, self.property_calculator) - self.turbine.calculator(p_high*(1-ploss)**3, - T_high, p_low*(1-ploss)**(-3), - self.property_calculator) - - # 假定低温回热器出口参数并进行迭代 + self.turbine.calculator( + p_high * (1 - ploss) ** 3, + T_high, + p_low * (1 - ploss) ** (-3), + self.property_calculator, + ) + mw = self.property_calculator.mw - T_hr_inlet = ((self.turbine.variables['outlet_state']['T'] - + main_compressor.variables['outlet_state']['T'])/2) - # T_lr_inlet = self.turbine.variables['outlet_state']['T'] * 1.01 + T_hr_inlet = ( + self.turbine.variables["outlet_state"]["T"] + + main_compressor.variables["outlet_state"]["T"] + ) / 2 + max_iter = 300 relax_fac = 0.4 for i in range(max_iter): - hr_inlet_state_hot_mol = self.property_calculator.calculate_properties(P=p_high*(1-ploss), T=T_hr_inlet) - # 单位转换成mol + hr_inlet_state_hot_mol = self.property_calculator.calculate_properties( + P=p_high * (1 - ploss), + T=T_hr_inlet, + ) hr_inlet_state_hot = { - 'P': hr_inlet_state_hot_mol['P'], - 'T': hr_inlet_state_hot_mol['T'], - 'h': hr_inlet_state_hot_mol['h']/mw, - 's': hr_inlet_state_hot_mol['s']/mw - } - - hrecuperator.calculator(hr_inlet_state_hot, - self.turbine.variables['outlet_state'], 0, - ploss, self.property_calculator) - lrecuperator.calculator(main_compressor.variables['outlet_state'], - hrecuperator.variables['hot_outlet_state'], 1, - ploss, self.property_calculator) - recompressor.calculator(p_in = p_low/(1-ploss), p_out = p_high*(1-ploss), - T_in = lrecuperator.variables['hot_outlet_state']['T'], - property_calculator = self.property_calculator) - self.concentrator.calculator(recompressor.variables['outlet_state'], - lrecuperator.variables['cold_outlet_state'], - self.x, self.property_calculator) - - Tc_outlet = self.concentrator.variables['outlet_state']['T'] + "P": hr_inlet_state_hot_mol["P"], + "T": hr_inlet_state_hot_mol["T"], + "h": hr_inlet_state_hot_mol["h"] / mw, + "s": hr_inlet_state_hot_mol["s"] / mw, + } + + hrecuperator.calculator( + hr_inlet_state_hot, + self.turbine.variables["outlet_state"], + 0, + ploss, + self.property_calculator, + ) + lrecuperator.calculator( + main_compressor.variables["outlet_state"], + hrecuperator.variables["hot_outlet_state"], + 1, + ploss, + self.property_calculator, + ) + recompressor.calculator( + p_in=p_low / (1 - ploss), + p_out=p_high * (1 - ploss), + T_in=lrecuperator.variables["hot_outlet_state"]["T"], + property_calculator=self.property_calculator, + ) + self.concentrator.calculator( + recompressor.variables["outlet_state"], + lrecuperator.variables["cold_outlet_state"], + self.x, + self.property_calculator, + ) + + Tc_outlet = self.concentrator.variables["outlet_state"]["T"] err = abs(Tc_outlet - T_hr_inlet) if err <= 1e-5: break if i == max_iter - 1: - raise ValueError(f"迭代次数超过范围，当前误差{err:.4f}") - + raise ValueError(f"Iteration limit exceeded; current error={err:.4f}") + T_hr_inlet = relax_fac * Tc_outlet + (1 - relax_fac) * T_hr_inlet - - hrecuperator.calculator(self.concentrator.variables['outlet_state'], - self.turbine.variables['outlet_state'], 0, - ploss, self.property_calculator) - lrecuperator.calculator(main_compressor.variables['outlet_state'], - hrecuperator.variables['hot_outlet_state'], 1, - ploss, self.property_calculator) - - main_compressor.variables['Wc'] *= (1-self.x) - recompressor.variables['Wc'] *= self.x + + hrecuperator.calculator( + self.concentrator.variables["outlet_state"], + self.turbine.variables["outlet_state"], + 0, + ploss, + self.property_calculator, + ) + lrecuperator.calculator( + main_compressor.variables["outlet_state"], + hrecuperator.variables["hot_outlet_state"], + 1, + ploss, + self.property_calculator, + ) + + main_compressor.variables["Wc"] *= 1 - self.x + recompressor.variables["Wc"] *= self.x self.recuperator.append(lrecuperator) self.recuperator.append(hrecuperator) self.compressor.append(main_compressor) self.compressor.append(recompressor) - self.condenser.calculator(lrecuperator.variables['hot_outlet_state'], - main_compressor.variables['inlet_state']) - self.condenser.variables['Q_out'] *= (1-self.x) - self.heater.calculator(hrecuperator.variables['cold_outlet_state'], - self.turbine.variables['inlet_state']) - # 计算循环效率 - cycle_eff = self.cycle_eff_calculator() - return cycle_eff + self.condenser.calculator( + lrecuperator.variables["hot_outlet_state"], + main_compressor.variables["inlet_state"], + ) + self.condenser.variables["Q_out"] *= 1 - self.x + self.heater.calculator( + hrecuperator.variables["cold_outlet_state"], + self.turbine.variables["inlet_state"], + ) - def base_params_single_optimize(self, fixed_var, base_params, target_var_name, bounds): - """ - 组件性能单变量优化器 - : 固定边界条件 - : 默认参数 - : 优化变量名称 - : 变量范围 - """ - - if target_var_name not in base_params: - raise ValueError(f"参数{target_var_name}不在参数字典中") - - def opt_fun(opt_var): - # 复制变量字典 - opt_param = base_params.copy() - # 修改要优化的变量为参数 - opt_param[target_var_name] = opt_var - # 带入循环参数计算 - try: - eff = self.RC( - T_low = fixed_var['T_low'], - T_high= fixed_var['T_high'], - p_low = fixed_var['p_low'], - p_high = fixed_var['p_high'], - ploss = fixed_var['ploss'], - param = opt_param) - return -eff - except Exception as e: - return 0.0 # 遇到物性计算崩溃时返回极差值 - res = minimize_scalar(opt_fun, bounds=bounds, method='bounded') - if res.success: - print(f"✅ 优化完成！") - print(f"👉 最佳 {target_var_name} = {res.x:.4f}") - print(f"👉 此时系统最高效率 = {-res.fun:.2%}\n") - else: - print("❌ 优化失败。") - - return res - - def fixed_params_single_optimize(self, fixed_params, target_var_name, params, bounds): - """ - 边界条件单变量优化器 - : 固定边界条件 - : 默认参数 - : 优化变量名称 - : 变量范围 - """ - if target_var_name not in fixed_params: - raise ValueError(f"参数{target_var_name}不在参数字典中") - - def opt_fun(opt_var): - # 复制变量 - opt_params = fixed_params.copy() - # 变量替换 - opt_params[target_var_name] = opt_var - # 带入循环 - try: - eff = self.RC(T_low = opt_params['T_low'], - T_high = opt_params['T_high'], - p_low = opt_params['p_low'], - p_high = opt_params['p_high'], - ploss = opt_params['ploss'], - param = params) - return -eff - except Exception as e: - return 0.0 - res = minimize_scalar(opt_fun, bounds=bounds, method='bounded') - if res.success: - print(f"✅ 优化完成！") - print(f"👉 最佳 {target_var_name} = {res.x:.4f}") - print(f"👉 此时系统最高效率 = {-res.fun:.2%}\n") - else: - print("❌ 优化失败。") - return res - - def plot_optimization_landscape(self, fixed_params, params, target_var_name, bounds, res, num_points=50): - """ - 绘制单变量优化地形图 - :param target_var_name: 要扫描和优化的变量名（如 'x'） - :param bounds: 扫描和优化的范围 (min, max) - :param num_points: 扫描的采样点数量，越大曲线越平滑，但计算越慢 - """ - print(f"开始对【{target_var_name}】进行区间扫描，共计算 {num_points} 个点...") - plt.rcParams['font.sans-serif'] = ['SimHei'] # Windows 用黑体 - plt.rcParams['axes.unicode_minus'] = False # 正常显示负号 - # 1. 生成扫描数组 - x_vals = np.linspace(bounds[0], bounds[1], num_points) - eff_vals = [] - valid_x = [] # 记录那些没有报错的 x - if target_var_name not in fixed_params: - # 2. 遍历计算曲线上的点 - for val in x_vals: - current_param = params.copy() - current_param[target_var_name] = val - try: - # 调用你的黑盒物理模型（注意：这里取正效率用于画图） - eff = self.RC( - T_low=fixed_params['T_low'], - T_high=fixed_params['T_high'], - p_low=fixed_params['p_low'], - p_high=fixed_params['p_high'], - ploss=fixed_params['ploss'], - param=current_param - ) - # 如果系统加了夹点校验且没通过，可能会返回 None 或者抛异常 - # 这里确保只有成功的点才画上去 - eff_vals.append(eff * 100) # 乘以 100 转换为百分比 - valid_x.append(val) - except Exception as e: - # 如果某个 x 导致计算崩溃，我们跳过这个点，不画它 - pass - else: - for val in x_vals: - current_param = fixed_params.copy() - current_param[target_var_name] = val - try: - # 调用你的黑盒物理模型（注意：这里取正效率用于画图） - eff = self.RC( - T_low=current_param['T_low'], - T_high=current_param['T_high'], - p_low=current_param['p_low'], - p_high=current_param['p_high'], - ploss=current_param['ploss'], - param=params - ) - # 如果系统加了夹点校验且没通过，可能会返回 None 或者抛异常 - # 这里确保只有成功的点才画上去 - eff_vals.append(eff * 100) # 乘以 100 转换为百分比 - valid_x.append(val) - except Exception as e: - # 如果某个 x 导致计算崩溃，我们跳过这个点，不画它 - pass - print("扫描完成！正在使用优化器寻找精确最高点...") - - # 4. 开始绘图 - plt.figure(figsize=(8, 6), dpi=120) # 设置画布大小和清晰度 - - # 画出目标函数曲线 - plt.plot(valid_x, eff_vals, linestyle='-', color='#1f77b4', linewidth=2, label='系统热效率曲线') - - # 如果优化成功，用醒目的红星标出最优点 - if res.success: - best_x = res.x - best_eff = -res.fun * 100 - plt.scatter(best_x, best_eff, color='red', marker='*', s=200, zorder=5, label=f'最优点 ({best_x:.4f}, {best_eff:.2f}%)') - - # 画辅助虚线对齐坐标轴 - plt.axvline(x=best_x, color='gray', linestyle='--', alpha=0.6) - plt.axhline(y=best_eff, color='gray', linestyle='--', alpha=0.6) - else: - print("优化结果未输入!") - # 设置图表装饰 - plt.title(f'系统热效率随 {target_var_name} 的变化趋势', fontsize=14) - plt.xlabel(f'优化变量: {target_var_name}', fontsize=12) - plt.ylabel('循环热效率 η (%)', fontsize=12) - plt.grid(True, linestyle=':', alpha=0.7) - plt.legend(fontsize=11) - - # 显示图像 - plt.tight_layout() - plt.show() - - # 调用测试 - # plot_optimization_landscape('x', bounds=(0.6, 0.95), num_points=40) - - - def sweep_and_optimize(self, fixed_params, params, sweep_var, sweep_bounds, opt_var, opt_bounds, num_points=50): - """ - 带内部动态优化的单变量扫描器 (极度通用版) - - :param fixed_params: 固定的边界条件字典 - :param params: 组件性能参数字典 - :param sweep_var: 你要扫描/遍历的变量名 (例如 'T_low') - :param sweep_bounds: 扫描变量的范围 (min, max) - :param opt_var: 在每个扫描点下，你需要动态寻找最优值的变量名 (例如 'x') - :param opt_bounds: 优化变量的搜索范围 (min, max) - :param num_points: 扫描点数 - """ - print(f"\n🚀 开始执行嵌套扫描：") - print(f" - 扫描变量 (X轴): 【{sweep_var}】 范围 {sweep_bounds}") - print(f" - 内部动态优化变量: 【{opt_var}】 范围 {opt_bounds}") - - # 1. 生成扫描节点 - sweep_vals = np.linspace(sweep_bounds[0], sweep_bounds[1], num_points) - - # 记录数据的列表 - valid_sweep_vals = [] - best_effs = [] - best_opt_vals = [] - - # 2. 开始逐点扫描 - for s_val in sweep_vals: - - # 【核心1：每次必须使用干净的字典副本】 - current_fixed = fixed_params.copy() - current_param = params.copy() - - # 判断扫描变量是属于 fixed_params 还是 params，并赋值 - if sweep_var in current_fixed: - current_fixed[sweep_var] = s_val - elif sweep_var in current_param: - current_param[sweep_var] = s_val - else: - raise ValueError(f"找不到扫描变量: {sweep_var}") - - # 3. 定义内部优化目标函数 (闭包) - def inner_objective(guess_val): - # 将优化器猜的值赋给 opt_var - current_param[opt_var] = guess_val - - try: - # 调用黑盒计算 - eff = self.RC( - T_low=current_fixed['T_low'], - T_high=current_fixed['T_high'], - p_low=current_fixed['p_low'], - p_high=current_fixed['p_high'], - ploss=current_fixed['ploss'], - param=current_param - ) - # 【预留口：此处可加入换热器内部夹点校验】 - # 取出两个换热器进行夹点校验 - for rec in self.recuperator: - if not rec.check_pinch_point(self.property_calculator): - return 0.0 # 核心！如果交叉了，直接返回 0 效率，强迫优化器换参数 - return -eff - except Exception: - return 0.0 - - # 4. 调用一维优化器 - res = minimize_scalar(inner_objective, bounds=opt_bounds, method='bounded') - - # 5. 结果校验与存储 - if res.success and -res.fun > 0: - best_eff = -res.fun * 100 - best_opt_val = res.x - - valid_sweep_vals.append(s_val) - best_effs.append(best_eff) - best_opt_vals.append(best_opt_val) - - print(f"✔️ {sweep_var} = {s_val:.2f} | 寻得最优 {opt_var} = {best_opt_val:.4f} | 最高效率 = {best_eff:.2f}%") - else: - print(f"❌ {sweep_var} = {s_val:.2f} | 优化失败或物理无解，已跳过") - - # 6. 调用画图方法 (将画图剥离，保持代码干净) - self._plot_results(sweep_var, valid_sweep_vals, best_effs, opt_var, best_opt_vals) - - return valid_sweep_vals, best_effs, best_opt_vals - - def _plot_results(self, sweep_var, x_data, y_eff_data, opt_var, y_opt_data): - """专门用来画图的内部方法，支持双Y轴""" - if not x_data: - print("没有有效数据可供绘制！") - return - - plt.rcParams['font.sans-serif'] = ['SimHei'] - plt.rcParams['axes.unicode_minus'] = False - - fig, ax1 = plt.subplots(figsize=(9, 6), dpi=120) - - # 画左Y轴：最高效率曲线 - color1 = '#1f77b4' - ax1.set_xlabel(f'扫描变量: {sweep_var}', fontsize=12) - ax1.set_ylabel('最优循环热效率 η (%)', color=color1, fontsize=12) - ax1.plot(x_data, y_eff_data, color=color1, linewidth=2.5, label='系统热效率') - ax1.tick_params(axis='y', labelcolor=color1) - ax1.grid(True, linestyle=':', alpha=0.6) - - # 画右Y轴：对应的最优分流量走势 - ax2 = ax1.twinx() - color2 = '#d62728' - ax2.set_ylabel(f'匹配的最优动态变量: {opt_var}', color=color2, fontsize=12) - ax2.plot(x_data, y_opt_data, color=color2, linestyle='--', linewidth=2, label=f'最优 {opt_var} 值') - ax2.tick_params(axis='y', labelcolor=color2) - - plt.title(f'系统最高效率及对应的最优 {opt_var} 随 {sweep_var} 的变化', fontsize=14) - fig.tight_layout() - plt.show() + return self.cycle_eff_calculator() diff --git a/brayton_cycle/optimization.py b/brayton_cycle/optimization.py new file mode 100644 index 0000000..96eb176 --- /dev/null +++ b/brayton_cycle/optimization.py @@ -0,0 +1,343 @@ +# -*- coding: utf-8 -*- +"""Optimization helpers for recompression Brayton cycle studies.""" + +import numpy as np +import matplotlib.pyplot as plt +from scipy.optimize import minimize_scalar + +from .cycles import BraytonCycle +from .sensitivity import RC_FIXED_KEYS, RC_PARAM_KEYS + + +INVALID_OBJECTIVE = 1.0e12 + + +def _require_keys(data, required_keys, data_name): + missing = [key for key in required_keys if key not in data] + if missing: + missing_text = ", ".join(missing) + raise ValueError(f"{data_name} missing required keys: {missing_text}") + + +def _resolve_variable_location(variable_name, fixed_params, params): + in_fixed = variable_name in fixed_params + in_params = variable_name in params + + if in_fixed and in_params: + raise ValueError( + f"{variable_name!r} exists in both fixed_params and params; " + "rename one of them or choose the target explicitly." + ) + if in_fixed: + return "fixed" + if in_params: + return "params" + raise ValueError(f"Unknown variable: {variable_name}") + + +def _with_updated_value(fixed_params, params, variable_name, value): + fixed = dict(fixed_params) + cycle_params = dict(params) + location = _resolve_variable_location(variable_name, fixed, cycle_params) + + if location == "fixed": + fixed[variable_name] = value + else: + cycle_params[variable_name] = value + + return fixed, cycle_params + + +def _validate_bounds(bounds, bounds_name): + if len(bounds) != 2: + raise ValueError(f"{bounds_name} must contain exactly two values") + if bounds[0] >= bounds[1]: + raise ValueError(f"{bounds_name} lower bound must be smaller than upper bound") + + +def _validate_num_points(num_points): + if num_points < 1: + raise ValueError("num_points must be at least 1") + + +def _evaluate_rc_cycle(fixed_params, params, refprop_path=None): + fixed = dict(fixed_params) + cycle_params = dict(params) + + _require_keys(fixed, RC_FIXED_KEYS, "fixed_params") + _require_keys(cycle_params, RC_PARAM_KEYS, "params") + + cycle_kwargs = {"name": "rc optimization evaluation"} + if refprop_path is not None: + cycle_kwargs["refprop_path"] = refprop_path + + cycle = BraytonCycle(**cycle_kwargs) + efficiency = cycle.RC( + T_low=fixed["T_low"], + T_high=fixed["T_high"], + p_low=fixed["p_low"], + p_high=fixed["p_high"], + ploss=fixed["ploss"], + param=cycle_params, + ) + return efficiency, cycle + + +def _pinch_points_ok(cycle): + if not cycle.recuperator: + return True + return all( + recuperator.check_pinch_point(cycle.property_calculator) + for recuperator in cycle.recuperator + ) + + +def _rc_objective(fixed_params, params, refprop_path=None, check_pinch=False): + try: + efficiency, cycle = _evaluate_rc_cycle(fixed_params, params, refprop_path) + if check_pinch and not _pinch_points_ok(cycle): + return INVALID_OBJECTIVE + return -efficiency + except Exception: + return INVALID_OBJECTIVE + + +def _annotate_result(result): + result.valid = bool( + result.success + and np.isfinite(result.fun) + and result.fun < INVALID_OBJECTIVE / 2 + ) + result.best_efficiency = -result.fun if result.valid else np.nan + return result + + +def optimize_rc_param( + fixed_params, + params, + target_var_name, + bounds, + refprop_path=None, + check_pinch=False, +): + """Optimize one RC component/cycle parameter for maximum efficiency.""" + if target_var_name not in params: + raise ValueError(f"{target_var_name!r} is not in params") + _validate_bounds(bounds, "bounds") + + def objective(value): + trial_params = dict(params) + trial_params[target_var_name] = value + return _rc_objective( + fixed_params, + trial_params, + refprop_path=refprop_path, + check_pinch=check_pinch, + ) + + result = minimize_scalar(objective, bounds=bounds, method="bounded") + return _annotate_result(result) + + +def optimize_rc_fixed_param( + fixed_params, + params, + target_var_name, + bounds, + refprop_path=None, + check_pinch=False, +): + """Optimize one RC boundary-condition parameter for maximum efficiency.""" + if target_var_name not in fixed_params: + raise ValueError(f"{target_var_name!r} is not in fixed_params") + _validate_bounds(bounds, "bounds") + + def objective(value): + trial_fixed = dict(fixed_params) + trial_fixed[target_var_name] = value + return _rc_objective( + trial_fixed, + params, + refprop_path=refprop_path, + check_pinch=check_pinch, + ) + + result = minimize_scalar(objective, bounds=bounds, method="bounded") + return _annotate_result(result) + + +def scan_rc_efficiency( + fixed_params, + params, + target_var_name, + bounds, + num_points=50, + refprop_path=None, +): + """Evaluate RC efficiency over a one-dimensional variable sweep.""" + _validate_bounds(bounds, "bounds") + _validate_num_points(num_points) + + x_values = [] + efficiencies = [] + errors = [] + + for value in np.linspace(bounds[0], bounds[1], num_points): + trial_fixed, trial_params = _with_updated_value( + fixed_params, params, target_var_name, value + ) + try: + efficiency, _ = _evaluate_rc_cycle( + trial_fixed, trial_params, refprop_path=refprop_path + ) + except Exception as exc: + errors.append((value, exc)) + continue + + x_values.append(value) + efficiencies.append(efficiency) + + return x_values, efficiencies, errors + + +def sweep_and_optimize_rc( + fixed_params, + params, + sweep_var, + sweep_bounds, + opt_var, + opt_bounds, + num_points=50, + refprop_path=None, + check_pinch=True, +): + """Sweep one variable and optimize another at each sweep point.""" + if sweep_var == opt_var: + raise ValueError("sweep_var and opt_var must be different variables") + _validate_bounds(sweep_bounds, "sweep_bounds") + _validate_bounds(opt_bounds, "opt_bounds") + _validate_num_points(num_points) + + valid_sweep_vals = [] + best_efficiencies = [] + best_opt_vals = [] + failures = [] + + for sweep_value in np.linspace(sweep_bounds[0], sweep_bounds[1], num_points): + current_fixed, current_params = _with_updated_value( + fixed_params, params, sweep_var, sweep_value + ) + + def objective(opt_value): + trial_fixed, trial_params = _with_updated_value( + current_fixed, current_params, opt_var, opt_value + ) + return _rc_objective( + trial_fixed, + trial_params, + refprop_path=refprop_path, + check_pinch=check_pinch, + ) + + result = minimize_scalar(objective, bounds=opt_bounds, method="bounded") + result = _annotate_result(result) + + if result.valid: + valid_sweep_vals.append(sweep_value) + best_efficiencies.append(result.best_efficiency) + best_opt_vals.append(result.x) + else: + failures.append((sweep_value, result)) + + return valid_sweep_vals, best_efficiencies, best_opt_vals, failures + + +def plot_optimization_landscape( + fixed_params, + params, + target_var_name, + bounds, + result=None, + num_points=50, + refprop_path=None, + show=True, +): + """Plot RC efficiency over a one-dimensional sweep.""" + x_values, efficiencies, errors = scan_rc_efficiency( + fixed_params, + params, + target_var_name, + bounds, + num_points=num_points, + refprop_path=refprop_path, + ) + efficiency_percent = [efficiency * 100 for efficiency in efficiencies] + + fig, ax = plt.subplots(figsize=(8, 6), dpi=120) + ax.plot(x_values, efficiency_percent, color="#1f77b4", linewidth=2) + ax.set_title(f"RC efficiency vs {target_var_name}") + ax.set_xlabel(target_var_name) + ax.set_ylabel("Cycle efficiency (%)") + ax.grid(True, linestyle=":", alpha=0.7) + + if result is not None and getattr(result, "valid", result.success): + best_efficiency = getattr(result, "best_efficiency", -result.fun) + ax.scatter( + result.x, + best_efficiency * 100, + color="red", + marker="*", + s=200, + zorder=5, + ) + ax.axvline(x=result.x, color="gray", linestyle="--", alpha=0.6) + ax.axhline(y=best_efficiency * 100, color="gray", linestyle="--", alpha=0.6) + + fig.tight_layout() + if show: + plt.show() + + return fig, ax, x_values, efficiency_percent, errors + + +def plot_sweep_optimization_results( + sweep_var, + sweep_values, + best_efficiencies, + opt_var, + best_opt_values, + show=True, +): + """Plot nested sweep and optimization results with two y axes.""" + if not sweep_values: + raise ValueError("No valid sweep data to plot") + + fig, ax1 = plt.subplots(figsize=(9, 6), dpi=120) + + color1 = "#1f77b4" + ax1.set_xlabel(sweep_var) + ax1.set_ylabel("Best cycle efficiency (%)", color=color1) + best_efficiency_percent = [ + efficiency * 100 for efficiency in best_efficiencies + ] + ax1.plot(sweep_values, best_efficiency_percent, color=color1, linewidth=2.5) + ax1.tick_params(axis="y", labelcolor=color1) + ax1.grid(True, linestyle=":", alpha=0.6) + + ax2 = ax1.twinx() + color2 = "#d62728" + ax2.set_ylabel(f"Best {opt_var}", color=color2) + ax2.plot( + sweep_values, + best_opt_values, + color=color2, + linestyle="--", + linewidth=2, + ) + ax2.tick_params(axis="y", labelcolor=color2) + + fig.tight_layout() + if show: + plt.show() + + return fig, ax1, ax2 diff --git a/brayton_cycle/sensitivity.py b/brayton_cycle/sensitivity.py new file mode 100644 index 0000000..8f80888 --- /dev/null +++ b/brayton_cycle/sensitivity.py @@ -0,0 +1,50 @@ +# -*- coding: utf-8 -*- +"""Evaluation helpers for cycle sensitivity analysis.""" + +from .cycles import BraytonCycle + + +RC_FIXED_KEYS = ("T_low", "T_high", "p_low", "p_high", "ploss") +RC_PARAM_KEYS = ( + "x", + "compressor_eff", + "recompressor_eff", + "turbine_eff", + "recuperator_eff", + "highT_recuperator_eff", +) + + +def _require_keys(data, required_keys, data_name): + missing = [key for key in required_keys if key not in data] + if missing: + missing_text = ", ".join(missing) + raise ValueError(f"{data_name} missing required keys: {missing_text}") + + +def evaluate_rc_efficiency(fixed_params, params, refprop_path=None): + """Return the recompression Brayton cycle thermal efficiency. + + The returned efficiency is a decimal value, for example 0.45 means 45%. + A new cycle instance is created for each evaluation so repeated sensitivity + runs do not reuse component state from previous cases. + """ + fixed = dict(fixed_params) + cycle_params = dict(params) + + _require_keys(fixed, RC_FIXED_KEYS, "fixed_params") + _require_keys(cycle_params, RC_PARAM_KEYS, "params") + + cycle_kwargs = {"name": "rc efficiency evaluation"} + if refprop_path is not None: + cycle_kwargs["refprop_path"] = refprop_path + + cycle = BraytonCycle(**cycle_kwargs) + return cycle.RC( + T_low=fixed["T_low"], + T_high=fixed["T_high"], + p_low=fixed["p_low"], + p_high=fixed["p_high"], + ploss=fixed["ploss"], + param=cycle_params, + ) diff --git a/brayton_cycle_test.py b/brayton_cycle_test.py index b3e0d60..b6714c3 100644 --- a/brayton_cycle_test.py +++ b/brayton_cycle_test.py @@ -1,7 +1,11 @@ # -*- coding: utf-8 -*- """Run an example recompression CO2 Brayton cycle sweep.""" -from brayton_cycle import BraytonCycle +from brayton_cycle import ( + BraytonCycle, + plot_sweep_optimization_results, + sweep_and_optimize_rc, +) def run_rc_sweep(): @@ -42,7 +46,7 @@ def run_rc_sweep(): print(Q_input + Wc - (Q_output + Wt)) print((Wt - Wc) / Q_input) - x_vals, eff_vals, opt_vals = brayton.sweep_and_optimize( + x_vals, eff_vals, opt_vals, failures = sweep_and_optimize_rc( fixed_params=fixed_var, params=param, sweep_var="T_low", @@ -51,6 +55,16 @@ def run_rc_sweep(): opt_bounds=(0.1, 0.5), num_points=30, ) + if failures: + print(f"Skipped {len(failures)} invalid sweep points.") + + plot_sweep_optimization_results( + "T_low", + x_vals, + eff_vals, + "x", + opt_vals, + ) return brayton, x_vals, eff_vals, opt_vals