Analyze power constraint violations to find near-feasible operating points.

import numpy as np

# Load results
results = np.load("../../sweep_output/operating_window_results.npz")

# Extract data
V = results["V"]
e = results["e"]
d_t = results["d_t"]
f = results["f"]
E_p = results["E_p"]
G_PW_OUT = results["G_PW_OUT"]  # Total power output (ARC002 naming)
power_error = results["power_error"]
converged = results["converged"]

print("=" * 70)
print("  Power Constraint Analysis")
print("=" * 70)

# Basic statistics
print(f"\nTotal power range: {G_PW_OUT.min() / 1e3:.1f} - {G_PW_OUT.max() / 1e3:.1f} kW")
print("Target power: 25.0 kW")
print(f"Power error range: {power_error.min():.3f} - {power_error.max():.3f}")

# Find points closest to target power (25 kW)
power_diff = np.abs(G_PW_OUT - 25000.0)
closest_indices = np.argsort(power_diff)[:10]

print("\n" + "-" * 70)
print("  Top 10 closest to 25 kW target:")
print("-" * 70)
header = (
    f"{'Rank':<6}{'V (kV)':<10}{'e (mm)':<10}{'d_t (mm)':<12}"
    f"{'f (kHz)':<12}{'E_p (mJ)':<12}{'Power (kW)':<12}"
)
print(header)
print("-" * 70)

for rank, idx in enumerate(closest_indices, 1):
    print(
        f"{rank:<6}{V[idx] / 1e3:<10.1f}{e[idx] * 1e3:<10.1f}{d_t[idx] * 1e3:<12.1f}"
        f"{f[idx] / 1e3:<12.1f}{E_p[idx] * 1e3:<12.1f}{G_PW_OUT[idx] / 1e3:<12.2f}"
    )

# Find minimum and maximum power points
min_idx = np.argmin(G_PW_OUT)
max_idx = np.argmax(G_PW_OUT)

print("\n" + "-" * 70)
print("  Minimum Power Point:")
print("-" * 70)
print(f"  V = {V[min_idx] / 1e3:.1f} kV")
print(f"  e = {e[min_idx] * 1e3:.1f} mm")
print(f"  d_t = {d_t[min_idx] * 1e3:.1f} mm")
print(f"  f = {f[min_idx] / 1e3:.1f} kHz")
print(f"  E_p = {E_p[min_idx] * 1e3:.1f} mJ")
print(f"  Total Power = {G_PW_OUT[min_idx] / 1e3:.2f} kW")

print("\n" + "-" * 70)
print("  Maximum Power Point:")
print("-" * 70)
print(f"  V = {V[max_idx] / 1e3:.1f} kV")
print(f"  e = {e[max_idx] * 1e3:.1f} mm")
print(f"  d_t = {d_t[max_idx] * 1e3:.1f} mm")
print(f"  f = {f[max_idx] / 1e3:.1f} kHz")
print(f"  E_p = {E_p[max_idx] * 1e3:.1f} mJ")
print(f"  Total Power = {G_PW_OUT[max_idx] / 1e3:.2f} kW")

# Analyze other constraints for best power matches
print("\n" + "-" * 70)
print("  Constraint Status for Top 3 Closest Points:")
print("-" * 70)

for rank, idx in enumerate(closest_indices[:3], 1):
    print(f"\nRank {rank}: Power = {G_PW_OUT[idx] / 1e3:.2f} kW")
    print(f"  Breakdown margin: {results['breakdown_margin'][idx]:.1f} V")
    print(f"  Coverage margin: {results['coverage_margin'][idx]:.3f}")
    print(f"  Density margin: {results['density_margin'][idx] / 1e6:.1f} MJ/m³")
    print(
        "  All other constraints: OK"
        if results["breakdown_margin"][idx] > 0
        and results["coverage_margin"][idx] > 0
        and results["density_margin"][idx] > 0
        else "  Some constraints violated"
    )

print("\n" + "=" * 70)
print("  Recommendation:")
print("=" * 70)
best_idx = closest_indices[0]
print("\nClosest operating point to 25 kW target:")
print(f"  V = {V[best_idx] / 1e3:.1f} kV")
print(f"  e = {e[best_idx] * 1e3:.1f} mm")
print(f"  d_t = {d_t[best_idx] * 1e3:.1f} mm")
print(f"  f = {f[best_idx] / 1e3:.1f} kHz")
print(f"  E_p = {E_p[best_idx] * 1e3:.1f} mJ")
print(f"  Actual Power = {G_PW_OUT[best_idx] / 1e3:.2f} kW")
print(f"  Power Error = {abs(G_PW_OUT[best_idx] - 25000) / 25000 * 100:.1f}%")
print()