Key Verbiage: Fault Tolerance
#!/usr/bin/env python3
"""
===============================================================================
BALANCED-ELECTRIC-5 — FULL SIX-PHASE MACHINE / LOAD / FAULT TEST
===============================================================================
Purpose
-------
A self-contained numerical A/B test between:
1) conventional balanced 3-phase
2) a six-phase machine with six independent phase branches
The test is deliberately conservative:
* same machine envelope
* same active copper volume
* same total copper resistance budget
* same DC/source voltage
* same conductor current limit
* same thermal limit
* same target mechanical output
* explicit electrical -> mechanical power accounting
* no voltage/current multiplication by simply adding waveforms
* no hidden free energy
* deterministic fault injection
The six-phase machine is NOT forced through the earlier A+F/B+E/C+D
pair-sum topology. That topology is retained only as an optional phasor
diagnostic. The main machine comparison preserves all six phases.
MODEL
-----
This is an engineering screening model, not an FEA package and not a
replacement for a manufactured-machine validation.
Each phase is represented by a sinusoidal winding coupled to an idealized
rotating rotor. The electromagnetic torque is obtained from instantaneous
phase current and the derivative of phase flux linkage with rotor angle:
e_k = - d(lambda_k)/dt
T_k = i_k * d(lambda_k)/d(theta)
The electrical model is:
V_k = R_k i_k + L_k di_k/dt + e_k
The mechanical model is:
J domega/dt = T_em - T_load - B omega
Thermal model:
Cth dT/dt = P_Cu - (T - Tambient)/Rth
The controller is an ideal current-regulated sinusoidal drive with a
voltage/current limiter. The same controller law and limits are used for
both machines.
For the six-phase machine, the six phase axes are:
0, 60, 120, 180, 240, 300 degrees
For the conventional machine:
0, 120, 240 degrees
A phase-fault is represented by opening a branch and commanding zero current.
The controller attempts to continue producing the requested torque with the
remaining healthy phases, subject to the same current/voltage/thermal limits.
The script reports both:
A) nominal equal-output operation
B) fault survival / derating
It also evaluates the earlier proposed pairing:
A+F, B+E, C+D
as a phasor-only diagnostic. It is NOT used as the six-phase machine's main
power path because simple voltage/current addition is not a physically valid
power-conserving machine connection by itself.
Dependencies
------------
Python 3.9+
numpy
matplotlib
Usage
-----
py balanced-electric5.py
Optional:
py balanced-electric5.py --seconds 0.60 --dt 2e-5
py balanced-electric5.py --seconds 1.00 --dt 1e-5 --fault-duration 0.20
py balanced-electric5.py --no-plots
py balanced-electric5.py --save-dir results
Outputs
-------
summary.txt
nominal.csv
faults.csv
harmonics.csv
energy_balance.csv
topology_pairing_diagnostic.csv
nominal_waveforms.png
torque_comparison.png
thermal_comparison.png
fault_survival.png
harmonic_spectrum.png
IMPORTANT
---------
Results are model results. Real motor performance requires electromagnetic
FEA, winding-factor/slotting data, inverter losses, iron-loss models,
mechanical measurements, insulation/protection design, and hardware tests.
===============================================================================
"""
from __future__ import annotations
import argparse
import csv
import math
import os
from dataclasses import dataclass
from pathlib import Path
from typing import Dict, List, Tuple
import numpy as np
try:
import matplotlib.pyplot as plt
except Exception:
plt = None
# =============================================================================
# CONFIGURATION
# =============================================================================
@dataclass
class Config:
# Simulation
seconds: float = 0.60
dt: float = 2.0e-5
sample_plot_every: int = 10
# Electrical source
f_e: float = 60.0
v_phase_rms: float = 120.0
current_limit: float = 18.0
voltage_limit: float = 240.0
# Machine target
pole_pairs: int = 2
target_speed_rpm: float = 1800.0
target_power_w: float = 5000.0
# Per-conductor / copper fairness
total_copper_volume: float = 1.0
copper_resistance_total_3ph: float = 0.90
copper_resistance_total_6ph: float = 0.90
# Winding parameters.
# These are screening values, scaled so the two systems have the same
# total copper budget.
inductance_phase_3ph: float = 0.012
inductance_phase_6ph: float = 0.006
flux_linkage_peak: float = 0.62
# Mechanical
inertia: float = 0.035
viscous_b: float = 0.003
# Thermal
ambient_c: float = 25.0
thermal_limit_c: float = 120.0
thermal_Rth_C_per_W: float = 0.045
thermal_Cth_J_per_C: float = 900.0
# Controller
current_loop_hz: float = 900.0
max_current_slew_a_per_s: float = 18000.0
speed_loop_hz: float = 25.0
speed_kp_torque: float = 0.045
speed_ki_torque: float = 0.65
# Fault testing
fault_duration_s: float = 0.20
fault_settle_s: float = 0.08
# Analysis
fft_harmonics_max: int = 30
rms_tail_fraction: float = 0.35
# =============================================================================
# BASIC UTILITIES
# =============================================================================
def rms(x: np.ndarray) -> float:
return float(np.sqrt(np.mean(np.square(x))))
def pct(x: float) -> float:
return 100.0 * x
def safe_div(a: float, b: float, default: float = 0.0) -> float:
if abs(b) < 1e-15:
return default
return a / b
def fmt(x: float, digits: int = 6) -> str:
if not math.isfinite(x):
return "nan"
return f"{x:.{digits}f}"
def phase_wrap_rad(x: float) -> float:
return (x + math.pi) % (2.0 * math.pi) - math.pi
def phasor_from_samples(x: np.ndarray, theta: np.ndarray) -> complex:
"""Fundamental coefficient against electrical angle theta."""
c = (2.0 / len(x)) * np.sum(x * np.exp(-1j * theta))
return complex(c)
def positive_sequence(ph: np.ndarray) -> complex:
"""abc positive sequence, x_a + a*x_b + a^2*x_c convention."""
a = np.exp(1j * 2.0 * math.pi / 3.0)
return (ph[0] + a * ph[1] + a * a * ph[2]) / 3.0
def negative_sequence(ph: np.ndarray) -> complex:
a = np.exp(1j * 2.0 * math.pi / 3.0)
return (ph[0] + a * a * ph[1] + a * ph[2]) / 3.0
def zero_sequence(ph: np.ndarray) -> complex:
return complex(np.mean(ph))
# =============================================================================
# MACHINE DEFINITION
# =============================================================================
@dataclass
class Machine:
name: str
axes: np.ndarray
resistance: np.ndarray
inductance: np.ndarray
flux_peak: float
current_limit: float
voltage_limit: float
pole_pairs: int
@property
def nph(self) -> int:
return len(self.axes)
@property
def copper_resistance_sum(self) -> float:
return float(np.sum(self.resistance))
@property
def copper_budget_relative(self) -> float:
return self.copper_resistance_sum
def make_machine(cfg: Config, nph: int) -> Machine:
axes = 2.0 * math.pi * np.arange(nph) / nph
# Equal total copper resistance budget.
# More phase branches => lower resistance per branch.
if nph == 3:
r = cfg.copper_resistance_total_3ph / nph
L = cfg.inductance_phase_3ph
else:
r = cfg.copper_resistance_total_6ph / nph
L = cfg.inductance_phase_6ph
return Machine(
name=f"{nph}-phase",
axes=axes,
resistance=np.full(nph, r, dtype=float),
inductance=np.full(nph, L, dtype=float),
flux_peak=cfg.flux_linkage_peak,
current_limit=cfg.current_limit,
voltage_limit=cfg.voltage_limit,
pole_pairs=cfg.pole_pairs,
)
# =============================================================================
# ELECTROMAGNETIC MODEL
# =============================================================================
def flux_linkage(machine: Machine, theta_m: float) -> np.ndarray:
"""
Sinusoidal phase flux linkage.
theta_m is mechanical rotor angle.
Electrical angle = pole_pairs * theta_m.
"""
theta_e = machine.pole_pairs * theta_m
return machine.flux_peak * np.cos(theta_e - machine.axes)
def dflux_dtheta(machine: Machine, theta_m: float) -> np.ndarray:
"""
d(lambda_k)/d(theta_m).
"""
theta_e = machine.pole_pairs * theta_m
return (
machine.flux_peak
* machine.pole_pairs
* (-np.sin(theta_e - machine.axes))
)
def back_emf(machine: Machine, theta_m: float, omega_m: float) -> np.ndarray:
# Terminal-voltage convention: v = R i + L di/dt + d(lambda)/dt.
# Therefore the source-side flux term is +d(lambda)/dt. The negative
# sign belongs to the separately defined induced/back-EMF quantity.
return dflux_dtheta(machine, theta_m) * omega_m
def electromagnetic_torque(machine: Machine, theta_m: float, current: np.ndarray) -> float:
return float(np.dot(current, dflux_dtheta(machine, theta_m)))
# =============================================================================
# CONTROL
# =============================================================================
@dataclass
class ControllerState:
speed_integral: float = 0.0
def desired_current_vector(
machine: Machine,
theta_m: float,
omega_m: float,
torque_cmd: float,
) -> np.ndarray:
"""
Minimum-norm current vector producing the requested instantaneous torque
under the sinusoidal flux model.
Since:
T = i dot d(lambda)/d(theta)
the minimum-norm solution is proportional to d(lambda)/d(theta).
"""
g = dflux_dtheta(machine, theta_m)
gg = float(np.dot(g, g))
if gg < 1e-20 or torque_cmd == 0.0:
return np.zeros(machine.nph)
i = (torque_cmd / gg) * g
# Enforce phase current limit.
peak = float(np.max(np.abs(i)))
if peak > machine.current_limit:
i *= machine.current_limit / peak
return i
def speed_controller(
cfg: Config,
controller: ControllerState,
omega_m: float,
target_omega: float,
base_torque: float,
dt: float,
torque_cap: float,
) -> float:
"""Speed regulator around the torque required by the commanded load.
The previous draft used only a speed-error torque command. That meant a
machine starting very near the target speed could initially receive almost
no torque while a 5 kW load was already applied. This version explicitly
supplies the steady-state load torque and uses the speed loop only as the
correction term.
"""
error = target_omega - omega_m
controller.speed_integral += error * dt
controller.speed_integral = float(
np.clip(controller.speed_integral, -torque_cap, torque_cap)
)
correction = (
cfg.speed_kp_torque * error
+ cfg.speed_ki_torque * controller.speed_integral
)
torque = base_torque + correction
return float(np.clip(torque, 0.0, torque_cap))
def load_torque_for_target_power(
cfg: Config,
omega_m: float,
) -> float:
if abs(omega_m) < 1e-9:
return 0.0
return cfg.target_power_w / abs(omega_m)
# =============================================================================
# PHYSICAL VOLTAGE LIMITING
# =============================================================================
def voltage_limited_current_command(
machine: Machine,
theta_m: float,
omega_m: float,
current_actual: np.ndarray,
current_target: np.ndarray,
dt: float,
) -> Tuple[np.ndarray, np.ndarray]:
"""Current controller with an explicit, power-consistent voltage limit.
The branch equation is evaluated as:
v = R*i_new + L*di/dt + d(lambda)/dt
with i_new = i + di*dt. The allowable di is solved directly from the
voltage inequality rather than applying a final voltage scaling that
would otherwise break the electrical energy balance.
"""
e = back_emf(machine, theta_m, omega_m)
tau = 1.0 / (2.0 * math.pi * 900.0)
raw_di = (current_target - current_actual) / max(tau, dt)
# v = R*(i + di*dt) + L*di + e
# = (R*i + e) + (L + R*dt)*di
base_v = machine.resistance * current_actual + e
coeff = machine.inductance + machine.resistance * dt
# Compute allowable derivative interval for each phase.
lo = (-machine.voltage_limit - base_v) / np.maximum(coeff, 1e-15)
hi = ( machine.voltage_limit - base_v) / np.maximum(coeff, 1e-15)
di = np.clip(raw_di, lo, hi)
i_new = current_actual + di * dt
i_new = np.clip(i_new, -machine.current_limit, machine.current_limit)
# Recompute di after the current hard limit so the reported voltage is
# exactly the voltage implied by the actual current transition.
di_actual = (i_new - current_actual) / max(dt, 1e-15)
v = machine.resistance * i_new + machine.inductance * di_actual + e
return i_new, v
# =============================================================================
# SIMULATION RESULT
# =============================================================================
@dataclass
class SimResult:
name: str
t: np.ndarray
theta: np.ndarray
omega: np.ndarray
speed_rpm: np.ndarray
current: np.ndarray
voltage: np.ndarray
torque: np.ndarray
torque_load: np.ndarray
copper_loss: np.ndarray
electrical_power: np.ndarray
mechanical_power: np.ndarray
thermal_c: np.ndarray
healthy: np.ndarray
fault_name: str
# =============================================================================
# MAIN SIMULATION
# =============================================================================
def simulate(
cfg: Config,
machine: Machine,
fault_phase: int | None = None,
fault_name: str = "none",
duration: float | None = None,
) -> SimResult:
if duration is None:
duration = cfg.seconds
n = int(math.ceil(duration / cfg.dt)) + 1
t = np.arange(n, dtype=float) * cfg.dt
theta = np.zeros(n)
omega = np.zeros(n)
speed_rpm = np.zeros(n)
current = np.zeros((n, machine.nph))
voltage = np.zeros((n, machine.nph))
torque = np.zeros(n)
torque_load = np.zeros(n)
copper_loss = np.zeros(n)
electrical_power = np.zeros(n)
mechanical_power = np.zeros(n)
thermal_c = np.zeros(n)
healthy = np.ones((n, machine.nph), dtype=bool)
target_omega = 2.0 * math.pi * cfg.target_speed_rpm / 60.0
# Start close to synchronous speed to make this an electrical-machine
# operating-point comparison rather than a startup contest.
omega[0] = target_omega * 0.985
thermal_c[0] = cfg.ambient_c
if fault_phase is not None:
healthy[:, fault_phase] = False
controller = ControllerState()
# Torque cap deliberately generous but finite.
torque_cap = max(
cfg.target_power_w / max(target_omega, 1e-9) * 2.5,
10.0,
)
for k in range(n - 1):
th = theta[k]
om = omega[k]
i_now = current[k]
# Load selected to represent target mechanical output at target speed.
# Add the modeled viscous-friction torque so the commanded torque
# corresponds to target shaft output rather than merely shaft+loss.
tl = load_torque_for_target_power(cfg, target_omega)
viscous_torque = cfg.viscous_b * target_omega
base_torque = tl + viscous_torque
torque_load[k] = tl
# Speed loop. The baseline is the load torque; the loop supplies only
# the correction needed to hold target speed.
t_cmd = speed_controller(
cfg,
controller,
om,
target_omega,
base_torque,
cfg.dt,
torque_cap,
)
# Fault: remove phase completely.
if fault_phase is not None:
# We do not simply delete the phase from the machine. Its branch
# remains electrically present but is opened, carrying zero current.
t_cmd = max(t_cmd, 0.0)
i_target = desired_current_vector(
machine,
th,
om,
t_cmd,
)
if fault_phase is not None:
i_target[fault_phase] = 0.0
i_new, v = voltage_limited_current_command(
machine,
th,
om,
i_now,
i_target,
cfg.dt,
)
if fault_phase is not None:
i_new[fault_phase] = 0.0
v[fault_phase] = 0.0
current[k + 1] = i_new
voltage[k] = v
# Recompute actual electromagnetic torque from actual current.
te = electromagnetic_torque(machine, th, i_new)
# If a branch is open, the controller cannot produce arbitrary torque.
# This is naturally reflected by the current vector.
torque[k] = te
# Mechanical dynamics.
accel = (
te
- tl
- cfg.viscous_b * om
) / cfg.inertia
omega[k + 1] = max(0.0, om + accel * cfg.dt)
theta[k + 1] = th + omega[k + 1] * cfg.dt
speed_rpm[k] = omega[k] * 60.0 / (2.0 * math.pi)
# Instantaneous electrical input.
pe = float(np.dot(v, i_new))
pc = float(np.sum(machine.resistance * np.square(i_new)))
pm = te * om
electrical_power[k] = pe
copper_loss[k] = pc
mechanical_power[k] = pm
# Lumped thermal model.
dT = (
pc
- (thermal_c[k] - cfg.ambient_c)
/ cfg.thermal_Rth_C_per_W
) / cfg.thermal_Cth_J_per_C
thermal_c[k + 1] = thermal_c[k] + dT * cfg.dt
speed_rpm[-1] = omega[-1] * 60.0 / (2.0 * math.pi)
torque_load[-1] = torque_load[-2]
torque[-1] = torque[-2]
voltage[-1] = voltage[-2]
copper_loss[-1] = copper_loss[-2]
electrical_power[-1] = electrical_power[-2]
mechanical_power[-1] = mechanical_power[-2]
return SimResult(
name=machine.name,
t=t,
theta=theta,
omega=omega,
speed_rpm=speed_rpm,
current=current,
voltage=voltage,
torque=torque,
torque_load=torque_load,
copper_loss=copper_loss,
electrical_power=electrical_power,
mechanical_power=mechanical_power,
thermal_c=thermal_c,
healthy=healthy,
fault_name=fault_name,
)
# =============================================================================
# STEADY-STATE METRICS
# =============================================================================
@dataclass
class Metrics:
name: str
speed_rpm: float
current_rms_total: float
current_rms_max: float
voltage_rms_total: float
voltage_peak_max: float
pin_w: float
pout_w: float
pcu_w: float
efficiency: float
pf: float
torque_mean_nm: float
torque_ripple_pct: float
thermal_c: float
thermal_margin_c: float
current_utilization_pct: float
voltage_utilization_pct: float
power_ripple_pct: float
energy_residual_pct: float
fault_name: str
def tail(x: np.ndarray, frac: float) -> np.ndarray:
start = int(len(x) * (1.0 - frac))
return x[max(start, 0):]
def calculate_metrics(
cfg: Config,
machine: Machine,
result: SimResult,
) -> Metrics:
sl = slice(int(len(result.t) * (1.0 - cfg.rms_tail_fraction)), None)
i_tail = result.current[sl]
v_tail = result.voltage[sl]
i_rms_each = np.sqrt(np.mean(i_tail * i_tail, axis=0))
v_rms_each = np.sqrt(np.mean(v_tail * v_tail, axis=0))
current_rms_total = float(np.sqrt(np.sum(i_rms_each ** 2)))
voltage_rms_total = float(np.sqrt(np.sum(v_rms_each ** 2)))
pin = float(np.mean(result.electrical_power[sl]))
pcu = float(np.mean(result.copper_loss[sl]))
pout = float(np.mean(result.mechanical_power[sl]))
apparent = voltage_rms_total * current_rms_total
pf = safe_div(pin, apparent)
tq = result.torque[sl]
tq_mean = float(np.mean(tq))
tq_pp = float(np.max(tq) - np.min(tq))
tq_ripple = safe_div(tq_pp, abs(tq_mean)) * 100.0
p = result.electrical_power[sl]
p_ripple = safe_div(
float(np.max(p) - np.min(p)),
max(abs(float(np.mean(p))), 1e-12),
) * 100.0
# Approximate conservation residual:
# Pin should approximately equal Pout + copper + mechanical viscous loss.
visc = float(np.mean(cfg.viscous_b * result.omega[sl] ** 2))
residual = pin - (pout + pcu + visc)
residual_pct = safe_div(abs(residual), max(abs(pin), 1.0)) * 100.0
thermal = float(np.max(result.thermal_c[sl]))
current_util = (
float(np.max(np.abs(i_tail))) / machine.current_limit * 100.0
)
voltage_util = (
float(np.max(np.abs(v_tail))) / machine.voltage_limit * 100.0
)
return Metrics(
name=machine.name,
speed_rpm=float(np.mean(result.speed_rpm[sl])),
current_rms_total=current_rms_total,
current_rms_max=float(np.max(i_rms_each)),
voltage_rms_total=voltage_rms_total,
voltage_peak_max=float(np.max(np.abs(v_tail))),
pin_w=pin,
pout_w=pout,
pcu_w=pcu,
efficiency=safe_div(pout, pin),
pf=pf,
torque_mean_nm=tq_mean,
torque_ripple_pct=tq_ripple,
thermal_c=thermal,
thermal_margin_c=cfg.thermal_limit_c - thermal,
current_utilization_pct=current_util,
voltage_utilization_pct=voltage_util,
power_ripple_pct=p_ripple,
energy_residual_pct=residual_pct,
fault_name=result.fault_name,
)
# =============================================================================
# HARMONIC ANALYSIS
# =============================================================================
def harmonic_amplitudes(
x: np.ndarray,
sample_rate: float,
fundamental_hz: float,
max_h: int,
) -> Dict[int, float]:
n = len(x)
window = np.hanning(n)
y = (x - np.mean(x)) * window
spec = np.fft.rfft(y)
freqs = np.fft.rfftfreq(n, 1.0 / sample_rate)
out = {}
for h in range(1, max_h + 1):
target = h * fundamental_hz
idx = int(np.argmin(np.abs(freqs - target)))
amp = 2.0 * abs(spec[idx]) / max(np.sum(window), 1e-12)
out[h] = float(amp)
return out
def harmonic_metrics(
cfg: Config,
result: SimResult,
) -> List[Tuple[str, int, float]]:
dt = cfg.dt
fs = 1.0 / dt
sl = slice(int(len(result.t) * (1.0 - cfg.rms_tail_fraction)), None)
i = result.current[sl]
rows = []
for phase in range(i.shape[1]):
amps = harmonic_amplitudes(
i[:, phase],
fs,
cfg.f_e,
cfg.fft_harmonics_max,
)
for h, amp in amps.items():
rows.append((result.name, phase, h, amp))
return rows
# =============================================================================
# SIX-PHASE / EARLIER PAIRING DIAGNOSTIC
# =============================================================================
def pairing_diagnostic() -> List[Dict[str, float | str]]:
"""
Reproduce the earlier geometric pairing without pretending that pairwise
voltage/current addition is a power-conserving machine connection.
Six source phases:
A=0 B=30 C=120 D=150 E=240 F=270
Proposed:
A+F, B+E, C+D
We calculate phasor amplitudes and sequence content only.
"""
names = ["A", "B", "C", "D", "E", "F"]
deg = np.array([0.0, 30.0, 120.0, 150.0, 240.0, 270.0])
z = np.exp(1j * np.deg2rad(deg))
pairs = [(0, 5), (1, 4), (2, 3)]
out = []
ph = np.array([z[a] + z[b] for a, b in pairs], dtype=complex)
# Normalize pair sum by sqrt(2), matching the earlier diagnostic.
ph /= math.sqrt(2.0)
p = positive_sequence(ph)
n = negative_sequence(ph)
z0 = zero_sequence(ph)
for idx, (a, b) in enumerate(pairs):
out.append({
"pair": f"{names[a]}+{names[b]}",
"amplitude": abs(ph[idx]),
"phase_deg": math.degrees(math.atan2(ph[idx].imag, ph[idx].real)) % 360.0,
"positive_sequence": abs(p),
"negative_sequence": abs(n),
"zero_sequence": abs(z0),
"negative_positive_pct": safe_div(abs(n), abs(p)) * 100.0,
})
return out
# =============================================================================
# FAULT TESTS
# =============================================================================
@dataclass
class FaultSummary:
architecture: str
fault: str
phase_index: int
nominal_pout_w: float
fault_pout_w: float
retained_output_pct: float
speed_rpm: float
torque_ripple_pct: float
max_current_util_pct: float
max_thermal_c: float
thermal_margin_c: float
protection_flag: str
def run_faults(
cfg: Config,
machine: Machine,
nominal_metrics: Metrics,
) -> List[FaultSummary]:
rows = []
for phase in range(machine.nph):
fault_name = f"{machine.name}_phase_{phase+1}_OPEN"
r = simulate(
cfg,
machine,
fault_phase=phase,
fault_name=fault_name,
duration=cfg.fault_duration_s,
)
m = calculate_metrics(cfg, machine, r)
retained = safe_div(m.pout_w, max(nominal_metrics.pout_w, 1e-9)) * 100.0
# Protection flag is descriptive, not a pass/fail claim about a real
# protection device.
if m.thermal_c >= cfg.thermal_limit_c:
protection = "THERMAL_LIMIT_REACHED"
elif m.current_utilization_pct >= 99.9:
protection = "CURRENT_LIMIT_ACTIVE"
elif m.voltage_utilization_pct >= 99.9:
protection = "VOLTAGE_LIMIT_ACTIVE"
else:
protection = "WITHIN_MODEL_LIMITS"
rows.append(
FaultSummary(
architecture=machine.name,
fault=fault_name,
phase_index=phase + 1,
nominal_pout_w=nominal_metrics.pout_w,
fault_pout_w=m.pout_w,
retained_output_pct=retained,
speed_rpm=m.speed_rpm,
torque_ripple_pct=m.torque_ripple_pct,
max_current_util_pct=m.current_utilization_pct,
max_thermal_c=m.thermal_c,
thermal_margin_c=m.thermal_margin_c,
protection_flag=protection,
)
)
return rows
# =============================================================================
# CSV OUTPUT
# =============================================================================
def write_csv(path: Path, rows: List[Dict | Tuple], headers: List[str] | None = None):
with path.open("w", newline="", encoding="utf-8") as f:
w = csv.writer(f)
if headers:
w.writerow(headers)
for row in rows:
if isinstance(row, dict):
w.writerow([row.get(h, "") for h in headers or row.keys()])
else:
w.writerow(list(row))
def write_nominal_csv(
path: Path,
results: List[SimResult],
every: int = 20,
):
rows = []
for r in results:
for k in range(0, len(r.t), every):
rows.append(
[
r.name,
r.t[k],
r.speed_rpm[k],
r.torque[k],
r.electrical_power[k],
r.mechanical_power[k],
r.copper_loss[k],
r.thermal_c[k],
*r.current[k].tolist(),
*r.voltage[k].tolist(),
]
)
max_ph = max(r.current.shape[1] for r in results)
headers = [
"architecture",
"time_s",
"speed_rpm",
"torque_Nm",
"electrical_power_W",
"mechanical_power_W",
"copper_loss_W",
"temperature_C",
]
for p in range(max_ph):
headers.append(f"I_phase_{p+1}_A")
for p in range(max_ph):
headers.append(f"V_phase_{p+1}_V")
# Pad shorter 3-phase rows.
padded = []
for row in rows:
base = row[:8]
vals = row[8:]
# row contains n current + n voltage.
# infer n from architecture.
nph = 3 if "3-phase" in row[0] else 6
ii = vals[:nph]
vv = vals[nph:]
padded.append(
base
+ ii + [""] * (max_ph - nph)
+ vv + [""] * (max_ph - nph)
)
with path.open("w", newline="", encoding="utf-8") as f:
w = csv.writer(f)
w.writerow(headers)
w.writerows(padded)
# =============================================================================
# REPORTING
# =============================================================================
def metrics_dict(m: Metrics) -> Dict[str, float | str]:
return {
"architecture": m.name,
"fault": m.fault_name,
"speed_rpm": m.speed_rpm,
"total_current_rms_A": m.current_rms_total,
"max_phase_current_rms_A": m.current_rms_max,
"total_voltage_rms_V": m.voltage_rms_total,
"max_voltage_V": m.voltage_peak_max,
"input_power_W": m.pin_w,
"output_power_W": m.pout_w,
"copper_loss_W": m.pcu_w,
"efficiency": m.efficiency,
"power_factor_model": m.pf,
"mean_torque_Nm": m.torque_mean_nm,
"torque_ripple_pct": m.torque_ripple_pct,
"thermal_C": m.thermal_c,
"thermal_margin_C": m.thermal_margin_c,
"current_utilization_pct": m.current_utilization_pct,
"voltage_utilization_pct": m.voltage_utilization_pct,
"power_ripple_pct": m.power_ripple_pct,
"energy_residual_pct": m.energy_residual_pct,
}
def print_metrics(m: Metrics):
print(f"\n[{m.name}] {m.fault_name}")
print(f" speed {m.speed_rpm:12.3f} rpm")
print(f" total RMS current {m.current_rms_total:12.4f} A")
print(f" max phase RMS current {m.current_rms_max:12.4f} A")
print(f" total RMS voltage {m.voltage_rms_total:12.4f} V")
print(f" max voltage {m.voltage_peak_max:12.4f} V")
print(f" input power {m.pin_w:12.4f} W")
print(f" mechanical output {m.pout_w:12.4f} W")
print(f" copper loss {m.pcu_w:12.4f} W")
print(f" efficiency {100*m.efficiency:12.4f} %")
print(f" model PF {m.pf:12.6f}")
print(f" mean torque {m.torque_mean_nm:12.4f} Nm")
print(f" torque ripple {m.torque_ripple_pct:12.4f} %")
print(f" max temperature {m.thermal_c:12.4f} C")
print(f" thermal margin {m.thermal_margin_c:12.4f} C")
print(f" current utilization {m.current_utilization_pct:12.4f} %")
print(f" voltage utilization {m.voltage_utilization_pct:12.4f} %")
print(f" input-power ripple {m.power_ripple_pct:12.4f} %")
print(f" energy residual {m.energy_residual_pct:12.6f} %")
def write_summary(
path: Path,
cfg: Config,
m3: Metrics,
m6: Metrics,
faults: List[FaultSummary],
pairing: List[Dict],
):
lines = []
lines.append("=" * 88)
lines.append("BALANCED-ELECTRIC-5 — FULL SIX-PHASE MACHINE / LOAD / FAULT TEST")
lines.append("=" * 88)
lines.append("")
lines.append("MODEL STATUS")
lines.append("This is a physically constrained screening model, not FEA or hardware validation.")
lines.append("The comparison preserves six independent phases rather than pair-summing them.")
lines.append("Voltage/current addition from the earlier A+F/B+E/C+D topology is reported separately.")
lines.append("")
lines.append("FAIRNESS CONSTRAINTS")
lines.append(f"Target mechanical power: {cfg.target_power_w:.3f} W")
lines.append(f"Target speed: {cfg.target_speed_rpm:.3f} rpm")
lines.append(f"Source phase RMS voltage: {cfg.v_phase_rms:.3f} V")
lines.append(f"Per-conductor current limit: {cfg.current_limit:.3f} A")
lines.append(f"Voltage limit: {cfg.voltage_limit:.3f} V")
lines.append(f"3-phase total copper R: {cfg.copper_resistance_total_3ph:.6f} ohm")
lines.append(f"6-phase total copper R: {cfg.copper_resistance_total_6ph:.6f} ohm")
lines.append(f"Thermal limit: {cfg.thermal_limit_c:.3f} C")
lines.append("")
lines.append("NOMINAL RESULTS")
lines.append("")
lines.append("Metric 3-phase 6-phase")
lines.append("-" * 70)
def line(label, a, b, suffix=""):
lines.append(
f"{label:<30} {a:>14.5f}{suffix:<8} {b:>14.5f}{suffix}"
)
line("Speed rpm", m3.speed_rpm, m6.speed_rpm)
line("Output W", m3.pout_w, m6.pout_w)
line("Input W", m3.pin_w, m6.pin_w)
line("Efficiency", 100*m3.efficiency, 100*m6.efficiency, "%")
line("Copper loss W", m3.pcu_w, m6.pcu_w)
line("Total RMS current A", m3.current_rms_total, m6.current_rms_total)
line("Max phase RMS current A", m3.current_rms_max, m6.current_rms_max)
line("Torque ripple", m3.torque_ripple_pct, m6.torque_ripple_pct, "%")
line("Power ripple", m3.power_ripple_pct, m6.power_ripple_pct, "%")
line("Max temperature C", m3.thermal_c, m6.thermal_c)
line("Thermal margin C", m3.thermal_margin_c, m6.thermal_margin_c)
line("Current utilization", m3.current_utilization_pct, m6.current_utilization_pct, "%")
line("Voltage utilization", m3.voltage_utilization_pct, m6.voltage_utilization_pct, "%")
line("Energy residual", m3.energy_residual_pct, m6.energy_residual_pct, "%")
lines.append("")
lines.append("FAULT SURVIVAL")
lines.append("-" * 88)
for f in faults:
lines.append(
f"{f.architecture:<10} {f.fault:<24} "
f"retained={f.retained_output_pct:8.3f}% "
f"speed={f.speed_rpm:8.2f} rpm "
f"T_ripple={f.torque_ripple_pct:9.3f}% "
f"Tmax={f.max_thermal_c:8.3f} C "
f"{f.protection_flag}"
)
lines.append("")
lines.append("EARLIER A+F / B+E / C+D PAIRING — PHASOR-ONLY DIAGNOSTIC")
lines.append("-" * 88)
for p in pairing:
lines.append(
f"{p['pair']:<8} amp={p['amplitude']:.9f} "
f"phase={p['phase_deg']:.6f} deg "
f"positive={p['positive_sequence']:.9f} "
f"negative={p['negative_sequence']:.9f} "
f"neg/pos={p['negative_positive_pct']:.3f}%"
)
lines.append("")
lines.append("INTERPRETATION RULES")
lines.append("-" * 88)
lines.append("1. No result here proves superiority of one machine technology.")
lines.append("2. Lower power ripple is not the same thing as higher efficiency.")
lines.append("3. A phase-fault result is meaningful only inside this model's limits.")
lines.append("4. The earlier pair-sum topology is not treated as a power-conserving machine connection.")
lines.append("5. Real claims require FEA, inverter-loss modeling, thermal validation and hardware tests.")
lines.append("6. The most useful comparison is equal output + equal copper + equal voltage/current/thermal limits.")
path.write_text("\n".join(lines), encoding="utf-8")
# =============================================================================
# PLOTS
# =============================================================================
def plot_nominal(
out: Path,
r3: SimResult,
r6: SimResult,
):
if plt is None:
return
sl3 = slice(int(len(r3.t) * 0.65), None)
sl6 = slice(int(len(r6.t) * 0.65), None)
fig, ax = plt.subplots(figsize=(12, 6))
ax.plot(r3.t[sl3], r3.torque[sl3], label="3-phase torque")
ax.plot(r6.t[sl6], r6.torque[sl6], label="6-phase torque")
ax.set_xlabel("Time (s)")
ax.set_ylabel("Electromagnetic torque (Nm)")
ax.set_title("Nominal Electromagnetic Torque")
ax.grid(True, alpha=0.25)
ax.legend()
fig.tight_layout()
fig.savefig(out / "torque_comparison.png", dpi=160)
plt.close(fig)
fig, ax = plt.subplots(figsize=(12, 6))
ax.plot(r3.t[sl3], r3.speed_rpm[sl3], label="3-phase speed")
ax.plot(r6.t[sl6], r6.speed_rpm[sl6], label="6-phase speed")
ax.set_xlabel("Time (s)")
ax.set_ylabel("Speed (rpm)")
ax.set_title("Nominal Speed")
ax.grid(True, alpha=0.25)
ax.legend()
fig.tight_layout()
fig.savefig(out / "nominal_speed.png", dpi=160)
plt.close(fig)
fig, ax = plt.subplots(figsize=(12, 6))
ax.plot(r3.t[sl3], r3.thermal_c[sl3], label="3-phase temperature")
ax.plot(r6.t[sl6], r6.thermal_c[sl6], label="6-phase temperature")
ax.set_xlabel("Time (s)")
ax.set_ylabel("Winding temperature (C)")
ax.set_title("Nominal Thermal Response")
ax.grid(True, alpha=0.25)
ax.legend()
fig.tight_layout()
fig.savefig(out / "thermal_comparison.png", dpi=160)
plt.close(fig)
fig, ax = plt.subplots(figsize=(12, 6))
ax.plot(r3.t[sl3], r3.electrical_power[sl3], label="3-phase input power")
ax.plot(r6.t[sl6], r6.electrical_power[sl6], label="6-phase input power")
ax.plot(r3.t[sl3], r3.mechanical_power[sl3], "--", label="3-phase mechanical")
ax.plot(r6.t[sl6], r6.mechanical_power[sl6], "--", label="6-phase mechanical")
ax.set_xlabel("Time (s)")
ax.set_ylabel("Power (W)")
ax.set_title("Electrical Input / Mechanical Output")
ax.grid(True, alpha=0.25)
ax.legend()
fig.tight_layout()
fig.savefig(out / "power_comparison.png", dpi=160)
plt.close(fig)
def plot_faults(
out: Path,
fault_rows: List[FaultSummary],
):
if plt is None or not fault_rows:
return
names = []
vals = []
groups = []
for r in fault_rows:
names.append(r.fault.replace(r.architecture + "_", ""))
vals.append(r.retained_output_pct)
groups.append(r.architecture)
x = np.arange(len(names))
fig, ax = plt.subplots(figsize=(13, 6))
for arch in sorted(set(groups)):
xx = [i for i, g in enumerate(groups) if g == arch]
yy = [vals[i] for i in xx]
ax.plot(xx, yy, "o-", label=arch)
ax.set_xticks(x)
ax.set_xticklabels(names, rotation=45, ha="right")
ax.set_ylabel("Retained output (%)")
ax.set_title("Single-Phase Open Fault Survival")
ax.grid(True, alpha=0.25)
ax.legend()
fig.tight_layout()
fig.savefig(out / "fault_survival.png", dpi=160)
plt.close(fig)
def plot_harmonics(
out: Path,
cfg: Config,
r3: SimResult,
r6: SimResult,
):
if plt is None:
return
fs = 1.0 / cfg.dt
fig, ax = plt.subplots(figsize=(12, 6))
for r in [r3, r6]:
sl = slice(int(len(r.t) * (1.0 - cfg.rms_tail_fraction)), None)
i = r.current[sl]
vals = []
for h in range(1, cfg.fft_harmonics_max + 1):
amps = []
for p in range(i.shape[1]):
a = harmonic_amplitudes(
i[:, p],
fs,
cfg.f_e,
cfg.fft_harmonics_max,
)
amps.append(a[h])
vals.append(np.sqrt(np.sum(np.square(amps))))
ax.stem(
np.arange(1, cfg.fft_harmonics_max + 1),
vals,
markerfmt=".",
basefmt=" ",
label=r.name,
)
ax.set_xlabel("Harmonic order")
ax.set_ylabel("Aggregate current harmonic amplitude")
ax.set_title("Current Harmonic Spectrum")
ax.grid(True, alpha=0.25)
ax.legend()
fig.tight_layout()
fig.savefig(out / "harmonic_spectrum.png", dpi=160)
plt.close(fig)
# =============================================================================
# MAIN
# =============================================================================
def parse_args() -> argparse.Namespace:
p = argparse.ArgumentParser(
description="Full six-phase vs conventional 3-phase machine/load/fault test."
)
p.add_argument("--seconds", type=float, default=0.60)
p.add_argument("--dt", type=float, default=2.0e-5)
p.add_argument("--fault-duration", type=float, default=0.20)
p.add_argument("--save-dir", type=str, default="balanced-electric5-results")
p.add_argument("--no-plots", action="store_true")
return p.parse_args()
def main():
args = parse_args()
cfg = Config(
seconds=args.seconds,
dt=args.dt,
fault_duration_s=args.fault_duration,
)
out = Path(args.save_dir)
out.mkdir(parents=True, exist_ok=True)
print("=" * 100)
print("BALANCED-ELECTRIC-5 — FULL SIX-PHASE MACHINE / LOAD / FAULT TEST")
print("=" * 100)
print()
print("This is a constrained engineering screening model.")
print("It preserves six independent phases; it does NOT pair-sum them in the")
print("main machine model. The earlier A+F/B+E/C+D connection is diagnostic only.")
print()
m3_machine = make_machine(cfg, 3)
m6_machine = make_machine(cfg, 6)
print("MACHINE CONFIGURATION")
print(f" 3-phase axes: {np.rad2deg(m3_machine.axes)}")
print(f" 6-phase axes: {np.rad2deg(m6_machine.axes)}")
print(f" total copper R, 3φ: {m3_machine.copper_resistance_sum:.6f} ohm")
print(f" total copper R, 6φ: {m6_machine.copper_resistance_sum:.6f} ohm")
print(f" target: {cfg.target_power_w:.1f} W @ {cfg.target_speed_rpm:.1f} rpm")
print()
print("RUNNING NOMINAL 3-PHASE...")
r3 = simulate(cfg, m3_machine)
m3 = calculate_metrics(cfg, m3_machine, r3)
print("RUNNING NOMINAL 6-PHASE...")
r6 = simulate(cfg, m6_machine)
m6 = calculate_metrics(cfg, m6_machine, r6)
print_metrics(m3)
print_metrics(m6)
print()
print("=" * 100)
print("DIRECT NOMINAL COMPARISON")
print("=" * 100)
def compare(label: str, a: float, b: float, unit: str = ""):
delta = b - a
rel = safe_div(delta, abs(a)) * 100.0
print(
f"{label:<28} 3φ={a:>12.5f}{unit:<5} "
f"6φ={b:>12.5f}{unit:<5} "
f"delta={delta:>12.5f}{unit:<5} ({rel:+.3f}%)"
)
compare("Output power", m3.pout_w, m6.pout_w, "W")
compare("Input power", m3.pin_w, m6.pin_w, "W")
compare("Efficiency", 100*m3.efficiency, 100*m6.efficiency, "%")
compare("Copper loss", m3.pcu_w, m6.pcu_w, "W")
compare("Total RMS current", m3.current_rms_total, m6.current_rms_total, "A")
compare("Max phase RMS current", m3.current_rms_max, m6.current_rms_max, "A")
compare("Torque ripple", m3.torque_ripple_pct, m6.torque_ripple_pct, "%")
compare("Power ripple", m3.power_ripple_pct, m6.power_ripple_pct, "%")
compare("Maximum temperature", m3.thermal_c, m6.thermal_c, "C")
compare("Current utilization", m3.current_utilization_pct, m6.current_utilization_pct, "%")
print()
print("=" * 100)
print("RUNNING SINGLE-PHASE OPEN FAULTS")
print("=" * 100)
fault_rows: List[FaultSummary] = []
fault_rows.extend(run_faults(cfg, m3_machine, m3))
fault_rows.extend(run_faults(cfg, m6_machine, m6))
for f in fault_rows:
print(
f"{f.architecture:>8} | {f.fault:<20} | "
f"retained={f.retained_output_pct:8.3f}% | "
f"speed={f.speed_rpm:8.2f} rpm | "
f"torque ripple={f.torque_ripple_pct:9.3f}% | "
f"Tmax={f.max_thermal_c:8.3f} C | "
f"{f.protection_flag}"
)
print()
print("=" * 100)
print("EARLIER A+F / B+E / C+D PAIRING DIAGNOSTIC")
print("=" * 100)
pairing = pairing_diagnostic()
for p in pairing:
print(
f"{p['pair']:<8} "
f"amp={p['amplitude']:.9f} "
f"phase={p['phase_deg']:.6f}° "
f"positive={p['positive_sequence']:.9f} "
f"negative={p['negative_sequence']:.9f} "
f"zero={p['zero_sequence']:.9f} "
f"neg/pos={p['negative_positive_pct']:.3f}%"
)
# CSVs
write_nominal_csv(out / "nominal.csv", [r3, r6], every=max(1, int(0.0002 / cfg.dt)))
metrics_rows = [
metrics_dict(m3),
metrics_dict(m6),
]
write_csv(
out / "nominal_metrics.csv",
metrics_rows,
list(metrics_rows[0].keys()),
)
fault_dict_rows = [vars(x) for x in fault_rows]
write_csv(
out / "faults.csv",
fault_dict_rows,
list(fault_dict_rows[0].keys()) if fault_dict_rows else None,
)
harm_rows = harmonic_metrics(cfg, r3) + harmonic_metrics(cfg, r6)
write_csv(
out / "harmonics.csv",
harm_rows,
["architecture", "phase", "harmonic", "amplitude"],
)
# Energy balance trace.
eb_rows = []
for r in [r3, r6]:
every = max(1, int(0.0002 / cfg.dt))
for k in range(0, len(r.t), every):
visc = cfg.viscous_b * r.omega[k] ** 2
residual = (
r.electrical_power[k]
- r.mechanical_power[k]
- r.copper_loss[k]
- visc
)
eb_rows.append([
r.name,
r.t[k],
r.electrical_power[k],
r.mechanical_power[k],
r.copper_loss[k],
visc,
residual,
])
write_csv(
out / "energy_balance.csv",
eb_rows,
[
"architecture",
"time_s",
"electrical_input_W",
"mechanical_output_W",
"copper_loss_W",
"viscous_loss_W",
"residual_W",
],
)
write_csv(
out / "topology_pairing_diagnostic.csv",
pairing,
[
"pair",
"amplitude",
"phase_deg",
"positive_sequence",
"negative_sequence",
"zero_sequence",
"negative_positive_pct",
],
)
write_summary(
out / "summary.txt",
cfg,
m3,
m6,
fault_rows,
pairing,
)
if not args.no_plots:
if plt is None:
print("\nmatplotlib unavailable: plots skipped.")
else:
plot_nominal(out, r3, r6)
plot_faults(out, fault_rows)
plot_harmonics(out, cfg, r3, r6)
print()
print("=" * 100)
print("FILES WRITTEN")
print("=" * 100)
for p in sorted(out.iterdir()):
print(f" {p}")
print()
print("=" * 100)
print("IMPORTANT INTERPRETATION")
print("=" * 100)
print("This test answers a much stronger question than the earlier waveform test:")
print("whether a six-phase machine can preserve useful output under matched")
print("copper/electrical/thermal constraints.")
print()
print("It does NOT establish commercial superiority. A real design still needs:")
print(" * electromagnetic FEA")
print(" * slot/winding-factor optimization")
print(" * iron-loss and saturation modeling")
print(" * inverter switching/conduction losses")
print(" * detailed thermal network")
print(" * insulation and fault-protection engineering")
print(" * acoustic/vibration testing")
print(" * hardware validation")
print()
print(f"Results directory: {out.resolve()}")
if __name__ == "__main__":
main()
#!/usr/bin/env python3
"""
BALANCED-ELECTRIC-6 — EXACT USER TOPOLOGY, POWER-CONSERVING REDUCED MODEL
===============================================================================
Compares conventional 3-phase against the exact proposed topology:
A=0°, B=30°, C=120°, D=150°, E=240°, F=270°
U=A+F, V=B+E, W=C+D
The six physical windings are retained. The pair connection is represented by
an incidence matrix C and reduced terminal variables:
i_branch = C i_pair
v_pair = C.T v_branch
This guarantees v_branch.T i_branch = v_pair.T i_pair.
Unlike the previous broken delivery, this file is ONLY the runnable program;
it never tries to write itself to /mnt/data. It creates output files next to
itself in balanced-electric6-output/.
This is an engineering screening model, not FEA or hardware certification.
It includes sinusoidal flux linkage, mutual inductance, copper loss, iron loss,
inverter loss, mechanical loss, thermal state, energy accounting, harmonics,
common-mode voltage, nominal operation, and open-pair/open-phase faults.
Requires: Python 3.10+, numpy. matplotlib is optional.
"""
from __future__ import annotations
import argparse, csv, math
from dataclasses import dataclass
from pathlib import Path
import numpy as np
try:
import matplotlib.pyplot as plt
except Exception:
plt = None
DEG=math.pi/180.0
TARGET_POWER=5000.0
TARGET_RPM=1800.0
OMEGA_TARGET=TARGET_RPM*2*math.pi/60
TOTAL_COPPER_R=0.90
CURRENT_LIMIT_RMS=18.0
TERMINAL_VOLTAGE_LIMIT_RMS=240.0
AMBIENT=25.0
THERMAL_LIMIT=120.0
THERMAL_RISE_PER_W=0.075
THERMAL_TAU=18.0
INERTIA=0.060
POLE_PAIRS=1
FLUX_PEAK=1.20
L3_SELF=0.012
L6_SELF=0.006
M6_NEAR=0.0012
M6_OPP=0.0005
IRON_BASE=38.0
IRON_SPEED=0.004
INV_FRAC=0.010
MECH_BASE=18.0
MECH_SPEED=0.0025
CONTROL_TAU=0.0008
DT_DEFAULT=4e-5
LOAD_TORQUE=TARGET_POWER/OMEGA_TARGET
AX3_DEG=np.array([0.,120.,240.])
AX6_DEG=np.array([0.,30.,120.,150.,240.,270.])
AX3=AX3_DEG*DEG
AX6=AX6_DEG*DEG
# [A B C D E F] = C [U V W]
C6=np.array([[1,0,0],[0,1,0],[0,0,1],[0,0,1],[0,1,0],[1,0,0]],float)
@dataclass
class System:
name:str
axes:np.ndarray
C:np.ndarray
R_branch:np.ndarray
L_branch:np.ndarray
R_term:np.ndarray
L_term:np.ndarray
flux_term_axes:np.ndarray
flux_term_coeff:np.ndarray
@dataclass
class Result:
system:System
fault:str
t:np.ndarray
speed:np.ndarray
theta:np.ndarray
i_branch:np.ndarray
v_branch:np.ndarray
i_term:np.ndarray
v_term:np.ndarray
torque:np.ndarray
pin:np.ndarray
pout:np.ndarray
pcu:np.ndarray
piron:np.ndarray
pinv:np.ndarray
pmechloss:np.ndarray
temp:np.ndarray
cmv:np.ndarray
summary:dict
def make_L(n,self_l,near,opp=0):
M=np.zeros((n,n))
for i in range(n):
for j in range(n):
if i==j: M[i,j]=self_l
else:
d=min((i-j)%n,(j-i)%n)
if d==1: M[i,j]=near
elif n==6 and d==3: M[i,j]=opp
else: M[i,j]=near*0.35
e=np.linalg.eigvalsh(M)
if e.min()<=1e-9: M+=np.eye(n)*(1e-9-e.min()+1e-10)
return M
def systems():
r3=TOTAL_COPPER_R/3
r6=TOTAL_COPPER_R/6
L3=make_L(3,L3_SELF,0.0015)
L6=make_L(6,L6_SELF,M6_NEAR,M6_OPP)
s3=System('3-phase',AX3,np.eye(3),np.full(3,r3),L3,np.eye(3)*r3,L3,AX3,np.ones(3))
# The pair terminal flux is C.T lambda_branch. It is generally not a
# single sinusoid; keeping the six components is important.
s6=System('USER-6',AX6,C6,np.full(6,r6),L6,C6.T@np.diag(np.full(6,r6))@C6,C6.T@L6@C6,AX6,np.ones(6))
return s3,s6
def branch_flux(s,theta):
return FLUX_PEAK*np.cos(POLE_PAIRS*theta-s.axes)
def branch_dflux(s,theta):
return -FLUX_PEAK*POLE_PAIRS*np.sin(POLE_PAIRS*theta-s.axes)
def branch_emf(s,theta,omega):
return branch_dflux(s,theta)*omega
def torque(s,theta,ib):
return float(np.dot(ib,branch_dflux(s,theta)))
def pair_current_reference(s,theta,torque_cmd):
g=s.C.T@branch_dflux(s,theta)
gg=float(g@g)
if gg<1e-14 or torque_cmd<=0: return np.zeros(3)
ip=(torque_cmd/gg)*g
# branch RMS current limit, because physical conductors carry pair current.
ib=s.C@ip
mx=float(np.max(np.abs(ib)))
# Instantaneous peak limit corresponding to RMS conductor limit.
lim=CURRENT_LIMIT_RMS*math.sqrt(2)
if mx>lim: ip*=lim/mx
return ip
def fft_thd(x,dt,f0):
x=np.asarray(x,float); n=len(x)
if n<128 or f0<=0:return 0.0
y=x-np.mean(x); w=np.hanning(n); sp=np.fft.rfft(y*w); fr=np.fft.rfftfreq(n,dt)
k=int(np.argmin(np.abs(fr-f0)))
if k<1 or abs(sp[k])<1e-14:return 0.0
bins=np.arange(2*k,len(sp),k); bins=bins[bins<len(sp)]
return float(np.sqrt(np.sum(np.abs(sp[bins])**2))/abs(sp[k]))
def simulate(s,seconds,dt,fault_pair=None,fault_start=None,fault_duration=0.0):
n=int(round(seconds/dt))+1; t=np.arange(n)*dt
speed=np.zeros(n); theta=np.zeros(n)
ib=np.zeros((n,6 if s.name=='USER-6' else 3)); vb=np.zeros_like(ib)
it=np.zeros((n,3)); vt=np.zeros((n,3)); te=np.zeros(n)
pin=np.zeros(n); pout=np.zeros(n); pcu=np.zeros(n); piron=np.zeros(n); pinv=np.zeros(n); pmech_loss_arr=np.zeros(n)
temp=np.full(n,AMBIENT); cm=np.zeros(n)
ip=np.zeros(3); omega=OMEGA_TARGET*0.985; th=0.0; T=AMBIENT
if fault_start is None:fault_start=seconds*0.45
fend=fault_start+fault_duration
target=OMEGA_TARGET
for k,tk in enumerate(t[:-1]):
active=fault_pair is not None and fault_start<=tk<fend
tl=LOAD_TORQUE
# PI-like correction, with explicit load-torque feed-forward.
corr=0.10*(target-omega)
torque_cmd=float(np.clip(tl+corr,0,2*tl))
ip_ref=pair_current_reference(s,th,torque_cmd)
if active: ip_ref[fault_pair]=0.0
eb=branch_emf(s,th,omega)
et=s.C.T@eb
# Reduced KVL: vt = Rterm*i + Lterm*di/dt + et.
raw=(ip_ref-ip)/CONTROL_TAU
base=s.R_term@ip+et
coeff=np.diag(s.L_term)+s.R_term*dt
# Diagonal approximation is used only for voltage limiter; actual
# coupled L is retained in the voltage calculation below.
# Conservative scalar scaling keeps the full vector within terminal limit.
di=raw.copy()
for _ in range(4):
trial=ip+di*dt
dtrial=(trial-ip)/dt
vtrial=s.R_term@trial+s.L_term@dtrial+et
mx=float(np.max(np.abs(vtrial)))
lim=TERMINAL_VOLTAGE_LIMIT_RMS*math.sqrt(2)
if mx>lim: di*=lim/mx
else: break
inew=ip+di*dt
# Terminal current is also conductor current; instantaneous limit.
limi=CURRENT_LIMIT_RMS*math.sqrt(2)
ibtrial=s.C@inew; imax=float(np.max(np.abs(ibtrial)))
if imax>limi:inew*=limi/imax
if active:inew[fault_pair]=0.0
di=(inew-ip)/dt
vt_now=s.R_term@inew+s.L_term@di+et
if active:vt_now[fault_pair]=0.0
ib_now=s.C@inew; vb_now=np.linalg.lstsq(s.C.T,vt_now,rcond=None)[0] if s.name=='USER-6' else vt_now.copy()
# For the series connection, branch voltages are not arbitrary: use the
# physical branch KVL. Pair power equals branch power exactly by C.
vb_phys=s.R_branch*ib_now+s.L_branch@ (s.C@di)+eb
if s.name=='USER-6':
# vb_phys sums to vt_now; numerical KVL defines physical branch volts.
vb_now=vb_phys
te_now=torque(s,th,ib_now)
p_e=float(np.dot(vt_now,inew))
p_c=float(np.sum(s.R_branch*ib_now**2))
pir=IRON_BASE+IRON_SPEED*omega**2
pml=MECH_BASE+MECH_SPEED*omega**2
piv=INV_FRAC*abs(p_e)
alpha=(te_now-tl-pml/max(omega,10))/INERTIA
omega=max(0,omega+alpha*dt); th+=omega*dt
heat=p_c+pir+piv
Ttarget=AMBIENT+THERMAL_RISE_PER_W*heat
T+=(Ttarget-T)*dt/THERMAL_TAU
speed[k]=omega*60/(2*math.pi);theta[k]=th;ib[k+1]=ib_now;vb[k]=vb_now
it[k+1]=inew;vt[k]=vt_now;te[k]=te_now;pin[k]=p_e;pout[k]=max(0,te_now*omega)
pcu[k]=p_c;piron[k]=pir;pinv[k]=piv;pmech_loss_arr[k]=pml;temp[k]=T;cm[k]=np.mean(vb_now);ip=inew
speed[-1]=omega*60/(2*math.pi);theta[-1]=th;te[-1]=te[-2];pin[-1]=pin[-2];pout[-1]=pout[-2];pcu[-1]=pcu[-2];piron[-1]=piron[-2];pinv[-1]=pinv[-2];pmech_loss_arr[-1]=pmech_loss_arr[-2];temp[-1]=temp[-2];vb[-1]=vb[-2];vt[-1]=vt[-2];ib[-1]=ib[-2];it[-1]=it[-2];cm[-1]=cm[-2]
cut=max(1,n//2); sl=slice(cut,None)
ir=np.sqrt(np.mean(ib[sl]**2,axis=0)); vr=np.sqrt(np.mean(vb[sl]**2,axis=0));
I=rms_vec(np.sqrt(np.sum(ib[sl]**2,axis=1))); V=rms_vec(np.sqrt(np.sum(vb[sl]**2,axis=1)))
Pin=float(np.mean(pin[sl])); Pout=float(np.mean(pout[sl])); tq=float(np.mean(te[sl]))
Ein=float(np.trapezoid(pin[sl],t[sl])); Eout=float(np.trapezoid(pout[sl],t[sl])); Ec=float(np.trapezoid(pcu[sl],t[sl])); Ei=float(np.trapezoid(piron[sl],t[sl])); Ev=float(np.trapezoid(pinv[sl],t[sl])); Em=float(np.trapezoid(pmech_loss_arr[sl],t[sl]))
w0=speed[cut]*2*math.pi/60;w1=speed[-1]*2*math.pi/60;dEk=.5*INERTIA*(w1*w1-w0*w0)
dEm=.5*float(it[-1]@s.L_term@it[-1])-0.5*float(it[cut]@s.L_term@it[cut])
# Electromagnetic energy closure: terminal electrical input must equal
# copper dissipation + electromagnetic mechanical conversion + change in
# magnetic stored energy. Iron/inverter/mechanical losses are modeled as
# separate system-level losses and are therefore not double-counted here.
Eem=float(np.trapezoid(pout[sl],t[sl]))
accounted=Eout+Ec+dEm
residual=100*(Ein-accounted)/max(abs(Ein),1e-12)
f0=max(speed[-1],1)/60
cmr=rms_vec(cm[sl])
summ={'speed_rpm':float(speed[-1]),'total_rms_current_A':I,'max_phase_rms_current_A':float(np.max(ir)),
'total_rms_voltage_V':V,'max_phase_rms_voltage_V':float(np.max(vr)),'input_power_W':Pin,'mechanical_output_W':Pout,
'copper_loss_W':float(np.mean(pcu[sl])),'iron_loss_W':float(np.mean(piron[sl])),'inverter_loss_W':float(np.mean(pinv[sl])),
'mechanical_loss_W':float(np.mean(pmech_loss_arr[sl])),'efficiency_pct':100*Pout/max(Pin,1e-12),
'model_pf':Pin/max(I*V,1e-12),'mean_torque_Nm':tq,'torque_ripple_pct':100*float(np.std(te[sl]))/max(abs(tq),1e-12),
'max_temperature_C':float(np.max(temp[sl])),'thermal_margin_C':THERMAL_LIMIT-float(np.max(temp[sl])),
'current_utilization_pct':100*float(np.max(ir))/CURRENT_LIMIT_RMS,'voltage_utilization_pct':100*float(np.max(vr))/TERMINAL_VOLTAGE_LIMIT_RMS,
'input_power_ripple_pct':100*float(np.std(pin[sl]))/max(abs(Pin),1e-12),
'current_thd_pct':100*float(np.mean([fft_thd(ib[sl,j],dt,f0) for j in range(ib.shape[1])])),
'voltage_thd_pct':100*float(np.mean([fft_thd(vb[sl,j],dt,f0) for j in range(vb.shape[1])])),
'common_mode_rms_V':cmr,'energy_residual_pct':residual,'electromagnetic_conversion_J':Eem,'magnetic_energy_delta_J':dEm}
return Result(s,'none' if fault_pair is None else f'{s.name}_pair_{fault_pair+1}_OPEN',t,speed,theta,ib,vb,it,vt,te,pin,pout,pcu,piron,pinv,pmech_loss_arr,temp,cm,summ)
def rms_vec(x):return float(np.sqrt(np.mean(np.square(x))))
def status(s):
if s['thermal_margin_C']<0:return 'THERMAL_LIMIT'
if s['max_phase_rms_current_A']>CURRENT_LIMIT_RMS*1.001:return 'CURRENT_LIMIT'
if s['max_phase_rms_voltage_V']>TERMINAL_VOLTAGE_LIMIT_RMS*1.001:return 'VOLTAGE_LIMIT'
return 'WITHIN_MODEL_LIMITS'
def print_summary(r):
s=r.summary;print(f'\n[{r.system.name}] {r.fault}')
rows=[('speed','speed_rpm','rpm'),('total RMS current','total_rms_current_A','A'),('max phase RMS current','max_phase_rms_current_A','A'),('total RMS voltage','total_rms_voltage_V','V'),('max phase RMS voltage','max_phase_rms_voltage_V','V'),('input power','input_power_W','W'),('mechanical output','mechanical_output_W','W'),('copper loss','copper_loss_W','W'),('iron loss','iron_loss_W','W'),('inverter loss','inverter_loss_W','W'),('mechanical loss','mechanical_loss_W','W'),('efficiency','efficiency_pct','%'),('model PF','model_pf',''),('mean torque','mean_torque_Nm','Nm'),('torque ripple','torque_ripple_pct','%'),('max temperature','max_temperature_C','C'),('thermal margin','thermal_margin_C','C'),('current utilization','current_utilization_pct','%'),('voltage utilization','voltage_utilization_pct','%'),('input-power ripple','input_power_ripple_pct','%'),('current THD','current_thd_pct','%'),('voltage THD','voltage_thd_pct','%'),('common-mode RMS','common_mode_rms_V','V'),('energy residual','energy_residual_pct','%')]
for a,b,c in rows:print(f' {a:24s} {s[b]:12.5f} {c}')
def compare(a,b):
print('\n'+'='*100+'\nDIRECT NOMINAL COMPARISON\n'+'='*100)
keys=[('Output power','mechanical_output_W','W'),('Input power','input_power_W','W'),('Efficiency','efficiency_pct','%'),('Copper loss','copper_loss_W','W'),('Iron loss','iron_loss_W','W'),('Inverter loss','inverter_loss_W','W'),('Total RMS current','total_rms_current_A','A'),('Max phase RMS current','max_phase_rms_current_A','A'),('Torque ripple','torque_ripple_pct','%'),('Power ripple','input_power_ripple_pct','%'),('Current THD','current_thd_pct','%'),('Voltage THD','voltage_thd_pct','%'),('Max temperature','max_temperature_C','C'),('Current utilization','current_utilization_pct','%'),('Voltage utilization','voltage_utilization_pct','%'),('Common-mode RMS','common_mode_rms_V','V'),('Energy residual','energy_residual_pct','%')]
for lab,k,u in keys:
x=a.summary[k];y=b.summary[k];d=y-x;p=100*d/max(abs(x),1e-12);print(f'{lab:26s} 3φ={x:12.5f}{u:2s} USER-6={y:12.5f}{u:2s} delta={d:12.5f}{u:2s} ({p:+8.3f}%)')
def diagnostic():
z=np.exp(1j*AX6);alpha=np.exp(1j*120*DEG);names=['A','B','C','D','E','F'];pairs=[(0,5),(1,4),(2,3)];out=[]
for name,(p,q) in zip(['A+F','B+E','C+D'],pairs):
x=np.zeros(3,dtype=complex);x[['A+F','B+E','C+D'].index(name)]=z[p]+z[q]
pos=abs((x[0]+alpha*x[1]+alpha**2*x[2])/3);neg=abs((x[0]+alpha**2*x[1]+alpha*x[2])/3);zero=abs(np.sum(x)/3)
out.append((name,abs(z[p]+z[q]),math.degrees(np.angle(z[p]+z[q]))%360,pos,neg,zero,100*neg/max(pos,1e-15)))
return out
def save_csv(r,out):
p=out/f'{r.system.name}_{r.fault}.csv'
with p.open('w',newline='',encoding='utf-8') as f:
w=csv.writer(f);w.writerow(['t_s','speed_rpm','theta_rad','torque_Nm','p_in_W','p_out_W','p_cu_W','p_iron_W','p_inv_W','p_mech_loss_W','temperature_C','common_mode_V'])
for k in range(len(r.t)):w.writerow([r.t[k],r.speed[k],r.theta[k],r.torque[k],r.pin[k],r.pout[k],r.pcu[k],r.piron[k],r.pinv[k],r.pmechloss[k],r.temp[k],r.cmv[k]])
def make_plots(a,b,out):
if plt is None:return
plots=[('speed.png','Speed (rpm)',[(a.speed,'3-phase'),(b.speed,'USER-6')]),('torque.png','Torque (N·m)',[(a.torque,'3-phase'),(b.torque,'USER-6')]),('power.png','Power (W)',[(a.pin,'3-phase input'),(b.pin,'USER-6 input'),(a.pout,'3-phase shaft'),(b.pout,'USER-6 shaft')]),('temperature.png','Temperature (°C)',[(a.temp,'3-phase'),(b.temp,'USER-6')])]
for fn,yl,series in plots:
fig,ax=plt.subplots(figsize=(11,6))
for y,l in series:ax.plot(a.t if len(y)==len(a.t) else b.t,y,label=l)
ax.set_xlabel('Time (s)');ax.set_ylabel(yl);ax.set_title('Balanced-Electric-6');ax.grid(True,alpha=.25);ax.legend();fig.tight_layout();fig.savefig(out/fn,dpi=160);plt.close(fig)
def main():
ap=argparse.ArgumentParser(description='Exact Balanced-Electric-6 topology screening model')
ap.add_argument('--seconds',type=float,default=.60);ap.add_argument('--dt',type=float,default=DT_DEFAULT);ap.add_argument('--fault-duration',type=float,default=.20);ap.add_argument('--fault-start',type=float,default=None);ap.add_argument('--no-plots',action='store_true')
a=ap.parse_args()
if a.seconds<=0 or a.dt<=0 or a.fault_duration<0:raise SystemExit('seconds/dt must be positive; fault-duration must be nonnegative')
out=Path(__file__).resolve().parent/'balanced-electric6-output';out.mkdir(parents=True,exist_ok=True)
s3,s6=systems()
print('='*100);print('BALANCED-ELECTRIC-6 — EXACT USER TOPOLOGY / POWER-CONSERVING MODEL');print('='*100)
print('\nEXACT USER TOPOLOGY');print(' A=0°, B=30°, C=120°, D=150°, E=240°, F=270°');print(' U=A+F, V=B+E, W=C+D')
print('\nPOWER-CONSERVING CONNECTION');print(' i_branch = C i_pair');print(' v_pair = C^T v_branch');print(' v_branch^T i_branch = v_pair^T i_pair')
print(f'\nFAIRNESS: total copper R = {TOTAL_COPPER_R:.6f} ohm in both systems');print(f' terminal voltage limit = {TERMINAL_VOLTAGE_LIMIT_RMS:.1f} V RMS');print(f' conductor current limit = {CURRENT_LIMIT_RMS:.1f} A RMS');print(f' target = {TARGET_POWER:.1f} W @ {TARGET_RPM:.1f} rpm')
print('\nRUNNING NOMINAL 3-PHASE...');n3=simulate(s3,a.seconds,a.dt)
print('RUNNING NOMINAL EXACT USER-6...');n6=simulate(s6,a.seconds,a.dt)
print_summary(n3);print_summary(n6);compare(n3,n6)
print('\n'+'='*100+'\nRUNNING SINGLE-PAIR / SINGLE-PHYSICAL-PHASE OPEN FAULTS\n'+'='*100)
faults=[]
for s,nom in [(s3,n3),(s6,n6)]:
for j in range(3):
r=simulate(s,a.seconds,a.dt,j,a.fault_start,a.fault_duration);faults.append(r)
retained=100*r.summary['mechanical_output_W']/max(nom.summary['mechanical_output_W'],1e-12)
if s.name=='USER-6':desc=f'pair {j+1} (two physical series branches)'
else:desc=f'phase {j+1}'
print(f'{s.name:8s} | {desc:31s} | retained={retained:8.3f}% | speed={r.summary["speed_rpm"]:8.2f} rpm | torque ripple={r.summary["torque_ripple_pct"]:8.3f}% | Tmax={r.summary["max_temperature_C"]:8.3f} C | {status(r.summary)}')
print('\n'+'='*100+'\nEARLIER A+F / B+E / C+D PHASOR DIAGNOSTIC\n'+'='*100)
for x in diagnostic():print(f'{x[0]:8s} amp={x[1]:.9f} phase={x[2]:.6f}° positive={x[3]:.9f} negative={x[4]:.9f} zero={x[5]:.9f} neg/pos={x[6]:.3f}%')
save_csv(n3,out);save_csv(n6,out)
for r in faults:save_csv(r,out)
report=out/'balanced-electric6-report.txt'
with report.open('w',encoding='utf-8') as f:
f.write('BALANCED-ELECTRIC-6 SCREENING REPORT\n'+'='*80+'\n')
f.write('Exact topology: A=0°, B=30°, C=120°, D=150°, E=240°, F=270°; U=A+F, V=B+E, W=C+D\n\n')
for name,r in [('3-phase',n3),('USER-6',n6)]:
f.write(name+'\n');[f.write(f' {k}: {v:.9f}\n') for k,v in r.summary.items()]
f.write('\nFAULTS\n')
for r in faults:f.write(f'{r.system.name:8s} {r.fault:24s} output={r.summary["mechanical_output_W"]:.3f}W speed={r.summary["speed_rpm"]:.2f}rpm ripple={r.summary["torque_ripple_pct"]:.3f}% {status(r.summary)}\n')
f.write('\nLIMITATIONS: lumped sinusoidal windings; no FEA, saturation, slotting, detailed winding factors, switching-device physics, insulation, arc/breaker physics, or certification.\n')
if not a.no_plots and plt is not None:make_plots(n3,n6,out);print(f'\nPlots written to: {out}')
print(f'Report written to: {report}');print(f'CSV/output directory: {out}')
print('\nThis is the standalone repaired program. It does not generate or overwrite itself.')
if __name__=='__main__':main()
balanced-6x6.zip (47.4 KB)
