Plot Generator Operating Window from VDEF Parameter Sweep.

This script loads results from the VDEF parameter sweep and visualizes the generator operating window with actual design points overlaid.

Usage:

python plot_generator_operating_window.py

Requirements:

• Run example_vdef_mdao.py first to generate parameter_sweep_cases.sql

Database file not found: parameter_sweep_cases.sql
Please run example_vdef_mdao.py first to generate the parameter sweep data.

To generate the data, run:
  python example_vdef_mdao.py

import sys
from pathlib import Path

import matplotlib.pyplot as plt
import numpy as np

try:
    import openmdao.api as om
except ImportError:
    print("Error: OpenMDAO is required. Install with: pip install openmdao")
    sys.exit(1)

# Add parent directory to path to import paroto
try:
    # When run as a script
    sys.path.insert(0, str(Path(__file__).parent.parent.parent))
except NameError:
    # When run by sphinx_gallery or in interactive mode
    sys.path.insert(0, str(Path.cwd().parent.parent))

from paroto.utils.plotting import add_design_points_to_operating_window, plot_operating_window


def load_parameter_sweep_results(db_file="parameter_sweep_cases.sql"):
    """Load results from OpenMDAO case recorder database.

    Parameters
    ----------
    db_file : str
        Path to SQLite database file

    Returns
    -------
    frequencies : np.ndarray
        Pulse frequencies (Hz)
    voltages : np.ndarray
        HV voltages (V)
    hv_constraints : np.ndarray
        HV operating window constraint values
    all_feasible : np.ndarray
        Boolean array indicating if all constraints are satisfied
    """
    db_path = Path(db_file)
    if not db_path.exists():
        raise FileNotFoundError(
            f"Database file not found: {db_file}\n"
            "Please run example_vdef_mdao.py first to generate the parameter sweep data."
        )

    # Load cases from database
    cr = om.CaseReader(db_file)
    cases = cr.list_cases("driver")

    frequencies = []
    voltages = []
    hv_constraints = []
    power_errors = []
    breakdown_margins = []

    print(f"Loading {len(cases)} cases from {db_file}...")

    for case_id in cases:
        try:
            case = cr.get_case(case_id)

            # Extract variables
            freq = case.get_val("pulse_frequency")[0]
            volt = case.get_val("hv_voltage")[0]
            hv_const = case.get_val("hv_operating_window_satisfied")[0]
            power_err = abs(case.get_val("power_check.power_error")[0])
            breakdown = case.get_val("breakdown_margin")[0]

            frequencies.append(freq)
            voltages.append(volt)
            hv_constraints.append(hv_const)
            power_errors.append(power_err)
            breakdown_margins.append(breakdown)

        except Exception as e:
            print(f"Warning: Failed to load case {case_id}: {e}")
            continue

    # Convert to arrays
    frequencies = np.array(frequencies)
    voltages = np.array(voltages)
    hv_constraints = np.array(hv_constraints)
    power_errors = np.array(power_errors)
    breakdown_margins = np.array(breakdown_margins)

    # Determine overall feasibility (all constraints satisfied)
    power_tolerance = 0.10
    all_feasible = (
        (power_errors < power_tolerance) & (breakdown_margins > 0) & (hv_constraints >= 1.0)
    )

    print(f"Loaded {len(frequencies)} cases successfully")
    print(f"  Feasible points (all constraints): {np.sum(all_feasible)}")
    print(f"  HV constraint violations: {np.sum(hv_constraints < 1.0)}")

    return frequencies, voltages, hv_constraints, all_feasible


def main():
    """Plot operating window with VDEF results."""
    # Load parameter sweep results
    try:
        frequencies, voltages, hv_constraints, all_feasible = load_parameter_sweep_results()
    except FileNotFoundError as e:
        print(f"\n{e}")
        print("\nTo generate the data, run:")
        print("  python example_vdef_mdao.py")
        return

    # Create figure with two subplots
    fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(16, 7))

    # ===== Plot 1: HV Constraint Only =====
    print("\nGenerating Plot 1: HV Constraint...")

    # Plot operating window
    f_boundary, V_boundary = plot_operating_window(
        ax1,
        Z=50.0,
        t_pulse=20e-9,
        P_max=1000.0,
        f_range=(0, 250000),
        V_range=(0, 85000),
        show_annotations=True,
    )

    # Overlay design points colored by HV constraint only
    add_design_points_to_operating_window(
        ax1, frequencies, voltages, hv_constraints, label_prefix="VDEF Design"
    )

    ax1.set_xlabel("Frequency (Hz)", fontsize=12)
    ax1.set_ylabel("HV Voltage (V)", fontsize=12)
    ax1.set_title(
        "Generator Operating Window\n(HV Constraint Only)", fontsize=13, fontweight="bold"
    )

    # ===== Plot 2: All Constraints =====
    print("Generating Plot 2: All Constraints...")

    # Plot operating window again
    f_boundary2, V_boundary2 = plot_operating_window(
        ax2,
        Z=50.0,
        t_pulse=20e-9,
        P_max=1000.0,
        f_range=(0, 250000),
        V_range=(0, 85000),
        show_annotations=True,
    )

    # Overlay design points colored by overall feasibility
    # Separate feasible and infeasible
    if np.any(all_feasible):
        ax2.scatter(
            frequencies[all_feasible],
            voltages[all_feasible],
            c="darkgreen",
            marker="o",
            s=50,
            alpha=0.7,
            edgecolors="green",
            linewidth=1.5,
            label="All Constraints Satisfied",
            zorder=5,
        )

    if np.any(~all_feasible):
        ax2.scatter(
            frequencies[~all_feasible],
            voltages[~all_feasible],
            c="darkred",
            marker="x",
            s=80,
            alpha=0.6,
            linewidth=2,
            label="At Least One Constraint Violated",
            zorder=5,
        )

    ax2.set_xlabel("Frequency (Hz)", fontsize=12)
    ax2.set_ylabel("HV Voltage (V)", fontsize=12)
    ax2.set_title(
        "Generator Operating Window\n(All Constraints: Power + Breakdown + HV)",
        fontsize=13,
        fontweight="bold",
    )
    ax2.legend(loc="upper right")

    plt.tight_layout()

    # Save figure
    output_file = "generator_operating_window.png"
    plt.savefig(output_file, dpi=150, bbox_inches="tight")
    print(f"\nFigure saved to: {output_file}")

    plt.show()

    print("\nDone!")


if __name__ == "__main__":
    main()