Brayton-Cycle-Optimization/brayton_cycle/sensitivity.py

# -*- 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",
)
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):
    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,
    )


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