Analysis routines

Higher-level analysis workflows that orchestrate repeated run_fault calls and answer parametric questions about an existing network — currently the inverse rho problem at the bus-grounding side.

Physical / modelling context

A fully built Network solves a forward problem: given the bus grounding impedances \(Z_{\text{B},i}(\rho_E, f)\), find the EPR at the fault bus. The natural inverse question — "how poor can the soil get before the EPR exceeds a touch-voltage limit \(u_{\max}\)?" — is recurring in safety-engineering workflows: it sets the worst-case specific earth resistance the network can tolerate while still satisfying the relevant standard (e.g. EN 50522) at the assumed fault current.

find_max_rho_scaling solves that problem by log-bisection over a uniform scaling factor \(c\) of \(\rho_E\) on a user-supplied list of selected buses. At every trial \(c\), the bus-grounding impedance is re-evaluated through the bus's own BusType.impedance_formula and run_fault is invoked. The bracket converges on the largest \(c\) for which \(|U_\text{EPR}(f)|_{\text{RMS}} \le u_{\max}\) at the fault bus. The sister routine find_max_rho_f_scaling (and the diagnostic helper evaluate_max_epr_under_k) extend the same idea to a frequency- dependent rho-f characteristic: rho is scaled while a separate factor \(k\) tunes the imaginary, frequency-coupling part of the formula, so the inverse problem can be posed against a parametric rho-f curve rather than a single scalar.

The shape of the rho-f curve at each bus is controlled entirely by the existing BusType.impedance_formula; the analysis routines only vary the scalar factors, so any user-defined parametric form (linear, square root, frequency-dependent, …) is supported transparently. Original rho values are restored via a finally block, so the network is left untouched even if the bisection fails.

Example

import groundinsight as gi

# Assume `net` is a built network with at least one fault and one
# source. Pick the buses whose specific_earth_resistance should be
# scaled jointly — typically all buses sharing a soil environment.
result = gi.find_max_rho_scaling(
    network=net,
    fault_name="fault1",
    bus_names=["bus_substation", "bus_fault"],
    u_max=200.0,            # touch-voltage limit in volts (RMS)
    c_bounds=(0.1, 100.0),  # search bracket on the scaling factor
    tol_rel=1e-3,
    max_iter=40,
)

# `result` is a dict with keys: c_max, u_epr_rms_at_c_max,
# rho_max (per-bus dict), iterations.
print(f"c_max = {result['c_max']:.3f}")
print(f"EPR at c_max = {result['u_epr_rms_at_c_max']:.1f} V")
for bus, rho_max in result["rho_max"].items():
    print(f"  {bus}: rho_max = {rho_max:.0f} Ω·m")
print(f"converged in {result['iterations']} iterations")

For the rho-f variant — useful when the bus impedance carries a frequency-coupling term whose magnitude should also be probed — see :func:find_max_rho_f_scaling.

API reference

analysis

Analysis subpackage.

Contains higher-level analysis routines on top of :mod:network_operations, such as inverse problems for the bus-grounding rho-f characteristic. These functions assume that the network is already fully built (buses, branches, faults, sources, paths) and orchestrate repeated :func:run_fault calls.

evaluate_max_epr_under_k

evaluate_max_epr_under_k(
    network: Network,
    bus_names: List[str],
    k: KVector,
    *,
    fault_scalings: Optional[Dict[float, float]] = None,
    run_fault_kwargs: Optional[Dict[str, Any]] = None
) -> Dict[str, float]

Evaluate the RMS EPR at each selected bus for a given k vector.

For every name in bus_names the bus grounding impedance is overwritten with the rho-f linear form Z(rho, f) = k1*rho + (k2+jk3)*f + (k4+jk5)*rho*f -- using the bus' own :attr:Bus.specific_earth_resistance for rho -- and the bus is swept as the active fault: an existing fault at that bus is reused, otherwise a temporary one is created with fault_scalings (or 1.0 at every simulation frequency by default) and removed at the end. The original bus impedances, the temporary faults and the previous active_fault are restored in a finally block, so the network is left exactly as it was on entry.

Parameters:

Name	Type	Description	Default
network	Network	The simulation network. Sources, branches, paths and the buses to be swept must already be configured.	required
bus_names	List[str]	Buses whose impedance is rewritten and which are swept as fault locations one by one. Must be non-empty and refer to buses in network.	required
k	KVector	5-tuple (k1, k2, k3, k4, k5) of real model parameters.	required
fault_scalings	Optional[Dict[float, float]]	Frequency-resolved scalings used for any fault that has to be created on the fly (no pre-existing fault at the swept bus). Defaults to {f: 1.0} for every simulation frequency. Existing faults are reused unmodified.	None
run_fault_kwargs	Optional[Dict[str, Any]]	Extra keyword arguments forwarded to :func:run_fault at every step.	None

Returns:

Type	Description
Dict[str, float]	Dict mapping each bus_name to its RMS EPR in volts at that
Dict[str, float]	bus when it is the active fault.

Raises:

Type	Description
ValueError	If bus_names is empty, k does not have length 5, or any name is not in the network.

Source code in src/groundinsight/analysis/inverse_rho_f.py

def evaluate_max_epr_under_k(
    network: Network,
    bus_names: List[str],
    k: KVector,
    *,
    fault_scalings: Optional[Dict[float, float]] = None,
    run_fault_kwargs: Optional[Dict[str, Any]] = None,
) -> Dict[str, float]:
    """Evaluate the RMS EPR at each selected bus for a given ``k`` vector.

    For every name in ``bus_names`` the bus grounding impedance is
    overwritten with the rho-f linear form ``Z(rho, f) = k1*rho +
    (k2+jk3)*f + (k4+jk5)*rho*f`` -- using the bus' own
    :attr:`Bus.specific_earth_resistance` for ``rho`` -- and the bus is
    swept as the active fault: an existing fault at that bus is reused,
    otherwise a temporary one is created with ``fault_scalings`` (or
    ``1.0`` at every simulation frequency by default) and removed at the
    end. The original bus impedances, the temporary faults and the
    previous ``active_fault`` are restored in a ``finally`` block, so the
    network is left exactly as it was on entry.

    Args:
        network: The simulation network. Sources, branches, paths and
            the buses to be swept must already be configured.
        bus_names: Buses whose impedance is rewritten and which are
            swept as fault locations one by one. Must be non-empty and
            refer to buses in ``network``.
        k: 5-tuple ``(k1, k2, k3, k4, k5)`` of real model parameters.
        fault_scalings: Frequency-resolved scalings used for any fault
            that has to be created on the fly (no pre-existing fault at
            the swept bus). Defaults to ``{f: 1.0}`` for every simulation
            frequency. Existing faults are reused unmodified.
        run_fault_kwargs: Extra keyword arguments forwarded to
            :func:`run_fault` at every step.

    Returns:
        Dict mapping each ``bus_name`` to its RMS EPR in volts at that
        bus when it is the active fault.

    Raises:
        ValueError: If ``bus_names`` is empty, ``k`` does not have
            length 5, or any name is not in the network.
    """
    if not bus_names:
        raise ValueError("bus_names must not be empty.")
    if len(k) != 5:
        raise ValueError(f"k must be a 5-tuple, got length {len(k)}.")
    missing = [b for b in bus_names if b not in network.buses]
    if missing:
        raise ValueError(f"Unknown bus(es) in network: {missing!r}.")

    rfk: Dict[str, Any] = dict(run_fault_kwargs) if run_fault_kwargs else {}
    if fault_scalings is None:
        fault_scalings = {float(f): 1.0 for f in network.frequencies}

    # Snapshot current bus impedances so the network can be restored.
    impedance_backup: Dict[str, Dict[float, ComplexNumber]] = {
        b: dict(network.buses[b].impedance) for b in bus_names
    }

    # Look up an existing fault per swept bus (if any).
    existing_faults_by_bus: Dict[str, str] = {}
    for fname, fault in network.faults.items():
        if fault.bus in bus_names and fault.bus not in existing_faults_by_bus:
            existing_faults_by_bus[fault.bus] = fname

    temp_faults_created: List[str] = []
    active_fault_backup = network.active_fault

    # Snapshot existing paths and clear them: ``run_fault`` only triggers
    # ``define_paths`` when ``network.paths`` is empty. We need fresh paths
    # for every (source, fault) combination -- including the temporary
    # faults we are about to create -- so we drop the cache for the
    # duration of the sweep and restore it at the end.
    paths_backup = dict(network.paths)
    network.paths.clear()

    def _temp_fault_name(b: str) -> str:
        base = f"_inv_rhof_{b}"
        candidate = base
        n = 0
        while candidate in network.faults:
            n += 1
            candidate = f"{base}_{n}"
        return candidate

    epr_per_bus: Dict[str, float] = {}

    try:
        # Overwrite every selected bus impedance with the k-form.
        for b in bus_names:
            bus = network.buses[b]
            rho = float(bus.specific_earth_resistance)
            new_imp: Dict[float, ComplexNumber] = {}
            for f in network.frequencies:
                z = _z_rho_f(k, rho, float(f))
                new_imp[float(f)] = ComplexNumber(real=z.real, imag=z.imag)
            bus.impedance = new_imp

        # Pre-create any temporary faults the sweep needs *before* the
        # first ``run_fault`` so ``define_paths`` enumerates paths for
        # all swept (source, fault) combinations in a single pass.
        sweep_fault_names: List[str] = []
        for b in bus_names:
            if b in existing_faults_by_bus:
                sweep_fault_names.append(existing_faults_by_bus[b])
            else:
                fname = _temp_fault_name(b)
                network.add_fault(
                    Fault(name=fname, bus=b, scalings=dict(fault_scalings))
                )
                temp_faults_created.append(fname)
                sweep_fault_names.append(fname)

        # Sweep: each fault once. The first call populates ``network.paths``
        # for every (source, fault) pair; subsequent calls reuse the cache.
        for b, fname in zip(bus_names, sweep_fault_names):
            run_fault(network, fault_name=fname, **rfk)
            result_bus = next(
                rb for rb in network.results[fname].buses if rb.name == b
            )
            epr_per_bus[b] = float(result_bus.uepr)
    finally:
        # Restore bus impedances regardless of success/failure.
        for b in bus_names:
            network.buses[b].impedance = impedance_backup[b]
        # Drop temporary faults and any results they produced.
        for fname in temp_faults_created:
            network.faults.pop(fname, None)
            network.results.pop(fname, None)
        # Restore the path cache (cleared at the start of the sweep).
        network.paths.clear()
        network.paths.update(paths_backup)
        # Restore active_fault best-effort.
        if active_fault_backup is not None:
            try:
                network.active_fault = active_fault_backup
            except Exception:  # pragma: no cover -- defensive
                logger.warning(
                    "Could not restore active_fault to %r.", active_fault_backup
                )

    return epr_per_bus

find_max_rho_f_scaling

find_max_rho_f_scaling(
    network: Network,
    bus_names: List[str],
    u_limit: float,
    k_ref: KVector,
    *,
    c_bounds: Tuple[float, float] = (0.001, 1000.0),
    tol_rel: float = 0.001,
    max_iter: int = 60,
    fault_scalings: Optional[Dict[float, float]] = None,
    run_fault_kwargs: Optional[Dict[str, Any]] = None
) -> Dict[str, Any]

Find the largest scaling factor c of a reference k-vector.

Performs a log-scale bisection on c so that k = c * k_ref yields a maximum RMS EPR (across all buses swept as faults) that is at or below u_limit. k_ref is typically a fit produced by :func:groundmeas.services.analytics.rho_f_model from measured impedance points; this function answers how much head-room the network has around that characteristic.

Parameters:

Name	Type	Description	Default
network	Network	The simulation network. Buses, branches, sources and paths must already be configured.	required
bus_names	List[str]	Buses whose impedance is parameterised by k and which are swept as fault locations.	required
u_limit	float	Upper bound on the RMS EPR (in volts) at any swept bus. Must be strictly positive.	required
k_ref	KVector	Reference 5-tuple (k1_ref, ..., k5_ref). Must have length 5 and not be the zero vector.	required
c_bounds	Tuple[float, float]	Search interval for c. Strictly positive, c_bounds[0] < c_bounds[1]. Defaults to (1e-3, 1e3).	(0.001, 1000.0)
tol_rel	float	Bisection tolerance on the relative bracket width (c_hi - c_lo) / c_lo. Defaults to 1e-3.	0.001
max_iter	int	Hard cap on bisection steps. Defaults to 60.	60
fault_scalings	Optional[Dict[float, float]]	See :func:evaluate_max_epr_under_k.	None
run_fault_kwargs	Optional[Dict[str, Any]]	Forwarded to :func:run_fault.	None

Returns:

Type	Description
Dict[str, Any]	Dict with keys
Dict[str, Any]	"c_max" (float): Largest scaling factor compatible with u_limit within the bracket.
Dict[str, Any]	"k_max" (Tuple[float, ...]): c_max * k_ref.
Dict[str, Any]	"max_epr_rms_at_c_max" (float): Maximum RMS EPR across all swept buses at c_max, in volts.
Dict[str, Any]	"epr_rms_per_bus_at_c_max" (Dict[str, float]): RMS EPR per swept bus at c_max.
Dict[str, Any]	"iterations" (int): Number of bisection steps taken (0 if the upper bracket bound is already admissible).

Raises:

Type	Description
ValueError	For invalid input or if c = c_bounds[0] already violates the EPR limit.

Examples:

>>> import groundinsight as gi
>>> from groundinsight.models.core_models import BusType, BranchType
>>> from groundinsight.analysis import find_max_rho_f_scaling
>>> bt = BusType(name="BT", system_type="Grounded",
...              voltage_level=20.0,
...              impedance_formula="rho * 0.01 + 0*f")
>>> brt = BranchType(name="BRT", grounding_conductor=True,
...                  self_impedance_formula="(0.25 + I*0.6)*l",
...                  mutual_impedance_formula="(0.0 + I*0.6)*l")
>>> net = gi.create_network(name="N", frequencies=[50])
>>> _ = gi.create_bus(name="b0", type=bt, network=net)
>>> _ = gi.create_bus(name="b1", type=bt, network=net)
>>> _ = gi.create_branch(name="br", type=brt, from_bus="b0",
...                      to_bus="b1", length=1.0, network=net)
>>> _ = gi.create_source(name="src", bus="b0",
...                      values={50: 100.0}, network=net)
>>> # k_ref reproducing the original Z = 0.01*rho on b0/b1.
>>> res = find_max_rho_f_scaling(
...     net, ["b0", "b1"], u_limit=200.0,
...     k_ref=(0.01, 0.0, 0.0, 0.0, 0.0),
... )
>>> isinstance(res["c_max"], float)
True

Source code in src/groundinsight/analysis/inverse_rho_f.py

def find_max_rho_f_scaling(
    network: Network,
    bus_names: List[str],
    u_limit: float,
    k_ref: KVector,
    *,
    c_bounds: Tuple[float, float] = (1e-3, 1e3),
    tol_rel: float = 1e-3,
    max_iter: int = 60,
    fault_scalings: Optional[Dict[float, float]] = None,
    run_fault_kwargs: Optional[Dict[str, Any]] = None,
) -> Dict[str, Any]:
    """Find the largest scaling factor ``c`` of a reference k-vector.

    Performs a log-scale bisection on ``c`` so that ``k = c * k_ref``
    yields a maximum RMS EPR (across all buses swept as faults) that is
    at or below ``u_limit``. ``k_ref`` is typically a fit produced by
    :func:`groundmeas.services.analytics.rho_f_model` from measured
    impedance points; this function answers *how much head-room the
    network has around that characteristic*.

    Args:
        network: The simulation network. Buses, branches, sources and
            paths must already be configured.
        bus_names: Buses whose impedance is parameterised by ``k`` *and*
            which are swept as fault locations.
        u_limit: Upper bound on the RMS EPR (in volts) at any swept bus.
            Must be strictly positive.
        k_ref: Reference 5-tuple ``(k1_ref, ..., k5_ref)``. Must have
            length 5 and not be the zero vector.
        c_bounds: Search interval for ``c``. Strictly positive,
            ``c_bounds[0] < c_bounds[1]``. Defaults to ``(1e-3, 1e3)``.
        tol_rel: Bisection tolerance on the relative bracket width
            ``(c_hi - c_lo) / c_lo``. Defaults to ``1e-3``.
        max_iter: Hard cap on bisection steps. Defaults to ``60``.
        fault_scalings: See :func:`evaluate_max_epr_under_k`.
        run_fault_kwargs: Forwarded to :func:`run_fault`.

    Returns:
        Dict with keys

        - ``"c_max"`` (float): Largest scaling factor compatible with
          ``u_limit`` within the bracket.
        - ``"k_max"`` (Tuple[float, ...]): ``c_max * k_ref``.
        - ``"max_epr_rms_at_c_max"`` (float): Maximum RMS EPR across all
          swept buses at ``c_max``, in volts.
        - ``"epr_rms_per_bus_at_c_max"`` (Dict[str, float]): RMS EPR
          per swept bus at ``c_max``.
        - ``"iterations"`` (int): Number of bisection steps taken
          (``0`` if the upper bracket bound is already admissible).

    Raises:
        ValueError: For invalid input or if ``c = c_bounds[0]`` already
            violates the EPR limit.

    Examples:
        >>> import groundinsight as gi
        >>> from groundinsight.models.core_models import BusType, BranchType
        >>> from groundinsight.analysis import find_max_rho_f_scaling
        >>> bt = BusType(name="BT", system_type="Grounded",
        ...              voltage_level=20.0,
        ...              impedance_formula="rho * 0.01 + 0*f")
        >>> brt = BranchType(name="BRT", grounding_conductor=True,
        ...                  self_impedance_formula="(0.25 + I*0.6)*l",
        ...                  mutual_impedance_formula="(0.0 + I*0.6)*l")
        >>> net = gi.create_network(name="N", frequencies=[50])
        >>> _ = gi.create_bus(name="b0", type=bt, network=net)
        >>> _ = gi.create_bus(name="b1", type=bt, network=net)
        >>> _ = gi.create_branch(name="br", type=brt, from_bus="b0",
        ...                      to_bus="b1", length=1.0, network=net)
        >>> _ = gi.create_source(name="src", bus="b0",
        ...                      values={50: 100.0}, network=net)
        >>> # k_ref reproducing the original Z = 0.01*rho on b0/b1.
        >>> res = find_max_rho_f_scaling(
        ...     net, ["b0", "b1"], u_limit=200.0,
        ...     k_ref=(0.01, 0.0, 0.0, 0.0, 0.0),
        ... )
        >>> isinstance(res["c_max"], float)
        True
    """
    if u_limit <= 0:
        raise ValueError(f"u_limit must be positive, got {u_limit!r}.")
    if not bus_names:
        raise ValueError("bus_names must not be empty.")
    if len(k_ref) != 5:
        raise ValueError(f"k_ref must be a 5-tuple, got length {len(k_ref)}.")
    if all(x == 0 for x in k_ref):
        raise ValueError("k_ref must not be the zero vector.")
    c_lo_init, c_hi_init = c_bounds
    if not (0 < c_lo_init < c_hi_init):
        raise ValueError(
            f"c_bounds must satisfy 0 < c_lo < c_hi, got {c_bounds!r}."
        )

    last_epr_per_bus: Dict[str, Dict[str, float]] = {"value": {}}

    def _max_epr_at(c: float) -> float:
        eprs = evaluate_max_epr_under_k(
            network, bus_names,
            k=tuple(c * x for x in k_ref),
            fault_scalings=fault_scalings,
            run_fault_kwargs=run_fault_kwargs,
        )
        last_epr_per_bus["value"] = eprs
        return max(eprs.values()) if eprs else 0.0

    epr_lo = _max_epr_at(c_lo_init)
    eprs_lo_snapshot: Dict[str, float] = dict(last_epr_per_bus["value"])
    if epr_lo > u_limit:
        raise ValueError(
            f"u_limit={u_limit:g} V is below the maximum EPR at "
            f"c={c_lo_init:g}: max EPR_RMS={epr_lo:g} V — no scaling "
            f"factor in the bracket {c_bounds!r} satisfies the constraint."
        )

    epr_hi = _max_epr_at(c_hi_init)
    iterations = 0
    if epr_hi <= u_limit:
        logger.info(
            "Bracket fully admissible: max EPR_RMS at c=%g is %g V <= "
            "u_limit=%g V. Consider widening c_bounds.",
            c_hi_init, epr_hi, u_limit,
        )
        c_max, epr_at = c_hi_init, epr_hi
        epr_per_bus_at = dict(last_epr_per_bus["value"])
    else:
        c_lo, c_hi = c_lo_init, c_hi_init
        while iterations < max_iter and (c_hi - c_lo) / c_lo > tol_rel:
            c_mid = math.sqrt(c_lo * c_hi)  # geometric mean -> log bisection
            epr_mid = _max_epr_at(c_mid)
            if epr_mid <= u_limit:
                c_lo, epr_lo = c_mid, epr_mid
                eprs_lo_snapshot = dict(last_epr_per_bus["value"])
            else:
                c_hi = c_mid
            iterations += 1
        c_max, epr_at = c_lo, epr_lo
        epr_per_bus_at = eprs_lo_snapshot

    return {
        "c_max": c_max,
        "k_max": tuple(c_max * x for x in k_ref),
        "max_epr_rms_at_c_max": epr_at,
        "epr_rms_per_bus_at_c_max": epr_per_bus_at,
        "iterations": iterations,
    }

find_max_rho_scaling

find_max_rho_scaling(
    network: Network,
    fault_name: str,
    bus_names: List[str],
    u_max: float,
    *,
    c_bounds: Tuple[float, float] = (0.001, 1000.0),
    tol_rel: float = 0.001,
    max_iter: int = 60,
    run_fault_kwargs: Dict[str, Any] = None
) -> Dict[str, Any]

Find the largest uniform rho-scaling factor compatible with an EPR limit.

Bisects the scalar c such that scaling every selected bus' specific earth resistance to c * rho_0 yields an RMS earth potential rise at the fault bus that just satisfies |u_EPR|_rms <= u_max. The bus impedance formula is re-evaluated through the existing :meth:Bus.calculate_impedance machinery, so any user-defined rho-f characteristic is honoured.

Parameters:

Name	Type	Description	Default
network	Network	The simulation network. Must already contain the named fault, the sources, the buses listed in bus_names and consistent paths from sources to fault.	required
fault_name	str	Name of the fault to evaluate. Used as active_fault during every bisection step.	required
bus_names	List[str]	Names of the buses whose specific earth resistance is uniformly scaled by the same factor c. Must be non-empty and refer to buses in network.	required
u_max	float	Upper bound on the RMS earth potential rise at the fault bus, in volts. Must be strictly positive. The RMS is taken over all simulation frequencies, matching :attr:ResultBus.uepr.	required
c_bounds	Tuple[float, float]	Search interval for the scaling factor c. Both bounds must be strictly positive, and c_bounds[0] < c_bounds[1]. Defaults to (1e-3, 1e3), i.e. six decades.	(0.001, 1000.0)
tol_rel	float	Relative tolerance on the bracket width (c_hi - c_lo) / c_lo at which the bisection terminates. Defaults to 1e-3.	0.001
max_iter	int	Hard cap on the number of bisection steps. Defaults to 60 (largely sufficient given the log bracketing and tol_rel).	60
run_fault_kwargs	Dict[str, Any]	Extra keyword arguments forwarded to :func:run_fault at every step (e.g. {"auto_parallel_coefficients": True}). Defaults to None.	None

Returns:

Type

Description

Dict[str, Any]

Dict[str, Any]: A dictionary with keys "c_max" (float): The largest scaling factor satisfying the EPR constraint within the bracket. If epr(c_hi) <= u_max already, the upper bound is returned and iterations is zero (the bracket should be widened). "u_epr_rms_at_c_max" (float): The RMS EPR at the fault bus evaluated at c_max, in volts. "rho_max" (Dict[str, float]): c_max * rho_0[bus] for every selected bus. "iterations" (int): Number of bisection steps taken.

Raises:

Type	Description
ValueError	If u_max is not positive, bus_names is empty, c_bounds is invalid, any name is not in the network, or the EPR at the lower bound c_bounds[0] already exceeds u_max (constraint cannot be satisfied within the bracket).

Examples:

>>> import groundinsight as gi
>>> from groundinsight.models.core_models import BusType, BranchType
>>> from groundinsight.analysis import find_max_rho_scaling
>>> bt = BusType(name="BT", system_type="Grounded",
...              voltage_level=20.0,
...              impedance_formula="rho * 0.01 + I * f * 0")
>>> brt = BranchType(name="BRT", grounding_conductor=True,
...                  self_impedance_formula="(0.25 + I*0.6)*l",
...                  mutual_impedance_formula="(0.0 + I*0.6)*l")
>>> net = gi.create_network(name="N", frequencies=[50])
>>> _ = gi.create_bus(name="b0", type=bt, network=net)
>>> _ = gi.create_bus(name="b1", type=bt, network=net)
>>> _ = gi.create_branch(name="br", type=brt, from_bus="b0",
...                      to_bus="b1", length=1.0, network=net)
>>> _ = gi.create_source(name="src", bus="b0",
...                      values={50: 100.0}, network=net)
>>> _ = gi.create_fault(name="flt", bus="b1",
...                     scalings={50: 1.0}, network=net)
>>> res = find_max_rho_scaling(net, "flt", ["b0", "b1"],
...                            u_max=200.0)
>>> isinstance(res["c_max"], float)
True

Source code in src/groundinsight/analysis/inverse_rho.py

def find_max_rho_scaling(
    network: Network,
    fault_name: str,
    bus_names: List[str],
    u_max: float,
    *,
    c_bounds: Tuple[float, float] = (1e-3, 1e3),
    tol_rel: float = 1e-3,
    max_iter: int = 60,
    run_fault_kwargs: Dict[str, Any] = None,
) -> Dict[str, Any]:
    """
    Find the largest uniform rho-scaling factor compatible with an EPR limit.

    Bisects the scalar ``c`` such that scaling every selected bus' specific
    earth resistance to ``c * rho_0`` yields an RMS earth potential rise at
    the fault bus that just satisfies ``|u_EPR|_rms <= u_max``. The bus
    impedance formula is re-evaluated through the existing
    :meth:`Bus.calculate_impedance` machinery, so any user-defined rho-f
    characteristic is honoured.

    Args:
        network (Network): The simulation network. Must already contain the
            named fault, the sources, the buses listed in ``bus_names`` and
            consistent paths from sources to fault.
        fault_name (str): Name of the fault to evaluate. Used as
            ``active_fault`` during every bisection step.
        bus_names (List[str]): Names of the buses whose specific earth
            resistance is uniformly scaled by the same factor ``c``. Must
            be non-empty and refer to buses in ``network``.
        u_max (float): Upper bound on the RMS earth potential rise at the
            fault bus, in volts. Must be strictly positive. The RMS is
            taken over all simulation frequencies, matching
            :attr:`ResultBus.uepr`.
        c_bounds (Tuple[float, float], optional): Search interval for the
            scaling factor ``c``. Both bounds must be strictly positive,
            and ``c_bounds[0] < c_bounds[1]``. Defaults to
            ``(1e-3, 1e3)``, i.e. six decades.
        tol_rel (float, optional): Relative tolerance on the bracket width
            ``(c_hi - c_lo) / c_lo`` at which the bisection terminates.
            Defaults to ``1e-3``.
        max_iter (int, optional): Hard cap on the number of bisection
            steps. Defaults to ``60`` (largely sufficient given the log
            bracketing and ``tol_rel``).
        run_fault_kwargs (Dict[str, Any], optional): Extra keyword
            arguments forwarded to :func:`run_fault` at every step (e.g.
            ``{"auto_parallel_coefficients": True}``). Defaults to
            ``None``.

    Returns:
        Dict[str, Any]: A dictionary with keys

            - ``"c_max"`` (float): The largest scaling factor satisfying
              the EPR constraint within the bracket. If ``epr(c_hi) <=
              u_max`` already, the upper bound is returned and
              ``iterations`` is zero (the bracket should be widened).
            - ``"u_epr_rms_at_c_max"`` (float): The RMS EPR at the fault
              bus evaluated at ``c_max``, in volts.
            - ``"rho_max"`` (Dict[str, float]): ``c_max * rho_0[bus]``
              for every selected bus.
            - ``"iterations"`` (int): Number of bisection steps taken.

    Raises:
        ValueError: If ``u_max`` is not positive, ``bus_names`` is empty,
            ``c_bounds`` is invalid, any name is not in the network, or
            the EPR at the lower bound ``c_bounds[0]`` already exceeds
            ``u_max`` (constraint cannot be satisfied within the bracket).

    Examples:
        >>> import groundinsight as gi
        >>> from groundinsight.models.core_models import BusType, BranchType
        >>> from groundinsight.analysis import find_max_rho_scaling
        >>> bt = BusType(name="BT", system_type="Grounded",
        ...              voltage_level=20.0,
        ...              impedance_formula="rho * 0.01 + I * f * 0")
        >>> brt = BranchType(name="BRT", grounding_conductor=True,
        ...                  self_impedance_formula="(0.25 + I*0.6)*l",
        ...                  mutual_impedance_formula="(0.0 + I*0.6)*l")
        >>> net = gi.create_network(name="N", frequencies=[50])
        >>> _ = gi.create_bus(name="b0", type=bt, network=net)
        >>> _ = gi.create_bus(name="b1", type=bt, network=net)
        >>> _ = gi.create_branch(name="br", type=brt, from_bus="b0",
        ...                      to_bus="b1", length=1.0, network=net)
        >>> _ = gi.create_source(name="src", bus="b0",
        ...                      values={50: 100.0}, network=net)
        >>> _ = gi.create_fault(name="flt", bus="b1",
        ...                     scalings={50: 1.0}, network=net)
        >>> res = find_max_rho_scaling(net, "flt", ["b0", "b1"],
        ...                            u_max=200.0)
        >>> isinstance(res["c_max"], float)
        True
    """
    if u_max <= 0:
        raise ValueError(f"u_max must be positive, got {u_max!r}.")
    if not bus_names:
        raise ValueError("bus_names must not be empty.")
    c_lo_init, c_hi_init = c_bounds
    if not (0 < c_lo_init < c_hi_init):
        raise ValueError(
            f"c_bounds must satisfy 0 < c_lo < c_hi, got {c_bounds!r}."
        )
    missing = [b for b in bus_names if b not in network.buses]
    if missing:
        raise ValueError(f"Unknown bus(es) in network: {missing!r}.")
    if fault_name not in network.faults:
        raise ValueError(f"Unknown fault {fault_name!r} in network.")

    rfk: Dict[str, Any] = dict(run_fault_kwargs) if run_fault_kwargs else {}

    fault_bus_name = network.faults[fault_name].bus

    # Snapshot original rhos so the network can be restored at the end.
    rho_0: Dict[str, float] = {
        b: float(network.buses[b].specific_earth_resistance) for b in bus_names
    }

    def _epr_rms_at(c: float) -> float:
        """Evaluate the RMS EPR at the fault bus for a given scaling factor."""
        for b in bus_names:
            bus = network.buses[b]
            bus.specific_earth_resistance = c * rho_0[b]
            bus.calculate_impedance(network.frequencies)
        run_fault(network, fault_name=fault_name, **rfk)
        result_bus = next(
            rb
            for rb in network.results[fault_name].buses
            if rb.name == fault_bus_name
        )
        return float(result_bus.uepr)

    iterations = 0
    try:
        c_lo, c_hi = c_lo_init, c_hi_init
        epr_lo = _epr_rms_at(c_lo)
        if epr_lo > u_max:
            raise ValueError(
                f"u_max={u_max:g} V is below the EPR at c_min={c_lo:g}: "
                f"|u_EPR|_rms(c_min)={epr_lo:g} V — no scaling factor in "
                f"the bracket {c_bounds!r} satisfies the constraint."
            )
        epr_hi = _epr_rms_at(c_hi)
        if epr_hi <= u_max:
            # Whole bracket is admissible. Return c_hi but flag with zero
            # iterations so the caller can recognise the boundary case.
            logger.info(
                "Bracket fully admissible: |u_EPR|_rms(c_hi=%g)=%g V <= "
                "u_max=%g V. Consider widening c_bounds.",
                c_hi, epr_hi, u_max,
            )
            c_max, epr_at = c_hi, epr_hi
        else:
            while iterations < max_iter and (c_hi - c_lo) / c_lo > tol_rel:
                c_mid = math.sqrt(c_lo * c_hi)  # geometric mean -> log bisection
                epr_mid = _epr_rms_at(c_mid)
                if epr_mid <= u_max:
                    c_lo, epr_lo = c_mid, epr_mid
                else:
                    c_hi, epr_hi = c_mid, epr_mid
                iterations += 1
            c_max, epr_at = c_lo, epr_lo
    finally:
        # Restore original rhos and recompute their impedances no matter what.
        for b in bus_names:
            bus = network.buses[b]
            bus.specific_earth_resistance = rho_0[b]
            bus.calculate_impedance(network.frequencies)

    return {
        "c_max": c_max,
        "u_epr_rms_at_c_max": epr_at,
        "rho_max": {b: c_max * rho_0[b] for b in bus_names},
        "iterations": iterations,
    }

select_rho_f_from_catalog

select_rho_f_from_catalog(
    network: Network,
    bus_names: List[str],
    u_limit: float,
    candidates: Dict[str, KVector],
    *,
    fault_scalings: Optional[Dict[float, float]] = None,
    run_fault_kwargs: Optional[Dict[str, Any]] = None,
    sort_by: Literal[
        "max_epr_asc", "max_epr_desc", "name"
    ] = "max_epr_asc"
) -> pl.DataFrame

Pick admissible rho-f characteristics from a user-provided catalog.

Evaluates every candidate k in the catalog with :func:evaluate_max_epr_under_k and reports, per candidate, the maximum RMS EPR observed across the bus sweep, the per-bus EPRs and whether the candidate satisfies the limit u_limit. The catalog is typically a hand-curated list of soil scenarios (e.g. dry sand, wet clay, permafrost) or rho-f fits from previous measurements.

Parameters:

Name	Type	Description	Default
network	Network	The simulation network. Sources, branches and the buses to be swept must already be configured.	required
bus_names	List[str]	Buses whose impedance is parameterised by every candidate k and which are swept as fault locations.	required
u_limit	float	Upper bound on the RMS EPR (in volts) at any swept bus. Must be strictly positive.	required
candidates	Dict[str, KVector]	Mapping {name: (k1, k2, k3, k4, k5)} of the candidate rho-f characteristics. Names must be unique (Python dict guarantees that), tuples must each have length 5. May be empty -- the result is then an empty DataFrame with the documented schema.	required
fault_scalings	Optional[Dict[float, float]]	Frequency-resolved scalings for any fault that has to be created on the fly. See :func:evaluate_max_epr_under_k.	None
run_fault_kwargs	Optional[Dict[str, Any]]	Forwarded to :func:run_fault.	None
sort_by	Literal['max_epr_asc', 'max_epr_desc', 'name']	How to sort the result rows. "max_epr_asc" (default): admissible candidates first, tightest-EPR-margin candidates at the top. "max_epr_desc": largest EPR first; useful to inspect the worst-case candidates. "name": lexicographic by candidate name.	'max_epr_asc'

Returns:

Type	Description
DataFrame	A Polars DataFrame with one row per candidate and the columns
DataFrame	"name" (str)
DataFrame	"k1", "k2", "k3", "k4", "k5" (float)
DataFrame	"max_epr_rms_V" (float): maximum RMS EPR across the bus sweep, in volts.
DataFrame	"admissible" (bool): True iff max_epr_rms_V <= u_limit.
DataFrame	one "epr_<bus>_V" column per swept bus, holding the RMS EPR at that bus when it is the active fault, in volts.

Raises:

Type	Description
ValueError	If u_limit is non-positive, bus_names is empty, or any candidate has a wrong-length k or refers to an unknown bus (the underlying helper validates this).

Examples:

>>> import groundinsight as gi
>>> from groundinsight.models.core_models import BusType, BranchType
>>> from groundinsight.analysis import select_rho_f_from_catalog
>>> bt = BusType(name="BT", system_type="Grounded",
...              voltage_level=20.0,
...              impedance_formula="rho * 0.01 + 0*f")
>>> brt = BranchType(name="BRT", grounding_conductor=True,
...                  self_impedance_formula="(0.25 + I*0.6)*l",
...                  mutual_impedance_formula="(0.0 + I*0.6)*l")
>>> net = gi.create_network(name="N", frequencies=[50])
>>> _ = gi.create_bus(name="b0", type=bt, network=net)
>>> _ = gi.create_bus(name="b1", type=bt, network=net)
>>> _ = gi.create_branch(name="br", type=brt, from_bus="b0",
...                      to_bus="b1", length=1.0, network=net)
>>> _ = gi.create_source(name="src", bus="b0",
...                      values={50: 100.0}, network=net)
>>> catalog = {
...     "low":   (0.005, 0.0, 0.0, 0.0, 0.0),
...     "med":   (0.01,  0.0, 0.0, 0.0, 0.0),
...     "high":  (0.05,  0.0, 0.0, 0.0, 0.0),
... }
>>> df = select_rho_f_from_catalog(
...     net, ["b0", "b1"], u_limit=20.0, candidates=catalog,
... )
>>> set(df.columns) >= {"name", "k1", "max_epr_rms_V", "admissible"}
True

Source code in src/groundinsight/analysis/inverse_rho_f.py

def select_rho_f_from_catalog(
    network: Network,
    bus_names: List[str],
    u_limit: float,
    candidates: Dict[str, KVector],
    *,
    fault_scalings: Optional[Dict[float, float]] = None,
    run_fault_kwargs: Optional[Dict[str, Any]] = None,
    sort_by: Literal["max_epr_asc", "max_epr_desc", "name"] = "max_epr_asc",
) -> pl.DataFrame:
    """Pick admissible rho-f characteristics from a user-provided catalog.

    Evaluates every candidate ``k`` in the catalog with
    :func:`evaluate_max_epr_under_k` and reports, per candidate, the
    maximum RMS EPR observed across the bus sweep, the per-bus EPRs and
    whether the candidate satisfies the limit ``u_limit``. The catalog is
    typically a hand-curated list of soil scenarios (e.g. dry sand, wet
    clay, permafrost) or rho-f fits from previous measurements.

    Args:
        network: The simulation network. Sources, branches and the buses
            to be swept must already be configured.
        bus_names: Buses whose impedance is parameterised by every
            candidate ``k`` and which are swept as fault locations.
        u_limit: Upper bound on the RMS EPR (in volts) at any swept bus.
            Must be strictly positive.
        candidates: Mapping ``{name: (k1, k2, k3, k4, k5)}`` of the
            candidate rho-f characteristics. Names must be unique
            (Python dict guarantees that), tuples must each have length
            5. May be empty -- the result is then an empty DataFrame
            with the documented schema.
        fault_scalings: Frequency-resolved scalings for any fault that
            has to be created on the fly. See
            :func:`evaluate_max_epr_under_k`.
        run_fault_kwargs: Forwarded to :func:`run_fault`.
        sort_by: How to sort the result rows.

            - ``"max_epr_asc"`` (default): admissible candidates first,
              tightest-EPR-margin candidates at the top.
            - ``"max_epr_desc"``: largest EPR first; useful to inspect
              the worst-case candidates.
            - ``"name"``: lexicographic by candidate name.

    Returns:
        A Polars DataFrame with one row per candidate and the columns

        - ``"name"`` (str)
        - ``"k1", "k2", "k3", "k4", "k5"`` (float)
        - ``"max_epr_rms_V"`` (float): maximum RMS EPR across the bus
          sweep, in volts.
        - ``"admissible"`` (bool): ``True`` iff
          ``max_epr_rms_V <= u_limit``.
        - one ``"epr_<bus>_V"`` column per swept bus, holding the RMS
          EPR at that bus when it is the active fault, in volts.

    Raises:
        ValueError: If ``u_limit`` is non-positive, ``bus_names`` is
            empty, or any candidate has a wrong-length ``k`` or refers
            to an unknown bus (the underlying helper validates this).

    Examples:
        >>> import groundinsight as gi
        >>> from groundinsight.models.core_models import BusType, BranchType
        >>> from groundinsight.analysis import select_rho_f_from_catalog
        >>> bt = BusType(name="BT", system_type="Grounded",
        ...              voltage_level=20.0,
        ...              impedance_formula="rho * 0.01 + 0*f")
        >>> brt = BranchType(name="BRT", grounding_conductor=True,
        ...                  self_impedance_formula="(0.25 + I*0.6)*l",
        ...                  mutual_impedance_formula="(0.0 + I*0.6)*l")
        >>> net = gi.create_network(name="N", frequencies=[50])
        >>> _ = gi.create_bus(name="b0", type=bt, network=net)
        >>> _ = gi.create_bus(name="b1", type=bt, network=net)
        >>> _ = gi.create_branch(name="br", type=brt, from_bus="b0",
        ...                      to_bus="b1", length=1.0, network=net)
        >>> _ = gi.create_source(name="src", bus="b0",
        ...                      values={50: 100.0}, network=net)
        >>> catalog = {
        ...     "low":   (0.005, 0.0, 0.0, 0.0, 0.0),
        ...     "med":   (0.01,  0.0, 0.0, 0.0, 0.0),
        ...     "high":  (0.05,  0.0, 0.0, 0.0, 0.0),
        ... }
        >>> df = select_rho_f_from_catalog(
        ...     net, ["b0", "b1"], u_limit=20.0, candidates=catalog,
        ... )
        >>> set(df.columns) >= {"name", "k1", "max_epr_rms_V", "admissible"}
        True
    """
    if u_limit <= 0:
        raise ValueError(f"u_limit must be positive, got {u_limit!r}.")
    if not bus_names:
        raise ValueError("bus_names must not be empty.")

    # Empty catalog -> return an empty DataFrame with the documented schema.
    if not candidates:
        schema = {
            "name": pl.Utf8,
            "k1": pl.Float64, "k2": pl.Float64, "k3": pl.Float64,
            "k4": pl.Float64, "k5": pl.Float64,
            "max_epr_rms_V": pl.Float64,
            "admissible": pl.Boolean,
        }
        for b in bus_names:
            schema[f"epr_{b}_V"] = pl.Float64
        return pl.DataFrame(schema=schema)

    rows: List[Dict[str, Any]] = []
    for name, k in candidates.items():
        if len(k) != 5:
            raise ValueError(
                f"Candidate {name!r}: k must be a 5-tuple, got length {len(k)}."
            )
        eprs = evaluate_max_epr_under_k(
            network, bus_names, k=tuple(k),
            fault_scalings=fault_scalings,
            run_fault_kwargs=run_fault_kwargs,
        )
        max_epr = max(eprs.values()) if eprs else 0.0
        row: Dict[str, Any] = {
            "name": name,
            "k1": float(k[0]),
            "k2": float(k[1]),
            "k3": float(k[2]),
            "k4": float(k[3]),
            "k5": float(k[4]),
            "max_epr_rms_V": float(max_epr),
            "admissible": bool(max_epr <= u_limit),
        }
        for b in bus_names:
            row[f"epr_{b}_V"] = float(eprs[b])
        rows.append(row)

    df = pl.DataFrame(rows)

    if sort_by == "max_epr_asc":
        # Admissible first (True -> 1 sorts after False -> 0 by default,
        # so descending on ``admissible`` puts True on top), then
        # ascending EPR within each block.
        df = df.sort(
            by=["admissible", "max_epr_rms_V"], descending=[True, False]
        )
    elif sort_by == "max_epr_desc":
        df = df.sort(by="max_epr_rms_V", descending=True)
    elif sort_by == "name":
        df = df.sort(by="name")
    else:  # pragma: no cover -- guarded by Literal type, defensive
        raise ValueError(
            f"sort_by must be one of "
            f"'max_epr_asc' | 'max_epr_desc' | 'name', got {sort_by!r}."
        )

    return df