循环单一敏感性分析与单变量优化demo

.gitignore

@@ -3,3 +3,4 @@
 __pycache__/
 *.pyc
 *.pyo
+examples/output/

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",

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,
+    )

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()

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()