Brayton-Cycle-Optimization/brayton_cycle/optimization.py

# -*- coding: utf-8 -*-
"""Optimization helpers for recompression Brayton cycle studies."""

import numpy as np
import matplotlib.pyplot as plt
from scipy.optimize import minimize_scalar

from .cycles import BraytonCycle
from .sensitivity import RC_FIXED_KEYS, RC_PARAM_KEYS


INVALID_OBJECTIVE = 1.0e12


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 _resolve_variable_location(variable_name, fixed_params, params):
    in_fixed = variable_name in fixed_params
    in_params = variable_name in params

    if in_fixed and in_params:
        raise ValueError(
            f"{variable_name!r} exists in both fixed_params and params; "
            "rename one of them or choose the target explicitly."
        )
    if in_fixed:
        return "fixed"
    if in_params:
        return "params"
    raise ValueError(f"Unknown variable: {variable_name}")


def _with_updated_value(fixed_params, params, variable_name, value):
    fixed = dict(fixed_params)
    cycle_params = dict(params)
    location = _resolve_variable_location(variable_name, fixed, cycle_params)

    if location == "fixed":
        fixed[variable_name] = value
    else:
        cycle_params[variable_name] = value

    return fixed, cycle_params


def _validate_bounds(bounds, bounds_name):
    if len(bounds) != 2:
        raise ValueError(f"{bounds_name} must contain exactly two values")
    if bounds[0] >= bounds[1]:
        raise ValueError(f"{bounds_name} lower bound must be smaller than upper bound")


def _validate_num_points(num_points):
    if num_points < 1:
        raise ValueError("num_points must be at least 1")


def _evaluate_rc_cycle(fixed_params, params, refprop_path=None):
    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 optimization evaluation"}
    if refprop_path is not None:
        cycle_kwargs["refprop_path"] = refprop_path

    cycle = BraytonCycle(**cycle_kwargs)
    efficiency = 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,
    )
    return efficiency, cycle


def _pinch_points_ok(cycle):
    if not cycle.recuperator:
        return True
    return all(
        recuperator.check_pinch_point(cycle.property_calculator)
        for recuperator in cycle.recuperator
    )


def _rc_objective(fixed_params, params, refprop_path=None, check_pinch=False):
    try:
        efficiency, cycle = _evaluate_rc_cycle(fixed_params, params, refprop_path)
        if check_pinch and not _pinch_points_ok(cycle):
            return INVALID_OBJECTIVE
        return -efficiency
    except Exception:
        return INVALID_OBJECTIVE


def _annotate_result(result):
    result.valid = bool(
        result.success
        and np.isfinite(result.fun)
        and result.fun < INVALID_OBJECTIVE / 2
    )
    result.best_efficiency = -result.fun if result.valid else np.nan
    return result


def optimize_rc_param(
    fixed_params,
    params,
    target_var_name,
    bounds,
    refprop_path=None,
    check_pinch=False,
):
    """Optimize one RC component/cycle parameter for maximum efficiency."""
    if target_var_name not in params:
        raise ValueError(f"{target_var_name!r} is not in params")
    _validate_bounds(bounds, "bounds")

    def objective(value):
        trial_params = dict(params)
        trial_params[target_var_name] = value
        return _rc_objective(
            fixed_params,
            trial_params,
            refprop_path=refprop_path,
            check_pinch=check_pinch,
        )

    result = minimize_scalar(objective, bounds=bounds, method="bounded")
    return _annotate_result(result)


def optimize_rc_fixed_param(
    fixed_params,
    params,
    target_var_name,
    bounds,
    refprop_path=None,
    check_pinch=False,
):
    """Optimize one RC boundary-condition parameter for maximum efficiency."""
    if target_var_name not in fixed_params:
        raise ValueError(f"{target_var_name!r} is not in fixed_params")
    _validate_bounds(bounds, "bounds")

    def objective(value):
        trial_fixed = dict(fixed_params)
        trial_fixed[target_var_name] = value
        return _rc_objective(
            trial_fixed,
            params,
            refprop_path=refprop_path,
            check_pinch=check_pinch,
        )

    result = minimize_scalar(objective, bounds=bounds, method="bounded")
    return _annotate_result(result)


def scan_rc_efficiency(
    fixed_params,
    params,
    target_var_name,
    bounds,
    num_points=50,
    refprop_path=None,
):
    """Evaluate RC efficiency over a one-dimensional variable sweep."""
    _validate_bounds(bounds, "bounds")
    _validate_num_points(num_points)

    x_values = []
    efficiencies = []
    errors = []

    for value in np.linspace(bounds[0], bounds[1], num_points):
        trial_fixed, trial_params = _with_updated_value(
            fixed_params, params, target_var_name, value
        )
        try:
            efficiency, _ = _evaluate_rc_cycle(
                trial_fixed, trial_params, refprop_path=refprop_path
            )
        except Exception as exc:
            errors.append((value, exc))
            continue

        x_values.append(value)
        efficiencies.append(efficiency)

    return x_values, efficiencies, errors


def sweep_and_optimize_rc(
    fixed_params,
    params,
    sweep_var,
    sweep_bounds,
    opt_var,
    opt_bounds,
    num_points=50,
    refprop_path=None,
    check_pinch=True,
):
    """Sweep one variable and optimize another at each sweep point."""
    if sweep_var == opt_var:
        raise ValueError("sweep_var and opt_var must be different variables")
    _validate_bounds(sweep_bounds, "sweep_bounds")
    _validate_bounds(opt_bounds, "opt_bounds")
    _validate_num_points(num_points)

    valid_sweep_vals = []
    best_efficiencies = []
    best_opt_vals = []
    failures = []

    for sweep_value in np.linspace(sweep_bounds[0], sweep_bounds[1], num_points):
        current_fixed, current_params = _with_updated_value(
            fixed_params, params, sweep_var, sweep_value
        )

        def objective(opt_value):
            trial_fixed, trial_params = _with_updated_value(
                current_fixed, current_params, opt_var, opt_value
            )
            return _rc_objective(
                trial_fixed,
                trial_params,
                refprop_path=refprop_path,
                check_pinch=check_pinch,
            )

        result = minimize_scalar(objective, bounds=opt_bounds, method="bounded")
        result = _annotate_result(result)

        if result.valid:
            valid_sweep_vals.append(sweep_value)
            best_efficiencies.append(result.best_efficiency)
            best_opt_vals.append(result.x)
        else:
            failures.append((sweep_value, result))

    return valid_sweep_vals, best_efficiencies, best_opt_vals, failures


def plot_optimization_landscape(
    fixed_params,
    params,
    target_var_name,
    bounds,
    result=None,
    num_points=50,
    refprop_path=None,
    show=True,
):
    """Plot RC efficiency over a one-dimensional sweep."""
    x_values, efficiencies, errors = scan_rc_efficiency(
        fixed_params,
        params,
        target_var_name,
        bounds,
        num_points=num_points,
        refprop_path=refprop_path,
    )
    efficiency_percent = [efficiency * 100 for efficiency in efficiencies]

    fig, ax = plt.subplots(figsize=(8, 6), dpi=120)
    ax.plot(x_values, efficiency_percent, color="#1f77b4", linewidth=2)
    ax.set_title(f"RC efficiency vs {target_var_name}")
    ax.set_xlabel(target_var_name)
    ax.set_ylabel("Cycle efficiency (%)")
    ax.grid(True, linestyle=":", alpha=0.7)

    if result is not None and getattr(result, "valid", result.success):
        best_efficiency = getattr(result, "best_efficiency", -result.fun)
        ax.scatter(
            result.x,
            best_efficiency * 100,
            color="red",
            marker="*",
            s=200,
            zorder=5,
        )
        ax.axvline(x=result.x, color="gray", linestyle="--", alpha=0.6)
        ax.axhline(y=best_efficiency * 100, color="gray", linestyle="--", alpha=0.6)

    fig.tight_layout()
    if show:
        plt.show()

    return fig, ax, x_values, efficiency_percent, errors


def plot_sweep_optimization_results(
    sweep_var,
    sweep_values,
    best_efficiencies,
    opt_var,
    best_opt_values,
    show=True,
):
    """Plot nested sweep and optimization results with two y axes."""
    if not sweep_values:
        raise ValueError("No valid sweep data to plot")

    fig, ax1 = plt.subplots(figsize=(9, 6), dpi=120)

    color1 = "#1f77b4"
    ax1.set_xlabel(sweep_var)
    ax1.set_ylabel("Best cycle efficiency (%)", color=color1)
    best_efficiency_percent = [
        efficiency * 100 for efficiency in best_efficiencies
    ]
    ax1.plot(sweep_values, best_efficiency_percent, color=color1, linewidth=2.5)
    ax1.tick_params(axis="y", labelcolor=color1)
    ax1.grid(True, linestyle=":", alpha=0.6)

    ax2 = ax1.twinx()
    color2 = "#d62728"
    ax2.set_ylabel(f"Best {opt_var}", color=color2)
    ax2.plot(
        sweep_values,
        best_opt_values,
        color=color2,
        linestyle="--",
        linewidth=2,
    )
    ax2.tick_params(axis="y", labelcolor=color2)

    fig.tight_layout()
    if show:
        plt.show()

    return fig, ax1, ax2