# -*- 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 from .components import ( Compressor, Concentrator, Condenser, Heater, Recuperator, Turbine, ) from .properties import CO2PropertyCalculator class BraytonCycle(): """循环计算""" 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 self.heater = None self.condenser = None self.concentrator = None self.refprop_path = refprop_path 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']) 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']) # 计算循环参数 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 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']) # 计算循环参数 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 = [] 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) # 假定低温回热器出口参数并进行迭代 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 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 = { '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}") 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 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 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()