simulation_engine.py

"""Simulation engine for charge/discharge cycles and termination handling."""

from __future__ import annotations

import csv
import json
import math
from datetime import datetime
from dataclasses import dataclass, field
from pathlib import Path
from typing import Any, Callable, List

from pathsim import Simulation


UPPER_VOLTAGE_CUTOFF_KEY = "Upper voltage cut-off [V]"


def round_up_to_0p01(value: float) -> float:
    """Round up to the nearest 0.01V."""
    return math.ceil(float(value) * 100.0) / 100.0


def derive_charge_voltage_limit(charge_voltage_target_v: float) -> float:
    """Derive controller voltage limit as +1% target rounded up to 0.01V."""
    return round_up_to_0p01(float(charge_voltage_target_v) * 1.01)


def derive_upper_voltage_cutoff(charge_voltage_target_v: float) -> float:
    """Derive PyBaMM upper voltage cut-off from charge target."""
    return float(charge_voltage_target_v) + 0.10


def apply_derived_upper_voltage_cutoff(
    parameter_overrides: dict[str, float] | None,
    charge_voltage_target_v: float,
) -> dict[str, float]:
    """Return overrides with derived upper cut-off injected when missing.

    Explicit override is always preserved.
    """
    out: dict[str, float] = {} if parameter_overrides is None else dict(parameter_overrides)
    if UPPER_VOLTAGE_CUTOFF_KEY not in out:
        out[UPPER_VOLTAGE_CUTOFF_KEY] = derive_upper_voltage_cutoff(charge_voltage_target_v)
    return out


@dataclass
class SimulationResult:
    """Container for traces and termination metadata."""

    time_trace: List[float] = field(default_factory=list)
    ia_trace: List[float] = field(default_factory=list)
    ib_trace: List[float] = field(default_factory=list)
    va_trace: List[float] = field(default_factory=list)
    vb_trace: List[float] = field(default_factory=list)
    soc_a_trace: List[float] = field(default_factory=list)
    soc_b_trace: List[float] = field(default_factory=list)
    cycle_index_trace: List[int] = field(default_factory=list)
    stop_reason: str = ""
    cycle_reasons: List[str] = field(default_factory=list)
    completed_cycles: int = 0
    full_cycle_count_requested: int = 0
    event_count_requested: int = 0
    event_count_executed: int = 0
    full_cycles_completed: int = 0
    cycle_summaries: List[dict[str, Any]] = field(default_factory=list)
    aging_summary: dict[str, Any] = field(default_factory=dict)
    aging_artifacts: dict[str, str] = field(default_factory=dict)


@dataclass
class SimulationConfig:
    """Runtime simulation configuration for charge/discharge control limits."""

    charge_voltage_target_v: float = 4.15
    charge_voltage_limit_v: float | None = None
    discharge_voltage_cutoff_v: float = 2.5
    discharge_voltage_limit_v: float | None = None
    cv_entry_v_per_cell: float | None = None
    safety_max_cell_v: float | None = None
    cv_taper_current_threshold_a_fraction: float = 0.05
    cv_safety_timeout_s: float = 43200.0
    charge_cv_min_current_a: float = 0.01
    charge_cv_margin_v: float | None = None
    cv_entry_abs_tol_v: float = 1e-9

    def __post_init__(self) -> None:
        self.charge_voltage_target_v = float(self.charge_voltage_target_v)
        if self.charge_voltage_limit_v is None:
            self.charge_voltage_limit_v = derive_charge_voltage_limit(self.charge_voltage_target_v)
        else:
            self.charge_voltage_limit_v = float(self.charge_voltage_limit_v)
        self.discharge_voltage_cutoff_v = float(self.discharge_voltage_cutoff_v)
        if self.safety_max_cell_v is None:
            self.safety_max_cell_v = float(self.charge_voltage_limit_v)
        else:
            self.safety_max_cell_v = float(self.safety_max_cell_v)
        self.cv_taper_current_threshold_a_fraction = float(self.cv_taper_current_threshold_a_fraction)
        self.cv_safety_timeout_s = float(self.cv_safety_timeout_s)
        self.charge_cv_min_current_a = float(self.charge_cv_min_current_a)
        self.cv_entry_abs_tol_v = float(self.cv_entry_abs_tol_v)

        if self.cv_entry_v_per_cell is None:
            self.cv_entry_v_per_cell = self.charge_voltage_target_v
        else:
            self.cv_entry_v_per_cell = float(self.cv_entry_v_per_cell)

        if self.discharge_voltage_limit_v is None:
            self.discharge_voltage_limit_v = self.discharge_voltage_cutoff_v
        else:
            self.discharge_voltage_limit_v = float(self.discharge_voltage_limit_v)

        if self.charge_cv_margin_v is None:
            self.charge_cv_margin_v = max(1e-6, self.safety_max_cell_v - self.charge_voltage_target_v)
        else:
            self.charge_cv_margin_v = max(1e-6, float(self.charge_cv_margin_v))


def _capacity_from_soc_window(delta_soc_pct: float, throughput_ah: float) -> float | None:
    """Estimate full-capacity equivalent from throughput and SOC window."""
    soc_window_fraction = abs(delta_soc_pct) / 100.0
    if soc_window_fraction <= 1e-12:
        return None
    return throughput_ah / soc_window_fraction


def _build_cycle_summaries(
    result: SimulationResult,
    pack_current: float,
    alternate_charge_discharge: bool,
    nominal_capacity_a: float,
    nominal_capacity_b: float,
) -> list[dict[str, Any]]:
    """Build cycle-level FCC/SOH summary rows from recorded traces."""
    rows: list[dict[str, Any]] = []
    if not result.time_trace or not result.cycle_index_trace:
        return rows

    nominal_capacity_a = max(float(nominal_capacity_a), 1e-12)
    nominal_capacity_b = max(float(nominal_capacity_b), 1e-12)
    nominal_capacity_pack = nominal_capacity_a + nominal_capacity_b

    max_cycle_num = max(int(x) for x in result.cycle_index_trace)
    for cycle_num in range(1, max_cycle_num + 1):
        indices = [i for i, c in enumerate(result.cycle_index_trace) if int(c) == cycle_num]
        if not indices:
            continue

        i_start = indices[0]
        i_end = indices[-1]

        start_soc_a = float(result.soc_a_trace[i_start])
        end_soc_a = float(result.soc_a_trace[i_end])
        start_soc_b = float(result.soc_b_trace[i_start])
        end_soc_b = float(result.soc_b_trace[i_end])

        delta_soc_a_pct = end_soc_a - start_soc_a
        delta_soc_b_pct = end_soc_b - start_soc_b

        throughput_a_ah = abs(delta_soc_a_pct) / 100.0 * nominal_capacity_a
        throughput_b_ah = abs(delta_soc_b_pct) / 100.0 * nominal_capacity_b
        throughput_pack_ah = throughput_a_ah + throughput_b_ah

        fcc_estimate_a_ah = _capacity_from_soc_window(delta_soc_a_pct, throughput_a_ah)
        fcc_estimate_b_ah = _capacity_from_soc_window(delta_soc_b_pct, throughput_b_ah)

        soc_window_pack_fraction = 0.0
        if nominal_capacity_pack > 0:
            soc_window_pack_fraction = (
                (abs(delta_soc_a_pct) / 100.0) * nominal_capacity_a
                + (abs(delta_soc_b_pct) / 100.0) * nominal_capacity_b
            ) / nominal_capacity_pack
        if soc_window_pack_fraction > 1e-12:
            fcc_estimate_pack_ah = throughput_pack_ah / soc_window_pack_fraction
        else:
            fcc_estimate_pack_ah = None

        soh_a_pct = None if fcc_estimate_a_ah is None else (fcc_estimate_a_ah / nominal_capacity_a) * 100.0
        soh_b_pct = None if fcc_estimate_b_ah is None else (fcc_estimate_b_ah / nominal_capacity_b) * 100.0
        soh_pack_pct = None
        if fcc_estimate_pack_ah is not None and nominal_capacity_pack > 0:
            soh_pack_pct = (fcc_estimate_pack_ah / nominal_capacity_pack) * 100.0

        signed_current = _current_for_cycle(pack_current, cycle_num - 1, alternate_charge_discharge)
        direction = "Discharge" if signed_current > 0.0 else "Charge"

        cycle_reason = ""
        if cycle_num - 1 < len(result.cycle_reasons):
            cycle_reason = result.cycle_reasons[cycle_num - 1]

        rows.append(
            {
                "cycle_index": cycle_num,
                "direction": direction,
                "time_start_s": float(result.time_trace[i_start]),
                "time_end_s": float(result.time_trace[i_end]),
                "duration_s": float(result.time_trace[i_end] - result.time_trace[i_start]),
                "start_soc_a_pct": start_soc_a,
                "end_soc_a_pct": end_soc_a,
                "delta_soc_a_pct": float(delta_soc_a_pct),
                "start_soc_b_pct": start_soc_b,
                "end_soc_b_pct": end_soc_b,
                "delta_soc_b_pct": float(delta_soc_b_pct),
                "throughput_a_ah": float(throughput_a_ah),
                "throughput_b_ah": float(throughput_b_ah),
                "throughput_pack_ah": float(throughput_pack_ah),
                "fcc_estimate_a_ah": fcc_estimate_a_ah,
                "fcc_estimate_b_ah": fcc_estimate_b_ah,
                "fcc_estimate_pack_ah": fcc_estimate_pack_ah,
                "soh_a_pct": soh_a_pct,
                "soh_b_pct": soh_b_pct,
                "soh_pack_pct": soh_pack_pct,
                "end_va_v": float(result.va_trace[i_end]),
                "end_vb_v": float(result.vb_trace[i_end]),
                "termination_reason": cycle_reason,
            }
        )

    return rows


def _build_aging_summary(cycle_rows: list[dict[str, Any]]) -> dict[str, Any]:
    """Build aggregate FCC/SOH reporting metrics for aging-study observability."""
    if not cycle_rows:
        return {
            "cycle_count": 0,
            "fcc_tracking_available": False,
            "soh_tracking_available": False,
            "notes": "No cycle data recorded.",
        }

    def _series(key: str) -> list[float]:
        vals: list[float] = []
        for row in cycle_rows:
            v = row.get(key)
            if isinstance(v, (int, float)):
                vals.append(float(v))
        return vals

    fcc_pack = _series("fcc_estimate_pack_ah")
    soh_pack = _series("soh_pack_pct")

    def _consistency(vals: list[float]) -> dict[str, Any]:
        if not vals:
            return {"count": 0, "min": None, "max": None, "span": None, "mean": None}
        vmin = min(vals)
        vmax = max(vals)
        return {
            "count": len(vals),
            "min": vmin,
            "max": vmax,
            "span": vmax - vmin,
            "mean": sum(vals) / len(vals),
        }

    return {
        "cycle_count": len(cycle_rows),
        "fcc_tracking_available": bool(fcc_pack),
        "soh_tracking_available": bool(soh_pack),
        "fcc_pack_stats": _consistency(fcc_pack),
        "soh_pack_stats": _consistency(soh_pack),
        "reporting_columns": [
            "cycle_index",
            "direction",
            "throughput_pack_ah",
            "fcc_estimate_pack_ah",
            "soh_pack_pct",
            "termination_reason",
        ],
    }


def _write_aging_artifacts(cycle_rows: list[dict[str, Any]], aging_summary: dict[str, Any]) -> dict[str, str]:
    """Write per-cycle FCC/SOH summary artifacts for aging-study reporting."""
    if not cycle_rows:
        return {}

    out_dir = Path("scripts/output") / f"aging_cycle_summary_{datetime.now().strftime('%Y%m%d_%H%M%S')}"
    out_dir.mkdir(parents=True, exist_ok=True)

    csv_path = out_dir / "cycle_summary.csv"
    json_path = out_dir / "aging_summary.json"

    fieldnames = list(cycle_rows[0].keys())
    with csv_path.open("w", newline="") as fp:
        writer = csv.DictWriter(fp, fieldnames=fieldnames)
        writer.writeheader()
        writer.writerows(cycle_rows)

    with json_path.open("w") as fp:
        json.dump({"aging_summary": aging_summary, "cycle_rows": cycle_rows}, fp, indent=2)

    return {
        "output_dir": str(out_dir),
        "cycle_summary_csv": str(csv_path),
        "aging_summary_json": str(json_path),
    }


def _record_step(
    sim: Simulation,
    pack_block,
    dt: float,
    elapsed: float,
    result: SimulationResult,
    cycle_index: int,
) -> float:
    """Advance one timestep and append outputs into trace arrays."""
    sim.timestep(dt=dt, adaptive=False)
    elapsed += dt

    ia = float(pack_block.outputs[0])
    ib = float(pack_block.outputs[1])
    va = float(pack_block.cell_a.get_voltage())
    vb = float(pack_block.cell_b.get_voltage())
    soc_a = float(pack_block.outputs[2])
    soc_b = float(pack_block.outputs[3])

    result.time_trace.append(elapsed)
    result.ia_trace.append(ia)
    result.ib_trace.append(ib)
    result.va_trace.append(va)
    result.vb_trace.append(vb)
    result.soc_a_trace.append(soc_a)
    result.soc_b_trace.append(soc_b)
    result.cycle_index_trace.append(int(cycle_index))

    return elapsed


# ===== Charge/CV Control Defaults (legacy constants) =====
# Runtime behavior now uses SimulationConfig passed to run_simulation.
CHARGE_VOLTAGE_TARGET = 4.15
CV_ENTRY_V_PER_CELL = CHARGE_VOLTAGE_TARGET
SAFETY_MAX_CELL_V = 4.2
CV_TAPER_CURRENT_THRESHOLD_A_FRACTION = 0.05
CV_SAFETY_TIMEOUT_S = 43200.0
_CHARGE_CV_MIN_CURRENT_A = 0.01
_CHARGE_CV_PER_CELL_MARGIN_V = max(1e-6, SAFETY_MAX_CELL_V - CHARGE_VOLTAGE_TARGET)
_DISCHARGE_VOLTAGE_CUTOFF_V = 2.5
_CV_ENTRY_ABS_TOL_V = 1e-9


def _at_or_above_with_tolerance(value: float, threshold: float, abs_tol: float) -> bool:
    """Return True when value is at/above threshold, including near-equality float noise."""
    return value >= (threshold - abs_tol) or math.isclose(value, threshold, rel_tol=0.0, abs_tol=abs_tol)


def _calculate_power_cv_charge_current(
    va: float,
    vb: float,
    cc_current: float,
    charge_voltage_target_v: float,
    charge_cv_margin_v: float,
    charge_cv_min_current_a: float,
    gamma: float = 1.5,
) -> float:
    """Calculate CV charge current with power-law shaping near target."""
    v_avg = (va + vb) / 2.0
    voltage_error = charge_voltage_target_v - v_avg
    voltage_window = charge_cv_margin_v

    if voltage_window <= 0:
        return 0.0

    ratio = voltage_error / voltage_window
    if ratio <= 0.0:
        i_cv = 0.0
    else:
        i_cv = abs(cc_current) * (min(1.0, ratio) ** gamma)

    return max(i_cv, charge_cv_min_current_a)


def _calculate_quadratic_cv_charge_current(
    va: float,
    vb: float,
    cc_current: float,
    charge_voltage_target_v: float,
    charge_cv_margin_v: float,
    charge_cv_min_current_a: float,
) -> float:
    """Calculate CV charge current with quadratic shaping near target."""
    return _calculate_power_cv_charge_current(
        va,
        vb,
        cc_current,
        charge_voltage_target_v,
        charge_cv_margin_v,
        charge_cv_min_current_a,
        gamma=2.0,
    )


def _calculate_cv_charge_current(
    va: float,
    vb: float,
    cc_current: float,
    capacity_ah: float,
    charge_voltage_target_v: float,
    charge_cv_margin_v: float,
    charge_cv_min_current_a: float,
) -> float:
    """
    Calculate proportional CV charge current (pack-level control).
    
    Uses target-based linear control so current tapers as V_avg approaches target:
    I_cv = I_cc * ((V_target - V_avg) / V_margin), clamped to [0, I_cc]
    
    Returns minimum _CHARGE_CV_MIN_CURRENT_A to maintain charging phase.
    """
    v_avg = (va + vb) / 2.0
    voltage_error = charge_voltage_target_v - v_avg
    voltage_window = charge_cv_margin_v
    
    if voltage_window <= 0:
        return 0.0
    
    ratio = voltage_error / voltage_window
    if ratio <= 0.0:
        i_cv = 0.0
    else:
        i_cv = abs(cc_current) * min(1.0, ratio)
    
    # Clamp to minimum to prevent oscillation
    i_cv = max(i_cv, charge_cv_min_current_a)
    
    return i_cv


def _evaluate_linear_cccv_policy(
    policy_input: dict[str, Any],
    policy_params: dict[str, Any],
    policy_state: dict[str, Any],
) -> tuple[float, dict[str, Any], dict[str, Any]]:
    """Default linear CCCV policy evaluator."""
    _ = policy_params
    target = _calculate_cv_charge_current(
        policy_input["va"],
        policy_input["vb"],
        policy_input["i_cc_magnitude_a"],
        policy_input["capacity_ah"],
        policy_input["charge_voltage_target_v"],
        policy_input["cv_margin_v"],
        policy_input["min_charge_current_a"],
    )
    telemetry = {
        "policy_mode": "normal",
        "voltage_error_v": policy_input["charge_voltage_target_v"] - policy_input["v_avg"],
        "effective_gain_a_per_v": abs(policy_input["i_cc_magnitude_a"]) / max(1e-12, policy_input["cv_margin_v"]),
    }
    return target, policy_state, telemetry


def _evaluate_power_cccv_policy(
    policy_input: dict[str, Any],
    policy_params: dict[str, Any],
    policy_state: dict[str, Any],
) -> tuple[float, dict[str, Any], dict[str, Any]]:
    """Power-law CCCV policy evaluator."""
    gamma = float(policy_params.get("gamma", 1.5))
    target = _calculate_power_cv_charge_current(
        policy_input["va"],
        policy_input["vb"],
        policy_input["i_cc_magnitude_a"],
        policy_input["charge_voltage_target_v"],
        policy_input["cv_margin_v"],
        policy_input["min_charge_current_a"],
        gamma=gamma,
    )
    telemetry = {
        "policy_mode": "normal",
        "voltage_error_v": policy_input["charge_voltage_target_v"] - policy_input["v_avg"],
        "effective_gain_a_per_v": "nonlinear",
        "gamma": gamma,
    }
    return target, policy_state, telemetry


def _evaluate_quadratic_cccv_policy(
    policy_input: dict[str, Any],
    policy_params: dict[str, Any],
    policy_state: dict[str, Any],
) -> tuple[float, dict[str, Any], dict[str, Any]]:
    """Quadratic CCCV policy evaluator."""
    _ = policy_params
    target = _calculate_quadratic_cv_charge_current(
        policy_input["va"],
        policy_input["vb"],
        policy_input["i_cc_magnitude_a"],
        policy_input["charge_voltage_target_v"],
        policy_input["cv_margin_v"],
        policy_input["min_charge_current_a"],
    )
    telemetry = {
        "policy_mode": "normal",
        "voltage_error_v": policy_input["charge_voltage_target_v"] - policy_input["v_avg"],
        "effective_gain_a_per_v": "nonlinear",
    }
    return target, policy_state, telemetry


CHARGE_POLICY_REGISTRY: dict[str, dict[str, Any]] = {
    "linear_cccv": {
        "policy_id": "linear_cccv",
        "policy_version": "1.0.0",
        "policy_capabilities": {
            "supports_cc": True,
            "supports_cv": True,
            "supports_step_charge": False,
            "supports_discharge": False,
        },
        "required_signals": ["time_s", "dt_s", "va", "vb", "v_avg", "v_max", "i_pack"],
        "optional_signals": ["temperature_c", "plating_risk_score", "diffusion_estimator", "anode_potential_v"],
        "evaluator": _evaluate_linear_cccv_policy,
    },
    "power_cccv": {
        "policy_id": "power_cccv",
        "policy_version": "1.0.0",
        "policy_capabilities": {
            "supports_cc": True,
            "supports_cv": True,
            "supports_step_charge": False,
            "supports_discharge": False,
        },
        "required_signals": ["time_s", "dt_s", "va", "vb", "v_avg", "v_max", "i_pack"],
        "optional_signals": ["temperature_c", "plating_risk_score", "diffusion_estimator", "anode_potential_v"],
        "evaluator": _evaluate_power_cccv_policy,
    },
    "quadratic_cccv": {
        "policy_id": "quadratic_cccv",
        "policy_version": "1.0.0",
        "policy_capabilities": {
            "supports_cc": True,
            "supports_cv": True,
            "supports_step_charge": False,
            "supports_discharge": False,
        },
        "required_signals": ["time_s", "dt_s", "va", "vb", "v_avg", "v_max", "i_pack"],
        "optional_signals": ["temperature_c", "plating_risk_score", "diffusion_estimator", "anode_potential_v"],
        "evaluator": _evaluate_quadratic_cccv_policy,
    },
}


def _get_charge_policy_spec(policy_id: str) -> dict[str, Any]:
    """Resolve policy metadata and evaluator from registry."""
    if policy_id not in CHARGE_POLICY_REGISTRY:
        available = ", ".join(sorted(CHARGE_POLICY_REGISTRY.keys()))
        raise ValueError(f"Unknown charge policy '{policy_id}'. Available: {available}")
    return CHARGE_POLICY_REGISTRY[policy_id]


def _validate_policy_signals(policy_spec: dict[str, Any], available_signals: set[str]) -> tuple[list[str], list[str]]:
    """Return missing required/optional signals for a policy."""
    missing_required = [s for s in policy_spec["required_signals"] if s not in available_signals]
    missing_optional = [s for s in policy_spec["optional_signals"] if s not in available_signals]
    return missing_required, missing_optional


def _build_stage3_policy_input(
    elapsed: float,
    dt: float,
    cycle_idx: int,
    va: float,
    vb: float,
    ia: float,
    ib: float,
    v_avg: float,
    v_max: float,
    i_pack: float,
    cc_current_magnitude: float,
    capacity_ah: float,
    charge_voltage_target_v: float,
    safety_max_cell_v: float,
    cv_margin_v: float,
    min_charge_current_a: float,
) -> dict[str, Any]:
    """Build canonical Stage 3 policy input payload."""
    return {
        "time_s": elapsed,
        "dt_s": dt,
        "cycle_index": cycle_idx + 1,
        "va": va,
        "vb": vb,
        "v_avg": v_avg,
        "v_max": v_max,
        "ia": ia,
        "ib": ib,
        "i_pack": i_pack,
        "i_cc_magnitude_a": cc_current_magnitude,
        "capacity_ah": capacity_ah,
        "charge_voltage_target_v": charge_voltage_target_v,
        "safety_max_cell_v": safety_max_cell_v,
        "cv_margin_v": cv_margin_v,
        "min_charge_current_a": min_charge_current_a,
        "temperature_c": None,
        "plating_risk_score": None,
        "diffusion_estimator": None,
        "anode_potential_v": None,
    }


def _smooth_cv_command(
    target: float,
    cc_current: float,
    cv_state: dict,
    dt: float,
) -> tuple[float, dict]:
    """
    Smooth CV command to prevent bang-bang oscillation (Slice 2).
    
    Uses previous command memory to enforce smooth ramp transitions.
    Prevents sudden jumps between high and low currents.
    
    Aggressive asymmetric ramp limiting:
    - Ramp DOWN: No limit (follow target immediately when decreasing)
    - Ramp UP: 5% per second (prevents overshoot oscillation)
    
    Strategy: Allow current to drop rapidly when proportional law commands low,
    but prevent rapid rise to reduce oscillatory feedback.
    
    Args:
        target: Target command from proportional control law (A)
        cc_current: Full CC current magnitude (A)
        cv_state: State dict with 'prev_command' and 'is_initialized'
        dt: Time step (s)
    
    Returns:
        (command, updated_cv_state): Smoothed command and updated state
    """
    # Initialize on first call
    if not cv_state.get('is_initialized', False):
        cv_state['prev_command'] = target
        cv_state['is_initialized'] = True
        return target, cv_state
    
    # On subsequent calls: smooth the transition with asymmetric ramp
    prev = cv_state['prev_command']
    delta = target - prev
    
    if delta > 0:
        # Ramping UP: slow (5% per second) to prevent overshoot
        max_ramp_per_step = cc_current * 0.05 * (dt / 1.0)
        if abs(delta) <= max_ramp_per_step:
            command = target
        else:
            command = prev + max_ramp_per_step
    else:
        # Ramping DOWN: no limit, follow immediately for taper convergence
        command = target
    
    # Update state for next step
    cv_state['prev_command'] = command
    
    return command, cv_state


def _calculate_per_cell_cv_currents(
    va: float,
    vb: float,
    ia_nominal: float,
    ib_nominal: float,
    v_limit: float = CHARGE_VOLTAGE_TARGET,
    margin: float = _CHARGE_CV_PER_CELL_MARGIN_V,
    per_cell_state: dict = None,
) -> tuple[float, float, dict]:
    """
    Apply per-cell voltage limiting during CV charging.
    
    When a cell exceeds v_limit, reduce its current proportionally.
    Proportional feedback margin: _CHARGE_CV_PER_CELL_MARGIN_V (50mV)
    
    Args:
        va, vb: Cell A and B voltages (V)
        ia_nominal, ib_nominal: Nominal charging currents (A)
        v_limit: Per-cell voltage limit (V)
        per_cell_state: State dict (unused for Slice 2 - kept for compatibility)
    
    Returns:
        (ia_adjusted, ib_adjusted, per_cell_state): Adjusted currents and state
    """
    if per_cell_state is None:
        per_cell_state = {}
    
    ia_adjusted = ia_nominal
    ib_adjusted = ib_nominal
    
    # Per-cell feedback: reduce current if cell exceeds limit
    if va > v_limit:
        # Proportional reduction: reaches 0 current at (v_limit + margin)
        over_voltage = va - v_limit
        if margin > 0:
            reduction_factor = max(0.0, 1.0 - (over_voltage / margin))
        else:
            reduction_factor = 0.0 if over_voltage > 0 else 1.0
        ia_adjusted = ia_nominal * reduction_factor
    
    if vb > v_limit:
        # Same proportional reduction for Cell B
        over_voltage = vb - v_limit
        if margin > 0:
            reduction_factor = max(0.0, 1.0 - (over_voltage / margin))
        else:
            reduction_factor = 0.0 if over_voltage > 0 else 1.0
        ib_adjusted = ib_nominal * reduction_factor
    
    return ia_adjusted, ib_adjusted, per_cell_state


def _check_cv_termination(
    va: float,
    vb: float,
    ia: float,
    ib: float,
    cv_duration_elapsed: float,
    capacity_ah: float,
    v_target: float,
    cv_safety_timeout_s: float = CV_SAFETY_TIMEOUT_S,
    cv_taper_current_threshold_a_fraction: float = CV_TAPER_CURRENT_THRESHOLD_A_FRACTION,
) -> str:
    """
    Check CV charging termination conditions.

    Normal success termination:
    - V_avg >= V_target AND I_avg_abs <= taper_current_threshold

    Safety fallback:
    - cv_elapsed > CV_SAFETY_TIMEOUT_S (logged as cv_safety_timeout)

    Returns termination reason or empty string to continue.
    """
    v_avg = (va + vb) / 2.0
    i_avg_abs = abs((ia + ib) / 2.0)
    taper_current_threshold = cv_taper_current_threshold_a_fraction * capacity_ah
    
    # Voltage must reach target first
    if v_avg >= v_target and i_avg_abs <= taper_current_threshold:
        return (
            "capacity_limit: cv_taper_current: "
            f"V_avg={v_avg:.3f}V >= {v_target:.2f}V "
            f"AND I_avg_abs={i_avg_abs:.4f}A <= {taper_current_threshold:.4f}A"
        )

    if cv_duration_elapsed > cv_safety_timeout_s:
        return (
            "time_limit: cv_safety_timeout: "
            f"cv_elapsed={cv_duration_elapsed:.0f}s > {cv_safety_timeout_s:.0f}s "
            f"(V_avg={v_avg:.3f}V, I_avg_abs={i_avg_abs:.4f}A)"
        )
    
    return ""


def _cycle_end_reason(
    is_discharge: bool,
    soc_a: float,
    soc_b: float,
    soc_lower_bound: float,
    soc_upper_bound: float,
    va: float = float("nan"),
    vb: float = float("nan"),
    discharge_voltage_only: bool = False,
    charge_voltage_only: bool = False,
    discharge_voltage_cutoff_v: float = _DISCHARGE_VOLTAGE_CUTOFF_V,
    charge_voltage_limit_v: float = 4.20,
) -> str:
    """Return cycle end reason, or empty string when cycle should continue.
    
    Added: voltage-limit-aware charge termination to prevent continued CC charge
    after cells reach the configured max voltage (PyBaMM DFN max voltage event surface).
    """
    import math
    if is_discharge:
        if not math.isnan(va) and not math.isnan(vb):
            if va <= discharge_voltage_cutoff_v or vb <= discharge_voltage_cutoff_v:
                v_min = min(va, vb)
                return (
                    "min_voltage_reached: Discharge end: at least one cell reached voltage cutoff "
                    f"({discharge_voltage_cutoff_v:.2f} V, min={v_min:.4f}V)"
                )
        if (not discharge_voltage_only) and (soc_a <= soc_lower_bound or soc_b <= soc_lower_bound):
            return (
                "capacity_limit: Discharge end: at least one cell reached lower SOC bound "
                f"({soc_lower_bound:.1f}%)"
            )
    if not is_discharge:
        # Charge: check voltage limit first (new behavior)
        if not math.isnan(va) and not math.isnan(vb):
            if va >= charge_voltage_limit_v or vb >= charge_voltage_limit_v:
                v_max = max(va, vb)
                return (
                    "max_voltage_reached: Charge end: at least one cell reached voltage limit "
                    f"({charge_voltage_limit_v:.2f} V, max={v_max:.4f}V)"
                )
        if (not charge_voltage_only) and (soc_a >= soc_upper_bound or soc_b >= soc_upper_bound):
            return (
                "capacity_limit: Charge end: at least one cell reached upper SOC bound "
                f"({soc_upper_bound:.1f}%)"
            )
    return ""


def _get_current_soc(pack_block) -> tuple[float, float]:
    """Read latest SOC values from block outputs."""
    soc_a = float(pack_block.outputs[2])
    soc_b = float(pack_block.outputs[3])
    return soc_a, soc_b


def _current_for_cycle(base_current: float, cycle_idx: int, alternate: bool) -> float:
    """Generate signed pack current for each cycle."""
    if not alternate:
        return float(base_current)
    if cycle_idx % 2 == 0:
        return float(base_current)
    return -float(base_current)


def _print_cycle_debug(cycle_idx: int, phase: str, soc_a: float, soc_b: float, direction: str) -> None:
    """Print cycle boundary SOC to verify state continuity across cycles."""
    print(
        f"[Cycle {cycle_idx + 1} {direction}] {phase}: "
        f"SOC_A={soc_a:.3f}%, SOC_B={soc_b:.3f}%"
    )


def _run_rest_phase(
    sim,
    current_source,
    parallel_pack,
    dt: float,
    rest_duration_s: float,
    elapsed: float,
    result: SimulationResult,
    cycle_idx: int,
) -> float:
    """Execute rest phase at 0A for specified duration."""
    rest_steps = int(rest_duration_s / dt)
    current_source.current = 0.0
    print(f"[Cycle {cycle_idx + 1}] Rest phase: {rest_duration_s:.0f}s ({rest_steps} steps)")
    
    for _ in range(rest_steps):
        elapsed = _record_step(sim, parallel_pack, dt, elapsed, result, cycle_index=cycle_idx + 1)
    
    print(f"[Cycle {cycle_idx + 1}] Rest phase complete")
    return elapsed


def run_simulation(
    blocks,
    connections,
    current_source,
    parallel_pack,
    dt: float,
    simulation_mode: str,
    sim_duration: float | None,
    max_cycle_time: float | None,
    pack_current: float,
    cycle_count: int = 1,
    alternate_charge_discharge: bool = False,
    soc_lower_bound: float = 0.0,
    soc_upper_bound: float = 100.0,
    charge_policy_id: str = "linear_cccv",
    charge_policy_params: dict[str, Any] | None = None,
    discharge_voltage_only: bool = False,
    charge_voltage_only: bool = False,
    rest_duration_s: float = 1800.0,
    sim_config: SimulationConfig | None = None,
) -> SimulationResult:
    """
    Execute pack simulation for fixed-time or cycle mode.

        Supports:
        - charge/discharge full-cycle termination
        - multiple event loop with alternating charge/discharge
        - termination conditions per event

        Semantics:
        - alternate_charge_discharge=False: cycle_count is event count
        - alternate_charge_discharge=True: cycle_count is full-cycle count,
            executed internally as 2 * cycle_count events
    """
    result = SimulationResult()
    runtime_config = SimulationConfig() if sim_config is None else sim_config

    sim = Simulation(blocks, connections=connections, dt=dt)
    sim.reset(0.0)

    elapsed = 0.0

    if simulation_mode == "Fixed Time Mode":
        if sim_duration is None:
            raise ValueError("Fixed Time Mode requires sim_duration.")
        n_steps = int(float(sim_duration) / float(dt))
        for _ in range(n_steps):
            elapsed = _record_step(sim, parallel_pack, dt, elapsed, result, cycle_index=1)
        result.stop_reason = "Fixed duration reached"
        result.completed_cycles = 1
        result.full_cycle_count_requested = 0
        result.event_count_requested = 0
        result.event_count_executed = 0
        result.full_cycles_completed = 0
        result.cycle_summaries = _build_cycle_summaries(
            result,
            pack_current=pack_current,
            alternate_charge_discharge=alternate_charge_discharge,
            nominal_capacity_a=parallel_pack.cell_a.capacity_ah,
            nominal_capacity_b=parallel_pack.cell_b.capacity_ah,
        )
        result.aging_summary = _build_aging_summary(result.cycle_summaries)
        result.aging_artifacts = _write_aging_artifacts(result.cycle_summaries, result.aging_summary)
        return result

    if abs(float(pack_current)) < 1e-12:
        raise ValueError("Full Cycle Mode requires non-zero pack current.")
    if max_cycle_time is None:
        raise ValueError("Full Cycle Mode requires max_cycle_time.")

    soc_lower_bound = max(0.0, min(100.0, float(soc_lower_bound)))
    soc_upper_bound = max(0.0, min(100.0, float(soc_upper_bound)))
    if soc_lower_bound >= soc_upper_bound:
        raise ValueError("SOC bounds must satisfy lower < upper.")

    max_steps = int(float(max_cycle_time) / float(dt))
    full_cycle_count_requested = max(1, int(cycle_count))
    if alternate_charge_discharge:
        event_count_requested = full_cycle_count_requested * 2
    else:
        event_count_requested = full_cycle_count_requested

    result.full_cycle_count_requested = full_cycle_count_requested
    result.event_count_requested = event_count_requested
    charge_policy_params = {} if charge_policy_params is None else dict(charge_policy_params)
    policy_spec = _get_charge_policy_spec(charge_policy_id)
    available_signals = {
        "time_s", "dt_s", "cycle_index",
        "va", "vb", "v_avg", "v_max", "ia", "ib", "i_pack",
        "i_cc_magnitude_a", "capacity_ah", "charge_voltage_target_v",
        "safety_max_cell_v", "cv_margin_v", "min_charge_current_a",
    }
    missing_required, missing_optional = _validate_policy_signals(policy_spec, available_signals)
    if missing_required:
        raise ValueError(
            f"Charge policy '{charge_policy_id}' missing required signals: {missing_required}"
        )
    if missing_optional:
        print(
            f"[Stage3Policy] optional signals unavailable for '{charge_policy_id}': {missing_optional}"
        )
    print(
        f"[Stage3Policy] selected='{policy_spec['policy_id']}' "
        f"version={policy_spec['policy_version']}"
    )

    for cycle_idx in range(event_count_requested):
        cycle_current = _current_for_cycle(pack_current, cycle_idx, alternate_charge_discharge)
        current_source.current = float(cycle_current)
        is_discharge = cycle_current > 0.0
        direction = "Discharge" if is_discharge else "Charge"

        # If a cycle starts exactly at the boundary, mark it complete and continue.
        start_soc_a, start_soc_b = _get_current_soc(parallel_pack)
        _print_cycle_debug(cycle_idx, "Cycle Start", start_soc_a, start_soc_b, direction)
        start_reason = _cycle_end_reason(
            is_discharge,
            start_soc_a,
            start_soc_b,
            soc_lower_bound,
            soc_upper_bound,
            discharge_voltage_only=discharge_voltage_only,
            charge_voltage_only=charge_voltage_only,
            discharge_voltage_cutoff_v=runtime_config.discharge_voltage_limit_v,
            charge_voltage_limit_v=runtime_config.charge_voltage_limit_v,
        )
        if start_reason:
            result.cycle_reasons.append(
                f"Cycle {cycle_idx + 1} ({direction}): start boundary reached ({start_reason})"
            )
            result.completed_cycles = cycle_idx + 1
            result.event_count_executed = result.completed_cycles
            result.full_cycles_completed = (
                result.completed_cycles // 2 if alternate_charge_discharge else result.completed_cycles
            )
            result.stop_reason = start_reason
            continue

        cycle_reason = ""
        
        # ===== CC/CV state tracking for charge cycles =====
        cv_mode_active = False
        cv_mode_duration = 0.0
        cc_current_magnitude = abs(cycle_current)
        cv_command_state = {}  # Slice 2: State for smooth CV command
        per_cell_state = {}  # Slice 2: State for per-cell limiting hysteresis
        policy_state: dict[str, Any] = {}
        
        for _ in range(max_steps):
            elapsed = _record_step(sim, parallel_pack, dt, elapsed, result, cycle_index=cycle_idx + 1)

            # ===== Stage 1: Read State =====
            # Gather all voltages, currents, SOC, and derived values
            soc_a = result.soc_a_trace[-1]
            soc_b = result.soc_b_trace[-1]
            va = result.va_trace[-1]
            vb = result.vb_trace[-1]
            ia = result.ia_trace[-1]
            ib = result.ib_trace[-1]
            v_avg = (va + vb) / 2.0
            v_max = max(va, vb)
            i_avg_abs = abs((ia + ib) / 2.0)
            i_pack = ia + ib

            # ===== Stage 2: Phase Evaluation =====
            # Determine CC/CV state and update CV latch
            if not is_discharge:  # Charge cycle only
                # CC->CV transition: target-based per-cell entry
                if cv_mode_active or _at_or_above_with_tolerance(
                    v_max,
                    runtime_config.cv_entry_v_per_cell,
                    runtime_config.cv_entry_abs_tol_v,
                ):
                    if not cv_mode_active:
                        cv_mode_active = True
                        cv_mode_duration = 0.0
                        print(
                            f"[Cycle {cycle_idx + 1}] CV entry t={elapsed:.1f}s "
                            f"V_avg={v_avg:.4f}V V_max={v_max:.4f}V "
                            f"I_pack={i_pack:.4f}A target={runtime_config.charge_voltage_target_v:.4f}V"
                        )
                    cv_mode_duration += dt
                else:
                    cv_mode_duration = 0.0

            # ===== Stage 3: Command Generation =====
            # Generate current command based on phase
            if not is_discharge:  # Charge cycle
                if cv_mode_active:
                    # CV mode: policy-driven target generation + existing smoothing chain
                    policy_input = _build_stage3_policy_input(
                        elapsed=elapsed,
                        dt=dt,
                        cycle_idx=cycle_idx,
                        va=va,
                        vb=vb,
                        ia=ia,
                        ib=ib,
                        v_avg=v_avg,
                        v_max=v_max,
                        i_pack=i_pack,
                        cc_current_magnitude=cc_current_magnitude,
                        capacity_ah=parallel_pack.cell_a.capacity_ah,
                        charge_voltage_target_v=runtime_config.charge_voltage_target_v,
                        safety_max_cell_v=runtime_config.safety_max_cell_v,
                        cv_margin_v=runtime_config.charge_cv_margin_v,
                        min_charge_current_a=runtime_config.charge_cv_min_current_a,
                    )
                    target, policy_state, _policy_telemetry = policy_spec["evaluator"](
                        policy_input,
                        charge_policy_params,
                        policy_state,
                    )
                    i_cv_pack, cv_command_state = _smooth_cv_command(target, cc_current_magnitude, cv_command_state, dt)
                    
                    # Calculate nominal per-cell currents using current divider
                    r_a = parallel_pack.cell_a.get_resistance()
                    r_b = parallel_pack.cell_b.get_resistance()
                    i_a_nominal = i_cv_pack * r_b / (r_a + r_b)
                    i_b_nominal = i_cv_pack * r_a / (r_a + r_b)
                    
                    # Apply per-cell feedback limiting
                    i_a_adjusted, i_b_adjusted, per_cell_state = _calculate_per_cell_cv_currents(
                        va,
                        vb,
                        i_a_nominal,
                        i_b_nominal,
                        runtime_config.charge_voltage_target_v,
                        runtime_config.charge_cv_margin_v,
                        per_cell_state,
                    )
                    
                    # Total pack current (sum of adjusted per-cell currents)
                    i_cv_adjusted = i_a_adjusted + i_b_adjusted
                    current_source.current = -i_cv_adjusted  # Negative for charging
                else:
                    # CC mode: constant current
                    current_source.current = float(cycle_current)

            # ===== Stage 4: Termination Candidate Evaluation =====
            # Check all termination conditions and build candidates
            cycle_reason = ""
            termination_logged = False

            # Candidate 0: wrapper error detected in Sequence A block update
            if bool(getattr(parallel_pack, "last_wrapper_err", False)):
                cycle_reason = "wrapper_error: Sequence A pack wrapper reported error"
                termination_logged = True
            
            # Candidate 1: Generic cycle boundary (discharge SOC/voltage, charge SOC/voltage)
            # Added: charge voltage limit to prevent continued CC after max voltage event
            if not cycle_reason:
                cycle_reason = _cycle_end_reason(
                    is_discharge,
                    soc_a,
                    soc_b,
                    soc_lower_bound,
                    soc_upper_bound,
                    va,
                    vb,
                    discharge_voltage_only=discharge_voltage_only,
                    charge_voltage_only=charge_voltage_only,
                    discharge_voltage_cutoff_v=runtime_config.discharge_voltage_limit_v,
                    charge_voltage_limit_v=runtime_config.charge_voltage_limit_v,
                )
            
            # Candidate 2: CV taper termination (charge only, when in CV)
            if not is_discharge and cv_mode_active:
                cv_termination_reason = _check_cv_termination(
                    va,
                    vb,
                    ia,
                    ib,
                    cv_mode_duration,
                    parallel_pack.cell_a.capacity_ah,
                    runtime_config.charge_voltage_target_v,
                    runtime_config.cv_safety_timeout_s,
                    runtime_config.cv_taper_current_threshold_a_fraction,
                )
                if cv_termination_reason:
                    cycle_reason = cv_termination_reason
                    print(
                        f"[Cycle {cycle_idx + 1}] Charge termination t={elapsed:.1f}s "
                        f"reason={cv_termination_reason} "
                        f"I_cell_avg_abs={i_avg_abs:.4f}A I_pack={i_pack:.4f}A"
                    )
                    termination_logged = True

            # Candidate 3: Safety termination (charge only, any phase)
            if not is_discharge and v_max >= runtime_config.safety_max_cell_v:
                cycle_reason = (
                    "max_voltage_reached: charge_safety_max_cell_v: "
                    f"V_max={v_max:.4f}V >= {runtime_config.safety_max_cell_v:.4f}V"
                )
                print(
                    f"[Cycle {cycle_idx + 1}] Charge termination t={elapsed:.1f}s "
                    f"reason={cycle_reason} I_pack={i_pack:.4f}A"
                )
                termination_logged = True

            # ===== Stage 5: Arbitration =====
            # Apply winning termination reason (already selected above via overwrite precedence)
            if cycle_reason and (not is_discharge) and (not termination_logged):
                print(
                    f"[Cycle {cycle_idx + 1}] Charge termination t={elapsed:.1f}s "
                    f"reason={cycle_reason} I_cell_avg_abs={i_avg_abs:.4f}A I_pack={i_pack:.4f}A"
                )
            
            # ===== Stage 6: Apply Result =====
            # Break cycle if termination reason exists, otherwise continue to next step
            if cycle_reason:
                break

        if not cycle_reason:
            cycle_reason = "time_limit: Max cycle time reached"

        end_soc_a, end_soc_b = _get_current_soc(parallel_pack)
        _print_cycle_debug(cycle_idx, "Cycle End", end_soc_a, end_soc_b, direction)
        result.cycle_reasons.append(f"Cycle {cycle_idx + 1} ({direction}): {cycle_reason}")

        # ===== Rest Phase =====
        # Insert rest period after each cycle (except on time_limit early exit)
        if cycle_reason != "time_limit: Max cycle time reached" and rest_duration_s > 0.0:
            if cycle_idx < event_count_requested - 1:
                elapsed = _run_rest_phase(sim, current_source, parallel_pack, dt, rest_duration_s, elapsed, result, cycle_idx)

        if cycle_reason == "time_limit: Max cycle time reached":
            result.completed_cycles = cycle_idx + 1
            result.event_count_executed = result.completed_cycles
            result.full_cycles_completed = (
                result.completed_cycles // 2 if alternate_charge_discharge else result.completed_cycles
            )
            result.stop_reason = cycle_reason
            break

        # B2 lifecycle transition: rebuild DFN simulations from final per-cell solutions.
        # This prevents reusing event-terminated solutions directly via step(...) in next cycle.
        if cycle_idx < event_count_requested - 1:
            for cell_model in (parallel_pack.cell_a, parallel_pack.cell_b):
                if hasattr(cell_model, "reinitialize_from_final_solution"):
                    try:
                        cell_model.reinitialize_from_final_solution()
                    except Exception as exc:
                        print(f"[SimulationEngine] lifecycle transition warning: {exc}")

        result.completed_cycles = cycle_idx + 1
        result.event_count_executed = result.completed_cycles
        result.full_cycles_completed = (
            result.completed_cycles // 2 if alternate_charge_discharge else result.completed_cycles
        )
        result.stop_reason = cycle_reason

    if event_count_requested > 1:
        if alternate_charge_discharge:
            result.stop_reason = (
                f"{result.full_cycles_completed}/{result.full_cycle_count_requested} full cycles completed "
                f"({result.event_count_executed}/{result.event_count_requested} events). "
                f"Last: {result.stop_reason}"
            )
        else:
            result.stop_reason = (
                f"{result.completed_cycles}/{event_count_requested} events completed. "
                f"Last: {result.stop_reason}"
            )

    result.cycle_summaries = _build_cycle_summaries(
        result,
        pack_current=pack_current,
        alternate_charge_discharge=alternate_charge_discharge,
        nominal_capacity_a=parallel_pack.cell_a.capacity_ah,
        nominal_capacity_b=parallel_pack.cell_b.capacity_ah,
    )
    result.aging_summary = _build_aging_summary(result.cycle_summaries)
    result.aging_artifacts = _write_aging_artifacts(result.cycle_summaries, result.aging_summary)

    return result