From 5407eebf2fe9190cb5867f4033874a40570452b1 Mon Sep 17 00:00:00 2001 From: ljz <425868052@qq.com> Date: Mon, 22 Jun 2026 15:11:17 +0800 Subject: [PATCH] =?UTF-8?q?=E5=BE=AA=E7=8E=AF=E5=8D=95=E4=B8=80=E6=95=8F?= =?UTF-8?q?=E6=84=9F=E6=80=A7=E5=88=86=E6=9E=90=E4=B8=8E=E5=8D=95=E5=8F=98?= =?UTF-8?q?=E9=87=8F=E4=BC=98=E5=8C=96demo?= MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit --- .gitignore | 1 + brayton_cycle/__init__.py | 8 +- brayton_cycle/sensitivity.py | 142 ++++++++++++++++++ ..._component_performance_sensitivity_demo.py | 135 +++++++++++++++++ examples/rc_design_sensitivity_demo.py | 136 +++++++++++++++++ 5 files changed, 421 insertions(+), 1 deletion(-) create mode 100644 examples/rc_component_performance_sensitivity_demo.py create mode 100644 examples/rc_design_sensitivity_demo.py diff --git a/.gitignore b/.gitignore index a1506c7..ac96d40 100644 --- a/.gitignore +++ b/.gitignore @@ -3,3 +3,4 @@ __pycache__/ *.pyc *.pyo +examples/output/ diff --git a/brayton_cycle/__init__.py b/brayton_cycle/__init__.py index fda55e3..f809fdf 100644 --- a/brayton_cycle/__init__.py +++ b/brayton_cycle/__init__.py @@ -18,7 +18,11 @@ from .optimization import ( sweep_and_optimize_rc, ) from .properties import CO2PropertyCalculator -from .sensitivity import evaluate_rc_efficiency +from .sensitivity import ( + evaluate_rc_efficiency, + local_rc_component_performance_sensitivity, + local_rc_design_sensitivity, +) __all__ = [ "BraytonCycle", @@ -30,6 +34,8 @@ __all__ = [ "Recuperator", "Turbine", "evaluate_rc_efficiency", + "local_rc_component_performance_sensitivity", + "local_rc_design_sensitivity", "optimize_rc_fixed_param", "optimize_rc_param", "plot_optimization_landscape", diff --git a/brayton_cycle/sensitivity.py b/brayton_cycle/sensitivity.py index 8f80888..3300060 100644 --- a/brayton_cycle/sensitivity.py +++ b/brayton_cycle/sensitivity.py @@ -13,6 +13,14 @@ RC_PARAM_KEYS = ( "recuperator_eff", "highT_recuperator_eff", ) +RC_DESIGN_VARIABLES = ("T_low", "T_high", "p_low", "p_high", "ploss", "x") +RC_COMPONENT_PERFORMANCE_VARIABLES = ( + "compressor_eff", + "recompressor_eff", + "turbine_eff", + "recuperator_eff", + "highT_recuperator_eff", +) def _require_keys(data, required_keys, data_name): @@ -48,3 +56,137 @@ def evaluate_rc_efficiency(fixed_params, params, refprop_path=None): ploss=fixed["ploss"], param=cycle_params, ) + + +def _set_rc_variable(fixed_params, params, variable_name, value): + fixed = dict(fixed_params) + cycle_params = dict(params) + + in_fixed = variable_name in fixed + in_params = variable_name in cycle_params + if in_fixed and in_params: + raise ValueError(f"{variable_name!r} exists in both fixed_params and params") + if not in_fixed and not in_params: + raise ValueError(f"Unknown RC variable: {variable_name}") + + if in_fixed: + fixed[variable_name] = value + else: + cycle_params[variable_name] = value + + return fixed, cycle_params + + +def local_rc_design_sensitivity( + fixed_params, + params, + variables=None, + relative_step=0.01, + absolute_steps=None, + refprop_path=None, +): + """Run one-at-a-time local sensitivity analysis for RC design variables. + + The returned rows use decimal efficiency values. For example, 0.46 means + 46%. The normalized sensitivity is: + + ((eff_plus - eff_minus) / eff_base) + / ((value_plus - value_minus) / value_base) + + so variables with different units can be compared directly. + """ + if relative_step <= 0: + raise ValueError("relative_step must be positive") + + variables = variables or RC_DESIGN_VARIABLES + absolute_steps = absolute_steps or {} + + base_fixed = dict(fixed_params) + base_params = dict(params) + base_efficiency = evaluate_rc_efficiency( + base_fixed, + base_params, + refprop_path=refprop_path, + ) + + rows = [] + for variable_name in variables: + if variable_name in base_fixed: + base_value = base_fixed[variable_name] + elif variable_name in base_params: + base_value = base_params[variable_name] + else: + raise ValueError(f"Unknown RC variable: {variable_name}") + + step = absolute_steps.get(variable_name) + if step is None: + step = abs(base_value) * relative_step + if step <= 0: + raise ValueError(f"Step for {variable_name!r} must be positive") + + minus_value = base_value - step + plus_value = base_value + step + minus_fixed, minus_params = _set_rc_variable( + base_fixed, + base_params, + variable_name, + minus_value, + ) + plus_fixed, plus_params = _set_rc_variable( + base_fixed, + base_params, + variable_name, + plus_value, + ) + + eff_minus = evaluate_rc_efficiency( + minus_fixed, + minus_params, + refprop_path=refprop_path, + ) + eff_plus = evaluate_rc_efficiency( + plus_fixed, + plus_params, + refprop_path=refprop_path, + ) + + derivative = (eff_plus - eff_minus) / (plus_value - minus_value) + if base_value == 0 or base_efficiency == 0: + normalized_sensitivity = None + else: + normalized_sensitivity = derivative * base_value / base_efficiency + + rows.append( + { + "variable": variable_name, + "base_value": base_value, + "minus_value": minus_value, + "plus_value": plus_value, + "base_efficiency": base_efficiency, + "minus_efficiency": eff_minus, + "plus_efficiency": eff_plus, + "derivative": derivative, + "normalized_sensitivity": normalized_sensitivity, + } + ) + + return rows + + +def local_rc_component_performance_sensitivity( + fixed_params, + params, + variables=None, + relative_step=0.01, + absolute_steps=None, + refprop_path=None, +): + """Run one-at-a-time local sensitivity analysis for RC component performance.""" + return local_rc_design_sensitivity( + fixed_params=fixed_params, + params=params, + variables=variables or RC_COMPONENT_PERFORMANCE_VARIABLES, + relative_step=relative_step, + absolute_steps=absolute_steps, + refprop_path=refprop_path, + ) diff --git a/examples/rc_component_performance_sensitivity_demo.py b/examples/rc_component_performance_sensitivity_demo.py new file mode 100644 index 0000000..3449d41 --- /dev/null +++ b/examples/rc_component_performance_sensitivity_demo.py @@ -0,0 +1,135 @@ +# -*- coding: utf-8 -*- +"""Demo: one-at-a-time component-performance sensitivity for an RC cycle.""" + +from pathlib import Path +import csv +import sys + +import matplotlib.pyplot as plt + + +PROJECT_ROOT = Path(__file__).resolve().parents[1] +if str(PROJECT_ROOT) not in sys.path: + sys.path.insert(0, str(PROJECT_ROOT)) + +from brayton_cycle import local_rc_component_performance_sensitivity # noqa: E402 + + +COMPONENT_PERFORMANCE_VARIABLES = ( + "compressor_eff", + "recompressor_eff", + "turbine_eff", + "recuperator_eff", + "highT_recuperator_eff", +) + + +def build_base_case(): + fixed_params = { + "T_high": 650 + 273.15, + "T_low": 42 + 273.15, + "p_high": 20.0e3, + "p_low": 9.09e3, + "ploss": 0.01, + } + params = { + "x": 0.279, + "compressor_eff": 0.9, + "recompressor_eff": 0.9, + "turbine_eff": 0.93, + "recuperator_eff": 0.94, + "highT_recuperator_eff": 0.96, + } + return fixed_params, params + + +def sort_by_importance(rows): + return sorted( + rows, + key=lambda row: abs(row["normalized_sensitivity"] or 0.0), + reverse=True, + ) + + +def print_summary(rows): + print("RC component-performance local sensitivity") + print("Efficiency values are decimals; 0.46 means 46%.") + print() + print( + f"{'variable':<24} {'base':>10} {'eff-':>12} " + f"{'eff+':>12} {'norm_sens':>14}" + ) + print("-" * 80) + for row in sort_by_importance(rows): + sensitivity = row["normalized_sensitivity"] + sensitivity_text = "nan" if sensitivity is None else f"{sensitivity: .6f}" + print( + f"{row['variable']:<24} " + f"{row['base_value']:>10.6g} " + f"{row['minus_efficiency']:>12.6f} " + f"{row['plus_efficiency']:>12.6f} " + f"{sensitivity_text:>14}" + ) + + +def save_csv(rows, output_path): + output_path.parent.mkdir(parents=True, exist_ok=True) + with output_path.open("w", newline="", encoding="utf-8") as file: + writer = csv.DictWriter(file, fieldnames=list(rows[0].keys())) + writer.writeheader() + writer.writerows(rows) + + +def plot_sensitivity(rows, output_path=None, show=False): + sorted_rows = sort_by_importance(rows) + variables = [row["variable"] for row in sorted_rows] + sensitivities = [ + row["normalized_sensitivity"] or 0.0 + for row in sorted_rows + ] + colors = ["#1f77b4" if value >= 0 else "#d62728" for value in sensitivities] + + fig, ax = plt.subplots(figsize=(8, 4.8), dpi=130) + ax.barh(variables, sensitivities, color=colors) + ax.axvline(0.0, color="black", linewidth=0.8) + ax.set_xlabel("Normalized sensitivity") + ax.set_title("RC component-performance sensitivity") + ax.invert_yaxis() + fig.tight_layout() + + if output_path is not None: + output_path.parent.mkdir(parents=True, exist_ok=True) + fig.savefig(output_path) + if show: + plt.show() + + return fig, ax + + +def main(save_outputs=True, show_plot=False): + fixed_params, params = build_base_case() + rows = local_rc_component_performance_sensitivity( + fixed_params, + params, + variables=COMPONENT_PERFORMANCE_VARIABLES, + relative_step=0.01, + ) + + print_summary(rows) + + if save_outputs: + output_dir = PROJECT_ROOT / "examples" / "output" + save_csv(rows, output_dir / "rc_component_performance_sensitivity.csv") + plot_sensitivity( + rows, + output_path=output_dir / "rc_component_performance_sensitivity.png", + show=show_plot, + ) + elif show_plot: + plot_sensitivity(rows, show=True) + + return rows + + +if __name__ == "__main__": + main() diff --git a/examples/rc_design_sensitivity_demo.py b/examples/rc_design_sensitivity_demo.py new file mode 100644 index 0000000..28e42c4 --- /dev/null +++ b/examples/rc_design_sensitivity_demo.py @@ -0,0 +1,136 @@ +# -*- coding: utf-8 -*- +"""Demo: one-at-a-time design-parameter sensitivity for an RC Brayton cycle.""" + +from pathlib import Path +import csv +import sys + +import matplotlib.pyplot as plt + + +PROJECT_ROOT = Path(__file__).resolve().parents[1] +if str(PROJECT_ROOT) not in sys.path: + sys.path.insert(0, str(PROJECT_ROOT)) + +from brayton_cycle import local_rc_design_sensitivity # noqa: E402 + + +DESIGN_VARIABLES = ( + "T_low", + "T_high", + "p_low", + "p_high", + "ploss", + "x", +) + + +def build_base_case(): + fixed_params = { + "T_high": 650 + 273.15, + "T_low": 42 + 273.15, + "p_high": 20.0e3, + "p_low": 9.09e3, + "ploss": 0.01, + } + params = { + "x": 0.279, + "compressor_eff": 0.9, + "recompressor_eff": 0.9, + "turbine_eff": 0.93, + "recuperator_eff": 0.94, + "highT_recuperator_eff": 0.96, + } + return fixed_params, params + + +def sort_by_importance(rows): + return sorted( + rows, + key=lambda row: abs(row["normalized_sensitivity"] or 0.0), + reverse=True, + ) + + +def print_summary(rows): + print("RC design-parameter local sensitivity") + print("Efficiency values are decimals; 0.46 means 46%.") + print() + print( + f"{'variable':<10} {'base':>12} {'eff-':>12} " + f"{'eff+':>12} {'norm_sens':>14}" + ) + print("-" * 66) + for row in sort_by_importance(rows): + sensitivity = row["normalized_sensitivity"] + sensitivity_text = "nan" if sensitivity is None else f"{sensitivity: .6f}" + print( + f"{row['variable']:<10} " + f"{row['base_value']:>12.6g} " + f"{row['minus_efficiency']:>12.6f} " + f"{row['plus_efficiency']:>12.6f} " + f"{sensitivity_text:>14}" + ) + + +def save_csv(rows, output_path): + output_path.parent.mkdir(parents=True, exist_ok=True) + with output_path.open("w", newline="", encoding="utf-8") as file: + writer = csv.DictWriter(file, fieldnames=list(rows[0].keys())) + writer.writeheader() + writer.writerows(rows) + + +def plot_sensitivity(rows, output_path=None, show=False): + sorted_rows = sort_by_importance(rows) + variables = [row["variable"] for row in sorted_rows] + sensitivities = [ + row["normalized_sensitivity"] or 0.0 + for row in sorted_rows + ] + colors = ["#1f77b4" if value >= 0 else "#d62728" for value in sensitivities] + + fig, ax = plt.subplots(figsize=(8, 4.8), dpi=130) + ax.barh(variables, sensitivities, color=colors) + ax.axvline(0.0, color="black", linewidth=0.8) + ax.set_xlabel("Normalized sensitivity") + ax.set_title("RC design-parameter sensitivity") + ax.invert_yaxis() + fig.tight_layout() + + if output_path is not None: + output_path.parent.mkdir(parents=True, exist_ok=True) + fig.savefig(output_path) + if show: + plt.show() + + return fig, ax + + +def main(save_outputs=True, show_plot=False): + fixed_params, params = build_base_case() + rows = local_rc_design_sensitivity( + fixed_params, + params, + variables=DESIGN_VARIABLES, + relative_step=0.01, + ) + + print_summary(rows) + + if save_outputs: + output_dir = PROJECT_ROOT / "examples" / "output" + save_csv(rows, output_dir / "rc_design_sensitivity.csv") + plot_sensitivity( + rows, + output_path=output_dir / "rc_design_sensitivity.png", + show=show_plot, + ) + elif show_plot: + plot_sensitivity(rows, show=True) + + return rows + + +if __name__ == "__main__": + main()