lujingze/Brayton-Cycle-Optimization
1
0 Fork 0
Code Issues Pull Requests Packages Projects Releases Wiki Activity

架构调整和敏感性分析代码附加

Browse Source
This commit is contained in:
ljz committed 2026-06-20 14:32:13 +08:00
1 parent 1d09599e8c
commit 469567c3db
5 files changed
+621 -429
Download Patch File Download Diff File Expand all files Collapse all files

No files matched your search

brayton_cycle/cycles.py
+196 -427
View File
@@ -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()