Note
Go to the end to download the full example code.
Find working VedEf parameters through systematic search.
This script attempts to find a working set of parameters by: 1. Starting with nominal values 2. Testing individual parameter variations 3. Diagnosing which constraints are violated
import importlib.util
from pathlib import Path
import numpy as np
# Import setup_vedef_problem from the numbered file
spec = importlib.util.spec_from_file_location(
"problem_module", Path(__file__).parent / "1_setup_vedef_problem.py"
)
problem_module = importlib.util.module_from_spec(spec)
spec.loader.exec_module(problem_module)
setup_vedef_problem = problem_module.setup_vedef_problem
def test_single_case(V, e, d_t, f, E_p, verbose=True):
"""Test a single parameter case and report detailed diagnostics.
Parameters
----------
V : float
High voltage (V)
e : float
Electrode gap distance (m)
d_t : float
Torch diameter (m)
f : float
Pulse frequency (Hz)
E_p : float
Energy per pulse (J)
verbose : bool
Print detailed diagnostics
Returns
-------
success : bool
True if case ran successfully
results : dict or None
Results dictionary if successful
"""
if verbose:
print(
f"\nTesting: V={V / 1e3:.1f}kV, e={e * 1e3:.1f}mm, d_t={d_t * 1e3:.1f}mm, "
f"f={f / 1e3:.1f}kHz, E_p={E_p * 1e3:.1f}mJ"
)
try:
# Create problem
prob = setup_vedef_problem(target_power=25000.0, max_energy_density=1e9, as_subsystem=False)
# Set fixed parameters
prob.set_val("pulse_duration", 1e-6, units="s")
prob.set_val("TP_QM", 160.0 / 86400.0, units="kg/s")
prob.set_val("ThermalDiameterModel.gas_density", 0.717, units="kg/m**3")
prob.set_val("ThermalDiameterModel.gas_heat_capacity", 2200.0, units="J/(kg*K)")
prob.set_val("ThermalDiameterModel.thermal_conductivity", 0.08, units="W/(m*K)")
prob.set_val("InitialTemperatureModel.gas_density", 0.717, units="kg/m**3")
prob.set_val("BreakdownModel.gas_properties_pressure", 101325.0, units="Pa")
# Set design variables (using ARC002 naming convention)
prob.set_val("G_UMAX_OUT", V, units="V")
prob.set_val("G_e", e, units="m")
prob.set_val("TP_D_OUT", d_t, units="m")
prob.set_val("G_F", f, units="Hz")
prob.set_val("G_Ep", E_p, units="J")
# Run model
prob.run_model()
# Extract results
results = {
"T_0": prob.get_val("T_0", units="K")[0],
"T_max": prob.get_val("T_max", units="K")[0],
"breakdown_voltage": prob.get_val("breakdown_voltage", units="V")[0],
"breakdown_margin": prob.get_val("breakdown_margin", units="V")[0],
"coverage_fraction": prob.get_val("coverage_fraction")[0],
"coverage_margin": prob.get_val("coverage_margin")[0],
"energy_density": prob.get_val("energy_density", units="J/m**3")[0],
"density_margin": prob.get_val("density_margin", units="J/m**3")[0],
"flow_velocity": prob.get_val("flow_velocity", units="m/s")[0],
"G_PW_OUT": prob.get_val("G_PW_OUT", units="W")[0],
"power_error": prob.get_val("PowerConstraint.power_error")[0],
}
# Check for NaN/Inf
has_nan = any(not np.isfinite(v) for k, v in results.items())
if verbose:
print(f" T_0 = {results['T_0']:.1f} K")
print(f" T_max = {results['T_max']:.1f} K")
print(f" Breakdown margin = {results['breakdown_margin'] / 1e3:.1f} kV")
print(f" Coverage margin = {results['coverage_margin']:.3f}")
print(f" Density margin = {results['density_margin'] / 1e6:.1f} MJ/m³")
print(f" Flow velocity = {results['flow_velocity']:.2f} m/s")
print(f" Total power = {results['G_PW_OUT'] / 1e3:.1f} kW")
print(f" Power error = {results['power_error'] * 100:.1f}%")
if has_nan:
print(" [X] FAILED: NaN/Inf in results")
return False, None
# Check constraints
constraints_ok = True
if results["breakdown_margin"] < 0:
print(" [X] Breakdown constraint violated")
constraints_ok = False
if results["coverage_margin"] < 0:
print(" [X] Coverage constraint violated")
constraints_ok = False
if results["density_margin"] < 0:
print(" [X] Density constraint violated")
constraints_ok = False
if results["power_error"] > 0.1:
print(" [X] Power constraint violated")
constraints_ok = False
if constraints_ok:
print(" [OK] SUCCESS: All constraints satisfied!")
else:
print(" [!] Partial: Model ran but constraints violated")
return not has_nan, results
except Exception as e:
if verbose:
print(f" [X] EXCEPTION: {type(e).__name__}: {e}")
return False, None
def main():
"""Search for working parameters."""
print("=" * 70)
print(" VedEf Parameter Search")
print("=" * 70)
# Test 1: Nominal values (from defaults in vedef_group.py)
print("\n" + "-" * 70)
print("Test 1: Nominal values from model defaults")
print("-" * 70)
success, results = test_single_case(
V=20000.0,
e=0.01,
d_t=0.02,
f=50000.0,
E_p=0.01,
)
if success and results:
print("\n[OK] Found working nominal case!")
return
# Test 2: Try different power levels (f * E_p = 25000 W)
print("\n" + "-" * 70)
print("Test 2: Varying frequency to match target power")
print("-" * 70)
for E_p in [0.001, 0.005, 0.010, 0.020, 0.030]:
f_required = 25000.0 / E_p
if f_required > 200e3: # Skip if frequency too high
continue
success, results = test_single_case(
V=20000.0,
e=0.01,
d_t=0.02,
f=f_required,
E_p=E_p,
)
if success and results and results["power_error"] < 0.1:
print("\n[OK] Found working case with matched power!")
return
# Test 3: Scan voltage
print("\n" + "-" * 70)
print("Test 3: Voltage scan at E_p=5mJ")
print("-" * 70)
E_p = 0.005
f = 25000.0 / E_p # 5 MHz
for V in np.linspace(10e3, 60e3, 6):
success, results = test_single_case(
V=V,
e=0.01,
d_t=0.02,
f=f,
E_p=E_p,
)
if success and results:
print(f"\n[OK] Working voltage found: {V / 1e3:.1f} kV")
# Test 4: Scan gap distance
print("\n" + "-" * 70)
print("Test 4: Gap distance scan")
print("-" * 70)
E_p = 0.005
f = 25000.0 / E_p
for e in np.linspace(0.005, 0.025, 5):
success, results = test_single_case(
V=30000.0,
e=e,
d_t=0.02,
f=f,
E_p=E_p,
)
if success and results:
print(f"\n[OK] Working gap found: {e * 1e3:.1f} mm")
print("\n" + "=" * 70)
print(" Search complete")
print("=" * 70)
if __name__ == "__main__":
main()