Battery Management System per Storage Grid-Scale

# Architettura software BMS - Stack completo

# Layer 1: Firmware (C/C++ su MCU)
# - Acquisizione dati ADC ad alta frequenza (1-1000 Hz)
# - Algoritmi real-time: SoC Coulomb counting, fault detection
# - Controllo attuatori: contattori, balancing, cooling

# Layer 2: Edge BMS Controller (Python/C++ su Linux SBC)
# - Algoritmi avanzati: EKF, SoH prediction, thermal model
# - Comunicazione con firmware via CAN/SPI
# - Buffer dati e aggregazione

# Layer 3: Rack/System BMS (Python su server industriale)
# - Coordinamento multi-rack
# - Interfaccia con PCS (Power Conversion System)
# - Comunicazione con EMS via Modbus TCP / IEC 61850

# Layer 4: Cloud/Edge Analytics
# - Fleet analytics su più siti
# - Training modelli ML per SoH prediction
# - Digital twin del sistema batterie

class BMSArchitecture:
    """
    Rappresentazione dell'architettura software BMS
    con responsabilità per ogni layer
    """

    LAYERS = {
        'firmware': {
            'language': 'C/C++',
            'os': 'FreeRTOS / Bare Metal',
            'cycle_time_ms': 1,          # 1 ms per fault detection
            'functions': [
                'cell_voltage_sampling',  # 1 kHz
                'current_integration',    # Coulomb counting
                'fault_detection',        # Hardware comparators
                'contactor_control',      # Fail-safe logic
                'passive_balancing'       # Resistive discharge
            ]
        },
        'bms_controller': {
            'language': 'Python / C++',
            'os': 'Linux (RT kernel)',
            'cycle_time_ms': 100,        # 100 ms per EKF update
            'functions': [
                'ekf_soc_estimation',
                'soh_prediction',
                'thermal_model',
                'active_balancing_control',
                'can_communication'
            ]
        },
        'system_bms': {
            'language': 'Python',
            'os': 'Linux',
            'cycle_time_ms': 1000,       # 1 s per grid commands
            'functions': [
                'multi_rack_coordination',
                'pcs_interface',          # Modbus TCP / IEC 61850
                'ems_interface',
                'scada_interface',
                'data_logging'
            ]
        }
    }

import numpy as np
from dataclasses import dataclass, field
from typing import Optional

@dataclass
class CoulombCounterSoC:
    """
    Stima SoC con Coulomb Counting.
    Semplice ma soggetto a deriva per errori di misura corrente.
    """
    capacity_ah: float           # capacità nominale [Ah]
    coulombic_efficiency: float  # Efficienza coulombica (tipica: 0.98-0.99 LFP, 0.995 NMC)
    soc: float = 1.0             # SoC iniziale [0-1]

    # Accumulo errori
    _accumulated_ah: float = field(default=0.0, repr=False)

    def update(self, current_a: float, dt_s: float) -> float:
        """
        Aggiorna SoC con misura di corrente.

        Args:
            current_a: Corrente [A]. Positivo = scarica, negativo = carica
            dt_s: Intervallo di campionamento [s]

        Returns:
            SoC aggiornato [0-1]
        """
        # Delta charge in Ah
        delta_ah = current_a * dt_s / 3600.0

        # Applica efficienza coulombica durante la carica
        if current_a < 0:  # Carica
            effective_delta = delta_ah * self.coulombic_efficiency
        else:              # Scarica
            effective_delta = delta_ah

        self._accumulated_ah += effective_delta

        # Aggiorna SoC
        new_soc = self.soc - (effective_delta / self.capacity_ah)

        # Clamp al range fisico [0, 1]
        self.soc = float(np.clip(new_soc, 0.0, 1.0))
        return self.soc

    def calibrate_from_ocv(self, ocv_v: float, ocv_soc_curve: dict) -> float:
        """
        Calibrazione SoC dalla curva OCV.
        Eseguita a riposo (corrente ~0 A per almeno 30 min).
        """
        # Interpolazione sulla curva OCV-SoC
        voltages = sorted(ocv_soc_curve.keys())
        socs = [ocv_soc_curve[v] for v in voltages]
        self.soc = float(np.interp(ocv_v, voltages, socs))
        self._accumulated_ah = 0.0
        return self.soc

import numpy as np
from scipy.interpolate import interp1d

class BatteryEKF:
    """
    Extended Kalman Filter per stima SoC con modello Thevenin 1° ordine.

    Stato: x = [SoC, V_RC]
    - SoC: State of Charge [0-1]
    - V_RC: Tensione sul circuito RC [V] (dinamica di polarizzazione)

    Riferimento: Plett, G.L. (2004) - "Extended Kalman Filtering for
    Battery Management Systems" - Journal of Power Sources
    """

    def __init__(self,
                 capacity_ah: float,
                 R0: float,
                 R1: float,
                 C1: float,
                 ocv_soc_table: tuple,
                 Q_noise: np.ndarray = None,
                 R_noise: float = None):
        """
        Args:
            capacity_ah: capacità nominale [Ah]
            R0: Resistenza ohmica [ohm]
            R1: Resistenza RC [ohm]
            C1: capacità RC [F]
            ocv_soc_table: (soc_array, ocv_array) per interpolazione
            Q_noise: Matrice covarianza rumore processo (2x2)
            R_noise: Varianza rumore misura tensione [V^2]
        """
        self.Q_batt = capacity_ah * 3600  # capacità in Coulomb
        self.R0 = R0
        self.R1 = R1
        self.C1 = C1

        # Curva OCV-SoC per interpolazione
        soc_pts, ocv_pts = ocv_soc_table
        self._ocv_func = interp1d(soc_pts, ocv_pts,
                                   kind='cubic',
                                   fill_value='extrapolate')
        # Derivata OCV rispetto a SoC (per linearizzazione)
        self._docv_dsoc = np.gradient(ocv_pts, soc_pts)
        self._docv_func = interp1d(soc_pts, self._docv_dsoc,
                                    kind='linear',
                                    fill_value='extrapolate')

        # Stato iniziale: x = [SoC, V_RC]
        self.x = np.array([1.0, 0.0])

        # Covarianza iniziale (incertezza elevata su SoC)
        self.P = np.diag([0.01, 0.001])  # [SoC^2, V^2]

        # Rumori
        self.Q = Q_noise if Q_noise is not None else np.diag([1e-6, 1e-8])
        self.R = R_noise if R_noise is not None else 1e-4  # 10 mV RMS

    def predict(self, current_a: float, dt_s: float) -> np.ndarray:
        """
        Step di predizione EKF.

        Modello dinamico (discretizzato con Eulero):
          SoC(k+1)  = SoC(k) - I*dt / Q_batt
          V_RC(k+1) = V_RC(k) * exp(-dt/(R1*C1)) + I * R1 * (1 - exp(-dt/(R1*C1)))

        Nota: corrente positiva = scarica (convenzione BMS)
        """
        soc, v_rc = self.x

        # Costante di tempo RC
        tau = self.R1 * self.C1
        exp_tau = np.exp(-dt_s / tau)

        # Predizione stato
        soc_pred = soc - (current_a * dt_s) / self.Q_batt
        v_rc_pred = v_rc * exp_tau + current_a * self.R1 * (1 - exp_tau)

        self.x = np.array([
            np.clip(soc_pred, 0.0, 1.0),
            v_rc_pred
        ])

        # Jacobiana del modello (matrice di transizione linearizzata)
        # F = d(f)/d(x)
        F = np.array([
            [1.0,  0.0],
            [0.0,  exp_tau]
        ])

        # Propagazione covarianza
        self.P = F @ self.P @ F.T + self.Q
        return self.x

    def update(self, v_terminal_measured: float, current_a: float) -> tuple:
        """
        Step di aggiornamento EKF con misura tensione terminale.

        Modello di osservazione:
          V_terminal = OCV(SoC) - I*R0 - V_RC

        Returns:
            (soc_estimate, covariance, innovation)
        """
        soc, v_rc = self.x

        # Tensione predetta dal modello
        v_oc = float(self._ocv_func(soc))
        v_predicted = v_oc - current_a * self.R0 - v_rc

        # Innovation (residuo)
        innovation = v_terminal_measured - v_predicted

        # Jacobiana dell'osservazione: H = d(h)/d(x)
        d_ocv_d_soc = float(self._docv_func(soc))
        H = np.array([[d_ocv_d_soc, -1.0]])

        # Covarianza dell'innovation
        S = H @ self.P @ H.T + self.R

        # Guadagno di Kalman
        K = self.P @ H.T / S

        # Aggiornamento stato e covarianza
        self.x = self.x + K.flatten() * innovation
        self.x[0] = np.clip(self.x[0], 0.0, 1.0)  # SoC in [0,1]

        I_matrix = np.eye(2)
        self.P = (I_matrix - np.outer(K.flatten(), H)) @ self.P

        return self.x[0], self.P[0, 0], innovation

    @property
    def soc(self) -> float:
        return float(self.x[0])

    @property
    def soc_uncertainty_1sigma(self) -> float:
        """Incertezza SoC a 1 sigma (68% confidence interval)"""
        return float(np.sqrt(self.P[0, 0]))


# ---- Esempio di utilizzo ----

def demo_ekf():
    # Curva OCV-SoC per LFP (LiFePO4) - valori tipici
    soc_pts = np.array([0.0, 0.1, 0.2, 0.3, 0.4, 0.5,
                        0.6, 0.7, 0.8, 0.9, 1.0])
    ocv_pts = np.array([3.0, 3.15, 3.22, 3.28, 3.30, 3.32,
                        3.33, 3.34, 3.35, 3.40, 3.65])  # Volt

    # Parametri batteria LFP grid-scale (es. CATL 280 Ah)
    bms_ekf = BatteryEKF(
        capacity_ah=280.0,
        R0=0.0002,     # 0.2 mohm - tipico LFP prismatico
        R1=0.0005,     # 0.5 mohm
        C1=5000.0,     # 5000 F = tau ~ 2.5 s
        ocv_soc_table=(soc_pts, ocv_pts),
        Q_noise=np.diag([1e-7, 1e-9]),
        R_noise=1e-5   # ~3.2 mV RMS tensione
    )

    # Simulazione: scarica a 0.5C per 100 step da 1 s
    dt = 1.0
    I_discharge = 140.0  # Ampere (0.5C su 280 Ah)

    results = []
    for step in range(100):
        # Predict
        bms_ekf.predict(I_discharge, dt)

        # Simula misura tensione con rumore
        true_ocv = float(np.interp(bms_ekf.soc, soc_pts, ocv_pts))
        v_meas = true_ocv - I_discharge * 0.0002 - bms_ekf.x[1]
        v_meas += np.random.normal(0, 0.003)  # 3 mV rumore

        # Update
        soc_est, variance, innov = bms_ekf.update(v_meas, I_discharge)

        results.append({
            'step': step,
            'soc': soc_est,
            'uncertainty': bms_ekf.soc_uncertainty_1sigma,
            'innovation_mv': innov * 1000
        })

    return results

if __name__ == '__main__':
    results = demo_ekf()
    print(f"SoC finale: {results[-1]['soc']:.3f}")
    print(f"Incertezza: ±{results[-1]['uncertainty']*100:.2f}% SoC")

import numpy as np
from dataclasses import dataclass

@dataclass
class BatteryDegradationModel:
    """
    Modello di degradazione batteria combinato calendar + cycle aging.

    Basato su: Wang et al. (2011) "Cycle-life model for graphite-LiFePO4 cells"
    Journal of Power Sources, adattato per applicazioni grid-scale.
    """

    # Parametri chimia LFP (calibrabili per NMC)
    # Calendar aging: Q_cal = B * exp(-Ea_cal / (R*T)) * sqrt(t)
    B_calendar: float = 14876      # Pre-exponential factor
    Ea_calendar: float = 24500.0   # Energia attivazione [J/mol] (LFP)

    # Cycle aging: Q_cyc = A * exp(Ea_cyc / (R*T)) * exp(b_dod * DoD) * N
    A_cycle: float = 7.543e5       # Pre-exponential factor
    Ea_cycle: float = -31700.0     # Energia attivazione [J/mol]
    b_dod: float = -0.836          # Coefficiente DoD

    R_gas: float = 8.314           # Costante gas [J/(mol*K)]

    def calendar_loss_fraction(self,
                                temp_k: float,
                                time_days: float,
                                avg_soc: float = 0.5) -> float:
        """
        Calcola la perdita di capacità per calendar aging.

        Args:
            temp_k: Temperatura media di stoccaggio [K]
            time_days: Tempo di calendario [giorni]
            avg_soc: SoC medio durante stoccaggio

        Returns:
            Frazione di capacità persa [0-1]
        """
        # Fattore temperatura (Arrhenius)
        k_cal = self.B_calendar * np.exp(-self.Ea_calendar / (self.R_gas * temp_k))

        # Fattore SoC (stress factor - più alto SoC = più degrado)
        soc_factor = 1 + 1.5 * (avg_soc - 0.5) ** 2

        # Legge di potenza radice quadrata (diffusione SEI)
        time_hours = time_days * 24
        loss = k_cal * soc_factor * np.sqrt(time_hours) / 100.0

        return float(np.clip(loss, 0.0, 1.0))

    def cycle_loss_fraction(self,
                             temp_k: float,
                             dod: float,
                             n_cycles: int,
                             avg_crate: float = 0.5) -> float:
        """
        Calcola la perdita di capacità per cycle aging.

        Args:
            temp_k: Temperatura operativa media [K]
            dod: Depth of Discharge [0-1]
            n_cycles: Numero di cicli equivalenti a piena profondità
            avg_crate: C-rate medio (1C = descarga in 1 ora)

        Returns:
            Frazione di capacità persa [0-1]
        """
        # Fattore temperatura
        k_cyc = self.A_cycle * np.exp(self.Ea_cycle / (self.R_gas * temp_k))

        # Fattore DoD (stress meccanico/chimico sugli elettrodi)
        dod_factor = np.exp(self.b_dod * (1 - dod))

        # Fattore C-rate (stress termico e meccanico)
        crate_factor = 1 + 0.1 * max(0, avg_crate - 0.5)

        loss = k_cyc * dod_factor * crate_factor * n_cycles / 100.0

        return float(np.clip(loss, 0.0, 1.0))

    def predict_soh(self,
                    temp_k: float,
                    time_days: float,
                    n_cycles: int,
                    dod: float = 0.8,
                    avg_soc: float = 0.5,
                    avg_crate: float = 0.3) -> dict:
        """
        Predice SoH combinando calendar e cycle aging.

        Returns:
            dict con SoH, perdita calendar, perdita cicli, RUL stimato
        """
        q_cal = self.calendar_loss_fraction(temp_k, time_days, avg_soc)
        q_cyc = self.cycle_loss_fraction(temp_k, dod, n_cycles, avg_crate)

        # Perdita totale (non lineare - interazione tra i meccanismi)
        q_total = q_cal + q_cyc - 0.3 * q_cal * q_cyc

        current_soh = 1.0 - q_total

        # Stima RUL (cicli rimanenti fino a SoH = 0.8)
        if q_cyc > 0 and n_cycles > 0:
            rate_per_cycle = q_cyc / n_cycles
            remaining_capacity_loss = max(0, current_soh - 0.8)
            rul_cycles = int(remaining_capacity_loss / rate_per_cycle) if rate_per_cycle > 0 else 0
        else:
            rul_cycles = 0

        return {
            'soh': float(np.clip(current_soh, 0.0, 1.0)),
            'soh_percent': float(np.clip(current_soh * 100, 0.0, 100.0)),
            'calendar_loss_pct': q_cal * 100,
            'cycle_loss_pct': q_cyc * 100,
            'total_loss_pct': q_total * 100,
            'rul_cycles': rul_cycles,
            'eol_threshold': 0.8
        }


# Esempio: BESS da 100 MWh in operazione per 2 anni
model = BatteryDegradationModel()

result = model.predict_soh(
    temp_k=298.15,     # 25°C
    time_days=730,     # 2 anni
    n_cycles=730,      # 1 ciclo/giorno
    dod=0.8,           # 80% DoD tipico per grid
    avg_soc=0.5,
    avg_crate=0.3      # 0.3C - tipico BESS 4h
)

print(f"SoH dopo 2 anni: {result['soh_percent']:.1f}%")
print(f"RUL stimato: {result['rul_cycles']} cicli rimanenti")
print(f"  - Calendar loss: {result['calendar_loss_pct']:.1f}%")
print(f"  - Cycle loss: {result['cycle_loss_pct']:.1f}%")

import numpy as np
from scipy.integrate import solve_ivp

class BatteryThermalModel:
    """
    Modello termico lumped-parameter per cella batterica.

    Modello a 2 nodi: core cella + superficie
    dT_core/dt = (Q_gen - (T_core - T_surf) / R_internal) / C_core
    dT_surf/dt = ((T_core - T_surf) / R_internal - (T_surf - T_amb) / R_external) / C_surf

    Riferimento: Bernardi et al. (1985) "A Mathematical Model of the
    Lithium/Polymer Battery" - J. Electrochem. Soc.
    """

    def __init__(self,
                 R_thermal_internal: float = 0.05,  # K/W - resistenza core-superficie
                 R_thermal_external: float = 0.5,   # K/W - resistenza superficie-ambiente
                 C_thermal_core: float = 100.0,      # J/K - capacità termica core
                 C_thermal_surf: float = 10.0,       # J/K - capacità termica superficie
                 cell_resistance_ohm: float = 0.001, # Ohm - resistenza interna per Q_gen
                 t_ambient_c: float = 25.0):

        self.R_int = R_thermal_internal
        self.R_ext = R_thermal_external
        self.C_core = C_thermal_core
        self.C_surf = C_thermal_surf
        self.R_cell = cell_resistance_ohm
        self.T_amb = t_ambient_c + 273.15  # Kelvin

        # Stato iniziale: T_core = T_surf = T_amb
        self.state = np.array([self.T_amb, self.T_amb])

    def _thermal_ode(self, t: float, T: np.ndarray, current_a: float) -> np.ndarray:
        """
        ODE del modello termico.
        T[0] = T_core, T[1] = T_surf (Kelvin)
        """
        T_core, T_surf = T

        # Generazione calore per effetto Joule: Q = I^2 * R
        # Per modello più accurato includere calore entropico
        Q_joule = (current_a ** 2) * self.R_cell
        Q_entropic = 0.0  # Semplificato (può essere negativo in carica LFP)
        Q_gen = Q_joule + Q_entropic

        # Flusso calore core -> superficie
        Q_cs = (T_core - T_surf) / self.R_int

        # Flusso calore superficie -> ambiente (cooling system)
        Q_sa = (T_surf - self.T_amb) / self.R_ext

        dT_core_dt = (Q_gen - Q_cs) / self.C_core
        dT_surf_dt = (Q_cs - Q_sa) / self.C_surf

        return np.array([dT_core_dt, dT_surf_dt])

    def simulate(self,
                 current_profile_a: np.ndarray,
                 dt_s: float = 1.0) -> dict:
        """
        Simula profilo termico per un profilo di corrente.

        Args:
            current_profile_a: Array correnti [A] nel tempo
            dt_s: Passo temporale [s]

        Returns:
            dict con temperature core, superficie e flag di allarme
        """
        n_steps = len(current_profile_a)
        T_core_arr = np.zeros(n_steps)
        T_surf_arr = np.zeros(n_steps)
        alarms = []

        current_state = self.state.copy()

        for i, current in enumerate(current_profile_a):
            # Integrazione numerica
            sol = solve_ivp(
                fun=lambda t, y: self._thermal_ode(t, y, current),
                t_span=(0, dt_s),
                y0=current_state,
                method='RK45',
                max_step=dt_s / 10
            )
            current_state = sol.y[:, -1]

            T_core_c = current_state[0] - 273.15
            T_surf_c = current_state[1] - 273.15

            T_core_arr[i] = T_core_c
            T_surf_arr[i] = T_surf_c

            # Fault detection termico
            alarm = self._check_thermal_faults(i, T_core_c, T_surf_c)
            if alarm:
                alarms.append(alarm)

        self.state = current_state

        return {
            'T_core_c': T_core_arr,
            'T_surf_c': T_surf_arr,
            'T_max_c': float(np.max(T_core_arr)),
            'alarms': alarms,
            'thermal_runaway_risk': float(np.max(T_core_arr)) > 60.0
        }

    def _check_thermal_faults(self,
                               step: int,
                               t_core_c: float,
                               t_surf_c: float) -> Optional[dict]:
        """Verifica soglie di allarme termico."""

        # LFP thresholds - più conservative di NMC
        WARN_TEMP = 45.0
        ALERT_TEMP = 55.0
        CRITICAL_TEMP = 65.0  # Pre-thermal runaway

        if t_core_c > CRITICAL_TEMP:
            return {'step': step, 'level': 'CRITICAL', 'T': t_core_c,
                    'action': 'EMERGENCY_DISCONNECT'}
        elif t_core_c > ALERT_TEMP:
            return {'step': step, 'level': 'ALERT', 'T': t_core_c,
                    'action': 'REDUCE_POWER_50PCT'}
        elif t_core_c > WARN_TEMP:
            return {'step': step, 'level': 'WARNING', 'T': t_core_c,
                    'action': 'INCREASE_COOLING'}
        return None


# Importazione Optional (mancante nell'esempio sopra per brevita)
from typing import Optional

# Simulazione: BESS in carica rapida (1C) a 25°C ambiente
model_thermal = BatteryThermalModel(
    R_thermal_external=0.3,   # Raffreddamento attivo a liquido
    cell_resistance_ohm=0.0005  # LFP 280Ah prismatica
)

# Profilo corrente: 1C charge (280A) per 30 minuti
current_profile = np.full(1800, -280.0)  # Negativo = carica

result = model_thermal.simulate(current_profile, dt_s=1.0)
print(f"Temperatura massima core: {result['T_max_c']:.1f}°C")
print(f"Allarmi generati: {len(result['alarms'])}")
print(f"Rischio thermal runaway: {result['thermal_runaway_risk']}")

from dataclasses import dataclass
from typing import List, Tuple
import numpy as np

@dataclass
class CellState:
    id: int
    voltage_v: float
    soc: float
    temperature_c: float
    capacity_ah: float

class ActiveBalancingController:
    """
    Controllore per active cell balancing con algoritmo SoC-based.
    Obiettivo: equalizzare SoC tra le celle, non la tensione.

    Nota: equalizzare SoC e più corretto di equalizzare tensione
    perchè celle con diverse capacità hanno curve OCV diverse.
    """

    def __init__(self,
                 balancing_current_a: float = 2.0,
                 soc_tolerance: float = 0.02,
                 min_imbalance_for_action: float = 0.03):
        """
        Args:
            balancing_current_a: Corrente di bilanciamento [A]
            soc_tolerance: Tolleranza SoC per considerare celle bilanciate
            min_imbalance_for_action: Imbalance minimo per avviare bilanciamento
        """
        self.I_bal = balancing_current_a
        self.soc_tol = soc_tolerance
        self.min_imbalance = min_imbalance_for_action

    def compute_balancing_plan(self,
                                cells: List[CellState],
                                dt_s: float = 10.0) -> List[dict]:
        """
        Calcola piano di bilanciamento ottimale.

        Algoritmo:
        1. Calcola SoC medio del pack
        2. Identifica celle con SoC sopra media (sorgenti) e sotto media (sink)
        3. Pianifica trasferimenti energia per minimizzare imbalance

        Returns:
            Lista di azioni: {'from_cell': id, 'to_cell': id,
                                  'duration_s': t, 'current_a': I}
        """
        if not cells:
            return []

        soc_values = np.array([c.soc for c in cells])
        soc_mean = np.mean(soc_values)
        max_imbalance = np.max(soc_values) - np.min(soc_values)

        # Non agire se imbalance e trascurabile
        if max_imbalance < self.min_imbalance:
            return []

        actions = []

        # Celle sopra media (sorgenti di energia)
        sources = [(c, c.soc - soc_mean)
                   for c in cells if c.soc - soc_mean > self.soc_tol]
        # Celle sotto media (riceventi di energia)
        sinks = [(c, soc_mean - c.soc)
                 for c in cells if soc_mean - c.soc > self.soc_tol]

        # Ordina per massimo imbalance
        sources.sort(key=lambda x: -x[1])
        sinks.sort(key=lambda x: -x[1])

        # Pianifica trasferimenti (algoritmo greedy)
        for src_cell, src_excess in sources:
            for snk_cell, snk_deficit in sinks:
                if src_excess < self.soc_tol or snk_deficit < self.soc_tol:
                    continue

                # Energia da trasferire (in termini di SoC)
                delta_soc = min(src_excess, snk_deficit) * 0.5  # Conservativo

                # Durata bilanciamento per trasferire delta_soc
                # delta_soc = I * t / (3600 * capacity_ah)
                avg_cap = (src_cell.capacity_ah + snk_cell.capacity_ah) / 2
                duration_s = (delta_soc * avg_cap * 3600) / self.I_bal

                if duration_s > dt_s:  # Azione significativa
                    actions.append({
                        'from_cell': src_cell.id,
                        'to_cell': snk_cell.id,
                        'current_a': self.I_bal,
                        'duration_s': min(duration_s, 300.0),  # Max 5 min per azione
                        'delta_soc': delta_soc
                    })

                src_excess -= delta_soc
                snk_deficit -= delta_soc

        return actions

    def estimate_balancing_time(self, cells: List[CellState]) -> float:
        """
        Stima il tempo necessario per bilanciare completamente il pack [ore].
        """
        soc_values = np.array([c.soc for c in cells])
        max_imbalance = np.max(soc_values) - np.min(soc_values)

        if max_imbalance < self.soc_tol:
            return 0.0

        avg_cap = np.mean([c.capacity_ah for c in cells])

        # Energia da trasferire (Ah) per una cella
        energy_to_transfer_ah = max_imbalance * avg_cap / 2

        # Ore necessarie a corrente I_bal
        hours = energy_to_transfer_ah / self.I_bal

        return float(hours)

from enum import Enum, auto
from dataclasses import dataclass, field
from typing import List, Optional, Callable
import time

class BMSState(Enum):
    """Stati del BMS - solo le transizioni permesse sono valide"""
    INIT           = auto()  # Avvio, auto-test
    STANDBY        = auto()  # Pronto, rete connessa, nessun flusso energetico
    PRECHARGE      = auto()  # Pre-carica condensatori (evita inrush)
    OPERATIONAL    = auto()  # Operativo normale (carica o scarica)
    CHARGING       = auto()  # In carica
    DISCHARGING    = auto()  # In scarica
    BALANCING      = auto()  # Cell balancing attivo
    FAULT_SOFT     = auto()  # Guasto recuperabile (overtemp warning, etc.)
    FAULT_HARD     = auto()  # Guasto critico - richiede reset manuale
    EMERGENCY_STOP = auto()  # Arresto di emergenza - disconnessione fisica
    SHUTDOWN       = auto()  # Spegnimento controllato

@dataclass
class FaultCode:
    code: str
    description: str
    severity: str   # 'WARNING', 'SOFT_FAULT', 'HARD_FAULT', 'EMERGENCY'
    recoverable: bool
    action: str

# Registry dei fault codes
FAULT_REGISTRY = {
    'OV_CELL':   FaultCode('OV_CELL',   'Cell overvoltage',        'HARD_FAULT', False, 'OPEN_CONTACTOR'),
    'UV_CELL':   FaultCode('UV_CELL',   'Cell undervoltage',       'HARD_FAULT', False, 'OPEN_CONTACTOR'),
    'OT_CELL':   FaultCode('OT_CELL',   'Cell overtemperature',    'HARD_FAULT', False, 'OPEN_CONTACTOR'),
    'OT_WARN':   FaultCode('OT_WARN',   'Temperature warning',     'WARNING',    True,  'REDUCE_POWER'),
    'OC_CHARGE': FaultCode('OC_CHARGE', 'Overcurrent in charge',   'HARD_FAULT', False, 'OPEN_CONTACTOR'),
    'OC_DISC':   FaultCode('OC_DISC',   'Overcurrent discharge',   'HARD_FAULT', False, 'OPEN_CONTACTOR'),
    'ISOLATION': FaultCode('ISOLATION', 'Isolation fault',         'EMERGENCY',  False, 'EMERGENCY_STOP'),
    'COMM_FAIL': FaultCode('COMM_FAIL', 'Communication failure',   'SOFT_FAULT', True,  'STANDBY'),
    'SOC_LOW':   FaultCode('SOC_LOW',   'SoC below minimum',       'SOFT_FAULT', True,  'STOP_DISCHARGE'),
    'SOC_HIGH':  FaultCode('SOC_HIGH',  'SoC above maximum',       'SOFT_FAULT', True,  'STOP_CHARGE'),
    'TR_DETECT': FaultCode('TR_DETECT', 'Thermal runaway detected', 'EMERGENCY', False, 'EMERGENCY_STOP'),
}

class BMSSafetyStateMachine:
    """
    Safety state machine per BMS grid-scale.
    Implementa le protezioni definite in IEC 62619 e IEEE 1625.
    """

    # Soglie di protezione (configurabili per chimica)
    PROTECTION_THRESHOLDS = {
        'LFP': {
            'cell_ov_v': 3.65,      # Overvoltage cutoff
            'cell_uv_v': 2.80,      # Undervoltage cutoff
            'cell_ov_warn_v': 3.60, # Overvoltage warning
            'cell_uv_warn_v': 2.90, # Undervoltage warning
            'temp_ot_c': 55.0,      # Overtemperature cutoff
            'temp_ot_warn_c': 45.0, # Overtemperature warning
            'temp_ut_c': -10.0,     # Undertemperature cutoff (charge only)
            'oc_charge_a': 1.0,     # Max C-rate charge (per unit)
            'oc_disc_a': 2.0,       # Max C-rate discharge (per unit)
        },
        'NMC': {
            'cell_ov_v': 4.20,
            'cell_uv_v': 3.00,
            'cell_ov_warn_v': 4.15,
            'cell_uv_warn_v': 3.10,
            'temp_ot_c': 50.0,
            'temp_ot_warn_c': 40.0,
            'temp_ut_c': 0.0,
            'oc_charge_a': 0.7,
            'oc_disc_a': 1.5,
        }
    }

    def __init__(self, chemistry: str = 'LFP',
                 on_state_change: Optional[Callable] = None):
        self.state = BMSState.INIT
        self.thresholds = self.PROTECTION_THRESHOLDS[chemistry]
        self.active_faults: List[FaultCode] = []
        self.fault_history: List[dict] = []
        self.on_state_change = on_state_change
        self._contactor_closed = False

    def check_and_transition(self, telemetry: dict) -> BMSState:
        """
        Verifica telemetria e aggiorna stato.

        Args:
            telemetry: dict con cell_voltages, cell_temps, pack_current, soc, etc.

        Returns:
            Nuovo stato BMS
        """
        faults = self._detect_faults(telemetry)

        if faults:
            self._handle_faults(faults)
        else:
            # Rimozione fault recuperabili se condizione normalizzata
            self._clear_recoverable_faults(telemetry)

        return self.state

    def _detect_faults(self, tel: dict) -> List[FaultCode]:
        """Rileva fault attivi nella telemetria."""
        detected = []
        thr = self.thresholds

        # 1. Cell Voltage Checks
        for v in tel.get('cell_voltages', []):
            if v > thr['cell_ov_v']:
                detected.append(FAULT_REGISTRY['OV_CELL'])
                break
            if v < thr['cell_uv_v']:
                detected.append(FAULT_REGISTRY['UV_CELL'])
                break

        # 2. Temperature Checks
        for t in tel.get('cell_temps', []):
            if t > thr['temp_ot_c']:
                detected.append(FAULT_REGISTRY['OT_CELL'])
                break
            if t > thr['temp_ot_warn_c']:
                detected.append(FAULT_REGISTRY['OT_WARN'])
                break

        # 3. Isolation Fault (impedance monitor)
        if tel.get('isolation_ok', True) == False:
            detected.append(FAULT_REGISTRY['ISOLATION'])

        # 4. Thermal runaway detection (gas sensor, rapid T rise)
        if self._detect_thermal_runaway(tel):
            detected.append(FAULT_REGISTRY['TR_DETECT'])

        return detected

    def _detect_thermal_runaway(self, tel: dict) -> bool:
        """
        Rilevamento precoce thermal runaway:
        - Rate of temperature rise > 1°C/s (early stage)
        - Gas sensor (CO, H2, elettrolita VOC) attivato
        - Calo improvviso tensione cella con surriscaldamento
        """
        t_rate = tel.get('temp_rate_c_per_s', 0.0)
        gas_alarm = tel.get('gas_sensor_alarm', False)

        return t_rate > 1.0 or gas_alarm

    def _handle_faults(self, faults: List[FaultCode]) -> None:
        """Gestisce i fault rilevati con transizione di stato appropriata."""

        # Priorità: EMERGENCY > HARD_FAULT > SOFT_FAULT > WARNING
        max_severity = max(faults, key=lambda f: [
            'WARNING', 'SOFT_FAULT', 'HARD_FAULT', 'EMERGENCY'
        ].index(f.severity))

        # Log fault
        for fault in faults:
            if fault not in self.active_faults:
                self.active_faults.append(fault)
                self.fault_history.append({
                    'timestamp': time.time(),
                    'code': fault.code,
                    'severity': fault.severity
                })

        # Transizione di stato
        if max_severity.severity == 'EMERGENCY':
            self._transition_to(BMSState.EMERGENCY_STOP)
            self._open_main_contactor(emergency=True)

        elif max_severity.severity == 'HARD_FAULT':
            self._transition_to(BMSState.FAULT_HARD)
            self._open_main_contactor(emergency=False)

        elif max_severity.severity == 'SOFT_FAULT':
            if self.state not in (BMSState.FAULT_HARD, BMSState.EMERGENCY_STOP):
                self._transition_to(BMSState.FAULT_SOFT)

    def _transition_to(self, new_state: BMSState) -> None:
        if new_state != self.state:
            old_state = self.state
            self.state = new_state
            if self.on_state_change:
                self.on_state_change(old_state, new_state)

    def _open_main_contactor(self, emergency: bool) -> None:
        """Apre il contattore principale (fisico)."""
        self._contactor_closed = False
        # In produzione: comando hardware GPIO/CAN al driver contattore

    def _clear_recoverable_faults(self, tel: dict) -> None:
        """Rimuove fault recuperabili se condizione normalizzata."""
        thr = self.thresholds

        self.active_faults = [
            f for f in self.active_faults
            if not f.recoverable
        ]

        if not self.active_faults and self.state == BMSState.FAULT_SOFT:
            self._transition_to(BMSState.STANDBY)

import numpy as np
from dataclasses import dataclass
from typing import Tuple

@dataclass
class ChargingOptimizer:
    """
    Ottimizzazione strategia di ricarica per massimizzare cicli vita.

    Strategie implementate:
    1. CC-CV (Constant Current - Constant Voltage) standard
    2. Multi-step CC (riduzione degrado agli elettrodi)
    3. Pulse charging (riduzione stress termico)
    """

    capacity_ah: float
    cell_voltage_max: float   # Es. 3.65V per LFP
    cell_voltage_min: float   # Es. 2.80V per LFP
    soc_max: float = 0.95     # Non caricare oltre 95%
    soc_min: float = 0.05     # Non scaricare sotto 5%

    def cc_cv_profile(self,
                       target_soc: float,
                       current_soc: float,
                       max_crate: float = 0.5) -> dict:
        """
        Genera profilo di ricarica CC-CV ottimizzato.

        Fase CC: ricarica a corrente costante fino a V_max - 50mV
        Fase CV: mantiene tensione costante, corrente decresce
        Terminazione: quando corrente scende sotto C/20
        """
        if current_soc >= target_soc:
            return {'phase': 'COMPLETE', 'current_a': 0, 'voltage_v': 0}

        soc_delta = target_soc - current_soc

        # Calcola corrente CC ottimale basata su soc_delta e temperatura
        # Riduzione corrente per SoC alto (vicino alla fine carica)
        if current_soc < 0.8:
            cc_crate = max_crate
        elif current_soc < 0.9:
            cc_crate = max_crate * 0.7  # Rallenta a 80%
        else:
            cc_crate = max_crate * 0.3  # CV-like region

        cc_current = cc_crate * self.capacity_ah

        # Tensione target (con margine di sicurezza per cell balancing)
        v_target = self.cell_voltage_max - 0.05  # 50 mV di margine

        return {
            'phase': 'CC' if current_soc < 0.9 else 'CV',
            'current_a': cc_current,
            'voltage_v': v_target,
            'estimated_minutes': (soc_delta * self.capacity_ah) / cc_current * 60
        }

    def calculate_optimal_dod(self,
                               daily_cycles: float,
                               target_years: float,
                               chemistry: str = 'LFP') -> dict:
        """
        Calcola il DoD ottimale per massimizzare l'energia totale throughput
        nell'arco della vita target.

        Trade-off: più DoD = più energia per ciclo ma meno cicli totali
        Ottimale = massimo di (DoD * cicli_a_quel_DoD)
        """
        # Modello empirico cicli vs DoD (semplificato)
        CYCLES_AT_DOD = {
            'LFP': {0.5: 8000, 0.6: 6500, 0.7: 5500, 0.8: 4500, 0.9: 3500, 1.0: 2500},
            'NMC': {0.5: 3000, 0.6: 2400, 0.7: 2000, 0.8: 1600, 0.9: 1200, 1.0: 900}
        }

        cycles_map = CYCLES_AT_DOD.get(chemistry, CYCLES_AT_DOD['LFP'])
        dod_values = sorted(cycles_map.keys())

        results = []
        required_cycles = daily_cycles * 365 * target_years

        for dod in dod_values:
            total_cycles = cycles_map[dod]
            energy_per_cycle_rel = dod  # Relativo
            total_energy_rel = total_cycles * energy_per_cycle_rel

            # La batteria regge i cicli richiesti?
            years_of_life = total_cycles / (daily_cycles * 365)

            results.append({
                'dod': dod,
                'total_cycles': total_cycles,
                'years_of_life': years_of_life,
                'total_energy_throughput': total_energy_rel,
                'meets_target': years_of_life >= target_years
            })

        # Ottimale: massimo energy throughput che soddisfa target anni
        valid = [r for r in results if r['meets_target']]
        optimal = max(valid, key=lambda x: x['total_energy_throughput']) if valid else results[-1]

        return {
            'optimal_dod': optimal['dod'],
            'expected_years': optimal['years_of_life'],
            'total_cycles_available': optimal['total_cycles'],
            'all_scenarios': results
        }


# Esempio
optimizer = ChargingOptimizer(
    capacity_ah=280.0,
    cell_voltage_max=3.65,
    cell_voltage_min=2.80
)

result = optimizer.calculate_optimal_dod(
    daily_cycles=1.5,   # 1.5 cicli/giorno (tipico arbitraggio + frequency reg)
    target_years=15.0,  # Vita target 15 anni
    chemistry='LFP'
)

print(f"DoD ottimale: {result['optimal_dod']*100:.0f}%")
print(f"Vita attesa: {result['expected_years']:.1f} anni")

from dataclasses import dataclass
from typing import List, Optional
import numpy as np

@dataclass
class GridDispatchCommand:
    """Comando di dispatch dalla rete o dall'EMS"""
    power_kw: float          # Positivo = scarica verso rete, negativo = carica da rete
    duration_s: int
    service_type: str        # 'FCR', 'aFRR', 'peak_shaving', 'arbitrage', 'ramp_control'
    priority: int            # 1 = massima priorità (safety), 10 = minima
    timestamp: float

class BESSGridDispatcher:
    """
    Dispatcher per BESS che gestisce comandi dalla rete
    rispettando i vincoli BMS (SoC, temperatura, fault state).

    Integra con EMS tramite Modbus TCP / IEC 61850 XCBR/MMXU
    """

    def __init__(self,
                 power_max_kw: float,
                 capacity_kwh: float,
                 soc_min: float = 0.1,
                 soc_max: float = 0.95,
                 ramp_rate_kw_per_s: float = None):

        self.P_max = power_max_kw
        self.E_total = capacity_kwh
        self.soc_min = soc_min
        self.soc_max = soc_max
        # Default ramp rate: full power in 1 secondo (tipico BESS moderno)
        self.ramp_rate = ramp_rate_kw_per_s or power_max_kw

        self._current_power_kw = 0.0
        self._current_soc = 0.5
        self._bms_state = 'OPERATIONAL'

    def execute_command(self,
                         cmd: GridDispatchCommand,
                         bms_telemetry: dict) -> dict:
        """
        Esegue un comando di dispatch verificando i vincoli BMS.

        Returns:
            {'executed_power_kw': float, 'curtailed': bool,
              'reason': str, 'available_energy_kwh': float}
        """
        self._current_soc = bms_telemetry.get('soc', self._current_soc)
        self._bms_state = bms_telemetry.get('state', self._bms_state)

        # 1. Verifica stato BMS
        if self._bms_state in ('FAULT_HARD', 'EMERGENCY_STOP'):
            return {
                'executed_power_kw': 0,
                'curtailed': True,
                'reason': f'BMS in stato {self._bms_state}',
                'available_energy_kwh': 0
            }

        # 2. Calcola potenza permessa con vincoli SoC
        requested_p = cmd.power_kw

        if requested_p > 0:  # Scarica
            # Energia disponibile sopra SoC minimo
            available_energy = max(0,
                (self._current_soc - self.soc_min) * self.E_total)

            # Potenza massima che non porta a SoC < soc_min nella durata
            p_max_soc = (available_energy / cmd.duration_s) * 3600

            max_discharge = min(self.P_max, p_max_soc)

            if requested_p > max_discharge:
                executed_p = max_discharge
                curtailed = True
                reason = f'SoC troppo basso: disponibili {available_energy:.1f} kWh'
            else:
                executed_p = requested_p
                curtailed = False
                reason = 'OK'

        else:  # Carica (potenza negativa)
            # capacità disponibile sotto SoC massimo
            available_cap = max(0,
                (self.soc_max - self._current_soc) * self.E_total)

            p_max_soc = (available_cap / cmd.duration_s) * 3600
            p_requested_abs = abs(requested_p)

            max_charge = min(self.P_max, p_max_soc)

            if p_requested_abs > max_charge:
                executed_p = -max_charge
                curtailed = True
                reason = f'SoC troppo alto: disponibili {available_cap:.1f} kWh'
            else:
                executed_p = requested_p
                curtailed = False
                reason = 'OK'

        # 3. Applica ramp rate limiting
        power_delta = executed_p - self._current_power_kw
        max_delta = self.ramp_rate * 0.1  # 100 ms step

        if abs(power_delta) > max_delta:
            executed_p = self._current_power_kw + np.sign(power_delta) * max_delta

        self._current_power_kw = executed_p

        return {
            'executed_power_kw': executed_p,
            'curtailed': curtailed,
            'reason': reason,
            'available_energy_kwh': abs(
                (self.soc_max - self._current_soc) * self.E_total
                if executed_p < 0
                else (self._current_soc - self.soc_min) * self.E_total
            )
        }

    def frequency_regulation_response(self,
                                        grid_freq_hz: float,
                                        nominal_freq_hz: float = 50.0,
                                        deadband_hz: float = 0.010) -> float:
        """
        Risposta automatica alla frequenza di rete (FCR - Frequency Containment Reserve).

        Regola europea: risposta proporzionale lineare tra ±200 mHz,
        potenza massima oltre ±200 mHz (IEC 61000-4-30).

        Returns:
            Potenza di risposta [kW] (positiva = iniezione in rete)
        """
        freq_deviation = grid_freq_hz - nominal_freq_hz

        # Deadband: nessuna risposta per deviazioni minori
        if abs(freq_deviation) <= deadband_hz:
            return 0.0

        # Risposta proporzionale (droop)
        effective_deviation = freq_deviation - np.sign(freq_deviation) * deadband_hz

        # Range proporzionale: ±200 mHz = ±100% potenza
        droop_range_hz = 0.200
        droop_response = np.clip(
            effective_deviation / droop_range_hz, -1.0, 1.0
        )

        # Inverti: frequenza bassa = rete ha bisogno di potenza = scarica BESS
        response_power = -droop_response * self.P_max

        return float(response_power)

# Stack tecnologico BMS grid-scale - 2025

BMS_TECH_STACK = {

    # LAYER 1: Cell Monitoring IC (Hardware)
    'cell_monitoring': {
        'vendors': ['Texas Instruments BQ76952', 'Analog Devices ADBMS6815',
                    'Renesas ISL94212', 'NXP MC33771'],
        'voltage_accuracy': '±0.5-2 mV',
        'current_integration': 'Shunt or Hall-effect sensor',
        'interface': 'SPI / isoSPI / CAN',
        'isolation': 'Galvanic (up to 1500V DC)'
    },

    # LAYER 2: Cell Controller MCU (Firmware)
    'cell_controller': {
        'hw': ['ST STM32H7', 'NXP S32K3', 'Renesas RH850'],
        'os': ['FreeRTOS', 'AUTOSAR CP', 'Bare Metal'],
        'language': 'C99/C11',
        'cycle_time': '1-10 ms',
        'standards': ['ISO 26262 (ASIL-D per EV)', 'IEC 61508 (SIL-2 per grid)']
    },

    # LAYER 3: BMS Controller (Edge Computing)
    'bms_controller': {
        'hw': ['Raspberry Pi CM4 Industrial', 'Kontron KBox A-202',
               'Beckhoff CX5200', 'NVIDIA Jetson (per ML)'],
        'os': 'Linux (PREEMPT-RT kernel)',
        'language': 'Python 3.11 + C extensions',
        'key_libs': ['NumPy', 'SciPy', 'filterpy (EKF)',
                     'scikit-learn', 'asyncio'],
        'comms': ['CANopen', 'Modbus RTU/TCP', 'EtherCAT'],
        'protocols': ['IEC 61850', 'OCPP 2.0.1 (per EV)']
    },

    # LAYER 4: System EMS (Server)
    'energy_management': {
        'platform': ['Python FastAPI', 'Node.js', 'Java Spring Boot'],
        'database': ['InfluxDB (timeseries)', 'PostgreSQL (config)',
                     'Redis (real-time cache)'],
        'message_broker': ['Apache Kafka', 'MQTT (EMQX)'],
        'grid_protocols': ['Modbus TCP', 'IEC 61850', 'DNP3', 'SunSpec'],
        'monitoring': ['Grafana', 'Prometheus', 'Victoria Metrics']
    },

    # LAYER 5: Cloud Analytics
    'cloud_analytics': {
        'platform': ['AWS IoT TwinMaker', 'Azure IoT Hub', 'GCP IoT Core'],
        'ml_platform': ['MLflow', 'Ray', 'TensorFlow Lite (edge inference)'],
        'analytics': ['Apache Spark (batch)', 'Apache Flink (streaming)'],
        'digital_twin': ['AWS IoT TwinMaker', 'Bentley iTwin', 'AVEVA PI']
    }
}

# Configurazione Modbus per comunicazione BMS-EMS
MODBUS_REGISTER_MAP = {
    # Input Registers (read-only)
    1000: ('soc_percent', 'uint16', 'x100'),        # SoC: 0-10000 = 0-100.00%
    1001: ('soh_percent', 'uint16', 'x100'),         # SoH: 0-10000 = 0-100.00%
    1002: ('pack_voltage', 'uint16', 'x10'),          # V: 0-65535 = 0-6553.5V
    1003: ('pack_current', 'int16', 'x10'),           # A: -32768-32767 = -3276.8 to 3276.7A
    1004: ('max_cell_temp', 'int16', 'x10'),          # °C: -500 to +1000 = -50.0 to 100.0°C
    1005: ('bms_state', 'uint16', 'enum'),            # 0=INIT, 1=STANDBY, ..., 9=EMERGENCY
    1006: ('active_faults_bitmask', 'uint32', 'bits'), # Bit per fault attivo
    1008: ('available_power_kw', 'int16', 'x1'),      # kW disponibile (pos=scarica, neg=carica)

    # Holding Registers (read-write)
    2000: ('power_setpoint_kw', 'int16', 'x1'),       # Setpoint potenza da EMS
    2001: ('charge_enable', 'uint16', 'bool'),         # 1 = abilita carica
    2002: ('discharge_enable', 'uint16', 'bool'),      # 1 = abilita scarica
    2003: ('soc_setpoint_percent', 'uint16', 'x100'), # SoC target per EMS
}