Parameter sweeps

Solve one fault once per parameter combination and stack the results into long-format frames that carry the parameters as columns. This is what the statistics below operate on — until such a frame exists there is nothing to summarise, because every accessor on Network reports a single solve.

rho_f_points builds the points from a catalogue of five-parameter rho-f vectors, which is the form groundfield exports.

sweep

Run one fault over a grid of parameter variations and collect it into one frame.

run_outage_study already answers "what changes when an element drops out". This module answers the other half: what changes when a parameter moves -- the rho-f characteristic of the faulted station, the soil resistivity, the harmonic content of the source, the fault location itself.

The unit of work is a :class:SweepPoint: a label plus the overrides that define it. Every point is applied, solved and restored, and the results are stacked into long-format Polars frames that carry the point's label and its parameters as ordinary columns. That last part is the point of the module -- until a frame like this exists there is nothing for a statistic to operate on, which is why :mod:groundinsight.analysis.statistics starts here.

Overrides are applied by writing directly onto Bus.impedance and friends, which reaches the solver untouched because run_fault does not recompute impedances (see :class:~groundinsight.electrical_network.ElectricalNetwork). Everything is restored in a finally block, so an exception in the middle of a sweep leaves the network exactly as it was found -- including the path cache and the active fault.

Example

Vary the rho-f characteristic of the faulted station over a grid and watch what it does to the potential rise and to the parallel impedances of the two feeder directions::

points = gi.rho_f_points(
    bus="Station_7",
    k_vectors={f"k1={k1:g}": (k1, 1e-4, 3e-4, 0.0, 0.0)
               for k1 in (0.01, 0.02, 0.05, 0.1)},
)
study = gi.run_sweep(
    net, fault="F_Station_7", points=points,
    cuts=[gi.Cut(name="left", branches=["C6_7"]),
          gi.Cut(name="right", branches=["C7_8"])],
)
study.impedances()   # Z_G and both reduction factors per point
study.cuts()         # Z_left / Z_right / r_left / r_right per point

SweepPoint

Bases: BaseModel

One parameter combination to solve.

Attributes:

Name	Type	Description
label	str	Identifies the point in every result frame. Must be unique in a sweep.
bus_impedance	dict of str to dict of float to complex, optional	Impedance tables written straight onto Bus.impedance. Bypasses the bus type's formula entirely -- this is how a rho-f characteristic measured or fitted elsewhere enters the study.
bus_rho	dict of str to float, optional	New specific_earth_resistance per bus, with the bus type's formula re-evaluated afterwards. Use this to vary the soil rather than the characteristic.
fault	(str, optional)	Solve a different fault at this point. Defaults to the sweep's fault.
fault_scalings	dict of float to float, optional	Replace the active fault's per-frequency scalings -- the harmonic content of the source.
parameters	dict of str to object, optional	Free-form values copied into every result row as columns, so a plot can be made against the physical parameter rather than against the label.

SweepResult

Bases: BaseModel

Everything a sweep collected, as long-format frames.

Attributes:

Name	Type	Description
fault	str	The fault the sweep was run for (points may override it individually).
labels	list of str	Point labels in the order they were solved.
failures	dict of str to str	Points that raised, mapped to the exception text. They are absent from the frames; the sweep does not abort on one bad point.

branches

branches() -> pl.DataFrame

Per-branch results of every point, stacked.

Source code in src/groundinsight/simulation/sweep.py

def branches(self) -> pl.DataFrame:
    """Per-branch results of every point, stacked."""
    return _or_empty(self._branches)

buses

buses() -> pl.DataFrame

Per-bus results of every point, stacked.

Source code in src/groundinsight/simulation/sweep.py

def buses(self) -> pl.DataFrame:
    """Per-bus results of every point, stacked."""
    return _or_empty(self._buses)

cuts

cuts() -> pl.DataFrame

Side impedances and side reduction factors, stacked. Empty when the sweep was run without cuts.

Source code in src/groundinsight/simulation/sweep.py

def cuts(self) -> pl.DataFrame:
    """Side impedances and side reduction factors, stacked. Empty when the
    sweep was run without cuts."""
    return _or_empty(self._cuts)

impedances

impedances() -> pl.DataFrame

Z_G and both reduction factors of every point, stacked.

Source code in src/groundinsight/simulation/sweep.py

def impedances(self) -> pl.DataFrame:
    """``Z_G`` and both reduction factors of every point, stacked."""
    return _or_empty(self._impedances)

rho_f_points

rho_f_points(
    *,
    bus: str,
    k_vectors: Dict[str, KVector],
    frequencies: Sequence[float],
    rho: float
) -> List[SweepPoint]

Build sweep points from a catalogue of rho-f parameter vectors.

Each vector is evaluated into an impedance table for bus over frequencies at the given rho, so the study varies the fitted characteristic of one station while the rest of the network stays put. k1 ... k5 and rho are copied into parameters, which puts them in the result frames as plottable columns.

Parameters:

Name	Type	Description	Default
bus	str	Name of the bus whose characteristic is varied.	required
k_vectors	dict of str to tuple	Label to (k1, k2, k3, k4, k5).	required
frequencies	sequence of float	Frequencies to evaluate the form at -- normally network.frequencies.	required
rho	float	Soil resistivity the form is evaluated at.	required

Returns:

Type	Description
list of SweepPoint

Raises:

Type	Description
ValueError	If a vector produces a non-positive real part at any frequency, which the solver rejects as non-passive. The offending label and frequency are named, because an unconstrained least-squares fit can land there and the failure is otherwise reported far from its cause.

Source code in src/groundinsight/simulation/sweep.py

def rho_f_points(
    *,
    bus: str,
    k_vectors: Dict[str, KVector],
    frequencies: Sequence[float],
    rho: float,
) -> List[SweepPoint]:
    """
    Build sweep points from a catalogue of rho-f parameter vectors.

    Each vector is evaluated into an impedance table for ``bus`` over
    ``frequencies`` at the given ``rho``, so the study varies the *fitted
    characteristic* of one station while the rest of the network stays put.
    ``k1`` ... ``k5`` and ``rho`` are copied into ``parameters``, which puts
    them in the result frames as plottable columns.

    Parameters
    ----------
    bus : str
        Name of the bus whose characteristic is varied.
    k_vectors : dict of str to tuple
        Label to ``(k1, k2, k3, k4, k5)``.
    frequencies : sequence of float
        Frequencies to evaluate the form at -- normally ``network.frequencies``.
    rho : float
        Soil resistivity the form is evaluated at.

    Returns
    -------
    list of SweepPoint

    Raises
    ------
    ValueError
        If a vector produces a non-positive real part at any frequency, which
        the solver rejects as non-passive. The offending label and frequency are
        named, because an unconstrained least-squares fit can land there and the
        failure is otherwise reported far from its cause.
    """
    points: List[SweepPoint] = []
    for label, k in k_vectors.items():
        if len(k) != 5:
            raise ValueError(
                f"rho-f vector '{label}' has {len(k)} entries; the standard "
                f"form takes exactly five, (k1, k2, k3, k4, k5)."
            )
        table: Dict[float, complex] = {}
        for freq in frequencies:
            z = _z_rho_f(tuple(k), rho, float(freq))
            if z.real <= 0.0:
                raise ValueError(
                    f"rho-f vector '{label}' gives Re(Z) = {z.real:.6g} Ohm at "
                    f"{freq} Hz with rho = {rho} Ohm*m, which is not a passive "
                    f"impedance and the solver will reject it. An unconstrained "
                    f"least-squares fit can produce this below the frequency "
                    f"range it was fitted on -- check the fit's validity range."
                )
            table[float(freq)] = z
        points.append(
            SweepPoint(
                label=label,
                bus_impedance={bus: table},
                parameters={
                    "bus": bus,
                    "rho_Ohm_m": rho,
                    "k1": k[0],
                    "k2": k[1],
                    "k3": k[2],
                    "k4": k[3],
                    "k5": k[4],
                },
            )
        )
    return points

run_sweep

run_sweep(
    network: Network,
    *,
    fault: str,
    points: Sequence[SweepPoint],
    cuts: Optional[Sequence[Cut]] = None,
    phase_current_mode: str = "auto",
    collect_branches: bool = False,
    on_error: str = "record"
) -> SweepResult

Solve one fault once per parameter point and stack the results.

Parameters:

Name	Type	Description	Default
network	Network	The network. Left exactly as found, whatever happens.	required
fault	str	Fault to solve, unless a point names its own.	required
points	sequence of SweepPoint	The parameter grid. Labels must be unique.	required
cuts	sequence of Cut	When given, :func:~groundinsight.analysis.analyze_cuts runs at every point and its frame is stacked into SweepResult.cuts().	None
phase_current_mode	('auto', 'paths')	Forwarded to :func:run_fault.	"auto"
collect_branches	bool	Branch results multiply the row count by the number of branches and are rarely what a parameter study plots, so they are off by default.	False
on_error	('record', 'raise')	"record" (default) notes a failing point in SweepResult.failures and carries on -- a single non-passive parameter combination should not throw away the rest of a long grid. "raise" propagates instead.	"record"

Returns:

Type	Description
SweepResult

Raises:

Type	Description
ValueError	If points is empty, if two points share a label, if the fault is unknown, or if on_error is not one of the two accepted values.

Source code in src/groundinsight/simulation/sweep.py

def run_sweep(
    network: Network,
    *,
    fault: str,
    points: Sequence[SweepPoint],
    cuts: Optional[Sequence[Cut]] = None,
    phase_current_mode: str = "auto",
    collect_branches: bool = False,
    on_error: str = "record",
) -> SweepResult:
    """
    Solve one fault once per parameter point and stack the results.

    Parameters
    ----------
    network : Network
        The network. Left exactly as found, whatever happens.
    fault : str
        Fault to solve, unless a point names its own.
    points : sequence of SweepPoint
        The parameter grid. Labels must be unique.
    cuts : sequence of Cut, optional
        When given, :func:`~groundinsight.analysis.analyze_cuts` runs at every
        point and its frame is stacked into ``SweepResult.cuts()``.
    phase_current_mode : {"auto", "paths"}, optional
        Forwarded to :func:`run_fault`.
    collect_branches : bool, optional
        Branch results multiply the row count by the number of branches and are
        rarely what a parameter study plots, so they are off by default.
    on_error : {"record", "raise"}, optional
        ``"record"`` (default) notes a failing point in ``SweepResult.failures``
        and carries on -- a single non-passive parameter combination should not
        throw away the rest of a long grid. ``"raise"`` propagates instead.

    Returns
    -------
    SweepResult

    Raises
    ------
    ValueError
        If ``points`` is empty, if two points share a label, if the fault is
        unknown, or if ``on_error`` is not one of the two accepted values.
    """
    if on_error not in ("record", "raise"):
        raise ValueError(
            f"on_error must be 'record' or 'raise', got {on_error!r}."
        )
    if not points:
        raise ValueError("A sweep needs at least one point.")
    labels = [p.label for p in points]
    duplicates = {label for label in labels if labels.count(label) > 1}
    if duplicates:
        raise ValueError(
            f"Sweep point label(s) {sorted(duplicates)} appear more than once. "
            f"Labels identify the rows of every result frame and have to be "
            f"unique."
        )
    if fault not in network.faults:
        raise ValueError(
            f"Fault '{fault}' does not exist in network '{network.name}'. "
            f"Available: {sorted(network.faults)}."
        )

    previous_active = network.active_fault
    bus_frames: List[pl.DataFrame] = []
    branch_frames: List[pl.DataFrame] = []
    impedance_frames: List[pl.DataFrame] = []
    cut_frames: List[pl.DataFrame] = []
    solved: List[str] = []
    failures: Dict[str, str] = {}

    for point in points:
        target = point.fault or fault
        if target not in network.faults:
            message = (
                f"point '{point.label}' names fault '{target}', which does not "
                f"exist in network '{network.name}'"
            )
            if on_error == "raise":
                raise ValueError(message)
            failures[point.label] = message
            continue

        saved = None
        try:
            saved = _apply(network, point)
            if point.fault_scalings is not None:
                fault_obj = network.faults[target]
                saved["scalings_fault"] = target
                saved["scalings"] = dict(fault_obj.scalings)
                fault_obj.scalings = dict(point.fault_scalings)

            run_fault(network, target, phase_current_mode=phase_current_mode)

            bus_frames.append(_tag(network.res_buses(fault=target), point))
            if collect_branches:
                branch_frames.append(
                    _tag(network.res_branches(fault=target), point)
                )
            impedance_frames.append(
                _tag(
                    network.res_all_impedances().filter(
                        pl.col("fault_name") == target
                    ),
                    point,
                )
            )
            if cuts:
                analysis = analyze_cuts(network, fault=target, cuts=cuts)
                cut_frames.append(_tag(analysis.to_polars(), point))
            solved.append(point.label)
        except Exception as exc:  # noqa: BLE001 -- recorded, then re-raised or not
            if on_error == "raise":
                raise
            logger.warning(
                "Sweep point '%s' failed and was skipped: %s", point.label, exc
            )
            failures[point.label] = str(exc)
        finally:
            if saved is not None:
                _restore(network, saved)

    if previous_active is not None and previous_active in network.faults:
        network.set_active_fault(previous_active)

    result = SweepResult(fault=fault, labels=solved, failures=failures)
    result._buses = pl.concat(bus_frames, how="diagonal") if bus_frames else None
    result._branches = (
        pl.concat(branch_frames, how="diagonal") if branch_frames else None
    )
    result._impedances = (
        pl.concat(impedance_frames, how="diagonal") if impedance_frames else None
    )
    result._cuts = pl.concat(cut_frames, how="diagonal") if cut_frames else None
    if failures:
        logger.warning(
            "%d of %d sweep points failed: %s",
            len(failures),
            len(points),
            ", ".join(sorted(failures)),
        )
    return result